mfalt
commented
5 years ago I could try to implement these features. AFAIK H\b is usually solved as in
https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/blob/master/src/eigenGeneral.jl#L28
using some simple iterations like
"""Assuming Hessenberg H"""
function ldiv!(H, rhs)
n = size(H,1)
for i = 1:n-1
Gi, r = givens(H, i, i+1, i)
lmul!(Gi, H) # Note: slow indexing
lmul!(Gi, rhs)
end
# H is now upper triangular
ldiv!(UpperTriangular(H), rhs)
return rhs
endHowever, the operation lmul!(Gi, H) modifies two rows of H, which can be quite slow compared to working on the columns, there are ways to get around this but they are quite ugly, for example:
"""Assuming Hessenberg H"""
function ldiv!(H::AbstractMatrix{T}, rhs) where {T}
n = size(H,1)
Ht = H'
rotations = LinearAlgebra.Givens{T}[] # Can be replaced with LinearAlgebra.Rotation
for i = n:-1:2
Gi, r = givens(Ht, i, i-1, i)
push!(rotations, Gi')
rmul!(H, Gi')
end
# H is now upper triangular
ldiv!(UpperTriangular(H), rhs)
# Should be possible to do lmul!(rotation::Rotation, rhs) here
for i in n-1:-1:1
lmul!(rotations[i], rhs)
end
return rhs
endwhere it can be cleaned up a bit if https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/givens.jl#L368 is relaxed to accept vectors (Rotation*Vector currently throws error). This second version is a bit faster, for the cost of readability, it uses the fact that
Hx=b <==>
H(G1*G2*...Gn)*(Gn'*Gn-1'*...*G1')x=b <==>
(Gn'*Gn-1'*...*G1')x = (H*(G1*G2*...Gn))\b
x = (Gn'*Gn-1'*...*G1')*(H*(G1*G2*...Gn))\bWhere the G1,G2,...Gn that are needed to make H*(G1*G2*...Gn) upper triangular are calculated as the transpose of givens on the transpose of H, starting with the last column, instead of first row.
Is it worthwhile using this approach or is there a good standard way of making the original version efficient from a memory-layout perspective?
Edit: the original discussion is found here https://discourse.julialang.org/t/how-to-optimize-simple-matrix-equation-solving/23035/16
stevengj
olof3
andreasnoack
I just noticed that the LinearAlgebra package only supports a very small number of specialized functions for
Hessenbergfactorization objects (#7965), which limits their practical utility. Some things that would be nice to have:H \ b: this can be done in O(n²) operations. Currently it gives aMethodError. Actually, as discussed below, we can even support fast(H-µI) \ b, where different shifts µ can be used without making copies. #31853H*scalarandH+µ*Ishould return anotherHessenbergobject; currently these give aMethodError. Ideally, as discussed below, the shift µ should be stored in theHessenbergobject itself, so that multiple shifts can be applied without allocating a whole new matrix. #31853H+µ*Iin O(n²) operations, e.g. using theorem 2.1 of this paper or maybe using Givens rotations to make H upper-triangular (though avoiding copying the whole matrix would be good). #31853hessenberg(::Hermitian)or real-symmetric — in this case the Hessenberg matrix is tridiagonal with associated speedups. #31853This came up on discourse, where someone wanted to solve
(A - λ*I) \ bproblems repeatedly for many differentλ— a nice way to do this is via a Hessenberg factorization ofA = QHQ'sinceA - λ*I = Q(H-λ*I)Q'and subsequent solves are O(n²). You could also use Schur for this sort of thing, but computing the Schur factorization is a lot more expensive (≈7×).