🔢 mathematics

The eigenvalue decomposition of normal matrices by the skew-symmetric part

This paper introduces a new, stable, and efficient method for computing the eigenvalue decomposition of dense real normal matrices by leveraging the decomposition of their skew-symmetric part, demonstrating competitive performance against existing algorithms and applicability to Riemannian barycenter computations on the special orthogonal group.

Original authors: Simon Mataigne, Kyle A. Gallivan

Published 2026-03-31

📖 4 min read🧠 Deep dive

CC BY 4.0

Original authors: Simon Mataigne, Kyle A. Gallivan

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 have a giant, complex puzzle made of numbers (a matrix). Your goal is to take this puzzle apart to see its fundamental building blocks (its eigenvalues and eigenvectors). This is a classic problem in mathematics, but for certain types of puzzles called Normal Matrices, there's a tricky catch: they are often messy and hard to solve directly.

This paper introduces a clever new shortcut. Instead of trying to solve the whole messy puzzle at once, the authors suggest looking at a specific, hidden "shadow" of the puzzle first.

Here is the breakdown of their method using simple analogies:

1. The Problem: The "Messy" Normal Matrix

Think of a Normal Matrix as a spinning top that is perfectly balanced but wobbles in a complex way. It has two types of movements:

Symmetric parts: Movements that look the same if you look in a mirror (like a perfect reflection).
Skew-symmetric parts: Movements that are the exact opposite of their mirror image (like a left hand vs. a right hand).

Traditionally, to understand how this top spins, mathematicians use a heavy, slow method called the QR algorithm. It's like trying to untangle a knot by pulling on every single string one by one. It works, but it takes a long time, especially if the top is huge.

2. The Insight: The "Skew-Symmetric Shadow"

The authors realized something brilliant: The "left-hand/right-hand" part (the skew-symmetric part) of the puzzle actually holds the secret to the whole thing.

Imagine the puzzle is a complex machine. The authors say, "Don't try to fix the whole machine. Just take out the gears that spin in opposite directions (the skew-symmetric part)."

In the world of math, this "opposite-spinning" part is actually much easier to solve than the whole machine. It's like finding a hidden key that fits a simpler lock.
Because the puzzle is "Normal," the structure of this simpler lock perfectly matches the structure of the complex machine.

3. The Method: The "Shadow Detective"

The new method works like a detective solving a crime by looking at a shadow:

Cast the Shadow: The computer first isolates the "skew-symmetric" part of the matrix. This is a quick, easy step.
Solve the Simple Puzzle: The computer solves this simpler puzzle. Because it's simpler, it can use very fast, pre-existing tools (like a specialized screwdriver) to get the answer almost instantly.
Reconstruct the Truth: The computer then uses the answer from the simple puzzle to figure out the answer for the complex original puzzle.
- The Catch: Sometimes, the shadow is a bit blurry (due to computer rounding errors). If the "spin speeds" of the puzzle are very similar, the shadow might get confused.
- The Fix: The authors built a "polishing step." If the shadow is blurry, they do a tiny bit of extra work to sharpen the image, ensuring the final answer is perfect.

4. Why This Matters: The "Special Orthogonal" Party

This method shines brightest when dealing with Orthogonal Matrices.

Analogy: Imagine a dance floor where everyone is spinning. In a "Random Orthogonal Matrix," almost everyone is spinning in pairs (complex numbers), and very few are standing still (real numbers).
The Advantage: Traditional methods try to handle the standing dancers and the spinning dancers all at once, which is slow. This new method realizes that the spinning dancers (the skew-symmetric part) tell you everything you need to know.
The Result: For these specific types of matrices, the new method is as fast as the fastest existing methods, but it's often much faster because it skips the heavy lifting of the traditional approach.

5. Real-World Application: Finding the "Average" Position

The paper ends with a cool application: The Riemannian Barycenter.

The Scenario: Imagine you have 100 different maps of the Earth, but each map is rotated slightly differently. You want to find the "average" map that represents the center of all these rotations.
The Bottleneck: To find this average, you have to constantly calculate the "logarithm" (a mathematical way to measure the difference) between rotations. This requires solving the matrix puzzle over and over again.
The Impact: By using this new "Shadow Detective" method, the computer can calculate these averages 2 to 4 times faster than before. This is huge for fields like robotics, computer vision, and statistics, where speed is everything.

Summary

The paper says: "Stop trying to solve the whole complex puzzle. Look at the 'opposite-spinning' shadow instead. It's easier to solve, and it reveals the solution to the whole thing almost instantly."

It's a smarter, faster way to untangle complex mathematical knots, especially when dealing with rotations and spins.

1. Problem Statement

The paper addresses the computation of the Eigenvalue Decomposition (EVD) or Real Schur Decomposition (RSD) for dense real normal matrices ( $A \in \mathbb{R}^{n \times n}$ where $AA^\top = A^\top A$ ).

Context: While efficient algorithms exist for Hermitian (symmetric) matrices and general matrices (via the QR algorithm), normal matrices present a unique challenge. Standard methods (like the shifted QR algorithm) reduce a normal matrix to an upper Hessenberg form and iterate until it becomes diagonal. However, this process ignores the inherent structure of normal matrices, leading to unnecessary computational costs.
Specific Challenge: The authors focus on cases where the matrix has few real eigenvalues (i.e., most eigenvalues are complex conjugate pairs), a property common in random orthogonal matrices (e.g., Haar-distributed). Existing methods often fail to exploit the specific structure of the skew-symmetric part of these matrices to achieve faster convergence or lower operation counts.

