Abstract

We study the direct and inverse spectral theory for a class of finite Hermitian banded matrices. Using the theory of matrix orthogonal polynomials, we provide an explicit procedure for reconstructing a banded matrix from a matrix-valued measure that encodes its spectral data. We establish necessary and sufficient conditions for a measure to be the spectral measure of a matrix in the examined class. We further analyze the connections between this spectral analysis, block tridiagonalization algorithms, and the Toda lattice evolution on banded matrices.

1 Introduction↩︎

Spectral analysis of Jacobi matrices plays a fundamental role in many areas of mathematics. In mathematical physics, inverse spectral theory is used to study integrable nonlinear systems such as the Toda lattice, where the equations of motion are expressed as an isospectral deformation of a Jacobi matrix. In numerical linear algebra, the bijection between Jacobi matrices and spectral measures provides a rigorous framework for analyzing Krylov methods, such as the Lanczos algorithm[1], [2] and conjugate gradient method[3]–[6], using the theory of orthogonal polynomials. In this paper, we aim to examine the spectral and inverse spectral analysis of a broader class of Hermitian banded matrices. We show that working with general bandwidths gives rise to spectral measures that are matrix-valued, connecting naturally to the theory of matrix orthogonal polynomials. Our objective is to fully characterize the connection between these objects by extending the classical strategies used in the tridiagonal case.

Explicitly, we study the spectral theory of finite \(N\times N\) matrices of the form \[\begin{align} \label{eq:BandedMatDef} \mathbf{J}=\begin{bmatrix} \mathbf{A}_0 & \mathbf{B}_0^* & & & & \\ \mathbf{B}_0 & \mathbf{A}_1 & \mathbf{B}_1^* & & & \\ & \mathbf{B}_1 & \mathbf{A}_2 & \mathbf{B}_2^* & & \\ & & \mathbf{B}_2 & \mathbf{A}_3 & \ddots & \\ & & & \ddots & \ddots & \mathbf{B}_{n-2}^* \\ & & & & \mathbf{B}_{n-2} & \mathbf{A}_{n-1} \end{bmatrix}, \end{align}\tag{1}\] where, for each \(j\), \(\mathbf{A}_j=\mathbf{A}_j(\mathbf{J})\) is Hermitian and each \(\mathbf{B}_j(\mathbf{J})\) has full rank. We assume that \[\label{eq:BlocksDim} \begin{align} \mathbf{A}_j(\mathbf{J}) \in \mathbb{C}^{k\times k} ~ \text{for} ~ j=0,1,\dots,n-2, \quad \mathbf{A}_{n-1}(\mathbf{J}) \in \mathbb{C}^{(k-\ell)\times(k-\ell)},\\ \mathbf{B}_j(\mathbf{J})\in \mathbb{C}^{k\times k} ~\text{for}~ j=0,1,\dots,n-3, \quad \mathbf{B}_{n-2}(\mathbf{J})\in \mathbb{C}^{(k-\ell) \times k}, \end{align}\tag{2}\] so that \(N = nk - \ell\) with \(0\leq\ell< k\). We further assume that \(\mathbf{B}_j(\mathbf{J})\) is in row echelon form with positive pivots. The precise class of matrices satisfying these conditions is formalized in Definition 2.

We define a spectral map that assigns to each \(\mathbf{J}\) a \(k\times k\) matrix-valued measure constructed from the eigenvalues of \(\mathbf{J}\) and the first \(k\) components of its normalized eigenvectors. Our main results provide a complete characterization of this map, showing that the associated (matrix) inner product is nondegenerate for polynomials of degree at most \(n-2\) and identifying the exact rank deficiency for those of degree \(n-1\) (Theorems 7 and 8). We further show that this spectral map is injective (Theorem 10) and give a reconstruction procedure for the inverse spectral problem using the theory of matrix orthogonal polynomials, proving that each banded matrix in this class is uniquely determined by its spectral measure (Theorem 11). This inverse spectral analysis also allows us to characterize the range of the spectral map and to formulate necessary and sufficient conditions (Corollary 2) for a measure to correspond to a matrix of the form 1 .

Although techniques and results for related classes of matrices have appeared in the literature (see Section 1.1), to our knowledge, no prior work has leveraged the connection between banded matrices and matrix orthogonal polynomials for finite matrices, particularly in cases where the final blocks are smaller than the preceding ones. An important objective of this work is to also explore the implications of these results in other applications, including the equivalence of block tridiagonalization algorithms and the evolution of the Toda lattice on banded matrices. These connections have not been thoroughly investigated and provide a natural generalization of the classical connections between Jacobi matrices, the Toda flow, and Krylov methods.

The remainder of this paper is organized as follows. In the rest of this section, we review the relevant literature, discuss the connection between banded matrices and block tridiagonalization algorithms such as block Lanczos and the Householder algorithms, and provide background on the Toda lattice, briefly highlighting its generalization to banded matrices. In Section 2, we introduce the notion of matrix-valued measures, examine their properties, and define matrix orthogonal polynomials along with their recurrence relations. Section 3 is devoted to the spectral analysis of matrices of the form 1 . We define the spectral map, characterize its range, and establish its injectivity. We also develop the inverse spectral theory and present an approach to recovering banded matrices from their spectral measure. Finally, Section 4 provides a detailed discussion of the Toda flow on banded matrices and analyzes the evolution of the spectral measure.

1.1 Related work↩︎

While the spectral analysis of Jacobi matrices [7]–[20] has been studied in great detail, techniques for analyzing the spectral theory of general banded matrices[9], [10], [21]–[31] are less developed, with most progress appearing only in recent years. In [9], [10], [32], [33], matrix orthogonal polynomials are used to establish a connection between matrix-valued measures and finite Hermitian block tridiagonal matrices with equal sized blocks, extending the classical spectral theory for Jacobi matrices. This approach was further developed in [28], where infinite matrices of this structure are analyzed thoroughly. To the best of our knowledge, this idea has not been applied to finite banded Hermitian matrices of the form 1 , where the last block may have a degenerate size. This generalization presents additional challenges, arising from the associated matrix-valued measure inducing only a quasi-inner product and the corresponding orthogonal polynomials potentially having degenerate norms. Adressing these difficulties forms the focus of the present work.

Alternative techniques have been developed to analyze the direct and inverse spectral theory of broader classes of banded matrices. In [25], [26], the authors study finite and infinite real symmetric banded matrices whose off-diagonal entries vanish beyond a certain index. The approach in [25], [26] is based on the linear interpolation theory for vector polynomials introduced in [34]. This interpolation theory enables the extension of the results in [22], [23] from pentadiagonal matrices to matrices of arbitrary bandwidth and leads to a unique reconstruction procedure for the inverse spectral problem, an approach that differs from the orthogonal polynomial methodology used in the present work.

A different approach is presented in [24], where the spectral theory of finite Hermitian block matrices is developed. This book introduces the notion of restricted spectral data on a completely extendable set, which refers to a subset of the spectral data that is sufficient to uniquely reconstruct the corresponding Hermitian matrix. Theorem 8 in [24] gives the necessary and sufficient conditions for such a set to exist and provides an reconstruction procedure. It is worth noting, however, that the results in [24] assume all blocks are of the same size.

The analysis in [27] addresses bounded banded operators that admit a positive bidiagonal factorization after an appropriate shift. This approach generalizes the spectral theorem beyond the setting of self-adjoint or normal operators and builds on ideas from the theory of oscillatory matrices. The methods in this work are closely related to multiple orthogonal polynomials.

In this paper, we extend the matrix orthogonal polynomial approach to finite Hermitian banded matrices whose structure is given by 1 . This class of matrices is covered by the approach in [24] only when all the diagonal and off-diagonal blocks have the same size, so our method provides a more general framework for cases where \(\ell>0\) in 2 . Additionally, our work gives a simplified and basic solution to the inverse spectral problem, drawing inspiration from the spectral analysis of Jacobi matrices and avoiding the complexities of linear interpolation theory in [25], [26] and the theory of multiple orthogonal polynomials in [27].

1.2 Equivalence of block tridiagonalization algorithms↩︎

The block Lanczos algorithm [35], [36] is an iterative procedure for constructing a block tridiagonal approximation of a Hermitian matrix. In its simplest form, it is given by Algorithm 1 in Appendix 5. Given a Hermitian matrix \(\mathbf{A}\in \mathbb{C}^{N\times N}\) and an initial block \(\mathbf{V}\in \mathbb{C}^{N\times k}\), the block Lanczos iteration at step \(n\leq \lceil N/k\rceil\) produces a block tridiagonal matrix \(\mathbf{J}_n\) with the structure of 1 and a sequence of matrices \(\mathbf{V}_1, \dots, \mathbf{V}_n \in \mathbb{C}^{N \times k}\) such that¹ \(\mathbf{V}_i^*\mathbf{V}_j = \delta_{ij} \mathbf{I}_k\) and \[\begin{align} \mathbf{J}_n = \mathbf{Q}_n^*\mathbf{A} \mathbf{Q}_n, \qquad \mathbf{Q}_n = \begin{bmatrix} \mathbf{V}_1 & \cdots & \mathbf{V}_n \end{bmatrix}. \end{align}\] The columns of \(\mathbf{Q}_n\) form an orthonormal basis for the degree \(n\) block Krylov subspace defined as \[\begin{align} \mathcal{K}_{n+1}(\mathbf{A},\mathbf{V}) := \mathrm{span}\left\{\mathbf{V}, \mathbf{A}\mathbf{V}, \dots, \mathbf{A}^n \mathbf{V} \right\}. \end{align}\] Here, the span is interpreted as the span of all columns of the matrices \(\mathbf{V},\mathbf{A}\mathbf{V},\dots,\mathbf{A}^n \mathbf{V}\). It is important to note that if \[\begin{align} \mathbf{K}(\mathbf{A}) := \begin{bmatrix} \mathbf{I}_{N\times k},& \mathbf{A} \mathbf{I}_{N\times k},& \dots,& \mathbf{A}^{\lceil N/k\rceil}\mathbf{I}_{N\times k} \end{bmatrix} \end{align}\] has full rank, then block Lanczos on \(\mathbf{A}\) with starting block \(\mathbf{I}_{N\times k}\) runs for \(\lceil N/k \rceil\) iterations and is said to run to completion. On the other hand, if \(n_0\) is the first index for which \(\mathcal{K}_{n_0}(\mathbf{A},\mathbf{V})=\mathcal{K}_{n_0+1}(\mathbf{A},\mathbf{V})\), then the algorithm terminates early at step \(n_0\), as further iterations no longer generate linearly independent vectors.

Block tridiagonalization can also be achieved by modifying the classical Householder tridiagonalization algorithm using a different elimination pattern. For a vector \(\mathbf{v}=[v_1^*,\dots,v_n^*]^*\in\mathbb{C}^n\), the associated Householder reflector is defined by \[\begin{align} \mathbf{H}(\mathbf{v}) = \mathbf{I}_n - 2 \mathbf{w} \mathbf{w}^*, \qquad \mathbf{w} = \frac{\mathbf{u}}{\|\mathbf{u}\|}, \qquad \mathbf{u} = \frac{|v_1|}{v_1}\|\mathbf{v}\|\mathbf{e}_1 + \mathbf{v}, \end{align}\] with the convention that \(|0|/0=1\). The Householder procedure uses these reflectors to construct a sequence of unitary matrices \(\{\mathbf{Q}_j\}_{j=1}^N\) that sequentially bring the leading \(n\times n\) principal subblocks of \[\begin{align} \label{eq:HouseholderTransf} \mathbf{Q}_{n}\cdots \mathbf{Q}_1 \mathbf{A} \mathbf{Q}_1^* \cdots \mathbf{Q}_{n}^* \end{align}\tag{3}\] to the desired block tridiagonal form. At step \(n\), the matrix \(\mathbf{Q}_n\) is constructed as \[\begin{align} \mathbf{Q}_n = \begin{bmatrix} \mathbf{I}_{k+n-1} & \mathbf{0} \\ \mathbf{0} & \mathbf{H}_n \end{bmatrix}, \qquad \mathbf{H}_n= \mathbf{H}(\mathbf{v}_n), \end{align}\] where \(\mathbf{v}_n\) is the vector formed by the last \(N-k-n+1\) entries of the \(n\)-th column of 3 . When \(\mathbf{Q}_n^*\) is applied on the left of 3 , it introduces zeros into rows \(k+n+1,\dots,N\) of the \(n\)-th column. Multiplication on the right by \(\mathbf{Q}_n\) preserves these zeros, and by symmetry the corresponding entries in the \(n\)-th row are also zero, yielding the desired block tridiagonal form.

A brief review of the Householder reduction to block tridiagonal form is provided in Appendix 5. The algorithm applies Householder reflectors sequentially, one transformation at a time. For improved computational efficiency, especially on modern architectures, several consecutive reflectors can be accumulated and applied simultaneously using block representations [37]–[39]. Such implementations rely on compact representations of products of Householder reflectors, such as the \(\mathbf{W}\mathbf{Y}^*\) and compact \(\mathbf{W}\mathbf{Y}^*\) representations[40], [41]. A related but conceptually distinct approach is based on block Householder reflectors, defined as matrices of the form \[\begin{align} \mathbf{H}(\mathbf{W}) = \mathbf{I}_n - 2 \mathbf{W} \mathbf{W}^* , \end{align}\] where \(\mathbf{W}\in \mathbb{C}^{N\times k}\) satisfies \(\mathbf{W}^* \mathbf{W} = \mathbf{I}_k\). These reflectors generalize classical Householder transformations, but they do not, in general, preserve the triangular structure of the off-diagonal blocks described in 1 . For further details, see [42].

The following result establishes an equivalence between block Lanczos and the Householder reduction. Although this equivalence is classical [43], we present a simple alternative argument based on the spectral theory of banded Hermitian matrices.

Theorem 1. If block Lanczos (Algorithm 1) applied to a Hermitian matrix \(\mathbf{A}\in \mathbb{C}^{N\times N}\) with starting block \(\mathbf{I}_{N\times k}\) runs to completion, then it produces the same block tridiagonal matrix as the Householder procedure (Algorithm 2).

Proof. Let \(\mathbf{J}_1\) and \(\mathbf{J}_2\) denote the block tridiagonal matrices produced by the block Lanczos algorithm and the Householder procedure, respectively, and write \[\begin{align} \label{eq:SimilarityEqsLanczosHouseholder} \mathbf{J}_1 = \mathbf{Q}_1^* \mathbf{A} \mathbf{Q}_1, \quad \mathbf{J}_2 = \mathbf{Q}_2^* \mathbf{A} \mathbf{Q}_2, \end{align}\tag{4}\] where \(\mathbf{Q}_1,\mathbf{Q}_2\) are unitary matrices whose first \(k\) columns are \(\mathbf{I}_{N\times k}\). Consider the spectral map \(\varphi\) defined in 24 below, and let \({\boldsymbol{\mu}}_1 = \varphi(\mathbf{J}_1)\) and \({\boldsymbol{\mu}}_2 = \varphi(\mathbf{J}_2)\). Since \(\mathbf{J}_1\) and \(\mathbf{J}_1\) have the same eigenvalues and identical first \(k\) eigenvector entries, it follows that \({\boldsymbol{\mu}}_1 = {\boldsymbol{\mu}}_2\). Injectivity of \(\varphi\) in Theorem 10 therefore gives \(\mathbf{J}_1=\mathbf{J}_2\). ◻

1.3 Toda flow↩︎

The finite Toda lattice, introduced by Morikazu Toda in 1967 [44], [45], is a completly integrable model for a nonlinear one-dimensional crystal. The integrability of the Toda lattice was established independently by Flaschka[46] and Manakov[47], who proved that the system can be written as \[\begin{align} \label{eq:toda} \partial_t \mathbf{X} = [\mathbf{X},\mathbf{B}(\mathbf{X})] = \mathbf{X}\mathbf{B}(\mathbf{X})-\mathbf{B}(\mathbf{X})\mathbf{X}, \quad \mathbf{X}(0)=\mathbf{X}_0, \end{align}\tag{5}\] where \(\mathbf{X}_0 \in \mathbb{C}^{N\times N}\) is a Jacobi matrix. The matrix \(\mathbf{B}(\mathbf{X})\) is defined by \[\begin{align} \label{eq:Bdef} \mathbf{B}(\mathbf{X}) = \mathbf{X}_- - \mathbf{X}_-^T, \end{align}\tag{6}\] with \(\mathbf{X}_-\) denoting the strictly lower triangular part of \(\mathbf{X}\).

Under the Toda flow, the eigenvalues of \(\mathbf{X}\) remain constant in time, and the first entries of the eigenvectors evolve in a simple way. Let \(\lambda_j(t)\) be the eigenvalues of \(\mathbf{X}(t)\) and \(\mathbf{v}_{1,j}(t)\) be the first entry of the \(j\)-th normalized eigenvector, and assume without loss of generality that ² \(\mathbf{v}_{1,j}(t)>0\). For each \(j = 1, \dots, N\), we have \[\begin{align} \label{eq:JacEvalEvecEvol} \lambda_j(t) = \lambda_j, \quad \text{and} \quad \mathbf{v}_{1,j}(t) = \frac{\mathbf{v}_{1,j}(0) e^{\lambda_j t}}{\left( \sum_{i=1}^N \mathbf{v}_{1,i}^2(0), e^{2\lambda_i t} \right)^{1/2}}, \end{align}\tag{7}\] where \(\lambda_1\geq \lambda_2\geq \dots\geq \lambda_N\) are the eigenvalues of \(\mathbf{X}_0\). This leads to a simple expression for the spectral measure of \(\mathbf{X}(t)\), given by \[\begin{align} \label{eq:JacSpectralMeasure} \mu_{\mathbf{X}(t)} = \sum_{j=1}^N w_j(t)\delta_{\lambda_j}, \qquad w_j(t) = |\mathbf{v}_{1,j}(t)|^2 = \frac{e^{2\lambda_j t} w_j(0)}{\sum_{j=1}^N e^{2\lambda_i t} w_i(0)} \quad \text{for }j=1,\dots,N. \end{align}\tag{8}\]

The evolution of \(\mu_{\mathbf{X}(t)}\) provides a direct procedure for solving the finite Toda lattice using inverse spectral methods. Starting from an initial Jacobi matrix \(\mathbf{X}_0\), one computes its eigenvalues and the first components of its eigenvectors, evolves the weights using an explicit exponential factor, and reconstructs the Jacobi matrix from the spectral measure. This solution process is remarkable because, despite the Toda lattice being nonlinear, its integrable structure provides spectral variables in which the evolution is explicit and simple.

There are limited results giving analogous formulae for the Toda lattice with banded matrix initial data. Results for pentadiagonal initial data appear in [48], where the evolution of the first two eigenvector components is derived. However, to the best of our knowledge, an explicit characterization of Toda lattice in terms of the evolution of a matrix-valued measure has not been obtained for general bandwidth. Using the spectral theory developed in Sections 3, we show that the Toda flow preserves the banded structure of matrices of the form 1 and induces a simple evolution on the associated matrix-valued spectral measure. Specifically, the spectral measure evolves as \[\begin{align} {\boldsymbol{\mu}}_{\mathbf{X}(t)} = \sum_{j=1}^N \mathbf{L}^{-1}(t) \big( e^{2\lambda_j t} \mathbf{v}_j(0) \mathbf{v}_j(0)^* \big) \mathbf{L}^{-*}(t) \, \delta_{\lambda_j}, \end{align}\] where \(\{\lambda_j\}_{j=1}^N\) are the eigenvalues of \(\mathbf{X}_0\), \(\mathbf{v}_j(0)\) denotes the first \(k\) components of the corresponding eigenvector, and \(k\) is the bandwidth or block size of \(\mathbf{X}_0\). The lower triangular matrix \(\mathbf{L}(t)\) has positive diagonal entries and satisfies \[\begin{align} \sum_{j=1}^N e^{2\lambda_j t} \mathbf{v}_j(0)\mathbf{v}_j(0)^* = \mathbf{L}(t) \mathbf{L}(t)^*, \end{align}\] and thus serves to give the Cholesky factorization of the sum. This shows that the Toda flow is also fully determined by the spectral data for banded matrices. This result is presented in Section 4.

