# CUDA Matrix inverse

**URL:** <https://discourse.julialang.org/t/cuda-matrix-inverse/53341>\
**Category:** GPU\
**Tags:** cuda, linearalgebra\
**Created:** [January 14, 2021, 4:23pm UTC](https://discourse.julialang.org/t/cuda-matrix-inverse/53341 "2021-01-14T16:23:11Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [January 14, 2021, 4:23pm UTC](https://discourse.julialang.org/t/cuda-matrix-inverse/53341/1 "2021-01-14T16:23:11Z")

</div>

Hello!

I try to inverse large matrix on GPU and try to use this code:

```julia
using CUDA, LinearAlgebra
function cuinv(m::Matrix{T}) where T
    if size(m, 1) != size(m, 2) throw(ArgumentError("Matrix not square.")) end
    A = CuArray(m)
    B = CuArray(Matrix{T}(I(size(A,1))))
    A, ipiv = CUDA.CUSOLVER.getrf!(A)
    Matrix{T}(CUDA.CUSOLVER.getrs!('N', A, ipiv, B))
end
m = rand(100, 100)
isapprox(cuinv(m), inv(m))
#true

```

Is it optimal solution?

I used this example for code above:

```julia
function Base.:\(_A, _B)
    A, B = copy(_A), copy(_B)
    A, ipiv = CUDA.CUSOLVER.getrf!(A)
    return CUDA.CUSOLVER.getrs!('N', A, ipiv, B)
end

```

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [January 14, 2021, 4:56pm UTC](https://discourse.julialang.org/t/cuda-matrix-inverse/53341/2 "2021-01-14T16:56:19Z")

</div>

This isn’t optimal because you are taking the inverse of a matrix. Don’t do this. Instead you should use a solver (if you only need the inverse once) or a factorization. This approach will be faster and or significantly more numerically stable.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 14, 2021, 5:16pm UTC](https://discourse.julialang.org/t/cuda-matrix-inverse/53341/3 "2021-01-14T17:16:26Z")

</div>

> [@Oscar\_Smith](#):
>
> This approach will be faster and or significantly more numerically stable.

Computing the matrix inverse explicitly is often wasteful, compared to computing the solution directly for your right-hand side(s), but it is not numerically unstable.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [January 14, 2021, 5:24pm UTC](https://discourse.julialang.org/t/cuda-matrix-inverse/53341/4 "2021-01-14T17:24:42Z")

</div>

[Why Shouldn't I Invert That Matrix?](http://gregorygundersen.com/blog/2020/12/09/matrix-inversion/) shows that inverting and multiplying has about 2x the forward error and 10x the backward error of a LU based solve (in the well conditioned case) and almost infinitely worse error for ill conditioned matrices. Furthermore the empirical results are in line with the theoretical ones.

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [January 14, 2021, 6:24pm UTC](https://discourse.julialang.org/t/cuda-matrix-inverse/53341/5 "2021-01-14T18:24:07Z")

</div>

my task is like this: mult x’ \* inv(m) \* x (where m - not changing, x - many different)

```julia
m = rand(10,10)
x = rand(10, 100)

function f1(m, x)
    lud = lu(m)
    for i = 1:100
        x' / lud
    end
end
function f2(m, x)
    im = inv(m)
    for i = 1:100
        x' * im
    end
end
function f3(m, x)
    lud = lu(m)
    for i = 1:100
        lud \ x
    end
end

@benchmark f1($m, $x)
@benchmark f2($m, $x)
@benchmark f3($m, $x)

```

result:

```julia
julia> @benchmark f1($m, $x)
BenchmarkTools.Trial:
  memory estimate: 794.78 KiB
  allocs estimate: 102
  --------------
  minimum time: 946.100 μs (0.00% GC)
  median time: 1.200 ms (0.00% GC)
  mean time: 1.646 ms (4.78% GC)
  maximum time: 12.848 ms (90.30% GC)
  --------------
  samples: 3037
  evals/sample: 1

julia> @benchmark f2($m, $x)
BenchmarkTools.Trial:
  memory estimate: 799.94 KiB
  allocs estimate: 104
  --------------
  minimum time: 448.300 μs (0.00% GC)
  median time: 560.500 μs (0.00% GC)
  mean time: 913.309 μs (6.89% GC)
  maximum time: 32.566 ms (0.00% GC)
  --------------
  samples: 5468
  evals/sample: 1
julia> 
julia> 

julia> @benchmark f3($m, $x)
BenchmarkTools.Trial:
  memory estimate: 794.78 KiB
  allocs estimate: 102
  --------------
  minimum time: 5.068 ms (0.00% GC)
  median time: 10.902 ms (0.00% GC)
  mean time: 10.238 ms (0.66% GC)
  maximum time: 28.456 ms (0.00% GC)
  --------------
  samples: 489
  evals/sample: 1

julia> 

```

variant 2 is fastest

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 14, 2021, 6:41pm UTC](https://discourse.julialang.org/t/cuda-matrix-inverse/53341/6 "2021-01-14T18:41:53Z")

</div>

> [@Oscar\_Smith](#):
>
> [Why Shouldn’t I Invert That Matrix?](http://gregorygundersen.com/blog/2020/12/09/matrix-inversion/) shows that inverting and multiplying has about 2x the forward error and 10x the backward error of a LU based solve

Being worse by a constant factor doesn’t mean it is “numerically unstable”. (The definition of numerical stability is insensitive to constant factors.) [This article](https://arxiv.org/abs/1201.6035), for example, argues that multiplication by the matrix inverse is ordinarily backwards stable.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 14, 2021, 6:57pm UTC](https://discourse.julialang.org/t/cuda-matrix-inverse/53341/7 "2021-01-14T18:57:20Z")

</div>

> [@PharmCat](#):
>
> my task is like this: mult `x’ * inv(m) * x` (where m - not changing, x - many different)

As a general rule of thumb, computing and re-using an explicit matrix inverse will be more efficient than re-using the LU factors for each solve if you have many more vectors than the dimension of your matrix, assuming the matrices are dense (not sparse or otherwise specially structured). (You should also consider pre-allocating result vectors and using `mul!` to work in-place.) This was the case in your little benchmark, but is it also true of your real application?

That being said, there might be more efficient methods if you are specifically interested in x^T M^{-1} x. For example, if you are computing the _average_ of x^T M^{-1} x for many random x, then there is an analytical formula for the result in terms of the trace of M^{-1}. More generally, the computation of x^T M^{-1} x is closely related to trace estimation. If M is a large sparse symmetric matrix, for example, there are efficient iterative methods called [“Lanczos quadrature”](https://locus.siam.org/doi/abs/10.1137/16M1104974?mobileUi=0) for computing x^T M^{-1} x _without_ computing M^{-1} x.

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [January 14, 2021, 7:02pm UTC](https://discourse.julialang.org/t/cuda-matrix-inverse/53341/8 "2021-01-14T19:02:59Z")

</div>

> [@Oscar\_Smith](#):
>
> [Why Shouldn't I Invert That Matrix?](http://gregorygundersen.com/blog/2020/12/09/matrix-inversion/) shows that inverting and multiplying has about 2x the forward error and 10x the backward error of a LU based solve (in the well conditioned case) and almost infinitely worse error for ill conditioned matrices. Furthermore the empirical results are in line with the theoretical ones.

[We had that discussion already](https://discourse.julialang.org/t/result-changes-depending-on-whether-matrix-is-sparse-or-not/30669/7).

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [January 14, 2021, 7:08pm UTC](https://discourse.julialang.org/t/cuda-matrix-inverse/53341/9 "2021-01-14T19:08:26Z")

</div>

I calculate REML function for big matrices. For small dimensions all works fine (I’m using sweep operator). (I use mul! or calculation with preallocations in code). But for large matrix I want to try calculate part on GPU.

logREML(\theta,\beta) = -\frac{N-p}{2} - \frac{1}{2}\sum\_{i=1}^nlog|V\_{i}|- -\frac{1}{2}log|\sum\_{i=1}^nX\_i'V\_i^{-1}X\_i|-\frac{1}{2}\sum\_{i=1}^n(y\_i - X\_{i}\beta)'V\_i^{-1}(y\_i - X\_{i}\beta)

---

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [January 14, 2021, 9:04pm UTC](https://discourse.julialang.org/t/cuda-matrix-inverse/53341/10 "2021-01-14T21:04:01Z")

</div>

Is there some reason you couldn’t write individual terms of that like this?

```julia
function individual_neg_log_likelihood(...)
      V = Symmetric(...) # assemble your covariance matrix
      Vf = cholesky!(Vf) # assuming positive definite
      logdet(Vf) + dot(X, Vf\X) + dot(y-X*beta, Vf\(y - X*beta)) # or something
end

neg_log_lik(...) = sum(individual_neg_log_likelihood, ...)

```

I’m not being careful with signs or anything (or the fact that the logdet is outside the sum in your second expression), but hopefully you see the point. The operations that you need to do with V\_i are certainly doable with a factorization. I’ve also never used a GPU, but I would be pretty shocked if it weren’t possible to compute a Cholesky factorization and do some solves on the GPU.

Quick edit here: If `X` is a matrix and not a vector, you should change the call to `dot` in the second term to something like `X'*(Vf\X)`, or something more thoughtful. Should have mentioned that I was assuming `X` was a vector, which doesn’t make sense.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [December 6, 2021, 12:41am UTC](https://discourse.julialang.org/t/cuda-matrix-inverse/53341/11 "2021-12-06T00:41:44Z")

</div>

Bumping this to ask if there is any overload of `inv` for `CuArray` nowadays? Can we just `inv(CuArray(A))` to use the GPU?

---

<div class="post-metadata">

**Author:** ![maleadt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maleadt/32/10097_2.png) [@maleadt](https://discourse.julialang.org/u/maleadt)\
**Post date:** [December 6, 2021, 7:54am UTC](https://discourse.julialang.org/t/cuda-matrix-inverse/53341/12 "2021-12-06T07:54:49Z")

</div>

No, only for diagonal matrices.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [December 6, 2021, 10:02pm UTC](https://discourse.julialang.org/t/cuda-matrix-inverse/53341/13 "2021-12-06T22:02:01Z")

</div>

Thank you @maleadt , is that because the functionality isn’t available in the underlying CUDA library or because the Julia package didn’t wrap it yet? 🙂

---

<div class="post-metadata">

**Author:** ![maleadt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maleadt/32/10097_2.png) [@maleadt](https://discourse.julialang.org/u/maleadt)\
**Post date:** [December 7, 2021, 6:47am UTC](https://discourse.julialang.org/t/cuda-matrix-inverse/53341/14 "2021-12-07T06:47:23Z")

</div>

CUDA provides getri, getrf and getrs, so IIUC it should be possible to implement that operation.