# Christoffel transformations for
matrix orthogonal polynomials in the real line and
the non-Abelian 2D Toda lattice hierarchy

Carlos Álvarez-Fernández  Departamento de Métodos Cuantitativos, Universidad
Pontificia Comillas, Calle de Alberto Aguilera 23, 28015-Madrid, Spain
<EMAIL_ADDRESS>,  Gerardo Ariznabarreta  Departamento de Física Teórica
II (Métodos Matemáticos de la Física), Universidad Complutense de Madrid,
Ciudad Universitaria, Plaza de Ciencias 1, 28040-Madrid, Spain
<EMAIL_ADDRESS>,  Juan Carlos García-Ardila  Departamento de Matemáticas,
Universidad Carlos III de Madrid, Avenida Universidad 30, 28911 Leganés, Spain
<EMAIL_ADDRESS>,  Manuel Mañas  Departamento de Física Teórica II
(Métodos Matemáticos de la Física), Universidad Complutense de Madrid, Ciudad
Universitaria, Plaza de Ciencias 1, 28040 Madrid, Spain<EMAIL_ADDRESS>and  Francisco Marcellán  Departamento de Matemáticas, Universidad Carlos III
de Madrid and Instituto de Ciencias Matemáticas (ICMAT), Avenida Universidad
30, 28911 Leganés, Spain<EMAIL_ADDRESS>
###### Abstract.

Given a matrix polynomial $W(x)$, matrix bi-orthogonal polynomials with
respect to the sesquilinear form

$\displaystyle\left\langle{P(x),Q(x)}\right\rangle_{W}=\int
P(x)W(x)\operatorname{d}\mu(x)(Q(x))^{\top},\quad P,Q\in\mathbb{R}^{p\times
p}[x],$

where $\mu(x)$ is a matrix of Borel measures supported in some infinite subset
of the real line, are considered. Connection formulas between the sequences of
matrix bi-orthogonal polynomials with respect to
$\left\langle{\cdot,\cdot}\right\rangle_{W}$ and matrix polynomials orthogonal
with respect to $\mu(x)$ are presented. In particular, for the case of
nonsingular leading coefficients of the perturbation matrix polynomial $W(x)$
we present a generalization of the Christoffel formula constructed in terms of
the Jordan chains of $W(x)$. For perturbations with a singular leading
coefficient several examples by Durán et al are revisited. Finally, we extend
these results to the non-Abelian 2D Toda lattice hierarchy.

###### Key words and phrases:

Matrix orthogonal polynomials, Block Jacobi matrices, Darboux–Christoffel
transformation, Block Cholesky decomposition, Block LU decomposition, quasi-
determinants, non-Abelian Toda hierarchy

###### 1991 Mathematics Subject Classification:

42C05,15A23

GA thanks financial support from the Universidad Complutense de Madrid Program
“Ayudas para Becas y Contratos Complutenses Predoctorales en España 2011"

MM & FM thanks financial support from the Spanish “Ministerio de Economía y
Competitividad" research project MTM2012-36732-C03-01, _Ortogonalidad y
aproximación; teoría y aplicaciones_

###### Contents

  1. 1 Introduction
    1. 1.1 Historical background and state of the art
    2. 1.2 Objectives, results and layout of the paper
    3. 1.3 On spectral theory of matrix polynomials
    4. 1.4 On orthogonal matrix polynomials
  2. 2 Connection formulas for Darboux transformations of Christoffel type
    1. 2.1 Connection formulas for bi-orthogonal polynomials
    2. 2.2 Connection formulas for the Christoffel–Darboux kernel
  3. 3 Monic matrix polynomial perturbations
    1. 3.1 The Christoffel formula for matrix bi-orthogonal polynomials
    2. 3.2 Degree one monic matrix polynomial perturbations
    3. 3.3 Examples
  4. 4 Singular leading coefficient matrix polynomial perturbations
  5. 5 Extension to non-Abelian 2D Toda hierarchies
    1. 5.1 Block Hankel moment matrices vs multi-component Toda hierarchies
    2. 5.2 The Christoffel transformation for the non-Abelian 2D Toda hierarchy

##  1\. Introduction

This paper is devoted to the extension of the Christoffel formula to the
Matrix Orthogonal Polynomials on the Real Line (MOPRL) and the non-Abelian 2D
Toda lattice hierarchy.

###  1.1. Historical background and state of the art

In 1858 the German mathematician Elwin Christoffel [32] was interested, in the
framework of Gaussian quadrature rules, in finding explicit formulas relating
the corresponding sequences of orthogonal polynomials with respect to two
measures $\operatorname{d}\mu$ (in the Christoffel’s discussion was just the
Lebesgue measure $\operatorname{d}\mu=\operatorname{d}x$ ) and
$d\hat{\mu}(x)=p(x)d\mu(x)$, with $p(x)=(x-q_{1})\cdots(x-q_{N})$ a signed
polynomial in the support of $\operatorname{d}\mu$, as well as the
distribution of their zeros as nodes in such quadrature rules, see [101]. The
so called Christoffel formula is a very elegant formula from a mathematical
point of view, and is a classical result which can be found in a number of
orthogonal polynomials textbooks, see for example [29, 95, 51]. Despite these
facts, we must mention that for computational and numerical purposes it is not
so practical, see [51]. These transformations have been extended from measures
to the more general setting of linear functionals. In the theory of orthogonal
polynomials with respect to a moment linear functional
$u\in\big{(}\mathbb{R}[x]\big{)}^{\prime}$, an element of the algebraic dual
(which coincides with the topological dual) of the linear space
$\mathbb{R}[x]$ of polynomials with real coefficients. Given a positive
definite linear moment functional, i.e.
$\Big{|}\big{(}\left\langle{u,x^{n+m}}\right\rangle\big{)}_{n,m=0}^{k}\Big{|}>0$,
$\forall k\in\mathbb{Z}_{+}:=\\{0,1.2,\dots\\}$, there exists a nontrivial
probability measure $\mu$ such that (see [9], [29], [95])
$\left\langle{u,x^{n}}\right\rangle=\int x^{m}\operatorname{d}\mu(x)$. Given a
moment linear functional $u$, its canonical or elementary Christoffel
transformation is a new moment functional given by $\hat{u}=(x-a)u$ with
$a\in\mathbb{R}$, see [24, 29, 97]. The right inverse of a Christoffel
transformation is called the Geronimus transformation. In other words, if you
have a moment linear functional $u$, its elementary or canonical Geronimus
transformation is a new moment linear functional $\check{u}$ such that
$(x-a)\check{u}=u$. Notice that in this case $\check{u}$ depends on a free
parameter, see [54, 78]. The right inverse of a general Christoffel
transformation is said to be a multiple Geronimus transformation, see [40].
All these transformations are refered as Darboux transformations, a name that
was first given in the context of integrable systems in [77]. In 1878 the
French mathematician Gaston Darboux, when studying the Sturm–Liouville theory
in [35], explicitly treated these transformations, which appeared for the
first time in [82]. In the framework of orthogonal polynomials on the real
line, such a factorization of Jacobi matrices has been studied in [24] and
[97]. They also play an important role in the analysis of bispectral problems,
see [58] and [57].

An important aspect of canonical Christoffel transformations is its relations
with $LU$ factorization (and its flipped version, an $UL$ factorization) of
the Jacobi matrix. A sequence of monic polynomials
$\\{P_{n}(x)\\}_{n=0}^{\infty}$ associated with a nontrivial probability
measure $\mu$ satisfies a three term recurrence relation (TTRR, in short)
$xP_{n}(x)=P_{n+1}(x)+b_{n}P_{n}(x)+a_{n}^{2}P_{n-1}(x),n\geq 0,$ with the
convention $P_{-1}(x)=0.$ If we denote by
$P(x)=[P_{0}(x),P_{1}(x),\dots]^{T}$, then the matrix representation of the
multiplication operator by $x$ is directly deduced from the TTRR and reads
$xP(x)=JP(x)$, where $J$ is a tridiagonal semi-infinite matrix such that the
entries in the upper diagonal are the unity. Assuming that $a$ is a real
number off the support of $\mu,$ then you have a factorization $J-aI=LU$,
where $L$ and $U$ are, respectively, lower unitriangular and upper triangular
matrices. The important observation is that the matrix $\hat{J}$ defined by
$\hat{J}-aI=UL$ is again a Jacobi matrix and the corresponding sequence of
monic polynomials $\\{\hat{P}_{n}(x)\\}_{n=0}^{\infty}$ associated with the
multiplication operator defined by $\hat{J}$ is orthogonal with respect to the
canonical Christoffel transformation of the measure $\mu$ defined as above.

For a moment linear functional $u$, the Stieltjes function
$S(x):=\sum_{n=0}^{\infty}\frac{\left\langle{u,x^{n}}\right\rangle}{x^{n+1}}$
plays an important role in the theory of orthogonal polynomials, due to its
close relation with the measure associated to $u$ as well as its (rational)
Padé Approximation, see [23, 63]. If you consider the canonical Christoffel
transformation $\hat{u}$ of the linear functional $u$, then its Stieltjes
function is $\hat{S}(x)=(x-a)S(x)-u_{0}$. This is a particular case of the
spectral linear transformations studied in [102].

Given a bilinear form
$L:\mathbb{R}[x]\times\mathbb{R}[x]\rightarrow\mathbb{R}$ one could consider
the following non-symmetric and symmetric bilinear perturbations

$\displaystyle\tilde{L}_{1}(p,q)=$ $\displaystyle L(wp,q),$
$\displaystyle\tilde{L}_{2}(p,q)=$ $\displaystyle L(p,wq),$
$\displaystyle\hat{L}(p,q)=$ $\displaystyle L(wp,wq),$

where $w(x)$ is a polynomial. The study of these perturbations can be found in
[25]. Taking into account the matrix representation of the multiplication
operator by $z$ is a Hessenberg matrix, the authors establish a relation
between the Hessenberg matrices associated with the initial and the perturbed
functional by using $LU$ and $QR$ factorization. They also give some algebraic
relations between the sequences of orthogonal polynomials associated with
those bilinear forms. The above perturbations can be seen as an extension of
the Christoffel transformation for bilinear forms. When the bilinear form is
defined by a nontrivial probability measure supported on the unit circle,
Christoffel transformations have been studied in [26] in the framework of CMV
matrices, i.e. the matrix representation of the multiplication operator by $z$
in terms of an orthonormal Laurent polynomial basis. Therein, the authors
state the explicit relation between the sequences of orthonormal Laurent
polynomials associated with a measure and its Christoffel transformation, as
well as its link with QR factorizations of such CMV matrices.

The theory of scalar orthogonal polynomials with respect to probability
measures supported either on the real line or the unit circle is a standard
and classic topic in approximation theory and it also has remarkable
applications in many domains as discrete mathematics, spectral theory of
linear differential operators, numerical integration, integrable systems,
among others. Some extensions of such a theory have been developed more
recently. One of the most exciting generalizations appears when you consider
non-negative Hermitian-valued matrix of measures of size $p\times p$ on a
$\sigma$-algebra of subsets of a space $\Omega$ such that each entry is
countably additive and you are interested in the analysis of the Hilbert space
of matrix valued functions of size $p\times p$ under the inner product
associated with such a matrix of measures. This question appears in the
framework of weakly stationary processes, see [89]. Notice that such an inner
product pays the penalty of the non-commutativity of matrices as well as the
existence of singular matrices with respect to the scalar case. By using the
standard Gram-Schmidt method for the canonical linear basis of the linear
space of polynomials with matrix coefficients a theory of matrix orthogonal
polynomials can be studied. The paper by M. G. Krein [68] is credited as the
first contribution in this topic. Despite they have been sporadically studied
during the last half century, there is an exhaustive bibliography focused on
inner products defined on the linear space of polynomials with matrix
coefficients as well as on the existence and analytic properties of the
corresponding sequences of matrix orthogonal polynomials in the real line (see
[43], [42], [80], [88], [94]) and their applications in Gaussian quadrature
for matrix-valued functions ([93]), scattering theory ([14], [53]) and system
theory ([50]). The work [33] constitutes an updated overview on these topics.

But, more recently, an intensive attention was paid to the spectral analysis
of second order linear differential operators with matrix polynomials as
coefficients. This work was motivated by the Bochner’s characterization of
classical orthogonal polynomials (Hermite, Laguerre and Jacobi) as
eigenfunctions of second order linear differential equations with polynomial
coefficients. The matrix case gives a more rich set of solutions. From the
pioneering work [44] some substantial progress has been done in the study of
families of matrix orthogonal polynomials associated to second order linear
differential operators as eigenfunctions and their structural properties (see
[43], [59], [60] as well as the survey [45]). Moreover, in [27] the authors
showed that there exist sequences of orthogonal polynomials satisfying a first
order linear matrix differential equation that constitutes a remarkable
difference with the scalar case where such a situation does not appear. The
spectral problem for second order linear difference operators with polynomial
coefficients has been considered in [13] as a first step in the general
approach. Therein four families of matrix orthogonal polynomials (as matrix
relatives of Charlier, Meixner, Krawtchouk scalar polynomials and another one
that seems not have any scalar relative) are obtained as illustrative examples
of the method described therein.

It is also a remarkable fact that matrix orthogonal polynomials appear in the
analysis of non standard inner products in the scalar case. Indeed, from the
study of higher order recurrence relations that some sequences of orthogonal
polynomials satisfy (see [44] where the corresponding inner product is
analyzed as an extension of the Favard’s theorem and [48], where the
connection with matrix orthogonal polynomials is stated), to the relation
between standard scalar polynomials associated with measures supported on
harmonic algebraic curves and matrix orthogonal polynomials deduced by a
splitting process of the first ones (see [74]) you get an extra motivation for
the study of matrix orthogonal polynomials. Matrix orthogonal polynomials
appear in the framework of orthogonal polynomials in several variables when
the lexicographical order is introduced. Notice that in such a case, the
moment matrix has a Hankel block matrix where each block is a Hankel matrix
i.e. it has a doubly Hankel structure, see [39].

Concerning spectral transformations, in [40] the authors show that the so
called multiple Geronimus transformations of a measure supported in the real
line yield a simple Geronimus transformation for a matrix of measures. This
approach is based on the analysis of general inner products
$\langle\cdot,\cdot\rangle$ such that the multiplication by a polynomial
operator $h$ is symmetric and satisfies an extra condition $\langle
h(x)p(x),q(x)\rangle=\int p(x)q(x)\operatorname{d}\mu(x),$ where $\mu$ is a
nontrivial probability measure supported on the real line. The connection
between the Jacobi matrix associated to the sequence of scalar polynomials
with respect to $\mu$ and the Hessenberg matrix associated with the
multiplication operator by $h$ is given in terms of the so called $UL$
factorizations. Notice that the connection between the Darboux process and the
noncommutative bispectral problem has been discussed in [56]. The analysis of
perturbations on the entries of the matrix of moments from the point of view
of the relations between the corresponding sequences of matrix orthogonal
polynomials was done in [30].

The seminal work of the Japanese mathematician Mikio Sato [91, 92] and later
on of the Kyoto school [36, 37, 38] settled the basis for a Grasmannian and
Lie group theoretical description of integrable hierarchies. Not much later
Motohico Mulase [83] gave a mathematical description of factorization
problems, dressing procedure, and linear systems as the keys for
integrability. It was not necessary to wait too long, in the development of
integrable systems theory, to find multicomponent versions of the integrable
Toda equations, [98, 99, 100] which later on played a prominent role in the
connection with orthogonal polynomials and differential geometry. The
multicomponent versions of the KP hierachy were analyzed in [21, 22] and [62,
70, 71] and in [72, 73] we can find a further study of the multi-component
Toda lattice hierarchy, block Hankel/Toeplitz reductions, discrete flows,
additional symmetries and dispersionless limits. For the relation with
multiple orthogonal polynomials see [8, 11].

The work of Mark Adler and Pierre van Moerbeke was fundamental to the
connection between integrable systems and orthogonal polynomials. They showed
that the Gauss–Borel factorization problem is the keystone for this
connection. In particular, their studies in the papers on the 2D Toda
hierarchy and what they called the discrete KP hierarchy [3, 4, 5, 6, 7]
clearly established –from a group-theoretical setup– why standard
orthogonality of polynomials and integrability of nonlinear equations of Toda
type where so close.

The relation of multicomponent Toda systems or non-Abelian versions of Toda
equations with matrix orthogonal polynomials was studied, for example, in [80,
11] (on the real line) and in [81, 17] (on the unit circle).

The approach to the Christoffel transformations in this paper, which is based
on the Gauss–Borel factorization problem, has been used before in different
contexts. It has been applied for the construction of discrete integrable
systems connected with orthogonal polynomials of diverse types,

  1. i)

The case of multiple orthogonal polynomials and multicomponent Toda was
analyzed in [12].

  2. ii)

In [15] we dealt with the case of matrix orthogonal Laurent polynomials on the
circle and CMV orderings.

  3. iii)

For orthogonal polynomials in several real variables see [16, 17] and [18] for
orthogonal polynomials on the unit torus and the multivariate extension of the
CMV ordering.

It is well known that there is a deep connection between discrete integrable
systems and Darboux transformations of continuous integrable systems, see for
example [41]. Finally, let us comment that, in the realm of several variables,
in [17, 18, 19] one can find extensions of the Christoffel formula to the
multivariate scenario with real variables and on the unit torus, respectively.

###  1.2. Objectives, results and layout of the paper

In this contribution, we focus our attention on the study of Christoffel
transformations (Darboux transformations in the language of integrable systems
[77], or Lévy transformations in the language of differential geometry [49])
for matrix sesquilinear forms. More precisely, given a matrix of measures
$\mu(x)$ and a matrix polynomial $W(x)$ we are going to deal with the
following matrix sesquilinear forms

(1) $\displaystyle\left\langle{P(x),Q(x)}\right\rangle_{W}=\int
P(x)W(x)\operatorname{d}\mu(x)(Q(x))^{\top}.$

We will first focus our attention on the existence of matrix bi-orthogonal
polynomials with respect to the sesquilinear form
$\left\langle{\cdot,\cdot}\right\rangle_{W}$ under some assumptions about the
matrix polynomial $W$. Once this is done, the next step will be to find an
explicit representation of such bi-orthogonal polynomials in terms of the
matrix orthogonal polynomials with respect to the matrix of measures
$\operatorname{d}\mu(x)$. We start with what we call connection formulas in
Proposition 18.

One of the main achievements of this paper is Theorem 2 where we extend the
Christoffel formula to MOPRL with a perturbation given by an arbitrary degree
monic matrix polynomial. For that aim we use the rich spectral theory
available today for these type of polynomials, in particular tools like root
polynomials and Jordan chains will be extremely useful, see [88, 76].
Following [59, 60, 27] some applications to the analysis of matrix orthogonal
polynomials which are eigenfunctions of second order linear differential
operators and related to polynomial perturbations of diagonal matrix of
measures $\operatorname{d}\mu(x)$ leading to sesquilinear forms as
$\left\langle{\cdot,\cdot}\right\rangle_{W}$ will be considered.

Next, to have a better understanding of the singular leading coefficient case
we concentrate on the study of some cases which generalize important examples
given by Alberto Grünbaum and Antonio Durán in [45, 46], in relation again
with second order linear differential operators.

Finally, we see that these Christoffel transformations extend to more general
scenarios in Integrable Systems Theory. In these cases we find the non-Abelian
Toda hierarchy which is relevant in string theory. In general, we have lost
the block Hankel condition, and we do not have anymore a matrix of measures
but only a sesquilinear form. We show that Theorem 2 also holds in this
general situation. At this point we must stress that for the non-Abelian Toda
equation we can find Darboux transformations (or Christoffel transformations)
in [90], see also [84], which contemplate only what it are called elementary
transformations and their iteration. Evidently, their constructions do not
cover by far what Theorem 2 does. There are many matrix polynomials that do
not factor in terms of linear matrix polynomials and, therefore, they cannot
be studied by means of the results in [90, 84]. We have been fortunate to have
at our disposal the spectral theory of [88, 76] that at the moment of the
publication of [90] was not so well known and under construction.

The layout of the paper is as follows. We continue this introduction with two
subsections that give the necessary background material regarding the spectral
theory of matrix polynomials and also of matrix orthogonal polynomials. Then,
in §2, we give the connection formulas for bi-orthogonal polynomials and for
the Christoffel–Darboux kernel, being this last result relevant to find the
dual polynomials in the family of bi-orthogonal polynomials. We continue in §3
discussing the non singular leading coefficient case, i.e., the monic matrix
polynomial perturbation. We find the Christoffel formula for matrix bi-
orthogonal polynomials and, as an example, we consider the degree one monic
matrix polynomial perturbations. We dedicate the rest of this section to
discuss some examples. In §4 we start the exploration of the singular leading
coefficient matrix polynomial perturbations and, despite we do not give a
general theory, we have been able to successfully discuss some relevant
examples. Finally, §5 is devoted to the study of the extension of the previous
results to the non-Abelian 2D Toda lattice hierarchy.

###  1.3. On spectral theory of matrix polynomials

Here we give some background material regarding matrix polynomials. For
further reading we refer the reader to [55]

######  Definition 1.

