🔢 mathematics

Preconditioned normal equations for solving discretised partial differential equations

This paper introduces a "normal" preconditioning strategy for solving non-symmetric linear systems from PDE discretizations by constructing preconditioners based on the associated normal PDE, demonstrating its effectiveness in achieving fast and stable convergence for convection-diffusion problems.

Original authors: Lorenzo Lazzarino, Yuji Nakatsukasa, Umberto Zerbinati

Published 2026-07-28

📖 3 min read🧠 Deep dive

CC BY 4.0

Original authors: Lorenzo Lazzarino, Yuji Nakatsukasa, Umberto Zerbinati

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 trying to solve a giant, tangled knot of equations that describe how heat spreads through a metal plate or how smoke drifts in a room. In the world of science and engineering, these are called Partial Differential Equations (PDEs). To solve them on a computer, scientists break the smooth, continuous world into a grid of tiny dots, turning the smooth physics into a massive list of numbers. This list is a "linear system," a giant puzzle where you have to find the right numbers to make everything balance.

Usually, these puzzles are tricky because the rules aren't symmetrical; the way heat moves forward isn't exactly the same as the way it would move backward if you rewound time. Because of this, the standard computer tools used to solve them can sometimes get stuck, spin their wheels, or take forever to find the answer. Scientists have developed a special toolbox called "Krylov subspace methods" to tackle these, but there's a debate about which tool is best. One popular tool, GMRES, is like a very thorough detective who checks every possible angle, but it can be slow and complicated to set up. Another tool, called CGNR, is like a sprinter who runs straight ahead, but it usually only works well if the puzzle is perfectly symmetrical. The big question is: Can we make the sprinter run fast and straight even when the puzzle is messy and asymmetrical?

This paper by Lorenzo Lazzarino, Yuji Nakatsukasa, and Umberto Zerbinati says, "Yes, but we have to change the track." They propose a clever trick called "preconditioned normal equations." Instead of trying to force the sprinter to run on the messy, asymmetrical track, they suggest transforming the problem into a new, perfectly symmetrical version of itself (called the "normal equation") and then running the sprinter on that. The magic lies in how they build the track. They introduce a new concept called a "normal preconditioner." Think of a preconditioner as a pair of special glasses that helps the computer see the problem clearly. The authors show that for this new symmetrical track, the "perfect glasses" aren't unique. In fact, there are many different pairs of glasses that work just as well, as long as they make the numbers in the puzzle look like they are all clustered around the number 1.

The researchers tested this idea on a classic problem: advection-diffusion, which is like tracking a puff of smoke being blown by the wind (advection) while also spreading out on its own (diffusion). They found that by looking at a "normal" version of the physics problem, they could build a preconditioner that made the computer solve the puzzle very fast in many cases. However, the results were nuanced: while the method achieved mesh-independent convergence for constant winds, the iteration counts did grow significantly when the wind was constant but the stabilization parameters were large, and the method struggled on coarser grids. Furthermore, when they tested it on complex, swirling winds (recirculating flows), the standard sparse version of their preconditioner degraded markedly. To fix this, they had to switch to a more computationally expensive "matrix-free" approach that used a global projection, which successfully handled the swirling winds but required more inner calculations. The paper suggests that this isn't just a theoretical curiosity; it's a practical, robust alternative to traditional tools, especially for difficult PDEs where fine-tuning classical preconditioners is hard. By treating the problem through this new "normal" lens, they turned a difficult, asymmetrical race into a smooth sprint, provided the right track adjustments are made for the specific wind conditions.

Technical Summary: Preconditioned Normal Equations for Solving Discretised Partial Differential Equations

Problem Statement
The paper addresses the challenge of solving large-scale, non-symmetric square linear systems ($Ax=b$) arising from the discretization of Partial Differential Equations (PDEs). While the Generalized Minimal Residual (GMRES) method is the standard choice for such systems, the authors argue for the viability of Conjugate Gradient (CG)-like methods applied to the normal equation ( $A^T A x = A^T b$ ). Specifically, they focus on the Conjugate Gradient Normal Residual (CGNR) and LSQR algorithms. A primary barrier to this approach has been the "squared condition number" of the normal equation ( $A^T A$ ) and the difficulty in constructing effective preconditioners that account for the non-symmetric nature of $A$ . Furthermore, standard preconditioning strategies for $Ax=b$ (aiming to cluster eigenvalues of $P^{-1}A$ ) do not necessarily translate to effective preconditioning for the normal equation.

Methodology and Theoretical Framework
The authors introduce a concept termed "normal preconditioning," which redefines what constitutes a "good" preconditioner in the context of the normal equation.

