On Wednesday, July 30, 2014 10:08:45 PM UTC-5, 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
>
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
