# Generalized Standard Triples for Algebraic Linearizations of Matrix
Polynomials††thanks: Received by the editors on Month/Day/Year. Accepted for
publication on Month/Day/Year. Handling Editor: Name of Handling Editor.
Corresponding Author: Name of Corresponding Author

Eunice Y. S. Chan Ontario Research Centre for Computer Algebra, School of
Mathematical and Statistical Sciences, Department of Applied Mathematics,
Western University<EMAIL_ADDRESS><EMAIL_ADDRESS>Robert M. Corless
Ontario Research Centre for Computer Algebra, School of Mathematical and
Statistical Sciences, Western University<EMAIL_ADDRESS>Leili Rafiee
Sevyeri33footnotemark: 3

###### Abstract

We define “generalized standard triples” $\bm{X}$, $z\bm{C}_{1}-\bm{C}_{0}$,
$\bm{Y}$ of regular matrix polynomials $\bm{P}(z)\in\mathbb{C}^{n\times n}$ in
order to use the representation $\bm{X}(z\bm{C}_{1}\leavevmode\nobreak\
-\leavevmode\nobreak\ \bm{C}_{0})^{-1}\bm{Y}\leavevmode\nobreak\
=\leavevmode\nobreak\ \bm{P}^{-1}(z)$ for $z\notin\Lambda(\bm{P}(z))$. This
representation can be used in constructing algebraic linearizations; for
example, for $\bm{H}(z)=z\bm{A}(z)\bm{B}(z)+\bm{C}\in\mathbb{C}^{n\times n}$
from linearizations for $\bm{A}(z)$ and $\bm{B}(z)$. This can be done even if
$\bm{A}(z)$ and $\bm{B}(z)$ are expressed in differing polynomial bases. Our
main theorem is that $\bm{X}$ can be expressed using the coefficients of the
expression $1=\sum_{k=0}^{\ell}e_{k}\phi_{k}(z)$ in terms of the relevant
polynomial basis. For convenience we tabulate generalized standard triples for
orthogonal polynomial bases, the monomial basis, and Newton interpolational
bases; for the Bernstein basis; for Lagrange interpolational bases; and for
Hermite interpolational bases. We account for the possibility of common
similarity transformations. We give explicit proofs for the less familiar
bases.

###### keywords:

Standard triple, regular matrix polynomial, polynomial bases, companion
matrix, colleague matrix, comrade matrix, algebraic linearization,
linearization of matrix polynomials.

65F15, 15A22, 65D05

##  1 Introduction

A matrix polynomial is normally defined using words such as “A matrix
polynomial $\bm{P}(\lambda)\in F^{m\times n}(\lambda)$ is a polynomial in the
variable $\lambda$ with coefficients that are $m$ by $n$ matrices with entries
from the field $F$.” Typically an expression in the monomial basis
$\lambda^{k}$ is given for $\bm{P}(\lambda)$, and only regular matrix
polynomials are considered, that is, with $m=n$ (we will use $n$ for the
dimension) and where $\det\bm{P}(\lambda)$ is not identically zero. Matrix
polynomials have many applications and their study is of both classic and
ongoing interest.

In this paper we wish to consider the case when any polynomial basis
$\phi_{k}(\lambda)$ is used, if the set $\\{\phi_{k}(\lambda)\\}$ for $0\leq
k\leq\ell$ forms a basis for polynomials of degree at most $\ell$. Thus, we
write our regular matrix polynomial as

(1.1) $\bm{P}(\lambda)=\sum_{k=0}^{\ell}\bm{P}_{k}\phi_{k}(\lambda)\>,$

where the matrices $\bm{P}_{k}\in F^{n\times n}$ are square, and the degree of
$\bm{P}(\lambda)$ is at most $\ell$. The upper bound $\ell$ on the degree is
also called the grade. The notion of “grade” is useful even for the monomial
basis, but it is especially useful if the basis is an interpolational basis or
the Bernstein basis, when the degree of the polynomial may not be clear from
the data.

Two matrix polynomials $\bm{P}_{1}(\lambda)$ and $\bm{P}_{2}(\lambda)$ are
called equivalent if there exist unimodular matrix polynomials (that is,
matrices with constant nonzero determinant) $\bm{E}(\lambda)$ and
$\bm{F}(\lambda)$ with
$\bm{P}_{1}(\lambda)=\bm{E}(\lambda)\bm{P}_{2}(\lambda)\bm{F}(\lambda)$. A
matrix pencil $\bm{L}(\lambda):=\lambda\bm{C}_{1}-\bm{C}_{0}$ is called a
linearization of the matrix polynomial $\bm{P}(\lambda)$ if both $\bm{C}_{1}$
and $\bm{C}_{0}$ are of dimension $N\geq n\ell$ and $\bm{L}(\lambda)$ is
equivalent to the block diagonal matrix
$\mathrm{diag}(\bm{P}(\lambda),\bm{I}_{N-n})$. Two linearizations
$\bm{L}_{m}(\lambda)$ and $\bm{L}_{\phi}(\lambda)$ are called strictly
equivalent if the corresponding matrices are equivalent:
$\bm{C}_{1,m}=\bm{E}\bm{C}_{1,\phi}\bm{F}$ and
$\bm{C}_{0,m}=\bm{E}\bm{C}_{0,\phi}\bm{F}$, with the same unimodular matrices
$\bm{E}$ and $\bm{F}$.

The reversal111This definition, which is standard, is particularly appropriate
for the monomial basis. The coefficients of the reversed matrix polynomial in
the monomial basis are simply the same matrices in reverse order. The notion
of a reversal, however, is independent of the basis used. In [8] we find a
slightly different definition, appropriate for computation in a Lagrange or
Hermite interpolational basis, which maps an arbitrary finite point to
infinity; this difference allows for greater numerical stability. of a matrix
polynomial of grade $\ell$ is the polynomial
$\mathrm{rev}\,\bm{P}(\lambda)=\lambda^{\ell}\bm{P}(\lambda^{-1})$. A
linearization $\bm{L}(\lambda)=\lambda\bm{C}_{1}-\bm{C}_{0}$ of $\bm{P}$ is
called a strong linearization if also $\mathrm{rev}\,\bm{L}(\lambda)$ is a
linearization of $\mathrm{rev}\,\bm{P}(\lambda)$.

For more information, consult [24]. See also [26], [18], and consult the
seminal book [17]. Linearizations using different polynomials bases were first
systematically studied in [1]. Some recent papers of interest include [4],
[25], [10], [14], and [29]; this is a very active area.

###  1.1 Organization of the Paper

In the remainder of this first section, we establish notation, give the
definition of a generalized standard triple, and give lemmas about common
similarity transformations222Transposition and flipping give altogether four
common variations of companion matrix pencils. Other variations are possible
(indeed, any similarity transformation will work) but these are the main
variations seen in the literature. We include these variations in enough
detail to help the reader with the bookkeeping.. At the end of this section,
we show how to use the generalized standard triple in the construction of
algebraic linearizations.

In section 2, we tabulate our results in detail for the generalized standard
triples for many common polynomial bases in sections 3.2, 3.3, and 3.4. In
section 4, we give the strict equivalence of generalized standard triples for
any polynomial basis. We give specific proofs in section 5.

###  1.2 Notation and Definition of a Generalized Standard Triple

If $\bm{L}(z)=z\bm{C}_{1}-\bm{C}_{0}\in\mathbb{C}^{N\times N}$ (usually
$N=n\ell$ but not always; for Lagrange and Hermite interpolational bases some
constructions use $N=(n+2)\ell$ and others $N=(n+1)\ell$) is a linearization
of $\bm{P}(z)$, then
$\mathrm{det}(\bm{P}(z))=\mathrm{det}(\bm{L}(z))=\mathrm{det}(z\bm{C}_{1}-\mathrm{\bm{C}_{0})}$.
The eigenvalues of $\bm{P}$ are thus computable from the generalized
eigenvalues of $\bm{L}$. For instance, this can be done with eig(C0,C1) in
Matlab, or Eigenvalues(C[0],C[1]) in Maple if the matrix variables are defined
appropriately (and are of complex floating-point type in Maple).

A standard pair $(\bm{X},\bm{T})$ for a regular matrix polynomial
$\bm{P}(\lambda)$ expressed in the monomial basis is defined in [17] or in
[24] as having the properties $\bm{X}$ has dimension $n\times n\ell$, $\bm{T}$
has dimension $n\ell\times n\ell$, and

(1.2) $\sum_{k=0}^{\ell}\bm{P}_{k}\bm{X}\bm{T}^{k}=0$

and that the $n\ell$ by $n\ell$ matrix

(1.3) $\bm{Q}=\begin{bmatrix}\bm{X}\\\ \bm{X}\bm{T}\\\ \vdots\\\
\bm{X}\bm{T}^{\ell-1}\end{bmatrix}$

is nonsingular. We can then define a third matrix

(1.4) $\bm{Y}=\bm{Q}^{-1}\begin{bmatrix}\bm{0}_{n}\\\ \vdots\\\ \bm{0}_{n}\\\
\bm{I}_{n}\end{bmatrix}$

and say that the triple $(\bm{X},\bm{T},\bm{Y})$ is a standard triple for a
monic $\bm{P}(\lambda)$. It is pointed out in [24] that monicity of
$\bm{P}(\lambda)$ is not required for many of the formulæ to do with standard
pairs (but is required for some). Theorem 2.6 of [17], which states that if
there are matrices $\bm{X}$, $\bm{T}$, and $\bm{Y}$ of dimension $n\times
n\ell$, $n\ell\times n\ell$, and $n\ell\times n$ for which

(1.5) $\bm{P}^{-1}(\lambda)=\bm{X}(\lambda\bm{I}_{n}-\bm{T})^{-1}\bm{Y}$

then $(\bm{X},\bm{T},\bm{Y})$ is a standard triple for $\bm{P}(\lambda)$. This
construction is clearly tied to the monomial basis, and we would like to
extend this to a form for other bases, and also to the non-monic case. In
particular we would like the following extension of Theorem 2.6 in [17] or
Theorem 12.1.4 in [16] to be available: If a matrix
$\bm{X}\in\mathbb{C}^{n\times N}$, the pencil $\bm{L}(z)$, and a matrix
$\bm{Y}\in\mathbb{C}^{N\times n}$ satisfies

(1.6) $\bm{P}^{-1}(z)=\bm{X}(z\bm{C}_{1}-\bm{C}_{0})^{-1}\bm{Y}$

for $z\notin\Lambda(\bm{P})$ (the set of polynomial eigenvalues of $\bm{P}$),
then $\bm{X}$, $\bm{L}(z)$, and $\bm{Y}$ form a generalized standard triple
for $\bm{P}(\lambda)$. This obviously requires regularity of $\bm{P}$.

This new definition would allow $N>n\ell$ and not just $N=n\ell$. Note that
the matrices $\bm{X}$ and $\bm{Y}$ do not depend on $z$, but the linearization
$\bm{L}(z)$ does, albeit only linearly; we could instead have chosen to use
the words “standard quadruple” to mean $(\bm{X},\bm{C}_{1},\bm{C}_{0},\bm{Y})$
where $z$ does not appear of any of these matrices, but this seems to be a
matter of aesthetics only.

In the case when the leading coefficient $\bm{A}_{\ell}$ of the matrix
polynomial is nonsingular, one can reduce to the monic case in any of several
ways. This is done in several places in the literature, and thus our
generalized standard triple is not very new in allowing for a non-identity
matrix coefficient of $\lambda$ in the linearization. For various reasons we
do not wish to do this here; for instance, if the polynomial basis is not
degree-graded, e.g. for the Bernstein basis, then it is not clear what one
should use for the leading coefficient. Similarly for the Lagrange and Hermite
interpolational bases.

All we will need for the purposes of this paper is that the representation
displayed in equation (1.6) holds (whether we say that it involves a “triple”
or a “quadruple”). The representation itself is what is useful in the
recursive construction of algebraic linearizations.

There is some risk of confusion with this definition: not all properties of
standard triples may hold for “generalized standard triples”, and so we should
take some care.

####  1.2.1 Similarity

If $\bm{X}$, $z\bm{C}_{1}-\bm{C}_{0}$, and $\bm{Y}$ form a generalized
standard triple, then so also do $\bm{X}\bm{S}$,
$\bm{S}^{-1}(z\bm{C}_{1}-\bm{C}_{0})\bm{S}$, and $\bm{S}^{-1}\bm{Y}$ for any
nonsingular matrix $\bm{S}$ of dimension $N$ by $N$.

######  Lemma 1.1.

If $\bm{S}$ is nonsingular333We only use similarities in this section, but as
a referee points out, we could also use equivalences, with different matrices
$\bm{S}$ on the left and right. and $\bm{B}_{1}=\bm{S}^{-1}\bm{C}_{1}\bm{S}$
and $\bm{B}_{0}=\bm{S}^{-1}\bm{C}_{0}\bm{S}$ so that the pencil
$z\bm{B}_{1}-\bm{B}_{0}$ has the same generalized eigenvalues as
$z\bm{C}_{1}-\bm{C}_{0}$, then another standard triple for $\bm{P}(z)$ is
$\widetilde{\bm{X}}$, $z\bm{B}_{1}-\bm{B}_{0}$, $\widetilde{\bm{Y}}$ where
$\widetilde{\bm{X}}=\bm{XS}$ and $\widetilde{\bm{Y}}=\bm{S}^{-1}\bm{Y}$.

######  Proof 1.2.

$\displaystyle\bm{P}^{-1}(z)$
$\displaystyle=\bm{X}(z\bm{C}_{1}-\bm{C}_{0})^{-1}\bm{Y}$
$\displaystyle=\bm{XSS}^{-1}(z\bm{C}_{1}-\bm{C}_{0})^{-1}\bm{SS}^{-1}\bm{Y}$
(1.7)
$\displaystyle=(\bm{XS})(\bm{S}^{-1}(z\bm{C}_{1}-\bm{C}_{0})\bm{S})^{-1}\bm{S}^{-1}\bm{Y}\>.$

Several similarities are used very frequently. For convenience we describe two
of the most common explicitly here.

######  Lemma 1.3 (Flipping).

Put $\bm{J}=\text{the}\;N\times N$ “anti-identity”, also called the sip
matrix, for standard involutory permutation, $\bm{J}_{i,j}=0$ unless $i+j=N+1$
when $\bm{J}_{i,N+1-i}\leavevmode\nobreak\ =\leavevmode\nobreak\ 1$. Then
$\bm{J}^{2}\leavevmode\nobreak\ =\leavevmode\nobreak\ \bm{I}$ and the
“flipped” linearization $\bm{L_{F}}(z)=\bm{J}(z\bm{C}_{1}-\bm{C}_{0})\bm{J}$
has in its generalized standard triple the matrices $\bm{X_{F}}=\bm{XJ}$ and
$\bm{Y_{F}}=\bm{JY}$. The paper [27] calls this matrix “$\bm{R}$”.

######  Proof 1.4.

Immediate.

######  Remark 1.5.

Flipping switches both the order of the equations and the order of the
variables. It obviously does not change eigenvalues. Flipping, transposition,
and flipping-with-transposition give four common equivalent linearizations
[31].

###  1.3 Algebraic Linearizations

