Merge pull request #17 from JuliaAlgebra/bl/hermitian

Add support for hermitian matrix

src/MultivariateMoments.jl

@@ -20,6 +20,7 @@ include("moment.jl")
 include("measure.jl")
 include("expectation.jl")
 include("symmatrix.jl")
+include("hermitian_matrix.jl")
 include("moment_matrix.jl")
 include("atomic.jl")
 include("extract.jl")

src/hermitian_matrix.jl

@@ -0,0 +1,93 @@
+export VectorizedHermitianMatrix
+
+"""
+    struct VectorizedHermitianMatrix{T} <: AbstractMatrix{T}
+        Q::Vector{T}
+        n::Int
+    end
+
+Hermitian ``n \\times n`` matrix storing the vectorized upper triangular
+real part of the matrix followed by the vectorized upper triangular
+imaginary part in the `Q` vector (this corresponds to the vectorized
+form of `ComplexOptInterface.HermitianPositiveSemidefiniteConeTriangle`).
+It implement the `AbstractMatrix` interface except for `setindex!` as it might
+break its symmetry. The [`symmetric_setindex!`](@ref) function should be used
+instead.
+"""
+struct VectorizedHermitianMatrix{T} <: AbstractMatrix{T}
+    Q::Vector{T}
+    n::Int
+end
+
+Base.copy(Q::VectorizedHermitianMatrix) = VectorizedHermitianMatrix(copy(Q.Q), Q.n)
+function Base.map(f::Function, Q::VectorizedHermitianMatrix)
+    if Q.n <= 1
+        return VectorizedHermitianMatrix(map(f, Q.Q), Q.n)
+    else
+        el = f(Q[1, 2])
+        ret = VectorizedHermitianMatrix(similar(Q.Q, typeof(real(el))), Q.n)
+        for j in 1:Q.n
+            for i in 1:j
+                if i == 1 && j == 2
+                    x = el
+                else
+                    x = f(Q[i, j])
+                end
+                symmetric_setindex!(ret, x, i, j)
+            end
+        end
+        return ret
+    end
+end
+
+imag_map(n, i, j) = trimap(n, n) + trimap(i, j - 1)
+imag_map(Q::VectorizedHermitianMatrix, i, j) = imag_map(Q.n, i, j)
+
+function trimat(::Type{T}, f, n, σ) where {T}
+    Q = Vector{T}(undef, N + trimap(n - 1, n - 1))
+    for i in 1:n
+        for j in 1:i
+            x = f(σ[i], σ[j])
+            Q[trimap(i, j)] = real(x)
+            if i != j
+                Q[imag_map(n, i, j)] = imag(x)
+            end
+        end
+    end
+    return VectorizedHermitianMatrix{T}(Q, n)
+end
+
+Base.size(Q::VectorizedHermitianMatrix) = (Q.n, Q.n)
+
+"""
+    symmetric_setindex!(Q::VectorizedHermitianMatrix, value, i::Integer, j::Integer)
+
+Set `Q[i, j]` to the value `value` and `Q[j, i]` to the value `-value`.
+"""
+function symmetric_setindex!(Q::VectorizedHermitianMatrix, value, i::Integer, j::Integer)
+    if i > j
+        symmetric_setindex!(Q, conj(value), j, i)
+    else
+        Q.Q[trimap(i, j)] = real(value)
+        if i == j
+            if !iszero(imag(value))
+                error("Cannot set diagonal entry ($i, $j) of hermitian matrix with a non-zero imaginary part $(imag(value))")
+            end
+        else
+            Q.Q[imag_map(Q, i, j)] = imag(value)
+        end
+    end
+end
+
+function Base.getindex(Q::VectorizedHermitianMatrix, i::Integer, j::Integer)
+    I, J = max(i, j), min(i, j)
+    r = Q.Q[trimap(I, J)]
+    if i == j
+        return Complex(r)
+    else
+        c = Q.Q[imag_map(Q, I, J)]
+        return Complex(r, i < j ? c : -c)
+    end
+end
+Base.getindex(Q::VectorizedHermitianMatrix, I::Tuple) = Q[I...]
+Base.getindex(Q::VectorizedHermitianMatrix, I::CartesianIndex) = Q[I.I]

src/moment_matrix.jl

@@ -6,17 +6,17 @@ using SemialgebraicSets
 abstract type AbstractMomentMatrix{T, B<:MB.AbstractPolynomialBasis} <: AbstractMeasureLike{T} end
 
 """
-    mutable struct MomentMatrix{T, B<:MultivariateBases.AbstractPolynomialBasis} <: AbstractMeasureLike{T}
-        Q::SymMatrix{T}
+    mutable struct MomentMatrix{T, B<:MultivariateBases.AbstractPolynomialBasis, MT<:AbstractMatrix{T}} <: AbstractMeasureLike{T}
+        Q::MT
         basis::B
         support::Union{Nothing, AlgebraicSet}
     end
 
 Measure ``\\nu`` represented by the moments of the monomial matrix ``x x^\\top`` in the symmetric matrix `Q`.
 The set of points that are zeros of all the polynomials ``p`` such that ``\\mathbb{E}_{\\nu}[p] = 0`` is stored in the field `support` when it is computed.
 """
