Storing the LU decomposition of a matrix

Specific Domains Numerics
1

Hi!
I have a calculation in which I need the determinant and possibly the inverse of a matrix (let’s call it X). I am, therefore, trying to store the LU decomposition of X and then use it later. To this end, I have the following function:

function logpdf!(X, Y)
           Y = lu(X)
           return logdet(Y)
end

However, it seems like Y is not modified. When I run the following code, I get:

julia> let
         X = rand(ComplexF64, 100, 100)
         Y = lu(X)
         X = rand(ComplexF64, 100, 100)
         logpdf!(X, Y)
         lu(X) == Y
       end
false

What am I doing wrong? Thanks for your help.

2

You’re not modifying the input Y; instead, you’re creating a new Y local to the function, and this does not affect the original Y outside of the function.

3

But why is this the case?

4

This section of the manual might help clear things up: Scope of Variables · The Julia Language

Inside your logpdf! function, Y = lu(x) rebinds Y to a new value lu(X), but only within the function’s local scope. This doesn’t affect the global Y because it’s a different Y that only exists inside your function. Your local Y shadows the Y passed as the argument.

And if you find all this a bit confusing, you’re not alone:

5

On top of @mthelm85 's explanation, if I interpret your intentions correctly you could (1) return both X and Y, (2) make logpdf! properly in-place by defining appropriate .=, or (3) define a lu!(Y, X) for dense matrices.

  1. Define logdetY, Y = logpdf(X) to return Y so you can keep it around. Note that the Y = lu(X) inside would still allocate, which you might not want.
  2. Inside your logpdf! define an in-place Y .= lu(X)
    to re-use Y. Even though lu factorization sometimes acts like array (e.g. \ or det), it’s a struct. There currently isn’t an in-place .= defined, perhaps because lu would allocate.
  3. If you really mean to keep Y in-place, you might need something like lu!(Y, X) which isn’t defined for dense matrices. There appears to be an kinda in-place Y = lu!(X) but AFAIK it saves space with X but still allocates Y. It looks like there are in-place methods for sparse matrices though.

Either (2) or (3) should only take a few lines of code. I’m not sure if they would be considered missing methods, and thus entertained as pull requests to LinearAlgebra.

BTW I think the preferred notation would be logpdf!(Y, X) with the in-place argument first.

6

I don’t think this is a question about scope, really. It’s a question about the difference between assignment and mutation, which is a perennial confusion (not unique to Julia), and has its own manual section: Assignment expressions and assignment versus mutation

2 Likes