2 Matrix-valued measures and polynomials↩︎

A (positive) matrix-valued measure \({\boldsymbol{\mu}}\) on \(\mathbb{R}\) is a countably-additive set function \[\begin{align} {\boldsymbol{\mu}}:\mathcal{B}(\mathbb{R}) \to \mathbb{C}^{k \times k}_+, \end{align}\] where \(\mathcal{B}(\mathbb{R})\) denotes the Borel \(\sigma\)-algebra on \(\mathbb{R}\) and \(\mathbb{C}^{k \times k}_+\) is the set of Hermitian positive semi-definite \(k\times k\) matrices. The measure \({\boldsymbol{\mu}}\) is said to be normalizable if \(\|{\boldsymbol{\mu}}(\mathbb{R})\| < \infty\) and \(\det {\boldsymbol{\mu}}(\mathbb{R}) \neq 0\). In this case, we assume³ \[\begin{align} {\boldsymbol{\mu}}(\mathbb{R}) = \mathbf{I}_k, \end{align}\] and we write \({\boldsymbol{\mu}}\in \mathcal{P}_k(\mathbb{R})\). A point \(x\) is said to be in the support of \({\boldsymbol{\mu}}\) if \({\boldsymbol{\mu}}(B) \neq \mathbf{0}\) for every open set \(B\) that contains \(x\). The set of all such points is denoted by \(\mathop{\mathrm{supp}}{\boldsymbol{\mu}}\).

A matrix-valued measure \({\boldsymbol{\mu}}\) induces both left and right quasi-inner products on matrix-valued Borel-measurable functions \(\mathbf{F},\mathbf{G}\colon\mathbb{R}\to\mathbb{C}^{k\times k}\), defined, respectively, by \[\begin{align} \langle \mathbf{F},\mathbf{G}\rangle_L := \int_{\mathbb{R}}\mathbf{F}(x){\boldsymbol{\mu}}(\mathrm d x)\mathbf{G}^*(x), \quad \langle \mathbf{F},\mathbf{G}\rangle_R := \int_{\mathbb{R}}\mathbf{F}^*(x){\boldsymbol{\mu}}(\mathrm d x)\mathbf{G}(x). \end{align}\] These products are related via \(\langle \mathbf{F}, \mathbf{G} \rangle_L = \langle \mathbf{F}^*, \mathbf{G}^* \rangle_R\). Throughout the paper, we work mainly with the right quasi-inner product and, with slight abuse of notation, write \(\langle \mathbf{F}, \mathbf{G} \rangle_{\boldsymbol{\mu}}\) in place of \(\langle \mathbf{F}, \mathbf{G} \rangle_R\). The right product satisfies the following properties:

\(\langle \mathbf{F}, \mathbf{G} \rangle_{\boldsymbol{\mu}}= \langle \mathbf{G}, \mathbf{F} \rangle_{\boldsymbol{\mu}}^*\),
For \(\mathbf{C}\in \mathbb{C}^{k\times k}\), \(\langle \mathbf{F},\mathbf{G} \mathbf{C}\rangle_{\boldsymbol{\mu}}= \langle \mathbf{F},\mathbf{G} \rangle_{\boldsymbol{\mu}}\mathbf{C}\) and \(\langle\mathbf{F} \mathbf{C},\mathbf{G}\rangle_{\boldsymbol{\mu}}= \mathbf{C}^* \langle \mathbf{F},\mathbf{G} \rangle_{\boldsymbol{\mu}}\),
\(\langle \mathbf{F}, \mathbf{F} \rangle_{\boldsymbol{\mu}}\in \mathbb{C}^{k\times k}_+\),
\(\| \mathbf{F} \| := (\mathop{\mathrm{tr}}\langle \mathbf{F}, \mathbf{F} \rangle_{\boldsymbol{\mu}})^{1/2}\) defines a semi-norm.

We note that \(\|\mathbf{F}\|\) cannot be expected to be a norm without additional assumptions. For example, consider \[\begin{align} {\boldsymbol{\mu}}= \sum_{i = 1}^n \mathbf{w}_i \mathbf{w}_i^* \delta_{x_i}, \end{align}\] for \(x_i < x_{i+1}\) and \(\sum_i \mathbf{w}_i \mathbf{w}_i^* = \mathbf{I}_k\). Then for any function \(\mathbf{F}\) satisfying \(\mathbf{F}(x_i) \mathbf{w}_i = \mathbf{0}\), we have \(\langle \mathbf{F}, \mathbf{F} \rangle_{\boldsymbol{\mu}}= \mathbf{0}\). This motivates the notion of \(n\)-definiteness, which guarantees that a sequence of \(n\) orthonormal polynomials is well-defined (see Section 2.2).

Definition 1. A measure \({\boldsymbol{\mu}}\in \mathcal{P}_k(\mathbb{R})\) is called \(n\)-definite if \(\|\mathbf{P}\| \neq 0\) for every matrix-valued polynomial \(\mathbf{P}\) of degree at most \(n-2\) satisfying \(\mathbf{P}(x) \neq \mathbf{0}\) for some \(x \in \mathbb{R}\).

2.1 Monic Orthogonal Polynomials↩︎

Consider a sequence of monic matrix-valued polynomials \(\mathbf{\Pi}_0(x),\mathbf{\Pi}_1(x),\mathbf{\Pi}_2(x),\dots\) defined by \[\begin{align} \label{eq:MonicDef} \mathbf{\Pi}_j(x) = \mathbf{\Pi}_j(x,{\boldsymbol{\mu}}) = \mathbf{I}_k~x^j + {\boldsymbol{\tau}}_j({\boldsymbol{\mu}}) x^{j-1} + \text{ lower order terms}, \quad {\boldsymbol{\tau}}_j:={\boldsymbol{\tau}}_j({\boldsymbol{\mu}})\in \mathbb{C}^{k\times k}. \end{align}\tag{9}\] These polynomials are said to be orthogonal with respect to the quasi-inner product \(\langle \cdot,\cdot \rangle_{\boldsymbol{\mu}}\) if⁴ \[\begin{align} \label{eq:MonicOrthogCond} \langle \mathbf{\Pi}_i(\diamond,{\boldsymbol{\mu}}),\mathbf{\Pi}_j(\diamond,{\boldsymbol{\mu}})\rangle_{\boldsymbol{\mu}}= {\boldsymbol{\gamma}}_j \delta_{ij}, \quad \text{where {\boldsymbol{\gamma}}_j:= {\boldsymbol{\gamma}}_j({\boldsymbol{\mu}})=\langle \mathbf{\Pi}_j(\diamond,{\boldsymbol{\mu}}),\mathbf{\Pi}_j(\diamond,{\boldsymbol{\mu}}) \rangle_{\boldsymbol{\mu}}\in \mathbb{C}_+^{k\times k}}. \end{align}\tag{10}\] We now show that, when \({\boldsymbol{\mu}}\) is \(n\)-definite, monic orthogonal polynomials of degree at most \(n-1\) are uniquely determined, and therefore \(\{{\boldsymbol{\gamma}}_j({\boldsymbol{\mu}})\}_{j=0}^{n-1}\) and \(\{\tau_j({\boldsymbol{\mu}})\}_{j=0}^{n-1}\) are well-defined. To establish this, we first note that any polynomial can be expressed as a linear combination of monic polynomials.

Lemma 1. Let \(\{\mathbf{\Pi}_j\}_{j\geq0}\) be a sequence of monic matrix-valued polynomials. For every polynomial \(\mathbf{P}\in \mathbb{C}^{k\times k}\) of degree \(p\), there exists unique matrices \(\{\Phi_j\}_{j=0}^{p}\) with \(\Phi_j \in \mathbb{C}^{k\times k}\) such that \(\mathbf{P}(x) =\sum_{j=0}^{p}\mathbf{\Pi}_j(x) \Phi_j\).

Theorem 2. If the \({\boldsymbol{\mu}}\) is \(n\)-definite, then the monic orthogonal polynomials \(\{\mathbf{\Pi}_j(\diamond,{\boldsymbol{\mu}})\}_{j=0}^{n-1}\) are uniquely defined and \(\det {\boldsymbol{\gamma}}_j({\boldsymbol{\mu}}) \neq 0\) for \(j \leq n-2\).

Proof. Existence of the monic orthogonal polynomials \(\mathbf{\Pi}_0,\dots, \mathbf{\Pi}_{n-1}\) follows recursively by applying the Gram-Schmidt process to the sequence (\(x\mapsto\mathbf{I}_k\), \(x\mapsto x\mathbf{I}_k\), \(\dots\), \(x\mapsto x^{n-1}\mathbf{I}_k\)), provided that each \({\boldsymbol{\gamma}}_j\) is nonsingular for \(1\leq j<n-1\).

To show \(\mathrm{det}~{\boldsymbol{\gamma}}_{j}\neq 0\) for \(j<n-1\), suppose for contradiction that \(j\) is the first index where \({\boldsymbol{\gamma}}_j\) is singular. In this case, there exists a nonzero vector \(\mathbf{v}\) such that \({\boldsymbol{\gamma}}_j \mathbf{v} = \mathbf{0}\). Define \(\mathbf{\Pi} = \mathbf{\Pi}_j\mathbf{v}\mathbf{v}^*\), then \(\mathbf{\Pi}\neq \mathbf{0}\) since its leading coefficient is \(\mathbf{v}\mathbf{v}^*\neq \mathbf{0}\). However, we have \(\langle \mathbf{\Pi},\mathbf{\Pi} \rangle_{\boldsymbol{\mu}}= \mathbf{v}\mathbf{v}^*{\boldsymbol{\gamma}}_j \mathbf{v}\mathbf{v}^* = \mathbf{0}\), which contradicts the fact that \({\boldsymbol{\mu}}\) is \(n\)-definite, and thus \(\det {\boldsymbol{\gamma}}_j \neq 0\).

To establish uniqueness, let \(\tilde{\mathbf{\Pi}}_j\) be another monic orthogonal polynomial of degree \(j\). Note that \(\mathbf{\Pi}_j-\tilde{\mathbf{\Pi}}_j\) is a polynomial of degree of \(j-1\), and hence \[\begin{align} \langle \mathbf{\Pi}_j-\tilde{\mathbf{\Pi}}_j, \mathbf{\Pi}_j-\tilde{\mathbf{\Pi}}_j\rangle_{\boldsymbol{\mu}}= \langle \mathbf{\Pi}_j , \mathbf{\Pi}_j-\tilde{\mathbf{\Pi}}_j \rangle_{\boldsymbol{\mu}}- \langle \tilde{\mathbf{\Pi}}_j , \mathbf{\Pi}_j-\tilde{\mathbf{\Pi}}_j \rangle_{\boldsymbol{\mu}}= \mathbf{0}, \end{align}\] where the second equality follows from Lemma 1. Finally, since \({\boldsymbol{\mu}}\) is \(n\)-definite, it follows that \(\tilde{\mathbf{\Pi}}_j = \mathbf{\Pi}_j\). ◻

The monic orthogonal polynomials satisfy a number of well-known properties (see, for instance, [28], [49]–[53]). We highlight below a fundamental result that will be important for our work.

Lemma 2. If \({\boldsymbol{\mu}}\) is \(n\)-definite and \(\mathbf{P}\) has degree \(p\leq n-2\), then the coefficients \(\{\Phi_j\}_{j=0}^p\) in Lemma 1 are given by \(\Phi_j = {\boldsymbol{\gamma}}_j^{-1}({\boldsymbol{\mu}}) \langle \mathbf{\Pi}_j,\mathbf{P} \rangle_{\boldsymbol{\mu}}\).

Proof. The statement follows from Lemma 1, the orthogonality relations in 10 , and Theorem 2, which ensures that \({\boldsymbol{\gamma}}_j\) is invertible for \(j\leq n-2\). ◻

Theorem 3. Assume that \({\boldsymbol{\mu}}\) is \(n\)-definite. Then the associated monic orthogonal matrix polynomials obey the following three-term recurrence relation, \[\begin{align} \label{eq:MonicRec} x\mathbf{\Pi}_j(x) = \mathbf{\Pi}_{j+1}(x) + \mathbf{\Pi}_j(x) \mathbf{C}_j + \mathbf{\Pi}_{j-1}(x) \mathbf{D}_j , \quad \text{for } j=0,1,\dots,n-2 \end{align}\tag{11}\] where \(\mathbf{\Pi}_{-1}\equiv 0\), \({\boldsymbol{\gamma}}_{-1} = \mathbf{I}_k\), \[\begin{align} \mathbf{C}_j({\boldsymbol{\mu}}) = {\boldsymbol{\gamma}}_{j}^{-1}\langle \diamond\mathbf{\Pi}_j(\diamond),\mathbf{\Pi}_j \rangle_{\boldsymbol{\mu}}= {\boldsymbol{\tau}}_j({\boldsymbol{\mu}})-{\boldsymbol{\tau}}_{j+1}({\boldsymbol{\mu}}), \quad \mathbf{D}_j({\boldsymbol{\mu}}) = {\boldsymbol{\gamma}}_{j-1}^{-1}({\boldsymbol{\mu}})\langle \diamond\mathbf{\Pi}_{j-1}(\diamond),\mathbf{\Pi}_j \rangle_{\boldsymbol{\mu}}= {\boldsymbol{\gamma}}_{j-1}^{-1}({\boldsymbol{\mu}}){\boldsymbol{\gamma}}_j({\boldsymbol{\mu}}). \end{align}\]

Proof. Since \(x\mathbf{\Pi}_j(x)-\mathbf{\Pi}_{j+1}(x)\) is a polynomial of degree \(j\) with leading coefficient \({\boldsymbol{\tau}}_j-{\boldsymbol{\tau}}_{j-1}\), Lemma 1 implies that \[\begin{align} x\mathbf{\Pi}_j(x)-\mathbf{\Pi}_{j+1}(x) = \sum_{i=0}^j \mathbf{\Pi}_i(x) \Phi_i, \quad \text{with} \quad \Phi_j = {\boldsymbol{\tau}}_j-{\boldsymbol{\tau}}_{j+1}. \end{align}\] Taking the inner product of both sides with \(\mathbf{\Pi}_\ell\) gives \[\begin{align} \langle \mathbf{\Pi}_\ell, \diamond\mathbf{\Pi}_j(\diamond)-\mathbf{\Pi}_{j+1} \rangle_{\boldsymbol{\mu}}= \sum_{i=0}^j \langle \mathbf{\Pi}_\ell , \mathbf{\Pi}_i\rangle_{\boldsymbol{\mu}}\Phi_i, = {\boldsymbol{\gamma}}_\ell \Phi_\ell. \end{align}\] However, the left-hand side can be rewritten as \[\begin{align} \langle \mathbf{\Pi}_\ell, \diamond\mathbf{\Pi}_j(\diamond)-\mathbf{\Pi}_{j+1} \rangle_{\boldsymbol{\mu}}= \langle \diamond\mathbf{\Pi}_\ell(\diamond), \mathbf{\Pi}_j \rangle_{\boldsymbol{\mu}}- \langle \mathbf{\Pi}_\ell,\mathbf{\Pi}_{j+1} \rangle_{\boldsymbol{\mu}}. \end{align}\] By the orthogonality condition 10 , both terms vanish when \(\ell\leq j-2\), and hence \(\Phi_\ell = \mathbf{0}\) in this case. Finally, we have \[\begin{align} \Phi_{j-1} = {\boldsymbol{\gamma}}_{j-1}^{-1} \langle \diamond\mathbf{\Pi}_{j-1}(\diamond), \mathbf{\Pi}_j \rangle_{\boldsymbol{\mu}}= {\boldsymbol{\gamma}}_{j-1}^{-1}\left( \langle \mathbf{\Pi}_{j}, \mathbf{\Pi}_j \rangle_{\boldsymbol{\mu}}- \langle \diamond\mathbf{\Pi}_{j-1}(\diamond) - \mathbf{\Pi}_j, \mathbf{\Pi}_j \rangle_{\boldsymbol{\mu}}\right) = {\boldsymbol{\gamma}}_{j-1}^{-1}{\boldsymbol{\gamma}}_j, \end{align}\] which completes our proof. ◻

2.2 Orthonormal Polynomials↩︎

If \({\boldsymbol{\mu}}\) is \(n\)-definite, then a family of orthonormal matrix polynomials \(\{\mathbf{P}_j(\diamond,{\boldsymbol{\mu}})\}_{j=0}^{n-2}\) is given by \[\begin{align} \label{eq:OpsNormalization} \mathbf{P}_j(x) := \mathbf{P}_j(x,{\boldsymbol{\mu}}) = \mathbf{\Pi}_j(x,{\boldsymbol{\mu}}) \mathbf{\boldsymbol{\gamma}}_j^{-1/2}({\boldsymbol{\mu}}) \mathbf{Q}_j, \end{align}\tag{12}\] where \(\mathbf{\Pi}_j(\diamond,{\boldsymbol{\mu}})\) and \({\boldsymbol{\gamma}}_j({\boldsymbol{\mu}})\) are defined in 10 , and \(\mathbf{Q}_j\) is an arbitrary unitary matrix. It is important to note that, for every choice of \(\{\mathbf{Q}_j\}_{j=0}^{n-2}\), we have \(\langle\mathbf{P}_i,\mathbf{P}_j\rangle_{\boldsymbol{\mu}}= \mathbf{I}_k\delta_{ij}\), so orthonormal polynomials are determined only up to a right multiplication by a unitary matrix. Throughout our work, we fix \(\mathbf{Q}_0 = \mathbf{I}_k\) so that \(\mathbf{P}_0 = \mathbf{I}_k\).

Theorem 3 implies that these polynomials follow a Hermitian recurrence relation \[\begin{align} \label{eq:OrthonormRec} x \mathbf{P}_j(x) = \mathbf{P}_{j+1}(x) \mathbf{B}_{j} + \mathbf{P}_j(x)\mathbf{A}_j + \mathbf{P}_{j-1}(x)\mathbf{B}^*_{j-1}, \quad j=0,\dots,n-3, \quad \mathbf{A}_j = \mathbf{A}_j^*, \quad \det \mathbf{B}_j\neq 0, \end{align}\tag{13}\] with the convention \(\mathbf{P}_{-1} \equiv \mathbf{0}\) and \(\mathbf{B}_{-1} \equiv \mathbf{I}_k\). The recurrence coefficients \(\mathbf{A}_j = \mathbf{A}_j({\boldsymbol{\mu}})\) and \(\mathbf{B}_j = \mathbf{B}_j({\boldsymbol{\mu}})\) are explicitly given by \[\begin{align} \label{eq:RecCoefForm} \mathbf{A}_j = \langle \mathbf{P}_j,\diamond \mathbf{P}_j(\diamond) \rangle_{\boldsymbol{\mu}}= \mathbf{Q}_j^*{\boldsymbol{\gamma}}_j^{1/2}\mathbf{C}_j{\boldsymbol{\gamma}}_j^{-1/2}\mathbf{Q}_j, \quad \text{and} \quad \mathbf{B}_j = \langle \mathbf{P}_{j+1}, \diamond \mathbf{P}_j(\diamond) \rangle_{\boldsymbol{\mu}}= \mathbf{Q}_{j+1}^*{\boldsymbol{\gamma}}_{j+1}^{-1/2}\mathbf{D}_{j+1}^* {\boldsymbol{\gamma}}_j^{1/2}\mathbf{Q}_j, \end{align}\tag{14}\] where \(\mathbf{C}_j\) and \(\mathbf{D}_j\) are defined in Theorem 3. Lemma 3 shows that there exists a unique choice of \(\{\mathbf{Q}_j\}_{j=0}^{n-2}\) that produces a three-term recurrence in which \(\mathbf{B}_j\), for \(j=0,\dots,n-3\), is upper triangular with positive diagonal entries.

