1 Introduction

In the direct eigenvalue problem for a given square matrix, the objective is to determine its spectrum, which includes its eigenvalues and eigenvectors. In contrast, an inverse eigenvalue problem (IEP) seeks to construct a matrix from a given set of eigenvalues. Solving an IEP is trivial when the matrix is free of structural restrictions. However, when the matrix's structure is constrained, the problem becomes more challenging. The book by Chu and Golub in [1] categorizes IEP using the structure of matrices and specified matrix entries.

The IEPs studied in this paper are part of a broader category that relies on specific matrix structures and entries. Some IEPs can be solved using the traditional linear algebra method. In contrast, others require more advanced techniques, such as those employed in Lax dynamical systems (e.g., see the paper by Chu in [2] and the paper by Chu and Colub in [3]). Previous works by Akaiwa et al. in [4], [5], [6], and [7] have demonstrated how discrete soliton theory (e.g., see the papers by Hirota in [8] and by Tsujimoto in [9]) can solve certain IEPs, particularly those involving band matrices with LR decomposition factored by bidiagonal matrices.

Received March 17th, 2025, Revised April 16th, 2025, Accepted for publication April 17th, 2025 Copyright © 2024 Published by ITB Institut for Research and Community Service, ISSN: 2337-5760, DOI: 10.5614/j.math.fund.sci.2024.56.3.5

In this paper, the matrix structure is constrained on the tridiagonal matrix called the Jacobi matrix or operator. The IEPs for the Jacobi matrix can describe various physical systems (see the book by Chu and Golub in [10]), such as a vibrating beam or rod, an oscillatory mass-spring system, a composite pendulum, and Strum–Liouville problems. The entries in the matrix obtained by solving IEPs correspond to different physical parameters, such as stress, mass, length, etc. One of the papers on the IEPs for the Jacobi matrix was presented by Boor and Golub in [11]. The Jacobi matrix is also closely related to the theory of orthogonal polynomials (see the book by Chihara in [12]).

Let us consider an × tridiagonal matrix , which is factored as = , using unit lower bidiagonal matrix and upper bidiagonal matrix . Let the matrices , and be expressed as

\[L = \begin{pmatrix} 1 & & & & \\ e_1 & 1 & & & \\ & \ddots & \ddots & \\ & & e_{m-1} & 1 \end{pmatrix}, \quad R = \begin{pmatrix} q_1 & 1 & & \\ & q_2 & \ddots & \\ & & \ddots & 1 \\ & & & q_m \end{pmatrix} \tag{1}\] and

\[A = \begin{pmatrix} x_1 & 1 & & & \\ y_1 & x_2 & \ddots & & \\ & \ddots & \ddots & 1 & \\ & & y_{m-1} & x_m \end{pmatrix}. \tag{2}\]

The form of = is called an LR decomposition. The recursion relations of the entries , , , are

\[x_k=q_k+e_{k-1},\quad k=1,2,\ldots,m, \qquad y_k=q_ke_k,\quad k=1,2,\ldots,m-1,\] where let \(e_0=0\).

Let 2−1 = for = 1,2, … , , and 2 = for = 1,2, … , − 1. We have data of the entries in as

\[U = \{u_1, u_2, \dots, u_{2m-1}\} = \{q_1, e_1, q_2, e_2, \dots, q_{m-1}, e_{m-1}, q_m\},\\] and split into two sets

\[U_{\text{specified}} = \{u_1, u_2, \dots, u_{m-1}\}, \quad U_{\text{unknown}} = \{u_m, u_{m+1}, \dots, u_{2m-1}\}.\]

The degrees of freedom of the entries , specified and unknown are 2 − 1, − 1 and . Suppose that all eigenvalues = {1, 2, … , } of are prescribed.

Note here that their order is free. If the entries \(U_{\text{specified}}\) are specified, then the remaining entries \(U_{\text{unknown}}\) should be uniquely determined from \(U_{\text{specified}}\) and \(\Lambda\) if they exist.

This paper aims to design an algorithm to solve this IEP using the techniques of discrete soliton theory. We derive recurrence relations between matrix entries and eigenvalues. The derived relations can compute unknown entries \(U_{\rm unknown}\) from \(U_{\rm specified}\) and \(\Lambda\).

2 Determinant Expression of the Eigenvalue Problem

We briefly show the techniques for expressing eigenproblems as determinants and moments in the discrete soliton theory. The book by Henrici in [13] provides proof of the following results.

Let \(A^{(0)} = A = LR\) and the eigenproblem of \(A^{(0)}\) be

