Computing Nearby Non-trivial Smith Forms

TL;DR

This paper introduces an optimization-based numerical method for computing the nearest matrix polynomial with a non-trivial Smith Normal Form, enabling error analysis and spectral property control.

Contribution

It presents a novel optimization approach to find nearby matrix polynomials with specified Smith form properties, including maximum McCoy rank, with practical implementation details.

Findings

Effective optimization technique for Smith form computation.

Demonstrated robustness of the method with Maple examples.

Generalization to maximum McCoy rank cases.

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.

Equations

$$\mathcal{S}=\left(\begin{array}{cccc}{s}_1& & & \\ & s_2& & \\ & & \ddots & \\ & & & s_n\end{array}\right)\in \mathbb{R}[t{]}^{n\times n},$$

$$\delta_{0}=1,~{}~{}\delta_{i}=\text{GCD}\biggl{\{}\text{all $i\times i$ minors of ${{\mathcal{A}}}$}\biggr{\}}\in{\mathbb{R}}[t].$$

$${{\mathcal{A}}}=\left(\begin{array}[]{cc}t^{2}-2t+1&\\ &t^{2}+2t+2\end{array}\right)=\operatorname{diag}(f,g)\in{\mathbb{R}}[t]^{2\times 2}.$$

$$\mathcal{K}=\bigl{\{}C_{0,0},\ldots,C_{0,k},tC_{1,0},\ldots,tC_{1,k},\ldots,t^{d}C_{d,0},\ldots,t^{d}C_{d,k}\bigr{\}},$$

$$p=\sum _{j=0}^d{p}_j{t}^t\in \mathbb{R}[t]\mapsto \text{mbox}vec(p)=(p_0,p_1,\dots ,p_d{)}^T\in {\mathbb{R}}^{(d+1)\times 1}$$

$$\mathcal{A}\in \mathbb{R}[t{]}^{m\times n}\mapsto \text{mbox}vec(\mathcal{A})=\left(\begin{array}{c}\text{mbox}vec({\mathcal{A}}_{11})\\ ⋮\\ \text{mbox}vec({\mathcal{A}}_{mn})\end{array}\right)\in {\mathbb{R}}^{mn(d+1)\times 1}.$$

$$\mathcal{A}\in \mathbb{R}[t{]}^{m\times n}\mapsto \text{mbox}pvec(\mathcal{A})=\left(\begin{array}{c}{\mathcal{A}}_{11}\\ ⋮\\ {\mathcal{A}}_{mn}\end{array}\right)\in \mathbb{R}[t{]}^{mn\times 1}.$$

$${\varphi }_{d_2}(a)=\left(\begin{array}{ccc}{a}_0& & \\ ⋮& \ddots & \\ a_{d_1}& & a_0\\ & \ddots & ⋮\\ & & a_{d_1}\end{array}\right)\in {\mathbb{R}}^{(d_1+d_2+1)\times (d_2+1)}.\text{mbox}Itfollowsthat\text{mbox}vec(ab)={\varphi }_{d_2}(a)\text{mbox}vec(b).$$

$$\text{mbox}vec(\mathcal{A}b)={\Phi }_{d_2}(\mathcal{A})\text{mbox}vec(b),$$

$${\Phi }_{d_2}(\mathcal{A})=\left(\begin{array}{ccc}{\varphi }_{d_2}({\mathcal{A}}_{11})& \cdots & {\varphi }_{d_2}({\mathcal{A}}_{1n})\\ ⋮& & ⋮\\ {\varphi }_{d_2}({\mathcal{A}}_{m1})& \cdots & {\varphi }_{d_2}({\mathcal{A}}_{mn})\end{array}\right)\in {\mathbb{R}}^{m(d_1+d_2+1)\times n(d_2+1)}.$$

$$\mathcal{A}\otimes \mathcal{B}=\left(\begin{array}{ccc}{\mathcal{A}}_{11}\mathcal{B}& \cdots & {\mathcal{A}}_{1n}\mathcal{B}\\ ⋮& & ⋮\\ {\mathcal{A}}_{m1}\mathcal{B}& \cdots & {\mathcal{A}}_{mn}\mathcal{B}\end{array}\right)\in \mathbb{R}[t{]}^{mk\times n\ell }.$$

$$\text{mbox}vec(AXB)=(B^T\otimes A)\text{mbox}vec(X).$$

$$\text{mbox}pvec(\mathcal{A}\mathcal{X}\mathcal{B})=({\mathcal{B}}^T\otimes \mathcal{A})\text{mbox}pvec(\mathcal{X}).$$

$$\text{mbox}vec(AX)=(X^T\otimes I)\text{mbox}vec(A)=(I\otimes A)\text{mbox}vec(X)=\text{mbox}vec(B).$$

$$\begin{gathered} J_{\det }=\frac{\partial \text{mbox}vec(\det (\mathcal{A}))}{\partial \text{mbox}vec(\mathcal{A})} \\ ={\Phi }_d(\text{mbox}pvec(\mathrm{Adj}(\mathcal{A}{)}^T{)}^T)\in {\mathbb{R}}^{(nd+1)\times n^2(d+1)}. \end{gathered}$$

$$\det (\mathcal{A}+\Delta \mathcal{A})=\det (\mathcal{A})+\text{mbox}pvec(\mathrm{Adj}(\mathcal{A}{)}^T{)}^T\text{mbox}pvec(\Delta \mathcal{A})+O({\Vert \Delta \mathcal{A}\Vert }_F^2),$$

$$\text{mbox}vec(\det (\mathcal{A}+\Delta \mathcal{A}))\approx \text{mbox}vec(\det (\mathcal{A}))+\text{mbox}vec(\text{mbox}pvec(\mathrm{Adj}(\mathcal{A}{)}^T{)}^T\text{mbox}pvec(\Delta \mathcal{A})).$$

$$\frac{\partial \text{mbox}vec(\det (\mathcal{A}))}{\partial \text{mbox}vec(\mathcal{A})}=\nabla (\det (\mathcal{A}))={\Phi }_d(\text{mbox}pvec(\mathrm{Adj}(\mathcal{A}{)}^T{)}^T).\text{qed}$$

$$J_{\mathrm{Adj}}={\left[{\Phi }_{(n-1)d}(I\otimes \mathcal{A})\right]}^{+}\left[{\Phi }_d(\text{mbox}pvec(I)\text{mbox}pvec(\mathrm{Adj}(\mathcal{A}{)}^T{)}^T)-{\Phi }_d(\mathrm{Adj}(\mathcal{A}{)}^T\otimes I)\right],$$

$$\text{mbox}pvec(\mathcal{A}\mathrm{Adj}(\mathcal{A}))=(\mathrm{Adj}(\mathcal{A}{)}^T\otimes I)\text{mbox}pvec(\mathcal{A})=(I\otimes \mathcal{A})\text{mbox}pvec(\mathrm{Adj}(\mathcal{A})),$$

$$\text{mbox}vec(\mathcal{A}\mathrm{Adj}(\mathcal{A}))={\Phi }_{(n-1)d}(I\otimes \mathcal{A})\text{mbox}vec(\mathrm{Adj}(\mathcal{A}))={\Phi }_d(\mathrm{Adj}(\mathcal{A}{)}^T\otimes I)\text{mbox}vec(\mathcal{A}).$$

$$\partial \text{mbox}vec(\mathcal{A}\mathrm{Adj}(\mathcal{A}))=(\partial {\Phi }_{(n-1)d}(I\otimes \mathcal{A}))\text{mbox}vec(\mathrm{Adj}(\mathcal{A}))+{\Phi }_{(n-1)d}(I\otimes \mathcal{A})\partial \text{mbox}vec(\mathrm{Adj}(\mathcal{A})).$$

$$\text{mbox}vec((\partial \mathcal{A})\mathrm{Adj}(\mathcal{A})+\mathcal{A}(\partial \mathrm{Adj}(\mathcal{A})))$$

$$\text{mbox}vec((\mathrm{Adj}(\mathcal{A}{)}^T\otimes I)\text{mbox}pvec(\partial \mathcal{A})+(I\otimes \mathcal{A})\text{mbox}pvec(\partial \mathrm{Adj}(\mathcal{A}))),$$

$${\Phi }_d((\mathrm{Adj}(\mathcal{A}{)}^T\otimes I))\text{mbox}vec(\partial \mathcal{A})+{\Phi }_{(n-1)d}(I\otimes A)\text{mbox}vec(\partial \mathrm{Adj}(\mathcal{A})).$$

$$\partial \text{mbox}vec(\det (\mathcal{A})I)=\text{mbox}vec(\partial \text{mbox}pvec(\det (\mathcal{A})I)),$$

$$\text{mbox}vec(\partial \text{mbox}pvec(\det (\mathcal{A})I))=\text{mbox}vec(\text{mbox}pvec(I)\text{mbox}pvec(\mathrm{Adj}(\mathcal{A}{)}^T{)}^T\text{mbox}pvec(\partial \mathcal{A})).$$

$$\text{mbox}vec(\text{mbox}pvec(I)\text{mbox}pvec(\mathrm{Adj}(\mathcal{A}{)}^T{)}^T)\text{mbox}pvec(\partial \mathcal{A})={\Phi }_d(\text{mbox}pvec(I)\text{mbox}pvec(\mathrm{Adj}(\mathcal{A}{)}^T{)}^T)\text{mbox}vec(\partial \mathcal{A}),$$

$${\Phi }_{(n-1)d}(I\otimes \mathcal{A})\frac{\partial \text{mbox}vec(\mathrm{Adj}(\mathcal{A}))}{\partial \text{mbox}vec(\mathcal{A})}={\Phi }_d(\text{mbox}pvec(I)\text{mbox}pvec(\mathrm{Adj}(\mathcal{A}{)}^T{)}^T)-{\Phi }_d(\mathrm{Adj}(\mathcal{A}{)}^T\otimes I).$$

$$J_{\mathrm{Adj}}={\left[{\Phi }_{(n-1)d}(I\otimes \mathcal{A})\right]}^{+}\left[{\Phi }_d(\text{mbox}pvec(I)\text{mbox}pvec(\mathrm{Adj}(\mathcal{A}{)}^T{)}^T)-{\Phi }_d(\mathrm{Adj}(\mathcal{A}{)}^T\otimes I)\right],$$

Full text

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.

url]https://cs.uwaterloo.ca/ mwg

url]https://cs.uwaterloo.ca/ jharalds

url]https://cs.uwaterloo.ca/ glabahn

1 Introduction

