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.1.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: 10
lin-sys: sparse-direct
nnz(A): 240, nnz(P): 60
Solving for lambda = 0
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 5.12e+00 8.58e-01 3.93e+01 2.40e+01 1.00e-01 7.29e-03
250| 7.51e-03 3.59e-05 2.58e-04 2.03e-04 3.19e-03 1.49e-02
500| 7.41e-03 1.60e-06 6.49e-05 -1.89e-05 9.33e-04 2.48e-02
750| 3.89e-03 1.04e-08 1.56e-06 -7.37e-07 9.33e-04 3.47e-02
800| 1.67e-08 3.76e-10 3.11e-12 -1.56e-12 9.33e-04 3.67e-02
------------------------------------------------------------------
status: solved
timings: total: 3.85e-02s = setup: 1.81e-03s + solve: 3.67e-02s
lin-sys: 3.24e-03s, cones: 3.16e-02s, accel: 3.59e-04s
------------------------------------------------------------------
objective = -0.000000
------------------------------------------------------------------
Error: 2.4700619951554748
Solving for lambda = 1e-06
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 3.04e-06 1.00e-06 1.76e-04 1.21e-04 9.33e-04 1.01e-04
250| 2.31e-02 9.45e-09 7.98e-07 4.11e-05 1.00e-06 1.03e-02
500| 7.27e-03 7.34e-10 4.85e-07 4.16e-05 1.00e-06 2.12e-02
750| 4.47e-03 3.69e-10 4.00e-07 4.18e-05 1.00e-06 3.06e-02
1000| 3.61e-03 1.61e-10 3.42e-07 4.19e-05 1.00e-06 3.78e-02
1250| 3.49e-03 1.59e-10 3.38e-07 4.19e-05 1.00e-06 4.49e-02
1500| 3.10e-03 9.65e-11 3.18e-07 4.20e-05 1.00e-06 5.18e-02
1750| 3.22e-03 3.06e-08 2.75e-07 4.20e-05 1.36e-05 5.90e-02
2000| 1.94e-03 2.37e-09 3.94e-07 4.22e-05 1.36e-05 6.64e-02
2250| 1.81e-03 5.11e-10 4.73e-07 4.23e-05 1.36e-05 7.35e-02
2500| 1.72e-03 5.65e-10 5.08e-07 4.23e-05 1.36e-05 8.07e-02
2750| 1.63e-03 6.25e-10 5.40e-07 4.23e-05 1.36e-05 8.83e-02
3000| 1.53e-03 7.00e-10 5.67e-07 4.24e-05 1.36e-05 9.56e-02
3250| 1.43e-03 4.32e-09 6.40e-07 4.24e-05 4.32e-05 1.03e-01
3500| 1.21e-03 6.07e-09 6.68e-07 4.25e-05 4.32e-05 1.09e-01
3750| 1.03e-03 9.64e-09 7.59e-07 4.27e-05 4.32e-05 1.16e-01
4000| 4.61e-04 2.64e-08 2.03e-07 4.31e-05 4.32e-05 1.23e-01
4250| 2.54e-04 2.09e-09 1.05e-07 4.30e-05 4.32e-05 1.29e-01
4500| 7.91e-05 1.08e-09 3.60e-08 4.31e-05 4.32e-05 1.36e-01
4750| 3.03e-05 1.24e-09 5.48e-09 4.31e-05 4.32e-05 1.42e-01
4775| 9.42e-06 5.65e-10 5.58e-09 4.31e-05 4.32e-05 1.43e-01
------------------------------------------------------------------
status: solved
timings: total: 1.43e-01s = setup: 2.17e-06s + solve: 1.43e-01s
lin-sys: 1.76e-02s, cones: 1.17e-01s, accel: 1.83e-03s
------------------------------------------------------------------
objective = 0.000043
------------------------------------------------------------------
Error : 0.039165998317111615
Solving for lambda = 5.011872336272725e-06
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 1.50e-01 4.95e-06 3.31e-05 1.92e-04 4.32e-05 8.42e-05
250| 3.43e-03 1.34e-07 1.58e-06 2.11e-04 4.32e-05 7.29e-03
500| 1.73e-03 2.11e-09 1.81e-06 2.11e-04 4.32e-05 1.47e-02
750| 1.59e-03 1.01e-09 1.90e-06 2.12e-04 4.32e-05 2.22e-02
1000| 1.52e-03 9.87e-10 1.95e-06 2.12e-04 4.32e-05 2.98e-02
1250| 1.34e-03 9.31e-09 2.10e-06 2.12e-04 1.38e-04 3.71e-02
1500| 1.01e-03 1.04e-08 2.13e-06 2.13e-04 1.38e-04 4.48e-02
1750| 6.80e-04 8.59e-09 1.73e-06 2.14e-04 1.38e-04 5.23e-02
2000| 1.45e+00 2.00e-04 1.53e-06 2.15e-04 1.38e-04 5.95e-02
2250| 4.87e-04 1.13e-08 1.40e-06 2.14e-04 1.38e-04 6.65e-02
2500| 7.75e-04 8.76e-08 1.49e-06 2.16e-04 1.38e-04 7.36e-02
2750| 6.69e-04 2.14e-08 2.14e-06 2.15e-04 1.38e-04 8.19e-02
3000| 2.18e-04 7.41e-09 7.61e-07 2.15e-04 1.38e-04 8.90e-02
3250| 2.17e-04 3.06e-08 1.10e-07 2.15e-04 1.38e-04 9.61e-02
3500| 2.40e-05 4.38e-10 7.45e-08 2.15e-04 1.38e-04 1.03e-01
3600| 4.83e-06 5.28e-10 4.07e-09 2.15e-04 1.38e-04 1.06e-01
------------------------------------------------------------------
status: solved
timings: total: 1.06e-01s = setup: 1.96e-06s + solve: 1.06e-01s
lin-sys: 1.29e-02s, cones: 8.62e-02s, accel: 1.09e-03s
------------------------------------------------------------------
objective = 0.000215
------------------------------------------------------------------
Error : 0.025334768895438196
Solving for lambda = 2.5118864315095822e-05
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 2.44e-01 2.49e-05 1.34e-04 9.44e-04 1.38e-04 9.16e-05
250| 1.61e-03 1.94e-07 5.77e-06 1.06e-03 1.38e-04 7.88e-03
500| 1.25e-03 5.66e-09 5.22e-06 1.06e-03 1.38e-04 1.55e-02
750| 1.08e-03 4.03e-09 4.49e-06 1.06e-03 1.38e-04 2.26e-02
1000| 8.60e-04 1.78e-09 3.69e-06 1.06e-03 1.38e-04 3.04e-02
1250| 7.46e-04 1.21e-08 3.55e-06 1.06e-03 4.46e-04 3.85e-02
1500| 5.59e-04 7.69e-09 2.78e-06 1.06e-03 4.46e-04 4.60e-02
1750| 4.57e-04 1.01e-08 2.73e-06 1.06e-03 4.46e-04 5.33e-02
2000| 2.15e-04 3.68e-09 1.35e-06 1.06e-03 4.46e-04 6.05e-02
2250| 9.98e-05 1.54e-09 6.55e-07 1.06e-03 4.46e-04 6.74e-02
2500| 7.51e-06 5.47e-10 4.88e-08 1.06e-03 4.46e-04 7.47e-02
------------------------------------------------------------------
status: solved
timings: total: 7.47e-02s = setup: 1.88e-06s + solve: 7.47e-02s
lin-sys: 9.02e-03s, cones: 6.11e-02s, accel: 8.98e-04s
------------------------------------------------------------------
objective = 0.001062
------------------------------------------------------------------
Error : 0.021416258779558488
Solving for lambda = 0.0001258925411794166
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 3.71e-01 1.17e-04 4.59e-04 4.61e-03 4.46e-04 8.17e-05
250| 1.44e-03 4.05e-07 1.33e-05 5.27e-03 4.46e-04 7.40e-03
500| 9.07e-04 4.06e-07 5.92e-06 5.28e-03 4.46e-04 1.52e-02
750| 1.14e-04 2.21e-09 1.42e-06 5.28e-03 4.46e-04 2.22e-02
1000| 1.70e-05 3.22e-10 1.49e-07 5.28e-03 4.46e-04 2.97e-02
1025| 9.37e-06 1.07e-10 1.06e-07 5.28e-03 4.46e-04 3.05e-02
------------------------------------------------------------------
status: solved
timings: total: 3.05e-02s = setup: 1.79e-06s + solve: 3.05e-02s
lin-sys: 3.54e-03s, cones: 2.50e-02s, accel: 5.09e-04s
------------------------------------------------------------------
objective = 0.005279
------------------------------------------------------------------
Error : 0.03355024082447967
Solving for lambda = 0.000630957344480193
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 2.03e+00 5.72e-04 7.00e-03 1.84e-02 4.46e-04 8.29e-05
250| 6.36e-03 1.41e-06 1.78e-04 2.61e-02 4.46e-04 7.31e-03
500| 2.35e-02 1.04e-05 5.01e-05 2.62e-02 4.46e-04 1.45e-02
750| 7.01e-04 1.89e-07 2.18e-05 2.62e-02 4.46e-04 2.16e-02
1000| 2.42e-05 1.05e-09 3.62e-07 2.62e-02 4.46e-04 2.87e-02
1100| 7.41e-06 1.86e-10 7.34e-08 2.62e-02 4.46e-04 3.16e-02
------------------------------------------------------------------
status: solved
timings: total: 3.16e-02s = setup: 1.74e-06s + solve: 3.16e-02s
lin-sys: 3.65e-03s, cones: 2.57e-02s, accel: 6.26e-04s
------------------------------------------------------------------
objective = 0.026224
------------------------------------------------------------------
Error : 0.06744664704036951
Solving for lambda = 0.0031622776601683794
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 1.04e+01 2.74e-03 2.66e-01 -2.34e-02 4.46e-04 9.11e-05
250| 8.77e-03 1.73e-05 2.41e-05 1.30e-01 2.01e-03 7.52e-03
500| 7.31e-05 3.12e-08 9.26e-07 1.30e-01 6.56e-03 1.47e-02
575| 1.47e-06 9.60e-10 4.27e-08 1.30e-01 6.56e-03 1.68e-02
------------------------------------------------------------------
status: solved
timings: total: 1.68e-02s = setup: 1.78e-06s + solve: 1.68e-02s
lin-sys: 1.92e-03s, cones: 1.39e-02s, accel: 1.96e-04s
------------------------------------------------------------------
objective = 0.129971
------------------------------------------------------------------
Error : 0.11835628778123863
Solving for lambda = 0.01584893192461111
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 3.47e+00 1.39e-02 3.91e-01 3.49e-01 6.56e-03 8.42e-05
250| 1.03e-04 1.62e-07 1.25e-05 6.43e-01 6.56e-03 7.30e-03
375| 1.71e-06 1.57e-09 1.05e-07 6.43e-01 6.56e-03 1.15e-02
------------------------------------------------------------------
status: solved
timings: total: 1.15e-02s = setup: 1.73e-06s + solve: 1.15e-02s
lin-sys: 1.25e-03s, cones: 9.60e-03s, accel: 1.12e-04s
------------------------------------------------------------------
objective = 0.642509
------------------------------------------------------------------
Error : 0.17415959268201173
Solving for lambda = 0.07943282347242805
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 1.82e+01 6.55e-02 1.19e+01 -3.25e+00 6.56e-03 8.98e-05
225| 4.87e-07 4.73e-09 1.27e-07 3.15e+00 3.61e-02 7.46e-03
------------------------------------------------------------------
status: solved
timings: total: 7.47e-03s = setup: 1.72e-06s + solve: 7.47e-03s
lin-sys: 7.51e-04s, cones: 6.26e-03s, accel: 5.70e-05s
------------------------------------------------------------------
objective = 3.149158
------------------------------------------------------------------
Error : 0.22262545190934194
Solving for lambda = 0.3981071705534969
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 1.65e+01 3.21e-01 5.33e+01 -1.46e+01 3.61e-02 8.96e-05
125| 7.27e-07 1.89e-08 8.50e-07 1.49e+01 1.83e-01 4.27e-03
------------------------------------------------------------------
status: solved
timings: total: 4.27e-03s = setup: 1.76e-06s + solve: 4.27e-03s
lin-sys: 4.41e-04s, cones: 3.58e-03s, accel: 2.71e-05s
------------------------------------------------------------------
objective = 14.885089
------------------------------------------------------------------
Error : 0.3057864439961101
Solving for lambda = 1.9952623149688788
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 1.67e+01 1.55e+00 2.52e+02 -9.34e+01 1.83e-01 8.80e-05
50| 2.65e-07 5.94e-09 4.41e-07 6.15e+01 1.83e-01 1.86e-03
------------------------------------------------------------------
status: solved
timings: total: 1.86e-03s = setup: 1.67e-06s + solve: 1.86e-03s
lin-sys: 1.67e-04s, cones: 1.56e-03s, accel: 7.00e-06s
------------------------------------------------------------------
objective = 61.454584
------------------------------------------------------------------
Error : 0.5124766513093888
Solving for lambda = 10.0
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 8.75e+01 7.97e+00 6.59e+03 -3.72e+03 1.83e-01 9.20e-05
50| 1.31e-06 2.45e-10 6.05e-08 1.32e+02 1.83e-01 1.81e-03
------------------------------------------------------------------
status: solved
timings: total: 1.81e-03s = setup: 1.71e-06s + solve: 1.81e-03s
lin-sys: 1.73e-04s, cones: 1.49e-03s, accel: 1.87e-05s
------------------------------------------------------------------
objective = 132.180037
------------------------------------------------------------------
Error : 0.9959215771443783