Thanks Chris your answer is spot on. I totally forgot about that stacking method which is also very useful to show that the 1st and 2nd normal equation are equivalent if one applies Tikhonov regularization.
I am already using iterative solvers for my problem. Interestingly the Kaczmarz method provides extremely good results for my matrix. Since it is challenging to choose lambda and the number of iterations in practice I want to look if a direct solver would work better. To you points: 1) If \ does not support sparse matrices this will be too slow (I am inherently assuming that the cholesky decomposition is also sparse...) 2) Can you point me to Julia code where a sparse QR is implemented? I know the argument about avoiding the normal equations quite well but in practice this does not always matter (e.g. if the matrix is an operator that you do not want to arrange explicitly) 3) Do you have a pointer to LSQR Julia code? Step 0) is already done. I get speedups of factor 100 and and this is quite crucial as the solver should be used in an online reconstruction framework. Cheers, Tobi Am Samstag, 10. Januar 2015 14:21:07 UTC+1 schrieb Chris Foster: > > 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] <javascript:>> 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 >
