Skip to content

Commit

Permalink
Preserve Hessenberg shape when possible (JuliaLang#40039)
Browse files Browse the repository at this point in the history 
* Preserve Hessenberg shape when possible

* Remove obsolete comments

* Add NEWS

* Support for unitful Hessenberg matrices

* Mark some tests for unitful matrices as broken
  • Loading branch information
@sostock @antoine-levitt
sostock authored and antoine-levitt committed May 9, 2021
1 parent e6da235 commit 6e52611
Show file tree
Hide file tree
Showing 4 changed files with 139 additions and 6 deletions.
1 change: 1 addition & 0 deletions NEWS.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,6 +93,7 @@ Standard library changes
* On aarch64, OpenBLAS now uses an ILP64 BLAS like all other 64-bit platforms. ([#39436])
* OpenBLAS is updated to 0.3.13. ([#39216])
* SuiteSparse is updated to 5.8.1. ([#39455])
* The shape of an `UpperHessenberg` matrix is preserved under certain arithmetic operations, e.g. when multiplying or dividing by an `UpperTriangular` matrix. ([#40039])
* `cis(A)` now supports matrix arguments ([#40194]).

#### Markdown
Expand Down
91 changes: 85 additions & 6 deletions stdlib/LinearAlgebra/src/hessenberg.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,17 +94,96 @@ Base.copy(A::Transpose{<:Any,<:UpperHessenberg}) = tril!(transpose!(similar(A.pa
rmul!(H::UpperHessenberg, x::Number) = (rmul!(H.data, x); H)
lmul!(x::Number, H::UpperHessenberg) = (lmul!(x, H.data); H)

# (future: we could also have specialized routines for UpperHessenberg * UpperTriangular)

fillstored!(H::UpperHessenberg, x) = (fillband!(H.data, x, -1, size(H,2)-1); H)

+(A::UpperHessenberg, B::UpperHessenberg) = UpperHessenberg(A.data+B.data)
-(A::UpperHessenberg, B::UpperHessenberg) = UpperHessenberg(A.data-B.data)
# (future: we could also have specialized routines for UpperHessenberg ± UpperTriangular)

# shift Hessenberg by λI
+(H::UpperHessenberg, J::UniformScaling) = UpperHessenberg(H.data + J)
-(J::UniformScaling, H::UpperHessenberg) = UpperHessenberg(J - H.data)
for T = (:UniformScaling, :Diagonal, :Bidiagonal, :Tridiagonal, :SymTridiagonal,
:UpperTriangular, :UnitUpperTriangular)
for op = (:+, :-)
@eval begin
$op(H::UpperHessenberg, x::$T) = UpperHessenberg($op(H.data, x))
$op(x::$T, H::UpperHessenberg) = UpperHessenberg($op(x, H.data))
end
end
end

for T = (:Number, :UniformScaling, :Diagonal)
@eval begin
*(H::UpperHessenberg, x::$T) = UpperHessenberg(H.data * x)
*(x::$T, H::UpperHessenberg) = UpperHessenberg(x * H.data)
/(H::UpperHessenberg, x::$T) = UpperHessenberg(H.data / x)
\(x::$T, H::UpperHessenberg) = UpperHessenberg(x \ H.data)
end
end

function *(H::UpperHessenberg, U::UpperOrUnitUpperTriangular)
T = typeof(oneunit(eltype(H))*oneunit(eltype(U)))
HH = similar(H.data, T, size(H))
copyto!(HH, H)
rmul!(HH, U)
UpperHessenberg(HH)
end
function *(U::UpperOrUnitUpperTriangular, H::UpperHessenberg)
T = typeof(oneunit(eltype(H))*oneunit(eltype(U)))
HH = similar(H.data, T, size(H))
copyto!(HH, H)
lmul!(U, HH)
UpperHessenberg(HH)
end

function /(H::UpperHessenberg, U::UpperTriangular)
T = typeof(oneunit(eltype(H))/oneunit(eltype(U)))
HH = similar(H.data, T, size(H))
copyto!(HH, H)
rdiv!(HH, U)
UpperHessenberg(HH)
end
function /(H::UpperHessenberg, U::UnitUpperTriangular)
T = typeof(oneunit(eltype(H))/oneunit(eltype(U)))
HH = similar(H.data, T, size(H))
copyto!(HH, H)
rdiv!(HH, U)
UpperHessenberg(HH)
end

function \(U::UpperTriangular, H::UpperHessenberg)
T = typeof(oneunit(eltype(U))\oneunit(eltype(H)))
HH = similar(H.data, T, size(H))
copyto!(HH, H)
ldiv!(U, HH)
UpperHessenberg(HH)
end
function \(U::UnitUpperTriangular, H::UpperHessenberg)
T = typeof(oneunit(eltype(U))\oneunit(eltype(H)))
HH = similar(H.data, T, size(H))
copyto!(HH, H)
ldiv!(U, HH)
UpperHessenberg(HH)
end

function *(H::UpperHessenberg, B::Bidiagonal)
TS = promote_op(matprod, eltype(H), eltype(B))
if B.uplo == 'U'
A_mul_B_td!(UpperHessenberg(zeros(TS, size(H)...)), H, B)
else
A_mul_B_td!(zeros(TS, size(H)...), H, B)
end
end
function *(B::Bidiagonal, H::UpperHessenberg)
TS = promote_op(matprod, eltype(B), eltype(H))
if B.uplo == 'U'
A_mul_B_td!(UpperHessenberg(zeros(TS, size(B)...)), B, H)
else
A_mul_B_td!(zeros(TS, size(B)...), B, H)
end
end

function /(H::UpperHessenberg, B::Bidiagonal)
A = Base.@invoke /(H::AbstractMatrix, B::Bidiagonal)
B.uplo == 'U' ? UpperHessenberg(A) : A
end

# Solving (H+µI)x = b: we can do this in O(m²) time and O(m) memory
# (in-place in x) by the RQ algorithm from:
Expand Down
52 changes: 52 additions & 0 deletions stdlib/LinearAlgebra/test/hessenberg.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,10 @@ module TestHessenberg

using Test, LinearAlgebra, Random

const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
isdefined(Main, :Furlongs) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "Furlongs.jl"))
using .Main.Furlongs

# for tuple tests below
≅(x,y) = all(p -> p[1] ≈ p[2], zip(x,y))

Expand Down Expand Up @@ -55,6 +59,54 @@ let n = 10
H = UpperHessenberg(Areal)
@test Array(Hc + H) == Array(Hc) + Array(H)
@test Array(Hc - H) == Array(Hc) - Array(H)
@testset "Preserve UpperHessenberg shape (issue #39388)" begin
for H = (UpperHessenberg(Areal), UpperHessenberg(Furlong.(Areal)))
if eltype(H) <: Furlong
A = Furlong.(rand(n,n))
d = Furlong.(rand(n))
dl = Furlong.(rand(n-1))
du = Furlong.(rand(n-1))
us = Furlong(1)*I
else
A = rand(n,n)
d = rand(n)
dl = rand(n-1)
du = rand(n-1)
us = 1*I
end
@testset "$op" for op = (+,-)
for x = (us, Diagonal(d), Bidiagonal(d,dl,:U), Bidiagonal(d,dl,:L),
Tridiagonal(dl,d,du), SymTridiagonal(d,dl),
UpperTriangular(A), UnitUpperTriangular(A))
@test op(H,x) == op(Array(H),x)
@test op(x,H) == op(x,Array(H))
@test op(H,x) isa UpperHessenberg
@test op(x,H) isa UpperHessenberg
end
end
A = randn(n,n)
d = randn(n)
dl = randn(n-1)
@testset "Multiplication/division" begin
for x = (5, 5I, Diagonal(d), Bidiagonal(d,dl,:U),
UpperTriangular(A), UnitUpperTriangular(A))
@test H*x == Array(H)*x broken = eltype(H) <: Furlong && x isa Bidiagonal
@test x*H == x*Array(H) broken = eltype(H) <: Furlong && x isa Bidiagonal
@test H/x == Array(H)/x broken = eltype(H) <: Furlong && x isa Union{Bidiagonal, Diagonal, UpperTriangular}
@test x\H == x\Array(H) broken = eltype(H) <: Furlong && x isa Union{Bidiagonal, Diagonal, UpperTriangular}
@test H*x isa UpperHessenberg broken = eltype(H) <: Furlong && x isa Bidiagonal
@test x*H isa UpperHessenberg broken = eltype(H) <: Furlong && x isa Bidiagonal
@test H/x isa UpperHessenberg broken = eltype(H) <: Furlong && x isa Union{Bidiagonal, Diagonal}
@test x\H isa UpperHessenberg broken = eltype(H) <: Furlong && x isa Union{Bidiagonal, Diagonal}
end
x = Bidiagonal(d, dl, :L)
@test H*x == Array(H)*x
@test x*H == x*Array(H)
@test H/x == Array(H)/x broken = eltype(H) <: Furlong
@test_broken x\H == x\Array(H) # issue 40037
end
end
end
end

@testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, Int), herm in (false, true)
Expand Down
1 change: 1 addition & 0 deletions test/testhelpers/Furlongs.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@ Base.floatmin(::Type{Furlong{p,T}}) where {p,T<:AbstractFloat} = Furlong{p}(floa
Base.floatmin(::Furlong{p,T}) where {p,T<:AbstractFloat} = floatmin(Furlong{p,T})
Base.floatmax(::Type{Furlong{p,T}}) where {p,T<:AbstractFloat} = Furlong{p}(floatmax(T))
Base.floatmax(::Furlong{p,T}) where {p,T<:AbstractFloat} = floatmax(Furlong{p,T})
Base.conj(x::Furlong{p,T}) where {p,T} = Furlong{p,T}(conj(x.val))

# convert Furlong exponent p to a canonical form. This
# is not type stable, but it doesn't matter since it is used
Expand Down

0 comments on commit 6e52611

Please sign in to comment.