🧮 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 that is 8 paths to track instead of 256.
This is the payoff of Bertini’s multihomogeneous start system: it tracks one path per actual solution, not one per total-degree phantom.
Setting up the system¶
We need 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('lam')
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('lam')
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('lam')
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()
good = solver.finite_solutions() # the n eigenpairs (lam is the last coord)
assert len(good) == n # one path per eigenvalue
return sorted(complex(s[len(s) - 1]).real for s in good)
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.
The same idea extends to matrices of variables (a NumPy object array of Variable), which is
where problems like rank conditions and regeneration are headed. Today the vector
equations are expanded into scalar polynomials; a future evaluation block will carry the
matrix structure all the way into the numerics for larger problems.
Complete example¶
The whole tutorial as one runnable script – assemble nothing, just run it:
"""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('lam')
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('lam')
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()
good = solver.finite_solutions() # the n eigenpairs (lam is the last coord)
assert len(good) == n # one path per eigenvalue
return sorted(complex(s[len(s) - 1]).real for s in good)
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()