Let $A_{0},A_{1}\cdots,A_{N}\in\mathbb{R}^{p\times p}$ be square matrices with
real entries. Then

(2) $\displaystyle W(x)=A_{N}x^{N}+A_{N-1}x^{N-1}+\cdots+A_{1}x+A_{0}$

is said to be a matrix polynomial of degree $N$, $\deg W=N$. The matrix
polynomial is said to be monic when $A_{N}=I_{p}$, where
$I_{p}\in\mathbb{R}^{p\times p}$ denotes the identity matrix. The linear space
of matrix polynomials with coefficients in $\mathbb{R}^{p\times p}$ will be
denoted by $\mathbb{R}^{p\times p}[x]$.

######  Definition 2.

We say that a matrix polynomial $W$ as in (2) is monic normalizable if $\det
A_{N}\neq 0$ and say that $\tilde{W}(x):=A_{N}^{-1}W(x)$ is its monic
normalization.

######  Definition 3.

The spectrum, or the set of eigenvalues, $\sigma(W)$ of a matrix polynomial
$W$ is the zero set of $\det W(x)$, i.e.

$\displaystyle\sigma(W):=Z(W)=\\{a\in\mathbb{C}:\det W(a)=0\\}.$

######  Proposition 1.

A monic normalizable matrix polynomial $W(x)$, $\deg W=N$, has $Np$ (counting
multiplicities) eigenvalues or zeros; i.e., we can write

$\displaystyle\det W(x)=\prod_{i=1}^{q}(x-x_{i})^{\alpha_{i}}$

with $Np=\alpha_{1}+\dots+\alpha_{q}$.

######  Remark 1.

In contrast with the scalar case, there exist matrix polynomials which do not
have a unique factorization in terms of degree one factors or even it could
happen that the factorization does not exist. For example, the matrix
polynomial

$\displaystyle W(x)=I_{2}x^{2}-\begin{bmatrix}2&0\\\ 0&2\end{bmatrix}x$

can be written as

$\displaystyle W(x)=\left(I_{2}x-\begin{bmatrix}1&-1\\\
-1&1\end{bmatrix}\right)\left(I_{2}x-\begin{bmatrix}1&1\\\
1&1\end{bmatrix}\right)\ \ \ \ \ \ \text{or}\ \ \ \ \ \
W(x)=\left(I_{2}x-\begin{bmatrix}2&0\\\ 0&2\end{bmatrix}\right)I_{2}x,$

but the polynomial

$W(x)=I_{2}x^{2}-\begin{bmatrix}0&0\\\ 1&0\end{bmatrix}$

cannot be factorized in terms of degree one matrix polynomials.

######  Definition 4.

  1. i)

Two matrix polynomials $W_{1},W_{2}\in\mathbb{R}^{m\times m}[x]$ are said to
be equivalent $W_{1}\sim W_{2}$ if there exist two matrix polynomials
$E,F\in\mathbb{R}^{m\times m}[x]$, with constant determinants (not depending
on $x$), such that $W_{1}(x)=E(x)W_{2}(x)F(x)$.

  2. ii)

A degree one matrix polynomial $I_{Np}x-A\in\mathbb{R}^{Np\times Np}$ is
called a linearization of a monic matrix polynomial $W\in\mathbb{R}^{p\times
p}[x]$ if

$\displaystyle I_{Np}x-A\sim\begin{bmatrix}W(x)&0\\\
0&I_{(N-1)p}\end{bmatrix}$

######  Definition 5.

Given a matrix polynomial $W(x)=I_{p}x^{N}+A_{N-1}x^{N-1}+\cdots+A_{0}$ its
companion matrix $C_{1}\in\mathbb{R}^{Np\times Np}$ is

$\displaystyle C_{1}:=\begin{bmatrix}0&I_{p}&0&\dots&0\\\
0&0&I_{p}&\ddots&0\\\ \vdots&\vdots&\ddots&\ddots&\\\ 0&0&0&&I_{p}\\\
-A_{0}&-A_{1}&-A_{2}&\dots&-A_{N-1}\end{bmatrix}.$

The companion matrix plays an important role in the study of the spectral
properties of a matrix polynomial $W(x)$, see for example [55, 75] and [76]

######  Proposition 2.

Given a monic matrix polynomial $W(x)=I_{p}x^{N}+A_{N-1}x^{N-1}+\cdots+A_{0}$
its companion matrix $C_{1}$ provides a linearization

$\displaystyle I_{Np}x-C_{1}\sim\begin{bmatrix}W(x)&0\\\
0&I_{(N-1)p}\end{bmatrix}$

where

$\displaystyle E(x)$
$\displaystyle=\begin{bmatrix}B_{N-1}(x)&B_{N-2}(x)&B_{N-3}(x)&\dots&B_{1}(x)&B_{0}(x)\\\
-I_{p}&0&0&\dots&0&0\\\ 0&-I_{p}&0&\dots&0&0\\\ 0&0&-I_{p}&&0&0\\\
\vdots&&&\ddots&\vdots\\\ 0&0&0&&-I_{p}&0\end{bmatrix},$ $\displaystyle F(x)$
$\displaystyle=\begin{bmatrix}I_{p}&0&0&\dots&0&0\\\
-I_{p}x&I_{p}&0&\dots&0&0\\\ 0&-I_{p}x&I_{p}&\ddots&0&0\\\
\vdots&\ddots&\ddots&\ddots&\vdots&\vdots\\\ 0&0&0&\ddots&I_{p}&0\\\
0&0&0&&-I_{p}x&I_{p}\\\ \end{bmatrix},$

with $B_{0}(x):=I_{p}$, $B_{r+1}(x)=xB_{r}(x)+A_{N-r-1}$, for
$r\in\\{0,1,\dots,N-2\\}$.

From here one deduces the important

######  Proposition 3.

The eigenvalues with multiplicities of a monic matrix polynomial coincide with
those of its companion matrix.

######  Proposition 4.

Any nonsingular matrix polynomial $W(x)\in\mathbb{C}^{m\times m}[x]$, $\det
W(x)\neq 0$, can be represented

$\displaystyle
W(x)=E(x_{0})\operatorname{diag}((x-x_{0})^{\kappa_{1}},\dots,(x-x_{0})^{\kappa_{m}})F(x_{0})$

at $x=x_{0}\in\mathbb{C}$, where $E(x_{0})$ and $F(x_{0})$ are nonsingular
matrices and $\kappa_{1}\leq\dots\leq\kappa_{m}$ are nonnegative integers.
Moreover, $\\{\kappa_{1},\dots,\kappa_{m}\\}$ are uniquely determined by $W$
and they are known as partial multiplicities of $W(x)$ at $x_{0}$.

######  Definition 6.

  1. i)

Given a monic matrix polynomial $W(x)\in\mathbb{R}^{p\times p}[x]$ with
eigenvalues and multiplicities $\\{x_{k},\alpha_{k}\\}_{k=1}^{q}$,
$Np=\alpha_{1}+\dots+\alpha_{q}$, a non-zero vector $v_{k,0}\in\mathbb{C}^{p}$
is said to be an eigenvector with eigenvalue $x_{k}$ whenever
$W(x_{k})v_{k,0}=0$, $v_{k,0}\in\operatorname{Ker}W(x_{k})\neq\\{0\\}$.

  2. ii)

A sequence of vectors $\\{v_{i,0},v_{i,1}\ldots,v_{1,m_{i}-1}\\}$ is said to
be a Jordan chain of length $m_{i}$ corresponding to $x_{i}\in\sigma(W)$ if
$v_{0,i}$ is an eigenvector of $W(x_{i})$ and

$\displaystyle\sum_{r=0}^{j}\frac{1}{r!}\frac{\operatorname{d}^{r}W}{\operatorname{d}x^{r}}\Big{|}_{x=x_{i}}v_{i,j-r}$
$\displaystyle=0,$ $\displaystyle j$ $\displaystyle=\\{0,\ldots,m_{i}-1\\}.$

  3. iii)

A root polynomial at an eigenvalue $x_{0}\in\sigma(W)$ of $W(x)$ is a non-zero
vector polynomial $v(x)\in\mathbb{C}^{p}[x]$ such that $W(x_{0})v(x_{0})=0$.
The multiplicity of this zero will be denoted by $\kappa$.

  4. iv)

The maximal length of a Jordan chain corresponding to the eigenvalue $x_{k}$
is called the multiplicity of the eigenvector $v_{0,k}$ and is denoted by
$m(v_{0,k})$.

The above definition generalizes the concept of Jordan chain for degree one
matrix polynomials [1].

######  Proposition 5.

The Taylor expansion of a root polynomial at a given eigenvalue
$x_{0}\in\sigma(W)$

$\displaystyle v(x)=\sum_{j=0}^{q}v_{j}(x-x_{0})^{j}$

provides a Jordan chain $\\{v_{0},v_{1},\dots,v_{\kappa-1}\\}$.

######  Proposition 6.

Given an eigenvalue $x_{0}\in\sigma(W)$, with multiplicity
$s=\dim\operatorname{Ker}W(x_{0})$, we can construct $s$ root polynomials

$\displaystyle v_{i}(x)=$
$\displaystyle\sum_{j=0}^{\kappa_{i}-1}v_{i,j}(x-x_{0})^{j},$ $\displaystyle
i$ $\displaystyle\in\\{1,\dots,s\\},$

where $v_{i}(x)$ is a root polynomial with the largest order $\kappa_{i}$
among all root polynomials whose eigenvector does not belong to
$\mathbb{C}\\{v_{1,0},\dots,v_{i-1,0}\\}$.

######  Definition 7.

A canonical set of Jordan chains of the monic matrix polynomial $W(x)$
corresponding to the eigenvalue $x_{0}\in\sigma(W)$ is, in terms of the root
polynomials described in Proposition 6, the set of vectors

$\displaystyle\\{v_{1,0}\dots,v_{1,\kappa_{1}-1},\dots,v_{s,0}\dots,v_{s,\kappa_{r}-1}\\}.$

######  Proposition 7.

For a monic matrix polynomial $W(x)$ the lengths
$\\{\kappa_{1},\dots,\kappa_{r}\\}$ of the Jordan chains in a canonical set of
Jordan chains of $W(x)$ corresponding to the eigenvalue $x_{0}$, see
Definition 7, are the nonzero partial multiplicities of $W(x)$ at $x=x_{0}$
described in Proposition 4.

######  Definition 8.

For each eigenvalue $x_{i}\in\sigma(W)$, with multiplicity $\alpha_{i}$ and
$s_{i}=\dim\operatorname{Ker}W(x_{i})$, we choose a canonical set of Jordan
chains

$\displaystyle\big{\\{}v_{j,0}^{(i)},\dots,v_{j,\kappa_{j}^{(i)}-1}^{(i)}\big{\\}},$
$\displaystyle j$ $\displaystyle=1,\dots,s_{i},$

and, consequently, with partial multiplicities satisfying
$\sum_{j=1}^{s_{i}}\kappa_{j}^{(i)}=\alpha_{i}$. Thus, we can consider the
following adapted root polynomials

(3) $\displaystyle
v_{j}^{(i)}(x)=\sum_{r=0}^{\kappa_{j}^{(i)}-1}v_{j,r}^{(i)}(x-x_{i})^{r}.$

######  Proposition 8.

Given a monic matrix polynomial $W(x)$ the adapted root polynomials given in
Definition 8 satisfy

$\displaystyle\frac{\operatorname{d}^{r}}{\operatorname{d}x^{r}}\Big{|}_{x=x_{i}}\big{(}W(x)v_{j}^{(i)}(x)\big{)}$
$\displaystyle=0,$ $\displaystyle r$
$\displaystyle=0,\dots,\kappa^{(i)}_{j}-1,$ $\displaystyle j$
$\displaystyle=1\dots,s_{i}.$

###  1.4. On orthogonal matrix polynomials

Recall that a sesquilinear form $\left\langle{\cdot,\cdot}\right\rangle$ on
the linear space $\mathbb{R}^{p\times p}[x]$ is a map

$\left\langle{\cdot,\cdot}\right\rangle:\mathbb{R}^{p\times
p}[x]\times\mathbb{R}^{p\times p}[x]\longrightarrow\mathbb{R}^{p\times p},$

such that for any triple $P,Q,R\in\mathbb{R}^{p\times p}[x]$ of matrix
polynomials we have

  1. i)

$\left\langle{AP(x)+BQ(x),R(x)}\right\rangle=A\left\langle{P(x),R(x)}\right\rangle+B\left\langle{Q(x),R(x)}\right\rangle$,
$\forall A,B\in\mathbb{R}^{p\times p}$.

  2. ii)

$\left\langle{P(x),AQ(x)+BR(x)}\right\rangle=\left\langle{P(x),Q(x)}\right\rangle
A^{\top}+\left\langle{P(x),R(x)}\right\rangle B^{\top}$, $\forall
A,B\in\mathbb{R}^{p\times p}$.

Here $A^{\top}$ denotes the transpose of $A$, an antiautomorphism of order two
in the ring of matrices.

######  Definition 9.

A sesquilinear form $\left\langle{\cdot,\cdot}\right\rangle$ is said to be non
degenerate if the leading principal sub-matrices of the corresponding Hankel
matrix of moments
$M:=\left(\left\langle{I_{p}x^{i}I,I_{p}x^{j}}\right\rangle\right)^{\infty}_{i,j=0}$
are nonsingular, and nontrivial if $\left\langle{\cdot,\cdot}\right\rangle$ is
a symmetric matrix sesquilinear form and
$\left\langle{P(x),P(x)}\right\rangle$ is a positive definite matrix for all
$P(x)\in\mathbb{R}^{p\times p}[x]$ with nonsingular leading coefficient.

Given a sesquilinear form $\left\langle{\cdot,\cdot}\right\rangle$, two
sequences of polynomials $\big{\\{}P_{n}^{[1]}(x)\big{\\}}_{n=0}^{\infty}$ and
$\big{\\{}P_{n}^{[2]}(x)\big{\\}}_{n=0}^{\infty}$ are said to be bi-orthogonal
with respect to $\left\langle{\cdot,\cdot}\right\rangle$ if

  1. i)

$\deg(P_{n}^{[1]})=\deg(P_{n}^{[2]})=n$ for all $n\in\mathbb{Z}_{+}$.

  2. ii)

$\left\langle{P_{n}^{[1]}(x),P_{m}^{[2]}(x)}\right\rangle=\delta_{n,m}H_{n}$
for all $n,m\in\mathbb{Z}_{+}$

where $H_{n}\neq 0$ and $\delta_{n,m}$ is the Kronecker delta. Here, it is
important to notice the order of the polynomials in the sesquilinear form;
i.e., if $n\neq m$ then
$\left\langle{P_{n}^{[2]}(x),P_{m}^{[1]}(x)}\right\rangle$ could be different
from $0$.

######  Remark 2.

Recall that if A is a positive semidefinite (resp. definite) matrix, then
there exists a unique positive semidefinite (resp. definite) matrix $B$ such
that $B^{2}=A$. $B$ is said to be the square root of $A$ (see [61], Theorem
7.2.6) and we denote it by $B=:A^{1/2}$. As in the scalar case, when
$\left\langle{\cdot,\cdot}\right\rangle$ is a sesquilinear form, we will write
the matrix $\left\langle{P,P}\right\rangle^{1/2}:=\|P\|$.

Let

$\displaystyle\mu=\begin{bmatrix}\mu_{1,1}&\dots&\mu_{1,p}\\\
\vdots&&\vdots\\\ \mu_{p,1}&\dots&\mu_{p,p}\end{bmatrix}$

be a $p\times p$ matrix of Borel measures in $\mathbb{R}$. Given any pair of
matrix polynomials $P(x),Q(x)\in\mathbb{R}^{p\times p}[x]$ we introduce the
following sesquilinear form

$\displaystyle\left\langle{P(x),Q(x)}\right\rangle=\int_{\mathbb{R}}P(x)\operatorname{d}\mu(x)(Q(x))^{\top}.$

In terms of the moments of the matrix of measures $\mu$ we define the matrix
moments as

$\displaystyle
m_{n}:=\int_{\mathbb{R}}x^{n}\operatorname{d}\mu(x)\in\mathbb{R}^{p\times p}$

and arrange them in the semi-infinite block matrix and its $k$-th truncation

$\displaystyle M$ $\displaystyle:=\begin{bmatrix}m_{0}&m_{1}&m_{2}&\cdots\\\
m_{1}&m_{2}&m_{3}&\cdots\\\ m_{2}&m_{3}&m_{4}&\cdots\\\
\vdots&\vdots&\vdots&\cdots\\\ \end{bmatrix},$ $\displaystyle M_{[k]}$
$\displaystyle:=\begin{bmatrix}m_{0}&\cdots&m_{k-1}\\\ \vdots&&\vdots\\\
m_{k-1}&\cdots&m_{2k-2}\end{bmatrix}.$

Following [15] we can prove

######  Proposition 9.

If $\det M_{[k]}\neq 0$ for $k\in\\{1,2,\dots\\}$, then there exists a unique
Gaussian factorization of the moment matrix $M$ given by

$\displaystyle M=S_{1}^{-1}H(S_{2})^{-\top},$

where $S_{1},S_{2}$ are lower unitriangular block matrices and $H$ is a
diagonal block matrix

$\displaystyle S_{i}$ $\displaystyle=\begin{bmatrix}I_{p}&0&0&\dots\\\
(S_{i})_{1,0}&I_{p}&0&\cdots\\\ (S_{i})_{2,0}&(S_{i})_{2,1}&I_{p}&\ddots\\\
&&&\ddots\end{bmatrix},$ $\displaystyle H$
$\displaystyle=\begin{bmatrix}H_{0}&0&0&\cdots\\\ 0&H_{1}&0&\ddots\\\
0&0&H_{2}&\ddots\\\ \vdots&\vdots&\ddots&\ddots\end{bmatrix},$ $\displaystyle
i$ $\displaystyle=1,2,$

with $(S_{i})_{n,m},H_{n}\in\mathbb{R}^{p\times p}$, $\forall
n,m\in\\{0,1,\dots\\}$. If $\mu=\mu^{\top}$ then we are dealing with a
Cholesky block factorization with $S_{1}=S_{2}$ and $H=H^{\top}$.

For $l\geq k$ we will also use the following bordered truncated moment matrix

$\displaystyle M_{[k,l]}^{[1]}$
$\displaystyle:=\left[\begin{array}[]{ccc}m_{0}&\cdots&m_{k-1}\\\
\vdots&&\vdots\\\ m_{k-2}&\cdots&m_{2k-3}\\\\[1.0pt] \hline\cr
m_{l}&\dots&m_{l+k-1}\end{array}\right],$

where we have replaced the last row of blocks,
$\begin{bmatrix}m_{k-1}&\dots&m_{2k-2}\end{bmatrix}$, of the truncated moment
matrix $M_{k}$ by the row of blocks
$\begin{bmatrix}m_{l}&\dots&m_{l+k-1}\end{bmatrix}$. We also need a similar
matrix but replacing the last block column of $M_{k}$ by a column of blocks as
indicated

$\displaystyle M_{[k,l]}^{[2]}$
$\displaystyle:=\left[\begin{array}[]{ccc|c}m_{0}&\cdots&m_{k-2}&m_{l}\\\
\vdots&&\vdots&\vdots\\\ m_{k-1}&\cdots&m_{2k-3}&m_{k+l-1}\end{array}\right].$

Using last quasi-determinants, see [52, 85], we find

######  Proposition 10.

If the last quasi-determinants of the truncated moment matrices are
nonsingular, i.e.,

$\displaystyle\det\Theta_{*}(M_{[k]})\neq$ $\displaystyle 0,$ $\displaystyle
k=1,2,\dots,$

then the Gauss–Borel factorization can be performed and the following
expressions

$\displaystyle H_{k}$ $\displaystyle=\Theta_{*}(M_{[k+1]}),$
$\displaystyle(S_{1}^{-1})_{k,l}$
$\displaystyle=\Theta_{*}(M^{[1]}_{[k,l+1]})\Theta_{*}(M_{[l+1]})^{-1},$
$\displaystyle(S_{2}^{-1})_{k,l}$
$\displaystyle=\big{(}\Theta_{*}(M_{[l+1]})^{-1}\Theta_{*}(M^{[2]}_{[k,l+1]})\big{)}^{\top},$

hold.

######  Definition 10.

We define $\chi(x):=[I_{p},I_{p}x,I_{p}x^{2},\dots]^{\top}$ and the vectors of
matrix polynomials $P^{[1]}=[P_{0}^{[1]},P_{1}^{[1]},\cdots]^{\top}$ and
$P^{[2]}=[P_{0}^{[2]},P_{1}^{[2]},\dots]^{\top}$, where

$\displaystyle P^{[1]}$ $\displaystyle:=S_{1}\chi(x),$ $\displaystyle P^{[2]}$
$\displaystyle:=S_{2}\chi(x).$

######  Proposition 11.

The matrix polynomials $P_{n}^{[i]}(x)$ are monic and $\deg P_{n}^{[i]}=n$,
$i=1,2$.

Observe that the moment matrix can be expressed as

$\displaystyle
M=\int_{\mathbb{R}}\chi(x)\operatorname{d}\mu(x)(\chi(x))^{\top}.$