Lemma 3. Let \({\boldsymbol{\mu}}\) be \(n\)-definite, then there exists a unique family of matrix-valued orthogonal polynomials \(\{\mathbf{P}_j(\diamond,{\boldsymbol{\mu}})\}_{j=0}^{n-2}\), i.e. a unique choice of unitary matrices \(\{\mathbf{Q}_j\}_{j=1}^{n-2}\) in 12 , such that the matrices \(\{\mathbf{A}_j({\boldsymbol{\mu}})\}_{j=0}^{n-3}\) are Hermitian and the matrices \(\{\mathbf{B}_j({\boldsymbol{\mu}})\}_{j=0}^{n-3}\) are upper triangular with positive diagonal entries.

Proof. Let \(\{\mathbf{\tilde{P}}_j\}_{j=0}^{n-2}\) be any family of matrix orthonormal polynomials associated with \({\boldsymbol{\mu}}\), with recurrence coefficients \(\{\mathbf{\tilde{A}}_j\}_{j=0}^{n-3}\) and \(\{\mathbf{\tilde{B}}_j\}_{j=0}^{n-3}\). We first construct a family \(\{\mathbf{P}_j\}\) inductively using unitary transformations to enforce the required conditions on the recurrence coefficients. Consider \[\begin{align} x\mathbf{\tilde{P}}_0(x) = \mathbf{\tilde{P}}_1(x) \mathbf{\tilde{B}}_0 + \mathbf{\tilde{P}}_0(x) \mathbf{\tilde{A}}_0, \end{align}\] and let \(\mathbf{\tilde{B}}_0 = \mathbf{Q}_1 \mathbf{B}_0\) be a QR decomposition, where \(\mathbf{Q}_1\) is unitary and \(\mathbf{B}_0\) is upper triangular with positive diagonal entries. Setting \[\begin{align} \mathbf{P}_0 = \mathbf{\tilde{P}}_0, \quad \mathbf{P}_1 = \mathbf{\tilde{P}}_1 \mathbf{Q}_1, \quad \text{and} \quad \mathbf{A}_0 = \mathbf{\tilde{A}}_0, \end{align}\] we find \[\begin{align} x\mathbf{P}_0(x) = \mathbf{P}_1(x) \mathbf{B}_0 + \mathbf{P}_0(x) \mathbf{A}_0. \end{align}\] Assume now that \(\{\mathbf{P}_j\}_{j=0}^{i}\) have been constructed so that the required conditions on \(\mathbf{A}_j\) and \(\mathbf{B}_j\) hold for all \(j<i\). Let \(\mathbf{Q}_1,\dots,\mathbf{Q}_i\) be the unitary matrices used in the previous steps, and by construction these matrices satisfy \[\begin{align} \mathbf{\tilde{B}}_j\mathbf{Q}_j = \mathbf{Q}_{j+1} \mathbf{B}_j, \quad \text{and} \quad \mathbf{A}_j = \mathbf{Q}_j^* \mathbf{\tilde{A}_j\mathbf{Q}_j} \quad \text{for j=1,\dots,i-1}. \end{align}\] Consider the recurrence for \(\mathbf{\tilde{P}}_{i+1}\), that is, \[\begin{align} x\mathbf{\tilde{P}}_i(x) = \mathbf{\tilde{P}}_{i+1}(x) \mathbf{\tilde{B}}_i + \mathbf{\tilde{P}}_i(x) \mathbf{\tilde{A}}_i + \mathbf{\tilde{P}}_{i-1}(x) \mathbf{\tilde{B}}_{i-1}. \end{align}\] Multiplying on the right by \(\mathbf{Q}_i\) gives \[\begin{align} x\mathbf{P}_i(x) = \mathbf{\tilde{P}}_{i+1}(x) \mathbf{\tilde{B}}_i \mathbf{Q}_{i} + \mathbf{P}_i(x) \mathbf{Q}_i^*\mathbf{\tilde{A}}_i \mathbf{Q}_i + \mathbf{\tilde{P}}_{i-1}\mathbf{Q}_{i-1}^*(x) \mathbf{\tilde{B}}_{i-1} \mathbf{Q}_i. \end{align}\] Using a QR decomposition \(\mathbf{\tilde{B}}_i \mathbf{Q}_i = \mathbf{Q}_{i+1}\mathbf{B}_i\), we define \(\mathbf{P}_{i+1}=\mathbf{\tilde{P}}_{i+1}\mathbf{Q}_{i+1}\) and \(\mathbf{A}_i = \mathbf{Q}_i^* \mathbf{A}_i \mathbf{Q}_i\). This guarantees that the \((i+1)\)-th recurrence has the desired form, establishing the existence of a normalization. Uniqueness follows from the uniqueness of the QR decomposition for invertible matrices. ◻

For the remainder of this paper, we fix the orthonormal polynomials \(\{\mathbf{P}_j(\diamond,{\boldsymbol{\mu}})\}_{j=0}^{n-2}\) to be the unique sequence given by Lemma 3. The \((n-1)\)-th orthonormal polynomial⁵ \(\mathbf{P}_{n-1}(\diamond,{\boldsymbol{\mu}})\) is not directly defined, since the normalization factor \({\boldsymbol{\gamma}}_{n-1}({\boldsymbol{\mu}})=\langle\mathbf{\Pi}_{n-1}(\diamond,{\boldsymbol{\mu}}),\mathbf{\Pi}_{n-1}(\diamond,{\boldsymbol{\mu}}) \rangle\) in 12 is potentially singular. In the following, we construct \(\mathbf{P}_{n-1}(\diamond,{\boldsymbol{\mu}})\) when \(\mathrm{rank}~{\boldsymbol{\gamma}}_{n-1}({\boldsymbol{\mu}}) = k-\ell\) for some \(0\leq \ell<k\).

Lemma 4. Let \(\mathbf{A} \in \mathbb{C}^{n\times n}\) be a positive semidefinite matrix of rank \(k\), then there exists a unique upper triangular matrix \(\mathbf{R}\in \mathbb{C}^{k\times n}\) in row echelon form and positive pivots such that \(\mathbf{A} = \mathbf{R}^* \mathbf{R}\).

Proof. Since \(\mathbf{A}\) is positive semidefinite of rank \(k\), there exists \(\mathbf{X} \in \mathbb{C}^{k\times n}\) such that \(\mathbf{A} = \mathbf{X}^* \mathbf{X}\). Because \(\mathbf{X}\) has full row rank, it can be factored as \(\mathbf{X} = \mathbf{Q} \mathbf{R}\) where \(\mathbf{Q} \in \mathbb{C}^{k\times k}\) is unitary and \(\mathbf{R} \in \mathbb{C}^{k\times n}\) is in row echelon form with positive pivots. Substituting into \(\mathbf{A} = \mathbf{X}^* \mathbf{X}\) gives \(\mathbf{A} = \mathbf{R}^* \mathbf{Q}^* \mathbf{Q} \mathbf{R} = \mathbf{R}^* \mathbf{R}\).

Suppose \(\mathbf{\tilde{R}}\) is another matrix in row echelon form with positive pivots such that \(\mathbf{A} = \mathbf{\tilde{R}}^* \mathbf{\tilde{R}}\). The pivot locations of \(\mathbf{R}\) and \(\mathbf{\tilde{R}}\) must coincide; otherwise, the submatrix of \(\mathbf{R}\) formed by the pivot columns is nonsingular while the corresponding submatrix of \(\mathbf{\tilde{R}}\) is singular. Moreover, defining⁶ \(\mathbf{U}:=\mathbf{R}\tilde{\mathbf{R}}^\dagger\in\mathbb{C}^{k\times k}\), we have \[\begin{align} \mathbf{U}^*\mathbf{U} = (\tilde{\mathbf{R}}^\dagger)^*\mathbf{R}^* \mathbf{R}\tilde{\mathbf{R}}^\dagger = (\tilde{\mathbf{R}} \tilde{\mathbf{R}}^\dagger)^* (\tilde{\mathbf{R}}\tilde{\mathbf{R}}^\dagger) = \mathbf{I}_k, \end{align}\] and \[\begin{align} \mathbf{U} \tilde{\mathbf{R}} = \mathbf{R} \tilde{\mathbf{R}}^\dagger \tilde{\mathbf{R}} = (\mathbf{R} \mathbf{R}^*)^{-1}(\mathbf{R}\mathbf{R}^*) \mathbf{R} \tilde{\mathbf{R}}^\dagger\tilde{\mathbf{R}} = (\mathbf{R}\mathbf{R}^*)^{-1}\mathbf{R} \tilde{\mathbf{R}}^* (\tilde{\mathbf{R}}\tilde{\mathbf{R}}^\dagger) \tilde{\mathbf{R}} = (\mathbf{R} \mathbf{R}^*)^{-1}(\mathbf{R} \mathbf{R}^*) \mathbf{R} = \mathbf{R}, \end{align}\] where we used the facts that \(\mathbf{R}^*\mathbf{R} = \tilde{\mathbf{R}}^*\tilde{\mathbf{R}}\) and \(\tilde{\mathbf{R}} \tilde{\mathbf{R}}^\dagger = \mathbf{I}_k\). We now claim that \(\mathbf{U}\) is upper triangular. Indeed, if \(\mathbf{U}_{ij}\neq 0\) for \(j<i\), the \(i\)-th row of \(\mathbf{R}\) is given by \(\mathbf{r}_i = \sum_{j=1}^k \mathbf{U}_{ij}\mathbf{\tilde{r}}_j\), and thus has a nonzero entry before its pivot. Since \(\mathbf{U}\) is unitary and all the pivots are positive, we conclude that \(\mathbf{U} = \mathbf{I}_k\) and therefore \(\mathbf{R}= \mathbf{\tilde{R}}\). ◻

Theorem 4. Assume \({\boldsymbol{\mu}}\) is \(n\)-definite and \(\mathrm{rank}~{\boldsymbol{\gamma}}_{n-1} = k-\ell\) with \(0\leq\ell<k\), where \({\boldsymbol{\gamma}}_{n-1}\) is defined in 10 . Let \(\{\mathbf{P}_j(\diamond,{\boldsymbol{\mu}})\}_{j=0}^{n-2}\) denote the orthonormal polynomials associated with \({\boldsymbol{\mu}}\), with recurrence coefficients \(\{\mathbf{A}_j({\boldsymbol{\mu}})\}_{j=0}^{n-3}\) and \(\{\mathbf{B}_j({\boldsymbol{\mu}})\}_{j=0}^{n-3}\). Consider the polynomial \[\begin{align} \label{eq:Pdef} \mathbf{P}(x):=\mathbf{P}(x,{\boldsymbol{\mu}}) = x\mathbf{P}_{n-2}(x,{\boldsymbol{\mu}})-\mathbf{P}_{n-2}(x,{\boldsymbol{\mu}})\mathbf{A}_{n-2}-\mathbf{P}_{n-3}(x,{\boldsymbol{\mu}})\mathbf{B}_{n-3}^*, \end{align}\tag{15}\] where \[\begin{align} \mathbf{A}_{n-2} := \mathbf{A}_{n-2}({\boldsymbol{\mu}}) = \langle \mathbf{P}_{n-2}(\diamond,{\boldsymbol{\mu}}),\diamond,\mathbf{P}_{n-2}(\diamond,{\boldsymbol{\mu}}) \rangle_{\boldsymbol{\mu}}, \end{align}\] and define \(\mathbf{B}_{n-2}:=\mathbf{B}_{n-2}({\boldsymbol{\mu}})\in \mathbb{C}^{(k-\ell)\times k}\) as the unique matrix in row echelon form with positive pivots satisfying \(\langle\mathbf{P},\mathbf{P} \rangle_{\boldsymbol{\mu}}=\mathbf{B}_{n-2}^*\mathbf{B}_{n-2}\). Then, the polynomial \[\begin{align} \label{eq:Pnm1Def} \mathbf{P}_{n-1}(x) := \mathbf{P}_{n-1}(x,{\boldsymbol{\mu}}) = \mathbf{P}(x,{\boldsymbol{\mu}}) \mathbf{B}_{n-2}^\dagger \in \mathbb{C}^{k\times (k-\ell)}, \end{align}\tag{16}\] satisfies \[\begin{align} \label{eq:Pnm1OrthogCond} \langle \mathbf{P}_{n-1},\mathbf{P}_j \rangle_{\boldsymbol{\mu}}= \mathbf{0} \quad \text{for }j<n-1, \quad \langle \mathbf{P}_{n-1},\mathbf{P}_{n-1}\rangle_{\boldsymbol{\mu}}= \mathbf{I}_{k-\ell}, \quad \text{and} \quad \langle \mathbf{P}_{n-1},\diamond \mathbf{P}_{n-2}(\diamond) \rangle_{\boldsymbol{\mu}}= \mathbf{B}_{n-2}. \end{align}\tag{17}\]

Proof. We first verify that \(\mathbf{P}_{n-1}\) is well-defined. From the recurrence of the orthonormal polynomials in 13 , the leading coefficient of \(\mathbf{P}\) is given by \(\mathbf{B}_0^{-1}\dots \mathbf{B}_{n-3}^{-1}\), which implies \[\begin{align} \mathrm{rank} \langle \mathbf{P},\mathbf{P}\rangle_{\boldsymbol{\mu}}= \mathrm{rank}\langle \mathbf{\Pi}_{n-1},\mathbf{\Pi}_{n-1}\rangle_{\boldsymbol{\mu}}= k-\ell. \end{align}\] Lemma 4 then guarantees the existence and uniqueness of \(\mathbf{B}_{n-2}\). Since \(\mathbf{B}_{n-2}\) has full row rank, it admits a right inverse, and therefore \(\mathbf{P}_{n-1}\) in 16 is well-defined. Now, recall that \[\begin{align} \mathbf{B}_{n-3} = \langle \mathbf{P}_{n-2}, \diamond\mathbf{P}_{n-3}(\diamond) \rangle_{\boldsymbol{\mu}}\quad \text{and} \quad \mathbf{A}_{n-2} = \langle \mathbf{P}_{n-2},\diamond\mathbf{P}_{n-2}(\diamond) \rangle_{\boldsymbol{\mu}}, \end{align}\] so \(\langle \mathbf{P},\mathbf{P}_j\rangle_{\boldsymbol{\mu}}= \mathbf{0}\) for \(j<n-1\), and the orthonormality conditions in 17 follow from the definition of \(\mathbf{B}_{n-2}\). Finally, the identity \(\mathbf{B}_{n-2} = \langle \mathbf{P}_{n-1}, \diamond \mathbf{P}_{n-2}(\diamond) \rangle_{\boldsymbol{\mu}}\) follows from \[\begin{align} \langle \mathbf{P},x\mathbf{P}_{n-2} \rangle_{\boldsymbol{\mu}}= \langle \mathbf{P},\mathbf{P} \rangle_{\boldsymbol{\mu}}= \mathbf{B}_{n-2}^*\mathbf{B}_{n-2}, \end{align}\] upon left multiplication by \((\mathbf{B}_{n-2}^\dagger)^*\). ◻

Remark 5. For consistency with 14 , we extend the definition of \(\mathbf{A}_j\) to \(j=n-1\) by setting \[\begin{align} \mathbf{A}_{n-1} := \mathbf{A}_{n-1}({\boldsymbol{\mu}}) = \langle \mathbf{P}_{n-1}(\diamond,{\boldsymbol{\mu}}),\diamond\,\mathbf{P}_{n-1}(\diamond,{\boldsymbol{\mu}})\rangle_{\boldsymbol{\mu}}\in \mathbb{C}^{(k-\ell)\times (k-\ell)} \end{align}\] This coefficient, together with \(\mathbf{A}_{n-2}\) and \(\mathbf{B}_{n-2}\) in Theorem 4, will be important for the inverse spectral analysis carried out in Section 3.4.

3 The spectral map and its inverse for Banded Hermitian matrices↩︎

3.1 The class \(\mathcal{J}_{k,N}\) and its spectral properties↩︎

In this section, we study the class of banded Hermitian matrices introduced in 1 and their connection to matrix-valued measures. We first formalize the set of all such matrices and introduce notation that will be used throughout the paper.

Definition 2. We denote by \(\mathcal{J}_{k,N}\) the set of all \(N\times N\) Hermitian matrices admitting a block tridiagonal representation of the form 1 with \(n=\lceil N/k\rceil\). For each \(\mathbf{J}\in \mathcal{J}_{k,N}\), we denote the diagonal and sub-diagonal blocks of \(\mathbf{J}\) respectively by \[\begin{align} \mathbf{A}_j:=\mathbf{A}_j(\mathbf{J}), \quad j=0,\dots,n-1, \quad \text{and} \quad \mathbf{B}_j = \mathbf{B}_j(\mathbf{J}), \quad j=0,\dotsm,n-2. \end{align}\] These blocks satisfy the following conditions, with \(\ell=nk-N\):

\(\mathbf{A}_0,\dots,\mathbf{A}_{n-2}\in \mathbb{C}^{k\times k}\) and \(\mathbf{A}_{n-1}\in \mathbb{C}^{(k-\ell)\times (k-\ell)}\) are Hermitian,
\(\mathbf{B}_0,\dots,\mathbf{B}_{n-3}\in \mathbb{C}^{k\times k}\) and \(\mathbf{B}_{n-2}\in \mathbb{C}^{(k-\ell)\times k}\) have full rank,
\(\mathbf{B}_0, \dots, \mathbf{B}_{n-3}\) are upper triangular with strictly positive diagonal entries,
\(\mathbf{B}_{n-2}\) is in row echelon form with positive leading entries.

The following result describes eigenvectors corresponding to repeated eigenvalues and the multiplicity of these eigenvalues.

Lemma 5. Let \(\mathbf{J} \in \mathcal{J}_{k,N}\) and suppose \(\mathbf{J} \mathbf{q}_j = \lambda_j \mathbf{q}_j\), \(j = 1,2,\ldots,N\), where \(\{\mathbf{q}_j\}_{j=1}^N\) is an orthonormal basis and \(\lambda_j \leq \lambda_{j+1}\). If \(\lambda_\ell = \lambda_j\) for \(\ell = j,\ldots,j+p\), then⁷ \[\begin{align} (\mathbf{q}_j)_{1:k}, (\mathbf{q}_{j+1})_{1:k}, \ldots, (\mathbf{q}_{j+p})_{1:k}, \end{align}\] are linearly independent. As a result, every eigenvalue of \(\mathbf{J}\) has multiplicity at most \(k\).

