Why Solve Mixed-Model Systems Iteratively?
In a genome-wide association analysis, the number of individuals can be large enough that repeatedly factoring a dense covariance matrix is expensive in both time and memory. Iterative methods offer another route: instead of explicitly calculating a matrix inverse, they use repeated matrix-vector products to approximate solutions to linear systems. Their value depends on more than the iteration count. The system must suit the solver, the matrix-vector products must be affordable, and the stopping rule must provide adequate accuracy for the intended analysis.
For a linear mixed model, write \(y=X\beta+u+\varepsilon\), where \(y\) is the vector of phenotypes, \(X\) is the fixed-effect design matrix, and \(u+\varepsilon\) has covariance matrix \(V\). In a common model, \(V=\sigma_g^2K+\sigma_e^2I\), with relationship matrix \(K\), identity matrix \(I\), and variance components \(\sigma_g^2\) and \(\sigma_e^2\). Given fixed variance components, generalized least squares uses \(V^{-1}y\) and \(V^{-1}X\). These quantities can be obtained by solving systems such as \(Vz=y\) and \(VW=X\), without forming \(V^{-1}\).
For the covariance system \(A=V\), standard conjugate gradients require \(V\) to be symmetric positive definite (SPD): \(V=V^\mathsf{T}\), and \(q^\mathsf{T}Vq>0\) for every nonzero vector \(q\). A symmetric relationship matrix with nonnegative eigenvalues and \(\sigma_e^2>0\) gives this property, since \(q^\mathsf{T}Vq=\sigma_g^2q^\mathsf{T}Kq+\sigma_e^2q^\mathsf{T}q>0\) when \(\sigma_g^2\geq0\). If the assumptions fail, ordinary conjugate gradients are not justified; a solver designed for the actual matrix structure is needed. In practice, numerical roundoff can also cause trouble if the system is nearly singular.
Conjugate Gradients and Stopping
Conjugate gradients (CG) begins with an initial guess \(\hat{x}_0\), forms its residual \(r_0=b-A\hat{x}_0\), and updates the solution along carefully chosen directions. For an SPD matrix, each new direction is conjugate to the earlier directions under the \(A\)-inner product: \(p_i^\mathsf{T}Ap_j=0\) for \(i\ne j\). This structure makes CG more efficient than repeatedly moving in the direction of the steepest local decrease. In exact arithmetic, it finds the solution in at most \(n\) iterations for an \(n\)-dimensional system. In floating-point arithmetic, convergence can take longer and is usually judged by a tolerance.
For \(k\geq0\), the standard CG direction update is \(p_{k+1}=r_{k+1}+\beta_kp_k\), with \(\beta_k=r_{k+1}^\mathsf{T}r_{k+1}/(r_k^\mathsf{T}r_k)\). These formulas assume the denominators are nonzero; if the residual is zero, the current iterate is already an exact solution. CG minimizes the error over an expanding sequence of Krylov subspaces, the spaces generated by \(r_0,Ar_0,A^2r_0,\ldots\). The matrix itself need not be stored if the product \(Aq\) can be computed when needed.
A practical stopping rule checks the relative residual \(\|r_k\|_2/\|b\|_2\), when \(b\ne0\), against a chosen tolerance. An absolute residual threshold can also be useful when the scale of \(b\) is meaningful. The tolerance should be set in light of the downstream calculation: a solve used repeatedly inside a statistical procedure may need tighter control than a preliminary computation. A small residual is evidence of a good fit to the equations, but does not always imply a small error in the solution.
Worked Example: Two Conjugate-Gradient Steps
Consider the hypothetical SPD system \(A x=b\), with \(A=\begin{pmatrix}4&1\\1&3\end{pmatrix}\), \(b=(1,2)^\mathsf{T}\), and initial guess \(\hat{x}_0=(0,0)^\mathsf{T}\). The matrix is symmetric, and its leading diagonal entries and determinant are positive (\(4>0\), \(12-1=11>0\)), so it is SPD. No genomic measurements are being represented by these illustrative numbers.
Initially, \(r_0=b-A\hat{x}_0=(1,2)^\mathsf{T}\) and \(p_0=r_0\). Thus \(r_0^\mathsf{T}r_0=5\) and \(Ap_0=(6,7)^\mathsf{T}\), so \(p_0^\mathsf{T}Ap_0=20\). The first step size is \(\alpha_0=5/20=0.25\), giving \(\hat{x}_1=(0.25,0.50)^\mathsf{T}\). Its residual is \(r_1=(1,2)^\mathsf{T}-0.25(6,7)^\mathsf{T}=(-0.50,0.25)^\mathsf{T}\).
Now \(r_1^\mathsf{T}r_1=0.3125\), so \(\beta_0=0.3125/5=0.0625\), and \(p_1=r_1+\beta_0p_0=(-0.4375,0.375)^\mathsf{T}\). Multiplication gives \(Ap_1=(-1.375,0.6875)^\mathsf{T}\), hence \(p_1^\mathsf{T}Ap_1=0.859375\) and \(\alpha_1=0.3125/0.859375=4/11\). Therefore \(\hat{x}_2=(0.25,0.50)^\mathsf{T}+(4/11)(-0.4375,0.375)^\mathsf{T}=(1/11,7/11)^\mathsf{T}\). Checking directly, \(A\hat{x}_2=(1,2)^\mathsf{T}=b\), so the residual is zero. In floating-point computation, the same check would typically give a small, rather than exactly zero, residual.
Residuals, Conditioning, and What They Do Not Guarantee
A useful accuracy result explains why a residual alone cannot fully characterize solution error. Let \(A\) be nonsingular, \(Ax=b\), and \(r=b-A\hat{x}\). When \(b\ne0\) (equivalently \(x\ne0\)), the relative error is bounded by the condition number of \(A\) times the relative residual:
To justify this, observe that \(r=A(x-\hat{x})\), so \(x-\hat{x}=A^{-1}r\) and \(\|x-\hat{x}\|_2\leq\|A^{-1}\|_2\|r\|_2\). Also \(b=Ax\), which gives \(\|b\|_2\leq\|A\|_2\|x\|_2\), or \(1/\|x\|_2\leq\|A\|_2/\|b\|_2\). Multiplying the two inequalities yields the stated bound. The nonzero condition matters: if \(b=0\), the exact solution is \(x=0\), so the relative-error denominators vanish. In that case, use absolute error or residual criteria instead.
The bound is a worst-case guarantee and can be loose. Still, it identifies a central diagnostic: when \(\kappa_2(A)\) is large, even a small relative residual may not ensure high relative accuracy in the solution. For SPD \(A\), CG also has an error bound in the \(A\)-norm, \(\|e_k\|_A\leq2[(\sqrt{\kappa_2(A)}-1)/(\sqrt{\kappa_2(A)}+1)]^k\|e_0\|_A\), where \(e_k=x-\hat{x}_k\) and \(\|v\|_A=\sqrt{v^\mathsf{T}Av}\). This motivates preconditioning: changing the system in a way that improves its conditioning while preserving a valid solution.
Worked Example: A Small Residual Can Give a Loose Error Bound
Take the hypothetical system \(A=\operatorname{diag}(1,100)\), \(b=(1,100)^\mathsf{T}\), and approximate solution \(\hat{x}=(1,0.99)^\mathsf{T}\). The exact solution is \(x=(1,1)^\mathsf{T}\). The matrix is SPD and has \(\kappa_2(A)=100\). Direct substitution gives \(A\hat{x}=(1,99)^\mathsf{T}\), so \(r=(0,1)^\mathsf{T}\).
Here \(\|r\|_2=1\) and \(\|b\|_2=\sqrt{10001}\approx100.005\), so the relative residual is about \(0.0099995\). The actual relative solution error is \(\|x-\hat{x}\|_2/\|x\|_2=0.01/\sqrt{2}\approx0.007071\). The bound gives at most \(100(0.0099995)\approx0.99995\). It is valid but very conservative: passing a relative-residual threshold near 0.01 does not, by itself, certify a correspondingly accurate solution for every system.
Preconditioning and Computational Diagnostics
A preconditioner is an SPD matrix \(M\) that is comparatively easy to apply and approximates \(A\) in a way that improves the effective conditioning. Preconditioned CG (PCG) is equivalent to applying CG to the symmetrically transformed system \(M^{-1/2}AM^{-1/2}z=M^{-1/2}b\), then recovering \(x=M^{-1/2}z\). The relevant condition number is that of \(M^{-1/2}AM^{-1/2}\), not simply the condition number of \(A\). A useful preconditioner reduces iteration work enough to compensate for its own construction and application costs. It must be SPD for standard PCG theory to apply.
Worked Example: How an Idealized Preconditioner Changes the System
For illustration, let \(A=\operatorname{diag}(1,100)\), \(b=(1,100)^\mathsf{T}\), and choose \(M=\operatorname{diag}(1,100)\). Both matrices are SPD. This is an idealized choice: it makes \(M=A\), so it demonstrates the algebra but is not generally a cheap or useful choice when the full matrix is difficult to handle.
The transformed matrix is \(M^{-1/2}AM^{-1/2}=I\), and the transformed right-hand side is \(M^{-1/2}b=(1,10)^\mathsf{T}\). The transformed solution is \(z=(1,10)^\mathsf{T}\), and \(x=M^{-1/2}z=(1,1)^\mathsf{T}\). The original residual is \(b-Ax=(0,0)^\mathsf{T}\). Here the condition number falls from 100 to 1, so CG solves the transformed problem in one iteration in exact arithmetic. Real preconditioners rarely remove all conditioning, and their benefit must be measured rather than assumed.
For GWAS mixed models, these solves may be needed for vectors such as \(y\), columns of \(X\), or variant data, depending on the association procedure. Check convergence for the original equations \(Az=b\), not only for a preconditioned residual: a transformed residual can look small while the original residual is not adequate. Useful checks include the relative and absolute residual, the number of iterations, whether the residual decreases consistently, and whether independently recomputed \(b-A\hat{x}\) agrees with the solver's reported residual. Recursive residual updates can drift through rounding, so occasional direct recomputation is prudent. In a statistical workflow, also check that the resulting GLS quantities and association statistics are stable under a tighter tolerance.
Confirm symmetry and positive definiteness for CG or PCG, and verify that the variance-component values and covariance operator correspond to the intended fitted model.
Record the tolerance, iteration count, and residual history. Investigate stalled or increasing residuals rather than treating a maximum iteration count as convergence.
Recompute the original-system residual and compare results at a tighter tolerance. Count matrix-vector products and include preconditioner setup and application costs when judging speed.
Choosing an Approach and Recognizing Failure
Direct factorization and iterative solution make different trade-offs. For a moderate dense SPD matrix, Cholesky factorization can be reliable and allows repeated solves with different right-hand sides, though its memory and computational costs grow rapidly with matrix size. The earlier material on eigenvalue decomposition describes an exact LMM scan strategy that can be efficient when its decomposition cost is justified and reused. Iterative methods avoid explicitly constructing a factorization, but may require many matrix-vector products, and each solve can have approximation error. The right choice depends on matrix size, structure, number of right-hand sides, reuse, and required accuracy.
A critical limitation is that iteration does not automatically make a costly matrix-vector product cheap. If applying \(V\) requires dense \(n\times n\) multiplication, one CG iteration still takes order \(n^2\) operations. An iterative approach is attractive when the covariance action can be computed efficiently and the system converges in a manageable number of iterations. The eigenvalue decomposition approach and iterative solvers also differ in what they provide: a decomposition can support repeated transformed calculations, while an iterative solve targets particular right-hand sides.
- Stopping on iteration count alone. A fixed number of iterations is not proof of convergence; inspect the original-system residual.
- Confusing residual with solution error. A small residual can coexist with substantial error when the system is ill-conditioned; compare tolerances and consider conditioning.
- Using CG on an unsuitable matrix. Nonsymmetry or failure of positive definiteness invalidates standard CG guarantees.
- Ignoring numerical or statistical consequences. If a tighter solve tolerance changes downstream estimates or test statistics materially, the original tolerance was not adequate for that analysis.
Check Your Understanding
Answer each question, then compare with the solutions below.
- What property of a covariance matrix is required by standard conjugate gradients, and what model condition can help ensure it for \(V=\sigma_g^2K+\sigma_e^2I\)?
- For an approximate solution \(\hat{x}\), why is a small residual not always a guarantee of a small relative error?
- Suppose \(b=0\) and \(A\) is nonsingular. Why should the relative-error bound above not be applied?
- What should be checked if a preconditioned residual is small but the original-system residual is not?
- Why might CG remain computationally expensive even if it converges in relatively few iterations?
Solutions. (1) Standard CG requires a symmetric positive-definite matrix. If \(K\) is symmetric positive semidefinite, \(\sigma_g^2\geq0\), and \(\sigma_e^2>0\), then \(V\) is positive definite. (2) The relative error is bounded in terms of the relative residual multiplied by \(\kappa_2(A)\); a large condition number can make that bound loose. (3) The exact solution is \(x=0\), so the relative-error denominators \(\|x\|_2\) and \(\|b\|_2\) are zero; use absolute criteria instead. (4) Recompute and assess \(b-A\hat{x}\) in the original system, and do not declare convergence until that residual meets the chosen criterion. (5) Each iteration requires a matrix-vector product and possibly preconditioner application; if those operations are expensive, even a modest iteration count can be costly.