๐งฎ Eigenvalues by homotopy continuationยถ
The eigenvalue problem \(A x = \lambda x\) is the cleanest place to see why numerical algebraic geometry cares about variable-group structure. Written as a polynomial system it is
a system that is degree 1 in the eigenvector coordinates \(x\) and degree 1 in the eigenvalue \(\lambda\). Treating \(x\) and \(\lambda\) as separate groups of variables, every equation has bidegree \((1,1)\), and the multihomogeneous Bรฉzout number is \(n\) โ exactly the number of eigenvalues. The naive total-degree count would be \(2^n\). For an \(8\times 8\) matrix, using a two-homogeneous variable group structure yields only 8 paths to track, instead of 256 if we use only a single variable group. This is the payoff of Bertiniโs multihomogeneous start system for this system: it tracks one path per actual solution.
Setting up the systemยถ
Weโll use numpy for the matrix and the cross-check, and bertini for the solve and
for writing the equations as actual linear algebra:
import numpy as np
import bertini as bertini
from bertini import ZeroDimSolver
Pick a small symmetric matrix (real, distinct eigenvalues make the check easy to read):
A = np.array([[2, 1, 0],
[1, 3, 1],
[0, 1, 4]])
n = A.shape[0]
With the linear-algebra layer the equations are linear algebra: make a vector of
variables for the eigenvector and a scalar for the eigenvalue, and \((A-\lambda I)x\)
is just A @ x - lam*x:
x = np.array(bertini.variables('x', n), dtype=object) # array([x0, x1, x2], dtype=object)
lam = bertini.Variable('๐')
Note
The matrix A here is integer, so it enters the function tree exactly.
bertini deliberately refuses plain Python float entries โ a 64-bit float
would cap the precision of the whole arbitrary-precision computation at ~16 digits. For
non-integer coefficients pass exact values
(bertini.coefficients([['5/2', '1'], ...]) accepts fractions, exact decimal
strings, and bertini multiprecision values).
Two ways to make it zero-dimensionalยถ
An eigenvector is only defined up to scale, so \((A-\lambda I)x = 0\) alone has a whole line of solutions for each eigenvalue (and the trivial \(x=0\)). There are two clean ways to cut that down to isolated points, and Bertini supports both.
Affine, with a normalization. Keep \(x\) an ordinary (affine) variable group and add one generic linear equation \(c\cdot x = 1\) to pin the scale:
def build_affine():
x = np.array(bertini.variables('x', n), dtype=object)
lam = bertini.Variable('๐')
c = np.array([5, 8, 3]) # any generic integer vector
sys = bertini.System()
sys.add_functions(A @ x - lam * x) # the rows of (A - lam I) x
sys.add_function(c @ x - 1) # fix the eigenvector scale
sys.add_variable_group(bertini.VariableGroup(list(x)))
sys.add_variable_group(bertini.VariableGroup([lam]))
return sys
Projective. Declare \(x\) a projective (homogeneous) variable group: the eigenvector then lives in \(\mathbb{P}^{n-1}\) natively, where scale is already quotiented out, so no normalization equation is needed:
def build_projective():
x = np.array(bertini.variables('x', n), dtype=object)
lam = bertini.Variable('๐')
sys = bertini.System()
sys.add_functions(A @ x - lam * x) # nothing else!
sys.add_hom_variable_group(bertini.VariableGroup(list(x))) # x in P^{n-1}
sys.add_variable_group(bertini.VariableGroup([lam]))
return sys
Both are multihomogeneous with the same Bรฉzout number \(n\), so both track exactly \(n\) paths โ the projective version just expresses the geometry directly instead of bolting on a normalization.
Solving, and reading off the eigenvaluesยถ
The solve is identical for either system. A small helper runs it and pulls out the
eigenvalues (lam is the last user coordinate of each solution; the rest are the
eigenvector):
def eigenvalues_of(system):
solver = ZeroDimSolver(system, mptype='adaptive', startsystem='mhom')
solver.solve()
solns = solver.finite_solutions() # the n eigenpairs (lam is the last coord)
assert len(solns) == n # one path per eigenvalue
return sorted(complex(s[len(s) - 1]).real for s in solns)
Both formulations recover numpyโs eigenvalues:
expected = sorted(np.linalg.eigvals(A).real)
for build in (build_affine, build_projective):
got = eigenvalues_of(build())
for g, e in zip(got, expected):
assert abs(g - e) < 1e-6
Comparing the twoยถ
Both track the same number of paths, but the projective system is leaner โ one fewer equation, and the eigenvector group needs no homogenizing coordinate โ so it is usually the faster of the two. Time them and see:
import time
for name, build in (("affine+normalization", build_affine),
("projective ", build_projective)):
t0 = time.perf_counter()
got = eigenvalues_of(build())
dt = time.perf_counter() - t0
print(f"{name}: {got} in {dt:.3f}s")
affine+normalization: [...] in ...s
projective : [...] in ...s
The projective system is the more direct statement of the problem and typically solves faster; the affine one is handy when you would rather not think projectively. Either way the solver returns the same eigenvalues โ pick whichever reads better for your problem.
Why this matters, and where it is goingยถ
Two things are worth noticing. First, the eigenvalues came out of path tracking, not a
dedicated eigensolver โ the same machinery solves any polynomial system, and here the
multihomogeneous structure made it efficient (n paths, not \(2^n\)). Second, we
never built \((A-\lambda I)x\) entry by entry: with \(x\) a vector of variables, the
condition reads A @ x - lam*x โ the linear algebra you would write on paper. That is
what bertiniโs linear-algebra support gives you, and this problem is its motivating
example. Coefficients on variables stay exact, so the precision of the solve is never
silently capped by a stray floating-point literal.
Complete exampleยถ
"""Eigenvalues by homotopy continuation -- Bertini 2 tutorial.
Solve the eigenvalue problem (A - lam I) x = 0 as a multihomogeneous polynomial
system, recovering the eigenvalues by path tracking.
Run: python eigenvalues_by_homotopy.py
"""
import time
import numpy as np
import bertini as bertini
from bertini import ZeroDimSolver
def setup_matrix():
# small symmetric matrix: real, distinct eigenvalues make the check easy to read
A = np.array([[2, 1, 0],
[1, 3, 1],
[0, 1, 4]])
n = A.shape[0]
return A, n
def build_affine(A, n):
x = np.array(bertini.variables('x', n), dtype=object)
lam = bertini.Variable('๐')
c = np.array([5, 8, 3]) # any generic integer vector
sys = bertini.System()
sys.add_functions(A @ x - lam * x) # the rows of (A - lam I) x
sys.add_function(c @ x - 1) # fix the eigenvector scale
sys.add_variable_group(bertini.VariableGroup(list(x)))
sys.add_variable_group(bertini.VariableGroup([lam]))
return sys
def build_projective(A, n):
x = np.array(bertini.variables('x', n), dtype=object)
lam = bertini.Variable('๐')
sys = bertini.System()
sys.add_functions(A @ x - lam * x) # nothing else!
sys.add_hom_variable_group(bertini.VariableGroup(list(x))) # x in P^{n-1}
sys.add_variable_group(bertini.VariableGroup([lam]))
return sys
def eigenvalues_of(system, n):
solver = ZeroDimSolver(system, mptype='adaptive', startsystem='mhom')
solver.solve()
solns = solver.finite_solutions() # the n eigenpairs (lam is the last coord)
assert len(solns) == n # one path per eigenvalue
return sorted(complex(s[len(s) - 1]).real for s in solns)
def check_both(A, n):
# both formulations recover numpy's eigenvalues
expected = sorted(np.linalg.eigvals(A).real)
for build in (build_affine, build_projective):
got = eigenvalues_of(build(A, n), n)
for g, e in zip(got, expected):
assert abs(g - e) < 1e-6
def compare_timing(A, n):
for name, build in (("affine+normalization", build_affine),
("projective ", build_projective)):
t0 = time.perf_counter()
got = eigenvalues_of(build(A, n), n)
dt = time.perf_counter() - t0
print(f"{name}: {got} in {dt:.3f}s")
def main():
A, n = setup_matrix()
check_both(A, n)
compare_timing(A, n)
if __name__ == '__main__':
main()