Add a UpperHessenberg row in LinearAlgebra.factorize?

[WalterMadelim]

My intuition is that there should also be a row about UpperHessenberg in this table (and what should be a proper factorization for it). Does my idea make sense?

[mikmoore]

Although UpperHessenberg is a factorization, in a sense, it is never returned by factorize (even when the input is UpperHessenberg) so it should not be in this table as an output structure. As for input structure, something with UpperHessenberg structure follows the already-listed rules:

julia> factorize(UpperHessenberg(randn(3, 3))) # LU in general
LU{Float64, Matrix{Float64}, Vector{Int64}}
L factor:
3×3 Matrix{Float64}:
  1.0       0.0       0.0
 -0.884411  1.0       0.0
 -0.0       0.419893  1.0
U factor:
3×3 Matrix{Float64}:
 -1.07354  -0.0216093  -1.83604
  0.0      -2.24899    -1.79199
  0.0       0.0         1.88132

julia> factorize(UpperHessenberg(triu!(randn(3, 3)))) # UpperTriangular if so
3×3 UpperTriangular{Float64, UpperHessenberg{Float64, Matrix{Float64}}}:
 0.472389   0.63327   1.01618
  ⋅        -1.58879  -0.534824
  ⋅          ⋅       -0.119193

julia> factorize(UpperHessenberg(diagm(randn(3)))) # Diagonal if so
3×3 Diagonal{Float64, Vector{Float64}}:
 0.0997544    ⋅         ⋅
  ⋅         -2.56062    ⋅
  ⋅           ⋅       -0.573808

One could collect a (partial) list of “fast-solving” types somewhere else, but factorize is not the place.

[stevengj]

@WalterMadelim

My intuition is that there should also be a row about UpperHessenberg in this table (and what should be a proper factorization for it). Does my idea make sense?

You usually don’t need to factorize an upper-Hessenberg matrix: if H is upper-Hessenberg, you can already solve Hx =b in linear time in the number of nonzeros (i.e. O(n^2) for n \times n upper-Hessenberg). This is already implemented in Julia: (H+­μI) \ x solvers for Hessenberg factorizations by stevengj · Pull Request #31853 · JuliaLang/julia · GitHub

In particular, not only can you do x = H \ b, but it also implements an in-place algorithm accessible via ldiv!(H, x .= b), including an in-place algorithm for shifted systems (H + \mu I)x = b without modifying H.

If you have an arbitrary matrix A, then there is an upper-Hessenberg factorization A = Q H Q^* that is already implemented in Julia. This is used internally as a first step for eigensolves, but is useful in its own right if you need to solve systems with A + \mu I = Q (H + \mu I) Q^* for many different shifts \mu, since you can re-use the same Hessenberg factorization H for all the shifts.

That being said, there is an open issue to support more functionality for Hessenberg matrices, and in particular fast QR, Schur, and eigen factorizations: more methods for Hessenberg factorizations · Issue #625 · JuliaLang/LinearAlgebra.jl · GitHub

[danielwe]

There is a specialized, fast ldiv!(::UpperHessenberg, _) method, but the dispatch isn’t set up such that you hit a fast path for \(::UpperHessenberg, _). The latter falls back to the generic \, which sends you through an LU factorization.

julia> @which UpperHessenberg([1.0 2.0; 3.0 4.0]) \ [5.0; 6.0]
\(A::AbstractMatrix, B::AbstractVecOrMat)
     @ LinearAlgebra ~/.julia/juliaup/julia-1.11.5+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:1118

=> LinearAlgebra.jl/src/generic.jl at release-1.11 · JuliaLang/LinearAlgebra.jl · GitHub

Seems like the implementation is focused on making Hessenberg factorization objects performant, rather than raw UpperHessenberg instances.

[stevengj]

@mikmoore

factorize(UpperHessenberg(randn(3, 3))) # LU in genera

Calling factorize on UpperHessenberg is currently a bad idea, because it ends up just calling an O(n^3) algorithm that doesn’t exploit the structure. The key point, however, is that factorization is not needed: you can use the matrix as-is in a solve.

@danielwe

There is a specialized, fast ldiv!(::UpperHessenberg, _) method, but the dispatch isn’t set up such that you hit a fast path for \(::UpperHessenberg, _). The latter falls back to the generic \, which sends you through an LU factorization.

Yikes, that’s an oversight, probably I thought \ dispatched to ldiv! at the time (the way it does for Factorization objects)? Should be an easy fix to add a specialized dispatch for this: add missing methods for division of Hessenberg matrices by stevengj · Pull Request #1322 · JuliaLang/LinearAlgebra.jl · GitHub

[WalterMadelim]

@stevengj

Calling factorize on UpperHessenberg is currently a bad idea

Does this also mean that calling factorize on a Tridiagonal is also undesired?
Note that Tridiagonal is simpler than UpperHessenberg.

julia> A
4×4 Tridiagonal{Int64, Vector{Int64}}:
 1  2  ⋅  ⋅
 1  2  3  ⋅
 ⋅  2  3  4
 ⋅  ⋅  3  4

julia> @which factorize(A)
factorize(A::Tridiagonal)
     @ LinearAlgebra K:\julia-1.11.4\share\julia\stdlib\v1.11\LinearAlgebra\src\lu.jl:656

julia> # in the place above, I see `factorize(A::Tridiagonal) = lu(A)`

julia> @which lu(A)
lu(A::AbstractMatrix{T}, args...; kwargs...) where T
     @ LinearAlgebra K:\julia-1.11.4\share\julia\stdlib\v1.11\LinearAlgebra\src\lu.jl:341

And this behavior is documented in the table (the image in my first post above).

[stevengj]

@WalterMadelim

Does this also mean that calling factorize on a Tridiagonal is also undesired?

No, because there is also a specialized factorization method for tridiagonal matrices. Although it’s true that you can solve tridiagonal systems directly, in linear time, without additional factorization.

[RoyiAvital]

@stevengj

In particular, not only can you do x = H \ b, but it also implements an in-place algorithm accessible via ldiv!(H, x .= b), including an in-place algorithm for shifted systems (H + \mu I)x = b without modifying H.

If you have an arbitrary matrix A, then there is an upper-Hessenberg factorization A = Q H Q^* that is already implemented in Julia. This is used internally as a first step for eigensolves, but is useful in its own right if you need to solve systems with A + \mu I = Q (H + \mu I) Q^* for many different shifts \mu, since you can re-use the same Hessenberg factorization H for all the shifts.

Is there a special syntax for that?

[stevengj]

@RoyiAvital

Is there a special syntax for that?

For the shifted solve? If you have a Hessenberg factorization object F, you can just do (F + μ*I) \ b, or ldiv!(F + μ*I, b), and it will compute (A + \mu I)^{-1} b without modifying the factorization. This is described in the hessenberg documentation. If you have an UpperHessenberg or SymTridiagonal matrix H, you can pass a shift=μ keyword parameter to ldiv!.