\[A^{(0)}\mathbf{v}_i^{(0)} = \lambda_i \mathbf{v}_i^{(0)}, \quad j = 1, 2, ..., m.\] (3)

Here, the m dimensional vector \(\mathbf{v}_{j}^{(0)}\) is the eigenvector of \(A^{(0)}\) for the eigenvalue \(\lambda_{j}\).

Suppose that all of \(\lambda_j\) is nonzero, and let \(\mathbf{v}^{(1)} = R\mathbf{v}^{(0)}/\lambda_j\) for j = 1, 2, ..., m. We substitute \(A^{(0)}\) in Eq. (3) by \(A^{(0)} = LR\) and use \(\mathbf{v}^{(1)} = R\mathbf{v}^{(0)}/\lambda_j\), we thus obtain the relations

\[R\mathbf{v}_{j}^{(0)} = \lambda_{j}\mathbf{v}^{(1)}, \ L\mathbf{v}_{j}^{(1)} = \mathbf{v}_{j}^{(0)}, \ j = 1, 2, ..., m.\] (4)

Let \(A^{(1)} = RL\). Using Eqs. (4), we obtain the eigenproblem of \(A^{(1)}\) by

\[A^{(1)}\mathbf{v}_{j}^{(1)} = \lambda_{j}\mathbf{v}_{j}^{(1)}, \quad j = 1, 2, ..., m.\]

Both eigenvalues of \(A^{(0)}\) and \(A^{(1)}\) are same \(\Lambda\).

Let the eigenvectors \(\mathbf{v}_{j}^{(n)}\) of \(A^{(n)}\) for n = 0.1 be expressed as

\[\mathbf{v}_{j}^{(n)} = \begin{pmatrix} v_{1,j}^{(n)} \\ v_{2,j}^{(n)} \\ \vdots \\ v_{m,j}^{(n)} \end{pmatrix}, \quad j = 1, 2, \dots, m, \quad n = 0, 1.\] (5)

We substitute Eq. (5) into Eqs. (4) and let \(v_{0,j}^{(n)}=0\), \(v_{m+1,j}^{(n)}=0\) for j=1,2,...,m and n=0,1. Recalling Eq. (1), we then have the relations,

\[q_i v_{i,j}^{(0)} + v_{i+1,j}^{(0)} = \lambda_j v_{i,j}^{(1)}, \quad e_{i-1} v_{i-1,j}^{(1)} + v_{i,j}^{(1)} = v_{i,j}^{(0)}\] (6)

for \(i=1,2,\ldots,m,\ j=1,2,\ldots,m\). Eliminating \(v_{i,j}^{(1)}\) and \(v_{i,j}^{(0)}\) in Eqs. (6), we obtain the recursion relations,

\[v_{i+1,j}^{(0)} = (\lambda_j - q_i - e_{i-1})v_{i,j}^{(0)} - q_{i-1}e_{i-1}v_{i-1,j}^{(0)}, \tag{7}\]

\[v_{i+1,j}^{(1)} = (\lambda_j - q_i - e_i)v_{i,j}^{(1)} - q_i e_{i-1}v_{i-1,j}^{(1)}\] (8)

for \(i=1,2,\ldots,m\), \(j=1,2,\ldots,m\), respectively. Using Eqs. (7), (8) for \(i=1,2,\ldots,m\) repeatedly, it turns out that \(v_{i+1,j}^{(n)}\) is expressed by the polynomial of \(\lambda_j\) of degree i with coefficients that consist of \(q_1,e_1,q_2,e_2,\ldots,q_i,e_i\) and arbitrary constant \(v_{1,j}^{(n)}\).

Let \(\phi_k^{(n)}(z)\) for k = 1, 2, ..., m and n = 0, 1 be monic polynomials of degree k for a variable \(z \in \mathbb{C}\) which satisfy the boundary conditions,

\[\phi_{-1}^{(n)}(z) = 0, \quad \phi_0^{(n)}(z) = 1, \quad \phi_m^{(n)}(z) = \det(zI_m - A^{(n)})\] (9)

for n = 0.1, and recursion relations,

\[z\phi_{k-1}^{(1)}(z) = \phi_k^{(0)}(z) + q_k\phi_{k-1}^{(0)}(z), \quad k = 1, 2, \dots, m,\] (10)

\[\phi_k^{(0)}(z) = \phi_k^{(1)}(z) + e_k \phi_{k-1}^{(1)}(z), \quad k = 0, 1, 2, ..., m.\] (11)

