I have a square sparse matrix
J of size
n can be around 1 million.
μ is a scalar.
I is the identity matrix.
v is a vector.
The system to be solved is
(J'J + μI) x = J'v
What I have now is the following:
H = J'*J + spdiagm(0=>fill(μ, n)) B = J'*v x = H\B
J is of type
n can be pretty large.
Is there some better way to solve this system?
Here is a coarse comparison between Julia and Matlab on my computer.
I am not sure whether this comparison is fair. But, I am very curious on how I could improve the Julia code to make things faster.
Note that although I am using a tridiagonal matrix here, the matrix
J may take other form in actual usage. It is a Jacobian matrix for a Fischer-Burmeister function.
using LinearAlgebra, SparseArrays function main() n = 1000000 J = spdiagm(0=>[9;fill(17,n-2);9], -1=>fill(-8,n-1), 1=>fill(-8,n-1)) v = fill(0.05,n) H = J'*J + spdiagm(0=>fill(0.05, n)) B = J'*v @time x = H\B return x end main() % 0.394470 seconds (48 allocations: 534.059 MiB, 1.99% gc time)
n = 1000000; J = spdiags([-8*ones(n,1) [9;17*ones(n-2,1);9] -8*ones(n,1)], [-1 0 1], n, n); v = 0.05*ones(n,1); H = J'*J + spdiags(0.05*ones(n,1), 0, n, n); B = J'*v; tic; x = H\B; toc; % Elapsed time is 0.101177 seconds.
Thanks for your time!
The suggestion from @ettersi really makes a difference.
using LinearAlgebra, SparseArrays using IterativeSolvers function main2() n = 1000000 J = spdiagm(0=>[9;fill(17,n-2);9], -1=>fill(-8,n-1), 1=>fill(-8,n-1)) v = fill(0.05,n) @time x = lsmr(J, v, λ=sqrt(0.05)) return x end main2() % 0.033991 seconds (34 allocations: 61.037 MiB)