Preconditioned Normal Equations

• Preconditioned normal equations are linear systems derived by applying preconditioning to the normal-equation formulation to enhance convergence and numerical stability.
• They play a key role in Krylov methods, randomized sketching, and mixed-precision computations, effectively reducing the condition number and spectral issues.
• Finite-precision analyses drive the use of iterative refinement and restarting strategies to achieve backward stability and maintain attainable accuracy.

Preconditioned normal equations are linear systems obtained by applying preconditioning to a normal-equation formulation of either a nonsymmetric square problem or a full-column-rank least-squares problem. For a square system $A x=b$, the basic normal equations are $A^T A x=A^T b$; for least squares, they arise from minimizing $\|A x-b\|_2$. The modern literature treats this formulation not merely as a reduction to a symmetric positive-definite system, but as a framework for Krylov methods, randomized sketching, mixed-precision computation, PDE-derived operators, and Schur-complement systems in optimization. A central theme is that preconditioning can recover fast convergence while finite-precision analysis determines when such methods are backward stable and when additional refinement is required (Epperly et al., 25 Feb 2025, Garrison et al., 17 Mar 2026, Lazzarino et al., 24 Feb 2025).

1. Algebraic forms and basic terminology

For a nonsingular square matrix $A\in\mathbb R^{n\times n}$, the unpreconditioned normal equations are

$A^T A\,x=A^T b.$

For a full-column-rank least-squares problem with $A\in\mathbb R^{m\times n}$, $m\ge n$, the same equation characterizes the minimizer of $\|A x-b\|_2$. Because $A^T A$ is symmetric positive definite in the full-rank case, it is natural to use CG-type methods; however, the squaring of the condition number makes preconditioning central (Ipsen, 24 Jul 2025, Lazzarino et al., 24 Feb 2025).

The recent least-squares literature distinguishes several preconditioned variants:

Variant	System	Characteristic
PNE	$A_p^T A_p\,y=A_p^T b,\quad R_s x=y$	Symmetrically preconditioned normal equations
HPNE	$A^T A x=A^T b$0	Half-preconditioned; generally nonsymmetric
SNE	$A^T A x=A^T b$1	Seminormal equations from thin QR
NNE	$A^T A x=A^T b$2	Uses any full-rank $A^T A x=A^T b$3 with the same column space as $A^T A x=A^T b$4

These formulations are analyzed for perturbation behavior in mixed precision and under randomized sketch-based preconditioning (Garrison et al., 17 Mar 2026).

For square nonsymmetric systems, right preconditioning is commonly written with $A^T A x=A^T b$5, so that

$A^T A x=A^T b$6

With both left and right preconditioners one may write

$A^T A x=A^T b$7

and in the notation of Epperly, Greenbaum, and Nakatsukasa one typically takes $A^T A x=A^T b$8, $A^T A x=A^T b$9, choosing $\|A x-b\|_2$0 to control $\|A x-b\|_2$1. Their basic assumptions include

$\|A x-b\|_2$2

hence $\|A x-b\|_2$3 and $\|A x-b\|_2$4 (Epperly et al., 25 Feb 2025).

2. Krylov formulations and preconditioning strategies

LSQR is algebraically equivalent to CG applied to the normal equations, but it is implemented through Lanczos bidiagonalization of the matrix or operator used in the least-squares formulation. In the right-preconditioned setting, one applies LSQR to

$\|A x-b\|_2$5

whose normal equations are

$\|A x-b\|_2$6

The matrix-vector operations are therefore

$\|A x-b\|_2$7

so each Krylov step uses four primitive solves or multiplications: $\|A x-b\|_2$8, $\|A x-b\|_2$9, $A\in\mathbb R^{n\times n}$0, and $A\in\mathbb R^{n\times n}$1. The Lanczos bidiagonalization recurrences are the standard LSQR recurrences with $A\in\mathbb R^{n\times n}$2 replaced by the composite operator $A\in\mathbb R^{n\times n}$3, and at iteration $A\in\mathbb R^{n\times n}$4 the method forms a $A\in\mathbb R^{n\times n}$5 bidiagonal matrix $A\in\mathbb R^{n\times n}$6 whose singular-value approximation yields $A\in\mathbb R^{n\times n}$7, hence $A\in\mathbb R^{n\times n}$8 (Epperly et al., 25 Feb 2025).

