Computing nearby non-trivial Smith forms
Introduction
For a given square matrix polynomial , one can find unimodular matrices such that is a diagonal matrix . Unimodular means that there is a polynomial inverse matrix, or equivalently, that the determinant is a nonzero constant from . The unimodular matrices encapsulate the polynomial row and column operations, respectively, needed for such a diagonalization. The best known diagonalization is the Smith Normal Form (SNF, or simply Smith form) of a matrix polynomial. Here where are monic and for . The Smith form always exists and is unique though the associated unimodular matrices , are not unique (see, e.g., (Kailath, 1980; Gohberg et al., 2009)). The diagonal entries are referred to as the invariant factors of .
Matrix polynomials and their Smith forms are used in many areas of computational algebra, control systems theory, differential equations and mechanics. The Smith form is important as it effectively reveals the structure of the polynomial lattice of rows and columns, as well as the effects of localizing at individual eigenvalues. That is, it characterizes how the rank decreases as the variable t is set to different values (especially eigenvalues, where the rank drops). The Smith form is closely related to the more general Smith-McMillan form for matrices of rational functions, a form that reveals the structure of the eigenvalue at infinity, or the infinite spectral structure. The eigenvalue at infinity is non-trivial if the leading coefficient matrix is rank deficient or equivalently, the determinant does not achieve the generic degree.
The algebra of matrix polynomials is typically described assuming that the coefficients are fixed and come from an exact arithmetic domain, usually the field of real or complex numbers. In this exact setting, computing the Smith form has been well studied, and very efficient procedures are available (see (Kaltofen and Storjohann, 2015) and the references therein). However, in some applications, particularly control theory and mechanics, the coefficients can come from measured data or contain some amount of uncertainty. Compounding this, for efficiency reasons such computations are usually performed using floating point to approximate computations in , introducing roundoff error. As such, algorithms must accommodate numerical inaccuracies and are prone to numerical instability.
Numerical methods to compute the Smith form of a matrix polynomial typically rely on linearization and orthogonal transformations (Van Dooren and Dewilde, 1983; Beelen and Van Dooren, 1988; Demmel and Kågström, 1993a, Demmel and Kågström, 1993b; Demmel and Edelman, 1995) to infer the Smith form of a nearby matrix polynomial via the Jordan blocks in the Kronecker canonical form (see (Kailath, 1980)). These linearization techniques are numerically backwards stable, and for many problems this is sufficient to ensure that the computed solutions are computationally useful when a problem is continuous.
Unfortunately, the eigenvalues of a matrix polynomial are not necessarily continuous functions of the coefficients of the matrix polynomial, and backwards stability is not always sufficient to ensure computed solutions are useful in the presence of discontinuities. Previous methods are also unstructured in the sense that the computed non-trivial Smith form may not be the Smith form of a matrix polynomial with a prescribed coefficient structure, for example, maintaining the degree of entries or not introducing additional non-zero entries or coefficients. In extreme instances, the unstructured backwards error can be arbitrarily small, while the structured distance to an interesting Smith form is relatively large. Finally, existing numerical methods can also fail to compute meaningful results on some problems due to uncertainties. Examples of such problems include nearly rank deficient matrix polynomials, repeated eigenvalues or eigenvalues that are close together and other ill-posed instances.
In this paper we consider the problem of computing a nearby matrix polynomial with a prescribed spectral structure, broadly speaking, the degrees and multiplicities of the invariant factors in the Smith form, or equivalently the structure and multiplicity of the eigenvalues of the matrix polynomial. The results presented in this paper extend those in the conference paper (Giesbrecht et al., 2018). This work is not so much about computing the Smith normal form of a matrix polynomial using floating point arithmetic, but rather our focus is on the computation of a nearby matrix polynomial with “an interesting” or non-generic Smith normal form. The emphasis in this work is on the finite spectral structure of a matrix polynomial, since the techniques described are easily generalized to handle the instance of the infinite spectral structure as a special case.
The Smith form of a matrix polynomial is not continuous with respect to the usual topology of the coefficients and the resulting limitations of backward stability is not the only issue that needs to be addressed when finding nearest objects in an approximate arithmetic environment. A second issue can be illustrated by recalling the well-known representation of the invariant factors of a matrix as ratios of the determinantal divisors , where In the case of matrix polynomials, computing the nearest non-trivial Smith form is thus equivalent to finding the nearest matrix polynomial whose polynomial entries have a non-trivial GCD. Recall that approximate GCD problems can have infima that are unattainable. That is, there are co-prime polynomials with nearby polynomials with a non-trivial GCD at distances arbitrarily approaching an infimum, while at the infimum itself the GCD is trivial (see, e.g., (Giesbrecht et al., 2019a)).
The issue of unattainable infima extends to the Smith normal form. As an example, consider If we look for nearby polynomials of degree at most 2 such that (i.e. a nontrivial gcd) at minimal distance for some , then it is shown in (Haraldson, 2015, Example 3.3.6) that this distance is . This distance has an infimum of 2 at . However at we have even though for all . For to have a non-trivial Smith form we must perturb such that they have a non-trivial GCD, and thus any such perturbation must be at a distance of at least 2. However, the perturbation of distance precisely 2 has a trivial Smith form. There is no merit to perturbing the off-diagonal entries of .
Our work indirectly involves measuring the sensitivity to the eigenvalues of and the determinant of . Thus we differ from most sensitivity and perturbation analysis (e.g., (Stewart, 1994; Ahmad and Alam, 2009)), since we also study how perturbations affect the invariant factors, instead of the roots of the determinant. Additionally our theory is able to support the instance of being rank deficient and having degree exceeding one. One may also approach the problem geometrically in the context of manifolds (Edelman et al., 1997, Edelman et al., 1999). We do not consider the manifold approach directly since it does not yield numerical algorithms.
Determining what it means for a matrix polynomial to have a non-trivial Smith form numerically and finding the distance from one matrix polynomial to another matrix polynomial having an interesting or non-trivial Smith form are formulated as finding solutions to continuous optimization problems. The main contributions of this paper are deciding when has an interesting Smith form, providing bounds on a “radius of triviality” around and a structured stability analysis on iterative methods to compute a structured matrix polynomial with desired spectral properties.
The remainder of the paper is organized as follows. In Section 2 we give the notation and terminology along with some needed background used in our work. Section 3 discusses the approximate Smith form computation as an optimization problem and provides some new bounds on the distance to non-triviality. We present an optimization algorithm in Section 4 with local stability properties and rapid local convergence to compute a nearby matrix polynomial with a non-trivial Smith form and discuss implementation details. A method to compute a matrix polynomial with a prescribed lower bound on the number of ones in the Smith form is discussed in Section 5. We discuss our implementation and include some examples in Section 6. The paper ends with a conclusion along with topics for future research.
A preliminary version of some of the results in this paper appears in the proceedings of the ISSAC 2018 conference (Giesbrecht et al., 2018).
Section snippets
Preliminaries
In this section we give the notation and formal definitions needed to precisely describe the problems summarized above. We also present some existing results used as building blocks for our work. In addition, we provide a basic description of matrix functions and their first-order derivatives (Jacobian matrices) which will be needed for the optimization work central to our results.
When does a numerical matrix polynomial have a trivial SNF?
In this section we consider the question of determining if a matrix polynomial has a non-trivial SNF, or rather how much do the coefficients need to be perturbed to have a non-trivial SNF. We provide a lower bound on this distance by analyzing the distance to a reduced-rank generalized Sylvester matrix.
Approximate SNF via optimization
In this section we formulate the approximate Smith form problem as the solution to a continuous constrained optimization problem. We assume that the solutions in question are attainable and develop a method with rapid local convergence. As the problem is non-convex, our convergence analysis will be local.
Lower McCoy rank approximation
In this section we describe how to find a nearby matrix polynomial of lower McCoy. Another way to formulate having a non-trivial SNF is to solve the minimization problem where must have the appropriate structure. Essentially this finds the smallest perturbation of with an eigenvalue that lowers the rank by at least 2. The auxiliary variables ω and B are used to enforce this constraint. Here is the conjugate
Implementation and examples
We have implemented our algorithms and techniques in the Maple computer algebra system.2 We use the variant of Levenberg-Marquardt discussed in Section 4 in several instances to solve the first-order necessary condition. All computations are done using hardware precision and measured in floating point operations, or FLOPs. The input size of our problem is measured in the dimension and degree of , which are n
Conclusion and topics for future research
In this paper we have shown that the problem of computing a nearby matrix polynomial with a non-trivial spectral structure can be solved by (mostly local) optimization techniques. Regularity conditions were shown to hold for most instances of the problems in question, ensuring that Lagrange multipliers exist under mild assumptions about the solutions. When Lagrange multipliers do not exist, alternative formulations that admit Lagrange multipliers have been proposed. Several of these algorithms
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgements
The authors would like to acknowledge the support of the Natural Sciences and Engineering Research Council (NSERC) Canada. Cette recherche a été financée par le Conseil de Recherches en Sciences Naturelles et en Génie du Canada (CRSNG). The second author would also like to acknowledge the support of the NSF (CCF-1708884) and NSA (H98230-18-1-0016).
References (37)
- et al.
Pseudospectra, critical points and multiple eigenvalues of matrix polynomials
Linear Algebra Appl.
(2009) - et al.
When are two numerical polynomials relatively prime?
J. Symb. Comput.
(1998) - et al.
An improved algorithm for the computation of Kronecker's canonical form of a singular pencil
Linear Algebra Appl.
(1988) - et al.
The dimension of matrices (matrix pencils) with given Jordan (Kronecker) canonical forms
Linear Algebra Appl.
(1995) Perturbation theory for rectangular matrix pencils
Linear Algebra Appl.
(1994)The computation of Kronecker's canonical form of a singular pencil
Linear Algebra Appl.
(1979)- et al.
The eigenstructure of an arbitrary polynomial matrix: computational aspects
Linear Algebra Appl.
(1983) Nonlinear Programming
(1999)- et al.
Computational complexity of μ calculation
IEEE Trans. Autom. Control
(1994) - et al.
The generalized Schur decomposition of an arbitrary pencil : robust software with error bounds and applications. Part I: theory and algorithms
ACM Trans. Math. Softw.
(1993)
The generalized Schur decomposition of an arbitrary pencil B: robust software with error bounds and applications. Part II: software and applications
ACM Trans. Math. Softw.
A geometric approach to perturbation theory of matrices and matrix pencils. Part I: versal deformations
SIAM J. Matrix Anal. Appl.
A geometric approach to perturbation theory of matrices and matrix pencils. Part II: a stratification-enhanced staircase algorithm
SIAM J. Matrix Anal. Appl.
On the quadratic convergence of the Levenberg-Marquardt method without nonsingularity assumption
Computing
Resultant properties of gcd of many polynomials and a factorization representation of gcd
Int. J. Control
Computing approximate greatest common right divisors of differential polynomials
Found. Comput. Math.
Computing the nearest rank-deficient matrix polynomial
Computing nearby non-trivial Smith forms
Cited by (3)
Computation of the nearest non-prime polynomial matrix: Structured low-rank approximation approach
2021, Linear Algebra and Its ApplicationsApproximate GCD by relaxed NewtonSLRA algorithm
2021, ACM Communications in Computer AlgebraRelaxed NewtonSLRA for Approximate GCD
2021, Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)