######  Proposition 12.

The families of monic matrix polynomials
$\big{\\{}P_{n}^{[1]}(x)\big{\\}}_{n=0}^{\infty}$ and
$\big{\\{}P_{n}^{[2]}(x)\big{\\}}_{n=0}^{\infty}$ are bi-orthogonal

$\displaystyle\left\langle{P^{[1]}_{n}(x),P^{[2]}_{m}(x)}\right\rangle$
$\displaystyle=\delta_{n,m}H_{n},$ $\displaystyle n,m$
$\displaystyle\in\mathbb{Z}_{+}.$

If $\mu=\mu^{\top}$ then $P_{n}^{[1]}=P_{n}^{[2]}=:P_{n}$ which in turn
conform an orthogonal set of monic matrix polynomials

$\displaystyle\left\langle{P_{n}(x),P_{m}(x)}\right\rangle$
$\displaystyle=\delta_{n,m}H_{n},$ $\displaystyle n,m$
$\displaystyle\in\mathbb{Z}_{+},$

and we can write $\|P_{n}\|=H_{n}^{1/2}$. These bi-orthogonal relations can be
recasted as

$\displaystyle\int_{\mathbb{R}}P_{n}^{[1]}(x)\operatorname{d}\mu(x)x^{m}$
$\displaystyle=\int_{\mathbb{R}}x^{m}\operatorname{d}\mu(x)(P_{n}^{[2]}(x))^{\top}=H_{n}\delta_{n,m},$
$\displaystyle m\leq n.$

###### Proof.

From the definition of the polynomials and the factorization problem we get

$\displaystyle\int_{\mathbb{R}}P^{[1]}(x)\operatorname{d}\mu(x)(P^{[2]}(x))^{\top}=\int_{\mathbb{R}}S_{1}\chi(x)\operatorname{d}\mu(x)\chi(x)^{T}{S_{2}}^{\top}=S_{1}M(S_{2})^{\top}=H.$

∎

######  Remark 3.

The matrix of measures $\operatorname{d}\mu(x)$ may undergo a similarity
transformation, $\operatorname{d}\mu(x)\mapsto\operatorname{d}\mu_{c}(x)$, and
be conjugate to $\operatorname{d}\mu(x)=B^{-1}\operatorname{d}\mu_{c}(x)B$,
where $B\in\mathbb{R}^{p\times p}$ is a nonsingular matrix. The relation
between the orthogonal polynomials given by these two measures is easily seen
to be

$\displaystyle\begin{aligned}
M&=B^{-1}M_{c}B,&S_{1}&=B^{-1}S_{c,1}B,&H&=B^{-1}H_{c}B,&(S_{2})^{\top}&=B^{-1}(S_{c,2})^{\top}B,\end{aligned}$
$\displaystyle\begin{aligned}
P^{[1]}_{n}&=B^{-1}P^{[1]}_{c,n}B,&(P^{[2]}_{n})^{\top}&=B^{-1}(P^{[2]}_{c,n})^{\top}B.\end{aligned}$

The shift matrix is the following semi-infinite block matrix

$\displaystyle\Lambda:=\begin{bmatrix}0&I_{p}&0&0&\dots\\\
0&0&I_{p}&0&\dots\\\ 0&0&0&I_{p}&\dots\\\
\vdots&\ddots&\ddots&\ddots&\ddots\\\ \end{bmatrix}$

which satisfies the important spectral property

$\displaystyle\Lambda\chi(x)=x\chi(x).$

######  Proposition 13.

The block Hankel symmetry of the moment matrix can be written as

$\displaystyle\Lambda M=M\Lambda^{\top}.$

Notice that this symmetry completely characterizes Hankel block matrices.

######  Proposition 14.

We have the following last quasi-determinantal expressions

$\displaystyle P^{[1]}_{n}(x)$
$\displaystyle=\Theta_{*}\begin{bmatrix}m_{0}&m_{1}&\cdots&m_{n-1}&I_{p}\\\
m_{1}&m_{2}&\cdots&m_{n}&I_{p}x\\\ \vdots&\vdots&&\vdots&\vdots\\\
m_{n-1}&m_{n}&\cdots&m_{2n-2}&I_{p}x^{n-1}\\\
m_{n}&m_{n+1}&\cdots&m_{2n-1}&I_{p}x^{n}\end{bmatrix},$
$\displaystyle(P^{[2]}_{n}(x))^{\top}$
$\displaystyle=\Theta_{*}\begin{bmatrix}m_{0}&m_{1}&\cdots&m_{n-1}&m_{n}\\\
m_{1}&m_{2}&\cdots&m_{n}&m_{n+1}\\\ \vdots&\vdots&&\vdots&\vdots\\\
m_{n-1}&m_{n}&\cdots&m_{2n-2}&m_{2n-1}\\\
I_{p}&I_{p}x&\cdots&I_{p}x^{n-1}&I_{p}x^{n}\end{bmatrix}.$

Given two sequences of matrix bi-orthonormal polynomials
$\\{P_{k}^{[1]}(x)\\}_{k=0}^{\infty}$ and
$\\{P_{k}^{[2]}(x)\\}_{k=0}^{\infty}$, with respect to
$\left\langle{\cdot,\cdot}\right\rangle$, we define the $n$-th
Christoffel–Darboux kernel matrix polynomial

(4) $\displaystyle
K_{n}(x,y):=\sum_{k=0}^{n}(P_{k}^{[2]}(y))^{\top}H_{k}^{-1}P^{[1]}_{k}(x).$

Named after [32, 34] see also [33].

######  Proposition 15 (ABC Theorem).

An Aitken–Berg–Collar type formula

$\displaystyle
K_{n}(x,y)=\big{[}I_{p},\dots,I_{p}x^{n}\big{]}(M_{n})^{-1}\begin{bmatrix}I_{p}\\\
\vdots\\\ I_{p}y^{n}\end{bmatrix},$

holds.

The scalar version was rediscovered and popularized by Berg [20], who found it
in a paper of Collar [31], who attributes it to his teacher, Aitken. As we are
inverting a Hankel block matrix we are dealing with a Hankel Bezoutian type
expression. This is connected with the following

######  Proposition 16 (Christoffel–Darboux formula).

The Christoffel–Darboux kernel satisfies

$(x-y)K_{n}(x,y)=(P^{[2]}_{n}(y))^{\top}(H_{n})^{-1}P^{[1]}_{n+1}(x)-(P^{[2]}_{n+1}(y))^{\top}(H_{n})^{-1}P^{[1]}_{n}(x).$

##  2\. Connection formulas for Darboux transformations of Christoffel type

Given a monic matrix polynomial $W(x)$ of degree $N$ we consider a new matrix
of measures of the form

$\displaystyle\operatorname{d}\mu(x)\mapsto\operatorname{d}\hat{\mu}(x):=W(x)\operatorname{d}\mu(x)$

with the corresponding perturbed sesquilinear form

$\displaystyle\left\langle{P(x),Q(x)}\right\rangle_{W}=\int
P(x)W(x)\operatorname{d}\mu(x)(Q(x))^{\top}.$

In the same way as above, the moment block matrix

$\displaystyle\hat{M}:=\int\chi(x)W(x)\operatorname{d}\mu(x)(\chi(x))^{\top}$

is introduced. Let us assume that the perturbed moment matrix admits a
Gaussian factorization

$\displaystyle\hat{M}=\hat{S}_{1}^{-1}\hat{H}(\hat{S}_{2})^{-\top},$

where $\hat{S}_{1},\hat{S}_{2}$ are lower unitriangular block matrices and
$\hat{H}$ is a diagonal block matrix

$\displaystyle\hat{S}_{i}$ $\displaystyle=\begin{bmatrix}I_{p}&0&0&\dots\\\
(\hat{S}_{i})_{1,0}&I_{p}&0&\cdots\\\
(\hat{S}_{i})_{2,0}&(\hat{S}_{i})_{2,1}&I_{p}&\ddots\\\
&&&\ddots\end{bmatrix},$ $\displaystyle\hat{H}$
$\displaystyle=\begin{bmatrix}\hat{H}_{0}&0&0&\cdots\\\
0&\hat{H}_{1}&0&\ddots\\\ 0&0&\hat{H}_{2}&\ddots\\\
\vdots&\vdots&\ddots&\ddots\end{bmatrix},$ $\displaystyle i$
$\displaystyle=1,2.$

Then, we have the corresponding perturbed bi-orthogonal matrix polynomials

$\displaystyle\hat{P}^{[i]}(x)$ $\displaystyle=\hat{S}_{i}\chi(x),$
$\displaystyle i$ $\displaystyle=1,2,$

with respect to the perturbed sesquilinear form
$\left\langle{\cdot,\cdot}\right\rangle_{W}$.

######  Remark 4.

The discussion for monic matrix polynomial perturbations and perturbations
with a matrix polynomial with non singular leading coefficients are
equivalent. Indeed, if instead of a monic matrix polynomial we have a matrix
polynomial $\tilde{W}(x)=A_{N}x^{N}+\dots+A_{0}$ with a nonsingular leading
coefficient, $\det A_{N}\neq 0$, then we could factor out $A_{N}$,
$\tilde{W}(x)=A_{N}W(x)$, where $W$ is monic. The moment matrices are related
by $\tilde{M}=A_{N}\hat{M}$ and, moreover,
$\tilde{S}_{1}=A_{N}\hat{S}_{1}(A_{N})^{-1}$, $\tilde{H}=A_{N}\hat{H}$,
$\tilde{S}_{2}=\hat{S}_{2}$, and
$\tilde{P}^{[1]}_{k}(x)=A_{N}P^{[1]}_{k}(x)(A_{N})^{-1}$ as well as
$\tilde{P}^{[2]}_{k}(x)=\hat{P}^{[2]}_{k}(x)$.

###  2.1. Connection formulas for bi-orthogonal polynomials

######  Proposition 17.

The moment matrix $M$ and the $W$ perturbed moment matrix $\hat{M}$ satisfy

$\displaystyle\hat{M}=W(\Lambda)M.$

######  Definition 11.

Let us introduce the following semi-infinite matrices

$\displaystyle\omega^{[1]}$ $\displaystyle:=\hat{S}_{1}W(\Lambda)S_{1}^{-1},$
$\displaystyle\omega^{[2]}$ $\displaystyle:=(S_{2}\hat{S}_{2}^{-1})^{\top},$

which we call resolvent or connection matrices.

######  Proposition 18 (Connection formulas).

Perturbed and non perturbed bi-orthogonal polynomials are subject to the
following linear connection formulas

(5) $\displaystyle\omega^{[1]}P^{[1]}(x)=$
$\displaystyle\hat{P}^{[1]}(x)W(x),$ (6) $\displaystyle P^{[2]}(x)=$
$\displaystyle(\omega^{[2]})^{\top}\hat{P}^{[2]}(x).$

######  Proposition 19.

The following relations hold

$\displaystyle\hat{H}\omega^{[2]}=\omega^{[1]}H.$

###### Proof.

From Proposition 17 and the $LU$ factorization we get

$\displaystyle\hat{S}_{1}^{-1}\hat{H}\hat{S}_{2}^{-\top}=W(\Lambda)S_{1}^{-1}HS_{2}^{-\top},$

so that

$\displaystyle\hat{H}(S_{2}\hat{S}_{2}^{-1})^{\top}=\hat{S}_{1}W(\Lambda)S_{1}^{-1}H$

and the result follows. ∎

From this result we easily get that

######  Proposition 20.

The resolvent matrix $\omega$ is a band upper triangular block matrix with all
the block superdiagonals above the $N$-th one equal to zero.

$\displaystyle\omega^{[1]}=\begin{bmatrix}\omega^{[1]}_{0,0}&\omega^{[1]}_{0,1}&\omega^{[1]}_{0,2}&\dots&\omega^{[1]}_{0,N-1}&I_{p}&0&0&\dots\\\
0&\omega^{[1]}_{1,1}&\omega^{[1]}_{1,2}&\dots&\omega^{[1]}_{1,N-1}&\omega^{[1]}_{1,N}&I_{p}&0&\dots\\\
0&0&\omega^{[1]}_{2,2}&\dots&\omega^{[1]}_{2,N-1}&\omega^{[1]}_{2,N}&\omega^{[1]}_{2,N+1}&I_{p}&\ddots\\\
&\ddots&\ddots&\ddots&&&&\ddots&\ddots\end{bmatrix}$

with

(7) $\displaystyle\omega^{[1]}_{k,k}=\hat{H}_{k}(H_{k})^{-1}.$

###  2.2. Connection formulas for the Christoffel–Darboux kernel

In order to relate the perturbed and non perturbed kernel matrix polynomials
let us introduce the following truncation of the connection matrix $\omega$.

######  Definition 12.

We introduce the lower unitriangular matrix
$\omega_{(n,N)}\in\mathbb{R}^{Np\times Np}$ with

$\displaystyle\omega_{(n,N)}:=\begin{cases}\begin{bmatrix}0&\dots&0&0&\dots&0\\\
\vdots&&&\vdots&\\\ 0&\dots&0&0&\dots&0\\\
\omega^{[1]}_{0,n+1}&\dots&\omega_{0,N-1}&I_{p}&\ddots&0\\\
\vdots&&&&\ddots&\vdots\\\
\omega^{[1]}_{n,n+1}&\cdots&&&\omega^{[1]}_{n,n+N-1}&I_{p}\end{bmatrix},&n<N,\\\\[4.0pt]
\begin{bmatrix}I_{p}&0&\dots&0&0\\\
\omega^{[1]}_{n-N+2,n+1}&I_{p}&\ddots&0&0\\\
\vdots&\ddots&\ddots&\ddots&\vdots\\\ \vdots&&\ddots&I_{p}&0\\\
\omega^{[1]}_{n,n+1}&\cdots&&\omega^{[1]}_{n,n+N-1}&I_{p}\end{bmatrix},&n\geq
N,\end{cases}$

and the diagonal block matrix

$\displaystyle\hat{H}_{n,N}=\operatorname{diag}(\hat{H}_{n-N+1},\dots,\hat{H}_{n}).$

Then, we can write the important

######  Theorem 1.

The perturbed and original Christoffel-Darboux kernels are related by the
following connection formulas

$\displaystyle
K_{n}(x,y)+\Big{[}(\hat{P}_{n-N+1}^{[2]}(x))^{\top},\dots,(\hat{P}_{n}^{[2]}(x))^{\top}\Big{]}(\hat{H}_{(n,N)})^{-1}\omega_{(n,N)}\begin{bmatrix}P^{[1]}_{n+1}(x)\\\
\vdots\\\ P^{[1]}_{n+N}(x)\end{bmatrix}=\hat{K}_{n}(x,y)W(x),$

by convention $\hat{P}_{j}^{[2]}=0$ whenever $j<0$.

###### Proof.

Consider the truncation

$\displaystyle(\omega^{[2]})_{[n+1]}:=\begin{bmatrix}I_{p}&\dots&\omega^{[2]}_{0,N-1}&\omega^{[2]}_{0,N}&0&0&\dots&0\\\
0&I_{p}&\dots&\omega^{[2]}_{1,N+1}&\omega^{[2]}_{1,N+2}&0&\dots&0\\\
\vdots&&\ddots&&&\ddots&\ddots\\\ &&&&&\\\
0&0&&0&I_{p}&\dots&&\omega^{[2]}_{n-N,n}\\\
\vdots&\vdots&&&\ddots&\ddots&&\vdots\\\ 0&0&&&&&I_{p}&\omega^{[2]}_{n-1,n}\\\
0&0&&&&&0&I_{p}\end{bmatrix}.$

Recalling (6) in the form
$(P^{[2]}(y))^{\top}=(\hat{P}^{[2]}(y))^{\top}\omega^{[2]}$ we see that
$\big{(}(\hat{P}^{[2]}(y))_{[n+1]}\big{)}^{\top}(\omega^{[2]})_{[n+1]}=\big{(}(P^{[2]}(y))_{[n+1]}\big{)}^{\top}$
holds for the $n$-th truncations of $P^{[2]}(y))$ and $\hat{P}^{[2]}(y)$.
Therefore,

$\displaystyle\big{(}(\hat{P}^{[2]}(y))_{[n+1]}\big{)}^{\top}(\omega^{[2]})_{[n+1]}(H_{[n+1]})^{-1}(P^{[1]}(x))_{[n+1]}=$
$\displaystyle\big{(}(P^{[2]}(y))_{[n+1]}\big{)}^{\top}(H_{[n+1]})^{-1}(P^{[1]}(x))_{[n+1]}$
$\displaystyle=$ $\displaystyle K_{n}(x,y).$

Now, we consider $(\omega^{[2]})_{[n+1]}(H_{[n+1]})^{-1}(P^{[1]}(x))_{[n+1]}$
and recall Proposition 19 in the form

$\displaystyle(\omega^{[2]})_{[n+1]}\big{(}H_{[n+1]}\big{)}^{-1}=(\hat{H}_{[n+1]})^{-1}(\omega^{[1]})_{[n+1]}$

which leads to

$\displaystyle(\omega^{[2]})_{n+1}(H_{n+1})^{-1}(P^{[1]}(x))_{n+1}=(\hat{H}_{[n+1]})^{-1}(\omega^{[1]})_{[n+1]}(P^{[1]}(x))_{[n+1]}.$

Observe also that

$\displaystyle(\omega^{[1]})_{[n+1]}(P^{[1]}(x))_{[n+1]}=(\omega^{[1]}P^{[1]}(x))_{[n+1]}-\begin{bmatrix}0_{(n-N)p\times
p}\\\ V_{N}(x)\end{bmatrix}$

with

$\displaystyle V_{N}(x)=\omega_{(n,N)}\begin{bmatrix}P^{[1]}_{n+1}(x)\\\
\vdots\\\ P^{[1]}_{n+N}(x)\end{bmatrix}.$

Hence, recalling (5) we get

$\displaystyle(\omega^{[1]})_{[n+1]}(P^{[1]}(x))_{[n+1]}=(\hat{P}^{[1]}(x))_{[n+1]}W(x)-\begin{bmatrix}0_{(n-N)p\times
p}\\\ V_{N}(x)\end{bmatrix},$

and consequently

$\big{(}(\hat{P}^{[2]}(y))_{[n+1]}\big{)}^{\top}(\omega^{[2]})_{[n+1]}(H_{[n+1]})^{-1}(P^{[1]}(x))_{[n+1]}\\\
\begin{aligned}
&=\big{(}(\hat{P}^{[2]}(y))_{[n+1]}\big{)}^{\top}(\hat{H}_{[n+1]})^{-1}(\hat{P}^{[1]}(x))_{[n+1]}W(x)-\big{(}(\hat{P}^{[2]}(y))_{[n+1]}\big{)}^{\top}(\hat{H}_{[n+1]})^{-1}\begin{bmatrix}0_{(n-N)p\times
p}\\\ V_{N}(x)\end{bmatrix}\\\
&=\hat{K}_{n+1}(x,y)W(y)-\big{(}(\hat{P}^{[2]}(y))_{[n+1]}\big{)}^{\top}(\hat{H}_{[n+1]})^{-1}\begin{bmatrix}0_{(n-N)p\times
p}\\\ V_{N}(x)\end{bmatrix}.\end{aligned}$

∎

##  3\. Monic matrix polynomial perturbations

In this section we study the case of perturbations by monic matrix polynomials
$W(x)$, which is equivalent to matrix polynomials with nonsingular leading
coefficients. Using the theory given in §1.3 we are able to extend the
celebrated Christoffel formula to this context.

###  3.1. The Christoffel formula for matrix bi-orthogonal polynomials

We are now ready to show how the perturbed set of matrix bi-ortogonal
polynomials $\\{\hat{P}_{n}^{[1]}(x),\hat{P}_{n}^{[2]}(x),\\}_{n=0}^{\infty}$
is related to the original set
$\\{P_{n}^{[1]}(x),P_{n}^{[2]}(x)\\}_{n=0}^{\infty}$

######  Proposition 21.

Let $v_{j}^{(i)}(x)$ be the adapted root polynomials of the monic matrix
polynomial $W(x)$ given in (3). Then, for each eigenvalue $x_{i}\in\sigma(W)$,
$i\in\\{1,\dots,q\\}$,

(8)
$\displaystyle\omega^{[1]}_{k,k}\frac{\operatorname{d}^{r}(P_{k}^{[1]}v_{j}^{(i)})}{\operatorname{d}x^{r}}\bigg{|}_{x=x_{i}}+\dots+\omega^{[1]}_{k,k+N-1}\frac{\operatorname{d}^{r}(P_{k+N-1}^{[1]}v_{j}^{(i)})}{\operatorname{d}x^{r}}\bigg{|}_{x=x_{i}}=-\frac{\operatorname{d}^{r}(P^{[1]}_{k+N}v_{j}^{(i)})}{\operatorname{d}x^{r}}\bigg{|}_{x=x_{i}},$

for $r=0,\dots,\kappa^{(i)}_{j}-1,$ and $j=1\dots,s_{i}$.

###### Proof.

From (5) we get