For normal equations arising from PDE discretization, CGNE and LSQR are often described as residual-minimizing methods in the original $A\in\mathbb R^{n\times n}$9-norm: $A^T A\,x=A^T b.$0 with the standard CG estimate

$A^T A\,x=A^T b.$1

This makes preconditioning the normal operator, rather than $A^T A\,x=A^T b.$2 directly, a natural design target when one can construct an effective SPD approximation to $A^T A\,x=A^T b.$3 or $A^T A\,x=A^T b.$4 (Lazzarino et al., 24 Feb 2025).

A recurrent comparison is with GMRES on the right-preconditioned nonsymmetric system

$A^T A\,x=A^T b.$5

For GMRES, the Saad–Schultz bound states that after $A^T A\,x=A^T b.$6 steps

$A^T A\,x=A^T b.$7

so $A^T A\,x=A^T b.$8 must approximate $A^T A\,x=A^T b.$9 on the spectrum of $A\in\mathbb R^{m\times n}$0. The contrast emphasized in the stability analysis is that LSQR-based methods only require $A\in\mathbb R^{m\times n}$1, whereas GMRES additionally benefits from a tightly clustered spectrum in the complex plane, not wrapping around the origin, and from well-conditioned eigenvectors (Epperly et al., 25 Feb 2025).

3. Finite-precision limits and attainable accuracy

A decisive issue is attainable accuracy in floating-point arithmetic. Even with a “near-perfect” preconditioner for which $A\in\mathbb R^{m\times n}$2, finite-precision LSQR can stagnate far above the backward-stable level

$A\in\mathbb R^{m\times n}$3

The reported experiments show that preconditioned LSQR may produce errors many orders of magnitude larger than classical direct methods, and that both the residual and the forward error can stall many orders of magnitude above $A\in\mathbb R^{m\times n}$4 (Epperly et al., 25 Feb 2025).

Two remedies are identified. The first is iterative refinement in the original system, with preconditioned LSQR used as the inner solver. In exact arithmetic, the refinement loop computes

$A\in\mathbb R^{m\times n}$5

approximately solves

$A\in\mathbb R^{m\times n}$6

by LSQR, and updates

$A\in\mathbb R^{m\times n}$7

The second is restarting preconditioned LSQR when the backward-error estimate

$A\in\mathbb R^{m\times n}$8

fails to decrease sufficiently over a fixed number of iterations. In the cited backward-error analysis, if $A\in\mathbb R^{m\times n}$9, each refinement step satisfies

$m\ge n$0

with $m\ge n$1. The resulting theorem states that if each inner LSQR solve takes

$m\ge n$2

iterations, then after

$m\ge n$3

refinement steps the overall procedure is backward stable, with only $m\ge n$4 calls to $m\ge n$5, $m\ge n$6, $m\ge n$7, and $m\ge n$8 (Epperly et al., 25 Feb 2025).

The same phenomenon appears for preconditioned CG on SPD normal equations. Periodic iterative refinement—apply PCG for $m\ge n$9 steps, recompute the residual in high precision, and repeat—restores

$\|A x-b\|_2$0

at essentially the same cost. This extends the attainable-accuracy picture beyond LSQR to the broader family of Krylov solvers on normal equations (Epperly et al., 25 Feb 2025).

A common source of confusion is the relation between convergence and accuracy. Small $\|A x-b\|_2$1 is sufficient for fast Krylov convergence, but it is not sufficient by itself for straightforward finite-precision LSQR to reach the backward-stable regime. By contrast, left-preconditioned LSQR appears empirically to attain a backward-stable solution without refinement when $\|A x-b\|_2$2 is small, “say $\|A x-b\|_2$3,” in

