Low-Rank Matrix Completion
This example demonstrates how to use the positive semidefinite cone, as well as reusing a cached workspace and using warm-starting.
Matrix completion is the problem of filling in missing data into a partially observed matrix where the measurements we are given have been corrupted by Gaussian noise. In low-rank matrix completion we have the additional prior knowledge that the matrix we are completing is low-rank. For simplicity we shall also assume that the matrix we are reconstructing is symmetric positive definite. A famous instance of this problem is the Netflix prize.
Concretely, we denote by \(\hat X \in \mathbf{R}^{n \times n}\) the true matrix corrupted by noise, and denote by \(\mathcal{I}\) the set of indices (row and column pairs) from which we receive noisy observations. We shall use the nuclear norm, denoted \(\| \cdot \|_*\), as a convex surrogate for rank, which we shall trade off against the observations using regularization parameter \(\lambda \geq 0\). The low-rank matrix completion problem is given by
\[\begin{split}\begin{array}{ll}
\mbox{minimize} & \frac{1}{2} \sum_{i,j \in \mathcal{I_n}} (X_{ij} - \hat X_{ij})^2 + \lambda \|X \|_* \\
\mbox{subject to} & X \succeq 0
\end{array}\end{split}\]
over variable \(X \in \mathbf{R}^{n \times n}\) (we use \(\cdot \succeq 0\) to indicate membership in the symmetric positive semidefinite cone).
We can convert this into a more standard form. First, let \(x = \mathrm{vec}(X)\) be the semidefinite vectorization of \(X\) described in cones (and concretely implemented in the code that follows). Further, let \(A\) be the linear operator that extracts the elements of \(x\) for which we have (noisy) observations, and let \(b = A \mathrm{vec}(\hat X)\). Since the nuclear norm of a positive semidefinite matrix is given by its trace we obtain
\[\begin{split}\begin{array}{ll}
\mbox{minimize} & \frac{1}{2} y^T y + \lambda \mathrm{vec}(I_n)^\top x \\
\mbox{subject to} & y = Ax - b \\
& \mathrm{mat}(x) \succeq 0
\end{array}\end{split}\]
over variable \(x \in \mathbf{R}^{n(n+1) /2}\) and \(y \in \mathbf{R}^{|\mathcal{I}|}\), where \(I_n\) is the \(n \times n\) identity matrix. From this formulation it is straightforward to convert it into the standard form accepted by SCS. The regularization parameter \(\lambda \geq 0\) trades off the rank of the solution and the quality of the fit, and so we solve the problem for many choices of \(\lambda\). Since \(\lambda\) enters only in the linear part of the objective function, we can reuse the matrix factorization and use warm starting to reduce the computation time.
Python code to solve this is below.
import scs
import numpy as np
import scipy as sp
from scipy import sparse
np.random.seed(1)
# The vec function as documented in api/cones
def vec(S):
n = S.shape[0]
S = np.copy(S)
S *= np.sqrt(2)
S[range(n), range(n)] /= np.sqrt(2)
return S[np.triu_indices(n)]
# The mat function as documented in api/cones
def mat(s):
n = int((np.sqrt(8 * len(s) + 1) - 1) / 2)
S = np.zeros((n, n))
S[np.triu_indices(n)] = s / np.sqrt(2)
S = S + S.T
S[range(n), range(n)] /= np.sqrt(2)
return S
dim = 15 # dim x dim matrix
vlen = int(dim * (dim + 1) / 2) # length of vector x = vec(X)
# Generate true matrix
rank = dim // 5 # low rank
X = np.random.randn(dim, rank)
X = X @ X.T
#############################################################################
# Let's first do some basic sanity checks to ensure that mat, vec are working:
# mat(vec( . )) should be identity
print(f"Should be ~ 0: {np.linalg.norm(X - mat(vec(X)))}")
# Trace( . ) should be vec(I)' vec( . )
print(f"Should be ~ 0: {np.trace(X) - vec(np.eye(dim)) @ vec(X)}")
#############################################################################
num_measurements = vlen // 2 # how many measurements are revealed
# Generate random measurement indices
measurement_idxs = np.random.choice(
np.arange(vlen), size=num_measurements, replace=False
)
# Create A matrix
Ad = np.zeros((num_measurements, vlen))
for i in range(num_measurements):
Ad[i, measurement_idxs[i]] = 1.0
# Noisy measurements of X
measurements = Ad @ vec(X) + 0.01 * np.random.randn(num_measurements) # + noise
# Auxiliary data
In = sparse.eye(vlen)
Im = sparse.eye(num_measurements)
On = sparse.csc_matrix((vlen, vlen))
Onm = sparse.csc_matrix((vlen, num_measurements))
# SCS data
P = sparse.block_diag([On, sparse.eye(num_measurements)], format="csc")
A = sparse.vstack(
[
# zero cone
sparse.hstack([Ad, -Im]),
# positive semidefinite cone
sparse.hstack([-In, Onm]),
],
format="csc",
)
b = np.hstack([measurements, np.zeros(vlen)])
c = np.hstack([np.zeros(vlen + num_measurements)])
data = dict(P=P, A=A, b=b, c=c)
cone = dict(z=num_measurements, s=dim)
# Setup workspace
solver = scs.SCS(data, cone, eps_abs=1e-6, eps_rel=1e-6)
print(f"Solving for lambda = 0")
sol = solver.solve() # lambda = 0
X_hat = mat(sol["x"][:vlen])
print(f"Error: {np.linalg.norm(X_hat - X) / np.linalg.norm(X)}")
# Solve for different values of lambda
lambdas = np.logspace(-6, 1, 11)
for lam in lambdas:
print(f"Solving for lambda = {lam}")
# Re-use workspace, just update the `c` vector
c_new = np.hstack([lam * vec(np.eye(dim)), np.zeros(num_measurements)])
solver.update(c=c_new)
# Solve updated problem
sol = solver.solve() # will warm-start automatically
X_hat = mat(sol["x"][:vlen])
# What is the norm error?
print(f"Error : {np.linalg.norm(X_hat - X) / np.linalg.norm(X)}")
After following the python install instructions, we can run the code yielding output:
Should be ~ 0: 3.0513465426178085e-15
Should be ~ 0: 0.0
------------------------------------------------------------------
SCS v3.3.1 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem: variables n: 180, constraints m: 180
cones: z: primal zero / dual free vars: 60
s: psd vars: 120, ssize: 1
settings: eps_abs: 1.0e-06, eps_rel: 1.0e-06, eps_infeas: 1.0e-07
alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
max_iters: 100000, normalize: 1, rho_x: 1.00e-06
acceleration_lookback: 10, acceleration_interval: 5
lin-sys: sparse-direct-mkl-pardiso
nnz(A): 240, nnz(P): 60
Solving for lambda = 0
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 4.93e+00 1.17e+00 4.92e+01 3.27e+01 1.00e-01 8.44e-04
250| 4.27e-03 4.14e-05 2.66e-04 1.84e-04 5.60e-04 9.65e-03
300| 4.53e-09 1.16e-10 4.61e-10 -2.31e-10 5.60e-04 1.15e-02
------------------------------------------------------------------
status: solved
timings: total: 1.15e-02s = setup: 5.79e-04s + solve: 1.10e-02s
lin-sys: 2.09e-03s, cones: 7.83e-03s, accel: 5.04e-04s
------------------------------------------------------------------
objective = -0.000000
------------------------------------------------------------------
Error: 0.5886099648031258
Solving for lambda = 1e-06
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 3.45e-05 1.00e-06 4.21e-05 5.30e-05 5.60e-04 1.08e-04
250| 5.64e-03 3.52e-08 8.59e-07 4.25e-05 1.49e-06 9.14e-03
500| 1.64e-03 9.16e-09 5.52e-07 4.24e-05 1.49e-06 1.86e-02
750| 1.78e-03 4.50e-09 5.20e-07 4.27e-05 1.49e-06 2.77e-02
1000| 1.35e-03 1.24e-08 3.61e-07 4.29e-05 1.49e-06 3.70e-02
1250| 4.69e-04 8.94e-09 1.36e-07 4.30e-05 1.49e-06 4.63e-02
1500| 2.09e-05 5.60e-10 2.20e-09 4.31e-05 1.49e-06 5.56e-02
1625| 8.68e-06 4.42e-10 4.90e-10 4.31e-05 1.49e-06 6.02e-02
------------------------------------------------------------------
status: solved
timings: total: 6.02e-02s = setup: 1.03e-06s + solve: 6.02e-02s
lin-sys: 1.10e-02s, cones: 4.47e-02s, accel: 1.78e-03s
------------------------------------------------------------------
objective = 0.000043
------------------------------------------------------------------
Error : 0.03896789533150445
Solving for lambda = 5.011872336272725e-06
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 3.47e+00 1.46e-05 1.74e-05 1.41e-04 1.49e-06 1.01e-04
250| 1.70e-03 3.19e-08 2.20e-06 2.12e-04 1.49e-06 1.03e-02
500| 1.46e-03 1.94e-08 1.93e-06 2.13e-04 1.49e-06 2.06e-02
750| 1.39e-03 3.11e-08 1.47e-06 2.13e-04 1.49e-06 3.01e-02
1000| 1.62e-03 1.04e-07 4.70e-07 2.14e-04 5.27e-06 3.98e-02
1250| 4.22e-04 8.03e-08 3.57e-07 2.14e-04 5.27e-06 4.95e-02
1500| 3.58e-05 2.62e-09 6.22e-10 2.15e-04 1.66e-06 5.91e-02
1750| 1.36e-05 1.04e-09 1.79e-09 2.15e-04 1.66e-06 6.87e-02
1850| 9.29e-06 5.50e-10 1.14e-09 2.15e-04 1.66e-06 7.34e-02
------------------------------------------------------------------
status: solved
timings: total: 7.34e-02s = setup: 1.29e-06s + solve: 7.34e-02s
lin-sys: 1.39e-02s, cones: 5.40e-02s, accel: 2.18e-03s
------------------------------------------------------------------
objective = 0.000215
------------------------------------------------------------------
Error : 0.025311368120136863
Solving for lambda = 2.5118864315095822e-05
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 1.84e+01 6.27e-05 1.44e-03 2.45e-05 1.66e-06 9.79e-05
250| 2.82e-03 2.93e-07 2.93e-06 1.06e-03 6.24e-06 1.35e-02
500| 4.06e-04 9.22e-08 1.85e-07 1.06e-03 6.24e-06 2.80e-02
750| 1.90e-04 4.92e-09 1.52e-08 1.06e-03 6.24e-06 4.18e-02
1000| 8.28e-05 3.80e-09 1.97e-09 1.06e-03 6.24e-06 5.37e-02
1200| 8.84e-06 5.40e-10 3.69e-09 1.06e-03 6.24e-06 6.14e-02
------------------------------------------------------------------
status: solved
timings: total: 6.14e-02s = setup: 1.06e-06s + solve: 6.14e-02s
lin-sys: 1.35e-02s, cones: 4.24e-02s, accel: 2.61e-03s
------------------------------------------------------------------
objective = 0.001062
------------------------------------------------------------------
Error : 0.021453797312837492
Solving for lambda = 0.0001258925411794166
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 4.50e+01 2.97e-04 1.31e-02 -3.04e-03 6.24e-06 9.88e-05
250| 6.03e-03 3.38e-07 2.00e-05 5.27e-03 6.24e-06 9.75e-03
500| 1.23e-03 7.09e-08 2.00e-06 5.28e-03 6.24e-06 2.31e-02
750| 6.22e-04 3.28e-07 1.41e-06 5.28e-03 2.08e-05 3.25e-02
1000| 6.93e-05 6.86e-08 6.16e-08 5.28e-03 2.08e-05 4.13e-02
1200| 7.61e-06 2.00e-08 6.94e-09 5.28e-03 2.08e-05 5.06e-02
------------------------------------------------------------------
status: solved
timings: total: 5.06e-02s = setup: 1.20e-06s + solve: 5.06e-02s
lin-sys: 9.91e-03s, cones: 3.62e-02s, accel: 2.08e-03s
------------------------------------------------------------------
objective = 0.005278
------------------------------------------------------------------
Error : 0.03356478285269981
Solving for lambda = 0.000630957344480193
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 1.46e+02 1.33e-03 1.40e-01 -5.12e-02 2.08e-05 1.08e-04
250| 1.01e-03 1.49e-07 7.90e-06 2.62e-02 2.08e-05 1.34e-02
500| 2.52e-05 1.05e-08 1.82e-07 2.62e-02 2.08e-05 2.39e-02
575| 2.64e-06 6.97e-10 1.73e-08 2.62e-02 2.08e-05 2.67e-02
------------------------------------------------------------------
status: solved
timings: total: 2.67e-02s = setup: 1.25e-06s + solve: 2.67e-02s
lin-sys: 4.85e-03s, cones: 1.92e-02s, accel: 1.37e-03s
------------------------------------------------------------------
objective = 0.026224
------------------------------------------------------------------
Error : 0.06744945233057856
Solving for lambda = 0.0031622776601683794
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 5.19e+02 3.80e-03 2.74e+00 -1.28e+00 2.08e-05 9.95e-05
250| 1.88e-04 4.34e-06 4.42e-06 1.30e-01 1.63e-04 9.74e-03
300| 6.66e-06 1.10e-07 4.37e-08 1.30e-01 1.63e-04 1.16e-02
------------------------------------------------------------------
status: solved
timings: total: 1.16e-02s = setup: 9.42e-07s + solve: 1.16e-02s
lin-sys: 2.01e-03s, cones: 8.57e-03s, accel: 4.76e-04s
------------------------------------------------------------------
objective = 0.129971
------------------------------------------------------------------
Error : 0.11835899695782559
Solving for lambda = 0.01584893192461111
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 1.97e+02 1.85e-02 4.97e+00 -2.09e+00 1.63e-04 1.01e-04
250| 1.40e-05 6.47e-08 5.30e-09 6.43e-01 1.63e-04 9.10e-03
275| 1.20e-07 6.40e-10 5.51e-09 6.43e-01 1.63e-04 1.02e-02
------------------------------------------------------------------
status: solved
timings: total: 1.02e-02s = setup: 9.82e-07s + solve: 1.02e-02s
lin-sys: 1.88e-03s, cones: 7.25e-03s, accel: 6.05e-04s
------------------------------------------------------------------
objective = 0.642509
------------------------------------------------------------------
Error : 0.17415994089940848
Solving for lambda = 0.07943282347242805
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 9.15e+02 6.42e-02 1.20e+02 -5.83e+01 1.63e-04 9.79e-05
225| 6.50e-06 2.55e-08 5.99e-08 3.15e+00 5.50e-04 8.59e-03
------------------------------------------------------------------
status: solved
timings: total: 8.59e-03s = setup: 8.81e-07s + solve: 8.59e-03s
lin-sys: 1.61e-03s, cones: 6.24e-03s, accel: 2.93e-04s
------------------------------------------------------------------
objective = 3.149158
------------------------------------------------------------------
Error : 0.2226257261519734
Solving for lambda = 0.3981071705534969
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 1.24e+03 3.17e-01 7.24e+02 -3.54e+02 5.50e-04 9.79e-05
150| 5.81e-06 7.49e-07 1.17e-08 1.49e+01 3.61e-03 5.76e-03
------------------------------------------------------------------
status: solved
timings: total: 5.76e-03s = setup: 8.71e-07s + solve: 5.76e-03s
lin-sys: 1.06e-03s, cones: 4.12e-03s, accel: 2.94e-04s
------------------------------------------------------------------
objective = 14.885088
------------------------------------------------------------------
Error : 0.3057863717395143
Solving for lambda = 1.9952623149688788
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 7.18e+02 1.58e+00 2.01e+03 -9.77e+02 3.61e-03 9.82e-05
100| 3.59e-09 1.46e-10 2.35e-09 6.15e+01 3.61e-03 4.12e-03
------------------------------------------------------------------
status: solved
timings: total: 4.13e-03s = setup: 8.11e-07s + solve: 4.12e-03s
lin-sys: 6.80e-04s, cones: 2.70e-03s, accel: 5.51e-04s
------------------------------------------------------------------
objective = 61.454584
------------------------------------------------------------------
Error : 0.5124766604261639
Solving for lambda = 10.0
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 3.61e+03 8.00e+00 5.12e+04 -2.57e+04 3.61e-03 1.15e-04
125| 5.03e-10 4.66e-08 1.65e-10 1.32e+02 5.64e-02 5.01e-03
------------------------------------------------------------------
status: solved
timings: total: 5.02e-03s = setup: 9.11e-07s + solve: 5.01e-03s
lin-sys: 8.88e-04s, cones: 3.27e-03s, accel: 5.84e-04s
------------------------------------------------------------------
objective = 132.180038
------------------------------------------------------------------
Error : 0.9959215533308903