Hi Tobi, I think cholfact(X'X + lambda I)\(X'y) should work: X'X is positive (semi-)definite, lambda*I is positive definite for lambda>0, and the sum of any two positive definite matrices is also positive definite.
Alternatively, if you're worried about the numerical conditioning of X'X, you can note that the regularized least squares minimization problem minimize ||Xb - y||^2 + λ||b||^2 is really just ordinary least squares with a larger system matrix: minimize ||Ab - z||^2 where A is the block matrix [X] [λI] and z is the vector made up blockwise as [y] [0] Now it just looks like ordinary least squares with a larger set of equations, so you can solve it with whatever sparse direct or iterative solver you can get your hands on. I like this way of writing regularization conditions - it makes it clear that you're just adjoining extra linear equations to the least squares system, which makes the task of translating it into numerical code easier. My advice is to try the following in order: 1) Try the builtin backdivide first (b = A\z). Looking at the source, I'm afraid julia may not yet have support for sparse direct least squares but I could be wrong. Maybe someone has built a package for this. Certainly matlab does well here by calling SuiteSparseQR with a preordering permutation to reduce fill in (possibly something like symamd(X'X) but I couldn't get it to tell me.). 2) Try a sparse direct solver, such as SuiteSparseQR with one of the preordering permutations. You need only R from the decomposition A = QR, from which it's then quite easy to solve the least squares problem as b = R\(R'\(A'z)). As far as I understand, this is significantly better conditioned than using the normal equations directly, since it avoids forming X'X. 3) If your problem is hitting memory limits, try an iterative solver like the LSQR method. In my meager experience this involves quite a bit more fiddling since you have to choose convergence thresholds, etc. It can also be quite a bit slower than a direct solver if you want high precision! Oh, I forgot step zero 0) Make sure it really has to be solved as a sparse problem. If the full matrix fits in memory, it's always worth trying to solve it directly just to see what happens ;-) (A direct full solution is also a good comparison point to determine whether any fancy sparse stuff is giving you the right answer.) By the way, if you want to know lots about solving least squares problems, the book "Numerical Methods for Least Squares Problems" by Ake Björck is just amazing. Cheers, ~Chris On Sat, Jan 10, 2015 at 7:27 PM, Tobias Knopp <[email protected]> wrote: > Hi, > > I have learned from another mailing list thread that I con solve the normal > equation > > X'X b = X'y > > with the sparse matrix X by > > Z = X' > b = cholfact(Z)\Zy > > However, I need to solve > > (X'X + lambda I) b = X'y > > i.e. the regularized normal equation. Does someone know a direct method for > this? IIUC cholfact cannot be used for this. Alternatively it would be > interesting if someone has code (or knows algorithms) for computing the SVD > of a sparse matrix. > > Thanks > > Tobi
