On Sat, Jan 10, 2015 at 11:46 PM, Tobias Knopp <tobias.kn...@googlemail.com> wrote: > 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.
That's interesting, I didn't know about the Kaczmarz method. Is it a standard method for your problem domain? > 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...) Correct. The docs suggest that there's no support for this yet in 0.4 dev: help?> \ Base.\(A, B) Matrix division using a polyalgorithm. For input matrices "A" and "B", the result "X" is such that "A*X == B" when "A" is square. The solver that is used depends upon the structure of "A". A direct solver is used for upper- or lower triangular "A". For Hermitian "A" (equivalent to symmetric "A" for non- complex "A") the "BunchKaufman" factorization is used. Otherwise an LU factorization is used. For rectangular "A" the result is the minimum-norm least squares solution computed by a pivoted QR factorization of "A" and a rank estimate of A based on the R factor. For sparse, square "A" the LU factorization (from UMFPACK) is used. To be honest, the last time I did this kind of stuff was in matlab (hence the comment about matlab's \). My experience was to try some iterative solvers first in anticipation that they'd be required, only to find that the matlab \ correctly detected a much better non-iterative method for my problem. A little depressing, but instructive :-) > 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) Well sometimes settling for less precision and faster code is the right thing anyway. I'd have no qualms about using the normal equations if they do the job for the problem at hand. I assume nobody has written a sparse direct QR in julia, especially given that suitesparseqr is available. Actually, there's an old issue about this: https://github.com/JuliaLang/julia/issues/3554 > 3) Do you have a pointer to LSQR Julia code? There appears to be an lsqr implementation in IterativeSolvers.jl. Cheers, ~Chris
