Add ADMM #1198
Add ADMM #1198
Conversation
odl/solvers/nonsmooth/admm.py
Outdated
| The functions ``f`` and ``g`` in the problem definition. They | ||
| need to implement the ``proximal`` method. | ||
| L : linear `Operator` | ||
| The linear operator that should be applied before ``f``. Its range |
odl/solvers/nonsmooth/admm.py
Outdated
| if callback is not None and not callable(callback): | ||
| raise TypeError('`callback` {} is not callable'.format(callback)) | ||
|
|
||
| # Initialize dual variables |
There was a problem hiding this comment.
Dual is probably the wrong name since we never look at convex conjugates, Just u and z should do.
examples/solvers/admm_tomography.py
Outdated
| # Combine functionals, order must correspond to the operator K | ||
| g = odl.solvers.SeparableSum(l2_norm, l1_norm) | ||
|
|
||
| # --- Select solver parameters and solve using PDHG --- # |
odl/solvers/nonsmooth/admm.py
Outdated
| f.proximal(tau)(x, out=x) | ||
|
|
||
| # L x^(k+1) + u^k | ||
| L(x, out=tmp_ran) |
There was a problem hiding this comment.
What is the difference between linearized ADMM and PDHG? In any case, I am pretty sure it is possible to rewrite your code calling "L" only once per iteration.
There was a problem hiding this comment.
The main difference is that ADMM uses the proximal of f in contrast to the proximal of g^* as in PDHG, so it's not a primal-dual method in that sense. But it's possible that one can still interpret it that way, dunno.
Anyway, you're right, I should be able to scrap one L evaluation, potentially at the cost of a new temporary.
odl/solvers/nonsmooth/admm.py
Outdated
| # x^k - (tau/sigma) L^*(Lx^k + u^k - z^k) | ||
| x.lincomb(1, x, -tau / sigma, tmp_dom) | ||
| # x^(k+1) = prox[tau*f](_) | ||
| f.proximal(tau)(x, out=x) |
There was a problem hiding this comment.
I like this line as it is very elegant. However, currently some of ODL's functionals do not support this (input and output are the same variable), see #1181.
There was a problem hiding this comment.
True, but this is not the place to implement workarounds. This code has to be able to assume that operators work with aliased input and output.
|
Obviously I should add the solver to the test suite. |
examples/solvers/admm_tomography.py
Outdated
| angle_partition = odl.uniform_partition(0, np.pi, 360) | ||
| # Detector: uniformly sampled, n = 512, min = -30, max = 30 | ||
| detector_partition = odl.uniform_partition(-30, 30, 512) | ||
| geometry = odl.tomo.Parallel2dGeometry(angle_partition, detector_partition) |
There was a problem hiding this comment.
I prefer if we use the parallel_beam_geometry helper in these small examples
examples/solvers/admm_tomography.py
Outdated
| l2_norm = odl.solvers.L2NormSquared(ray_trafo.range).translated(data) | ||
|
|
||
| # Isotropic TV-regularization i.e. the l1-norm | ||
| l1_norm = 0.015 * odl.solvers.L1Norm(gradient.range) |
examples/solvers/admm_tomography.py
Outdated
| @@ -0,0 +1,87 @@ | |||
| """Total variation tomography using ADMM. | |||
examples/solvers/admm_tomography.py
Outdated
| Where ``A`` is a parallel beam ray transform, ``grad`` the spatial | ||
| gradient and ``g`` is given noisy data. | ||
|
|
||
| See the documentation of the `admm` method for further details. |
examples/solvers/admm_tomography.py
Outdated
| # --- Select solver parameters and solve using PDHG --- # | ||
|
|
||
| # Estimated operator norm, add 10 percent to ensure ||K||_2^2 * sigma * tau < 1 | ||
| op_norm = 1.1 * odl.power_method_opnorm(op) |
There was a problem hiding this comment.
Overall we should move towards op.norm() in these cases now that it's implemented.
There was a problem hiding this comment.
Don't blame me for all the above, it's literally copied from the PDHG example (except for the admm typos).
There was a problem hiding this comment.
I guess we could make a "mother of all PRs" to fix that
odl/solvers/nonsmooth/admm.py
Outdated
| # tmp_ran <- L x^(k+1) | ||
| L(x, out=tmp_ran) | ||
| # z^(k+1) <- prox[sigma*g](L x^(k+1) + u^k) | ||
| g.proximal(sigma)(tmp_ran + u, out=z) # 1 copy here |
There was a problem hiding this comment.
why not tmp_ran += u on the line above?
There was a problem hiding this comment.
I need L(x) later to add the new u. The alternative would be to use another temporary, but that seems overkill to me.
| assert x == x0 | ||
|
|
||
| # Check that a provided callback is actually called | ||
| class CallbackTest(Callback): |
There was a problem hiding this comment.
class CallbackTest(Callback):
was_called = False
def __call__(self, *args, **kwargs):
self.was_called = True
|
Done fixing, will merge after CI |
|
In particular, I'll keep the roles of |
Depending on the outcome of #1197,
fandgmay change roles.