# How to solve this Ax=b faster?
**URL:** https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693
**Category:** Numerics
**Tags:** linearalgebra, linearsolve
**Created:** [May 26, 2022, 5:34am UTC](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693 "2022-05-26T05:34:42Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@Minimum-Pollution-96](https://discourse.julialang.org/u/Minimum-Pollution-96)
#### Post date: [May 26, 2022, 5:34am UTC](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693/1 "2022-05-26T05:34:42Z")
I am trying to solve a sparse stiffness matrix (positive definite) of size `1986360x1986360 SparseMatrixCSC{Float64, Int64} with 18327153 stored entries`.
Here’s a sample code.
```julia
using SparseArrays
using LinearAlgebra
A = sprandn(1986360, 1986360, 2e-6).*20;
A = A+A'+100I
b = zeros(size(a,1));
b[[1:331060;].*6] .= 200;
x = A\b
```
The problem is that the last step takes so much time as it uses only one core (have never actually finished it). I tried running it via julia REPL normally as well as by launching julia with `julia --threads=16`.
I also tried solving it as `x=cholesky(A)\b` but that also uses only one core. It also throws `SuiteSparse.CHOLMOD.CHOLMODException("out of memory")` after sometime even though 23 gb of memory(out of 32) is available.
Is there anyway to make julia use all the cores? Also, I am trying solve such a large matrix for the first time. So, please do suggest if there is something I can do to solve this faster.
---
### Author:  [@Jae-Mo\_Lihm](https://discourse.julialang.org/u/Jae-Mo_Lihm)
#### Post date: [May 26, 2022, 8:19am UTC](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693/2 "2022-05-26T08:19:12Z")
Have you tried iterative solvers?
There is Krylov.jl and many others.
[https://juliasmoothoptimizers.github.io/Krylov.jl/stable/solvers/spd/](https://juliasmoothoptimizers.github.io/Krylov.jl/stable/solvers/spd/)
---
### Author:  [@greatpet](https://discourse.julialang.org/u/greatpet)
#### Post date: [May 26, 2022, 8:29am UTC](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693/3 "2022-05-26T08:29:37Z")
Random sparse systems should be solved by iterative solvers, since the default direct solver causes fill-in during Gaussian elimination and blows up the memory usage.
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [May 26, 2022, 8:35am UTC](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693/4 "2022-05-26T08:35:37Z")
> [@Minimum-Pollution-96](#):
>
> `[1:331060;].*6`
Side note: Don’t do this. It allocates an array - no! _two_ arrays - needlessly. It is of course not an effect of any significance here, but it’s a good idea to mind unnecessary allocations in general.
Instead, do `(1:331060) .* 6` or `6:6:(331060*6)`.
---
### Author:  [@Minimum-Pollution-96](https://discourse.julialang.org/u/Minimum-Pollution-96)
#### Post date: [May 26, 2022, 9:01am UTC](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693/5 "2022-05-26T09:01:44Z")
Thank you! I’ll check this out.
---
### Author:  [@Minimum-Pollution-96](https://discourse.julialang.org/u/Minimum-Pollution-96)
#### Post date: [May 26, 2022, 9:01am UTC](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693/6 "2022-05-26T09:01:56Z")
Thanks. Noted.
---
### Author:  [@Minimum-Pollution-96](https://discourse.julialang.org/u/Minimum-Pollution-96)
#### Post date: [May 26, 2022, 9:05am UTC](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693/7 "2022-05-26T09:05:40Z")
Could you explain this a little more? I thought that sparse systems have separate solvers that act on them. [Sparse Linear Algebra · The Julia Language](https://docs.julialang.org/en/v1/stdlib/SuiteSparse/)
But I do face out of memory error while trying `cholesky(A)`.
---
### Author:  [@greatpet](https://discourse.julialang.org/u/greatpet)
#### Post date: [May 26, 2022, 9:09am UTC](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693/8 "2022-05-26T09:09:27Z")
Sparse matrices encountered in many real-world problems have some structures, which makes it possible to do Gaussian elimination without turning the matrix into a dense matrix in intermediate steps (called “fill-in”). If your sparse matrix is completely random, then Gaussian elimination doesn’t work well, and you should switch to a package that implements iterative methods which are insensitive to the fill-in problem.
---
### Author:  [@Andrea\_Vigliotti](https://discourse.julialang.org/u/Andrea_Vigliotti)
#### Post date: [May 26, 2022, 9:57am UTC](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693/9 "2022-05-26T09:57:06Z")
you might want to look at this [reply](https://discourse.julialang.org/t/sparse-solve-with-sparse-rhs/81595/3) by @stevengj, I haven’t tried it yet myself, but he mentions two parallel solvers [MUMPS.jl](https://github.com/JuliaSmoothOptimizers/MUMPS.jl) and [Pardiso.jl](https://github.com/JuliaSparse/Pardiso.jl) that might be of help
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 26, 2022, 12:32pm UTC](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693/10 "2022-05-26T12:32:13Z")
Just use LinearSolve.jl
[http://linearsolve.sciml.ai/dev/](http://linearsolve.sciml.ai/dev/)
And then run the gambit.
```julia
for alg in [KrylovJL_GMRES(), KLUFactorization(), KrylovKitJL_GMRES(), IterativeSolversJL_GMRES(), UMFPACKFactorization(), MKLPardisoFactorize(), MKLPardisoIterate()]
sol = solve(prob,alg)
end
```
Likely MKLPardiso is best without a preconditioner. For the GMRES ones, you may want to define a preconditioner like:
[http://linearsolve.sciml.ai/dev/basics/Preconditioners/](http://linearsolve.sciml.ai/dev/basics/Preconditioners/)
---
### Author:  [@Minimum-Pollution-96](https://discourse.julialang.org/u/Minimum-Pollution-96)
#### Post date: [May 27, 2022, 10:17am UTC](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693/11 "2022-05-27T10:17:26Z")
> [@Andrea\_Vigliotti](#):
>
> you might want to look at this [reply](https://discourse.julialang.org/t/sparse-solve-with-sparse-rhs/81595/3) by @stevengj, I haven’t tried it yet myself, but he mentions two parallel solvers [MUMPS.jl](https://github.com/JuliaSmoothOptimizers/MUMPS.jl) and [Pardiso.jl](https://github.com/JuliaSparse/Pardiso.jl) that might be of help
It seems like Pardiso is not free for academic users anymore. Lets me sign up but received a mail saying academic license discontinued from Jan 2022.
> [@ChrisRackauckas](#):
>
> Just use LinearSolve.jl
>
> [http://linearsolve.sciml.ai/dev/](http://linearsolve.sciml.ai/dev/)
>
> And then run the gambit.
>
> ```julia
> for alg in [KrylovJL_GMRES(), KLUFactorization(), KrylovKitJL_GMRES(), IterativeSolversJL_GMRES(), UMFPACKFactorization(), MKLPardisoFactorize(), MKLPardisoIterate()]
> sol = solve(prob,alg)
> end
>
> ```
>
> Likely MKLPardiso is best without a preconditioner. For the GMRES ones, you may want to define a preconditioner like:
>
> [http://linearsolve.sciml.ai/dev/basics/Preconditioners/](http://linearsolve.sciml.ai/dev/basics/Preconditioners/)
Thank you. I was able to solve it with Krylov. But I’ll check this out as well!
---
### Author:  [@photor](https://discourse.julialang.org/u/photor)
#### Post date: [May 27, 2022, 10:57am UTC](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693/12 "2022-05-27T10:57:37Z")
Learned an important point from you and @Jae-Mo_Lihm
---
### Author:  [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [May 27, 2022, 11:18am UTC](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693/13 "2022-05-27T11:18:01Z")
If this problem is occuring as part of an optimization problem, I can recommend COSMO.jl
---
### Author:  [@Andrea\_Vigliotti](https://discourse.julialang.org/u/Andrea_Vigliotti)
#### Post date: [May 28, 2022, 7:52am UTC](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693/14 "2022-05-28T07:52:28Z")
I tried LinearSolve on a 3D Finite Element Problem with 68’170 elements and 82’424 nodes of the type
\mathbf{A}\mathbf{u} = \mathbf{r}
where \mathbf{A} is a sparse SPD matrix with 162’372 rows with 10’889’190 non zero elements (on average 67 non zero entries per row). I tried the following
```julia
using using LinearSolve
prob = LinearProblem(A,r);
for alg in [KrylovJL_GMRES(), KLUFactorization(), KrylovKitJL_GMRES(), IterativeSolversJL_GMRES(), UMFPACKFactorization()]
@time sol = LinearSolve.solve(prob,alg)
end
180.594400 seconds (2.56 M allocations: 2.356 GiB, 3.02% gc time, 0.45% compilation time)
217.641759 seconds (4.28 M allocations: 3.479 GiB, 0.43% gc time)
181.348235 seconds (244.91 k allocations: 2.233 GiB, 3.01% gc time, 0.05% compilation time)
954.111518 seconds (4.03 M allocations: 450.389 MiB, 0.05% gc time, 0.17% compilation time)
13.303307 seconds (121.40 k allocations: 2.719 GiB, 5.24% gc time, 0.96% compilation time)
```
the MKL solvers gave and `UndefVarError` so I couldnt try them.
I also tried `Krylov.jl`, and got
```julia
using Krylov
@time x, stats = Krylov.cg(A,r)
80.857184 seconds (1.28 M allocations: 77.709 MiB, 0.88% compilation time)
```
if now I compare these performances with the standard `LinearAlgebra.jl` package, this is what I get
```julia
@time A\r;
3.062161 seconds (661 allocations: 1.390 GiB, 1.07% gc time, 0.13% compilation time)
```
it seems that the `LinearAlgebra` solver, on this problem, outperforms all the others, should that be the case or am I not getting the best speed from the other solvers for some reason? or is it the size of the problem that influences the performances.
(please, do repost this as a new thread if think it is appropriate)
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 28, 2022, 11:42am UTC](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693/15 "2022-05-28T11:42:13Z")
> [@Andrea\_Vigliotti](#):
>
> it seems that the `LinearAlgebra` solver, on this problem, outperforms all the others, should that be the case or am I not getting the best speed from the other solvers for some reason? or is it the size of the problem that influences the performances.
And it’s not just compile time? The LinearAlgebra default is just `UMFPACKFactorization` so it should be the same time. Do you have a good way for me to reproduce this?
For the Krylov methods, did you try adding any preconditioners? Algebraic multigrid is probably a good idea here.
---
### Author:  [@Andrea\_Vigliotti](https://discourse.julialang.org/u/Andrea_Vigliotti)
#### Post date: [May 28, 2022, 4:07pm UTC](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693/16 "2022-05-28T16:07:28Z")
The matrix \mathbf{A} is the result of a FE mesh with hexahedral elements of a solid mechanics problem, with 3 degrees of freedom (DoFs) per node. I have Dirichlet boundary conditions on the boundary and free internal nodes, thus the problem is constructed as follows
\mathbf{A} = \mathbf{Kt}[\textrm{iDoFs},\textrm{iDoFs}]\\ \mathbf{r} = \mathbf{Kt}[\textrm{iDoFs},\textrm{bDoFs}]\mathbf{s}[\textrm{bDoFs}]
where \mathbf{Kt} is the stiffness matrix of the entire domain, \textrm{bDoFs} are the indices of the boundary DoFs and \textrm{iDoFs} are the indices of the internal DoFs, \mathbf{s} is the entire DoFs vector. \mathbf{A} is Hermitian, with `size(A)=(162372, 162372)`
Depending on the position of the element, whether internal or on the boundary, each node can be surrounded from 1 up to 8 elements, which decides the number of of non zero entries per each row of the matrix, all diagonal entries are non zero. It is a quite structured matrix, which is difficult to obtain by random generation, if it helps can share the matrix and the vectors as `JLD2` files.
This is what I tried
```julia
u = @btime A\r
@show norm(r-A*u);
3.843 s (52 allocations: 1.39 GiB)
norm(r - A * u) = 5.766109468498555e-12
```
```julia
u = @btime lu(A)\r
@show norm(r-A*u);
10.340 s (72 allocations: 2.39 GiB)
norm(r - A * u) = 4.227052608619862e-12
```
then with `LinearSolve.jl`
```julia
prob = LinearProblem(A,r);
sol = @btime LinearSolve.solve(prob, UMFPACKFactorization())
@show norm(r-A*sol.u);
10.984 s (94 allocations: 2.71 GiB)
norm(r - A * sol.u) = 4.227052608619862e-12
```
then using precodintioners from `AlgebraicMultigrid.jl` with Krylov methods I get
```julia
m = AlgebraicMultigrid.ruge_stuben(A)
p = AlgebraicMultigrid.aspreconditioner(m)
sol = @btime LinearSolve.solve(prob, LinearSolve.KrylovKitJL_GMRES(), Pl=p)
@show norm(r-A*sol.u);
51.581 s (603 allocations: 471.36 MiB)
norm(r - A * sol.u) = 7.902120848944424e-5
```
and
```julia
m = AlgebraicMultigrid.smoothed_aggregation(A)
p = AlgebraicMultigrid.aspreconditioner(m)
sol = @btime LinearSolve.solve(prob, LinearSolve.KrylovKitJL_GMRES(), Pl=p)
@show norm(r-A*sol.u);
46.302 s (613 allocations: 477.56 MiB)
norm(r - A * sol.u) = 0.0001005812834686952
```
all of these have run on `Jupyter-lab` notebook with a single thread. Different runs of `@btime` differ of up to 1 sec.
Please do let me know if you have any question or need more information.
side note: `lu` and `UMFPACKFactorization()` terminate on the same solution:
```julia
sol = LinearSolve.solve(prob, UMFPACKFactorization())
u = lu(A)\r
@show norm(u-sol.u);
norm(u - sol.u) = 0.0
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 28, 2022, 4:39pm UTC](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693/17 "2022-05-28T16:39:39Z")
What LinearSolve.jl is doing is equivalent to the `lu(A)\b` route, which will default to UMFPACK in Base Julia, but I did not know that there was a specialization of `\` that was faster. Interesting, I wonder what shortcut it takes. @rayegun could you take a look at this?
---
### Author:  [@Andrea\_Vigliotti](https://discourse.julialang.org/u/Andrea_Vigliotti)
#### Post date: [May 28, 2022, 5:15pm UTC](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693/18 "2022-05-28T17:15:57Z")
this is what I get from `@which A\r` :
`~/julia/julia-1.6.4/share/julia/stdlib/v1.6/SparseArrays/src/linalg.jl:1551`
---
### Author:  [@rayegun](https://discourse.julialang.org/u/rayegun)
#### Post date: [May 28, 2022, 6:43pm UTC](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693/19 "2022-05-28T18:43:44Z")
I’ll take a look
---
### Author:  [@rayegun](https://discourse.julialang.org/u/rayegun)
#### Post date: [May 28, 2022, 6:48pm UTC](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693/20 "2022-05-28T18:48:48Z")
Unfortunately that’s the dispatch method, which doesn’t give much info, you said A is Hermitian?
Can you send me the JLD2 file? If you need to send it privately you can transmit it to kimmerer@mit.edu somehow.
It looks like it’s calling CHOLMOD, by dispatching on `ishermitian`. Should we duplicate this codepath @ChrisRackauckas? I’m not even sure if we wrap CHOLMOD very well in LinearSolve.jl yet. `julia-1.7.2+0~x64/share/julia/stdlib/v1.7/SuiteSparse/src/cholmod.jl:1573`. I am working on familiarizing myself with the CHOLMOD codebase anyway right now.
[Next page](https://discourse.julialang.org/t/how-to-solve-this-ax-b-faster/81693.md?page=2)