# Non-allocating matrix inversion

**URL:** <https://discourse.julialang.org/t/non-allocating-matrix-inversion/62264>\
**Category:** Performance\
**Tags:** matrices\
**Created:** [June 2, 2021, 11:48am UTC](https://discourse.julialang.org/t/non-allocating-matrix-inversion/62264 "2021-06-02T11:48:48Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![nandoconde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nandoconde/32/19497_2.png) [@nandoconde](https://discourse.julialang.org/u/nandoconde)\
**Post date:** [June 2, 2021, 11:48am UTC](https://discourse.julialang.org/t/non-allocating-matrix-inversion/62264/1 "2021-06-02T11:48:48Z")

</div>

Is there any way to perform matrix inversion?

Say I have a `Vector` of `Matrix` where I want to perform the inverse, and then save them in the same spot where they came from in the array. I had thought of this:

```julia
function bucle!(Zd)
    @simd for i in eachindex(Zd)
        @inbounds Zd[i] = inv(Zd[i])
    end
    return Zd
end

```

But it allocates:

```julia
julia> @btime bucle!($Zd)
  1.125 ms (4004 allocations: 2.34 MiB)

```

Thanks!

---

<div class="post-metadata">

**Author:** ![oheil](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oheil/32/220745_2.png) [@oheil](https://discourse.julialang.org/u/oheil)\
**Post date:** [June 2, 2021, 12:13pm UTC](https://discourse.julialang.org/t/non-allocating-matrix-inversion/62264/2 "2021-06-02T12:13:03Z")

</div>

For a single matrix this could be:

```julia
julia> function inv!(m)
          m .= inv(m)
       end
inv! (generic function with 1 method)

```

But `inv` already allocates:

```julia
julia> @btime inv($m)
  838.983 ns (4 allocations: 1.95 KiB)

julia> @btime inv!($m)
  812.727 ns (4 allocations: 1.95 KiB)

```

Your function would be than:

```julia
julia> function bucle!(Zd)
           @simd for i in eachindex(Zd)
               @inbounds inv!(Zd[i])
           end
           return Zd
       end

julia> Zd=[rand(Float64,3,3), rand(Float64,3,3), rand(Float64,3,3)]

julia> @btime bucle!($Zd)
  2.511 μs (12 allocations: 5.86 KiB)

```

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [June 2, 2021, 12:30pm UTC](https://discourse.julialang.org/t/non-allocating-matrix-inversion/62264/3 "2021-06-02T12:30:05Z")

</div>

If your matrices are small (i.e. probably less than 10x10), then you can use [https://github.com/JuliaArrays/StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) to make your original function faster and non-allocating with no changes at all to the function itself:

```julia
julia> using BenchmarkTools

julia> using StaticArrays

julia> function bucle!(Zd)
           @simd for i in eachindex(Zd)
               @inbounds Zd[i] = inv(Zd[i])
           end
           return Zd
       end
bucle! (generic function with 1 method)

julia> Zd = [rand(SMatrix{3, 3, Float64}) for _ in 1:3];

julia> @btime bucle!($Zd);
  25.224 ns (0 allocations: 0 bytes)

```

On my machine, using `SMatrix` as the element of `Zd` results in a 50X speedup over plain `Matrix` for a 3x3 input (expect this improvement to decrease as the matrices get larger and eventually go away completely for large enough matrices).

---

<div class="post-metadata">

**Author:** ![giordano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/giordano/32/2166_2.png) [@giordano](https://discourse.julialang.org/u/giordano)\
**Post date:** [June 2, 2021, 1:09pm UTC](https://discourse.julialang.org/t/non-allocating-matrix-inversion/62264/4 "2021-06-02T13:09:13Z")

</div>

> [@nandoconde](#):
>
> Say I have a `Vector` of `Matrix` where I want to perform the inverse

Am I missing something or your code is doing the element-wise reciprocal, and not matrix inversion, that you could do allocation-free with

```julia
julia> bucle!(m) = m .= inv.(m)
bucle! (generic function with 1 method)

julia> @btime bucle!(m) setup = (m = collect(reshape(1.0:100.0, 10, 10)))
  368.077 ns (0 allocations: 0 bytes)
10×10 Matrix{Float64}:
 1.0 0.0909091 0.047619 0.0322581 0.0243902 0.0196078 0.0163934 0.0140845 0.0123457 0.010989
 0.5 0.0833333 0.0454545 0.03125 0.0238095 0.0192308 0.016129 0.0138889 0.0121951 0.0108696
 0.333333 0.0769231 0.0434783 0.030303 0.0232558 0.0188679 0.015873 0.0136986 0.0120482 0.0107527
 0.25 0.0714286 0.0416667 0.0294118 0.0227273 0.0185185 0.015625 0.0135135 0.0119048 0.0106383
 0.2 0.0666667 0.04 0.0285714 0.0222222 0.0181818 0.0153846 0.0133333 0.0117647 0.0105263
 0.166667 0.0625 0.0384615 0.0277778 0.0217391 0.0178571 0.0151515 0.0131579 0.0116279 0.0104167
 0.142857 0.0588235 0.037037 0.027027 0.0212766 0.0175439 0.0149254 0.012987 0.0114943 0.0103093
 0.125 0.0555556 0.0357143 0.0263158 0.0208333 0.0172414 0.0147059 0.0128205 0.0113636 0.0102041
 0.111111 0.0526316 0.0344828 0.025641 0.0204082 0.0169492 0.0144928 0.0126582 0.011236 0.010101
 0.1 0.05 0.0333333 0.025 0.02 0.0166667 0.0142857 0.0125 0.0111111 0.01

```

?

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [June 2, 2021, 1:26pm UTC](https://discourse.julialang.org/t/non-allocating-matrix-inversion/62264/5 "2021-06-02T13:26:24Z")

</div>

Numerically speaking, you never actually want to calculate matrix inverses. The only thing a matrix inverse is good for is solving `Ax=b` for `x`, which can be done more efficiently and stably by QR decomposition and gauss elimination (and in specific cases even better, specialized techniques).

The inverse is only really useful for symbolic reasoning.

If in doubt, do `x = A\b` rather than `x = inv(A)*b`

---

<div class="post-metadata">

**Author:** ![ranocha](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ranocha/32/35588_2.png) [@ranocha](https://discourse.julialang.org/u/ranocha)\
**Post date:** [June 2, 2021, 1:28pm UTC](https://discourse.julialang.org/t/non-allocating-matrix-inversion/62264/6 "2021-06-02T13:28:16Z")

</div>

Expanding that a bit: If you want to solve small linear systems, use StaticArrays.jl. If you want to solve bigger linear systems, look for `factorize!`. Otherwise, it would be good to know more about your application to help you.

---

<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:** [June 2, 2021, 1:50pm UTC](https://discourse.julialang.org/t/non-allocating-matrix-inversion/62264/7 "2021-06-02T13:50:00Z")

</div>

> [@oheil](#):
>
> But `inv` already allocates:

There is actually an undocumented `LinearAlgebra.inv!` function that acts in-place on an LU factorization, which can be constructed in-place with `lu!`:

```julia-auto
inv!(A) = LinearAlgebra.inv!(lu!(A))

```

However, it still allocates some small 1d arrays internally:

```julia-auto
julia> A = rand(1000,1000);

julia> @btime inv($A);
  30.062 ms (5 allocations: 8.13 MiB)

julia> @btime inv!($A);
  28.390 ms (3 allocations: 508.09 KiB)

```

To completely eliminate all allocations, you need to pre-allocate the low-level workspace vectors and call lower-level LAPACK functions (based on code for [`getrf!`](https://github.com/JuliaLang/julia/blob/fd905fa2f0ab5624646404bf4f807f122a09ded1/stdlib/LinearAlgebra/src/lapack.jl#L562-L582) and [`getri!`](https://github.com/JuliaLang/julia/blob/fd905fa2f0ab5624646404bf4f807f122a09ded1/stdlib/LinearAlgebra/src/lapack.jl#L1033-L1062)).

```julia-auto
import Base: require_one_based_indexing
import LinearAlgebra: checksquare, chkstride1
import LinearAlgebra.BLAS: @blasfunc, BlasInt
import LinearAlgebra.LAPACK: chkargsok, chklapackerror, liblapack, getrf!, getri!

for (getrf, getri, elty) in ((:dgetrf_,:dgetri_,:Float64), (:sgetrf_,:sgetri_,:Float32), (:zgetrf_,:zgetri_,:ComplexF64), (:cgetrf_,:cgetri_,:ComplexF32))
    @eval function getrf!(A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt})
        require_one_based_indexing(A)
        chkstride1(A)
        chkstride1(ipiv)
        m, n = size(A)
        lda = max(1,stride(A, 2))
        length(ipiv) ≥ min(m,n) || throw(DimensionMismatch())
        info = Ref{BlasInt}()
        ccall((@blasfunc($getrf), liblapack), Cvoid,
              (Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty},
               Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}),
              m, n, A, lda, ipiv, info)
        chkargsok(info[])
        A, ipiv, info[] #Error code is stored in LU factorization type
    end
    @eval function getri!(A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}, work::AbstractVector{$elty})
        require_one_based_indexing(A, ipiv)
        chkstride1(A, ipiv)
        n = checksquare(A)
        if n != length(ipiv)
            throw(DimensionMismatch("ipiv has length $(length(ipiv)), but needs $n"))
        end
        lda = max(1,stride(A, 2))
        lwork = BlasInt(length(work))
        lwork ≥ n || throw(DimensionMismatch())
        info = Ref{BlasInt}()
        ccall((@blasfunc($getri), liblapack), Cvoid,
              (Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt},
               Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}),
              n, A, lda, ipiv, work, lwork, info)
        chklapackerror(info[])
        A
    end
end

struct InvWorkspace{T}
    ipiv::Vector{BlasInt}
    work::Vector{T}
end

function InvWorkspace(A::AbstractMatrix{T}) where {T}
    n = checksquare(A)
    ipiv = Vector{BlasInt}(undef, n)
    work = Vector{T}(undef, 64*n)
    return InvWorkspace{T}(ipiv, work)
end

function inv!(A::AbstractMatrix{T}, W::InvWorkspace{T}) where {T}
    getrf!(A, W.ipiv)
    return getri!(A, W.ipiv, W.work)
end

```

which gives

```julia-auto
julia> W = InvWorkspace(A);

julia> @btime inv!($A, $W);
  28.423 ms (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [June 2, 2021, 1:56pm UTC](https://discourse.julialang.org/t/non-allocating-matrix-inversion/62264/8 "2021-06-02T13:56:48Z")

</div>

```julia
using LinearAlgebra, BenchmarkTools, BenchmarkPlots, StatsPlots
bg = BenchmarkGroup()
bg[:inv] = @benchmarkable A*inv(B)*b setup=(A=rand(100,100);B=rand(100,100);b=rand(100))
bg[:backslash] = @benchmarkable A*(B\b) setup=(A=rand(100,100);B=rand(100,100);b=rand(100))
bg[:slash] = @benchmarkable (A/B)*b setup=(A=rand(100,100);B=rand(100,100);b=rand(100))
bg[:qr] = @benchmarkable A*(qr(B)\b) setup=(A=rand(100,100);B=rand(100,100);b=rand(100))
bg[:lu] = @benchmarkable A*(lu(B)\b) setup=(A=rand(100,100);B=rand(100,100);b=rand(100))
bg[:inv_pre] = @benchmarkable A*Q*b setup=(A=rand(100,100);B=rand(100,100);Q=inv(B);b=rand(100))
bg[:qr_pre] = @benchmarkable A*(Q\b) setup=(A=rand(100,100);B=rand(100,100);Q=qr(B);b=rand(100))
bg[:lu_pre] = @benchmarkable A*(Q\b) setup=(A=rand(100,100);B=rand(100,100);Q=lu(B);b=rand(100))
trial = run(bg)
plot(trial; yscale=:log10)

```

![image](https://global.discourse-cdn.com/julialang/original/3X/7/f/7f1c3cf13385d2244e930209a06315d324597692.png)

Even ignoring the numerical instability of the inversion process, you can see that `A*(B\b)` is faster than inversion, and that `qr(B)` is faster when used as a precomputed auxiliary.

Edit: Added in `LU` decomposition, because it’s the preferred algorithm.

---

<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:** [June 2, 2021, 1:59pm UTC](https://discourse.julialang.org/t/non-allocating-matrix-inversion/62264/9 "2021-06-02T13:59:03Z")

</div>

> [@gustaphe](#):
>
> Even ignoring the numerical instability of the inversion process, you can see that `A*(B\b)` is faster than inversion, and that `qr(B)` is faster when used as a precomputed auxiliary.

Yes, the standard caveat applies that you are normally better off not using matrix inverses at all, but rather computing the factorization object (typically LU factorization for solving linear systems).

---

<div class="post-metadata">

**Author:** ![nandoconde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nandoconde/32/19497_2.png) [@nandoconde](https://discourse.julialang.org/u/nandoconde)\
**Post date:** [June 2, 2021, 2:59pm UTC](https://discourse.julialang.org/t/non-allocating-matrix-inversion/62264/10 "2021-06-02T14:59:43Z")

</div>

Yep, I think that using `StaticArrays` for my problem will be the wisest. All of my matrices are 2x2 or 4x4.

But I have a question: why do you use `SMatrix` instead of `MMatrix`? Isn’t it mutable, since you are reusing the space?

---

<div class="post-metadata">

**Author:** ![nandoconde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nandoconde/32/19497_2.png) [@nandoconde](https://discourse.julialang.org/u/nandoconde)\
**Post date:** [June 2, 2021, 3:01pm UTC](https://discourse.julialang.org/t/non-allocating-matrix-inversion/62264/11 "2021-06-02T15:01:26Z")

</div>

Thanks for the insight. I am aware of it, but unfortunately there is no way around this. I need sometimes to plot each of the elements of the inverse matrix.

There is no multiplication afterwards, so I get no benefit from that. 😉

---

<div class="post-metadata">

**Author:** ![nandoconde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nandoconde/32/19497_2.png) [@nandoconde](https://discourse.julialang.org/u/nandoconde)\
**Post date:** [June 2, 2021, 3:02pm UTC](https://discourse.julialang.org/t/non-allocating-matrix-inversion/62264/12 "2021-06-02T15:02:25Z")

</div>

Whoa, this is deep level stuff haha. Thanks for the info!

---

<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:** [June 2, 2021, 4:14pm UTC](https://discourse.julialang.org/t/non-allocating-matrix-inversion/62264/13 "2021-06-02T16:14:52Z")

</div>

> [@nandoconde](#):
>
> But I have a question: why do you use `SMatrix` instead of `MMatrix` ? Isn’t it mutable, since you are reusing the space?

You don’t need to overwrite individual elements of the matrix — you are overwriting the whole matrix with its inverse. So just `SMatrix` is fine. (The temporary intermediate calculations are done in registers or on the stack, so no need to worry about allocations there.)

> I think that using StaticArrays for my problem will be the wisest. All of my matrices are 2x2 or 4x4.

Yes, in this case they will be an order of magnitude (at least!) faster than calling out to LAPACK.

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [April 23, 2022, 9:56pm UTC](https://discourse.julialang.org/t/non-allocating-matrix-inversion/62264/14 "2022-04-23T21:56:43Z")

</div>

Perhaps it’s worth adding to this thread that if your matrix is symmetric/hermitian positive definite, you can compute an in-place inverse with zero (pre-)allocations through a Cholesky decomposition. The undocumented `LinearAlgebra.inv!(cholesky!(A))` gets you most of the way there, but to get every last ounce of performance you need to use the LAPACK API as shown below (in both cases you get a substantial speedup compared to the naive out-of-place `inv(A)` just from the fact that you’re using a Cholesky decomposition).

Here’s how you do this through the LAPACK API. Note that this also works for `CuArray`s because `CUDA.jl` adds `CUSOLVER`-wrapping methods to `LAPACK.potr{f,s,i}!`.

```julia
using LinearAlgebra
using LinearAlgebra: BlasReal, BlasFloat

function inv!(A::Union{Symmetric{<:BlasReal},Hermitian{<:BlasFloat}})
    _, info = LAPACK.potrf!(A.uplo, A.data)
    (info == 0) || throw(PosDefException(info))
    LAPACK.potri!(A.uplo, A.data)
    return A
end

```

Note that this function relies on `A` being `<: Union{Symmetric,Hermitian}` such that only one triangle of `A.data`, indicated by `A.uplo`, is considered data both for input and output. The LAPACK functions will use the other triangle as workspace and leave garbage there upon returning. If you need `A.data` to contain a valid inverse at the end you need to copy data over from the active triangle. This is what `LinearAlgebra.inv!(cholesky!(A))` does for you and is the reason it takes a little more time.

Apart from this, `LinearAlgebra.inv!(cholesky!(A))` and `inv!(A)` are basically equivalent; they make the same LAPACK calls under the hood, and while `cholesky!` has a type instability and hence a dynamic lookup and two tiny allocations, the performance hit from this is entirely insignificant (about 0.4 μs when I isolate the difference in benchmarks, compared to \>100 μs for computing a 100x100 inverse on my laptop).

Benchmarks at 100x100:

```julia
julia> using BenchmarkTools

julia> function randposdef(N)
           A = randn(N, N)
           return Symmetric(A * A' + I)
       end
randposdef (generic function with 1 method)

julia> @btime inv(A) setup=(A = randposdef(100)) evals=1;
  177.223 μs (5 allocations: 129.17 KiB)

julia> @btime inv(cholesky!(A)) setup=(A = randposdef(100)) evals=1;
  128.572 μs (4 allocations: 78.22 KiB)

julia> @btime LinearAlgebra.inv!(cholesky!(A)) setup=(A = randposdef(100)) evals=1;
  121.586 μs (2 allocations: 48 bytes)

julia> @btime inv!(A) setup=(A = randposdef(100)) evals=1;
  114.315 μs (0 allocations: 0 bytes)

```