Implemented Theorem 7 in Miya14

Project.toml

@@ -7,6 +7,7 @@ version = "0.1.6"
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OpenBLASConsistentFPCSR_jll = "6cdc7f73-28fd-5e50-80fb-958a8875b1af"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -23,4 +24,4 @@ julia = "1.6"
 
 [extras]
 IntervalConstraintProgramming = "138f1668-1576-5ad7-91b9-7425abbf3153"
-LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043"
+LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043"

docs/src/api/eigenvalues.md

@@ -20,3 +20,10 @@ Pages=["verify_eigs.jl"]
 Private=false
 ```
 
+## Singular Values Verification
+
+```@autodocs
+Modules=[IntervalLinearAlgebra]
+Pages=["miyajima_svd.jl"]
+Private=false
+```

docs/src/references.md

@@ -104,6 +104,34 @@ L. Jaulin and B. Desrochers, [*Introduction to the algebra of separators with ap
 ```
 ---
 
+
+#### [MIYA14]
+
+```@raw html
+<ul><li>
+```
+S. Miyajima, *Verified bounds for all the singular values of matrix*, CJapan J. Indust. Appl. Math. DOI 10.1007/s13160-014-0145-5 
+```@raw html
+<li style="list-style: none"><details>
+<summary>bibtex</summary>
+```
+```
+@article{miyajima2014svd,
+  title={Verified bounds for all the singular values of matrix},
+  author={Miyajima, Shinya},
+  journal={Japan Journal of Industrial and Applied Mathematics},
+  volume={31},
+  pages={513--539},
+  year={2014},
+  publisher={Springer}
+}
+```
+```@raw html
+</details></li></ul>
+```
+---
+
+
 #### [NEU90]
 
 ```@raw html

src/IntervalLinearAlgebra.jl

@@ -41,12 +41,49 @@ include("pils/pils_solvers.jl")
 
 include("eigenvalues/interval_eigenvalues.jl")
 include("eigenvalues/verify_eigs.jl")
+include("eigenvalues/miyajima_svd.jl")
+
+include("numerical_test/multithread.jl")
+
+using LinearAlgebra
+
+if Sys.ARCH == :x86_64
+    using OpenBLASConsistentFPCSR_jll
+else
+    @warn "The behaviour of multithreaded OpenBlas on this architecture is unclear,
+    we will fallback to single threaded OpenBLAS"
+end
 
 function  __init__()
     @require IntervalConstraintProgramming = "138f1668-1576-5ad7-91b9-7425abbf3153" include("linear_systems/oettli_nonlinear.jl")
     @require LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043" include("linear_systems/oettli_linear.jl")
+    if Sys.ARCH == :x86_64
+        @info "Switching to OpenBLAS with ConsistentFPCSR = 1 flag enabled, guarantees
+        correct floating point rounding mode over all threads."
+        BLAS.lbt_forward(OpenBLASConsistentFPCSR_jll.libopenblas_path; verbose =  true)
+
+        N = BLAS.get_num_threads()
+        K = 1024
+        if NumericalTest.rounding_test(N, K)
+            @info "OpenBLAS is giving correct rounding on a ($K,$K) test matrix on $N threads"
+        else
+            @warn "OpenBLAS is not rounding correctly on the test matrix"
+            @warn "The number of BLAS threads was set to 1 to ensure rounding mode is consistent"    
+            if !NumericalTest.rounding_test(1, K)
+                @warn "The rounding test failed on 1 thread"
+            end
+        end
+    else
+        BLAS.set_num_threads(1)
+        @warn "The number of BLAS threads was set to 1 to ensure rounding mode is consistent"
+        if !NumericalTest.rounding_test(1, 1024)
+            @warn "The rounding test failed on 1 thread"
+        end
+    end
 end
 
 set_multiplication_mode(config[:multiplication])
 
+using LinearAlgebra
+
 end

src/eigenvalues/miyajima_svd.jl

@@ -0,0 +1,180 @@
+# In this file we implement some of the algorithms from
+#
+# Shinya Miyajima
+#
+# Verified bounds for all the singular values of matrix
+# Japan J. Indust. Appl. Math.
+# DOI 10.1007/s13160-014-0145-5 
+
+abstract type AbstractIntervalSVDSolver end
+
+struct M1 <: AbstractIntervalSVDSolver end
+struct R1 <: AbstractIntervalSVDSolver end
+
+# The method is called svdbox to mirror eigenbox, even if it would be best to call it 
+# svd interval
+
+using LinearAlgebra
+
+function _power_iteration_singularvalue(A, max_iter)
+    n = size(A, 2)
+    xp = rand(n)
+    @inbounds for _ in 1:max_iter
+    xp .= A'*(A*xp)
+    xp ./= norm(xp)
+    end
+    return xp
+end
+
+# we use this function to bound the $2$ norm of a matrix,
+# as remarked in Rump, Verified bounds for singular values,
+# Equation 3.3
+function _bound_perron_frobenius_singularvalue(M, max_iter=10)
+
+    if any(M.<0)
+        @warn "The matrix has nonpositive entries; M'*M could be nonnegative, but careful"
+    end
+
+    size(M) == (1, 1) && return M[1]
+    xpf = IA.Interval.(_power_iteration_singularvalue(mid.(M), max_iter))
+    Mxpf = M'* (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
+
+export svdbox
+"""
+    svdbox(A[, method = R1()])
+
+Return an enclosure for all the singular values of `A`.
+
+# Input 
+-   `A` -- interval matrix
+- `method` -- the method used to compute the enclosure
+    Possible values are
+        - M1 Algorithm from Theorem 7 in [[MIYA14]](@ref)
+        - R1 Rump Algorithm from Theorem 5 in [[MIYA14]](@ref)
+
+### Algorithm 
+
+The algorithms used by the function are described in [[MIYA14]](@ref).
+
+### Example
+The matrix encloses a matrix that has singular values 
+``3, √5, 2, 0``.
+
+```julia
+julia> A = [0.9..1.1 0 0 0 2; 0 0 3 0 0; 0 0 0 0 0; 0 2 0 0 0]
+4×5 Matrix{Interval{Float64}}:
+     [0.899999, 1.10001]  [0, 0]  [0, 0]  [0, 0]  [2, 2]
+ [0, 0]                   [0, 0]  [3, 3]  [0, 0]  [0, 0]
+ [0, 0]                   [0, 0]  [0, 0]  [0, 0]  [0, 0]
+ [0, 0]                   [2, 2]  [0, 0]  [0, 0]  [0, 0]
+
+ julia> IntervalLinearAlgebra.svdbox(A, IntervalLinearAlgebra.M1())
+ 4-element Vector{Interval{Float64}}:
+  [2.89999, 3.10001]
+  [2.13606, 2.33607]
+  [1.89999, 2.10001]
+ [-0.100001, 0.100001]
+```
+
+"""
+svdbox(A) = svdbox(A, R1())
+
+function svdbox(A::AbstractMatrix{Interval{T}}, ::M1) where T
+    mA = mid.(A)
+    svdA  = svd(mA)
+    U = Interval{T}.(svdA.U)
+    Vt = Interval{T}.(svdA.Vt)
+    Σ = diagm(Interval{T}.(svdA.S))
+    V = Interval{T}.(svdA.V)
+
+    E = U*Σ*Vt - A
+    normE = sqrt(_bound_perron_frobenius_singularvalue(abs.(E)))
+
+    F = Vt*V-I
+    normF = sqrt(_bound_perron_frobenius_singularvalue(abs.(F)))
+
+    G = U'*U-I
+    normG = sqrt(_bound_perron_frobenius_singularvalue(abs.(G)))
+
+    @assert normF < 1 "It is not possible to verify the singular values with this precision"
+    @assert normG < 1 "It is not possible to verify the singular values with this precision"
+
+    svdbounds = [hull(σ*sqrt((1-normF)*(1-normG))-normE, 
+                     σ*sqrt((1+normF)*(1+normG))+normE)
+                for σ in diag(Σ)]
+
+    return svdbounds
+end
+
+function certifysvd(A::AbstractMatrix{Interval{T}}, U, V, Vt) where {T}
+    IU = Interval{T}.(U)
+    IV = Interval{T}.(V)
+    IVt = Interval{T}.(Vt)
+
+    F = IVt*IV-I
+    G = (IU)'*IU-I
+    normF = sqrt(_bound_perron_frobenius_singularvalue(abs.(F)))
+    normG = sqrt(_bound_perron_frobenius_singularvalue(abs.(G)))
+
+    @assert normF < 1 "It is not possible to verify the singular values with this precision"
+    @assert normG < 1 "It is not possible to verify the singular values with this precision"
+
+    M = IU'*A*IV
+    D = diag(M)
+    E = M-diagm(D)
+
+    normE = sqrt(_bound_perron_frobenius_singularvalue(abs.(E)))
+
+    svdbounds = [hull((abs(σ)-normE)/sqrt((1+normF)*(1+normG)), 
+                     (abs(σ)+normE)/sqrt((1-normF)*(1-normG)))
+                for σ in D]
+
+    return sort!(svdbounds, rev = true)
+end
+
+function certifysvd(A::AbstractMatrix{Complex{Interval{T}}}, U, V, Vt) where {T}
+    IU = Interval{T}.(real(U))+im*Interval{T}.(imag(U))
+    IV = Interval{T}.(real(V))+im*Interval{T}.(imag(V))
+    IVt = Interval{T}.(real(Vt))+im*Interval{T}.(imag(Vt))
+
+    F = IVt*IV-I
+    G = (IU)'*IU-I
+    normF = sqrt(_bound_perron_frobenius_singularvalue(abs.(F)))
+    normG = sqrt(_bound_perron_frobenius_singularvalue(abs.(G)))
+
+    @assert normF < 1 "It is not possible to verify the singular values with this precision"
+    @assert normG < 1 "It is not possible to verify the singular values with this precision"
+
+    M = IU'*A*IV
+    D = diag(M)
+    E = M-diagm(D)
+
+    normE = sqrt(_bound_perron_frobenius_singularvalue(abs.(E)))
+
+    svdbounds = [hull((abs(σ)-normE)/sqrt((1+normF)*(1+normG)), 
+                     (abs(σ)+normE)/sqrt((1-normF)*(1-normG)))
+                for σ in D]
+
+    return sort!(svdbounds, rev = true)
+end
+
+
+function svdbox(A::AbstractMatrix{Interval{T}}, ::R1) where T
+    mA = mid.(A)
+    svdA  = svd(mA)
+    return certifysvd(A, svdA.U, svdA.V, svdA.Vt)
+end
+
+function svdbox(A::AbstractMatrix{Complex{Interval{T}}}, ::R1) where T
+    mA = mid.(A)
+    svdA  = svd(mA)
+    return certifysvd(A, svdA.U, svdA.V, svdA.Vt)
+end

src/multiplication.jl

@@ -26,18 +26,45 @@ Sets the algorithm used to perform matrix multiplication with interval matrices.
 """
 function set_multiplication_mode(multype)
     type = MultiplicationType{multype}()
+
+    # real matrices
     @eval *(A::AbstractMatrix{Interval{T}}, B::AbstractMatrix{Interval{T}}) where T =
         *($type, A, B)
 
     @eval *(A::AbstractMatrix{T}, B::AbstractMatrix{Interval{T}}) where T = *($type, A, B)
-
+    
     @eval *(A::AbstractMatrix{Interval{T}}, B::AbstractMatrix{T}) where T = *($type, A, B)
-
+    
     @eval *(A::Diagonal, B::AbstractMatrix{Interval{T}}) where T = *($type, A, B)
+
     @eval *(A::AbstractMatrix{Interval{T}}, B::Diagonal) where T = *($type, A, B)
+
     config[:multiplication] = multype
 end
 
+function *(A::AbstractMatrix{Complex{Interval{T}}}, B::AbstractMatrix) where T
+    return real(A)*B+im*imag(A)*B
+end
+
+function *(A::AbstractMatrix, B::AbstractMatrix{Complex{Interval{T}}}) where T
+    return A*real(B)+im*A*imag(B)
+end
+
+function *(A::AbstractMatrix{Complex{Interval{T}}}, B::AbstractMatrix{Complex{Interval{T}}}) where T
+    rA, iA = real(A), imag(A)
+    rB, iB = real(B), imag(B)
+    return rA*rB-iA*iB+im*(iA*rB+rA*iB)
+end
+
+function *(A::AbstractMatrix{Complex{T}}, B::AbstractMatrix{Interval{T}}) where {T}
+    return real(A)*B+im*imag(A)*B
+end
+
+function *(A::AbstractMatrix{Interval{T}}, B::AbstractMatrix{Complex{T}}) where {T}
+    return A*real(B)+im*A*imag(B)
+end
+
+
 function *(::MultiplicationType{:slow}, A, B)
     TS = promote_type(eltype(A), eltype(B))
     return mul!(similar(B, TS, (size(A,1), size(B,2))), A, B)

src/numerical_test/multithread.jl

@@ -0,0 +1,38 @@
+module NumericalTest
+
+export rounding_test
+
+function test_matrix(k)
+    A = zeros(Float64, (k, k))
+    A[:, end] = fill(2^(-53), k)
+    for i in 1:k-1
+        A[i,i] = 1.0
+    end
+    return A
+end
+
+using LinearAlgebra
+"""
+    rounding_test(n, k)
+
+Let `u=fill(2^(-53), k-1)` and let A be the matrix
+[I u;
+0 2^(-53)]
+
+This test checks the result of A*A' in different rounding modes,
+running BLAS on `n` threads
+"""
+function rounding_test(n,k)
+
+
+    BLAS.set_num_threads( n )
+    A = test_matrix( k )
+    B = setrounding(Float64, RoundUp) do
+        BLAS.gemm('N', 'T', 1.0, A, A)
+    end
+
+    return all([B[i,i]==nextfloat(1.0) for i in 1:k-1])
+end
+
+
+end

test/runtests.jl

@@ -7,9 +7,10 @@ include("test_classify.jl")
 include("test_multiplication.jl")
 include("test_utils.jl")
 
-
+include("test_numerical_test/test_numerical_test.jl")
 include("test_eigenvalues/test_interval_eigenvalues.jl")
 include("test_eigenvalues/test_verify_eigs.jl")
+include("test_eigenvalues/test_miyajima_svd.jl")
 include("test_solvers/test_enclosures.jl")
 include("test_solvers/test_epsilon_inflation.jl")
 include("test_solvers/test_precondition.jl")