Repository navigation
Conversation
Benchmark Results (Julia v1)Time benchmarks
Memory benchmarks
|
|
@OlivierHnt kind reminder :) |
|
I have not forgotten but I do not know the math behind these algorithms so I need to find time to look into that. Btw I put your fix to disambiguate |
|
It's already the end of October 🤦 Ok I will put on an alarm on this. |
|
Ok so computing the extreme eigenvalues with the interval version of CI is failing pretty badly now though because the export list in Dealing with LinearAlgebra.jl is always a mess in my experience because there are so many special type matrices and they keep changing the internal logic and dispatch... Some comments
I also feel like there is an issue with |
- 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
|
CI is passing now and it should have addressed your comments |
| ### 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()) | ||
| ``` |
There was a problem hiding this comment.
Should this be in the Interfaces part of the documentation?
There was a problem hiding this comment.
Also something is a little misleading. Gershgorin and the contraction argument are two different things.
It seems like IntervalEigen for eigvals means Gershgorin, or did I misunderstand?
| function LinearAlgebra.eigvals!(A::AbstractMatrix{<:Interval}; 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) |
There was a problem hiding this comment.
you could use <:RealOrComplexI directly instead of also defining an other method below for complex interval, no?
| 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) |
There was a problem hiding this comment.
Rohn is always better than Gershgorin?
Should we take an intersection?
| 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) |
There was a problem hiding this comment.
but here we can guarantee that the eigenvalues are real actually
| IntervalEigenRohn | ||
|
|
||
| Rohn eigenvalue enclosure (fast, O(n³)). Returns a single interval enclosing all | ||
| eigenvalues and [`nai`](@ref) eigenvectors. Valid for any interval radius. |
| # 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")) |
There was a problem hiding this comment.
Are these really needed? I did not try it but are these not already throwing an error?
- 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
As discussed during the meeting, the current implementation in the LinearAlgebra extension is incorrect for large intervals while the implementation in IntervalLinearAlgebra works for large interval but is its computation of
eigvalsfor the matrix of centers and radius is just using floating point computation. By combining both, we have a fully rigorous eigen function for matrices of intervals.This was done by Claude Code but I have done a pass too, it looks good.
algCloses #753