Skip to content

README.md

Latest commit

 

History

History
189 lines (140 loc) · 6.04 KB

README.md

File metadata and controls

189 lines (140 loc) · 6.04 KB

IPsolve – Python

A Python implementation of the interior-point solver for piecewise linear-quadratic (PLQ) composite problems, based on the JMLR paper:

Aravkin, Burke, Pillonetto. Sparse/Robust Estimation and Kalman Smoothing with Nonsmooth Log-Concave Densities: Modeling, Computation, and Theory. JMLR 14(82):2689–2728, 2013.

Quick start

cd python
pip install -e ".[test]"           # editable install with test deps
pytest                             # run 12 tests (verified vs CVXPY)
from ipsolve import solve, huber, l1, l2, hinge, qreg

# Least squares
x = solve(H, z)

# Robust regression
x = solve(H, z, meas=huber(kappa=1.0))

# Lasso
x = solve(H, z, proc=l1(lam=0.5))

# Huber + L1 (robust sparse)
x = solve(H, z, meas=huber(kappa=1.0), proc=l1(lam=0.5))

# SVM  (was: proc_mMult=1.0/lam ... now just proc=l2(lam=lam))
w = solve(yX, np.ones(m), meas=hinge(), proc=l2(lam=0.01))

# Quantile regression
x = solve(H, z, meas=qreg(tau=0.25))

# Box constraints
x = solve(H, z, bounds=(-1, 1))

What it solves

$$ \min_{x} ;; \ell^{!\top} x ;+; \rho_{\text{meas}}!\bigl(z - Hx\bigr) ;+; \rho_{\text{proc}}!\bigl(Kx - k\bigr) \quad\text{s.t.}\quad Ax \le a $$

where $\rho_{\text{meas}}$ and $\rho_{\text{proc}}$ are PLQ penalties (ℓ₁, ℓ₂, Huber, hinge, quantile, Vapnik, …) represented in their conjugate/dual form via matrices $(M, B, C, c, b)$.

The solver uses a barrier interior-point method with Schur-complement elimination and optional Mehrotra predictor-corrector.

Design

Module Purpose
ipsolve.plq PLQ dataclass, compose(), ensure_sparse helpers
ipsolve.penalties Penalty library (10 penalties -- see table below)
ipsolve.kkt SchurOperator for mat-vec/dense Schur complement; KKT residual & solve
ipsolve.solver Interior-point loop with extracted _max_step / _armijo_search helpers
ipsolve.api High-level solve(H, z, meas, proc, ...) entry point

Penalties

Each penalty is a factory function. Call without a dimension for solve():

meas=huber(kappa=1.0)   # lazy spec -- dimension filled in by solve()
proc=l1(lam=0.5)

Call with a dimension to get the raw PLQ object:

plq = huber(100, kappa=1.0)   # returns PLQ directly
Name Formula Key parameter
l2 (lam/2) |v|^2 lam (or mMult)
l1 lam * |v|_1 lam
huber Huber loss kappa
hinge lam * sum pos(v) lam
vapnik lam * sum max(|v| - eps, 0) lam, eps
qreg check / pinball loss tau, lam
qhuber quantile Huber (smooth quantile) tau, kappa
infnorm lam * |v|_inf lam
logreg sum log(1 + exp(v)) alpha
hybrid sum(sqrt(1 + v^2/s) - 1) scale

Examples

All examples produce figures saved to figures/. Run from the repo root:

python examples/lasso.py
python examples/robust_regression.py
python examples/svm.py
python examples/huber_l1.py
python examples/quantile_regression.py
python examples/cvar_portfolio.py
python examples/infnorm.py
python examples/kalman_demo.py
Example Description MATLAB equivalent
lasso.py Sparse signal recovery via ℓ₂ + ℓ₁ lassoAndEnet.m
robust_regression.py Huber loss vs LS with outliers RobustRegression.m
svm.py Binary classification with hinge + ℓ₂ SVMdemo.m
huber_l1.py Robust sparse regression, verified vs CVXPY HuberL1Example.m
quantile_regression.py Multi-quantile curves with pinball loss —
cvar_portfolio.py CVaR portfolio with simplex constraints CvarDemo.m
infnorm.py ℓ∞-norm regularisation (simplex dual) infNorm.m
kalman_demo.py Kalman smoothing: LS, robust, constrained KalmanDemo.m

Tests

12 tests verify all penalties against CVXPY ground truth:

pytest -v

Tests cover: ℓ₂, Huber, lasso, Huber+ℓ₁, SVM, ℓ₁ fitting, quantile regression, ℓ₁+ℓ₁, Vapnik, box constraints, ℓ∞ norm, and Kalman smoothing.

Future work

  • Nonlinear extension -- Gauss-Newton outer loop for rho(F(x)) where F is nonlinear
  • Logistic regression / hybrid examples

Inexact (CG) Schur solve

The SchurOperator class in kkt.py encapsulates the Schur complement Omega = B' T^{-1} B as either a dense matrix (for direct LU) or a matrix-free LinearOperator (for conjugate-gradient).

Pass inexact=True to solve() to use the CG path. This avoids forming and factoring the dense (N x N) Schur complement, reducing memory from O(N^2) to O(nnz) and scaling better for overdetermined problems (m >> n).

Scaling benchmark

python benchmarks/scaling_demo.py

Part 1 — Direct (LU) vs CG comparison

Lasso (ℓ₂ + ℓ₁, underdetermined m < n):

m n Direct (s) IP CG CG (s) IP CG Rel diff
100 200 0.08 17 0 0.22 17 1960 3.4e-8
200 500 0.66 20 0 1.15 20 4588 6.9e-8
500 1000 3.89 20 0 5.26 20 4499 5.1e-5

Huber regression (overdetermined m ≫ n) — CG wins:

m n Direct (s) IP CG CG (s) IP CG Speedup Rel diff
500 200 0.09 10 0 0.08 10 233 1.2× 4.5e-9
1000 500 1.09 13 0 0.55 13 413 2.0× 1.7e-9
2000 1000 10.14 16 0 2.84 16 522 3.6× 5.4e-10

Part 2 — CG-only scaling (large problems)

Lasso (CG-only):

m n Time (s) IP iters CG iters Objective
500 1000 5.4 20 4499 24.59
1000 2000 29.4 23 6209 45.15
2000 5000 369.5 36 16221 120.27

Huber (CG-only):

m n Time (s) IP iters CG iters
2000 1000 2.5 16 522
5000 2000 9.0 14 351
10000 5000 70.6 16 536