I have been looking at a similar problem, so here is an example of how this can be solved quite a bit quicker
H = randn(100,100)
b = randn(100)
λs = 1.0:1:1000
@btime testslow(H, λs, b) #677.874 ms (7001 allocations: 154.48 MiB)
@btime testfast(H, λs, b) #48.968 ms (3018 allocations: 670.38 KiB)
with the code
using LinearAlgebra
# Solve H\rhs and store result in rhs, where H is on upper hessenberg form
function backsolve_hess!(H::AbstractMatrix{T}, rhs::AbstractArray{T}) where {T}
# Make H UpperTriangular and apply givens to rhs
n = size(H,1)
for i = 1:n-1
Gi, r = givens(H, i, i+1, i)
# Note: slow indexing
lmul!(Gi, H)
lmul!(Gi, rhs)
end
# H is now upper triangular
ldiv!(UpperTriangular(H), rhs)
return rhs
end
# Solve tmp = (Q*H*Q'+λI)\rhs, hiven hessenberg H, unitary Q, and b = Q'*rhs
function hessenberg_solve!(tmp, H, Q, λ, b, tmpH)
tmpH .= H
# tmp = Q'*rhs
tmp .= b
for i = 1:size(tmpH,1)
tmpH[i,i] += λ
end
# tmp = (H+λI)⁻¹Q'*rhs
backsolve_hess!(tmpH, tmp)
# tmp = Q*(H+λI)⁻¹Q'*rhs
lmul!(Q, tmp)
return tmp
end
function testslow(A, λs, b)
s = zero(eltype(A))
for λ in λs
tmp = (A + λ*I)\b
s += norm(tmp)
end
return s
end
function testfast(A, λs, b)
# A = Q*H*Q'
HF = hessenberg(A)
# Type insability here
H = HF.H::typeof(A)
Q = HF.Q::LinearAlgebra.HessenbergQ{eltype(A), typeof(A)}
tmpH = similar(A)
tmp = similar(b)
# b2 = Q⁻¹b=Q'*b
b2 = Q'*b
s = zero(eltype(HF))
for λ in λs
hessenberg_solve!(tmp, H, Q, λ, b2, tmpH)
s += norm(tmp)
end
return s
end
It is noteworthy that the line lmul!(Gi, H) takes most of the time. It actually accesses H roughly like H[k,:], so it is probably possible to make the code even faster by working with the transpose some way.
.