# Large scale triangular least squares

I have to solve the following least squares problem: $$\| \left[ \begin{smallmatrix} \mathbf{L} \\ \mathbf{I} \end{smallmatrix} \right]\mathbf{x} - \mathbf{b} \|_2^2$$ where $\mathbf{L} \in \mathbb{R}^{n\times n}$ is $O(n)$ sparse lower triangular matrix, $\mathbf{I} \in \mathbb{R}^{n \times n}$ is the identity and $\mathbf{b} = \left[ \begin{smallmatrix} \mathbf{b}_1 \\ \mathbf{b}_2 \end{smallmatrix} \right] \in \mathbb{R}^{2n}$.

Hence, solving individual system $\mathbf{Lx} = \mathbf{b}_1$ is of $O(n)$ complexity, by the forward substitution algorithm, but the least squares fit is expensive.

I am open to any suggestion, including fast approximate stochastic solvers etc. Of course, it would be perfect if one is aware of a direct method that exploits this kind of structure.

• The normal equation is like $(L^tL + I)x = L^tb_1 + b_2$. If the conditioning is not too bad, maybe you can try to solve it with conjugate gradient, this just requires a function to compute $y = (L^tL + I)x$, which can be decomposed in $y_1 = Lx; y= L^t y_1 + x$ if you want to avoid assembling $L^tL$. You may also want to compute $Diag(L^tL + I)$ for preconditioning it with Jacobi. – BrunoLevy Oct 1 '15 at 15:03
• The problem is that $\mathbf{L}$ comes from discretization of a PDE, so it is (increasingly) ill-conditioned. Forming the normal equations would square the condition number, of course. Shortly, we tried CG with different preconditioners, and it does work up to some scale. I was hoping more for a method that scales linearly with dimension (e.g. algebraic multigrid - like). – Srdan Oct 2 '15 at 9:09

• Indeed, it is in a way Tikhonov regularization, and we are actually using LSMR right now, to avoid building the normal equations. By "large scale" I mean "potentially huge scale" (right now, about $O(10^7)$, but this is actually quite modest). – Srdan Oct 2 '15 at 9:12