use Symmetric() to get rid of complex eigen() results

TuringLang · Nov 18, 2021 · fde53e4 · fde53e4
1 parent c202784
commit fde53e4
Show file tree

Hide file tree

Showing 2 changed files with 10 additions and 1 deletion.
diff --git 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
@@ -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)