$\|A x-b\|_2$4

iterations, but the cited work states explicitly that no proof is yet available and that a complete analysis remains open (Epperly et al., 25 Feb 2025).

4. Randomized and mixed-precision preconditioners

A prominent recent direction computes the preconditioner from a randomized sketch. For full-column-rank $\|A x-b\|_2$5, one forms a “fast-JL” sketch

$\|A x-b\|_2$6

where $\|A x-b\|_2$7 is a random sign diagonal, $\|A x-b\|_2$8 is a real or complex Fourier or Hadamard transform, and $\|A x-b\|_2$9 samples $A^T A$0 rows from the identity. In lower precision $A^T A$1, one computes the thin QR factorization

$A^T A$2

promotes $A^T A$3 to working precision, and sets $A^T A$4. If

$A^T A$5

then with probability $A^T A$6,

$A^T A$7

This construction underlies both symmetrically preconditioned normal equations and half-preconditioned normal equations in mixed precision (Garrison et al., 17 Mar 2026).

The associated perturbation bounds are formulated in terms of the relative normal-equations residual

$A^T A$8

For PNE, the computed solution satisfies a bound scaled by $A^T A$9, growing linearly in $A_p^T A_p\,y=A_p^T b,\quad R_s x=y$0, with the max-term equal to $A_p^T A_p\,y=A_p^T b,\quad R_s x=y$1 when $A_p^T A_p\,y=A_p^T b,\quad R_s x=y$2 and growing like $A_p^T A_p\,y=A_p^T b,\quad R_s x=y$3 when $A_p^T A_p\,y=A_p^T b,\quad R_s x=y$4. For HPNE, the qualitative dependence is the same: weak on preconditioner error, but growing with large $A_p^T A_p\,y=A_p^T b,\quad R_s x=y$5. The seminormal equations inherit the $A_p^T A_p\,y=A_p^T b,\quad R_s x=y$6 sensitivity, while the not-normal equations also display residual-dependent growth. A significant conclusion is that the conditioning depends only mildly on the quality of the preconditioner, but it does depend on the size of the least-squares residual—even if the normal equations do not originate from a least-squares problem (Garrison et al., 17 Mar 2026).

The mixed-precision implementation includes an automatic precision-selection rule. In single precision one estimates $A_p^T A_p\,y=A_p^T b,\quad R_s x=y$7 via a fast $A_p^T A_p\,y=A_p^T b,\quad R_s x=y$8-norm condition estimator; the sketch precision is then chosen as half if $A_p^T A_p\,y=A_p^T b,\quad R_s x=y$9, single if $A^T A x=A^T b$00, and double if $A^T A x=A^T b$01 or overflow. The PNE or HPNE solve is performed in double precision. In experiments with $A^T A x=A^T b$02, $A^T A x=A^T b$03, $A^T A x=A^T b$04, the new bounds track the observed $A^T A x=A^T b$05 versus $A^T A x=A^T b$06, and for $A^T A x=A^T b$07, PNE and HPNE attain the same accuracy as MATLAB’s backslash. In mixed precision with $A^T A x=A^T b$08, $A^T A x=A^T b$09 or $A^T A x=A^T b$10, $A^T A x=A^T b$11, and a single-precision preconditioner, the new bounds remain informative and for $A^T A x=A^T b$12, mixed-precision PNE and HPNE again match mldivide’s accuracy. On NVIDIA H100 GPUs, mixed precision delivers up to a $A^T A x=A^T b$13–$A^T A x=A^T b$14 speedup over double-PNE and often out-runs the QR solver for $A^T A x=A^T b$15 up to a few thousand (Garrison et al., 17 Mar 2026).

A closely related analysis of randomized preconditioned normal equations emphasizes that, with an effective preconditioner, the solution accuracy is almost as accurate as the QR-based MATLAB backslash even for highly ill-conditioned matrices, and that the perturbation bound reduces to the standard least-squares form when $A^T A x=A^T b$16 and $A^T A x=A^T b$17. Numerical experiments reported there extend to matrices with $A^T A x=A^T b$18, with observed relative errors tracking the QR-based solution and plateauing at $A^T A x=A^T b$19 (Ipsen, 24 Jul 2025).