Proof. Suppose, for contradiction, that the vectors \((\mathbf{q}_j)_{1:k}, \ldots, (\mathbf{q}_{j+p})_{1:k}\) are linearly dependent, then there exists a linear combination that gives an eigenvector \(\mathbf{v}\) for \(\lambda_j\) that has its first \(k\) entries being all zeros, i.e. \[\begin{align} \mathbf{v} = \begin{bmatrix} \mathbf{0}^T & \mathbf{v}_1^T& \cdots~ & \mathbf{v}_{n-1}^T \end{bmatrix}^T \in \mathbb{C}^N. \end{align}\] Using the recurrence implied by \((\mathbf{J} - \lambda_j) \mathbf{v} = \mathbf{0}\), we find \[\begin{align} \mathbf{A}_0 \mathbf{0} + \mathbf{B}^*_0 \mathbf{v}_1 &= \mathbf{0} \quad \Rightarrow \quad \mathbf{v}_1 = \mathbf{0},\\ &\vdots\\ \mathbf{B}_{j-1} \mathbf{0} + \mathbf{A}_j \mathbf{0} + \mathbf{B}_j^* \mathbf{v}_{j+1} &= \mathbf{0} \quad \Rightarrow \quad \mathbf{v}_{j+1} = \mathbf{0},\\ &\vdots\\ \mathbf{B}_{n-3} \mathbf{0} + \mathbf{A}_{n-2} \mathbf{0} + \mathbf{B}_{n-2}^* \mathbf{v}_{n-1} &= \mathbf{0} \quad \Rightarrow \quad \mathbf{v}_{n-1} = \mathbf{0}. \end{align}\] The last equality follows from the fact that \(\mathbf{B}_{n-2}^*\) has a left inverse. This implies that \[\begin{align} \mathbf{q}_j, \mathbf{q}_{j+1}, \ldots, \mathbf{q}_{j+p}, \end{align}\] are linearly dependent, a contradiction. As a direct consequence, if an eigenvalue \(\lambda_j\) had multiplicity \(p > k\), the corresponding vectors \((\mathbf{q}_j)_{1:k}, \ldots, (\mathbf{q}_{j+p-1})_{1:k}\) would be linearly dependent, which is impossible. Therefore, the multiplicity of any eigenvalue cannot exceed \(k\). ◻

3.2 Matrix polynomial maps↩︎

Let \(\mathbf{P}(x) = \sum_{j=0}^p \Phi_j\,x^j\) be a matrix polynomial with coefficients \(\Phi_j \in \mathbb{C}^{k\times k}\). For \(\mathbf{J}\in \mathbb{C}^{N\times N}\) and \(\mathbf{E}\in\mathbb{C}^{N\times k}\), we define the matrix polynomial map \[\begin{align} \label{eq:PolMapDef} \mathbf{P}(\mathbf{J})\circ\mathbf{E} := \sum_{j=0}^p \mathbf{J}^j\,\mathbf{E}\,\Phi_j \in \mathbb{C}^{N\times k}. \end{align}\tag{18}\] This notation was first introduced in [54] and has been used in [55]–[59] to analyze block Krylov subspaces. An analogous definition applies to vector-valued polynomials, i.e., when \(\Phi_j \in \mathbb{C}^k\). We collect the elementary properties of the polynomial map in the following lemma.

Lemma 6. Let \(\mathbf{P},\mathbf{Q}\), and \(\mathbf{R}\) be matrix polynomials and let \(j\) be a non-negative integer. The polynomial map, defined in 18 , satisfies the following relations:

if \(\mathbf{P}(x) = \mathbf{Q}(x)+\mathbf{R}(x)\), then \(\mathbf{P}(\mathbf{J})\circ \mathbf{E} = \left(\mathbf{Q}(\mathbf{J})\circ \mathbf{E}\right) + \left(\mathbf{R}(\mathbf{J}) \circ \mathbf{E}\right)\),
if \(\mathbf{P}(x) = x^j \,\mathbf{Q}(x)\), then \(\mathbf{P}(\mathbf{J})\circ \mathbf{E} = \mathbf{J}^j\,\left(\mathbf{Q}(\mathbf{J})\circ \mathbf{E}\right)\),
if \(\mathbf{P}(x) = \mathbf{Q}(x)\, \mathbf{C}\) for some \(\mathbf{C}\in \mathbb{C}^{k\times k}\), then \(\mathbf{P}(\mathbf{J})\circ \mathbf{E} = \left(\mathbf{Q}(\mathbf{J}) \circ \mathbf{E}\right)\,\mathbf{C}\).

Throughout most of this document, \(\mathbf{E}\) is chosen as a block selection matrix \(\mathbf{E}_j\) defined by \[\begin{align} \label{eq:SelMat} \mathbf{E}_j := \mathbf{E}_{j}^{N,k} = \begin{bmatrix} \mathbf{0}_{k\times k} &\cdots & \mathbf{0}_{k\times k}& \mathbf{I}_k & \mathbf{0}_{k\times k} & \cdots & \mathbf{0}_{k\times (k-\ell)} \end{bmatrix}^T \in \mathbb{C}^{N\times k} \quad j = 1,\dots,n-1, \end{align}\tag{19}\] where \(\ell = nk-N\) and \(n=\lceil N/k\rceil\). In other words, all blocks of \(\mathbf{E}_j\) are zero except for the \(j\)-th block of size \(k\times k\), which is the identity. For simplicity, the superscripts will be omitted when the dimensions are clear from context.

Lemma 7. Consider \(\mathbf{J} \in \mathcal{J}_{k,N}\) and let \(\{\mathbf{B}_j\}_{j=0}^{n-2}\) denote its subdiagonal blocks. For any integer \(1 \le j \le \lceil N/k\rceil-1\), we have \[\begin{align} \label{eq:PowersJ} \mathbf{E}_{j+1}^* \mathbf{J}^j \mathbf{E}_1 = \mathbf{B}_{j-1} \mathbf{B}_{j-2} \cdots \mathbf{B}_0. \end{align}\tag{20}\]

Proof. For \(j=1\), the statement is trivially true. Now, suppose that 20 holds for \(j-1\), then \[\begin{align} \mathbf{J}^j \mathbf{E}_1 = \mathbf{J} \mathbf{J}^{j-1} \mathbf{E}_1 = \mathbf{J} \begin{bmatrix} \mathbf{X}_1^T & \mathbf{X}_2^T & \dots & \mathbf{X}_{j}^T & \mathbf{0}^T & \cdots & \mathbf{0}^T \end{bmatrix}^T, \quad \text{with} \quad \mathbf{X}_j = \mathbf{B}_{j-2} \cdots \mathbf{B}_0. \end{align}\] This implies that \(\mathbf{E}_{j+1}^* \mathbf{J}^j \mathbf{E}_1 = \mathbf{B}_{j-1} \mathbf{X}_{j}\). ◻

Lemma 8. Consider \(\mathbf{J} \in \mathcal{J}_{k,N}\) with associated blocks \(\{\mathbf{A}_j\}_{j=0}^{n-1}\) and \(\{\mathbf{B}_j\}_{j=0}^{n-2}\) as in Definition 2. Let \(\mathbf{P}_{-1}(x)= \mathbf{0}_k\) and \(\mathbf{P}_0(x)=\mathbf{I}_k\), and define a sequence of matrix polynomials \(\{\mathbf{P}_j\}_{j=1}^{n-1}\) recursively by \[\begin{align} \label{eq:PolRec} \mathbf{P}_{j}(x) = \mathbf{Q}_{j}(x) \mathbf{B}_{j-1}^\dagger, \qquad \mathbf{Q}_{j}(x) = x \mathbf{P}_{j-1}(x) - \mathbf{P}_{j-1}(x)\mathbf{A}_{j-1} - \mathbf{P}_{j-2}(x)\mathbf{B}^*_{j-2}, \qquad j=1,\dots,n-1, \end{align}\tag{21}\] where \(\mathbf{B}_{j-1}^\dagger\) denotes the right inverse of \(\mathbf{B}_{j-1}\) and \(\mathbf{B}_{-1}=\mathbf{I}_k\). Then, for each \(j=0,\dots,n-1\), we have \[\begin{align} \label{eq:OpsFromBandedProp} \mathbf{P}_j(\mathbf{J})\circ \mathbf{E}_1 = \mathbf{E}_{j+1}. \end{align}\tag{22}\]

Proof. For \(j=0\), the claim is trivially satisfied. Next, for \(j=1,\dots,n-1\), using Lemma 6 together with 21 , we find \[\begin{align} \mathbf{P}_j(\mathbf{J})\circ\mathbf{E}_1 = \mathbf{J} \left( \mathbf{P}_{j-1}(\mathbf{J})\circ\mathbf{E}_1 \right)\mathbf{B}_{j-1}^\dagger - \left( \mathbf{P}_{j-1}(\mathbf{J})\circ\mathbf{E}_1 \right)\mathbf{A}_{j-1}\mathbf{B}_{j-1}^\dagger - \left( \mathbf{P}_{j-2}(\mathbf{J})\circ\mathbf{E}_1 \right)\mathbf{B}_{j-2}^* \mathbf{B}_{j-1}^\dagger. \end{align}\] From the definition of \(\mathbf{J}\) in 1 , we have \[\label{eq:ERec} \begin{align} \mathbf{J} \mathbf{E}_{j} = \mathbf{E}_{j-1}\mathbf{B}_{j-2}^*+\mathbf{E}_{j}\mathbf{A}_{j-1} + \mathbf{E}_{j+1}\mathbf{B}_{j-1}, \quad 1\leq j\leq n-1. \end{align}\tag{23}\] Multiplying 23 on the right by \(\mathbf{B}_{j-1}^\dagger\), we see that \(\mathbf{E}_{j+1}\) satisfies the same recurrence as \(\mathbf{P}_j(\mathbf{J})\circ\mathbf{E}_1\), which completes the proof. ◻

3.3 The spectral map↩︎

We define the spectral map \(\varphi: \mathcal{J}_{k,N}\to \mathcal{P}_k(\mathbb{R})\) by \[\begin{align} \label{eq:SpectMapDef1} \varphi(\mathbf{J}) = \sum_{j=1}^N \mathbf{v}_j \mathbf{v}_j^* \delta_{\lambda_j}, \end{align}\tag{24}\] where \(\{\lambda_j\}_{j=1}^N\) are the eigenvalues of \(\mathbf{J}\) and \(\mathbf{v}_j\) is the first \(k\) components of the \(j\)th normalized eigenvector. Equivalently, \(\varphi(\mathbf{J})\) can be written as \[\begin{align} \label{eq:SpectMapDef2} \varphi(\mathbf{J}) = \sum_{j=1}^m \mathbf{V}_j\mathbf{V}_j^* \delta_{x_j}, \end{align}\tag{25}\] where \(\{x_j\}_{j=1}^m\) are the distinct eigenvalues of \(\mathbf{J}\), and \(\mathbf{V}_j\) is the matrix formed by the first \(k\) rows of the eigenvector matrix associated with \(x_j\). This definition does not depend on the choice of eigenbasis. If \(\mathbf{\tilde{V}}_j\) corresponds to a different choice of eigenvectors for \(x_j\), then there exists unitary matrix \(\mathbf{Q}_j\) such that \(\tilde{\mathbf{V}}_j = \mathbf{V}_j \mathbf{Q}_j\), and therefore \(\mathbf{V}_j \mathbf{V}_j^* = \mathbf{\tilde{V}}_j \mathbf{\tilde{V}}_j^*\). It is also worth noting that Lemma 5 implies each \(\mathbf{V}_j \mathbf{V}_j^*\) has rank \(n_j \le k\), and the total rank satisfies \(\sum_{j} n_j = N\) and \(\sum_{j} \mathbf{V}_j \mathbf{V}_j^* = \mathbf{I}_k\).

Remark 6. By Herglotz’s representation theorem[60], the spectral map \(\varphi(\mathbf{J})={\boldsymbol{\mu}}\) admits an equivalent an equivalent characterization as the unique probability measure satisfying \[\begin{align} \mathbf{I}_{N\times k}^*\, (\mathbf{J}-z)^{-1} \,\mathbf{I}_{N\times k} = \int_{\mathbb{R}}\frac{{\boldsymbol{\mu}}(\mathrm{d}x)}{x-z} \quad \text{for }\mathrm{Im}\,z>0. \end{align}\] This formulation is standard in the spectral theory of infinite matrices and operators, and reduces to 24 in the finite case.

Lemma 9. If \(\mathbf{J} \in \mathcal{J}_{k,N}\), then the moments of \({\boldsymbol{\mu}}=\varphi(\mathbf{J})\) satisfy \[\begin{align} \int_{\mathbb{R}} x^i {\boldsymbol{\mu}}(\mathrm d x) = {\mathbf{E}_1}^*\mathbf{J}^i\mathbf{E}_1\quad \text{for }i\geq 0. \end{align}\] Moreover, for polynomials \(\mathbf{P},\mathbf{Q}\), we have \[\begin{align} \langle \mathbf{P}, \mathbf{Q} \rangle_{\boldsymbol{\mu}}= (\mathbf{P}(\mathbf{J})\circ\mathbf{E}_1)^*(\mathbf{Q}(\mathbf{J})\circ\mathbf{E}_1). \end{align}\]

Proof. The first identity follows directly from the definition of the spectral measure. Specifically, let \(\mathbf{J} = \mathbf{U} \Lambda \mathbf{U}^*\) denote the eigendecomposition of the banded matrix, then \[\begin{align} \int_{\mathbb{R}} x^i {\boldsymbol{\mu}}(\mathrm d x) = \sum_{j=1}^N \lambda_j^i \mathbf{v}_j \mathbf{v}_j^* = \mathbf{E}_1^* \mathbf{U} \Lambda^i\mathbf{U}^* \mathbf{E}_1 = \mathbf{E}_1^*\mathbf{J}^i\mathbf{E}_1. \end{align}\] Now, suppose that \(\mathbf{P}(x) = \sum_{i=0}^p \mathbf{C}_ix^i\) and \(\mathbf{Q}(x) = \sum_{i=0}^q \mathbf{D}_i x^i\), then \[\begin{align} (\mathbf{P}(\mathbf{J})\circ\mathbf{E}_1)^* (\mathbf{Q}(\mathbf{J})\circ\mathbf{E}_1) = \sum_{i=0}^p\sum_{j=0}^q \mathbf{C}_i^* \left( {\mathbf{E}_1}^*\mathbf{J}^{i+j}\mathbf{E}_1 \right)\mathbf{D}_j = \sum_{i=0}^p\sum_{j=0}^q \mathbf{C}_i^* \left( \int_{\mathbb{R}}x^{i+j}{\boldsymbol{\mu}}(\mathrm dx) \right)\mathbf{D}_j = \langle \mathbf{P}, \mathbf{Q} \rangle_{\boldsymbol{\mu}}. \end{align}\] ◻

Theorem 7 shows that the quasi-inner product induced by the measure \({\boldsymbol{\mu}}= \varphi(\mathbf{J})\), where \(\mathbf{J} \in \mathcal{J}_{k,N}\) with \(N = kn - \ell\), defines a norm on the space of polynomials of degree at most \(n-2\), implying that \({\boldsymbol{\mu}}\) is \(n\)-definite. Moreover, there exist exactly \(\ell\) independent vector polynomials of degree \(n-1\) with zero norm. These results are fundamental for the identifying the range of the spectral map in Section 3.4.

Theorem 7. Let \({\boldsymbol{\mu}}= \varphi(\mathbf{J})\), where \(\mathbf{J} \in \mathcal{J}_{k,N}\) with \(N = kn - \ell\), \(0 \leq \ell < k\), and \[\begin{align} \label{eq:NdDef} \mathcal{N}_d = \left\{ \mathbf{p}(x) = \sum_{j=0}^d \mathbf{c}_j x^j,~ \mathbf{c}_j \in \mathbb{C}^k: \langle \mathbf{p}, \mathbf{p} \rangle_{\boldsymbol{\mu}}= 0 \right\}. \end{align}\tag{26}\] Then the dimension of \(\mathcal{N}_d\) satisfies \[\begin{align} \label{eq:NdDimEq} \dim \mathcal{N}_d = \begin{cases} 0, & \text{if } d < n - 1, \\ \ell, & \text{if } d = n-1. \end{cases} \end{align}\tag{27}\]

Proof. Suppose first that \(d < n - 1\), and let \(\mathbf{p}(x) = \sum_{j=0}^d \mathbf{c}_j x^j\) with \(\mathbf{c}_d \neq \mathbf{0}\). Assume \(\langle \mathbf{p}, \mathbf{p} \rangle_{\boldsymbol{\mu}}= 0\), then by Lemma 9 we have \((\mathbf{p}(\mathbf{J})\circ\mathbf{E}_1)^*(\mathbf{p}(\mathbf{J})\circ\mathbf{E}_1) = 0\), which implies \[\begin{align} \label{eq:ContradicEqForDef} \mathbf{J}^d \mathbf{E}_1 \mathbf{c}_d + \sum_{j=0}^{d-1} \mathbf{J}^j \mathbf{E}_1 \mathbf{c}_j = \mathbf{0}. \end{align}\tag{28}\] Multiplying 28 on the left by \(\mathbf{E}_{d+1}^*\) and using Lemma 7, we get \(\mathbf{B}_{p-1}\dots\mathbf{B}_0 \mathbf{c}_d = \mathbf{0}\), which contradicts the fact that \(\mathbf{B}_{p-1}\dots\mathbf{B}_0\) is nonsingular.

Next, assume \(\dim \mathcal{N}_{n-1} = r > \ell\). Then there exist linearly independent polynomials \(\{ \mathbf{p}_i(x) \}_{i=1}^r\) such that \[\begin{align} \mathbf{p}_i(x) = \sum_{j=1}^{n-1} \mathbf{c}_j^{(i)} x^j \in \mathbb{C}^k,\quad \mathbf{c}_{n-1}^{(i)} \neq \mathbf{0}, \quad \langle \mathbf{p}_i, \mathbf{p}_i \rangle_{\boldsymbol{\mu}}= 0. \end{align}\] We note that \(\langle\mathbf{p}_i,\mathbf{p}_i \rangle_{\boldsymbol{\mu}}= 0\) is equivalent to \[\begin{align} \label{eq:ZeroNormEquiv} \mathbf{V}_j^* \mathbf{p}_i(x_j) = \mathbf{0} \qquad \text{for } j=1,\dots,m, \end{align}\tag{29}\] where \({\boldsymbol{\mu}}= \sum_{j=1}^m \mathbf{V}_j \mathbf{V}_j^* \delta_{x_j}\). We also claim that the leading coefficients \(\{\mathbf{c}_{n-1}^{(i)}\}_{i=1}^r\) are linearly independent. Indeed, if there exists nonzero constants \(\alpha_1,\dots,\alpha_r\) such that \(\sum_{i=1}^r \alpha_i \mathbf{c}_{n-1}^{(i)} = \mathbf{0}\), then the polynomial \(\mathbf{p}(x) = \sum_i\alpha_i \mathbf{p}_i(x)\) is of degree \(n-2\) and satisfies \(\langle \mathbf{p},\mathbf{p} \rangle_{\boldsymbol{\mu}}= 0\) by 29 , contradicting \(\mathrm{dim}~ \mathcal{N}_{n-2}=0\) because \(\mathbf{p} \neq \mathbf{0}\) by linear independence. Now, using Lemmas 7 and 9 again, we have \(\mathbf{B}_{n-2}\dots\mathbf{B}_0 \mathbf{c}_{n-1}^{(i)} = \mathbf{0}\) for \(i = 1,\dots,r\). Since \(\mathbf{B}_{n-1}\dots\mathbf{B}_0\) is nonsingular, we find \(\dim \ker(\mathbf{B}_{n-2}) = r > \ell\), which contradicts the assumption that \(\mathbf{B}_{n-2} \in \mathbb{C}^{(k-\ell)\times k}\) has full rank, and therefore \(r \leq \ell\).

