reduce lineardiff convex problem to 2 * order variables via thin QR #237

lineardiff solves, per window of length $N$,

$$\min_{\mathbf{z}} \Vert \mathbf{y} - \mathbf{M}\mathbf{z}\Vert_2^2 + \Vert \boldsymbol{\lambda} \odot \mathbf{z}\Vert_1, \qquad \mathbf{M} = [\textstyle\int \mathbf{Y};\ \mathbf{B}]^\top \in \mathbb{R}^{N\times 2k},\quad \mathbf{z} = [\mathbf{a}; \mathbf{c}],\quad \boldsymbol{\lambda} = [10^{-6}\mathbf{1};\ \gamma\mathbf{1}]$$

where $\odot$ is elementwise multiplication. CVXPY thus sees an $N$-length residual despite only $2k$ unknowns ($k$ = order). It reduces exactly to a $2k$-sized problem.

Derivation

Thin QR: $\mathbf{M} = \mathbf{Q}\mathbf{R}$, $\mathbf{Q} \in \mathbb{R}^{N\times 2k}$ with orthonormal columns, $\mathbf{R} \in \mathbb{R}^{2k\times 2k}$ upper triangular. Decompose $\mathbf{y} = \mathbf{Q}\mathbf{Q}^\top \mathbf{y} + (\mathbf{I} - \mathbf{Q}\mathbf{Q}^\top)\mathbf{y}$. The fit $\mathbf{M}\mathbf{z} = \mathbf{Q}\mathbf{R}\mathbf{z}$ lies in $\text{range}(\mathbf{Q})$, orthogonal to $(\mathbf{I} - \mathbf{Q}\mathbf{Q}^\top)\mathbf{y}$, so by Pythagoras

$$\Vert \mathbf{y} - \mathbf{M}\mathbf{z}\Vert_2^2 = \Vert \mathbf{Q}\mathbf{Q}^\top \mathbf{y} - \mathbf{Q}\mathbf{R}\mathbf{z}\Vert_2^2 + \Vert(\mathbf{I} - \mathbf{Q}\mathbf{Q}^\top)\mathbf{y}\Vert_2^2 = \Vert \mathbf{Q}^\top \mathbf{y} - \mathbf{R}\mathbf{z}\Vert_2^2 + \text{const},$$

using $\Vert \mathbf{Q}\mathbf{v}\Vert_2 = \Vert \mathbf{v}\Vert_2$. The constant is the part of $\mathbf{y}$ no $\mathbf{z}$ can fit. With $\mathbf{q} = \mathbf{Q}^\top \mathbf{y} \in \mathbb{R}^{2k}$:

$$\min_{\mathbf{z}} \Vert \mathbf{R}\mathbf{z} - \mathbf{q}\Vert_2^2 + \Vert \boldsymbol{\lambda} \odot \mathbf{z}\Vert_1$$

The data enter only through $(\mathbf{R}, \mathbf{q})$, as sufficient statistics. Missing samples, masked by *obs, zero whole rows of $\mathbf{M}$ and $\mathbf{y}$ and contribute to neither term.

Why not the Gram form. Expanding gives $\mathbf{z}^\top \mathbf{M}^\top \mathbf{M} \mathbf{z} - 2(\mathbf{M}^\top \mathbf{y})^\top \mathbf{z} + \text{const}$, equally $N$-independent, with $\mathbf{M}^\top \mathbf{M} = \mathbf{R}^\top \mathbf{R}$. But $\text{cond}(\mathbf{M}^\top \mathbf{M}) = \text{cond}(\mathbf{M})^2$, whereas Householder QR yields $\mathbf{R}$ directly with $\text{cond}(\mathbf{R}) = \text{cond}(\mathbf{M})$.

Column scaling. CVXPY canonicalizes $\Vert \mathbf{R}\mathbf{z} - \mathbf{q}\Vert_2^2$ into the quadratic $\mathbf{z}^\top \mathbf{R}^\top \mathbf{R} \mathbf{z}$ regardless, so the squaring reappears inside the solver. The columns of $\mathbf{M}$ (repeated integrals, $\tau^j/j!$) differ by orders of magnitude even after the $\tau$ scaling, so $\text{cond}(\mathbf{R}) \sim 10^4\text{–}10^5$ (4e4 at order 3, irregular steps). The unreduced problem tolerates that, but squared to ~1e9 Clarabel stalls at user_limit with wrong coefficients. Let $\mathbf{D} = \text{diag}(\mathbf{d})$, $d_i = \Vert \mathbf{R}_{:,i}\Vert_2$, and substitute $\mathbf{z} = \mathbf{D}^{-1}\mathbf{w}$:

$$\min_{\mathbf{w}} \Vert(\mathbf{R}\mathbf{D}^{-1})\mathbf{w} - \mathbf{q}\Vert_2^2 + \Vert (\boldsymbol{\lambda} \oslash \mathbf{d}) \odot \mathbf{w}\Vert_1, \qquad \mathbf{z} = \mathbf{w} \oslash \mathbf{d}$$

where $\oslash$ is elementwise division. It is still a weighted $\ell_1$ problem and exactly equivalent. $\mathbf{R}\mathbf{D}^{-1}$ has unit-norm columns (Jacobi preconditioning), and since $\Vert \mathbf{R}_{:,i}\Vert_2 = \Vert \mathbf{Q}\mathbf{R}_{:,i}\Vert_2 = \Vert \mathbf{M}_{:,i}\Vert_2$, this is column normalization of $\mathbf{M}$, computed from $\mathbf{R}$. It brought cond from ~9e4 to ~3e2, so ~1e5 after squaring.

The total is still $O(Nk^2)$ for the QR, but everything CVXPY and the solver touch is $N$-independent.

Consequences

• Caching: the problem shape depends only on order, so _PROBLEM_CACHE can be keyed by order alone instead of also depending on window length. The 64-entry clear and the separate B_p parameter both go away. Currently truncated edge windows add ~11 keys per run, at 10–25 ms of canonicalization each, also incurred for every window_size tried in optimize.
• Prototype: tests pass, including equivariance; outputs agree to ~1e-9; one run is ~1.3× faster. The gain is small because CVXPY's fixed per-call overhead (~1.7 ms) dominates.

Filed by Claude (Claude Code) on behalf of @pavelkomarov.