# Storing the LU decomposition of a matrix

**URL:** <https://discourse.julialang.org/t/storing-the-lu-decomposition-of-a-matrix/122138>\
**Category:** Numerics\
**Tags:** question\
**Created:** [November 1, 2024, 2:42pm UTC](https://discourse.julialang.org/t/storing-the-lu-decomposition-of-a-matrix/122138 "2024-11-01T14:42:53Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Gattu\_Mytraya](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gattu_mytraya/32/45429_2.png) [@Gattu\_Mytraya](https://discourse.julialang.org/u/Gattu_Mytraya)\
**Post date:** [November 1, 2024, 2:42pm UTC](https://discourse.julialang.org/t/storing-the-lu-decomposition-of-a-matrix/122138/1 "2024-11-01T14:42:53Z")

</div>

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:

```julia
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
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.

---

<div class="post-metadata">

**Author:** ![mthelm85](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mthelm85/32/224164_2.png) [@mthelm85](https://discourse.julialang.org/u/mthelm85)\
**Post date:** [November 1, 2024, 3:33pm UTC](https://discourse.julialang.org/t/storing-the-lu-decomposition-of-a-matrix/122138/2 "2024-11-01T15:33:14Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Gattu\_Mytraya](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gattu_mytraya/32/45429_2.png) [@Gattu\_Mytraya](https://discourse.julialang.org/u/Gattu_Mytraya)\
**Post date:** [November 1, 2024, 3:42pm UTC](https://discourse.julialang.org/t/storing-the-lu-decomposition-of-a-matrix/122138/3 "2024-11-01T15:42:42Z")

</div>

But why is this the case?

---

<div class="post-metadata">

**Author:** ![mthelm85](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mthelm85/32/224164_2.png) [@mthelm85](https://discourse.julialang.org/u/mthelm85)\
**Post date:** [November 1, 2024, 4:11pm UTC](https://discourse.julialang.org/t/storing-the-lu-decomposition-of-a-matrix/122138/4 "2024-11-01T16:11:59Z")

</div>

This section of the manual might help clear things up: [Scope of Variables · The Julia Language](https://docs.julialang.org/en/v1/manual/variables-and-scoping/#scope-of-variables)

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](https://en.wikipedia.org/wiki/Variable_shadowing) the `Y` passed as the argument.

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

> [@The manual's section on variable scope sucks](https://discourse.julialang.org/t/the-manuals-section-on-variable-scope-sucks/119522):
>
> I know it’s not for the lack of trying, but the manual’s section on scope is unclear if not outright wrong. Issue 1 The manual says: If you assign to an existing local, it always [emphasis not mine] updates that existing local: you can only shadow a local by explicitly declaring a new local in a nested scope with the local keyword. a few lines below, it says: When x = occurs in a local scope, Julia applies the following rules to decide what the expression means based on where the assign…

---

<div class="post-metadata">

**Author:** ![apo383](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/apo383/32/11272_2.png) [@apo383](https://discourse.julialang.org/u/apo383)\
**Post date:** [November 1, 2024, 7:27pm UTC](https://discourse.julialang.org/t/storing-the-lu-decomposition-of-a-matrix/122138/5 "2024-11-01T19:27:31Z")

</div>

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.

---

<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:** [November 1, 2024, 7:57pm UTC](https://discourse.julialang.org/t/storing-the-lu-decomposition-of-a-matrix/122138/6 "2024-11-01T19:57:46Z")

</div>

> [@mthelm85](#):
>
> This section of the manual might help clear things up: [Scope of Variables · The Julia Language](https://docs.julialang.org/en/v1/manual/variables-and-scoping/#scope-of-variables)

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](https://docs.julialang.org/en/v1/manual/variables/#man-assignment-expressions)