{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "brsPGTNRhIQD"
   },
   "source": [
    "# Problem 3\n",
    "The infinite matrix $ A $ has entries $ a_{11} = 1, a_{12} = \\frac{1}{2}, a_{21} = \\frac{1}{3}, a_{13} = \\frac{1}{4}, a_{22} = \\frac{1}{5}, a_{31} = \\frac{1}{6} $ ... is a bounded operator on $ \\ell_2$. What is $\\Vert{A}\\Vert$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "l1Jzx5S4hIQG"
   },
   "source": [
    "### Finding a general formula for the matrix's terms\n",
    "We have that\n",
    "\n",
    "$ A = \\begin{bmatrix} \n",
    "1   & 1/2 & 1/4 & 1/7 & \\cdots \\\\\n",
    "1/3 & 1/5 & 1/8 & \\cdots & \\vdots \\\\\n",
    "1/6 & 1/9 & \\cdots & \\ddots & \\vdots \\\\\n",
    "1/10 & \\cdots & \\cdots & \\ddots & \\vdots \\\\\n",
    "\\vdots & \\vdots & \\vdots & \\vdots &\\vdots & \\\\\n",
    "\\end{bmatrix} $\n",
    "\n",
    "Consider the matrix with all of the reciprocals of the elements in A:\n",
    "\n",
    "\n",
    "$ B = \\begin{bmatrix}\n",
    "1 & 2 & 4 & 7 & 11 & 16 & \\cdots \\\\\n",
    "3 & 5 & 8 & 12 & 17 & \\cdots \\\\\n",
    "6 & 9 & 13 & 18 & \\cdots \\\\\n",
    "10 & 14 & 19 & \\cdots \\\\\n",
    "15 & 20 \\cdots \\\\\n",
    "21 \\cdots \\\\\n",
    "\\end {bmatrix} $\n",
    "\n",
    "Noting that the first column of $B$ is just the triangular numbers, we guess that the general term $B_{ij}$ can be represented as a quadratic polynomial in both i and j, that is $ B_{ij} = a i^2 + b ij + c j^2 + d i + e j + f $ with $i,j\\geq0$.\n",
    "\n",
    "We can pick six elements from $B$ and solve the resulting linear system of equations to get that\n",
    "$$ B_{ij} = \\frac{1}{2} i^2 + ij + \\frac{1}{2}j^2 + \\frac{3}{2} i + \\frac{1}{2} j + 1, \\; \\; A_{ij} = \\frac{1}{B_{ij}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EghWW2YFhIQH"
   },
   "source": [
    "### Brute force\n",
    "We first find the spectral norm of the $n \\times n$  upper left quadrant of $A$ for $n=1,2,4,...,4096$ and see if we can get converging answers to get some accurate digits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 672
    },
    "executionInfo": {
     "elapsed": 52249,
     "status": "ok",
     "timestamp": 1602316767832,
     "user": {
      "displayName": "Husnain Raza",
      "photoUrl": "",
      "userId": "16146658556003895350"
     },
     "user_tz": 300
    },
    "id": "oamnOnUXhIQJ",
    "outputId": "4096c28f-5786-4513-c21e-9f1e6aa60fb9"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 1.0\n",
      "2 1.1833501765516565\n",
      "4 1.2525373975167995\n",
      "8 1.2700463058540765\n",
      "16 1.2735252154501302\n",
      "32 1.2741181443691536\n",
      "64 1.274209131297663\n",
      "128 1.2742221200377764\n",
      "256 1.2742238859485546\n",
      "512 1.2742241184580805\n",
      "1024 1.274224148449696\n",
      "2048 1.2742241522691704\n"
     ]
    }
   ],
   "source": [
    "from numpy.linalg import norm\n",
    "A = lambda i,j: 1/(0.5*i*i + i*j + 0.5*j*j + 1.5*i + 0.5*j + 1)\n",
    "for n in [2**i for i in range(12)]:\n",
    "    print(n, norm([[A(i,j) for i in range(n)] for j in range(n)],2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "hY-BUgffhIQV"
   },
   "source": [
    "From the above, we note that:\n",
    " - the limit $ \\lim_{N \\to \\infty} \\Vert [A_{ij}]_{0\\le i,j \\le N} \\Vert$ seems to converge\n",
    " - it seems that the norm of A is around `1.2742241`\n",
    "\n",
    "We can confirm these digits by calculating the spectral norm of a larger block of the matrix directly"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-CAHa2j0hIQd"
   },
   "source": [
    "Recall: $ \\Vert{A}\\Vert_2 = \\sqrt{\\lambda_{\\max}(A^T A)} $\n",
    "where $\\lambda_{\\max}$ is the largest eigenvalue of $A$\n",
    "\n",
    "Calculating this on the upper-left $4096\\times4096$ quadrant of $A$ gives us the answer to 10 digits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "id": "fOuHjD9IhIQi"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 7.29 s, sys: 2.09 s, total: 9.38 s\n",
      "Wall time: 4.76 s\n",
      "CPU times: user 2min 39s, sys: 43.9 s, total: 3min 23s\n",
      "Wall time: 1min 23s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(1.2742241527518152+0j)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy\n",
    "N = 4096\n",
    "xx, yy = np.meshgrid(range(N), range(N))\n",
    "M = np.vectorize(A)(yy, xx)\n",
    "%time N = M @ M.T\n",
    "%time max_eigval = max(scipy.linalg.eigvals(N))\n",
    "max_eigval**0.5"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "Problem 3.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}