# Slow lower triangular solves when compared to full Cholesky

**URL:** <https://discourse.julialang.org/t/slow-lower-triangular-solves-when-compared-to-full-cholesky/68122>\
**Category:** Numerics\
**Tags:** numerics\
**Created:** [September 13, 2021, 10:12pm UTC](https://discourse.julialang.org/t/slow-lower-triangular-solves-when-compared-to-full-cholesky/68122 "2021-09-13T22:12:59Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [September 13, 2021, 10:12pm UTC](https://discourse.julialang.org/t/slow-lower-triangular-solves-when-compared-to-full-cholesky/68122/1 "2021-09-13T22:12:59Z")

</div>

Can somebody here help me understand why in this code solving with the L in A = LL^T is an order of magnitude slower than just solving using A? Here is a MWE:

```julia
using BenchmarkTools

if !(@isdefined x) # because I lazily re-run in a REPL.
  const A = randn(5_000, 5_000)
  const Mf = cholesky!(Symmetric(A'A))
end

# With:
# BLAS.get_num_threads() == 4
# BLAS.vendor() == :mkl
# versioninfo() == (v1.6.2, official binary, etc)

print("Full solve:") # 5.050 ms, 0 alloc.
@btime ldiv!($Mf, x) (setup=x=randn(5_000))
print("Cholesky factor solve:") # 59.8 ms, 3 alloc, 190 MiB.
@btime ldiv!($Mf.L, x) (setup=x=randn(5_000))

```

At larger matrix sizes, the difference is even more pronounced. For a 22k x 22k matrix, the differences was 0.1 seconds vs 1.5 seconds. Looking at the `which(LinearAlgebra.ldiv!, ...)` outputs, it looks like with the triangular `Mf.L` Julia is calling `trtrs!`, whereas for `Mf` it is calling `potrs!`. I am sure that `potrs!` is doing something smarter than just two sequential triangular solves, but it seems a little surprising that this `ldiv!` makes big allocations and is that much slower. Intuitively, shouldn’t the solve with a triangular factor be faster, if anything?

If it matters, I have an intel CPU that does use AVX-512, and like the snippet above says, I use MKL with 4 BLAS threads. From looking around, I see other posts on the forum that say that `trtrs!` is slower by a factor of two or three from other triangular solvers, but I haven’t seen this particular issue. Apologies if this is a rehash of a known issue or discussion.

EDIT:  
Okay, there is definitely something confusing happening here. I deleted this topic because at first I think I had just made a mistake with my BenchmarkTools interpolation. When I did

```julia
@btime ldiv!($(Mf.L), x) (setup=x=randn(5000))

```

then the times looked reasonable. But the issue persists when I just do

```julia
const _x = randn(5000)
@time Mf.L\x # 0.05 seconds
@time Mf\x # 0.005 seconds

```

So it’s still an order of magnitude difference in this less careful benchmark. I honestly am not sure anymore if this is just a weird benchmarking issue or if `trtrs!` vs `potrs!` is actually causing a problem.

So many edits. Sorry everyone.

---

<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:** [September 14, 2021, 12:50am UTC](https://discourse.julialang.org/t/slow-lower-triangular-solves-when-compared-to-full-cholesky/68122/2 "2021-09-14T00:50:31Z")

</div>

Doesn’t `Mf.L` first have to copy the triangular factor into a newly allocated matrix? This has O(n^2) complexity similar to the solve itself, so it is quite a significant cost.

When you interpolate `$(Mf.L)` into `@btime`, you are precomputing this copy operation, removing it from the benchmark time.

---

<div class="post-metadata">

**Author:** ![artkuo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/artkuo/32/28577_2.png) [@artkuo](https://discourse.julialang.org/u/artkuo)\
**Post date:** [September 14, 2021, 1:33am UTC](https://discourse.julialang.org/t/slow-lower-triangular-solves-when-compared-to-full-cholesky/68122/3 "2021-09-14T01:33:16Z")

</div>

Indeed I get big differences as suggested by @stevengj:

```julia
julia> @btime ldiv!($Mf.L, x) (setup=x=randn(5000));
  169.805 ms (3 allocations: 190.73 MiB)

julia> @btime Mf.L;
  167.914 ms (3 allocations: 190.73 MiB)

```

Inspection reveals that `Mf.L` calls `getproperty`, which does a `copy` and explains almost all the time and allocations. So let’s do the copy in setup:

```julia
julia> @btime ldiv!(MfL,x) setup=(MfL=Mf.L;x=randn(5000));
  5.993 ms (0 allocations: 0 bytes)

```

I believe you’re not intended do use `.L` performantly, and it’s intended for most people to just do a solve.

---

<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:** [September 14, 2021, 1:42am UTC](https://discourse.julialang.org/t/slow-lower-triangular-solves-when-compared-to-full-cholesky/68122/4 "2021-09-14T01:42:47Z")

</div>

> [@artkuo](#):
>
> I believe you’re not intended do use `.L` performantly, and it’s intended for most people to just do a solve.

In some applications you need it. But in that case you can just retrieve the L factor once, when the factorization is constructed — it’s much cheaper than computing the Cholesky factorization in the first place.

---

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [September 14, 2021, 12:48pm UTC](https://discourse.julialang.org/t/slow-lower-triangular-solves-when-compared-to-full-cholesky/68122/5 "2021-09-14T12:48:12Z")

</div>

Thank you both! I should have thought to look in the source about the difference in `getfield` for `:L` and `:U`. I was under the impression that `getfield` with `:L` was giving me a slightly fancier or wrapped version of `getfield(Mf, :U)'`, but it actually does make the [eager copy](https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/cholesky.jl#L506).

It’s sort of an interesting performance “gotcha”, because spamming `@time Mf.U'\x` in the REPL in the same setting shows that even that is faster than the `@time Mf\x`. It makes sense, because again looking at the source and examining `Mf.factors`, it’s clear that LAPACK stores the actual output in the U matrix, and there is no `uplo` option for `cholesky!(M::Matrix, [...])`. I’m a bit curious about the eager copy. I would have expected a function like `lowerfactor(...)` or whatever to be eager and make a copy and a `getproperty` to return the lazy `getproperty(..., :U)'`. Is there a simple explanation for that choice?

In any case, I appreciate the sanity check about this (and the explanation about why the `@btime` interpolation was changing things). Thanks again!

---

<div class="post-metadata">

**Author:** ![artkuo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/artkuo/32/28577_2.png) [@artkuo](https://discourse.julialang.org/u/artkuo)\
**Post date:** [September 14, 2021, 6:15pm UTC](https://discourse.julialang.org/t/slow-lower-triangular-solves-when-compared-to-full-cholesky/68122/6 "2021-09-14T18:15:27Z")

</div>

Actually I think you are right, there may be some performance left on the table. There is no field `U` or `L` but rather `factors` which can store `U` and/or `L` depending on field `uplo`. In this case, it’s set to `U` so you are right, `Mf.U` is blazingly fast. But if you ask for the wrong triangular, it makes a curious choice:

```julia
    if d === :U
        return UpperTriangular(sym_uplo(Cuplo) == d ? Cfactors : copy(Cfactors'))

```

This says that asking for `L` (when `U` is stored) copies the adjoint and makes a view of it, whereas it seems like it could have just made a view and taken the adjoint, which should be lazy and not need a copy. (I’m not an expert, so please correct me if I’m way off base.)

Though in my mind, the problem isn’t just the eagerness of the copy, but also the fact that you wanted to look at it. It’s like Schroedinger’s cat, which we all know allocates memory once you try to observe its state :).

I wonder whether you really need `L` and what you plan to do with it. When you do a solve, it’s not like there’s any need to explicitly make L or U. It’s known that there’s a square `factors` of which only a triangular part is needed. So if `getproperty` can’t be improved, maybe your code can access `factors` directly and do the right thing with it. (EDIT: missing apostrophe)

---

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [September 14, 2021, 8:12pm UTC](https://discourse.julialang.org/t/slow-lower-triangular-solves-when-compared-to-full-cholesky/68122/7 "2021-09-14T20:12:25Z")

</div>

I hear you. In general I don’t need a specific factor. But for some applications of Hutchinson-type stochastic trace estimation, having a symmetric factor can be enormously beneficial. I shouldn’t let the discussion in this (solved) thread drift too much, but I often need to estimate traces like \text{tr}(A^{-1} B), where A = L L^T. And as it turns out, estimating the trace with

\sum\_{j=1}^M \textbf{u}\_j L^{-1} B L^{-T} \textbf{u}\_j

instead of

\sum\_{j=1}^M \textbf{u}\_j A^{-1} B \textbf{u}\_j

(where the \textbf{u}\_j vectors have iid Rademacher components) has substantial benefits in variance reduction. For the specific application in Gaussian processes that I work on, there is a brief discussion that the symmetrized version is at least as good [here](https://www.mcs.anl.gov/~anitescu/PUBLICATIONS/2014/steinChenAnitescu-AOAS14.pdf), although the theorem itself pertains to the non-symmetrized estimator.

In any case, I’ll just work with U instead and everything is fine. At this point I’m just complaining because I personally would prefer the behavior of `getfield(M, :L)` to give me the lazy `getfield(M, :U)'`. But life goes on.