Ai-Long Zheng

New progress in real and complex polynomial root-finding

Matrix methods are increasingly popular for polynomial root-finding. The idea is to approximate the roots as the eigenvalues of the companion or generalized companion matrix associated with an input polynomial. The algorithms also solve secular equation. QR algorithm is the most customary method for eigen-solving, but we explore the inverse Rayleigh quotient iteration instead, which turns out to be competitive with the most popular root-finders because of its excellence in exploiting matrix structure. To advance the iteration we preprocess the matrix and incorporate Newtons linearization, repeated squaring, homotopy continuation techniques, and some heuristics. The resulting algorithms accelerate the known numerical root-finders for univariate polynomial and secular equations, and are particularly well suited for the acceleration by using parallel processing. Furthermore, even on serial computers the acceleration is dramatic for numerical approximation of the real roots in the typical case where they are much less numerous than all complex roots.

Superfast algorithms for Cauchy-like matrix computations and extensions

Abstract An effective algorithm of [M. Morf, Ph.D. Thesis, Department of Electrical Engineering, Stanford University, Stanford, CA, 1974; in: Proceedings of the IEEE International Conference on ASSP, IEEE Computer Society Press, Silver Spring, MD, 1980, pp. 954–959; R.R. Bitmead and B.D.O. Anderson, Linear Algebra Appl. 34 (1980) 103–116] computes the solution x → =T −1 b → to a strongly nonsingular Toeplitz or Toeplitz-like linear system T x → = b → , a short displacement generator for the inverse T −1 of T , and det T . We extend this algorithm to the similar computations with n×n Cauchy and Cauchy-like matrices. Recursive triangular factorization of such a matrix can be computed by our algorithm at the cost of executing O (nr 2 log 3 n) arithmetic operations, where r is the scaling rank of the input Cauchy-like matrix C (r=1 if C is a Cauchy matrix). Consequently, the same cost bound applies to the computation of the determinant of C , a short scaling generator of C −1 , and the solution to a nonsingular linear system of n equations with such a matrix C . (Our algorithm does not use the reduction to Toeplitz-like computations.) We also relax the assumptions of strong nonsingularity and even nonsingularity of the input not only for the computations in the field of complex or real numbers, but even, where the algorithm runs in an arbitrary field. We achieve this by using randomization, and we also show a certain improvement of the respective algorithm by Kaltofen for Toeplitz-like computations in an arbitrary field. Our subject has close correlation to rational tangential (matrix) interpolation under passivity condition (e.g., to Nevanlinna–Pick tangential interpolation problems) and has further impact on the decoding of algebraic codes.

Newton's Iteration for Inversion of Cauchy-Like and Other Structured Matrices

We specify some initial assumptions that guarantee rapid refinement of a rough initial approximation to the inverse of a Cauchy-like matrix, by means of our new modification of Newtons iteration, where the input, output, and all the auxiliary matrices are represented with their short generators defined by the associated scaling operators. The computations are performed fast since they are confined to operations with short generators of the given and computed matrices. Because of the known correlations among various structured matrices, the algorithm is immediately extended to rapid refinement of rough initial approximations to the inverses of Vandermonde-like, Chebyshev?Vandermonde-like, and Toeplitz-like matrices, where again the computations are confined to operations with short generators of the involved matrices.

Randomized preconditioning of the MBA algorithm

MBA algorithm inverts a structured matrix in nearly linear arithmetic time but requires a serious restriction on the input class. We remove this restriction by means of randomization and extend the progress to some fundamental computations with polynomials, e.g., computing their GCDs and AGCDs, where most effective known algorithms rely on computations with matrices having Toeplitz-like structure. Furthermore, our randomized algorithms fix rank deficiency and ill conditioning of general and structured matrices. At the end we comment on a wide range of other natural extensions of our progress and underlying ideas.

Root-finding by expansion with independent constraints

