I have a square sparse matrix `J`

of size `n`

* `n`

where `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 `SparseMatrixCSC`

and `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.

Julia:

```
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)
```

Matlab:

```
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)
```