# Solve the System of Equations (J' J + μI) x = J'v

**URL:** <https://discourse.julialang.org/t/solve-the-system-of-equations-j-j-i-x-jv/42766>\
**Category:** Performance\
**Tags:** linearalgebra\
**Created:** [July 9, 2020, 3:01am UTC](https://discourse.julialang.org/t/solve-the-system-of-equations-j-j-i-x-jv/42766 "2020-07-09T03:01:14Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Norman](https://avatars.discourse-cdn.com/v4/letter/n/97f17d/32.png) [@Norman](https://discourse.julialang.org/u/Norman)\
**Post date:** [July 9, 2020, 3:01am UTC](https://discourse.julialang.org/t/solve-the-system-of-equations-j-j-i-x-jv/42766/1 "2020-07-09T03:01:14Z")

</div>

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)

```

---

<div class="post-metadata">

**Author:** ![lucas711642](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lucas711642/32/12051_2.png) [@lucas711642](https://discourse.julialang.org/u/lucas711642)\
**Post date:** [July 9, 2020, 3:45am UTC](https://discourse.julialang.org/t/solve-the-system-of-equations-j-j-i-x-jv/42766/2 "2020-07-09T03:45:13Z")

</div>

If μ is positive, then you can use the `shift` keyword of `cholesky`, replacing

```julia
(J' * J + spdiagm(0 => fill(μ, n))) \ B

```

by

```julia
cholesky(J' * J; shift = μ) \ B

```

Note that you don’t need to construct the diagonal matrix in your code, and you could have used `H = J' * J + μ * I` instead.

---

<div class="post-metadata">

**Author:** ![ettersi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ettersi/32/6829_2.png) [@ettersi](https://discourse.julialang.org/u/ettersi)\
**Post date:** [July 9, 2020, 4:01am UTC](https://discourse.julialang.org/t/solve-the-system-of-equations-j-j-i-x-jv/42766/3 "2020-07-09T04:01:46Z")

</div>

Your linear system are the normal equations of the least-squares problem

\min\_x \|\hat J x - \hat v\|\_2

where

\hat J = \begin{pmatrix} J \\ \sqrt{\mu} \, I\_n \end{pmatrix} ,\qquad \hat v = \begin{pmatrix} v \\ 0\_n \end{pmatrix}

with I\_n, 0\_n the n \times n identity and zero matrices, respectively. It is usually recommended to tackle the original least squares problem rather than the normal equations since the normal equations are much more ill-conditioned.

Since your least squares is very large and sparse, I recommend you try solving it using iterative methods like [LSMR](https://juliamath.github.io/IterativeSolvers.jl/dev/linear_systems/lsmr) or [LSQR](https://juliamath.github.io/IterativeSolvers.jl/dev/linear_systems/lsqr) which are implemented in IterativeSolvers.jl.

---

<div class="post-metadata">

**Author:** ![Joris\_Pinkse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joris_pinkse/32/216398_2.png) [@Joris\_Pinkse](https://discourse.julialang.org/u/Joris_Pinkse)\
**Post date:** [July 9, 2020, 4:04am UTC](https://discourse.julialang.org/t/solve-the-system-of-equations-j-j-i-x-jv/42766/4 "2020-07-09T04:04:46Z")

</div>

Julia compiles things before it runs. This means that if you only run some code once, then it can take longer. So only time on the second run (add another main call) to make the comparison fair. Julia speed is moreover sensitive to the way you programs things. I don’t know enough about the sparse array functions to comment. When I run your code, Julia appears to only be using one thread, whereas Matlab may be using more. There is a limit on the number of threads that OpenBLAS uses (8?), which is probably not a constraint here, but which you can change if you compile from source. There are other libraries that you could use like MKL. Finally, look into PackageCompiler: it’s my best friend.

---

<div class="post-metadata">

**Author:** ![jlperla](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlperla/32/34332_2.png) [@jlperla](https://discourse.julialang.org/u/jlperla)\
**Post date:** [July 9, 2020, 4:50am UTC](https://discourse.julialang.org/t/solve-the-system-of-equations-j-j-i-x-jv/42766/5 "2020-07-09T04:50:47Z")

</div>

> [@ettersi](#):
>
> Since your least squares is very large and sparse, I recommend you try solving it using iterative methods like [LSMR](https://juliamath.github.io/IterativeSolvers.jl/dev/linear_systems/lsmr) or [LSQR](https://juliamath.github.io/IterativeSolvers.jl/dev/linear_systems/lsqr) which are implemented in IterativeSolvers.jl.

yes. To reiterate, those algorithms implement the L2 regularized least squares directly, which I believe is the exact form of your problem.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [July 9, 2020, 9:31am UTC](https://discourse.julialang.org/t/solve-the-system-of-equations-j-j-i-x-jv/42766/6 "2020-07-09T09:31:40Z")

</div>

> [@ettersi](#):
>
> Your linear system are the normal equations of the least-squares problem
> 
> minx∥^Jx−^v∥2 \min\_x |\hat J x - \hat v|\_2
> 
> where ^J\hat J denotes JJ extended by an extra row of √μ\sqrt{\mu} and ^v\hat v denotes vv extended by an extra 0, i.e.
> 
> ^J=(J√μ),^v=(v0).

I don’t think this reformulation is correct. The normal equations for the least squares problem you wrote are (J' J + \mu \textbf{1}) x = J' v. Notice the \mu \textbf{1} not the original \mu I.

Edit: the equations actually describe the optimality conditions of minimizing || J x - v ||\_2^2 + \mu ||x||\_2^2.  
Edit 2: if \mu is positive then J' J + \mu I is positive definite, so you can use `IterativeSolvers.cg`. You may not need to form J' J + \mu I at all, just use the matrix-free API from [https://juliamath.github.io/IterativeSolvers.jl/dev/getting\_started/#matrixfree-1](https://juliamath.github.io/IterativeSolvers.jl/dev/getting_started/#matrixfree-1).

---

<div class="post-metadata">

**Author:** ![ettersi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ettersi/32/6829_2.png) [@ettersi](https://discourse.julialang.org/u/ettersi)\
**Post date:** [July 10, 2020, 2:36pm UTC](https://discourse.julialang.org/t/solve-the-system-of-equations-j-j-i-x-jv/42766/7 "2020-07-10T14:36:44Z")

</div>

Correct, my original formulation was wrong. I’ve updated my answer.  
@Norman’s solution is correct anyway since he used the \lambda parameter of `lsmr` rather than the augmented linear system.

Applying CG to the normal equations does not solve the ill-conditioning problem since the matvec J^TJ \, v wipes out information in v regardless of how you evaluate the product. This is the magic of `lsqr`, it implicitly runs CG on J^TJ + \mu I, but it does it in a clever way to avoid this loss of information.