2. Methodology

The core insight of the paper is that for a normal matrix $A$ , the skew-symmetric part $\Omega = \frac{A - A^\top}{2}$ and the symmetric part $H = \frac{A + A^\top}{2}$ commute ( $H\Omega = \Omega H$ ). Consequently, they share the same invariant subspaces. The proposed method leverages the structure of $\Omega$ to recover the EVD of $A$ without performing a full Hessenberg reduction on $A$ .

The algorithm proceeds in three main phases:

A. Decomposition of the Skew-Symmetric Part

Tridiagonalization: Compute the skew-symmetric tridiagonal form of $\Omega$ using Householder reflectors.
Permutation: Apply an "even-odd" permutation to transform the tridiagonal skew-symmetric matrix into a block form containing a bidiagonal matrix $B$ .
SVD: Compute the Singular Value Decomposition (SVD) of the bidiagonal matrix $B$ . This yields the singular values (which correspond to the magnitudes of the imaginary parts of $A$ 's eigenvalues) and the singular vectors.
Initial Schur Vectors: Construct an initial matrix of Schur vectors $\tilde{Q}$ from the SVD components.

B. Classification and Refinement of Invariant Subspaces

The method categorizes the eigenvalues of $A$ into three sets based on the singular values of $\Omega$ :

Set 1 (Distinct Non-zero Imaginary Parts): If the imaginary parts of the eigenvalues are distinct, the invariant subspaces of $\Omega$ are identical to those of $A$ . No further computation is needed for these vectors.
Set 2 (Repeated Imaginary Parts): If multiple eigenvalues share the same imaginary part (a cluster), the invariant subspace of $\Omega$ is larger than that of $A$ . The algorithm extracts the corresponding block from $\tilde{Q}$ , computes the projection $M = V^\top A V$ , and solves a smaller symmetric skew-Hamiltonian Eigenvalue Problem (or a general normal EVP) to refine the vectors.
Set 3 (Real Eigenvalues): If $\Omega$ has a null space (zero singular values), these correspond to the real eigenvalues of $A$ . The algorithm extracts the null space vectors, projects $A$ onto this subspace ( $H = V^\top A V$ ), and solves a small symmetric EVP to find the real eigenvalues and refine the vectors.

C. Floating-Point Stability and Correction

The authors analyze the sensitivity of these subspaces to perturbations. They introduce Algorithm 4.1, which includes:

Clustering: Identifying clusters of singular values in $\Omega$ to determine which subspaces require refinement.
Correction Steps: If the separation between eigenvalues is small (leading to potential loss of accuracy), the algorithm performs a refinement step (e.g., using a Jacobi method) on the specific invariant subspaces to restore user-specified accuracy.

3. Key Contributions

Novel Algorithm: A new method (Algorithm 3.1/4.1) that computes the EVD of real normal matrices by decomposing the skew-symmetric part rather than the full matrix.
Theoretical Foundation: Rigorous proofs establishing the relationship between the invariance groups of the skew-symmetric part and the original normal matrix, specifically handling cases with repeated imaginary parts and real eigenvalues.
Complexity Analysis:
- Best Case: For matrices with distinct imaginary parts (common in random orthogonal matrices), the method requires $\sim \frac{14}{3}n^3$ floating-point operations (flops). This is equivalent to the cost of a Hessenberg factorization but avoids the expensive iterative QR steps required to diagonalize the Hessenberg form.
- Comparison: It is significantly faster than standard LAPACK routines (gees) for large $n$ , offering speedups of up to 10x for $n \approx 100$ and 2x for $n=10,000$ .
Stability Analysis: A detailed sensitivity analysis showing that while the method is generally stable, specific eigenvalue configurations (clustered imaginary parts) can degrade accuracy. The paper provides theoretical bounds and practical correction steps to mitigate this.

4. Numerical Results

Accuracy: Experiments on random normal matrices and orthogonal matrices show that the method achieves machine precision accuracy ( $\sim 10^{-15}$ to $10^{-16}$ in double precision) for the Schur form and eigenvalues.
Performance:
- The method outperforms LAPACK's gees (Schur decomposition) and hessenberg routines in runtime for large matrices.
- It is particularly efficient for Haar-distributed orthogonal matrices, where the probability of having distinct imaginary parts is 1, allowing the algorithm to skip the most expensive refinement steps.
Application: The method was applied to compute Riemannian barycenters on the special orthogonal group $SO(n)$. In this iterative application, the speedup was significant (factors of 1.7 to 4.4), as the matrix logarithm (requiring EVD) is the bottleneck.

5. Significance

Efficiency for Normal Matrices: This work provides the first computationally efficient method specifically tailored for dense real normal matrices that exploits their structure to bypass the standard Hessenberg reduction bottleneck.
Manifold Optimization: The method is highly relevant for optimization on matrix manifolds (like $SO(n)$ and Stiefel manifolds), where computing matrix logarithms and barycenters is frequent. The speedup in these iterative algorithms can drastically reduce the time required for statistical analysis on manifolds.
Robustness: By providing a framework to detect and correct for eigenvalue clustering, the method ensures high accuracy even in worst-case scenarios, making it a viable alternative to general-purpose solvers.

In summary, the paper proposes a specialized, high-performance algorithm for normal matrices that transforms the problem into a bidiagonal SVD and small symmetric EVPs, offering superior speed and comparable accuracy to state-of-the-art general solvers, particularly for orthogonal and random normal matrices.