---
title: Preconditioned Normal Equations
url: https://www.emergentmind.com/topics/preconditioned-normal-equations
type: topic
---

# Preconditioned Normal Equations

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 [2502.17767][2603.16644][2502.17626].

## 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 [2507.18466][2502.17626].

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_p^T A\,x=A_p^T b\) | Half-preconditioned; generally nonsymmetric |
| SNE | \(R^T R\,x=A^T b\) | Seminormal equations from thin QR |
| NNE | \(B^T A\,x=B^T b\) | Uses any full-rank \(B\) with the same column space as \(A\) |

These formulations are analyzed for perturbation behavior in mixed precision and under randomized sketch-based preconditioning [2603.16644].

For square nonsymmetric systems, right preconditioning is commonly written with \(y=M_2 x\), so that
\[
A^T A\,M_2^{-1}y=A^T b.
\]
With both left and right preconditioners one may write
\[
M_1^{-1}(A^T A)M_2^{-1}y=M_1^{-1}A^T b,
\]
and in the notation of Epperly, Greenbaum, and Nakatsukasa one typically takes \(M_1=P^T\), \(M_2=P\), choosing \(P\) to control \(\mathrm{cond}(A P^{-1})\). Their basic assumptions include
\[
\mathrm{cond}(A P^{-1})\le C,\qquad \|A P^{-1}\|>1,
\]
hence \(1\le \|A P^{-1}\|\le \mathrm{cond}(A P^{-1})\) and \(\|P^{-1}\|\le \mathrm{cond}(A P^{-1})/\|A\|\) [2502.17767].

## 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
\[
\min_z \|b-A P^{-1} z\|_2,
\]
whose normal equations are
\[
(P^{-T}A^T A P^{-1})\,y=P^{-T}A^T b,\qquad x=P^{-1}y.
\]
The matrix-vector operations are therefore
\[
v\mapsto A(P^{-1}v),\qquad w\mapsto A^T(P^{-T}w),
\]
so each Krylov step uses four primitive solves or multiplications: \(A\), \(A^T\), \(P^{-1}\), and \(P^{-T}\). The Lanczos bidiagonalization recurrences are the standard LSQR recurrences with \(A\) replaced by the composite operator \(A P^{-1}\), and at iteration \(k\) the method forms a \(k\times k\) bidiagonal matrix \(B_k\) whose singular-value approximation yields \(y_k\), hence \(x_k=P^{-1}y_k\) [2502.17767].

For normal equations arising from PDE discretization, CGNE and LSQR are often described as residual-minimizing methods in the original \(2\)-norm:
\[
x_k=\arg\min_{x\in\mathcal K_k(B,A^T b)} \|A x-b\|_2,\qquad B=A^T A,
\]
with the standard CG estimate
\[
\frac{\|x_k-x_*\|_B}{\|x_0-x_*\|_B}
\le
2\left(\frac{\kappa(B)^{1/2}-1}{\kappa(B)^{1/2}+1}\right)^k
=
2\left(\frac{\kappa(A)-1}{\kappa(A)+1}\right)^k.
\]
This makes preconditioning the normal operator, rather than \(A\) directly, a natural design target when one can construct an effective SPD approximation to \(A^T A\) or \(A^T T A\) [2502.17626].

A recurrent comparison is with GMRES on the right-preconditioned nonsymmetric system
\[
A P^{-1} y=b,\qquad x=P^{-1}y.
\]
For GMRES, the Saad–Schultz bound states that after \(k\) steps
\[
\|b-A P^{-1}y_k\|
\le
\min_{p(0)=1,\deg p\le k}\|p(A P^{-1})\|\,\|b\|,
\]
so \(p\) must approximate \(1/z\) on the spectrum of \(A P^{-1}\). The contrast emphasized in the stability analysis is that LSQR-based methods only require \(\mathrm{cond}(A P^{-1})=O(1)\), whereas GMRES additionally benefits from a tightly clustered spectrum in the complex plane, not wrapping around the origin, and from well-conditioned eigenvectors [2502.17767].

