constrained_lasso: the exact path with p > n, under constraints, and what a ridge costs

preprints/constrained_lasso/constrained_lasso.tex

@@ -43,6 +43,7 @@
 \input{figures/gaines.tex}
 \input{figures/monotone.tex}
 \input{figures/stress_numbers.tex}
+\input{figures/wide_numbers.tex}
 
 \newtheorem{proposition}{Proposition}
 \theoremstyle{remark}
@@ -74,7 +75,10 @@
 On caps, signs, equality and inequality rows, monotone fits and zero sums, on simulated
 and real data, up to 500 predictors and on problems built to break a path algorithm,
 the paths reproduce LARS and an interior-point solver and meet the LASSO's optimality
-conditions to rounding, or the CLA stops with an error. 
+conditions to rounding, or the CLA stops with an error. With more predictors than
+observations the paths stay exact under bounds, signs, inequality rows and zero sums.
+The path algorithms for general constraints add a ridge there, and the ridge moves the
+end of the path.
 
 \medskip\noindent
 \textbf{Keywords:} constrained LASSO; critical line algorithm; solution path;
@@ -193,7 +197,9 @@ \section{A constrained regression}
 (Section~\ref{sec:limits}). For $C\beta=0$
 alone, the Python package \texttt{c-lasso} \cite{classo} has a path algorithm with a
 small ridge term added. Neither was derived with portfolios in mind, so
-Section~\ref{sec:evidence} uses both as independent checks. Hobson, Butcher and Lentati
+Section~\ref{sec:evidence} uses both as independent checks. For zero sums on blocks of
+coefficients, the comlasso algorithm of Jeon et al.\ \cite{jeon2020} traces the exact
+path, also with $p>n$; it takes no bounds or inequality rows. Hobson, Butcher and Lentati
 \cite{hobson2026} revisit Gaines et al.\ with an algorithm of their own, CLASSO. The
 linear system of Gaines et al.\ is singular whenever a binding row vanishes on the
 active coefficients, as a sign constraint on a zero coefficient does, so CLASSO takes
@@ -220,6 +226,17 @@ \section{A constrained regression}
 Each is a few lines on top of a released CLA. Readers who know the CLA can skip
 Section~\ref{sec:cla}.
 
+\paragraph{More predictors than observations} With $p>n$ the encodings part from the
+other path algorithms for general constraints. Gaines et al., CLASSO and
+\texttt{c-lasso} add a ridge to $X\trans X$ there, and so trace a nearby problem. With
+$p>n$ that problem is not close at its end: the ridge decides which interpolating fit the
+path runs to. As shipped, \texttt{c-lasso} ends up to \clWideGapHi{} times the largest
+coefficient away from the exact path on a zero-sum problem with $p>n$
+(Section~\ref{sec:evidence}). The CLA needs no ridge, and traces the exact path under
+bounds, signs, inequality rows and zero sums alike (Table~\ref{tab:wide}). Microbiome
+studies often have more taxa than samples, so this is the case a user of the zero-sum
+LASSO meets first.
+
 \paragraph{Notation} The design is $X\in\R^{n\times p}$, with $n$ observations and $p$
 predictors, and the response is $y\in\R^n$. Write $G=X\trans X$ and $c=X\trans y$. The
 constrained LASSO is
@@ -623,6 +640,19 @@ \section{Numerical evidence}
 near $10^{16}$ it has the largest residual in the table, \stKKT{}
 (Section~\ref{sec:limits}).
 
+\paragraph{More predictors than observations} \emph{Wide, zero-sum} is the only
+problem so far with both $p>n$ and a constraint, and a zero sum alone has an exact path
+in comlasso \cite{jeon2020} too. Table~\ref{tab:wide} puts on \emph{Wide}, again with
+$n=\wdN$, $p=\wdP$ and \stSeeds{} seeds, the constraints comlasso does not take:
+nonnegativity; a box at half the largest coefficient of the least-squares fit of
+smallest norm; two inequality rows that cap the net effect of three predictors, as in
+\emph{Cholesterol}, and bind; the row of \emph{Pinned}; and a zero sum with the signs
+of \emph{Zero-sum, signs} or with the caps of the box. \texttt{cvxcla} declines none of
+them. Every path meets the optimality conditions to \wdKKT{} or better along every
+segment, and the interior-point solver agrees to \wdGap{}. Its objective is never lower
+than the path's by more than \wdIpmBetter{}, relative, and that only near the end of a
+path, where the fit almost interpolates and the objective is small.
+
 \paragraph{Two independent path algorithms} Two path algorithms derived for the LASSO,
 with no reference to portfolios, can be run against the CLA. Both trace the same path.
 
