# Lasso¶

Lasso is a well known technique for sparse linear regression. It is obtained by adding an $$\ell_1$$ regularization term in the objective,

$\begin{array}{ll} \mbox{minimize} & \frac{1}{2} \| Ax - b \|_2^2 + \gamma \| x \|_1 \end{array}$

where $$x \in \mathbf{R}^{n}$$ is the vector of parameters, $$A \in \mathbf{R}^{m \times n}$$ is the data matrix, and $$\gamma > 0$$ is the weighting parameter. The problem has the following equivalent form,

$\begin{split}\begin{array}{ll} \mbox{minimize} & \frac{1}{2} y^T y + \gamma \boldsymbol{1}^T t \\ \mbox{subject to} & y = Ax - b \\ & -t \le x \le t \end{array}\end{split}$

In order to get a good trade-off between sparsity of the solution and quality of the linear fit, we solve the problem for varying weighting parameter $$\gamma$$. Since $$\gamma$$ enters only in the linear part of the objective function, we can reuse the matrix factorization and enable warm starting to reduce the computation time.

## Python¶

import osqp
import numpy as np
import scipy as sp
from scipy import sparse

# Generate problem data
sp.random.seed(1)
n = 10
m = 1000
x_true = np.multiply((np.random.rand(n) > 0.8).astype(float),
np.random.randn(n)) / np.sqrt(n)
gammas = np.linspace(1, 10, 11)

# Auxiliary data
In = sparse.eye(n)
Im = sparse.eye(m)
On = sparse.csc_matrix((n, n))
Onm = sparse.csc_matrix((n, m))

# OSQP data
P = sparse.block_diag([On, sparse.eye(m), On], format='csc')
q = np.zeros(2*n + m)
sparse.hstack([In, Onm, -In]),
sparse.hstack([In, Onm, In])], format='csc')
l = np.hstack([b, -np.inf * np.ones(n), np.zeros(n)])
u = np.hstack([b, np.zeros(n), np.inf * np.ones(n)])

# Create an OSQP object
prob = osqp.OSQP()

# Setup workspace
prob.setup(P, q, A, l, u, warm_start=True)

# Solve problem for different values of gamma parameter
for gamma in gammas:
# Update linear cost
q_new = np.hstack([np.zeros(n+m), gamma*np.ones(n)])
prob.update(q=q_new)

# Solve
res = prob.solve()


## Matlab¶

% Generate problem data
rng(1)
n = 10;
m = 1000;
x_true = (randn(n, 1) > 0.8) .* randn(n, 1) / sqrt(n);
b = Ad * x_true + 0.5 * randn(m, 1);
gammas = linspace(1, 10, 11);

% OSQP data
P = blkdiag(sparse(n, n), speye(m), sparse(n, n));
q = zeros(2*n+m, 1);
speye(n), sparse(n, m), -speye(n);
speye(n), sparse(n, m), speye(n);];
l = [b; -inf*ones(n, 1); zeros(n, 1)];
u = [b; zeros(n, 1); inf*ones(n, 1)];

% Create an OSQP object
prob = osqp;

% Setup workspace
prob.setup(P, q, A, l, u, 'warm_start', true);

% Solve problem for different values of gamma parameter
for i = 1 : length(gammas)
% Update linear cost
gamma = gammas(i);
q_new = [zeros(n+m,1); gamma*ones(n,1)];
prob.update('q', q_new);

% Solve
res = prob.solve();
end


## CVXPY¶

from cvxpy import *
import numpy as np
import scipy as sp
from scipy import sparse

# Generate problem data
sp.random.seed(1)
n = 10
m = 1000
A = sparse.random(m, n, density=0.5)
x_true = np.multiply((np.random.rand(n) > 0.8).astype(float),
np.random.randn(n)) / np.sqrt(n)
b = A.dot(x_true) + 0.5*np.random.randn(m)
gammas = np.linspace(1, 10, 11)

# Define problem
x = Variable(n)
gamma = Parameter(nonneg=True)
objective = 0.5*sum_squares(A*x - b) + gamma*norm1(x)
prob = Problem(Minimize(objective))

# Solve problem for different values of gamma parameter
for gamma_val in gammas:
gamma.value = gamma_val
prob.solve(solver=OSQP, warm_start=True)


## YALMIP¶

% Generate problem data
rng(1)
n = 10;
m = 1000;
A = sprandn(m, n, 0.5);
x_true = (randn(n, 1) > 0.8) .* randn(n, 1) / sqrt(n);
b = A * x_true + 0.5 * randn(m, 1);
gammas = linspace(1, 10, 11);

% Define problem
x = sdpvar(n, 1);
gamma = sdpvar;
objective = 0.5*norm(A*x - b)^2 + gamma*norm(x,1);

% Solve with OSQP
options = sdpsettings('solver', 'osqp');
x_opt = optimizer([], objective, options, gamma, x);

% Solve problem for different values of gamma parameter
for i = 1 : length(gammas)
x_opt(gammas(i));
end