Computing Lower Rank Approximations of Matrix Polynomials

Given an input matrix polynomial whose coefficients are floating point numbers, we consider the problem of finding the nearest matrix polynomial which has rank at most a specified value. This generalizes the problem of finding a nearest matrix polynomial that is algebraically singular with a prescribed lower bound on the dimension given in a previous paper by the authors. In this paper we prove that such lower rank matrices at minimal distance always exist, satisfy regularity conditions, and are all isolated and surrounded by a basin of attraction of non-minimal solutions. In addition, we present an iterative algorithm which, on given input sufficiently close to a rank-at-most matrix, produces that matrix. The algorithm is efficient and is proven to converge quadratically given a sufficiently good starting point. An implementation demonstrates the effectiveness and numerical robustness of our algorithm in practice.


Introduction
Matrix polynomials appear in many areas of computational algebra, control systems theory, differential equations, and mechanics. The algebra of matrix polynomials is typically described assuming that the individual polynomial coefficients come from an exact arithmetic domain. However, in the case of applications these coefficients typically have numeric coefficients, usually real or complex numbers. As such, arithmetic can have numerical errors and algorithms are prone to numerical instability.
Numerical errors have an impact, for example, in determining the rank of a matrix polynomial with floating point coefficients. In an exact setting determining the rank or determinant of a matrix polynomial is straightforward, and efficient procedures are available, for example from Storjohann and Villard (2005). However, in a numeric environment, a matrix polynomial may appear to have full or high rank while at the same time being close to one having lower rank. Here "close" is defined naturally under the Frobenius norm on the underlying coefficient matrices of the matrix polynomial. Rather than computing the rank of the given matrix polynomial exactly, one can ask how far away it is from one that is rank-deficient, and then to find one at that distance. In the case of matrices with constant entries this is a problem solved via the Singular Value Decomposition (SVD). However, in the case of matrix polynomials no equivalent rank revealing factorization has thus far been available.
In this paper we consider the problem of computing the nearest matrix polynomial to an input matrix polynomial in R[t] m×n having rank at most a specified value r. More precisely, given an integer r and an A ∈ R[t] m×n of full rank, we want to compute ∆A ∈ R[t] m×n with deg(∆A ij ) ≤ deg A ij (or similar degree constraints to be specified later), such that A + ∆A has rank at most n − r and where ∆A is minimized. In the case where n − r is one less than the row or column size then this is the problem of finding the nearest matrix polynomial which is singular.
A reasonable metric for measuring closeness on the space of matrix polynomials over the reals is the Frobenius norm. For a matrix polynomial A ∈ R[t] m×n , with (i, j) entry A ij ∈ R[t], the Frobenius norm is given by where, for a polynomial a ∈ R[t], the coefficient 2-norm is defined by a = 0≤i≤deg a a i t i , a 2 = a 2 2 = 0≤i≤deg a a 2 i . (1. 2) The main results in this paper center on the characterization of the geometry of minimal solutions. We show that minimal solutions exist, that is, for a given r there exists a ∆A ∈ R[t] m×n of minimal norm such that A +∆A has rank at most n − r and meets the required degree constraints on perturbed coefficients. In addition, we show that minimal solutions are isolated and are surrounded by a non-trivial open neighbourhood of non-minimal solutions. Also regularity and second-order sufficiency conditions are generically satisfied and a restricted version of the problem always satisfies these conditions. Finally we show that we can also generalize our results to the lower rank approximation instance of matrix polynomials generated by an affine structure 1 , and so generalize to low-rank approximations of structured matrices by taking the degree to be zero.
We demonstrate efficient algorithms for computing our minimal lower rank approximants. That is, for an input matrix polynomial A ∈ R[t] m×n (with prescribed affine structure) sufficiently close to a singular matrix polynomial, we give an iterative scheme which converges to a rank at most matrix polynomial at minimal distance, at a provably quadratic rate of convergence. We further generalize the iterative scheme so that it converges to a matrix polynomial with a kernel of dimension at least r, at a minimal distance and a provable quadratic rate of convergence. Finally, we also discuss a Maple implementation which demonstrates the convergence and numerical robustness of our iterative scheme.

