Method for obtaining least-square solution of sparse overdetermined linear system?

[Peter Simon]

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

[Douglas Bates]

Sparse least squares is more frequently performed by solving the sparse, symmetric positive-definite system

 X'X b = X'y

than by, say, taking a sparse QR decomposition of X.  However, it is important to use a sparse Cholesky factorization of X'X rather than, say, an LU factorization.  The sparse Cholesky is the best understood and best developed of all direct sparse solvers. 

The cholfact method for sparse matrices does work with rectangular matrices in the formo Z = X'.  That is

b = cholfact(Z)\Zy

solves

ZZ' b = Z y

Here is an example using the Julia sparse version of the matrix afiro defined in test/suitesparse.jl of the Julia sources

julia> tt = sparse(afiro)

27x51 sparse matrix with 102 Float64 entries:

[3 ,  1]  =  1.0

[4 ,  2]  =  1.0

[7 ,  3]  =  1.0

[8 ,  4]  =  1.0

[9 ,  5]  =  1.0

[10,  6]  =  1.0

[13,  7]  =  1.0

⋮

[21, 47]  =  2.279

[14, 48]  =  1.4

[16, 48]  =  -1.0

[16, 49]  =  1.0

[25, 49]  =  -1.0

[15, 50]  =  1.0

[27, 50]  =  1.0

[16, 51]  =  1.0

julia> z = sparse(afiro)

27x51 sparse matrix with 102 Float64 entries:

[3 ,  1]  =  1.0

[4 ,  2]  =  1.0

[7 ,  3]  =  1.0

[8 ,  4]  =  1.0

[9 ,  5]  =  1.0

[10,  6]  =  1.0

[13,  7]  =  1.0

⋮

[21, 47]  =  2.279

[14, 48]  =  1.4

[16, 48]  =  -1.0

[16, 49]  =  1.0

[25, 49]  =  -1.0

[15, 50]  =  1.0

[27, 50]  =  1.0

[16, 51]  =  1.0

julia> srand(1234321)

julia> y = rand(size(z,2))

51-element Array{Float64,1}:

 0.0944218

 0.936611 

 0.258327 

 0.930924 

 0.555283 

 0.87151  

 0.041553 

 0.968779 

 0.653566 

 0.458101 

 ⋮        

 0.0239988

 0.926619 

 0.522457 

 0.110692 

 0.487905 

 0.669727 

 0.33555  

 0.212379 

 0.849889 

julia> b = cholfact(z)\z*y

27-element Array{Float64,1}:

  0.809255 

 -0.923453 

  0.250299 

  0.2811   

  0.146158 

 -0.276562 

  0.620823 

  0.772582 

  0.434435 

  0.606513 

  ⋮        

  0.674793 

  0.235333 

  0.495153 

  0.338174 

  0.394105 

  0.602457 

  0.0148562

  1.21874  

  0.566282 

julia> sol = z'b

51-element Array{Float64,1}:

 0.250299

 0.2811  

 0.620823

 0.772582

 0.434435

 0.606513

 0.36439 

 0.658524

 0.356608

 0.443731

 ⋮       

 0.406332

 0.728272

 0.655013

 0.438805

 0.89312 

 0.776574

 0.130504

 0.610959

 0.14536 

julia> norm(z * (y - sol))

5.9349937888025225e-15

[Dominique Orban]

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.

[Peter Simon]

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