To establish \(r=\ell\), consider the polynomial \(\mathbf{Q}_{n-1}\) defined in Lemma 8 and note that \[\begin{align} \label{eq:RankPTilde} \langle \mathbf{Q}_{n-1},\mathbf{Q}_{n-1} \rangle_{\boldsymbol{\mu}}&= \langle \diamond\mathbf{P}_{n-2}(\diamond)-\mathbf{P}_{n-2}\mathbf{A}_{n-2}-\mathbf{P}_{n-3}\mathbf{B}_{n-3}^*, \diamond\mathbf{P}_{n-2}(\diamond)-\mathbf{P}_{n-2}\mathbf{A}_{n-2}-\mathbf{P}_{n-3}\mathbf{B}_{n-3}^* \rangle_{\boldsymbol{\mu}}= \mathbf{B}_{n-2}^*\mathbf{B}_{n-2}, \end{align}\tag{30}\] where the second equality follows from Lemma 8 and Lemma 9. Hence \(\mathrm{rank}~\langle \mathbf{Q}_{n-1}, \mathbf{Q}_{n-1} \rangle_{\boldsymbol{\mu}}= k-\ell\), and there exists linearly independent vectors \(\{\mathbf{v}_i\}_{i=1}^{\ell}\) such that \[\begin{align} \langle \mathbf{p}_i, \mathbf{p}_i \rangle = 0, \qquad \mathbf{p}_i(x) := \mathbf{Q}_{n-1}(x)\mathbf{v}_i, \quad \text{for } i=1,\dots \ell. \end{align}\] Since the leading coefficient of \(\mathbf{Q}_{n-1}\), given by \(\mathbf{B}_0^{-1}\cdots \mathbf{B}_{n-3}^{-1}\), is nonsingular, the vector polynomials \(\{\mathbf{p}_i\}_{i=1}^\ell\) are linearly independent and have nonzero leading coefficients. We conclude that \(\{\mathbf{p}_i\}_{i=1}^\ell \subset \mathcal{N}_{n-1}\) which implies \(\mathrm{dim} ~ \mathcal{N}_{n-1} \geq \ell\), and therefore \(\mathrm{dim} ~ \mathcal{N}_{n-1} = \ell\). ◻

Theorem 7 admits an equivalent formulation in terms of matrix polynomials and the structure of the spectral measure, as presented in Theorem 8. Statement [item:Nd-rank] in Theorem 8 shows that when \(N = kn - \ell\), the norms of the matrix orthogonal polynomials are of full rank for all polynomials of degree up to \(n-2\), and of rank \(k - \ell\) for those of degree \(n-1\). Statement [item:Nd-meas] describes the conditions on the support and weights of a measure required for the existence of such matrix orthogonal polynomials.

Theorem 8. Let \(\mathbf{J} \in \mathcal{J}_{k,N}\) with \(N = kn - \ell\) and \(0 \leq \ell < k\), and let \({\boldsymbol{\mu}}= \varphi(\mathbf{J})=\sum_{j=1}^m \mathbf{V}_j\mathbf{V}_j^* \delta_{x_j}\) be the associated spectral measure, as defined in 25 . The following statements are equivalent:

\(\mathcal{N}_d\) defined in 26 satisfies 27 .
For every monic matrix polynomial \(\mathbf{\Pi}\) of degree \(d\), we have \[\begin{align} \mathrm{rank}~\langle \mathbf{\Pi},\mathbf{\Pi}\rangle_{{\boldsymbol{\mu}}} \geq \begin{cases} k, & d < n-1,\\ k-\ell, & d = n-1, \end{cases} \end{align}\] and, in particular, \[\begin{align} \mathrm{rank}~\langle \mathbf{\Pi}_{n-1}, \mathbf{\Pi}_{n-1}\rangle_{\boldsymbol{\mu}}= k-\ell \end{align}\] where \(\mathbf{\Pi}_{n-1}=\mathbf{\Pi}_{n-1}(\diamond,{\boldsymbol{\mu}})\) denotes the monic orthogonal polynomial of degree \(n-1\).
Let \[\begin{align} \mathbf{X} = \mathrm{diag}(\underbrace{x_1,\dots,x_1}_{n_1},\underbrace{x_2,\dots,x_2}_{n_2},\ldots,\underbrace{x_m,\dots,x_m}_{n_m}), \quad \mathbf{V}^* = \begin{bmatrix} \mathbf{V}_1 & \cdots & \mathbf{V}_m \end{bmatrix}, \quad \mathbf{V} \in \mathbb{C}^{N \times k}, \end{align}\] where \(n_j = \mathrm{rank}~\mathbf{V}_j\), then the matrix \[\begin{align} \mathbf{M}_d(\mathbf{V}, \mathbf{X}) := \begin{bmatrix} \mathbf{V} & \mathbf{X} \mathbf{V} & \cdots & \mathbf{X}^d \mathbf{V} \end{bmatrix} \in \mathbb{C}^{N\times (d+1)k} \end{align}\] has full rank, that is, \[\begin{align} \mathrm{rank}~\mathbf{M}_d(\mathbf{V},\mathbf{X}) = \min\{(d+1) k, N\}. \end{align}\]

Proof. The proof proceeds in two steps. We first show that [item:Nd-dim] and [item:Nd-rank] are equivalent, and then establish the equivalence between [item:Nd-dim] and [item:Nd-meas].

Step 1: Equivalence of [item:Nd-dim] and [item:Nd-rank]. First, suppose that [item:Nd-dim] holds. If \(\mathbf{\Pi}\) is a monic polynomial of degree \(d\), then \(\mathrm{rank}\,\langle \mathbf{\Pi}, \mathbf{\Pi} \rangle_{\boldsymbol{\mu}}= k\) for \(d < n-1\); otherwise, there would exist a nonzero vector \(\mathbf{v}\) such that \(\mathbf{p}(x) = \mathbf{\Pi}(x)\mathbf{v}\) satisfies \(\langle \mathbf{p}, \mathbf{p} \rangle_{\boldsymbol{\mu}}= 0\), contradicting [item:Nd-dim]. Next, assume that \(\mathrm{rank}\,\langle \mathbf{\Pi}, \mathbf{\Pi} \rangle_{\boldsymbol{\mu}}= k - r\) for \(d = n-1\) with \(r > \ell\), then there exists linearly independent vectors \(\{\mathbf{v}_i\}_{i=1}^r\) such that \[\begin{align} \langle \mathbf{p}_i, \mathbf{p}_i \rangle_{\boldsymbol{\mu}}= 0, \quad \mathbf{p}_i(x) := \mathbf{\Pi}(x)\mathbf{v}_i, \quad i = 1,\dots,r. \end{align}\] Since \(\mathbf{\Pi}\) is monic, each \(\mathbf{p}_i\) is a vector polynomial of degree \(n-1\) and the set \(\{\mathbf{p}_i\}_{i=1}^r\) is linearly independent, implying \(\dim \mathcal{N}_{n-1} = r \geq \ell\), again contradicting [item:Nd-dim]. To show that the rank condition holds for \(\mathbf{\Pi}_{n-1}\), consider the polynomial \(\mathbf{Q}_{n-1}\) defined in Lemma 8. Since \(\mathbf{\Pi}_{n-1}(x) = \mathbf{Q}_{n-1}(x)\mathbf{B}_{n-3}\cdots\mathbf{B}_{0}\), we have \[\begin{align} \mathrm{rank}~\langle \mathbf{\Pi}_{n-1}, \mathbf{\Pi}_{n-1} \rangle_{\boldsymbol{\mu}} = \mathrm{rank}~\langle \mathbf{Q}_{n-1}, \mathbf{Q}_{n-1} \rangle_{\boldsymbol{\mu}} = k - \ell, \end{align}\] where the second equality follows from 30 .

Now suppose that [item:Nd-rank] holds, and assume that there exists a non-trivial polynomial \(\mathbf{p}\in \mathbb{C}^k\) of degree \(d\) with \(d<n-1\) such that \(\langle \mathbf{p}, \mathbf{p} \rangle_{\boldsymbol{\mu}}= 0\). Let \(\mathbf{P} \in\mathbb{C}^{k\times k}\) be a matrix polynomial of degree \(d\) whose leading coefficient \(\mathbf{C}_d\) is invertible and whose first column is \(\mathbf{p}\). Define \(\mathbf{\Pi} = \mathbf{P} \mathbf{C}_d^{-1}\) and \(\mathbf{v} = \mathbf{C}_d\mathbf{e}_1\), then \[\begin{align} \mathbf{v}^* \langle \mathbf{\Pi}, \mathbf{\Pi} \rangle_{\boldsymbol{\mu}}\mathbf{v} = \langle \mathbf{p}, \mathbf{p} \rangle_{\boldsymbol{\mu}}= 0. \end{align}\] Hence \(\langle \mathbf{\Pi}, \mathbf{\Pi} \rangle_{\boldsymbol{\mu}}\) is rank deficient, which leads to a contradiction. Since \(\mathrm{rank}\,\langle \mathbf{\Pi}_{n-1}, \mathbf{\Pi}_{n-1} \rangle_{\boldsymbol{\mu}}= k - \ell\), we have \(\dim \mathcal{N}_{n-1} \ge \ell\), and it remains to prove that \(\dim \mathcal{N}_{n-1} \le \ell\). If this were false, there would exist \(r > \ell\) linearly independent vector polynomials \(\{\mathbf{p}_i\}_{i=1}^r\) of degree \(n-1\) such that \(\langle \mathbf{p}_i, \mathbf{p}_i \rangle_{\boldsymbol{\mu}}= 0\) for \(i = 1, \dots, r\). Since \(\mathrm{dim}~\mathcal{N}_{n-2}=0\), the leading coefficients of \(\{\mathbf{p}_i\}_{i=1}^r\) are linearly independent. As a result, there exists a matrix polynomial \(\mathbf{P}\) of degree \(d\) whose leading coefficient \(\mathbf{C}_{n-1}\) is invertible and whose first \(r\) columns are given by \(\mathbf{p}_1,\dots, \mathbf{p}_r\). Setting \(\mathbf{\Pi} = \mathbf{P} \mathbf{C}_d^{-1}\) and \(\mathbf{v}_i = \mathbf{C}_{n-1} \mathbf{e}_i\) for \(i=1,\dots,r\), we find \[\begin{align} \mathbf{v}_i^* \langle \mathbf{\Pi}, \mathbf{\Pi} \rangle_{\boldsymbol{\mu}}\mathbf{v}_i = \langle \mathbf{p}_i, \mathbf{p}_i \rangle_{\boldsymbol{\mu}} = 0, \qquad i = 1, \dots, r, \end{align}\] so \(\mathrm{rank}~\langle \mathbf{\Pi}, \mathbf{\Pi} \rangle_{\boldsymbol{\mu}}<k-\ell\), contradicting [item:Nd-rank].

Step 2: Equivalence of [item:Nd-dim] and [item:Nd-meas]. For a polynomial \(\mathbf{p}(x) = \sum_{j=0}^{d} \mathbf{c}_j x^j\), we have \[\begin{align} \begin{bmatrix} \mathbf{V}_1^* \mathbf{p}(x_1) \\ \mathbf{V}_2^* \mathbf{p}(x_2) \\ \vdots \\ \mathbf{V}_m^* \mathbf{p}(x_m) \end{bmatrix} = \begin{bmatrix} \mathbf{V}_1^* & x_1 \mathbf{V}_1^* & \cdots & x_1^{d} \mathbf{V}_1^* \\ \mathbf{V}_2^* & x_2 \mathbf{V}_2^* & \cdots & x_2^{d} \mathbf{V}_2^* \\ \vdots & \vdots && \vdots \\ \mathbf{V}_m^* & x_m \mathbf{V}_m^* & \cdots & x_m^{d} \mathbf{V}_m^* \end{bmatrix}\begin{bmatrix} \mathbf{c}_0 \\ \mathbf{c}_1 \\ \vdots \\ \mathbf{c}_d \end{bmatrix} = \mathbf{M}_d(\mathbf{V}, \mathbf{X}) \begin{bmatrix} \mathbf{c}_0 \\ \mathbf{c}_1 \\ \vdots \\ \mathbf{c}_d \end{bmatrix}, \end{align}\] which implies \(\mathrm{dim}~\mathcal{N}_d = \mathrm{dim}~\mathrm{ker} ~\mathbf{M}_d(\mathbf{V},\mathbf{X})\). Thus, the matrix \(\mathbf{M}_d(\mathbf{V},\mathbf{X})\) has full rank if and only if [item:Nd-dim] is satisfied, completing our proof. ◻

Remark 9. For the spectral measure \({\boldsymbol{\mu}}\) defined in 25 , the matrices \(\mathbf{V}_j\) are defined only up to right multiplication by a unitary matrix. However, for any unitary matrices \(\mathbf{Q}_j\), setting \(\mathbf{\tilde{V}}_j = \mathbf{V}_j \mathbf{Q}_j\) gives \[\begin{align} \mathbf{M}_d(\mathbf{\tilde{V}}, \mathbf{X}) = \mathrm{diag}(\mathbf{Q}_1^*,\dots,\mathbf{Q}_m^*)\, \mathbf{M}_d(\mathbf{V}, \mathbf{X}). \end{align}\] Therefore, if the condition in Theorem 8 [item:Nd-meas] is satisfied for one choice \(\{\mathbf{V}_j\}_{j=1}^m\), it is satisfied for any other choice \(\{\mathbf{\tilde{V}}_j\}_{j=1}^m\).

The results concerning the null space \(\mathcal{N}_d\) and the rank of matrix polynomials for higher degrees can be easily obtained once the conditions in Theorem 8 are satisfied, as stated in the next corollary.

Corollary 1. Let \({\boldsymbol{\mu}}\) be a spectral measure associated with a matrix in \(\mathcal{J}_{k,N}\) with \(N = kn-\ell\). For every \(d \ge n\), there exists a monic matrix polynomial \(\mathbf{\Pi}\) of degree \(d\) such that \(\langle \mathbf{\Pi}, \mathbf{\Pi} \rangle_{{\boldsymbol{\mu}}} = \mathbf{0}\). As a result, we have \[\begin{align} \mathrm{dim}~\mathcal{N}_d = k(d+1)-N, \qquad \text{and} \qquad \langle \mathbf{\Pi}_d,\mathbf{\Pi}_d \rangle_{\boldsymbol{\mu}}= \mathbf{0}, \end{align}\] where \(\mathbf{\Pi}_d=\mathbf{\Pi}_d(\diamond,{\boldsymbol{\mu}})\) is any monic polynomial of degree \(d\) satisfying \(\langle \mathbf{\Pi}_d,\mathbf{P}\rangle_{\boldsymbol{\mu}}=\mathbf{0}\) for all matrix polynomials \(\mathbf{P}\) of degree \(p<d\).

Proof. Theorems 7 and Theorem 8 [item:Nd-meas] guarantee that, for \(d \geq n\), the matrix \(\mathbf{M}_{d-1}(\mathbf{V}, \mathbf{X})\) has linearly independent rows and therefore admits a right inverse \(\mathbf{M}_{d-1}^{\dagger}(\mathbf{V}, \mathbf{X})\). Define \(\{\mathbf{C}_j\}_{j=0}^{d-1}\) such that \(\mathbf{C}_j \in \mathbb{C}^{k\times k}\) and \[\begin{align} \begin{bmatrix} \mathbf{C}_0 \\ \mathbf{C}_1\\ \vdots \\ \mathbf{C}_{d-1} \end{bmatrix} = \mathbf{M}_{d-1}^\dagger(\mathbf{V},\mathbf{X}) \begin{bmatrix} -\mathbf{V}_1^* x_1^d \\ -\mathbf{V}_2^* x_2^d \\ \vdots \\ -\mathbf{V}_m^* x_m^d \end{bmatrix}, \end{align}\] and consider \(\mathbf{\Pi}(x) = \mathbf{I}_k x^d+\sum_{j=0}^{d-1}\mathbf{C}_j x^j\). By construction, we have \[\begin{align} \begin{bmatrix} \mathbf{V}_1^* \mathbf{\Pi}(x_1)\\ \mathbf{V}_2^* \mathbf{\Pi}(x_2)\\ \vdots\\ \mathbf{V}_m^* \mathbf{\Pi}(x_m) \end{bmatrix} = \mathbf{M}_d(\mathbf{V},\mathbf{X}) \begin{bmatrix} \mathbf{C}_0\\ \vdots \\ \mathbf{C}_{d-1}\\ \mathbf{I}_k \end{bmatrix} = \mathbf{M}_{d-1}(\mathbf{V},\mathbf{X}) \begin{bmatrix} \mathbf{C}_0 \\ \mathbf{C}_1 \\ \vdots \\ \mathbf{C}_{d-1} \end{bmatrix} + \begin{bmatrix} \mathbf{V}_1^* x_1^d \\ \mathbf{V}_2^* x_2^d \\ \vdots \\ \mathbf{V}_m^* x_m^d \end{bmatrix}, \end{align}\] and therefore \[\begin{align} \label{eq:ZeroCond} \mathbf{V}_j^* \mathbf{\Pi}(x_j) = \mathbf{0} \quad \text{for}~ j=1,\dots,m, \end{align}\tag{31}\] which implies \(\langle \mathbf{\Pi}, \mathbf{\Pi} \rangle_{\boldsymbol{\mu}}= \mathbf{0}\).

Since \(\dim \mathcal{N}_{n-1} = \ell\) and each degree increment gives at most \(k\) dimensions, we have \[\begin{align} \mathrm{dim}~\mathcal{N}_d \leq \ell+k(d-n+1) = k(d+1)-N. \end{align}\] The columns of \(\mathbf{\Pi}\) are linearly independent elements of \(\mathcal{N}_d\) for all \(d \ge n\), so the above inequality is in fact an equality. Finally, since \(\mathbf{\Pi}_d\) is a monic orthogonal polynomial of degree \(d\), the difference \(\mathbf{\Pi}_d - \mathbf{\Pi}\) has degree at most \(d-1\). Using Lemma 1, together with the orthogonality of \(\mathbf{\Pi}_d\) and 31 , we have \[\begin{align} \langle \mathbf{\Pi}_n-\mathbf{ \Pi},\mathbf{\Pi}_n-\mathbf{ \Pi} \rangle_{\boldsymbol{\mu}}= \langle \mathbf{\Pi}_n ,\mathbf{\Pi}_n-\mathbf{ \Pi} \rangle_{\boldsymbol{\mu}}- \langle \mathbf{\Pi},\mathbf{\Pi}_n-\mathbf{ \Pi}\rangle_{\boldsymbol{\mu}}= \mathbf{0}. \end{align}\] This implies that \(\mathbf{V}_j^* \left(\mathbf{\Pi}_n(x_j)-\mathbf{ \Pi}(x_j)\right) = \mathbf{0}\) for \(j=1,\dots,m\), and hence \(\langle \mathbf{\Pi}_n, \mathbf{\Pi}_n \rangle_{\boldsymbol{\mu}}= \mathbf{0}\). ◻

In the remainder of this section, we show that each spectral measure corresponds to a unique matrix \(\mathbf{J}\) within the class \(\mathcal{J}_{k,N}\). Hence, the spectral map is invertible and the explicit construction of its inverse will be addressed in the next section.

Theorem 10. The spectral map \(\varphi: \mathcal{J}_{k,N}\to \mathcal{P}_k(\mathbb{R})\), defined in 24 , is injective.

