mainalgorithm¶
Main algorithms for recovering weights as eigenvalues from mixed moment matrices.
Example: Three circles on a line¶
Imports:
sage: from momentproblems import moment_functionals
sage: from momentproblems.mainalgorithm import *
sage: from momentproblems.mainalgorithm import _check_weights
sage: tol = 1e-7
sage: dim = 2
Construction of moment matrices of components and mixtures. Note that we have more pencils than needed (s > r):
sage: functionals = [moment_functionals.mk_moment_ellipse(center=xy, radius=(1,1)) for xy in [(0,0), (-3,3/2), (2,-1)]]
sage: weights = [[1, 1, 1], [2, 3, 6], [-1, 3, 5], [2, 1, -2]]
sage: Hs = [moment_functionals.hankel(f, dim=dim, d=6) for f in functionals]
sage: Ms = [Hs[0].parent().sum(weights[i][j] * Hs[j] for j in range(len(Hs))) for i in range(len(weights))]
Symbolically¶
sage: λ_ij = weights_from_moment_matrices_fixed_δ(Ms, d=4, δ=2, dim=dim, tol=tol)
sage: assert λ_ij.base_ring().is_exact()
sage: _check_weights(weights, λ_ij, tol)
True
sage: λ_ij = weights_from_moment_matrices(Ms, d=4, dim=dim, tol=tol)
sage: assert λ_ij.base_ring().is_exact()
sage: _check_weights(weights, λ_ij, tol)
True
In the case that δ=0, we find more common eigenvalues than desired:
sage: Hs2 = [moment_functionals.hankel(f, dim=dim, d=4) for f in functionals]
sage: Ms2 = [Hs2[0].parent().sum(weights[i][j] * Hs2[j] for j in range(len(Hs2))) for i in range(len(weights))]
sage: λ_ij = weights_from_moment_matrices_fixed_δ(Ms2, d=4, δ=0, dim=dim, tol=tol)
sage: λ_ij.ncols() > len(Hs2)
True
We can compute the eigenvector corresponding to the extraneous eigenvalue which exists due to the symmetry structure between the components:
sage: extraneous = ProjectiveSpace(len(Ms2)-1)(*sum(matrix(weights).columns()))
sage: from momentproblems.commoneigenspaces import common_eigenspaces_symbolic
sage: eigvals, eigmat = common_eigenspaces_symbolic(Ms2)
sage: monoms = vector(moment_functionals.monomial_basis(dim=dim, d=4))
sage: p = eigmat.column(eigvals.index(extraneous)) * monoms; p / p.lc()
y0^4 + 8*y0^3*y1 + 24*y0^2*y1^2 + 32*y0*y1^3 + 16*y1^4 - 5*y0^2 - 20*y0*y1 - 20*y1^2 + 25/8
For δ=1, the number of common eigenvalues is correct and we recover the
defining ideals:
sage: Hs2 = [moment_functionals.hankel(f, dim=dim, d=5) for f in functionals]
sage: Ms2 = [Hs2[0].parent().sum(weights[i][j] * Hs2[j] for j in range(len(Hs2))) for i in range(len(weights))]
sage: λ_ij = weights_from_moment_matrices_fixed_δ(Ms2, d=4, δ=1, dim=dim, tol=tol)
sage: λ_ij.ncols() == len(Hs2)
True
sage: _check_weights(weights, λ_ij, tol)
True
sage: Hs_rec = [Hj/Hj[0,0] for Hj in unmixed_moment_matrices(Ms2, λ_ij)]
sage: H0_rec, = [Hj for Hj in Hs_rec if Hj[0,1] == 0]
sage: monoms = vector(moment_functionals.monomial_basis(dim=dim, d=5))
sage: ker = (H0_rec.right_kernel_matrix() * monoms); ker.column()
[ -y0^4 - 2*y0^2*y1^2 - y1^4 + 1]
[-y0^5 - 2*y0^3*y1^2 - y0*y1^4 + y0]
[-y0^4*y1 - 2*y0^2*y1^3 - y1^5 + y1]
[ -y0^4 - y0^2*y1^2 + y0^2]
[ -y0^3*y1 - y0*y1^3 + y0*y1]
[ -y0^2*y1^2 - y1^4 + y1^2]
[ -y0^5 - y0^3*y1^2 + y0^3]
[ -y0^4*y1 - y0^2*y1^3 + y0^2*y1]
[ -y0^3*y1^2 - y0*y1^4 + y0*y1^2]
[ -y0^2*y1^3 - y1^5 + y1^3]
sage: ideal(ker.list()).groebner_basis()
[y0^2 + y1^2 - 1]
Numerically¶
sage: Hs = [H.change_ring(RDF) for H in Hs]
sage: Ms = [Hs[0].parent().sum(weights[i][j] * Hs[j] for j in range(len(Hs))) for i in range(len(weights))]
Reconstruction of weights via two different algorithms.
Ellipses are generated by polynomials of degree δ=2:
sage: λ_ij = weights_from_moment_matrices_fixed_δ(Ms, d=4, δ=2, dim=dim, tol=tol)
sage: _check_weights(weights, λ_ij, tol)
True
We have r+1=3 components:
sage: λ_ij = weights_from_moment_matrices(Ms, d=4, dim=dim, tol=tol)
sage: _check_weights(weights, λ_ij, tol)
True
Convex case¶
Reconstruction using information about convex hulls in case of positive-semidefinite matrices:
sage: weights = [[4, 2, 1], [2, 3, 6], [1, 3, 5], [2, 1, 2]]
sage: Hs = [moment_functionals.hankel(f, dim=dim, d=6) for f in functionals]
sage: Ms = [Hs[0].parent().sum(weights[i][j] * Hs[j] for j in range(len(Hs))) for i in range(len(weights))]
Check positive-semidefiniteness, so all eigenvalues are non-negative (Sage 9.3
has a method is_positive_semidefinite for this):
sage: is_positive_semidefinite = lambda M: all(abs(e.imag()) < tol and e.real() > -tol for e in M.change_ring(CDF).eigenvalues())
sage: all(is_positive_semidefinite(Mi) for Mi in Ms)
True
The reconstruction (symbolically):
sage: λ_ij = weights_from_moment_matrices_convex(Ms, tol=tol)
sage: _check_weights(weights, λ_ij, tol)
True
Numerically:
sage: Hs = [H.change_ring(RDF) for H in Hs]
sage: Ms = [Hs[0].parent().sum(weights[i][j] * Hs[j] for j in range(len(Hs))) for i in range(len(weights))]
sage: all(is_positive_semidefinite(Mi) for Mi in Ms)
True
sage: λ_ij = weights_from_moment_matrices_convex(Ms, tol=tol)
sage: _check_weights(weights, λ_ij, tol) # not tested (numeric convex hull not fully implemented)
True
With complex weights:
sage: weights = [[4, CDF(0,2), CDF(1,1)], [CDF(2,-2), 3, CDF(6,-1)], [CDF(1,1), CDF(3,4), 5], [CDF(2,-2), CDF(-1,-4), CDF(-2,2)]]
sage: Ms = [Hs[0].parent().sum(weights[i][j] * Hs[j] for j in range(len(Hs))) for i in range(len(weights))]
This requires a linear combination of the positive-semidefinite matrices Hs
with positive weights, so that the coordinates for the eigenvalues can be
choosen with respect to this matrix:
sage: H = sum(Hs).change_ring(CDF)
sage: λ_ij = weights_from_moment_matrices_convex(Ms, H, tol=tol)
sage: _check_weights(weights, λ_ij, tol) # not tested (numeric convex hull not fully implemented)
True
Example: straight lines on the torus¶
The components are of degree (1,1), (1,1) and (2,1). Therefore, max-degree 3 is enough to find common eigenspaces. As the components are of max-degree 2, an upper bound on δ is 2, so d+δ ≤ 5:
sage: functionals = [moment_functionals.mk_moment_torus_line(ratio=q, translate=c) for q, c in
....: [((1, 1), e^(i*pi/4)), ((-1, 1), 1), ((-2, 1), 1)]] #, ((2, 1), e^(i*1))]]
sage: D = 5
sage: Ts = [moment_functionals.toeplitz(f, dim=dim, d=D).change_ring(CDF) for f in functionals]
sage: weights = [[1, 1, 1], [2, 3, 6], [-1, 3, 5]]
sage: Ms = [Ts[0].parent().sum(weights[i][j] * Ts[j] for j in range(len(Ts))) for i in range(len(weights))]
sage: δ = 2
sage: λ_ij = weights_from_moment_matrices_fixed_δ(Ms, d=D-δ, δ=δ, dim=dim, tol=tol, torus=True)
sage: _check_weights(weights, λ_ij, tol)
True
sage: λ_ij = weights_from_moment_matrices(Ms, d=3, dim=dim, tol=tol, torus=True)
sage: _check_weights(weights, λ_ij, tol)
True
sage: λ_ij.ncols() == len(Ts)
True
From the weights λ_ij, we recover the original Toeplitz matrices (up to permutation), after normalizing them:
sage: Ts_rec = [Tj/Tj[0,0] for Tj in unmixed_moment_matrices(Ms, λ_ij)]
sage: all(any((Tj_rec - Tl).norm() < tol for Tl in Ts) for Tj_rec in Ts_rec)
True
In the case that δ=0, we find more common eigenvalues than desired:
sage: λ_ij = weights_from_moment_matrices_fixed_δ(Ms, d=D, δ=0, dim=dim, tol=tol, torus=True)
sage: λ_ij.ncols() > len(Ts)
True
Example: affine line segments¶
sage: functionals = [moment_functionals.mk_moment_affine_line_segment(start, end) for start, end in
....: [[(2,3,1), (3,4,0)], [(1,-3,-5), (0,2,0)], [(-3,-7,2), (3,-4,-4)]]]
sage: dim, D = 3, 3
sage: Hs = [moment_functionals.hankel(f, dim=dim, d=D).change_ring(RDF) for f in functionals]
sage: weights = [[1, 1, 1], [2, 3, 6], [-1, 3, 5]]
sage: Ms = [Hs[0].parent().sum(weights[i][j] * Hs[j] for j in range(len(Hs))) for i in range(len(weights))]
Each of the three components lives in a variety of degree 1, so total degree 2 is enough to find common eigenspaces (polynomials that vanish on two of three components). An upper bound on δ is 1, so d+δ ≤ 3:
sage: δ = 1
sage: λ_ij = weights_from_moment_matrices_fixed_δ(Ms, d=D-δ, δ=δ, dim=dim, tol=tol)
sage: _check_weights(weights, λ_ij, tol)
True
sage: λ_ij = weights_from_moment_matrices(Ms, d=2, dim=dim, tol=tol)
sage: _check_weights(weights, λ_ij, tol)
True
sage: λ_ij.ncols() == len(Hs)
True
In this case, we already find the correct weights when δ=0:
sage: λ_ij = weights_from_moment_matrices_fixed_δ(Ms, d=D, δ=0, dim=dim, tol=tol)
sage: _check_weights(weights, λ_ij, tol)
True
Note that picking d larger than necessary can introduce numerical errors that requires increasing the tolerance:
sage: D = 4
sage: Hs = [moment_functionals.hankel(f, dim=dim, d=D).change_ring(RDF) for f in functionals]
sage: Ms = [Hs[0].parent().sum(weights[i][j] * Hs[j] for j in range(len(Hs))) for i in range(len(weights))]
sage: λ_ij = weights_from_moment_matrices_fixed_δ(Ms, d=3, δ=1, dim=dim, tol=tol)
sage: _check_weights(weights, λ_ij, tol)
False
sage: λ_ij = weights_from_moment_matrices_fixed_δ(Ms, d=3, δ=1, dim=dim, tol=1e2*tol)
sage: _check_weights(weights, λ_ij, 1e2*tol)
True
Example: random functionals on some trigonometric varieties¶
These functionals will not be positive-semidefinite. Recovery works nonetheless. The components are of bidegree (1,2), (2,1), (1,1). Thus, max-degree 3 is enough to find common eigenspaces, i.e. there is a polynomial of degree (3,3) vanishing on the first two components, but not on the third. As the varieties are defined by some unstructured polynomials, we expect δ=0 to be sufficient for recovery.
sage: L.<x,y> = LaurentPolynomialRing(QQ)
sage: random_fun = lambda: ZZ.random_element(1, 30)
sage: functionals = [moment_functionals.random_moment_functional(L.ideal([p]), random_fun=random_fun)
....: for p in [y^2 - x*y + 1, x^2 + y + x*y + y, x*y + 3*x + 1]]
sage: D = 3
sage: dim = L.ngens()
sage: Ts = [moment_functionals.toeplitz(f, dim=dim, d=D).change_ring(CDF) for f in functionals]
sage: weights = [[1, 1, 1], [2, 3, 6], [-1, 3, 5]]
sage: Ms = [Ts[0].parent().sum(weights[i][j] * Ts[j] for j in range(len(Ts))) for i in range(len(weights))]
sage: δ = 0
sage: λ_ij = weights_from_moment_matrices_fixed_δ(Ms, d=D-δ, δ=δ, dim=dim, tol=tol, torus=True)
sage: _check_weights(weights, λ_ij, tol)
True
sage: λ_ij = weights_from_moment_matrices(Ms, d=D, dim=dim, tol=tol, torus=True)
sage: _check_weights(weights, λ_ij, tol)
True
Example: tensor decomposition¶
We demonstrate how to compute a tensor decomposition by computing the eigenvalues of its slices, in case the rank of the tensor is at most one of the dimensions.
sage: K = QQ
sage: n = 7
sage: r = n-2
sage: s = r+3
sage: assert r <= n and r <= s
sage: x = [vector(K, [ZZ.random_element(1,30) for _ in range(s)]) for _ in range(r)]
sage: y = [vector(K, [ZZ.random_element(1,30) for _ in range(n)]) for _ in range(r)]
sage: z = [vector(K, [ZZ.random_element(1,30) for _ in range(n)]) for _ in range(r)]
sage: λ = [K(ZZ.random_element(1,30)) for _ in range(r)]
sage: import numpy as np
sage: from momentproblems.commoneigenspaces import common_eigenspaces_symbolic
sage: t = np.einsum('ja,jb,jc->abc', (matrix.diagonal(λ) * matrix(x)).numpy(int), matrix(y).numpy(int), matrix(z).numpy(int))
sage: Xs = [matrix(K, np.ascontiguousarray(t[i,:,:])) for i in range(s)]
sage: Vx, ev_x = eigenspaces_lifted(Xs, tol=None)[:2]
sage: YZs = unmixed_moment_matrices(Xs, matrix.column(K, ev_x))
sage: def check(YZ): # the rank-1 matrices are what we expect
....: u, = [w[0] for e,w,m in YZ.eigenvectors_right() if e != 0]
....: v, = [w[0] for e,w,m in YZ.eigenvectors_left() if e != 0]
....: return any(matrix([u, y[j]]).rank() == 1 == matrix([v, z[j]]).rank() == 1 for j in range(r))
sage: _check_weights(matrix.column(K, x), matrix.column(K, ev_x), tol=None)
True
sage: all(check(YZ) for YZ in YZs)
True
- mainalgorithm.shift_invariant_eigenvalue_indices(Ur, eigspaces, d, δ, dim, tol, basis_fun, deg_fun)¶
- mainalgorithm.eigenspaces_lifted(Ms, M=None, *, tol, convex=False, **kwds)¶
- mainalgorithm.eigenvalues_from_indices(eigvals, indices, *, exact=False)¶
- mainalgorithm.convex_hull_vertices(eigvals, *, tol, exact)¶
- mainalgorithm.weights_from_moment_matrices_fixed_δ(Ms, d, δ, dim, tol, torus=False, **kwds)¶
- mainalgorithm.weights_from_moment_matrices_convex(Ms, M=None, *, tol, **kwds)¶
- mainalgorithm.eigenvalues_dimension(eigvals, *, exact, tol)¶
Return the dimension of the projective subspace spanned by the eigenvalues
- mainalgorithm.weights_from_moment_matrices(Ms0, d, dim, tol, torus=False, δ_MAX=10, **kwds)¶
- mainalgorithm.unmixed_moment_matrices(Ms, λ_ij)¶
