(H+­μI) \ x solvers for Hessenberg factorizations (#31853)

NEWS.md

@@ -36,6 +36,8 @@ Standard library changes
 * The BLAS submodule no longer exports `dot`, which conflicts with that in LinearAlgebra ([#31838]).
 * `diagm` and `spdiagm` now accept optional `m,n` initial arguments to specify a size ([#31654]).
 
+* `Hessenberg` factorizations `H` now support efficient shifted solves `(H+µI) \ b` and determinants, and use a specialized tridiagonal factorization for Hermitian matrices. There is also a new `UpperHessenberg` matrix type ([#31853]).
+
 #### SparseArrays
 
 

stdlib/LinearAlgebra/docs/src/index.md

@@ -168,8 +168,9 @@ as well as whether hooks to various optimized methods for them in LAPACK are ava
 | [`Hermitian`](@ref)           | [Hermitian matrix](https://en.wikipedia.org/wiki/Hermitian_matrix)                            |
 | [`UpperTriangular`](@ref)     | Upper [triangular matrix](https://en.wikipedia.org/wiki/Triangular_matrix)                    |
 | [`UnitUpperTriangular`](@ref) | Upper [triangular matrix](https://en.wikipedia.org/wiki/Triangular_matrix) with unit diagonal |
-| [`LowerTriangular`](@ref)     | Lower [triangular matrix](https://en.wikipedia.org/wiki/Triangular_matrix)                    |
+| [`LowerTriangular`](@ref)     | Lower [triangular matrix](https://en.wikipedia.org/wiki/Triangular_matrix)                    |     |
 | [`UnitLowerTriangular`](@ref) | Lower [triangular matrix](https://en.wikipedia.org/wiki/Triangular_matrix) with unit diagonal |
+| [`UpperHessenberg`](@ref)     | Upper [Hessenberg matrix](https://en.wikipedia.org/wiki/Hessenberg_matrix)
 | [`Tridiagonal`](@ref)         | [Tridiagonal matrix](https://en.wikipedia.org/wiki/Tridiagonal_matrix)                        |
 | [`SymTridiagonal`](@ref)      | Symmetric tridiagonal matrix                                                                  |
 | [`Bidiagonal`](@ref)          | Upper/lower [bidiagonal matrix](https://en.wikipedia.org/wiki/Bidiagonal_matrix)              |
@@ -186,6 +187,7 @@ as well as whether hooks to various optimized methods for them in LAPACK are ava
 | [`UnitUpperTriangular`](@ref) |     |     | MV  | MV  | [`inv`](@ref), [`det`](@ref)                                |
 | [`LowerTriangular`](@ref)     |     |     | MV  | MV  | [`inv`](@ref), [`det`](@ref)                                |
 | [`UnitLowerTriangular`](@ref) |     |     | MV  | MV  | [`inv`](@ref), [`det`](@ref)                                |
+| [`UpperHessenberg`](@ref)     |     |     |     | MM  | [`inv`](@ref), [`det`](@ref)                                |
 | [`SymTridiagonal`](@ref)      | M   | M   | MS  | MV  | [`eigmax`](@ref), [`eigmin`](@ref)                          |
 | [`Tridiagonal`](@ref)         | M   | M   | MS  | MV  |                                                             |
 | [`Bidiagonal`](@ref)          | M   | M   | MS  | MV  |                                                             |
@@ -269,6 +271,12 @@ Stacktrace:
 [...]
 ```
 
+If you need to solve many systems of the form `(A+μI)x = b` for the same `A` and different `μ`, it might be beneficial
+to first compute the Hessenberg factorization `F` of `A` via the [`hessenberg`](@ref) function.
+Given `F`, Julia employs an efficient algorithm for `(F+μ*I) \ b` (equivalent to `(A+μ*I)x \ b`) and related
+operations like determinants.
+
+
 ## [Matrix factorizations](@id man-linalg-factorizations)
 
 [Matrix factorizations (a.k.a. matrix decompositions)](https://en.wikipedia.org/wiki/Matrix_decomposition)
@@ -319,6 +327,7 @@ LinearAlgebra.LowerTriangular
 LinearAlgebra.UpperTriangular
 LinearAlgebra.UnitLowerTriangular
 LinearAlgebra.UnitUpperTriangular
+LinearAlgebra.UpperHessenberg
 LinearAlgebra.UniformScaling
 LinearAlgebra.lu
 LinearAlgebra.lu!

stdlib/LinearAlgebra/src/LinearAlgebra.jl

@@ -52,6 +52,7 @@ export
     UpperTriangular,
     UnitLowerTriangular,
     UnitUpperTriangular,
+    UpperHessenberg,
     Diagonal,
     UniformScaling,
 
@@ -356,7 +357,6 @@ include("triangular.jl")
 
 include("factorization.jl")
 include("qr.jl")
-include("hessenberg.jl")
 include("lq.jl")
 include("eigen.jl")
 include("svd.jl")
@@ -367,6 +367,7 @@ include("bunchkaufman.jl")
 include("diagonal.jl")
 include("bidiag.jl")
 include("uniformscaling.jl")
+include("hessenberg.jl")
 include("givens.jl")
 include("special.jl")
 include("bitarray.jl")

stdlib/LinearAlgebra/src/adjtrans.jl

@@ -273,7 +273,7 @@ pinv(v::TransposeAbsVec, tol::Real = 0) = pinv(conj(v.parent)).parent
 \(u::AdjOrTransAbsVec, v::AdjOrTransAbsVec) = pinv(u) * v
 
 
-## right-division \
+## right-division /
 /(u::AdjointAbsVec, A::AbstractMatrix) = adjoint(adjoint(A) \ u.parent)
 /(u::TransposeAbsVec, A::AbstractMatrix) = transpose(transpose(A) \ u.parent)
 /(u::AdjointAbsVec, A::Transpose{<:Any,<:AbstractMatrix}) = adjoint(conj(A.parent) \ u.parent) # technically should be adjoint(copy(adjoint(copy(A))) \ u.parent)

stdlib/LinearAlgebra/src/factorization.jl

@@ -71,13 +71,18 @@ function Base.show(io::IO, ::MIME"text/plain", x::Transpose{<:Any,<:Factorizatio
 end
 
 # With a real lhs and complex rhs with the same precision, we can reinterpret
-# the complex rhs as a real rhs with twice the number of columns
+# the complex rhs as a real rhs with twice the number of columns or rows
 function (\)(F::Factorization{T}, B::VecOrMat{Complex{T}}) where T<:BlasReal
     require_one_based_indexing(B)
     c2r = reshape(copy(transpose(reinterpret(T, reshape(B, (1, length(B)))))), size(B, 1), 2*size(B, 2))
     x = ldiv!(F, c2r)
     return reshape(copy(reinterpret(Complex{T}, copy(transpose(reshape(x, div(length(x), 2), 2))))), _ret_size(F, B))
 end
+function (/)(B::VecOrMat{Complex{T}}, F::Factorization{T}) where T<:BlasReal
+    require_one_based_indexing(B)
+    x = rdiv!(copy(reinterpret(T, B)), F)
+    return copy(reinterpret(Complex{T}, x))
+end
 
 function \(F::Factorization, B::AbstractVecOrMat)
     require_one_based_indexing(B)
@@ -95,6 +100,24 @@ function \(adjF::Adjoint{<:Any,<:Factorization}, B::AbstractVecOrMat)
     ldiv!(adjoint(F), BB)
 end
 
+function /(B::AbstractMatrix, F::Factorization)
+    require_one_based_indexing(B)
+    TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F)))
+    BB = similar(B, TFB, size(B))
+    copyto!(BB, B)
+    rdiv!(BB, F)
+end
+function /(B::AbstractMatrix, adjF::Adjoint{<:Any,<:Factorization})
+    require_one_based_indexing(B)
+    F = adjF.parent
+    TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F)))
+    BB = similar(B, TFB, size(B))
+    copyto!(BB, B)
+    rdiv!(BB, adjoint(F))
+end
+/(adjB::AdjointAbsVec, adjF::Adjoint{<:Any,<:Factorization}) = adjoint(adjF.parent \ adjB.parent)
+/(B::TransposeAbsVec, adjF::Adjoint{<:Any,<:Factorization}) = adjoint(adjF.parent \ adjoint(B))
+
 # support the same 3-arg idiom as in our other in-place A_*_B functions:
 function ldiv!(Y::AbstractVecOrMat, A::Factorization, B::AbstractVecOrMat)
     require_one_based_indexing(Y, B)
@@ -120,3 +143,10 @@ end
 # fallback methods for transposed solves
 \(F::Transpose{<:Any,<:Factorization{<:Real}}, B::AbstractVecOrMat) = adjoint(F.parent) \ B
 \(F::Transpose{<:Any,<:Factorization}, B::AbstractVecOrMat) = conj.(adjoint(F.parent) \ conj.(B))
+
+/(B::AbstractMatrix, F::Transpose{<:Any,<:Factorization{<:Real}}) = B / adjoint(F.parent)
+/(B::AbstractMatrix, F::Transpose{<:Any,<:Factorization}) = conj.(conj.(B) / adjoint(F.parent))
+/(B::AdjointAbsVec, F::Transpose{<:Any,<:Factorization{<:Real}}) = B / adjoint(F.parent)
+/(B::TransposeAbsVec, F::Transpose{<:Any,<:Factorization{<:Real}}) = B / adjoint(F.parent)
+/(B::AdjointAbsVec, F::Transpose{<:Any,<:Factorization}) = conj.(conj.(B) / adjoint(F.parent))
+/(B::TransposeAbsVec, F::Transpose{<:Any,<:Factorization}) = conj.(conj.(B) / adjoint(F.parent))

stdlib/LinearAlgebra/src/givens.jl

@@ -248,6 +248,8 @@ function givensAlgorithm(f::Complex{T}, g::Complex{T}) where T<:AbstractFloat
     return cs, sn, r
 end
 
+givensAlgorithm(f, g) = givensAlgorithm(promote(float(f), float(g))...)
+
 """
 
     givens(f::T, g::T, i1::Integer, i2::Integer) where {T} -> (G::Givens, r::T)