rigorous computation of eigenvalues of matrices (#68)

* rigorous computation of eigenvalues of symmetric matrices

* specialized algorithm for simple eigenvalues

* added docstrings to docs

* fix merging conflicts

* updated algorithm to work with generic matrices

* fix eigenvectors for general matrices

docs/src/api/eigenvalues.md

@@ -1,10 +1,22 @@
+# Eigenvalues computations
+
 ```@index
 Pages = ["eigenvalues.md"]
 ```
 
+## Interval matrices eigenvalues
+
 ```@autodocs
 Modules=[IntervalLinearAlgebra]
 Pages=["interval_eigenvalues.jl"]
 Private=false
 ```
 
+## Floating point eigenvalues verification
+
+```@autodocs
+Modules=[IntervalLinearAlgebra]
+Pages=["verify_eigs.jl"]
+Private=false
+```
+

docs/src/references.md

@@ -211,6 +211,33 @@ S.M. Rump, [*Verification methods: Rigorous results using floating-point arithme
 ```
 ---
 
+#### [RUM01]
+
+```@raw html
+<ul><li>
+```
+Rump, Siegfried M. [*Computational error bounds for multiple or nearly multiple eigenvalues*](https://www.tuhh.de/ti3/paper/rump/Ru99c.pdf), Linear algebra and its applications 324.1-3 (2001): 209-226.
+```@raw html
+<li style="list-style: none"><details>
+<summary>bibtex</summary>
+```
+```
+@article{rump2001computational,
+  title={Computational error bounds for multiple or nearly multiple eigenvalues},
+  author={Rump, Siegfried M},
+  journal={Linear algebra and its applications},
+  volume={324},
+  number={1-3},
+  pages={209--226},
+  year={2001},
+  publisher={Elsevier}
+}
+```
+```@raw html
+</details></li></ul>
+```
+---
+
 #### [RUM99]
 
 ```@raw html

src/IntervalLinearAlgebra.jl

@@ -16,6 +16,8 @@ end
 
 @reexport using LinearAlgebra, IntervalArithmetic
 
+const IA = IntervalArithmetic
+
 export
     set_multiplication_mode, get_multiplication_mode,
     LinearKrawczyk, Jacobi, GaussSeidel, GaussianElimination, HansenBliekRohn, NonLinearOettliPrager, LinearOettliPrager,
@@ -24,7 +26,7 @@ export
     comparison_matrix, interval_norm, interval_isapprox, Orthants,
     is_H_matrix, is_strongly_regular, is_strictly_diagonally_dominant, is_Z_matrix, is_M_matrix,
     rref,
-    eigenbox
+    eigenbox, verify_eigen, bound_perron_frobenius_eigenvalue
 
 
 include("linear_systems/enclosures.jl")
@@ -38,4 +40,5 @@ include("classify.jl")
 include("rref.jl")
 
 include("eigenvalues/interval_eigenvalues.jl")
+include("eigenvalues/verify_eigs.jl")
 end

src/eigenvalues/verify_eigs.jl

@@ -0,0 +1,170 @@
+"""
+    verify_eigen(A[, λ, X0]; w=0.1, ϵ=1e-16, maxiter=10)
+
+Finds a rigorous bound for the eigenvalues and eigenvectors of `A`. Eigenvalues are treated
+as simple.
+
+### Input
+
+- `A`       -- matrix
+- `λ`       -- (optional) approximate value for an eigenvalue of `A`
+- `X0`      -- (optional) eigenvector associated to `λ`
+- `w`       -- relative inflation parameter
+- `ϵ`       -- absolute inflation parameter
+- `maxiter` -- maximum number of iterations
+
+### Output
+
+- Interval bounds on eigenvalues and eigenvectors.
+- A boolean certificate (or a vector of booleans if all eigenvalues are computed) `cert`.
+  If `cert[i]==true`, then the bounds for the ith eigenvalue and eigenvectore are rigorous,
+  otherwise not.
+
+### Algorithm
+
+The algorithm for this function is described in [[RUM01]](@ref).
+
+### Example
+
+```julia
+julia> A = Symmetric([1 2;2 3])
+2×2 Symmetric{Int64, Matrix{Int64}}:
+ 1  2
+ 2  3
+
+julia> evals, evecs, cert = verify_eigen(A);
+
+julia> evals
+2-element Vector{Interval{Float64}}:
+ [-0.236068, -0.236067]
+  [4.23606, 4.23607]
+
+julia> evecs
+2×2 Matrix{Interval{Float64}}:
+ [-0.850651, -0.85065]  [0.525731, 0.525732]
+  [0.525731, 0.525732]  [0.85065, 0.850651]
+
+julia> cert
+2-element Vector{Bool}:
+ 1
+ 1
+```
+"""
+function verify_eigen(A; kwargs...)
+    evals, evecs = eigen(mid.(A))
+
+    T = interval_eigtype(A, evals[1])
+    evalues = similar(evals, T)
+    evectors = similar(evecs, T)
+
+    cert = Vector{Bool}(undef, length(evals))
+    @inbounds for (i, λ₀) in enumerate(evals)
+        λ, v, flag = verify_eigen(A, λ₀, view(evecs, :,i); kwargs...)
+        evalues[i] = λ
+        evectors[:, i] .= v
+        cert[i] = flag
+
+    end
+    return evalues, evectors, cert
+end
+
+function verify_eigen(A, λ, X0; kwargs...)
+    ρ, X, cert = _verify_eigen(A, λ, X0; kwargs...)
+    return (real(λ) ± ρ) + (imag(λ) ± ρ) * im, X0 + X, cert
+end
+
+function verify_eigen(A::Symmetric, λ, X0; kwargs...)
+    ρ, X, cert = _verify_eigen(A, λ, X0; kwargs...)
+    return λ ± ρ, X0 + real.(X), cert
+end
+
+function _verify_eigen(A, λ::Number, X0::AbstractVector;
+                      w=0.1, ϵ=floatmin(), maxiter=10)
+
+    _, v = findmax(abs.(X0))
+
+    R = mid.(A) - λ * I
+    R[:, v] .= -X0
+    R = inv(R)
+    C = IA.Interval.(A) - λ * I
+    Z = -R * (C * X0)
+    C[:, v] .= -X0
+    C = I - R * C
+    Zinfl = w * IA.Interval.(-mag.(Z), mag.(Z)) .+ IA.Interval(-ϵ, ϵ)
+
+    X = Complex.(Z)
+    cert = false
+    @inbounds for _ in 1:maxiter
+        Y = (real.(X) + Zinfl) + (imag.(X) + Zinfl) * im
+
+        Ytmp = Y * Y[v]
+        Ytmp[v] = 0
+
+        X = Z + C * Y + R * Ytmp
+        cert = all(X .⊂ Y)
+        cert && break
+    end
+
+    ρ = mag(X[v])
+    X[v] = 0
+
+    return ρ, X, cert
+end
+
+
+"""
+    bound_perron_frobenius_eigenvalue(A, max_iter=10)
+
+Finds an upper bound for the Perron-Frobenius eigenvalue of the **non-negative** matrix `A`.
+
+### Input
+
+- `A` -- square real non-negative matrix
+- `max_iter` -- maximum number of iterations of the power method used internally to compute
+    an initial approximation of the Perron-Frobenius eigenvector
+
+### Example
+
+```julia-repl
+julia> A = [1 2;3 4]
+2×2 Matrix{Int64}:
+ 1  2
+ 3  4
+
+julia> bound_perron_frobenius_eigenvalue(A)
+5.372281323275249
+```
+"""
+function bound_perron_frobenius_eigenvalue(A::AbstractMatrix{T}, max_iter=10) where {T<:Real}
+    any(A .< 0) && throw(ArgumentError("Matrix contains negative entries"))
+    return _bound_perron_frobenius_eigenvalue(A, max_iter)
+end
+
+function _bound_perron_frobenius_eigenvalue(M, max_iter=10)
+
+    size(M, 1) == 1 && return M[1]
+    xpf = IA.Interval.(_power_iteration(M, max_iter))
+    Mxpf = M * xpf
+    ρ = zero(eltype(M))
+    @inbounds for (i, xi) in enumerate(xpf)
+        iszero(xi) && continue
+        tmp = Mxpf[i] / xi
+        ρ = max(ρ, tmp.hi)
+    end
+    return ρ
+end
+
+function _power_iteration(A, max_iter)
+    n = size(A,1)
+    xp = rand(n)
+    @inbounds for _ in 1:max_iter
+    xp .= A*xp
+    xp ./= norm(xp)
+    end
+    return xp
+end
+
+
+interval_eigtype(::Symmetric, ::T) where {T<:Real} = Interval{T}
+interval_eigtype(::AbstractMatrix, ::T) where {T<:Real} = Complex{Interval{T}}
+interval_eigtype(::AbstractMatrix, ::Complex{T}) where {T<:Real} = Complex{Interval{T}}

src/utils.jl

@@ -1,4 +1,5 @@
 mid(x::Real) = x
+mid(x::Complex) = x
 
 """
     interval_isapprox(a::Interval, b::Interval; kwargs)
@@ -172,3 +173,7 @@ end
 
 Base.firstindex(O::Orthants) = 1
 Base.lastindex(O::Orthants) = length(O)
+
+
+_unchecked_interval(x::Real) = Interval(x)
+_unchecked_interval(x::Complex) = Interval(real(x)) + Interval(imag(x)) * im

test/runtests.jl

@@ -9,8 +9,8 @@ include("test_utils.jl")
 
 
 include("test_eigenvalues/test_interval_eigenvalues.jl")
-
+include("test_eigenvalues/test_verify_eigs.jl")
 include("test_solvers/test_enclosures.jl")
 include("test_solvers/test_epsilon_inflation.jl")
-include("test_solvers/test_oettli_prager.jl")
 include("test_solvers/test_precondition.jl")
+include("test_solvers/test_oettli_prager.jl")

test/test_eigenvalues/test_verify_eigs.jl

@@ -0,0 +1,37 @@
+@testset "verify perron frobenius " begin
+    A = [1 2;3 4]
+    ρ = bound_perron_frobenius_eigenvalue(A)
+    @test  (5+sqrt(big(5)))/2 ≤ ρ
+end
+
+@testset "verified eigenvalues" begin
+    n = 5 # matrix size
+
+    # symmetric case
+    ev = sort(randn(n))
+    D = Diagonal(ev)
+    Q, _ = qr(rand(n, n))
+    A = Symmetric(IA.Interval.(Q) * D * IA.Interval.(Q)')
+
+    evals, evecs, cert = verify_eigen(A)
+    @test all(cert)
+    @test all(ev .∈ evals)
+
+
+    # real eigenvalues case
+    P = rand(n, n)
+    Pinv, _ = epsilon_inflation(P, Diagonal(ones(n)))
+    A = IA.Interval.(P) * D * Pinv
+
+    evals, evecs, cert = verify_eigen(A)
+    @test all(cert)
+    @test all(ev .∈ evals)
+
+    # test complex eigenvalues
+    ev = sort(rand(Complex{Float64}, n), by = x -> (real(x), imag(x)))
+    A = IA.Interval.(P) * Matrix(Diagonal(ev)) * Pinv
+
+    evals, evecs, cert = verify_eigen(A)
+    @test all(cert)
+    @test all(ev .∈ evals)
+end