From d9b47348a246f0e8a774f2289b24bc12ee089fab Mon Sep 17 00:00:00 2001 From: "Documenter.jl" Date: Wed, 4 Oct 2023 12:37:23 +0000 Subject: [PATCH] build based on 3c5f4ba --- previews/PR73/.documenter-siteinfo.json | 2 +- previews/PR73/atoms/index.html | 54 ++++++++++++------------- previews/PR73/index.html | 2 +- previews/PR73/moments/index.html | 6 +-- 4 files changed, 32 insertions(+), 32 deletions(-) diff --git a/previews/PR73/.documenter-siteinfo.json b/previews/PR73/.documenter-siteinfo.json index fe58334..2c9d79e 100644 --- a/previews/PR73/.documenter-siteinfo.json +++ b/previews/PR73/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.9.3","generation_timestamp":"2023-10-04T12:36:57","documenter_version":"1.1.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.9.3","generation_timestamp":"2023-10-04T12:37:16","documenter_version":"1.1.0"}} \ No newline at end of file diff --git a/previews/PR73/atoms/index.html b/previews/PR73/atoms/index.html index af49899..cbfc14b 100644 --- a/previews/PR73/atoms/index.html +++ b/previews/PR73/atoms/index.html @@ -2,37 +2,37 @@ Atoms extraction · MultivariateMoments
GitHub

Atoms extration

Vectorized matrix

MultivariateMoments.SymMatrix — Type
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.

source
MultivariateMoments.VectorizedHermitianMatrix — Type
struct VectorizedHermitianMatrix{T} <: AbstractMatrix{Complex{T}}
+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.

source
MultivariateMoments.VectorizedHermitianMatrix — Type
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.

source
MultivariateMoments.square_getindex — Function
square_getindex!(Q::SymMatrix, I)

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

source
square_getindex!(Q::VectorizedHermitianMatrix, I)

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

source
MultivariateMoments.symmetric_setindex! — Function
symmetric_setindex!(Q::SymMatrix, value, i::Integer, j::Integer)

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

source
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.

source

Moment matrix

MultivariateMoments.MomentMatrix — Type
mutable struct MomentMatrix{T, B<:MultivariateBases.AbstractPolynomialBasis, MT<:AbstractMatrix{T}} <: AbstractMeasureLike{T}
+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.

source
MultivariateMoments.square_getindex — Function
square_getindex!(Q::SymMatrix, I)

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

source
square_getindex!(Q::VectorizedHermitianMatrix, I)

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

source
MultivariateMoments.symmetric_setindex! — Function
symmetric_setindex!(Q::SymMatrix, value, i::Integer, j::Integer)

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

source
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.

source

Moment matrix

MultivariateMoments.MomentMatrix — Type
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.

source

Atomic measure

MultivariateMoments.WeightedDiracMeasure — Type
struct WeightedDiracMeasure{T}
+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.

source

Atomic measure

MultivariateMoments.AtomicMeasure — Type
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.

source

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! — Function
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.

source
MultivariateMoments.atomic_measure — Function
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
+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.

source

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! — Function
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.

source
MultivariateMoments.atomic_measure — Function
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.

source

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 — Type
struct FlatExtension{
+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.

source

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 — Type
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.

source

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 — Type
struct MacaulayNullspace{T,MT<:AbstractMatrix{T},BT}
+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.

source

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 — Type
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.

source
MultivariateMoments.ImageSpaceSolver — Type
struct ImageSpaceSolver{A<:LowRankLDLTAlgorithm,S<:MacaulayNullspaceSolver}
+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.

source
MultivariateMoments.ImageSpaceSolver — Type
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.

source

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.low_rank_ldlt — Function
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.

source

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.low_rank_ldlt — Function
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.

source
MultivariateMoments.LowRankLDLT — Type
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.

source

The choice of cutoff between the significant and negligeable eigen/singular values is parametrized by the following interface:

MultivariateMoments.accuracy — Function
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.

source
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.

source
MultivariateMoments.doubt — Function
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.

source

The rank check can be chosen among the following:

MultivariateMoments.UserRank — Type
struct UserRank <: RankCheck
+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.

source

The choice of cutoff between the significant and negligeable eigen/singular values is parametrized by the following interface:

MultivariateMoments.accuracy — Function
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.

source
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.

source
MultivariateMoments.doubt — Function
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.

source

The rank check can be chosen among the following:

MultivariateMoments.UserRank — Type
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:
@@ -41,10 +41,10 @@
  > 0.05
    1.0e-5
    5.0e-6
-3
source
MultivariateMoments.FixedRank — Type
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
source
MultivariateMoments.FixedRanks — Type
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])
@@ -54,27 +54,27 @@
 2
 
 julia> rank_from_singular_values([1, 1e-1, 5e-2, 1e-5, 5e-6], check)