An algebraic linearization, as referred to in the title of this present note,
is defined in [6] as a linearization of a matrix polynomial
$\bm{H}(\lambda)=z\bm{A}(\lambda)\bm{B}(\lambda)+\bm{C}$ constructed
recursively from linearizations of the lower-degree component matrix
polynomials $\bm{A}(\lambda)$ and $\bm{B}(\lambda)$, together with a constant
matrix $\bm{C}$. Algebraic linearizations are typically strong linearizations,
which not only preserve eigenvalues but also their partial multiplicities,
even at infinity [6]. Algebraic linearizations offer a new, potentially more
numerically stable, class of linearizations. The recursive construction of
algebraic linearizations relies on the generalized standard triples of each of
the component matrix polynomials, and (as does the unrelated paper [29])
allows different polynomial bases to be used for each component. This present
note provides some explicit formulas for generalized standard triples in
various bases, for reference. As one reviewer points out, these formulas could
simply be obtained by reading the proofs that these linearizations are indeed
linearizations; one purpose of this paper is simply convenience.

If $\bm{A}(z)$ and $\bm{B}(z)$ have the generalized standard triple
representations
$\bm{A}^{-1}(z)=\bm{X}_{A}(z\bm{D}_{A}-\bm{E}_{A})^{-1}\bm{Y}_{A}$ and
$\bm{B}^{-1}(z)=\bm{X}_{B}(z\bm{D}_{B}-\bm{E}_{B})^{-1}\bm{Y}_{B}$, then the
pencil $z\bm{D}_{H}-\bm{E}_{H}$ is a linearization of
$\bm{H}(z)=z\bm{A}(z)\bm{B}(z)+\bm{C}$, where the matrices $\bm{D}_{H}$ and
$\bm{E}_{H}$ are given as follows:

(1.8) $\bm{D}_{H}=\left[\begin{array}[]{ccc}\bm{D}_{A}&&\\\ &\bm{I}_{n}&\\\
&&\bm{D}_{B}\end{array}\right]$

and

(1.9)
$\bm{E}_{H}=\left[\begin{array}[]{ccc}\bm{A}&\bm{0}_{N_{A},n}&-\bm{Y}_{A}\bm{C}\bm{X}_{B}\\\
-\bm{X}_{A}&\bm{0}_{n}&\bm{0}_{n,N_{B}}\\\
\bm{0}_{N_{B},N_{A}}&-\bm{Y}_{B}&\bm{B}\end{array}\right]$

For a proof and some examples, see [6].

##  2 Representation of Matrix Polynomials using Generalized Standard Triples

Since the $\phi_{k}(x)$, $0\leq k\leq\ell$ form a basis, we may express the
polynomial $1$ in that basis: then $1=\sum_{k=0}^{n-1}e_{k}\phi_{k}(z)$
defines the coefficients $e_{k}$ uniquely. Putting

(2.10)
$\bm{X}=\begin{bmatrix}e_{n-1}&e_{n-2}&\cdots&e_{1}&e_{0}\end{bmatrix}\otimes\bm{I}$

always gives our generalized standard triple
$\bm{P}^{-1}(z)=\bm{X}(z\bm{C}_{1}-\bm{C}_{0})^{-1}\bm{Y}$ with
$\bm{Y}=\begin{bmatrix}\bm{I}_{n}&\bm{0}_{n}&\bm{0}_{n}&\cdots&\bm{0}_{n}\end{bmatrix}^{T}$.
The proof is simple and (with appropriate modifications for the bases)
universal: the change of basis matrix $\bm{\Phi}$ has as its last row of its
inverse the coefficients of $1$ in that basis. This is exactly the translation
of the monomial standard triple element $\bm{X}$ into the new basis.

$\displaystyle\bm{P}^{-1}(z)$
$\displaystyle=\bm{X}\left(z\bm{C}_{1}-\bm{C}_{0}\right)^{-1}\bm{Y}$
$\displaystyle=\bm{X}\left(z\bm{E}\bm{C}_{1,m}\bm{F}-\bm{E}\bm{C}_{0,m}\bm{F}\right)^{-1}\bm{Y}$
(2.11)
$\displaystyle=\bm{X}\bm{\Phi}^{-1}\left(z\bm{C}_{1,m}-\bm{C}_{0,m}\right)^{-1}\bm{E}^{-1}\bm{Y}\>.$

Thus the new $\bm{X}$ is as claimed. It is a separate matter to show that
$\bm{E}^{-1}[\bm{I},0,\ldots,0]^{T}=[\bm{I},0,\ldots,0]^{T}$ again, but it
always is for the bases that we consider.

######  Remark 2.1.

There are linearizations not explicitly considered in this paper; for
instance, a referee has pointed out that when a matrix polynomial is expressed
in a basis where the elements satisfy a linear recurrence, then there is an
automatic way to build what is called a CORK linearization. See [18] and [32]
for details. Whatever the linearization, though, the previous universal
theorem shows how to construct the generalized standard triple.

That proof sketch may not be convincing, in part because details are omitted.
In what follows we examine specific cases in detail and supply specific proofs
for each basis. Indeed, much of the utility of this paper is simply writing
down those details, which will allow easier programming for the uses of these
generalized standard triples.

##  3 Examples of generalized standard triples

In this section, we tabulate generalized standard triples for four classes of
linearizations. We do so by examples of companion pencils, leaving the reader
to do the necessary tensor products to produce linearizations. This saves some
space in the presentation. In contrast, in section 5 where we gave proofs, we
do so in full generality.

###  3.1 Companion matrices

In the special case $n=1$, a linearization is usually called a “companion
pencil” or Frobenius pencil444The Frobenius form of a matrix is related, but
different: see for instance [30].. Thus finding roots of a scalar polynomial
can be done by finding generalized eigenvalues of the companion pencil. In the
monic case, $\bm{C}_{1}$ becomes the identity matrix and the generalized
eigenproblem becomes a standard eigenproblem. Kublanovskaya calls these
“accompanying pencils” in [22]. For bases other than the monomial, the
unfortunate nomenclature “colleague matrix” or “comrade matrix” is also used.
This nomenclature hinders citation search and we prefer “generalized
companion”, if a distinction is needed. See [26].

Construction of a linearization from a companion pencil is a simple matter of
the Kronecker (tensor) product: given $\bm{C}_{1}$,
$\bm{C}_{0}\in\mathbb{C}^{n\times n}$, take
$\widetilde{\bm{C}_{1}}=\bm{C}_{1}\otimes\bm{I}_{n}$ and then replace each
block $p_{k}\bm{I}_{n}$ with the corresponding matrix coefficient
$\bm{P}_{k}\in\mathbb{C}^{r\times r}$ (the first $p_{k}$, in
$p_{k}\bm{I}_{n}$, is the symbolic coefficient from
$p(z)=\sum_{k=0}^{\ell}p_{k}\phi_{k}(z)$; the matrix coefficient
$\bm{P}_{k}\in\mathbb{C}^{r\times r}$ is from
$\bm{P}(z)=\sum_{k=0}^{\ell}\bm{P}_{k}\phi_{k}(z).)$ This will be clearer by
example.

###  3.2 Bases with three-term recurrence relations

The monomial basis, the shifted monomial basis, the Taylor basis, the Newton
interpolational bases, and many common orthogonal polynomial bases all have
three-term recurrence relations that, except for initial cases, can be written

(3.12)
$z\phi_{k}(z)=\alpha_{k}\phi_{k+1}(z)+\beta_{k}\phi_{k}(z)+\gamma_{k}\phi_{k-1}(z)\>.$

We give a selection in table 1, and refer the reader to Section 18.9 of the
Digital Library of Mathematical Functions (dlmf.nist.gov) for more. See also
[15].

$\phi_{k}(z)$ | Name | $\alpha_{k}$ | $\beta_{k}$ | $\gamma_{k}$ | $\phi_{0}$ | $\phi_{1}$
---|---|---|---|---|---|---
$z^{k}$ | monomial | $1$ | $0$ | $0$ | $1$ | $z$
$(z-a)^{k}$ | shifted monomial | $1$ | $a$ | $0$ | $1$ | $z-a$
$\nicefrac{{(z-a)^{k}}}{{k!}}$ | Taylor | $n+1$ | $a$ | $0$ | $1$ | $z-a$
$\prod_{j=0}^{k-1}(z-\tau_{j})$ | Newton interpolational | $1$ | $\tau_{n}$ | $0$ | $1$ | $z-\tau_{0}$
$T_{k}(z)=\cos\left(k\cos^{-1}(z)\right)$ | Chebyshev | $\nicefrac{{1}}{{2}}$ | $0$ | $\nicefrac{{1}}{{2}}$ | $1$ | $z$
$P_{k}(z)$ | Legendre | $\nicefrac{{(k+1)}}{{(2k+1)}}$ | $0$ | $\nicefrac{{k}}{{(2k+1)}}$ | $1$ | $z$
Table 1: A short list of three-term recurrence relations for some important
polynomial bases. For a more comprehensive list, see The Digital Library of
Mathematical Functions. These relations and others are coded in Walter
Gautschi’s packages OPQ and SOPQ [15] and in the MatrixPolynomialObject
implementation package in Maple (see [19]).

For all such bases, we have the companion pencil555For exposition, we follow
Peter Lancaster’s dictum, namely that the $5\times 5$ case almost always gives
the idea.

(3.13) $\displaystyle\bm{C}_{1}$
$\displaystyle=\begin{bmatrix}\dfrac{p_{5}}{\alpha_{4}}&&&&\\\ &1&&&\\\
&&1&&\\\ &&&1&\\\ &&&&1\end{bmatrix}$ (3.14) $\displaystyle\bm{C}_{0}$
$\displaystyle=\begin{bmatrix}-p_{4}+\dfrac{\beta_{4}}{\alpha_{4}}p_{5}&-p_{3}+\dfrac{\gamma_{4}}{\alpha_{4}}p_{5}&-p_{2}&-p_{1}&-p_{0}\\\
\alpha_{3}&\beta_{3}&\gamma_{3}&&\\\ &\alpha_{2}&\beta_{2}&\gamma_{2}&\\\
&&\alpha_{1}&\beta_{1}&\gamma_{1}\\\ &&&\alpha_{0}&\beta_{0}\end{bmatrix}$

and

(3.15) $\displaystyle\bm{X}$
$\displaystyle=\begin{bmatrix}0&0&0&0&1\end{bmatrix}$ (3.16)
$\displaystyle\bm{Y}$
$\displaystyle=\begin{bmatrix}1&0&0&0&0\end{bmatrix}^{\mathrm{T}}\>.$

For instance, a flipped and transposed pencil of this class for the Chebyshev
case is666For the matrix polynomial case, each $\bm{P}_{k}$ would be
transposed.

(3.17) $\bm{L}(z)=\begin{bmatrix}z&-\dfrac{1}{2}&&&p_{0}\\\
-1&z&-\dfrac{1}{2}&&p_{1}\\\ &-\dfrac{1}{2}&z&-\dfrac{1}{2}&p_{2}\\\
&&-\dfrac{1}{2}&z&p_{3}+p_{5}\\\ &&&-\dfrac{1}{2}&2zp_{5}+p_{4}\end{bmatrix}$

has flipped and transposed $\bm{X}=\begin{bmatrix}0&0&0&0&1\end{bmatrix}$,
$\bm{Y}=\begin{bmatrix}1&0&0&0&0\end{bmatrix}^{\mathrm{T}}$. As another
instance, a Newton interpolational basis on the nodes $\tau_{0}$, $\tau_{1}$,
$\ldots$, $\tau_{5}$ has a companion pencil

(3.18) $z\begin{bmatrix}p_{5}&&&&\\\ &1&&&\\\ &&1&&\\\ &&&1&\\\
&&&&1\end{bmatrix}-\begin{bmatrix}-p_{4}+\tau_{4}p_{5}&-p_{3}&-p_{2}&-p_{1}&-p_{0}\\\
1&\tau_{3}&&&\\\ &1&\tau_{2}&&\\\ &&1&\tau_{1}&\\\
&&&1&\tau_{0}\end{bmatrix}\>.$

The corresponding linearization is

(3.19) $z\begin{bmatrix}\bm{P}_{5}&&&&\\\ &\bm{I}_{n}&&&\\\ &&\bm{I}_{n}&&\\\
&&&\bm{I}_{n}&\\\
&&&&\bm{I}_{n}\end{bmatrix}-\begin{bmatrix}-\bm{P}_{4}+\tau_{4}\bm{P}_{5}&-\bm{P}_{3}&-\bm{P}_{2}&-\bm{P}_{1}&-\bm{P}_{0}\\\
\bm{I}_{n}&\tau_{3}\bm{I}_{n}&&&\\\ &\bm{I}_{n}&\tau_{2}\bm{I}_{n}&&\\\
&&\bm{I}_{n}&\tau_{1}\bm{I}_{n}&\\\
&&&\bm{I}_{n}&\tau_{0}\bm{I}_{n}\end{bmatrix}\>.$

###  3.3 The Bernstein basis

The set of polynomials $\\{B_{k}^{\ell}(z)\\}_{k=0}^{\ell}$ is a set of
$\ell+1$ polynomials each of exact degree $\ell$ that forms a basis for
polynomials of degree at most $\ell$. They have many applications, for example
in Computer Aided Geometric Design (CAGD), and many important properties
including that of optimal condition number over all bases positive on
$\left[0,1\right]$. They do not satisfy a simple three term recurrence
relation of the form discussed in section 3.2. See [13], [12], and [11] for
more details of Bernstein bases.

A companion pencil for $p_{5}(z)=\sum_{k=0}^{5}p_{k}B_{k}^{5}(z)$ is

(3.20) $\displaystyle\bm{C}_{1}$
$\displaystyle=\begin{bmatrix}-p_{4}+\dfrac{1}{5}p_{5}&-p_{3}&-p_{2}&-p_{1}&-p_{0}\\\
1&\dfrac{2}{4}&&&\\\ &1&\dfrac{3}{3}&&\\\ &&1&\dfrac{4}{2}&\\\
&&&1&\dfrac{5}{1}\end{bmatrix}$ (3.21) $\displaystyle\bm{C}_{0}$
$\displaystyle=\begin{bmatrix}-p_{4}&-p_{3}&-p_{2}&-p_{1}&-p_{0}\\\ 1&0&&&\\\
&1&0&&\\\ &&1&0&\\\ &&&1&0\end{bmatrix}$ (3.22) $\displaystyle\bm{X}$
$\displaystyle=\begin{bmatrix}\dfrac{1}{5}&\dfrac{2}{5}&\dfrac{3}{5}&\dfrac{4}{5}&\dfrac{5}{5}\end{bmatrix}$
(3.23) $\displaystyle\bm{Y}$
$\displaystyle=\begin{bmatrix}1&0&0&0&0\end{bmatrix}^{\mathrm{T}}\>.$

We have $p^{-1}(z)=\bm{X}(z\bm{C}_{1}-\bm{C}_{0})^{-1}\bm{Y}$ if $p(z)\neq 0$.
This pencil was first analyzed in [20] and [21]. One of the present authors
independently invented and implemented a version of this linearization in
Maple (except using $\bm{P}^{\mathrm{T}}(z)$, and reversed from the above
form) in about $2004$. For a review of Bernstein linearization, see [27]. For
a proof of its numerical stability, see the original thesis [20]. We supply a
proof in section 5. The standard triple is, we believe, new to this paper.

###  3.4 The Lagrange interpolational basis

There are by now several Lagrange basis pencils and linearizations. The use of
barycentric forms means that Lagrange interpolation is efficient and
numerically stable. For many sets of nodes (Chebyshev nodes on $[-1,1]$, or
roots of unity on the unit disk) the resulting interpolant is also well-
conditioned, and can even be “better than optimal” [9], see also [5]. The
linearization we use here is “too large” and has (numerically harmless in our
experience) spurious roots at infinity777This numerical harmlessness needs
some explanation. In brief, Lagrange basis matrix polynomial eigenvalues will
be well-conditioned only in a compact region determined by the interpolation
nodes, and are increasingly ill-conditioned towards infinity; in practice this
means only small changes in the data are needed to perturb large finite ill-
conditioned eigenvalues out to infinity. Any eigenvalues produced numerically
that are well outside the region determined by the interpolation nodes are
likely easily perturbed all the way to infinity, and can be safely ignored.;
for alternative formulations see [32], [28]. Then pencil is
$z\bm{C}_{1}-\bm{C}_{0}$ where