## 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 \(\mathrm{cond}(A P^{-1})\ll \mathrm{cond}(A)\), finite-precision LSQR can stagnate far above the backward-stable level
\[
\|b-A x_k\|\le O(\|A\|\,\|x_k\|\,u).
\]
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 \(u\) [2502.17767].

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
\[
r_i=b-Ax_i,\qquad c_i=P^{-1}A^T r_i,
\]
approximately solves
\[
(P^{-T}A^T A P^{-1})\,\Delta y_i=c_i
\]
by LSQR, and updates
\[
x_{i+1}=x_i+P^{-1}\Delta y_i.
\]
The second is restarting preconditioned LSQR when the backward-error estimate
\[
\bar b_i=\frac{\|b-Ax_i\|}{\|A\|\,\|x_i\|}
\]
fails to decrease sufficiently over a fixed number of iterations. In the cited backward-error analysis, if \(\|A P^{-1}\|u\ll 1\), each refinement step satisfies
\[
\|b-Ax_{i+1}\|\le \eta\,\|b-Ax_i\|+O(\|A\|\,\|x_i\|\,u),
\]
with \(\eta=O(\mathrm{cond}(A P^{-1})u)<1\). The resulting theorem states that if each inner LSQR solve takes
\[
q=O(\log(1/(\mathrm{cond}(A P^{-1})u)))
\]
iterations, then after
\[
t=O\!\left(\frac{\log(1/u)}{\log(1/(\mathrm{cond}(A P^{-1})u))}\right)
\]
refinement steps the overall procedure is backward stable, with only \(O(\log(1/u))\) calls to \(A\), \(A^T\), \(P^{-1}\), and \(P^{-T}\) [2502.17767].

The same phenomenon appears for preconditioned CG on SPD normal equations. Periodic iterative refinement—apply PCG for \(q\) steps, recompute the residual in high precision, and repeat—restores
\[
\|b-Ax_t\|\le O(\|A\|\,\|x_t\|\,u)
\]
at essentially the same cost. This extends the attainable-accuracy picture beyond LSQR to the broader family of Krylov solvers on normal equations [2502.17767].

