osmlem

odl.solvers.iterative.statistical.osmlem(op, x, data, niter, callback=None, **kwargs)[source]

Ordered Subsets Maximum Likelihood Expectation Maximation algorithm.

This solver attempts to solve:

max_x L(x | data)

where L(x, | data) is the likelihood of x given data. The likelihood depends on the forward operators op[0], ..., op[n-1] such that (approximately):

op[i](x) = data[i]
Parameters:
opsequence of Operator

Forward operators in the inverse problem.

xop.domain element

Vector to which the result is written. Its initial value is used as starting point of the iteration, and its values are updated in each iteration step. The initial value of x should be non-negative.

datasequence of op.range element-like

Right-hand sides of the equation defining the inverse problem.

niterint

Number of iterations.

callbackcallable, optional

Function called with the current iterate after each iteration.

Other Parameters:
sensitivitiesfloat or op.domain element-like, optional

The algorithm contains an A^T 1 term, if this parameter is given, it is replaced by it. Default: op[i].adjoint(op[i].range.one())

See also

mlem

Ordinary MLEM algorithm without subsets.

loglikelihood

Function for calculating the logarithm of the likelihood

Notes

Given forward models A_i, and data g_i, i = 1, ..., M, the algorithm attempts to find an x that maximizes:

\prod_{i=1}^M P(g_i | g_i \text{ is } Poisson(A_i(x)) \text{ distributed}).

The algorithm is explicitly given by partial updates:

x_{n + m/M} = \frac{x_{n + (m - 1)/M}}{A_i^* 1} A_i^* (g_i / A_i(x_{n + (m - 1)/M}))

for m = 1, ..., M and x_{n+1} = x_{n + M/M}.

The algorithm is not guaranteed to converge, but works for many practical problems.

References

Natterer, F. Mathematical Methods in Image Reconstruction, section 5.3.2.