Atoms extration

struct SymMatrix{T} <: AbstractMatrix{T}
    Q::Vector{T}
    n::Int
end

Symmetric $n \times n$ matrix storing the vectorized upper triangular part of the matrix in the Q vector (this corresponds to the vectorized form of MathOptInterface.PositiveSemidefiniteConeTriangle). It implement the AbstractMatrix interface except for setindex! as it might break its symmetry. The symmetric_setindex! function should be used instead.

struct VectorizedHermitianMatrix{T} <: AbstractMatrix{Complex{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! function should be used instead.

square_getindex!(Q::SymMatrix, I)

Return the SymMatrix corresponding to Q[I, I].

square_getindex!(Q::VectorizedHermitianMatrix, I)

Return the VectorizedHermitianMatrix corresponding to Q[I, I].

symmetric_setindex!(Q::SymMatrix, value, i::Integer, j::Integer)

Set Q[i, j] and Q[j, i] to the value value.

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.

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.

moment_matrix(μ::Measure, x)

Creates a matrix the moment matrix for the moment matrix $x x^\top$ using the moments of μ.

struct WeightedDiracMeasure{T}
    center::Vector{T}
    weight::T
end

Represents a weighted dirac measure with centered at center and with weight weight.

struct AtomicMeasure{T, AT, V} <: AbstractMeasureLike{T}
    variables::V                           # Vector/Tuple of variables
    atoms::Vector{WeightedDiracMeasure{T, AT}} # Atoms of the measure
end

An measure is said to be atomic if it is a sum of weighted dirac measures. For instance, $\eta = 2 \delta_{(1, 0)} + 3 \delta_{(1/2, 1/2)}$ is an atomic measure since it is a sum of the diracs centered at $(1, 0)$ and $(1/2, 1/2)$ and weighted respectively by 2 and 3. That is, $\mathbb{E}_{\eta}[p] = 2p(1, 0) + 3p(1/2, 1/2)$.

The AtomicMeasure struct stores an atomic measure where variables contains the variables and atoms contains atoms of the measure.

function compute_support!(ν::MomentMatrix, rank_check, solver) end

Computes the support field of ν wth solver using a low-rank decomposition with the rank evaluated with rank_check.

atomic_measure(ν::MomentMatrix, rank_check::RankCheck, [dec::LowRankLDLTAlgorithm], [solver::SemialgebraicSets.AbstractAlgebraicSolver])

Return an AtomicMeasure with the atoms of ν if it is atomic or nothing if ν is not atomic. The rank_check and dec parameters are passed as is to the low_rank_ldlt function. By default, dec is an instance of SVDLDLT. The extraction relies on the solution of a system of algebraic equations. using solver. For instance, given a MomentMatrix, μ, the following extract atoms using a rank_check of 1e-4 for the low-rank decomposition and homotopy continuation to solve the obtained system of algebraic equations.

using HomotopyContinuation
solver = SemialgebraicSetsHCSolver(; compile = false)
atoms = atomic_measure(ν, 1e-4, solver)

If no solver is given, the default solver of SemialgebraicSets is used which currently computes the Gröbner basis, then the multiplication matrices and then the Schur decomposition of a random combination of these matrices. For floating point arithmetics, homotopy continuation is recommended as it is more numerically stable than Gröbner basis computation.

struct FlatExtension{
    MMS<:SemialgebraicSets.AbstractMultiplicationMatricesSolver,
}
    multiplication_matrices_solver::MMS
end

Given a moment matrix satisfying the flat extension property described in [L09, Section 5.3], computes the multiplication matrices using the formula given in [L09, Lemma 6.21] or [LLR08, Section 4.4.4]. Given the multiplication matrices, the solutions are obtained with multiplication_matrices_solver.

[L09] Laurent, Monique. Sums of squares, moment matrices and optimization over polynomials. Emerging applications of algebraic geometry (2009): 157-270.

[LLR08] Lasserre, Jean Bernard and Laurent, Monique, and Rostalski, Philipp. Semidefinite characterization and computation of zero-dimensional real radical ideals. Foundations of Computational Mathematics 8 (2008): 607-647.

struct MacaulayNullspace{T,MT<:AbstractMatrix{T},BT}
    matrix::MT
    basis::BT
    accuracy::T
end

This matrix with rows indexed by basis can either be interpreted as the right null space of a Macaulay matrix with columns indexed by basis (resp. or the image space of a moment matrix with rows and columns indexed by basis). The value of matrix[i, j] should be interpreted as the value of the ith element of basis for the jth generator of the null space (resp. image) space.

struct ImageSpaceSolver{A<:LowRankLDLTAlgorithm,S<:MacaulayNullspaceSolver}
    ldlt::A
    null::S
end

Computes the image space of the moment matrix using ldlt and then solve it with null.

LowRankLDLTAlgorithm

Method for computing a $r \times n$ matrix U of a $n \times n$ rank $r$ psd matrix Q such that Q = U'U.

ShiftCholeskyLDLT <: LowRankLDLTAlgorithm

Shift the matrix by shift times the identity matrix before cholesky.

SVDLDLT <: LowRankLDLTAlgorithm

Use SVD decomposition.

MultivariateMoments.low_rank_ldlt(Q::AbstractMatrix, dec::LowRankLDLTAlgorithm, ranktol)

Returns a $r \times n$ matrix $U$ of a $n \times n$ rank $r$ positive semidefinite matrix $Q$ such that $Q = U^\top U$. The rank of $Q$ is the number of singular values larger than ranktol${} \cdot \sigma_1$ where $\sigma_1$ is the largest singular value.

struct LowRankLDLT{T,Tr,C<:RankCheck}
    L::Matrix{T}
    singular_values::Vector{Tr}
    rank_check::C
end

Low-Rank cholesky decomposition L * Diagonal(singular_values) * L' of size (n, r) of a n x n matrix with singular values singular_values[1] > ... > singular_values[n]. The rank was chosen given singular_values using rank_check via the rank_from_singular_values function.

abstract type RankCheck end

Method for computing the rank with rank_from_singular_values based on a list of singular values.

rank_from_singular_values(σ, check::RankCheck)

Return the rank of a matrix with singular values σ (in decreasing order) using check.

accuracy(σ, r, check::RankCheck)

Returns a value measuring the accuracy of the rank check check returning rank r. This is used by Echelon to determine the accuracy to use for the Gaussian elimination.

accuracy(chol::LowRankLDLT)

Return the ratio rtol = σ[r+1]/σ[1] where σ is the vector of singular values of the matrix for which chol is the rank-r Cholesky decomposition. This is a good relative tolerance to use with the matrix as σ[r+1] is the first singular value that was discarded.

doubt(σ, check::RankCheck)

Returns a value measuring the doubt of the rank check check. Lower values means more doubt so less certainty. This is used by FallbackRank for deciding whether the fallback should be used.

struct UserRank <: RankCheck
    pagesize::Int
end

The user chooses the rank given the singular values in a REPL.TerminalMenus.RadioMenu of page size pagesize.

Example

julia> rank_from_singular_values([1, 1e-1, 5e-2, 1e-5, 5e-6], UserRank())
Choose the last significant singular value:
   1.0
   0.1
 > 0.05
   1.0e-5
   5.0e-6
3

struct FixedRank <: RankCheck
    r::Int
end

The rank is hardcoded to r, independently of the singular values.

Example

julia> rank_from_singular_values([1, 1e-1, 5e-2, 1e-5, 5e-6], FixedRank(3))
3

mutable struct FixedRanks <: RankCheck
    r::Vector{Int}
    current::Int
end

The ith rank is hardcoded to r[i], independently of the singular values. The field current indicates how many ranks have already been asked. When current is length(r), no rank can be asked anymore.

Example

julia> check = FixedRanks([2, 3])
FixedRanks([2, 3], 0)

julia> rank_from_singular_values([1, 1e-1, 5e-5, 1e-5, 5e-6], check)
2

julia> rank_from_singular_values([1, 1e-1, 5e-2, 1e-5, 5e-6], check)
3

struct AbsoluteRankTol{T} <: RankCheck
    tol::T
end

The rank is the number of singular values that are strictly larger than tol.

Example

julia> rank_from_singular_values([1, 1e-1, 5e-2, 1e-5, 5e-6], AbsoluteRankTol(1e-4))
3

struct LeadingRelativeRankTol{T} <: RankCheck
    tol::T
end

The rank is the number of singular values that are strictly larger than tol * maximum(σ) where maximum(σ) is the largest singular value.

Example

When the matrix is obtained from a homogeneous problem where the scaling is irrelevant, LeadingRelativeRankTol may be preferable to AbsoluteRankTol as shown below

julia> rank_from_singular_values(1e6 * [1, 1e-1, 5e-2, 1e-5, 5e-6], AbsoluteRankTol(1e-4))
5

julia> rank_from_singular_values(1e6 * [1, 1e-1, 5e-2, 1e-5, 5e-6], LeadingRelativeRankTol(1e-4))
3

struct DifferentialRankTol{T} <: RankCheck
    tol::T
end

The rank is the number of singular values before a singular value (not included) is tol times the previous one (included).

Example

It is sometimes difficult to figure out the tolerance to use in LeadingRelativeRankTol. For instance, choosing 1e-3 will consider 1e-3 in the example below as not part of the rank while DifferentialRankTol would include it because it is close to the previous singular value.

julia> rank_from_singular_values([1, 1e-1, 5e-2, 1e-2, 5e-3, 1e-3, 1e-6, 5e-7], LeadingRelativeRankTol(1e-3))
5

julia> rank_from_singular_values([1, 1e-1, 5e-2, 1e-2, 5e-3, 1e-3, 1e-6, 5e-7], DifferentialRankTol(1e-2))
6

struct LargestDifferentialRank <: RankCheck
end

The rank is the number of singular values until the singular value that has the largest ratio with the next singular value.

Example

It is sometimes difficult to figure out the tolerance to use in DifferentialRankTol. For instance, choosing 1e-2 will consider 1e-2, 5e-2 and 1e-3 in the example below as not part of the rank while LargestDifferentialRank would include them because there is a largest gap later.

julia> rank_from_singular_values([1, 1e-2, 5e-2, 1e-3, 1e-6, 5e-7], DifferentialRankTol(1e-2))
1

julia> rank_from_singular_values([1, 1e-2, 5e-2, 1e-3, 1e-6, 5e-7], LargestDifferentialRank())
4

struct FallbackRank{T,D,F} <: RankCheck
    tol::T
    default::D
    fallback::F
end

Defaults to checking the rank with default and falls back to fallback if the doubt is strictly larger than tol. By default, fallback is UserRank.

Example

The advantage of UserRank is that the user get to see if the rank check is ambiguous and act accordingly. The downside is that it might be cumbersome if there are many rank checks to do. With FallbackRank, the user only has to sort out the tricky ones. In the example below, the first example is handled by LargestDifferentialRank. For the second one, the user sees that this is a tricky choice and can manually choose one of the two ranks, then see the result of the rest of his code using this value of the code and then choose the other rank and see the impact of this different choice.

julia> check = FallbackRank(1e-1, LargestDifferentialRank())
FallbackRank{Float64, LargestDifferentialRank, UserRank}(0.1, LargestDifferentialRank(), UserRank(8))

julia> rank_from_singular_values([1, 1e-2, 5e-2, 1e-3, 1e-6, 5e-7], check)
4

julia> rank_from_singular_values([1, 1e-2, 5e-2, 1e-3, 5e-6, 5e-7], check)
Choose the last significant singular value:
 > 1.0
   0.01
   0.05
   0.001
   5.0e-6
   5.0e-7
1

julia> rank_from_singular_values([1, 1e-2, 5e-2, 1e-3, 5e-6, 5e-7], check)
Choose the last significant singular value:
   1.0
   0.01
   0.05
 > 0.001
   5.0e-6
   5.0e-7
4

struct ShiftNullspace{D,C<:RankCheck} <: MacaulayNullspaceSolver
    check::C
end

From the MacaulayNullspace, computes multiplication matrices by exploiting the shift property of the rows [DBD12]. The rank check check is used to determine the standard monomials among the row indices of the null space.

[DBD12] Dreesen, Philippe, Batselier, Kim, and De Moor, Bart. Back to the roots: Polynomial system solving, linear algebra, systems theory. IFAC Proceedings Volumes 45.16 (2012): 1203-1208.

struct Echelon{D,S<:Union{Nothing,SemialgebraicSets.AbstractAlgebraicSolver}} <:
    MacaulayNullspaceSolver
    solver::S
end

Given a MacaulayNullspace, computes its echelon form (corresponding to the Canonical Null Space of [D13]) with Gaussian elimination. From this echelon form, the left null space can easily be computed using using [HL05, (8)]. This left null space forms a system of polynomial equations.

Let B be the set of monomials corresponding to the rows of the pivots of this echelon form. If the moment matrix satisfies the flat extension property described in [L09, Section 5.3] then all monomials of the border of B (as defined in [LLR08, (2.3)]) will correspond to a row of of the matrix. In that case, the polynomial of the system obtained by [HL05, (8)] form a rewriting family for B [L09, (2.16)] a.k.a. a B-border prebasis (as defined in [LLR08, (2.4)]). Therefore, they can be used to compute multiplication matrices.

[HL05] Henrion, D. & Lasserre, J-B. Detecting Global Optimality and Extracting Solutions of GloptiPoly 2005

[D13] Dreesen, Philippe. Back to the Roots: Polynomial System Solving Using Linear Algebra Ph.D. thesis (2013)

[L09] Laurent, Monique. Sums of squares, moment matrices and optimization over polynomials. Emerging applications of algebraic geometry (2009): 157-270.

[LLR08] Lasserre, Jean Bernard and Laurent, Monique, and Rostalski, Philipp. Semidefinite characterization and computation of zero-dimensional real radical ideals. Foundations of Computational Mathematics 8 (2008): 607-647.

[S04] Stetter, Hans J. Numerical polynomial algebra. Society for Industrial and Applied Mathematics, 2004.

@enum LinearDependence INDEPENDENT TRIVIAL DEPENDENT

Linear dependence of the element of a basis representing the indices of the rows of a [MacaulayNullspace]. DEPENDENT indicates that it is linearly dependent to rows that appear earlier in the basis. TRIVIAL indicates that the element was not in the original basis so it is trivially independent.

struct StaircaseDependence
    standard_or_corner::Bool
    linear::LinearDependence
end

Dependence of the element of a basis representing the indices of the rows of a [MacaulayNullspace]. If the field standard_or_corner is true then the elements is either standard or is a corner (depending on the linear dependence encoded in the linear field). Otherwise, it is a border that is not a corner or it is not even a border. See LinearDependence for the linear field.

struct BasisDependence{D,B<:MB.AbstractPolynomialBasis}
    basis::B
    dependence::Vector{D}
end

The dependence between the elements of a basis.

struct BorderBasis{D,T,MT<:AbstractMatrix{T},B}
    dependence::BasisDependence{D,B}
    matrix::MT
end

This matrix with rows indexed by standard and columns indexed by border a standard-border basis of the ideal border .- matrix' * standard [LLR08, Section 2.5]. For solving this with a multiplication matrix solver, it is necessary for the basis border to be a superset of the set of corners of standard and it is sufficient for it to be the border of standard.

[LLR08] Lasserre, Jean Bernard and Laurent, Monique, and Rostalski, Philipp. Semidefinite characterization and computation of zero-dimensional real radical ideals. Foundations of Computational Mathematics 8 (2008): 607-647.

Base.@kwdef struct AlgebraicBorderSolver{
    D,
    A<:Union{Nothing,SS.AbstractGröbnerBasisAlgorithm},
    S<:Union{Nothing,SS.AbstractAlgebraicSolver},
}
    algorithm::A = nothing
    solver::S = nothing
end

Solve a border basis using algorithm and solver by first converting it to aBorderBasis{D}`.

struct AlgebraicFallbackBorderSolver{
    M<:Union{Nothing,SS.AbstractMultiplicationMatricesSolver},
    S<:AlgebraicBorderSolver,
}
    multiplication_matrices_solver::M
    algebraic_solver::S
end

Solve with multiplication_matrices_solver and if it fails, falls back to solving the algebraic system formed by the border basis with algebraic_solver.

struct MomentMatrixWeightSolver
    rtol::T
    atol::T
end

Given a moment matrix ν and the atom centers, determine the weights by solving a linear system over all the moments of the moment matrix, keeping duplicates (e.g., entries corresponding to the same monomial).

If the moment values corresponding to the same monomials are known to be equal prefer MomentVectorWeightSolver instead.

struct MomentVectorWeightSolver{T}
    rtol::T
    atol::T
end

Given a moment matrix ν and the atom centers, first convert the moment matrix to a vector of moments, using measure(ν; rtol=rtol, atol=atol) and then determine the weights by solving a linear system over the monomials obtained.

If the moment values corresponding to the same monomials can have small differences, measure can throw an error if rtol and atol are not small enough. Alternatively to tuning these tolerances MomentVectorWeightSolver can be used instead.