🔢 mathematics

Computing matrix functions associated with a Hermitian--definite pencil

This paper investigates the numerical evaluation of $Af(A^{-1}B)$ for Hermitian-definite pencils by analyzing problem conditioning and comparing algorithms based on Schur decomposition with either matrix square roots or Cholesky factorization, ultimately demonstrating that the Cholesky-based approach offers superior accuracy and efficiency.

Original authors: Dario A. Bini, Massimiliano Fasi, Bruno Iannazzo

Published 2026-05-25

📖 4 min read🧠 Deep dive

CC BY 4.0

Original authors: Dario A. Bini, Massimiliano Fasi, Bruno Iannazzo

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are a master chef trying to create a very specific, complex dish. You have two main ingredients: Ingredient A (which is perfectly stable and solid, like a firm block of cheese) and Ingredient B (which is a bit more variable, like a delicate spice blend).

Your goal is to mix them together in a very precise way to create a final sauce. The recipe calls for you to take Ingredient A, mix it with a special transformation of Ingredient B, and then apply a "flavor function" (let's call it $f$ ) to the result. In math-speak, this is written as $A \cdot f(A^{-1}B)$ .

The problem is that mixing these ingredients directly is messy. If you try to melt the cheese (invert A) and mix it with the spices first, you might lose the texture or introduce errors. The paper by Bini, Fasi, and Iannazzo is essentially a guide on how to mix these ingredients most efficiently and accurately without ruining the dish.

Here is the breakdown of their findings using simple analogies:

1. The Problem: The "Messy Kitchen"

The authors look at a mathematical operation used in many fields (like physics and engineering) where you need to combine two matrices (grids of numbers).

The Naive Approach: The easiest way to cook this is to just dump everything into a pot, stir it, and hope for the best. In math, this means calculating $A^{-1}B$ first. The paper shows this is like trying to bake a cake by mixing the flour and eggs after they've already been burnt. It's slow, and if your ingredients aren't perfect (a concept called "conditioning"), the final cake tastes terrible (high error).
The "Square Root" Approach: One way to fix this is to find the "square root" of the cheese block (Ingredient A) first. This is like grating the cheese into a fine powder before mixing. It works better than the naive approach, but grating cheese is hard work (computationally expensive) and can still be a bit shaky if the cheese is old.

2. The Solution: The "Cholesky" Shortcut

The authors propose a new, smarter way to cook. Instead of grating the cheese (finding the square root), they suggest using a Cholesky factorization.

The Metaphor: Imagine instead of grating the cheese, you simply slice it into perfect, neat triangular wedges. This is much faster to do and keeps the structure of the cheese intact.
The Result: They combine this "slicing" method with a standard mixing technique (Schur decomposition). They call this the Cholesky–Schur algorithm.

3. The Comparison: Who Wins?

The authors ran a "cooking competition" with five different methods:

Naive: The messy, direct mix. (Slowest and least accurate).
Square Root + Schur: The grating method. (Better, but still heavy lifting).
Double Schur: A method that slices the ingredients twice. (Good, but not the best).
Cholesky + Schur (The Winner): The slicing method.
Cholesky + Schur (for two solid ingredients): A variation if both ingredients are stable.

The Verdict:

Speed: The Cholesky methods were the fastest. They required about half the "effort" (computational steps) of the Square Root method.
Accuracy: The Cholesky methods were also the most precise. They produced a "sauce" that was much closer to the perfect theoretical taste, even when the ingredients were imperfect or "ill-conditioned" (like old, crumbly cheese).
Why? The math analysis showed that the "slicing" (Cholesky) method introduces fewer errors during the mixing process than the "grating" (Square Root) method.

4. What They Didn't Do

It is important to note what this paper doesn't say.

They did not test this on massive, sparse datasets (like a giant spreadsheet with mostly empty cells). Their methods are for "moderate-sized" matrices.
They did not apply this to specific real-world problems like medical imaging or climate modeling in this specific text. They focused purely on the mathematical mechanics of the mixing process itself.
They did not solve the problem of applying this to a single vector (a single line of data) rather than the whole matrix. That is a different challenge they mention for future work.

Summary

Think of this paper as a culinary guide for mathematicians. They took a complex recipe ( $Af(A^{-1}B)$ ), tested five different ways to prepare it, and proved that slicing the ingredients neatly (Cholesky factorization) is superior to grating them (Matrix Square Root). It's faster, it's cheaper to do, and it results in a much tastier (more accurate) final product.

Problem Statement
The paper addresses the numerical evaluation of the quantity $\phi(A, B) := Af(A^{-1}B)$ , where $A \in \mathbb{C}^{n \times n}$ is Hermitian positive definite, $B \in \mathbb{C}^{n \times n}$ is Hermitian, and $f$ is a function defined on the spectrum of $A^{-1}B$ . This problem arises in the context of Hermitian–definite matrix pencils $B - \lambda A$ and is central to computing various matrix means (e.g., arithmetic, harmonic, geometric) and solving systems of ordinary differential equations. The authors focus on Problem P1: computing the full matrix $\phi(A, B)$ for matrices of moderate size. A key challenge is to avoid the explicit formation of $A^{-1}B$ , which can be numerically unstable, and to minimize computational cost while maintaining numerical stability.

Methodology and Algorithms
The authors analyze the conditioning of the problem and propose several algorithms that leverage the spectral properties of the pencil. They utilize the equivalence $\phi(A, B) = A^{1/2}f(A^{-1/2}BA^{-1/2})A^{1/2}$ to transform the problem into evaluating a function of a Hermitian matrix. The study compares five specific implementations:

Naive: Computes $A^{-1}B$ explicitly and applies $f$ .
Sqrt Schur (Algorithm 1): Computes the matrix square root $A^{1/2}$ , forms $A^{-1/2}BA^{-1/2}$ , performs a Schur decomposition, applies $f$ to the eigenvalues, and reconstructs the result.
Schur Schur (Algorithm 2): Uses two Schur decompositions (one for $A$ and one for the transformed matrix) without computing a matrix square root explicitly.
Chol Schur (Algorithm 3): Replaces the matrix square root and inversion with the Cholesky factorization $A = R^*R$ . It computes $R^{-1}BR^{-*}$ , performs a Schur decomposition, and reconstructs the result.
Chol Schur PD (Algorithm 4): An extension of Algorithm 3 for the case where $B$ is also positive definite, utilizing Cholesky factorizations for both $A$ and $B$ .

The authors provide a detailed rounding error analysis in floating-point arithmetic for the square root-based approach (Algorithm 1) and the Cholesky-based approach (Algorithm 3). They derive explicit bounds for the forward error, decomposing the total error into a perturbation error (due to the initial factorization) and a roundoff error accumulated during the subsequent steps.

Key Contributions and Results

Conditioning Analysis: The authors derive explicit expressions for the Fréchet derivative of $\phi(A, B)$ and establish bounds for its condition number. These bounds depend on the condition numbers of $A$ and $B$ and the divided differences of the function $f$ and its dual.
Computational Complexity: Theoretical flop counts are provided. Algorithm 1 (Sqrt Schur) requires approximately $27n^3$ flops. Algorithm 3 (Chol Schur) reduces this to approximately $13n^3$ flops (or $(12 + 2/3)n^3$ with further optimizations), effectively halving the computational cost by avoiding the expensive matrix square root.
Error Analysis: The rounding error analysis reveals that the error bound for Algorithm 3 is theoretically superior to that of Algorithm 1. Specifically, the perturbation error in Algorithm 1 scales with $n^3$ (due to the matrix square root), whereas in Algorithm 3, it scales linearly with $n$ (due to Cholesky factorization). Furthermore, the analysis suggests that the error in Algorithm 3 is less sensitive to the condition number of $A$ compared to Algorithm 1.
Numerical Experiments: Experiments conducted on random positive definite matrices confirm the theoretical predictions.
- Accuracy: The Cholesky-based algorithms (Algorithms 3 and 4) demonstrate significantly higher accuracy, particularly when $A$ is ill-conditioned. Their forward errors grow more moderately with the condition number of $A$ compared to the square root-based and naive methods.
- Performance: Algorithms 3 and 4 are the fastest, outperforming the naive approach by 1–2 orders of magnitude and the square root-based approach by a noticeable margin.
- Comparison: The "Chol Schur" and "Chol Schur PD" variants exhibit nearly identical performance and accuracy, suggesting that the additional Cholesky factorization of $B$ in Algorithm 4 offers a minor theoretical improvement but is practically comparable to Algorithm 3.

Significance and Claims
The paper claims that algorithms based on the Cholesky factorization are both more accurate and more efficient than those based on the matrix square root for computing $\phi(A, B)$ . The authors assert that their numerical experiments confirm the theoretical analysis, showing that the Cholesky-based methods achieve a better balance of computational cost and numerical stability. They conclude that for moderate-sized matrices, the Cholesky–Schur approach (Algorithm 3) is the preferred method. The paper explicitly notes that it does not address the case of large, sparse matrices (Problem P2), where the resulting dense matrix cannot be stored, and suggests that such cases require different techniques like Krylov subspace methods, which are left for future work.

1. The Problem: The "Messy Kitchen"

2. The Solution: The "Cholesky" Shortcut

3. The Comparison: Who Wins?

4. What They Didn't Do

Summary

More like this