Elimination methods are highly effective for the solution of linear and nonlinear systems of equations, but reversal of the elimination principle can be beneficial as well: competent incorporation of additional independent constraints and variables or more generally immersion of the original computational problem into a larger task, defined by a larger number of independent constraints and variables can improve global convergence of iterative algorithms, that is their convergence from the start. A well known example is the dual linear and nonlinear programming, which enhances the power of optimization algorithms. We believe that this is just an ad hoc application of general Principle of Expansion with Independent Constraints; it should be explored systematically for devising iterative algorithms for the solution of equations and systems of equations and for optimization. At the end of this paper we comment on other applications and extensions of this principle. Presently we show it at work for the approximation of a single zero of a univariate polynomial p of a degree n. Empirical global convergence of the known algorithms for this task is much weaker than that of the algorithms for all n zeros, such as Weierstrass-Durand-Kerners root-finder, which reduces its root-finding task to Vietes (Vietas) system of n polynomial equations with n unknowns. We adjust this root-finder to the approximation of a single zero of p, preserve its fast global convergence and decrease the number of arithmetic operations per iteration from quadratic to linear. Together with computing a zero of a polynomial p, the algorithm deflates this polynomial as by-product, and then could be reapplied to the quotient to approximate the next zero of p. Alternatively by using m processors that exchange no data, one can concurrently approximate up to m zeros of p. Our tests confirm the efficiency of the proposed algorithms. Technically our root-finding boils down to computations with structured matrices, polynomials and partial fraction decompositions. Our study of these links can be of independent interest; e.g., as by-product we express the inverse of a Sylvester matrix via its last column, thus extending the celebrated result of Gohberg and Sementsul (1972) [22] from Toeplitz to Sylvester matrix inverses.

On Isolation of Real and Nearly Real Zeros of a Univariate Polynomial and Its Splitting into Factors

We show two simple algorithms for isolation of the real and nearly real zeros of a univariate polynomial, as well as of those zeros that lie on or near a fixed circle on the complex plane. We also simplify slightly approximation of complex zeros of a polynomial with real coefficients.

Advancing matrix computations with randomized preprocessing

The known algorithms for linear systems of equations perform significantly slower where the input matrix is ill conditioned, that is lies near a matrix of a smaller rank. The known methods counter this problem only for some important but special input classes, but our novel randomized augmentation techniques serve as a remedy for a typical ill conditioned input and similarly facilitates computations with rank deficient input matrices. The resulting acceleration is dramatic, both in terms of the proved bit-operation cost bounds and the actual CPU time observed in our tests. Our methods can be effectively applied to various other fundamental matrix and polynomial computations as well.

Eigen-solving via reduction to DPR1 matrices

Highly effective polynomial root-finders have been recently designed based on eigen-solving for DPR1 (that is diagonal + rank-one) matrices. We extend these algorithms to eigen-solving for the general matrix by reducing the problem to the case of the DPR1 input via intermediate transition to a TPR1 (that is triangular + rank-one) matrix. Our transforms use substantially fewer arithmetic operations than the QR classical algorithms but employ non-unitary similarity transforms of a TPR1 matrix, whose representation tends to be numerically unstable. We, however, operate with TPR1 matrices implicitly, as with the inverses of Hessenberg matrices. In this way our transform of an input matrix into a similar DPR1 matrix partly avoids numerical stability problems and still substantially decreases arithmetic cost versus the QR algorithm.

Real and complex polynomial root-finding with eigen-solving and preprocessing

Recent progress on root-finding for polynomial and secular equations largely relied on eigen-solving for the associated companion and diagonal plus rank-one generalized companion matrices. By applying to them Rayleigh quotient iteration, we could have already competed with the current best polynomial root-finders, but we achieve further speedup by applying additive preprocessing. Moreover our novel rational maps of the input matrix enables us to direct the iteration to approximating only real roots, so that we dramatically accelerate their numerical computation in the important case where they are much less numerous than all complex roots.

Preconditioning, randomization, solving linear systems, eigen-solving, and root-finding

We propose novel randomized preprocessing techniques for solving linear systems of equations and eigen-solving with extensions to the solution of polynomial and secular equations. According to our formal study and extensive experiments, the approach turns out to be effective, particularly in the case of structured input matrices.