Inplace version for LU decomposition in Julia

Hi there!

I have to compute the lu decomposition of a matrix that is dependent - I was thinking to create a buffer first and then inplace compute the lu decomposition when needed, but after googling and searching through some threads, I still have not found the correct function to use with the LinearAlgebra.jl package.

Below is a MWE of my use case. I would like to substitute with something like! and find out how exactly to assign the buffer that would work with this.

using Distributions, LinearAlgebra, Random

rng = Xoshiro(1)
d = InverseWishart(10, [1.5 0.8 0.5 ; 0.8  2. .5 ; .5 .5 1.] )

for t in Base.OneTo(T)
    mat = rand(d)
    mat_lu =
    # do some computations with mat_lu
1 Like

Does the existing lu! fit your purpose?


Thanks! That is a bit embarassing :D.

Does this only work with Sparse Matrices or also out of the box with standard LU decomposed matrices as buffer? I.e., the following would not work (continuation from MWE):

buffer = rand(d) )
mat = rand(d)!(buffer, mat) # MethodError: no method matching lu!(::LU{Float64, Matrix{Float64}, Vector{Int64}}, ::Matrix{Float64})

I couldn’t link to it but I think the second docstring for lu! on this page is more generic. From what I understand, both the L and U parts can be stored in the matrix itself

I think the 2 argument variant only works for sparse matrices and is intended for factorizing another sparse matrix with the same sparsity structure as the previous matrix.

In your case I think you need to keep the reference to mat like so:

mat = rand(d)
buffer =
mat = rand!(mat)
buffer =!(mat)

You could maybe wrap this logic into a struct to make it more ergonomic to use.


Thanks a lot! This still seems to be allocating though, i.e.:

using BenchmarkTools

mat = rand(d)
buffer =
mat = rand!(mat)
buffer =!(mat)

@btime!($mat) #186.438 ns (1 allocation: 80 bytes)

I assume this is because mat is a standard matrix, while the! output is a LU{Float64, Matrix{Float64}, Vector{Int64}} output.

Is there a way around this? For a single argument function, I assume I would already need to provide a LU buffer, but then I would already need to destructure the time dependent covariance matrix into a LU decomposition, which is what I want to avoid in the first place.

Judging by the source code of lu!, I don’t think there is any way around allocation of ipiv without manual hacking.
Note that strictly speaking, the ! doesn’t guarantee the absence of allocation: it just warns the user that some inputs might be overwritten.

1 Like

FastLapackInterface.jl provides ways to separate memory allocation from (a subset of) LAPACK computations. It allows to create a “workspace” which can be reused for several calls to the same functions on similar matrices.

The documentation contains an example of this for LU factorization (i.e. re-using the same “workspace” for multiple calls to getrf!):

RecursiveFactorizations.jl has a non-allocating lu! where you can pass ipiv.