A package implementing algorithms and screening rules for the quadratic LASSO problem
Project description
qlasso
This package provides algorithms and screening rules to solve the quadratic lasso problem
\underset{x}{\min}\ \|Ax-c\|^2 + \lambda \|x\|_1^2
or, more generally, the quadratic group-lasso
\underset{X}{\min}\ \|AX-K\|^2 + \lambda \left( \sum_{i=1}^p \|x_i\|_2 \right)^2,
where $A\in\mathbb{R}^{n\times p}$, $c\in\mathbb{R}^n$, $K\in\mathbb{R}^{n\times r}$ and $x_i\in\mathbb{R}^r$ denotes the $i$th row of the matrix $X\in\mathbb{R}^{p \times r}$. It implements the techniques developped in a paper recently submitted by G. Sagnol and L. Pronzato, "Fast Screening Rules for Optimal Design via Quadratic Lasso Reformulation", which proposes new screening rules and a homotopy algorithm for solving Bayes optimal design problems. This work is a follow-up of the article Sagnol & Pauwels (2019). Statistical Papers 60(2):215--234, where it was shown that the quadratic lasso is equivalent to the problem of Bayes $c$-optimal design:
\min_{w\geq 0, \sum_{i=1}^p w_i=1}\quad c^T \left( \sum_{i=1}^p w_i a_i a_i^T + \lambda I_n \right)^{-1} c
and the quadratic group-lasso is equivalent to the problem of Bayes $L$-optimal design:
\min_{w\geq 0, \sum_{i=1}^p w_i=1}\quad \operatorname{trace}\ K^T \left( \sum_{i=1}^p w_i a_i a_i^T + \lambda I_n \right)^{-1} K,
with $a_i\in\mathbb{R}^n$ the $i$th column of $A$.
usage
Load the modules qlasso_instance
and qlasso_solvers
import qlasso.qlasso_instance as qli
import qlasso.qlasso_solvers as qls
# load a small random instance
rand_instance = qli.SLassoInstance.random_instance(p=50,n=8)
# load the mnist instance considered in the paper
mnist_instance = qli.SLassoInstance.mnist_instance(sample_size=600, resize=1)
#inspect size of `A` and `c`
print('MNIST. Shape of A:', mnist_instance.A.shape, ' -- Shape of c:', mnist_instance.c.shape)
print('RAND. Shape of A:', rand_instance.A.shape, ' -- Shape of c:', rand_instance.c.shape)
Solve the MNIST instance with the CD solver, and D1-screening rule run every 10 iterations
Acceleration should become clearly visible after roughly 50 iterations
solver = qls.CDSolver(mnist_instance,
maxit=1000,
print_frequency=1,
screening_rules=['D1'],
screen_frequency=10,
apply_screening=True,
tol=1e-6)
#solve the problem with the value lambda=0.1 for the regularization parameter
solver.solve(lbda = 0.1)
#print optimal solution of the quadratic lasso
print('optimal x-solution:')
solver.print_x()
#print the corresponding c-optimal design
print()
print('optimal design:')
solver.print_w()
Solve the small random instance with a specific solver and compare the effect of different screening rules:
# you can replace ['MWU'] with the solver of your choice
algo = {'fista': qls.FistaSolver, # FISTA
'CD': qls.CDSolver, # (Block-)Coordinate Descent
'FW': qls.FWSolver, # Frank-Wolfe
'MWU': qls.MultiplicativeSolver, # Multiplicative Weight Update
}['MWU']
solver = algo(rand_instance,
maxit=1000,
print_frequency=1,
screening_rules=['B1','B2','B3','D1','D2'],
screen_frequency=5,
apply_screening=False, #this is required to obtain the true rejection rate of each screening rule
#turn this option to True to obtain a CPU speed-up
tol=1e-6)
solver.solve(lbda = 0.1)
NB: The rule B4
considered in the paper is also implemented, but it only works for $A$-optimal design problems (i.e., with $r=n$, $K=I$).
Now, we can draw a plot to visualize the number of design points eliminated by each rule, as a function of the iteration count:
# plot the rejection rate by different screening screening_rules
import matplotlib.pyplot as plt
for rule in solver.screening_rejection:
plt.plot(solver.screen_iter, solver.screening_rejection[rule], label=rule)
plt.legend()
plt.show()
Solve a quadratic lasso instance using an adaptation of the Homotopy algorithm (LARS)
The orginal algorithm for the standard (with non-squared penalty) lasso was described in
Osborne, Presnell & Turlach (2000). IMA Journal of Numerical Analysis, 20(3):389--403.
and
Efron, Hastie, Johnstone & Tibshirani (2004). The Annals of Statistics, 32(2):407--499.
lars_solver = qls.LarsSolver(mnist_instance)
lars_solver.solve(lbda=0.1)
#print the c-optimal design
print()
print('optimal design:')
lars_solver.print_w()
Solve a quadratic lasso instance with a SOCP solver
You can replace 'cvxopt' with the name of a commercial solver installed on your system, such as gurobi
or mosek
socp_solver = qls.PicosSolver(mnist_instance)
socp_solver.solve(lbda=0.1,verbosity=1,solver='cvxopt',tol=1e-4) # solves a SOCP reformulation of the Quadratic Lasso problem
Project details
Release history Release notifications | RSS feed
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.