A common source of confusion is the relation between convergence and accuracy. Small \(\mathrm{cond}(A P^{-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 \(\mathrm{cond}(P^{-1}A)\) is small, “say \(\lesssim 10\),” in
\[
O(\mathrm{cond}(P^{-1}A)^{1/2}\,\log_{10}(\|A\|/u))
\]
iterations, but the cited work states explicitly that no proof is yet available and that a complete analysis remains open [2502.17767].

## 4. Randomized and mixed-precision preconditioners

A prominent recent direction computes the preconditioner from a randomized sketch. For full-column-rank \(A\), one forms a “fast-JL” sketch
\[
\Omega=SFD\in\mathbb R^{d\times m},
\]
where \(D\) is a random sign diagonal, \(F\) is a real or complex Fourier or Hadamard transform, and \(S\) samples \(d\) rows from the identity. In lower precision \(u_p\), one computes the thin QR factorization
\[
\Omega A=Q_s R_s,
\]
promotes \(R_s\) to working precision, and sets \(A_p=A R_s^{-1}\). If
\[
d \ge 2\,m\,\mu\,(1+\tfrac{\epsilon}{3})\,\frac{\ln(n/\delta)}{\epsilon^2},
\]
then with probability \(1-\delta\),
\[
\kappa(A_p)\le \sqrt{\frac{1+\epsilon}{1-\epsilon}}.
\]
This construction underlies both symmetrically preconditioned normal equations and half-preconditioned normal equations in mixed precision [2603.16644].

The associated perturbation bounds are formulated in terms of the relative normal-equations residual
\[
\rho=\frac{\|A\hat x-b\|}{\|A\|\,\|\hat x\|}.
\]
For PNE, the computed solution satisfies a bound scaled by \(\kappa(R_s)\kappa(A_p)\approx \kappa(A)\), growing linearly in \(\epsilon_A\), with the max-term equal to \(1\) when \(\rho\ll 1\) and growing like \(\rho\) when \(\rho\gg 1\). For HPNE, the qualitative dependence is the same: weak on preconditioner error, but growing with large \(\rho\). The seminormal equations inherit the \(\kappa(A)^2\) 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 [2603.16644].

The mixed-precision implementation includes an automatic precision-selection rule. In single precision one estimates \(\kappa(A)\approx 10^{\kappa_0}\) via a fast \(1\)-norm condition estimator; the sketch precision is then chosen as half if \(\kappa_0<4\), single if \(4\le \kappa_0<8\), and double if \(\kappa_0\ge 8\) or overflow. The PNE or HPNE solve is performed in double precision. In experiments with \(m=6000\), \(n=100\), \(\kappa(A)=10^4\), the new bounds track the observed \(\|x_*-\hat x\|/\|\hat x\|\) versus \(\rho\), and for \(\rho\gtrsim 10^{-4}\), PNE and HPNE attain the same accuracy as MATLAB’s backslash. In mixed precision with \(m=6000\), \(n=100\) or \(1000\), \(\kappa(A)=10^8\), and a single-precision preconditioner, the new bounds remain informative and for \(\rho\gtrsim 10^{-6}\), mixed-precision PNE and HPNE again match mldivide’s accuracy. On NVIDIA H100 GPUs, mixed precision delivers up to a \(10\)–\(20\%\) speedup over double-PNE and often out-runs the QR solver for \(n\) up to a few thousand [2603.16644].

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 \(\kappa(A_p)\approx 1\) and \(\kappa(R_s)\approx \kappa(A)\). Numerical experiments reported there extend to matrices with \(\kappa(A)=10^{12}\), with observed relative errors tracking the QR-based solution and plateauing at \(\kappa(A)\epsilon\) [2507.18466].

## 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 \(\mathcal L:V\to V'\), the “normal PDE” is based on the composition \(\mathcal L^*\mathcal L\), which is self-adjoint and elliptic. After discretization, one solves
\[
A^T T A\,x=A^T T b
\]
with a preconditioner
\[
G=P^T T P\approx A^T T A,
\]
where \(T\) is a Riesz map associated with the chosen inner product. The spectral rationale is that \(G^{-1}A^T T A\) is spectrally equivalent to the Gram matrix of \(A P^{-1}\) in the \(T\)-inner product, so its eigenvalues are the squares of the \(T\)-singular values of \(A P^{-1}\). An ideal normal preconditioner is therefore any \(P\) for which \(A P^{-1}\) is \(T\)-orthogonal [2502.17626].

For the space-time fractional advection-diffusion equation, Zhao, Jin, and Lin derive a nonsymmetric Toeplitz-like matrix \(I+A\) and form the normal equations
\[
(I+A)^T(I+A)\,u=(I+A)^T b.
\]
Their preconditioner is \(M=P_\ell^T P_\ell\), where \(P_\ell=I+A_\ell\) is obtained by band truncation of the Toeplitz blocks. Because \(\|A-A_\ell\|/\|A\|=O(\ell^{-\beta})\), the eigenvalues of the preconditioned normal matrix are clustered about \(1\), with
\[
\kappa(M^{-1}N)\le \frac{1+\epsilon}{1-\epsilon}.
\]
The reported PCGNR iteration counts for Example 4.1, with bandwidth \(\ell=8\), range from \(3.1\) to \(5.1\) as \(m=n\) increases from \(16\) to \(256\), while unpreconditioned CGNR grows from \(12.4\) to \(1060.2\). The per-iteration cost is \(O(N\log N+N\ell)\) [1510.05089].

For all-at-once block Toeplitz systems from evolutionary PDEs, Hon et al. symmetrize the problem and solve the normal equations
\[
\T_H^T \T_H\,\mathbf u=\T_H^T\mathbf f_H
\]
with a discrete-sine-transform-diagonalizable SPD preconditioner derived from the symbol \(|g_\theta|^2\). Their analysis proves low-rank structure in \(M^K-P^K\), from which they obtain an eigenvalue clustering result for \(P^{-1}M\) and, in the CGNE setting, the statement that at least \(mn-m\) eigenvalues are exactly \(1\). In exact arithmetic, PCGNE therefore terminates in at most \(m+1\) steps, and in practice the number of iterations is independent of the number of time steps \(n\) [2201.10062].

For large sparse least-squares problems, Al Daas and Grigori develop two-level Schwarz preconditioners for the normal equations matrix \(C=A^T A\). The local splittings are algebraic SPSD splittings, the coarse space is built from local generalized eigenproblems, and the main condition-number estimate is
\[
\kappa(M_{\rm add}^{-1}C)\le (k_c+1)\Bigl(2+(2k_c+1)\frac{k_m}{\tau}\Bigr).
\]
The upper bound is independent of the number of subdomains \(N\), adjustable by the threshold parameter \(\tau\), and implemented on top of PETSc using only \(150\) lines of Fortran, C, or Python code. In a strong-scaling test on an \(8\) million \(\times\,7.5\) million problem, setup and solve times scaled well from \(16\) to \(512\) subdomains, and a single global choice \(\tau=0.6\) worked well across all reported problems [2107.09006].

## 6. Sampling, sparsification, and optimization-derived systems

Sampling-based preconditioners exploit the additive structure of \(A^T A=\sum_i a_i^T a_i\). Chen and Wu sample rows with probabilities
\[
p_i=\frac{\|a_i\|_2^2}{\|A\|_F^2},
\]
which gives an unbiased estimator of \(A^T A\) and minimizes the Frobenius-norm variance among such choices. From the sampled matrix \(A_s=S A\), they build \(M=A_s^T A_s\) and define an SPD preconditioner through \(t\) symmetric Gauss–Seidel sweeps approximating solves with \(M\). The resulting PCG+RS method preserves sparsity because \(S\) merely copies and rescales rows of \(A\). In the reported sparse sprand example with \(m=90\,000\), \(n=300\), iteration counts drop from \(98\) to \(21\) when \(\kappa(A^T A)=3.9\times 10^3\), and from \(296\) (fail) to \(70\) when \(\kappa(A^T A)=7.4\times 10^5\) [1806.02968].

In interior-point and proximal methods of multipliers for linear and convex quadratic programming, the Newton step produces a regularized normal-equations matrix
\[
M_{NE,k}=A(Q+\Theta_k^{-1}+\rho_k I_n)^{-1}A^T+\delta_k I_m.
\]
The preconditioning strategy of Frangioni and coauthors sparsifies the diagonal weighting by thresholding small entries and forms
\[
P_{NE,k}=A E_k A^T+\delta_k I_m.
\]
Their spectral analysis shows that exactly \(m-\mathrm{rank}(A)\) eigenvalues of \(P_{NE,k}^{-1}\tilde M_{NE,k}\) are equal to \(1\), while the remaining eigenvalues satisfy
\[
1\le \lambda \le 1+\frac{C_{E,k}\mu_k}{\delta_k}\sigma_{\max}^2(A).
\]
On the Netlib LP test set, IP–PMM+PCG on the normal equations solves \(100\%\) of the \(96\) reported problems in \(141.3\) seconds, using \(2\,907\) IP steps and \(101\,382\) CG iterations [1912.10064].

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 \(P_{NE,(k_c,k_r)}\), at least
\[
\max\{m-(2k_r+k_c),0\}
\]
eigenvalues are exactly \(1\), 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 [2107.06822].

Indefinite least-squares problems provide a different perspective. Li and Meng transform the normal equations for
\[
\min_{x\in\mathbb R^n}(b-Ax)^T J (b-Ax),\qquad J=\mathrm{diag}(I_p,-I_q),
\]
into a sparse block \(3\times 3\) linear system and use GMRES with an accelerated block preconditioner
\[
\mathcal P=
\begin{pmatrix}
I & A_1 & 0\\
A_1^T & \alpha I & -A_2^T\\
0 & A_2 & I
\end{pmatrix}.
\]
All eigenvalues of \(\mathcal P^{-1}\mathcal A\) are real and cluster around \(1\) as \(\alpha\to 0_+\). In the reported tests, \(\mathcal P\) requires only \(2\) GMRES iterations across the tabulated examples, while competing preconditioners require \(3\)–\(12\) or \(4\) iterations depending on the problem class [2505.17504].

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 [2502.17767][2603.16644].

Source: https://www.emergentmind.com/topics/preconditioned-normal-equations