#!/usr/bin/env python
# coding: utf-8

# # Problem 3
# 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$?

# ### Finding a general formula for the matrix's terms
# 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}}$$

# ### Brute force
# 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

# In[26]:


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))


# From the above, we note that:
#  - the limit $ \lim_{N \to \infty} \Vert [A_{ij}]_{0\le i,j \le N} \Vert$ seems to converge
#  - it seems that the norm of A is around `1.2742241`
# 
# We 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:

# In[24]:


import numpy as np
import scipy
N = 4096
xx, yy = np.meshgrid(range(N), range(N))
M = np.vectorize(A)(yy, xx)
get_ipython().run_line_magic('time', 'N = M @ M.T')
get_ipython().run_line_magic('time', 'max_eigval = max(scipy.linalg.eigvals(N))')
max_eigval**0.5