-3
source
MultivariateMoments.AbsoluteRankTol — Type
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
source
MultivariateMoments.LeadingRelativeRankTol — Type
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
source
MultivariateMoments.DifferentialRankTol — Type
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
source
MultivariateMoments.LargestDifferentialRank — Type
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
source

Given the MacaulayNullspace, there are two approaches implemented to obtain the moment matrices:

Given the MacaulayNullspace, there are two approaches implemented to obtain the moment matrices:

MultivariateMoments.ShiftNullspace — Type
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.

source
MultivariateMoments.Echelon — Type
struct Echelon{D,S<:Union{Nothing,SemialgebraicSets.AbstractAlgebraicSolver}} <:
+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.

source
MultivariateMoments.Echelon — Type
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.

Note

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.

source

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 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.

source
MultivariateMoments.StaircaseDependence — Type
struct StaircaseDependence
+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.

Note

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.

source

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 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.

source
MultivariateMoments.StaircaseDependence — Type
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.

source
MultivariateMoments.BasisDependence — Type
struct BasisDependence{D,B<:MB.AbstractPolynomialBasis}
+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.

source
MultivariateMoments.BasisDependence — Type
struct BasisDependence{D,B<:MB.AbstractPolynomialBasis}
     basis::B
     dependence::Vector{D}
-end

The dependence between the elements of a basis.

Tip

If the number of variables is 2 or 3, it can be visualized with Plots.jl.

source

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 — Type
struct BorderBasis{D,T,MT<:AbstractMatrix{T},B}
+end

The dependence between the elements of a basis.

Tip

If the number of variables is 2 or 3, it can be visualized with Plots.jl.

source

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 — Type
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.

