diff --git a/Project.toml b/Project.toml index 02d875c8..bb5d28d4 100644 --- a/Project.toml +++ b/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" \ No newline at end of file +LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043" diff --git a/docs/src/api/eigenvalues.md b/docs/src/api/eigenvalues.md index 424dcadb..887a1250 100644 --- a/docs/src/api/eigenvalues.md +++ b/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 +``` \ No newline at end of file diff --git a/docs/src/references.md b/docs/src/references.md index 1415abab..69f9c5b6 100644 --- a/docs/src/references.md +++ b/docs/src/references.md @@ -104,6 +104,34 @@ L. Jaulin and B. Desrochers, [*Introduction to the algebra of separators with ap ``` --- + +#### [MIYA14] + +```@raw html + +``` +--- + + #### [NEU90] ```@raw html diff --git a/src/IntervalLinearAlgebra.jl b/src/IntervalLinearAlgebra.jl index a3b2a156..bcb1e5ba 100644 --- a/src/IntervalLinearAlgebra.jl +++ b/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 diff --git a/src/eigenvalues/miyajima_svd.jl b/src/eigenvalues/miyajima_svd.jl new file mode 100644 index 00000000..9b7aaa16 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/multiplication.jl b/src/multiplication.jl index b4342be6..637d9c9c 100644 --- a/src/multiplication.jl +++ b/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) diff --git a/src/numerical_test/multithread.jl b/src/numerical_test/multithread.jl new file mode 100644 index 00000000..5c01e024 --- /dev/null +++ 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 \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index de035dc7..2fba2a32 100644 --- a/test/runtests.jl +++ b/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") diff --git a/test/test_eigenvalues/test_miyajima_svd.jl b/test/test_eigenvalues/test_miyajima_svd.jl new file mode 100644 index 00000000..daf2535d --- /dev/null +++ b/test/test_eigenvalues/test_miyajima_svd.jl @@ -0,0 +1,20 @@ +@testset "verify singular value perron frobenius " begin + A = [0.5 0.5;0.3 0.7] + ρ = IntervalLinearAlgebra._bound_perron_frobenius_singularvalue(A) + @test 1 ≤ ρ +end + +@testset "verified svd" begin + A = [0.9..1.1 0 0 0 2; 0 0 3 0 0; 0 0 0 0 0; 0 2 0 0 0] + Σ = IntervalLinearAlgebra.svdbox(A, IntervalLinearAlgebra.M1()) + + @test all([3 ∈ Σ[1], 5 ∈ Σ[2]^2, 2 ∈ Σ[3], 0 ∈ Σ[4]]) + + Σ = IntervalLinearAlgebra.svdbox(A, IntervalLinearAlgebra.R1()) + @test all([3 ∈ Σ[1], 5 ∈ Σ[2]^2, 2 ∈ Σ[3], 0 ∈ Σ[4]]) + + A = im*A + + Σ = IntervalLinearAlgebra.svdbox(A, IntervalLinearAlgebra.R1()) + @test all([3 ∈ Σ[1], 5 ∈ Σ[2]^2, 2 ∈ Σ[3], 0 ∈ Σ[4]]) +end diff --git a/test/test_multiplication.jl b/test/test_multiplication.jl index d7c9b9ee..caea00ae 100644 --- a/test/test_multiplication.jl +++ b/test/test_multiplication.jl @@ -4,14 +4,22 @@ @test get_multiplication_mode() == Dict(:multiplication => :fast) A = [2..4 -2..1; -1..2 2..4] + imA = im*A + set_multiplication_mode(:slow) @test A * A == [0..18 -16..8; -8..16 0..18] + @test imA * imA == -1*[0..18 -16..8; -8..16 0..18] set_multiplication_mode(:rank1) @test A * A == [0..18 -16..8; -8..16 0..18] + @test imA * imA == -1*[0..18 -16..8; -8..16 0..18] set_multiplication_mode(:fast) @test A * A == [-2..19.5 -16..10; -10..16 -2..19.5] @test A * mid.(A) == [5..12.5 -8..2; -2..8 5..12.5] @test mid.(A) * A == [5..12.5 -8..2; -2..8 5..12.5] + + @test imA * imA == -1*[-2..19.5 -16..10; -10..16 -2..19.5] + @test mid.(A) * imA == im*[5..12.5 -8..2; -2..8 5..12.5] + @test imA * mid.(A) == im*[5..12.5 -8..2; -2..8 5..12.5] end diff --git a/test/test_numerical_test/test_numerical_test.jl b/test/test_numerical_test/test_numerical_test.jl new file mode 100644 index 00000000..d5e79107 --- /dev/null +++ b/test/test_numerical_test/test_numerical_test.jl @@ -0,0 +1,11 @@ +@testset "Test Numerical Test" begin + A = IntervalLinearAlgebra.NumericalTest.test_matrix(4) + @test A == [1.0 0 0 2^(-53); + 0 1.0 0 2^(-53); + 0 0 1.0 2^(-53); + 0 0 0 2^(-53)] + + # we test the singlethread version, so that CI points + # out if setting rounding modes is broken + @test IntervalLinearAlgebra.NumericalTest.rounding_test(1, 2) == true +end \ No newline at end of file