5
$\begingroup$

I have to solve the following least squares problem: \begin{equation} \| \left[ \begin{smallmatrix} \mathbf{L} \\ \mathbf{I} \end{smallmatrix} \right]\mathbf{x} - \mathbf{b} \|_2^2 \end{equation} 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.

Improve this question
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – BrunoLevy
    Oct 1, 2015 at 15:03
  • $\begingroup$ 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). $\endgroup$
    – Srdan
    Oct 2, 2015 at 9:09

1 Answer 1

Reset to default
4
$\begingroup$

This looks isomorphic to Tikhonov regularization (a.k.a. ridge regression). Some googling brought up LSQR and the newer LSMR. Those links both have implementations in a number of languages. For large scale problems, Petsc has KSPLSQR built in.

Depending on what you mean by 'large scale', mlpack may also work. Mlpack has a tutorial covering both its command line program as well as its C++ api.

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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). $\endgroup$
    – Srdan
    Oct 2, 2015 at 9:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.