🔢 mathematics

Finding roots of complex analytic functions via generalized colleague matrices

This paper presents a numerical scheme for finding all roots of complex analytic functions in a square domain by generalizing the classical colleague matrix approach through the construction of well-conditioned polynomial bases satisfying three-term recurrences and introducing a specialized QR algorithm to compute the eigenvalues of the resulting generalized colleague matrices.

Original authors: Hanwen Zhang, Vladimir Rokhlin

Published 2026-04-13

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Hanwen Zhang, Vladimir Rokhlin

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 treasure hunter trying to find hidden islands (roots) in a vast, foggy ocean (the complex plane). The ocean is filled with mysterious functions that twist and turn, making it impossible to see where the islands are just by looking.

This paper presents a new, highly efficient map-making technique to find all the islands in a specific square-shaped region of the ocean, no matter how tricky the water gets.

Here is the breakdown of their method using simple analogies:

1. The Problem: Finding Roots in the Fog

In math, a "root" is a spot where a function equals zero. On a simple straight line (like a real number line), we have old, reliable tools to find these spots. But in the complex plane (which has two dimensions: left/right and up/down), the water is choppy. Traditional tools often get lost or break down because the "waves" (singularities) can be just outside the area you are looking at, messing up your calculations.

2. The Solution: Approximating with a "Magic Lattice"

Instead of trying to solve the complex function directly, the authors suggest building a simpler, fake version of it—a polynomial—that behaves almost exactly like the real thing inside your square domain.

The Old Way: Imagine trying to draw a perfect circle using only straight lines (monomials). It's possible, but you need thousands of lines, and the drawing gets wobbly and inaccurate.
The New Way: They use a special set of "building blocks" (a polynomial basis) that fit together perfectly, like a high-quality jigsaw puzzle. These blocks are defined by a simple rule called a three-term recurrence. Think of it like a family tree where every new generation is just a mix of the two parents before it. This makes the math very stable and easy to handle.

The Secret Sauce: Usually, these building blocks only work well on straight lines. The authors discovered a trick: if you assign random weights to your measurements (like giving random importance to different points on the map), these blocks suddenly work beautifully in the complex plane too. It's like discovering that if you shake the dice just right, you can build a stable tower out of sand.

3. The "Generalized Colleague Matrix": The Treasure Map

Once they have built this perfect polynomial approximation, they need to find its roots.

