commoneigenspaces

commoneigenspaces.common_eigenspaces_symbolic(As, *, extend=False, max_tries=30, multiplicities=False, B=None)

TESTS:

sage: from momentproblems.commoneigenspaces import common_eigenspaces_symbolic
sage: K = GF(101)
sage: As = [
....:     matrix.diagonal(K, [2,2,3,3,4]),
....:     matrix.diagonal(K, [3,2,3,3,5]),
....:     matrix.identity(K, 5),
....:     matrix.diagonal(K, [5,2,4,4,1])]
sage: eigvals, eigmat = common_eigenspaces_symbolic(As, extend=False)
sage: eigmat.dimensions()
(5, 5)
sage: Γ = sorted(eigvals); Γ
[(1 : 1 : 51 : 1),
 (4 : 5 : 1 : 1),
 (26 : 26 : 76 : 1),
 (26 : 26 : 76 : 1),
 (61 : 41 : 81 : 1)]
sage: assert all([all(γ[j] * As[l] * v == γ[l] * As[j] * v
....:                 for j in range(1,len(γ)) for l in range(j))
....:             for γ, v in zip(eigvals, eigmat.columns())])

sage: P, Q = [matrix.random(K, As[0].ncols(), algorithm='unimodular') for _ in range(2)]
sage: Bs = [P * A * Q for A in As]
sage: eigvals2, eigmat2 = common_eigenspaces_symbolic(Bs, extend=False)
sage: sorted(eigvals2) == Γ
True
commoneigenspaces.common_eigenspaces_numeric(Δs, *, tol, B=None, algorithm='fast')

TESTS:

Generate input:

sage: from momentproblems import *
sage: from momentproblems.commoneigenspaces import common_eigenspaces_numeric
sage: tol = 1e-7
sage: functionals = [moment_functionals.mk_moment_ellipse(center=xy, radius=(1,1)) for xy in [(0,0), (-2,1), (2,-1)]]
sage: weights = [[4, 2, 1], [2, 3, 6], [1, 3, 5], [2, 1, 2]]
sage: Hs = [moment_functionals.hankel(f, dim=2, 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))]
sage: ker_M = intersections.null_space_intersection(*Ms, tol=tol)
sage: Vr = intersections.null_space(ker_M.H, tol=tol)
sage: VrH = Vr.H
sage: Δs = [VrH * Mi * Vr for Mi in Ms]  # now, Δ_i is a regular pencil

This matrix is non-singular, as it is a positive linear combination of positive-semidefinite matrices forming a regular pencil:

sage: Δ = sum(Δs)

Check that we have computed eigenvectors and eigenvalues with coordinates choosen with respect to Δ:

sage: eigspaces = list(common_eigenspaces_numeric(Δs, tol=tol, B=Δ, algorithm='fast'))
sage: len(eigspaces)
4
sage: for Q, Z, ee in eigspaces:
....:     ee = matrix(ee)
....:     for i in range(len(Δs)):
....:         for k in range(ee.ncols()):
....:             assert ((Δs[i] - ee[i,k] * Δ) * Z[:,k]).norm() < tol

sage: eigspaces = list(common_eigenspaces_numeric(Δs, tol=tol, B=Δ, algorithm='cluster'))
sage: len(eigspaces)
4
sage: for Q, Z, ee in eigspaces:
....:     ee = matrix(ee)
....:     for i in range(len(Δs)):
....:         for k in range(ee.ncols()):
....:             assert ((Δs[i] - ee[i,k] * Δ) * Z[:,k]).norm() < tol

With unspecified choice of coordinates:

sage: eigspaces = list(common_eigenspaces_numeric(Δs, tol=tol))
sage: len(eigspaces)
4
sage: for Q, Z, ee in eigspaces:
....:     ee = matrix(ee)
....:     for i in range(len(Δs)):
....:       for j in range(i):
....:         for k in range(ee.ncols()):
....:             scale = 0.1 / (abs(ee[j,k]) + abs(ee[i,k]))  # TODO an attempt to rescue the accuracy
....:             assert (scale * (ee[j,k] * Δs[i] - ee[i,k] * Δs[j]) * Z[:,k]).norm() < tol
sage: As = [
....:     matrix.diagonal(RDF, [2,2,3,3,4], sparse=False),
....:     matrix.diagonal(RDF, [3,2,3,3,5], sparse=False),
....:     matrix.diagonal(RDF, [1,1,1,1,1], sparse=False),
....:     matrix.diagonal(RDF, [5,2,4,4,1], sparse=False)]
sage: eigspaces = list(common_eigenspaces_numeric(As, B=matrix.identity(RDF, 5), tol=1e-13, algorithm='cluster'))
sage: list(map(matrix.column, sorted([matrix(ee).columns() for Q,Z,ee in eigspaces])))
[
[2.0]  [2.0]  [3.0 3.0]  [4.0]
[2.0]  [3.0]  [3.0 3.0]  [5.0]
[1.0]  [1.0]  [1.0 1.0]  [1.0]
[2.0], [5.0], [4.0 4.0], [1.0]
]