source
MultivariateMoments.AlgebraicBorderSolver — Type
Base.@kwdef struct AlgebraicBorderSolver{
+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.

source
MultivariateMoments.AlgebraicBorderSolver — Type
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}`.

source
MultivariateMoments.AlgebraicFallbackBorderSolver — Type
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.

source

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.

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 — Type
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.

source
MultivariateMoments.MomentVectorWeightSolver — Type
struct MomentVectorWeightSolver{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.

source
MultivariateMoments.MomentVectorWeightSolver — Type
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.

source
+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.

source diff --git a/previews/PR73/index.html b/previews/PR73/index.html index 5616b0a..662fc69 100644 --- a/previews/PR73/index.html +++ b/previews/PR73/index.html @@ -1,2 +1,2 @@ -Index · MultivariateMoments
GitHub
+Index · MultivariateMoments
GitHub
diff --git a/previews/PR73/moments/index.html b/previews/PR73/moments/index.html index e53a428..380057e 100644 --- a/previews/PR73/moments/index.html +++ b/previews/PR73/moments/index.html @@ -1,4 +1,4 @@ -Moments and expectation · MultivariateMoments
GitHub

Moments and expectation

Moment

Given a measure $\mu$ and a monomial $m$, the moment $m$ of the measure is defined by the expectation $\mathbb{E}_\mu[m]$. Given a monomial and a value for the moment, a moment can be created using the moment function

The moment function returns an AbstractMoment which is a subtype of AbstractMomentLike. An AbstractMomentLike is a type that can act like an AbstractMoment (it is similar to MultivariatePolynomials' AbstractMonomialLike, AbstractTermLike and AbstractPolynomialLike), that is, it implements the following two functions

Measure

Given a monomials and a values for the moments, a "measure" can be created using the measure function

MultivariateMoments.measure — Function
measure(a::AbstractVector{T}, X::AbstractVector{<:AbstractMonomial}; rtol=Base.rtoldefault(T), atol=zero(T))

Creates a measure with moments moment(a[i], X[i]) for each i. An error is thrown if there exists i and j such that X[i] == X[j] but !isapprox(a[i], a[j]; rtol=rtol, atol=atol).

source

The measure function returns an AbstractMeasure which is a subtype of AbstractMeasureLike. Note that it does not actually compute the probability density function of a measure having these moments, it simply stores a vector of moments belonging to a hypothetical measure. However, it acts like a measure when taking its scalar product with a polynomial.

An AbstractMeasureLike is a type that can act like an AbstractMeasure, that is, it implements the following two functions

MultivariatePolynomials.variables — Method
variables(μ::AbstractMeasureLike)

Returns the variables of μ in decreasing order. Just like in MultivariatePolynomials, it could contain variables of zero degree in every monomial.

source

The moments of the dirac measure for a vector of monomials can be obtained by the dirac function

MultivariateMoments.dirac — Function
dirac(X::AbstractVector{<:AbstractMoment}, s::AbstractSubstitution...)

Creates the dirac measure by evaluating the moments of X using s.

Examples

Calling dirac([x*y, x*y^2], x=>3, y=>2) should the measure with moment x*y of value 6 and moment x*y^2 of value 12.

source

Expectation

The expectation of polynomial with respect to a measure can be computed either using MultivariateMoments.expectation or using the Base.dot scalar product.

MultivariateMoments.expectation — Function
MultivariateMoments.expectation(μ::AbstractMeasureLike, p::AbstractPolynomialLike)
-MultivariateMoments.expectation(p::AbstractPolynomialLike, μ::AbstractMeasureLike)

Computes the expectation $\mathbb{E}_{\mu}[p]$.

source
+Moments and expectation · MultivariateMoments
GitHub

Moments and expectation

Moment

Given a measure $\mu$ and a monomial $m$, the moment $m$ of the measure is defined by the expectation $\mathbb{E}_\mu[m]$. Given a monomial and a value for the moment, a moment can be created using the moment function

The moment function returns an AbstractMoment which is a subtype of AbstractMomentLike. An AbstractMomentLike is a type that can act like an AbstractMoment (it is similar to MultivariatePolynomials' AbstractMonomialLike, AbstractTermLike and AbstractPolynomialLike), that is, it implements the following two functions

Measure

Given a monomials and a values for the moments, a "measure" can be created using the measure function

MultivariateMoments.measure — Function
measure(a::AbstractVector{T}, X::AbstractVector{<:AbstractMonomial}; rtol=Base.rtoldefault(T), atol=zero(T))

Creates a measure with moments moment(a[i], X[i]) for each i. An error is thrown if there exists i and j such that X[i] == X[j] but !isapprox(a[i], a[j]; rtol=rtol, atol=atol).

source

The measure function returns an AbstractMeasure which is a subtype of AbstractMeasureLike. Note that it does not actually compute the probability density function of a measure having these moments, it simply stores a vector of moments belonging to a hypothetical measure. However, it acts like a measure when taking its scalar product with a polynomial.

An AbstractMeasureLike is a type that can act like an AbstractMeasure, that is, it implements the following two functions

MultivariatePolynomials.variables — Method
variables(μ::AbstractMeasureLike)

Returns the variables of μ in decreasing order. Just like in MultivariatePolynomials, it could contain variables of zero degree in every monomial.

source

The moments of the dirac measure for a vector of monomials can be obtained by the dirac function

MultivariateMoments.dirac — Function
dirac(X::AbstractVector{<:AbstractMoment}, s::AbstractSubstitution...)

Creates the dirac measure by evaluating the moments of X using s.

Examples

Calling dirac([x*y, x*y^2], x=>3, y=>2) should the measure with moment x*y of value 6 and moment x*y^2 of value 12.

source

Expectation

The expectation of polynomial with respect to a measure can be computed either using MultivariateMoments.expectation or using the Base.dot scalar product.

MultivariateMoments.expectation — Function
MultivariateMoments.expectation(μ::AbstractMeasureLike, p::AbstractPolynomialLike)
+MultivariateMoments.expectation(p::AbstractPolynomialLike, μ::AbstractMeasureLike)

Computes the expectation $\mathbb{E}_{\mu}[p]$.

source