Atoms extration
Vectorized matrix
MultivariateMoments.SymMatrix — Typestruct SymMatrix{T} <: AbstractMatrix{T}
+ Q::Vector{T}
+ n::Int
+endSymmetric $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.
MultivariateMoments.VectorizedHermitianMatrix — Typestruct VectorizedHermitianMatrix{T} <: AbstractMatrix{Complex{T}}
+ Q::Vector{T}
+ n::Int
+endHermitian $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.
MultivariateMoments.square_getindex — Functionsquare_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].
MultivariateMoments.symmetric_setindex! — Functionsymmetric_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.
Moment matrix
MultivariateMoments.MomentMatrix — Typemutable struct MomentMatrix{T, B<:MultivariateBases.AbstractPolynomialBasis, MT<:AbstractMatrix{T}} <: AbstractMeasureLike{T}
+ Q::MT
+ basis::B
+ support::Union{Nothing, AlgebraicSet}
+endMeasure $\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.
MultivariateMoments.moment_matrix — Functionmoment_matrix(μ::Measure, x)Creates a matrix the moment matrix for the moment matrix $x x^\top$ using the moments of μ.
Atomic measure
MultivariateMoments.WeightedDiracMeasure — Typestruct WeightedDiracMeasure{T}
+ center::Vector{T}
+ weight::T
+endRepresents a weighted dirac measure with centered at center and with weight weight.
MultivariateMoments.AtomicMeasure — Typestruct AtomicMeasure{T, AT, V} <: AbstractMeasureLike{T}
+ variables::V # Vector/Tuple of variables
+ atoms::Vector{WeightedDiracMeasure{T, AT}} # Atoms of the measure
+endAn 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.
Atoms extraction
Given a MomentMatrix containing the moments of an atomic measure, atomic_measure attempts to recover the dirac centers and weights by first computing an algebraic system with the atom centers as solution with compute_support!.
MultivariateMoments.compute_support! — Functionfunction compute_support!(ν::MomentMatrix, rank_check, solver) endComputes the support field of ν wth solver using a low-rank decomposition with the rank evaluated with rank_check.
MultivariateMoments.atomic_measure — Functionatomic_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.
For the first step of compute_support!, there are two approaches. The first one is to exploit the flat extension to directly obtain the multiplication matrices.
MultivariateMoments.FlatExtension — Typestruct FlatExtension{
+ MMS<:SemialgebraicSets.AbstractMultiplicationMatricesSolver,
+}
+ multiplication_matrices_solver::MMS
+endGiven 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.
The second approach is to first obtain the image space of the moment matrix, represented as a MacaulayNullspace and to then compute the multiplication matrices from this image space. This is implemented by the ImageSpaceSolver.
MultivariateMoments.MacaulayNullspace — Typestruct MacaulayNullspace{T,MT<:AbstractMatrix{T},BT}
+ matrix::MT
+ basis::BT
+ accuracy::T
+endThis 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.
MultivariateMoments.ImageSpaceSolver — Typestruct ImageSpaceSolver{A<:LowRankLDLTAlgorithm,S<:MacaulayNullspaceSolver}
+ ldlt::A
+ null::S
+endComputes the image space of the moment matrix using ldlt and then solve it with null.
This image space is obtained from a low rank LDLT decomposition of the moment matrix. This decomposition can either be obtained by a cholesky or SVD decomposition from which we remove the rows corresponding to the negligeable eigen/singular values.
MultivariateMoments.LowRankLDLTAlgorithm — TypeLowRankLDLTAlgorithmMethod for computing a $r \times n$ matrix U of a $n \times n$ rank $r$ psd matrix Q such that Q = U'U.
MultivariateMoments.ShiftCholeskyLDLT — TypeShiftCholeskyLDLT <: LowRankLDLTAlgorithmShift the matrix by shift times the identity matrix before cholesky.
MultivariateMoments.SVDLDLT — TypeSVDLDLT <: LowRankLDLTAlgorithmUse SVD decomposition.
MultivariateMoments.low_rank_ldlt — FunctionMultivariateMoments.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.
MultivariateMoments.LowRankLDLT — Typestruct LowRankLDLT{T,Tr,C<:RankCheck}
+ L::Matrix{T}
+ singular_values::Vector{Tr}
+ rank_check::C
+endLow-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.
The choice of cutoff between the significant and negligeable eigen/singular values is parametrized by the following interface:
MultivariateMoments.RankCheck — Typeabstract type RankCheck endMethod for computing the rank with rank_from_singular_values based on a list of singular values.
MultivariateMoments.rank_from_singular_values — Functionrank_from_singular_values(σ, check::RankCheck)Return the rank of a matrix with singular values σ (in decreasing order) using check.
MultivariateMoments.accuracy — Functionaccuracy(σ, 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.
MultivariateMoments.doubt — Functiondoubt(σ, 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.
The rank check can be chosen among the following:
MultivariateMoments.UserRank — Typestruct UserRank <: RankCheck
+ pagesize::Int
+endThe 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
+3MultivariateMoments.FixedRank — Typestruct FixedRank <: RankCheck
+ r::Int
+endThe 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))
+3MultivariateMoments.FixedRanks — Typemutable struct FixedRanks <: RankCheck
+ r::Vector{Int}
+ current::Int
+endThe 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)
+3MultivariateMoments.AbsoluteRankTol — Typestruct AbsoluteRankTol{T} <: RankCheck
+ tol::T
+endThe 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))
+3MultivariateMoments.LeadingRelativeRankTol — Typestruct LeadingRelativeRankTol{T} <: RankCheck
+ tol::T
+endThe 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))
+3MultivariateMoments.DifferentialRankTol — Typestruct DifferentialRankTol{T} <: RankCheck
+ tol::T
+endThe 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))
+6MultivariateMoments.LargestDifferentialRank — Typestruct LargestDifferentialRank <: RankCheck
+endThe 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())
+4MultivariateMoments.FallbackRank — Typestruct FallbackRank{T,D,F} <: RankCheck
+ tol::T
+ default::D
+ fallback::F
+endDefaults 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
+4Given the MacaulayNullspace, there are two approaches implemented to obtain the moment matrices:
MultivariateMoments.ShiftNullspace — Typestruct ShiftNullspace{D,C<:RankCheck} <: MacaulayNullspaceSolver
+ check::C
+endFrom 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.
MultivariateMoments.Echelon — Typestruct Echelon{D,S<:Union{Nothing,SemialgebraicSets.AbstractAlgebraicSolver}} <:
+ MacaulayNullspaceSolver
+ solver::S
+endGiven 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.
In the context of compute_support!, if the MacaulayNullspace was obtained using SVDLDLT, the left null space can easily be obtained from the singular vectors corresponding to the negligeable singular values that were removed. However, as mentioned [LLR08, Section 4.4.5], these will usually give an overdetermined bases. As shown in [S04, Section 10.8.2], it is desirable to avoid overdetermined bases because it could lead to inconsistencies in the basis for numerical reasons. For this reason, this method computes the left null space from the echelon form of the significant singular values instead.
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.
The Echelon uses the RowEchelon package to determine the independent rows (which is not numerically stable) while the ShiftNullspace uses RankChecks with the singular values so it should have better numerical behavior. They can either simply distinguish the dependency of rows with LinearDependence or use a sieve with StaircaseDependence to save some the computation of the singular values for some submatrices:
MultivariateMoments.LinearDependence — Type@enum LinearDependence INDEPENDENT TRIVIAL DEPENDENTLinear 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.
MultivariateMoments.StaircaseDependence — Typestruct StaircaseDependence
+ standard_or_corner::Bool
+ linear::LinearDependence
+endDependence 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.
MultivariateMoments.BasisDependence — Typestruct BasisDependence{D,B<:MB.AbstractPolynomialBasis}
+ basis::B
+ dependence::Vector{D}
+endThe dependence between the elements of a basis.
If the number of variables is 2 or 3, it can be visualized with Plots.jl.
The relationship between the dependent and the independent rows are then stored in a BorderBasis. By default, calling solve with only a BorderBasisas argument or providing aSemialgebraicSets.AbstractMultiplicationMatricesSolveras second argument will try to construct the moment matrices from these and find the solutions from these moment matrices. Alternatively, give an [AlgebraicBorderSolver](@ref) or a [AlgebraicFallbackBorderSolver`](@ref) as second argument to solve the system formed by these relations with these instead.
MultivariateMoments.BorderBasis — Typestruct BorderBasis{D,T,MT<:AbstractMatrix{T},B}
+ dependence::BasisDependence{D,B}
+ matrix::MT
+endThis 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.
MultivariateMoments.AlgebraicBorderSolver — TypeBase.@kwdef struct AlgebraicBorderSolver{
+ D,
+ A<:Union{Nothing,SS.AbstractGröbnerBasisAlgorithm},
+ S<:Union{Nothing,SS.AbstractAlgebraicSolver},
+}
+ algorithm::A = nothing
+ solver::S = nothing
+endSolve a border basis using algorithm and solver by first converting it to aBorderBasis{D}`.
MultivariateMoments.AlgebraicFallbackBorderSolver — Typestruct AlgebraicFallbackBorderSolver{
+ M<:Union{Nothing,SS.AbstractMultiplicationMatricesSolver},
+ S<:AlgebraicBorderSolver,
+}
+ multiplication_matrices_solver::M
+ algebraic_solver::S
+endSolve with multiplication_matrices_solver and if it fails, falls back to solving the algebraic system formed by the border basis with algebraic_solver.
Once the center of the atoms are determined, a linear system is solved to determine the weights corresponding to each dirac. By default, MomentMatrixWeightSolver is used by atomic_measure so that if there are small differences between moment values corresponding to the same monomial in the matrix (which can happen if these moments were computed numerically by a semidefinite proramming solvers, e.g., with SumOfSquares), the linear system handles that automatically.
MultivariateMoments.MomentMatrixWeightSolver — Typestruct MomentMatrixWeightSolver
+ rtol::T
+ atol::T
+endGiven 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.
MultivariateMoments.MomentVectorWeightSolver — Typestruct MomentVectorWeightSolver{T}
+ rtol::T
+ atol::T
+endGiven 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.