Here, \(I_m\) is an identity matrix of size m. Using techniques in the theory of orthogonal polynomials (see the paper by Spiridonov and Zhedanov in [14]), the sequences of polynomials \(\{\phi_k^{(0)} \mid k=0,1,...\}\) and \(\{\phi_k^{(1)} \mid k=0,1,...\}\) can be proved to be orthogonal under some linear integration, respectively. The relations (10) and (11) transform between \(\{\phi_k^{(0)}\}\) and \(\{\phi_k^{(1)}\}\) each other, and they are called Christoffel and Geronimus transformations, respectively.

Suppose that \(v_{i,j}^{(n)}\) satisfy

\[v_{1,j}^{(n)} \neq 0, \quad j = 1, 2, ..., m,\] (12)

\[v_{i,j}^{(n)} = v_{1,j}^{(n)} \phi_{i-1}^{(n)} (\lambda_j), \quad i = 0,1,...,m, m+1, \quad j = 1,2,...,m\] (13)

for n = 0,1. We substitute Eq. (13) at i = 0,1 and m + 1 into \(v_{0,j}^{(n)} = 0\), \(v_{1,j}^{(n)} = v_{1,j}^{(n)}\) and \(v_{m+1,j}^{(n)} = 0\). We then have Eqs. (9), because \(\phi_m^{(n)}(\lambda_j) = \det(\lambda_j I_m - A^{(n)}) = 0\) for j = 1,2,...,m. We substitute Eq. (13) into Eqs. (6), thus we obtain Eqs. (10), (11). Therefore, the assumption (13) is proved. If \(v_{1,j}^{(n)} = 0\) then \(\mathbf{v}_j^{(n)} = \mathbf{0}\). We thus need Eq. (12).

Let the characteristic polynomial \(p(z) = \det(zI_m - A) = \prod_{j=1}^m (z - \lambda_j)\) of A be expanded by

\[p(z) = z^m + a_1 z^{m-1} + a_2 z^{m-2} + \dots + a_{m-1} z + a_m.\] (14)

Using the coefficients, \(\{a_1, a_2, ..., a_m\}\), of p(z) in Eq. (14), we introduce the moments \(F = \{f_0, f_1, ..., f_{2m+1}\}\) which satisfy the linear equation,

\[f_i + a_1 f_{i-1} + a_2 f_{i-2} + \dots + a_{m-1} f_{i-m-1} + a_m f_{i-m} = 0\] (15)

for \(i=m,m+1,\ldots,2m+1\) with some initial values \(F_{\rm init}=\{f_0,f_1,\ldots,f_{m-1}\}\). Suppose that \(f_0=1\). It follows from \(\Lambda,F_{\rm init}\) and \(f_0=1\) that the degree of freedom of the moments F is 2m-1.

We introduce the Hankel determinants \(\tau_k^{(n)}\) of size k for the moments \(\{f_n\}\) defined by \(\tau_{-1}^{(n)}=0\), \(\tau_0^{(n)}=1\) and

\[\tau_{k}^{(n)} = \begin{vmatrix} f_{n} & f_{n+1} & \cdots & f_{n+k-1} \\ f_{n+1} & f_{n+2} & \cdots & f_{n+k} \\ \vdots & \vdots & \ddots & \vdots \\ f_{n+k-1} & f_{n+k} & \cdots & f_{n+2k-2} \end{vmatrix}, \quad k = 1, 2, \dots, m, m + 1\] (16)

for n = 0,1. At k = m + 1 in Eq. (16), it follows from Eq. (15) that \(\tau_{m+1}^{(n)} = 0\) for n = 0,1.

We introduce the Hadamard polynomials \(H_k^{(n)}(z)\) of degree k for the moments \(\{f_n\}\) and a variable \(z \in \mathbb{C}\) defined by \(H_{-1}^{(n)}(z) = 0\), \(H_0^{(n)}(z) = 1\) and

\[H_k^{(n)}(z) = \frac{T_k^{(n)}(z)}{\tau_k^{(n)}}, \quad k = 1, 2, ..., m, \quad n = 0, 1,\] (17)

where

\[T_k^{(n)}(z) = \begin{vmatrix} f_n & f_{n+1} & \cdots & f_{n+k-1} & f_{n+k} \\ f_{n+1} & f_{n+2} & \cdots & f_{n+k} & f_{n+k+1} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ f_{n+k-1} & f_{n+k} & \cdots & f_{n+2k-2} & f_{n+2k-1} \\ 1 & z & \cdots & z^{k-1} & z^k \end{vmatrix}.\](18)

