# Large scale ridge regression

I'm trying to solve a problem of the form

$\min_x \frac{1}{2}||Ax-b||^2_2 + \frac{\rho}{2}||x-z||^2_F$

where both $x$ and $b$ are high dimensional, and $b$ is much higher dimensional than $x$. The solution is given by $x^* = (A^T A+\rho I)^{-1}(A^T b + z)$, but the problem is so large that even inverting $A^T A + \rho I$ is infeasible. However, due to structure in the problem we can efficiently multiply by $A$ and $A^T$. Basically this is large scale linear ridge regression. What would be the ideal algorithm for efficiently implementing this minimization? Would something like biconjugate gradient work?

• How large is large? Also, when calculating $x^*$, we don't actually invert the matrix. Typical methods are Cholesky decomposition (when $A^TA + \rho I$ is well-conditioned), QR decomposition, and SVD. Cholesky is the fastest, and since a large part of the point of the formulation is to make $A^TA$ well-conditioned, may well work for you. – jbowman Oct 28 '13 at 20:33
• $x$ can easily be tens of thousands of dimensions, $b$ can be hundreds of millions. Working with video data. – David Pfau Oct 28 '13 at 20:35
• Is it sparse data? If so, you may be you find this articule interesting. statslab.cam.ac.uk/~rds37/papers/regres6l.pdf -It's about a min-hash technique for very large ( x and b) sparse data. – Manuel Oct 28 '13 at 21:24

## 1 Answer

I've found that LSQR is ideal for problems like this - I've used it successfully for operators of about 3e5 * 1e6 or so. Check http://www.stanford.edu/group/SOL/software/lsqr.html for details. I've used Friedlander's (I think) C port and the python port, which I have (hastily and sloppily) ported to R.

• This looks like the right thing for me. I'm using Matlab, and there is a built-in implementation for the unregularized case. I can't find a Matlab implementation of the regularized case, however. – David Pfau Oct 28 '13 at 22:45
• Did you check stanford.edu/group/SOL/software/lsqr/matlab? – alex Oct 28 '13 at 23:55
• Yes, that version has regularization built in. Also, I realized you could just fold the regularization term into the $Ax-b$ term. Either way, it works. – David Pfau Oct 29 '13 at 1:10