# Linear system problems

**URL:** <https://discourse.julialang.org/t/linear-system-problems/54519>\
**Category:** GPU\
**Created:** [February 3, 2021, 8:50am UTC](https://discourse.julialang.org/t/linear-system-problems/54519 "2021-02-03T08:50:45Z")\
**Posts on this page:** 16\
**Page:** 2

<div class="post-metadata">

**Author:** ![Ribeiro](https://avatars.discourse-cdn.com/v4/letter/r/d9b06d/32.png) [@Ribeiro](https://discourse.julialang.org/u/Ribeiro)\
**Post date:** [February 3, 2021, 7:59pm UTC](https://discourse.julialang.org/t/linear-system-problems/54519/21 "2021-02-03T19:59:58Z")

</div>

> [@fph](#):
>
> It seems very peculiar that `inv(A) * b` is more accurate than `A \ b`

That seems to have been a mistake on my side. A\b is only inaccurate on the GPU for Float32 (maybe on the CPU for Float32 too. I haven’t tested that).

> [@rveltz](#):
>
> change `A .= I .+ A ./N` and everything should work, debugging will be easier. Now, if your matrix is really rand(N,N), i’d rather do the math than the computations…

Sorry, I don’t understand. That command doesn’t work. The random matrices I’m using in these examples are just placeholders to demonstrate my issue. My real matrix is quite different.

> [@ctkelley](#):
>
> This looks like a sign that your matrix is very ill-conditioned. Have you tried running `cond` to look into this?

`cond(rand(1000,1000)) > 20000`  
Again, I’m using these random matrices as tests that I can easily share. Not sure what that number means.

> [@rveltz](#):
>
> you should print the residual of `x2=linsolve(A,b)`

I looked closely and linsolve is complaining that it didn’t converge, but if it didn’t converge after all that time, then it’s not nearly as good as A\b.

Thanks, everyone, for chiming in!

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [February 3, 2021, 8:06pm UTC](https://discourse.julialang.org/t/linear-system-problems/54519/22 "2021-02-03T20:06:48Z")

</div>

> [@Ribeiro](#):
>
> That command doesn’t work. The random matrices I’m using in these examples are just placeholders to demonstrate my issue.

Exactly but `A=rand(N,N)` is a bad choice for gmres. Try `A = I + rand(N,N) ./N` instead.

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [February 3, 2021, 8:13pm UTC](https://discourse.julialang.org/t/linear-system-problems/54519/23 "2021-02-03T20:13:20Z")

</div>

> [@Ribeiro](#):
>
> cond(rand(1000,1000)) \> 20000

That tells me that you’re likely to lose 4–5 figures in the computation so it would not be a surprise to see `norm(A*x - b, Inf)` as large as 10^{-12}. If you use the L^2 norm, the you could also scale by the square root of the dimension, get up to 10^{-10}, which is what @BLI got in an post upstream.

---

<div class="post-metadata">

**Author:** ![Ribeiro](https://avatars.discourse-cdn.com/v4/letter/r/d9b06d/32.png) [@Ribeiro](https://discourse.julialang.org/u/Ribeiro)\
**Post date:** [February 3, 2021, 8:50pm UTC](https://discourse.julialang.org/t/linear-system-problems/54519/24 "2021-02-03T20:50:25Z")

</div>

> [@rveltz](#):
>
> Exactly but `A=rand(N,N)` is a bad choice for gmres. Try `A = I + rand(N,N) ./N` instead.

I see. You’re right, the matrix itself plays a huge role, of course.

```julia
N=10000;
A2=I+rand(N,N)/N;
x=rand(N);
b2=A2*x;
@time x1=A2\b2; # 1.4 sec, 1e-11 accuracy
@time x2=linsolve(A2,b2); # 0.11 sec, 1e-11 accuracy
@time x3=gmres(A2,b2); # 0.005 sec, 1e-5 accuracy

```

This is on the CPU. So GMRES is really fast, but (using default parameters), not nearly as accurate as the others (probably good enough, though) and linsolve is 10x faster than \ with the same accuracy (for this case), which makes sense!  
So I tested all of the options on my real matrix (a small one, though) and got this:

```julia
N=1914
invert and multiply: 0.06 sec
\: 0.03 sec, 0.07 sec GPU
linsolve: 0.004 sec, 0.04 GPU
gmres: 0.002 sec

```

All of them were equally accurate (1.5e-3, which is ok). So on a CPU, for my specific problem it looks like GMRES works really well.  
So now I have to nag you about the GPU again. linsolve is clearly not running well on the GPU (I still get the warnings) and gmres doesn’t run at all. You said you use them on GPUs. What’s your trick?  
(Really appreciate the help!)

> [@ctkelley](#):
>
> That tells me that you’re likely to lose 4–5 figures in the computation so it would not be a surprise to see `norm(A*x - b, Inf)` as large as 10−1210^{-12} . If you use the L2L^2 norm, the you could also scale by the square root of the dimension, get up to 10−1010^{-10} , which is what @BLI got in an post upstream.

Sorry, but that is going way over my head. I looked at my real (smallish) case and on my 1914x1914 matrix I get cond=2.06. What does that tell you? Thanks a lot!

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [February 3, 2021, 8:57pm UTC](https://discourse.julialang.org/t/linear-system-problems/54519/25 "2021-02-03T20:57:32Z")

</div>

> [@Ribeiro](#):
>
> linsolve is clearly not running well on the GPU (I still get the warnings)

did you do `linsolve(x->Agpu*x, bgpu)`?

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [February 3, 2021, 8:58pm UTC](https://discourse.julialang.org/t/linear-system-problems/54519/26 "2021-02-03T20:58:55Z")

</div>

> [@Ribeiro](#):
>
> GMRES is really fast, but (using default parameters), not nearly as accurate as the others (probably good enough, though)

linsolve is a gmres… you need to tune the options of IterativeSolvers to match KrylovKit

---

<div class="post-metadata">

**Author:** ![Ribeiro](https://avatars.discourse-cdn.com/v4/letter/r/d9b06d/32.png) [@Ribeiro](https://discourse.julialang.org/u/Ribeiro)\
**Post date:** [February 3, 2021, 9:02pm UTC](https://discourse.julialang.org/t/linear-system-problems/54519/27 "2021-02-03T21:02:38Z")

</div>

> [@rveltz](#):
>
> did you do `linsolve(x->Agpu*x, bgpu)` ?

That is 4x slower than `linsolve(A,b)`;

> [@rveltz](#):
>
> linsolve is a gmres… you need to tune the options of IterativeSolvers to match KrylovKit

Oh. My bad. Thanks for the info.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [February 3, 2021, 9:47pm UTC](https://discourse.julialang.org/t/linear-system-problems/54519/28 "2021-02-03T21:47:30Z")

</div>

> [@Ribeiro](#):
>
> That is 4x slower than `linsolve(A,b)` ;

but do you get the warnings?

---

<div class="post-metadata">

**Author:** ![Ribeiro](https://avatars.discourse-cdn.com/v4/letter/r/d9b06d/32.png) [@Ribeiro](https://discourse.julialang.org/u/Ribeiro)\
**Post date:** [February 3, 2021, 10:00pm UTC](https://discourse.julialang.org/t/linear-system-problems/54519/29 "2021-02-03T22:00:39Z")

</div>

> [@rveltz](#):
>
> but do you get the warnings?

Yes. You don’t?

```julia
A=CUDA.rand(3,3).+1;
x=CUDA.rand(3);
b=A*x;
x1=linsolve(z->A*z, b)
┌ Warning: Performing scalar operations on GPU arrays: This is very slow, consider disallowing these operations with `allowscalar(false)`
└ @ GPUArrays ~/.julia/packages/GPUArrays/WV76E/src/host/indexing.jl:43

```

You’re able to run linsolve if you set `CUDA.allowscalar(false)`?

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [February 3, 2021, 10:39pm UTC](https://discourse.julialang.org/t/linear-system-problems/54519/30 "2021-02-03T22:39:10Z")

</div>

> [@Ribeiro](#):
>
> I get cond=2.06. What does that tell you? Thanks a lot!

That condition number is completely benign and a dimension of 2000 is not going to do much damage. Is this the problem where `A\b` failed to give a decent residual on the GPU in single precision? If so, that is a real mystery. There’s theory that says that should not happen. Does the same thing happen in Float32 on your CPU?

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [February 3, 2021, 11:57pm UTC](https://discourse.julialang.org/t/linear-system-problems/54519/31 "2021-02-03T23:57:17Z")

</div>

> [@Ribeiro](#):
>
> You’re able to run linsolve if you set `CUDA.allowscalar(false)` ?

Yes. Something is off and I dont have a gpu to test.

---

<div class="post-metadata">

**Author:** ![Ribeiro](https://avatars.discourse-cdn.com/v4/letter/r/d9b06d/32.png) [@Ribeiro](https://discourse.julialang.org/u/Ribeiro)\
**Post date:** [February 4, 2021, 9:12am UTC](https://discourse.julialang.org/t/linear-system-problems/54519/32 "2021-02-04T09:12:00Z")

</div>

> [@ctkelley](#):
>
> That condition number is completely benign and a dimension of 2000 is not going to do much damage. Is this the problem where `A\b` failed to give a decent residual on the GPU in single precision? If so, that is a real mystery. There’s theory that says that should not happen. Does the same thing happen in Float32 on your CPU?

I just tested it and this matrix (cond=2) is fine on the GPU using Float32. So the accuracy problem was strongly linked to the poorly conditioned random matrices I was using for testing. Thanks!

> [@rveltz](#):
>
> Yes. Something is off and I dont have a gpu to test.

Huh. Odd. There is a high probability that I’m doing something stupid =].  
If you at some point get to test on a GPU and see that this:

```julia
using KrylovKit
using CUDA
CUDA.allowscalar(false);
A=CUDA.rand(3,3).+1;
x=CUDA.rand(3);
b=A*x;
x1=linsolve(z->A*z, b)

```

works, then please let me know. Maybe you load some other stuff that is needed or have some CUDA settings that I don’t (or linsolve settings!).  
Thanks a lot!

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [February 4, 2021, 9:47am UTC](https://discourse.julialang.org/t/linear-system-problems/54519/33 "2021-02-04T09:47:31Z")

</div>

> [@Ribeiro](#):
>
> `x1=linsolve(z->A*z, b)`

You need to know that KK works for very general vectors, so you have to define [this](https://jutho.github.io/KrylovKit.jl/stable/) for it to work. I am surprised that it does not work here though. Maybe @juthohaegeman has an idea…

The reason I am surprised is because it used to be required to export the right `mul, axpy` but not anymore as you can see [here](https://rveltz.github.io/BifurcationKit.jl/dev/tutorials2b/#The-Swift-Hohenberg-equation-(non-local)-on-the-GPU-(Advanced)-1). In this tutorial, you are shown how to solve linear equations on GPU (among other things).

Ok, the culprit is `mul!` at `src/linsolve/gmres.jl:11`. You can fix it with (with type piracy 🤭?)

```julia
using KrylovKit, LinearAlgebra
using CUDA
CUDA.allowscalar(false);
A=CUDA.rand(3,3).+1;
x=CUDA.rand(3);
b=A*x;

LinearAlgebra.mul!(x::CuArray{T,1},y::CuArray{T,1},α) where T = x .= α .* y

julia> x1=linsolve(z->A*z, b, b)
(Float32[0.72570944, 0.23375773, 0.08379406], ConvergenceInfo: one converged value after 3 iterations and 11 applications of the linear map;
norms of residuals are given by (0.0f0,).
)

```

---

<div class="post-metadata">

**Author:** ![Ribeiro](https://avatars.discourse-cdn.com/v4/letter/r/d9b06d/32.png) [@Ribeiro](https://discourse.julialang.org/u/Ribeiro)\
**Post date:** [February 4, 2021, 10:31am UTC](https://discourse.julialang.org/t/linear-system-problems/54519/34 "2021-02-04T10:31:17Z")

</div>

@rveltz Thanks for all the tips. A lot of it goes over my head, but it helps!  
Your hacks work and linsolve now runs on the GPU! It is slightly slower and less reliable than \ though.  
On my real problem, \ took 0.07 sec, while linsolve took 0.1 sec, both with good accuracy.  
On the big random matrix, linsolve diverged, while \ worked great.  
So it seems that \ still beats the other options on the GPU, which is surprising to me, but fine for now.  
Thanks again!

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [February 4, 2021, 10:47am UTC](https://discourse.julialang.org/t/linear-system-problems/54519/35 "2021-02-04T10:47:00Z")

</div>

> [@Ribeiro](#):
>
> Your hacks work and linsolve now runs on the GPU! It is slightly slower and less reliable than \ though.

Again, I disagree! Iterative methods should work best for `A = I + rand(N,N) ./ N` with N large enough.

This thing to remember is that you cannot use gmres as `\`. It is not a black box: you need to know some properties of A (perturbation of identity, preconditioner) for it to work well.

---

<div class="post-metadata">

**Author:** ![Ribeiro](https://avatars.discourse-cdn.com/v4/letter/r/d9b06d/32.png) [@Ribeiro](https://discourse.julialang.org/u/Ribeiro)\
**Post date:** [February 4, 2021, 10:57am UTC](https://discourse.julialang.org/t/linear-system-problems/54519/36 "2021-02-04T10:57:31Z")

</div>

@rveltz I believe you. I’m not versed in GMRES yet, so I am indeed using it as a black box. I might have to revisit this in the near future, but will need to educate myself a bit before conducting more tests.  
Thanks!

[Previous page](https://discourse.julialang.org/t/linear-system-problems/54519.md?page=1)