Note here that both \(H_k^{(n)}(z)\) and \(T_k^{(n)}(z)\) are polynomials of degree k, and \(H_k^{(n)}(z)\) is monic but \(T_k^{(n)}(z)\) is not. At k=m in Eqs. (17), (18), it follows from Eqs. (14), (15) that

\[H_m^{(n)}(\lambda_j) = 0, \quad j = 1, 2, ..., m\] (19)

for n = 0,1.

The determinants \(\tau_k^{(n)}\) in Eq. (16) and \(T_k^{(n)}(z)\) in Eq. (18) satisfy some identities such as

\[z\tau_k^{(n)}T_{k-1}^{(n+1)}(z) = \tau_{k-1}^{(n+1)}T_k^{(n)}(z) + \tau_k^{(n+1)}T_{k-1}^{(n)}(z), \tag{20}\]

\[\tau_k^{(n+1)} T_k^{(n)}(z) = \tau_k^{(n)} T_k^{(n+1)}(z) + \tau_{k+1}^{(n)} T_{k-1}^{(n+1)}(z), \tag{21}\] which are called Jacobi's identity and Plücker relation, respectively. The proofs of the determinants identities (20), (21) are written in [13] precisely. Suppose that the Hankel determinants \(\tau_k^{(n)}\) satisfy

\[\tau_k^{(n)} \neq 0, \quad k = 1, 2, ..., m, \quad n = 0, 1.\] (22)

Dividing Eqs. (20), (21) by \(\tau_k^{(n)} \tau_{k-1}^{(n+1)}\) or \(\tau_k^{(n)} \tau_k^{(n+1)}\) and using Eq. (17), we obtain the relations of \(H_k^{(n)}(z)\) by

\[zH_{k-1}^{(n+1)}(z) = H_k^{(n)}(z) + \frac{\tau_{k-1}^{(n)}\tau_k^{(n+1)}}{\tau_k^{(n)}\tau_{k-1}^{(n+1)}}H_{k-1}^{(n)}(z), \tag{23}\]

\[H_k^{(n)}(z) = H_k^{(n+1)}(z) + \frac{\tau_{k+1}^{(n)} \tau_{k-1}^{(n+1)}}{\tau_k^{(n)} \tau_k^{(n+1)}} H_{k-1}^{(n+1)}(z). \tag{24}\]

Comparing Eqs. (9), (10), (11) with Eqs. (19), (23), (24), we thus obtain the determinant expression of \(q_k\), \(e_k\) and \(\phi_k^{(n)}(z)\) by

\[q_k = \frac{\tau_{k-1}^{(0)} \tau_k^{(1)}}{\tau_k^{(0)} \tau_{k-1}^{(1)}}, \quad k = 1, 2, \dots, m,\] (25)

\[e_k = \frac{\tau_{k+1}^{(0)} \tau_{k-1}^{(1)}}{\tau_k^{(0)} \tau_k^{(0)}}, \quad k = 0, 1, \dots, m,\] (26)

and

\[\phi_k^{(n)}(z) = H_k^{(n)}(z), \quad k = -1, 0, 1, \dots, m - 1, m\] (27)

for n = 0,1.

The mentioned relations are the one-step mapping from the eigenvector \(\mathbf{v}^{(0)}\) of \(A^{(0)}\) to \(\mathbf{v}^{(1)}\) of \(A^{(1)}\). We employ them to connect \(\Lambda\), \(U_{\text{specified}}\) and \(U_{\text{unknown}}\) to derive the algorithm to solve the IEP. In this paper, we do not use the mappings,

\[A^{(n)} = L^{(n)}R^{(n)} \to A^{(n+1)} = R^{(n)}L^{(n)} = L^{(n+1)}R^{(n+1)}, \quad n = 0,1,...\]

explicitly. The matrices \(L^{(n)}\), \(R^{(n)}\) and \(A^{(n)}\) are same as Eqs. (1), (2) with replacing the entries \(x_k\), \(y_k\), \(q_k\), \(e_k\) by \(x_k^{(n)}\), \(y_k^{(n)}\), \(q_k^{(n)}\), \(e_k^{(n)}\). The mapping \(L^{(n+1)}R^{(n+1)}=R^{(n)}L^{(n)}\) is called an LR transformation, and the relations of the entries are written by

\[q_k^{(n+1)} + e_{k-1}^{(n+1)} = q_k^{(n)} + e_k^{(n)}, \quad k = 1, 2, ..., m,\] (28)

\[q_k^{(n+1)}e_k^{(n+1)} = q_{k+1}^{(n)}e_k^{(n)}, \quad k = 1, 2, ..., m-1.\] (29)

