Triangular solve and inversions have poor performance

[jpdoane]

Solving and inverting triangular matrices seems to be implemented inefficiently:

n=8
X = randn(n, n)
L = LowerTriangular(randn(n, n))
U = L'

# These are mathematically equivalent but differ in performance
@assert (X'/U)' ≈ L\X  
@btime $L \ $X          # 2.492 μs (2 allocations: 592 bytes)
@btime ($X' / $U)'      # 409.151 ns (2 allocations: 592 bytes)

# This also occurs when inverting triangular matrics
I_F64 = Float64.(collect(I(n)))

@btime inv($L)          # 2.551 μs (2 allocations: 592 bytes)
@btime inv($U)'         # 2.492 μs (2 allocations: 592 bytes)
@btime $L \ I           # 2.427 μs (2 allocations: 592 bytes)
@btime (I / $U)'        # 2.496 μs (2 allocations: 592 bytes)
@btime $L \ $I_F64      # 2.418 μs (2 allocations: 592 bytes)
@btime ($I_F64 / $U)'   # 351.958 ns (2 allocations: 592 bytes)


julia> versioninfo()
Julia Version 1.11.5
Commit 760b2e5b739 (2025-04-14 06:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 112 × Intel(R) Xeon(R) Gold 6348 CPU @ 2.60GHz
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, icelake-server)
Threads: 32 default, 0 interactive, 16 GC (on 112 virtual cores)

This seems to be due to different BLAS functions being called under the hood: LinearAlgebra.generic_mattridiv! calls trtrs, while LinearAlgebra.generic_trimatdiv! call trsm, which appears to be much slower for some reason. Oddly, this Stack Overflow comment suggests that trtrs typically wraps trsm in most implementations and thus shouldn't be any faster, so thats a bit odd...

Likelrelated to issue 636

Suggested fix is replacing calls to trsm with trtrs

[dkarrasch]

I think there is a typo: @btime $L / $X # 4.894 μs (8 allocations: 1.86 KiB) should be @btime $L \ $X # 4.894 μs (8 allocations: 1.86 KiB), otherwise you're solving with the dense X. Differences in right- and left-solves should be due to memory access patterns.

[jpdoane]

Argh, thank you, that is indeed have a typo in my MWE code. I edited to original post to correct

Issue still stands, but represents ~5x performance difference rather than 10x difference

@assert (X'/U)' ≈ L\X  # These are mathematically equivalent but differ in performance
@btime $L \ $X          # 2.492 μs (2 allocations: 592 bytes)
@btime ($X' / $U)'      # 409.151 ns (2 allocations: 592 bytes)

[dkarrasch]

Have you tried calling the other BLAS function manually and benchmark? We have tried to change something, but then it got reverted: #1194... The change was simple, but I didn't have the stamina to go down the rabbit hole.

[dkarrasch]

The LAPACK one is literally just a "wrapper" of the BLAS one, adding zero pivot handling but nothing else. But who knows if there is some subtle threading issue behind?

[jpdoane]

Yes, BLAS.trsm configured as an ldiv seems to be faster than LAPACK.trtrs

function ldiv_trtrs!(C::Matrix{T}, L::LowerTriangular{T}, X::AbstractMatrix{T}) where T
    uplo = LinearAlgebra.uplo_char(L)
    tfunc  = LinearAlgebra.wrapperop(parent(L)) === identity ? 'N' : LinearAlgebra.wrapperop(parent(L)) === transpose ? 'T' : 'C'
    isunitc = LinearAlgebra.isunit_char(L)
    LAPACK.trtrs!(uplo,  tfunc, isunitc, LinearAlgebra._unwrap_at(parent(L)), copy!(C, X))
end

function ldiv_trsm!(C::Matrix{T}, L::LowerTriangular{T}, X::AbstractMatrix{T}) where T
    uplo = LinearAlgebra.uplo_char(L)
    tfunc  = LinearAlgebra.wrapperop(parent(L)) === identity ? 'N' : LinearAlgebra.wrapperop(parent(L)) === transpose ? 'T' : 'C'
    isunitc = LinearAlgebra.isunit_char(L)
    BLAS.trsm!('L', uplo, tfunc, isunitc, one(T), LinearAlgebra._unwrap_at(parent(L)), copy!(C, X))
end

n=8
X = randn(n, n)
L = LowerTriangular(randn(n, n))
C=copy(X)

@assert L\X ≈ ldiv_trtrs!(C,L,X)
@assert L\X ≈ ldiv_trsm!(C,L,X)

@btime $L \ $X                  # 2.375 μs (2 allocations: 592 bytes)
@btime ldiv_trtrs!($C,$L,$X)    # 2.234 μs (0 allocations: 0 bytes)
@btime ldiv_trsm!($C,$L,$X)     # 395.090 ns (0 allocations: 0 bytes)

[jpdoane]

Aha - it does appear to be threading related!

# default (8 threads)
BLAS.set_num_threads(8)
@btime $L \ $X                  # 2.375 μs (2 allocations: 592 bytes)
@btime ldiv_trtrs!($C,$L,$X)    # 2.234 μs (0 allocations: 0 bytes)
@btime ldiv_trsm!($C,$L,$X)     # 395.090 ns (0 allocations: 0 bytes)

# single threaded
BLAS.set_num_threads(1)
@btime $L \ $X                  # 441.056 ns (2 allocations: 592 bytes)
@btime ldiv_trtrs!($C,$L,$X)    # 403.565 ns (0 allocations: 0 bytes)
@btime ldiv_trsm!($C,$L,$X)     # 396.313 ns (0 allocations: 0 bytes)

[dkarrasch]

That seems to confirm what @amontoison claimed/observed, which was the reason for #1194. So, perhaps we try to pick it up again? We would need to figure out what went wrong last time and then we could re-land #1194. There were more reasons to do it, actually.

[dkarrasch]

BTW, what's the difference on larger matrices? Does the difference become negligible?

[jpdoane]

Yes, for bigger matrices the difference between trtrs and trsm becomes negligible for both single and multithreaded, and multithreaded performance is faster. Crossover point seems to be about n=25 or so

[stevengj]

See also https://discourse.julialang.org/t/ldiv-for-cholesky-is-slower-than-two-substitutions/130092/13?u=stevengj — apparently for OpenBLAS it is faster to do trtrs for a single rhs, and it automatically calls trsm otherwise?