Lasso
This example demonstrates quadratic objectives, as well as reusing a cached workspace and using warm-starting.
In the lasso the goal is to find a sparse vector that fits some measurements. The \(\ell_1\) norm is used as a convex surrogate for sparsity, and a regularization parameter \(\lambda \geq 0\) trades off sparsity and quality of fit. Concretely the lasso solves
\[\begin{array}{ll}
\mbox{minimize} & \frac{1}{2} \| Ax - b \|_2^2 + \lambda \| x \|_1
\end{array}\]
over variable \(x \in \mathbf{R}^{n}\), with data \(A \in \mathbf{R}^{m \times n}\) and \(b \in \mathbf{R}^n\). The problem has the following equivalent form,
\[\begin{split} \begin{array}{ll}
\mbox{minimize} & \frac{1}{2} y^T y + \lambda {\bf 1}^T t \\
\mbox{subject to} & y = Ax - b \\
& -t \le x \le t
\end{array}\end{split}\]
over variables \(x \in \mathbf{R}^{n}\), \(t \in \mathbf{R}^{n}\), \(y \in \mathbf{R}^{m}\). 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 sparsity 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
# Generate problem data
sp.random.seed(1)
np.random.seed(1)
n = 200 # Variables
m = 100 # Measurements
Ad = sparse.random(m, n, density=0.5) # Measurement matrix
x_true = sparse.random(n, 1, density=0.1) # True sparse vector
x_true = np.array(x_true.todense()).squeeze()
measurements = Ad @ x_true + 0.1 * np.random.randn(m)
measurements = np.array(measurements).squeeze()
# The smallest value of lambda with solution all-zeros
lambda_max = np.linalg.norm(Ad.T @ measurements, np.inf)
# Auxiliary data
In = sparse.eye(n)
Im = sparse.eye(m)
On = sparse.csc_matrix((n, n))
Onm = sparse.csc_matrix((n, m))
# SCS data
P = sparse.block_diag([On, sparse.eye(m), On], format="csc")
q = np.zeros(2 * n + m)
A = sparse.vstack(
[
# zero cone
sparse.hstack([Ad, -Im, Onm.T]),
# positive cones
sparse.hstack([In, Onm, -In]),
sparse.hstack([-In, Onm, -In]),
],
format="csc",
)
b = np.hstack([measurements, np.zeros(n), np.zeros(n)])
c = np.zeros(2 * n + m)
data = dict(P=P, A=A, b=b, c=c)
cone = dict(z=m, l=2 * n)
print(f"Solving for lambda = 0")
# Setup workspace
solver = scs.SCS(
data,
cone,
eps_abs=1e-6,
eps_rel=1e-6,
)
sol = solver.solve() # lambda = 0
x = sol["x"][:n]
print(f"Error : {np.linalg.norm(x_true - x) / np.linalg.norm(x_true)}")
# Solve for different values of lambda
lambdas = np.logspace(-2, np.log10(lambda_max), 11)
for lam in lambdas:
print(f"Solving for lambda = {lam}")
# Re-use workspace, just update the `c` vector
c_new = np.hstack([np.zeros(n + m), lam * np.ones(n)])
solver.update(c=c_new)
# Solve updated problem
sol = solver.solve() # will warm-start automatically
x = sol["x"][:n]
# What is the norm error?
print(f"Error : {np.linalg.norm(x_true - x) / np.linalg.norm(x_true)}")
After following the python install instructions, we can run the code yielding output:
Solving for lambda = 0
------------------------------------------------------------------
SCS v3.1.1 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem: variables n: 500, constraints m: 500
cones: z: primal zero / dual free vars: 100
l: linear vars: 400
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): 10900, nnz(P): 100
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 4.54e+00 7.55e-02 9.57e-01 4.89e-01 1.00e-01 1.68e-04
25| 1.56e-07 5.84e-08 3.28e-07 1.64e-07 1.00e-01 1.36e-03
------------------------------------------------------------------
status: solved
timings: total: 7.39e-03s = setup: 6.02e-03s + solve: 1.37e-03s
lin-sys: 1.06e-03s, cones: 3.43e-05s, accel: 6.73e-06s
------------------------------------------------------------------
objective = 0.000000
------------------------------------------------------------------
Error : 0.7256063655877213
Solving for lambda = 0.01
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 1.41e-01 1.00e-02 7.69e-02 3.85e-02 1.00e-01 1.48e-04
250| 1.34e-04 1.56e-06 1.40e-06 1.25e-01 1.69e-02 1.14e-02
500| 5.72e-05 8.48e-07 4.82e-07 1.25e-01 1.69e-02 2.19e-02
750| 2.01e-05 1.20e-07 5.32e-07 1.25e-01 1.69e-02 3.25e-02
875| 4.30e-06 1.16e-07 2.10e-07 1.25e-01 1.69e-02 3.79e-02
------------------------------------------------------------------
status: solved
timings: total: 3.80e-02s = setup: 4.35e-06s + solve: 3.80e-02s
lin-sys: 3.20e-02s, cones: 7.18e-04s, accel: 6.56e-04s
------------------------------------------------------------------
objective = 0.125039
------------------------------------------------------------------
Error : 0.2259107853405316
Solving for lambda = 0.02497591197515687
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 1.26e+00 1.50e-02 6.09e+00 -2.92e+00 1.69e-02 1.59e-04
250| 8.05e-04 1.01e-05 1.23e-04 3.10e-01 1.69e-02 1.37e-02
500| 2.98e-04 1.20e-06 1.85e-05 3.10e-01 1.69e-02 2.76e-02
750| 1.20e-04 1.01e-06 9.39e-06 3.10e-01 1.69e-02 4.12e-02
1000| 9.60e-05 7.37e-07 3.64e-06 3.10e-01 1.69e-02 5.46e-02
1250| 4.72e-05 8.34e-07 6.33e-06 3.10e-01 1.69e-02 6.76e-02
1500| 2.25e-05 3.04e-07 2.25e-06 3.10e-01 1.69e-02 7.83e-02
1750| 5.62e-06 2.28e-08 4.66e-07 3.10e-01 1.69e-02 8.91e-02
1950| 4.71e-06 6.17e-08 1.22e-07 3.10e-01 1.69e-02 9.78e-02
------------------------------------------------------------------
status: solved
timings: total: 9.78e-02s = setup: 4.45e-06s + solve: 9.78e-02s
lin-sys: 8.29e-02s, cones: 2.08e-03s, accel: 2.57e-03s
------------------------------------------------------------------
objective = 0.309833
------------------------------------------------------------------
Error : 0.21880188151191785
Solving for lambda = 0.06237961789907842
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 3.14e+00 3.74e-02 3.87e+01 -1.90e+01 1.69e-02 1.61e-04
250| 2.30e-03 2.35e-05 4.98e-04 7.59e-01 1.69e-02 1.09e-02
500| 6.38e-04 8.53e-06 1.13e-04 7.59e-01 1.69e-02 2.16e-02
750| 8.74e-04 9.97e-06 1.78e-05 7.59e-01 1.69e-02 3.20e-02
1000| 1.71e-04 3.89e-07 1.42e-05 7.59e-01 1.69e-02 4.27e-02
1250| 8.82e-05 2.24e-07 6.85e-06 7.59e-01 1.69e-02 5.30e-02
1500| 8.80e-06 1.58e-08 3.94e-07 7.59e-01 1.69e-02 6.34e-02
1675| 5.44e-06 8.79e-09 1.93e-07 7.59e-01 1.69e-02 7.04e-02
------------------------------------------------------------------
status: solved
timings: total: 7.04e-02s = setup: 4.46e-06s + solve: 7.04e-02s
lin-sys: 6.07e-02s, cones: 1.32e-03s, accel: 1.40e-03s
------------------------------------------------------------------
objective = 0.759146
------------------------------------------------------------------
Error : 0.19472289935884313
Solving for lambda = 0.15579878456913032
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 7.83e+00 9.34e-02 2.43e+02 -1.21e+02 1.69e-02 1.52e-04
250| 6.59e-03 2.83e-05 1.68e-03 1.81e+00 1.69e-02 1.07e-02
500| 5.37e-03 4.13e-05 8.56e-04 1.82e+00 1.69e-02 2.15e-02
750| 2.12e-03 7.02e-06 3.70e-04 1.82e+00 1.69e-02 3.22e-02
1000| 9.68e-04 1.98e-06 1.49e-04 1.82e+00 1.69e-02 4.30e-02
1250| 9.76e+00 1.16e-01 1.21e-01 1.66e+00 1.69e-02 5.40e-02
1500| 6.24e-05 2.01e-07 3.74e-06 1.82e+00 1.69e-02 6.48e-02
1750| 1.50e-05 1.46e-07 6.33e-07 1.82e+00 1.69e-02 7.51e-02
1900| 4.08e-06 1.56e-08 4.53e-08 1.82e+00 1.69e-02 8.15e-02
------------------------------------------------------------------
status: solved
timings: total: 8.15e-02s = setup: 4.07e-06s + solve: 8.15e-02s
lin-sys: 7.04e-02s, cones: 1.47e-03s, accel: 1.75e-03s
------------------------------------------------------------------
objective = 1.818814
------------------------------------------------------------------
Error : 0.16652260857219253
Solving for lambda = 0.3891216729235027
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 1.96e+01 2.33e-01 1.52e+03 -7.58e+02 1.69e-02 1.54e-04
250| 1.60e-02 1.13e-04 2.43e-03 4.22e+00 1.69e-02 1.15e-02
500| 3.91e-03 6.14e-05 1.38e-03 4.26e+00 5.49e-02 2.35e-02
750| 1.10e-03 8.02e-06 8.38e-04 4.27e+00 1.81e-01 3.44e-02
1000| 6.18e-06 3.09e-07 1.36e-06 4.27e+00 1.81e-01 4.43e-02
1025| 4.23e-06 1.61e-07 5.60e-07 4.27e+00 1.81e-01 4.53e-02
------------------------------------------------------------------
status: solved
timings: total: 4.53e-02s = setup: 3.94e-06s + solve: 4.53e-02s
lin-sys: 3.80e-02s, cones: 8.57e-04s, accel: 6.28e-04s
------------------------------------------------------------------
objective = 4.269514
------------------------------------------------------------------
Error : 0.16508930825747362
Solving for lambda = 0.9718668650563185
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 4.55e+00 5.83e-01 8.78e+02 -4.35e+02 1.81e-01 1.58e-04
250| 4.54e-03 3.43e-05 3.36e-03 1.01e+01 1.81e-01 1.07e-02
500| 6.93e-04 3.79e-06 4.76e-04 1.01e+01 1.81e-01 2.15e-02
725| 1.15e-07 5.49e-08 1.98e-10 1.01e+01 5.85e-01 3.20e-02
------------------------------------------------------------------
status: solved
timings: total: 3.20e-02s = setup: 3.95e-06s + solve: 3.20e-02s
lin-sys: 2.73e-02s, cones: 5.78e-04s, accel: 5.07e-04s
------------------------------------------------------------------
objective = 10.143419
------------------------------------------------------------------
Error : 0.17887188557683356
Solving for lambda = 2.427326127321828
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 3.52e+00 1.46e+00 1.69e+03 -8.37e+02 5.85e-01 5.24e-04
250| 9.19e-05 3.91e-06 1.18e-04 2.44e+01 5.85e-01 1.16e-02
350| 3.72e-08 1.71e-08 1.53e-07 2.44e+01 5.85e-01 1.63e-02
------------------------------------------------------------------
status: solved
timings: total: 1.63e-02s = setup: 3.98e-06s + solve: 1.63e-02s
lin-sys: 1.38e-02s, cones: 2.89e-04s, accel: 2.45e-04s
------------------------------------------------------------------
objective = 24.402281
------------------------------------------------------------------
Error : 0.2605486204939154
Solving for lambda = 6.062468369098836
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 8.79e+00 3.64e+00 1.06e+04 -5.29e+03 5.85e-01 1.55e-04
250| 8.09e-05 8.61e-05 3.41e-05 5.78e+01 5.85e-01 1.09e-02
325| 2.57e-06 1.02e-07 5.29e-06 5.78e+01 5.85e-01 1.46e-02
------------------------------------------------------------------
status: solved
timings: total: 1.46e-02s = setup: 4.46e-06s + solve: 1.46e-02s
lin-sys: 1.25e-02s, cones: 2.66e-04s, accel: 2.45e-04s
------------------------------------------------------------------
objective = 57.753523
------------------------------------------------------------------
Error : 0.4805364286021535
Solving for lambda = 15.141567633878541
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 2.20e+01 9.08e+00 6.64e+04 -3.31e+04 5.85e-01 1.54e-04
150| 1.74e-07 2.20e-06 8.87e-06 1.31e+02 6.39e+00 8.10e-03
------------------------------------------------------------------
status: solved
timings: total: 8.11e-03s = setup: 4.38e-06s + solve: 8.10e-03s
lin-sys: 6.28e-03s, cones: 1.46e-04s, accel: 1.09e-04s
------------------------------------------------------------------
objective = 130.731598
------------------------------------------------------------------
Error : 0.728897685138138
Solving for lambda = 37.81744603896349
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 5.02e+00 2.27e+01 3.78e+04 -1.88e+04 6.39e+00 3.06e-04
125| 1.11e-07 4.49e-07 3.88e-05 2.66e+02 6.39e+00 9.53e-03
------------------------------------------------------------------
status: solved
timings: total: 9.54e-03s = setup: 6.99e-06s + solve: 9.54e-03s
lin-sys: 7.88e-03s, cones: 2.06e-04s, accel: 1.81e-04s
------------------------------------------------------------------
objective = 265.601214
------------------------------------------------------------------
Error : 0.8689013367110517
Solving for lambda = 94.45252033943963
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 1.25e+01 5.66e+01 2.37e+05 -1.18e+05 6.39e+00 2.22e-04
250| 1.32e-06 1.19e-12 4.39e-03 3.80e+02 1.60e+01 1.32e-02
325| 1.47e-14 1.61e-08 2.15e-11 3.80e+02 6.73e+04 1.69e-02
------------------------------------------------------------------
status: solved
timings: total: 1.69e-02s = setup: 6.50e-06s + solve: 1.69e-02s
lin-sys: 1.28e-02s, cones: 2.90e-04s, accel: 2.29e-04s
------------------------------------------------------------------
objective = 379.601832
------------------------------------------------------------------
Error : 1.0000000000000042