The mapping (28)-(29) is called the qd algorithm in numerical analysis and the discrete Toda equation in the soliton theory. Rutishauser presented the origin of the qd algorithm in [15]. The origin of discrete soliton equations was analyzed by Hirota (e.g., see the paper by Hirota et al. in [16]).

3 Algorithm to Construct the Matrix

We derive an algorithm to construct A with the prescribed \(\Lambda\) and the specified entries \(U_{\text{specified}}\).

Let \(\psi_{2k-1} = \phi_k^{(0)}\) for k = 1, 2, ..., m and \(\psi_{2k} = \phi_k^{(1)}\) for k = 1, 2, ..., m - 1. We introduce the sequence of the polynomials \(\psi_i\) defined by

\[\begin{split} \Psi &= \{\psi_1, \psi_2, \psi_3, \psi_4, \dots, \psi_{2m-3}, \psi_{2m-2}, \psi_{2m-1}\} \\ &= \{\phi_1^{(0)}, \phi_1^{(1)}, \phi_2^{(0)}, \phi_2^{(1)}, \dots, \phi_{m-1}^{(0)}, \phi_{m-1}^{(1)}, \phi_m^{(0)}\}. \end{split}\]

Let m-1 of polynomials of \(\Psi\) be \(\Psi_{\text{specified}} = \{\psi_1, \psi_2, ..., \psi_{m-1}\}\) and m of the remaining polynomials of \(\Psi\) be \(\Psi_{\text{unknown}} = \{\psi_m, \psi_{m+1}, ..., \psi_{2m-1}\}\). From Eqs. (10), (11), we obtain the following Algorithm 1.

Algorithm 1. According to the following procedure and using \(U_{\text{specified}}\), we can compute \(\Psi_{\text{specified}}\).

• 1. Let \(\phi_0^{(0)}(z) = 1\), \(\phi_0^{(1)}(z) = 1\).
• 2. foreach i = 1, 2, ..., m 1 do
• 3. if i is odd then
• 4. Compute \(\psi_i(z) = \phi_k^{(0)}(z) = z\phi_{k-1}^{(1)}(z) q_k\phi_{k-1}^{(0)}(z)\) with \(k = \frac{i+1}{2}\).
• 5. else
• 6. Compute \(\psi_i(z) = \phi_k^{(1)}(z) = \phi_k^{(0)}(z) e_k \phi_{k-1}^{(1)}(z)\) with \(k = \frac{i}{2}\).
• 7. endif
• 8. enddo

We expand the polynomials \(\phi_k^{(n)}(z)\) and express them by

\[\phi_k^{(n)}(z) = z^k + \sum_{j=1}^k b_{k,j}^{(n)} z^{k-j}, \quad k = 1,2,...,m\] for n = 0,1. From Eqs. (17), (18) and (27), we expand the determinants in Eq. (18) by using cofactor expansion along the last row. We then obtain the coefficients \(b_{k,j}^{(n)}\) of \(\phi_k^{(n)}\) by

\[b_{k,j}^{(n)} = \frac{\Delta_{k-j}^{(n)}}{\tau_k^{(n)}}, \quad j = 1, 2, \dots, k, \quad k = 1, 2, \dots, m\] (30)

for n = 0,1, where \(\Delta_{k-j}^{(n)}\) is the (k+1,k-j+1) cofactor of the determinant in Eq. (18). The coefficient (30) is the same as the solution of the linear equation

\[\begin{pmatrix} f_{n} & f_{n+1} & \cdots & f_{n+k-1} & f_{n+k} \\ f_{n+1} & f_{n+2} & \cdots & f_{n+k} & f_{n+k+1} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ f_{n+k-1} & f_{n+k} & \cdots & f_{n+2k-2} & f_{n+2k-1} \\ 0 & 0 & \cdots & 0 & 1 \end{pmatrix} \begin{pmatrix} b_{k,k}^{(n)} \\ b_{k,k-1}^{(n)} \\ \vdots \\ b_{k,1}^{(n)} \\ 1 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ 1 \end{pmatrix}\](31)

by using the Cramer formula and Eq. (16). The linear equation (31) is rewritten as

\[f_{n+k+i} + \sum_{j=1}^{k} b_{k,j}^{(n)} f_{n+k+i-j} = 0, \quad i = 0,1,2,...,k-1\] (32)

for n = 0,1, k = 1,2,...,m. At (n,i) = (0,k-1),(1,k-1) in Eq. (32), we obtain the relations

\[f_{2k-1} = -\sum_{i=1}^{k} b_{k,i}^{(0)} f_{2k-1-i}, \quad k = 1, 2, ..., m,\] (33)

\[f_{2k} = -\sum_{j=1}^{k} b_{k,j}^{(1)} f_{2k-j}, \quad k = 1, 2, ..., m.\] (34)

Let m-1 of the moments of F be \(F_{\text{specified}}=\{f_1,\ldots,f_{m-1}\}\). From Eqs. (33), (34), we obtain the following Algorithm 2.

Algorithm 2. According to the following procedure and using the coefficients \(b_{k,i}^{(n)}\) of \(\phi_k^{(n)}\) in \(\Psi_{\text{specified}}\), we can compute \(F_{\text{specified}}\).

• 1. Let \(f_0 = 1\).
• 2. foreach i = 1, 2, ..., m 1 do
• 3. Compute \(f_i = -\sum_{j=1}^k b_{k,j}^{(0)} f_{i-j}\) with \(k = \frac{i+1}{2}\) if i is odd.
• 4. Compute \(f_i = -\sum_{j=1}^k b_{k,j}^{(1)} f_{i-j}\) with \(k = \frac{i}{2}\) if i is even.
• 5. enddo

Let m of the moments of F be \(F_{\text{unknown}} = \{f_m, f_{m+1}, \dots, f_{2m-1}\}\). From Eq. (15), we obtain the following Algorithm 3.

Algorithm 3. According to the following procedure and using \(f_0 = 1\), \(F_{\text{specified}}\) and the coefficients \(\{a_1, a_2, ..., a_m\}\) in Eq. (14) form \(\Lambda\), we can compute \(F_{\text{unknown}}\).

• 1. foreach i = m, m + 1, ..., 2m 1 do
• 2. Compute \(f_i = -\sum_{j=1}^m a_j f_{i-j}\).
• 3. enddo

We recall that \(\tau_0^{(0)} = \tau_0^{(1)} = 1\), \(\tau_1^{(0)} = f_0 = 1\), and \(\tau_{m+1}^{(0)} = \tau_{m+1}^{(1)} = 0\). Let \(\sigma_{2k-2} = \tau_k^{(0)}\) and \(\sigma_{2k-1} = \tau_k^{(1)}\) for k = 0, 1, ..., m+1. We introduce the sequence of the Hankel determinants \(\tau_k^{(n)}\) defined by

\[\begin{split} \Sigma &= \{\sigma_{-2}, \sigma_{-1}, \sigma_0, \sigma_1, \sigma_2, \sigma_3, \dots, \sigma_{2m-2}, \sigma_{2m-1}, \sigma_{2m}, \sigma_{2m+1}\} \\ &= \{1, 1, 1, \tau_1^{(1)}, \tau_2^{(0)}, \tau_2^{(1)}, \dots, \tau_m^{(0)}, \tau_m^{(1)}, 0, 0\}. \end{split}\]

The Hankel determinant \(\sigma_i\) can be computed by Eq. (16) from \(\{f_0, f_1, ..., f_i\}\) for i=1,2,...,2m+1. We substitute \(\tau_k^{(0)}=\sigma_{2k-2}\) and \(\tau_k^{(1)}=\sigma_{2k-1}\) in Eqs. (25), (26). From \(u_{2k-1}=q_k, u_{2k}=e_k\), we thus have relations

\[u_i = u_{2k-1} = q_k = \frac{\sigma_i \sigma_{i-3}}{\sigma_{i-1} \sigma_{i-2}}, \quad i = 2k-1, \quad k = 1, 2, ..., m,\] (35)

\[u_i = u_{2k} = e_k = \frac{\sigma_i \sigma_{i-3}}{\sigma_{i-1} \sigma_{i-2}}, \quad i = 2k, \quad k = 1, 2, \dots, m-1.\] (36)

The members \(\sigma_i\) of \(\Sigma\) to express Eqs. (35), (36) for all members of \(U_{\text{specified}}\) are \(\Sigma_{\text{specified}} = \{\sigma_{-2}, \sigma_{-1}, \sigma_0, \sigma_1, \dots, \sigma_{m-1}\}\). Similarly, the members for \(U_{\text{unknown}}\) are \(\Sigma_{\text{unknown}} = \{\sigma_{m-3}, \sigma_{m-2}, \dots, \sigma_{2m-1}\}\). From Eqs. (16), (35)-(36), we obtain the following Algorithm 4.

Algorithm 4. According to the following procedure and using \(f_0 = 1\), \(F_{\text{specified}}\) and \(F_{\text{unknown}}\), we can compute \(U_{\text{unknown}}\) if the procedure does not abort.