@@ -649,11 +679,15 @@ \section{Numerical evidence}
 its start can differ from the one measured here.
 
 \texttt{c-lasso} \cite{classo} has a path algorithm for $C\beta=0$ alone, without bounds
-or signs, so it runs on \emph{Zero-sum} and \emph{sCD14}. It adds a ridge term
-$10^{-3}I$ to $X\trans X$, fixed in the package, and that term is the whole of its
-distance to the CLA: its \clPoints{} grid points lie up to \clGap{} from the CLA path as
-shipped, relative to the largest coefficient, and \clGapUnridged{} with the ridge set
-to $10^{-12}$ by hand.
+or signs, so it runs on \emph{Zero-sum}, \emph{sCD14} and \emph{Wide, zero-sum}. It
+adds a ridge term $10^{-3}I$ to $X\trans X$, fixed in the package, and that term is the
+whole of its distance to the CLA: on the first two its \clPoints{} grid points lie up to
+\clGap{} from the CLA path as shipped, relative to the largest coefficient, and
+\clGapUnridged{} with the ridge set to $10^{-12}$ by hand. With $p>n$ the same ridge is
+not small. On \emph{Wide, zero-sum}, with its grid carried down to \clWideLamin{} times
+its largest $\lambda$ so that it reaches the end of the path, \texttt{c-lasso} as
+shipped lies \clWideGapLo{} to \clWideGapHi{} times the largest coefficient from the CLA
+over the \stSeeds{} seeds, and \clWideGapUnridged{} with the ridge at $10^{-12}$.
 
 Two further packages fit on a grid of $\lambda$ rather than tracing the path, so they
 cannot check the knots and are not run: \texttt{zeroSum} \cite{zerosum}, for zero-sum
@@ -689,6 +723,21 @@ \section{Numerical evidence}
 \label{tab:stress}
 \end{table}
 
+\begin{table}[t]
+\centering
+\small
+\begin{tabular}{lrcc}
+\toprule
+\emph{Wide} with & knots & gap & KKT \\
+\midrule
+\input{figures/wide.tex}
+\end{tabular}
+\caption{More predictors than observations, $n=\wdN$ and $p=\wdP$, under constraints
+beyond a zero sum, through the checks of Table~\ref{tab:evidence} against the
+interior-point solver, the worst over \stSeeds{} seeds. None is declined.}
+\label{tab:wide}
+\end{table}
+
 
 \section{Where the two forms stop}
 \label{sec:limits}
@@ -711,12 +760,18 @@ \section{Where the two forms stop}
 The path ends where the fit first interpolates, the limit $\lambda\to0^+$: the
 interpolating fit of smallest $\ell_1$ norm under the constraints. Past it the CLA would
 run on through interpolating fits of larger norm, which are not on the LASSO path, so
-\texttt{cvxcla} stops there. \emph{Wide, zero-sum} (Table~\ref{tab:stress}) tests this.
-The other path algorithms tested here do not give this path: Gaines et al.\ add a ridge
-of $10^{-4}$ when $n<p$, and \texttt{c-lasso} always adds one of $10^{-3}$. Nor does
-CLASSO \cite{hobson2026}, not run here, which adds one of $10^{-4}\max_j G_{jj}$. It is the
-case a user of the zero-sum LASSO meets first, since microbiome data often have more
-taxa than samples.
+\texttt{cvxcla} stops there. \emph{Wide, zero-sum} (Table~\ref{tab:stress}) tests this,
+and Table~\ref{tab:wide} tests it under bounds, signs and inequality rows. The other path
+algorithms for general constraints do not give this path. Gaines et al.\ add a ridge of
+$10^{-4}$ when $n<p$, \texttt{c-lasso} always adds one of $10^{-3}$, and CLASSO
+\cite{hobson2026}, not run here, adds one of $10^{-4}\max_j G_{jj}$. Each then traces a
+nearby problem whose path ends elsewhere. On \emph{Wide, zero-sum} the ridge of
+\texttt{c-lasso} moves its path by up to \clWideGapHi{} times the largest coefficient
+(Section~\ref{sec:evidence}). On their own problems with $p>n$, Hobson et al.\ report
+objectives up to $4\times10^{-2}$, relative, above the unregularised ones at their
+default ridge
+\cite[Table~3]{hobson2026}. For a zero sum alone, comlasso \cite{jeon2020} traces the
+exact path with $p>n$ and no ridge, and its authors make the same case against one.
 
 \paragraph{Numerically singular free blocks} The bordered system \eqref{eq:border}
 needs the columns on the free set to be independent. When two columns nearly coincide
