6.9 Efficient linear solvers🔗⧉

A quadratic restoration model balances agreement with the measurements against a penalty on the image:

Here $D$ is the identity for a small-value penalty, or a spatial difference operator for a smoothness penalty. Setting the gradient to zero gives the normal equations

It is enormous. For an image of $N = W\cdot H\cdot C$ values the operator $A^\top A$ is $N\times N$, tens of millions on a side, and yet solving this is routine on a laptop. This chapter is how. The trick is not a faster inverse but refusing to form the matrix at all, and the same machinery reaches well past deblurring: gradient-domain HDR, Chapter 10.8, Poisson image editing, colorization, and matting are each this same solve with a different operator $A$.

6.9.1 The problem, and the one trick🔗⧉

Start from the unregularized objective $f(x) = \tfrac12\lVert Ax - y\rVert^2$ (plus a regularizer, held aside for a paragraph). To locate its minimum, set the gradient to zero. The gradient of $\tfrac12\lVert Ax-y\rVert^2$ is

Setting the gradient to zero gives $A^\top Ax=A^\top y$. At an unregularized optimum, the data residual is orthogonal to the column space of $A$. Adding the quadratic penalty contributes $\lambda D^\top Dx$, yielding $(A^\top A+\lambda D^\top D)x=A^\top y$.

For a large system, implement the normal operator through its components: $Qv=A^\top(Av)+\lambda D^\top(Dv)$. There is no need to store $Q$. For the special case of periodic convolution, the data contribution is

That leaves one character to explain, $A^\top$. For periodic boundaries, the transpose of a blur is another convolution. The transpose of "convolve by $k$" is "convolve by the flipped kernel" $\tilde k(\mathbf u) = k(-\mathbf u)$, the kernel mirrored through its center, which turns convolution into correlation. So the two operations any solver ever asks for are

for periodic boundaries, each one a call to the convolution routine: apply_blur(x) = conv(x, k) for $A$, apply_blur_T(y) = conv(y, flip(k)) for $A^\top$. No matrix is formed, stored, or inverted anywhere. This is what matrix-free means, and it is the reason a problem with tens of millions of unknowns fits in memory and runs on a laptop.

For periodic convolution, the adjoint uses the flipped kernel (conjugate-flipped for complex data). Finite boundaries require more care. If $A=CKE$ extends the input with $E$, convolves with $K$, then crops with $C$, its transpose is $A^\top=E^\top K^\top C^\top$. The transpose of cropping inserts zeros; the transpose of reflection or replication accumulates repeated contributions. Reusing the forward boundary rule with a flipped kernel need not give the adjoint.

For a fixed warp or downsampler, scatter each output value back to the input locations with exactly the weights used by the forward gather. This is back-projection, not inverse warping. The dot-product test below checks the complete implementation, including boundaries.

For example, let a forward resampler form one output from two input pixels with weights one-quarter and three-quarters. Its transpose sends an output residual back to those same two input positions with those same weights. If several outputs used one input, their back-projected contributions add there. Running the geometric warp backward usually performs a new interpolation instead, and therefore need not implement the transpose. This distinction matters because an iterative solver uses the adjoint to turn output errors into changes of the input estimate. A plausible-looking back-projection is not enough; check the dot-product identity.

6.9.2 Direct methods, and why they wall out🔗⧉

A direct solver factors a matrix and back-substitutes. Cholesky applies when the system matrix is symmetric positive-definite; QR is useful for least squares and avoids explicitly forming normal equations. Sparse factorizations reorder variables to reduce fill-in. These methods can be effective for moderate images or when the same system is solved repeatedly with different measurements.

Dense storage is prohibitive at image scale. Sparse factors may also require much more memory than the original operator because elimination creates new nonzeros, called fill-in. Whether a direct solve is practical depends on dimensions, sparsity pattern, ordering, and reuse. Large systems often favor matrix-free iteration.

There is one exception, and it is so fast that entire workflows are built around it when it applies. When the blur is shift-invariant and the image is treated as wrapping around (periodic boundaries), $A$ is circulant, and a circulant operator is diagonalized by the Fourier transform L6.8show. A Laplacian regularizer diagonalizes in the same basis. Once everything is diagonal there is no coupled system left to solve at all: we transform, divide frequency by frequency, and transform back. Written out, the regularized solution is