In the old days, mathematicians built a "Companion Matrix" (a giant grid of numbers) where the roots were hidden as "eigenvalues" (special numbers that reveal the grid's secrets).
The authors built a "Generalized Colleague Matrix." Imagine this as a special, slightly distorted treasure map. It's not a full, messy grid; it's a streamlined, structured map (mostly empty space with a few key numbers) that holds the location of every root.

4. The "Complex Orthogonal QR" Algorithm: The Scavenger Hunt

Now that they have the map, they need to read it. They use a special algorithm called the QR algorithm to extract the roots.

Think of this algorithm as a highly organized scavenger hunt. It takes the messy map and performs a series of precise rotations (like turning a Rubik's cube) to line up the numbers until the roots pop out clearly.
Because the map has a special structure, this hunt is incredibly fast and doesn't lose accuracy, even when the numbers get huge or tiny. It's like having a GPS that never loses signal, even in a storm.

5. The Adaptive Strategy: Zooming In

Sometimes, the ocean is too rough for one big map. The authors also created an adaptive version of their method.

Imagine you are looking at a coastline. If the coast is smooth, you can draw it on a small piece of paper. But if there's a jagged cliff or a cluster of islands, you zoom in and draw that specific area on a larger, more detailed map.
Their algorithm automatically detects where the function is "rough" (hard to approximate) and splits the square domain into smaller squares, solving the problem piece by piece. This ensures they never miss a root, even if they are clustered tightly together.

Why This Matters

This method is a game-changer because:

It finds everything: It guarantees finding all roots in the area, not just the easy ones.
It's fast: It scales efficiently, meaning it can handle very complex functions without taking forever.
It's robust: It works even when there are dangerous "storms" (singularities) just outside the area you are studying.

In a nutshell: The authors invented a way to turn a chaotic, 2D mathematical puzzle into a structured, solvable grid by using random weights to create stable building blocks, and then using a specialized "scavenger hunt" algorithm to find the hidden answers. It's a bit like using a random shuffle to organize a deck of cards so perfectly that you can instantly find the Ace of Spades every time.

1. Problem Statement

The paper addresses the numerical challenge of finding all roots of an analytic function $f(z)$ within a square domain in the complex plane.

Context: Finding roots of polynomials on the real line is a well-established problem, often solved by approximating the function with Chebyshev polynomials and computing the eigenvalues of a "colleague matrix" (a structured matrix similar to a companion matrix).
Challenge: Extending this approach to the complex plane is difficult because:
1. There are no orthogonal polynomials on complex domains that satisfy the three-term recurrence relations required for efficient colleague matrix construction. The self-adjointness of multiplication by $x$ (which guarantees three-term recurrences in real intervals) does not hold for complex $z$ .
2. Attempts to use standard "complex inner products" (omitting conjugation) lead to highly unstable bases (e.g., the integral of $z^2$ over a square centered at the origin is zero, causing breakdown).
3. Standard QR algorithms for finding eigenvalues may suffer from numerical instability when applied to the specific structures arising in the complex domain.

2. Methodology

The authors propose a scheme that generalizes the classical real-line approach to the complex plane through three main components:

A. Construction of Special Polynomial Bases

Instead of seeking classical orthogonal polynomials, the authors construct a basis $\{P_j(z)\}$ that satisfies a three-term recurrence relation:
$zP_j(z) = \beta_j P_{j-1}(z) + \alpha_{j+1} P_j(z) + \beta_{j+1} P_{j+1}(z)$

Mechanism: They utilize complex orthogonalization with respect to a complex inner product defined by random weights.
- Given a set of points $\{z_i\}$ on the boundary of the domain and random complex weights $\{w_i\}$ , they define an inner product $[u, v]_w = \sum w_i u_i v_i$ (without conjugation).
- By applying a Gram-Schmidt-like process to the vector $b=(1, \dots, 1)^T$ using a diagonal matrix $Z = \text{diag}(z_i)$ , they generate vectors $q_j$ that correspond to polynomial values $P_j(z_i)$ .
Key Observation: While the resulting basis is not orthogonal in the classical sense, using random weights prevents the basis from becoming ill-conditioned. The condition number of the basis grows almost linearly with the order $n$ , making it suitable for high-degree approximations.

B. Generalized Colleague Matrices

Once the function $f(z)$ is approximated by a polynomial $p(z) = \sum c_j P_j(z)$ in this basis, the roots of $p(z)$ are found by computing the eigenvalues of a Generalized Colleague Matrix $C$ :
$C = A + e_n q^T$

$A$ is a complex symmetric tridiagonal matrix constructed from the recurrence coefficients $\alpha$ and $\beta$ .
$e_n q^T$ is a rank-1 update derived from the expansion coefficients $c_j$ .
Unlike classical colleague matrices (which are Hermitian), $C$ is complex symmetric plus a rank-1 update. This structure allows the matrix to be represented compactly by a small number of vectors (generators) rather than a dense matrix.

C. Complex Orthogonal QR Algorithm

To find the eigenvalues of $C$ , the authors develop a specialized Complex Orthogonal QR algorithm.

Structure Preservation: The algorithm maintains the structure of $C$ (complex symmetric + rank-1) throughout the iterations.
Transformations: Instead of unitary rotations (SU(2)) used in the Hermitian case, the algorithm uses complex orthogonal matrices ( $Q^T = Q^{-1}$ ) derived from Lemma 5 in the paper.
Stability Correction: A critical contribution is a correction step to handle cases where the rank-1 update vector $q$ has very large norms (which can occur when approximating analytic functions). The algorithm explicitly corrects rounding errors that would otherwise destroy the precision of the eigenvalues, ensuring structured backward stability.
Complexity: The algorithm finds all $n$ roots in $O(n^2)$ operations.

D. Adaptive Rootfinding

The paper presents two algorithms:

Non-adaptive: Uses a fixed polynomial order $n$ over the entire square.
Adaptive: Recursively subdivides the square domain into smaller squares. If the polynomial expansion of order $n_{exp}$ does not achieve the desired accuracy on a sub-square, that square is divided further. This handles functions with nearby singularities or rapid oscillations efficiently.

3. Key Contributions

Generalization of Colleague Matrices: Successfully extended the colleague matrix concept from real intervals to complex square domains by constructing polynomial bases via complex orthogonalization with random weights.
Stable Basis Construction: Demonstrated empirically that random weights in the complex inner product yield well-conditioned bases satisfying three-term recurrences, overcoming the theoretical barrier of non-existence of orthogonal polynomials in the complex plane.
Special-Purpose QR Algorithm: Developed a structured QR algorithm for complex symmetric matrices with rank-1 updates that is $O(n^2)$ and numerically stable, including specific corrections for large rounding errors.
Adaptive Framework: Provided a robust adaptive strategy for rootfinding that automatically refines the domain near singularities or clustering roots.

4. Numerical Results

The authors tested the algorithms on several functions:

$f_{cosh}$ : A function with a pole near the square. The method successfully found roots despite the singularity, with coefficients decaying geometrically.
$f_{poly}$ : A polynomial with simple roots. The method achieved full machine precision (double and extended) regardless of the magnitude of the vector $q$ or the polynomial order.
$f_{mult}$ : A polynomial with multiple roots. The error scaled as $u^{1/m}$ (where $m$ is multiplicity), consistent with theoretical expectations for multiple roots.
$f_{clust}$ : A function with clustering zeros near a singularity. The adaptive algorithm successfully resolved the roots by subdividing the domain, demonstrating that higher expansion orders reduce the number of subdivisions but may slightly reduce root accuracy due to the Maximum Principle (larger domains $\to$ larger function values $\to$ larger absolute errors).
$f_{entire}$ : An entire function over a large domain. The method found all 150 roots efficiently.

Performance: The algorithms are fast (milliseconds to seconds for $n \approx 100$ ) and exhibit structured backward stability in practice, even though a formal proof of stability for the complex orthogonal case is not yet provided.

5. Significance

This work provides a robust, efficient, and general-purpose tool for rootfinding of analytic functions in the complex plane.

Theoretical Impact: It bridges a gap in numerical analysis by showing how to construct useful polynomial bases in the complex plane without classical orthogonality.
Practical Impact: The method is superior to naive approaches (like discretizing the boundary and solving non-linear systems) because it finds all roots simultaneously with high precision and provides error estimates.
Future Work: The authors note that while the method is currently demonstrated on squares, the underlying theory (complex orthogonalization) suggests it can be generalized to arbitrary compact domains. The primary open problems are the theoretical understanding of why random weights yield well-conditioned bases and the formal proof of the QR algorithm's stability.

The code for these algorithms is made available publicly, facilitating further research and application in fields requiring complex rootfinding (e.g., fluid dynamics, quantum mechanics, and control theory).