Proof. Let \(\mathbf{J} \in \mathcal{J}_{k,N}\) and suppose that \({\boldsymbol{\mu}}= \varphi(\mathbf{J})\). Consider the sequences of polynomials \(\{\mathbf{Q}_j\}_{j=0}^{n-1}\) and \(\{\mathbf{P}_j\}_{j=0}^{n-1}\) associated with \(\mathbf{J}\) as defined in Lemma 8, then \[\begin{align} \langle \mathbf{P}_i, \mathbf{P}_j \rangle_{\boldsymbol{\mu}} = (\mathbf{P}_i(\mathbf{J})\circ \mathbf{E}_1)^* \, (\mathbf{P}_j(\mathbf{J}) \circ\mathbf{E}_1) = \mathbf{E}_{i+1}^* \mathbf{E}_{j+1} = \mathbf{I} \delta_{ij}, \end{align}\] where the first equality follows from Lemma 9 and the second from Lemma 8. Thus, \(\{\mathbf{P}_j\}_{j=0}^{n-1}\) are orthonormal polynomials with respect to \({\boldsymbol{\mu}}\). Now let \(\mathbf{\hat{J}} \in \mathcal{J}_{k,N}\) be another banded matrix such that \(\varphi(\mathbf{\hat{J}}) = {\boldsymbol{\mu}}\), and let \(\{\mathbf{Q}_j\}_{j=0}^{n-1}\) and \(\{\mathbf{\hat{P}}_j\}_{j=0}^{n-1}\) be the associated polynomials from Lemma 8. By the same argument, \(\{\mathbf{\hat{P}}_j\}_{j=0}^{n-1}\) also form a set of orthonormal polynomials with respect to \({\boldsymbol{\mu}}\).

We now prove that the two families of polynomials coincide. Let \(\{\mathbf{A}_j\}_{j=0}^{n-1}\), \(\{\mathbf{B}_j\}_{j=0}^{n-2}\) and \(\{\hat{\mathbf{A}}_j\}_{j=0}^{n-1}\), \(\{\hat{\mathbf{B}}_j\}_{j=0}^{n-2}\) denote the diagonal and off-diagonal blocks of \(\mathbf{J}\) and \(\hat{\mathbf{J}}\), respectively, and define \[\begin{align} \mathbf{C}_j = \mathbf{B}_{j-1}\dots \mathbf{B}_0, \quad \hat{\mathbf{C}}_j = \hat{\mathbf{B}}_{j-1}\dots \hat{\mathbf{B}}_0, \quad \text{for }j=0,\dots,n-1. \end{align}\] The recurrence relation in Lemma 8 implies that \(\mathbf{\Pi}_j := \mathbf{Q}_j \mathbf{C}_{j-1}\) and \(\hat{\mathbf{\Pi}}_j := \hat{\mathbf{Q}}_j \hat{\mathbf{C}}_{j-1}\) are monic orthogonal polynomials of degree \(j\). By Theorem 2 and Theorem 7, monic orthogonal polynomials of degree at most \(n-1\) are unique, so \(\mathbf{\Pi}_j = \hat{\mathbf{\Pi}}_j\) and \[\begin{align} \label{eq:EqMonicNorms} \langle \mathbf{\Pi}_j,\mathbf{\Pi}_j \rangle_{\boldsymbol{\mu}}= \langle \hat{\mathbf{\Pi}}_j, \hat{\mathbf{\Pi}}_j \rangle_{\boldsymbol{\mu}}, \quad \text{for }j=0,\dots,n-1. \end{align}\tag{32}\] Using Lemma 8 and Lemma 9, we have \[\begin{align} \label{eq:QjNorm} \langle \mathbf{Q}_{j}, \mathbf{Q}_{j} \rangle_{\boldsymbol{\mu}}= \langle \diamond\mathbf{P}_{j-1}(\diamond)-\mathbf{P}_{j-1}\mathbf{A}_{j-1}-\mathbf{P}_{j-2}\mathbf{B}_{j-2}^*, \diamond\mathbf{P}_{j-1}(\diamond)-\mathbf{P}_{j-1}\mathbf{A}_{j-1}-\mathbf{P}_{j-2}\mathbf{B}_{j-2}^* \rangle_{\boldsymbol{\mu}}= \mathbf{B}_{j-1}^*\mathbf{B}_{j-1}, \end{align}\tag{33}\] and similarly \(\langle \hat{\mathbf{Q}}_j,\hat{\mathbf{Q}}_j \rangle_{{\boldsymbol{\mu}}} = \hat{\mathbf{B}}_{j-1}^* \hat{\mathbf{B}}_{j-1}\), for \(j=1,\dots, n-1\). Combining 32 and 33 gives \[\begin{align} \mathbf{C}_j^*\mathbf{C}_j = \hat{\mathbf{C}}_j^*\hat{\mathbf{C}}_j, \quad \text{for } j=0,\dots,n-1. \end{align}\] Lemma 4 shows that \(\mathbf{C}_j = \hat{\mathbf{C}}_j\) and thus \(\mathbf{Q}_j = \hat{\mathbf{Q}}_j\) which implies \(\mathbf{P}_j = \hat{\mathbf{P}}_j\). Finally, by another application of Lemma 8 and Lemma 9, we have \[\begin{align} \mathbf{A}_j = (\mathbf{P}_j(\mathbf{J})\circ\mathbf{E}_1)^* \mathbf{J} (\mathbf{P}_j(\mathbf{J}) \circ\mathbf{E}_1) = \langle \mathbf{P}_j,\diamond\mathbf{P}_j(\diamond) \rangle_{\boldsymbol{\mu}}= \langle \mathbf{\hat{P}}_j, \diamond\mathbf{\hat{P}}_j(\diamond)\rangle_{\boldsymbol{\mu}}= (\mathbf{\hat{P}}_j(\mathbf{\hat{J}})\circ\mathbf{E}_1)^* \mathbf{\hat{J}}(\mathbf{\hat{P}}_j(\mathbf{\hat{J}})\circ\mathbf{E}_1)= \mathbf{\hat{A}}_j, \end{align}\] and \[\begin{align} \mathbf{B}_j = (\mathbf{P}_{j+1}(\mathbf{J})\circ\mathbf{E}_1)^* \mathbf{J} (\mathbf{P}_j(\mathbf{J})\circ\mathbf{E}_1) = \langle \mathbf{P}_{j+1},\diamond\mathbf{P}_j(\diamond) \rangle_{\boldsymbol{\mu}}= \langle \mathbf{\hat{P}}_{j+1}, \diamond\mathbf{\hat{P}}_j(\diamond)\rangle_{\boldsymbol{\mu}}= (\mathbf{\hat{P}}_{j+1}(\mathbf{\hat{J}})\circ\mathbf{E}_1)^* \mathbf{\hat{J}} (\mathbf{\hat{P}}_j(\mathbf{\hat{J}})\circ\mathbf{E}_1)= \mathbf{\hat{B}}_j, \end{align}\] and therefore \(\mathbf{J} = \mathbf{\hat{J}}\). ◻

3.4 Inverse spectral map↩︎

Given that the spectral map \(\varphi\) is injective, it is natural to determine its range and describe the procedure for constructing its inverse. Motivated by Theorems 7 and 8, we expect the range of \(\varphi\) to consist of the measures belonging to the set \[\begin{align} \label{eq:Mkn} \mathcal{M}_{k,N} := \left\{ \sum_{j=1}^m \mathbf{W}_j\delta_{x_j} : \sum_{j=1}^m\mathbf{W}_j=\mathbf{I}_k,\; \mathbf{W}_j\!\ge\!0,\; \textstyle\sum_{j=1}^m\mathrm{rank}\,\mathbf{W}_j=N,\; \text{\eqref{eq:NdDimEq} holds for } N=kn-\ell,\;0\le\ell<k \right\}. \end{align}\tag{34}\] In this section, we prove that the range of \(\varphi\) coincides exactly with \(\mathcal{M}_{N,k}\) and define the inverse spectral map, which is closely related to the construction of orthonormal matrix polynomials.

Let \({\boldsymbol{\mu}}\in \mathcal{M}_{k,N}\), then by Lemma 8 the sequences of monic orthogonal polynomials \(\{\mathbf{\Pi}_j(\diamond,{\boldsymbol{\mu}})\}_{j=0}^{n-1}\) and orthonormal matrix polynomials \(\{\mathbf{P}_j(\diamond,{\boldsymbol{\mu}})\}_{j=0}^{n-2}\) are well defined and are given by the recurrence relations in 11 and 13 . Moreover, the \((n-1)\)-th orthonormal polynomial \(\mathbf{P}_{n-1}\) is uniquely defined, as established in Theorem 4. The recurrence coefficients \(\{\mathbf{A}_j\}_{j=0}^{n-1}\) and \(\{\mathbf{B}_j\}_{j=0}^{n-2}\) associated with the orthonormal polynomials \(\{\mathbf{P}_{j}\}_{j=0}^{n-1}\) define a Hermitian, block tridiagonal matrix \[\begin{align} \label{eq:BlockJac} \mathbf{J} = \mathbf{J}({\boldsymbol{\mu}}) = \begin{bmatrix} \mathbf{A}_0 & \mathbf{B}_0^* & & & & \\ \mathbf{B}_0 & \mathbf{A}_1 & \mathbf{B}_1^* & & & \\ & \mathbf{B}_1 & \mathbf{A}_2 & \mathbf{B}_2^* & & \\ & & \mathbf{B}_2 & \mathbf{A}_3 & \ddots & \\ & & & \ddots & \ddots & \mathbf{B}_{n-2}^* \\ & & & & \mathbf{B}_{n-2} & \mathbf{A}_{n-1} \end{bmatrix}\in \mathbb{C}^{N\times N}, \quad N = kn-\ell. \end{align}\tag{35}\] By enforcing the normalization from Lemma 3, the matrix \(\mathbf{J}({\boldsymbol{\mu}})\) belongs to \(\mathcal{J}_{k,N}\) which shows that the correspondence \[\begin{align} \label{eq:InvSpecDef} {\boldsymbol{\mu}}\in \mathcal{M}_{k,N} ~ \overset{\psi}{\mapsto} ~\mathbf{J}({\boldsymbol{\mu}}) \in \mathcal{J}_{k,N}, \end{align}\tag{36}\] is a well-defined map.

Lemma 10. For a measure \({\boldsymbol{\mu}}\in \mathcal{M}_{k,N}\) with \(N = kn-\ell\), denote by \(\mathbf{J} = \psi({\boldsymbol{\mu}})\) the associated banded Hermitian matrix from \(\mathcal{J}_{k,N}\) and by \(\{\mathbf{P}_j(\diamond,{\boldsymbol{\mu}})\}_{j=0}^{n-1}\) the corresponding orthonormal polynomials defined in 12 and 16 . Define \[\begin{align} \mathcal{P}^*(x):=\mathcal{P}^*(x,{\boldsymbol{\mu}}) = \begin{bmatrix} \mathbf{P}_0(x,{\boldsymbol{\mu}}) & \mathbf{P}_1(x,{\boldsymbol{\mu}}) & \dots & \mathbf{P}_{n-2}(x,{\boldsymbol{\mu}}) & \mathbf{P}_{n-1}(x,{\boldsymbol{\mu}}) \end{bmatrix}, \end{align}\] then \(\mathcal{P}(x)\) satisfies \[\begin{align} \label{eq:RecMatForm} x \,\mathcal{P}(x) = \mathbf{J} \,\mathcal{P}(x) + \mathcal{R}(x), \end{align}\tag{37}\] where \(\mathcal{R}^*(x) := \mathcal{R}(x,{\boldsymbol{\mu}})^*\) is given by \[\begin{align} \mathcal{R}(x,{\boldsymbol{\mu}})^*= \begin{bmatrix} \mathbf{0}_{k\times k} & \dots & \mathbf{0}_{k\times k} & \mathbf{R}_{n-2}(x,{\boldsymbol{\mu}}) & \mathbf{R}_{n-1}(x,{\boldsymbol{\mu}}) \end{bmatrix}, \quad \mathbf{R}_{n-2}(x)\in \mathbb{C}^{k\times k}, \quad \mathbf{R}_{n-1}(x) \in \mathbb{C}^{k\times \ell}, \end{align}\] and satisfies \[\begin{align} \langle \mathbf{R}_{n-2}(\diamond,{\boldsymbol{\mu}}) ,\mathbf{R}_{n-2}(\diamond,{\boldsymbol{\mu}}) \rangle_{\boldsymbol{\mu}}= \mathbf{0}, \qquad \langle \mathbf{R}_{n-1}(\diamond,{\boldsymbol{\mu}}) ,\mathbf{R}_{n-1}(\diamond,{\boldsymbol{\mu}}) \rangle_{\boldsymbol{\mu}}= \mathbf{0}. \end{align}\]

Proof. The recurrence relations satisfied by the orthonormal polynomial \(\{\mathbf{P}_j\}_{j=0}^{n-2}\) in 13 imply that the first \(n-2\) blocks of \(\mathcal{R}\) are zero. Moreover, from the structure of \(\mathbf{J}\), we have \[\begin{align} \mathbf{R}_{n-2}(x) = x\, \mathbf{P}_{n-2}(x)-\mathbf{P}_{n-3}(x) \mathbf{B}_{n-3}^*-\mathbf{P}_{n-2}(x)\mathbf{A}_{n-2}-\mathbf{P}_{n-1}(x) \mathbf{B}_{n-2}, \end{align}\] and \[\begin{align} \mathbf{R}_{n-1}(x) = x\,\mathbf{P}_{n-1}(x) - \mathbf{P}_{n-2}(x)\mathbf{B}_{n-2}^*-\mathbf{P}_{n-1}(x)\mathbf{A}_{n-1}. \end{align}\] Using the definition of \(\mathbf{P}(x)\) in Theorem 4, we find \[\begin{align} \mathbf{R}_{n-2}(x) = \mathbf{P}(x) - \mathbf{P}_{n-1}(x)\mathbf{B}_{n-2} = \mathbf{P}(x)\left(\mathbf{I}_k-\mathbf{B}_{n-2}^\dagger\mathbf{B}_{n-2} \right), \end{align}\] where the last equality follows from \(\mathbf{P}_{n-1}(x) = \mathbf{P}(x) \mathbf{B}_{n-2}^\dagger\). Finally, since \(\langle \mathbf{P},\mathbf{P}\rangle_{\boldsymbol{\mu}}= \mathbf{B}_{n-2}^*\mathbf{B}_{n-2}\), we have \[\begin{align} \langle \mathbf{R}_{n-2}, \mathbf{R}_{n-2} \rangle_{\boldsymbol{\mu}}= (\mathbf{I}_k-\mathbf{B}_{n-2}^\dagger \mathbf{B}_{n-2})^*\langle \mathbf{P},\mathbf{P} \rangle_{\boldsymbol{\mu}}(\mathbf{I}_k-\mathbf{B}_{n-2}^\dagger \mathbf{B}_{n-2}) = \mathbf{0}. \end{align}\]

It remains to show that \(\mathbf{R}_{n-1}(x)\) has zero norm. Recall from Corollary 1 that \(\langle\mathbf{\Pi}_n,\mathbf{\Pi}_n \rangle_{\boldsymbol{\mu}}=\mathbf{0}\), where \(\mathbf{\Pi}_n(x)\) is the \(n\)-th monic orthogonal polynomial. Substituting the three-term recurrence from Theorem 3 gives \[\begin{align} \label{eq:NormxPi} \langle \diamond\mathbf{\Pi}_{n-1}(\diamond),\diamond\mathbf{\Pi}_{n-1}(\diamond) \rangle_{\boldsymbol{\mu}}= {\boldsymbol{\gamma}}_{n-1}\mathbf{C}_{n-1}^2 + \mathbf{D}_{n-1}^* {\boldsymbol{\gamma}}_{n-2}\mathbf{D}_{n-1}, \end{align}\tag{38}\] where \({\boldsymbol{\gamma}}_j\) is defined in 10 and \(\mathbf{C}_{n-1}, \mathbf{D}_{n-1}\) are as stated in Theorem 3. To connect \(\mathbf{\Pi}_{n-1}\) and \(\mathbf{P}_{n-1}\), we decompose \({\boldsymbol{\gamma}}_{n-1}\) as \[\begin{align} {\boldsymbol{\gamma}}_{n-1} = \left(\mathbf{U}_{n-1} \mathbf{\Lambda}_{n-1}^{1/2}\right)\left(\mathbf{U}_{n-1} \mathbf{\Lambda}_{n-1}^{1/2}\right)^*, \end{align}\] where \(\mathbf{\Lambda}_{n-1}\) is a diagonal matrix containing the nonzero eigenvalues of \({\boldsymbol{\gamma}}_{n-1}\), and \(\mathbf{U}_{n-1}\) is a matrix whose columns consist of the corresponding eigenvectors. This gives the normalization \[\begin{align} \label{eq:LastNormalization} \mathbf{P}_{n-1}(x) = \mathbf{\Pi}_{n-1}(x) \mathbf{U}_{n-1}\mathbf{\Lambda}_{n-1}^{-1/2}\mathbf{Q}_{n-1}, \end{align}\tag{39}\] which follows the same structure as the normalizations in 12 , for instance, \[\begin{align} \mathbf{P}_{n-2}(x) = \mathbf{\Pi}_{n-2}(x) {\boldsymbol{\gamma}}_{n-2}^{-1/2}\mathbf{Q}_{n-2}. \end{align}\] Here, \(\mathbf{Q}_{n-2}\) is the unique unitary matrix determined by Lemma 3 and \(\mathbf{Q}_{n-1}\) is the unitary matrix that guarantees that \(\mathbf{B}_{n-2}\) is in row echelon form with positive pivots. In other words, \(\mathbf{Q}_{n-2}\) and \(\mathbf{Q}_{n-1}\) are unitary matrices that satisfy \[\label{eq:ABrel} \begin{align} \mathbf{A}_{n-1} &= \langle \diamond\mathbf{P}_{n-1}(\diamond),\mathbf{P}_{n-1} \rangle_{\boldsymbol{\mu}}\\ &= \mathbf{Q}_{n-1}^* \mathbf{\Lambda}_{n-1}^{-1/2}\mathbf{U}_{n-1}^* \langle \diamond\mathbf{\Pi}_{n-1}(\diamond),\mathbf{\Pi}_{n-1}\rangle_{\boldsymbol{\mu}}\mathbf{U}_{n-1}\mathbf{\Lambda}_{n-1}^{-1/2} \mathbf{Q}_{n-1} \\ &= \mathbf{Q}_{n-1}^* \mathbf{\Lambda}_{n-1}^{1/2}\mathbf{U}_{n-1}^*\mathbf{C}_{n-1} \mathbf{U}_{n-1}\mathbf{\Lambda}_{n-1}^{-1/2} \mathbf{Q}_{n-1},\\ \mathbf{B}_{n-2} &= \langle \diamond \mathbf{P}_{n-1}(\diamond),\mathbf{P}_{n-2} \rangle_{\boldsymbol{\mu}}\\ &= \mathbf{Q}_{n-1}^* \mathbf{\Lambda}_{n-1}^{-1/2}\mathbf{U}_{n-1}^*\langle \diamond\mathbf{\Pi}_{n-1}(\diamond),\mathbf{\Pi}_{n-2} \rangle_{\boldsymbol{\mu}}{\boldsymbol{\gamma}}_{n-2}^{-1/2} \mathbf{Q}_{n-2} \\ &= \mathbf{Q}_{n-1}^* \mathbf{\Lambda}_{n-1}^{-1/2}\mathbf{U}_{n-1}^*\mathbf{D}_{n-1}^* {\boldsymbol{\gamma}}_{n-2}^{1/2}\mathbf{Q}_{n-2}. \end{align}\tag{40}\] Combining these identities with 38 and 39 , gives \[\label{eq:xPnm1} \begin{align} \langle \diamond\mathbf{P}_{n-1}(\diamond), \diamond \mathbf{P}_{n-1}(\diamond) \rangle_{\boldsymbol{\mu}}&= \mathbf{Q}_{n-1}^* \mathbf{\Lambda}_{n-1}^{-1/2}\mathbf{U}_{n-1}^*\langle \diamond\mathbf{\Pi}_{n-1}(\diamond),\diamond\mathbf{\Pi}_{n-1}(\diamond)\rangle_{\boldsymbol{\mu}}\mathbf{U}_{n-1}\mathbf{\Lambda}_{n-1}^{-1/2} \mathbf{Q}_{n-1}\\ &= \mathbf{Q}_{n-1}^* \mathbf{\Lambda}_{n-1}^{-1/2}\mathbf{U}_{n-1}^* \left( {\boldsymbol{\gamma}}_{n-1}\mathbf{C}_{n-1}^2 + \mathbf{D}_{n-1}^* {\boldsymbol{\gamma}}_{n-2}\mathbf{D}_{n-1} \right) \mathbf{U}_{n-1}\mathbf{\Lambda}_{n-1}^{-1/2} \mathbf{Q}_{n-1} \\ &= \mathbf{A}_{n-1}^2 + \mathbf{B}_{n-2} \mathbf{B}_{n-2}^*. \end{align}\tag{41}\] A direct computation using 41 gives \(\langle \mathbf{R}_{n-1},\mathbf{R}_{n-1} \rangle_{\boldsymbol{\mu}}= \mathbf{0}\), which concludes the proof. ◻

