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$?
We have that
$ A = \begin{bmatrix} 1 & 1/2 & 1/4 & 1/7 & \cdots \\ 1/3 & 1/5 & 1/8 & \cdots & \vdots \\ 1/6 & 1/9 & \cdots & \ddots & \vdots \\ 1/10 & \cdots & \cdots & \ddots & \vdots \\ \vdots & \vdots & \vdots & \vdots &\vdots & \\ \end{bmatrix} $
Consider the matrix with all of the reciprocals of the elements in A:
$ B = \begin{bmatrix} 1 & 2 & 4 & 7 & 11 & 16 & \cdots \\ 3 & 5 & 8 & 12 & 17 & \cdots \\ 6 & 9 & 13 & 18 & \cdots \\ 10 & 14 & 19 & \cdots \\ 15 & 20 \cdots \\ 21 \cdots \\ \end {bmatrix} $
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$.
We can pick six elements from $B$ and solve the resulting linear system of equations to get that $$ 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}}$$
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
from numpy.linalg import norm
A = lambda i,j: 1/(0.5*i*i + i*j + 0.5*j*j + 1.5*i + 0.5*j + 1)
for n in [2**i for i in range(12)]:
print(n, norm([[A(i,j) for i in range(n)] for j in range(n)],2))
1 1.0 2 1.1833501765516565 4 1.2525373975167995 8 1.2700463058540765 16 1.2735252154501302 32 1.2741181443691536 64 1.274209131297663 128 1.2742221200377764 256 1.2742238859485546 512 1.2742241184580805 1024 1.274224148449696 2048 1.2742241522691704
From the above, we note that:
1.2742241We can confirm these digits by calculating the spectral norm of a larger block of the matrix directly
Recall: $ \Vert{A}\Vert_2 = \sqrt{\lambda_{\max}(A^T A)} $ where $\lambda_{\max}$ is the largest eigenvalue of $A$
Calculating this on the upper-left $4096\times4096$ quadrant of $A$ gives us the answer to 10 digits:
import numpy as np
import scipy
N = 4096
xx, yy = np.meshgrid(range(N), range(N))
M = np.vectorize(A)(yy, xx)
%time N = M @ M.T
%time max_eigval = max(scipy.linalg.eigvals(N))
max_eigval**0.5
CPU times: user 7.29 s, sys: 2.09 s, total: 9.38 s Wall time: 4.76 s CPU times: user 2min 39s, sys: 43.9 s, total: 3min 23s Wall time: 1min 23s
(1.2742241527518152+0j)