bartorch.nlop.GaussNewton

class bartorch.nlop.GaussNewton(image_shape, pattern=None, trajectory=None, *, iterations=1, redu=2.0, alpha_min=0.0, batch=1, cg_maxiter=30, cg_tol=0.0, cg_lambda=0.0, kspace_shape=None, coil_shape=None, weights=None, basis=None, mask=None, sobolev=(220.0, 32.0), c=1.0, real=False, sos=False, toeplitz=True)

One or more Gauss-Newton steps of nlinv, as an operator that differentiates.

Takes (y, xn, x0, alpha) and returns the iterate:

y

The data as the step takes it – coil images, which prepare() makes from k-space and a sampling pattern.

xn

The iterate: the image and the coil coefficients in one flat vector. start() makes the one BART starts from, and decompose() reads one back.

x0

The centre the step regularises towards, which nlinv leaves at the starting point.

alpha

The Tikhonov weight. A number is taken as one; with iterations above one it is the initial weight, divided by redu after each step and never taken below alpha_min, which is what IRGNM does.

Every one of those carries a gradient, which is what makes this an unrolled network’s cell rather than a solver’s inside. A denoiser between two of these – or a term in the regularisation centre – is NLINV-Net.

Parameters:

• image_shape (tuple of int) – Coil-image shape, (coils, *spatial), as the linear encodings take it.

• pattern (tensor, optional) – Only to give the shape BART should expect; the pattern itself is an argument of prepare(), because a network is handed one per call. Ones over the image by default.

• trajectory (tensor, optional) – Off the grid. Its shape is what is kept, for the same reason.

• iterations (int) – Steps. One is a step; more is that many, with the weight decaying, wrapped in BART’s own gradient checkpointing.

• redu (float) – What the weight is divided by after each step.

• alpha_min (float) – What the weight decays towards.

• batch (int) – Independent copies of the model, stacked on BART’s batch axis. This is the leading axis of every argument.

• cg_maxiter (int, float, float) – The conjugate gradients inside each step: iter_conjgrad_conf’s maxiter, tol and l2lambda.

• cg_tol (int, float, float) – The conjugate gradients inside each step: iter_conjgrad_conf’s maxiter, tol and l2lambda.

• cg_lambda (int, float, float) – The conjugate gradients inside each step: iter_conjgrad_conf’s maxiter, tol and l2lambda.

• weights (tensor, optional) – As NonlinearSense takes them. These are held by the model.

• basis (tensor, optional) – As NonlinearSense takes them. These are held by the model.

• mask (tensor, optional) – As NonlinearSense takes them. These are held by the model.

• sobolev (tuple of float) – (a, b) of the coil weighting (1 + a |k|^2)^(-b/2).

• real (bool) – Constrain the image to be real (nlinv -R).

• sos (bool) – Normalise the coils by their root sum of squares.

• toeplitz (bool) – Off the grid, use a point-spread convolution for the normal operator.

Notes

The step is the library’s, and so is what it costs. Each one runs a conjugate-gradient solve whose backward pass is another, by norm_inv’s implicit differentiation rather than by unrolling – norm_inv_der_src and norm_inv_adj_src in nlops/norm_inv.c. So the memory of a K-step unroll is K cells and not K times the inner iterations, and BART’s nlop_checkpoint_create_F takes it down further.

The coil weighting is worth knowing about before reading a gradient. BART weights the coil half of the state by (1 + a |k|^2)^(-b/2), and its default b = 32 is a sixteenth power: over the state of a small fit the gradient of that half spans tens of decades and its tail runs below float32’s smallest normal number. Below that edge the arithmetic is the platform’s business rather than the library’s – a right-hand side whose norm is no longer a normal number is one BART’s checkeps declines to iterate on, and the solve comes back untouched, with Warning: data corrupted in the log and a gradient of zeros. Forward, none of this matters and the default is what nlinv reconstructs with; a gradient that has to be meaningful in the coil coefficients wants a gentler weighting, which is what sobolev=(220.0, 8.0) is.

__init__(image_shape, pattern=None, trajectory=None, *, iterations=1, redu=2.0, alpha_min=0.0, batch=1, cg_maxiter=30, cg_tol=0.0, cg_lambda=0.0, kspace_shape=None, coil_shape=None, weights=None, basis=None, mask=None, sobolev=(220.0, 32.0), c=1.0, real=False, sos=False, toeplitz=True)

Methods

__init__(image_shape[, pattern, trajectory, ...])
adjoint(dy)	DF(x)^H dy at the last evaluated point, for one input and one output.
chain(other, *[, output, input])	Output output of self into input input of other.
combine(other)	self and other side by side, sharing nothing.
decompose()	x -> (image, sensitivities), through the model's transforms.
del_out([output])	Drop an output, and everything computed only for it.
derivative(dx)	DF(x) dx at the last evaluated point, for one input and one output.
dup([a, b])	Make two inputs of the same shape one input, kept at a.
flatten([inputs_only])	Every input as one flat vector, and every output as another.
forward(*xs)	(y, xn, x0, alpha) -> the iterate.
jacobian([output, input])	DF/dx_input of one output, as a LinearOperator.
join()	(image, coefficients) -> x: split() the other way.
linearize(*xs)	The derivative at x, as a LinearOperator.
link([output, input])	Tie an output back into an input; both arguments go away.
permute_inputs(perm)	Reorder the inputs: the new input i is the old perm[i].
permute_outputs(perm)	Reorder the outputs: the new output o is the old perm[o].
pin(input, value)	Fix one input to value; the input goes away.
prepare()	(kspace, pattern) -> the data a step takes.
reshape_input(input, shape)	This operator with one input's shape written differently.
reshape_output(output, shape)	This operator with one output's shape written differently.
shift_input(new, old)	Move one input to another position, the rest closing up behind it.
shift_output(new, old)	Move one output to another position, the rest closing up behind it.
split()	x -> (image, coefficients), with no transform applied.
stack_inputs(a, b, axis)	Make two inputs one, concatenated along a C-order axis.
stack_outputs(a, b, axis)	Make two outputs one, concatenated along a C-order axis.
start([batch, device])	The iterate nlinv starts from: an image of ones, and no coils.
weight(alpha[, device])	alpha as the operator takes it.

Attributes

data_shape	coil images, batch first.
device
ishape	The domain, for an operator with one input.
ishapes	The shape of each argument, outputs and inputs kept apart.
oshape	The codomain, for an operator with one output.
oshapes
output_shapes	The same for what comes back.
shapes	What each argument is, without the axes BART leaves empty.
state_shape	the image and the coil coefficients laid end to end, batch first.