The "Cross" Effect: The paper establishes that for a preconditioner $G = P^T P$ to be effective for the normal equation (specifically for left-preconditioned CGNR), the matrix $AP^{-1}$ must have clustered singular values, not necessarily clustered eigenvalues. This creates a "cross" effect: a good left preconditioner for CGNR requires a good right preconditioner for the original system $Ax=b$ (in terms of singular values).
Functional Analytic Perspective: The authors highlight a functional inconsistency in the standard normal equation $A^T A x = A^T b$ when derived from PDEs. The operator $A$ maps a discrete space $V_h$ to its dual $V'_h$ , while $A^T$ maps $V_h$ to $V'_h$ as well. Consequently, $A^T A$ maps $V_h \to V'_h \to V_h$ , but the intermediate spaces are inconsistent with the continuous problem's function spaces.
The Modified Normal Equation: To resolve this, the authors propose a modified normal equation:
$A^T T A x = A^T T b$
where $T$ is a discrete Riesz map associated with a specific inner product on the function space (e.g., the $H^1_0$ inner product). This ensures the operator $A^T T A$ maps $V_h \to V_h$ consistently.
Deriving the "Normal" PDE: By selecting an appropriate Riesz map $T$ , the authors show that the operator $A^T T A$ corresponds to the discretization of a specific continuous PDE, termed the "normal PDE." For the advection-diffusion equation in the advection-dominated regime, choosing the $H^1_0$ inner product yields a normal PDE that is a reaction-diffusion equation.

Key Contributions

Definition of Normal Preconditioners: The paper formalizes the notion that an ideal preconditioner for the normal equation is not unique (unlike the inverse of $A$ ). Instead, any $P$ that makes $AP^{-1}$ orthogonal (or close to orthogonal) with respect to the inner product induced by $T$ is an ideal preconditioner. This significantly expands the search space for effective preconditioners.
Construction Strategy: A general framework is proposed:
1. Identify the correct Riesz map $\tau$ for the continuous PDE.
2. Derive the associated "normal PDE" from the operator $A^T T A$ .
3. Discretize this normal PDE to form a preconditioner $G = P^T T P$ .
Practical Implementation: The authors demonstrate that the resulting normal PDE (reaction-diffusion) is symmetric positive definite (SPD) and can be efficiently inverted using standard algebraic multigrid (AMG) methods, such as PETSc GAMG, avoiding the need for complex, problem-specific non-symmetric preconditioners.

Numerical Results
The paper presents numerical experiments on advection-diffusion problems with varying viscosity ( $\nu$ ) and mesh sizes, comparing the proposed preconditioned CGNR against FGMRES preconditioned with Approximate Ideal Restriction (AIR).

Mesh Independence: The proposed method achieves mesh-independent iteration counts for both constant and recirculating wind fields. As the mesh is refined, the iteration count remains stable or decreases.
Robustness: For highly singularly perturbed problems (small $\nu$ ), the preconditioned CGNR converges within a small number of iterations (typically under 20, often under 10), whereas standard FGMRES with AIR often stagnates or fails to converge within 40 iterations on coarse meshes.
Stabilization Consistency: The experiments show that the preconditioner must be stabilized consistently with the original operator (e.g., using the same Streamline-Upwind/Petrov-Galerkin (SUPG) terms) to maintain spectral equivalence.
Double-Glazing Benchmark: For a recirculating flow problem (Elman, Silvester, and Wathen benchmark), the authors employ a matrix-free implementation of the projected operator. The outer CGNR iterations converge in as few as 2 iterations on fine meshes, demonstrating the method's effectiveness even for complex, non-constant wind fields.

Significance and Claims
The paper positions itself as a "proof of concept" rather than a claim of universal superiority over all existing methods. The authors assert that:

Feasibility: Solving non-symmetric PDE systems via the preconditioned normal equation is a viable and competitive alternative to classical Krylov methods like GMRES.
Simplicity: The approach allows practitioners to leverage simpler, robust SPD preconditioning strategies (like AMG) for difficult non-symmetric problems, bypassing the need for intricate, ad-hoc non-symmetric preconditioner tuning.
Non-Uniqueness: The non-uniqueness of ideal normal preconditioners provides a rich design space for constructing efficient solvers.
Stability: When combined with iterative refinement, preconditioned CGNR/LSQR can provide backward stable solutions.

The authors acknowledge that while their method is effective, it is not necessarily "better" than a perfectly tuned geometric or algebraic preconditioner for the original system. However, it offers a robust, systematic framework for problems where such tuning is difficult or infeasible. Future work is planned to extend this framework to other singularly perturbed problems (Oseen, Helmholtz) and to explore mixed reformulations to handle the density of the normal operator.

More like this