@@ -834,8 +889,10 @@ \section{Conclusion}
 exact knots or stops with an error when its arithmetic cannot be trusted. Where an
 encoding starts degenerate, \texttt{cvxcla} relaxes the constraint, extrapolates and
 certifies the result, so the checks reach monotone fits, caps and rows under a zero sum
-and weights that sum to one; with a singular Gram they reach zero sums with $p>n$. Only
-rows that exclude $\beta=0$ stay out. The checks also made the CLA better at
+and weights that sum to one. With $p>n$ they stay exact, with no ridge, under bounds,
+signs, inequality rows and zero sums, where the other path algorithms for general
+constraints trace a ridged problem that ends elsewhere. Only rows that exclude $\beta=0$
+stay out. The checks also made the CLA better at
 portfolios, where it now takes tied forecasts and low-rank covariances. The practical
 consequence is simple: \texttt{lars\_path}
 cannot take even a cap on the coefficients, and a CLA that takes general linear rows
@@ -933,6 +990,12 @@ \section*{Acknowledgements}
 2020.
 \newblock \href{https://doi.org/10.1080/01621459.2019.1609970}{\texttt{doi:10.1080/01621459.2019.1609970}}.
 
+\bibitem{jeon2020}
+J.-J. Jeon, Y.~Kim, S.~Won, and H.~Choi.
+\newblock Primal path algorithm for compositional data analysis.
+\newblock \emph{Computational Statistics \& Data Analysis}, 148:106958, 2020.
+\newblock \href{https://doi.org/10.1016/j.csda.2020.106958}{\texttt{doi:10.1016/j.csda.2020.106958}}.
+
 
 \bibitem{lawson1974}
 C.~L. Lawson and R.~J. Hanson.

preprints/constrained_lasso/python/classo_check.py

@@ -16,7 +16,8 @@
 LASSO under ``C beta = 0``. Its last release, 1.0.11, uses ``np.infty``, which NumPy 2
 removed, so it runs here in its own environment, pinned to NumPy 1.26, and talks to
 ``make_figures.py`` through a file: that script writes the exact paths of the Zero-sum
-problem and of the sCD14 data to ``figures/zero_sum_paths.npz``, and this one runs
+problem, of the sCD14 data and of Wide, zero-sum (p > n, one per seed) to
+``figures/zero_sum_paths.npz``, and this one runs
 c-lasso's ``Path-Alg`` on the same problems, measures how far its answers lie from the
 exact path, and writes the numbers the note quotes to ``figures/classo.tex``. c-lasso
 also imports pandas and matplotlib without declaring them, hence the two pins above.
@@ -25,6 +26,9 @@
 ``eps_L2 * I``, with ``eps_L2 = 1e-3``, to ``X^T X``, a default of ``parameters_for_update``
 that no public setting reaches. That ridge term is all of its distance to the exact path,
 so the script runs it twice, as shipped and with the ridge lowered to ``RIDGE`` by hand.
+With n >= p the ridge moves the path only slightly. With p > n it decides where the path
+ends: the ridged path heads for a different interpolating fit. There c-lasso's grid is
+carried down to ``WIDE_LAMIN`` times its largest lambda, so that it reaches the end.
 
 The line ``# run: after figures`` above makes ``make figures`` run this script after
 ``make_figures.py``, here and in CI.
@@ -46,9 +50,10 @@
 POINTS = 100  # the grid c-lasso reports its path on
 REPEATS = 5  # each run timed this often, and the fastest counts
 RIDGE = 1e-12  # eps_L2 for the second run, in place of the shipped 1e-3
+WIDE_LAMIN = 1e-9  # c-lasso's smallest lambda, relative, for p > n; the shipped 1e-3 stops short of the end
 
 
-def classo_path(x: np.ndarray, y: np.ndarray) -> tuple[np.ndarray, np.ndarray, float]:
+def classo_path(x: np.ndarray, y: np.ndarray, lamin: float = 1e-3) -> tuple[np.ndarray, np.ndarray, float]:
     """c-lasso's Path-Alg on the zero-sum LASSO: lambda in the note's scaling, beta, seconds.
 
     c-lasso's R1 problem is ``||y - X b||^2 + lam ||b||_1`` under ``C b = 0``, without
@@ -64,6 +69,7 @@ def classo_path(x: np.ndarray, y: np.ndarray) -> tuple[np.ndarray, np.ndarray, f
         selection.PATH, selection.CV, selection.StabSel, selection.LAMfixed = True, False, False, False
         selection.PATHparameters.numerical_method = "Path-Alg"
         selection.PATHparameters.Nlam = POINTS
+        selection.PATHparameters.lamin = lamin
         start = time.perf_counter()
         problem.solve()
         best = min(best, time.perf_counter() - start)
@@ -101,22 +107,32 @@ def main() -> None:
     shipped, unridged = {}, {}
     for name in sorted({key.split(".")[0] for key in data.files}):
         x, y, knots, exact = data[f"{name}.x"], data[f"{name}.y"], data[f"{name}.lam"], data[f"{name}.beta"]
-        lam, beta, seconds = classo_path(x, y)
+        lamin = WIDE_LAMIN if name.startswith("wide") else 1e-3
+        lam, beta, seconds = classo_path(x, y, lamin)
         shipped[name] = gap(lam, beta, knots, exact)
         with_ridge(RIDGE)
-        lam, beta, _ = classo_path(x, y)
+        lam, beta, _ = classo_path(x, y, lamin)
         unridged[name] = gap(lam, beta, knots, exact)
         with_ridge(1e-3)
         print(f"{name}: c-lasso Path-Alg, {len(lam)} grid points in {seconds * 1e3:.1f} ms, largest gap to "
               f"the exact path {shipped[name]:.1e}, and {unridged[name]:.1e} with the ridge at {RIDGE:.0e}")
-        # As shipped c-lasso is close to the path but off it by its ridge; without the
-        # ridge it is on it. Flag it if either ever changes.
-        assert 1e-9 < shipped[name] < 1e-2, shipped[name]
+        # As shipped c-lasso is close to the path but off it by its ridge, and with p > n
+        # far off it; without the ridge it is on it. Flag it if any of that ever changes.
+        if name.startswith("wide"):
+            assert shipped[name] > 1e-1, shipped[name]
+        else:
+            assert 1e-9 < shipped[name] < 1e-2, shipped[name]
         assert unridged[name] < 1e-10, unridged[name]
+    square = [name for name in shipped if not name.startswith("wide")]
+    wide = [name for name in shipped if name.startswith("wide")]
     (FIGURES / "classo.tex").write_text(
         rf"\newcommand{{\clPoints}}{{{POINTS}}}" + "\n"
-        + rf"\newcommand{{\clGap}}{{{sci(max(shipped.values()))}}}" + "\n"
-        + rf"\newcommand{{\clGapUnridged}}{{{sci(max(unridged.values()))}}}" + "\n"
+        + rf"\newcommand{{\clGap}}{{{sci(max(shipped[n] for n in square))}}}" + "\n"
+        + rf"\newcommand{{\clGapUnridged}}{{{sci(max(unridged[n] for n in square))}}}" + "\n"
+        + rf"\newcommand{{\clWideGapLo}}{{{min(shipped[n] for n in wide):.2f}}}" + "\n"
+        + rf"\newcommand{{\clWideGapHi}}{{{max(shipped[n] for n in wide):.2f}}}" + "\n"
+        + rf"\newcommand{{\clWideGapUnridged}}{{{sci(max(unridged[n] for n in wide))}}}" + "\n"
+        + rf"\newcommand{{\clWideLamin}}{{{sci(WIDE_LAMIN)}}}" + "\n"
     )