Add simple eigenvalue computations and determinant

src/matmul.jl

@@ -65,6 +65,73 @@ end
 
 
 
+# matrix eigenvalues
+
+function LinearAlgebra.eigvals!(A::AbstractMatrix{<:Interval}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby)
+    # note: this function does not overwrite `A`
+    v = _eigvals(A, permute, scale, sortby)
+    isreal(v) && return v
+    _fold_conjugate!(v)
+    isreal(v) && return real(v)
+    return v
+end
+
+LinearAlgebra.eigvals!(A::AbstractMatrix{<:Complex{<:Interval}}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
+    # note: this function does not overwrite `A`
+    _eigvals(A, permute, scale, sortby)
+
+function _eigvals(A, permute, scale, sortby)
+    # Gershgorin circle theorem
+    B = _similarity_transform(A, permute, scale, sortby)
+    v = LinearAlgebra.diag(B)
+    T = eltype(v)
+    for j ∈ axes(B, 1)
+        r = zero(T)
+        for i ∈ axes(B, 2)
+            if i ≠ j
+                r += abs(B[i,j])
+            end
+        end
+        v[j] = interval(v[j], r; format = :midpoint)
+    end
+    return v
+end
+
+function _similarity_transform(A, permute, scale, sortby)
+    mA = mid.(A)
+    mλ, mV = LinearAlgebra.eigen(mA; permute = permute, scale = scale, sortby = sortby)
+    mλ .+= LinearAlgebra.diag(mV \ (mA * mV - mV * LinearAlgebra.Diagonal(mλ)))
+    Λ = LinearAlgebra.Diagonal(interval(mλ))
+    V = interval(mV)
+    V .= Λ .+ inv(V) * (A * V - V * Λ)
+    return V
+end
+
+function _fold_conjugate!(v)
+    for i ∈ eachindex(v)
+        vᵢ = v[i]
+        idxs = findall(j -> (j ≠ i) & !isdisjoint_interval(conj(vᵢ), v[j]), eachindex(v))
+        if isempty(idxs)
+            v[i] = real(vᵢ)
+        else
+            w = view(v, idxs)
+            z = conj(intersect_interval(conj(vᵢ), reduce(intersect_interval, w)))
+            z = complex(setdecoration(real(z), min(decoration(real(vᵢ)), minimum(decoration ∘ real, w))), setdecoration(imag(z), min(decoration(imag(vᵢ)), minimum(decoration ∘ imag, w))))
+            v[i] = z
+        end
+    end
+    return v
+end
+
+
+
+# matrix determinant
+
+LinearAlgebra.det(A::AbstractMatrix{<:Interval}) = real(reduce(*, LinearAlgebra.eigvals(A)))
+LinearAlgebra.det(A::AbstractMatrix{<:Complex{<:Interval}}) = reduce(*, LinearAlgebra.eigvals(A))
+
+
+
 # matrix multiplication
 
 """