Thanks to both Doug and Dominique for your very informative replies. For now, I've implemented the cholfact method suggested by Doug. It's providing a very nice speedup relative to my original naive implementation of sparse backslash on the normal equations.
--Peter On Thursday, July 31, 2014 9:26:31 AM UTC-7, Dominique Orban wrote: > > On Wednesday, July 30, 2014 11:08:45 PM UTC-4, Peter Simon wrote: > >> Is there such a method in Julia? Right now I'm multiplying both sides by >> the transpose of the coefficient matrix to form and solve the normal >> equations using \. I suspect that this is not very efficient compared to a >> sparse least squares solver. >> >> Thanks, >> --Peter >> > > It depends on how large/sparse your overdetermined matrix A is. As > factorizations go, you have a choice between performing the Cholesky > factorization of A'A or the QR factorization of A. Some smart Cholesky > Note that if A = QR, then A'A = R'R and R' is the Cholesky factor of A'A. > However, it's often more stable to compute the QR factorization. The > trouble with it is that if your A is large, Q is also large and is > typically dense. In Matlab, there's a way to request the Q-less QR > factorization: > > R = qr(A) > > but I'm not seeing that in the Julia documentation. It might be there and > not be apparent (to me). With R alone, you can solve the least-squares > problem, but it's recommended to perform one step of iterative refinement. > > If factorization isn't an option, the first thing to try is to run the > conjugate gradient method on the normal equations (don't form them, just > encode the appropriate matrix-vecor product). You may want to look at my > linear operator library: https://github.com/dpo/linop.jl > > There are specialized iterative methods for least-squares problems, that > are typically more stable than CG (in the same way that QR is more stable > than Cholesky). One of them is LSQR and is available in Julia here: > https://github.com/JuliaLang/IterativeSolvers.jl/blob/master/src/lsqr.jl > (I've never tested it). There are others. > > I hope this helps. > >