Here $D$ must also be diagonalized by the same Fourier basis; for a multi-component gradient penalty, sum the squared responses of its components. A denominator can still vanish if both the data operator and penalty annihilate the same frequency. Periodic boundaries give a circular problem. Some symmetric Neumann-boundary operators instead admit a cosine transform; Neumann boundaries are not periodic. Padding an inverse problem or windowing its observations does not generally make its solution identical to the original finite-boundary problem. Fast methods such as Krishnan and Fergus use spectral solves inside a broader optimization.

For a general operator or boundary model, apply the normal-equation operator without forming it:

We solve $Qx=b$. This notation keeps the regularizer present in every iterative update.

6.9.3 Basic iterative solvers🔗⧉

Gradient descent improves an initial estimate by stepping opposite the gradient. For the regularized quadratic above, the gradient is $Qx-b$, so the update is

where $\eta$ is the step size. Applying $Q$ includes both the forward/adjoint pair and the penalty operators:

x = initial_estimate
for t in 0 .. max_iterations-1:
    g = Q(x) - b
    if norm(g) <= atol + rtol * norm(b): break
    x = x - eta * g

For positive-definite $Q$, a fixed step requires $0<\eta<2/\lambda_{\max}(Q)$. Use an estimate of the largest eigenvalue or a line search; an arbitrary step can diverge.

The classical stationary iterations (Jacobi, Gauss-Seidel, Richardson) are variations on this same theme, differing in how they precondition or order the update. Their convergence depends on the operator and update parameters. They can be slow on an elongated quadratic. The least-squares bowl is badly elongated: the blur squashes high frequencies hard and low frequencies barely, so the surface is far steeper in some directions than others, a long narrow trough rather than a round basin. Plain gradient descent, forever stepping straight downhill, zig-zags across such a trough, taking thousands of timid steps to inch along its length (Figure 6.9.3⛶ fullscreen). How elongated the bowl is, and so how slow the solver, is measured by the condition number of $A^\top A$; the worse the conditioning, the more the ill-posedness of the problem shows up as sluggish convergence.

6.9.4 Conjugate gradient🔗⧉

Conjugate gradient solves a symmetric positive-definite system using matrix-vector products. With $\lambda>0$, $Q=A^\top A+\lambda D^\top D$ is positive-definite exactly when $\ker A\cap\ker D=\{0\}$. An identity penalty ensures this, but a gradient penalty leaves constants unpenalized: the data or another constraint must determine them.

CG constructs directions conjugate with respect to $Q$, so a line minimization does not undo previous ones in exact arithmetic. In that ideal arithmetic it terminates in at most $N$ steps; finite precision and the eigenvalue distribution determine practical convergence. The objective decreases, but the Euclidean residual norm need not decrease at every step. Shewchuk gives a geometric derivation.

x = initial_estimate
r = b - Q(x)
d = r
rho = dot(r, r)
tol = atol + rtol * norm(b)
for t in 0 .. max_iterations-1:
    if sqrt(rho) <= tol: break
    q = Q(d)
    denom = dot(d, q)
    if denom <= 0: report_failure("operator is not positive definite"); stop
    alpha = rho / denom
    x = x + alpha * d
    r = r - alpha * q
    rho_new = dot(r, r)
    d = r + (rho_new / rho) * d
    rho = rho_new
# Verify convergence with the true residual, not only the recurrence.
converged = norm(b - Q(x)) <= tol
return x, converged

LSQR and LSMR apply $A$ and $A^\top$ to solve rectangular least squares without explicitly forming $A^\top A$. For the regularized problem, use the augmented operator $[A;\sqrt\lambda D]$ and target $[y;0]$. This avoids the rounding error and storage of an explicitly formed normal matrix; it does not remove the inverse problem’s intrinsic conditioning. MINRES handles symmetric indefinite systems, and GMRES general nonsymmetric systems.

