Regularized least squares#
The mathematics behind Solvers and
Regularizers: the estimation problem, the
properties of the encoding that make it ill-posed, the proximal operators and
the convergence conditions of each iteration. The code is
src/bartorch/optim/_linear.py (the solvers), src/bartorch/optim/_blocks.py
(one step of each), src/bartorch/priors/ (the terms) and
src/bartorch/apps/_pics.py (the choice of iteration); the tests are
tests/test_optim_iterators.py, tests/test_solve.py, tests/test_priors.py
and tests/test_apps.py.
Measurement model#
The data \(y \in \mathbb{C}^M\) and the image \(x \in \mathbb{C}^N\) are related by \(y = Ax + \varepsilon\), with \(A\) the encoding operator and \(\varepsilon\) noise that is complex Gaussian with independent, identically distributed entries after prewhitening. \(A\) is applied without being formed as a matrix. The adjoint \(A^H\) is defined by \(\langle Ax, y\rangle = \langle x, A^H y\rangle\), and the normal operator \(A^HA\) is self-adjoint and positive semidefinite. Every solver needs only applications of \(A\) and \(A^H\), or of \(A^HA\). \(A^Hy\) is the simplest estimate; for Cartesian SENSE it is the coil combination of the zero-filled inverse FFTs and retains the aliasing.
Least squares and conditioning#
The least-squares estimate solves the normal equations \(A^HA\hat{x} = A^Hy\). With the singular value decomposition \(A = \sum_i \sigma_i u_i v_i^H\), the minimizer of least norm is
Property |
Condition |
Consequence |
|---|---|---|
Rank deficiency |
some \(x_0 \ne 0\) satisfy \(Ax_0 = 0\) |
\(x^\dagger + x_0\) fits the data equally well; the data carry no information about \(x_0\) |
Ill-conditioning |
full column rank, small \(\sigma_{\min}\) |
unique minimizer, with the noise along \(v_i\) amplified by \(1/\sigma_i\) |
An undersampled parallel-imaging acquisition measures fewer k-space positions than the image has voxels, but each with several coils. Where the coils separate the voxels that alias onto one another, \(A\) has full column rank and the problem is ill-conditioned rather than underdetermined; the spatially resolved noise amplification is the g-factor of SENSE.[1] Beyond the acceleration the coil geometry supports, and outside the region the coils cover, \(A\) is rank deficient. The eigenvalues of \(A^HA\) are the \(\sigma_i^2\), so \(\kappa(A^HA) = \kappa(A)^2\).
Conjugate gradients and early stopping#
CG on the normal equations, started from \(x_0 = 0\),
keeps its iterates in the range of \(A^H\) and so converges to \(x^\dagger\); a
nonzero \(x_0\) keeps its null-space component. On noisy data the iteration
shows semi-convergence: early iterates recover the components with large
\(\sigma_i\), later ones fit the noise along those with small \(\sigma_i\). The
iteration count acts as a regularization parameter,[2] as in iterative
SENSE,[3] with a best value that depends on the data and the
noise level.
Regularization functionals#
Regularization |
\(R(x)\) |
Effect |
|---|---|---|
Tikhonov |
\(\tfrac12 \lVert x \rVert_2^2\) |
Solution \((A^HA + \lambda I)^{-1}A^Hy\): the component along \(v_i\) is weighted by \(\sigma_i^2/(\sigma_i^2+\lambda)\). Not a spatial smoothness prior. |
Generalized Tikhonov |
\(\tfrac12 \lVert Gx - b \rVert_2^2\) |
Shrinkage toward a reference \(b\), or, with \(G\) a finite difference, a quadratic smoothness prior |
Sparsity in a transform |
\(\lVert \Psi x \rVert_1\) |
Few significant coefficients, for example of a wavelet transform[4] |
Total variation |
\(\sum_r \lVert (\nabla x)_r \rVert_2\) |
Piecewise-constant images[5] |
Locally low rank |
\(\sum_b \lVert B_b x \rVert_*\) |
Image blocks, as matrices over a contrast or coefficient axis, of low rank |
The quadratic functionals keep the problem linear least squares, which
CG solves with Tikhonov
terms stacked beneath the encoding. The others are not differentiable.
Proximal operators#
The data term \(f(x) = \tfrac12\lVert Ax-y\rVert_2^2\) has gradient \(A^H(Ax-y)\), Lipschitz continuous with \(L = \lambda_{\max}(A^HA)\). A nonsmooth \(g\) enters through
defined for proper, closed, convex \(g\) and \(\tau > 0\). For the \(\ell_1\) norm
of a complex array it is soft thresholding of the modulus with the phase kept,
zero where \(v_i = 0\); thresholding a group jointly (the gradient components of
isotropic total variation, or a term’s joint_axes) replaces \(\lvert v_i\rvert\)
by the group’s \(\ell_2\) norm. For the indicator function of a convex set it is
the projection onto the set.
A term of bartorch.priors is \(g(Gx)\). For unitary \(G\),
\(\operatorname{prox}_{\tau g\circ G}(v) = G^H\operatorname{prox}_{\tau g}(Gv)\),
and BART’s wavelet and locally low-rank terms apply their transform inside the
proximal operator and have \(G = I\). This is exact for the Haar wavelet on even
lengths, which is orthogonal, and for non-overlapping blocks. BART’s Daubechies
and CDF 4/4 transforms in symmetric extension produce more coefficients than
samples and are not orthogonal (test_synthesis_is_not_the_adjoint_of_a_redundant_transform
in tests/test_wavelet_thresh.py); with them, with overlapping blocks, and with
the random shifts both terms draw by default, the step approximates that
proximal operator. For a general \(G\) the solver has to apply \(G\) separately.
Iterations#
Algorithm |
Problem and splitting |
Transform \(G\) |
Per iteration |
Convergence conditions |
|---|---|---|---|---|
CG |
\(\min \tfrac12\lVert Ax-y\rVert^2 + \tfrac{\lambda}{2}\lVert x\rVert^2 + \sum_i \tfrac{w_i}{2}\lVert G_ix-b_i\rVert^2\) |
Any, inside the quadratic terms |
One normal operator |
Positive semidefinite normal operator |
IST |
\(x \leftarrow \operatorname{prox}_{\tau g}(x - \tau\nabla f(x))\) |
\(I\) only |
One normal operator, one proximal operator |
Convex \(f\), \(g\); \(0 < \tau \le 1/L\); objective error \(O(1/k)\) |
FISTA |
IST with Nesterov momentum |
\(I\) only |
As IST |
As IST; objective error \(O(1/k^2)\)[6] |
ADMM |
\(f(x) + \sum_j g_j(z_j)\) subject to \(z_j = G_jx\) |
Any |
A CG solve with \(A^HA + \rho\sum_j G_j^HG_j\), one proximal operator per term, a dual update |
Convex terms, any \(\rho > 0\)[7] |
Primal-dual (Chambolle–Pock) |
Saddle-point form with a dual variable per \(g_j(G_jx)\) |
Any |
\(A\), \(G_j\), their adjoints and proximal operators; no inner solve |
Convex terms, \(\sigma\tau\lVert K\rVert^2 < 1\) for the stacked operator \(K\)[8] |
POCS |
Repeated projection onto convex sets |
Inside the projections |
One projection per set |
Closed convex sets with nonempty intersection |
IST, FISTA and the primal-dual iteration take BART’s default step of 0.95,
not divided by \(L\) unless eigen=True. ADMM’s maxiter is BART’s budget of
inner conjugate-gradient iterations over the whole run, not a count of outer
steps. apps.pics chooses the iteration as BART’s italgo_choose does: CG
without terms or with an \(\ell_2\) term only, NIHT for the hard-thresholding
terms, ADMM when a total-variation-type or imaginary-part term appears, and
otherwise FISTA for a first term and ADMM once there is a second.
References#
Further reading: Parikh N, Boyd S. Proximal algorithms. Found Trends Optim 1(3):127–239 (2014). doi:10.1561/2400000003