# LU decomposition in Julia

**URL:** <https://discourse.julialang.org/t/lu-decomposition-in-julia/77664>\
**Category:** New to Julia\
**Created:** [March 10, 2022, 2:20am UTC](https://discourse.julialang.org/t/lu-decomposition-in-julia/77664 "2022-03-10T02:20:39Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![lbc546](https://avatars.discourse-cdn.com/v4/letter/l/71e660/32.png) [@lbc546](https://discourse.julialang.org/u/lbc546)\
**Post date:** [March 10, 2022, 2:20am UTC](https://discourse.julialang.org/t/lu-decomposition-in-julia/77664/1 "2022-03-10T02:20:39Z")

</div>

Trying to rewrite the `lu_nopivot` from this answer [matrix - Perform LU decomposition without pivoting in MATLAB - Stack Overflow](https://stackoverflow.com/a/41151228) into JULIA and use only one loop.  
The julia code I wrote

```julia
using LinearAlgebra
function lu_nopivot(A)
    n = size(A, 1)
    L = Matrix{eltype(A)}(I, n, n)
    U = copy(A)
    for k = 1:n
        L[k+1:n,k] = U[k+1:n,k] / U[k,k]
        U[k+1:n,:] = U[k+1:n,:] - L[k+1:n,k]*U[k,:]  
    end
    return L, U
end

```

But calling the function  
`L, U = lu_nopivot(A)` gives an error `MethodError: no method matching *(::Vector{Float64}, ::Vector{Float64})` on `L[k+1:n,k]*U[k,:]` Tried the matlab version and it worked fine. What could be the reason that make this failed?

---

<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:** [March 10, 2022, 2:29am UTC](https://discourse.julialang.org/t/lu-decomposition-in-julia/77664/2 "2022-03-10T02:29:44Z")

</div>

You want `.*` for broadcasted multiplication or `dot` for a dot product.

That said, you could just call `lu(A, NoPivot())`

---

<div class="post-metadata">

**Author:** ![lbc546](https://avatars.discourse-cdn.com/v4/letter/l/71e660/32.png) [@lbc546](https://discourse.julialang.org/u/lbc546)\
**Post date:** [March 10, 2022, 2:52am UTC](https://discourse.julialang.org/t/lu-decomposition-in-julia/77664/4 "2022-03-10T02:52:56Z")

</div>

Think I am looking to multiply two matrices here, isn’t `*` the way to go?

---

<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:** [March 10, 2022, 3:03am UTC](https://discourse.julialang.org/t/lu-decomposition-in-julia/77664/5 "2022-03-10T03:03:19Z")

</div>

The things you are multiplying aren’t matrices, (based on the error saying you can’t use `*` on 2 `Vector{Float64}`s). You might be running an old version of the code, what happens if you restart your repl?

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [March 10, 2022, 3:07am UTC](https://discourse.julialang.org/t/lu-decomposition-in-julia/77664/6 "2022-03-10T03:07:49Z")

</div>

You want a dot product here.

```julia
julia> using LinearAlgebra

julia> x = [1.0; 2.0; 3.0]; y = [0.0; 1.0; 2.0];

julia> x*y
ERROR: MethodError: no method matching *(::Vector{Float64}, ::Vector{Float64})
Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...) at ~/packages/julia-1.7.0/share/julia/base/operators.jl:655
  *(::StridedMatrix{T}, ::StridedVector{S}) where {T<:Union{Float32, Float64, ComplexF32, ComplexF64}, S<:Real} at ~/packages/julia-1.7.0/share/julia/stdlib/v1.7/LinearAlgebra/src/matmul.jl:44
  *(::StridedVecOrMat, ::Adjoint{<:Any, <:LinearAlgebra.LQPackedQ}) at ~/packages/julia-1.7.0/share/julia/stdlib/v1.7/LinearAlgebra/src/lq.jl:266
  ...
Stacktrace:
 [1] top-level scope
   @ REPL[15]:1

julia> dot(x,y)
8.0

```

`x*y` works in Matlab because Matlab doesn’t distinguish scalars, vectors, and matrices (e.g. a Matlab scalar is really a 1 x 1 complex matrix). Julia is more particular about types in order to achieve efficiency and generality.

EDIT: I see Oscar\_Smith pointed to `dot` above.

---

<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:** [March 10, 2022, 3:16am UTC](https://discourse.julialang.org/t/lu-decomposition-in-julia/77664/7 "2022-03-10T03:16:20Z")

</div>

> [@lbc546](#):
>
> `U[k,:] `

I think the confusion here is that, even though this is selecting the k-th row of `U`, the result of a 1d slice in Julia is always treated as a “column vector” (i.e. a 1d slice is a 1d array, no matter what dimension is sliced).

So, I think what you want is:

```julia
U[k+1:n,:] = U[k+1:n,:] - L[k+1:n,k]*transpose(U[k,:])

```

to make sure that `U[k,:]` is treated as a row vector for your rank-1 update here.

See also [Unintuitive Julia result: selecting a row from a matrix](https://discourse.julialang.org/t/unintuitive-julia-result-selecting-a-row-from-a-matrix/66222) and [Problem: extracting a row from an array, returns a column](https://discourse.julialang.org/t/problem-extracting-a-row-from-an-array-returns-a-column/37331) and [Trouble Understanding Slicing](https://discourse.julialang.org/t/trouble-understanding-slicing/36821) and [RFC: Drop dimensions indexed by scalars by mbauman · Pull Request #13612 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/pull/13612) (which first implemented this) and [Taking vector transposes seriously · Issue #4774 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/issues/4774#issuecomment-59422003) (where it was proposed as “[APL](https://en.wikipedia.org/wiki/APL_(programming_language))-style indexing”).

---

<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:** [March 10, 2022, 3:32am UTC](https://discourse.julialang.org/t/lu-decomposition-in-julia/77664/8 "2022-03-10T03:32:38Z")

</div>

I would again like to point out that `lu(A, NoPivot())` will be way way way faster.  
For a 1000x1000 matrix, `lu_nopivot` takes 2.4 seconds, while `lu(A, NoPivot())` takes 185 ms (12x faster).

Also, `lu` allocates 7mb vs `lu_nopivot` which allocates 11 GB of ram (over 1000x less)

---

<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:** [March 10, 2022, 3:32am UTC](https://discourse.julialang.org/t/lu-decomposition-in-julia/77664/9 "2022-03-10T03:32:46Z")

</div>

> [@stevengj](#):
>
> So, I think what you want is:

With this fix, your code works:

```julia
julia> A = rand(3,3)
3×3 Matrix{Float64}:
 0.244342 0.6094 0.281167
 0.0340127 0.130047 0.941087
 0.7212 0.69946 0.656365

julia> L, U = lu_nopivot(A);
([1.0 0.0 0.0; 0.1392011090670886 1.0 0.0; 2.9515958101582 -24.309759573913595 1.0], [0.24434249240692318 0.609400324922777 0.2811674033872621; 0.0 0.04521819259713014 0.9019484309838302; 0.0 0.0 21.75262222865455])

julia> L
3×3 Matrix{Float64}:
 1.0 0.0 0.0
 0.139201 1.0 0.0
 2.9516 -24.3098 1.0

julia> U
3×3 Matrix{Float64}:
 0.244342 0.6094 0.281167
 0.0 0.0452182 0.901948
 0.0 0.0 21.7526

julia> F = lu(A, NoPivot())
LU{Float64, Matrix{Float64}}
L factor:
3×3 Matrix{Float64}:
 1.0 0.0 0.0
 0.139201 1.0 0.0
 2.9516 -24.3098 1.0
U factor:
3×3 Matrix{Float64}:
 0.244342 0.6094 0.281167
 0.0 0.0452182 0.901948
 0.0 0.0 21.7526

```

> [@Oscar\_Smith](#):
>
> I would again like to point out that `lu(A, NoPivot())` will be way way way faster.

Definitely, but anyone who is doing LU without pivoting (which is numerically unstable) doesn’t care about performance anyway — it’s only used for learning exercises.