Nonlinear forward models

Nonlinear forward models#

TL;DR

  • Joint estimation of the image and the coil sensitivities (nlinv) and model-based parameter estimation (moba) have forward operators that are nonlinear in the unknowns.

  • Both are solved by iteratively regularized Gauss-Newton: each step solves a linearized least-squares problem with a Tikhonov term whose weight decreases geometrically, and the number of steps acts as a regularization parameter.

  • nlinv resolves the ambiguity between image and sensitivities with a Sobolev weighting of the coils; the data term of a signal model is nonconvex in its parameters, so the result depends on the starting point.

  • A model-based reconstruction estimates parameter maps directly; a subspace reconstruction keeps the forward model linear and fits the parameters afterwards.

Inverse problems and their solvers assumes a known, linear forward operator. Two common MRI reconstructions do not satisfy that assumption: the coil sensitivities can be unknowns alongside the image, and the image can be a function of physical parameters. Both lead to a nonlinear operator \(F\) and are solved by the same Gauss-Newton method.

Two nonlinear models#

Joint estimation of image and sensitivities. When the acquisition has no autocalibration (ACS) region large enough for a separate calibration, or the coil sensitivity maps are to be estimated from all of the data, both are unknowns:

\[ y_c = P F \left( S_c \cdot x \right). \]

The operator is bilinear: linear in \(x\) for fixed \(S\), linear in \(S\) for fixed \(x\), and nonlinear in the pair. This is nonlinear inversion, BART’s nlinv,[1] which bartorch.tools.nlinv() runs; its sensitivities are the calibration for a subsequent linear reconstruction.

Model-based parameter estimation. A quantitative experiment acquires a series of contrasts whose dependence on tissue parameters \(\theta\) is given by a signal model \(M\) — a mono-exponential decay in \(T_2\), an inversion recovery in \(T_1\), or a Bloch simulation of the sequence:

\[ y_{c,e} = P_e F \left( S_c \cdot M_e(\theta) \right), \]

with \(e\) the contrast index. The unknowns are the parameter maps. This is model-based reconstruction, BART’s moba,[2][3][4] which bartorch.apps.moba() performs over a TorchSim signal model; bartorch.apps.mobafit() fits the same models voxel by voxel to reconstructed images.

A Gauss-Newton solver requires the value \(F(x)\), the derivative \(DF_x\) at a point as a linear operator, and its adjoint \(DF_x^H\). A NonlinearOperator provides the three, and F.linearize(x) returns \(DF_x\) as a LinearOperator.

Iteratively regularized Gauss-Newton#

Each Gauss-Newton step replaces \(F\) by its linearization at the current iterate \(x_k\) and solves the resulting linear problem with a Tikhonov term centred on a reference \(x_{\mathrm{ref}}\):

\[ x_{k+1} = \arg\min_x \; \left\lVert DF_{x_k}(x - x_k) - \left(y - F(x_k)\right) \right\rVert_2^2 + \alpha_k \lVert x - x_{\mathrm{ref}} \rVert_2^2 , \]

that is, \(x_{k+1} = x_k + (DF^H DF + \alpha_k I)^{-1}\left[DF^H (y - F(x_k)) - \alpha_k (x_k - x_{\mathrm{ref}})\right]\) with \(DF = DF_{x_k}\). The weight decreases geometrically, \(\alpha_{k+1} = (\alpha_k - \alpha_{\min})/q + \alpha_{\min}\) with \(q = 2\) by default.[5] The early steps are strongly regularized and change the iterate only along well-determined directions; later steps admit the poorly determined ones. The number of steps is a regularization parameter: stopping early gives a smoother estimate, and continuing eventually fits the noise.

IRGNM implements this loop and IRGNMBlock one step of it. Without an inner solver the linear problem is solved by conjugate gradients inside BART, as in nlinv; with inner= it is passed to a solver of bartorch.optim, and that solver’s regularization terms add a penalty to the step.

Identifiability#

The bilinear model has a symmetry: \(S_c \cdot x\) is unchanged by \(S_c \mapsto \gamma S_c\) and \(x \mapsto x / \gamma\) for any nonzero function \(\gamma(r)\), so the data do not determine how structure is shared between the two factors. nlinv resolves this with the prior knowledge that sensitivities are smooth, built into the model: the coil unknown is a k-space representation \(\hat{s}\), and the sensitivities are

\[ S = \mathcal{F}^{-1}\!\left[ (1 + a \lvert k \rvert^2)^{-b/2}\, \hat{s} \right], \]

a Sobolev weighting that attenuates the high spatial frequencies of the coil estimate. The IRGNM penalty on \(\hat{s}\) is then a penalty on the Sobolev norm of \(S\), and the image receives the high-frequency structure. The product \(x \cdot \sqrt{\sum_c \lvert S_c \rvert^2}\) is invariant to \(\gamma\) up to its phase, which is why nlinv reports the image multiplied by the root sum of squares of the estimated sensitivities. NonlinearSense is this model; CoilSense() is the unweighted product in front of any linear encoding, which leaves the regularization of the coils to the caller.

Signal models have no such symmetry — the model fixes the meaning of each map — but the data term is nonconvex in \(\theta\), so the result depends on the starting point. The models of bartorch.nlop impose bounds by solving for a transformed variable, so that every iterate stays within the bounds; initial() builds a starting point from parameter values and split() converts a solution back to named maps in physical units.

Approaches to parameter mapping#

Approach

Unknown

Forward model

Cross-contrast information in the reconstruction

Problem

Output

Reconstruction, then voxel-wise fit

One image per contrast

Linear, \(P_e F S\)

Only through a joint regularizer, if one is used; none if the contrasts are reconstructed separately

Linear reconstruction, then a nonlinear fit per voxel

Contrast images, then parameter maps

Subspace reconstruction

Coefficient maps of a low-dimensional basis \(\Phi\)

Linear, \(P_e F S\, \Phi\)

The signal evolution is restricted to the span of \(\Phi\)[6][7]

Linear; convex with a convex regularizer

Coefficient maps, then parameter maps by dictionary matching or fitting

Model-based reconstruction

Parameter maps \(\theta\)

Nonlinear, \(P_e F S\, M(\theta)\)

The signal model couples every contrast to the same parameters

Nonlinear, nonconvex

Parameter maps

A model-based reconstruction estimates fewer unknowns from the same data — for a multi-echo experiment, a few maps rather than one image per echo — and applies regularization to the maps. This improves the estimate when the model describes the signal; a signal the model does not describe, such as partial volume of two tissues or an imperfect refocusing, becomes a bias in the maps. A subspace reconstruction restricts the signal to a basis derived from a dictionary of simulated signals, keeps a linear forward model, and defers the nonlinear estimation of the parameters to a separate step. Subspace-constrained T1 mapping and Parameter maps straight from k-space apply the second and third routes.

Signal models and their derivatives#

The signal models — InversionRecovery(), MultiEcho(), Bloch() — are TorchSim simulators presented as BART nonlinear operators. A signal model is voxel-wise: the signal of a voxel depends only on that voxel’s parameters, so one forward-mode automatic differentiation pass gives the derivative for the whole volume and no Jacobian matrix is formed. TorchOperator does the same for any differentiable PyTorch function. Composition with @ places a model in front of an encoding; the derivative of E @ M at \(\theta\) is \(E\, DM_\theta\). Differentiation through reconstruction describes how gradients pass through the Gauss-Newton steps themselves.

References#