-mutable struct MomentMatrix{T, B<:MB.AbstractPolynomialBasis} <: AbstractMomentMatrix{T, B}
-    Q::SymMatrix{T}
+mutable struct MomentMatrix{T, B<:MB.AbstractPolynomialBasis, MT<:AbstractMatrix{T}} <: AbstractMomentMatrix{T, B}
+    Q::MT
     basis::B
     support::Union{Nothing, AlgebraicSet}
 end
@@ -63,11 +63,9 @@ function MomentMatrix(Q::AbstractMatrix, monos::AbstractVector)
 end
 moment_matrix(Q::AbstractMatrix, monos) = MomentMatrix(Q, monos)
 
-function getmat(μ::MomentMatrix)
-    Matrix(μ.Q)
-end
+getmat(μ::MomentMatrix) = Matrix(μ.Q)
 
-function measure(ν::MomentMatrix{T, <:MB.MonomialBasis}) where T
+function measure(ν::MomentMatrix{T, <:MB.MonomialBasis, SymMatrix{T}}) where T
     n = length(ν.basis)
     monos = ν.basis.monomials
     measure(ν.Q.Q, [monos[i] * monos[j] for i in 1:n for j in 1:i])

src/symmatrix.jl

@@ -1,3 +1,5 @@
+export SymMatrix
+
 """
     struct SymMatrix{T} <: AbstractMatrix{T}
         Q::Vector{T}
@@ -26,7 +28,7 @@ Base.map(f::Function, Q::SymMatrix) = SymMatrix(map(f, Q.Q), Q.n)
 
 # j <= i
 function trimap(i, j)
-    div((i-1)*i, 2) + j
+    div((i - 1) * i, 2) + j
 end
 
 function trimat(::Type{T}, f, n, σ) where {T}

test/hermitian_matrix.jl

@@ -0,0 +1,28 @@
+using Test
+using MultivariateMoments
+
+@testset "VectorizedHermitianMatrix" begin
+    Q = VectorizedHermitianMatrix([1, 2, 3, -1], 2)
+    @test 1 == @inferred Q[1, 1]
+    @test 2 - im == @inferred Q[1, 2]
+    @test 2 + im == @inferred Q[2, 1]
+    @test 3 == @inferred Q[2, 2]
+    symmetric_setindex!(Q, 4, 1, 1)
+    @test Q.Q == [4, 2, 3, -1]
+    symmetric_setindex!(Q, 5 - 2im, 1, 2)
+    @test Q.Q == [4, 5, 3, -2]
+    P = copy(Q)
+    @test P.n == 2
+    symmetric_setindex!(P, 6, 2, 2)
+    @test Q.Q == [4, 5, 3, -2]
+    @test P.n == 2
+    @test P.Q == [4, 5, 6, -2]
+    R = map(i -> i - 1, P)
+    @test R.n == 2
+    @test R.Q == [3, 4, 5, -2]
+    @test P.Q == [4, 5, 6, -2]
+    S = map(conj, R)
+    @test S.n == 2
+    @test S.Q == [3, 4, 5, 2]
+    @test_throws ErrorException symmetric_setindex!(S, im, 1, 1)
+end

test/moment_matrix.jl

@@ -1,20 +1,3 @@
-@testset "SymMatrix" begin
-    Q = SymMatrix([1, 2, 3], 2)
-    symmetric_setindex!(Q, 4, 1, 1)
-    @test Q.Q == [4, 2, 3]
-    symmetric_setindex!(Q, 5, 1, 2)
-    @test Q.Q == [4, 5, 3]
-    P = copy(Q)
-    @test P.n == 2
-    symmetric_setindex!(P, 6, 2, 2)
-    @test Q.Q == [4, 5, 3]
-    @test P.n == 2
-    @test P.Q == [4, 5, 6]
-    R = map(i -> i - 1, P)
-    @test R.n == 2
-    @test R.Q == [3, 4, 5]
-end
-
 @testset "MomentMatrix" begin
     Mod.@polyvar x y
     @test_throws ArgumentError moment_matrix(measure([1], [x]), [y])

test/runtests.jl

@@ -7,6 +7,9 @@ using SemialgebraicSets
 using Test
 using LinearAlgebra
 
+include("symmatrix.jl")
+include("hermitian_matrix.jl")
+
 # Taken from JuMP/test/solvers.jl
 function try_import(name::Symbol)
     try

test/symmatrix.jl

@@ -0,0 +1,19 @@
+using Test
+using MultivariateMoments
+
+@testset "SymMatrix" begin
+    Q = SymMatrix([1, 2, 3], 2)
+    symmetric_setindex!(Q, 4, 1, 1)
+    @test Q.Q == [4, 2, 3]
+    symmetric_setindex!(Q, 5, 1, 2)
+    @test Q.Q == [4, 5, 3]
+    P = copy(Q)
+    @test P.n == 2
+    symmetric_setindex!(P, 6, 2, 2)
+    @test Q.Q == [4, 5, 3]
+    @test P.n == 2
+    @test P.Q == [4, 5, 6]
+    R = map(i -> i - 1, P)
+    @test R.n == 2
+    @test R.Q == [3, 4, 5]
+end