For a given square matrix polynomial ${{\mathcal{A}}}\in{\mathbb{R}}[t]^{n\times n}$, one can find unimodular matrices ${\mathcal{U}},{\mathcal{V}}\in{\mathbb{R}}[t]^{n\times n}$ such that ${\mathcal{U}}{{\mathcal{A}}}{\mathcal{V}}$ is a diagonal matrix ${\mathcal{S}}$. Unimodular means that there is a polynomial inverse matrix, or equivalently, that the determinant is a nonzero constant from ${\mathbb{R}}$. The unimodular matrices ${\mathcal{U}},{\mathcal{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. Here

[TABLE]

where $s_{1},\ldots,s_{n}\in{\mathbb{R}}[t]$ are monic and $s_{i}{\mskip 3.0mu|\mskip 3.0mu}s_{i+1}$ for $1\leq i<n$. The Smith form always exists and is unique though the associated unimodular matrices ${\mathcal{U}}$, ${\mathcal{V}}$ are not unique (see, e.g., (Kailath, 1980; Gohberg et al., 2009)). The diagonal entries $s_{1},\ldots,s_{n}$ are referred to as the invariant factors of ${{\mathcal{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 ${\mathbb{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, b; 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, Haraldson, and Labahn, 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 $s_{1},\ldots,s_{n}$ of a matrix ${{\mathcal{A}}}\in{\mathbb{R}}[t]^{n\times n}$ as ratios $s_{i}=\delta_{i}/\delta_{i-1}$ of the determinantal divisors $\delta_{0},\delta_{1},\ldots,\delta_{n}\in{\mathbb{R}}[t]$, where

[TABLE]

In the case of $2\times 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, Haraldson, and Kaltofen, 2019a)).

The issue of unattainable infima extends to the Smith normal form. As an example, consider

[TABLE]

If we look for nearby polynomials $\widetilde{f},\widetilde{g}\in{\mathbb{R}}[t]$ of degree at most $2$ such that $\gcd(\widetilde{f},\widetilde{g})=\gamma t+1$ (i.e. a nontrivial gcd) at minimal distance ${\|f-\widetilde{f}\|}_{2}^{2}+{\|g-\widetilde{g}\|}_{2}^{2}$ for some $\gamma\in{\mathbb{R}}$, then it is shown in (Haraldson, 2015, Example 3.3.6) that this distance is $(5\gamma^{4}-4\gamma^{3}+14\gamma^{2}+2)/(\gamma^{4}+\gamma^{2}+1)$. This distance has an infimum of $2$ at $\gamma=0$. However at $\gamma=0$ we have $\gcd(\widetilde{f},\widetilde{g})=1$ even though $\deg(\gcd(\widetilde{f},\widetilde{g}))>0$ for all $\gamma\neq 0$. For ${{\mathcal{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 ${{\mathcal{A}}}$.

Our work indirectly involves measuring the sensitivity to the eigenvalues of ${{\mathcal{A}}}$ and the determinant of ${{\mathcal{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 ${{\mathcal{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, 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 ${{\mathcal{A}}}$ has an interesting Smith form, providing bounds on a “radius of triviality” around ${{\mathcal{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, Haraldson, and Labahn, 2018).

2 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.

2.1 Notation and Terminology

We make extensive use of the following terminology and definitions. A matrix polynomial ${{\mathcal{A}}}\in{\mathbb{R}}[t]^{n\times n}$ is an ${n\times n}$ matrix whose entries are polynomials. Typically we also work with matrices whose entries have degree bound $d$ and use the notation ${\mathbb{R}}_{\leq d}[t]^{n\times n}$ to describe this set. Alternatively, we may express matrix polynomials as ${{\mathcal{A}}}=\sum_{1\leq j\leq d}{A}_{j}t^{j}$ where ${A}_{j}\in{\mathbb{R}}^{n\times n}$. The degree of a matrix polynomial $d$ is defined to be the degree of the highest-order non-zero entry of ${{\mathcal{A}}}$, or the largest index $j$ such that ${A}_{j}\neq 0$. We say that ${{\mathcal{A}}}$ has full rank or is regular if $\det({{\mathcal{A}}})\nequiv 0$. As noted in the introduction, ${{\mathcal{A}}}$ is said to be unimodular if $\det({{\mathcal{A}}})\in{\mathbb{R}}\backslash\{0\}$. The (finite) eigenvalues are the roots of $\det(A)\in{\mathbb{R}}[t]$. The norm of a polynomial $a\in{\mathbb{R}}[t]$ is defined as ${\|a\|}={\|a\|}_{2}={\|(a_{0},a_{1},\ldots,a_{d},0,\ldots,0)\|}_{2}$. For matrix polynomials we define ${\|{{\mathcal{A}}}\|}={\|{{\mathcal{A}}}\|}_{F}=\sqrt{\sum_{i,j}{\|{\mathcal{A_{i,j}}}\|}_{2}^{2}}$. Our choice of norm is a distributed coefficient norm, sometimes known as the Frobenius norm.

Definition 2.1** (SVD – Golub and Van Loan 2012).**

The Singular Value Decomposition (SVD) of $A\in{\mathbb{R}}^{n\times n}$ is given by $U^{T}\Sigma V$, where $U,V\in{\mathbb{R}}^{n\times n}$ satisfy $U^{T}U=I$, $V^{T}V=I$ and $\Sigma=\operatorname{diag}(\sigma_{1},\ldots,\sigma_{n})$ is a diagonal matrix with non-negative real entries in descending order of magnitude, the singular values of $A$. The distance to the nearest (unstructured) matrix of rank $m<n$ is $\sigma_{m+1}(A)$.

For scalar matrices we frequently write ${\|\cdot\|}_{2}$ for the largest singular value, and $\sigma_{\min}(\cdot)$ for the smallest singular value.

In this paper we are mainly concerned with coefficient structures that preserve the zero coefficient structure of a matrix polynomial, that is, we generally do not change zero coefficients to non-zero, or increase the degrees of matrix entries.

Definition 2.2** (Affine/Linear Structure).**

A non-zero matrix polynomial ${{\mathcal{A}}}\in{\mathsf{R}}[t]^{n\times n}$ of degree at most $d$ has a linear structure from a set $\mathcal{K}$ if ${{\mathcal{A}}}\in\mbox{{span}}(\mathcal{K})$ as a vector space over ${\mathsf{R}}$, where

[TABLE]

where $C_{l,j}\in{\mathsf{R}}^{n\times n}$ for $0\leq j\leq k$, where $k>0$ is a finite index variable. If ${{\mathcal{A}}}=\mathcal{C}_{0}+\mathcal{C}_{1}$, where $\mathcal{C}_{0}\in{\mathsf{R}}[t]^{n\times n}$ is fixed and $\mathcal{C}_{1}\in\mbox{{span}}(\mathcal{K})$, then ${{\mathcal{A}}}$ is said to have an affine structure from the set $\mathcal{K}$.

Linearly and affine linearly structured matrices are best thought of as imposing linear equality constraints on the entries. Examples of matrices with a linear structure include matrices with prescribed zero entries or coefficients, Toeplitz/Hankel matrices, Sylvester matrices, resultant-like matrices, Ruppert matrices and several other matrices appearing in symbolic-numeric computation. Matrices with an affine structure include all matrices with a linear structure and, in addition, matrices having prescribed non-zero constant entries/coefficients, for example monic matrix polynomials.

Recall that the rank of a matrix polynomial is the maximum number of linearly independent rows or columns as a vector space over ${\mathbb{R}}(t)$. This is the same as the rank of the matrix ${{\mathcal{A}}}(\omega)$ for any $\omega\in{\mathbb{C}}$ that is not an eigenvalue of ${{\mathcal{A}}}(t)$. If we allow evaluation at eigenvalues, then the McCoy rank is the lowest rank when ${{\mathcal{A}}}$ is evaluated at an eigenvalue.

Definition 2.3** (McCoy Rank and Non-Trivial SNF).**

The McCoy rank of ${{\mathcal{A}}}$ is $\min_{\omega\in{\mathbb{C}}}\{\operatorname{rank}{{\mathcal{A}}}(\omega)\}$, the lowest rank possible when ${{\mathcal{A}}}$ is evaluated at any $\omega\in{\mathbb{C}}$. Note that the rank of ${{\mathcal{A}}}$ only drops at all if it is evaluated at an eigenvalue $\omega\in{\mathbb{C}}$. The McCoy rank is also the number of ones in the Smith form. Equivalently, if ${{\mathcal{A}}}$ has $r$ non-trivial invariant factors, then the McCoy rank of ${{\mathcal{A}}}$ is $n-r$. The matrix polynomial ${{\mathcal{A}}}$ is said to have a non-trivial or interesting Smith normal form if the McCoy rank is at most $n-2$, or equivalently, if it has two or more invariant factors of non-zero degree.

Problem 2.4** (Approximate SNF Problem).**

Given a matrix polynomial ${{\mathcal{A}}}\in{\mathbb{R}}[t]^{n\times n}$, find the distance to a non-trivial SNF. That is, find a matrix polynomial $\widehat{{{\mathcal{A}}}}\in{\mathbb{R}}[t]^{n\times n}$ of prescribed coefficient structure that has a prescribed McCoy rank of at most $n-r$, for $r\geq 2$, such that ${\|{{\mathcal{A}}}-{\mathcal{\widehat{A}}}\|}$ is minimized under ${\|\cdot\|}$, if such an $\widehat{{{\mathcal{A}}}}$ exists.

We will consider ${\|\cdot\|}={\|\cdot\|}_{F}$ to be the Frobenius norm. The nearest matrix of McCoy rank at most $n-2$, if it exists, is called the approximate SNF.

Problem 2.5** (Lower McCoy Rank Approximation Problem).**

Compute the nearest matrix with McCoy rank $n-r$ matrix, if it exists, for $r\geq 3$.

In a generic sense, the nearest matrix polynomial with an interesting SNF will have McCoy rank $n-2$ with probability one, but many matrices arising from applications are expected to have more interesting (i.e. the invariant factors have a richer or non-generic multiplicity structure) Smith forms nearby.

As described in the introduction, it is possible that the distance to a non-trivial SNF is not attainable. That is, there is a solution that is approached asymptotically, but where the Smith form is trivial at the infimum. Fortunately, in most instances of interest, solutions will generally be attainable. We will also later discuss how to identify and compute unattainable solutions. Problem 2.4 and Problem 2.5 admit the nearest rank $n-1$ or rank $n-2$ matrix polynomial as a special case. However, the computational challenges are fundamentally different for non-trivial instances.

2.2 Embedding into Scalar Domains

In our study of nearest non-trivial Smith forms we often make use of the representation of the diagonal elements as ratios of GCDs of sub-determinants. When the coefficients are polynomials with numeric coefficients it is helpful to embed the arithmetic operations of polynomial multiplication and polynomial GCD into a matrix problem having numeric coefficients (i.e., from ${\mathbb{R}}$). Such an embedding allows us to employ a number of tools, including condition numbers and perturbations of matrix functions.

We start with some basic notation for mapping matrices and polynomials to vectors.

Definition 2.6** (Vec Operator).**

We define the operator ${\mbox{vec}}:{\mathbb{R}}[t]\to{\mathbb{R}}^{(d+1)\times 1}$ as follows:

[TABLE]

The vec operator ${\mbox{vec}}(\cdot)$ is extended to map a matrix ${\mathbb{R}}[t]^{m\times n}$ to a single vector in ${\mathbb{R}}^{mn(d+1)\times 1}$ by stacking columns of (padded) coefficient vectors on top of each other.

[TABLE]

It is sometimes useful to reduce matrix polynomials to vectors of polynomials in ${\mathbb{R}}[t]$ rather than vectors over ${\mathbb{R}}$.

Definition 2.7** (Polynomial Vec Operator).**

The pvec operator maps ${\mathbb{R}}[t]^{m\times n}$ to a vector ${\mathbb{R}}[t]^{nm\times 1}$ as

[TABLE]

We define the vectorization of matrix polynomials in this somewhat non-standard way so that we can later facilitate the computation of derivatives of matrix polynomial valued functions.

To describe polynomial multiplication in terms of linear maps over scalars we have:

Definition 2.8** (Convolution Matrix).**

Polynomial multiplication between polynomials $a,b\in{\mathbb{R}}[t]$, of degrees $d_{1}$ and $d_{2}$, respectively may be expressed as a Toeplitz-matrix-vector product. We define

[TABLE]

When $a$ is non-zero, we can also define division through pseudo-inversion or linear least squares. In a similar manner, we can define the product of matrix polynomials through a Toeplitz-block matrix.

Definition 2.9** (Block Convolution Matrix).**

We can express multiplication of a matrix and vector of polynomials, ${{\mathcal{A}}}\in{\mathbb{R}}[t]^{m\times n}$ and ${\mathcal{b}}\in{\mathbb{R}}[t]^{n\times 1}$, of degrees at most $d_{1}$ and $d_{2}$ respectively, as a scalar linear system

[TABLE]

where

[TABLE]

The block convolution matrix is sometimes referred to as a “Sylvester matrix” associated with ${{\mathcal{A}}}$. However, we reserve the term “Sylvester matrix” for the more standard linear system appearing from the GCD of two (or more) polynomials. The block convolution matrix is a scalar matrix whose entries have a linear (Toeplitz-block) structure.

Definition 2.10** (Kronecker Product).**

The Kronecker product of ${{\mathcal{A}}}\in{\mathbb{R}}[t]^{m\times n}$ and ${{\mathcal{B}}}\in{\mathbb{R}}[t]^{k\times\ell}$ denoted as ${{\mathcal{A}}}\otimes{{\mathcal{B}}}$ is the $mk\times n\ell$ matrix over ${\mathbb{R}}[t]$ defined as

[TABLE]

This definition of Kronecker product, sometimes referred to as the “outer product”, also holds for scalar matrices (and vectors).

Lemma 2.11**.**

For scalar matrices of compatible dimension $A,X$ and $B$ over ${\mathbb{R}}$, we can

[TABLE]

Likewise, for matrix polynomials ${{\mathcal{A}}},{\mathcal{X}}$ and ${{\mathcal{B}}}$ of compatible dimension over ${\mathbb{R}}[t]$, we have

[TABLE]

The Kronecker product can also be used to re-write matrix equations of the form $AX=B$, for matrices $A$, $B$ and $X$ of compatible dimensions, to

[TABLE]

2.3 Derivatives of Matrix Polynomial Valued Functions

In this paper we will need to compute derivatives of some important matrix polynomial valued functions, namely the determinant and adjoint. This problem is approached in the context of computing the Jacobian matrix of a vector valued function. The analysis in this section will be useful for showing that Lagrange multipliers typically exist in the optimization problems encountered. The quantities computed can also be used to derive first-order perturbation bounds for these matrix polynomial valued functions with respect to ${\|\cdot\|}_{F}$.

Recall that the adjoint of a matrix polynomial ${{\mathcal{A}}}\in{\mathbb{R}}[t]^{n\times n}$, denoted by $\operatorname{Adj}({{\mathcal{A}}})\in{\mathbb{R}}[t]^{n\times n}$, is the transpose of the cofactor matrix. Thus $\operatorname{Adj}({{\mathcal{A}}})_{ij}=(-1)^{i+j}\det({{\mathcal{A}}}[j|i])$ where ${{\mathcal{A}}}[j|i]$ is the $(j,i)$ cofactor of ${{\mathcal{A}}}$, that is, the matrix formed by removing row $j$ and column $i$ from ${{\mathcal{A}}}$. When ${{\mathcal{A}}}$ has full rank, ${{\mathcal{A}}}$ satisfies ${{\mathcal{A}}}\operatorname{Adj}({{\mathcal{A}}})=\det({{\mathcal{A}}})I$.

The determinant of an $n\times n$ matrix polynomial having entries of degree at most $d$ can be viewed as a mapping from ${\mathbb{R}}^{n^{2}(d+1)}\to{\mathbb{R}}^{nd+1}$, since the determinant has degree at most $nd$. With this same viewpoint, we can view the adjoint of a matrix polynomial as a mapping from ${\mathbb{R}}^{n^{2}(d+1)}\to{\mathbb{R}}^{n^{2}((n-1)d+1)}$, since the degree of the entries of the adjoint are at most $(n-1)d$. Our notation for computing derivatives of vector valued functions follows that of (Magnus and Neudecker, 1988).

It is not surprising that the determinant of a matrix polynomial has a similar identity (Magnus and Neudecker, 1988) to the well-known scalar identity $\nabla\det(A)={\mbox{vec}}((\operatorname{Adj}(A)^{T})^{T}).$

Theorem 2.12**.**

Let ${{\mathcal{A}}}\in{\mathbb{R}}[t]^{n\times n}$ have degree at most $d$, then

[TABLE]

Proof.

We note that from generalizing the scalar identity $\nabla\det(\cdot)={\mbox{vec}}(\operatorname{Adj}(\cdot)^{T})^{T}$, we can write a first-order expansion of the determinant as

[TABLE]

and ignoring higher-order terms we obtain the scalar expression

[TABLE]

The Jacobian can be extracted by (padding with zero coefficient entries as necessary) writing ${\mbox{vec}}({\mbox{pvec}}(\operatorname{Adj}({{\mathcal{A}}})^{T})^{T}{\mbox{pvec}}(\Delta{{\mathcal{A}}}))=J_{\det}{\mbox{vec}}(\Delta{{\mathcal{A}}})$ as a matrix-vector product. Thus, using block-convolution matrices we have

[TABLE]

Now that we have a closed-form expression for the derivative of the determinant, it is useful to derive a closed-form expression for the adjoint matrix. The closed-form expression reveals rank information, and is independently useful for optimization algorithms requiring derivatives. The rank information is useful to obtain insights about the existence of Lagrange multipliers. If $J_{\det}$ has full or locally constant (row) rank then constraint qualifications will hold for several constrained optimization problems involving the determinant. If $\operatorname{Adj}({{\mathcal{A}}})$ is non-zero then one can often infer the existence of Lagrange multipliers for other problems as well.

Theorem 2.13**.**

Let ${{\mathcal{A}}}\in{\mathbb{R}}[t]^{n\times n}$ have degree at most $d$ and rank $n$. The Jacobian of $\operatorname{Adj}({{\mathcal{A}}})$ is $J_{\operatorname{Adj}}\in{\mathbb{R}}^{(n^{2}((n-1)d+1))\times n^{2}(d+1)}$ with

[TABLE]

where $I$ is understood to be the $n\times n$ identity matrix and for a scalar matrix $A$ of full rank, $A^{+}$ is the Moore-Penrose pseudo-inverse arising from the SVD.

Proof.

First recall that if ${{\mathcal{A}}}$ has full rank, then ${{\mathcal{A}}}\operatorname{Adj}({{\mathcal{A}}})=\operatorname{Adj}({{\mathcal{A}}}){{\mathcal{A}}}=\det({{\mathcal{A}}})I$. This expression defines the adjoint matrix when ${{\mathcal{A}}}$ has full rank. We can write

[TABLE]

thus converting to a linear system over ${\mathbb{R}}$ produces

[TABLE]

Applying the product rule yields

[TABLE]

Next we observe that (1) has the same coefficients as the expression

[TABLE]

which is equivalent to

[TABLE]

which reduces to

[TABLE]

We now have the derivative of the left hand side the expression ${{\mathcal{A}}}\operatorname{Adj}({{\mathcal{A}}})=\det({{\mathcal{A}}})I$. Differentiation of the right hand side yields

[TABLE]

which is equivalent to the expression

[TABLE]

Converting (3) into a linear system over ${\mathbb{R}}$ leads to

[TABLE]

which is the derivative of the right-hand side.

Combining (2) and (4) we have

[TABLE]

Assuming that ${{\mathcal{A}}}$ has full rank so $\Phi_{(n-1)d}({\mbox{pvec}}(I\otimes{{\mathcal{A}}}))$ is pseudo-invertible, we can write

[TABLE]

which completes the proof. ∎

An observation that is important later is that the derivative of the adjoint has a Toeplitz-block structure. More importantly, the bandwidth is $O(d)$, and we only need to compute $O(n^{2})$ columns instead of $O(n^{2}d)$. We also note that $J_{\operatorname{Adj}}$ may be padded with zeros, since ${{\mathcal{A}}}$ may not have generic degrees.

Corollary 2.14**.**

If ${{\mathcal{A}}}$ has full rank then $J_{\operatorname{Adj}}$ has full rank.

Proof.

The matrix $\Phi_{(n-1)d}(I\otimes{{\mathcal{A}}})$ has full rank since $I\otimes{{\mathcal{A}}}$ has full rank. The matrix

[TABLE]

is a rank one update to a matrix polynomial. By evaluating (5) at a complex number $\omega$ that is not an eigenvalue of ${{\mathcal{A}}}$ we can show that (5) has full rank. Let $A={{\mathcal{A}}}(\omega)$, so $A\in{\mathbb{R}}^{n\times n}$ has full rank.

Using the Sherman-Morrison formula  (Higham, 2002, pg. 487) for rank 1 updates to a matrix, we need to verify that

[TABLE]

in order to ensure that (5) has full rank. We have that

[TABLE]

thus (5) has full rank. Note we used the identities for matrices $X,Y$ and $Z$ of appropriate dimension, that ${\mbox{vec}}(XYZ)=(Z^{T}\otimes X){\mbox{vec}}(Y)$ and ${\mbox{vec}}(X^{T})^{T}{\mbox{vec}}(Y)={\mbox{Tr}}(XY)$. Again, we have that

[TABLE]

has full rank, thus $J_{\operatorname{Adj}}$ is a product of two matrices of full rank, so $J_{\operatorname{Adj}}$ must also have full rank. ∎

Corollary 2.14 implies that Lagrange multipliers will exist to several optimization problems involving the adjoint matrix as a constraint, since the Jacobian matrix of the adjoint has full rank. The linear independent constraint qualification or the constant rank constraint qualification will hold for several optimization problems of the form

[TABLE]

for some reasonably prescribed ${\mathcal{F}}\in{\mathbb{R}}[t]^{n\times n}$.

Remark 2.15**.**

If ${{\mathcal{A}}}$ is rank deficient, then the derivative is still defined, but not necessarily by Theorem 2.13. If $\operatorname{rank}({{\mathcal{A}}})\leq n-3$ then $J_{\operatorname{Adj}}=0$, since all $(n-3)\times(n-3)$ minors vanish ($J_{\operatorname{Adj}}$ consists of the coefficients of these minors). If $\operatorname{rank}({{\mathcal{A}}})=n-1$ or $\operatorname{rank}({{\mathcal{A}}})=n-2$ then $J_{\operatorname{Adj}}$ is still defined and in both cases $J_{\operatorname{Adj}}\neq 0$. However $J_{\operatorname{Adj}}$ is not necessarily described by Theorem 2.13.

For several affine or linear perturbation structures (such as ones that preserve the degree of entries or the support of entries), Theorem 2.13 and the associated Corollary 2.14 will hold (after deleting some extraneous rows or columns).

3 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.

3.1 Embeddings into generalized Sylvester matrices and approximate GCDs

In the introduction we demonstrated that some nearby non-trivial Smith Forms are unattainable. In this subsection we investigate why these unattainable values occur. We first review some basic results needed to analyze the topology of the approximate Smith form problem.

For a matrix ${{\mathcal{A}}}\in{\mathbb{R}}[x]^{n\times n}$, we know that $s_{n}=\delta_{n}/\delta_{n-1}$, the quotient of the determinant and the GCD of all $(n-1)\times(n-1)$ minors. Since these minors are precisely the entries of the adjoint matrix, it follows that ${{\mathcal{A}}}$ has a non-trivial Smith form if and only if the GCD of all entries of the adjoint is non-trivial, that is, $\deg(\gcd(\{\operatorname{Adj}({{\mathcal{A}}})_{ij}\}))\geq 1$. In order to obtain bounds on the distance to a matrix having a non-trivial Smith form, we consider an approximate GCD problem of the form

[TABLE]

If this was a classical approximate GCD problem, then the use of Sylvester-like matrices would be sufficient. However, in our problem the degrees of the entries of the adjoint may change under perturbations. In order to perform an analysis, we need to study a family of generalized Sylvester matrices that allow higher-degree zero coefficients to be perturbed.

The computation of the GCD of many polynomials is typically embedded into a scalar matrix problem using the classical Sylvester matrix. However, in our case we want to look at GCDs of nearby polynomials but with the added wrinkle that the degrees of the entries of the individual polynomials may change under perturbations. In order to perform such an analysis, we need to study a family of generalized Sylvester matrices that allow higher-degree zero coefficients to be perturbed.

Let ${\bf{f}}=(f_{1},\ldots,$ $f_{k})\in{\mathbb{R}}[t]^{k}$ be a vector of polynomials with degrees ${\bf{d}}=(d_{1},\ldots,d_{k})$ ordered as $d_{j}\geq d_{j+1}$ for $1\leq j\leq k-1$. Set $d=d_{1}$ and $\ell=\max(d_{2},\ldots,d_{k})$ and suppose that for each $i\in\{2,\ldots,k\}$ we have $f_{i}=\sum_{1\leq j\leq\ell}f_{ij}t^{j}$.

Definition 3.1** (Generalized Sylvester Matrix).**

The generalized Sylvester matrix of ${\bf{f}}$ is defined as

[TABLE]

Some authors, e.g., (Fatouros and Karcanias, 2003; Vardulakis and Stoyle, 1978), refer to such a matrix as an expanded Sylvester matrix or generalized resultant matrix. The generalized Sylvester matrix has many useful properties pertaining to the Bézout coefficients. However, we are only concerned with the well known result that $\gcd({\bf{f}})=\gcd(f_{1},\ldots,f_{k})=1$ if and only if $\operatorname{\mbox{{Syl}}}_{{\bf{d}}}({\bf{f}})$ has full rank (Vardulakis and Stoyle, 1978).

Sometimes treatreating a polynomial of degree $d$ as one of larger degree is useful. This can be accomplished by constructing a similar matrix and padding rows and columns with zero entries. The generalized Sylvester matrix of degree at most ${\bf{d}}^{\prime}\geq{\bf{d}}$ (component-wise) of ${\bf{f}}$ is defined analogously as $\operatorname{\mbox{{Syl}}}_{{\bf{d}}^{\prime}}({\bf{f}})$, taking $d$ to be the largest degree entry and $\ell$ to be the largest degree of the remaining entries of ${\bf{d}}^{\prime}$. Note that $\ell=d$ is possible and typical. If the entries of ${\bf{f}}$ have a non-trivial GCD (that is possibly unattainable) under a perturbation structure $\Delta{\bf{f}}$, then it is necessary that $\operatorname{\mbox{{Syl}}}_{{\bf{d}}^{\prime}}({\bf{f}})$ is rank deficient, and often this will be sufficient.

If we view the entries of ${\bf{f}}$ as polynomials of degree ${\bf{d}}^{\prime}$ and $d^{\prime}_{i}>d_{i}$ for all $i$, then the entries of ${\bf{f}}$ have an unattainable GCD of distance zero, typically of the form $1+\varepsilon t\sim t+\varepsilon^{-1}$. In other words, the underlying approximate GCD problem is ill-posed in a sense that the solution is unattainable. In order to study the theory of unattainable GCD’s, sometimes referred to as GCD’s at infinity, we need to study the notion of a degree reversed polynomial.

Lemma 3.2**.**

If $\max({\bf{d}})=\max({\bf{d}}^{\prime})$ then $\operatorname{\mbox{{Syl}}}_{{\bf{d}}}({\bf{f}})$ has full rank if and only if and $\operatorname{\mbox{{Syl}}}_{{\bf{d}}^{\prime}}({\bf{f}})$ has full rank.

Proof.

Let $d$ and $\ell$ be the largest and second largest entries of ${\bf{d}}$ and $\ell^{\prime}$ be the second largest entry of ${\bf{d}}^{\prime}$. The result follows from the main theorem of Vardulakis and Stoyle (1978) by considering the case of $\ell^{\prime}=d$. ∎

This lemma characterizes the (generic) case when elements of maximal degree of ${\bf{f}}$ do not change under perturbations, in which case the generalized Sylvester matrix still meaningfully encodes GCD information. However, it is possible that $\operatorname{\mbox{{Syl}}}_{{\bf{d}}}({\bf{f}})$ has full rank and $\operatorname{\mbox{{Syl}}}_{{\bf{d}}^{\prime}}({\bf{f}})$ is rank deficient but the distance to a non-trivial GCD is not zero. This can occur when $d_{j}=d^{\prime}_{j}$ for some $j$ and ${\bf{d}}^{\prime}\geq{\bf{d}}$. To understand the most general case, we need to look at generalized Sylvester matrices of the reversal of several polynomials.

Definition 3.3**.**

The degree $d$ reversal of $f\in{\mathbb{R}}[t]$ of degree at most $d$ is defined as $\mbox{{rev}}_{d}(f)=t^{d}f(t^{-1})$. For a vector of polynomials ${\bf{f}}\in{\mathbb{R}}[t]^{k}$ of degrees at most ${\bf{d}}=(d_{1},\ldots,d_{k})$ the degree ${\bf{d}}$ reversal of ${\bf{f}}$ is the vector $\mbox{{rev}}_{{\bf{d}}}({\bf{f}})=(\mbox{{rev}}_{d_{1}}(f_{1}),\ldots,\mbox{{rev}}_{d_{k}}(f_{k}))$.

The following theorem enables us to determine if unattainable solutions are occurring in an approximate GCD problem with an arbitrary (possibly non-linear) structure on the coefficients.

Theorem 3.4**.**

Let ${\bf{f}}$ be a vector of non-zero polynomials of degree at most $d$. Suppose that $\operatorname{\mbox{{Syl}}}_{\bf{d}}({\bf{f}})$ has full rank and $\operatorname{\mbox{{Syl}}}_{{\bf{d}}^{\prime}}({\bf{f}})$ is rank deficient, where the perturbations $\Delta{\bf{f}}$ have degrees at most ${\bf{d}}^{\prime}$ and the entries of ${\bf{f}}$ have degrees ${\bf{d}}$. Then ${\bf{f}}$ has an unattainable non-trivial GCD of distance zero under the perturbation structure $\Delta{\bf{f}}$ if and only if $\operatorname{\mbox{{Syl}}}(\mbox{{rev}}_{{\bf{d}}^{\prime}}({\bf{f}}))$ is rank deficient.

Proof.

Suppose that $\operatorname{\mbox{{Syl}}}(\mbox{{rev}}_{{\bf{d}}^{\prime}}({\bf{f}}))$ has full rank. Then $\gcd(\mbox{{rev}}_{{\bf{d}}^{\prime}}({\bf{f}}))=1$, hence ${\bf{f}}$ does not have an unattainable non-trivial GCD, since $\gcd({\bf{f}})=1$. Conversely, suppose that $\operatorname{\mbox{{Syl}}}(\mbox{{rev}}_{{\bf{d}}^{\prime}}({\bf{f}}))$ is rank deficient. Then, $t$ is a factor of $\gcd(\mbox{{rev}}_{{\bf{d}}^{\prime}}({\bf{f}}))$ but $t$ is not a factor of $\gcd(\mbox{{rev}}_{\bf{d}}({\bf{f}}))$. Accordingly, all entries of ${\bf{f}}+\Delta{\bf{f}}$ may increase in degree and so the distance of ${\bf{f}}$ having a non-trivial GCD is zero, and so is unattainable. ∎

If the generalized Sylvester matrix of ${\bf{f}}$ has full rank, but the generalized Sylvester matrix that encodes the perturbations ${\bf{f}}+\Delta{\bf{f}}$ is rank deficient, then either there is an unattainable GCD, or the generalized Sylvester matrix is rank deficient due to over-padding with zeros. Theorem 3.4 provides a reliable way to detect this over-padding.

Definition 3.5**.**

We say that ${{\mathcal{A}}}$ has an unattainable non-trivial Smith form if $\gcd(\operatorname{Adj}({{\mathcal{A}}}))=1$ and $\gcd(\operatorname{Adj}({{\mathcal{A}}}+\widetilde{\Delta}{{\mathcal{A}}}))\neq 1$ for an arbitrarily small perturbation $\widetilde{\Delta}{{\mathcal{A}}}=\widetilde{\Delta}(\Delta{{\mathcal{A}}})$ of some prescribed affine structure.

Note that $\widetilde{\Delta}{{\mathcal{A}}}$ just means that perturbations to ${{\mathcal{A}}}$ are structured as an affine function of $\Delta{{\mathcal{A}}}$. It is important to carefully consider structured perturbations, because some matrix polynomials have an unattainable non-trivial SNF under unstructured perturbations, but have an attainable non-trivial SNF under structured perturbations (perturbations that preserve the degree of entries or support of entries are structured). Solutions that cannot be attained correspond to an eigenvalue at infinity of ${{\mathcal{A}}}$ with a non-trivial spectral structure. Such examples are easily constructed when $\det({{\mathcal{A}}})$ or $\operatorname{Adj}({{\mathcal{A}}})$ have non-generic degrees.

Example 3.6**.**

Let

[TABLE]

Then ${\mathcal{C}}$ has an unattainable non-trivial Smith form if all perturbations to ${{\mathcal{A}}}$ are support or degree preserving (i.e. they preserve zero entries or do not increase the degree of each entry), both linear structures. Note that ${{\mathcal{A}}}$ and ${\mathcal{C}}$ are both unimodular. However small perturbations to the non-zero coefficients of ${{\mathcal{A}}}$ make ${{\mathcal{A}}}+\Delta{{\mathcal{A}}}$ non-unimodular.

The Smith form of $\mbox{{rev}}({\mathcal{C}})=t{\mathcal{C}}|_{t=t^{-1}}$ is

[TABLE]

which implies that the eigenvalue at infinity of ${{\mathcal{A}}}$ has a non-trivial spectral structure. The eigenvalue at infinity having a non-trivial spectral structure implies that the SNF of ${\mathcal{C}}$ is unattainable. Note that this is equivalent to saying that ${\mathcal{C}}$ has a non-trivial Smith-McMillan form.

These examples are non-generic. Generically, the degree of all entries in the adjoint will be $(n-1)d$ and will remain unchanged locally under perturbations to the coefficients. Computing the distance to the nearest matrix polynomial with a non-trivial Smith form under a prescribed perturbation structure can be formulated as finding the nearest rank deficient (structured) generalized Sylvester matrix of the adjoint or the ${\bf{d}}^{\prime}$ reversal of the adjoint.

3.2 Nearest Rank Deficient Structured Generalized Sylvester Matrix

Suppose that ${{\mathcal{A}}}\in{\mathbb{R}}[t]^{n\times n}$ of degree at most $d$ has a trivial Smith form and does not have an unattainable non-trivial Smith form. Then one method to compute a lower bound on the distance the entries of ${{\mathcal{A}}}$ need to be perturbed to have an attainable or unattainable non-trivial Smith form is to solve

[TABLE]

Here ${\bf{d}}^{\prime}$ is the vector of the largest possible degrees of each entry of $\operatorname{Adj}({{\mathcal{A}}}+\widetilde{\Delta}{{\mathcal{A}}})$, and $\widetilde{\Delta}{{\mathcal{A}}})$ is a prescribed linear or affine perturbation structure.

It is sufficient to compute $\max({\bf{d}}^{\prime})$, a quantity which will generically be $(n-1)d$. For non-generic instances we require the computation of ${\bf{d}}^{\prime}$. This optimization problem is non-convex, but multi-linear in each coefficient of $\Delta{{\mathcal{A}}}$.

We do not attempt to solve this problem directly via numerical techniques, since it enforces a necessary condition that is often sufficient. Instead we use it to develop a theory of solutions which can be exploited by faster and more robust numerical methods.

Lemma 3.7**.**

Let ${\bf{f}}$ be a vector of polynomials with degrees ${\bf{d}}$ and admissible perturbations $\Delta{\bf{f}}$ of degrees ${\bf{d}}^{\prime}$ where $\max({\bf{d}})\leq\max({\bf{d}}^{\prime})$. Then the family of generalized Sylvester matrices $\operatorname{\mbox{{Syl}}}_{{\bf{d}}^{\prime}}({\bf{f}})$ of rank at least $e$ form an open set under the perturbations $\Delta{\bf{f}}$.

Proof.

By the degree assumption on $\Delta{\bf{f}}$ we have that for an infinitesimal $\Delta{\bf{f}}$ that $\operatorname{\mbox{{Syl}}}_{{\bf{d}}^{\prime}}({\bf{f}})$ and $\operatorname{\mbox{{Syl}}}_{{\bf{d}}^{\prime}}(\Delta{\bf{f}})$ have the same dimension. Accordingly, let us suppose that $\operatorname{\mbox{{Syl}}}_{{\bf{d}}^{\prime}}({\bf{f}})$ has rank at least $e$. Then the Sylvester matrix in question must have rank at least $e$ in an open-neighborhood around it. In particular, when ${\|\operatorname{\mbox{{Syl}}}_{{\bf{d}}^{\prime}}(\Delta{\bf{f}})\|}_{2}<\sigma_{e}(\operatorname{\mbox{{Syl}}}_{{\bf{d}}^{\prime}}({\bf{f}}))$ then $\operatorname{rank}(\operatorname{\mbox{{Syl}}}_{{\bf{d}}^{\prime}}({\bf{f}}+\Delta{\bf{f}}))\geq\operatorname{rank}(\operatorname{\mbox{{Syl}}}_{{\bf{d}}^{\prime}}({\bf{f}}))$ and the result follows. ∎

Theorem 3.8**.**

The optimization problem (6) has an attainable global minimum under linear perturbation structures.

Proof.

Let $\mathcal{S}$ be the set of all rank at most $e-1$ generalized Sylvester matrices of prescribed shape by ${\bf{d}}^{\prime}$ and $\operatorname{Adj}({{\mathcal{A}}})$. Lemma 3.7 implies that $\mathcal{S}$ is topologically closed.

Let $\mathcal{R}=\{\operatorname{\mbox{{Syl}}}_{{\bf{d}}^{\prime}}(\operatorname{Adj}({\mathcal{C}}))\text{ subject to }{\|{\mathcal{C}}\|}\leq{\|{{\mathcal{A}}}\|}\}$, where the generalized Sylvester matrices are padded with zeros to have the appropriate dimension if required. Since $\Delta{{\mathcal{A}}}$ has a linear perturbation structure, a feasible point is always ${C}=-{{\mathcal{A}}}$. By inspection $\mathcal{R}$ is seen to be a non-empty set that is bounded and closed.

The functional ${\|\cdot\|}$ is continuous over the non-empty closed and bounded set $\mathcal{S}\cap\mathcal{R}$. Let ${{\mathcal{B}}}\in\mathcal{S}\cap\mathcal{R}$. By Weierstrass’s theorem ${\|{{\mathcal{A}}}-{\mathcal{B}}\|}$ has an attainable global minimum over $\mathcal{S}\cap\mathcal{R}$. ∎

Note that if a feasible point exists under an affine perturbation structure, then a solution to the optimization problem exists as well. What this result says is that computing the distance to non-triviality is generally a well-posed problem, even though computing a matrix polynomial of minimum distance may be ill-posed (the solution is unattainable). The same results also hold when working over the ${\bf{d}}^{\prime}$ reversed coefficients. A similar argument is employed by (Kaltofen et al., 2007, Theorem 2.1).

3.3 Bounds on the Distance to non-triviality

Suppose that ${{\mathcal{A}}}\in{\mathbb{R}}[t]^{n\times n}$, of degree at most $d$, has a trivial Smith form and does not have an unattainable non-trivial Smith form. This section provides some basic bounds on the distance coefficients of ${{\mathcal{A}}}$ need to be perturbed to have a non-trivial Smith form. The bounds we derive are unstructured, although they can be generalized to several perturbation structures (such as ones that preserve the degree or support of entries) in a straight forward manner.

If we consider the mapping $\operatorname{Adj}(\cdot)$ as a vector-valued function from ${\mathbb{R}}^{n^{2}(d+1)}\to{\mathbb{R}}^{n^{2}((n-1)d+1)}$ (with some coordinates possibly fixed to zero), then we note that the mapping is locally Lipschitz. More precisely, there exists $c>0$ such that for a sufficiently small $\Delta{{\mathcal{A}}}$,

[TABLE]

The quantity $c$ can be approximately bounded above by the (scalar) Jacobian matrix $\nabla\operatorname{Adj}(\cdot)$ evaluated at ${{\mathcal{A}}}$. A local upper bound for $c$ is approximately ${\|\nabla\operatorname{Adj}({{\mathcal{A}}})\|}_{2}$. We can invoke Theorem 2.13 if ${{\mathcal{A}}}$ has full rank. By considering $\hat{c}={\left\|\left[\Phi_{(n-1)d}(I\otimes{{\mathcal{A}}})\right]^{+}\right\|}_{2}$, we obtain the (absolute) first-order approximate perturbation bound

[TABLE]

The entries of $\nabla\operatorname{Adj}({{\mathcal{A}}})$ consist of the coefficients of the $(n-2)\times(n-2)$ minors of ${{\mathcal{A}}}$. This follows because $\operatorname{Adj}(\cdot)$ is a multi-linear vector mapping and the derivative of each entry is a coefficient of the leading coefficient with respect to the variable of differentiation. The size of each minor can be bounded above (albeit poorly) by Hadamard’s inequality (Goldstein-Graham variant, see (Lossers, 1974)). As such, we have the sequence of bounds

[TABLE]

where ${\|{{\mathcal{A}}}\|}_{\infty}$ is understood to be a vector norm and ${\|\nabla\operatorname{Adj}({{\mathcal{A}}})\|}_{\infty}$ is understood to be a matrix norm. The bound in question can be used in conjunction with the SVD to obtain a lower bound on the distance to a matrix polynomial with a non-trivial Smith form.

Theorem 3.9**.**

Suppose that ${\bf{d}}^{\prime}=(\gamma,\gamma\ldots,\gamma)$ and $\operatorname{\mbox{{Syl}}}_{{\bf{d}}^{\prime}}(\operatorname{Adj}({{\mathcal{A}}}))$ has rank $e$. Then an approximate lower bound on the distance to non-triviality is

[TABLE]

Proof.

We note that for polynomials ${\bf{f}}$ with degrees ${\bf{d}}^{\prime}$ that ${\|\operatorname{\mbox{{Syl}}}_{{\bf{d}}^{\prime}}({\bf{f}})\|}_{F}=\gamma{\|{\bf{f}}\|}_{F}$. Accordingly, if $\Delta{{\mathcal{A}}}$ is a minimal perturbation to non-triviality, then

[TABLE]

and the theorem follows by a simple rearrangement. Note that ${\|\cdot\|}_{2}\leq{\|\cdot\|}_{F}$. ∎

If ${\bf{d}}^{\prime}$ has different entries, then $\ell{\|{\bf{f}}\|}_{F}\leq{\|\operatorname{\mbox{{Syl}}}_{{\bf{d}}^{\prime}}({\bf{f}})\|}_{F}\leq\gamma{\|{\bf{f}}\|}_{F},$ where $\gamma$ and $\ell$ are the largest and second-largest entries of ${\bf{d}}^{\prime}$. The lower bound provided can also be improved using the Karmarkar-Lakshman distance (Karmarkar and Lakshman, 1996) in lieu of the smallest singular value of the generalized Sylvester matrix, the ${\bf{d}}^{\prime}$ reversal of the adjoint or other approximate GCD lower bounds (e.g., (Beckermann and Labahn, 1998)).

4 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.

4.1 Constrained Optimization Formulation

An equivalent statement to ${{\mathcal{A}}}$ having a non-trivial attainable Smith form is that $\operatorname{Adj}({{\mathcal{A}}})={\mathcal{F}}^{*}h$ where ${\mathcal{F}}^{*}$ is a vector (or matrix) of scalar polynomials and $h$ is a divisor of $\gcd(\operatorname{Adj}({{\mathcal{A}}}))$. This directly leads to the following optimization problem:

[TABLE]

This is a multi-linearly structured approximate GCD problem which is a non-convex optimization problem. Instead of finding a rank deficient Sylvester matrix, we directly enforce that the entries of $\operatorname{Adj}({\mathcal{A}})$ have a non-trivial GCD. The normalization requirement that $\mathcal{N}_{h}{\mbox{vec}}(h)=1$ is chosen to force $h$ to have a non-zero degree, so that $h$ is not a scalar. One useful normalization is to define $\mathcal{N}_{h}$ such that $\operatorname{{lcoeff}}(h)=1$ (that is $\operatorname{{lcoeff}}(\cdot)$ is the leading coefficient of a polynomial). Explicitly, we assume the degree of the approximate GCD is known and make it monic. Of course, other valid normalizations also exist.

Since we are working over ${\mathbb{R}}[t]$, there will always be a quadratic, linear or zero factor of attainable solutions. If $h=0$ then the approximate SNF of ${\mathcal{A}}$ is rank deficient and computing approximate SNF reduces to the nearest rank at-most $n-1$ or $n-2$ matrix polynomial problems, both of which are well-understood (Giesbrecht, Haraldson, and Labahn, 2017, 2019b). Assuming that we are now working in the nonzero case, we can assume generically that $\deg(h)=1$ or $\deg(h)=2$.

4.2 Lagrange Multipliers and Optimality Conditions

In order to solve our problem we will employ the method of Lagrange multipliers. The Lagrangian is defined as

[TABLE]

where $\lambda$ is a vector of Lagrange multipliers.

A necessary first-order condition (KKT condition, e.g. (Bertsekas, 1999)) for a tuple $z^{\star}=z^{\star}(\Delta{\mathcal{A}},{\mathcal{F}}^{*},h,\lambda)$ to be a regular (attainable) minimizer is that the gradient of $L$ vanishes, that is,

[TABLE]

Let $J$ be the Jacobian matrix of the constraints defined as

[TABLE]

The second-order sufficiency condition for optimality at a local minimizer $z^{\star}$ is that

[TABLE]

that is, the Hessian with respect to $x=x(\Delta{{\mathcal{A}}},{\mathcal{F}}^{*},h)$ is positive definite over the kernel of the Jacobian of the constraints. The vector $x$ corresponds to the variables in the affine structure of $\Delta{{\mathcal{A}}}$,${\mathcal{F}}^{*}$, and $h$. If (8) and (9) both hold, then $z^{\star}$ is necessarily a local minimizer of (7). Of course, it is also necessary that $\ker(J(z^{\star}))^{T}\nabla^{2}_{xx}L(z^{\star})\ker(J(z^{\star}))\succeq 0$ at a minimizer, which is the second-order necessary condition. Our strategy for computing a local solution is to solve $\nabla L=0$ using a Newton-like method.

4.3 An Implementation with Local Quadratic

Convergence

A problem with Newton-like methods is that when the Hessian is rank deficient or ill-conditioned, then the Newton step becomes ill-defined or the rate of convergence degrades. The proposed formulation of our problem can encounter a rank deficient Hessian (this is due to over padding some vectors with zero entries or redundant constraints). Despite this we are still able to obtain a method with rapid local convergence under a very weak normalization assumption.

In order to obtain rapid convergence we make use of the Levenberg-Marquart (LM) algorithm. If $H=\nabla^{2}L$, then the LM iteration is defined as repeatedly solving for $z^{(k+1)}=z^{(k)}+\Delta z^{(k)}$ by

[TABLE]

for some $\ell>0$ while using ${\|\nabla L\|}_{2}$ as a merit function. The speed of convergence depends on the choice of $\nu_{k}>0$. Note that since LM is essentially a regularized Gauss-Newton method, when the Hessian is rank deficient then we may converge to a stationary point of the merit function. If convergence to a stationary point of the merit function is detected, then the method of Wright (2005) can be used to replace LM in several instances.

Yamashita and Fukushima (2001) show that, under a local-error bound condition, a system of non-linear equations $g(z)=0$ approximated by LM will converge quadratically to a solution with a suitable initial guess. Essentially, what this says is that to obtain rapid convergence it is sufficient for regularity ($J$ having full rank) to hold or second-order sufficiency, but it is not necessary to satisfy both. Note that we assume Lagrange multipliers exist. However, unlike the case when $J$ has full rank, the multipliers need not be unique. The advantage of LM over other Newton-like methods is that this method is globalized111Here “globalized” means that the method will converge to a stationary point of the merit function, not a local extremum of the problem. in exchange for an extra matrix multiplication, as $H^{T}H+\nu_{k}I$ is always positive definite, and hence always a descent direction for the merit function. We make the choice of $\nu_{k}\approx{\|g(z)\|}_{2}$ based on the results of Fan and Yuan (2005).

Definition 4.1** (Local Error Bound).**

Let $Z^{\star}$ be the set of all solutions to $g(z)=0$ and $X$ be a subset of ${\mathbb{R}}^{\ell}$ such that $X\cap Z^{\star}\neq\emptyset$. We say that ${\|g(z)\|}$ provides a local error bound on $g(z)=0$ if there exists a positive constant $c$ such that $c\cdot\mbox{{dist}}(z,Z^{\star})\leq{\|g(z)\|}$ for all $z\in X$, where $\mbox{{dist}}(\cdot)$ is the distance between a point and a set.

In this section it is useful to consider $g(z)=\nabla L(z)$, as we need the local error bounds to estimate $\nabla L(z)=0$.

Theorem 4.2**.**

If the second-order sufficiency condition (9) holds at an attainable solution to (7), then the local error-bound property holds.

Proof.

This result follows immediately from Section 3 of Wright (2005) and the references therein. ∎

The bounds of Wright can be used to infer when quadratic convergence occurs for Newton-like methods. In this problem, perturbations to $x$ are important in understanding how the problem behaves locally.

Remark 4.3**.**

Let $z=z(x,\lambda)$ where $x$ is a vector of variables and $\lambda$ is a vector of Lagrange multipliers, and define $g(z)=\nabla L(z)$. First suppose that both the second-order sufficiency condition (9) and first-order necessary condition (8) hold at the point $z^{\star}$. We can write the first-order expansion

[TABLE]

noting that $g(z^{\star})=0$. It is useful to observe that

[TABLE]

If $\Delta x=0$ then the error-bound from Hoffman (1952) (main theorem) applies and we have that there exists $c_{hof}>0$ such that $c_{hof}{\|\Delta\lambda\|}\leq{\|g(x,\lambda+\Delta\lambda)\|}$. If $\Delta x\neq 0$ then ${\left\|\begin{pmatrix}H_{xx}(z^{\star})\\ J(z^{\star})\end{pmatrix}\Delta x\right\|}\approx{\|g(x+\Delta x,\lambda)\|}$ and (9) implies that $H(z^{\star})(\Delta z)=0\implies\Delta x=0$, so

[TABLE]

so there exists $c_{\sigma_{\min}}>0$ when ${\|\Delta x\|}$ is sufficiently small such that $c_{\sigma_{\min}}{\|\Delta x\|}\leq{\|g(x+\Delta x,\lambda)\|}$. Note that $c_{\sigma_{\min}}\approx\sigma_{\min}\begin{pmatrix}H_{xx}(z^{\star})\\ J(z^{\star})\end{pmatrix}.$

The first-order approximation implies that when ${\|\Delta z\|}$ is sufficiently small that

[TABLE]

The key idea is to separate the problem into the cases of $\Delta x=0$ and $\Delta x\neq 0$, and then derive error bounds for each case. The important part of the discussion is that if one can estimate $c_{\sigma_{\min}}$ then one can often infer when quadratic convergence occurs.

The second-order sufficiency assumption is not necessary to derive error bounds bounds. It is straightforward to show the local error bound property holds if $J(z^{\star})$ has full rank, as the Lagrange multipliers will be (locally) unique, hence the solution is (locally) unique. Alternatively, if $J$ had constant rank in a non-trivial open neighborhood around a solution, then a similar argument could be made about the local error-bound property.

Theorem 4.4**.**

The second-order sufficiency condition holds at minimal solutions with Lagrange multipliers of minimal norm if $h$ is of maximal degree and monic and the minimal structured perturbation ${\|\Delta{{\mathcal{A}}}^{\star}\|}$ is sufficiently small.

Proof.

The Hessian of $L$ with respect to $x=x(\Delta A,{\mathcal{F}}^{*},h)$ is

[TABLE]

where $F$ is a square matrix with zero diagonal whose entries are a multi-linear polynomial in $\lambda$ and $\Delta{{\mathcal{A}}}$ and $E^{T}$ is a matrix whose entries are homogeneous linear functions in $\lambda$.

If $\Delta{{\mathcal{A}}}^{\star}=0$ then $\lambda^{\star}=0$. Hence both $E=0$ and $F=0$ and so, if $y\in\ker(H_{xx})\cap\ker(J)$ then $y=\begin{pmatrix}0&y_{2}&y_{3}\end{pmatrix}^{T}$. The Jacobian of the constraints may be written (up to permutation) as

[TABLE]

where $*$ are blocks corresponding to differentiating with respect to variables in $\Delta{{\mathcal{A}}}$ and the blocks $\mathcal{C}_{{\mathcal{F}}^{*}}$ and $\mathcal{C}_{h}$ consist of block convolution and convolution matrices that correspond to multiplication by ${\mathcal{F}}^{*}$ and $h$, respectively. The block $\mathcal{N}_{h}$ contains a normalization vector to ensure that $h$ has the appropriate degree. $Jy=0$ implies that there exists a vector of polynomials $v$ and a polynomial $u$ with the same degrees as ${\mathcal{F}}^{*}$ and $h$ such that ${\mathcal{F}}^{*}u+vh=0$ and $\mathcal{N}_{h}{\mbox{vec}}(u)=0$.

We have that $h$ is a factor of both ${\mathcal{F}}^{*}u$ and $vh$. Since $\gcd({\mathcal{F}}^{*},h)=1$ it must be that $h$ is a factor of $u$. It follows that $\deg(u)=\deg(h)$, so there exists some $\alpha\neq 0$ such that $\alpha u=h$. Since $h$ is monic, we have that $\mathcal{N}_{h}{\mbox{vec}}(h)=1$ but $\mathcal{N}_{h}{\mbox{vec}}(u)=0$, which implies that $\alpha=0$, and so $u=0$. We have that $vh=0$ and this implies $v=0$. Hence $\ker(J)\cap\ker(H_{xx})=0$ and second-order sufficiency holds when ${\|\Delta{{\mathcal{A}}}^{*}\|}=0$.

If ${\|\Delta{{\mathcal{A}}}^{*}\|}$ is sufficiently small, then ${\|F\|}$ will be sufficiently small so that $F+2I$ has full rank. Accordingly, we have that

[TABLE]

We remark that the techniques in the proof are very similar to those of Zeng and Dayton (2004) and Giesbrecht, Haraldson, and Kaltofen (2019a) to show that a Jacobian matrix appearing in approximate GCD computations of two (or more) polynomials has full rank. If we over-estimated the degrees of ${\mathcal{F}}^{*}$ then $H_{xx}$ would have some columns and rows consisting of zero (the block-convolution matrices would be padded with extra zero entries).

In the proof of Theorem 4.4 we note that

[TABLE]

The matrix $F=\nabla_{\Delta{{\mathcal{A}}}}\lambda^{T}J_{\operatorname{Adj}}({{\mathcal{A}}}+\Delta{{\mathcal{A}}})$ will consist of coefficients of the $(n-3)\times(n-3)$ minors of ${{\mathcal{A}}}+\Delta{{\mathcal{A}}}$ scaled by entries of $\lambda$. Accordingly, $F$ will generally not have $-2$ as an eigenvalue.

Remark 4.5**.**

Thus far we have assumed that Lagrange multipliers exist at the current solutions of interest, which are attainable solutions that have full rank. Corollary 2.14 and the proof of Theorem 4.4 imply that Lagrange multipliers generally exist under these assumptions for several perturbation structures, since we need to solve

[TABLE]

of which $J$ generally has constant or full rank. Of course if the solution was unattainable then the GCD constraints would break down as there is a “solution at infinity” in a sense that ${\|h\|}\to\infty$ as $\Delta{{\mathcal{A}}}\to\Delta{{\mathcal{A}}}^{\star}$.

The implication of the local-error bound property holding is that one can reasonably approximate when quadratic convergence occurs by estimating $\sigma_{\min}\left(\left[H_{xx}~{}|~{}J^{T}\right]\right)$ and $c_{hof}$. In particular, these quantities act as a structured condition number on the system. A structured backwards-error analysis of existing techniques can be performed using these quantities. Additionally, it is somewhat generic that $F+2I$ has full rank, hence the local error-bound will hold for most instances of the approximate SNF problem with an attainable solution. It is also important to note that we did not explicitly use the adjoint matrix. Indeed the result remains valid if we replace the adjoint with minors of prescribed dimension. Likewise, if ${{\mathcal{A}}}$ is an ill-posed instance of lower McCoy rank or approximate SNF without an attainable global minimum, then optimizing over a reversal of each entry of $\operatorname{Adj}({{\mathcal{A}}}+\Delta{{\mathcal{A}}})$ would yield a non-trivial answer and the same stability properties would hold. Thus, poorly posed problems also remain poorly posed if slightly perturbed.

Corollary 4.6**.**

The LM algorithm for solving $\nabla L=0$ has quadratic convergence under the assumptions of Theorem 4.2 and using $\nu_{k}={\|\nabla(L(z^{k}))\|}_{2}$.

Proof.

The quantity $\nabla L$ is a multivariate polynomial, hence it is locally Lipschitz. Second-order sufficiency holds, thus we have the local error bound property is satisfied. The method converges rapidly with a suitable initial guess. ∎

Note that for several perturbation structures if the adjoint has generic degrees, then the Jacobian of the constraints will have full rank, and a standard Newton iteration is also well-defined, and will converge quadratically as well.

In the next section we discuss a technique that possibly forgoes rapid local convergence, but has a polynomial per iteration cost to compute a low McCoy rank approximation.

4.4 Computational Challenges and Initial Guesses

The most glaring problem in deriving a fast iterative algorithm for the approximate Smith form problem is that the matrix $\operatorname{Adj}({\mathcal{A}}+\Delta{\mathcal{A}})$ has exponentially many coefficients as a multivariate polynomial in $\Delta{\mathcal{A}}$. This means computing the adjoint matrix symbolically as an ansatz is not feasible. In order to solve (8) we instead approximate the derivatives of the coefficients of the adjoint numerically.

To compute an initial guess, we can use $\Delta{{\mathcal{A}}}_{init}=0$ and take ${\mathcal{F}}^{*}$ and $h$ to be a reasonable approximation to an approximate GCD of $\operatorname{Adj}({{\mathcal{A}}})$, which will often be valid as per Theorem 4.2. To make sure the point is feasible, one can use a variant of Newton’s method to project to a feasible point. Corollary 2.14 implies that with a suitable initial guess, reasonable variants of Newton’s method (such as LM) will converge quadratically to a feasible point, assuming one exists.

Another technique is to take two rows or columns of ${{\mathcal{A}}}$ and perturb them so that the $2n$ entries have a non-trivial GCD. To find the best guess with this technique, $O(n^{2})$ approximate GCD computations on $O(n)$ polynomials of degree $d$ need to be performed. In the next section we will discuss more sophisticated techniques.

4.5 Attaining Unattainable Solutions

If a solution is unattainable then the degrees of all the entries of the adjoint matrix may change in an open neighborhood around a solution. If $\Delta{{\mathcal{A}}}^{\star}$ is an unattainable solution (of full rank) to (7) then $h(t)=t$ is clearly not a solution since $h(t)=t$ being a solution implies that such a solution would be attainable. Let $d_{\operatorname{Adj}}$ be the generic degree of $\operatorname{Adj}({{\mathcal{A}}}+\Delta{{\mathcal{A}}})$, then $t$ is a factor of $\gcd(\mbox{{rev}}_{d_{\operatorname{Adj}}}(\operatorname{Adj}({{\mathcal{A}}}+\Delta{{\mathcal{A}}}^{\star})))$. The reversed adjoint has no GCD at infinity by assumption, as such a GCD at infinity would be an attainable solution to the original problem. Accordingly, we note that Theorem 4.4 applies after some straightforward modifications, since

[TABLE]

are essentially (block) permutations of each other.

Since $\mbox{{rev}}_{d_{\operatorname{Adj}}}(\operatorname{Adj}({{\mathcal{A}}}+\Delta{{\mathcal{A}}})))$ achieves the generic degree, Lagrange multipliers should exist as we can apply Corollary 2.14 on $\nabla{\mbox{vec}}(\mbox{{rev}}_{d_{\operatorname{Adj}}}(\operatorname{Adj}({{\mathcal{A}}}+\Delta{{\mathcal{A}}})))$ by permuting entries, and the underlying approximate GCD problem is well-posed. Thus the problem will also typically admit Lagrange multipliers.

The essential ingredient in Theorem 4.4 is the normalization of the underlying approximate GCD problem. This means that “backwards stable” algorithms will compute the exact SNF of a nearby matrix polynomial that has no meaning in the context of computation. This generally occurs because the radius of uncertainty, usually proportional to unit rounding errors, contains infinitely many matrix polynomials with a non-trivial SNF. The backwards stability is not meaningful in this context, because the instance of the problem is not continuous. In such instances, computing the SNF is most likely the wrong

problem to be considering. Instead, computing the spectral structure of eigenvalues at infinity is most likely the appropriate problem. However there exist instances where both problems could be simultaneously poorly conditioned.

If the reversed problem has a radius of stability with respect to Theorem 4.4, then the original problem has a radius of instability, meaning that the iterates will converge to a point where ${\|h\|}$ is excessively large. In other words, if an instance of a problem is ill-posed, then it cannot be regularized — the finite and infinite eigenvalues and their spectral structure is indistinguishable in floating point arithmetic — in the context of the QZ decomposition, GUPTRI (Demmel and Kågström, 1993a, b) or similar algorithms. There are some instances where attempting to compute the SNF numerically is not possible and should not be attempted. In the context of an optimization problem, we can of course regularize the problem as we have just described. Van Dooren (1979) suggests that ill-posed problems should be formulated as an optimization problem as a means of regularization to overcome some of the numerical difficulties.

5 Lower McCoy Rank Approximation

In this section we describe how to find a nearby matrix polynomial of lower McCoy. Another way to formulate ${{\mathcal{A}}}$ having a non-trivial SNF is to solve the minimization problem

[TABLE]

where $\Delta{{\mathcal{A}}}$ must have the appropriate structure. Essentially this finds the smallest perturbation of ${{\mathcal{A}}}$ with an eigenvalue that lowers the rank by at least $2$. The auxiliary variables $\omega$ and $B$ are used to enforce this constraint. Here $B^{*}$ is the conjugate transpose of $B$, and $B^{*}B=I_{2}$ ensures that the kernel vectors are linearly independent and do not tend towards zero.

The optimization is unstable if $\omega$ is reasonably large, since the largest terms appearing are proportional to $O((d+1){\|{{\mathcal{A}}}\|}_{\infty}|\omega|^{d})$. To remedy this, if we assume that a solution to the optimization problem (10) exists and has full rank, then we may transform ${{\mathcal{A}}}+\Delta{{\mathcal{A}}}$ into a degree-one matrix polynomial (also known as a matrix pencil) with the same spectral properties, known as a linearization. If there is no full-rank solution one can simply take a lower-rank approximation (Giesbrecht, Haraldson, and Labahn, 2019b) and extract a square matrix polynomial of full rank that may be linearized. Alternatively, one may forgo the linearization and work directly with a problem that is more poorly conditioned. For the rest of this section we will assume, without loss of generality, that ${{\mathcal{A}}}$ and the solutions to the low McCoy rank problem have full rank.

We can encode the spectral structure and SNF of ${{\mathcal{A}}}$ as the following degree-one matrix polynomial (sometimes referred to as the companion linearization (Gohberg et al., 2009)) of the form ${\mathcal{P}}\in{\mathbb{R}}[t]^{nd\times nd}$, defined as

[TABLE]

This particular linearization encodes the SNF of ${{\mathcal{A}}}$, as $\operatorname{SNF}({\mathcal{P}})=\operatorname{diag}(I,I,\ldots,I,\operatorname{SNF}({{\mathcal{A}}}))$. It follows that ${{\mathcal{A}}}$ has a non-trivial SNF if and only if ${\mathcal{P}}$ has a non-trivial SNF. If we preserve the affine structure of ${\mathcal{P}}$ and only perturb blocks corresponding to ${{\mathcal{A}}}$, then the reduction to a pencil will be sufficient. Other linearizations are possible as well. The pencil is generally better behaved numerically since the largest entry upon evaluation at a $\omega\in{\mathbb{C}}$ is proportional to $O(d{\|{{\mathcal{A}}}\|}_{\infty}|\omega|)$ rather than $O({\|A\|}_{\infty}|\omega|^{d})$, albeit with matrices that are $d$ times larger.

5.1 Fast Low McCoy Rank via Optimization

One way to approach the lower McCoy rank approximation problem is to study all the minors (or sufficiently many) of a matrix polynomial. This method immediately generalizes from the previous section, however is not practical for computational purposes since the number of minors grows exponentially in the dimension. Instead, we can approach the problem by formulating it as an optimization problem, one that is remarkably similar to structured lower rank approximation of scalar matrices. This similarity facilitates computing an initial guess for the following optimization problem using the SVD.

The lower McCoy rank approximation problem may be formulated as the following real optimization problem: to find the nearest matrix polynomial to ${{\mathcal{A}}}\in{\mathbb{R}}[t]^{n\times n}$ with McCoy rank $n-r$, find the perturbation $\Delta{{\mathcal{A}}}\in{\mathbb{R}}[t]^{n\times n}$ which minimizes

[TABLE]

Note that the perturbation $\Delta{{\mathcal{A}}}$ is real valued in this problem. The unitary constraint on $B$ ensures that $\operatorname{rank}(B)=r$ and each column of $B$ remains away from zero. Accordingly, $\omega\in{\mathbb{C}}$ will be an eigenvalue of $({{\mathcal{P}}}+\Delta{{\mathcal{P}}})(\omega)$ since $\operatorname{rank}(({{\mathcal{P}}}+\Delta{{\mathcal{P}}})(\omega))\leq nd-r$, and thus the McCoy rank of ${{\mathcal{A}}}+\Delta{{\mathcal{A}}}$ is at-most $n-r$.

Real matrix polynomials can have complex eigenvalues and so complex numbers must necessarily appear in the constraints. The constraints arising from the complex numbers may be divided into real parts and imaginary parts, denoted as $\Re(\cdot)$ and $\Im(\cdot)$, respectively. By dividing the constraint into real and imaginary parts, we are able to solve an equivalent optimization problem completely with real variables. This ensures that $\Im(\Delta{{\mathcal{A}}})=0$, that is, the perturbations are real. Since ${{\mathcal{A}}}+\Delta{{\mathcal{A}}}$ may have complex eigenvalues (but entries with real coefficients), we require that $\operatorname{SNF}({{\mathcal{A}}}+\Delta{{\mathcal{A}}})$ has entries from ${\mathbb{R}}[t]$. Accordingly, we need to interpret the auxiliary variable $\omega$. The instance of $\Im(\omega)=0$ corresponds to $t-\omega$ as an invariant factor, while $\Im(\omega)\neq 0$ corresponds to the real irreducible quadratic $(t-\omega)(t-\overline{\omega})$. Thus at a solution, we are able to recover a real invariant factor regardless if $\omega$ has a non-zero imaginary part.

In order to approach the problem using the method of Lagrange multipliers we define the Lagrangian as

[TABLE]

and proceed to solve $\nabla L=0$. In our implementation we again make use of the LM method, although given the relatively cheap gradient cost, a first-order method will often be sufficient and faster. The problem is essentially tri-linear, and structurally similar to affinely structured low rank approximation, of which Lagrange multipliers will exist for most instances.

It is important to note that an attainable solution to this problem is not guaranteed, as it is possible for ${\|\omega\|}\to\infty$ as $\Delta{{\mathcal{A}}}\to\Delta{{\mathcal{A}}}^{\star}$. Such an instance is an unattainable solution in the context of Section 4.5. These solutions behave like an infinite eigenvalue and can be handed by specifically considering the eigenvalue $t=0$ of the reversed matrix polynomial.

5.2 Computing an Initial Guess

In order to compute an initial guess to (11) we exploit the pseudo tri-linearity of the problem. If two of $\Delta{{\mathcal{A}}}$, $\omega$ and $B$ are fixed then the problem is linear (or a linear surrogate can be solved) in the other variable. Despite the unitary constraint on $B$ being non-linear, it is not challenging to handle. Any full rank $B$ is suitable for an initial guess, since we may orthonormalize $B$ to satisfy the constraint that $B^{*}B=I_{r}$.

First we approximate the determinant of ${{\mathcal{A}}}$ and consider initial guesses where $\sigma_{n-r}({{\mathcal{A}}}(\omega^{init}))$ is reasonably small. If $\sigma_{n-r}({{\mathcal{A}}}(\omega^{init}))$ is reasonably small, then $\omega^{init}$ is (approximately) an eigenvalue of a nearby matrix polynomial of reduced McCoy rank. The zeros and local extrema of $\det({{\mathcal{A}}})$ are suitable candidates for computing an initial guess for $\omega$. The kernel $B^{init}$ can be approximated from the smallest $r$ singular vectors of ${{\mathcal{A}}}(\omega^{init})$. This ensures that $B^{init}$ is unitary and spans the kernel of a nearby rank deficient (scalar) matrix.

To compute an initial guess for $\Delta{{\mathcal{A}}}$ we can take $\Delta{{\mathcal{A}}}^{init}=0$, or solve a linear least squares problem where $B$ and $\omega$ are fixed. Alternatively, one may project to a feasible point by using a variant of Newton’s method, using $\Delta{{\mathcal{A}}}^{init}=0$, $\omega^{init}$ and $B^{init}$ as an initial guess for the Newton iteration to solve $({{\mathcal{A}}}+\Delta{{\mathcal{A}}})(\omega)B=0$ and $B^{*}B=I_{r}$. A feasible point computed by Newton’s method tends not to perturb $\Delta{{\mathcal{A}}}$ very much, whereas the least squares approximation may perturb ${{\mathcal{A}}}$ by an unnecessarily large amount.

5.3 About Global Optimization Methods

The problems previously discussed are NP hard to solve exactly and to approximate with coefficients from ${\mathbb{Q}}$. This follows since affinely structured low rank approximation (Braatz et al., 1994; Poljak and Rohn, 1993) is a special case. If we consider a matrix polynomial of degree zero, then this is a scalar matrix with an affine structure. The approximate SNF will be a matrix of rank at most $n-2$, and finding the nearest affinely structured singular matrix is NP hard.

Despite the problem being intractable in the worst case, not all instances are necessarily hard. The formulation (11) is multi-linear and polynomial, hence amenable to the sum of squares hierarchy. Lasserre’s sum of squares hierarchy (Lasserre, 2001) is a global framework for polynomial optimization that asymptotically approximates a lower bound. Accordingly, if ${\|\omega^{opt}\|}$ is bounded, then sum of squares techniques should yield insight into the problem.

6 Implementation and Examples

We have implemented our algorithms and techniques in the Maple computer algebra system222Sample code is at https://www.scg.uwaterloo.ca/software/GHL2018jsc-code-2018-11-28.tgz. . 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 ${{\mathcal{A}}}$, which are $n$ and $d$ respectively. The cost of most quasi-Newton methods is roughly proportional to inverting the Hessian matrix, which is $O(\ell^{3}),$ where $\ell$ is the number of variables in the problem.

The method of Section 4 requires approximately $O((n^{3}d)^{3})=O(n^{9}d^{3})$ FLOPs per iteration in an asymptotically optimal implementation with cubic matrix inversion, which is the cost of inverting the Hessian. Computing the Hessian costs roughly ${O\tilde{\phantom{\imath}}}(n^{4}d^{2}\times(n^{2})^{2})={O\tilde{\phantom{\imath}}}(n^{8}d^{2})$ FLOPs333For $\phi,\psi\colon\mathbb{R}\to\mathbb{R}$, $\phi={O\tilde{\phantom{\imath}}}(\psi)$ iff $\phi=O(\psi(\log|\psi|)^{c})$ for some absolute constant $c\geq 0$, i.e., we ignore log factors. using a blocking procedure, assuming the adjoint computation runs in ${O\tilde{\phantom{\imath}}}(n^{4}d)$ FLOPs (which can be done via interpolation in a straightforward manner)There are $O(n^{3}d)$ Lagrange multipliers since the adjoint has degree at most $(n-1)d$. Using reverse-mode automatic differentiation to compute $\nabla^{2}L$, this can be accomplished in ${O\tilde{\phantom{\imath}}}(n^{4}d\times n^{3}d)={O\tilde{\phantom{\imath}}}(n^{7}d^{2})$ FLOPs.

The method of Section 5 has a Hessian matrix of size $O(n^{2}d^{2})\times O(n^{2}d^{2})$ in the case of a rank zero McCoy rank approximation. Accordingly, the per iteration cost is roughly $O(n^{6}d^{6})$ FLOPs. If the linearization is not performed, then the per-iteration cost is $O(n^{6}d^{3})$ FLOPs. Given the lack of expensive adjoint computation, a first-order method will typically require several orders of magnitude fewer FLOPs per iteration (ignoring the initial setup cost), with local linear convergence.

Example 6.1** (Nearest Interesting SNF).**

Consider the matrix polynomial ${{\mathcal{A}}}$ with a trivial SNF

[TABLE]

of the form $\operatorname{diag}(1,\ldots,1,\det({{\mathcal{A}}}))$.

If we prescribe the perturbations to leave zero coefficients unchanged, then using the methods of Section 4 and Section 5 results in a local minimizer ${{\mathcal{A}}}+\Delta{{\mathcal{A}}}_{opt}$ given by

[TABLE]

with ${\|\Delta{{\mathcal{A}}}_{opt}\|}\approx.164813183138322$. The SNF of ${{\mathcal{A}}}+\Delta{{\mathcal{A}}}_{opt}$ is approximately

[TABLE]

where $s_{1}\approx t^{2}+0.0632934647739423t+0.960572576466186$. The factor $s_{1}$ corresponds to $\omega_{opt}\approx-0.0316467323869714-0.979576980535687i$.

The method discussed in Section 4 converges to approximately $14$ decimal points of accuracy444$\nabla L=0$ is solved to $14$ digits of accuracy; the extracted quantities are accurate to approximately the same amount. after $69$ iterations and the method of Section 5 converges to the same precision after approximately $34$ iterations. The initial guess used in both instances was $\Delta{{\mathcal{A}}}_{init}=0$. The initial guesses of ${\mathcal{F}}^{*}$ and $h$ were computed by an approximate GCD routine. For the initial guess of $\omega$ we chose a root or local extrema of $\det({{\mathcal{A}}})$ that minimized the second-smallest singular value of ${{\mathcal{A}}}(\omega)$, one of which is $\omega_{init}\approx-.12793-1.0223i$.

Example 6.2** (Lowest McCoy Rank Approximation).**

Let ${{\mathcal{A}}}$ be as in the previous example and consider the [math]-McCoy rank approximation problem with the same prescribed perturbation structure.

In this case we compute a local minimizer ${{\mathcal{A}}}+\Delta{{\mathcal{A}}}_{opt}$ given by

[TABLE]

with ${\|\Delta{{\mathcal{A}}}_{opt}\|}\approx.824645447014665$ after $34$ iterations to $14$ decimal points of accuracy. We compute $\omega_{opt}\approx-1.18536618732372i$ which corresponds to the single invariant factor $s_{1}\approx t^{2}+1.4051$. The SNF of ${{\mathcal{A}}}+\Delta{{\mathcal{A}}}_{opt}$ is of the form $(s_{1},s_{1},s_{1},s_{1})$.

7 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 are shown to be theoretically robust with a suitable initial guess. In general, reasonable quasi-Newton methods will have rapid local convergence under normalization assumptions for all the problems considered.

There are a number of problems that remain open for future work. In particular in the case of nearby nontrivial Smith forms there is the question of obtaining such forms via polynomial row and column operations, that is, finding the unimodular matrix multipliers that will produce our nearest Smith form. Preliminary work on this topic, including the formulation as an optimization problem and the proving of the existence of Lagrange multipliers for the optimization can be found in the thesis of Haraldson (2019). In some cases it may be practical to prescribe the degree structure, also called the structural supports, of the eigenvalues or the invariant factors of a nearby matrix polynomial. In this case, rather than look for a closest non-trivial SNF one would be interested in a closest SNF having a particular degree structure. As before this can be formulated as an optimization problem with early results available in (Haraldson, 2019).