# coding: utf-8

# # Rank-Revealing QR

# **Note:** `scipy.linalg`, not `numpy.linalg`!

# In[6]:

import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as pt


# ## Obtain a low-rank matrix

# In[10]:

n = 100
A0 = np.random.randn(n, n)
U0, sigma0, VT0 = la.svd(A0)
print(la.norm((U0*sigma0).dot(VT0) - A0))

sigma = np.exp(-np.arange(n))

A = (U0 * sigma).dot(VT0)


# ## Run the factorization

# Compute the QR factorization with `pivoting=True`.

# In[26]:

Q, R, perm = la.qr(A, pivoting=True)


# First of all, check that we've obtained a valid factorization

# In[27]:

la.norm(A[:, perm] - Q.dot(R), 2)


# In[28]:

la.norm(Q.dot(Q.T) - np.eye(n))


# Next, examine $R$:

# In[29]:

pt.imshow(np.log10(1e-15+np.abs(R)))
pt.colorbar()


# Specifically, recall that the diagonal of $R$ in QR contains column norms:

# In[30]:

pt.semilogy(np.abs(np.diag(R)))


# * In the case of `scipy`'s transform, diagonal entries of $R$ are guaranteed non-increasing.
# * But there is a whole science to how to choose the permutations (or other source vectors)
#     * and what promises one is able to make as a result of that