Solving Multivariate Polynomial Systems and Rectangular Multiparameter Eigenvalue Problems with MacaulayLab

Abstract: We present the Matlab toolbox MacaulayLab, which implements numerical linear algebra algorithms for solving multivariate polynomial systems and rectangular multiparameter eigenvalue problems. Its structure and functionality are the result of several years of research and algorithmic development. We demonstrate how the software works and compare its performance with other software packages, such as PNLA, PHCpack, and MultiParEig. Some core features of MacaulayLab are the fact that it solves two key problems via one common approach, works independently of the chosen polynomial basis and monomial order, and is capable of dealing with positive-dimensional solution sets at infinity. The toolbox (including its future updates) and a large collection of test problems are freely available online.

• The paper presents MacaulayLab, a unified Macaulay-matrix toolbox that solves zero-dimensional multivariate polynomial systems and rectangular multiparameter eigenvalue problems through numerical linear algebra.
• The paper demonstrates modular features including basis and monomial-order independence, solutions-at-infinity deflation, shared Schur decompositions, optional clustering, and residuals based on smallest singular values.
• The paper reports a database of 290 polynomial systems and 30 RMEPs, with sparse block-wise null-space methods running about 20–27 times faster than a column-space configuration and outperforming MultiParEig on the tested RMEPs.

MacaulayLab is a MATLAB toolbox that solves two superficially distinct problems—finding the common roots of multivariate polynomial systems and computing the eigenvalues of rectangular multiparameter eigenvalue problems (RMEPs)—through a single numerical linear algebra methodology built on the (block) Macaulay matrix (2605.20884). The paper describes the toolbox's architecture, its distinguishing features (basis and monomial order independence, handling of solutions at infinity), an accompanying database of test problems, and benchmark comparisons against PNLA, PHCpack, and MultiParEig.

A unified problem formulation

The unifying observation is that both problems are special cases of a system of matrix polynomial equations $\mpoly[j](\mathbf{x})\,\mathbf{y} = \mathbf{0}$ with $\|\mathbf{y}\| = 1$. When the coefficient matrices are scalar ($k = l = 1$), this reduces to the classical multivariate polynomial root-finding problem. When $s = 1$ and $k \geq l + m - 1$, the problem becomes an RMEP: finding the points $\mathbf{x}^* \in \mathbb{C}^m$ where the polynomial matrix $\mpoly(\mathbf{x}^*)$ drops below its normal rank, together with associated eigenvectors. RMEPs arise in the higher-order Heine–Stieltjes problem, cryptanalysis (minimum rank problems), globally optimal identification of misfit–latency models, and SISO model order reduction. The toolbox is described as one of the first dedicated to rectangular (as opposed to square) multiparameter eigenvalue problems.

The block Macaulay matrix approach

Both problems are solved identically. The matrix equations are multiplied by monomials up to total degree $d$, and the coefficients are arranged into a (block) Macaulay matrix $\mt{M}(d)$. For sufficiently large $d$, the nullity of $\|\mathbf{y}\| = 1$0 equals the number of solutions of a zero-dimensional problem, and the null space has a shift-invariant structure: its basis consists of (block) multivariate Vandermonde vectors evaluated at the solutions. Selecting rows corresponding to standard monomials yields $\|\mathbf{y}\| = 1$1 generalized multiplication matrices whose joint eigenvalues are the solutions. The toolbox offers two variants:

• Null space based: a basis $\|\mathbf{y}\| = 1$2 of the null space is computed (iteratively, recursively, or via sparse constructions), and multiplication maps are formed after a column compression that deflates solutions at infinity.
• Column space based: a backward QR decomposition of the column-permuted (block) Macaulay matrix yields the blocks $\|\mathbf{y}\| = 1$3 and $\|\mathbf{y}\| = 1$4 directly, and the deflation of solutions at infinity happens implicitly—avoiding explicit null space computation.

Two structural caveats govern correctness. First, the method requires the nullity to stabilize at the number of solutions, so the solution set must be zero-dimensional; the notable exception is a positive-dimensional solution set located entirely at infinity, which the column compression can still remove. Second, multiple solutions degrade accuracy; an optional clustering step, based on the eigenvalues of a random linear shift polynomial (following Corless et al.'s reordered Schur idea), groups and refines them via geometric means.

Software structure and key features

The solver proceeds in six steps: enlarging the solution subspace, checking its rank structure (block-wise by default, which is more robust and efficient than row-wise checks), column compression, constructing and solving the shift problems, optional clustering, and residual computation. The shift problems use a random linear polynomial plus the $\|\mathbf{y}\| = 1$5 coordinate functions; a single Schur decomposition is reused across all $\|\mathbf{y}\| = 1$6 problems so that solution components align on the diagonal without matching.

Two design choices stand out. The first is independence of polynomial basis and monomial order: all shift and evaluation operations are delegated to pluggable <basis> and <order> functions, with the monomial basis and graded reverse lexicographic order as defaults, and a pre-implemented Chebyshev basis—valuable since orthogonal bases have superior numerical properties. The second is the residual computation via the smallest singular value of the evaluated (matrix) polynomial, which uniformly handles both problem types.

Database and configuration comparison

An accompanying database of 290 multivariate polynomial systems and 30 RMEPs is distributed with the toolbox, including problems in non-monomial bases and problems with positive-dimensional solution sets at infinity (e.g., cyclic5, katsura7, h2fourdisk). Benchmarking configuration choices on the noon4 system and the cube RMEP shows that the sparse, block-wise null space approach is roughly 20–27× faster than the column space approach with column-wise rank checks, and the sparse variant additionally avoids materializing the often-large Macaulay matrix. Enabling the positive-dimensional flag roughly doubles runtime due to extra rank checks.

Comparison with other solvers

Performance profiles over 50 polynomial systems show PHClab (homotopy continuation) is fastest on more than half the problems, but MacaulayLab is competitive at larger ratios $\|\mathbf{y}\| = 1$7 and solves a slightly larger subset—consistent with the authors' position that homotopy methods, while efficient, are limited to square systems and can suffer ill-conditioning. For the 30 RMEPs, MacaulayLab outperforms both MultiParEig approaches (operator determinant and randomized sketching), which require linearization; the authors concede that MultiParEig remains faster for two- or three-parameter problems with large coefficient matrices, while MacaulayLab gains the advantage as the number of spectral parameters grows.

Limitations and open questions

The paper is candid about restrictions: the core method requires zero-dimensional affine solution sets (with the infinity-only exception), accuracy degrades with multiple solutions, and rank decisions are delicate and error-prone. Planned developments include exploiting matrix sparsity, support-based matrix size reduction, projective shifts to avoid rank checks, simultaneous triangularization or tensor formulations of the shift problems, and an interactive solver mode. Whether projective shifts can robustly eliminate rank checks, and how block column approaches affect the robustness of the column space variant, remain open questions the authors explicitly defer to future research.

Conclusion

MacaulayLab consolidates several years of algorithmic development into a modular, freely available toolbox that treats polynomial system solving and rectangular multiparameter eigenvalue computation as instances of one Macaulay-matrix realization problem. Its basis- and order-independence, handling of solutions at infinity, and the accompanying problem database make it a practical complement to homotopy-continuation and square-MEP solvers, particularly for overdetermined polynomial systems and higher-dimensional RMEPs.