Matrix completion using trace norm regularization
This shows the use of Douglas-Rachford's algorithm to perform matrix completion with trace norm (nuclear norm) regularization.
import numpy as np
import matplotlib.pyplot as plt
We consider $n \times n$ matrices.
n = 100
Generate a rank $r$ random matrix.
r = 10
x0 = np.dot(np.random.randn(n,r),np.random.randn(r,n))
Display the singular values.
plt.plot(np.linalg.svd(x0)[1], '.-')
[<matplotlib.lines.Line2D at 0x7f2a0b413d90>]
Random selection of $p$ indices.
def indices(n,p): return np.random.permutation(n**2)[0:p]
Associated forward (measurement) operator $\Phi$ and its adjoint $\Phi^\top$.
def phi(x,I): return x.flatten()[I]
def phi_adj(y,I):
x = np.zeros(n**2)
x[I] = y
return x.reshape((n,n))
def dotp(x,y): return np.sum( x.flatten()*y.flatten() )
Check adjointness of $\Phi$ and $\Phi^\top$.
p = 12
x = np.random.randn(n,n)
y = np.random.randn(p)
I = indices(n,p)
# must be 0
print( dotp(y,phi(x,I)) - dotp(phi_adj(y,I),x) )
-8.881784197001252e-16
Define the proximal operator of $F(x) = \iota_{\Phi \cdot = y}(x)$ (ie the ortho-projector).
def prox_F(x,y,I): return x+phi_adj(y-phi(x,I),I)
Check that prox$_F \circ$ prox$_F$ = prox$_F$ (since it is a projector).
x1 = prox_F(x,y,I)
print( np.linalg.norm(x1-prox_F(x1,y,I)) )
2.1677797545656835e-16
Compute the proximal operator of $G(x):=\|x\|_*$ the nuclear norm, which is the soft-thesholding of the singular values.
def soft_thresholding(x, u): return np.minimum(x+u, np.maximum(0,x-u))
Displays the soft thresholding.
t = np.linspace(-3,3,1000)
plt.plot(t, soft_thresholding(t,1.1));
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-1-9a2d8f83ceca> in <module>() ----> 1 t = np.linspace(-3,3,1000) 2 plt.plot(t, soft_thresholding(t,1.1)); NameError: name 'np' is not defined
def prox_G(x, gamma):
u,s,vh = np.linalg.svd(x)
soft_s=soft_thresholding(s,gamma)
return (u*soft_s)@vh
_,a,_ = np.linalg.svd(x)
_,b,_ = np.linalg.svd(prox_G(x,6))
plt.plot(a)
plt.plot(b);
Compute the symmetrized proximal operator (rprox).
def rprox_F(x,gamma,y,I):
return 2*prox_F(x,y,I)-x
def rprox_G(x,gamma):
return 2*prox_G(x,gamma)-x
def n_norm(x): return np.linalg.norm(x,'nuc')
Douglas-Rachford's algorithm.
def DR(y, mu, gamma,n_iter,I):
x_t = np.zeros((n,n))
x_list = []
x_t_list = []
nucl_norm=[]
for i in range(n_iter):
x_t_list.append(x_t.copy())
x_t=(1-mu/2)*x_t+mu/2*rprox_G(rprox_F(x_t,gamma,y,I),gamma)
x=prox_F(x_t,y,I)
x_list.append(x.copy())
nucl_norm.append(n_norm(x))
return x, nucl_norm
Test the DR algorithm.
mu = 1
gamma = 1
p = int( .5*n*n )
I = indices(n,p)
x,nucl_norm = DR(phi(x0,I), mu, gamma, 500, I )
plt.plot(np.log10( nucl_norm-np.min(nucl_norm)+1e-20 ))
[<matplotlib.lines.Line2D at 0x7f2a0ae17c90>]
Test with varying number $p$ of measurements.
r = 4
x0 = np.random.randn(n,r) @ np.random.randn(r,n)
# we varie p
n_test = 10 # number of tests
p_list = np.round( np.linspace(.02,.5,n_test)*n**2 )
err = []
n_iter = 500
import progressbar
for i in progressbar.progressbar(range(n_test)):
p = int( p_list[i] )
I = indices(n,p)
y = phi(x0,I)
x,nucl_norm = DR(y, mu, gamma, n_iter, I )
err.append( np.linalg.norm(x-x0) )
100% (10 of 10) |########################| Elapsed Time: 0:00:27 Time: 0:00:27
plt.plot(p_list/(n**2), err/np.linalg.norm(x0));
[<matplotlib.lines.Line2D at 0x7f2a0ac19b90>]