Theorem 11. Let \({\boldsymbol{\mu}}\in \mathcal{M}_{k,N}\), then \({\boldsymbol{\mu}}= \varphi(\psi({\boldsymbol{\mu}}))\).

Proof. Let \({\boldsymbol{\mu}}\in \mathcal{M}_{k,N}\) and write \({\boldsymbol{\mu}}= \sum_{j=1}^m \mathbf{V}_j \mathbf{V}_j^* \delta_{x_j}\). Using Lemma 10 and evaluating 37 at \(x_j\), followed by a multiplication on the right by \(\mathbf{V}_j\), we find for \(\mathbf{J} = \psi({\boldsymbol{\mu}})\), \[\begin{align} x_j \mathcal{P}(x_j)\mathbf{V}_j = \mathbf{J} \mathcal{P}(x_j) \mathbf{V}_j \qquad j=1,\dots,m. \end{align}\] Since \(\mathbf{V}_j\) has full rank \(n_j\), the resulting vectors \(\{\mathcal{P} \, \mathbf{V}_j \, \mathbf{e}_i\}_{i=1}^{n_j}\) are linearly independent, proving that the eigenvalues of \(\mathbf{J} = \psi({\boldsymbol{\mu}})\) coincide with \(x_j\) and that the first \(k\) rows of the corresponding eigenvectors are exactly \(\mathbf{V}_j\). These quantities uniquely determine the measure associated with \(\mathbf{J}\) via the mapping \(\varphi\), and this measure coincides with \({\boldsymbol{\mu}}\), i.e. \({\boldsymbol{\mu}}= \varphi(\mathbf{J})\). ◻

Collecting the results of this section, Theorems 7 and 8 show that \(\mathrm{Ran}~\varphi \subset \mathcal{M}_{k,N}\), where \(\mathcal{M}_{k,N}\) is precisely the domain on which \(\psi\) is well defined. Theorem 11 further shows that \(\varphi\) is surjective onto \(\mathcal{M}_{k,N}\), i.e. \(\mathrm{Ran}~\varphi = \mathcal{M}_{k,N}\), and that \(\varphi\) is a left-inverse of \(\psi\). Theorem 10 establishes that \(\varphi\) is injective, and together these properties imply that \(\varphi\) is also a right-inverse. Indeed, for any \(\mathbf{J} \in \mathcal{J}_{k,N}\), applying Theorem 11 with \({\boldsymbol{\mu}}= \varphi(\mathbf{J})\) gives \(\varphi(\mathbf{J}) = \varphi\left(\psi(\varphi(\mathbf{J}))\right)\). By injectivity of \(\varphi\), it follows that \(\psi(\varphi(\mathbf{J})) = \mathbf{J}\). Therefore, we conclude that \(\varphi\colon \mathcal{J}_{k,N}\to \mathcal{M}_{k,N}\) is a bijection and \(\psi = \varphi^{-1}\). As a consequence of Theorems 8 [item:Nd-meas] and 11, we also obtain a full characterization of measures associated with matrices in \(\mathcal{J}_{k,N}\).

Corollary 2. Let \({\boldsymbol{\mu}}= \sum_{j=1}^m \mathbf{V}_j\mathbf{V}_j^* \delta_{x_j}\) where the points \(\{x_j\}_{j=1}^m\) are distinct and each \(\mathbf{V}_j\in \mathbb{C}^{k}\times n_j\) with \(n_j\leq k\) and \(\sum_j n_j = N\). Define \[\begin{align} \mathbf{X} = \mathrm{diag}(\underbrace{x_1,\dots,x_1}_{n_1},\underbrace{x_2,\dots,x_2}_{n_2},\ldots,\underbrace{x_m,\dots,x_m}_{n_m}), \quad \mathbf{V}^* = \begin{bmatrix} \mathbf{V}_1 & \cdots & \mathbf{V}_m \end{bmatrix}, \quad \mathbf{V} \in \mathbb{C}^{N \times k}. \end{align}\] The measure \({\boldsymbol{\mu}}\) is the spectral measure of a matrix in \(\mathcal{J}_{k,N}\) iff the matrix \[\begin{align} \begin{bmatrix} \mathbf{V} & \mathbf{X} \mathbf{V} & \dots & \mathbf{X}^d \mathbf{V} \end{bmatrix} \in \mathbb{C}^{N\times (d+1)k}, \end{align}\] is full rank.

4 Toda Flow on Banded Hermitian Matrices↩︎

In this section, we extend the Toda flow from Section 1.3 to the class of banded Hermitian matrices \(\mathcal{J}_{k,N}\) defined in [def:JkN]. We show that the fundamental properties of the classical Toda lattice persist in this broader setting, and that the corresponding matrix-valued spectral measure admits a similar evolution to equations 8 and 8 .

First, we show that the solution of the Toda flow 5 can be expressed in terms of a matrix exponential involving the initial condition. For completeness, we include a proof following the approach described in [61]. We begin by stating a lemma that sets the notation and introduces an important element of the proof.

Lemma 11. Every square matrix \(\mathbf{A}\) admits a unique decomposition of the form \[\begin{align} \mathbf{A} = \pi_S(\mathbf{A}) + \pi_U(\mathbf{A}), \end{align}\] where \(\pi_S(\mathbf{A})\) is skew-Hermitian and \(\pi_U(\mathbf{A})\) is upper triangular.

Proof. The existence argument follows directly by taking \[\begin{align} (\pi_S(\mathbf{A}))_{ij} = \begin{cases} \mathbf{A}_{ij} & \text{if } j > i, \\ 0 & \text{if } i = j, \\ -\mathbf{A}_{ji}^* & \text{if } j < i, \end{cases} \quad \text{and} \quad (\pi_U(\mathbf{A}))_{ij} = \begin{cases} 0 & \text{if } j > i, \\ \mathbf{A}_{ij} & \text{if } i = j, \\ \mathbf{A}_{ij} + \mathbf{A}_{ji}^* & \text{if } j < i. \end{cases} \end{align}\] For uniqueness, suppose \(\mathbf{A} = \mathbf{S}_1 + \mathbf{U}_1 = \mathbf{S}_2 + \mathbf{U}_2\), where \(\mathbf{S}_1, \mathbf{S}_2\) are skew-Hermitian and \(\mathbf{U}_1, \mathbf{U}_2\) are upper triangular. It follows that \[\begin{align} \mathbf{S}_1 - \mathbf{S}_2 = \mathbf{U}_2 - \mathbf{U}_1. \end{align}\] The left-hand side is skew-Hermitian, while the right-hand side is upper triangular. The only matrix that is both skew-Hermitian and upper triangular is the zero matrix, which shows that the decomposition is unique. ◻

Theorem 12. Let \(\mathbf{X}_0\in \mathbb{C}^{N\times N}\) be a Hermitian matrix, and let \(\mathbf{X}(t)\in\mathbb{C}^{N\times N}\) denote the solution to the Toda flow 5 with initial condition \(\mathbf{X}_0\). The solution \(\mathbf{X}(t)\) can be expressed as \[\begin{align} \label{eq:todaSol} \mathbf{X}(t) = \mathbf{Q}^*(t) \mathbf{X}_0 \mathbf{Q}(t) = \mathbf{R}(t) \mathbf{X}_0 \mathbf{R}^{-1}(t), \end{align}\tag{42}\] where \(\mathbf{Q}(t) \mathbf{R}(t)\) is the QR decomposition of \(\exp(t\mathbf{X}_0)\), with \(\mathbf{Q}(t)\) unitary and \(\mathbf{R}(t)\) upper triangular with positive diagonal entries.

Proof. We first establish the equality \(\mathbf{X} = \mathbf{Q}^* \mathbf{X}_0 \mathbf{Q}\). Since \(\exp(t\mathbf{X}_0)\) depends smoothly on \(t\), its QR-decomposition is differentiable [62], [63]. Differentiating \(\exp(t\mathbf{X}_0) = \mathbf{Q}\mathbf{R}\) gives \[\begin{align} \label{eq:DerivQRToda} \partial_t \mathbf{Q} \mathbf{R} + \mathbf{Q} \partial_t \mathbf{R} = \mathbf{X}_0 \mathbf{Q} \mathbf{R}. \end{align}\tag{43}\] Multiplying on the left by \(\mathbf{Q}^*\) and on the right by \(\mathbf{R}^{-1}\), we get \[\begin{align} \mathbf{Q}^* \partial_t \mathbf{Q} + \partial_t \mathbf{R} \mathbf{R}^{-1} = \mathbf{Q}^* \mathbf{X}_0 \mathbf{Q}. \end{align}\] Note that \(\mathbf{Q}^*\mathbf{Q}=\mathbf{I}_N\), which implies \[\begin{align} \partial_t \mathbf{Q}^*\mathbf{Q} = \mathbf{Q}^*\partial_t \mathbf{Q} = 0, \end{align}\] so \(\mathbf{Q}^*\partial_t \mathbf{Q}\) is skew-Hermitian and thus \[\begin{align} \pi_{S}\left(\mathbf{Q}^*\mathbf{X}_0\mathbf{Q}\right) = \mathbf{Q}^*\partial_t \mathbf{Q}, \quad \text{and} \quad \pi_{U}\left(\mathbf{Q}^*\mathbf{X}_0\mathbf{Q}\right) = \partial_t \mathbf{R}\mathbf{R}^{-1}. \end{align}\]

Now, define \(\tilde{\mathbf{X}} = \mathbf{Q}^*\mathbf{X}_0\mathbf{Q}\) and note that \(\tilde{\mathbf{X}}(0) = \mathbf{X}_0\) and \[\begin{align} \partial_t\tilde{\mathbf{X}} &= \partial_t \mathbf{Q}^* \mathbf{X}_0 \mathbf{Q}+\mathbf{Q}^*\mathbf{X}_0\partial_t \mathbf{Q}\\ &= \left(\partial_t \mathbf{Q}^*\mathbf{Q}\right)\mathbf{Q}^*\mathbf{X}_0\mathbf{Q} + \mathbf{Q}^*\mathbf{X}_0\mathbf{Q}\left(\mathbf{Q}^*\partial_t \mathbf{Q}\right)\\ &= \mathbf{Q}^*\mathbf{X}_0\mathbf{Q}\left(\mathbf{Q}^*\partial_t \mathbf{Q}\right) -\left(\mathbf{Q}^*\partial_t \mathbf{Q}\right)\mathbf{Q}^*\mathbf{X}_0 \mathbf{Q}\\ &= [\tilde{\mathbf{X}},\pi_S(\tilde{\mathbf{X}})]. \end{align}\] On the other hand, we have \(\mathbf{X} = (\mathbf{X}-\mathbf{B}(\mathbf{X})) + \mathbf{B}(\mathbf{X})\) which implies \(\pi_U(\mathbf{X}) = \mathbf{X}-\mathbf{B}(\mathbf{X})\) and \(\pi_S(\mathbf{X}) = \mathbf{B}(\mathbf{X})\). Substituting into 5 gives \[\begin{align} \partial_t \mathbf{X} = [\mathbf{X},\pi_S(\mathbf{X})]. \end{align}\] Thus, both \(\mathbf{X}\) and \(\tilde{\mathbf{X}}\) satisfy the same differential equation with identical initial condition which establishes the first equality in 42 . The second equality follows from the fact that \(\mathbf{X}_0\) commutes with \(\exp(t\mathbf{X}_0) = \mathbf{Q}\mathbf{R}\). Specifically, we have \(\mathbf{X}_0 \mathbf{Q} \mathbf{R} = \mathbf{Q}\mathbf{R}\mathbf{X}_0\) which gives \(\mathbf{Q}^*\mathbf{X}_0\mathbf{Q} = \mathbf{R}\mathbf{X}_0\mathbf{R}^{-1},\) thereby completing the proof. ◻

Corollary 3. The Toda flow defines an isospectral flow on matrices in \(\mathcal{J}_{k,N}\). Moreover, it preserves the band size and the structure of \(\mathbf{X}_0\), that is \(\mathbf{X}(t)\in \mathcal{J}_{k,N}\) for all \(t\).

Proof. The first equality \(\mathbf{X} = \mathbf{Q}^* \mathbf{X}_0 \mathbf{Q}\) in Theorem 12 shows that \(\mathbf{X}\) and \(\mathbf{X}_0\) are similar, so 5 is an isospectral flow. The structure of the solution also guarantees that \(\mathbf{X}\) remains Hermitian, and the second equality in 42 implies that \(\mathbf{X}_{ij} = 0\) whenever \(i - j > k\). Together, these properties show that \(X(t)\) is banded with bandwidth \(k\) for all \(t\).

A closer look at \(\mathbf{X} = \mathbf{R} \mathbf{X}_0\mathbf{R}^{-1}\) shows that the off-diagonal blocks \(\{\mathbf{B}_j\}_{j=0}^{n-2}\) of \(\mathbf{X}\) satisfy \(\mathbf{B}_j = \mathbf{R}_j \mathbf{B}_j(0)\mathbf{U}_j\) where \(\mathbf{R}_j\) and \(\mathbf{U}_j\) are upper triangular matrices with positive diagonal entries. Let \(p_i\) denote the pivot column of row \(i\) in \(\mathbf{B}_j(0)\), so that \((\mathbf{B}_j(0))_{i,\ell} = 0\) for \(\ell < p_i\) and \(p_1 < p_2 < \cdots\). For \(\ell < p_i\), we have \[\begin{align} (\mathbf{R}_j \mathbf{B}_j(0))_{i,\ell} = \sum_{k \geq i} (\mathbf{R}_j)_{i,k} (\mathbf{B}_j(0))_{k,\ell} = 0, \end{align}\] since \((\mathbf{B}_j(0))_{k,\ell} = 0\) for all \(k \geq i\), as \(p_k \geq p_i > \ell\). A similar argument shows that right multiplication by \(\mathbf{U}_j\) preserves zeros to the left of each pivot. Thus \(\mathbf{B}_j\) has the same pivot structure as \(\mathbf{B}_j(0)\), and the pivot entries satisfy \[\begin{align} (\mathbf{B}_j)_{i,\,p_i} = (\mathbf{R}_j)_{i,i}\,(\mathbf{B}_j(0))_{i,\,p_i}\,(\mathbf{U}_j)_{p_i,\,p_i} > 0. \end{align}\] We conclude that \(\mathbf{X} \in \mathcal{J}_{k,N}\), completing the proof. ◻

The next result describes the time evolution of the spectral measure associated with \(\mathbf{X}(t)\). It shows that the unnormalized weights evolve by an exponential scaling of the initial weights, after which normalization is achieved through the matrix \(\mathbf{L}(t)\).

Theorem 13. Suppose that \(\mathbf{X}(t)\) is a solution to the Toda flow, \(\mathbf{X}(0) \in \mathcal{J}_{N,k}\), and denote the spectral measure of \(\mathbf{X}(t)\) by \[\begin{align} {\boldsymbol{\mu}}_{\mathbf{X}(t)} = \sum_{j=1}^m \mathbf{V}_j(t) \mathbf{V}_j^*(t)\delta_{\lambda_j}, \quad \mathbf{V}_j(t)\in \mathbb{R}^{k\times \ell_j}. \end{align}\] Then \[\begin{align} \mathbf{V}_j(t) \mathbf{V}_j^*(t) = \mathbf{L}^{-1}(t)\left( e^{2\lambda_j t} \mathbf{V}_j(0)\mathbf{V}_j(0)^*\right) \mathbf{L}^{-*}(t), \end{align}\] where \(\mathbf{L}^{-1}(t) = \mathbf{I}_{N\times k}^* \mathbf{R}^{-*}(t)\mathbf{I}_{N\times k}\) and \(\mathbf{R}(t)\) is defined in Theorem 12. Equivalently, \(\mathbf{L}(t)\) is the lower triangular matrix satisfying \(\sum_{j} e^{2\lambda_j t} \mathbf{V}_j(0)\mathbf{V}_j(0)^* = \mathbf{L}(t) \mathbf{L}(t)^*\).

Proof. Let \(\mathbf{X} = \mathbf{U} \mathbf{\Lambda} \mathbf{U}^*\) be the eigendecomposition of \(\mathbf{X}\), and let \(\mathbf{S}_j\in \mathbb{R}^{N\times \ell_j}\) be a column selection matrix that extracts the columns of \(\mathbf{U}\) corresponding to the eigenvalue \(\lambda_j\). In other words, \(\mathbf{S}_j\) satisfies \[\begin{align} \mathbf{X} \mathbf{U} \mathbf{S}_j = \lambda_j \mathbf{U} \mathbf{S}_j. \end{align}\] The evolution of the eigenvectors follows directly from equation 5 . In particular, \[\begin{align} \mathbf{U} = \mathbf{Q}^* \mathbf{U}(0) \mathbf{O}, \end{align}\] where \(\mathbf{O}(t)\) is a unitary matrix such that \(\mathbf{O}_j(t) = \mathbf{S}_j^T\mathbf{O}(t)\mathbf{S}_j\) is also unitary. Since \(\mathbf{V}_j\) is given by the first \(k\) rows of the corresponding eigenvectors, i.e. \(\mathbf{V}_j = \mathbf{I}_{N\times k}^* \mathbf{U} \mathbf{S}_j\), we have \[\begin{align} \label{eq:BandMeasPfEq1} \mathbf{V}_j = \mathbf{I}_{N\times k}^*\mathbf{Q}^* \mathbf{U}(0) \mathbf{O} \mathbf{S}_j. \end{align}\tag{44}\]

Recall that \(\exp(t\mathbf{X}_0) = \mathbf{Q} \mathbf{R}\) where \(\mathbf{Q}\) is unitary and \(\mathbf{R}\) is upper triangular, so \(\mathbf{Q}^* = \mathbf{R}^{-*}\exp(t\mathbf{X}_0)\) which implies \[\begin{align} \label{eq:BandMeasPfEq2} \mathbf{I}_{N\times k}^* \mathbf{Q}^* = \mathbf{I}_{N\times k}^* \mathbf{R}^{-*}\exp(t\mathbf{X}_0) = \mathbf{L}^{-1} \mathbf{I}_{N\times k}^* \mathbf{U}(0)\exp(t\mathbf{\Lambda})\mathbf{U}(0)^*, \quad \text{where} \quad \mathbf{L}^{-1} = \mathbf{I}_{N\times k}^*\mathbf{R}^{-*}\mathbf{I}_{N\times k}. \end{align}\tag{45}\] Combining both 44 and 45 with the fact that \(\mathbf{O} \mathbf{S}_j = \mathbf{S}_j \mathbf{O}_j\), we find \[\begin{align} \mathbf{V}_j = \mathbf{L}^{-1} \mathbf{I}_{N\times k}^* \mathbf{U}(0)\exp(t\mathbf{\Lambda}) \mathbf{S}_j \mathbf{O}_j = \mathbf{L}^{-1} \left(e^{\lambda_j t}\mathbf{U}_k(0)\mathbf{S}_j\right) \mathbf{O}_j = \mathbf{L}^{-1} \left(e^{\lambda_j t }\mathbf{V}_j(0)\right) \mathbf{O}_j, \end{align}\] and therefore \[\begin{align} \mathbf{V}_j \mathbf{V}_j^* =\mathbf{L}^{-1} \left( e^{2\lambda_j t} \mathbf{V}_j(0)\mathbf{V}_j^*(0)\right) \mathbf{L}^{-*}, \end{align}\] which completes our proof. ◻

