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

# # Interpolative Decomposition

# In[1]:


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


# ## Obtain a low-rank matrix

# In[2]:


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

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

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


# ## Run the factorization

# In[3]:


import scipy.linalg.interpolative as sli


# Compute a fixed-rank factorization:
# 
# (There's also an adaptive, fixed-precision mode.)

# In[4]:


k = 20
idx, proj = sli.interp_decomp(A, k)


# Examine `idx`:

# In[5]:


idx


# What does `numpy.argsort` do?

# In[6]:


sort_idx = np.argsort(idx)


# In[7]:


idx[sort_idx]


# Reconstruct the matrix:

# In[8]:


B = A[:,idx[:k]]
P = np.hstack([np.eye(k), proj])[:,np.argsort(idx)]
Aapprox = B@P


# In[9]:


la.norm(A - Aapprox, 2)


# What's the structure of $P$?
# 
# (ignoring the column permuation)

# In[10]:


pt.imshow(np.hstack([np.eye(k), proj]))


# In[ ]: