Krylov and core transformation algorithms for an inverse eigenvalue problem to compute recurrences of multiple orthogonal polynomials

This paper develops and analyzes two numerical algorithms based on inverse eigenvalue problems and linear algebra techniques to compute recurrence coefficients for multiple orthogonal polynomials on the step-line, demonstrating their accuracy and stability through experiments on both ill-conditioned and well-conditioned examples.

The Gist

Here is an explanation of the paper, translated into everyday language with creative analogies.

The Big Picture: Reversing the Recipe

Imagine you are a master chef. You have a secret recipe (a set of rules) that tells you how to bake a perfect cake. If you follow the recipe, you get a specific cake.

In the world of mathematics, Orthogonal Polynomials are like those cakes. They are special mathematical functions used to solve complex problems in physics, engineering, and statistics. Usually, mathematicians know the "ingredients" (the rules or weights) and use them to bake the cake (compute the polynomial).

This paper is about doing the opposite.

Imagine someone hands you a finished cake and says, "I don't know the recipe, but I know exactly what this cake tastes like and how it looks. Can you figure out the secret recipe that created it?"

This is called an Inverse Eigenvalue Problem. The "cake" is the data (specific points and weights), and the "recipe" is a special grid of numbers (a matrix) that generates the polynomials. The authors of this paper developed two new, high-tech kitchen tools to reverse-engineer this recipe.

---

The Cast of Characters

1. Multiple Orthogonal Polynomials (MOPs):
   Think of standard polynomials as a solo singer hitting a perfect note. MOPs are like a choir. They have to sing in harmony with multiple different conductors (measures) at the same time. It's much harder to keep everyone in tune with multiple conductors than with just one.

2. The "Step-Line" Recurrence:
   To keep the choir in sync, you need a conductor's baton that follows a specific pattern. In math, this pattern is called a "recurrence relation." It's a set of instructions like, "To get the next note, take the previous note, add a little bit of this, and subtract a bit of that." The paper focuses on a specific, zig-zagging pattern called the "step-line."

3. The Goal:
   The authors want to find the exact numbers (the coefficients) for that baton's instructions, given only the final result (the nodes and weights).

---

The Two New Tools (Algorithms)

The paper proposes two different ways to solve this "reverse recipe" puzzle.

Tool 1: The "Krylov" Method (The Echo Chamber)

• The Analogy: Imagine you are in a large, empty hall (a Krylov subspace). You shout a sound (a starting vector), and it bounces off the walls. You listen to the echo, shout again, and listen to the new echo.
• How it works: This method uses a process called Biorthogonal Lanczos. It's like having two choirs singing back and forth. One choir sings a note, and the other choir listens and adjusts its pitch to be perfectly "orthogonal" (at a right angle, or perfectly distinct) to the first.
• The Magic: By carefully listening to these echoes and adjusting the pitch step-by-step, the algorithm builds the "baton instructions" (the matrix) piece by piece.
• The Catch: In a noisy room (computer with limited precision), echoes can get messy. The authors found that if they didn't clean up the noise (re-orthogonalize), the choir would eventually lose its tune, and the recipe would be wrong.

Tool 2: The "Core Transformation" Method (The Sculptor)

• The Analogy: Imagine you have a block of marble (a simple diagonal matrix of numbers) and you want to carve it into a specific, intricate statue (the banded Hessenberg matrix).
• How it works: This method uses Gaussian Elimination as a chisel. It starts with a simple, flat block of data. Then, it applies a series of tiny, precise cuts (transformations) to move numbers around.
• The "Chasing" Step: As the sculptor cuts away the excess stone, a "bulge" (an unwanted number) might pop up in the wrong place. The algorithm has a special move called "bulge chasing" to push that unwanted number down the line until it disappears off the edge of the statue.
• The Result: After a sequence of cuts and chases, the simple block of marble is transformed into the complex, structured statue (the recurrence matrix) needed to generate the polynomials.

---

The Experiment: Did the Tools Work?

The authors tested these tools on two famous types of "cakes" (mathematical problems based on Kravchuk and Hahn polynomials).

• The Problem: These specific cakes are notoriously difficult. They are like trying to balance a house of cards in a hurricane. The math is ill-conditioned, meaning a tiny mistake in the ingredients leads to a massive error in the final cake.
• The Results:
  • The "Echo Chamber" (Krylov) without cleaning: Failed miserably. The noise overwhelmed the signal, and the recipe was garbage.
  • The "Echo Chamber" with cleaning (Reorthogonalization): Did a great job! It kept the choir in tune even in the hurricane.
  • The "Sculptor" (Core Transformation): Also did a great job. It was very stable and accurate.

The Verdict: Both the "cleaned-up Echo Chamber" and the "Sculptor" are excellent tools. However, the Sculptor is a bit more robust for the messy, difficult problems, while the Echo Chamber (with cleaning) is very accurate but requires more computing power to keep the noise down.

Why Should You Care?

You might not be baking mathematical cakes, but these tools are used everywhere:

• Random Matrix Theory: Understanding how complex systems (like the stock market or quantum particles) behave.
• Signal Processing: Cleaning up noisy data.
• Approximation: Predicting future trends based on past data.

By figuring out how to reverse-engineer these complex mathematical recipes, the authors have given scientists better tools to understand and predict the behavior of complex systems, even when the data is messy and the math is tricky. They turned a chaotic puzzle into a solvable recipe.

Technical Summary

Here is a detailed technical summary of the paper "Krylov and Core Transformation Algorithms for an Inverse Eigenvalue Problem to Compute Recurrences of Multiple Orthogonal Polynomials."

1. Problem Statement

The paper addresses the computation of recurrence coefficients for Multiple Orthogonal Polynomials (MOPs) defined on a "step-line" with respect to a system of discrete positive measures.

• Context: Classical orthogonal polynomials ($r=1$) satisfy a three-term recurrence relation linked to a symmetric tridiagonal (Jacobi) matrix. MOPs ($r > 1$) satisfy a more complex $(r+2)$-term recurrence relation, represented by a banded upper Hessenberg matrix ($H_N$).
• The Inverse Eigenvalue Problem (IEP): The authors formulate the problem as an IEP. Given:
  1. A set of distinct nodes $\{z_i\}_{i=1}^N$ (eigenvalues).
  2. Weight vectors $\{\alpha_{j,i}\}$ associated with $r$ discrete measures.
  • Goal: Construct the banded Hessenberg matrix $H_N$ containing the recurrence coefficients such that its eigenvalues are the nodes $z_i$, and its eigenvectors relate to the weights and the orthogonal polynomials.
• Specific Constraints: The problem requires finding a matrix $H_N$ and biorthogonal bases $W$ and $V$ such that:
  • $W^T V = I_N$ (Biorthogonality).
  • $W^T Z V = H_N$ (where $Z$ is the diagonal matrix of nodes).
  • $H_N$ has a specific banded structure with ones on the first subdiagonal (corresponding to monic Type II polynomials).

2. Methodology

The authors propose and analyze two distinct numerical approaches to solve this IEP, both leveraging the link between MOPs and Krylov subspaces.

A. Biorthogonal Lanczos Approach (Krylov-based)

This method exploits the connection between the recurrence relations of MOPs and block Krylov subspaces.

• Mechanism: It uses a specialized biorthogonal Lanczos algorithm with multiple starting vectors ($w_1, w_2$ for Type I and $v_1$ for Type II).
• Algorithm (IEP KRYL):
  • Generates nested biorthogonal bases for the Krylov subspaces $K_m(Z, v_1)$ and $K^\square_m(Z, [w_1, w_2])$.
  • Computes the recurrence coefficients iteratively via short recurrences.
• Numerical Stability Enhancements:
  • Reorthogonalization: To combat the loss of biorthogonality in finite precision, the authors implement full reorthogonalization (Algorithm 4.2: IEP KRYLREORTH(full)). This increases complexity from $O(N^2)$ to $O(N^3)$ but significantly improves accuracy.
  • Normalization: The algorithm handles the scaling of basis vectors to ensure the subdiagonal entries of $H_N$ are exactly 1, a requirement for monic polynomials.

B. Core Transformation Approach

This method generalizes techniques used for standard orthogonal polynomials to the MOP setting using Gaussian eliminations (core transformations).

• Mechanism: It applies a sequence of non-unitary similarity transformations to the diagonal matrix $Z$ to transform it into the desired banded Hessenberg form $H_N$.
• Algorithm (IEP CORE):
  • Elimination: Uses lower/upper triangular eliminators (2x2 active blocks) to zero out specific entries in the weight and node vectors ($w_1, w_2, v_1$).
  • Biorthogonalization: Enforces the condition $W^T V = I$ via LU factorizations during the process.
  • Chasing: Performs "bulge chasing" steps to move and eliminate unwanted non-zero elements outside the prescribed bandwidth, restoring the banded structure.
  • Scaling: A final diagonal scaling ensures the subdiagonal entries are unity.

3. Key Contributions

1. Formulation of MOP Recurrence as IEP: The paper rigorously establishes the link between the step-line recurrence of MOPs and an inverse eigenvalue problem involving a banded Hessenberg matrix and biorthogonal eigenvectors.
2. Algorithm Development:
   • Adaptation of the block biorthogonal Lanczos method for MOPs.
   • Development of a core transformation algorithm specifically tailored for the multi-measure case ($r=2$ in the experiments, but generalizable).
3. Numerical Analysis:
   • Detailed analysis of breakdowns (near-zero pivots) and strategies to handle them (e.g., look-ahead, though found less critical than conditioning issues).
   • Investigation of conditioning: The authors demonstrate that the problem is inherently ill-conditioned for classical MOP families (Kravchuk, Hahn) but well-conditioned for random inputs.
4. Comparative Evaluation: A comprehensive benchmarking of the algorithms against analytical solutions and high-precision quadruple-precision baselines.

4. Results and Experiments

The authors tested the algorithms on two types of problems:

• Ill-Conditioned Cases: Multiple Kravchuk and Hahn polynomials (where analytical solutions exist).
• Well-Conditioned Cases: Random weights and Chebyshev nodes.

Key Findings:

• Accuracy:
  • For ill-conditioned problems (Kravchuk/Hahn), both IEP CORE and IEP KRYLREORTH(full) perform similarly and significantly better than the standard IEP KRYL (without reorthogonalization).
  • For well-conditioned problems, IEP KRYLREORTH(full) achieves the highest accuracy, maintaining relative errors around $10^{-17}$for moderate$N$, whereas IEP CORE showed higher sensitivity to random weight variations.
• Stability:
  • Backward Stability: Both CORE and KRYLREORTH(full) are backward stable for $N < 25$. Instability appears for larger $N$ in the ill-conditioned cases due to the problem's inherent conditioning, not the algorithm itself.
  • Biorthogonality: The standard Lanczos method (IEP KRYL) suffers from exponential loss of biorthogonality ($W^T V \neq I$) as $N$ increases. Reorthogonalization and Core methods maintain biorthogonality much longer.
• Complexity:
  • IEP KRYL: $O(N^2)$.
  • IEP KRYLREORTH(full): $O(N^3)$ (due to full reorthogonalization).
  • IEP CORE: $O(N^2)$ (similar to standard core transformations).

5. Significance and Conclusion

• Theoretical Insight: The paper successfully bridges the gap between the theory of Multiple Orthogonal Polynomials and modern numerical linear algebra techniques (Krylov subspaces and core transformations).
• Practical Utility: The algorithms provide a robust way to compute recurrence coefficients for MOPs, which are essential for applications in rational approximation, random matrix theory, and spectral theory of non-symmetric operators.
• Recommendation:
  • For ill-conditioned problems (common in classical MOP families), IEP CORE and IEP KRYLREORTH(full) are the preferred methods.
  • For well-conditioned problems where maximum accuracy is required, IEP KRYLREORTH(full) is superior despite the higher computational cost.
• Future Work: The authors suggest extending these methods to other recurrence relations (e.g., nearest-neighbor), exploring vector continued fractions, and developing updating/downdating algorithms to add/remove nodes efficiently.

In summary, the paper provides a rigorous numerical framework for solving the inverse eigenvalue problem for MOPs, demonstrating that while the problem is often ill-conditioned, specific Krylov and Core transformation strategies can yield stable and accurate results.