(3.24) $\displaystyle\bm{C}_{1}$ $\displaystyle=\begin{bmatrix}0&&&&\\\
&1&&&\\\ &&1&&\\\ &&&1&\\\ &&&&1\end{bmatrix}$ (3.25)
$\displaystyle\bm{C}_{0}$
$\displaystyle=\begin{bmatrix}0&-\rho_{0}&-\rho_{1}&-\rho_{2}&-\rho_{3}&-\rho_{4}\\\
\beta_{0}&\tau_{0}&&&&\\\ \beta_{1}&&\tau_{1}&&&\\\ \beta_{2}&&&\tau_{2}&&\\\
\beta_{3}&&&&\tau_{3}&\\\ \beta_{4}&&&&&\tau_{4}\end{bmatrix}\>.$

Then $\mathrm{det}(\tau_{k}\bm{C}_{1}-\bm{C}_{0})=\rho_{k}$, $0\leq k\leq 4$
and $\mathrm{deg}(z\bm{C}_{1}-\bm{C}_{0})\leq 4$. Thus,
$p(z)=\mathrm{det}(z\bm{C}_{1}-\bm{C}_{0})$ interpolates the given data,
assuming the $\tau_{k}$ are distinct. Here the barycentric weights $\beta_{k}$
are found by partial fraction expansion of $\omega(z)^{-1}$ where

(3.26)
$\omega(z)=(z-\tau_{0})(z-\tau_{1})(z-\tau_{2})(z-\tau_{3})(z-\tau_{4})$

is the node polynomial. Explicitly,

(3.27) $\dfrac{1}{\omega(z)}=\sum_{k=0}^{5}\dfrac{\beta_{k}}{z-\tau_{k}}$

so

(3.28) $\beta_{k}=\prod_{\begin{subarray}{c}j=0\\\ j\neq
k\end{subarray}}^{5}(\tau_{k}-\tau_{j})^{-1}\>.$

The $\bm{X}$ and $\bm{Y}$ for the standard triple are

(3.29) $\displaystyle\bm{X}$
$\displaystyle=\begin{bmatrix}0&1&1&1&1&1\end{bmatrix}\>,$ (3.30)
$\displaystyle\bm{Y}$
$\displaystyle=\begin{bmatrix}1&0&0&0&0&0\end{bmatrix}^{\mathrm{T}}$

Notice in this case that for the linearization $N=(\ell+2)n$ while
$\mathrm{deg}\ p\leq\ell$, and therefore there are at least $2n$ eigenvalues
at infinity. This can be inconvenient if $n$ is at all large.

###  3.5 Hermite interpolational basis

The companion pencil of the previous section has been extended to Hermite
interpolational bases, where some of the nodes have “flowed together,”
collapsing to fewer distinct nodes888A formal definition can be found in [7],
for instance. The essential idea is that given two distinct pieces of data
$(\tau_{k},p(\tau_{k}))$ and $(\tau_{k+1},p(\tau_{k+1}))$, we also know the
forward difference $(p_{k+1}-p_{k})/(\tau_{k+1}-\tau_{k})$. In the limit as
one node approaches (flows towards) the other, we still know two pieces of
information: $p(\tau_{k})$ and $p^{\prime}(\tau_{k})$. Hermite interpolation
captures this idea.. We suppose that at each node $\tau_{i}$, there are now
$s_{i}\geq 1$ consecutive pieces of information known, namely
$\bm{P}(\tau_{i})$, $\bm{P}^{\prime}(\tau_{i})/1!$,
$\bm{P}^{\prime\prime}(\tau_{i})/2!$, and so on up to the last one, the value
of the $s_{i}-1$-th derivative at $z=\tau_{i}$, namely
$\bm{P}^{(s_{i}-1)}(\tau_{i})/(s_{i}-1)!$. The integer $s_{i}$ is called the
confluency of the node. The known pieces of information are the local Taylor
coefficients of the polynomial fitting the data:

(3.31) $\rho_{i,j}=\dfrac{f^{(j)}(\tau_{i})}{j!}\>,\quad 0\leq j\leq
s_{i}-1\>.$

Note that the derivative $\bm{P}^{\prime}(z)$ of a matrix polynomial is a
straightforward extension to matrices of the ordinary derivative. It is
isomorphic to the matrix with entries that are the ordinary derivatives of the
original matrix.

The companion pencil (that is, the scalar case) of the previous section
changes to the following elegant form. The matrix $\bm{C}_{1}$ is unchanged,

(3.32) $\bm{C}_{1}=\begin{bmatrix}0&&&&\\\ &1&&&\\\ &&\ddots&&\\\ &&&1&\\\
&&&&1\end{bmatrix}\>,$

being $(\ell+2)$ by $(\ell+2)$ as before, although now

(3.33) $\ell=-1+\sum_{i=0}^{m}s_{i}$

is the grade of the resulting polynomial. The matrix $\bm{C}_{0}$ changes,
picking up transposed Jordan-like blocks for each distinct node. For instance,
suppose we have two distinct nodes, $\tau_{0}$ and $\tau_{1}$. Suppose further
that $\tau_{0}$ has confluency $s_{0}=3$ while $\tau_{1}$ has confluency
$s_{1}=2$. This means that we know $f(\tau_{0})$,
$\nicefrac{{f^{\prime}(\tau_{0})}}{{1!}}$,
$\nicefrac{{f^{\prime\prime}(\tau_{0})}}{{2!}}$, $f(\tau_{1})$ and
$\nicefrac{{f^{\prime}(\tau_{1})}}{{1!}}$. Then,

(3.34)
$\bm{C}_{0}=\left[\begin{array}[]{cccccc}0&-\nicefrac{{f^{\prime\prime}(\tau_{0})}}{{2!}}&-\nicefrac{{f^{\prime}(\tau_{0})}}{{1!}}&-f(\tau_{0})&-\nicefrac{{f^{\prime}(\tau_{1})}}{{1!}}&-f(\tau_{1})\\\
\cdashline{2-4}\beta_{02}&\lx@intercol\hfil\tau_{0}\hfil\lx@intercol&&\lx@intercol\hfil\hfil\lx@intercol&&\\\
\beta_{01}&\lx@intercol\hfil
1\hfil\lx@intercol&\tau_{0}&\lx@intercol\hfil\hfil\lx@intercol&&\\\
\beta_{00}&\lx@intercol\hfil\hfil\lx@intercol&1&\lx@intercol\hfil\tau_{0}\hfil\lx@intercol&&\\\
\cdashline{2-6}\beta_{11}&&&&\lx@intercol\hfil\tau_{1}\hfil\lx@intercol&\lx@intercol\hfil\hfil\lx@intercol\\\
\beta_{10}&&&&\lx@intercol\hfil
1\hfil\lx@intercol&\lx@intercol\hfil\tau_{1}\hfil\lx@intercol\\\
\cdashline{5-6}\end{array}\right]$

Note the reverse ordering of the derivative values in this formulation. The
barycentric weights $\beta_{ij}$ again come from the partial fraction
expansion of the reciprocal of the node polynomial

(3.35) $\omega(z)=\prod_{i=0}^{m}(z-\tau_{i})^{s_{i}}\>.$

That is,

(3.36)
$\dfrac{1}{\omega(z)}=\sum_{i=0}^{m}\sum_{j=0}^{s_{i}-1}\dfrac{\beta_{ij}}{(z-\tau_{i})^{j+1}}\>.$

For the standard triple, take in the scalar case

(3.37) $\bm{Y}=\begin{bmatrix}1&0&\cdots&0\end{bmatrix}^{\mathrm{T}}$

but for $\bm{X}$ take the coefficients of the expansion of the polynomial $1$
in this particular Hermite interpolational basis: it is equal to $1$ at each
node but has all derivatives zero at each node. That is, put

(3.38) $\begin{cases}\rho_{ij}=1\quad&\text{if }j=0\>,\\\
0\quad&\text{otherwise}\>,\end{cases}$

and sort them in order:

(3.39)
$\bm{X}=\begin{bmatrix}0&\rho_{0,s_{0}-1}&\rho_{0,s_{0}-2}&\cdots&\rho_{0,0}&\rho_{1,s_{1}-1}&\cdots&\rho_{n,0}\end{bmatrix}\>.$