• 1. foreach i = m 3, m 2, ..., 2m 1 do
• 2. Compute \(\sigma_i = \tau_k^{(1)}\) with \(k = \frac{i+1}{2}\) by Eq. (16) if i is odd.
• 3. Compute \(\sigma_i = \tau_k^{(0)}\) with \(k = \frac{i+2}{2}\) by Eq. (16) if i is even.

• 4. if = 0 then this procedure abort.
• 5. enddo
• 6. foreach = , + 1, … ,2 − 1 do
• 7. Compute = with = +1 2 by Eq. (35) if is odd.
• 8. Compute = with = 2 by Eq. (36) if is even.
• 9. enddo

Algorithms 1-3 can finish successfully under any , specified. However, Algorithm 4 has two cases in which it either finishes successfully or aborts because of the division by 0. In Algorithm 4, we need the conditions (22) for expressing all of , as Eqs. (25), (26). Therefore, we need the condition of all ≠ 0 and ≠ 0.

As a result, we proved the following Theorem 1.

Theorem 1. If Algorithm 4 finishes successfully, then Algorithms 1-4 can compute = under the condition of all ≠ 0, ≠ 0 from the prescribed eigenvalues and the specified entries specified. If Algorithm 4 aborts, then there does not exist = under the condition of all ≠ 0, ≠ 0.

4 Examples

We demonstrate the proposed algorithm's efficiency by showing some examples, including cases with real, complex, and multiple eigenvalues. We also include a case for the algorithm that does not solve the IEP.

We consider 4 × 4 matrix = be expressed by

\[L = \begin{pmatrix} 1 & 0 & 0 & 0 \\ e_1 & 1 & 0 & 0 \\ 0 & e_2 & 1 & 0 \\ 0 & 0 & e_3 & 1 \end{pmatrix}, \quad R = \begin{pmatrix} q_1 & 1 & 0 & 0 \\ 0 & q_2 & 1 & 0 \\ 0 & 0 & q_3 & 1 \\ 0 & 0 & 0 & q_4 \end{pmatrix},\] \[A = LR = \begin{pmatrix} q_1 & 1 & 0 & 0 \\ q_1e_1 & q_2 + e_1 & 1 & 0 \\ 0 & q_2e_2 & q_3 + e_2 & 1 \\ 0 & 0 & q_3e_3 & q_4 + e_3 \end{pmatrix}.\]

Using the proposed Algorithms 1-4, we solve the IEP with the prescribed eigenvalues \(\Lambda = \{\lambda_1, \lambda_2, \lambda_3, \lambda_4\}\) of A and the specified entries \(U_{\text{specified}} = \{u_1, u_2, u_3\} = \{q_1, e_1, q_2\}\) to obtain unknown entries

\[U_{\text{unknown}} = \{u_4, u_5, u_6, u_7\} = \{e_2, q_3, e_3, q_4\}.\]

Let m = 4 and

\[\begin{split} \psi_1 &= \phi_1^{(0)} = z + b_{1,1}^{(0)}, \quad \psi_2 = \phi_1^{(1)} = z + b_{1,1}^{(1)}, \\ \psi_3 &= \phi_2^{(0)} = z^2 + b_{2,1}^{(0)} z + b_{2,2}^{(0)}, \\ F_{\text{specified}} &= \{f_1, f_2, f_3\}, \quad F_{\text{unknown}} = \{f_4, f_5, f_6, f_7\}, \\ p(z) &= (z - \lambda_1) \cdots (z - \lambda_4) = z^4 + a_1 z^3 + a_2 z^2 + a_3 z + a_4, \\ \Sigma_{\text{unknown}} &= \{\sigma_1, \sigma_2, \dots, \sigma_7\} = \{\tau_1^{(1)}, \tau_2^{(0)}, \tau_2^{(1)}, \tau_3^{(0)}, \tau_3^{(1)}, \tau_4^{(0)}, \tau_4^{(1)}\}. \end{split}\]

The results of examples are as follows.

4.1 Case 1: Real and Single Eigenvalues

Inputs: \[\Lambda = \{1,2,3,4\}, U_{\text{specified}} = \{2,1,3\}.\]

Algorithm 1: \[\psi_1 = z - 2\], \(\psi_2 = z - 3\) and \(\psi_3 = z^2 - 6z + 6\).

Algorithm 2: \[f_0 = 1\], \(F_{\text{specified}} = \{2,6,24\}\).

Algorithm 3:

\[p(z) = z^4 - 10z^3 + 35z^2 - 50z + 24\],
\(F_{\text{unknown}} = \{106,472,2066,8864\}\).