Previous research
Much of the work in this area has often been done under the heading of matrix pencils. See Gohberg et al. (2009) for an excellent overview. Nonsingular (full rank) square matrix polynomials are sometimes referred to as regular matrix polynomials.
In the case of finding the nearest singular matrix pencil this problem was solved by the present authors in Giesbrecht et al. (2017). Previous to that this problem was posed for linear matrix pencils in Byers and Nichols (1993) and followed up in Byers et al. (1998). The nearest singular matrix polynomial relates to the stability of polynomial eigenvalue problems, linear time invariant systems and differential-algebraic equations studied subsequently in (Kressner and Voigt, 2015;Guglielmi et al., 2017). For non-linear matrix polynomials/pencils, previous works rely on embedding a non-linear (degree greater than 1) matrix polynomial into a linear matrix polynomial of much higher order. Theorem 1.1 and Section 7.2 of Gohberg et al. (2009) shows that any regular A ∈ R[t] n×n of degree d, is equivalent to a linear matrix polynomial B = B 0 + tB 1 , for B 0 , B 1 ∈ R nd×nd . However, this equivalence is (obviously) not an isomorphism, nor is it distance preserving 2 . Hence a nearby singular matrix polynomial to B ∈ R[t] nd×nd (even when constrained to a degree one perturbation) almost certainly does not correspond to a 1 A matrix A ∈ F m×n has an affine structure over a ring F if it can be written as . . , B L } ⊆ F m×n and c i ∈ F. If B 0 is the zero matrix, then the structure is said to be linear. Examples of linear structures include symmetric and hermitian matrices while matrices with an affine structure include entries that are fixed non-zero coefficients, such as monic matrix polynomials.
2 The equivalence mapping is not surjective.
nearby singular matrix polynomial to A ∈ R[t] n×n . Moreover, even if one was to perturb to a rank-reduced matrix within the image of the linearization, the inverse image would not necessarily have reduced rank. In Lawrence and Corless (2015) a more sophisticated linearization with an eye towards ameliorating this is explored.
In the context of computer algebra the notion of symbolic-numeric algorithms for polynomials has been an active area of research for a number of years, and the general framework of finding nearby instances with a desired algebraic property is being thoroughly explored. Closest to our work here is work on approximate Greatest Common Divisors (GCD) Corless et al. (1995); Beckermann and Labahn (1998b,a), multivariate polynomial factorizations Kaltofen et al. (2008), and especially the optimization-based approaches employing the Structured Total Least Norm algorithm Li et al. (2005); Kaltofen et al. (2005Kaltofen et al. ( , 2006; Zhi (2007) and Riemannian SVD Botting et al. (2005). More recently, we have explored computing the approximate GCRD of (non-commutative) differential polynomials (Giesbrecht and Haraldson, 2014;Giesbrecht et al., 2016) and resolve similar issues.
The computer algebra community has made impressive progress on fast, exact algorithms for matrix polynomials, including nearly optimal algorithms for computing ranks, factorizations and various normal forms; see Kaltofen and Storjohann (2015) and references therein for a recent overview. Part of our goal in this current paper is establish a basis for extending the reach of these symbolic techniques to matrices of polynomials with floating point coefficients.
In a more general setting our problem can be formulated as a Structured Low Rank Approximation (SLRA) problem. A popular method to solve SLRA problems is the Structured Total Least Norm (STLN) approach (Rosen et al., 1996(Rosen et al., , 1998. These are iterative methods and in general their convergence to stationary points is linear (first order), rather than quadratic, unless additional assumptions are made. In the event STLN converges to a solution, there may be other solutions arbitrarily nearby, as second order sufficient conditions may not hold. The SLRA problem is a non-linear least squares problem and accordingly other techniques such as the Restricted and Riemannian SVD (De Moor, 1993, 1994, 1995 provide general tools for solving such problems. Other heuristic tools applicable to our problem include variable projection Pereyra, 1973, 2003) and Newton's method (Abatzoglou et al., 1991). We would expect these methods to perform very poorly in our case, as one can expect problems with large residuals to perform poorly and the rational function arising from variable projection can be too costly to deal with for modestly sized problems. The problem may also be considered as optimization on a manifold (Absil et al., 2009), however we do not explicitly consider this approach. For a detailed survey of affinely structured low-rank approximation, see (Markovsky, 2008(Markovsky, , 2011. Other methods for structured low-rank approximation involve the family of lift and project algorithms, with the best known being Cadzow's algorithm (Cadzow, 1988). More recently Schost and Spaenlehauer (2016) gives a sequence of alternating projections that provably converge quadratically to a fixed point. However, lift and project algorithms do not generally satisfy necessary first order (see (Bertsekas, 1999)) optimality conditions, and while they may converge (quickly) to a fixed point, there is no guarantee that the fixed point is an optimal solution, though it is usually quite good. In any case, for specific problems such as ours, understanding the geometry of the minimal solutions (and hence the well-posedness of the problem) is key to effective algorithms for their computation.
SLRA problems are in general NP-hard to solve, see for example (Poljak and Rohn, 1993;Braatz et al., 1994). They are also hard to approximate under affinely structured matrices over Q. In general the hardness stems from determining if a bilinear system of equations admits a non-trivial solution. In the instance of classical matrix polynomials it is trivial to construct feasible points since the underlying scalar matrix problem is linearly structured.
All of our contributions apply to matrix polynomials with an affine structure provided that feasible points exist, that is, singular matrix polynomials with a prescribed structure exist, which is NP-hard in general. In particular, in the degree zero case our algorithms and techniques apply to affine SLRA problems. Thus, computing the nearest (affinely structured) matrix polynomial is equivalent to SLRA problems with an affine structure.
While the contributions in this paper focus on local properties of SLRA, the local properties also imply global results. The Sum of Squares (SOS) hierarchy is a global framework for studying polynomial optimization problems subject to polynomial constraints Lasserre (2001). The SOS optimization tools have found experimental success in computing structured distances to singularity and extracting minimizers when the solutions are locally unique, see for example Henrion and Lasserre (2006). In general the SOS hierarchy converges for an infinite order of relaxations, but for several problems the relaxations converge after a finite order. The finite convergence is in polynomial time with respect to the input and the number of relaxations. In particular, this finite convergence was observed for affine SLRA problems in Henrion and Lasserre (2006) but little theory was provided to indicate the reason why. The later work of Nie (2014) shows that, under regularity and second-order sufficiency conditions, finite convergence always occurs and that it is possible to extract a minimal solution. In our contributions we prove that second-order sufficiency and regularity conditions hold generically (and if they do not, then they will hold on a restricted subset of the problem). The corollary to this is that the SOS hierarchy will have finite convergence for affine SLRA problems if a solution exists (such as computing the distance of the nearest rank-deficient matrix polynomial) and if the embedding is minimal then a minimizer may be extracted as well. Another useful feature of the SOS hierarchy is even if convergence cannot be certified, a structured lower-bound is obtained.

Outline
In Sections 2 and 3 we describe tools needed for our constructions and then explore the geometry of our problem. We show that the problem is locally well-posed. One cannot expect the nearest rank at most matrix polynomial to be unique. However under weak normalization assumptions, we show that solutions are locally unique in a closed-ball around them. To complement the separation of solutions, we also show that for an equivalent problem, solutions corresponding to a different closed ball are separated by at least a constant amount independent of the dimension of the space.
In Section 4 we give an equality constrained variant of Newtons' method for computing via post-refinement the nearest rank at most matrix polynomial. The main idea is to compute an initial guess with a suitable first order or lift-and project method. We are able to prove that, with a suitable initial guess and regularity assumptions, our algorithm generally has local quadratic convergence except for degenerate cases. This is done by deriving closed-form expressions for the Jacobian of the constraints and the Hessian of the Lagrangian. When we refer to the speed of convergence, we refer to quotient rates as is typical in the nomenclature.
In Section 5 we describe our prototype implementation, including heuristics for starting points and other improvements. We discuss the numerical performance of the algorithm and give examples demonstrating convergence. results for a low-rank approximation of matrix polynomials. The paper ends with a conclusion and topics for future research.

Preliminaries and Geometry
In this section we will introduce some basic definitions and explore the numerical geometry of our lower rank problem. Canonically we will let be a matrix polynomial, with coefficients A 0 , . . . , A d ∈ R n×n . In the case of rectangular matrix polynomials we are able to pad the matrix with zeros, thus embedding the problem into one with square matrix polynomials. Thus we will let is a Toeplitz matrix. Such matrices are conveniently used to describe polynomial multiplication in the sense that if c = a · b with a of degree d and c ∈ R[t] of degree at most µ − 1, then vec(c) = φ(a) · vec(b) where vec(p) is the vector of coefficients of a polynomial.
. This property is central to our work in the coming sections.
For ease of notation we will take N = n(µ + d) = n 2 d + n(d + 1), M = nµ = n 2 d + n and R ≥ 1 when dealing with R-embeddings in subsequent sections. We note that A is a block-Toeplitz matrix, and as such one method of understanding the lower rank problem is to find close by structured rank deficient block-Toeplitz matrices, a typical structured low rank approximation problem. Some authors refer to such embeddings as a (permuted) Sylvester matrix associated with A. We avoid this terminology as it is ambiguous when considering Sylvester matrices occurring in (approximate) GCD computations. Unlike the standard linearizations in (Gohberg et al., 2009, Section 7.2) used to turn arbitrary degree matrix pencils into linear pencils, this Rembedding is kernel preserving for matrix polynomials of arbitrary degree.
Note that this is minimizing a convex objective function subject to nonconvex constraints. However, the equality constraints are linear in each argument. It is still not clear that Problem 2.2 is well-posed in the current form. We will prove that solutions exist, that is, there is an attainable global minimum value and not an infimum.
Proof. Suppose that A has rank s < n. By permuting rows and columns we may assume without loss of generality that the leading s × s submatrix of A is non-singular. There is a unique vector of the form Multiplying through by γ, we find that b = γc satisfies the requirements of the lemma.
Lemma 2.4. A is singular if and only if A does not have full column rank.
A has a non-trivial kernel and, b ∈ ker A by construction. Conversely, suppose that A has full rank. Then for all b ∈ R[t] n×1 we have Ab = 0 which implies that A b = 0 or ker A is trivial.
We recall the Singular Value Decomposition as the primary tool for finding the distance to the nearest unstructured rank deficient matrix over R or C.
is a diagonal matrix consisting of the singular values of C in descending order of magnitude. See (Golub and Van Loan, 2012).
The following fact is a standard motivation for the SVD.
Fact 2.6 (Eckart and Young (1936)). Suppose C = QΣP T ∈ R N ×M as above has full column rank, with N ≥ M . Then ∆C = Q diag(0, . . . , 0, −σ M )P T is such that C + ∆C has column rank at most M − 1, ∆C F = σ M , and ∆C is a perturbation of minimal Frobenius norm which reduces the column rank of C.
Lemma 2.7. Given a non-singular A ∈ R[t] n×n , and ∆A ∈ R[t] n×n such that B = A + ∆A is singular, it is the case that ∆A ≥ σ nµ ( A).
Proof. By Lemma 2.4 above, B is not of full column rank. Thus, by Fact 2.6 ∆A F ≥ σ nµ (A).
Corollary 2.8. The set of rank r > 1 matrices over R[t] n×n of degree at most d is open, or equivalently, the set of all matrices of rank at-most n − r over R[t] n×n of degree at most d is closed.
S is the intersection of a closed and bounded set and a closed set, hence S is closed and bounded. S is isomorphic to some closed and bounded subset of Euclidean space, hence by the Heine-Borel theorem, S is compact. To show the set is non-empty, we note that, by the degree assumption on ∆A, ∆A = −A is a feasible point independent of rank. Let C ∈ S then A − C 2 F = ∆A 2 F is a continuous function over a compact set. By Weierstrass' theorem it has an attainable global minimum.
It is important not to over-constrain the problem with a choice of ∆A, since otherwise the feasible set might be empty. Another reasonable choice of ∆A which we can handle, is that the perturbation has the same coefficient structure/support as A, that is, zero terms in polynomial entries are preserved.
We note that this result says nothing about uniqueness or separation of solutions or any local properties. All that has been shown is that if the perturbations are in the same space as the input, and one seeks a rank atmost approximation, then there is an attainable global minimum value, i.e. not an infimum. If one wants a minimal solution with the rank being exactly n − r, then there is no guarantee that there is an attainable global minimum to Problem 2.2.

Rank Factorization
A natural formulation of the problem that encompasses the rank implicitly is to perform a rank factorization and write A + ∆A = UV for U ∈ R[t] n×(n−r) and V ∈ R[t] (n−r)×m . Here UV is subject to some constraints that preserve the structure of ∆A (i.e., that we do not perturb any coefficients we are not allowed to, typically that deg ∆A ij ≤ deg A ij , but possibly also preserving the zero coefficients and not introducing a larger support). This is a non-linear least squares problem. However solutions are not unique. Indeed, if Z ∈ R[t] (n−r)×(n−r) is unimodular (i.e., det(Z) ∈ R * ), then UZ, Z −1 V is another rank n − r factorization, and we obtain an infinite family. While normalizing over matrix polynomial rank-factorizations is difficult, it is much easier to exploit the quasi-distance preserving property of · F and look at rank-factorizations of A, that do not necessarily correspond to U and V.

Embedded Rank Factorization
Definition 3.1. Let N = (µ+d)n, M = nµ and R > 0. A rank factorization of A + ∆A is given by writing A + ∆A = U V where U ∈ R N ×R and V ∈ R R×M are arbitrary (unstructured) matrices over R.
Our goal is to find U, V with appropriate dimensions which minimize and such that ∆ A has the correct Toeplitz-block structure (i.e., it is an Rembedding of a matrix polynomial). This is a problem with a non-convex objective function (that is convex in each argument) and non-convex constraints. We note that U , V have no direct connection with U, V ∈ R[t] n×n .
One may always write A + ∆A this way via the SVD for fixed A and ∆A, so in particular the optimal solution can be written as a rank factorization. The problem min A − U V 2 such that U V has the same structure as ∆A is generally ill-posed and needs to be constrained to do any meaningful analysis, as there are numerous degrees of freedom. At first glance, optimizing over rank factorizations appears to be a harder problem than the original. However it is helpful to perform analysis on this formulation. In particular, we are able to prove that optimal values of ∆A that satisfy first order conditions (which contains all useful perturbations) are separated by a constant amount, and that equivalence classes of solutions are isolated. Additionally, this formulation of the problem is convex in each argument (but not jointly convex) and is amenable to block coordinate descent methods.
We next need to demonstrate that the condition that the matrix ∆ A = A − U V is the R-embedding of some matrix polynomial ∆A ∈ R[t] n×n can be phrased as a single polynomial being zero. A is generated by a linear We note that there exist c i such that Γ(∆ A) = 0 if and only if ∆ A is an R-embedding of a matrix polynomial.
Problem 3.2. With A, U, V as above, the constrained R-embedded rank factorization problem consists of computing min A − U V 2 F subject to the constraints that U T U − I = 0 and Γ(U, V ) = 0. If R = M − 1, then this encodes all rank deficient matrix polynomials.
It is still not clear that Problem 3.2 is well-posed, as there are many degrees of freedom in V , and this matrix can have arbitrary rank. The enforcement of U as an orthogonal matrix (U T U − I = 0) is allowed for without loss of generality. Informally then we are looking at all rank factorizations where where U is orthogonal and Γ(U, V ) = 0, that is, the product satisfies the block-Toeplitz structure on ∆A.
We employ the machinery of non-linear optimization to describe the geometry of the minimal solutions, and hence the nearest appropriately structured matrices. See (Bertsekas, 1999) for an excellent overview.  (Bertsekas (1999, Section 3.1.1)). For a sufficiently large ρ > 0, one has that 3 Problem 3.2 is equivalent to computing a solution to the unconstrained optimization problem All the interesting solutions to the minimization of Φ(U, V ) occur at stationary points. The first-order necessary condition (on V ) of gradients vanishing gives us If we assume that the constraints are active, that is U is orthogonal and that Γ(U, V ) = 0, then we have U T A − V = 0. Of course, there is the other first order necessary condition requiring that However, we do not need to employ this explicitly in the following.
Theorem 3.4 (Strong Separation of Objective). Suppose ∆A and ∆A are distinct (local) optimal solutions to Problem 3.2 that satisfy first order necessary conditions. Then ∆A − ∆A 2 ≥ σ min ( A), where σ min (·) is the smallest non-trivial singular value.
Proof. From the previously discussed necessary first order condition we have that there exists U ∈ R N ×R , V ∈ R R×M and U ∈ R N ×R and V ∈ R R ×M such that Note that R and R need not be the same. From this we obtain the sequence of lower bounds The symmetric matrix W = U T U U T U is a product of matrices whose non-zero eigenvalues have magnitude 1. Symmetric matrices have real eigenvalues, and the non-zero eigenvalues of W will be ±1, since U and U are orthogonal. Thus W 2 = 1.
W must have at least one negative eigenvalue or zero eigenvalue by the orthogonality assumption, since W = I. Since W is symmetric, we can diagonalize W as a matrix with ±1 and 0 entries on the diagonal. It follows that I − W 2 ≥ 1 and the theorem follows.
While the separation bound exploited properties of the rank factorization, these bounds hold for all formulations of the problem.
Corollary 3.5. All locally optimal solutions satisfying first order necessary conditions are isolated modulo equivalence classes.
Proof. Suppose the contrary, that is that (U, V ) is a solution corresponding to ∆A and (U , V ) is a solution corresponding to ∆A . The objective function and constraints are locally Lipschitz continuous, so let s > 0 be a Lipschitz constant with respect to · F in some open neighborhood.
If we take 0 < ε < σ min ( A) s then we have which is a contradiction to Theorem 3.4.
Implicitly the matrix V parametrizes the kernel of A. If we normalize the kernel of A to contain R-embeddings of primitive kernel vectors then the matrix V can be made locally unique, although we do not employ this in the rank-factorization formulation directly.
Corollary 3.6. Under a suitable choice of · we have that minimal solutions are separated. In particular, separation holds for · 1 .
The proof follows immediately from equivalence of matrix norms, as norms are equivalent in a finite dimensional space.
While there are too many degrees of freedom to easily obtain a (locally) quadratically convergent minimization over the rank factorization, the rank factorization does yield non-trivial insights into the geometry of the solution space. In particular, the isolation of solutions indicates first order (gradient) methods will perform well on the problem. In the next section we will introduce a locally quadratically convergent algorithm for an equivalent form of Problem 2.2 that reduces each equivalence class of solutions to a single solution.

An Iterative Algorithm for Lower Rank Approximation
In this section we propose an iterative algorithm to solve Problem 2.2 based on Newton's method for constrained optimization. Sufficient conditions for quadratic convergence are that the second-order sufficiency holds (Wright, 2005) and local Lipschitz continuity of the objective and constraints. We ensure these conditions hold for non-degenerate problems by working on a restricted space of minimal R-embeddings that remove degrees of freedom.

Minimal System of Equations
In order to compute a nearby rank n − r approximation we want to solve the non-convex optimization problem min ∆A 2 F subject to (A + ∆A)B = 0 rank(B) = r. (4.1) In the instance of (structured) scalar matrices the rank constraint can be enforced by ensuring that B has orthogonal columns 4 or is in a column reduced echelon form. In the instance of matrix polynomials this is not sufficient, since polynomial multiples of the same vector will have linearly independent combined coefficient vectors. In order to apply these normalizations on the coefficient vectors of B we require that the columns be represented with a minimal number of equations with respect to B.
, a dimension 1 subspace). We say that b is minimally degree R-embedded in A if (1) b is minimally R-embedded in A and (2) b corresponds to a primitive kernel vector b, that is gcd(b 1 , . . . , b n ) = 1.
We note that this definition ensures minimally R-embedded vectors are unique (up to scaling a factor), or that ( A j + ∆A j ) B[ * , j] = 0 has a (locally) unique solution for fixed ∆A. In the minimal embedding, we will assume, without loss of generality, that redundant or equations known in advance, such as 0 = 0, ∆ A ij = 0 or B ij = 0 corresponding to known entries are removed for some indices of i and j. Some of these trivial equations occur because of the CREF assumption, while others occur from over-estimating degrees of entries.
This allows us to reformulate (A + ∆A)B = 0 as a (bi-linear) system of equations where the j th column of B is minimally degree embedded in the system ( A j + ∆ A j ). We also note that assuming B is in a column-reduced echelon form essentially requires us to guess the pivots in advance of the optimal solution, which is only possible with a good initial guess. The benefit of this approach is that if the pivots are not guessed correctly, we are still able to compute a n − r approximation of A.
In order to exclude trivial solutions, we can assume that the pivot elements of B have a norm bounded away from zero. Let N( b i ) be a normalization vector such that N( b i ) T b i = 1 which implies that the CREF pivots are bounded away from zero. For example, take the pivot to have unit norm. Note that other normalization vectors are possible, such as N( b i ) = b i (which corresponds to each column having a unit norm) if the initial guess is adequately close, or we could take the pivot element to be a monic polynomial. Of course there are several other permissible normalizations.
Define the matrix A i to have the column b i = B[1..n, i] minimally degree embedded. We can express (4.2) in a vector-matrix form as follows.
has a (locally) unique solution for fixed ∆A.

Lagrange Multipliers and Optimality Conditions
In order to solve (4.1) we will use the method of Lagrange multipliers (Bertsekas, 1999).
Let M (∆A, B) be the vector of residuals corresponding to (4.3), then the Lagrangian is defined as where λ = (λ 1 , . . . , λ # residuals ) T is a vector of Lagrange multipliers. We will find it convenient to define x = vec(∆A) vec(B) to be the combined vector of unknowns corresponding to A and B. Let ∇ 2 xx L denote the Hessian matrix of L with respect to x and J be the Jacobian of the residuals of the constraints, i.e. J = ∇ x M (∆A, B). Necessary optimality conditions at a point (x * , λ * ) (Bertsekas, 1999)  These conditions are known as the second-order sufficiency conditions Bertsekas (1999). We note that (4.6) implies that minimal solutions are locally unique, and will fail to hold if minimal solutions are not locally unique. The idea is to show that (4.6) holds in the minimal embedding, which allows us to construct an algorithm with rapid local convergence. We will adopt that notation that ψ( b i ) corresponds to ψ( b i )vec( A i + ∆ A i ) = 0. Here we use the bi-linearity of (4.3) to write the same system using a matrix with entries from B instead of vec(A + ∆A).

The Jacobian
The closed-form expression for the Jacobian of the residuals (up to permutation) in (4.3) is given by (4.7) Unlike the case of a single kernel vector in (Giesbrecht et al., 2017), J may be rank deficient since some equations corresponding to low (high) index entries may be redundant at the solution. The Lagrange multipliers will not be unique in this particular scenario and the rate of convergence may degrade if Newton's method is used. In the instance of r = 1 then we present the following result (Giesbrecht et al., 2017).
Theorem 4.4. Suppose that r = 1 and b 1 is minimally degree R-embedded in A 1 , then J has full rank when (4.5) holds.
Proof. We show that J has full row rank by contradiction. If this matrix was rank deficient, then one row is a linear combination of the others. This means that one of the equations in the constraints is trivial or the solution is not regular (see (Bertsekas, 1999, Section 3.1)). As we are only concerned about regular solutions, this contradicts the minimal degree R-embedding.
The corollary to this is that in the minimal embedding regularity conditions hold and it is straight forward to obtain rapid local convergence.

The Hessian
The Hessian matrix, ∇ 2 L is straight forward to compute as The following theorem shows that second-order sufficiency holds for the instance of r = 1. The case of r > 1 follows immediately by induction. This is in contrast to Theorem 4.4, which does not always hold for r > 1.
Theorem 4.5 (Second Order Sufficiency Holds). Suppose that A + ∆A has a minimally degree R-embedded kernel vector b, i.e. r = 1 in (4.4), then at a minimal solution, the second order sufficiency condition (4.6) holds in the minimal embedding of b.
Proof. If ∆A = 0 at the local minimizer (x * , λ * ) then We have that for y ∈ span(K) such that Jy = 0 implies that Ay = 0 and N( b) T y = 0. It follows that ker A = span( b), thus we have y = b or y = 0 via the minimal degree R-embedding, thus y = 0 as b / ∈ span(K). Hence, second-order sufficiency holds, as ker J ∩ K = 0.
If ∆A = 0 then we have that The matrix E is linear in λ, however the precise tensor decomposition is irrelevant to the proof. If E has full rank, then ∇ 2 xx L has full rank and we are done, so suppose that E is rank deficient. If E is rank deficient, then one can eliminate a row of E and column of E T without affecting H via symmetric row and column updates. We observe that ker(H + E) ⊆ ker H and the result follows.
Proof. The proof is almost the same as Theorem 4.5 and follows by induction on r since each block is decoupled.
We now have all of the ingredients for an iterative method with rapid local convergence.

Iterative Post-Refinement
Newton's method for equality constrained minimization problems can be interpreted as solving the non-linear system of equations ∇L = 0. Newton's method is based on the iterative update scheme (4.8) If r = 1 then ∇ 2 L has full rank and the iteration is well defined by matrix inversion. If r > 1 then we consider the quasi-Newton method defined as for a suitably chosen parameter µ k . Taking µ k = ∇L(x k , λ k ) 1 one has provably quadratic convergence (Wright, 2005, Theorem 4.2) with x k and λ k chosen sufficiently close to the optimal solution.
We now have a method to compute a nearby rank deficient matrix polynomial with a rate of convergence that is quadratic, provided that the initial values of x are chosen to be sufficiently close to the optimal solution.

Implementation, Description and Comparison
In this section we discuss implementation details and demonstrate our implementation for computing the nearest rank deficient matrix polynomial. All algorithms are implemented in Maple 2016. All experiments are done using quad precision floating point arithmetic, with about 35 decimal digits of accuracy. We compare some degree one examples to the recent results of (Guglielmi et al., 2017).
To compute an approximate kernel vector, first we use the SVD to compute an approximate kernel of an R-embedded (nearly) rank deficient matrix polynomial. Next we use structured orthogonal elimination RQ (LQ) decomposition to produce a minimally (degree) R-embedded vector from the kernel. In the case of several kernel vectors we use a modified Gaussian elimination on an embedding of an approximate kernel obtained by the SVD and approximate GCD to find nearby approximate kernel vectors that are primitive.

Description of Algorithm
We now formally describe an algorithm for computing the nearest matrix polynomial of a prescribed rank. The algorithm has no global convergence guarantees, however a globally convergent (although not necessarily optimal) algorithm can be developed in a straight forward manner via augmenting our second order algorithm with a first order one, and removing content from kernel vectors if necessary.
The size of ∇ 2 L is O(r 2 n 4 d 2 ) and accordingly each iteration has a cost of O(r 6 n 12 d 6 ) flops using standard matrix multiplication, where r is the dimension of the kernel.

Algorithm 1 : Iterative Kernel Post-Refinement
Require: • Full rank matrix polynomial A ∈ R[t] n×n • (Approximately) Rank deficient matrix polynomial C ∈ R[t] n]×n • Approximate kernel vectors c 1 , . . . , c r ∈ R[t] n×1 of the desired degree/displacement structure • Displacement structure matrix ∆A to optimize over Ensure: • Singular matrix A + ∆A with B ⊂ ker(A + ∆A) or an indication of failure.

4: Compute
x + ∆x λ + ∆λ by solving (4.9) until ∆x ∆λ 2 is sufficiently small or divergence is detected. 5: Return the locally optimal ∆A and B or an indication of failure.

Nearest Rank Deficient Linearly and Affinely Structured Matrix
In this section we consider Examples 2.10, 2.11 and 2.12 from Guglielmi et al. (2017), where we compare our results to real perturbations. Note that complex perturbations are a straight-forward generalization of the theory presented here, and can be re-formulated as a problem over R.
The technique of Guglielmi et al. (2017) poses computing a nearby rankdeficient linear matrix pencil by verifying that sufficiently many images of the matrix polynomial are singular, so that det(A + ∆A) ≡ 0. The problem is then posed as a solution to a system of Ordinary Differential Equations (ODE), assuming that certain genericity conditions on the eigenvalues of the solution hold 5 . They consider the instances of computing A 0 and A 1 with a common kernel vector, and the instance where A 0 and A 1 do not have a common kernel. Additionally, perturbations affecting only one of A 0 and A 1 are considered. We note that the solutions to the ODEs do not necessarily satisfy necessary optimality conditions (4.5), and accordingly will generally not be local minimizers.

Nearest Affinely Structured Examples I
Consider first the matrix polynomial coming from Examples 2.10 and 2.12 of Guglielmi et al. (2017) Example 5.1. If we assume that A 1 is constant, then this is finding the (locally) nearest matrix polynomial with an affine structure since A 1 has nonzero fixed constants. First let's assume that zero entries are preserved, this is a linear structure on A 0 .
To compute an initial guess for b we use the SVD on A and extract a guess from the smallest singular vector. This gives us For an initial guess on A we take A init = A. Note that we do not need an initial guess that is singular, it just needs to be "sufficiently close" to a singular matrix polynomial.
If we do not allow perturbations to zero-coefficients, that is, A 0 [1, 1] and A 0 [2, 3] may not be perturbed, then after five iterations of plain Newton's method (see (Giesbrecht et al., 2017)) we compute  Guglielmi et al. (2017) report an upper-bound on the distance to singularity allowing complex perturbations, that is ∆A ∈ C[t] n×n of ∆ C A F ≈ 0.1357 in Example 2.10. In Example 2.12, Guglielmi et al. (2017) report an upper-bound on the distance to singularity allowing real perturbations, ∆ R A F ≈ 0.1366. Although we only consider real perturbations, both bounds are improved. We conjecture that the complex bound can be improved further.
If we allow perturbations to A 0 and A 1 , then this is some form of finding the nearest rank deficient matrix polynomial. The question is whether to allow degree or support preserving perturbations. Again, we will use the same initial guesses as the previous example.
Example 5.3. In the instance of degree preserving perturbations we compute after five iterations of Newton's method Example 5.4. In the instance of support preserving we compute after five iterations of Newton's method, Guglielmi et al. (2017) report an upper-bound on the distance to singularity of ∆ deg A F ≈ 0.1193 in Example 2.12. This bound is larger than the one computed in Example 5.3.

Nearest Affinely Structured Examples II
Example 5.5. Next we consider the the matrix polynomial A in Example 2.11 of (Guglielmi et al., 2017)  To compute an initial guess for we take A init = A and take b init is computed from the smallest singular vector of A. We note that this initial guess does not attempt to find a nearby singular matrix polynomial for the initial guess, all that is needed is ∇L(x init , λ init ) is reasonably small to obtain convergence.
Using a globalized variant of Newton's method based on Levenberg-Marquardt we compute The previously noted small quadratic terms were at roughly machine precision (the computation is done with 35 digits of precision) and truncated. Guglielmi et al. (2017) obtain a result on this past example that produces an upper bound on the distance to singularity of 0.9438619. Their computation is accurate to seven decimal points, and accordingly our post-refinement has an improvement of about 0.000297. This is not surprising, since we solve the necessary conditions (4.5) directly with a reasonable initial guess.

Lower Rank Approximation of a 4 × 4 Matrix
In this following example we consider computing a lower-rank approximation to a given matrix polynomial. Consider the 4 × 4 matrix polynomial A, defined as Example 5.6. We will consider a displacement structure on the kernel as well in this example, where higher-order zero terms are not perturbed from the initial guess. For the entries of ∆A we preserve higher-order zero terms, and allow low-order terms to be perturbed. This is a linearly structured problem, on both the main variable ∆A and the auxiliary kernel variable B.
To ensure the rank constraint holds, we will additionally assume that the kernel, B is in a CREF (while B is obviously not) and the columns have unit norm. This normalization is (locally) equivalent to the ones discussed in Section 4.2. Having B in a CREF ensures that the two kernel vectors are locally linearly independent during the iteration. Of course perturbing both pivots to zero is possible (although this is sub-optimal). In such a scenario linear independence can no longer be guaranteed, and the iteration would need to be re-ininitialized.
For the initial guess we use A init = A and take B init as A natural question is what happens if we change the displacement structure on the kernel? To investigate this behavior, we consider an equivalent representation of the previously used kernel, except that B is in a CREF directly.
Example 5.7. If we change the kernel B init to be used in the initialization of the previous example, then we compute a perturbation with ∆A ≈ 0.0008408.
In either case, we obtain comparable answers that are a reasonable lowerrank approximation, and can likely be improved by relaxing restrictions on the displacement structure on B or B. It is important to note that relaxing the degree bounds to be (n − r)d in general on all non-zero entries (where entries are zero if they are in the same row as a CREF pivot) will likely lead to a better approximation, however one may lose quadratic convergence if doing so, since iterates may no longer have primitive kernel vectors, and (4.6) will no longer hold. As discussed in Section 4, it is generally difficult to determine the CREF pivots of the kernel unless the initial guess is very accurate.
The structure of the kernel is an important consideration when deciding upon an initial guess. It is preferable to restrict fewer coefficients, however the iteration requires a better initialization due to the increased number of possible descent directions. In such scenarios for maximum flexibility, a globalized variant of Newton's method is required. Like-wise, the structure for ∆A is also an important choice. Restricting which terms can be changed has a large influence on the (approximate) distance to singularity (of prescribed kernel dimension).
Another way to approach the lower-rank approximation problem is to use alternating projections or alternating directions of descent (since the objective is bi-linear with bi-linear constraints, it is convex in each argument) on the rank factorization in Section 3. Since solutions in one coordinate, ∆A are isolated, one can expect linear convergence with a reasonable algorithm. The lack-of normalization required overcomes the difficulty of choosing a suitable kernel displacement structure, however convergence would be linear at best and determining the dimensions of U and V is another problem to be discussed. It is also worth noting that Algorithm 5.1 requires more computational resources per iteration as r increases, however a rank factorization requires fewer computational resources per iteration as r increases.

Conclusions and Future Work
We have shown that finding lower-rank approximations of matrix polynomials can be established as a numerically well-posed problem and is amenable to first and second order optimization methods. The existence and isolation of solutions is established along with an algorithm exploiting affine structures to obtain locally quadratic convergence under mild normalization assumptions.
Along with considering the lower-rank approximation of matrix polynomials, we present a generalization of the theory to matrix polynomials with an arbitrary affine structure. We provide examples of how the structure of permissible perturbations and prescribed kernel structure impacts the distance to solutions.
We also regard this current paper as a first step towards a formally robust approach to non-linear matrix polynomials, in the spirit of recent work with symbolic-numeric algorithms for polynomials. Problems such as approximate matrix polynomial division, GCRD and factorization all have applications which can benefit from these modern tools.