For the earlier instance (two nodes, of confluency 3 and 2, respectively,

(3.40)
$\bm{X}=\begin{bmatrix}0&\smash[b]{\underbrace{\begin{matrix}0&0&1\end{matrix}}_{\text{for
$\tau_{0}$}}}&\smash[b]{\underbrace{\begin{matrix}0&1\end{matrix}}_{\text{for
$\tau_{1}$}}}\end{bmatrix}\>.$

Then

(3.41) $p^{-1}(z)=\bm{X}(z\bm{C}_{1}-\bm{C}_{0})^{-1}\bm{Y}\>.$

######  Remark 3.1.

We may re-order the nodes in any fashion we like, and each ordering generates
its own companion pencil (both Hermite and Lagrange). We may also find a
pencil where the confluent data is ordered $p(\tau_{i})$,
$\nicefrac{{p^{\prime}(\tau_{i})}}{{1!}}$,
$\nicefrac{{p^{\prime\prime}(\tau_{i})}}{{2!}}$, etc., although we have not
done so.

If there is just one node of confluency $\ell$, we recover the standard
Frobenius companion form (plus two infinite roots):

(3.42) $\begin{bmatrix}0&&&&\\\ &1&&&\\\ &&\ddots&&\\\ &&&1&\\\
&&&&1\end{bmatrix}\quad,\quad\begin{bmatrix}0&-p_{\ell-1}&-p_{\ell-2}&\cdots&-p_{1}&-p_{0}\\\
1&\tau_{0}&&&&\\\ 0&1&\tau_{0}&&&\\\ 0&&1&\ddots&&\\\
\vdots&&&\ddots&\tau_{0}&\\\ 0&&&&1&\tau_{0}\end{bmatrix}\>.$

Here $p_{k}=\nicefrac{{p^{(k)}(\tau_{0})}}{{k!}}$ is the ordinary coefficient
in the expansion $p(z)=\sum_{k=0}^{\ell}p_{k}(z-\tau_{0})^{k}$. The numerical
stability of these Hermite interpolational companions has been studied briefly
[23] but much remains unknown. We confine ourselves in this paper to the study
of the standard triple.

To make a linearization out of these companion pencils, take the Kronecker
tensor product with $\bm{I}_{n}$, and insert the appropriate matrix polynomial
values and derivative values.

######  Remark 3.2.

The modified linearizations of [32] also have standard triples that can be
used for algebraic linearization, and arguably should be tabled here as well.
They have the advantage of including fewer eigenvalues at infinity, or no
spurious eigenvalues at infinity, which may lead to better algebraic
linearizations. However, they are more involved, and we have less numerical
experience with them. In particular we do not understand their dependence on
the ordering of the nodes, and so we leave their analysis to a future study.

##  4 Strict Equivalence of Generalized Standard Triples for any polynomial
basis

######  Theorem 4.1.

If $\phi_{k}(z)$ for $0\leq k\leq\ell$ is one of the degree-graded polynomial
bases (e.g. Chebyshev, Newton, Jacobi $P$) or a Bernstein basis although they
are not degree-graded, then the companion pencil for a polynomial $p(z)$
expressed in that basis is strictly equivalent to the second companion pencil
for same polynomial expressed in the monomial basis. That is, there exist
unimodular matrices $\bm{E}$ and $\bm{F}$ for which
$\bm{C}_{1,m}=\bm{E}\bm{C}_{1,\phi}\bm{F}$ and
$\bm{C}_{0,m}=\bm{E}\bm{C}_{0,\phi}\bm{F}$. The matrix $\bm{E}$ will depend on
the given polynomial $p(z)$. Here the subscript $m$ is short for “monomial".

######  Proof 4.2.

Denote the change-of-bases matrix for polynomials up to degree $\ell-1$ by
$\bm{\Phi}$. This matrix is $\ell$ by $\ell$. Then we have

(4.43) $\begin{bmatrix}\phi_{\ell-1}(z)\\\ \phi_{\ell-2}(z)\\\ \vdots\\\
\phi_{1}(z)\\\ \phi_{0}(z)\end{bmatrix}=\bm{\Phi}\begin{bmatrix}z^{\ell-1}\\\
\vdots\\\ z^{2}\\\ z\\\ 1\end{bmatrix}\>.$

In all cases considered here, the companion matrix pencil for a polynomial
$p(z)$ of exact degree $\ell$ has null vectors of the form

(4.44) $\bm{N}=\begin{bmatrix}\phi_{\ell-1}(\lambda)\\\
\phi_{\ell-2}(\lambda)\\\ \vdots\\\ \phi_{1}(\lambda)\\\
\phi_{0}(\lambda)\end{bmatrix}$

where $\lambda$ is a root of $p(z)$. That is,

(4.45) $\left(\lambda\bm{C}_{1,\phi}-\bm{C}_{0,\phi}\right)\bm{N}=\bm{0}\>.$

Using the $\bm{\Phi}$ formula above, we have $\bm{F}=\bm{\Phi}$. By direct
computation, we find that
$\bm{E}=\bm{C}_{0,m}\bm{\Phi}^{-1}\bm{C}_{0,\phi}^{-1}$ necessarily giving
$\bm{E}\bm{C}_{0,\phi}\bm{F}=\bm{C}_{0,m}$. Since $\bm{F}=\bm{\Phi}$ is
unimodular and nonsingular, all that remains is to show that $\bm{E}$ is
unimodular and nonsingular, and that it satisfies
$\bm{E}\bm{C}_{1,\phi}\bm{F}=\bm{C}_{1,m}$ as well. Since $\bm{C}_{1,\phi}$ is
the identity matrix except for the $1,1$ entry which is $a_{\ell}\neq 0$, the
leading coefficient of the polynomial, this last is straightforward. Indeed,
the action of premultiplying by $\bm{C}_{0,m}$ and postmultiplying by
$\bm{C}_{0,\phi}^{-1}$ cancels the coefficient $a_{\ell}$ in the $1,1$ entry
of $\bm{E}$, and thus by continuity this construction works even if
$a_{\ell}=0$, that is for polynomials of grade $\ell$. Then since $\bm{\Phi}$
is upper triangular for degree-graded matrices and lower triangular for
Bernstein matrices, and $\bm{E}$ has in that case zeros in the first column
below the diagonal entry, $\bm{E}$ is unimodular.

To make this work for regular matrix polynomials of dimension $n$, we must use
the tensor product $\bm{\Phi}\otimes\bm{I}$.

We illustrate this proof with a four by four example in the Bernstein basis.
If
$p(z)=a_{0}B_{0}^{4}(z)+a_{1}B_{1}^{4}(z)+a_{2}B_{2}^{4}(z)+a_{3}B_{3}^{4}(z)+a_{4}B_{4}^{4}(z)$,
then

(4.46)
$\bm{E}=\left[\begin{array}[]{cccc}1&1/4\,a_{{0}}-a_{{1}}+3/2\,a_{{2}}&-1/6\,a_{{0}}+2/3\,a_{{1}}&1/4\,a_{{0}}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&1/4&0&0\\\ \vskip 6.0pt plus 2.0pt
minus 2.0pt\cr 0&1/4&1/6&0\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0&1/4&1/3&1/4\end{array}\right]$

and

(4.47) $\bm{F}=\left[\begin{array}[]{cccc}4&0&0&0\\\ \vskip 6.0pt plus 2.0pt
minus 2.0pt\cr-6&6&0&0\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 4&-8&4&0\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr-1&3&-3&1\end{array}\right]\>.$

Direct computation shows that both matrices are nonsingular irrespective of
the values of the $a_{k}$, and that these transform the Bernstein companion
pencil

(4.48)
$\bm{C}_{0,\mathrm{Bernstein}}=\left[\begin{array}[]{cccc}-a_{{3}}&-a_{{2}}&-a_{{1}}&-a_{{0}}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr 1&0&0&0\\\ \vskip 6.0pt plus 2.0pt
minus 2.0pt\cr 0&1&0&0\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0&0&1&0\end{array}\right]$

and

(4.49)
$\bm{C}_{1,\mathrm{Bernstein}}=\left[\begin{array}[]{cccc}\dfrac{a_{{4}}}{4}-a_{{3}}&-a_{{2}}&-a_{{1}}&-a_{{0}}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr 1&2/3&0&0\\\ \vskip 6.0pt plus 2.0pt
minus 2.0pt\cr 0&1&3/2&0\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0&0&1&4\end{array}\right]$

to

(4.50)
$\bm{C}_{0,\mathrm{monomial}}=\left[\begin{array}[]{cccc}-b_{3}&-b_{2}&-b_{1}&-b_{{0}}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr 1&0&0&0\\\ \vskip 6.0pt plus 2.0pt
minus 2.0pt\cr 0&1&0&0\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0&0&1&0\end{array}\right]$

with $-b_{3}=4\,a_{{0}}-12\,a_{{1}}+12\,a_{{2}}-4\,a_{{3}}$,
$-b_{2}=-6\,a_{{0}}+12\,a_{{1}}-6\,a_{{2}}$, $-b_{1}=4\,a_{{0}}-4\,a_{{1}}$,
and $-b_{0}=-a_{0}$. Also, $\bm{E}\bm{C}_{1,\mathrm{Bernstein}}\bm{F}$ becomes
the identity matrix except the $1,1$ entry is
$b_{4}=a_{{0}}-4\,a_{{1}}+6\,a_{{2}}-4\,a_{{3}}+a_{{4}}$. These are the
correct coefficients of the same polynomial expressed in the monomial basis.

###  4.1 The Lagrange interpolational case

As is usual we denote a Lagrange basis element on the distinct nodes
$[\tau_{0},\tau_{1},\ldots,\tau_{\ell}]$ by $\ell_{k}(z)=\beta_{k}\prod_{j\neq
k}(z-\tau_{j})$. The use of the symbol $\ell$ by itself denotes an integer,
namely the grade of the polynomial; the distinction is that $\ell$ with a
subscript and a variable $\ell_{j}(z)$ denotes a Lagrange basis polynomial.
This should not cause confusion.

######  Remark 4.3.

Many people think of “interpolation” as meaning the construction of a monomial
basis polynomial $p(z)=a_{\ell}z^{\ell}+\cdots+a_{0}$ that fits the given data
$p(\tau_{k})=\rho_{k}$ for $0\leq k\leq\ell$. This is naive. Interpolation
truly means constructing a polynomial in any basis that we may use to evaluate
$p(z)$ for $z$ different to the values at the nodes. The most stable and
convenient way to do this is by the barycentric form of Lagrange interpolants.
Constructing an interpolant in a Newton basis by using divided differences or
the monomial basis by using a Vandermonde matrix is changing the basis.
Changing bases can have condition number exponential in the degree, and is
usually a bad idea. In practice, we use the barycentric form [3]. For the
purposes of proof of equivalence, we here occasionally use the Vandermonde
matrix, and we think about the explicit construction of the monomial basis.
This is not used in numerical practice.

######  Theorem 4.4.

If $\phi_{k}(z)$ for $0\leq k\leq\ell$ is a Lagrange basis on distinct nodes
$\tau_{0}$, $\tau_{1}$, $\ldots$, $\tau_{\ell}$, then the $\ell+2$ by $\ell+2$
companion pencil for a polynomial $p(z)$ expressed in that basis is strictly
equivalent to the second companion pencil for same polynomial expressed in the
monomial basis but regarded as having grade $\ell+2$ (i.e. with zero
coefficients padding the terms $z^{\ell+2}$ and $z^{\ell+1}$. That is, there
exist unimodular matrices $\bm{E}$ and $\bm{F}$ for which
$\bm{C}_{1,\phi}=\bm{E}\bm{C}_{1,m}\bm{F}$ and
$\bm{C}_{0,\phi}=\bm{E}\bm{C}_{0,m}\bm{F}$, where now the second companion
matrices in the monomial basis have dimension larger by two than needed for
the exact degree. The matrix $\bm{E}$ will depend on the given polynomial
$p(z)$.

######  Proof 4.5.

For the Lagrange basis companion pencil, the right null vector is of the form

(4.51) $\bm{N}=\begin{bmatrix}w(\lambda)\\\ \ell_{0}(\lambda)\\\
\ell_{1}(\lambda)\\\ \vdots\\\ \ell_{\ell}(\lambda)\end{bmatrix}\>.$

Here $w(z)=\prod_{k=0}^{\ell}(z-\tau_{k})$ is of degree $\ell+1$ and all the
other entries, being elements of the Lagrange basis on $\ell+1$ nodes, are of
degree $\ell$. Thus $\bm{\Phi}$ is dimension $\ell+2$ by $\ell+2$ and has
first column $e_{1}$; that is, $1$ in the first entry and zeros below it. The
rest of the first row contains the coefficients of $w(z)$ expanded in the
monomial basis. The remaining rows of $\bm{{\Phi}}$ contain the coefficients
of the monomial expansions of the Lagrange basis polynomial; that is, the
inverse of the transposed Vandermonde matrix, which relates the Lagrange
interpolation basis to the monomial basis. Call that block $\hat{\bm{\Phi}}$.
Explicitly,

(4.52)
$\hat{}\bm{\Phi}^{-1}=\begin{bmatrix}\tau_{0}^{\ell}&\tau_{1}^{\ell}&\cdots&\tau_{\ell}^{\ell}\\\
\tau_{0}^{\ell-1}&\tau_{1}^{\ell-1}&\cdots&\tau_{\ell}^{\ell-1}\\\
\vdots&&&\vdots\\\ \tau_{0}&\tau_{1}&\cdots&\tau_{\ell}\\\ 1&1&\cdots&1\\\
\end{bmatrix}\>.$

Then $\bm{E}$ is straightforwardly seen to be
$\mathrm{diag}(1,\hat{}\bm{\Phi}^{-1})$ and as in the degree-graded case
$\bm{F}=\bm{\Phi}$. Further, $\bm{E}\bm{C}_{0,\phi}\bm{F}$ has as its first
row $[0,-\rho\hat{}\bm{\Phi}]$. But $\rho\hat{}\bm{\Phi}$ is by the
Vandermonde matrix simply the negative of the vector of monomial coefficients,
$-[0,a_{\ell-1},a_{\ell-2},\ldots,a_{0}]$. More, the block underneath, namely
$\hat{}\bm{\Phi}^{-1}\mathrm{diag}(\tau_{0},\tau_{1},\ldots,\tau_{\ell})\hat{}\bm{\Phi}$
turns out to be simple, because

(4.53)
$\hat{}\bm{\Phi}^{-1}\mathrm{diag}(\tau_{0},\tau_{1},\ldots,\tau_{\ell})=\begin{bmatrix}\tau_{0}^{\ell+1}&\tau_{1}^{\ell+1}&\cdots&\tau_{\ell}^{\ell+1}\\\
\tau_{0}^{\ell}&\tau_{1}^{\ell}&\cdots&\tau_{\ell}^{\ell}\\\
\vdots&&&\vdots\\\ \tau_{0}^{2}&\tau_{1}^{2}&\cdots&\tau_{\ell}^{2}\\\
\tau_{0}&\tau_{1}&\cdots&\tau_{\ell}\\\ \end{bmatrix}$

Multiplying this by $\hat{}\bm{\Phi}$ shifts the identity matrix down one
diagonal, giving the correct form for the second companion matrix
$\bm{C}_{0,m}$.

As in the previous theorem, to construct $\bm{E}$ and $\bm{F}$ for
$n$-dimensional regular matrix polynomials, we must take the tensor product
$\bm{\Phi}\otimes\bm{I}$.

Again we illustrate this proof with an example, this time of interpolation at
the four points $[-1,-1/2,1/2,1]$. This will give rise to a polynomial of
degree at most $3$. If the values this polynomial takes at these four points
are $\rho_{0}$, $\rho_{1}$, $\rho_{2}$, and $\rho_{3}$, then the equivalent
polynomial expressed in the monomial basis has coefficients

$\displaystyle a_{0}$
$\displaystyle=-1/6\,\rho_{{0}}+2/3\,\rho_{{1}}+2/3\,\rho_{{2}}-1/6\,\rho_{{3}}$
$\displaystyle a_{1}$
$\displaystyle=1/6\,\rho_{{0}}-4/3\,\rho_{{1}}+4/3\,\rho_{{2}}-1/6\,\rho_{{3}}$
$\displaystyle a_{2}$
$\displaystyle=2/3\,\rho_{{0}}-2/3\,\rho_{{1}}-2/3\,\rho_{{2}}+2/3\,\rho_{{3}}$
(4.54) $\displaystyle a_{3}$
$\displaystyle=2/3\,\rho_{{3}}-4/3\,\rho_{{2}}+4/3\,\rho_{{1}}-2/3\,\rho_{{0}}\>.$

Expressing this as a polynomial of grade $5$, that is $p(z)=0\cdot
z^{5}+0\cdot z^{4}+a_{3}z^{3}+a_{2}z^{2}+a_{1}z+a_{0}$, we get the second
companion pencil

(4.55)
$\bm{C}_{0,\mathrm{monomial}}=\left[\begin{array}[]{ccccc}0&-a_{{3}}&-a_{{2}}&-a_{{1}}&-a_{{0}}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr 1&0&0&0&0\\\ \vskip 6.0pt plus 2.0pt
minus 2.0pt\cr 0&1&0&0&0\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0&0&1&0&0\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0&0&0&1&0\end{array}\right]$

and

(4.56) $\bm{C}_{1,\mathrm{monomial}}=\left[\begin{array}[]{ccccc}0&0&0&0&0\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&1&0&0&0\\\ \vskip 6.0pt plus 2.0pt
minus 2.0pt\cr 0&0&1&0&0\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0&0&0&1&0\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0&0&0&0&1\end{array}\right]\>.$

The Lagrange basis companion matrix of [9], which admittedly has two extra
infinite eigenvalues, is

(4.57)
$\bm{C}_{0,\mathrm{Lagrange}}=\left[\begin{array}[]{ccccc}0&-\rho_{{3}}&-\rho_{{2}}&-\rho_{{1}}&-\rho_{{0}}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr 2/3&1&0&0&0\\\ \vskip 6.0pt plus 2.0pt
minus 2.0pt\cr-4/3&0&1/2&0&0\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
4/3&0&0&-1/2&0\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr-2/3&0&0&0&-1\end{array}\right]$

and

(4.58) $\bm{C}_{1,\mathrm{Lagrange}}=\left[\begin{array}[]{ccccc}0&0&0&0&0\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&1&0&0&0\\\ \vskip 6.0pt plus 2.0pt
minus 2.0pt\cr 0&0&1&0&0\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0&0&0&1&0\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0&0&0&0&1\end{array}\right]\>.$

For these interpolation nodes, direct computation shows

(4.59) $\bm{E}=\left[\begin{array}[]{ccccc}1&0&0&0&0\\\ \vskip 6.0pt plus
2.0pt minus 2.0pt\cr 0&1&1/8&-1/8&-1\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0&1&1/4&1/4&1\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&1&1/2&-1/2&-1\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&1&1&1&1\end{array}\right]$

and

(4.60) $\bm{F}=\left[\begin{array}[]{ccccc}1&0&-5/4&0&1/4\\\ \vskip 6.0pt plus
2.0pt minus 2.0pt\cr 0&2/3&2/3&-1/6&-1/6\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr 0&-4/3&-2/3&4/3&2/3\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0&4/3&-2/3&-4/3&2/3\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0&-2/3&2/3&1/6&-1/6\end{array}\right]\>.$

##  5 Individual Proofs

In section 2 we give a short universal proof of all the theorems in this
section. Each individual proof in this section is therefore redundant. We
include them here both for surety (giving two proofs of each theorem) and
because they give insight and may be relevant to any numerical analysis. We
will use the Schur Complement, in the following form: assuming a matrix
$\bm{R}$ is partitioned into

(5.61) $\displaystyle\bm{R}=\begin{bmatrix}\bm{A}&\bm{B}\\\
\bm{C}&\bm{D}\end{bmatrix}$

where $\bm{A}\in\mathbb{C}^{{}_{\ell}\times_{\ell}}$,
$\bm{B}\in\mathbb{C}^{r\times(N-r)}$, $\bm{C}\in\mathbb{C}^{(N-r)\times r}$
and $\bm{D}\in\mathbb{C}^{(N-r)\times(N-r)}$ is assumed invertible, then

(5.62) $\displaystyle\bm{R}=\begin{bmatrix}\bm{I}_{n}&\bm{BD}^{-1}\\\
0&\bm{I}\end{bmatrix}\begin{bmatrix}\bm{A}-\bm{BD}^{-1}\bm{C}&0\\\
\bm{C}&\bm{D}\end{bmatrix}.$

If further the Schur Complement $\bm{A}-\bm{BD}^{-1}\bm{C}$ is invertible,
then

(5.63)
$\displaystyle\bm{R}^{-1}=\begin{bmatrix}(\bm{A}-\bm{BD}^{-1}\bm{C})^{-1}&-(\bm{A}-\bm{BD}^{-1}\bm{C})^{-1}\bm{BD}^{-1}\\\
&\\\
-\bm{D}^{-1}\bm{C}(\bm{A}-\bm{BD}^{-1}\bm{C})^{-1}&\bm{D}^{-1}+\bm{D}^{-1}\bm{C}(\bm{A}-\bm{BD}^{-1}\bm{C})^{-1}\bm{BD}^{-1}\end{bmatrix}$

as can be verified by block multiplication of $\bm{R}$ or by $\bm{R}$. We will
use $\bm{S}$ for the Schur Complement $\bm{S}=\bm{A}-\bm{BD}^{-1}\bm{C}$. We
will take $\bm{R}=z\bm{C}_{1}-\bm{C}_{0}$. We may already use this to
establish for each of the four classes of linearizations that

(5.64)
$\displaystyle\det\bm{R}=\det(z\bm{C}_{1}-\bm{C}_{0})=\det(\bm{A}-\bm{BD}^{-1}\bm{C})\det\bm{D}=\det\bm{P}(z)\>.$

Notice that the coefficients of $\bm{P}$ do not appear in the $\bm{D}$ block
(in any of our linearizations). Thus the Schur Complement carries all the
information particular to $\bm{P}(z)$. The computations verifying (5.64) are
not obvious but in each case $\bm{D}^{-1}$ plays an important role. We will
see that generically $\bm{D}^{-1}$ exists, except for isolated values of $z$,
which we can safely ignore and recover later by continuity.
We take each case in turn.

######  Theorem 5.1.

If
$\bm{C}_{1}=\mathrm{diag}\begin{bmatrix}\dfrac{1}{\alpha_{\ell-1}}\bm{P}_{\ell}&\bm{I}_{n}&\bm{I}_{n}&\cdots&\bm{I}_{n}\end{bmatrix}$
and

(5.65)
$\displaystyle\bm{C}_{0}=\begin{bmatrix}\dfrac{\beta_{\ell-1}}{\alpha_{\ell-1}}\bm{P}_{\ell}-\bm{P}_{\ell-1}&\dfrac{\gamma_{\ell-1}}{\alpha_{\ell-1}}\bm{P}_{\ell}-\bm{P}_{\ell-2}&-\bm{P}_{\ell-1}&\cdots&-\bm{P}_{0}\\\
\alpha_{\ell-2}\bm{I}_{\ell}&\beta_{\ell-2}\bm{I}_{n}&\gamma_{\ell-2}\bm{I}_{n}&&\\\
&\alpha_{\ell-3}\bm{I}_{n}&\beta_{\ell-3}\bm{I}_{n}&\gamma_{\ell-3}\bm{I}_{n}&\\\
&&\ddots&\ddots&\gamma_{1}\bm{I}_{n}\\\
&&&\alpha_{0}\bm{I}_{n}&\beta_{0}\bm{I}_{n}\end{bmatrix}$

and $\bm{X}=\begin{bmatrix}0&0&\cdots&0&\bm{I}_{n}\end{bmatrix}$ and
$\bm{Y}=\begin{bmatrix}\bm{I}_{n}&0&0&\cdots&0\end{bmatrix}$ then
$\bm{X}(z\bm{C}_{1}-\bm{C}_{0})^{-1}\bm{Y}=\bm{P}^{-1}(z)$ where
$\bm{P}(z)=\sum_{k=0}^{\ell}\bm{P}_{k}\phi_{k}(z)$ except for such $z$ that
$\det\bm{P}(z)=0$. As in section 3.2 the polynomials $\phi_{k}(z)$ satisfy
$z\phi_{k}=\alpha_{k}\phi_{k+1}+\beta_{k}\phi_{k}+\gamma_{k}\phi_{k-1}$,
$\phi_{-1}=0$, $\phi_{0}=1$,
$\phi_{1}=\nicefrac{{(z-\beta_{0})}}{{\alpha_{0}}}$. In this theorem,
$\ell\geq 2$ and $N=\ell n$, and if $\bm{P}_{\ell}\neq 0_{n}$ then degree
$\bm{P}=\ell$.

That this is a linearization is well-known; see e.g. [2]. We only prove
$\bm{P}^{-1}(z)=\bm{XR}^{-1}\bm{Y}$, here.

######  Proof 5.2.

We use the first block column of Schur Complement inverse formula

(5.66) $\displaystyle\bm{R}^{-1}=\begin{bmatrix}\bm{S}^{-1}&*\enspace\\\
-\bm{D}^{-1}\bm{CS}^{-1}&*\enspace\end{bmatrix}\>.$

Here

(5.67)
$\displaystyle\bm{D}=\begin{bmatrix}(z-\beta_{\ell-2})\bm{I}_{n}&-\gamma_{\ell-2}\bm{I}_{n}&&&\\\
-\alpha_{\ell-3}\bm{I}_{n}&(z-\beta_{\ell-3})\bm{I}_{n}&-\gamma_{\ell-3}\bm{I}_{n}&&\\\
&-\alpha_{\ell-4}\bm{I}_{n}&&\ddots&\\\
&&\ddots&\ddots&-\gamma_{1}\bm{I}_{n}\\\
&&&-\alpha_{0}\bm{I}_{n}&(z-\beta_{0})\bm{I}\end{bmatrix}$

is block tridiagonal, and

(5.68) $\displaystyle\bm{C}=\begin{bmatrix}-\alpha_{\ell-2}\bm{I}_{n}\\\ 0\\\
0\\\ \vdots\\\ 0\end{bmatrix}\>.$

By inspection $\bm{V}=-\bm{D}^{-1}\bm{C}$ is

(5.69) $\displaystyle\bm{V}=q\begin{bmatrix}\phi_{\ell-2}(z)\bm{I}_{n}\\\
\vdots\\\ \phi_{2}(z)\bm{I}_{n}\\\ \phi_{1}(z)\bm{I}_{n}\\\
\phi_{0}(z)\bm{I}_{n}\end{bmatrix}$

for some constant $q$, because

(5.70)
$\displaystyle-\alpha_{k}\phi_{k+1}(z)+(z-\beta_{k})\phi_{k}(z)-\gamma_{k}\phi_{k-1}(z)=0$

for $k=0,1,\cdots,\ell-3$. The constant $q$ is obtained from

(5.71) $\displaystyle
q\cdot(z-\beta_{\ell-2})\phi_{\ell-2}(z)-q\cdot\gamma_{\ell-2}\phi_{\ell-3}(z)=+\alpha_{\ell-2}$

or

(5.72) $\displaystyle q\cdot[\phi_{\ell-1}(z)]=+1$

So

(5.73) $\displaystyle q=\dfrac{+1}{\phi_{\ell-1}(z)}\>.$

It follows that

$\displaystyle\bm{S}=\dfrac{z-\beta_{\ell-1}}{\alpha_{\ell-1}}\bm{P}_{\ell}+\bm{P}_{\ell-1}+\begin{bmatrix}\dfrac{-\gamma_{\ell-1}}{\alpha_{\ell-1}}\bm{P}_{\ell}+\bm{P}_{\ell-2}&\bm{P}_{\ell-3}&\cdots&\bm{P}_{0}\end{bmatrix}\begin{bmatrix}\phi_{\ell-2}(z)\\\
\phi_{\ell-3}(z)\\\ \vdots\\\ \phi_{0}(z)\\\
\end{bmatrix}\cdot\dfrac{1}{\phi_{\ell-1}(z)}$
$\displaystyle=\dfrac{\dfrac{z-\beta_{\ell-1}}{\alpha_{\ell-1}}\phi_{\ell-1}(z)\bm{P}_{\ell}+\phi_{\ell-1}(z)\bm{P}_{\ell-1}-\dfrac{\gamma_{\ell-1}}{\alpha_{\ell-1}}\phi_{\ell-2}(z)\bm{P}_{\ell}+\phi_{\ell-2}(z)\bm{P}_{\ell-2}+\cdots+\phi_{0}(z)\bm{P}_{0}}{\phi_{\ell-1}(z)}$
(5.74)
$\displaystyle=\dfrac{\sum_{k=0}^{\ell}\phi_{k}(z)\bm{P}_{k}}{\phi_{\ell-1}(z)}=\dfrac{\bm{P}(z)}{\phi_{\ell-1}(z)}\>.$

Thus

(5.75)
$\displaystyle-\bm{D}^{-1}\bm{CS}^{-1}=\begin{bmatrix}\phi_{\ell-2}(z)\bm{I}_{n}\\\
\vdots\\\ \phi_{0}(z)\bm{I}_{n}\end{bmatrix}\bm{P}^{-1}(z)$

because $\dfrac{1}{\phi_{\ell-1}(z)}\bm{S}^{-1}=\bm{P}^{-1}(z)$. Finally,
$\phi_{0}(z)=1$, so the bottom block is $\bm{P}^{-1}(z)$, establishing that

(5.76)
$\displaystyle\bm{X}=\begin{bmatrix}0&0&\cdots&0&\bm{I}_{n}\end{bmatrix}$
(5.77)
$\displaystyle\bm{Y}=\begin{bmatrix}\bm{I}_{n}&0&\cdots&0&0\end{bmatrix}^{T}$

will produce $\bm{XR}^{-1}\bm{Y}=\bm{P}^{-1}(z)$.

######  Theorem 5.3.

Put

(5.78)
$\displaystyle\bm{C}_{1}=\begin{bmatrix}\dfrac{1}{\ell}\bm{P}_{\ell}-\bm{P}_{\ell-1}&-\bm{P}_{\ell-2}&\cdots&-\bm{P}_{1}&-\bm{P}_{0}\\\
\bm{I}_{n}&\dfrac{2}{\ell-1}\bm{I}_{n}&&&\\\
&\bm{I}_{n}&\dfrac{3}{\ell-2}\bm{I}_{n}&&\\\ &&\ddots&\ddots&\\\
&&&\bm{I}_{n}&\dfrac{\ell}{1}\bm{I}_{n}\end{bmatrix}$

and

(5.79)
$\displaystyle\bm{C}_{0}=\begin{bmatrix}-\bm{P}_{\ell-1}&-\bm{P}_{\ell-2}&\cdots&-\bm{P}_{1}&-\bm{P}_{0}\\\
\bm{I}_{n}&0&&&\\\ &\bm{I}_{n}&0&&\\\ &&\ddots&\ddots&\\\
&&&\bm{I}_{n}&0\end{bmatrix}$

and $\bm{Y}=\begin{bmatrix}\bm{I}_{n}&0&\cdots&0&0\end{bmatrix}^{T}$ with
$\bm{X}=\begin{bmatrix}\dfrac{1}{\ell}\bm{I}_{n}&\dfrac{2}{\ell}\bm{I}_{n}&\dfrac{3}{\ell}\bm{I}_{n}&\cdots&\dfrac{\ell}{\ell}\bm{I}_{n}\end{bmatrix}$.
Then $\bm{X}(z\bm{C}_{1}-\bm{C}_{0})^{-1}\bm{Y}=\bm{P}^{-1}(z)$, unless
$z\in\Lambda(\bm{P})$, and
$\det\bm{P}(z)=\det\bm{R}(z)=\det(z\bm{C}_{1}-\bm{C}_{0})$.

######  Proof 5.4.

This linearization in proved e.g. in [27], but for convenience we supply one
here as well. The Schur factoring is

(5.80) $\displaystyle\bm{R}=\begin{bmatrix}\bm{I}_{n}&\bm{BD}^{-1}\\\
0&\bm{I}_{N-r}\end{bmatrix}\begin{bmatrix}\bm{S}&0\\\
\bm{C}&\bm{D}\end{bmatrix}$

where $\bm{S}=\bm{A}-\bm{BD}^{-1}\bm{C}$ is the Schur Complement. Here

(5.81) $\displaystyle A=\dfrac{z}{\ell}\bm{P}_{\ell}+(1-z)\bm{P}_{\ell-1}$
(5.82)
$\displaystyle\bm{B}=\begin{bmatrix}(1-z)\bm{P}_{\ell-2}&(1-z)\bm{P}_{\ell-3}&\cdots&(1-z)\bm{P}_{0}\end{bmatrix}$
(5.83) $\displaystyle\bm{C}=\begin{bmatrix}(z-1)\bm{I}_{n}\\\ 0\\\ 0\\\
\vdots\\\ 0\end{bmatrix}$

and

(5.84) $\displaystyle\bm{D}=\begin{bmatrix}\dfrac{2}{\ell-1}z\bm{I}_{n}&&&\\\
(z-1)\bm{I}_{n}&\dfrac{3}{\ell-2}z\bm{I}_{n}&&\\\ &(z-1)\bm{I}_{n}&\ddots&\\\
&&\ddots&\\\ &&(z-1)\bm{I}_{n}&\dfrac{\ell}{1}z\bm{I}_{n}\end{bmatrix}$

Therefore $\bm{V}=\bm{D}^{-1}\bm{C}$ satisfies

(5.85) $\displaystyle\begin{bmatrix}\dfrac{2}{\ell-1}z\bm{I}_{n}&&&\\\
(z-1)\bm{I}_{n}&\dfrac{3}{\ell-2}z\bm{I}_{n}&&\\\
&(z-1)\bm{I}_{n}&\dfrac{4}{\ell-3}z\bm{I}_{n}&\\\ &&\ddots&\\\ &&\ddots&\\\
&&(z-1)\bm{I}_{n}&\dfrac{\ell}{1}z\bm{I}_{n}\end{bmatrix}\begin{bmatrix}\bm{v}_{1}\\\
\bm{v}_{2}\\\ \vdots\\\
\bm{v}_{\ell-1}\end{bmatrix}=\begin{bmatrix}(z-1)\bm{I}_{n}\\\ 0\\\ \vdots\\\
0\end{bmatrix}$

So

(5.86)
$\displaystyle\bm{v}_{1}=\dfrac{\ell-1}{2}\left(\dfrac{z-1}{z}\right)\bm{I}_{n}=-\dfrac{\ell-1}{2}\left(\dfrac{1-z}{z}\right)\bm{I}_{n}$
(5.87)
$\displaystyle\bm{v}_{2}=-\dfrac{\ell-2}{3}\cdot\bm{v}_{1}=-\dfrac{\ell-2}{3}\cdot\dfrac{\ell-1}{2}\cdot\left(\dfrac{1-z}{z}\right)^{2}\bm{I}_{n}$
(5.88)
$\displaystyle\bm{v}_{3}=-\dfrac{\ell-3}{4}\cdot\dfrac{\ell-2}{3}\cdot\dfrac{\ell-1}{2}\left(\dfrac{1-z}{z}\right)^{3}\bm{I}_{n}$

and so on; by inspection, confirmed by a formal induction not given here,

(5.89)
$\displaystyle\bm{v}_{k}=-\dfrac{(\ell-1)!}{(\ell-k-1)!(k+1)!}\left(\dfrac{1-z}{z}\right)^{k}\bm{I}_{n}=-\dfrac{1}{\ell}{\ell\choose
k+1}\left(\dfrac{1-z}{z}\right)^{k}\bm{I}_{n}$

for $k=1,\cdots,\ell-1$. Thus

$\displaystyle\bm{S}$
$\displaystyle=\dfrac{z}{\ell}\bm{P}_{\ell}+(1-z)\bm{P}_{\ell-1}+(1-z)\begin{bmatrix}\bm{P}_{\ell-2}&\bm{P}_{\ell-3}&\cdots&\bm{P}_{0}\end{bmatrix}\begin{bmatrix}\dfrac{1}{\ell}{\ell\choose
2}\left(\dfrac{1-z}{z}\right)\bm{I}_{n}\\\ \dfrac{1}{\ell}{\ell\choose
3}\left(\dfrac{1-z}{z}\right)^{2}\bm{I}_{n}\\\ \vdots\\\
\dfrac{1}{\ell}{\ell\choose\ell}\left(\dfrac{1-z}{z}\right)^{\ell-1}\bm{I}_{n}\\\
\end{bmatrix}$ $\displaystyle=\dfrac{1}{\ell
z^{\ell-1}}\cdot\left[z^{\ell}\bm{P}_{\ell}+\ell
z^{\ell-1}(1-z)\bm{P}_{\ell-1}+{\ell\choose
2}z^{\ell-2}(1-z)^{2}\bm{P}_{\ell-2}+\cdots+{\ell\choose\ell}(1-z)^{\ell}\bm{P}_{0}\right]$
(5.90) $\displaystyle=\dfrac{\bm{P}(z)}{\ell z^{\ell-1}}\>.$

Hence

$\displaystyle\det\bm{R}$ $\displaystyle=\det\bm{S}\det\bm{D}$
$\displaystyle=\dfrac{\det\bm{P}(z)}{(\ell
z^{\ell-1})^{n}}\cdot\left(\dfrac{2}{\ell-1}\cdot\dfrac{3}{\ell-2}\cdots\dfrac{\ell-1}{2}\cdot\ell\cdot
z\right)^{n}$ (5.91) $\displaystyle=\det\bm{P}(z)\>.$

This establishes the linearization. Moreover,

(5.92) $\displaystyle\bm{S}^{-1}=\ell z^{\ell-1}\bm{P}^{-1}(z)$

and the first column of $\bm{R}^{-1}$ is

(5.93) $\displaystyle\begin{bmatrix}\bm{S}^{-1}\\\
-\bm{D}^{-1}\bm{CS}^{-1}\end{bmatrix}=\begin{bmatrix}\ell
z^{\ell-1}\bm{P}^{-1}\\\ \ell
z^{\ell-1}\cdot\dfrac{1}{\ell}\leavevmode\resizebox{}{14.0pt}{$\displaystyle\binom{\ell}{2}$}\left(\dfrac{1-z}{z}\right)\bm{P}^{-1}\\\
\ell
z^{\ell-1}\cdot\dfrac{1}{\ell}\leavevmode\resizebox{}{14.0pt}{$\displaystyle\binom{\ell}{3}$}\left(\dfrac{1-z}{z}\right)^{2}\bm{P}^{-1}\\\
\ell
z^{\ell-1}\cdot\dfrac{1}{\ell}\leavevmode\resizebox{}{14.0pt}{$\displaystyle\binom{\ell}{4}$}\left(\dfrac{1-z}{z}\right)^{3}\bm{P}^{-1}\\\
\vdots\\\ \ell
z^{\ell-1}\cdot\dfrac{1}{\ell}\leavevmode\resizebox{}{14.0pt}{$\displaystyle\binom{\ell}{\ell}$}\left(\dfrac{1-z}{z}\right)^{\ell-1}\bm{P}^{-1}\\\
\end{bmatrix}=\begin{bmatrix}nz^{n-1}\bm{P}^{-1}\\\ {\ell\choose
2}z^{\ell-2}(1-z)\bm{P}^{-1}\\\ {\ell\choose
3}z^{\ell-3}(1-z)^{2}\bm{P}^{-1}\\\ \vdots\\\
{\ell\choose\ell}z^{0}(1-z)^{\ell-1}\bm{P}^{-1}\end{bmatrix}$

We now notice that $1$, expressed as a linear combination of

(5.94) ${\ell\choose 1}z^{\ell-1},{\ell\choose
2}z^{\ell-2}(1-z),\cdots,{\ell\choose\ell}z^{0}(1-z)^{\ell-1}$

is

$\displaystyle 1$ $\displaystyle=\dfrac{1}{\ell}\cdot{\ell\choose
1}z^{\ell-1}+\dfrac{2}{\ell}\cdot{\ell\choose
2}z^{\ell-2}(1-z)+\cdots+\dfrac{\ell}{\ell}\cdot{\ell\choose\ell}z^{0}(1-z)^{\ell-1}$
$\displaystyle={\ell-1\choose 0}z^{\ell-1}(1-z)^{0}+{\ell-1\choose
1}z^{\ell-2}(1-z)^{1}+\cdots+{\ell-1\choose\ell-1}z^{0}(1-z)^{\ell-1}$ (5.95)
$\displaystyle=(z+1-z)^{\ell-1}\>.$

Indeed we use a degree-reduced Bernstein bases here, ${\ell-1\choose
k}z^{k}(1-z)^{\ell-1-k}$, to express $1$.
In any case, the coefficients of $1$ give us our $\bm{X}$ vector:
$\bm{XR}^{-1}\bm{Y}=\bm{P}^{-1}(z)$.

######  Theorem 5.5 (Lagrange Basis).

If $\bm{P}(z)\in\mathbb{C}^{n\times n}$ is of degree at most $\ell$, and takes
the values $\bm{\rho}_{k}\in\mathbb{C}^{n\times n}$ at the $\ell+1$ distinct
nodes $z=\tau_{k}$, $0\leq k\leq\ell$, i.e
$\bm{P}(\tau_{k})=\bm{\rho}_{k}\in\mathbb{C}^{n\times n}$, and the reciprocal
of the node polynomial $\omega(z)=\prod_{k=0}^{\ell}(z-\tau_{k})$ has partial
fraction expansion

(5.96)
$\displaystyle\dfrac{1}{\omega(z)}=\sum_{k=0}^{\ell}\dfrac{\beta_{k}}{z-\tau_{k}}$

then a linearization for $\bm{P}(z)$ is $z\bm{C}_{1}-\bm{C}_{0}$ where
$\bm{C}_{1}=\mathrm{diag}(\bm{0}_{n},\bm{I}_{n},\bm{I}_{n},\cdots,\bm{I}_{n})$
with $\ell+2$ diagonal blocks, so $N=(\ell+2)r$, and

(5.97)
$\displaystyle\bm{C}_{0}=\begin{bmatrix}0&-\bm{\rho}_{0}&-\bm{\rho}_{1}&-\bm{\rho}_{2}&\cdots&-\bm{\rho}_{\ell}\\\
\beta_{0}\bm{I}_{n}&\tau_{0}\bm{I}_{n}&&&&\\\
\beta_{1}\bm{I}_{n}&&\tau_{1}\bm{I}_{n}&&&\\\
\beta_{2}\bm{I}_{n}&&&\tau_{2}\bm{I}_{n}&&\\\ \vdots&&&&\ddots&\\\
\beta_{\ell}\bm{I}_{n}&&&&&\tau_{\ell}\bm{I}\end{bmatrix}\>.$

Moreover, if $\bm{Y}=\begin{bmatrix}\bm{I}_{n}&0&0&\cdots&0\end{bmatrix}^{T}$
and
$\bm{X}=\begin{bmatrix}0_{n}&\bm{I}_{n}&\bm{I}_{n}&\cdots&\bm{I}_{n}\end{bmatrix}$
then $\bm{X}(z\bm{C}_{1}-\bm{C}_{0})^{-1}\bm{Y}=\bm{P}^{-1}(z)$ where
$z\in\Lambda(\bm{P})$.

######  Proof 5.6.

Again we use the Schur complement: $\bm{S}=\bm{A}-\bm{BD}^{-1}\bm{C}$ where
here

(5.98) $\displaystyle\bm{A}$ $\displaystyle=\bm{0}_{n}$ (5.99)
$\displaystyle\bm{B}$
$\displaystyle=-\begin{bmatrix}\bm{\rho}_{0}&\bm{\rho}_{1}&\cdots&\bm{\rho}_{\ell}\end{bmatrix}$
(5.100) $\displaystyle\bm{D}^{-1}$
$\displaystyle=\mathrm{diag}\left(\dfrac{1}{z-\tau_{0}}\bm{I}_{n},\dfrac{1}{z-\tau_{1}}\bm{I}_{n},\cdots,\dfrac{1}{z-\tau_{\ell}}\bm{I}_{n}\right)$
(5.101) $\displaystyle\bm{C}$
$\displaystyle=\begin{bmatrix}\beta_{0}\bm{I}_{n}\\\ \beta_{1}\bm{I}_{n}\\\
\vdots\\\ \beta_{\ell}\bm{I}_{n}\end{bmatrix}$

So

(5.102)
$\displaystyle\bm{S}=\sum_{k=0}^{\ell}\dfrac{\beta_{k}}{z-\tau_{k}}\bm{\rho}_{k}=\omega(z)^{-1}\bm{P}(z)$

from the first barycentric formula [3].
Note the first column of $\bm{R}^{-1}(z)$ is $\begin{bmatrix}\bm{S}^{-1}\\\
-\bm{CD}^{-1}\bm{S}^{-1}\end{bmatrix}$ or

(5.103) $\displaystyle\begin{bmatrix}\omega(z)\bm{P}^{-1}(z)\\\
\left(\dfrac{\beta_{0}}{z-\tau_{0}}\right)\omega(z)\bm{P}^{-1}(z)\\\
\left(\dfrac{\beta_{1}}{z-\tau_{1}}\right)\omega(z)\bm{P}^{-1}(z)\\\ \vdots\\\
\left(\dfrac{\beta_{\ell}}{z-\tau_{\ell}}\right)\omega(z)\bm{P}^{-1}(z)\\\
\end{bmatrix}$

Note that
$\sum_{k=0}^{\ell}\dfrac{\beta_{k}}{z-\tau_{k}}=\dfrac{1}{\omega(z)}$, so

(5.104)
$\displaystyle\begin{bmatrix}0&\bm{I}_{n}&\bm{I}_{n}&\cdots&\bm{I}_{n}\end{bmatrix}\cdot\bm{R}^{-1}\begin{bmatrix}\bm{I}_{n}\\\
0\\\ \vdots\\\
0\end{bmatrix}=\left(\sum_{k=0}^{\ell}\dfrac{\beta_{k}}{z-\tau_{k}}\right)\omega(z)\bm{P}^{-1}(z)=\bm{P}^{-1}(z)$

######  Theorem 5.7.

In the Hermite interpolational bases on $m+1$ nodes each with coefficiency
$s_{i}$, so the degree $\ell$ is at most $-1+\sum_{k=0}^{m}s_{k}$, the
barycentric weights are

(5.105)
$\displaystyle\dfrac{1}{\omega(z)}=\sum_{i=0}^{m}\sum_{j=0}^{s_{i}-1}\dfrac{\beta_{ij}}{(z-\tau_{i})^{j+1}}$

As in the Lagrange case,
$\bm{C}_{1}=\mathrm{diag}(0,\bm{I}_{n},\cdots,\bm{I}_{n})$. $\bm{C}_{0}$ is as
below:

(5.106)
$\displaystyle\bm{C}_{0}=\begin{bmatrix}0&-\hat{\bm{\rho}}_{0}&-\hat{\bm{\rho}}_{1}&\cdots&-\hat{\bm{\rho}}_{m}\\\
\beta_{0,s_{0}-1}\bm{I}_{n}&\bm{J}^{T}_{0}&&&\\\
\beta_{0,s_{0}-2}\bm{I}_{n}&&\bm{J}^{T}_{1}&&\\\ \vdots&&&\ddots&\\\
\beta_{m,s_{m}-1}&&&&\bm{J}^{T}_{m}\end{bmatrix}$

where each block per node of data is collected in the $n\times m\ell$ block
matrix

(5.107)
$\displaystyle\hat{\bm{\rho}}_{i}=\begin{bmatrix}\bm{\rho}_{i,s_{i}-1}&\bm{\rho}_{i,s_{i}-2}&\cdots&\bm{\rho}_{i,0}\end{bmatrix}\>.$

Each diagonal node block is a tensor product of a transposed Jordan block:

(5.108) $\displaystyle\bm{J}_{i}=\begin{bmatrix}\tau_{i}\bm{I}_{n}&&&&\\\
\bm{I}_{n}&\tau_{i}\bm{I}_{n}&&&\\\ &\bm{I}_{n}&\tau_{i}\bm{I}_{n}&&\\\
&&\ddots&\ddots&\\\ &&&\bm{I}_{n}&\tau_{i}\bm{I}_{n}\end{bmatrix}\>.$

This form arises naturally on letting distinct Lagrange nodes flow together in
a limit.

Express $1$ as a polynomial in this basis. Then
$1\longleftrightarrow\rho_{00}=1,\rho_{10}=1,\cdots,\rho_{n0}=1$ and all other
components are zero. Put

(5.109)
$\bm{X}=\begin{bmatrix}0&\smash[b]{\underbrace{\begin{matrix}0&\cdots&0&1\end{matrix}}_{\text{$S_{0}$
entries}}}&\smash[b]{\underbrace{\begin{matrix}0&\cdots&0&1\end{matrix}}_{\text{$S_{1}$
entries}}}&\cdots&1\end{bmatrix}\otimes\bm{I}_{n}$

and $\bm{Y}=\begin{bmatrix}\bm{I}_{n}&0&0&\cdots&0\end{bmatrix}$.
A similar but more involved computation than in theorem 5.5 gives

(5.110)
$\bm{S}=\dfrac{1}{\omega(z)}\bm{P}(z)=\sum_{i=0}^{m}\sum_{j=0}^{s_{i}-1}\sum_{k=0}^{j}\beta_{ij}\bm{\rho}_{ik}(z-\tau_{i})^{k-j-1}$

and $\bm{D}^{-1}\bm{C}$ contains just the correct powers of $(z-\tau_{i})$
divided into $\beta_{ij}$ to make the sums come out right; the inverse of the
block

(5.111) $\begin{bmatrix}(z-\tau_{0})\bm{I}_{n}&&&\\\
-\bm{I}_{n}&(z-\tau_{0})\bm{I}_{n}&&&\\\ &-\bm{I}_{n}&\ddots&&\\\
&&\ddots&&\\\ &&&-\bm{I}_{n}&(z-\tau_{0})\bm{I}_{n}\end{bmatrix}$

is

(5.112) $\begin{bmatrix}\dfrac{1}{z-\tau_{0}}\bm{I}_{n}&&&\\\
\dfrac{1}{(z-\tau_{0})^{2}}\bm{I}_{n}&\dfrac{1}{z-\tau_{0}}\bm{I}_{n}&&\\\
\dfrac{1}{(z-\tau_{0})^{3}}\bm{I}_{n}&\dfrac{1}{(z-\tau_{0})^{2}}\bm{I}_{n}&\dfrac{1}{z-\tau_{0}}\bm{I}_{n}&\\\
\vdots&&\ddots&\\\
\dfrac{1}{(z-\tau_{0})^{s_{0}}}\bm{I}_{n}&&&\dfrac{1}{z-\tau_{0}}\bm{I}_{n}\end{bmatrix}\>.$

and thus each block is reminiscent of theorem 5.1, in fact.

######  Remark 5.8.

In every case $\bm{X}=[\text{coefficients of}\;1]\otimes\bm{I}$,
$\bm{Y}=[1,0,\cdots,0]\otimes\bm{I}$. This is in agreement with our universal
proof in section 2.

##  6 Examples

In this section, we will show some experiments done in Maple 2017 to
demonstrate that standard triples introduced in section works for the
different bases. We wrote our own code for constructing the linearizations
rather than using Maple’s built-in CompanionMatrix function since the result
of the built-in function is the flipped and transposed version of the
companion matrices compared to the structure in this paper.

For the following examples, we check the correctness of the standard triple
for each of the following examples by rearranging the resolvent form

$\displaystyle\bm{P}^{-1}(z)$
$\displaystyle=\bm{X}\left(z\bm{C}_{1}-\bm{C}_{0}\right)^{-1}\bm{Y}$
$\displaystyle\bm{I}_{n}$
$\displaystyle=\bm{X}\left(z\bm{C}_{1}-\bm{C}_{0}\right)^{-1}\bm{Y}\bm{P}(z)\>.$

Since these computations are done exactly, the result will exactly equal the
identity matrix. For the Lagrange basis example, since we construct our
companion matrices using $\tau$ and $\rho$ instead of the matrix polynomial
itself, $\bm{P}(z)$ is constructed using the barycentric Lagrange
interpolation formula, which can be derived from equation (5.102). The Hermite
interpolational basis examples are handled similarly to the Lagrange case,
where $\bm{P}(z)$ is the Hermite interpolation polynomial, which can be
derived from equation (5.110).

###  6.1 Bases with three-term recurrence relations

######  Example 6.1 (Chebyshev basis of the first kind).

(6.113) $\displaystyle\bm{P}(z)$
$\displaystyle=\begin{bmatrix}\nicefrac{{1}}{{5}}&\nicefrac{{7}}{{100}}\\\
-\nicefrac{{93}}{{200}}&-\nicefrac{{29}}{{200}}\end{bmatrix}T_{0}(z)+\begin{bmatrix}\nicefrac{{53}}{{300}}&\nicefrac{{7}}{{60}}\\\
\nicefrac{{2}}{{25}}&\nicefrac{{3}}{{50}}\end{bmatrix}T_{1}(z)$
$\displaystyle+\begin{bmatrix}-\nicefrac{{9}}{{80}}&-\nicefrac{{13}}{{80}}\\\
\nicefrac{{57}}{{400}}&-\nicefrac{{47}}{{400}}\end{bmatrix}T_{2}(z)+\begin{bmatrix}-\nicefrac{{3}}{{250}}&-\nicefrac{{31}}{{500}}\\\
-\nicefrac{{77}}{{500}}&\nicefrac{{27}}{{250}}\end{bmatrix}T_{3}(z)\>.$

The standard triple for equation (6.113) is

$\displaystyle\bm{C}_{0}=\begin{bmatrix}-{\nicefrac{{21}}{{200}}}&-{\nicefrac{{29}}{{400}}}&-{\nicefrac{{107}}{{750}}}&-{\nicefrac{{17}}{{50}}}&-{\nicefrac{{73}}{{200}}}&\nicefrac{{1}}{{25}}\\\
-{\nicefrac{{13}}{{400}}}&-{\nicefrac{{1}}{{400}}}&-{\nicefrac{{49}}{{300}}}&-{\nicefrac{{283}}{{750}}}&-{\nicefrac{{27}}{{200}}}&{\nicefrac{{9}}{{25}}}\\\
\nicefrac{{1}}{{2}}&0&0&0&\nicefrac{{1}}{{2}}&0\\\
0&\nicefrac{{1}}{{2}}&0&0&0&\nicefrac{{1}}{{2}}\\\ 0&0&1&0&0&0\\\
0&0&0&1&0&0\end{bmatrix}\quad\bm{C}_{1}=\begin{bmatrix}-{\nicefrac{{33}}{{250}}}&-{\nicefrac{{9}}{{25}}}&0&0&0&0\\\
{\nicefrac{{3}}{{25}}}&-{\nicefrac{{37}}{{250}}}&0&0&0&0\\\ 0&0&1&0&0&0\\\
0&0&0&1&0&0\\\ 0&0&0&0&1&0\\\ 0&0&0&0&0&1\end{bmatrix}$
$\displaystyle\bm{X}=\begin{bmatrix}0&0&0&0&1&0\\\
0&0&0&0&0&1\end{bmatrix}\quad\bm{Y}=\begin{bmatrix}1&0\\\ 0&1\\\ 0&0\\\ 0&0\\\
0&0\\\ 0&0\end{bmatrix}\>.$

Then,

$\bm{X}\left(z\bm{C}_{1}-\bm{C}_{0}\right)^{-1}\bm{Y}\bm{P}(z)=\bm{I}_{2}\>.$

which indicates that the standard triple is correct.

######  Example 6.2 (Newton Interpolational Basis).

(6.114) $\tau=\left[\mathrm{seq}\left(\cos\left(\frac{\pi\cdot
k}{3}\right),k=0..3\right)\right]=\left[1,\nicefrac{{1}}{{2}},-\nicefrac{{1}}{{2}},-1\right]$
(6.115) $\displaystyle\bm{P}(z)$ $\displaystyle=\begin{bmatrix}6&25\\\
-1&5\end{bmatrix}\prod_{j=0}^{0}(z-\tau_{j})+\begin{bmatrix}-{\nicefrac{{80}}{{3}}}&{\nicefrac{{25}}{{3}}}\\\
{\nicefrac{{43}}{{3}}}&{\nicefrac{{94}}{{3}}}\end{bmatrix}\prod_{j=0}^{1}(z-\tau_{j})$
$\displaystyle+\begin{bmatrix}{\nicefrac{{77}}{{4}}}&{\nicefrac{{31}}{{4}}}\\\
\nicefrac{{9}}{{4}}&-{\nicefrac{{25}}{{2}}}\end{bmatrix}\prod_{j=0}^{2}(z-\tau_{j})+\begin{bmatrix}{\nicefrac{{86}}{{5}}}&-{\nicefrac{{61}}{{5}}}\\\
4&-{\nicefrac{{48}}{{5}}}\end{bmatrix}\prod_{j=0}^{3}(z-\tau_{j})$

The standard triple for equation (6.115) is

$\displaystyle\bm{C}_{0}=\begin{bmatrix}-{\nicefrac{{557}}{{20}}}&-{\nicefrac{{33}}{{20}}}&{\nicefrac{{80}}{{3}}}&-{\nicefrac{{25}}{{3}}}&-6&-25\\\
-{\nicefrac{{17}}{{4}}}&{\nicefrac{{173}}{{10}}}&-{\nicefrac{{43}}{{3}}}&-{\nicefrac{{94}}{{3}}}&1&-5\\\
1&0&1/2&0&0&0\\\ 0&1&0&1/2&0&0\\\ 0&0&1&0&1&0\\\
0&0&0&1&0&1\end{bmatrix}\quad\bm{C}_{1}=\begin{bmatrix}{\frac{86}{5}}&-{\frac{61}{5}}&0&0&0&0\\\
4&-{\frac{48}{5}}&0&0&0&0\\\ 0&0&1&0&0&0\\\ 0&0&0&1&0&0\\\ 0&0&0&0&1&0\\\
0&0&0&0&0&1\end{bmatrix}$ $\displaystyle\bm{X}=\begin{bmatrix}0&0&0&0&1&0\\\
0&0&0&0&0&1\end{bmatrix}\quad\bm{Y}=\begin{bmatrix}1&0\\\ 0&1\\\ 0&0\\\ 0&0\\\
0&0\\\ 0&0\end{bmatrix}\>.$

###  6.2 Bernstein Basis

######  Example 6.3 (Non-singular leading coefficient case).

(6.116) $\displaystyle\bm{P}(z)$
$\displaystyle=\begin{bmatrix}\nicefrac{{4}}{{25}}&\nicefrac{{99}}{{100}}\\\
\nicefrac{{9}}{{100}}&\nicefrac{{3}}{{5}}\end{bmatrix}B^{3}_{0}(z)+\begin{bmatrix}-\nicefrac{{17}}{{25}}&\nicefrac{{11}}{{50}}\\\
-\nicefrac{{67}}{{100}}&\nicefrac{{7}}{{50}}\end{bmatrix}B^{3}_{1}(z)$
$\displaystyle+\begin{bmatrix}-\nicefrac{{59}}{{100}}&-\nicefrac{{31}}{{50}}\\\
\nicefrac{{3}}{{25}}&-\nicefrac{{33}}{{100}}\end{bmatrix}B^{3}_{2}(z)+\begin{bmatrix}\nicefrac{{41}}{{50}}&\nicefrac{{21}}{{50}}\\\
\nicefrac{{18}}{{25}}&\nicefrac{{9}}{{50}}\end{bmatrix}B^{3}_{3}(z)\>.$

The standard triple for equation (6.116) is

$\displaystyle\bm{C}_{0}=\begin{bmatrix}{\nicefrac{{59}}{{100}}}&{\nicefrac{{31}}{{50}}}&{\nicefrac{{17}}{{25}}}&-{\nicefrac{{11}}{{50}}}&-{\nicefrac{{4}}{{25}}}&-{\nicefrac{{99}}{{100}}}\\\
-{\nicefrac{{3}}{{25}}}&{\nicefrac{{33}}{{100}}}&{\nicefrac{{67}}{{100}}}&-{\nicefrac{{7}}{{50}}}&-{\nicefrac{{9}}{{100}}}&-\nicefrac{{3}}{{5}}\\\
1&0&0&0&0&0\\\ 0&1&0&0&0&0\\\ 0&0&1&0&0&0\\\
0&0&0&1&0&0\end{bmatrix}\quad\bm{C}_{1}=\begin{bmatrix}{\nicefrac{{259}}{{300}}}&{\nicefrac{{19}}{{25}}}&{\nicefrac{{17}}{{25}}}&-{\nicefrac{{11}}{{50}}}&-{\nicefrac{{4}}{{25}}}&-{\nicefrac{{99}}{{100}}}\\\
{\nicefrac{{3}}{{25}}}&{\nicefrac{{39}}{{100}}}&{\nicefrac{{67}}{{100}}}&-{\nicefrac{{7}}{{50}}}&-{\nicefrac{{9}}{{100}}}&-\nicefrac{{3}}{{5}}\\\
1&0&1&0&0&0\\\ 0&1&0&1&0&0\\\ 0&0&1&0&3&0\\\ 0&0&0&1&0&3\end{bmatrix}$
$\displaystyle\bm{X}=\begin{bmatrix}\nicefrac{{1}}{{3}}&0&\nicefrac{{2}}{{3}}&0&1&0\\\
0&\nicefrac{{1}}{{3}}&0&\nicefrac{{2}}{{3}}&0&1\end{bmatrix}\quad\bm{Y}=\begin{bmatrix}1&0\\\
0&1\\\ 0&0\\\ 0&0\\\ 0&0\\\ 0&0\end{bmatrix}\>.$

Then,

$\bm{X}\left(z\bm{C}_{1}-\bm{C}_{0}\right)^{-1}\bm{Y}\bm{P}(z)=\bm{I}_{2}\>.$

######  Example 6.4 (Singular leading coefficient case).

(6.117) $\displaystyle\bm{P}(z)$
$\displaystyle=\begin{bmatrix}\nicefrac{{29}}{{100}}&-\nicefrac{{8}}{{25}}\\\
\nicefrac{{7}}{{10}}&-\nicefrac{{1}}{{100}}\end{bmatrix}B^{3}_{0}(z)+\begin{bmatrix}-\nicefrac{{41}}{{50}}&\nicefrac{{41}}{{100}}\\\
-\nicefrac{{7}}{{10}}&\nicefrac{{91}}{{100}}\end{bmatrix}B^{3}_{1}(z)$
$\displaystyle+\begin{bmatrix}\nicefrac{{9}}{{10}}&\nicefrac{{19}}{{100}}\\\
\nicefrac{{4}}{{5}}&\nicefrac{{22}}{{25}}\end{bmatrix}B^{3}_{2}(z)+\begin{bmatrix}1&1\\\
\nicefrac{{9851}}{{1980}}&0\end{bmatrix}B^{3}_{3}(z)\>.$

Expressing equation (6.117) into the monomial basis, we have

$\bm{P}(z)=\begin{bmatrix}{\nicefrac{{29}}{{100}}}&-{\nicefrac{{8}}{{25}}}\\\
{\nicefrac{{7}}{{10}}}&-{\nicefrac{{1}}{{100}}}\end{bmatrix}+\begin{bmatrix}{\nicefrac{{29}}{{100}}}&-{\nicefrac{{8}}{{25}}}\\\
{\nicefrac{{7}}{{10}}}&-{\nicefrac{{1}}{{100}}}\end{bmatrix}z+\begin{bmatrix}{\nicefrac{{849}}{{100}}}&-{\nicefrac{{57}}{{20}}}\\\
{\nicefrac{{87}}{{10}}}&-{\nicefrac{{57}}{{20}}}\end{bmatrix}z^{2}+\begin{bmatrix}-{\nicefrac{{89}}{{20}}}&{\nicefrac{{99}}{{50}}}\\\
-{\nicefrac{{89}}{{396}}}&\nicefrac{{1}}{{10}}\end{bmatrix}z^{3}\>.$

Taking the determinant of the leading coefficient

$\det\left(\begin{bmatrix}-{\nicefrac{{89}}{{20}}}&{\nicefrac{{99}}{{50}}}\\\
-{\nicefrac{{89}}{{396}}}&\nicefrac{{1}}{{10}}\end{bmatrix}\right)=\left(-{\nicefrac{{89}}{{20}}}\right)\left(\nicefrac{{1}}{{10}}\right)-\left({\nicefrac{{99}}{{50}}}\right)\left(-{\nicefrac{{89}}{{396}}}\right)=0\>,$

we can observe that leading coefficient is singular, and thus, this matrix
polynomial is non-monic. The standard triple for equation (6.117) is

$\displaystyle\bm{C}_{0}=\begin{bmatrix}-{\nicefrac{{9}}{{10}}}&-{\nicefrac{{19}}{{100}}}&{\nicefrac{{41}}{{50}}}&-{\nicefrac{{41}}{{100}}}&-{\nicefrac{{29}}{{100}}}&{\nicefrac{{8}}{{25}}}\\\
-\nicefrac{{4}}{{5}}&-{\nicefrac{{22}}{{25}}}&{\nicefrac{{7}}{{10}}}&-{\nicefrac{{91}}{{100}}}&-{\nicefrac{{7}}{{10}}}&{\nicefrac{{1}}{{100}}}\\\
1&0&0&0&0&0\\\ 0&1&0&0&0&0\\\ 0&0&1&0&0&0\\\
0&0&0&1&0&0\end{bmatrix}\quad\bm{C}_{1}=\begin{bmatrix}-{\nicefrac{{17}}{{30}}}&{\nicefrac{{43}}{{300}}}&{\nicefrac{{41}}{{50}}}&-{\nicefrac{{41}}{{100}}}&-{\nicefrac{{29}}{{100}}}&{\nicefrac{{8}}{{25}}}\\\
{\nicefrac{{5099}}{{5940}}}&-{\nicefrac{{22}}{{25}}}&{\nicefrac{{7}}{{10}}}&-{\nicefrac{{91}}{{100}}}&-{\nicefrac{{7}}{{10}}}&{\nicefrac{{1}}{{100}}}\\\
1&0&1&0&0&0\\\ 0&1&0&1&0&0\\\ 0&0&1&0&3&0\\\ 0&0&0&1&0&3\end{bmatrix}$
$\displaystyle\bm{X}=\begin{bmatrix}\nicefrac{{1}}{{3}}&0&\nicefrac{{2}}{{3}}&0&1&0\\\
0&\nicefrac{{1}}{{3}}&0&\nicefrac{{2}}{{3}}&0&1\end{bmatrix}\quad\bm{Y}=\begin{bmatrix}1&0\\\
0&1\\\ 0&0\\\ 0&0\\\ 0&0\\\ 0&0\end{bmatrix}\>.$

Then,

$\bm{X}\left(z\bm{C}_{1}-\bm{C}_{0}\right)^{-1}\bm{Y}\bm{P}(z)=\bm{I}_{2}\>.$

###  6.3 Lagrange Basis

######  Example 6.5.

(6.118) $\displaystyle\tau$
$\displaystyle=\left[\mathrm{seq}\left(\cos\left(\frac{\pi\cdot
k}{2}\right),k=0..2\right)\right]=\left[1,0,-1\right]$ $\displaystyle\rho$
$\displaystyle=\left[\bm{I}_{2},\bm{I}_{2},\bm{I}_{2}\right]=\left[\begin{bmatrix}1&0\\\
0&1\end{bmatrix},\begin{bmatrix}1&0\\\ 0&1\end{bmatrix},\begin{bmatrix}1&0\\\
0&1\end{bmatrix}\right]$

$\displaystyle\bm{C}_{0}=\begin{bmatrix}0&0&-1&0&-1&0&-1&0\\\
0&0&0&-1&0&-1&0&-1\\\ \nicefrac{{1}}{{2}}&0&1&0&0&0&0&0\\\
0&\nicefrac{{1}}{{2}}&0&1&0&0&0&0\\\ -1&0&0&0&0&0&0&0\\\ 0&-1&0&0&0&0&0&0\\\
2&0&0&0&0&0&-1&0\\\
0&\nicefrac{{1}}{{2}}&0&0&0&0&0&-1\end{bmatrix}\quad\bm{C}_{1}=\begin{bmatrix}0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0\\\ 0&0&1&0&0&0&0&0\\\ 0&0&0&1&0&0&0&0\\\ 0&0&0&0&1&0&0&0\\\
0&0&0&0&0&1&0&0\\\ 0&0&0&0&0&0&1&0\\\ 0&0&0&0&0&0&0&1\end{bmatrix}$
$\displaystyle\bm{X}=\begin{bmatrix}0&0&1&0&1&0&1&0\\\
0&0&0&1&0&1&0&1\end{bmatrix}\quad\bm{Y}=\begin{bmatrix}1&0\\\ 0&1\\\ 0&0\\\
0&0\\\ 0&0\\\ 0&0\end{bmatrix}\>.$

Using the barycentric Lagrange interpolation formula, we construct our matrix
polynomial

$\displaystyle\bm{P}(z)$
$\displaystyle=\begin{bmatrix}\left(z-1\right)z\left(z+1\right)\left(\frac{1}{2\left(z-1\right)}-\frac{1}{z}+\frac{1}{2\left(z+1\right)}\right)&0\\\
0&\left(z-1\right)z\left(z+1\right)\left(\frac{1}{2\left(z-1\right)}-\frac{1}{z}+\frac{1}{2\left(z+1\right)}\right)\end{bmatrix}$
$\displaystyle=\begin{bmatrix}1&0\\\ 0&1\end{bmatrix}$

that corresponds to the given $\tau$ and $\rho$ from equation 6.118.
Therefore, $\bm{P}^{-1}(z)$

$\bm{X}\left(z\bm{C}_{1}-\bm{C}_{0}\right)^{-1}\bm{Y}=\bm{I}_{2}$

###  6.4 Hermite Interpolational Basis

######  Example 6.6 (Polynomial case).

Let

(6.119) $\tau=\left[-1,-\frac{1}{2},\frac{1}{2},1\right]$

and

$z$ | $P(z)$ | $P^{\prime}(z)$ | $P^{\prime\prime}(z)$
---|---|---|---
$\>\>\>\tau_{0}=-1$ | $1$ | $0$ | $0$
$\phantom{-}\tau_{1}=-\frac{1}{2}$ | $1$ |  |
$\tau_{2}=\frac{1}{2}$ | $1$ |  |
$\tau_{3}=1$ | $1$ | $0$ |

Note that this polynomial is identically $1$: its values at all nodes are $1$,
and all derivatives at all nodes are $0$. This demonstrates explicitly that
the degree of the polynomial is not necessarily revealed by the grade, which
here is $\ell=5$. The standard triple is then

$\displaystyle\bm{C}_{0}=\begin{bmatrix}0&0&-1&-1&-1&0&0&-1\\\
\nicefrac{{1}}{{6}}&1&0&0&0&0&0&0\\\ -{\nicefrac{{25}}{{36}}}&1&1&0&0&0&0&0\\\
{\nicefrac{{32}}{{27}}}&0&0&\nicefrac{{1}}{{2}}&0&0&0&0\\\
-{\nicefrac{{32}}{{9}}}&0&0&0&-\nicefrac{{1}}{{2}}&0&0&0\\\
\nicefrac{{1}}{{3}}&0&0&0&0&-1&0&0\\\ {\nicefrac{{11}}{{9}}}&0&0&0&0&1&-1&0\\\
{\nicefrac{{331}}{{108}}}&0&0&0&0&0&1&-1\end{bmatrix}\quad\bm{C}_{1}=\begin{bmatrix}0&0&0&0&0&0&0&0\\\
0&1&0&0&0&0&0&0\\\ 0&0&1&0&0&0&0&0\\\ 0&0&0&1&0&0&0&0\\\ 0&0&0&0&1&0&0&0\\\
0&0&0&0&0&1&0&0\\\ 0&0&0&0&0&0&1&0\\\ 0&0&0&0&0&0&0&1\end{bmatrix}$
$\displaystyle\bm{X}=\begin{bmatrix}0&0&1&1&1&0&0&1\end{bmatrix}\quad\bm{Y}=\begin{bmatrix}1\\\
0\\\ 0\\\ 0\\\ 0\\\ 0\\\ 0\\\ 0\end{bmatrix}\>.$

Using the first barycentric representation, the Hermite interpolation
polynomial of the given data is

$\displaystyle P(z)$
$\displaystyle=\left(z+1\right)^{3}\left(z+\frac{1}{2}\right)\left(z-\frac{1}{2}\right)\left(z-1\right)^{2}\cdot$
$\displaystyle\qquad\qquad\left({\frac{331}{108z+108}}+{\frac{11}{9\left(z+1\right)^{2}}}+\frac{1}{3}\left(z+1\right)^{-3}-{\frac{32}{9z+\frac{9}{2}}}+{\frac{32}{27z-{\frac{27}{2}}}}-{\frac{25}{36z-36}}+\frac{1}{6}\left(z-1\right)^{-2}\right)$
$\displaystyle=1\>,$

as discussed. Therefore the companion pencil has no finite eigenvalues, in
exact arithmetic. Numerically, it can be expected to have eigenvalues around
$O(1/\mu)^{1/7}$ where $\mu$ is the unit roundoff; here the exponent is $7$,
two more than the grade because only two of the spurious eigenvalues at
infinity are detected and removed precisely [23]. Indeed that is what occurs
(calculations not shown here). Returning to the example, calculating the
resolvent form gives

$\bm{X}\left(z\bm{C}_{1}-\bm{C}_{0}\right)^{-1}\bm{Y}=1\>,$

and therefore (due to multiplying this by $P(z)=1$), this shows that the
standard triple for the Hermite interpolating basis is correct.

######  Example 6.7 (Matrix polynomial case).

Let

$\tau=\left[0,1\right]$

and

$z$ | $\bm{P}(z)$ | $\bm{P}^{\prime}(z)$
---|---|---
$\tau_{0}=0$ | $\begin{bmatrix}-1&0\\\ -1&1\end{bmatrix}$ |
$\tau_{1}=1$ | $\begin{bmatrix}0&\phantom{-}1\\\ 1&-1\end{bmatrix}$ | $\begin{bmatrix}\phantom{-}1&-1\\\ -1&\phantom{-}0\end{bmatrix}$

Then, the standard triple is

$\displaystyle\bm{C}_{0}=\begin{bmatrix}0&0&-1&1&0&-1&1&0\\\
0&0&1&0&-1&1&1&-1\\\ 1&0&1&0&0&0&0&0\\\ 0&1&0&1&0&0&0&0\\\ -1&0&1&0&1&0&0&0\\\
0&-1&0&1&0&1&0&0\\\ 1&0&0&0&0&0&0&0\\\
0&1&0&0&0&0&0&0\end{bmatrix}\quad\bm{C}_{1}=\begin{bmatrix}0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0\\\ 0&0&1&0&0&0&0&0\\\ 0&0&0&1&0&0&0&0\\\ 0&0&0&0&1&0&0&0\\\
0&0&0&0&0&1&0&0\\\ 0&0&0&0&0&0&1&0\\\ 0&0&0&0&0&0&0&1\end{bmatrix}$
$\displaystyle\bm{X}=\begin{bmatrix}0&0&0&0&1&0&1&0\\\
0&0&0&0&0&1&0&1\end{bmatrix}\quad\bm{Y}=\begin{bmatrix}1&0\\\ 0&1\\\ 0&0\\\
0&0\\\ 0&0\\\ 0&0\\\ 0&0\\\ 0&0\end{bmatrix}\>.$

The Hermite interpolating polynomial is

$\bm{P}(z)=\begin{bmatrix}z-1&-2{z}^{2}+3z\\\
-3{z}^{2}+5z-1&2{z}^{2}-4z+1\end{bmatrix}$

and the resolvent form is

$\bm{X}\left(z\bm{C}_{1}-\bm{C}_{0}\right)^{-1}\bm{Y}=\begin{bmatrix}{\dfrac{-2{z}^{2}+4z-1}{6{z}^{4}-21{z}^{3}+23{z}^{2}-8z+1}}&{\dfrac{-2{z}^{2}+3z}{6{z}^{4}-21{z}^{3}+23{z}^{2}-8z+1}}\\\
{\dfrac{-3{z}^{2}+5z-1}{6{z}^{4}-21{z}^{3}+23{z}^{2}-8z+1}}&{\dfrac{-z+1}{6{z}^{4}-21{z}^{3}+23{z}^{2}-8z+1}}\end{bmatrix}\>.$

Then,

$\bm{X}\left(z\bm{C}_{1}-\bm{C}_{0}\right)^{-1}\bm{Y}\bm{P}(z)=\bm{I}_{2}\>,$

which indicates that the standard triples is correct.

##  7 Concluding remarks

The generalized standard triple (or standard quadruple, if you prefer) that we
propose in this paper for convenience in algebraic linearization may have
other uses. As pointed out on p. 28 of [17] many of the properties stated in
that work for monic polynomials are valid for non-monic polynomials with the
appropriate changes made. We have not attempted a comprehensive categorization
of those changes for other purposes.

We have established that these generalized standard triples allow the
resolvent representation for the matrix polynomial, equation (1.6), that is
useful for algebraic linearization.

## Acknowledgments

We acknowledge the support of Western University, The National Science and
Engineering Research Council of Canada, the Ontario Graduate Scholarship (OGS)
program, the University of Alcalá, the Ontario Research Centre of Computer
Algebra, and the Rotman Institute of Philosophy. Part of this work was
developed while RMC was visiting the University of Alcalá, in the frame of the
project Giner de los Rios. The authors would also like to thank Peter
Lancaster for teaching RMC long ago the value of the $5\times 5$ example.
Similarly we thank John C. Butcher for the proper usage of the word
“interpolational”. In addition, thank Françoise Tisseur for her thorough
comments on an earlier version of this paper.

## References

  * [1] Amir Amiraslani, Robert M. Corless, and Peter Lancaster.  Linearization of matrix polynomials expressed in polynomial bases.  IMA Journal of Numerical Analysis, 29(1):141–157, 2008.
  * [2] Stephen Barnett.  Some applications of the comrade matrix.  International Journal of Control, 21(5):849–855, 1975.
  * [3] Jean-Paul Berrut and Lloyd N. Trefethen.  Barycentric Lagrange interpolation.  SIAM Review, 46(3):501–517, 2004.
  * [4] Timo Betcke, Nicholas J. Higham, Volker Mehrmann, Christian Schröder, and Françoise Tisseur.  NLEVP: A collection of nonlinear eigenvalue problems.  ACM Transactions on Mathematical Software (TOMS), 39(2):7, 2013\.
  * [5] J. M. Carnicer, Y. Khiar, and J. M. Peña.  Optimal stability of the Lagrange formula and conditioning of the Newton formula.  Journal of Approximation Theory, 2017.
  * [6] Eunice Y. S. Chan, Robert M. Corless, Laureano Gonzalez-Vega, J. Rafael Sendra, and Juana Sendra.  Algebraic linearizations for matrix polynomials.  Linear Algebra and its Applications, 563:373–399, 2019.
  * [7] Robert M. Corless and Nicolas Fillion.  Polynomial and rational interpolation.  In A Graduate Introduction to Numerical Methods, pages 331–401. Springer, 2013.
  * [8] Robert M Corless, Nargol Rezvani, and Amirhossein Amiraslani.  Pseudospectra of matrix polynomials that are expressed in alternative bases.  Mathematics in Computer Science, 1(2):353–374, 2007.
  * [9] Robert M. Corless and Stephen M. Watt.  Bernstein bases are optimal, but, sometimes, Lagrange bases are better.  In Proceedings of SYNASC, Timisoara, pages 141–153. MIRTON Press, 2004.
  * [10] Froilán M Dopico, Piers W. Lawrence, Javier Pérez, and Paul Van Dooren.  Block Kronecker linearizations of matrix polynomials and their backward errors.  Numerische Mathematik, 140(2):373–426, 2018.
  * [11] Rida T. Farouki.  The Bernstein polynomial basis: A centennial retrospective.  Computer Aided Geometric Design, 29(6):379–419, 2012.
  * [12] Rida T. Farouki and T. Goodman.  On the optimal stability of the Bernstein basis.  Mathematics of Computation of the American Mathematical Society, 65(216):1553–1566, 1996.
  * [13] Rida T. Farouki and V. T. Rajan.  On the numerical condition of polynomials in Bernstein form.  Computer Aided Geometric Design, 4(3):191–216, 1987.
  * [14] Heike Faßbender and Philip Saltenberger.  On vector spaces of linearizations for matrix polynomials in orthogonal bases.  Linear Algebra and its Applications, 525:59–83, 2017.
  * [15] Walter Gautschi.  Orthogonal polynomials in MATLAB: Exercises and Solutions, volume 26.  SIAM, 2016.
  * [16] Israel Gohberg, Peter Lancaster, and Leiba Rodman.  Indefinite Linear Algebra and Applications.  Springer, 2005.
  * [17] Israel Gohberg, Peter Lancaster, and Leiba Rodman.  Matrix Polynomials.  SIAM Classics in Applied Mathematics, 2009.
  * [18] Stefan Güttel and Françoise Tisseur.  The nonlinear eigenvalue problem.  Acta Numerica, 26:1–94, 2017.
  * [19] David J. Jeffrey and Robert M. Corless.  Linear algebra in maple.  In Leslie Hogben, editor, Handbook of Linear Algebra, chapter 89. Chapman and Hall/CRC, 2013.
  * [20] Guðbjörn F. Jónsson.  Eigenvalue methods for accurate solution of polynomial equations. PhD thesis, Center for Applied Mathematics, Cornell University, Ithaca, NY, 2001.
  * [21] Guðbjörn F. Jónsson and Stephen Vavasis.  Solving polynomials with small leading coefficients.  SIAM Journal on Matrix Analysis and Applications, 26(2):400–414, 2004.
  * [22] V. N. Kublanovskaya.  Methods and algorithms of solving spectral problems for polynomial and rational matrices.  Journal of Mathematical Sciences, 96(3):3085–3287, 1999.
  * [23] Piers W. Lawrence and Robert M. Corless.  Numerical stability of barycentric Hermite root-finding.  In Proceedings of the 2011 International Workshop on Symbolic-Numeric Computation, pages 147–148. ACM, 2012.
  * [24] Jörg Liesen and Christian Mehl.  Matrix polynomials.  In Leslie Hogben, editor, Handbook of Linear Algebra, chapter 18. Chapman and Hall/CRC, 2013.
  * [25] D. Steven Mackey, Niloufer Mackey, Christian Mehl, and Volker Mehrmann.  Vector spaces of linearizations for matrix polynomials.  SIAM Journal on Matrix Analysis and Applications, 28(4):971–1004, 2006.
  * [26] D. Steven Mackey, Niloufer Mackey, and Françoise Tisseur.  Polynomial eigenvalue problems: Theory, computation, and structure.  In Numerical Algebra, Matrix Theory, Differential-Algebraic Equations and Control Theory, pages 319–348. Springer, 2015.
  * [27] D. Steven Mackey and Vasilije Perović.  Linearizations of matrix polynomials in Bernstein bases.  Linear Algebra and its Applications, 501:162–197, 2016.
  * [28] Yuji Nakatsukasa, Vanni Noferini, and Alex Townsend.  Vector spaces of linearizations for matrix polynomials: a bivariate polynomial approach.  SIAM Journal on Matrix Analysis and Applications, 38(1):1–29, 2017\.
  * [29] Leonardo Robol, Raf Vandebril, and Paul Van Dooren.  A framework for structured linearizations of matrix polynomials in various bases.  SIAM Journal on Matrix Analysis and Applications, 38(1):188–216, 2017.
  * [30] Arne Storjohann.  An ${O}(n^{3})$ algorithm for the Frobenius normal form.  In Proceedings of the 1998 international symposium on Symbolic and algebraic computation, pages 101–105. ACM, 1998.
  * [31] Olga Taussky and Hans Zassenhaus.  On the similarity transformation between a matrix and its transpose.  Pacific Journal of Mathematics, 9(3):893–896, 1959.
  * [32] Roel Van Beeumen, Wim Michiels, and Karl Meerbergen.  Linearization of Lagrange and Hermite interpolating matrix polynomials.  IMA Journal of Numerical Analysis, 35(2):909–930, 2015.