Algorithm 4:

\[\Sigma_{\text{unknown}} = \{2,2,12,-4,-104,-32,-768\},\\] \[U_{\text{unknown}} = \left\{-\frac{1}{3},\frac{13}{3},-\frac{12}{13},\frac{12}{13}\right\}.\]

Output: Successfully, the tridiagonal matrix A = LR with \(\Lambda\) and \(U_{\text{specified}}\) is computed by

\[A = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 0 & -\frac{1}{3} & 1 & 0 \\ 0 & 0 & -\frac{12}{13} & 1 \end{pmatrix} \begin{pmatrix} 2 & 1 & 0 & 0 \\ 0 & 3 & 1 & 0 \\ 0 & 0 & \frac{13}{3} & 1 \\ 0 & 0 & 0 & \frac{12}{13} \end{pmatrix} = \begin{pmatrix} 2 & 1 & 0 & 0 \\ 2 & 4 & 1 & 0 \\ 0 & -1 & 4 & 1 \\ 0 & 0 & -4 & 0 \end{pmatrix}.\]

4.2 Case 2: Complex and Multiple Eigenvalues

Inputs: \(\Lambda = \{1 + i, 1 + i, 1 - i, 1 - i\}, U_{\text{specified}} = \{2,1,3\}.\)

Algorithm 1: \(\psi_1 = z - 2\), \(\psi_2 = z - 3\) and \(\psi_3 = z^2 - 6z + 6\).

Algorithm 2: \(f_0 = 1\), \(F_{\text{specified}} = \{2,6,24\}\).

Algorithm 3:

\[p(z) = z^4 - 4z^3 + 8z^2 - 8z + 4\], \(F_{\text{unknown}} = \{60,88,40,-160\}\).

Algorithm 4:

\[\Sigma_{\text{unknown}} = \{2, 2, 12, -96, -2688, 25600, 102400\},\] \[U_{\text{unknown}} = \left\{-8, \frac{14}{3}, \frac{25}{21}, \frac{1}{7}\right\}.\]

Output: Successfully, the tridiagonal matrix A = LR with \(\Lambda\) and \(U_{\text{specified}}\) is computed by

\[A = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 0 & -8 & 1 & 0 \\ 0 & 0 & \frac{25}{21} & 1 \end{pmatrix} \begin{pmatrix} 2 & 1 & 0 & 0 \\ 0 & 3 & 1 & 0 \\ 0 & 0 & \frac{14}{3} & 1 \\ 0 & 0 & 0 & \frac{1}{7} \end{pmatrix} = \begin{pmatrix} 2 & 1 & 0 & 0 \\ 2 & 4 & 1 & 0 \\ 0 & -24 & -\frac{10}{3} & 1 \\ 0 & 0 & \frac{50}{9} & \frac{4}{3} \end{pmatrix}.\]

4.3 Case 3: Nonexistent

Inputs: \(\Lambda = \{1,2,3,4\}, U_{\text{specified}} = \{1,-1,2\}.\)

Algorithm 1: \(\psi_1 = z - 1\), \(\psi_2 = z\) and \(\psi_3 = z^2 - 2z + 2\).

Algorithm 2: \(f_0 = 1\), \(F_{\text{specified}} = \{1,0,-2\}\).

Algorithm 3: \(p(z) = z^4 - 10z^3 + 35z^2 - 50z + 24\),

\[F_{\text{unknown}} = \{6,106,750,4138\}.\]

Algorithm 4:

\[\Sigma_{\text{unknown}} = \{1, -1, -2, -10, -240, 0, 0\}, \quad U_{\text{unknown}} = \{-5, 12, 0, *\}.\]

The symbol '*' indicates failure to compute ⁴ = 7⁴ 6⁵ = 0⋅(−10) 0⋅(−240) .

Output: The tridiagonal matrix = with and specified does not exist under the condition of all ≠ 0 and ≠ 0.

5 Conclusion

This paper presents an algorithm for constructing a tridiagonal matrix factored by bidiagonal matrices based on a given set of eigenvalues and specified matrix entries. We derive recurrence relations connecting the eigenvalues to the unknown matrix entries using the discrete soliton theory. The algorithm was tested on examples, including real, complex, and multiple eigenvalue cases, demonstrating its efficiency and applicability to different IEP scenarios. Future work could explore extending this approach to other classes of band matrices and further generalizing the problem.

Acknowledgement

Japan Society for the Promotion of Science, Grants-in-Aid for Scientific Research, Grant Number JP23K03223 supported this work.