diff --git a/src/bounds/ellipsoid.jl b/src/bounds/ellipsoid.jl index 6cfac229..3ff7d13f 100644 --- a/src/bounds/ellipsoid.jl +++ b/src/bounds/ellipsoid.jl @@ -38,7 +38,9 @@ volume(ell::Ellipsoid) = ell.volume # Returns the principal axes axes(ell::Ellipsoid) = ell.axes -function decompose(A::AbstractMatrix) +decompose(A::AbstractMatrix) = decompose(Symmetric(A)) # ensure that eigen() always returns real values + +function decompose(A::Symmetric) E = eigen(A) axlens = @. 1 / sqrt(E.values) axes = E.vectors * Diagonal(axlens) diff --git a/test/bounds/bounds.jl b/test/bounds/bounds.jl index 5d625a9f..d6bb2bdd 100644 --- a/test/bounds/bounds.jl +++ b/test/bounds/bounds.jl @@ -12,6 +12,13 @@ const BOUNDST = [ Bounds.MultiEllipsoid ] +@testset "pathological cases" begin + A_almost_symmetric = [2.6081830533175096e8 -5.4107420917559285e6 -1.9314298704966028e9 -2.360066561768968e9; -5.410742091755895e6 379882.440454782 6.715028007245775e7 2.0195280814040575e7; -1.931429870496611e9 6.715028007245693e7 9.811342987452753e10 -4.6579127705367036e7; -2.3600665617689605e9 2.0195280814042665e7 -4.6579127705418006e7 9.80946804720486e10] + # shouldn't fail: + ell = Bounds.Ellipsoid(zeros(4), A_almost_symmetric) + Bounds.volume(ell) +end + @testset "interface - $B, $T, D=$D" for B in BOUNDST, T in [Float32, Float64], D in 1:20 # creation, inspection bound = B(T, D) diff --git a/test/bounds/ellipsoids.jl b/test/bounds/ellipsoids.jl index b1a997b6..81b0b759 100644 --- a/test/bounds/ellipsoids.jl +++ b/test/bounds/ellipsoids.jl @@ -11,7 +11,8 @@ const NMAX = 20 @test volume(ell) ≈ volume_prefactor(N) * scale^N axs, axlens = decompose(ell) @test axlens ≈ fill(scale, N) - @test axs ≈ Bounds.axes(ell) ≈ scale * Matrix(I, N, N) + @test axs ≈ Bounds.axes(ell) + @test norm.(eachcol(axs)) == fill(scale, N) end @testset "Scaling" begin