Elsevier

Journal of Symbolic Computation

Volume 102, January–February 2021, Pages 304-327
Journal of Symbolic Computation

Computing nearby non-trivial Smith forms

https://doi.org/10.1016/j.jsc.2019.10.019Get rights and content

Abstract

We consider the problem of computing the nearest matrix polynomial with a non-trivial Smith Normal Form. We show that computing the Smith form of a matrix polynomial is amenable to numeric computation as an optimization problem. Furthermore, we describe an effective optimization technique to find a nearby matrix polynomial with a non-trivial Smith form. The results are then generalized to include the computation of a matrix polynomial having a maximum specified number of ones in the Smith Form (i.e., with a maximum specified McCoy rank).

We discuss the geometry and existence of solutions and how our results can be used for an error analysis. We develop an optimization-based approach and demonstrate an iterative numerical method for computing a nearby matrix polynomial with the desired spectral properties. We also describe an implementation of our algorithms and demonstrate the robustness with examples in Maple.

Introduction

For a given square matrix polynomial AR[t]n×n, one can find unimodular matrices U,VR[t]n×n such that UAV is a diagonal matrix S. Unimodular means that there is a polynomial inverse matrix, or equivalently, that the determinant is a nonzero constant from R. The unimodular matrices U,V 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. HereS=(s1s2sn)R[t]n×n, where s1,,snR[t] are monic and si|si+1 for 1i<n. The Smith form always exists and is unique though the associated unimodular matrices U, V are not unique (see, e.g., (Kailath, 1980; Gohberg et al., 2009)). The diagonal entries s1,,sn are referred to as the invariant factors of A.

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 R, 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 s1,,sn of a matrix AR[t]n×n as ratios si=δi/δi1 of the determinantal divisors δ0,δ1,,δnR[t], whereδ0=1,δi=GCD{all i×i minors of A}R[t]. In the case of 2×2 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, considerA=(t22t+1t2+2t+2)=diag(f,g)R[t]2×2. If we look for nearby polynomials f˜,g˜R[t] of degree at most 2 such that gcd(f˜,g˜)=γt+1 (i.e. a nontrivial gcd) at minimal distance ff˜22+gg˜22 for some γR, then it is shown in (Haraldson, 2015, Example 3.3.6) that this distance is (5γ44γ3+14γ2+2)/(γ4+γ2+1). This distance has an infimum of 2 at γ=0. However at γ=0 we have gcd(f˜,g˜)=1 even though deg(gcd(f˜,g˜))>0 for all γ0. For A to have a non-trivial Smith form we must perturb f,g 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 A.

Our work indirectly involves measuring the sensitivity to the eigenvalues of A and the determinant of A. 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 A 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 A has an interesting Smith form, providing bounds on a “radius of triviality” around A 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 A having a non-trivial SNF is to solve the minimization problemminΔAF2subject to(A(ω)+ΔA(ω))B=0 and BB=I2,for some ωC and BCn×2, where ΔA must have the appropriate structure. Essentially this finds the smallest perturbation of A with an eigenvalue that lowers the rank by at least 2. The auxiliary variables ω and B are used to enforce this constraint. Here B 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 A, 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)

  • J. Demmel et al.

    The generalized Schur decomposition of an arbitrary pencil AλB: robust software with error bounds and applications. Part II: software and applications

    ACM Trans. Math. Softw.

    (1993)
  • A. Edelman et al.

    A geometric approach to perturbation theory of matrices and matrix pencils. Part I: versal deformations

    SIAM J. Matrix Anal. Appl.

    (1997)
  • A. Edelman et al.

    A geometric approach to perturbation theory of matrices and matrix pencils. Part II: a stratification-enhanced staircase algorithm

    SIAM J. Matrix Anal. Appl.

    (1999)
  • J.-Y. Fan et al.

    On the quadratic convergence of the Levenberg-Marquardt method without nonsingularity assumption

    Computing

    (2005)
  • S. Fatouros et al.

    Resultant properties of gcd of many polynomials and a factorization representation of gcd

    Int. J. Control

    (2003)
  • M. Giesbrecht et al.

    Computing approximate greatest common right divisors of differential polynomials

    Found. Comput. Math.

    (2019)
  • M. Giesbrecht et al.

    Computing the nearest rank-deficient matrix polynomial

  • M. Giesbrecht et al.

    Computing nearby non-trivial Smith forms

  • Cited by (3)

    • Approximate GCD by relaxed NewtonSLRA algorithm

      2021, ACM Communications in Computer Algebra
    • Relaxed NewtonSLRA for Approximate GCD

      2021, Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)
    View full text