5. PDE-derived and structure-exploiting preconditioners

In PDE settings, preconditioned normal equations are often constructed from the underlying operator rather than from the matrix alone. For a differential operator $A^T A x=A^T b$20, the “normal PDE” is based on the composition $A^T A x=A^T b$21, which is self-adjoint and elliptic. After discretization, one solves

$A^T A x=A^T b$22

with a preconditioner

$A^T A x=A^T b$23

where $A^T A x=A^T b$24 is a Riesz map associated with the chosen inner product. The spectral rationale is that $A^T A x=A^T b$25 is spectrally equivalent to the Gram matrix of $A^T A x=A^T b$26 in the $A^T A x=A^T b$27-inner product, so its eigenvalues are the squares of the $A^T A x=A^T b$28-singular values of $A^T A x=A^T b$29. An ideal normal preconditioner is therefore any $A^T A x=A^T b$30 for which $A^T A x=A^T b$31 is $A^T A x=A^T b$32-orthogonal (Lazzarino et al., 24 Feb 2025).

For the space-time fractional advection-diffusion equation, Zhao, Jin, and Lin derive a nonsymmetric Toeplitz-like matrix $A^T A x=A^T b$33 and form the normal equations

$A^T A x=A^T b$34

Their preconditioner is $A^T A x=A^T b$35, where $A^T A x=A^T b$36 is obtained by band truncation of the Toeplitz blocks. Because $A^T A x=A^T b$37, the eigenvalues of the preconditioned normal matrix are clustered about $A^T A x=A^T b$38, with

$A^T A x=A^T b$39

The reported PCGNR iteration counts for Example 4.1, with bandwidth $A^T A x=A^T b$40, range from $A^T A x=A^T b$41 to $A^T A x=A^T b$42 as $A^T A x=A^T b$43 increases from $A^T A x=A^T b$44 to $A^T A x=A^T b$45, while unpreconditioned CGNR grows from $A^T A x=A^T b$46 to $A^T A x=A^T b$47. The per-iteration cost is $A^T A x=A^T b$48 (Zhao et al., 2015).

For all-at-once block Toeplitz systems from evolutionary PDEs, Hon et al. symmetrize the problem and solve the normal equations

$A^T A x=A^T b$49

with a discrete-sine-transform-diagonalizable SPD preconditioner derived from the symbol $A^T A x=A^T b$50. Their analysis proves low-rank structure in $A^T A x=A^T b$51, from which they obtain an eigenvalue clustering result for $A^T A x=A^T b$52 and, in the CGNE setting, the statement that at least $A^T A x=A^T b$53 eigenvalues are exactly $A^T A x=A^T b$54. In exact arithmetic, PCGNE therefore terminates in at most $A^T A x=A^T b$55 steps, and in practice the number of iterations is independent of the number of time steps $A^T A x=A^T b$56 (Hon et al., 2022).

For large sparse least-squares problems, Al Daas and Grigori develop two-level Schwarz preconditioners for the normal equations matrix $A^T A x=A^T b$57. The local splittings are algebraic SPSD splittings, the coarse space is built from local generalized eigenproblems, and the main condition-number estimate is

$A^T A x=A^T b$58

The upper bound is independent of the number of subdomains $A^T A x=A^T b$59, adjustable by the threshold parameter $A^T A x=A^T b$60, and implemented on top of PETSc using only $A^T A x=A^T b$61 lines of Fortran, C, or Python code. In a strong-scaling test on an $A^T A x=A^T b$62 million $A^T A x=A^T b$63 million problem, setup and solve times scaled well from $A^T A x=A^T b$64 to $A^T A x=A^T b$65 subdomains, and a single global choice $A^T A x=A^T b$66 worked well across all reported problems (Daas et al., 2021).

6. Sampling, sparsification, and optimization-derived systems