The geometry is easiest to picture with elongated elliptical contours of a quadratic objective. Steepest descent points across the narrow direction of the valley, so successive line searches can zigzag while making little progress along it. Conjugate gradient chooses directions that are orthogonal under the system matrix's quadratic form. In exact arithmetic, minimizing along one of these directions does not undo the minimizations along the earlier ones. This gives the finite-dimensional termination result for a symmetric positive-definite system. Roundoff and conditioning affect the practical iteration count, so the useful implementation still checks residuals and can benefit greatly from preconditioning.

6.9.5 Preconditioning🔗⧉

A preconditioner approximates the system in a cheaper form. Let $M\approx Q^{-1}$ denote the operation applied to a residual. In preconditioned CG, set $z=Mr$, use $r^\top z$ in place of $r^\top r$, and form search directions from $z$. The preconditioner must preserve the required symmetry and positive definiteness. Diagonal scaling, incomplete Cholesky, an approximate Fourier solve, and suitable multigrid cycles are useful choices. They reduce solver effort; they do not add information to the measurements.

The choices approximate different parts of the difficulty. Diagonal scaling rescales coordinates with very different local gains. Incomplete Cholesky approximates a factorization while limiting fill-in. A Fourier preconditioner solves a nearby shift-invariant problem cheaply even when the actual operator has spatial variation or different boundaries. A multigrid cycle approximates corrections across several spatial scales. Applying a more accurate approximation costs more, so compare total solution time rather than iteration count alone.

6.9.6 Multigrid🔗⧉

For elliptic systems such as Poisson equations, local relaxation damps oscillatory error quickly and smooth error slowly. Multigrid transfers the residual to a coarse system, solves for a correction there, prolongs that correction, and smooths again Brandt 1977. The coarse operator must be consistent with the transfer operators, for example $Q_c=RQP$ for restriction $R$ and prolongation $P$. These are solver operators, not automatically a Gaussian pyramid’s reduce and expand. A well-designed multigrid method can achieve nearly linear work for suitable problems; arbitrary deblurring systems need not have the same smoothing behavior.

Consider an error that varies gradually from one side of the image to the other. A local relaxation step can make neighboring pixels agree better while barely changing that broad drift. On a coarser grid the same drift spans fewer samples, so a correction can act across it more effectively. Restrict the residual to the coarse grid, solve an approximate correction problem there, prolongate that correction, and relax again on the fine grid. The hierarchy changes which error scales are easy to remove. The transfer operators and coarse problem must still match the discretization; downsampling the current image alone is not a multigrid solver.

Descending to coarser grids and returning with corrections gives the V-cycle its name. Its purpose differs from the image pyramid used for blending: the coarse levels solve error-correction equations, and the fine estimate is updated by the prolonged correction. The connection to Chapter 6.6, Pyramids and wavelets is the usefulness of scale, not the interchangeability of every reduce and expand routine.

6.9.7 Which solver when, and the edges🔗⧉

Step back and survey the toolbox, because the best solver depends on the problem. The premise never changes: $A^\top A$ is huge but sparse and structured, so we never go near a direct dense inverse. What varies is the shortcut the structure allows.

Stop according to the residual of the system being solved, $\lVert Qx-b\rVert\le\mathrm{atol}+\mathrm{rtol}\lVert b\rVert$, with an iteration limit and a reported convergence status. The data mismatch $\lVert Ax-y\rVert$ may remain nonzero at the optimum. Choosing regularization from the noise level and choosing solver tolerance are different decisions. Warm starts can help when successive problems are similar.

One edge marks the end of the linear toolkit. Everything here solves a quadratic objective, the least-squares fit with a Gaussian (Tikhonov) prior, whose minimizer is a single linear system. A sparsity or total-variation prior (the L1 penalties of Chapter 7.5, Image priors) makes the objective non-quadratic, and no single linear solve reaches its minimum. The standard move is to wrap a linear solve in an outer loop, iteratively reweighted least squares, ADMM, or a proximal method, some of which use these linear solvers as subproblems. A proximal-gradient step may instead need only operator products and a proximal map. So even for the nonlinear priors, this chapter is the inner engine; the outer loop is Chapter 7.5, Image priors' concern.

Before solving, check the adjoint on several random vectors and on vectors supported near the boundaries:

Choose $\epsilon$ for the arithmetic precision and operation size. Also check linearity and the symmetry of $Q$. A wrong adjoint can cause a wrong solution, loss of symmetry, or solver failure.