$\displaystyle\omega^{[1]}_{k,k}P_{k}^{[1]}(x)+\dots+\omega^{[1]}_{k,k+N-1}P_{k+N-1}^{[1]}(x)+P_{k+N}(x)=\hat{P}^{[1]}_{k}(x)W(x).$

Now, according to Proposition 8 we have

(9) $\displaystyle\begin{aligned}
\frac{\operatorname{d}^{r}}{\operatorname{d}x^{r}}\bigg{|}_{x=x_{i}}\big{(}\hat{P}^{[1]}_{k}Wv_{j}^{(i)}\big{)}&=\sum_{s=0}^{r}\binom{r}{s}\frac{\operatorname{d}^{r-s}\hat{P}^{[1]}_{k}}{\operatorname{d}x^{r-s}}\bigg{|}_{x=x_{i}}\frac{\operatorname{d}^{s}(Wv_{j}^{(i)})}{\operatorname{d}x^{s}}\bigg{|}_{x=x_{i}}\\\
&=0,\end{aligned}$

for $r=0,\dots,\kappa^{(i)}_{j}-1$ and $j=1\dots,s_{i}$. ∎

Recall that $\sum_{j=1}^{s_{i}}\kappa_{j}^{(i)}=\alpha_{i}$ and that the sum
of all multiplicities $\alpha_{i}$ is $Np=\sum_{i=1}^{q}\alpha_{i}$,
$q=\\#\sigma(W)$.

######  Definition 13.

  1. i)

For each eigenvalue $x_{i}\in\sigma(W)$, in terms of the adapted root
polynomials $v_{j}^{(i)}(x)$ of the monic matrix polynomial $W(x)$ given in
(3), we introduce the vectors

$\displaystyle\pi_{j,k}^{(r),(i)}:=\frac{\operatorname{d}^{r}(P_{k}^{[1]}v_{j}^{(i)})}{\operatorname{d}x^{r}}\bigg{|}_{x=x_{i}}\in\mathbb{C}^{p}$

and arrange them in the partial row matrices
$\pi_{k}^{(i)}\in\mathbb{C}^{p\times\alpha_{i}p}$ given by

$\displaystyle\pi_{k}^{(i)}=\Big{[}\pi^{(0),(i)}_{1,k},\dots,\pi^{(\kappa^{(i)}_{1}-1),(i)}_{1,k},\dots,\pi^{(0),(i)}_{s_{i},k},\dots,\pi^{(\kappa_{s_{i}}^{(i)}-1),(i)}_{s_{i},k}\Big{]}.$

We collect all them as

$\displaystyle\pi_{k}:=\big{[}\pi_{k}^{(1)},\dots,\pi_{k}^{(q)}\big{]}\in\mathbb{C}^{p\times
Np}.$

Finally, we have

$\displaystyle\Pi_{k,N}:=\begin{bmatrix}\pi_{k}\\\ \vdots\\\
\pi_{k+N-1}\end{bmatrix}\in\mathbb{C}^{Np\times Np}.$

  2. ii)

In a similar manner, we define

$\displaystyle\gamma_{j,n}^{(r),(i)}(y)$
$\displaystyle:=\frac{\operatorname{d}^{r}\big{(}K_{n}(x,y)v_{j}^{(i)}(x)\big{)}}{\operatorname{d}x^{r}}\bigg{|}_{x=x_{i}}\in\mathbb{C}^{p}[y],$
$\displaystyle\gamma_{n}^{(i)}(y)$
$\displaystyle:=\big{[}\gamma^{(0),(i)}_{1,n}(y),\dots,\gamma^{(\kappa^{(i)}_{1}-1),(i)}_{1,n}(y),\dots,\gamma^{(0),(i)}_{s_{i},n}(y),\dots,\gamma^{(\kappa_{s_{i}}^{(i)}-1),(i)}_{s_{i},n}(y)\big{]}\in\mathbb{C}^{p\times\alpha_{i}}[y],$

and, as above, collect all of them in

$\displaystyle\gamma_{n}(y)=\big{[}\gamma_{n}^{(1)},\dots,\gamma_{n}^{(q)}\big{]}\in\mathbb{C}^{p\times
Np}[y].$

######  Theorem 2 (The Christoffel formula for matrix bi-orthogonal
polynomials).

The perturbed set of matrix bi-orthogonal polynomials
$\\{\hat{P}_{k}^{[1]}(x),\hat{P}_{k}^{[2]}(x)\\}_{k=0}^{\infty}$, whenever
$\det\Pi_{k,N}\neq 0$, can be written as the following last quasideterminant

(12) $\displaystyle\hat{P}_{k}^{[1]}(x)W(x)$
$\displaystyle=\Theta_{*}\left[\begin{array}[]{c|c}\Pi_{k,N}&\begin{matrix}P^{[1]}_{k}(x)\\\
\vdots\\\ P^{[1]}_{k+N-1}(x)\end{matrix}\\\
\hline\cr\pi_{k+N}&P^{[1]}_{k+N}(x)\end{array}\right],$ (15)
$\displaystyle(\hat{P}_{k}^{[2]}(x))^{\top}(\hat{H}_{k})^{-1}$
$\displaystyle=\Theta_{*}\left[\begin{array}[]{c|c}\Pi_{k+1,N}&\begin{matrix}0\\\
\vdots\\\ 0\\\ I_{p}\end{matrix}\\\
\hline\cr\gamma_{k}(x)&0\end{array}\right].$

Moreover, the new matrix squared norms are

(18) $\displaystyle\hat{H}_{k}$
$\displaystyle=\Theta_{*}\left[\begin{array}[]{c|c}\Pi_{k,N}&\begin{matrix}H_{k}\\\
0\\\ \vdots\\\ 0\end{matrix}\\\ \hline\cr\pi_{k+N}&0\end{array}\right].$

###### Proof.

We assume that $P_{j}^{[2]}=0$ whenever $j<0$. To prove (12) notice that from
(8) one deduces for the rows of the connection matrix that

(19)
$\displaystyle[\omega^{[1]}_{k,k},\dots,\omega^{[1]}_{k,k+n-1}]=-\pi_{k+N}(\Pi_{k,N})^{-1}.$

Now, using (5) we get

$\displaystyle[\omega^{[1]}_{k,k},\dots,\omega^{[1]}_{k,k+N-1}]\begin{bmatrix}P^{[1]}_{k}(x)\\\
\vdots\\\
P^{[1]}_{k+N-1}(x)\end{bmatrix}+P^{[1]}_{k+N}(x)=\hat{P}^{[1]}_{k}(x)W(x).$

and (12) follows immediately.

To deduce (15) for $k\geq N$ notice that Theorem 1 together with (9) yields

$\displaystyle\gamma^{(r),(i)}_{j,k}(x)+\Big{[}(\hat{P}_{k-N+1}^{[2]}(x))^{\top},\dots,(\hat{P}_{k}^{[2]}(x))^{\top}\Big{]}(\hat{H}_{(k,N)})^{-1}\omega_{(k,N)}\begin{bmatrix}\pi^{(r),(i)}_{j,k+1}\\\
\vdots\\\ \pi^{(r),(i)}_{j,k+N}\end{bmatrix}=0$

for $r=0,\dots,\kappa^{(i)}_{j}-1$ and $j=1\dots,s_{i}$. We arrange these
equations in a matrix form to get

$\displaystyle\gamma_{k}(x)+\Big{[}(\hat{P}_{k-N+1}^{[2]}(x))^{\top},\dots,(\hat{P}_{k}^{[2]}(x))^{\top}\Big{]}(\hat{H}_{(k,N)})^{-1}\omega_{(k,N)}\Pi_{k+1,N}=0.$

Therefore, assuming that $\det\Pi_{k+1,N}\neq 0$, we get

$\displaystyle\Big{[}(\hat{P}_{k-N+1}^{[2]}(x))^{\top},\dots,(\hat{P}_{k}^{[2]}(x))^{\top}\Big{]}(\hat{H}_{(k,N)})^{-1}\omega_{(k,N)}=-\gamma_{k}(x)(\Pi_{k+1,N})^{-1},$

which, in particular, gives

$\displaystyle(\hat{P}_{k}^{[2]}(x))^{\top}(\hat{H}_{k})^{-1}=-\gamma_{k}(x)(\Pi_{k+1,N})^{-1}\begin{bmatrix}0\\\
\vdots\\\ 0\\\ I_{p}\end{bmatrix}.$

Finally, (18) is a consequence of (7) and (19). ∎

###  3.2. Degree one monic matrix polynomial perturbations

Let us illustrate the situation with the most simple case of a perturbation of
degree one monic polynomial matrix

$\displaystyle W(x)=I_{p}x-A.$

The spectrum $\sigma(I_{p}x-A)=\sigma(A)=\\{x_{1},\dots,x_{q}\\}$ is
determined by the zeroes of the characteristic polynomial of $A$

$\displaystyle\det(I_{p}x-A)=(x-x_{1})^{\alpha_{1}}\cdots(x-x_{q})^{\alpha_{q}},$

and for each eigenvalue let $s_{i}=\dim\operatorname{Ker}(I_{p}x_{i}-A)$ be
the corresponding geometric multiplicity, and $\kappa^{(i)}_{j}$,
$j=1,\dots,s_{i}$, its partial multiplicities, so that
$\alpha_{i}=\sum_{j=1}^{s_{i}}\kappa_{j}^{(i)}$ is the algebraic multiplicity
(the order of the eigenvalue as a zero of the characteristic polynomial of
$A$). After a similarity transformation of $A$ we will get its canonical
Jordan form. With no lack of generality we assume that $A$ is already given in
Jordan canonical form