Sampling-based preconditioners exploit the additive structure of $A^T A x=A^T b$67. Chen and Wu sample rows with probabilities

$A^T A x=A^T b$68

which gives an unbiased estimator of $A^T A x=A^T b$69 and minimizes the Frobenius-norm variance among such choices. From the sampled matrix $A^T A x=A^T b$70, they build $A^T A x=A^T b$71 and define an SPD preconditioner through $A^T A x=A^T b$72 symmetric Gauss–Seidel sweeps approximating solves with $A^T A x=A^T b$73. The resulting PCG+RS method preserves sparsity because $A^T A x=A^T b$74 merely copies and rescales rows of $A^T A x=A^T b$75. In the reported sparse sprand example with $A^T A x=A^T b$76, $A^T A x=A^T b$77, iteration counts drop from $A^T A x=A^T b$78 to $A^T A x=A^T b$79 when $A^T A x=A^T b$80, and from $A^T A x=A^T b$81 (fail) to $A^T A x=A^T b$82 when $A^T A x=A^T b$83 (Chen et al., 2018).

In interior-point and proximal methods of multipliers for linear and convex quadratic programming, the Newton step produces a regularized normal-equations matrix

$A^T A x=A^T b$84

The preconditioning strategy of Frangioni and coauthors sparsifies the diagonal weighting by thresholding small entries and forms

$A^T A x=A^T b$85

Their spectral analysis shows that exactly $A^T A x=A^T b$86 eigenvalues of $A^T A x=A^T b$87 are equal to $A^T A x=A^T b$88, while the remaining eigenvalues satisfy

$A^T A x=A^T b$89

On the Netlib LP test set, IP–PMM+PCG on the normal equations solves $A^T A x=A^T b$90 of the $A^T A x=A^T b$91 reported problems in $A^T A x=A^T b$92 seconds, using $A^T A x=A^T b$93 IP steps and $A^T A x=A^T b$94 CG iterations (Bergamaschi et al., 2019).

A related general-purpose preconditioning framework for regularized interior-point methods proposes positive-definite normal-equation preconditioners obtained by dropping only complete columns and complete rows so that symmetry is preserved and the eigenvalues remain real. For the block-diagonal preconditioner $A^T A x=A^T b$95, at least

$A^T A x=A^T b$96

eigenvalues are exactly $A^T A x=A^T b$97, and the remaining eigenvalues lie in an explicit real interval bounded away from zero. The restriction to whole-column and whole-row sparsification is motivated precisely by the goal of avoiding complex-conjugate eigenpairs in the preconditioned matrix (Gondzio et al., 2021).

Indefinite least-squares problems provide a different perspective. Li and Meng transform the normal equations for

$A^T A x=A^T b$98

into a sparse block $A^T A x=A^T b$99 linear system and use GMRES with an accelerated block preconditioner

$\|A x-b\|_2$00

All eigenvalues of $\|A x-b\|_2$01 are real and cluster around $\|A x-b\|_2$02 as $\|A x-b\|_2$03. In the reported tests, $\|A x-b\|_2$04 requires only $\|A x-b\|_2$05 GMRES iterations across the tabulated examples, while competing preconditioners require $\|A x-b\|_2$06–$\|A x-b\|_2$07 or $\|A x-b\|_2$08 iterations depending on the problem class (Li et al., 23 May 2025).

Taken together, these results show that preconditioned normal equations form a broad numerical paradigm rather than a single algorithmic recipe. The common structure is the replacement of a difficult system by a symmetric positive-definite or otherwise structured surrogate whose spectrum can be controlled by sampling, operator-based modeling, algebraic decomposition, or optimization-specific sparsification. The principal qualification is numerical: fast convergence of the Krylov iteration does not automatically imply full attainable accuracy, and the most recent analyses treat iterative refinement, restarting, residual dependence, and mixed precision as intrinsic parts of the method rather than as secondary implementation details (Epperly et al., 25 Feb 2025, Garrison et al., 17 Mar 2026).