A numerically stable, BLAS/LAPACK-compatible $O(m^2)$ QR factorization for Cauchy-like matrices and its integration into Fastgcd for approximate polynomial GCD

Fastgcd, a leading algorithm for the approximate greatest common divisor (GCD) of two univariate polynomials with inexact coefficients, attains its low $O(m^2)$ per-iteration cost by factoring displacement-structured Cauchy-like matrices (obtained from the Toeplitz-like Sylvester, subresultant and Jacobian matrices through an FFT-based conversion) with a fast LU factorization known as the modified GKO method. This LU route does not preserve row-wise norms, and in benchmarks it loses to the highly optimized, though asymptotically cubic, Householder QR routines of BLAS/LAPACK. Nagasaka's 2021 comparative study of approximate-GCD algorithms therefore posed an explicit open problem: devise a numerically stable QR factorization of complexity $O(m^2)$ for Cauchy-like matrices converted from Toeplitz-like matrices, admitting a BLAS/LAPACK-compatible (i.e. practically fast) implementation, whose substitution for the fast LU inside Fastgcd would improve the algorithm in both theory and practice. The fast structured-QR methods available at the time (Lanczos-type and rank-structured approaches of 2003-2010) had been tried in exactly this setting and found numerically unstable. No published algorithm is known to meet the three joint requirements of quadratic complexity, demonstrated numerical stability for Cauchy-like matrices arising from Toeplitz-like ones, and practical BLAS/LAPACK-level speed. The problem remains open.

Background