$\displaystyle\left[\begin{array}[]{ccccccc}\leavevmode\hbox to6.67pt{\vbox
to6.67pt{\pgfpicture\makeatletter\hbox{\hskip 3.33301pt\lower-3.33301pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\mathcal{J}_{\kappa^{(1)}_{1}}(x_{1})&&0&&&&\\\
&\ddots&&&&0&\\\ 0&&\mathcal{J}_{\kappa^{(1)}_{s_{1}}}(x_{1})\leavevmode\hbox
to6.67pt{\vbox to6.67pt{\pgfpicture\makeatletter\hbox{\hskip
3.33301pt\lower-3.33301pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\leavevmode\hbox
to4.8pt{\vbox to17.8pt{\pgfpicture\makeatletter\hbox{\hskip
2.4pt\lower-8.40001pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{ { {}{}{}{}{}}{}{}{}{{}}{} {
{}{}{}{}{}}{}{}{}{{}}{}{}{}{}{{}}{}\pgfsys@beginscope\pgfsys@invoke{
}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{-2.0pt}{9.00002pt}\pgfsys@moveto{-2.0pt}{9.00002pt}\pgfsys@lineto{-2.0pt}{-8.00002pt}\pgfsys@lineto{2.0pt}{-8.00002pt}\pgfsys@lineto{2.0pt}{9.00002pt}\pgfsys@closepath\pgfsys@moveto{2.0pt}{-8.00002pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}&&&&\\\
&&&\ddots&&&\\\ &&&&\leavevmode\hbox to6.67pt{\vbox
to6.67pt{\pgfpicture\makeatletter\hbox{\hskip 3.33301pt\lower-3.33301pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\mathcal{J}_{\kappa^{(q)}_{1}}(x_{q})&&0\\\
&0&&&&\ddots&\\\
&&&&0&&\mathcal{J}_{\kappa^{(q)}_{s_{q}}}(x_{q})\leavevmode\hbox
to6.67pt{\vbox to6.67pt{\pgfpicture\makeatletter\hbox{\hskip
3.33301pt\lower-3.33301pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\leavevmode\hbox
to4.8pt{\vbox to17.8pt{\pgfpicture\makeatletter\hbox{\hskip
2.4pt\lower-8.40001pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{ { {}{}{}{}{}}{}{}{}{{}}{} {
{}{}{}{}{}}{}{}{}{{}}{}{}{}{}{{}}{}\pgfsys@beginscope\pgfsys@invoke{
}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{-2.0pt}{9.00002pt}\pgfsys@moveto{-2.0pt}{9.00002pt}\pgfsys@lineto{-2.0pt}{-8.00002pt}\pgfsys@lineto{2.0pt}{-8.00002pt}\pgfsys@lineto{2.0pt}{9.00002pt}\pgfsys@closepath\pgfsys@moveto{2.0pt}{-8.00002pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{array}\right],$

where the Jordan blocks corresponding to each eigenvalue are given by

$\displaystyle\mathcal{J}_{\kappa_{j}^{(i)}}(x_{i})$
$\displaystyle:=\begin{bmatrix}x_{i}&1&0&&\\\ 0&x_{i}&1&&\\\
&&\ddots&\ddots&\\\ &&&x_{i}&1\\\
&&&0&x_{i}\end{bmatrix}\in\mathbb{R}^{\kappa_{j}^{(i)}\times\kappa_{j}^{(i)}},$
$\displaystyle j$ $\displaystyle=1,\dots,s_{i}.$

For each eigenvalue $x_{i}$ we pick a basis
$\big{\\{}v^{(i)}_{0,j}\big{\\}}_{j=1}^{s_{i}}$ of
$\operatorname{Ker}(I_{p}x_{i}-A)$, then look for vectors
$\\{v_{r,j}^{(i)}\\}_{r=1}^{\kappa_{j}^{(i)}-1}$ such that

$\displaystyle(I_{p}x_{i}-A)v_{r,j}^{(i)}=$ $\displaystyle-v_{r-1,j}^{(i)},$
$\displaystyle r$ $\displaystyle=1,\dots,\kappa_{j}^{(i)}-1,$

so that $\\{v^{(i)}_{r,j}\\}_{r=0}^{\kappa^{(i)}_{j}}$ is a Jordan chain. As
we are dealing with $A$ in its canonical form the vectors $v^{(i)}_{r,j}$ can
be identified with those of the canonical basis $\\{e_{i}\\}_{i=1}^{p}$ of
$\mathbb{R}^{p}$ with $e_{i}=(0,\cdots 1_{i},\cdots 0)_{p}^{\top}$. Indeed, we
have

$\displaystyle
v^{(i)}_{r,j}=e_{\alpha_{1}+\dots+\alpha_{i-1}+\kappa^{(i)}_{1}+\dots+\kappa^{(i)}_{j-1}+r+1}.$

Then, consider the polynomial vectors

$\displaystyle v^{(i)}_{j}(x)$
$\displaystyle=\sum_{r=0}^{\kappa^{(i)}_{j}-1}v^{(i)}_{r,j}(x-x_{i})^{r},$
$\displaystyle j$ $\displaystyle=1,\dots,s_{i},$ $\displaystyle i$
$\displaystyle=1\dots,q,$

which satisfy

$\displaystyle\frac{\operatorname{d}^{r}}{\operatorname{d}x^{r}}\bigg{|}_{x=x_{i}}\Big{(}(I_{p}x-A)v^{(i)}_{j}\Big{)}$
$\displaystyle=0,$ $\displaystyle r$
$\displaystyle=0,\dots,\kappa^{(i)}_{j}-1,$ $\displaystyle j$
$\displaystyle=1,\dots,s_{i},$ $\displaystyle i$ $\displaystyle=1,\dots,q.$

Consequently,

(20)
$\displaystyle\frac{\operatorname{d}^{r}}{\operatorname{d}x^{r}}\bigg{|}_{x=x_{i}}\Big{(}\hat{P}^{[1]}(x)(I_{p}x-A)v^{(i)}_{j}\Big{)}$
$\displaystyle=0,$ $\displaystyle r$
$\displaystyle=0,\dots,\kappa^{(i)}_{j}-1,$ $\displaystyle j$
$\displaystyle=1,\dots,s_{i},$ $\displaystyle i$ $\displaystyle=1,\dots,q.$

Now, let us notice the following simple fact

######  Proposition 22.

For a given matrix $A\in\mathbb{R}^{p\times p}$ any matrix polynomial
$P(x)=\sum_{k=0}^{n}P_{k}x^{k}\in\mathbb{R}^{p\times p}[x]$, $\deg P=n$, can
be written as

$\displaystyle P$ $\displaystyle=\sum_{k=0}^{n}P_{k}^{(A)}(I_{p}x-A)^{k}.$

In particular, we have $P_{0}^{(A)}=P(A):=\sum_{k=0}^{n}P_{k}A^{k}$.

######  Proposition 23.

We can write

$\displaystyle\frac{1}{r!}\frac{\operatorname{d}^{r}(Pv^{(i)}_{j})}{\operatorname{d}x^{r}}\bigg{|}_{x=x_{i}}$
$\displaystyle=P(A)v^{(i)}_{r,j},$ $\displaystyle r$
$\displaystyle=0,\dots,\kappa^{(i)}_{j}-1,$ $\displaystyle j$
$\displaystyle=1,\dots,s_{i},$ $\displaystyle i$ $\displaystyle=1,\dots,q.$

###### Proof.

Observe that

$\displaystyle\frac{\operatorname{d}^{s}P}{\operatorname{d}x^{s}}\bigg{|}_{x=x_{i}}v^{(i)}_{r-s,j}$
$\displaystyle=\sum_{l=r}^{k}\frac{l!}{(l-s)!}P^{(A)}_{l}(I_{p}x_{i}-A)^{l-s}v^{(i)}_{r-s,j}$
$\displaystyle=\sum_{l=s}^{r}(-1)^{l-s}\frac{l!}{(l-s)!}P^{(A)}_{l}v^{(i)}_{r-l,j},$
$\displaystyle s$ $\displaystyle=1,\dots,r,$ $\displaystyle j$
$\displaystyle=1,\dots,s_{i},$ $\displaystyle i=1,\dots,q,$

and, consequently,

$\displaystyle\frac{\operatorname{d}^{r}(Pv^{(i)}_{j})}{\operatorname{d}x^{r}}\bigg{|}_{x=x_{i}}$
$\displaystyle=\sum_{s=0}^{r}\frac{r!}{s!}\sum_{l=s}^{r}(-1)^{l-s}\frac{l!}{(l-s)!}P^{(A)}_{l}v^{(i)}_{r-l,j}$
$\displaystyle=r!\sum_{s=0}^{r}\sum_{l=s}^{r}(-1)^{l-s}\binom{l}{s}P^{(A)}_{l}v^{(i)}_{r-l,j}$
$\displaystyle=r!\sum_{m=0}^{r}(-1)^{r-m}\Big{[}\sum_{s=0}^{r-m}(-1)^{s}\binom{r-m}{s}\Big{]}P^{(A)}_{r-m}v^{(i)}_{m,j}$
$\displaystyle=r!P(A)v^{(i)}_{r,j},$

for $r=0,\dots,\kappa^{(i)}_{j}-1$, $j=1,\dots,s_{i}$ and $i=1,\dots,q$. ∎

In the following we’ll use

$\displaystyle
K_{n}(y,A):=\sum_{m=0}^{n}P^{[2]}_{m}(y)H_{m}^{-1}P^{[1]}_{m}(A).$

######  Proposition 24 (Degree one Christoffel formula).

If $W(x)=I_{p}x-A$ and $\det P_{n}^{[1]}(A)\neq 0$ for $n\in\mathbb{Z}_{+}$,
then the Christoffel formulas can be written as

$\displaystyle\hat{P}^{[1]}_{n}(x)(I_{p}x-A)$
$\displaystyle=\Theta_{*}\begin{bmatrix}P_{n}^{[1]}(A)&P_{n}^{[1]}(x)\\\
P_{n+1}^{[1]}(A)&P_{n+1}^{[1]}(x)\end{bmatrix}$
$\displaystyle(\hat{P}_{n}^{[2]}(y))^{\top}$
$\displaystyle=\Theta_{*}\begin{bmatrix}P_{n+1}^{[1]}(A)&I\\\
K_{n}(y,A)&0\end{bmatrix}$
$\displaystyle=P_{n+1}^{[1]}(x)-P_{n+1}^{[1]}(A)\big{[}P_{n}^{[1]}(A)\big{]}^{-1}P_{n}^{[1]}(x),$
$\displaystyle=-K_{n}(y,A)\big{[}P_{n}^{[1]}(A)\big{]}^{-1}.$

For the perturbed matrix squared norms we have

$\displaystyle\hat{H}_{n}$
$\displaystyle=\Theta_{*}\begin{bmatrix}P^{[1]}(A)&H_{n}\\\
P^{[n+1]}(A)&0\end{bmatrix}$
$\displaystyle=-P_{n+1}^{[1]}(A)\big{[}P_{n}^{[1]}(A)\big{]}^{-1}H_{n}.$

###### Proof.

According to (5) and Theorem 1

$\displaystyle\omega_{n,n}P_{n}^{[1]}(x)+P_{n+1}^{[1]}(x)$
$\displaystyle=\hat{P}_{n}^{[1]}(x)(I_{p}x-A),$ $\displaystyle
K_{n}(x,y)+(\hat{P}_{n}^{[2]}(y))^{\top}P^{[1]}_{n+1}(x)$
$\displaystyle=\hat{K}_{n}(x,y)(I_{p}x-A),$

and using (20) we conclude

$\displaystyle\omega_{n,n}\frac{\operatorname{d}^{r}(P_{k}^{[1]}v^{(i)}_{j})}{\operatorname{d}x^{r}}\bigg{|}_{x=x_{i}}+\frac{\operatorname{d}^{r}(P_{n+1}^{[1]}v^{(i)}_{j})}{\operatorname{d}x^{r}}\bigg{|}_{x=x_{i}}$
$\displaystyle=0,$
$\displaystyle\frac{\operatorname{d}^{r}}{\operatorname{d}x^{r}}\bigg{|}_{x=x_{i}}\big{(}(K_{n}(x,y)v^{(i)}_{j}(x)\big{)}+(\hat{P}_{n}^{[2]}(y))^{\top}\frac{\operatorname{d}^{r}(P_{n+1}^{[1]}v^{(i)}_{j})}{\operatorname{d}x^{r}}\bigg{|}_{x=x_{i}}$
$\displaystyle=0,$

for $r=0,\dots,\kappa^{(i)}_{j}-1$, $j=1,\dots,s_{i}$ and $i=1,\dots,q$. From
Proposition 23 and the fact that the ordered arrangement of the Jordan chain
vectors $v^{(i)}_{r,j}$ gives the identity matrix $I_{p}$ we conclude

$\displaystyle\omega_{n,n}P^{[1]}_{k}(A)+P^{[1]}_{n+1}(A)$ $\displaystyle=0,$
$\displaystyle K_{n}(y,A)+(\hat{P}_{n}^{[2]}(y))^{\top}P^{[1]}_{n+1}(A)$
$\displaystyle=0.$

∎

We now illustrate the Christoffel formula in the matrix orthogonal polynomial
context with a simple case. We will study what is the effect of the
Christoffel transformation on a positive Borel scalar measure
$\operatorname{d}\mu(x)$, thus the perturbed matrix of measures is
$(I_{2}x-A)\operatorname{d}\mu(x)$. The perturbed monic orthogonal polynomials
will be expressed, see Proposition 24,

(21) $\displaystyle\hat{P}^{[1]}_{n}(x)$
$\displaystyle=\Big{(}I_{2}p_{n+1}(x)-p_{n+1}(A)\big{(}p_{n}(A)\big{)}^{-1}p_{n}(x)\Big{)}(I_{2}x-A)^{-1},$
(22) $\displaystyle\hat{P}^{[2]}_{n}(x)$
$\displaystyle=-K_{n}(x,A)\big{[}P_{n}^{[1]}(A)\big{]}^{-1}$

where $p_{n}(x),K_{n}(x,y)$ are the scalar orthogonal polynomials and kernel
polynomials associated with the original scalar positive Borel measure
$\operatorname{d}\mu(x)$. Observe that despite starting with a set of
orthogonal polynomials the perturbation generates a set of bi-orthogonal
matrix polynomials. As the original measure is scalar, if we ensure that
$A=A^{\top}$ is symmetric, we will get
$\hat{P}_{n}(x):=P^{[1]}_{n}(x)=P^{[2]}_{n}(x)$, a new set of orthogonal
matrix polynomials.

However, this will be a very trivial situation as we have

######  Proposition 25.

The matrix orthogonal polynomials $\\{\hat{P}_{n}(x)\\}_{n=0}^{\infty}$ of the
matrix of measures $(I_{p}x-A)\operatorname{d}\mu(x)$, where $A=A^{\top}$ is
symmetric and $\operatorname{d}\mu$ is a positive Borel scalar measure are
similar to diagonal matrix orthogonal polynomials.

###### Proof.

Being the matrix $A$ symmetric it will be always diagonalizable

$\displaystyle A=QDQ^{\top},$

where $Q$ is an orthogonal matrix $Q^{\top}=Q^{-1}$ and
$D=\operatorname{diag}(x_{1},\dots,x_{p})$, is a diagonal matrix that collects
the eigenvalues, not necessarily different, of $A$.

At the end, the new orthogonal polynomials will be

(23) $\displaystyle\hat{P}_{n}(x)$
$\displaystyle=Q\begin{bmatrix}\frac{p_{n+1}(x)-\frac{p_{n+1}(x_{1})}{p_{n}(x_{1})}p_{n}(x)}{x-x_{1}}&0&\dots&0\\\
0&\frac{p_{n+1}(x)-\frac{p_{n+1}(x_{2})}{p_{n}(x_{2})}p_{n}(x)}{x-x_{2}}&&0\\\
\vdots&\ddots&\ddots&\vdots\\\ \\\
0&0&&\frac{p_{n+1}(x)-\frac{p_{n+1}(x_{p})}{p_{n}(x_{p})}p_{n}(x)}{x-x_{p}}\end{bmatrix}Q^{\top},$

and the result is proven. ∎

Thus, we have a diagonal bunch of elementary Darboux transformations of the
original scalar orthogonal polynomials associated to the scalar measure
$\operatorname{d}\mu$. This situation reappears even when the matrix is not
symmetric but diagonalizable, as we will have again that the perturbed matrix
orthogonal polynomials will be similar to a similar bunch of elementary
Darboux transformations of the original scalar orthogonal polynomials.

$\displaystyle\hat{P}^{[1]}_{n}(x)$
$\displaystyle=Q\begin{bmatrix}\frac{p_{n+1}(x)-\frac{p_{n+1}(x_{1})}{p_{n}(x_{1})}p_{n}(x)}{x-x_{1}}&0&\dots&0\\\
0&\frac{p_{n+1}(x)-\frac{p_{n+1}(x_{2})}{p_{n}(x_{2})}p_{n}(x)}{x-x_{2}}&&0\\\
\vdots&&\ddots&\vdots\\\ \\\
0&0&&\frac{p_{n+1}(x)-\frac{p_{n+1}(x_{p})}{p_{n}(x_{p})}p_{n}(x)}{x-x_{p}}\end{bmatrix}Q^{-1},$
$\displaystyle\hat{P}^{[2]}_{n}(x)$
$\displaystyle=-Q\begin{bmatrix}\frac{K_{n}(x,x_{1})}{p_{n}(x_{1})}&0&\dots&0\\\
0&\frac{K_{n}(x,x_{2})}{p_{n}(x_{2})}&&0\\\ \vdots&&\ddots&\vdots\\\ \\\
0&0&&\frac{K_{n}(x,x_{p})}{p_{n}(x_{p})}\end{bmatrix}Q^{-1},$

where $Q$ does not need to be an orthogonal matrix.

If the matrix is not diagonalizable and has nontrivial Jordan blocks the
situation is different. Let us explore this case when $p=2$. Hence, we
consider

$\displaystyle W(x)$ $\displaystyle=I_{2}x-A,$

with

$\displaystyle A$ $\displaystyle=M\begin{bmatrix}x_{1}&1\\\
0&x_{1}\end{bmatrix}M^{-1}.$

Now we have only one eigenvalue $\sigma(A)=\\{x_{1}\\}$, with a length 2
Jordan chain. Thus, there is a linear basis
$\\{v_{1},v_{2}\\}\subset\mathbb{R}^{2}$, $(A-x_{1})v_{1}=v_{2}$,
$(A-x_{1})v_{2}=0$, with $v_{i}=[v_{i,1}\quad v_{i,2}]^{\top}$,
$i\in\\{1,2\\}$, such that

$\displaystyle M=\begin{bmatrix}v_{1,1}&v_{2,1}\\\
v_{1,2}&v_{2,2}\end{bmatrix}.$

Therefore,

$\displaystyle p_{n+1}(A)\big{(}p_{n}(A)\big{)}^{-1}$
$\displaystyle=M\begin{bmatrix}\frac{p_{n+1}(x_{1})}{p_{n}(x_{1})}&\frac{W(p_{n},p_{n+1})(x_{1})}{(p_{n}(x_{1}))^{2}}\\\
0&\frac{p_{n+1}(x_{1})}{p_{n}(x_{1})}\end{bmatrix}M^{-1},$
$\displaystyle(I_{2}x-A)^{-1}$
$\displaystyle=M\begin{bmatrix}\frac{1}{x-x_{1}}&\frac{1}{(x-x_{1})^{2}}\\\
0&\frac{1}{x-x_{1}}\end{bmatrix}M^{-1},$

where

$\displaystyle
W(p_{n},p_{n+1})(x)=p_{n}(x)p^{\prime}_{n+1}-p_{n+1}(x)p^{\prime}_{n}(x)$

is the Wronskian of two consecutive orthogonal polynomials.

Hence,

(24) $\displaystyle\hat{P}^{[1]}_{n}(x)$
$\displaystyle=M\begin{bmatrix}\frac{p_{n+1}(x)-\frac{p_{n+1}(x_{1})}{p_{n}(x_{1})}p_{n}(x)}{x-x_{1}}&\frac{p_{n+1}(x)-\frac{p_{n+1}(x_{1})}{p_{n}(x_{1})}p_{n}(x)-\frac{W(p_{n},p_{n+1})(x_{1})}{(p_{n}(x_{1}))^{2}}(x-x_{1})p_{n}(x)}{(x-x_{1})^{2}}\\\
0&\frac{p_{n+1}(x)-\frac{p_{n+1}(x_{1})}{p_{n}(x_{1})}p_{n}(x)}{x-x_{1}}\end{bmatrix}M^{-1},$
$\displaystyle\hat{P}^{[2]}(x)$
$\displaystyle=-M\begin{bmatrix}\frac{K_{n}(x,x_{1})}{p_{n}(x_{1})}&-\frac{K_{n}(x,x_{1})}{p_{n}(x_{1})}p_{n}^{\prime}(x_{1})+\frac{1}{p_{n}(x_{1})}\frac{\partial
K_{n}(x,y)}{\partial y}\big{|}_{y=x_{1}}\\\
0&\frac{K_{n}(x,x_{1})}{p_{n}(x_{1})}\end{bmatrix}M^{-1}$

Observe that the polynomials

$\displaystyle p_{n+1}(x)-\frac{p_{n+1}(x_{1})}{p_{n}(x_{1})}p_{n}(x),$
$\displaystyle
p_{n+1}(x)-\frac{p_{n+1}(x_{1})}{p_{n}(x_{1})}p_{n}(x)-\frac{W(p_{n},p_{n+1})(x_{1})}{(p_{n}(x_{1}))^{2}}(x-x_{1})p_{n}(x)$

have a zero at $x=x_{1}$ of order $1$ and $2$, respectively.

###  3.3. Examples

######  Example 1.

In [59] the authors define the notion of a classical pair $\\{w(x),D\\}$,
where $w(x)$ is a symmetric matrix valued weight function and $D$ is a second
order linear ordinary differential operator. In that paper a weight function
is said to be classical if there exists a second order linear ordinary
differential operator $D$ with matrix valued polynomial coefficients
$A_{j}(t)$, $\deg A_{j}\leq j$, of the form
$D=A_{2}(x)\frac{\operatorname{d}^{2}}{\operatorname{d}x^{2}}+A_{1}(x)\frac{\operatorname{d}}{\operatorname{d}x}+A_{0}(x)$,
such that $\left\langle{DP,Q}\right\rangle=\left\langle{P,DQ}\right\rangle$
for all matrix valued polynomial functions $P(x)$ and $Q(x)$. Then, the pair
$\\{w,D\\}$ is called a classical pair. In example 5.1 in [59] they present a
family of Jacobi type classical pairs that contains, up to equivalence, all
classical pairs of size two where $w(x)=x^{\alpha}(1-x)^{\beta}F(x)$, with
$\alpha,\beta>-1$ and $0<x<1$, and such that $F(x)$ is of degree one and which
are irreducible (in the sense that they are not equivalent to a direct sum of
classical pairs of size one). As we will show they are a direct sum of
orthogonal polynomials of size 1 produced by two degree one Christoffel
transformations of the scalar Jacobi polynomials with zeros at $x=0,1$. Thus,
we are faced with two scalar monic Jacobi polynomials with each of the two
parameters $\alpha$ and $\beta$ shifted by one, respectively. In [96] an
analysis of the reducibility of matrix weights is given. In particular, in
Example 2.4 they consider the case $\alpha=\beta$. We must stress that, as was
pointed in [59], reducibility of the matrix of weights $w(x)$ do not imply the
reducibility of the classical pair $\\{w(x),D\\}$. Indeed, despite that the
matrix of weights in this example is reducible the corresponding second order
linear differential operator is not.

The classical pair $\\{w(x)=x^{\alpha}(1-x)^{\beta}F(x),D\\}$ is given by

$\displaystyle F(x)$ $\displaystyle=F_{1}x+F_{0},$ $\displaystyle F_{1}$
$\displaystyle=\begin{bmatrix}0&-a\\\ \
-a&\frac{\beta-\alpha}{\alpha+1}a\end{bmatrix},$ $\displaystyle F_{0}$
$\displaystyle=\begin{bmatrix}1&1\\\ 1&1\end{bmatrix},$ $\displaystyle\text{
with }a$ $\displaystyle=\frac{\alpha+\beta+2}{\alpha+1},$

and a second order matrix linear ordinary differential operator

$\displaystyle
D=x(1-x)\frac{\operatorname{d}^{2}}{\operatorname{d}x^{2}}+(X-xU)\frac{\operatorname{d}}{\operatorname{d}x}+V$

where $U,V,X$ are constant matrices depending on a parameter $u$. The sequence
of orthogonal polynomials
$\big{\\{}\tilde{\hat{P}}^{(\alpha,\beta)}_{n}(x)\big{\\}}_{n=0}^{\infty}$
associated with the classical pair is not given in [59]. Here an explicit
representation of $\tilde{\hat{P}}^{(\alpha,\beta)}_{n}(x)$ using Darboux
transformations is deduced. In order to do it we consider that we have an
initial alternative Jacobi measure
$\operatorname{d}\mu(x)=x^{\alpha}(1-x)^{\beta}I_{2}$, with $\alpha,\beta>-1$
and $0<x<1$, which is perturbed by a degree one matrix polynomial $F$. This
matrix polynomial is not monic but its leading coefficient is non singular and
we can write

$\displaystyle F(x)$ $\displaystyle=F_{1}W(x),$ $\displaystyle W(x)$
$\displaystyle=I_{2}x-A,$ $\displaystyle A$
$\displaystyle:=-F_{1}^{-1}F_{0}=\frac{1}{a}\begin{bmatrix}\frac{\beta+1}{\alpha+1}&\frac{\beta+1}{\alpha+1}\\\
1&1\end{bmatrix},$

in terms of a degree one monic matrix polynomial $W(x)$. We have that $A$ has
two different eigenvalues $\sigma(A)=\\{0,1\\}$ with corresponding
eigenvectors $[1,-1]^{\top}$ and
$\Big{[}\frac{\beta+1}{\alpha+1},1\Big{]}^{\top}$, the matrix
$M:=\begin{bmatrix}1&\frac{\beta+1}{\alpha+1}\\\ -1&1\end{bmatrix}$ allows to
write $A=M\operatorname{diag}(0,1)M^{-1}$.

Remember, as was noticed in Remark 4, that from the monic orthogonal
polynomials $\hat{P}_{n}^{(\alpha,\beta),[1]}(x)$ with respect to $W$, we get

$\displaystyle\tilde{\hat{P}}_{n}^{(\alpha,\beta)}(x)=F_{1}\hat{P}_{n}^{(\alpha,\beta),[1]}(x)F_{1}^{-1},$

which are the monic orthogonal polynomials with respect to $w(x)$. As the
matrix of measures $F(x)\operatorname{d}\mu(x)$ is symmetric, the bi-
orthogonality collapses to orthogonality and the super-indexes $[1]$ and $[2]$
can be omitted. We will do the same with
$\hat{P}_{n}^{(\alpha,\beta),[1]}=\hat{P}_{n}^{(\alpha,\beta)}$.

Following [28, 29] we conclude that the set of monic matrix orthogonal
polynomials $\\{P_{n}^{(\alpha,\beta)}(x)\\}_{n=0}^{\infty}$ with respect to
$\operatorname{d}\mu(x)$ are
$P_{n}^{(\alpha,\beta)}(x)=p_{n}^{(\alpha,\beta)}(x)I_{2}$,111We must be
careful at this point with the notation. This is not the scalar standard
Jacobi polynomial usually denoted by the same symbol. In fact, if
$\mathcal{P}^{(\alpha,\beta)}(z)$ denotes the standard Jacobi polynomials,
then
$p^{(\alpha,\beta)}(x)=\frac{2^{n}}{S_{n}(\alpha,\beta)}\mathcal{P}^{(\beta,\alpha)}(2x-1)$,
notice the interchange between the parameters $\alpha\rightleftharpoons\beta$
and the linear transformation of the independent variable $x$. with the
alternative Jacobi polynomials $p_{n}^{(\alpha,\beta)}(x)$ given by

$\displaystyle p_{n}^{(\alpha,\beta)}(x)$
$\displaystyle=\frac{1}{S_{n}(\alpha,\beta)}\sum_{k=0}^{n}\binom{n+\alpha}{n-k}\binom{n+\beta}{k}x^{n-k}(x-1)^{k},$
with $\displaystyle S_{n}(\alpha,\beta)$
$\displaystyle=\binom{2n+\beta+\alpha}{n}.$

We easily see that

$\displaystyle p^{(\alpha,\beta)}_{n}(0)$
$\displaystyle=\frac{1}{S_{n}(\alpha,\beta)}\binom{n+\alpha}{n},$
$\displaystyle p^{(\alpha,\beta)}_{n}(1)$
$\displaystyle=\frac{1}{S_{n}(\alpha,\beta)}\binom{n+\beta}{n},$

so that

$\displaystyle\frac{p^{(\alpha,\beta)}_{n+1}(0)}{p^{(\alpha,\beta)}_{n}(0)}$
$\displaystyle=(n+1+\alpha)\rho^{(\alpha,\beta)}_{n},$
$\displaystyle\frac{p^{(\alpha,\beta)}_{n+1}(1)}{p^{(\alpha,\beta)}_{n}(1)}$
$\displaystyle=(n+1+\beta)\rho^{(\alpha,\beta)}_{n},$

where

$\displaystyle\rho^{(\alpha,\beta)}_{n}$
$\displaystyle:=\frac{(n+\beta+\alpha+1)}{(2n+\beta+\alpha+2)(2n+\beta+\alpha+1)}.$

From (23) we conclude

$\displaystyle\hat{P}_{n}^{(\alpha,\beta)}(x)$
$\displaystyle=M\begin{bmatrix}\frac{p_{n+1}^{(\alpha,\beta)}(x)-(n+1+\alpha)\rho^{(\alpha,\beta)}_{n}p_{n}^{(\alpha,\beta)}(x)}{x}&0\\\
0&\frac{p_{n+1}^{(\alpha,\beta)}(x)-(n+1+\beta)\rho^{(\alpha,\beta)}_{n}p_{n}^{(\alpha,\beta)}(x)}{x-1}\end{bmatrix}M^{-1}.$

However, we must notice that these two Darboux transformations correspond to
the following transformations of the Jacobi measure

$\displaystyle x^{\alpha}(x-1)^{\beta}\mapsto x(x^{\alpha}(x-1)^{\beta})$
$\displaystyle=x^{\alpha+1}(x-1)^{\beta},$ $\displaystyle
x^{\alpha}(x-1)^{\beta}\mapsto(x-1)(x^{\alpha}(x-1)^{\beta})=x^{\alpha}(x-1)^{\beta+1},$

i.e. the transformations correspond to the shifts $\alpha\mapsto\alpha+1$ and
$\beta\mapsto\beta+1$, respectively. Consequently,

$\displaystyle\hat{P}_{n}^{(\alpha,\beta)}(x)$
$\displaystyle=M\begin{bmatrix}p_{n}^{(\alpha+1,\beta)}(x)&0\\\
0&p_{n}^{(\alpha,\beta+1)}(x)\end{bmatrix}M^{-1}.$

With the matrix

$\displaystyle\tilde{M}:=\begin{bmatrix}-1&1\\\
\frac{1+\beta}{1+\alpha}&1\end{bmatrix}$

we can write $F_{1}M=-a\tilde{M}$. We finally get the monic matrix orthogonal
polynomials

$\displaystyle\tilde{\hat{P}}_{n}^{(\alpha,\beta)}(x)=\tilde{M}\begin{bmatrix}p_{n}^{(\alpha+1,\beta)}(x)&0\\\
0&p_{n}^{(\alpha,\beta+1)}(x)\end{bmatrix}\tilde{M}^{-1},$

for the matrix of measures $\tilde{W}(x)\operatorname{d}\mu(x)$ in example 5.1
of [59] which are explicitly expressed in terms of scalar Jacobi polynomials
as follows

$\displaystyle\tilde{\hat{P}}_{n}^{(\alpha,\beta)}(x)$
$\displaystyle=\frac{1}{2+\alpha+\beta}\begin{bmatrix}(\alpha+1)p_{n}^{(\alpha+1,\beta)}(x)+(\beta+1)p_{n}^{(\alpha,\beta+1)}(x)&-(\alpha+1)(p_{n}^{(\alpha+1,\beta)}(x)-p_{n}^{(\alpha,\beta+1)}(x))\\\
-(\beta+1)(p_{n}^{(\alpha+1,\beta)}(x)-p_{n}^{(\alpha,\beta+1)}(x))&(\beta+1)p_{n}^{(\alpha+1,\beta)}(x)+(\alpha+1)p_{n}^{(\alpha,\beta+1)}(x)\end{bmatrix}.$

To conclude with this example let us mention that in [27] it was found that
these matrix orthogonal polynomials also obey a first order ordinary
differential equation. From our point of view this is just a consequence of a
remarkable fact regarding the Darboux transformations
$p^{(\alpha+1,\beta)}(x),p^{(\alpha,\beta+1)}(x)$ of the original alternative
Jacobi polynomials. Under the hypergeometric function description of the
Jacobi polynomials one gets recurrences for the Jacobi polynomials. In
particular, from the Gauss’ contiguous relations one gets the first order
differential relations

$\displaystyle\Big{(}x\frac{\operatorname{d}}{\operatorname{d}x}+\alpha+1\Big{)}p_{n}^{(\alpha+1,\beta)}(x)$
$\displaystyle=(\alpha+1+n)p_{n}^{(\alpha,\beta+1)}(x),$
$\displaystyle\Big{(}(x-1)\frac{\operatorname{d}}{\operatorname{d}x}+\beta+1\Big{)}p_{n}^{(\alpha,\beta+1)}(x)$
$\displaystyle=(\beta+1+n)p_{n}^{(\alpha+1,\beta)}(x).$

This first order linear ordinary differential system can be re-casted as a
matrix linear differential equation as follows

(25) $\bigg{(}\begin{bmatrix}0&x-1\\\
x&0\end{bmatrix}\frac{\operatorname{d}}{\operatorname{d}x}+\begin{bmatrix}a_{1}&\beta+1\\\
\alpha+1&a_{2}\end{bmatrix}\bigg{)}\begin{bmatrix}p_{n}^{(\alpha+1,\beta)}(x)&0\\\
0&p_{n}^{(\alpha,\beta+1)}(x)\end{bmatrix}\\\
=\begin{bmatrix}p_{n}^{(\alpha+1,\beta)}(x)&0\\\
0&p_{n}^{(\alpha,\beta+1)}(x)\end{bmatrix}\begin{bmatrix}a_{1}&n+\beta+1\\\
n+\alpha+1&a_{2}\end{bmatrix},$

where $a_{1},a_{2}\in\mathbb{R}$. This equation is invariant under
multiplication on the right and on the left hand side by arbitrary diagonal
matrices

$\begin{bmatrix}l_{1}&0\\\ 0&l_{2}\end{bmatrix}\bigg{(}\begin{bmatrix}0&x-1\\\
x&0\end{bmatrix}\frac{\operatorname{d}}{\operatorname{d}x}+\begin{bmatrix}a_{1}&\beta+1\\\
\alpha+1&a_{2}\end{bmatrix}\bigg{)}\begin{bmatrix}r_{1}&0\\\
0&r_{2}\end{bmatrix}\begin{bmatrix}p_{n}^{(\alpha+1,\beta)}(x)&0\\\
0&p_{n}^{(\alpha,\beta+1)}(x)\end{bmatrix}\\\
=\begin{bmatrix}p_{n}^{(\alpha+1,\beta)}(x)&0\\\
0&p_{n}^{(\alpha,\beta+1)}(x)\end{bmatrix}\begin{bmatrix}l_{1}&0\\\
0&l_{2}\end{bmatrix}\begin{bmatrix}a_{1}&n+\beta+1\\\
n+\alpha+1&a_{2}\end{bmatrix}\begin{bmatrix}r_{1}&0\\\ 0&r_{2}\end{bmatrix}.$

After the similarity transformation,
$A\mapsto\tilde{M}^{-\top}A\tilde{M}^{\top}$ we find out that the orthogonal
polynomial $\tilde{\hat{P}}_{n}$ satisfies

$\displaystyle\Big{(}A_{1}(x)\frac{\operatorname{d}}{\operatorname{d}x}+A_{0}\Big{)}\big{(}\tilde{\hat{P}}_{n}\big{)}^{\top}=\big{(}\tilde{\hat{P}}_{n}\big{)}^{\top}\Lambda_{n}$

where

$\displaystyle A_{1}(x)$
$\displaystyle=\frac{\alpha+1}{\alpha+\beta+2}\bigg{(}\begin{bmatrix}-\delta&\delta_{-}\\\
\delta_{+}&\delta\end{bmatrix}x+d\begin{bmatrix}1&1\\\
-1&-1\end{bmatrix}\bigg{)},$ $\displaystyle A_{0}$
$\displaystyle=\frac{\alpha+1}{\alpha+\beta+2}\begin{bmatrix}C_{+}-\Delta&\frac{\beta+1}{\alpha+1}(C+\Delta_{-})\\\
C+\Delta_{+}&C_{-}+\Delta\end{bmatrix},$ $\displaystyle\Lambda_{n}$
$\displaystyle=\frac{\alpha+1}{\alpha+\beta+2}\begin{bmatrix}-\delta&\delta_{-}\\\
\delta_{+}&\delta\end{bmatrix}n+A_{0},$

with $d=l_{1}r_{2}$ and

$\displaystyle\delta$
$\displaystyle=l_{1}r_{2}+\frac{\beta+1}{\alpha+1}l_{2}r_{1},$
$\displaystyle\delta_{-}$
$\displaystyle=-l_{1}r_{2}+\Big{(}\frac{\beta+1}{\alpha+1}\Big{)}^{2}l_{2}r_{1},$
$\displaystyle\delta_{+}$ $\displaystyle=l_{1}r_{2}-l_{2}r_{1},$
$\displaystyle\Delta$ $\displaystyle=(\beta+1)(l_{1}r_{2}+l_{2}r_{1}),$
$\displaystyle\Delta_{-}$
$\displaystyle=-l_{1}r_{2}(\alpha+1)+l_{2}r_{1}(\beta+1),$
$\displaystyle\Delta_{+}$
$\displaystyle=l_{1}r_{2}(\beta+1)-l_{2}r_{1}(\alpha+1),$ $\displaystyle C$
$\displaystyle=-l_{1}r_{1}a_{1}+l_{2}r_{2}a_{2},$ $\displaystyle C_{-}$
$\displaystyle=\frac{\beta+1}{\alpha+1}l_{1}r_{1}a_{1}+l_{2}r_{2}a_{2},$
$\displaystyle C_{+}$
$\displaystyle=l_{1}r_{1}a_{1}+\frac{\beta+1}{\alpha+1}l_{2}r_{2}a_{2}.$

When we take $l_{1}=l_{2}=-1$, $r_{1}=r_{2}=1$ and $a_{1}=\beta+1$ and
$a_{2}=\alpha+1$ we get the first order ordinary differential system in §4 of
[27].

######  Remark 5.

The discussion in this example, regarding the Jacobi polynomials
$p^{(\alpha+1,\beta)}(x)$ and $p^{(\alpha,\beta+1)}(x)$ and the use of the
Gauss’ contiguous relations, connects with the results in [66], Remark 2.8.,
see also [64, 65].

######  Example 2.

Here we analyze the Chebyshev example taken from [27], that gives an example
of a family of matrix orthogonal polynomials which satisfy a first order
linear ordinary differential equation. In §3 of [27] we find a set of MOPS
related with the measure $\tilde{W}(x)\operatorname{d}\mu(x)$ where

$\displaystyle\tilde{W}(x)$ $\displaystyle:=\begin{bmatrix}1&x\\\
x&1\end{bmatrix},$ $\displaystyle\operatorname{d}\mu(x)$
$\displaystyle=\frac{1}{\sqrt{1-x^{2}}}.$

We have a nonsingular leading coefficient $\begin{bmatrix}0&1\\\
1&0\end{bmatrix}$ so that

$\displaystyle\tilde{W}(x)$ $\displaystyle=\begin{bmatrix}0&1\\\
1&0\end{bmatrix}W(x),$ $\displaystyle W(x):=I_{2}x-\begin{bmatrix}0&-1\\\
-1&0\end{bmatrix}.$

Following Remark 4 we shall analyze the Darboux transformations
$\operatorname{d}\mu(x)\mapsto W(x)\operatorname{d}\mu(x)$.

Thus, using (23) we can write the perturbed monic matrix orthogonal
polynomials as follows

$\displaystyle\hat{P}_{n}(x)$
$\displaystyle=Q\begin{bmatrix}\frac{t_{n+1}(x)-\frac{t_{n+1}(-1)}{t_{n}(-1)}t_{n}(x)}{x+1}&0\\\
0&\frac{t_{n+1}(x)-\frac{t_{n+1}(1)}{t_{n}(1)}t_{n}(x)}{x-1}\end{bmatrix}Q^{\top},$
$\displaystyle Q$ $\displaystyle:=\frac{1}{\sqrt{2}}\begin{bmatrix}1&1\\\
1&-1\end{bmatrix},$

where $\\{t_{n}(x)\\}_{n=0}^{\infty}$ are the monic Chebyshev polynomials of
first kind, i.e., $t_{n}(x)=2^{-n+1}T_{n}(x)$ with $T_{n}$ the first kind
Chebyshev polynomial of degree $n$. Therefore, recalling that $T_{n}(\pm
1)=(\pm 1)^{n}$ we get

$\displaystyle\frac{t_{n+1}(\mp 1)}{t_{n}(\mp 1)}$
$\displaystyle=\mp\frac{1}{2},$ $\displaystyle t_{n+1}(x)-\frac{t_{n+1}(\mp
1)}{t_{n}(\mp 1)}t_{n}(x)$ $\displaystyle=\frac{1}{2^{n}}(T_{n+1}(x)\pm
T_{n}(x)).$

Now, recalling the mutual recurrence relation satisfied by Chebyshev
polynomials of the first and second kind, denoted these last ones by $U_{n}$,

$\displaystyle T_{n+1}(x)$ $\displaystyle=xT_{n}(x)-(1-x^{2})U_{n-1}(x),$
$\displaystyle T_{n}(x)$ $\displaystyle=U_{n}(x)-xU_{n-1}(x),$

which implies $T_{n+1}(x)=xU_{n}(x)-U_{n-1}(x)$, we deduce

$\displaystyle T_{n+1}(x)\pm T_{n}(x)=(x\pm 1)(U_{n}(x)\mp U_{n-1}(x)).$

Consequently,

$\displaystyle\hat{P}_{n}(x)$
$\displaystyle=\frac{1}{2^{n}}Q\begin{bmatrix}U_{n}(x)-U_{n-1}(x)&0\\\
0&U_{n}(x)+U_{n-1}(x)\end{bmatrix}Q^{\top},$ $\displaystyle Q$
$\displaystyle:=\frac{1}{\sqrt{2}}\begin{bmatrix}1&1\\\ 1&-1\end{bmatrix}.$

The matrix orthogonal polynomials associated with the original measure
$\tilde{W}(x)\operatorname{d}\mu(x)$ can be recovered from this by a
similarity transformation with $\begin{bmatrix}0&1\\\ 1&0\end{bmatrix}$ so
that

$\displaystyle\tilde{\hat{P}}_{n}(x)$
$\displaystyle=\frac{1}{2^{n+1}}\begin{bmatrix}1&-1\\\
1&1\end{bmatrix}\begin{bmatrix}U_{n}(x)-U_{n-1}(x)&0\\\
0&U_{n}(x)+U_{n-1}(x)\end{bmatrix}\begin{bmatrix}1&1\\\ -1&1\end{bmatrix}$
$\displaystyle=\frac{1}{2^{n}}\begin{bmatrix}U_{n}(x)&-U_{n-1}(x)\\\
-U_{n-1}(x)&U_{n}(x)\end{bmatrix}.$

######  Remark 6.

The polynomials $\\{U_{n}\mp U_{n-1}\\}_{n=0}^{\infty}$ with $U_{-1}=0$, which
are orthogonal with respect to the measures $\frac{x\pm 1}{\sqrt{1-x^{2}}}$,
are the well known Chebyshev polynomials of the third and fourth kind,
respectively.

######  Remark 7.

The symmetric structure of the MOPS can be encoded in the equation

$\displaystyle\begin{bmatrix}0&1\\\
1&0\end{bmatrix}\tilde{\hat{P}}_{n}(x)=\tilde{\hat{P}}_{n}(x)\begin{bmatrix}0&1\\\
1&0\end{bmatrix}.$

As in the Jacobi case, the new two scalar families of orthogonal polynomials
are related through

(26) $\displaystyle\Big{(}(x\mp
1)\frac{\operatorname{d}}{\operatorname{d}x}+\frac{1}{2}\Big{)}\big{(}U_{n}(x)\pm
U_{n-1}(x)\big{)}=(n+\frac{1}{2})\big{(}U_{n}(x)\mp U_{n-1}(x)\big{)}.$

This follows from

$\displaystyle(x\mp 1)(U^{\prime}_{n}(x)\pm U^{\prime}_{n-1}(x))$
$\displaystyle=(x\mp
1)\frac{\operatorname{d}}{\operatorname{d}x}\Big{(}\frac{T_{n+1}(x)\mp
T_{n}(x)}{x\mp 1}\Big{)}$ $\displaystyle=T^{\prime}_{n+1}(x)\mp
T^{\prime}(x)-\frac{T_{n+1}(x)\mp T_{n}(x)}{x\mp 1}$
$\displaystyle=(n+1)U_{n}(x)-nU_{n-1}(x)-U_{n}(x)\mp U_{n-1}(x).$

Here we have used that $T_{n}^{\prime}=nU_{n-1}$.

Differential equation (26) can be written in matrix form

(27) $\bigg{(}\begin{bmatrix}0&x-1\\\
x+1&0\end{bmatrix}\frac{\operatorname{d}}{\operatorname{d}x}+\begin{bmatrix}a_{1}&\frac{1}{2}\\\
\frac{1}{2}&a_{2}\end{bmatrix}\bigg{)}\begin{bmatrix}U_{n}(x)-U_{n-1}(x)&0\\\
0&U_{n}(x)+U_{n-1}(x)\end{bmatrix}\\\ =\begin{bmatrix}U_{n}(x)-U_{n-1}(x)&0\\\
0&U_{n}(x)+U_{n-1}(x)\end{bmatrix}\begin{bmatrix}a_{1}&n+\frac{1}{2}\\\
n+\frac{1}{2}&a_{2}\end{bmatrix}$

where $a_{1},a_{2}\in\mathbb{R}$ are arbitrary constants. Notice also that
this matrix equation is invariant under multiplication on the right and on the
left hand sides by arbitrary diagonal matrices
$L=\operatorname{diag}(l_{1},l_{2})$ and $R=\operatorname{diag}(r_{1},r_{2})$,

$\begin{bmatrix}l_{1}&0\\\ 0&l_{2}\end{bmatrix}\bigg{(}\begin{bmatrix}0&x-1\\\
x+1&0\end{bmatrix}\frac{\operatorname{d}}{\operatorname{d}x}+\begin{bmatrix}a_{1}&\frac{1}{2}\\\
\frac{1}{2}&a_{2}\end{bmatrix}\bigg{)}\begin{bmatrix}r_{1}&0\\\
0&r_{2}\end{bmatrix}\begin{bmatrix}U_{n}(x)-U_{n-1}(x)&0\\\
0&U_{n}(x)+U_{n-1}(x)\end{bmatrix}\\\ =\begin{bmatrix}U_{n}(x)-U_{n-1}(x)&0\\\
0&U_{n}(x)+U_{n-1}(x)\end{bmatrix}\begin{bmatrix}l_{1}&0\\\
0&l_{2}\end{bmatrix}\begin{bmatrix}a&n+\frac{1}{2}\\\
n+\frac{1}{2}&b\end{bmatrix}\begin{bmatrix}r_{1}&0\\\ 0&r_{2}\end{bmatrix}.$

After the similarity transformation we find out that the orthogonal polynomial
$\tilde{\hat{P}}_{n}$ satisfies

$\displaystyle\Big{(}A_{1}(x)\frac{\operatorname{d}}{\operatorname{d}x}+A_{0}\Big{)}\tilde{\hat{P}}_{n}(x)=\tilde{\hat{P}}_{n}(x)\Lambda_{n},$

where

$\displaystyle A_{1}(x)$
$\displaystyle=\begin{bmatrix}-\delta_{+}x+\delta_{-}&\delta_{-}x-\delta_{+}\\\
-\delta_{-}x+\delta_{+}&\delta_{+}x-\delta_{-}\end{bmatrix},$
$\displaystyle\Lambda_{n}$
$\displaystyle=\begin{bmatrix}A_{+}-\delta_{+}\Big{(}n+\frac{1}{2}\Big{)}&A_{-}+\delta_{-}\Big{(}n+\frac{1}{2}\Big{)}\\\
A_{-}-\delta_{-}\Big{(}n+\frac{1}{2}\Big{)}&A_{+}+\delta_{+}\Big{(}n+\frac{1}{2}\Big{)}\end{bmatrix},$

and $A_{0}=\Lambda_{n=0}$ with

$\displaystyle\delta_{\pm}$ $\displaystyle=l_{1}r_{2}\pm l_{2}r_{1},$
$\displaystyle A_{\pm}$ $\displaystyle=l_{1}r_{1}a_{1}\pm l_{2}r_{2}a_{2}.$

Equations (3.1) and (3.2) of §3 of [27] can be recovered choosing
$\Big{(}\delta_{+}=A_{+}=0,\delta_{-}=1,A_{-}=-\frac{1}{2}\Big{)}$ and
$\Big{(}\delta_{-}=A_{-}=0,\delta_{+}=-1,A_{+}=\frac{1}{2}\Big{)}$,
respectively.

However, they are all equivalent to (27), another form of writing (26). It is
in fact this last equation (26) a quite interesting one. Indeed, we have two
families of Darboux transformed orthogonal polynomials interconnected by two
first order differential equations. Moreover, we conclude

$\displaystyle\Big{(}(x^{2}-1)\frac{\operatorname{d}^{2}}{\operatorname{d}x^{2}}+(2x\mp
1)\frac{\operatorname{d}}{\operatorname{d}x}+\frac{1}{4}\Big{)}\big{(}U_{n}(x)\mp
U_{n-1}(x)\big{)}=\Big{(}n+\frac{1}{2}\Big{)}^{2}\big{(}U_{n}(x)\mp
U_{n-1}(x)\big{)}.$

#### Example 3

Here we comment on the matrix Gegenbauer matrix valued polynomials discussed
in [66]. In this case the matrix of weights is a symmetric matrix,
$W^{(\nu)}:[-1,1]\to\mathbb{R}^{N\times N}$, with matrix coefficients of the
form

$\displaystyle(W^{(\nu)}(x))_{i,j}$
$\displaystyle:=(1-x^{2})^{\nu-1/2}\sum_{k=\operatorname{max}(0,i+j+1-N)}\alpha^{(\nu)}_{k}(i,j)C^{(\nu)}_{i+j-2k}(x),$
$\displaystyle i\geq j,$

where $\alpha^{(\nu)}_{t}(i,j)$ are some coefficients and $C^{(\nu)}_{n}(x)$
stands for the Gegenbauer or ultraspherical polynomials. Erik Koelink and
Pablo Román kindly communicated us a nice feature of the matrix Christoffel
transformation discussed in this paper when acting on this reach family of
MOPS: two families of Gegenbauer MOPS associated with matrices of weights
$W^{(\nu_{1})}(x)$ and $W^{(\nu_{2})}(x)$, such that
$\nu_{1}-\nu_{2}=m\in\mathbb{Z}$, are linked by a matrix Christoffel
transformation. Now, the perturbing polynomial $W(x)$ has $\deg W=2m$. These
examples are, in general, reducible to two irreducible blocks of sizes $N/2$,
for $N$ even, and $(N+1)/2$ and $(N-1)/2$ for odd $N$. For a discussion on the
orthogonal and non orthogonal reducibility of these examples see [66, 67].

##  4\. Singular leading coefficient matrix polynomial perturbations

After studying some examples that the literature provides us with, one may
realize that, even thought it is generic to assume the perturbing matrix
polynomial ${W}(x)$ to have a nonsingular leading coefficient, many examples
do have a singular matrix as its leading coefficient. This situation is a
special feature of the matrix case setting since in the scalar case, having a
singular leading term would mean that this coefficient is just zero
(affecting, of course, to the degree of the polynomial). For this reason, when
dealing with this kind of matrix polynomials talking about their degree should
make no sense. The effect that this fact has on our reasoning is that since
$\deg[\det{W}(x)]\leq Np$ the information encoded in the zeroes (and
corresponding adapted polynomials) of $\det{W}(x)$ is no longer enough to make
the matrices $\Pi_{kN}$ of the needed size. Therefore, there will be no way to
express the perturbed polynomials just in terms of the initial ones evaluated
at the zeroes of $\det{W}(x)$ and the method to find a Christoffel type
formula fails. However, the information that seems to be missing in these
cases may actually not be necessary due to the singular character of the
leading coefficient of the perturbing polynomial. Let us consider the
following example to take a glimpse of this scenario.

Let us pick up some scalar measure $\operatorname{d}\mu(x)$ and its associated
monic OPs $\big{\\{}p_{k}(x)\big{\\}}_{k=0}^{\infty}$ together with their
norms and three term recurrence relation

$\displaystyle h_{k}\delta_{kj}$ $\displaystyle:=\langle p_{k},p_{j}\rangle,$
$\displaystyle xp_{k}(x)$
$\displaystyle=J_{k,k-1}p_{k-1}(x)+J_{k,k}p_{k}(x)+p_{k+1}(x),$ $\displaystyle
J_{k,k-1}$ $\displaystyle=\frac{h_{k}}{h_{k-1}}>0.$

Now, consider its $2q\times 2q$ matrix diagonal extension
$\in\mathbb{R}^{2q\times 2q}[x]$

$\displaystyle P_{k}(x)$ $\displaystyle:=p_{k}(x)I_{2q},$ $\displaystyle
H_{k}$ $\displaystyle:=h_{k}I_{2q}.$

Our aim is to consider the following matrix polynomial (with singular leading
coefficient)

$\displaystyle{W}(x)$
$\displaystyle:=\begin{bmatrix}I_{q}+AA^{\top}x^{2}&Ax\\\
A^{\top}x&I_{q}\end{bmatrix},$ $\displaystyle A\in\mathbb{R}^{q\times q},$

which is inspired by the $q=1$ case $\begin{bmatrix}1+a^{2}x^{2}&ax\\\
ax&1\end{bmatrix}$ (see [45], and references therein) and study the
corresponding perturbations of our initial scalar measure; i.e.,
$\operatorname{d}\hat{\mu}(x):={W}(x)\operatorname{d}\mu(x)$ in order to
obtain the transformed matrix orthogonal polynomials

$\displaystyle\hat{P}_{k}(x)$
$\displaystyle:=\begin{bmatrix}(\hat{P}_{k})_{1,1}&(\hat{P}_{k})_{1,2}\\\
(\hat{P}_{k})_{2,1}&(\hat{P}_{k})_{2,2}\end{bmatrix},$
$\displaystyle\hat{P}_{k}(x)\in\mathbb{R}^{2q\times
2q}[x],(\hat{P}_{k})_{i,j}\in\mathbb{R}^{q\times q}[x],$
$\displaystyle\langle\hat{P}_{k},\hat{P}_{j}\rangle_{W}$
$\displaystyle:=\delta_{k,j}\hat{H}_{k}=\delta_{k,j}\begin{bmatrix}(\hat{H}_{k})_{1,1}&(\hat{H}_{k})_{1,2}\\\
(\hat{H}_{k})_{2,1}&(\hat{H}_{k})_{2,2}\end{bmatrix},$
$\displaystyle\hat{H}_{k}\in\mathbb{R}^{2q\times
2q},(\hat{H}_{k})_{i,j}\in\mathbb{R}^{q\times q}.$

We have splitted them up this way for computational purposes. Notice that
since ${W}(x)={W}(x)^{\top}$ we have
$\hat{M}=\hat{M}^{\top}:=\hat{S}^{-1}\hat{H}[\hat{S}^{-1}]^{\top}$ and,
therefore, $\hat{P}^{[1]}=\hat{P}^{[2]}:=\hat{P}$ and
$\hat{H}_{k}=(\hat{H}_{k})^{\top}$.

Let us point out that

$\displaystyle{W}(x)$ $\displaystyle=\mathcal{W}(x)\mathcal{W}(x)^{\top},$
$\displaystyle\mathcal{W}$ $\displaystyle:=\begin{bmatrix}I_{q}&Ax\\\
0&I_{q}\end{bmatrix},$
$\displaystyle\mathcal{W}^{-1}=\begin{bmatrix}I_{q}&-Ax\\\
0&I_{q}\end{bmatrix}.$

This implies that $\det W=\det\mathcal{W}=1$ and, consequently, there is no
spectral analysis to perform as there are non eigenvalues at all. Thus, the
relation between the original and perturbed measures and moment matrices is

$\displaystyle[\mathcal{W}(x)]^{-1}\operatorname{d}\hat{\mu}$
$\displaystyle=\operatorname{d}\mu[\mathcal{W}(x)]^{\top},$
$\displaystyle[\mathcal{W}(\Lambda)]^{-1}\hat{M}$
$\displaystyle=M[\mathcal{W}(\Lambda)]^{\top}.$

######  Definition 14.

We introduce the resolvent or connection matrix

$\displaystyle\omega:=\hat{S}\mathcal{W}(\Lambda)S^{-1}.$

######  Proposition 26.

The matrix $\omega$ is block tridiagonal, having only its diagonal and first
superdiagonal and subdiagonal nonzero, and satisfies

$\displaystyle\omega^{-1}$ $\displaystyle=H[\omega]^{\top}\hat{H}^{-1}.$

Moreover, we have the important connection formula

$\displaystyle\omega P=\hat{P}\mathcal{W}(x).$

###### Proof.

The first relation is a consequence of the $LU$ factorization of the moment
matrices and the connection formula is a straightforward consequence of the
definition of $\omega$. ∎

######  Proposition 27.

  1. i)

The matrices

$\displaystyle\rho_{k+1}$
$\displaystyle:=\left(I_{q}+J_{k+1,k}A^{\top}A\right)^{-1},$ $\displaystyle k$
$\displaystyle\in\\{-1,0,1,\dots\\},$

exist.

  2. ii)

The perturbed MOPS can be writen in terms of the original OPS as follows

$\displaystyle\hat{P}_{1,k+1}(x)\mathcal{W}(x)$
$\displaystyle=-\begin{bmatrix}J_{k+1,k}J_{k+1,k+1}A\rho_{k+1}A^{\top}&0\\\
J_{k+1,k}\rho_{k+1}A^{\top}&0\end{bmatrix}p_{k}(x)+\begin{bmatrix}I_{q}&J_{k+1,k+1}A\rho_{k+1}\\\
0&\rho_{k+1}\end{bmatrix}p_{k+1}(x)+\begin{bmatrix}0&A\\\
0&0\end{bmatrix}p_{k+2}(x),$

for $k\in\\{-1,0,1,\dots\\}.$

###### Proof.

From the $(k+1)$-th row of the connection formula we have that

$\displaystyle\omega_{k+1,k}p_{k}(x)+\omega_{k+1,k+1}p_{k+1}(x)+\omega_{k+1,k+2}p_{k+2}(x)=\hat{P}_{k+1}(x)\mathcal{W}(x),$

but from the Definition 14 and Proposition 26 one realizes that the previous
expression reads

$\displaystyle\hat{P}_{k+1}(x)\mathcal{W}(x)=\hat{H}_{k+1}\begin{bmatrix}0&-A\\\
0&0\end{bmatrix}^{\top}h_{k}^{-1}p_{k}(x)+\begin{bmatrix}(\omega_{k+1,k+2})_{11}&(\omega_{k+1,k+2})_{12}\\\
(\omega_{k+1,k+2})_{21}&(\omega_{k+1,k+2})_{22}\end{bmatrix}p_{k+1}(x)+\begin{bmatrix}0&A\\\
0&0\end{bmatrix}p_{k+2}(x).$

Now, taking into account that both $(\hat{P}_{k+1})_{11},(\hat{P}_{k+1})_{22}$
are monic $q\times q$ polynomials of degree $k+1$, while
$(\hat{P}_{k+1})_{12},(\hat{P}_{k+1})_{21}$ are $q\times q$ polynomials of
degree less than $k+1$, it is not hard to see (after using the three term
recurrence relation of the initial polynomials) that

$\displaystyle(\omega_{k+1,k+2})_{11}$ $\displaystyle=I_{q},$
$\displaystyle(\omega_{k+1,k+2})_{12}$
$\displaystyle=J_{k+1,k+1}A-h_{k}^{-1}(\hat{H}_{k+1})_{12}A^{\top}A,$
$\displaystyle(\omega_{k+1,k+2})_{21}$ $\displaystyle=0,$
$\displaystyle(\omega_{k+1,k+2})_{22}$
$\displaystyle=I_{q}-h_{k}^{-1}(\hat{H}_{k+1})_{22}A^{\top}A.$

Hence, we have every coefficient that appears in the connection formula in
terms of the still unknown norms of the MOPS. Therefore, we just need to
compute the second block column of the following integral

$\int\left[\omega_{k+1,k}p_{k}(x)+\omega_{k+1,k+1}p_{k+1}(x)+\omega_{k+1,k+2}p_{k+2}(x)\right]\left[(\mathcal{W}(x))^{\top}x^{k+1}\right]\operatorname{d}\mu(x)\\\
\begin{aligned}
&=\int\hat{P}_{k+1}(x)\mathcal{W}(x)\left[(\mathcal{W}(x))^{\top}x^{k+1}\right]\operatorname{d}\mu(x)\\\
&=\int\hat{P}_{k+1}(x)\operatorname{d}\hat{\mu}(x)x^{k+1}\\\
&=\hat{H}_{k+1},\end{aligned}$

which yields

$\displaystyle(\hat{H}_{k+1})_{12}$
$\displaystyle=J_{k+1,k+1}h_{k+1}A\left(I_{q}+J_{k+1,k}A^{\top}A\right)^{-1},$
$\displaystyle(\hat{H}_{k+1})_{22}$
$\displaystyle=h_{k+1}\left(I_{q}+J_{k+1,k}A^{\top}A\right)^{-1}.$

Therefore,

$\displaystyle(\omega_{k+1,k+2})_{12}$
$\displaystyle=J_{k+1,k+1}A\left(I_{q}+J_{k+1,k}A^{\top}A\right)^{-1},$
$\displaystyle(\omega_{k+1,k+2})_{22}$
$\displaystyle=\left(I_{q}+J_{k+1,k}A^{\top}A\right)^{-1},$

and the stated result follows. ∎

For $q=1$ and the classical measures we have, see [47],

######  Corollary 1.

Starting from the classical measures

  1. i)

Hermite monic polynomials
$\big{\\{}\mathcal{H}_{k}(x)\big{\\}}_{k=0}^{\infty}$ with norm
$h_{k}=\pi^{\frac{1}{2}}\frac{k!}{2^{k}}$

$\displaystyle J_{k+1,k}$ $\displaystyle=\frac{k+1}{2},$ $\displaystyle
J_{k+1,k+1}$ $\displaystyle=0,$
$\displaystyle\rho_{k+1}:=\frac{2}{2+a^{2}(k+1)}.$

  2. ii)

Laguerre monic polynomials
$\big{\\{}\mathcal{L}^{\alpha}_{k}(x)\big{\\}}_{k=0}^{\infty}$ with norm
$h_{k}=k!\Gamma(k+1+\alpha)$

$\displaystyle J_{k+1,k}$ $\displaystyle=(k+1)(k+\alpha+2),$ $\displaystyle
J_{k+1,k+1}$ $\displaystyle=(2k+\alpha+3),$
$\displaystyle\rho_{k+1}:=\frac{1}{1+a^{2}(k+1)(k+1+\alpha)}.$

and perturbing them by the matrix polynomial

$\displaystyle{W}(x)$
$\displaystyle=\mathcal{W}(x)\left(\mathcal{W}(x)\right)^{\top},$
$\displaystyle\mathcal{W}(x)=\begin{pmatrix}1&ax\\\ 0&1\end{pmatrix},$

one obtains the perturbed MOPS related to the classical OPS as follows

$\displaystyle\hat{\mathcal{H}}_{k+1}(x)\mathcal{W}(x)$
$\displaystyle=\begin{pmatrix}0&0\\\
\frac{\rho_{k+1}-1}{a}&0\end{pmatrix}\mathcal{H}_{k}(x)+\begin{pmatrix}1&0\\\
0&\rho_{k+1}\end{pmatrix}\mathcal{H}_{k+1}(x)+\begin{pmatrix}0&a\\\
0&0\end{pmatrix}\mathcal{H}_{k+2}(x),$
$\displaystyle\hat{\mathcal{L}}^{\alpha}_{k+1}(x)\mathcal{W}(x)$
$\displaystyle=\begin{pmatrix}(\rho_{k+1}-1)(2k+3+\alpha)&0\\\
\frac{\rho_{k+1}-1}{a}&0\end{pmatrix}\mathcal{L}^{\alpha}_{k}(x)+\begin{pmatrix}1&a(2k+3+\alpha)\rho_{k+1}\\\
0&\rho_{k+1}\end{pmatrix}\mathcal{L}^{\alpha}_{k+1}(x)+\begin{pmatrix}0&a\\\
0&0\end{pmatrix}\mathcal{L}^{\alpha}_{k+2}(x).$

##  5\. Extension to non-Abelian 2D Toda hierarchies

Matrix orthogonal polynomials are connected with non-Abelian Toda lattices,
see [80, 11].

###  5.1. Block Hankel moment matrices vs multi-component Toda hierarchies

Let us take $M=(m_{i,j})_{i,j=0}^{\infty}$, $m_{i,j}\in\mathbb{R}^{p\times p}$
a semi-infinite block matrix having a Gaussian factorization

$\displaystyle M=(S_{1})^{-1}H(S_{2})^{-\top},$

where $S_{1},S_{2}$ are lower uni-triangular block matrices and $H$ is block
diagonal. Notice that conditions for this factorization to hold were given in
Proposition 10.

######  Definition 15.

We introduce some continuous flows or perturbations of this semi-infinite
matrix. For that aim we first consider the diagonal matrices

$\displaystyle t_{i,j}=$
$\displaystyle\operatorname{diag}(t_{i,j,1},\dots,t_{i,j,p})\in\mathbb{R}^{p\times
p},$ $\displaystyle i$ $\displaystyle=1,2,$ $\displaystyle j$
$\displaystyle\in\mathbb{Z}_{+},$

the semi-infinite undressed wave matrices

$\displaystyle V_{i}^{(0)}$
$\displaystyle:=\exp\Big{(}\sum_{j=0}^{\infty}t_{i,j}\Lambda^{j}\Big{)},$
$\displaystyle i$ $\displaystyle=1,2,$

and the perturbed matrix $M(t)$, $t=(t_{1},t_{2})$,
$t_{i}=\\{t_{i,j,a}\\}_{\begin{subarray}{c}j\in\mathbb{Z}_{+}\\\
a\in\\{1,\dots,p\\}\end{subarray}}$

$\displaystyle M(t)=$ $\displaystyle
V_{1}^{(0)}(t_{1})M\big{(}V_{2}^{(0)}(t_{2})\big{)}^{-\top}.$

Observe that we do not require any Hankel form for the matrix $M$, modelled by
$\Lambda M=M\Lambda^{\top}$. However, if $M(0)$ is a Hankel matrix $M(t)$ is
also a Hankel matrix taking into account $\Lambda M(t)=M(t)\Lambda^{\top}$.
Hence, if $\operatorname{d}\mu(x)$ is the initial matrix of measures, then the
new matrix of measures $\operatorname{d}\mu(x,t)$ will be

$\displaystyle\operatorname{d}\mu(x,t)=\exp\Big{(}\sum_{j=0}^{\infty}t_{1,j}x^{j}\Big{)}\operatorname{d}\mu(x)\exp\Big{(}-\sum_{j=0}^{\infty}t_{2,j}x^{j}\Big{)}.$

Here $M(t)$ will be the moment matrix of the matrix of measures. Moreover, if
at any time the matrix of measures is block Hankel then it was and it will be
a Hankel block matrix at any time. If we assume that we can perform the
Gaussian factorization again we can write

$\displaystyle M(t)$ $\displaystyle=(S_{1}(t))^{-1}H(t)(S_{2}(t))^{-\top}.$

As we know, for the block Hankel case we are dealing with bi-orthogonal or
orthogonal polynomials with respect to the associated matrix of measures. What
happens in the general case? Following [3] and [71] we can understand the
Gaussian factorization also as a bi-orthogonality condition. The semi-infinite
vectors of polynomials will be

$\displaystyle P^{[1]}(x)$ $\displaystyle:=S_{1}(t)\chi(x),$ $\displaystyle
P^{[2]}(x)$ $\displaystyle:=S_{2}(t)\chi(x),$

and we consider a sesquilinear form in $\mathbb{R}^{p\times p}[x]$, see §1.4,
that for any couple of matrix polynomials $P=\sum_{k=0}^{\deg P}p_{k}x^{k}$
and $Q(x)=\sum_{l=0}^{\deg Q}q_{l}x^{l}$ is defined by

$\displaystyle\left\langle{P(x),Q(x)}\right\rangle=\sum_{\begin{subarray}{c}k=1,\dots,\deg
P\\\ l=1,\dots,\deg Q\end{subarray}}p_{k}M_{k,l}(q_{l})^{\top},$

where we can interpret

$\displaystyle M_{k,l}=\left\langle{x^{k}1_{p},x^{l}1_{p}}\right\rangle$

as the Gram matrix of the sesquilinear form. With respect to this sesquilinear
form we have the bi-orthogonality condition

$\displaystyle\left\langle{P_{k}^{[1]}(x),P_{l}^{[2]}(x)}\right\rangle=H_{k}\delta_{k,l}.$

For a block Hankel initial condition this sesquilinear form is just the
sequilinear product associated with a linear functional of a measure. In [10]
different examples are discussed for the matrix orthogonal polynomials
scenario. For example, multigraded Hankel matrices $M$ fulfilling

$\displaystyle\Big{(}\sum_{a=1}^{p}\Lambda^{n_{a}}E_{a,a}\Big{)}M=M\Big{(}\sum_{a=1}^{p}\big{(}\Lambda^{\top}\big{)}^{m_{a}}E_{a,a}\Big{)},$

where $n_{1},\dots,n_{p},m_{1},\dots,m_{p}$ are positive integers, can be
realized as

$\displaystyle M_{k,l}=\int x^{k}\operatorname{d}\mu^{(l)}(x)$

in terms of matrices of measures $\operatorname{d}\mu^{(l)}(x)$ which satisfy
the following periodicity condition

(28)
$\displaystyle\operatorname{d}\mu^{(l+m_{a})}_{a,b}(x)=x^{n_{a}}\operatorname{d}\mu^{(l)}_{a,b}(x).$

Therefore, given the measures
$\operatorname{d}\mu^{(0)}_{a,b},\dots,\operatorname{d}\mu^{(m_{b}-1)}_{a,b}$
we can recover all the others from (28). In this case, we have generalized
orthogonality conditions like

$\displaystyle\int P^{[1]}_{k}(x)\operatorname{d}\mu^{(l)}(x)$
$\displaystyle=0,$ $\displaystyle l=0,\dots,k-1.$

Coming back to the Gaussian factorization, we consider the wave matrices

$\displaystyle V_{1}(t):=$ $\displaystyle S_{1}(t)V_{1}^{(0)}(t_{1}),$
$\displaystyle\tilde{V}_{2}(t):=$
$\displaystyle\tilde{S}_{2}(t)(V^{(0)}_{2}(t_{2}))^{\top},$

where $\tilde{S}_{2}(t):=H(t)(S_{2}(t))^{-\top}$.

######  Proposition 28.

The wave matrices satisfy

(29) $\displaystyle\big{(}V_{1}(t)\big{)}^{-1}\tilde{V}_{2}(t)=M.$

###### Proof.

It is a consequence of the Gaussian factorization. ∎

Given a semi-infinite matrix $A$ we have unique splitting $A=A_{+}+A_{-}$
where $A_{+}$ is an upper triangular block matrix while is $A_{-}$ a strictly
lower triangular block matrix.

######  Proposition 29.

The following equations hold

$\displaystyle\frac{\partial S_{1}}{\partial t_{1,j,a}}(S_{1})^{-1}$
$\displaystyle=-\Big{(}S_{1}E_{a,a}\Lambda^{j}(S_{1})^{-1}\Big{)}_{-},$
$\displaystyle\frac{\partial S_{1}}{\partial t_{2,j,a}}(S_{1})^{-1}$
$\displaystyle=\Big{(}\tilde{S}_{2}E_{a,a}\big{(}\Lambda^{\top}\big{)}^{j}(\tilde{S}_{2})^{-1}\Big{)}_{-},$
$\displaystyle\frac{\partial\tilde{S}_{2}}{\partial
t_{1,j,a}}(\tilde{S}_{2})^{-1}$
$\displaystyle=\Big{(}S_{1}E_{a,a}\Lambda^{j}(S_{1})^{-1}\Big{)}_{+},$
$\displaystyle\frac{\partial\tilde{S}_{2}}{\partial
t_{2,j,a}}(\tilde{S}_{2})^{-1}$
$\displaystyle=-\Big{(}\tilde{S}_{2}E_{a,a}\big{(}\Lambda^{\top}\big{)}^{j}(\tilde{S}_{2})^{-1}\Big{)}_{+}.$

###### Proof.

Taking right derivatives of (29) yields

$\displaystyle\frac{\partial V_{1}}{\partial t_{i,j,a}}(V_{1})^{-1}$
$\displaystyle=\frac{\partial\tilde{V}_{2}}{\partial
t_{i,j,a}}(\tilde{V}_{2})^{-1},$ $\displaystyle i$
$\displaystyle\in\\{1,2\\},$ $\displaystyle j$
$\displaystyle\in\mathbb{Z}_{+},$

where

$\displaystyle\frac{\partial V_{1}}{\partial t_{1,j,a}}(V_{1})^{-1}$
$\displaystyle=\frac{\partial S_{1}}{\partial
t_{1,j,a}}(S_{1})^{-1}+S_{1}E_{a,a}\Lambda^{j}(S_{1})^{-1},$
$\displaystyle\frac{\partial V_{1}}{\partial t_{2,j,a}}(V_{1})^{-1}$
$\displaystyle=\frac{\partial S_{1}}{\partial t_{2,j,a}}(S_{1})^{-1},$
$\displaystyle\frac{\partial\tilde{V}_{2}}{\partial
t_{1,j,a}}(\tilde{V}_{2})^{-1}$
$\displaystyle=\frac{\partial\tilde{S}_{2}}{\partial
t_{1,j,a}}(\tilde{S}_{2})^{-1},$
$\displaystyle\frac{\partial\tilde{V}_{2}}{\partial
t_{2,j,a}}(\tilde{V}_{2})^{-1}$
$\displaystyle=\frac{\partial\tilde{S}_{2}}{\partial
t_{2,j,a}}(\tilde{S}_{2})^{-1}+\tilde{S}_{2}E_{a,a}\big{(}\Lambda^{\top}\big{)}^{j}(\tilde{S}_{2})^{-1},$

and the result follows immediately. ∎

As a consequence, we derive

######  Proposition 30.

The multicomponent 2D Toda lattice equations

$\displaystyle\frac{\partial}{\partial t_{2,1,b}}\Big{(}\frac{\partial
H_{k}}{\partial
t_{1,1,a}}(H_{k})^{-1}\Big{)}+E_{a,a}H_{k+1}E_{b,b}(H_{k})^{-1}-H_{k}E_{b,b}(H_{k-1})^{-1}E_{a,a}=0.$

hold.

###### Proof.

From Proposition 29 we get

$\displaystyle\frac{\partial H_{k}}{\partial t_{1,1,a}}(H_{k})^{-1}$
$\displaystyle=\beta_{k}E_{a,a}-E_{a,a}\beta_{k+1},$
$\displaystyle\frac{\partial\beta_{k}}{\partial t_{2,1,b}}$
$\displaystyle=H_{k}E_{b,b}(H_{k-1})^{-1},$

where $\beta_{k}\in\mathbb{R}^{p\times p}$, $k=1,2,\dots$, are the first
subdiagonal coefficients in $S_{1}$. ∎

The multi-component Toda and KP hierarchies were introduced in [100]. In [70,
71] its relevance in integrable aspects of differential geometry was
emphasized, and in [62] a representation approach was developed, while in [2,
11] it was used in relation with multiple orthogonality. A comprehensive
approach to multi-component 2D Toda hierarchy with applications in
dispersionless integrability or generalized orthogonal polynomials can be
found in [72, 73, 10].

If we introduce the total flows given by the derivatives

$\displaystyle\partial_{i,j}:=\sum_{a=1}^{p}\frac{\partial}{\partial
t_{i,j,a}}$

we get the non-Abelian 2D Toda lattice

$\displaystyle\partial_{2,1}\big{(}\partial_{1,1}(H_{k})\cdot(H_{k})^{-1}\big{)}+H_{k+1}(H_{k})^{-1}-H_{k}(H_{k-1})^{-1}=0.$

The non-Abelian Toda lattice was introduced in the context of string theory by
Polyakov, [86, 87], and then studied under the inverse spectral transform by
Mikhailov [79] and Riemann surface theory by Krichever [69]. The Darboux
transformations were considered in [90] and later in [84].

The non-Abelian 2D Toda lattice hierarchy is a reduction of the multicomponent
hierarchy by taking the diagonal time matrices
$t_{i,j}=\operatorname{diag}(t_{i,j,1},\dots,t_{i,j,p})$ proportional to the
identity; i.e.,

$\displaystyle t_{i,j}$ $\displaystyle\mapsto t_{i,j}I_{p},$ $\displaystyle
t_{i,j}$ $\displaystyle\in\mathbb{R}.$

These equations are just the first members of an infinite set of nonlinear
partial differential equations, an integrable hierarchy. Its elements are
given by

######  Definition 16.

The partial, Lax and Zakharov–Shabat matrices are given by

$\displaystyle\Pi_{1,a}$ $\displaystyle:=S_{1}E_{a,a}(S_{1})^{-1},$
$\displaystyle\Pi_{2,a}$
$\displaystyle:=\tilde{S}_{2}E_{a,a}(\tilde{S}_{2})^{-1},$ $\displaystyle
L_{1}$ $\displaystyle:=S_{1}\Lambda(S_{1})^{-1},$ $\displaystyle L_{2}$
$\displaystyle:=\tilde{S}_{2}\Lambda^{\top}(\tilde{S}_{2})^{-1},$
$\displaystyle B_{1,j,a}$
$\displaystyle:=\big{(}\Pi_{1,a}(L_{1})^{j}\big{)}_{+},$ $\displaystyle
B_{2,j,a}$ $\displaystyle:=\big{(}\Pi_{2,a}(L_{2})^{j}\big{)}_{-}.$

######  Proposition 31 (The integrable hierarchy).

The wave matrices obey the evolutionary linear systems

$\displaystyle\frac{\partial V_{1}}{\partial t_{1,j,a}}$
$\displaystyle=B_{1,j,a}V_{1},$ $\displaystyle\frac{\partial V_{1}}{\partial
t_{2,j,a}}$ $\displaystyle=B_{2,j,a}V_{1},$
$\displaystyle\frac{\partial\tilde{V}_{2}}{\partial t_{1,j,a}}$
$\displaystyle=B_{1,j,a}\tilde{V}_{2},$
$\displaystyle\frac{\partial\tilde{V}_{2}}{\partial t_{2,j,a}}$
$\displaystyle=B_{2,j,a}\tilde{V}_{2},$

the partial and Lax matrices are subject to the following _Lax equations_

$\displaystyle\frac{\partial\Pi_{i^{\prime},a^{\prime}}}{\partial t_{i,j,a}}$
$\displaystyle=\Big{[}B_{i,j,a},\Pi_{i^{\prime},a^{\prime}}\Big{]},$
$\displaystyle\frac{\partial L_{i^{\prime}}}{\partial t_{i,j,a}}$
$\displaystyle=\Big{[}B_{i,j,a},L_{i^{\prime}}\Big{]},$

and Zakharov–Sabat matrices fulfill the following _Zakharov–Shabat equations_

$\displaystyle\frac{\partial B_{i^{\prime},j^{\prime},a^{\prime}}}{\partial
t_{i,j,a}}-\frac{\partial B_{i,j,a}}{\partial
t_{i^{\prime},j^{\prime},a^{\prime}}}+\big{[}B_{i,j,a},B_{i^{\prime},j^{\prime},a^{\prime}}\big{]}=0.$

###### Proof.

Follows from Proposition 29. ∎

Given two semi-infinite block matrices $A,B$ the notation $[A,B]=AB-BA$ stands
for the usual commutator of matrices.

A crucial observation, regarding orthogonal polynomials, must be pointed out.
When orthogonal polynomials are involved, and the matrices to factorize are
block Hankel, equivalently $\Lambda M=M\Lambda^{\top}$, we get
$L_{1}=S_{1}\Lambda
S_{1}^{-1}=\tilde{S}_{2}\Lambda^{\top}\tilde{S}_{2}^{-1}=L_{2}$. As the reader
may have noticed the Lax matrices $L_{1}$ and $L_{2}$ are, by construction,
lower and upper Hessenberg block matrices, respectively. However, when the
Hankel property holds both Lax matrices are equal,

$\displaystyle L_{1}=L_{2},$

and, therefore, we are faced to a tridiagonal block matrix; i.e., a Jacobi
block matrix. Moreover, this Hankel condition implies an invariance property
under the flows introduced, as we have that $M(t)=V_{1}^{(0)}(t_{1}-t_{2})M$,
i.e., there are only one type of flows. This condition also implies that for
the total flows we have

$\displaystyle(\partial_{1,j}+\partial_{2,j})V_{1}$
$\displaystyle=V_{1}\Lambda^{j},$
$\displaystyle(\partial_{1,j}+\partial_{2,j})\tilde{V}_{2}$
$\displaystyle=\tilde{V}_{2}(\Lambda^{\top})^{j},$
$\displaystyle(\partial_{1,j}+\partial_{2,j})L_{1}$ $\displaystyle=0,$
$\displaystyle(\partial_{1,j}+\partial_{2,j})L_{2}$ $\displaystyle=0.$

Therefore, in the block Hankel case we are dealing with the multicomponent 1D
Toda hierarchy.

###  5.2. The Christoffel transformation for the non-Abelian 2D Toda hierarchy

The idea is to follow what we did in §2.1 and consider an initial condition
$\hat{M}$ at $t=0$, this is

$\displaystyle\hat{M}=W(\Lambda)M$

for a matrix polynomial $W(x)\in\mathbb{R}^{p\times p}[x]$. Observe that using
the scalar times $t_{i,j}\in\mathbb{R}$ of the non-Abelian flows determined by

$\displaystyle V_{i}^{(0)}$
$\displaystyle:=\exp\Big{(}\sum_{j=0}^{\infty}t_{i,j}\Lambda^{j}\Big{)},$
$\displaystyle i$ $\displaystyle=1,2,$

the perturbed matrix is given by

$\displaystyle\hat{M}(t)=$ $\displaystyle
V_{1}^{(0)}(t_{1})\hat{M}\big{(}V_{2}^{(0)}(t_{2})\big{)}^{-\top}$
$\displaystyle=$ $\displaystyle W(\Lambda)M(t).$

Here we have used that $[W(\Lambda),V^{(0)}_{1}(t)]=0$, $\forall
t_{1,j}\in\mathbb{R}$. Let us stress that we could request only $t_{1,j}$ to
be scalars an let $t_{2,j}$ to be diagonal matrices. Despite this is a more
general situation, we prefer to show how the method works in this simpler
scenario.

Assuming that the block Gauss factorization hold, we proceed as in §2.1 and
introduce the resolvents

$\displaystyle\omega^{[1]}(t)$
$\displaystyle:=\hat{S}_{1}(t)W(\Lambda)\big{(}S_{1}(t)\big{)}^{-1},$
$\displaystyle\omega^{[2]}(t)$
$\displaystyle:=(S_{2}(t)\big{(}\hat{S}_{2}(t)\big{)}^{-1})^{\top}.$

From the $LU$ factorization we get

$\displaystyle\big{(}\hat{S}_{1}(t)\big{)}^{-1}\hat{H}(t)\big{(}\hat{S}_{2}(t)\big{)}^{-\top}=W(\Lambda)\big{(}S_{1}(t)\big{)}^{-1}H(t)\big{(}S_{2}(t)\big{)}^{-\top},$

so that

$\displaystyle\hat{H}(t)(S_{2}(t)\big{(}\hat{S}_{2}(t)\big{)}^{-1})^{\top}=\hat{S}_{1}(t)W(\Lambda)\big{(}S_{1}(t)\big{)}^{-1}H(t),$

and, consequently,

$\displaystyle\hat{H}(t)\omega^{[2]}(t)=\omega^{[1]}(t)H(t)$

holds. Hence, as in the static case where the variable $t$ does not appear, we
have that this $t$-dependent resolvent matrix has a band block upper
triangular structure

$\displaystyle\omega^{[1]}=\begin{bmatrix}\omega^{[1]}_{0,0}&\omega^{[1]}_{0,1}&\omega^{[1]}_{0,2}&\dots&\omega^{[1]}_{0,N-1}&I_{p}&0&0&\dots\\\
0&\omega^{[1]}_{1,1}&\omega^{[1]}_{1,2}&\dots&\omega^{[1]}_{1,N-1}&\omega^{[1]}_{1,N}&I_{p}&0&\dots\\\
0&0&\omega^{[1]}_{2,2}&\dots&\omega^{[1]}_{2,N-1}&\omega^{[1]}_{2,N}&\omega^{[1]}_{2,N+1}&I_{p}&\ddots\\\
&\ddots&\ddots&\ddots&&&&\ddots&\ddots\end{bmatrix}$

with

$\displaystyle\hat{H}_{k}(t)=\omega^{[1]}_{k,k}(t)H_{k}(t),$

and the connection formulas described in Proposition 19 hold in this wider
context.

Moreover, if $W(x)$ is a monic polynomial we can ensure that the Christoffel
formula is also fulfilled for the non-Abelian 2D Toda and Theorem 2 remains
valid also in this scenario. Formulas (12) and (18) hold directly and need no
further explanation. However, (15) needs the following brief discussion. The
Christoffel–Darboux kernel is defined exactly as we did in (4) but very
probably there is no such a formula as the CD formula given in Proposition 16
is present in this scenario. However, as was shown in [12], there are cases,
such as the multigraded reductions, where one has a generalized CD formula.

## References

  * [1] K. Abadir and J. Magnus, _Matrix Algebra,_ Cambridge University Press, Cambridge 2005.
  * [2] M. Adler and P. van Moerbeke, _Group factorization, moment matrices and Toda lattices_ , Int. Math. Res. Notices 12 (1997) 556-572.
  * [3] M. Adler and P. van Moerbeke, _Generalized orthogonal polynomials, discrete KP and Riemann–Hilbert problems_ , Commun. Math. Phys. 207 (1999) 589-620.
  * [4] M. Adler and P. van Moerbeke, _Vertex operator solutions to the discrete KP hierarchy_ , Commun. Math. Phys. 203 (1999) 185-210.
  * [5] M. Adler and P. van Moerbekee, _The spectrum of coupled random matrices_ , Ann. Math. 149 (1999) 921-976.
  * [6] M. Adler and P. van Moerbeke, _Hermitian, symmetric and symplectic random ensembles: PDEs for the distribution of the spectrum_ , Ann. Math. 153 (2001) 149-189
  * [7] M. Adler and P. van Moerbeke, _Darboux transforms on band matrices, weights and associated polynomials_ , Int. Math. Res. Notices 18 (2001) 935-984.
  * [8] M. Adler, P. van Moerbeke, and P. Vanhaecke, _Moment matrices and multi- component KP, with applications to random matrix theory_ , Commun. Math. Phys. 286 (2009) 1-38.
  * [9] N. I. Akhiezer, _Classical moment problem and some related questions in analysis,_ Translated by N. Kemmer Hafner Publishing Co., New York 1965.
  * [10] C. Álvarez-Fernández, M. Mañas and U. Fidalgo Prieto, _The multicomponent 2D Toda hierarchy: generalized matrix orthogonal polynomials, multiple orthogonal polynomials and Riemann–Hilbert problems_ , Inverse Probl. 26 (2010) 055009 (15pp).
  * [11] C. Álvarez-Fernández, U. Fidalgo Prieto, and M. Mañas, _Multiple orthogonal polynomials of mixed type: Gauss-Borel factorization and the multi-component 2D Toda hierarchy_ , Adv. Math. 227 (2011) 1451-1525.
  * [12] C. Álvarez-Fernández and M. Mañas, _On the Christoffel–Darboux formula for generalized matrix orthogonal polynomials_ , J. Math. Anal. Appl. 418 (2014) 238-247.
  * [13] R. Álvarez-Nodarse, A. J. Durán, and A. M. de los Ríos, _Orthogonal matrix polynomials satisfying second order difference equation,_ J. Approx. Theory 169 (2013) 40-55.
  * [14] A. I. Aptekarev and E. M. Nikishin, _The scattering problem for a discrete Sturm-Liouville operator_ , Mat. Sb 121 (163): 327-358 (1983); Math. USSR Sb. 49 (1984) 325-355.
  * [15] G. Ariznabarreta and M. Mañas, _Matrix orthogonal Laurent polynomials on the unit circle and Toda type integrable systems,_ Adv. Math. 264 (2014) 396-463.
  * [16] G. Ariznabarreta and M. Mañas, _Multivariate orthogonal polynomials and integrable systems_ , arXiv:1409.0570
  * [17] G. Ariznabarreta and M. Mañas, _Darboux transformations for multivariate orthogonal polynomials,_ arXiv:1503.04786.
  * [18] G. Ariznabarreta and M. Mañas, _Multivariate orthogonal Laurent polynomials and integrable systems,_ arXiv:1506.08708.
  * [19] G. Ariznabarreta and M. Mañas, _Linear spectral transformations for multivariate orthogonal polynomials and multispectral Toda hierarchies_ arXiv:1511.09129.
  * [20] C. Berg, _Fibonacci numbers and orthogonal polynomials_ , J. Comput. Appl. Math. 17 (2011) 75–88.
  * [21] M. J. Bergvelt and A. P. E. ten Kroode, _Tau-functions and zero-curvature equations of Toda-AKNS type_ , J. Math. Phys. 29 (1988) 1308-1320.
  * [22] M. J. Bergvelt and A. P. E. ten Kroode, _Partitions, Vertex Operators Constructions and Multi-Component KP Equations_ , Pac. J. Math. 171 (1995) 23-88.
  * [23] C. Brezinski, _Padé-type approximation and general orthogonal polynomials_. Int. S. Num. M. 50. Birkhäuser Verlag, Basel-Boston, Mass. 1980.
  * [24] M. I. Bueno and F. Marcellán, _Darboux transformation and perturbation of linear functionals_ , Linear Algebra Appl. 384 (2004) 215-242.
  * [25] M. I. Bueno and F. Marcellán, _Polynomial perturbations of bilinear functionals and Hessenberg matrices_ , Linear Algebra Appl.  414 (2006) 64-83.
  * [26] M. J. Cantero, F. Marcellán, L. Moral, and L. Velázquez, _Darboux transformations for CMV matrices_ , arXiv:1503.05003.
  * [27] M. M. Castro and F. A. Grünbaum, _Orthogonal matrix polynomials satisfying first order differential equations: a collection of instructive examples,_ J. Nonlinear Math. Phy. 4 (2005) 63-76.
  * [28] V. S. Chelyshkov, _Alternative Jacobi polynomials and orthogonal exponentials_ , arXiv:1105.1838.
  * [29] T. S. Chihara, _An Introduction to Orthogonal Polynomials_. In :_Mathematics and its Applications Series_ , Vol. 13. Gordon and Breach Science Publishers, New York-London-Paris, 1978.
  * [30] A. E. Choque-Rivero and L. E. Garza, _Moment perturbation of matrix polynomials,_ Integr. Transf. Spec. F. 26 (2015) 177-191.
  * [31] A. R. Collar, _On the reciprocation of certain matrices_ , P. Roy. Soc. Edinb. A 59 (1939) 195-206.
  * [32] E. B. Christoffel, _Über die Gaussische Quadratur und eine Verallgemeinerung derselben_ , J. Reine Angew. Math. 55 (1858) 61–82.
  * [33] D. Damanik, A. Pushnitski, and B. Simon, _The analytic theory of matrix orthogonal polynomials_ , Surveys in Approximation Theory 4 (2008) 1-85.
  * [34] G. Darboux, _Mémoire sur l’approximation des fonctions de très-grands nombres, et sur une classe étendue de développements en série_ , J. Math. Pur. Appl. (3) 4 (1878) 5–56; 377–416.
  * [35] G. Darboux,_Sur une proposition relative aux équations linéaires_ , Comptes Rendus hebd. séances Acad. Sci. Paris 94 (1882) 1456-1459.
  * [36] E. Date, M. Jimbo, M. Kashiwara, and T. Miwa, _Operator approach to the Kadomtsev-Petviashvili equation. Transformation groups for soliton equations. III_ , J. Phys. Soc. Jpn. 50 (1981) 3806-3812.
  * [37] E. Date, M. Jimbo, M. Kashiwara, and T. Miwa, _Transformation groups for soliton equations. Euclidean Lie algebras and reduction of the KP hierarchy_ , Publ. Res. I. Math. Sci. 18 (1982) 1077-1110.
  * [38] E. Date, M. Jimbo, M. Kashiwara, and T. Miwa, _Transformation groups for soliton equations_ in _Nonlinear Integrable Systems-Classical Theory and Quantum Theory_ M. Jimbo and T. Miwa (eds.). World Scientific, Singapore, 1983.
  * [39] A. M. Delgado, J. S. Geronimo, P. Iliev, and F. Marcellán, _Two variable orthogonal polynomials and structured matrices,_ SIAM J. Matrix Anal. A. 28 (2006) 118-147.
  * [40] M. Derevyagin, J. C. García-Ardila, and F. Marcellán, _Multiple Geronimus transformations_ , Linear Algebra Appl. 454 (2014) 158-183.
  * [41] A. Doliwa, P. M. Santini, and M. Mañas, _Transformations of Quadrilateral Lattices_ , J. Math. Phys. 41 (2000) 944-990.
  * [42] A. J. Durán _Markov’s theorem for orthogonal matrix polynomials,_ Canad. J. Math. 48 (1996) 1180-1195.
  * [43] A. J. Durán, _On orthogonal polynomials with respect to a definite positive matrix of measures_ , Canad. J. Math. 47 (1995) 88-112.
  * [44] A. J. Durán, _Matrix inner product having a matrix symmetric second order differential operator,_ Rocky Mt. J. Math. 27 (1997) 585-600.
  * [45] A. J. Durán and F. A. Grünbaum, _A survey on orthogonal matrix polynomials satisfying second order differential equations,_ J. Comput. Appl. Math. 178 (2005) 169-190.
  * [46] A. J. Durán and F. A. Grünbaum, _Structural formulas for orthogonal matrix polynomials satisfying second order differential equations,_ Constr. Approx. 22 (2005) 255-271.