Remark 14. The result in Theorem 13 holds for any Hermitian initial data, however its importance is when \(\mathbf{X}(0) \in \mathcal{J}_{k,N}\), since in that case \({\boldsymbol{\mu}}_{\mathbf{X}}\) can be used to recover the solution \(\mathbf{X}\) using the inverse spectral theory developed in Section 3.4.

Remark 15. When \(k=1\), Theorem 13 reduces to the classical evolution of the spectral measure for Jacobi matrices. In this case, each \(\mathbf{V}_j(t)\) becomes a scalar \(\mathbf{v}_{1,j}(t)\), and \(\mathbf{L}(t)\) simplifies to \[\begin{align} L(t) = \Bigg(\sum_{i=1}^N e^{2\lambda_i t} v_{1,j}^2(0)\Bigg)^{1/2}. \end{align}\] Substituting into the formula of Theorem 13, we find the standard expression in 8 .

Acknowledgements↩︎

The authors thank Maxim Yattselev for pointing out several relevant references. This material is based upon work supported by NSF DMS-2306438 (TT). Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author and do not necessarily reflect the views of the National Science Foundation.

5 Block Lanczos and Householder↩︎

Figure 2: Householder reduction to block tridiagonal form

References↩︎

[1]

T. Chen and T. Trogdon, “Stability of the Lanczos algorithm on matrices with regular spectral distributions,” Linear Algebra and its Applications, vol. 682, pp. 191–237, 2024, doi: 10.1016/j.laa.2023.11.006.

[2]

C. A. Younes, X. Ding, and T. Trogdon, arXiv preprint 2504.03066, 2025“A Lanczos-Based Algorithmic Approach for Spike Detection in Large Sample Covariance Matrices.” arXiv, 2025, doi: 10.48550/arXiv.2504.03066.

[3]

X. Ding and T. Trogdon, “The conjugate gradient algorithm on a general class of spiked covariance matrices,” Quarterly of Applied Mathematics, vol. 80, no. 1, pp. 99–155, 2022, doi: 10.1090/qam/1605.

[4]

P. Deift and T. Trogdon, “The conjugate gradient algorithm on well-conditioned Wishart matrices is almost deterministic,” Quarterly of Applied Mathematics, vol. 79, no. 1, pp. 125–161, 2021, doi: 10.1090/qam/1574.

[5]

E. Paquette and T. Trogdon, “Universality for the Conjugate Gradient and MINRES Algorithms on Sample Covariance Matrices,” Communications on Pure and Applied Mathematics, vol. 76, no. 5, pp. 1085–1136, 2023, doi: 10.1002/cpa.22081.

[6]

X. Ding and T. Trogdon, “A Riemann–Hilbert Approach to the Perturbation Theory for Orthogonal Polynomials: Applications to Numerical Linear Algebra and Random Matrix Theory,” International Mathematics Research Notices, vol. 2024, no. 5, pp. 3975–4061, 2024, doi: 10.1093/imrn/rnad142.

[7]

P. Deift, Orthogonal Polynomials and Random Matrices: A Riemann-Hilbert Approach: A Riemann-Hilbert Approach. American Mathematical Soc., 2000.

[8]

F. Gesztesy and B. Simon, “M-Functions and inverse spectral analysis for finite and semi-infinite Jacobi matrices,” Journal d’Analyse Mathamétique, vol. 73, no. 1, pp. 267–297, 1997, doi: 10.1007/BF02788147.

[9]

Yu. M. Berezanskiı̆, Expansions in eigenfunctions of selfadjoint operators, vol. 17. American Mathematical Society, 1968.

[10]

F. V. Atkinson, “Discrete and Continuous Boundary Problems.” Journal of Applied Mathematics and Mechanics, vol. 46, no. 5, pp. 327–327, 1966, doi: 10.1002/zamm.19660460520.

[11]

A. J. Antony and M. Krishna, “Inverse spectral theory for Jacobi matrices and their almost periodicity,” Proceedings - Mathematical Sciences, vol. 104, no. 4, pp. 777–818, 1994, doi: 10.1007/BF02830803.

[12]

F. Gesztesy, M. Krishna, and G. Teschl, “On isospectral sets of Jacobi operators,” Communications in Mathematical Physics, vol. 181, no. 3, pp. 631–645, 1996, doi: 10.1007/BF02101290.

[13]

H. Hochstadt, “On the construction of a Jacobi matrix from spectral data,” Linear Algebra and its Applications, vol. 8, no. 5, pp. 435–446, 1974, doi: 10.1016/0024-3795(74)90077-9.

[14]

H. Hochstadt, “On the construction of a Jacobi matrix from mixed given data,” Linear Algebra and its Applications, vol. 28, pp. 113–115, 1979, doi: 10.1016/0024-3795(79)90124-1.

[15]

P. Deift and T. Nanda, “On the determination of a tridiagonal matrix from its spectrum and a submatrix,” Linear Algebra and its Applications, vol. 60, pp. 43–55, 1984, doi: 10.1016/0024-3795(84)90069-7.

[16]

D. R. Masson and J. Repka, “Spectral Theory of Jacobi Matrices in \(\ell^2(\mathbb{Z})\) and the \(su(1,1)\) Lie Algebra,” SIAM Journal on Mathematical Analysis, vol. 22, no. 4, pp. 1131–1146, 1991, doi: 10.1137/0522073.

[17]

W. E. Ferguson, “The construction of Jacobi and periodic Jacobi matrices with prescribed spectra,” Mathematics of Computation, vol. 35, no. 152, pp. 1203–1220, 1980, doi: 10.1090/S0025-5718-1980-0583498-3.

[18]

C. de Boor and G. H. Golub, “The numerically stable reconstruction of a Jacobi matrix from spectral data,” Linear Algebra and its Applications, vol. 21, no. 3, pp. 245–260, 1978, doi: 10.1016/0024-3795(78)90086-1.

[19]

G. Teschl, “Trace Formulas and Inverse Spectral Theory for Jacobi Operators,” Communications in Mathematical Physics, vol. 196, no. 1, pp. 175–202, 1998, doi: 10.1007/s002200050419.

[20]

N. I. Akhiezer, The Classical Moment Problem and Some Related Questions in Analysis. Society for Industrial; Applied Mathematics, 2020.

[21]

S. M. Zagorodnyuk, “The direct and inverse spectral problems for (2N+1)-diagonal complex transposition-antisymmetric matrices,” Methods of Functional Analysis and Topology, vol. 14, no. 2, pp. 124–131, 2008, Accessed: Nov. 13, 2025. [Online]. Available: https://mfat.imath.kiev.ua/article/?id=450.

[22]

M. A. Kudryavtsev, “The direct and the inverse problem of spectral analysis for five-diagonal symmetric matrices. I.” Mat. Fiz. Anal. Geom., vol. 5, no. 3–4, pp. 182–202, 1998.

[23]

M. A. Kudryavtsev, “The direct and the inverse problem of spectral analysis for five-diagonal symmetric matrices. II.” Mat. Fiz. Anal. Geom., vol. 6, no. 1–2, pp. 55–80, 1999.

[24]

V. Marchenko and V. Slavin, Inverse Problems in the Theory of Small Oscillations. American Mathematical Soc., 2018.

[25]

M. Kudryavtsev, S. Palafox, and L. O. Silva, arXiv preprint 1409.3868, 2017“Inverse spectral analysis for a class of finite band symmetric matrices.” arXiv, 2017, doi: 10.48550/arXiv.1409.3868.

[26]

M. Kudryavtsev, S. Palafox, and L. O. Silva, “Inverse spectral analysis for a class of infinite band symmetric matrices,” Journal of Mathematical Analysis and Applications, vol. 445, no. 1, pp. 762–783, 2017, doi: 10.1016/j.jmaa.2016.07.057.

[27]

A. Branquinho, A. Foulquié-Moreno, and M. Mañas, “Spectral theory for bounded banded matrices with positive bidiagonal factorization and mixed multiple orthogonal polynomials,” Advances in Mathematics, vol. 434, p. 109313, 2023, doi: 10.1016/j.aim.2023.109313.

[28]

D. Damanik, A. Pushnitski, and B. Simon, arXiv preprint 0711.2703, 2008“The Analytic Theory of Matrix Orthogonal Polynomials.” arXiv, 2014, doi: 10.48550/arXiv.0711.2703.

[29]

F. W. Biegler-König, “Construction of band matrices from spectral data,” Linear Algebra and its Applications, vol. 40, pp. 79–87, 1981, doi: 10.1016/0024-3795(81)90141-5.

[30]

B. Beckermann and A. Osipov, “Some Spectral Properties of Infinite Band Matrices,” Numerical Algorithms, vol. 34, no. 2, pp. 173–185, 2003, doi: 10.1023/B:NUMA.0000005361.17723.a4.

[31]

M. P. Mattis and H. Hochstadt, “On the construction of band matrices from spectral data,” Linear Algebra and its Applications, vol. 38, pp. 109–119, 1981, doi: 10.1016/0024-3795(81)90012-4.

[32]

A. J. Duran, “A Generalization of Favard’s Theorem for Polynomials Satisfying a Recurrence Relation,” Journal of Approximation Theory, vol. 74, no. 1, pp. 83–109, 1993, doi: 10.1006/jath.1993.1055.

[33]

H. Dette and W. J. Studden, “Matrix measures, moment spaces and Favard’s theorem for the interval [0,1] and [0,\({\infty}\)),” Linear Algebra and its Applications, vol. 345, no. 1, pp. 169–193, 2002, doi: 10.1016/S0024-3795(01)00493-1.

[34]

M. Kudryavtsev, S. Palafox, and L. O. Silva, “On a linear interpolation problem for \(n\)-dimensional vector polynomials,” Journal of Approximation Theory, vol. 199, pp. 45–62, 2015, doi: 10.1016/j.jat.2015.06.006.

[35]

G. H. Golub and R. Underwood, “The Block Lanczos Method for Computing Eigenvalues,” in Mathematical Software, Academic Press, 1977, pp. 361–377.

[36]

R. R. Underwood, PhD thesis, ProQuest Publication No. AAI7525622“An iterative block lanczos method for the solution of large sparse symmetric eigenproblems,” PhD thesis, Stanford University, Stanford, CA, USA, 1975.

[37]

G. H. Golub and C. F. V. Loan, Matrix Computations. JHU Press, 2013.

[38]

C. Bischof and X. Sun, “On orthogonal block elimination,” Argonne National Lab., IL (United States). Mathematics; Computer Science Div., MCS-P–450-0794, 1996. Accessed: Dec. 19, 2025. [Online]. Available: https://www.osti.gov/biblio/438466.

[39]

J. J. Dongarra, D. C. Sorensen, and S. J. Hammarling, “Block reduction of matrices to condensed forms for eigenvalue computations,” Journal of Computational and Applied Mathematics, vol. 27, no. 1, pp. 215–227, 1989, doi: 10.1016/0377-0427(89)90367-1.

[40]

C. Bischof and C. Van Loan, “The WY Representation for Products of Householder Matrices,” SIAM Journal on Scientific and Statistical Computing, vol. 8, no. 1, pp. s2–s13, 1987, doi: 10.1137/0908009.

[41]

R. Schreiber and C. Van Loan, “A Storage-efficient \(WY\) Representation for Products of Householder Transformations,” SIAM Journal on Scientific and Statistical Computing, vol. 10, no. 1, pp. 53–57, 1989, doi: 10.1137/0910005.

[42]

R. Schreiber and B. Parlett, “Block Reflectors: Theory and Computation,” SIAM Journal on Numerical Analysis, vol. 25, no. 1, pp. 189–205, 1988, doi: 10.1137/0725014.

[43]

A. Bloemendal and B. Virág, “Limits of spiked random matrices II,” The Annals of Probability, vol. 44, no. 4, pp. 2726–2769, 2016, doi: 10.1214/15-AOP1033.

[44]

M. Toda, “Vibration of a Chain with Nonlinear Interaction,” Journal of the Physical Society of Japan, vol. 22, no. 2, pp. 431–436, 1967, doi: 10.1143/JPSJ.22.431.

[45]

M. Toda, Theory of Nonlinear Lattices, vol. 20. Berlin, Heidelberg: Springer, 1989.

[46]

H. Flaschka, “The Toda lattice. II. Existence of integrals,” Physical Review B, vol. 9, no. 4, pp. 1924–1925, 1974, doi: 10.1103/PhysRevB.9.1924.

[47]

S. V. Manakov, “Complete integrability and stochastization of discrete dynamical systems,” Soviet Journal of Experimental and Theoretical Physics, vol. 40, no. 2, pp. 269–274, 1975, Accessed: Nov. 27, 2025. [Online]. Available: https://ui.adsabs.harvard.edu/abs/1975JETP...40..269M/abstract.

[48]

T. Nanda, “Isospectral flow on band matrices,” Ph.D. thesis, New York University, 1982.

[49]

A. Sinap and W. Van Assche, “Orthogonal matrix polynomials and applications,” Journal of Computational and Applied Mathematics, vol. 66, no. 1, pp. 27–52, 1996, doi: 10.1016/0377-0427(95)00193-X.

[50]

M. G. Kreı̆n, “Fundamental aspects of the representation theory of Hermitian operators with deficiency index \((n_+, n_-)\),” in American mathematical society translations, series 2, vol. 97, Providence, Rhode Island: American Mathematical Society, 1971, pp. 75–143.

[51]

L. Miranian, “Matrix-valued orthogonal polynomials on the real line: Some extensions of the classical theory,” Journal of Physics A: Mathematical and General, vol. 38, no. 25, p. 5731, 2005, doi: 10.1088/0305-4470/38/25/009.

[52]

A. J. Duran, “On Orthogonal Polynomials With Respect to a Positive Definite Matrix of Measures,” Canadian Journal of Mathematics, vol. 47, no. 1, pp. 88–112, 1995, doi: 10.4153/CJM-1995-005-8.

[53]

A. Sinap and W. Van Assche, “Polynomial interpolation and Gaussian quadrature for matrix-valued functions,” Linear Algebra and its Applications, vol. 207, pp. 71–114, 1994, doi: 10.1016/0024-3795(94)90005-1.

[54]

M. D. Kent, “Chebyshev, Krylov, Lanczos: Matrix relationships and computations,” Ph.{D}., Stanford University, United States – California, 1989.

[55]

A. Casulli and L. Robol, “An Efficient Block Rational Krylov Solver for Sylvester Equations with Adaptive Pole Selection,” SIAM Journal on Scientific Computing, vol. 46, no. 2, pp. A798–A824, 2024, doi: 10.1137/23M1548463.

[56]

A. Frommer, K. Lund, and D. B. Szyld, “Block Krylov Subspace Methods for Functions of Matrices II: Modified Block FOM,” SIAM Journal on Matrix Analysis and Applications, vol. 41, no. 2, pp. 804–837, 2020, doi: 10.1137/19M1255847.

[57]

V. Simoncini and E. Gallopoulos, “Convergence properties of block GMRES and matrix polynomials,” Linear Algebra and its Applications, vol. 247, pp. 97–119, 1996, doi: 10.1016/0024-3795(95)00093-3.

[58]

V. Simoncini, “Ritz and Pseudo-Ritz values using matrix polynomials,” Linear Algebra and its Applications, vol. 241–243, pp. 787–801, 1996, doi: 10.1016/0024-3795(95)00682-6.

[59]

S. Elsworth and S. Güttel, “The Block Rational Arnoldi Method,” SIAM Journal on Matrix Analysis and Applications, vol. 41, no. 2, pp. 365–388, 2020, doi: 10.1137/19M1245505.

[60]

F. Gesztesy and E. Tsekanovskii, “On matrix-valued herglotz functions,” Mathematische Nachrichten, vol. 218, no. 1, pp. 61–138, 2000, doi: https://doi.org/10.1002/1522-2616(200010)218:1<61::AID-MANA61>3.0.CO;2-D.

[61]

P. Deift, G. Dubach, C. Tomei, and T. Trogdon, The Toda lattice and universality for the computation of the eigenvalues of a random matrix. Cambridge University Press, 2025.

[62]

L. Dieci, R. D. Russell, and E. S. Van Vleck, “On the Compuation of Lyapunov Exponents for Continuous Dynamical Systems,” SIAM Journal on Numerical Analysis, vol. 34, no. 1, pp. 402–423, 1997, doi: 10.1137/S0036142993247311.

[63]

L. Dieci and T. Eirola, “On Smooth Decompositions of Matrices,” SIAM Journal on Matrix Analysis and Applications, vol. 20, no. 3, pp. 800–819, 1999, doi: 10.1137/S0895479897330182.

For integers \(m,n\), the matrix \(\mathbf{I}_{m\times n}\) denotes the \(m\times n\) identity matrix, and we abbreviate \(\mathbf{I}_n := \mathbf{I}_{n\times n}\).↩︎
For a Jacobi matrix, all eigenvalues are simple, and the first component of each eigenvector is nonzero.↩︎
The normalized measure \({\boldsymbol{\nu}}\) is given by \({\boldsymbol{\nu}}(x) = {\boldsymbol{\mu}}(\mathbb{R})^{-1/2} {\boldsymbol{\mu}}(x) {\boldsymbol{\mu}}(\mathbb{R})^{-1/2}\). Equivalently, if \({\boldsymbol{\mu}}(\mathbb{R}) = \mathbf{L} \mathbf{L}^*\) is the Cholesky decomposition of the total mass, then \({\boldsymbol{\nu}}\) can be expressed as \({\boldsymbol{\nu}}(x) = \mathbf{L}^{-1} {\boldsymbol{\mu}}(x) \mathbf{L}^{-*}\).↩︎
We use \(\diamond\) as a placeholder for the variable of a function. For instance, \(\frac{1}{\diamond}\) denotes the function \(f(x) = \frac{1}{x}\).↩︎
In this paper, we define the \((n-1)\)-th orthonormal polynomial as the unique \(k\times (k-\ell)\) polynomial \(\mathbf{P}_{n-1}\), where \(\mathrm{rank}~\gamma_{n-1}=k-\ell\), satisfying \(\langle \mathbf{P}_i,\mathbf{P}_j \rangle_{\boldsymbol{\mu}}= \mathbf{0}_{k\times(k-\ell)}\) for \(i\neq j\) and \(\langle \mathbf{P}_i,\mathbf{P}_i \rangle_{\boldsymbol{\mu}}= \mathbf{I}_{k-\ell}\) for \(i=0,\dots,n-1\), with \(\langle \mathbf{P}_{n-1},\diamond\mathbf{P}_{n-2}(\diamond) \rangle_{\boldsymbol{\mu}}\) is in row echelon form with positive pivots.↩︎
Here \(\mathbf{A}^\dagger\) denotes the Moore-Penrose pseudoinverse of \(\mathbf{A}\).↩︎
We use the notation \(\mathbf{x}_{i:j}\) to denote the entries of \(\mathbf{x}\) from index \(i\) through \(j\).↩︎

Banded Hermitian Matrices, Matrix Orthogonal Polynomials, and the Toda Lattice