In symbolic-numeric computation, the approximate GCD problem asks: given univariate polynomials $f$ and $g$ whose coefficients carry a priori errors, find perturbations $\Delta f, \Delta g$ of minimal norm such that $f+\Delta f$ and $g+\Delta g$ have a nontrivial common divisor, and compute that divisor. The problem goes back to numerical variants of the Euclidean algorithm (Schönhage's quasi-GCD, Dunaway, Noda-Sasaki) and received its modern formulation through Emiris-Galligo-Lombardi; it underlies applications from blind image deconvolution to control theory.

Nearly all leading algorithms reduce the problem to matrix computations on the Sylvester matrix $\mathrm{Syl}(f,g)$, the coefficient matrix of the polynomial identity $f\,u = g\,v$ (a matrix whose submatrices, the subresultant matrices, become rank-deficient exactly when a common factor of a given degree exists). QRGCD applies pivoted QR factorizations to these matrices; UVGCD combines QR factorizations with a Gauss-Newton refinement; structured total least norm (STLN) methods and GPGCD treat the task as structured low-rank approximation or optimization; Fastgcd, introduced by Bini and Boito, exploits matrix structure to reduce the asymptotic cost. Nagasaka's 2021 study disassembled these algorithms into interchangeable subroutines (degree search, candidate factors, refinement) and benchmarked all reassembled combinations under uniform implementations, concluding that UVGCD gives the smallest perturbations while Fastgcd has the best theoretical complexity but disappointing practical speed.

The structure Fastgcd exploits is displacement structure. A matrix $S$ has displacement rank at most $r$ if it satisfies $O_L S - S O_R = V H$ for thin matrices $V \in \mathbb{C}^{m\times r}$ and $H \in \mathbb{C}^{r\times n}$ called generators, with fixed displacement operators $O_L, O_R$. When the operators are shift matrices, $S$ is Toeplitz-like; the Sylvester, subresultant and Jacobian matrices of a polynomial pair have displacement rank 2, 2 and 5, so they are compressible to generator form. When the operators are diagonal, $S$ is Cauchy-like, and any Toeplitz-like matrix can be converted to Cauchy-like form $C = F_m S D^{-1} F_n^{H}$ (with $F$ denoting Fourier matrices) by transforming only the generators with FFTs. For Cauchy-like matrices there is a fast LU factorization with pivoting operating on generators, the GKO method of Gohberg-Kailath-Olshevsky in the numerically more robust variant of Gu, of complexity $O(m^2)$; this is the modified GKO method that Fastgcd uses for degree detection and for the least-squares solves in its Gauss-Newton refinement.

The difficulty is stability and constants. Householder QR is norm-wise backward stable and preserves row-wise norms, which is why QR-based subroutines calling the highly optimized kernels of BLAS/LAPACK (the standard libraries of tuned dense linear-algebra operations; ATLAS, OpenBLAS and MKL are widely used optimized distributions) dominate in practice despite $O(m^3)$ complexity, while the fast LU does not preserve row-wise norms and suffers in comparison. Fast QR factorizations of Cauchy-like matrices do exist: a Lanczos-type algorithm of Fasino and Gemignani, a QR-based solver for rank structured matrices of Delvaux and Van Barel, the fast QR of Gemignani-Van Barel-Delvaux, and the rank-structured approach of Delvaux-Gemignani-Van Barel, which converts the structured QR into operations on quasiseparable (rank-structured) representations. However, when Nagasaka implemented the last of these for Cauchy-like matrices converted from Toeplitz-like matrices, the result was numerically unstable, and none of these methods is organized around the blocked, cache-friendly kernels that make BLAS/LAPACK fast. A further constraint is that factors of the converted Cauchy-like matrix cannot be transformed back into factors of the original Toeplitz-like matrix (the FFT destroys triangularity), so any fast QR must be carried out and used entirely in generator form for rank detection, null-space vectors, and least-squares solves. This is the gap the source paper identifies: a stable, practical, $O(m^2)$ QR factorization for Cauchy-like matrices is missing, and its absence is what prevents Fastgcd from combining its asymptotic advantage with QR-level robustness and speed.

Problem Statement

As posed in the concluding section of Nagasaka (2021), the problem is to devise a QR factorization algorithm for Cauchy-like matrices - specified by their displacement generators, in particular the Cauchy-like matrices obtained by converting the Toeplitz-like Sylvester, subresultant and Jacobian matrices that arise in approximate-GCD computations - that simultaneously satisfies the following joint requirements:

1. Complexity. The factorization (or exactly the information Fastgcd draws from it: numerical rank for degree detection, least-squares solutions, and null-space vectors for candidate and refined factors) is computed in $O(m^2)$ floating-point operations for $m \times n$ matrices with $m \approx n$ and small fixed displacement rank, operating only on generators (FFT-based generator conversions with $O(rm\log m)$ cost are admissible; forming dense factors is not the intent, since factors of the converted matrix cannot be transformed back into factors of the original Toeplitz-like matrix).

2. Numerical stability. The algorithm is numerically stable, at minimum with backward-error behavior comparable to dense Householder QR, and this is demonstrably true for the class that matters here: Cauchy-like matrices converted from Toeplitz-like matrices, the class in which the rank-structured method of Delvaux-Gemignani-Van Barel (2010) was found to be unstable.

3. Practicality. The algorithm admits a BLAS/LAPACK-compatible implementation, meaning that its dominant kernels are the highly optimized routines of those libraries (or comparably vectorized, cache-friendly operations), so that the asymptotic $O(m^2)$ advantage is realized as a wall-clock advantage over the $O(m^3)$ dense LAPACK QR at the polynomial degrees occurring in practice.

In addition, the source states as its motivation that replacing the LU factorization with such a QR factorization may make Fastgcd much better: substituting the QR factorization for the modified GKO fast LU factorization in Fastgcd's degree-search and refinement steps is expected to improve Fastgcd's practical accuracy and/or speed. This expected consequence is the motivation for the problem, not an additional requirement on top of the three joint requirements above.

References

1. Kosaku Nagasaka, "Toward the best algorithm for approximate GCD of univariate polynomials", Journal of Symbolic Computation 105 (2021) 4-27. DOI: 10.1016/j.jsc.2019.08.004. URL: https://doi.org/10.1016/j.jsc.2019.08.004 (manuscript: https://da.lib.kobe-u.ac.jp/da/kernel/90007781/90007781.pdf)
2. Dario A. Bini, Paola Boito, "Structured matrix-based methods for polynomial ε-gcd: analysis and comparisons", Proceedings of the 2007 International Symposium on Symbolic and Algebraic Computation (ISSAC 2007), ACM (2007) 9-16. DOI: 10.1145/1277548.1277551. URL: https://doi.org/10.1145/1277548.1277551
3. Zhonggang Zeng, "The numerical greatest common divisor of univariate polynomials", in Randomization, Relaxation, and Complexity in Polynomial Equation Solving, Contemporary Mathematics 556, American Mathematical Society (2011) 187-217. DOI: 10.1090/conm/556/11014. URL: https://doi.org/10.1090/conm/556/11014
4. R. M. Corless, S. M. Watt, L. Zhi, "QR Factoring to Compute the GCD of Univariate Approximate Polynomials", IEEE Transactions on Signal Processing 52 (2004) 3394-3402. DOI: 10.1109/TSP.2004.837413. URL: https://doi.org/10.1109/TSP.2004.837413
5. Akira Terui, "GPGCD: an iterative method for calculating approximate GCD of univariate polynomials", Theoretical Computer Science 479 (2013) 127-149. DOI: 10.1016/j.tcs.2012.10.023. URL: https://doi.org/10.1016/j.tcs.2012.10.023
6. Mark Giesbrecht, Joseph Haraldson, George Labahn, "Computing the nearest rank-deficient matrix polynomial", Proceedings of the 2017 ACM International Symposium on Symbolic and Algebraic Computation (ISSAC 2017), ACM (2017) 181-188. DOI: 10.1145/3087604.3087648. URL: https://doi.org/10.1145/3087604.3087648
7. Ming Gu, "Stable and efficient algorithms for structured systems of linear equations", SIAM Journal on Matrix Analysis and Applications 19 (1998) 279-306. DOI: 10.1137/S0895479895291273. URL: https://doi.org/10.1137/S0895479895291273
8. Dario Fasino, Luca Gemignani, "A Lanczos-type algorithm for the QR factorization of Cauchy-like matrices", in Fast Algorithms for Structured Matrices: Theory and Applications, Contemporary Mathematics 323, American Mathematical Society (2003) 91-104. DOI: 10.1090/conm/323/05698. URL: https://doi.org/10.1090/conm/323/05698
9. Steven Delvaux, Marc Van Barel, "A QR-based solver for rank structured matrices", SIAM Journal on Matrix Analysis and Applications 30 (2008) 464-490. DOI: 10.1137/060654979. URL: https://doi.org/10.1137/060654979
10. Luca Gemignani, Marc Van Barel, Steven Delvaux, "Fast QR factorization of Cauchy-like matrices", Linear Algebra and its Applications 428 (2008) 697-711. DOI: 10.1016/j.laa.2007.08.038. URL: https://doi.org/10.1016/j.laa.2007.08.038
11. Steven Delvaux, Luca Gemignani, Marc Van Barel, "QR-factorization of displacement structured matrices using a rank structured matrix approach", in Numerical Methods for Structured Matrices and Applications, Operator Theory: Advances and Applications 199, Birkhauser (2010) 229-254. DOI: 10.1007/978-3-7643-8996-3_11. URL: https://doi.org/10.1007/978-3-7643-8996-3_11
12. Masaru Sanuki, "Challenge to fast and stable computation of approximate univariate GCD, based on displacement structures", Proceedings of the 2011 International Workshop on Symbolic-Numeric Computation (SNC 2011), ACM (2012) 178-186. DOI: 10.1145/2331684.2331713. URL: https://doi.org/10.1145/2331684.2331713
13. Konstantin Usevich, Ivan Markovsky, "Variable projection methods for approximate (greatest) common divisor computations", Theoretical Computer Science 681 (2017) 176-198. DOI: 10.1016/j.tcs.2017.03.028. URL: https://doi.org/10.1016/j.tcs.2017.03.028
14. Antonio Fazzi, Nicola Guglielmi, Ivan Markovsky, "An ODE-based method for computing the approximate greatest common divisor of polynomials", Numerical Algorithms 81 (2019) 719-740. DOI: 10.1007/s11075-018-0569-0. URL: https://doi.org/10.1007/s11075-018-0569-0
15. Nithin Govindarajan, Raphael Widdershoven, Shivkumar Chandrasekaran, Lieven De Lathauwer, "A fast algorithm for computing Macaulay null spaces of bivariate polynomial systems", SIAM Journal on Matrix Analysis and Applications 45 (2024) 368-396. DOI: 10.1137/23M1550414. URL: https://doi.org/10.1137/23M1550414
16. Jianlin Xia, Yuanzhe Xi, Ming Gu, "A superfast structured solver for Toeplitz linear systems via randomized sampling", SIAM Journal on Matrix Analysis and Applications 33 (2012) 837-858. DOI: 10.1137/110831982. URL: https://doi.org/10.1137/110831982
17. Bernhard Beckermann, Daniel Kressner, Heather Wilber, "Compression properties for large Toeplitz-like matrices", Numerical Algorithms 100 (2025) 2041-2067. DOI: 10.1007/s11075-025-02185-8. URL: https://doi.org/10.1007/s11075-025-02185-8 (preprint arXiv:2502.09823, https://arxiv.org/abs/2502.09823)
18. Sara Khichane, Vincent Neiger, "Matrices with displacement structure: a deterministic approach for linear systems and nullspace bases", arXiv:2603.02425 (2026). URL: https://arxiv.org/abs/2603.02425
19. Kosaku Nagasaka, "Approximate GCD of several multivariate sparse polynomials based on SLRA interpolation", Journal of Symbolic Computation 127 (2025) 102368. DOI: 10.1016/j.jsc.2024.102368. URL: https://doi.org/10.1016/j.jsc.2024.102368