A preprint describes a symbolic-numeric workflow—a mix of algebraic manipulation and numerical matrix calculations—that uses joint eigenvectors, or common directions extracted from several related matrices, to find solutions of polynomial equations, including repeated solutions, and to decompose symmetric tensors. The authors illustrate it with one polynomial-system computation and one tensor computation.
The calculation behind a repeated solution
The workflow first constructs a truncated normal form, which packages the system into a finite algebraic calculation, and selects a basis. It then builds multiplication matrices, calculates the eigenvalues of a random combination, and uses joint diagonalization for simple eigenvalues or joint Schur factorization for multiple ones. A default clustering threshold of 1.e-5 groups nearby values to assign multiplicities—the number of times a solution is repeated.
At the multiple solution (0,0), the calculation also returned an inverse-system basis, an algebraic way to describe local structure around a repeated solution. One element, labeled D[3], generated the same span as the local-algebra basis B0.
In the polynomial example, the system had two equations in two variables. The joint Schur calculation produced solution columns corresponding to (0,0) and (-2,2), with reported multiplicities of 3 and 1, respectively.
A second algebraic route
The team then used a resultant/Macaulay construction, a second way to organize the polynomial relations. It produced a 6-by-10 matrix. Its truncated normal form used the same basis B as the Gröbner-basis construction and gave root accuracy of the same order, although the paper reports no numerical error value or systematic benchmark.
The tensor example
For the tensor example, a symmetric tensor was represented by a homogeneous polynomial of degree 5. Apolar duality and Hankel constructions at degrees 2 and 3 extracted a truncated normal form using the basis (1, x1, x2, x1x2); the displayed normal form was a 4-by-10 matrix followed by two multiplication matrices.
The displayed projective point columns were (1,-2,2) and (1,0,0), and the reported degree vector was k=(0,2). A least-squares solve of a Vandermonde-like system returned one constant complex weight and the quadratic expression x_2^2 + 2x_1x_2 + x_1^2 + x_0x_2.
For this worked example, the reconstructed tensor had a reported apolar-norm error of 9.342121376789858e-14, which the authors describe as about 10^-13. It is a numerical diagnostic for that example, not a formal error bound or evidence that the same result will hold for other tensors.
What the examples leave open
The evidence remains narrow. The paper reports only same-order root accuracy for the Macaulay route, not a numerical error value or a systematic benchmark. The tensor error is likewise a single-example diagnostic, so it does not establish performance across other tensors or applications.
Numerical sensitivity is another open question. The paper specifies a default clustering threshold of 1.e-5, but reports no test of how assignments or joint factorizations change with the threshold, numerical precision or the random combination used in the eigensolve.
The document is an arXiv v1 preprint dated 20 Aug 2026. It identifies AlgebraicSolvers.jl as the package used for the illustrated techniques and points readers to its documentation.
Paper data and sources
Original title: Eigensolvers for polynomial roots and tensor decomposition
Authors: Enrica Barrilli, Bernard Mourrain
Journal/Repository: arXiv
Status: Preprint, not yet peer-reviewed
First online: 2026-08-20
DOI: Not available
Original paper · Full text