From f23633e1465795501291d9c54b9a69363b17d01f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 15 Apr 2026 17:34:18 +0200 Subject: [PATCH 01/10] Import eigenbox from IntervalLinearAlgebra --- ext/IntervalArithmeticLinearAlgebraExt.jl | 107 +++++++++++++++++++++- src/IntervalArithmetic.jl | 11 +++ 2 files changed, 117 insertions(+), 1 deletion(-) diff --git a/ext/IntervalArithmeticLinearAlgebraExt.jl b/ext/IntervalArithmeticLinearAlgebraExt.jl index 7dae8fb1..07ed8580 100644 --- a/ext/IntervalArithmeticLinearAlgebraExt.jl +++ b/ext/IntervalArithmeticLinearAlgebraExt.jl @@ -151,7 +151,9 @@ function LinearAlgebra.eigen!(A::AbstractMatrix{<:RealOrComplexI}; permute::Bool _fold_conjugate!(eltype(A), true_λ, true_v) foreach(j -> true_v[:,j] .*= interval(ref_scale[j]), 1:n) else - true_λ = fill(nai(eltype(λ_bar)), n) + # contraction mapping failed; fall back to eigenbox for eigenvalue bounds + box = IntervalArithmetic.eigenbox(A) + true_λ = fill(box, n) true_v = fill(nai(eltype(v_bar)), n, n) end _ensure_ng_flag!(true_v, all(isguaranteed, A)) @@ -177,6 +179,109 @@ function _fold_conjugate!(::Type{<:Interval}, λ, v) return λ, v end +# eigenvalue enclosure (eigenbox) + +struct Orthants + n::Int +end + +Base.eltype(::Type{Orthants}) = Vector{Int} +Base.length(O::Orthants) = 2^(O.n) + +function Base.iterate(O::Orthants, state=1) + state > 2 ^ O.n && return nothing + vec = -2*digits(state-1, base=2, pad=O.n) .+ 1 + return (vec, state+1) +end + +function Base.getindex(O::Orthants, i::Int) + 1 <= i <= length(O) || throw(BoundsError(O, i)) + return -2*digits(i-1, base=2, pad=O.n) .+ 1 +end + +Base.firstindex(O::Orthants) = 1 +Base.lastindex(O::Orthants) = length(O) + +# rigorous eigmax/eigmin via interval eigvals + +function _interval_eigmax(A::LinearAlgebra.Symmetric{<:Real, <:AbstractMatrix{<:Real}}) + λs = LinearAlgebra.eigvals(LinearAlgebra.Symmetric(interval.(A))) + return sup(maximum(real.(λs))) +end + +function _interval_eigmin(A::LinearAlgebra.Symmetric{<:Real, <:AbstractMatrix{<:Real}}) + λs = LinearAlgebra.eigvals(LinearAlgebra.Symmetric(interval.(A))) + return inf(minimum(real.(λs))) +end + +function IntervalArithmetic.eigenbox(A::LinearAlgebra.Symmetric{Interval{T}, Matrix{Interval{T}}}, ::IntervalArithmetic.Rohn) where {T} + AΔ = LinearAlgebra.Symmetric(IntervalArithmetic.radius.(A)) + Ac = LinearAlgebra.Symmetric(mid.(A)) + + ρ = _interval_eigmax(AΔ) + λmax = _interval_eigmax(Ac) + λmin = _interval_eigmin(Ac) + return interval(λmin - ρ, λmax + ρ) +end + +function IntervalArithmetic.eigenbox(A::LinearAlgebra.Symmetric{Interval{T}, Matrix{Interval{T}}}, ::IntervalArithmetic.Hertz) where {T} + n = LinearAlgebra.checksquare(A) + Amax = Matrix{T}(undef, n, n) + Amin = Matrix{T}(undef, n, n) + + λmin = T(Inf) + λmax = T(-Inf) + @inbounds for z in Orthants(n) + first(z) < 0 && continue + for j in 1:n + for i in 1:j + if z[i] == z[j] + Amax[i, j] = sup(A[i, j]) + Amin[i, j] = inf(A[i, j]) + else + Amax[i, j] = inf(A[i, j]) + Amin[i, j] = sup(A[i, j]) + end + end + end + + candmax = _interval_eigmax(LinearAlgebra.Symmetric(Amax)) + candmin = _interval_eigmin(LinearAlgebra.Symmetric(Amin)) + λmin = min(λmin, candmin) + λmax = max(λmax, candmax) + end + return interval(λmin, λmax) +end + +function IntervalArithmetic.eigenbox(A::AbstractMatrix{Interval{T}}, + method::IntervalArithmetic.AbstractIntervalEigenSolver) where {T} + λ = IntervalArithmetic.eigenbox(LinearAlgebra.Symmetric(interval.(T, 0.5*(A + A'))), method) + + n = size(A, 1) + S = 0.5*(A - A') + μ = IntervalArithmetic.eigenbox(LinearAlgebra.Symmetric(interval.(T, [zeros(Interval{T}, n, n) S'; S zeros(Interval{T}, n, n)])), method) + + return λ + μ*im +end + +function IntervalArithmetic.eigenbox(M::AbstractMatrix{Complex{Interval{T}}}, + method::IntervalArithmetic.AbstractIntervalEigenSolver) where {T} + A = real.(M) + B = imag.(M) + λ = IntervalArithmetic.eigenbox(LinearAlgebra.Symmetric(interval.(T, 0.5*[A+A' B'-B; B-B' A+A'])), method) + μ = IntervalArithmetic.eigenbox(LinearAlgebra.Symmetric(interval.(T, 0.5*[B+B' A-A'; A'-A B+B'])), method) + return λ + μ*im +end + +function IntervalArithmetic.eigenbox(M::LinearAlgebra.Hermitian{Complex{Interval{T}}, Matrix{Complex{Interval{T}}}}, + method::IntervalArithmetic.AbstractIntervalEigenSolver) where {T} + A = real(M) + B = imag(M) + return IntervalArithmetic.eigenbox(LinearAlgebra.Symmetric([A B'; B A]), method) +end + +IntervalArithmetic.eigenbox(A) = IntervalArithmetic.eigenbox(A, IntervalArithmetic.Rohn()) + # matrix inversion # note: use the contraction mapping theorem, only works when the entries of A have small radii diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl index 7978085b..3bb77923 100644 --- a/src/IntervalArithmetic.jl +++ b/src/IntervalArithmetic.jl @@ -293,4 +293,15 @@ bareinterval(::Type{BigFloat}, a::AbstractIrrational) = _unsafe_bareinterval(Big return :($x) # set body of the function to return the precomputed result end +# eigenvalue enclosure types and functions (implementations in LinearAlgebra extension) + +abstract type AbstractIntervalEigenSolver end + +struct Rohn <: AbstractIntervalEigenSolver end +struct Hertz <: AbstractIntervalEigenSolver end + +function eigenbox end + + export eigenbox, AbstractIntervalEigenSolver, Rohn, Hertz + end From bbe8d1dfeba2d7c69188349891e297007010cefb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 15 Apr 2026 17:39:27 +0200 Subject: [PATCH 02/10] Attribution --- ext/IntervalArithmeticLinearAlgebraExt.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/ext/IntervalArithmeticLinearAlgebraExt.jl b/ext/IntervalArithmeticLinearAlgebraExt.jl index 07ed8580..145a311a 100644 --- a/ext/IntervalArithmeticLinearAlgebraExt.jl +++ b/ext/IntervalArithmeticLinearAlgebraExt.jl @@ -180,6 +180,12 @@ function _fold_conjugate!(::Type{<:Interval}, λ, v) end # eigenvalue enclosure (eigenbox) +# The Orthants iterator, eigenbox (Rohn and Hertz methods), and the reduction of +# general/complex/Hermitian matrices to the symmetric case are derived from +# IntervalLinearAlgebra.jl (https://github.com/JuliaIntervals/IntervalLinearAlgebra.jl) +# Copyright (c) 2021 Luca Ferranti, MIT License. +# See: Hladík, Daney, Tsigaridas, "Bounds on real eigenvalues and singular values +# of interval matrices", APNUM 2013 (https://doi.org/10.1016/j.apnum.2012.09.003). struct Orthants n::Int From 37db3eb242d34d2d45ba76a552f1724399dc8474 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 15 Apr 2026 17:49:13 +0200 Subject: [PATCH 03/10] Rename to eigen --- ext/IntervalArithmeticLinearAlgebraExt.jl | 60 +++++++++++------------ src/IntervalArithmetic.jl | 11 ----- 2 files changed, 28 insertions(+), 43 deletions(-) diff --git a/ext/IntervalArithmeticLinearAlgebraExt.jl b/ext/IntervalArithmeticLinearAlgebraExt.jl index 145a311a..75201c3a 100644 --- a/ext/IntervalArithmeticLinearAlgebraExt.jl +++ b/ext/IntervalArithmeticLinearAlgebraExt.jl @@ -151,8 +151,8 @@ function LinearAlgebra.eigen!(A::AbstractMatrix{<:RealOrComplexI}; permute::Bool _fold_conjugate!(eltype(A), true_λ, true_v) foreach(j -> true_v[:,j] .*= interval(ref_scale[j]), 1:n) else - # contraction mapping failed; fall back to eigenbox for eigenvalue bounds - box = IntervalArithmetic.eigenbox(A) + # contraction mapping failed; fall back to eigenvalue enclosure + box = _eigbox(A) true_λ = fill(box, n) true_v = fill(nai(eltype(v_bar)), n, n) end @@ -179,35 +179,26 @@ function _fold_conjugate!(::Type{<:Interval}, λ, v) return λ, v end -# eigenvalue enclosure (eigenbox) -# The Orthants iterator, eigenbox (Rohn and Hertz methods), and the reduction of +# eigenvalue enclosure for the eigen fallback +# The Rohn method, Hertz method, Orthants iterator, and the reduction of # general/complex/Hermitian matrices to the symmetric case are derived from # IntervalLinearAlgebra.jl (https://github.com/JuliaIntervals/IntervalLinearAlgebra.jl) # Copyright (c) 2021 Luca Ferranti, MIT License. # See: Hladík, Daney, Tsigaridas, "Bounds on real eigenvalues and singular values # of interval matrices", APNUM 2013 (https://doi.org/10.1016/j.apnum.2012.09.003). -struct Orthants +struct _Orthants n::Int end -Base.eltype(::Type{Orthants}) = Vector{Int} -Base.length(O::Orthants) = 2^(O.n) +Base.length(O::_Orthants) = 2^(O.n) -function Base.iterate(O::Orthants, state=1) +function Base.iterate(O::_Orthants, state=1) state > 2 ^ O.n && return nothing vec = -2*digits(state-1, base=2, pad=O.n) .+ 1 return (vec, state+1) end -function Base.getindex(O::Orthants, i::Int) - 1 <= i <= length(O) || throw(BoundsError(O, i)) - return -2*digits(i-1, base=2, pad=O.n) .+ 1 -end - -Base.firstindex(O::Orthants) = 1 -Base.lastindex(O::Orthants) = length(O) - # rigorous eigmax/eigmin via interval eigvals function _interval_eigmax(A::LinearAlgebra.Symmetric{<:Real, <:AbstractMatrix{<:Real}}) @@ -220,7 +211,9 @@ function _interval_eigmin(A::LinearAlgebra.Symmetric{<:Real, <:AbstractMatrix{<: return inf(minimum(real.(λs))) end -function IntervalArithmetic.eigenbox(A::LinearAlgebra.Symmetric{Interval{T}, Matrix{Interval{T}}}, ::IntervalArithmetic.Rohn) where {T} +# Rohn: fast eigenvalue enclosure for symmetric interval matrices + +function _eigbox_rohn(A::LinearAlgebra.Symmetric{Interval{T}, Matrix{Interval{T}}}) where {T} AΔ = LinearAlgebra.Symmetric(IntervalArithmetic.radius.(A)) Ac = LinearAlgebra.Symmetric(mid.(A)) @@ -230,14 +223,16 @@ function IntervalArithmetic.eigenbox(A::LinearAlgebra.Symmetric{Interval{T}, Mat return interval(λmin - ρ, λmax + ρ) end -function IntervalArithmetic.eigenbox(A::LinearAlgebra.Symmetric{Interval{T}, Matrix{Interval{T}}}, ::IntervalArithmetic.Hertz) where {T} +# Hertz: exact hull for symmetric interval matrices (exponential complexity) + +function _eigbox_hertz(A::LinearAlgebra.Symmetric{Interval{T}, Matrix{Interval{T}}}) where {T} n = LinearAlgebra.checksquare(A) Amax = Matrix{T}(undef, n, n) Amin = Matrix{T}(undef, n, n) λmin = T(Inf) λmax = T(-Inf) - @inbounds for z in Orthants(n) + @inbounds for z in _Orthants(n) first(z) < 0 && continue for j in 1:n for i in 1:j @@ -259,35 +254,36 @@ function IntervalArithmetic.eigenbox(A::LinearAlgebra.Symmetric{Interval{T}, Mat return interval(λmin, λmax) end -function IntervalArithmetic.eigenbox(A::AbstractMatrix{Interval{T}}, - method::IntervalArithmetic.AbstractIntervalEigenSolver) where {T} - λ = IntervalArithmetic.eigenbox(LinearAlgebra.Symmetric(interval.(T, 0.5*(A + A'))), method) +# eigenvalue enclosure for general and complex interval matrices (reduced to symmetric case) + +function _eigbox(A::LinearAlgebra.Symmetric{Interval{T}, Matrix{Interval{T}}}) where {T} + return _eigbox_rohn(A) +end + +function _eigbox(A::AbstractMatrix{Interval{T}}) where {T} + λ = _eigbox(LinearAlgebra.Symmetric(interval.(T, 0.5*(A + A')))) n = size(A, 1) S = 0.5*(A - A') - μ = IntervalArithmetic.eigenbox(LinearAlgebra.Symmetric(interval.(T, [zeros(Interval{T}, n, n) S'; S zeros(Interval{T}, n, n)])), method) + μ = _eigbox(LinearAlgebra.Symmetric(interval.(T, [zeros(Interval{T}, n, n) S'; S zeros(Interval{T}, n, n)]))) return λ + μ*im end -function IntervalArithmetic.eigenbox(M::AbstractMatrix{Complex{Interval{T}}}, - method::IntervalArithmetic.AbstractIntervalEigenSolver) where {T} +function _eigbox(M::AbstractMatrix{Complex{Interval{T}}}) where {T} A = real.(M) B = imag.(M) - λ = IntervalArithmetic.eigenbox(LinearAlgebra.Symmetric(interval.(T, 0.5*[A+A' B'-B; B-B' A+A'])), method) - μ = IntervalArithmetic.eigenbox(LinearAlgebra.Symmetric(interval.(T, 0.5*[B+B' A-A'; A'-A B+B'])), method) + λ = _eigbox(LinearAlgebra.Symmetric(interval.(T, 0.5*[A+A' B'-B; B-B' A+A']))) + μ = _eigbox(LinearAlgebra.Symmetric(interval.(T, 0.5*[B+B' A-A'; A'-A B+B']))) return λ + μ*im end -function IntervalArithmetic.eigenbox(M::LinearAlgebra.Hermitian{Complex{Interval{T}}, Matrix{Complex{Interval{T}}}}, - method::IntervalArithmetic.AbstractIntervalEigenSolver) where {T} +function _eigbox(M::LinearAlgebra.Hermitian{Complex{Interval{T}}, Matrix{Complex{Interval{T}}}}) where {T} A = real(M) B = imag(M) - return IntervalArithmetic.eigenbox(LinearAlgebra.Symmetric([A B'; B A]), method) + return _eigbox(LinearAlgebra.Symmetric([A B'; B A])) end -IntervalArithmetic.eigenbox(A) = IntervalArithmetic.eigenbox(A, IntervalArithmetic.Rohn()) - # matrix inversion # note: use the contraction mapping theorem, only works when the entries of A have small radii diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl index 3bb77923..7978085b 100644 --- a/src/IntervalArithmetic.jl +++ b/src/IntervalArithmetic.jl @@ -293,15 +293,4 @@ bareinterval(::Type{BigFloat}, a::AbstractIrrational) = _unsafe_bareinterval(Big return :($x) # set body of the function to return the precomputed result end -# eigenvalue enclosure types and functions (implementations in LinearAlgebra extension) - -abstract type AbstractIntervalEigenSolver end - -struct Rohn <: AbstractIntervalEigenSolver end -struct Hertz <: AbstractIntervalEigenSolver end - -function eigenbox end - - export eigenbox, AbstractIntervalEigenSolver, Rohn, Hertz - end From c277c6931fb7d28462fe25df58d5ce5438396ed5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 16 Apr 2026 10:22:37 +0200 Subject: [PATCH 04/10] Fix --- ext/IntervalArithmeticLinearAlgebraExt.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/ext/IntervalArithmeticLinearAlgebraExt.jl b/ext/IntervalArithmeticLinearAlgebraExt.jl index 75201c3a..4963cd08 100644 --- a/ext/IntervalArithmeticLinearAlgebraExt.jl +++ b/ext/IntervalArithmeticLinearAlgebraExt.jl @@ -55,8 +55,9 @@ end # matrix eigenvalues -function LinearAlgebra.eigvals!(A::AbstractMatrix{<:Interval}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) +function LinearAlgebra.eigvals!(A::AbstractMatrix{<:Interval}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby, kwargs...) # note: this function does not overwrite `A` + # kwargs absorbs e.g. `alg` passed by Symmetric eigvals in Julia ≥ 1.12 λ = _eigvals(A, permute, scale, sortby) isreal(λ) && return real(λ) _fold_conjugate!(λ) @@ -64,7 +65,7 @@ function LinearAlgebra.eigvals!(A::AbstractMatrix{<:Interval}; permute::Bool=tru return λ end -LinearAlgebra.eigvals!(A::AbstractMatrix{<:Complex{<:Interval}}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = +LinearAlgebra.eigvals!(A::AbstractMatrix{<:Complex{<:Interval}}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby, kwargs...) = # note: this function does not overwrite `A` _eigvals(A, permute, scale, sortby) From 6c595ade477438b59dd57bca3c4ed675ff60d65f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 16 Apr 2026 10:30:20 +0200 Subject: [PATCH 05/10] Fix --- ext/IntervalArithmeticLinearAlgebraExt.jl | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/ext/IntervalArithmeticLinearAlgebraExt.jl b/ext/IntervalArithmeticLinearAlgebraExt.jl index 4963cd08..591c1fe2 100644 --- a/ext/IntervalArithmeticLinearAlgebraExt.jl +++ b/ext/IntervalArithmeticLinearAlgebraExt.jl @@ -55,9 +55,24 @@ end # matrix eigenvalues -function LinearAlgebra.eigvals!(A::AbstractMatrix{<:Interval}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby, kwargs...) +if isdefined(LinearAlgebra, :Algorithm) + struct IntervalEigen <: LinearAlgebra.Algorithm end +else + struct IntervalEigen end +end + +# Julia ≥ 1.12: override default_eigen_alg so that the Symmetric/Hermitian eigvals/eigen +# paths pass IntervalEigen() instead of a LAPACK algorithm. +if isdefined(LinearAlgebra, :default_eigen_alg) + LinearAlgebra.default_eigen_alg(::Union{ + LinearAlgebra.Symmetric{<:RealIntervalType}, + LinearAlgebra.Symmetric{<:Complex{<:RealIntervalType}}, + LinearAlgebra.Hermitian{<:RealIntervalType}, + LinearAlgebra.Hermitian{<:Complex{<:RealIntervalType}}}) = IntervalEigen() +end + +function LinearAlgebra.eigvals!(A::AbstractMatrix{<:Interval}; alg::IntervalEigen=IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) # note: this function does not overwrite `A` - # kwargs absorbs e.g. `alg` passed by Symmetric eigvals in Julia ≥ 1.12 λ = _eigvals(A, permute, scale, sortby) isreal(λ) && return real(λ) _fold_conjugate!(λ) @@ -65,7 +80,7 @@ function LinearAlgebra.eigvals!(A::AbstractMatrix{<:Interval}; permute::Bool=tru return λ end -LinearAlgebra.eigvals!(A::AbstractMatrix{<:Complex{<:Interval}}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby, kwargs...) = +LinearAlgebra.eigvals!(A::AbstractMatrix{<:Complex{<:Interval}}; alg::IntervalEigen=IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = # note: this function does not overwrite `A` _eigvals(A, permute, scale, sortby) From b007161e3eb2c34a3b9b9528760b32ff124a7b76 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 17 Apr 2026 09:09:24 +0200 Subject: [PATCH 06/10] Add explicit algorithms --- ext/IntervalArithmeticLinearAlgebraExt.jl | 120 ++++++++++++------ src/IntervalArithmetic.jl | 47 +++++++ .../ext/IntervalArithmeticLinearAlgebraExt.jl | 66 ++++++++++ 3 files changed, 195 insertions(+), 38 deletions(-) diff --git a/ext/IntervalArithmeticLinearAlgebraExt.jl b/ext/IntervalArithmeticLinearAlgebraExt.jl index 591c1fe2..64f148ee 100644 --- a/ext/IntervalArithmeticLinearAlgebraExt.jl +++ b/ext/IntervalArithmeticLinearAlgebraExt.jl @@ -53,25 +53,12 @@ function LinearAlgebra.opnormInf(A::AbstractMatrix{T}) where {T<:RealOrComplexI} return convert(Tnorm, nrm) end -# matrix eigenvalues +# matrix eigenvalues – algorithm types +# Types are defined in IntervalArithmetic; we use them here for dispatch. -if isdefined(LinearAlgebra, :Algorithm) - struct IntervalEigen <: LinearAlgebra.Algorithm end -else - struct IntervalEigen end -end - -# Julia ≥ 1.12: override default_eigen_alg so that the Symmetric/Hermitian eigvals/eigen -# paths pass IntervalEigen() instead of a LAPACK algorithm. -if isdefined(LinearAlgebra, :default_eigen_alg) - LinearAlgebra.default_eigen_alg(::Union{ - LinearAlgebra.Symmetric{<:RealIntervalType}, - LinearAlgebra.Symmetric{<:Complex{<:RealIntervalType}}, - LinearAlgebra.Hermitian{<:RealIntervalType}, - LinearAlgebra.Hermitian{<:Complex{<:RealIntervalType}}}) = IntervalEigen() -end +# Algorithm types (IntervalArithmetic.IntervalEigen etc.) are used with explicit prefix below. -function LinearAlgebra.eigvals!(A::AbstractMatrix{<:Interval}; alg::IntervalEigen=IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) +function LinearAlgebra.eigvals!(A::AbstractMatrix{<:Interval}; alg::IntervalArithmetic.AbstractIntervalEigenAlg=IntervalArithmetic.IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) # note: this function does not overwrite `A` λ = _eigvals(A, permute, scale, sortby) isreal(λ) && return real(λ) @@ -80,10 +67,28 @@ function LinearAlgebra.eigvals!(A::AbstractMatrix{<:Interval}; alg::IntervalEige return λ end -LinearAlgebra.eigvals!(A::AbstractMatrix{<:Complex{<:Interval}}; alg::IntervalEigen=IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = +LinearAlgebra.eigvals!(A::AbstractMatrix{<:Complex{<:Interval}}; alg::IntervalArithmetic.AbstractIntervalEigenAlg=IntervalArithmetic.IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = # note: this function does not overwrite `A` _eigvals(A, permute, scale, sortby) +# More specific methods for Symmetric/Hermitian wrappers to take priority over +# LinearAlgebra's methods on Julia ≥ 1.12 (which require alg::Algorithm). +const _SymHermInterval = Union{ + LinearAlgebra.Symmetric{<:RealIntervalType}, + LinearAlgebra.Symmetric{<:Complex{<:RealIntervalType}}, + LinearAlgebra.Hermitian{<:RealIntervalType}, + LinearAlgebra.Hermitian{<:Complex{<:RealIntervalType}}} + +function LinearAlgebra.eigvals(A::_SymHermInterval; alg::IntervalArithmetic.AbstractIntervalEigenAlg=IntervalArithmetic.IntervalEigen(), sortby::Union{Function,Nothing}=nothing) + sb = sortby === nothing ? LinearAlgebra.eigsortby : sortby + return LinearAlgebra.eigvals!(Matrix(A); alg, permute=true, scale=true, sortby=sb) +end + +function LinearAlgebra.eigvals!(A::_SymHermInterval; alg::IntervalArithmetic.AbstractIntervalEigenAlg=IntervalArithmetic.IntervalEigen(), sortby::Union{Function,Nothing}=nothing) + sb = sortby === nothing ? LinearAlgebra.eigsortby : sortby + return LinearAlgebra.eigvals!(Matrix(A); alg, permute=true, scale=true, sortby=sb) +end + function _eigvals(A, permute, scale, sortby) # Gershgorin circle theorem B = _similarity_transform(A, permute, scale, sortby) @@ -133,13 +138,45 @@ LinearAlgebra.det(A::AbstractMatrix{<:Interval}) = real(reduce(*, LinearAlgebra. LinearAlgebra.det(A::AbstractMatrix{<:Complex{<:Interval}}) = reduce(*, LinearAlgebra.eigvals(A)) # matrix eigendecomposition -# note: use the contraction mapping theorem, only works when the entries of A have small radii, and A has simple eigenvalues -LinearAlgebra.eigen(A::AbstractMatrix{<:RealOrComplexI}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = - LinearAlgebra.eigen!(A; permute, scale, sortby) +LinearAlgebra.eigen(A::AbstractMatrix{<:RealOrComplexI}; alg::IntervalArithmetic.AbstractIntervalEigenAlg=IntervalArithmetic.IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = + LinearAlgebra.eigen!(A; alg, permute, scale, sortby) -function LinearAlgebra.eigen!(A::AbstractMatrix{<:RealOrComplexI}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) +LinearAlgebra.eigen(A::_SymHermInterval; alg::IntervalArithmetic.AbstractIntervalEigenAlg=IntervalArithmetic.IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = + LinearAlgebra.eigen!(A; alg, permute, scale, sortby) + +function LinearAlgebra.eigen!(A::AbstractMatrix{<:RealOrComplexI}; alg::IntervalArithmetic.AbstractIntervalEigenAlg=IntervalArithmetic.IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) # note: this function does not overwrite `A` + true_λ, true_v = _eigen_dispatch(alg, A, permute, scale, sortby) + _ensure_ng_flag!(true_v, all(isguaranteed, A)) + return LinearAlgebra.Eigen(true_λ, true_v) +end + +# IntervalEigen (default): try contraction mapping, fall back to Rohn on failure +function _eigen_dispatch(::IntervalArithmetic.IntervalEigen, A, permute, scale, sortby) + contraction_ok, true_λ, true_v = _eigen_contraction(A, permute, scale, sortby) + contraction_ok && return true_λ, true_v + return _eigen_eigbox(A, _eigbox_rohn) +end + +# IntervalEigenContraction: contraction mapping only (nai if it fails) +function _eigen_dispatch(::IntervalArithmetic.IntervalEigenContraction, A, permute, scale, sortby) + contraction_ok, true_λ, true_v = _eigen_contraction(A, permute, scale, sortby) + contraction_ok && return true_λ, true_v + n = LinearAlgebra.checksquare(A) + T = eltype(A) + return fill(nai(T), n), fill(nai(T), n, n) +end + +# IntervalEigenRohn: Rohn eigenvalue enclosure +_eigen_dispatch(::IntervalArithmetic.IntervalEigenRohn, A, permute, scale, sortby) = _eigen_eigbox(A, _eigbox_rohn) + +# IntervalEigenHertz: Hertz exact hull +_eigen_dispatch(::IntervalArithmetic.IntervalEigenHertz, A, permute, scale, sortby) = _eigen_eigbox(A, _eigbox_hertz) + +# shared: contraction mapping theorem (tight bounds when radii are small and eigenvalues are simple) + +function _eigen_contraction(A, permute, scale, sortby) λ, v = LinearAlgebra.eigen!(mid.(A); permute, scale, sortby) n = length(λ) inds = [argmax(i -> abs(v[i,j]), 1:n) for j ∈ 1:n] @@ -166,14 +203,21 @@ function LinearAlgebra.eigen!(A::AbstractMatrix{<:RealOrComplexI}; permute::Bool true_v = interval.(v, r; format = :midpoint) _fold_conjugate!(eltype(A), true_λ, true_v) foreach(j -> true_v[:,j] .*= interval(ref_scale[j]), 1:n) + return true, true_λ, true_v else - # contraction mapping failed; fall back to eigenvalue enclosure - box = _eigbox(A) - true_λ = fill(box, n) - true_v = fill(nai(eltype(v_bar)), n, n) + return false, nothing, nothing end - _ensure_ng_flag!(true_v, all(isguaranteed, A)) - return LinearAlgebra.Eigen(true_λ, true_v) +end + +# shared: eigenvalue enclosure via eigbox (returns single enclosing interval for all eigenvalues, nai eigenvectors) + +function _eigen_eigbox(A, eigbox_method) + n = LinearAlgebra.checksquare(A) + box = _eigbox(A, eigbox_method) + v_bar = interval(mid.(A)) # just need the element type + true_λ = fill(box, n) + true_v = fill(nai(eltype(v_bar)), n, n) + return true_λ, true_v end _fold_conjugate!(::Type{<:ComplexI}, λ, v) = (λ, v) @@ -272,32 +316,32 @@ end # eigenvalue enclosure for general and complex interval matrices (reduced to symmetric case) -function _eigbox(A::LinearAlgebra.Symmetric{Interval{T}, Matrix{Interval{T}}}) where {T} - return _eigbox_rohn(A) +function _eigbox(A::LinearAlgebra.Symmetric{Interval{T}, Matrix{Interval{T}}}, method=_eigbox_rohn) where {T} + return method(A) end -function _eigbox(A::AbstractMatrix{Interval{T}}) where {T} - λ = _eigbox(LinearAlgebra.Symmetric(interval.(T, 0.5*(A + A')))) +function _eigbox(A::AbstractMatrix{Interval{T}}, method=_eigbox_rohn) where {T} + λ = _eigbox(LinearAlgebra.Symmetric(interval.(T, 0.5*(A + A'))), method) n = size(A, 1) S = 0.5*(A - A') - μ = _eigbox(LinearAlgebra.Symmetric(interval.(T, [zeros(Interval{T}, n, n) S'; S zeros(Interval{T}, n, n)]))) + μ = _eigbox(LinearAlgebra.Symmetric(interval.(T, [zeros(Interval{T}, n, n) S'; S zeros(Interval{T}, n, n)])), method) return λ + μ*im end -function _eigbox(M::AbstractMatrix{Complex{Interval{T}}}) where {T} +function _eigbox(M::AbstractMatrix{Complex{Interval{T}}}, method=_eigbox_rohn) where {T} A = real.(M) B = imag.(M) - λ = _eigbox(LinearAlgebra.Symmetric(interval.(T, 0.5*[A+A' B'-B; B-B' A+A']))) - μ = _eigbox(LinearAlgebra.Symmetric(interval.(T, 0.5*[B+B' A-A'; A'-A B+B']))) + λ = _eigbox(LinearAlgebra.Symmetric(interval.(T, 0.5*[A+A' B'-B; B-B' A+A'])), method) + μ = _eigbox(LinearAlgebra.Symmetric(interval.(T, 0.5*[B+B' A-A'; A'-A B+B'])), method) return λ + μ*im end -function _eigbox(M::LinearAlgebra.Hermitian{Complex{Interval{T}}, Matrix{Complex{Interval{T}}}}) where {T} +function _eigbox(M::LinearAlgebra.Hermitian{Complex{Interval{T}}, Matrix{Complex{Interval{T}}}}, method=_eigbox_rohn) where {T} A = real(M) B = imag(M) - return _eigbox(LinearAlgebra.Symmetric([A B'; B A])) + return _eigbox(LinearAlgebra.Symmetric([A B'; B A]), method) end # matrix inversion diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl index 7978085b..b56352ff 100644 --- a/src/IntervalArithmetic.jl +++ b/src/IntervalArithmetic.jl @@ -122,6 +122,53 @@ function configure_power(power::Symbol) return power end +# eigendecomposition algorithm types + +""" + AbstractIntervalEigenAlg + +Abstract supertype for interval eigendecomposition algorithms. + +See also: [`IntervalEigen`](@ref), [`IntervalEigenContraction`](@ref), +[`IntervalEigenRohn`](@ref), [`IntervalEigenHertz`](@ref). +""" +abstract type AbstractIntervalEigenAlg end + +""" + IntervalEigen + +Default eigendecomposition algorithm for interval matrices. Tries the contraction +mapping theorem first (tight bounds when radii are small and eigenvalues are simple), +and falls back to the Rohn eigenvalue enclosure if it fails. +""" +struct IntervalEigen <: AbstractIntervalEigenAlg end + +""" + IntervalEigenContraction + +Eigendecomposition via the contraction mapping theorem only. Returns `nai` if the +contraction fails (e.g. when interval radii are too large). +""" +struct IntervalEigenContraction <: AbstractIntervalEigenAlg end + +""" + IntervalEigenRohn + +Rohn eigenvalue enclosure (fast, O(n³)). Returns a single interval enclosing all +eigenvalues and `nai` eigenvectors. Valid for any interval radius. +""" +struct IntervalEigenRohn <: AbstractIntervalEigenAlg end + +""" + IntervalEigenHertz + +Hertz exact hull of eigenvalues (exponential O(2ⁿ) complexity). Returns a single +interval enclosing all eigenvalues and `nai` eigenvectors. +""" +struct IntervalEigenHertz <: AbstractIntervalEigenAlg end + +export AbstractIntervalEigenAlg, IntervalEigen, IntervalEigenContraction, IntervalEigenRohn, IntervalEigenHertz + """ MatMulMode{T} diff --git a/test/ext/IntervalArithmeticLinearAlgebraExt.jl b/test/ext/IntervalArithmeticLinearAlgebraExt.jl index 73335594..330fcf61 100644 --- a/test/ext/IntervalArithmeticLinearAlgebraExt.jl +++ b/test/ext/IntervalArithmeticLinearAlgebraExt.jl @@ -113,6 +113,72 @@ end @test isequal_interval(v3[1, 1], complex(interval(1), interval(0))) end +@testset "Eigendecomposition" begin + exact_eigs = LinearAlgebra.eigvals(LinearAlgebra.Symmetric([1.0 0.5; 0.5 2.0])) + + # small radius: contraction mapping succeeds, tight individual bounds + A_small = interval.([1.0 0.5; 0.5 2.0]) + r_small = LinearAlgebra.eigen(A_small) + @test in_interval(exact_eigs[1], real(r_small.values[1])) + @test in_interval(exact_eigs[2], real(r_small.values[2])) + @test !any(isnai, r_small.vectors) + + # moderate radius: contraction mapping may fail, fallback provides valid enclosure + A_mid = interval.([1.0 0.5; 0.5 2.0], 0.3; format = :midpoint) + r_mid = LinearAlgebra.eigen(A_mid) + @test in_interval(exact_eigs[1], real(r_mid.values[1])) + @test in_interval(exact_eigs[2], real(r_mid.values[2])) + + # large radius: contraction mapping fails, eigenvalue enclosure fallback + A_large = interval.([1.0 0.5; 0.5 2.0], 1.0; format = :midpoint) + r_large = LinearAlgebra.eigen(A_large) + @test in_interval(exact_eigs[1], real(r_large.values[1])) + @test in_interval(exact_eigs[2], real(r_large.values[2])) + @test all(isnai, r_large.vectors) + + # general (non-symmetric) matrix with large radius + M = [0.0 -1.0; 2.0 -0.5] + mid_eigs = LinearAlgebra.eigvals(M) + A_gen = interval.(M, 0.5; format = :midpoint) + r_gen = LinearAlgebra.eigen(A_gen) + for λ in mid_eigs + @test in_interval(real(λ), real(r_gen.values[1])) + @test in_interval(imag(λ), imag(r_gen.values[1])) + end + + # 3×3 symmetric with large radius + S3_mid = [2.0 1.0 0.0; 1.0 3.0 1.0; 0.0 1.0 4.0] + exact3 = LinearAlgebra.eigvals(LinearAlgebra.Symmetric(S3_mid)) + A3 = interval.(S3_mid, 0.5; format = :midpoint) + r3 = LinearAlgebra.eigen(A3) + for λ in exact3 + @test in_interval(λ, real(r3.values[1])) + end + + # explicit algorithm selection: IntervalEigenContraction on small-radius matrix + r_contraction = LinearAlgebra.eigen(A_small; alg = IntervalEigenContraction()) + @test in_interval(exact_eigs[1], real(r_contraction.values[1])) + @test in_interval(exact_eigs[2], real(r_contraction.values[2])) + @test !any(isnai, r_contraction.vectors) + + # explicit algorithm selection: IntervalEigenContraction on large-radius matrix → nai + r_contraction_large = LinearAlgebra.eigen(A_large; alg = IntervalEigenContraction()) + @test all(isnai, r_contraction_large.values) + @test all(isnai, r_contraction_large.vectors) + + # explicit algorithm selection: IntervalEigenRohn + r_rohn = LinearAlgebra.eigen(A_large; alg = IntervalEigenRohn()) + @test in_interval(exact_eigs[1], real(r_rohn.values[1])) + @test in_interval(exact_eigs[2], real(r_rohn.values[2])) + @test all(isnai, r_rohn.vectors) + + # explicit algorithm selection: IntervalEigenHertz + r_hertz = LinearAlgebra.eigen(A_large; alg = IntervalEigenHertz()) + @test in_interval(exact_eigs[1], real(r_hertz.values[1])) + @test in_interval(exact_eigs[2], real(r_hertz.values[2])) + @test all(isnai, r_hertz.vectors) +end + @testset "inv" begin A = interval.([1.0 2.0; 3.0 4.0]) iA = inv(A) From b1059b5da4c9dfbf96c3385655df624325aad420 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 17 Apr 2026 11:57:50 +0200 Subject: [PATCH 07/10] Improve docs --- docs/src/manual/usage.md | 22 ++++++++++++++++++++++ src/IntervalArithmetic.jl | 8 ++++---- 2 files changed, 26 insertions(+), 4 deletions(-) diff --git a/docs/src/manual/usage.md b/docs/src/manual/usage.md index d3de3e68..b228428d 100644 --- a/docs/src/manual/usage.md +++ b/docs/src/manual/usage.md @@ -181,6 +181,28 @@ and deal properly with the decorations of intervals. +### Linear algebra + +Interval matrices support standard linear algebra operations from the `LinearAlgebra` +standard library. The eigendecomposition (`eigen`, `eigvals`) uses several algorithms +that can be selected via the `alg` keyword argument: + +- [`IntervalEigen`](@ref IntervalArithmetic.IntervalEigen) (default): tries the contraction mapping theorem first, falls back to the Rohn enclosure. +- [`IntervalEigenContraction`](@ref IntervalArithmetic.IntervalEigenContraction): contraction mapping only. +- [`IntervalEigenRohn`](@ref IntervalArithmetic.IntervalEigenRohn): Rohn eigenvalue enclosure. +- [`IntervalEigenHertz`](@ref IntervalArithmetic.IntervalEigenHertz): Hertz exact hull (exponential complexity). + +All algorithms are subtypes of [`AbstractIntervalEigenAlg`](@ref IntervalArithmetic.AbstractIntervalEigenAlg). + +```@repl usage +using LinearAlgebra +A = interval.([1.0 0.5; 0.5 2.0]) +eigen(A) +eigen(A; alg = IntervalEigenRohn()) +``` + + + ### Custom interval bounds type A `BareInterval{T}` or `Interval{T}` have the restriction `T <: Union{Rational,AbstractFloat}` which is the parametric type for the bounds of the interval. Supposing one wishes to use their own numeric type `MyNumType <: Union{Rational,AbstractFloat}`, they must provide their own arithmetic operations (with correct rounding!). diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl index b56352ff..d43f56c1 100644 --- a/src/IntervalArithmetic.jl +++ b/src/IntervalArithmetic.jl @@ -146,8 +146,8 @@ struct IntervalEigen <: AbstractIntervalEigenAlg end """ IntervalEigenContraction -Eigendecomposition via the contraction mapping theorem only. Returns `nai` if the -contraction fails (e.g. when interval radii are too large). +Eigendecomposition via the contraction mapping theorem only. Returns [`nai`](@ref) +if the contraction fails (e.g. when interval radii are too large). """ struct IntervalEigenContraction <: AbstractIntervalEigenAlg end @@ -155,7 +155,7 @@ struct IntervalEigenContraction <: AbstractIntervalEigenAlg end IntervalEigenRohn Rohn eigenvalue enclosure (fast, O(n³)). Returns a single interval enclosing all -eigenvalues and `nai` eigenvectors. Valid for any interval radius. +eigenvalues and [`nai`](@ref) eigenvectors. Valid for any interval radius. """ struct IntervalEigenRohn <: AbstractIntervalEigenAlg end @@ -163,7 +163,7 @@ struct IntervalEigenRohn <: AbstractIntervalEigenAlg end IntervalEigenHertz Hertz exact hull of eigenvalues (exponential O(2ⁿ) complexity). Returns a single -interval enclosing all eigenvalues and `nai` eigenvectors. +interval enclosing all eigenvalues and [`nai`](@ref) eigenvectors. """ struct IntervalEigenHertz <: AbstractIntervalEigenAlg end From 682c3811acd57501c0caee284d10fa635ba302d4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sat, 2 May 2026 19:57:16 +0200 Subject: [PATCH 08/10] Simplify --- ext/IntervalArithmeticLinearAlgebraExt.jl | 33 +++++++++++------------ 1 file changed, 15 insertions(+), 18 deletions(-) diff --git a/ext/IntervalArithmeticLinearAlgebraExt.jl b/ext/IntervalArithmeticLinearAlgebraExt.jl index 64f148ee..8d75d155 100644 --- a/ext/IntervalArithmeticLinearAlgebraExt.jl +++ b/ext/IntervalArithmeticLinearAlgebraExt.jl @@ -53,10 +53,7 @@ function LinearAlgebra.opnormInf(A::AbstractMatrix{T}) where {T<:RealOrComplexI} return convert(Tnorm, nrm) end -# matrix eigenvalues – algorithm types -# Types are defined in IntervalArithmetic; we use them here for dispatch. - -# Algorithm types (IntervalArithmetic.IntervalEigen etc.) are used with explicit prefix below. +# matrix eigenvalues function LinearAlgebra.eigvals!(A::AbstractMatrix{<:Interval}; alg::IntervalArithmetic.AbstractIntervalEigenAlg=IntervalArithmetic.IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) # note: this function does not overwrite `A` @@ -147,20 +144,20 @@ LinearAlgebra.eigen(A::_SymHermInterval; alg::IntervalArithmetic.AbstractInterva function LinearAlgebra.eigen!(A::AbstractMatrix{<:RealOrComplexI}; alg::IntervalArithmetic.AbstractIntervalEigenAlg=IntervalArithmetic.IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) # note: this function does not overwrite `A` - true_λ, true_v = _eigen_dispatch(alg, A, permute, scale, sortby) + true_λ, true_v = _eigen(alg, A, permute, scale, sortby) _ensure_ng_flag!(true_v, all(isguaranteed, A)) return LinearAlgebra.Eigen(true_λ, true_v) end # IntervalEigen (default): try contraction mapping, fall back to Rohn on failure -function _eigen_dispatch(::IntervalArithmetic.IntervalEigen, A, permute, scale, sortby) +function _eigen(::IntervalArithmetic.IntervalEigen, A, permute, scale, sortby) contraction_ok, true_λ, true_v = _eigen_contraction(A, permute, scale, sortby) contraction_ok && return true_λ, true_v - return _eigen_eigbox(A, _eigbox_rohn) + return _eigen(IntervalArithmetic.IntervalEigenRohn(), A, permute, scale, sortby) end # IntervalEigenContraction: contraction mapping only (nai if it fails) -function _eigen_dispatch(::IntervalArithmetic.IntervalEigenContraction, A, permute, scale, sortby) +function _eigen(::IntervalArithmetic.IntervalEigenContraction, A, permute, scale, sortby) contraction_ok, true_λ, true_v = _eigen_contraction(A, permute, scale, sortby) contraction_ok && return true_λ, true_v n = LinearAlgebra.checksquare(A) @@ -168,12 +165,6 @@ function _eigen_dispatch(::IntervalArithmetic.IntervalEigenContraction, A, permu return fill(nai(T), n), fill(nai(T), n, n) end -# IntervalEigenRohn: Rohn eigenvalue enclosure -_eigen_dispatch(::IntervalArithmetic.IntervalEigenRohn, A, permute, scale, sortby) = _eigen_eigbox(A, _eigbox_rohn) - -# IntervalEigenHertz: Hertz exact hull -_eigen_dispatch(::IntervalArithmetic.IntervalEigenHertz, A, permute, scale, sortby) = _eigen_eigbox(A, _eigbox_hertz) - # shared: contraction mapping theorem (tight bounds when radii are small and eigenvalues are simple) function _eigen_contraction(A, permute, scale, sortby) @@ -211,9 +202,15 @@ end # shared: eigenvalue enclosure via eigbox (returns single enclosing interval for all eigenvalues, nai eigenvectors) -function _eigen_eigbox(A, eigbox_method) +function _eigen( + alg::Union{IntervalArithmetic.IntervalEigenHertz,IntervalArithmetic.IntervalEigenRohn}, + A, + permute, + scale, + sortby, +) n = LinearAlgebra.checksquare(A) - box = _eigbox(A, eigbox_method) + box = _eigbox(A, alg) v_bar = interval(mid.(A)) # just need the element type true_λ = fill(box, n) true_v = fill(nai(eltype(v_bar)), n, n) @@ -273,7 +270,7 @@ end # Rohn: fast eigenvalue enclosure for symmetric interval matrices -function _eigbox_rohn(A::LinearAlgebra.Symmetric{Interval{T}, Matrix{Interval{T}}}) where {T} +function _eigbox(A::LinearAlgebra.Symmetric{Interval{T}, Matrix{Interval{T}}}, ::IntervalArithmetic.IntervalEigenRohn) where {T} AΔ = LinearAlgebra.Symmetric(IntervalArithmetic.radius.(A)) Ac = LinearAlgebra.Symmetric(mid.(A)) @@ -285,7 +282,7 @@ end # Hertz: exact hull for symmetric interval matrices (exponential complexity) -function _eigbox_hertz(A::LinearAlgebra.Symmetric{Interval{T}, Matrix{Interval{T}}}) where {T} +function _eigbox(A::LinearAlgebra.Symmetric{Interval{T}, Matrix{Interval{T}}}, ::IntervalArithmetic.IntervalEigenHertz) where {T} n = LinearAlgebra.checksquare(A) Amax = Matrix{T}(undef, n, n) Amin = Matrix{T}(undef, n, n) From c6790a87f8e2326b0d365ff4fd9a33f118c970f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sun, 1 Nov 2026 11:31:26 +0100 Subject: [PATCH 09/10] Address review comments MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Rohn bounds rounded outward - propagate NG flag to eigenvalues, use complex(λ, μ) - convert 1//2 to interval instead of 0.5 - block non-rigorous LinearAlgebra matrix functions of Symmetric/Hermitian interval matrices - run the contraction on Matrix(A) for Symmetric/Hermitian wrappers - fall back to Rohn when mid.(DF) is singular (repeated eigenvalues) - return NaI for unbounded entries - eigvals uses the alg keyword - move tests to test/ext and update the export list --- ext/IntervalArithmeticLinearAlgebraExt.jl | 219 +++++++++++------- src/IntervalArithmetic.jl | 7 +- test/IntervalArithmetic.jl | 2 + .../ext/IntervalArithmeticLinearAlgebraExt.jl | 54 ++++- 4 files changed, 190 insertions(+), 92 deletions(-) diff --git a/ext/IntervalArithmeticLinearAlgebraExt.jl b/ext/IntervalArithmeticLinearAlgebraExt.jl index 8d75d155..b313e5bb 100644 --- a/ext/IntervalArithmeticLinearAlgebraExt.jl +++ b/ext/IntervalArithmeticLinearAlgebraExt.jl @@ -55,38 +55,69 @@ end # matrix eigenvalues -function LinearAlgebra.eigvals!(A::AbstractMatrix{<:Interval}; alg::IntervalArithmetic.AbstractIntervalEigenAlg=IntervalArithmetic.IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) +const _EigenAlg = IntervalArithmetic.AbstractIntervalEigenAlg + +function LinearAlgebra.eigvals!(A::AbstractMatrix{<:Interval}; alg::_EigenAlg=IntervalArithmetic.IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) # note: this function does not overwrite `A` - λ = _eigvals(A, permute, scale, sortby) - isreal(λ) && return real(λ) - _fold_conjugate!(λ) - isreal(λ) && return real(λ) + λ = _eigvals(alg, A, permute, scale, sortby) + _ensure_ng_flag!(λ, all(isguaranteed, A)) return λ end -LinearAlgebra.eigvals!(A::AbstractMatrix{<:Complex{<:Interval}}; alg::IntervalArithmetic.AbstractIntervalEigenAlg=IntervalArithmetic.IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = +function LinearAlgebra.eigvals!(A::AbstractMatrix{<:Complex{<:Interval}}; alg::_EigenAlg=IntervalArithmetic.IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) # note: this function does not overwrite `A` - _eigvals(A, permute, scale, sortby) + λ = _eigvals(alg, A, permute, scale, sortby) + _ensure_ng_flag!(λ, all(isguaranteed, A)) + return λ +end -# More specific methods for Symmetric/Hermitian wrappers to take priority over -# LinearAlgebra's methods on Julia ≥ 1.12 (which require alg::Algorithm). +# more specific methods for `Symmetric` and `Hermitian` wrappers to take priority over +# the LinearAlgebra methods, which require `alg::LinearAlgebra.Algorithm` on Julia ≥ 1.12 const _SymHermInterval = Union{ - LinearAlgebra.Symmetric{<:RealIntervalType}, - LinearAlgebra.Symmetric{<:Complex{<:RealIntervalType}}, - LinearAlgebra.Hermitian{<:RealIntervalType}, - LinearAlgebra.Hermitian{<:Complex{<:RealIntervalType}}} + LinearAlgebra.Symmetric{<:Interval}, + LinearAlgebra.Symmetric{<:Complex{<:Interval}}, + LinearAlgebra.Hermitian{<:Interval}, + LinearAlgebra.Hermitian{<:Complex{<:Interval}}} -function LinearAlgebra.eigvals(A::_SymHermInterval; alg::IntervalArithmetic.AbstractIntervalEigenAlg=IntervalArithmetic.IntervalEigen(), sortby::Union{Function,Nothing}=nothing) - sb = sortby === nothing ? LinearAlgebra.eigsortby : sortby - return LinearAlgebra.eigvals!(Matrix(A); alg, permute=true, scale=true, sortby=sb) +LinearAlgebra.eigvals(A::_SymHermInterval; alg::_EigenAlg=IntervalArithmetic.IntervalEigen(), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = + LinearAlgebra.eigvals!(Matrix(A); alg, sortby) + +LinearAlgebra.eigvals!(A::_SymHermInterval; alg::_EigenAlg=IntervalArithmetic.IntervalEigen(), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = + LinearAlgebra.eigvals!(Matrix(A); alg, sortby) + +# IntervalEigen (default): Gershgorin circle theorem, falls back to Rohn on failure +function _eigvals(::IntervalArithmetic.IntervalEigen, A, permute, scale, sortby) + λ = _eigvals_gershgorin(A, permute, scale, sortby) + all(isbounded, λ) && return λ + return _eigvals(IntervalArithmetic.IntervalEigenRohn(), A, permute, scale, sortby) end -function LinearAlgebra.eigvals!(A::_SymHermInterval; alg::IntervalArithmetic.AbstractIntervalEigenAlg=IntervalArithmetic.IntervalEigen(), sortby::Union{Function,Nothing}=nothing) - sb = sortby === nothing ? LinearAlgebra.eigsortby : sortby - return LinearAlgebra.eigvals!(Matrix(A); alg, permute=true, scale=true, sortby=sb) +_eigvals(alg::IntervalArithmetic.IntervalEigenContraction, A, permute, scale, sortby) = + first(_eigen(alg, A, permute, scale, sortby)) + +function _eigvals(alg::Union{IntervalArithmetic.IntervalEigenRohn,IntervalArithmetic.IntervalEigenHertz}, A, permute, scale, sortby) + n = LinearAlgebra.checksquare(A) + _isbounded_matrix(A) || return fill(nai(complex(eltype(A))), n) + return fill(_eigbox(A, alg), n) end -function _eigvals(A, permute, scale, sortby) +function _eigvals_gershgorin(A::AbstractMatrix{<:Interval}, permute, scale, sortby) + _isbounded_matrix(A) || return fill(nai(eltype(A)), LinearAlgebra.checksquare(A)) + λ = _gershgorin(A, permute, scale, sortby) + isreal(λ) && return real(λ) + _fold_conjugate!(λ) + isreal(λ) && return real(λ) + return λ +end + +function _eigvals_gershgorin(A::AbstractMatrix{<:Complex{<:Interval}}, permute, scale, sortby) + _isbounded_matrix(A) || return fill(nai(eltype(A)), LinearAlgebra.checksquare(A)) + return _gershgorin(A, permute, scale, sortby) +end + +_isbounded_matrix(A) = all(x -> isbounded(real(x)) & isbounded(imag(x)), A) + +function _gershgorin(A, permute, scale, sortby) # Gershgorin circle theorem B = _similarity_transform(A, permute, scale, sortby) v = LinearAlgebra.diag(B) @@ -136,36 +167,48 @@ LinearAlgebra.det(A::AbstractMatrix{<:Complex{<:Interval}}) = reduce(*, LinearAl # matrix eigendecomposition -LinearAlgebra.eigen(A::AbstractMatrix{<:RealOrComplexI}; alg::IntervalArithmetic.AbstractIntervalEigenAlg=IntervalArithmetic.IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = +LinearAlgebra.eigen(A::AbstractMatrix{<:RealOrComplexI}; alg::_EigenAlg=IntervalArithmetic.IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = LinearAlgebra.eigen!(A; alg, permute, scale, sortby) -LinearAlgebra.eigen(A::_SymHermInterval; alg::IntervalArithmetic.AbstractIntervalEigenAlg=IntervalArithmetic.IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = - LinearAlgebra.eigen!(A; alg, permute, scale, sortby) +LinearAlgebra.eigen(A::_SymHermInterval; alg::_EigenAlg=IntervalArithmetic.IntervalEigen(), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = + LinearAlgebra.eigen!(Matrix(A); alg, sortby) + +LinearAlgebra.eigen!(A::_SymHermInterval; alg::_EigenAlg=IntervalArithmetic.IntervalEigen(), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = + LinearAlgebra.eigen!(Matrix(A); alg, sortby) -function LinearAlgebra.eigen!(A::AbstractMatrix{<:RealOrComplexI}; alg::IntervalArithmetic.AbstractIntervalEigenAlg=IntervalArithmetic.IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) +function LinearAlgebra.eigen!(A::AbstractMatrix{<:RealOrComplexI}; alg::_EigenAlg=IntervalArithmetic.IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) # note: this function does not overwrite `A` true_λ, true_v = _eigen(alg, A, permute, scale, sortby) - _ensure_ng_flag!(true_v, all(isguaranteed, A)) + ng_flag = all(isguaranteed, A) + _ensure_ng_flag!(true_λ, ng_flag) + _ensure_ng_flag!(true_v, ng_flag) return LinearAlgebra.Eigen(true_λ, true_v) end -# IntervalEigen (default): try contraction mapping, fall back to Rohn on failure +# IntervalEigen (default): contraction mapping theorem, falls back to Rohn on failure function _eigen(::IntervalArithmetic.IntervalEigen, A, permute, scale, sortby) + _isbounded_matrix(A) || return _eigen_nai(A) contraction_ok, true_λ, true_v = _eigen_contraction(A, permute, scale, sortby) contraction_ok && return true_λ, true_v return _eigen(IntervalArithmetic.IntervalEigenRohn(), A, permute, scale, sortby) end -# IntervalEigenContraction: contraction mapping only (nai if it fails) +# IntervalEigenContraction: contraction mapping theorem only (NaI if it fails) function _eigen(::IntervalArithmetic.IntervalEigenContraction, A, permute, scale, sortby) + _isbounded_matrix(A) || return _eigen_nai(A) contraction_ok, true_λ, true_v = _eigen_contraction(A, permute, scale, sortby) contraction_ok && return true_λ, true_v + return _eigen_nai(A) +end + +function _eigen_nai(A) n = LinearAlgebra.checksquare(A) - T = eltype(A) + T = complex(eltype(A)) return fill(nai(T), n), fill(nai(T), n, n) end -# shared: contraction mapping theorem (tight bounds when radii are small and eigenvalues are simple) +# contraction mapping theorem: tight bounds when the entries of A have small radii, and +# A has simple eigenvalues function _eigen_contraction(A, permute, scale, sortby) λ, v = LinearAlgebra.eigen!(mid.(A); permute, scale, sortby) @@ -183,7 +226,10 @@ function _eigen_contraction(A, permute, scale, sortby) DF[n+(j-1)*n+1:n+j*n,n+(j-1)*n+1:n+j*n] .= A - LinearAlgebra.UniformScaling(λ_bar[j]) end - approx_DF⁻¹ = interval(inv(mid.(DF))) + # `mid.(DF)` is singular when `mid.(A)` has repeated eigenvalues + lu_mid_DF = LinearAlgebra.lu(mid.(DF); check = false) + LinearAlgebra.issuccess(lu_mid_DF) || return false, nothing, nothing + approx_DF⁻¹ = interval(inv(lu_mid_DF)) Y = LinearAlgebra.norm(approx_DF⁻¹ * F, Inf) Z₁ = LinearAlgebra.opnorm(approx_DF⁻¹ * DF - interval(LinearAlgebra.I), Inf) Z₂ = interval(2) * LinearAlgebra.opnorm(approx_DF⁻¹, Inf) @@ -200,23 +246,6 @@ function _eigen_contraction(A, permute, scale, sortby) end end -# shared: eigenvalue enclosure via eigbox (returns single enclosing interval for all eigenvalues, nai eigenvectors) - -function _eigen( - alg::Union{IntervalArithmetic.IntervalEigenHertz,IntervalArithmetic.IntervalEigenRohn}, - A, - permute, - scale, - sortby, -) - n = LinearAlgebra.checksquare(A) - box = _eigbox(A, alg) - v_bar = interval(mid.(A)) # just need the element type - true_λ = fill(box, n) - true_v = fill(nai(eltype(v_bar)), n, n) - return true_λ, true_v -end - _fold_conjugate!(::Type{<:ComplexI}, λ, v) = (λ, v) function _fold_conjugate!(::Type{<:Interval}, λ, v) @@ -236,7 +265,14 @@ function _fold_conjugate!(::Type{<:Interval}, λ, v) return λ, v end -# eigenvalue enclosure for the eigen fallback +# eigenvalue enclosure: a single interval enclosing all the eigenvalues, NaI eigenvectors + +function _eigen(alg::Union{IntervalArithmetic.IntervalEigenRohn,IntervalArithmetic.IntervalEigenHertz}, A, permute, scale, sortby) + λ = _eigvals(alg, A, permute, scale, sortby) + n = length(λ) + return λ, fill(nai(complex(eltype(A))), n, n) +end + # The Rohn method, Hertz method, Orthants iterator, and the reduction of # general/complex/Hermitian matrices to the symmetric case are derived from # IntervalLinearAlgebra.jl (https://github.com/JuliaIntervals/IntervalLinearAlgebra.jl) @@ -256,33 +292,26 @@ function Base.iterate(O::_Orthants, state=1) return (vec, state+1) end -# rigorous eigmax/eigmin via interval eigvals +# rigorous bounds on the extreme eigenvalues of a symmetric matrix of floats +# note: call the Gershgorin enclosure directly to avoid the fallback to Rohn -function _interval_eigmax(A::LinearAlgebra.Symmetric{<:Real, <:AbstractMatrix{<:Real}}) - λs = LinearAlgebra.eigvals(LinearAlgebra.Symmetric(interval.(A))) - return sup(maximum(real.(λs))) -end +_eigmax(A::AbstractMatrix{<:AbstractFloat}) = + maximum(x -> sup(real(x)), _eigvals_gershgorin(interval(A), true, true, LinearAlgebra.eigsortby)) -function _interval_eigmin(A::LinearAlgebra.Symmetric{<:Real, <:AbstractMatrix{<:Real}}) - λs = LinearAlgebra.eigvals(LinearAlgebra.Symmetric(interval.(A))) - return inf(minimum(real.(λs))) -end +_eigmin(A::AbstractMatrix{<:AbstractFloat}) = + minimum(x -> inf(real(x)), _eigvals_gershgorin(interval(A), true, true, LinearAlgebra.eigsortby)) # Rohn: fast eigenvalue enclosure for symmetric interval matrices -function _eigbox(A::LinearAlgebra.Symmetric{Interval{T}, Matrix{Interval{T}}}, ::IntervalArithmetic.IntervalEigenRohn) where {T} - AΔ = LinearAlgebra.Symmetric(IntervalArithmetic.radius.(A)) - Ac = LinearAlgebra.Symmetric(mid.(A)) - - ρ = _interval_eigmax(AΔ) - λmax = _interval_eigmax(Ac) - λmin = _interval_eigmin(Ac) - return interval(λmin - ρ, λmax + ρ) +function _eigbox(A::LinearAlgebra.Symmetric{Interval{T},Matrix{Interval{T}}}, ::IntervalArithmetic.IntervalEigenRohn) where {T} + Ac = mid.(Matrix(A)) + ρ = _eigmax(radius.(Matrix(A))) + return interval(_eigmin(Ac), _eigmax(Ac)) + interval(-ρ, ρ) end # Hertz: exact hull for symmetric interval matrices (exponential complexity) -function _eigbox(A::LinearAlgebra.Symmetric{Interval{T}, Matrix{Interval{T}}}, ::IntervalArithmetic.IntervalEigenHertz) where {T} +function _eigbox(A::LinearAlgebra.Symmetric{Interval{T},Matrix{Interval{T}}}, ::IntervalArithmetic.IntervalEigenHertz) where {T} n = LinearAlgebra.checksquare(A) Amax = Matrix{T}(undef, n, n) Amin = Matrix{T}(undef, n, n) @@ -292,7 +321,7 @@ function _eigbox(A::LinearAlgebra.Symmetric{Interval{T}, Matrix{Interval{T}}}, : @inbounds for z in _Orthants(n) first(z) < 0 && continue for j in 1:n - for i in 1:j + for i in 1:n if z[i] == z[j] Amax[i, j] = sup(A[i, j]) Amin[i, j] = inf(A[i, j]) @@ -303,44 +332,58 @@ function _eigbox(A::LinearAlgebra.Symmetric{Interval{T}, Matrix{Interval{T}}}, : end end - candmax = _interval_eigmax(LinearAlgebra.Symmetric(Amax)) - candmin = _interval_eigmin(LinearAlgebra.Symmetric(Amin)) - λmin = min(λmin, candmin) - λmax = max(λmax, candmax) + λmin = min(λmin, _eigmin(Amin)) + λmax = max(λmax, _eigmax(Amax)) end return interval(λmin, λmax) end # eigenvalue enclosure for general and complex interval matrices (reduced to symmetric case) -function _eigbox(A::LinearAlgebra.Symmetric{Interval{T}, Matrix{Interval{T}}}, method=_eigbox_rohn) where {T} - return method(A) -end - -function _eigbox(A::AbstractMatrix{Interval{T}}, method=_eigbox_rohn) where {T} - λ = _eigbox(LinearAlgebra.Symmetric(interval.(T, 0.5*(A + A'))), method) +function _eigbox(A::AbstractMatrix{Interval{T}}, alg) where {T} + all(isequal_interval.(A, transpose(A))) && return _eigbox(LinearAlgebra.Symmetric(A), alg) + half = interval(T, 1//2) + λ = _eigbox(LinearAlgebra.Symmetric(half * (A + A')), alg) n = size(A, 1) - S = 0.5*(A - A') - μ = _eigbox(LinearAlgebra.Symmetric(interval.(T, [zeros(Interval{T}, n, n) S'; S zeros(Interval{T}, n, n)])), method) + S = half * (A - A') + Z = zeros(Interval{T}, n, n) + μ = _eigbox(LinearAlgebra.Symmetric([Z S'; S Z]), alg) - return λ + μ*im + return complex(λ, μ) end -function _eigbox(M::AbstractMatrix{Complex{Interval{T}}}, method=_eigbox_rohn) where {T} +function _eigbox(M::AbstractMatrix{Complex{Interval{T}}}, alg) where {T} A = real.(M) B = imag.(M) - λ = _eigbox(LinearAlgebra.Symmetric(interval.(T, 0.5*[A+A' B'-B; B-B' A+A'])), method) - μ = _eigbox(LinearAlgebra.Symmetric(interval.(T, 0.5*[B+B' A-A'; A'-A B+B'])), method) - return λ + μ*im + # Hermitian: the eigenvalues of `M` are those of the real symmetric matrix `[A B'; B A]` + all(isequal_interval.(M, adjoint(M))) && return _eigbox(LinearAlgebra.Symmetric([A B'; B A]), alg) + half = interval(T, 1//2) + λ = _eigbox(LinearAlgebra.Symmetric(half * [A+A' B'-B; B-B' A+A']), alg) + μ = _eigbox(LinearAlgebra.Symmetric(half * [B+B' A-A'; A'-A B+B']), alg) + return complex(λ, μ) end -function _eigbox(M::LinearAlgebra.Hermitian{Complex{Interval{T}}, Matrix{Complex{Interval{T}}}}, method=_eigbox_rohn) where {T} - A = real(M) - B = imag(M) - return _eigbox(LinearAlgebra.Symmetric([A B'; B A]), method) +# matrix functions of `Symmetric` and `Hermitian` interval matrices +# note: the LinearAlgebra methods rely on `eigen`, but assume that the eigenvectors are +# orthonormal and that the eigenvalues are accurate up to roundoff errors (e.g. `sqrt` +# rounds slightly negative eigenvalues to zero), which is not rigorous + +Base.exp(A::_SymHermInterval) = exp(Matrix(A)) +Base.log(A::_SymHermInterval) = log(Matrix(A)) + +for f ∈ (:sqrt, :cbrt, :cis, :cos, :sin, :tan, :cosh, :sinh, :tanh, :acos, :asin, :atan, :acosh, :asinh, :atanh, :sincos) + @eval Base.$f(::_SymHermInterval; kwargs...) = _throw_matrix_function($f) end +Base.:^(A::_SymHermInterval, p::Integer) = Matrix(A)^p +Base.:^(A::_SymHermInterval, p::Real) = isinteger(p) ? Matrix(A)^Integer(p) : _throw_matrix_function(^) + +LinearAlgebra.svd(::_SymHermInterval; full::Bool=false) = _throw_matrix_function(LinearAlgebra.svd) +LinearAlgebra.svdvals!(::_SymHermInterval) = _throw_matrix_function(LinearAlgebra.svdvals!) + +_throw_matrix_function(f) = throw(ArgumentError("`$f` is not supported for `Symmetric` or `Hermitian` interval matrices")) + # matrix inversion # note: use the contraction mapping theorem, only works when the entries of A have small radii diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl index d43f56c1..b94733ad 100644 --- a/src/IntervalArithmetic.jl +++ b/src/IntervalArithmetic.jl @@ -137,9 +137,10 @@ abstract type AbstractIntervalEigenAlg end """ IntervalEigen -Default eigendecomposition algorithm for interval matrices. Tries the contraction -mapping theorem first (tight bounds when radii are small and eigenvalues are simple), -and falls back to the Rohn eigenvalue enclosure if it fails. +Default eigendecomposition algorithm for interval matrices. For `eigen`, tries the +contraction mapping theorem first (tight bounds when radii are small and eigenvalues are +simple). For `eigvals`, uses the Gershgorin circle theorem after a similarity transform. +In both cases, falls back to the Rohn eigenvalue enclosure if it fails. """ struct IntervalEigen <: AbstractIntervalEigenAlg end diff --git a/test/IntervalArithmetic.jl b/test/IntervalArithmetic.jl index 52025f91..a0e2adae 100644 --- a/test/IntervalArithmetic.jl +++ b/test/IntervalArithmetic.jl @@ -32,6 +32,8 @@ end :mag, :mid, :midradius, :mig, :mince, :mince!, :nai, :numtype, :overlap, :pieces, :pow, :pown, :precedes, :radius, :rootn, :sample, :setdisplay, :strictprecedes, :sup, :trv, :union_interval] + append!(expected, [:AbstractIntervalEigenAlg, :IntervalEigen, :IntervalEigenContraction, + :IntervalEigenHertz, :IntervalEigenRohn]) @test sort(names(IntervalArithmetic)) == sort(expected) @test haskey(Base.Docs.meta(IntervalArithmetic), diff --git a/test/ext/IntervalArithmeticLinearAlgebraExt.jl b/test/ext/IntervalArithmeticLinearAlgebraExt.jl index 330fcf61..e439d6be 100644 --- a/test/ext/IntervalArithmeticLinearAlgebraExt.jl +++ b/test/ext/IntervalArithmeticLinearAlgebraExt.jl @@ -94,7 +94,7 @@ end @test all(B .=== A) W = A .+ interval(-1, 1) - EW = eigen(W) + EW = eigen(W; alg = IntervalEigenContraction()) @test all(isnai, EW.values) @test all(isnai, EW.vectors) @@ -177,6 +177,58 @@ end @test in_interval(exact_eigs[1], real(r_hertz.values[1])) @test in_interval(exact_eigs[2], real(r_hertz.values[2])) @test all(isnai, r_hertz.vectors) + + # `eigvals` uses the algorithm + λ_rohn = eigvals(A_large; alg = IntervalEigenRohn()) + @test all(λ -> issubset_interval(interval(-1.2, 4.2), λ), λ_rohn) + λ_hertz = eigvals(A_large; alg = IntervalEigenHertz()) + @test all(λ -> issubset_interval(λ, λ_rohn[1]), λ_hertz) + @test all(isnai, eigvals(A_large; alg = IntervalEigenContraction())) + @test all(isguaranteed, eigvals(A_large)) + + # Rohn bounds are rigorous + D = interval.([1.0 0.0; 0.0 1.0], 0.1; format = :midpoint) + @test all(λ -> in_interval(0.9, λ) && in_interval(1.1, λ), eigvals(D; alg = IntervalEigenRohn())) + + # NG flag + @test !any(isguaranteed, eigen(A_large .+ 1; alg = IntervalEigenRohn()).values) + @test !any(isguaranteed, eigvals(A_large .+ 1; alg = IntervalEigenRohn())) + @test !any(isguaranteed, eigen(LinearAlgebra.Symmetric(A_large .+ 1)).values) + @test all(isguaranteed, eigen(A_gen; alg = IntervalEigenRohn()).values) + + # `Symmetric` and `Hermitian` wrappers with non-thin diagonal + r_sym = eigen(LinearAlgebra.Symmetric(A_mid)) + @test in_interval(exact_eigs[1], real(r_sym.values[1])) + @test in_interval(exact_eigs[2], real(r_sym.values[2])) + r_herm = eigen(LinearAlgebra.Hermitian(complex.(A_mid))) + @test in_interval(exact_eigs[1], real(r_herm.values[1])) + @test in_interval(exact_eigs[2], real(r_herm.values[2])) + @test in_interval(exact_eigs[1], eigvals(LinearAlgebra.Symmetric(A_mid))[1]) + @test in_interval(exact_eigs[1], real(eigvals(LinearAlgebra.Hermitian(complex.(A_mid)))[1])) + + # repeated eigenvalues fall back to Rohn + r_rep = eigen(interval.([1.0 0.0; 0.0 1.0])) + @test all(λ -> isequal_interval(λ, interval(1)), r_rep.values) + @test all(isnai, eigen(interval.([1.0 0.0; 0.0 1.0]); alg = IntervalEigenContraction()).values) + + # unbounded entries + U = [interval(1) interval(0, Inf); interval(0) interval(1)] + for alg ∈ (IntervalEigen(), IntervalEigenContraction(), IntervalEigenRohn(), IntervalEigenHertz()) + @test all(isnai, eigen(U; alg).values) + @test all(isnai, eigvals(U; alg)) + end + + # LinearAlgebra matrix functions of `Symmetric` and `Hermitian` matrices are not rigorous + S = interval.([2.0 1.0 0.0; 1.0 2.0 0.0; 0.0 0.0 -2.0^-50]) + @test_throws ArgumentError sqrt(S) + @test_throws ArgumentError sqrt(LinearAlgebra.Symmetric(S)) + @test_throws ArgumentError LinearAlgebra.Symmetric(S)^0.5 + @test_throws ArgumentError LinearAlgebra.svd(LinearAlgebra.Symmetric(S)) + @test_throws ArgumentError LinearAlgebra.svdvals!(LinearAlgebra.Symmetric(S)) + @test all(isequal_interval.(LinearAlgebra.Symmetric(S)^2, S * S)) + @test all(isequal_interval.(LinearAlgebra.Symmetric(S)^2.0, S * S)) + @test all(isequal_interval.(exp(LinearAlgebra.Symmetric(A_small)), exp(A_small))) + @test all(isequal_interval.(log(LinearAlgebra.Symmetric(A_small)), log(A_small))) end @testset "inv" begin From 4d76209940593eae134867f0621716bc6411419e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 4 Nov 2026 11:25:47 +0100 Subject: [PATCH 10/10] Address second round of review comments - Intersect Gershgorin with Rohn in the default `eigvals`, as neither is always tighter, and use it as the fallback of the default `eigen` - Fix Gershgorin: the discs are complex even for real matrices, and overlapping discs are merged into the hull of their connected component - Only use the symmetric reduction of Rohn and Hertz for `Symmetric` and `Hermitian` wrappers, as a plain matrix with symmetric entries contains non-symmetric matrices - Return real eigenvalues for real symmetric and Hermitian wrappers - Fix the similarity transform for complex Hermitian midpoint matrices - Merge the real and complex `eigvals!` methods - Remove the `svdvals!` override, which is rigorous since it relies on `eigvals` - Move the linear algebra docs to a new Interfaces page --- docs/make.jl | 3 +- docs/src/interfaces/linearalgebra.md | 61 ++++++ docs/src/manual/usage.md | 22 --- ext/IntervalArithmeticLinearAlgebraExt.jl | 182 ++++++++++++------ src/IntervalArithmetic.jl | 20 +- .../ext/IntervalArithmeticLinearAlgebraExt.jl | 42 +++- 6 files changed, 234 insertions(+), 96 deletions(-) create mode 100644 docs/src/interfaces/linearalgebra.md diff --git a/docs/make.jl b/docs/make.jl index d4d350ea..809874e8 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -30,7 +30,8 @@ makedocs(; "API" => "manual/api.md" ], "Interfaces" => [ - "Arblib.jl" => "interfaces/arblib.md" + "Arblib.jl" => "interfaces/arblib.md", + "LinearAlgebra" => "interfaces/linearalgebra.md" ], "Examples" => ["Rigorous computation of ``\\pi``" => "examples/pi.md"], ], diff --git a/docs/src/interfaces/linearalgebra.md b/docs/src/interfaces/linearalgebra.md new file mode 100644 index 00000000..db97945c --- /dev/null +++ b/docs/src/interfaces/linearalgebra.md @@ -0,0 +1,61 @@ +# LinearAlgebra + +Interval matrices support standard linear algebra operations from the `LinearAlgebra` +standard library. + +```@repl linearalgebra +using IntervalArithmetic, LinearAlgebra +A = interval.([1.0 0.5; 0.5 2.0], 0.1; format = :midpoint) +``` + + + +## Eigenvalues and eigenvectors + +The eigenvalues (`eigvals`) and the eigendecomposition (`eigen`) are computed by one of +several algorithms, all subtypes of +[`AbstractIntervalEigenAlg`](@ref IntervalArithmetic.AbstractIntervalEigenAlg), which +can be selected via the `alg` keyword argument: + +| Algorithm | `eigvals` | `eigen` | +| :-------- | :-------- | :------ | +| [`IntervalEigen`](@ref IntervalArithmetic.IntervalEigen) (default) | Gershgorin circle theorem intersected with Rohn | contraction mapping theorem, falls back to `eigvals` and NaI eigenvectors | +| [`IntervalEigenContraction`](@ref IntervalArithmetic.IntervalEigenContraction) | contraction mapping theorem | contraction mapping theorem | +| [`IntervalEigenRohn`](@ref IntervalArithmetic.IntervalEigenRohn) | Rohn enclosure | Rohn enclosure, NaI eigenvectors | +| [`IntervalEigenHertz`](@ref IntervalArithmetic.IntervalEigenHertz) | Hertz exact hull (exponential complexity) | Hertz exact hull, NaI eigenvectors | + +The Gershgorin circle theorem and the contraction mapping theorem enclose each +eigenvalue individually, but only succeed when the interval radii are small enough. +The Rohn and Hertz methods return a single interval enclosing all the eigenvalues, and +are valid for any interval radius. As neither the Gershgorin circle theorem nor the Rohn +enclosure is always tighter than the other, the default takes their intersection. + +```@repl linearalgebra +eigvals(A) +eigen(A) +eigvals(Symmetric(A); alg = IntervalEigenRohn()) +eigvals(Symmetric(A); alg = IntervalEigenHertz()) +``` + +When the algorithm fails, the eigenvalues and eigenvectors are NaIs. + +```@repl linearalgebra +eigen(interval.([1.0 0.5; 0.5 2.0], 1.0; format = :midpoint); alg = IntervalEigenContraction()) +``` + +A matrix of intervals contains non-symmetric matrices, even if its entries are +symmetric. Use the `Symmetric` and `Hermitian` wrappers to only consider the symmetric +and Hermitian matrices: their eigenvalues are real, and the Rohn and Hertz enclosures +are tighter. + +```@repl linearalgebra +eigvals(A; alg = IntervalEigenRohn()) +eigvals(Symmetric(A)) +``` + +!!! note + The matrix functions of `LinearAlgebra` for `Symmetric` and `Hermitian` matrices (e.g. + `sqrt`, `cos`, `svd`) assume that the eigenvectors are orthonormal and that the + eigenvalues are accurate up to roundoff errors. They throw an `ArgumentError` for + interval matrices, except `exp`, `log` and integer powers, which have rigorous + implementations. diff --git a/docs/src/manual/usage.md b/docs/src/manual/usage.md index b228428d..d3de3e68 100644 --- a/docs/src/manual/usage.md +++ b/docs/src/manual/usage.md @@ -181,28 +181,6 @@ and deal properly with the decorations of intervals. -### Linear algebra - -Interval matrices support standard linear algebra operations from the `LinearAlgebra` -standard library. The eigendecomposition (`eigen`, `eigvals`) uses several algorithms -that can be selected via the `alg` keyword argument: - -- [`IntervalEigen`](@ref IntervalArithmetic.IntervalEigen) (default): tries the contraction mapping theorem first, falls back to the Rohn enclosure. -- [`IntervalEigenContraction`](@ref IntervalArithmetic.IntervalEigenContraction): contraction mapping only. -- [`IntervalEigenRohn`](@ref IntervalArithmetic.IntervalEigenRohn): Rohn eigenvalue enclosure. -- [`IntervalEigenHertz`](@ref IntervalArithmetic.IntervalEigenHertz): Hertz exact hull (exponential complexity). - -All algorithms are subtypes of [`AbstractIntervalEigenAlg`](@ref IntervalArithmetic.AbstractIntervalEigenAlg). - -```@repl usage -using LinearAlgebra -A = interval.([1.0 0.5; 0.5 2.0]) -eigen(A) -eigen(A; alg = IntervalEigenRohn()) -``` - - - ### Custom interval bounds type A `BareInterval{T}` or `Interval{T}` have the restriction `T <: Union{Rational,AbstractFloat}` which is the parametric type for the bounds of the interval. Supposing one wishes to use their own numeric type `MyNumType <: Union{Rational,AbstractFloat}`, they must provide their own arithmetic operations (with correct rounding!). diff --git a/ext/IntervalArithmeticLinearAlgebraExt.jl b/ext/IntervalArithmeticLinearAlgebraExt.jl index b313e5bb..2dfd39b1 100644 --- a/ext/IntervalArithmeticLinearAlgebraExt.jl +++ b/ext/IntervalArithmeticLinearAlgebraExt.jl @@ -57,14 +57,7 @@ end const _EigenAlg = IntervalArithmetic.AbstractIntervalEigenAlg -function LinearAlgebra.eigvals!(A::AbstractMatrix{<:Interval}; alg::_EigenAlg=IntervalArithmetic.IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) - # note: this function does not overwrite `A` - λ = _eigvals(alg, A, permute, scale, sortby) - _ensure_ng_flag!(λ, all(isguaranteed, A)) - return λ -end - -function LinearAlgebra.eigvals!(A::AbstractMatrix{<:Complex{<:Interval}}; alg::_EigenAlg=IntervalArithmetic.IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) +function LinearAlgebra.eigvals!(A::AbstractMatrix{<:RealOrComplexI}; alg::_EigenAlg=IntervalArithmetic.IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) # note: this function does not overwrite `A` λ = _eigvals(alg, A, permute, scale, sortby) _ensure_ng_flag!(λ, all(isguaranteed, A)) @@ -79,19 +72,50 @@ const _SymHermInterval = Union{ LinearAlgebra.Hermitian{<:Interval}, LinearAlgebra.Hermitian{<:Complex{<:Interval}}} +# the eigenvalues of Hermitian matrices are real, and so are the eigenvectors of real +# symmetric matrices +const _RealSymInterval = Union{LinearAlgebra.Symmetric{<:Interval},LinearAlgebra.Hermitian{<:Interval}} + +_eigenvalues(::Union{_RealSymInterval,LinearAlgebra.Hermitian{<:Complex{<:Interval}}}, λ) = real.(λ) +_eigenvalues(::_SymHermInterval, λ) = λ + +_eigenvectors(::_RealSymInterval, v) = real.(v) +_eigenvectors(::_SymHermInterval, v) = v + LinearAlgebra.eigvals(A::_SymHermInterval; alg::_EigenAlg=IntervalArithmetic.IntervalEigen(), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = - LinearAlgebra.eigvals!(Matrix(A); alg, sortby) + LinearAlgebra.eigvals!(A; alg, sortby) + +function LinearAlgebra.eigvals!(A::_SymHermInterval; alg::_EigenAlg=IntervalArithmetic.IntervalEigen(), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) + # note: this function does not overwrite `A` + λ = _eigenvalues(A, _eigvals(alg, A, true, true, sortby)) + _ensure_ng_flag!(λ, all(isguaranteed, A)) + return λ +end -LinearAlgebra.eigvals!(A::_SymHermInterval; alg::_EigenAlg=IntervalArithmetic.IntervalEigen(), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = - LinearAlgebra.eigvals!(Matrix(A); alg, sortby) +# the Gershgorin circle theorem and the contraction mapping theorem are applied to the +# dense matrix, as the LinearAlgebra methods for the wrappers may throw for intervals +_dense(A) = A +_dense(A::_SymHermInterval) = Matrix(A) -# IntervalEigen (default): Gershgorin circle theorem, falls back to Rohn on failure +# IntervalEigen (default): Gershgorin circle theorem intersected with Rohn, as neither of +# the two is always tighter than the other function _eigvals(::IntervalArithmetic.IntervalEigen, A, permute, scale, sortby) - λ = _eigvals_gershgorin(A, permute, scale, sortby) - all(isbounded, λ) && return λ - return _eigvals(IntervalArithmetic.IntervalEigenRohn(), A, permute, scale, sortby) + box = _eigvals(IntervalArithmetic.IntervalEigenRohn(), A, permute, scale, sortby) + λ = _eigvals_gershgorin(_dense(A), permute, scale, sortby) + all(isbounded, λ) || return box + return _intersect_eigbox.(λ, box) end +_intersect_eigbox(λ::Interval, box::Interval) = _intersect_decorated(λ, box) +_intersect_eigbox(λ::Interval, box::Complex{<:Interval}) = _intersect_decorated(λ, real(box)) +# `box` is real only if the matrix is real symmetric or Hermitian +_intersect_eigbox(λ::Complex{<:Interval}, box::Interval) = complex(_intersect_decorated(real(λ), box)) +_intersect_eigbox(λ::Complex{<:Interval}, box::Complex{<:Interval}) = + complex(_intersect_decorated(real(λ), real(box)), _intersect_decorated(imag(λ), imag(box))) + +_intersect_decorated(x::Interval, y::Interval) = + IntervalArithmetic.setdecoration(intersect_interval(x, y), min(decoration(x), decoration(y))) + _eigvals(alg::IntervalArithmetic.IntervalEigenContraction, A, permute, scale, sortby) = first(_eigen(alg, A, permute, scale, sortby)) @@ -103,25 +127,26 @@ end function _eigvals_gershgorin(A::AbstractMatrix{<:Interval}, permute, scale, sortby) _isbounded_matrix(A) || return fill(nai(eltype(A)), LinearAlgebra.checksquare(A)) - λ = _gershgorin(A, permute, scale, sortby) - isreal(λ) && return real(λ) - _fold_conjugate!(λ) - isreal(λ) && return real(λ) + λ, component = _gershgorin(A, permute, scale, sortby) + _fold_conjugate!(λ, component) + all(isthinzero ∘ imag, λ) && return real(λ) return λ end function _eigvals_gershgorin(A::AbstractMatrix{<:Complex{<:Interval}}, permute, scale, sortby) _isbounded_matrix(A) || return fill(nai(eltype(A)), LinearAlgebra.checksquare(A)) - return _gershgorin(A, permute, scale, sortby) + return first(_gershgorin(A, permute, scale, sortby)) end _isbounded_matrix(A) = all(x -> isbounded(real(x)) & isbounded(imag(x)), A) function _gershgorin(A, permute, scale, sortby) - # Gershgorin circle theorem + # Gershgorin circle theorem: a connected component of `k` discs contains exactly `k` + # eigenvalues, so each disc is replaced by the hull of its connected component B = _similarity_transform(A, permute, scale, sortby) - v = LinearAlgebra.diag(B) - T = eltype(v) + c = LinearAlgebra.diag(B) + T = real(eltype(c)) + v = Vector{Complex{T}}(undef, length(c)) for j ∈ axes(B, 1) r = zero(T) for i ∈ axes(B, 2) @@ -129,32 +154,50 @@ function _gershgorin(A, permute, scale, sortby) r += abs(B[i,j]) end end - v[j] = interval(v[j], r; format = :midpoint) + # the disc is in the complex plane even if the matrix is real + v[j] = complex(interval(real(c[j]), r; format = :midpoint), interval(imag(c[j]), r; format = :midpoint)) + end + component = collect(eachindex(v)) + for i ∈ eachindex(v), j ∈ eachindex(v) + if (component[i] ≠ component[j]) & !isdisjoint_interval(v[i], v[j]) + replace!(component, component[j] => component[i]) + end end - return v + for k ∈ unique(component) + idxs = findall(==(k), component) + v[idxs] .= Ref(reduce(_hull_decorated, view(v, idxs))) + end + return v, component end +_hull_decorated(x::Interval, y::Interval) = + IntervalArithmetic.setdecoration(hull(x, y), min(decoration(x), decoration(y))) +_hull_decorated(x::Complex{<:Interval}, y::Complex{<:Interval}) = + complex(_hull_decorated(real(x), real(y)), _hull_decorated(imag(x), imag(y))) + function _similarity_transform(A, permute, scale, sortby) mA = mid.(A) mλ, mV = LinearAlgebra.eigen(mA; permute = permute, scale = scale, sortby = sortby) - mλ .+= LinearAlgebra.diag(mV \ (mA * mV - mV * LinearAlgebra.Diagonal(mλ))) + # note: `mλ` is real for Hermitian `mA`, while the correction may be complex + mλ = mλ + LinearAlgebra.diag(mV \ (mA * mV - mV * LinearAlgebra.Diagonal(mλ))) Λ = LinearAlgebra.Diagonal(interval(mλ)) V = interval(mV) V .= Λ .+ inv(V) * (A * V - V * Λ) return V end -function _fold_conjugate!(λ) +# the eigenvalues of a real matrix come in conjugate pairs: the conjugate of an eigenvalue +# in the component of `λ[i]` lies in one of the components intersecting `conj(λ[i])` +function _fold_conjugate!(λ, component) + μ = copy(λ) for i ∈ eachindex(λ) - λᵢ = λ[i] - idxs = findall(j -> (j ≠ i) & !isdisjoint_interval(conj(λᵢ), λ[j]), eachindex(λ)) - if isempty(idxs) - λ[i] = real(λᵢ) - else - w = view(λ, idxs) - z = conj(intersect_interval(conj(λᵢ), reduce(intersect_interval, w))) - z = complex(IntervalArithmetic.setdecoration(real(z), min(decoration(real(λᵢ)), minimum(decoration ∘ real, w))), IntervalArithmetic.setdecoration(imag(z), min(decoration(imag(λᵢ)), minimum(decoration ∘ imag, w)))) - λ[i] = z + idxs = findall(j -> !isdisjoint_interval(conj(μ[i]), μ[j]), eachindex(μ)) + if idxs == [i] + # the only eigenvalue of the component is its own conjugate + λ[i] = complex(real(μ[i])) + elseif !isempty(idxs) + z = conj(reduce(_hull_decorated, view(μ, idxs))) + λ[i] = complex(_intersect_decorated(real(μ[i]), real(z)), _intersect_decorated(imag(μ[i]), imag(z))) end end return λ @@ -171,10 +214,16 @@ LinearAlgebra.eigen(A::AbstractMatrix{<:RealOrComplexI}; alg::_EigenAlg=Interval LinearAlgebra.eigen!(A; alg, permute, scale, sortby) LinearAlgebra.eigen(A::_SymHermInterval; alg::_EigenAlg=IntervalArithmetic.IntervalEigen(), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = - LinearAlgebra.eigen!(Matrix(A); alg, sortby) + LinearAlgebra.eigen!(A; alg, sortby) -LinearAlgebra.eigen!(A::_SymHermInterval; alg::_EigenAlg=IntervalArithmetic.IntervalEigen(), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = - LinearAlgebra.eigen!(Matrix(A); alg, sortby) +function LinearAlgebra.eigen!(A::_SymHermInterval; alg::_EigenAlg=IntervalArithmetic.IntervalEigen(), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) + # note: this function does not overwrite `A` + true_λ, true_v = _eigen(alg, A, true, true, sortby) + ng_flag = all(isguaranteed, A) + true_λ = _ensure_ng_flag!(_eigenvalues(A, true_λ), ng_flag) + true_v = _ensure_ng_flag!(_eigenvectors(A, true_v), ng_flag) + return LinearAlgebra.Eigen(true_λ, true_v) +end function LinearAlgebra.eigen!(A::AbstractMatrix{<:RealOrComplexI}; alg::_EigenAlg=IntervalArithmetic.IntervalEigen(), permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) # note: this function does not overwrite `A` @@ -185,18 +234,20 @@ function LinearAlgebra.eigen!(A::AbstractMatrix{<:RealOrComplexI}; alg::_EigenAl return LinearAlgebra.Eigen(true_λ, true_v) end -# IntervalEigen (default): contraction mapping theorem, falls back to Rohn on failure -function _eigen(::IntervalArithmetic.IntervalEigen, A, permute, scale, sortby) +# IntervalEigen (default): contraction mapping theorem, falls back to the eigenvalues of +# `IntervalEigen` and NaI eigenvectors on failure +function _eigen(alg::IntervalArithmetic.IntervalEigen, A, permute, scale, sortby) _isbounded_matrix(A) || return _eigen_nai(A) - contraction_ok, true_λ, true_v = _eigen_contraction(A, permute, scale, sortby) + contraction_ok, true_λ, true_v = _eigen_contraction(_dense(A), permute, scale, sortby) contraction_ok && return true_λ, true_v - return _eigen(IntervalArithmetic.IntervalEigenRohn(), A, permute, scale, sortby) + λ = _eigvals(alg, A, permute, scale, sortby) + return λ, fill(nai(complex(eltype(A))), length(λ), length(λ)) end # IntervalEigenContraction: contraction mapping theorem only (NaI if it fails) function _eigen(::IntervalArithmetic.IntervalEigenContraction, A, permute, scale, sortby) _isbounded_matrix(A) || return _eigen_nai(A) - contraction_ok, true_λ, true_v = _eigen_contraction(A, permute, scale, sortby) + contraction_ok, true_λ, true_v = _eigen_contraction(_dense(A), permute, scale, sortby) contraction_ok && return true_λ, true_v return _eigen_nai(A) end @@ -293,25 +344,24 @@ function Base.iterate(O::_Orthants, state=1) end # rigorous bounds on the extreme eigenvalues of a symmetric matrix of floats -# note: call the Gershgorin enclosure directly to avoid the fallback to Rohn - -_eigmax(A::AbstractMatrix{<:AbstractFloat}) = - maximum(x -> sup(real(x)), _eigvals_gershgorin(interval(A), true, true, LinearAlgebra.eigsortby)) +# note: call the Gershgorin enclosure directly to avoid the intersection with Rohn -_eigmin(A::AbstractMatrix{<:AbstractFloat}) = - minimum(x -> inf(real(x)), _eigvals_gershgorin(interval(A), true, true, LinearAlgebra.eigsortby)) +function _eigrange(A::AbstractMatrix{<:AbstractFloat}) + λ = _eigvals_gershgorin(interval(A), true, true, LinearAlgebra.eigsortby) + return minimum(x -> inf(real(x)), λ), maximum(x -> sup(real(x)), λ) +end # Rohn: fast eigenvalue enclosure for symmetric interval matrices -function _eigbox(A::LinearAlgebra.Symmetric{Interval{T},Matrix{Interval{T}}}, ::IntervalArithmetic.IntervalEigenRohn) where {T} - Ac = mid.(Matrix(A)) - ρ = _eigmax(radius.(Matrix(A))) - return interval(_eigmin(Ac), _eigmax(Ac)) + interval(-ρ, ρ) +function _eigbox(A::LinearAlgebra.Symmetric{Interval{T}}, ::IntervalArithmetic.IntervalEigenRohn) where {T} + λmin, λmax = _eigrange(mid.(Matrix(A))) + ρ = last(_eigrange(radius.(Matrix(A)))) + return interval(λmin, λmax) + interval(-ρ, ρ) end # Hertz: exact hull for symmetric interval matrices (exponential complexity) -function _eigbox(A::LinearAlgebra.Symmetric{Interval{T},Matrix{Interval{T}}}, ::IntervalArithmetic.IntervalEigenHertz) where {T} +function _eigbox(A::LinearAlgebra.Symmetric{Interval{T}}, ::IntervalArithmetic.IntervalEigenHertz) where {T} n = LinearAlgebra.checksquare(A) Amax = Matrix{T}(undef, n, n) Amin = Matrix{T}(undef, n, n) @@ -332,16 +382,28 @@ function _eigbox(A::LinearAlgebra.Symmetric{Interval{T},Matrix{Interval{T}}}, :: end end - λmin = min(λmin, _eigmin(Amin)) - λmax = max(λmax, _eigmax(Amax)) + λmin = min(λmin, first(_eigrange(Amin))) + λmax = max(λmax, last(_eigrange(Amax))) end return interval(λmin, λmax) end # eigenvalue enclosure for general and complex interval matrices (reduced to symmetric case) +# note: a plain matrix contains non-symmetric matrices even if its entries are symmetric, +# only the `Symmetric` and `Hermitian` wrappers restrict to the symmetric ones + +_eigbox(A::LinearAlgebra.Hermitian{<:Interval}, alg) = _eigbox(LinearAlgebra.Symmetric(Matrix(A)), alg) + +function _eigbox(M::LinearAlgebra.Hermitian{<:Complex{<:Interval}}, alg) + # the eigenvalues of `M` are those of the real symmetric matrix `[A B'; B A]` + A = real.(Matrix(M)) + B = imag.(Matrix(M)) + return _eigbox(LinearAlgebra.Symmetric([A B'; B A]), alg) +end + +_eigbox(A::LinearAlgebra.Symmetric{<:Complex{<:Interval}}, alg) = _eigbox(Matrix(A), alg) function _eigbox(A::AbstractMatrix{Interval{T}}, alg) where {T} - all(isequal_interval.(A, transpose(A))) && return _eigbox(LinearAlgebra.Symmetric(A), alg) half = interval(T, 1//2) λ = _eigbox(LinearAlgebra.Symmetric(half * (A + A')), alg) @@ -356,8 +418,6 @@ end function _eigbox(M::AbstractMatrix{Complex{Interval{T}}}, alg) where {T} A = real.(M) B = imag.(M) - # Hermitian: the eigenvalues of `M` are those of the real symmetric matrix `[A B'; B A]` - all(isequal_interval.(M, adjoint(M))) && return _eigbox(LinearAlgebra.Symmetric([A B'; B A]), alg) half = interval(T, 1//2) λ = _eigbox(LinearAlgebra.Symmetric(half * [A+A' B'-B; B-B' A+A']), alg) μ = _eigbox(LinearAlgebra.Symmetric(half * [B+B' A-A'; A'-A B+B']), alg) @@ -379,8 +439,8 @@ end Base.:^(A::_SymHermInterval, p::Integer) = Matrix(A)^p Base.:^(A::_SymHermInterval, p::Real) = isinteger(p) ? Matrix(A)^Integer(p) : _throw_matrix_function(^) +# note: `svdvals` is rigorous as it only relies on `eigvals` LinearAlgebra.svd(::_SymHermInterval; full::Bool=false) = _throw_matrix_function(LinearAlgebra.svd) -LinearAlgebra.svdvals!(::_SymHermInterval) = _throw_matrix_function(LinearAlgebra.svdvals!) _throw_matrix_function(f) = throw(ArgumentError("`$f` is not supported for `Symmetric` or `Hermitian` interval matrices")) diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl index b94733ad..4a5ecadd 100644 --- a/src/IntervalArithmetic.jl +++ b/src/IntervalArithmetic.jl @@ -137,17 +137,21 @@ abstract type AbstractIntervalEigenAlg end """ IntervalEigen -Default eigendecomposition algorithm for interval matrices. For `eigen`, tries the -contraction mapping theorem first (tight bounds when radii are small and eigenvalues are -simple). For `eigvals`, uses the Gershgorin circle theorem after a similarity transform. -In both cases, falls back to the Rohn eigenvalue enclosure if it fails. +Default eigendecomposition algorithm for interval matrices: +- `eigvals` uses the Gershgorin circle theorem, after a similarity transform by the + approximate eigenvectors of the midpoint matrix, intersected with + [`IntervalEigenRohn`](@ref). When the Gershgorin circle theorem fails (e.g. when + interval radii are too large), it returns the enclosure of [`IntervalEigenRohn`](@ref). +- `eigen` uses the contraction mapping theorem, as [`IntervalEigenContraction`](@ref). + When it fails, it returns the eigenvalues of `eigvals` and [`nai`](@ref) eigenvectors. """ struct IntervalEigen <: AbstractIntervalEigenAlg end """ IntervalEigenContraction -Eigendecomposition via the contraction mapping theorem only. Returns [`nai`](@ref) +Eigendecomposition via the contraction mapping theorem only, which gives tight bounds +when the interval radii are small and the eigenvalues are simple. Returns [`nai`](@ref) if the contraction fails (e.g. when interval radii are too large). """ struct IntervalEigenContraction <: AbstractIntervalEigenAlg end @@ -156,7 +160,8 @@ struct IntervalEigenContraction <: AbstractIntervalEigenAlg end IntervalEigenRohn Rohn eigenvalue enclosure (fast, O(n³)). Returns a single interval enclosing all -eigenvalues and [`nai`](@ref) eigenvectors. Valid for any interval radius. +eigenvalues and, when called with `eigen`, [`nai`](@ref) eigenvectors. Valid for any +interval radius. """ struct IntervalEigenRohn <: AbstractIntervalEigenAlg end @@ -164,7 +169,8 @@ struct IntervalEigenRohn <: AbstractIntervalEigenAlg end IntervalEigenHertz Hertz exact hull of eigenvalues (exponential O(2ⁿ) complexity). Returns a single -interval enclosing all eigenvalues and [`nai`](@ref) eigenvectors. +interval enclosing all eigenvalues and, when called with `eigen`, [`nai`](@ref) +eigenvectors. """ struct IntervalEigenHertz <: AbstractIntervalEigenAlg end diff --git a/test/ext/IntervalArithmeticLinearAlgebraExt.jl b/test/ext/IntervalArithmeticLinearAlgebraExt.jl index e439d6be..96adc032 100644 --- a/test/ext/IntervalArithmeticLinearAlgebraExt.jl +++ b/test/ext/IntervalArithmeticLinearAlgebraExt.jl @@ -142,8 +142,7 @@ end A_gen = interval.(M, 0.5; format = :midpoint) r_gen = LinearAlgebra.eigen(A_gen) for λ in mid_eigs - @test in_interval(real(λ), real(r_gen.values[1])) - @test in_interval(imag(λ), imag(r_gen.values[1])) + @test any(μ -> in_interval(real(λ), real(μ)) && in_interval(imag(λ), imag(μ)), r_gen.values) end # 3×3 symmetric with large radius @@ -179,10 +178,27 @@ end @test all(isnai, r_hertz.vectors) # `eigvals` uses the algorithm - λ_rohn = eigvals(A_large; alg = IntervalEigenRohn()) + λ_rohn = eigvals(LinearAlgebra.Symmetric(A_large); alg = IntervalEigenRohn()) @test all(λ -> issubset_interval(interval(-1.2, 4.2), λ), λ_rohn) - λ_hertz = eigvals(A_large; alg = IntervalEigenHertz()) + λ_hertz = eigvals(LinearAlgebra.Symmetric(A_large); alg = IntervalEigenHertz()) @test all(λ -> issubset_interval(λ, λ_rohn[1]), λ_hertz) + + # default intersects Gershgorin with Rohn, which is tighter than both + A_01 = interval.([1.0 0.5; 0.5 2.0], 0.1; format = :midpoint) + λ_default = eigvals(LinearAlgebra.Symmetric(A_01)) + λ_box = eigvals(LinearAlgebra.Symmetric(A_01); alg = IntervalEigenRohn()) + @test all(λ -> issubset_interval(λ, λ_box[1]), λ_default) + @test !isequal_interval(λ_default[1], λ_box[1]) + @test isdisjoint_interval(λ_default[1], λ_default[2]) + @test in_interval(exact_eigs[1], λ_default[1]) && in_interval(exact_eigs[2], λ_default[2]) + + # a plain matrix with symmetric entries also contains non-symmetric matrices + # e.g. `[0 1; -1 0]` with eigenvalues ±i + K = [interval(0) interval(-1, 1); interval(-1, 1) interval(0)] + for alg ∈ (IntervalEigen(), IntervalEigenRohn(), IntervalEigenHertz()) + @test all(λ -> in_interval(1, imag(λ)) && in_interval(-1, imag(λ)), eigvals(K; alg)) + @test all(λ -> isthinzero(imag(λ)), eigvals(LinearAlgebra.Symmetric(K); alg)) + end @test all(isnai, eigvals(A_large; alg = IntervalEigenContraction())) @test all(isguaranteed, eigvals(A_large)) @@ -206,6 +222,21 @@ end @test in_interval(exact_eigs[1], eigvals(LinearAlgebra.Symmetric(A_mid))[1]) @test in_interval(exact_eigs[1], real(eigvals(LinearAlgebra.Hermitian(complex.(A_mid)))[1])) + # eigenvalues of `Symmetric` and `Hermitian` matrices are real + C = interval.([1.0 1e-8; 1e-8 1.0]) + for alg ∈ (IntervalEigen(), IntervalEigenContraction(), IntervalEigenRohn(), IntervalEigenHertz()) + @test eigvals(LinearAlgebra.Symmetric(C); alg) isa Vector{Interval{Float64}} + @test eigvals(LinearAlgebra.Hermitian(C); alg) isa Vector{Interval{Float64}} + @test eigvals(LinearAlgebra.Hermitian(complex.(C)); alg) isa Vector{Interval{Float64}} + @test eigen(LinearAlgebra.Symmetric(C); alg).values isa Vector{Interval{Float64}} + @test eigen(LinearAlgebra.Symmetric(C); alg).vectors isa Matrix{Interval{Float64}} + @test eigen(LinearAlgebra.Hermitian(complex.(C)); alg).values isa Vector{Interval{Float64}} + @test eigen(LinearAlgebra.Hermitian(complex.(C)); alg).vectors isa Matrix{Complex{Interval{Float64}}} + end + # but not those of complex symmetric matrices + CS = interval.(complex.([1.0 1.0; 1.0 2.0], [1.0 0.0; 0.0 0.0])) + @test all(λ -> !in_interval(0, imag(λ)), eigvals(LinearAlgebra.Symmetric(CS))) + # repeated eigenvalues fall back to Rohn r_rep = eigen(interval.([1.0 0.0; 0.0 1.0])) @test all(λ -> isequal_interval(λ, interval(1)), r_rep.values) @@ -224,7 +255,8 @@ end @test_throws ArgumentError sqrt(LinearAlgebra.Symmetric(S)) @test_throws ArgumentError LinearAlgebra.Symmetric(S)^0.5 @test_throws ArgumentError LinearAlgebra.svd(LinearAlgebra.Symmetric(S)) - @test_throws ArgumentError LinearAlgebra.svdvals!(LinearAlgebra.Symmetric(S)) + σ = LinearAlgebra.svdvals(LinearAlgebra.Symmetric(S)) + @test in_interval(3, σ[1]) && in_interval(1, σ[2]) && in_interval(2.0^-50, σ[3]) @test all(isequal_interval.(LinearAlgebra.Symmetric(S)^2, S * S)) @test all(isequal_interval.(LinearAlgebra.Symmetric(S)^2.0, S * S)) @test all(isequal_interval.(exp(LinearAlgebra.Symmetric(A_small)), exp(A_small)))