ENH: add nonlinear version of ADMM, closes #1201 #1202
Conversation
kohr-h
force-pushed
the
issue-1201__nonlinear_admm
branch
from
d0c81ad
to
7c4f701
Compare
|
Ready for review. I also found the bug that was causing headaches for aliased input and output, see #1181. It was in |
| # Define number of energy bins on both sides as well as the source spectrum | ||
| num_bins_reco = 4 | ||
| num_bins_data = 2 | ||
| assert num_bins_reco >= num_bins_data |
There was a problem hiding this comment.
Not at all sure we want asserts in examples
| # tau^k <- opnorm_factor / (delta * ||dL(x^k)||^2) | ||
| A = L.derivative(x) | ||
| xstart = A.domain.one() if i == 0 else x | ||
| A_norm = power_method_opnorm(A, xstart, maxiter=opnorm_maxiter) |
There was a problem hiding this comment.
power method in each iterate is going to be expensive not to say un-usable. Is there no way around this?
| @@ -0,0 +1,169 @@ | |||
| """Total variation spectral tomography using preconditioned nonlinear ADMM. | |||
There was a problem hiding this comment.
Examples like this truly need to go into a "complicated" folder so to say, this is getting out of hand.
Nothing against the example per se, just that we want beginners to find them "in the right order".
| group_size = num_bins_reco // num_bins_data | ||
| # Note: the spectrum (i.e. intensity per channel) is additive on the data | ||
| # side, i.e., adding more entries adds to the total signal in the data bins | ||
| source_spectrum = [8e4, 1e5, 8e4, 6e4] |
There was a problem hiding this comment.
I'll have to work a bit on this example anyway. Currently it's a bit silly since for n x n problems, you can actually linearize.
| @@ -159,3 +170,212 @@ def admm_linearized_simple(x, f, g, L, tau, sigma, niter, **kwargs): | |||
| u = L(x) + u - z | |||
| if callback is not None: | |||
| callback(x) | |||
|
|
|||
|
|
|||
| def admm_precon_nonlinear(x, f, g, L, delta, niter, sigma=None, **kwargs): | |||
There was a problem hiding this comment.
Not sure that "precon" should have its own function, IMO it should just be optional and this can be called "admm_nonlinear"
There was a problem hiding this comment.
Also, are we sure that this is the admm_nonlinear method, or can we get a naming collision here in the future?
There was a problem hiding this comment.
As I understand it, the preconditioning is part of the method, and the way it is used makes the method explicit, i.e., everything can be expressed in terms of proximals. The ADMM method has implicit steps in the form of some argmin which is not a proximal.
| if not 0 < sigma < 1 / delta: | ||
| raise ValueError( | ||
| '`sigma` must lie strictly between 0 and 1/delta`, got ' | ||
| 'sigma={:.4} and (1/delta)={:.4}'.format(sigma_in, 1 / delta)) |
There was a problem hiding this comment.
4 digits is bad practice for errors, just print all of them.
| y = L.range.zero() | ||
| # mu = [mu_old +] delta * (L(x) [- y]) | ||
| mu = delta * L(x) | ||
| # mubar = 2 * mu [- mu_old] |
| out.lincomb(1.0 / (1 + 2 * sig * lam), x) | ||
| else: | ||
| out.lincomb(1.0 / (1 + 2 * sig * lam), x, | ||
| 2 * sig * lam / (1 + 2 * sig * lam), g) |
There was a problem hiding this comment.
mixing 1.0 and 1 here, lets go for one of them
| @@ -72,5 +74,65 @@ def test_admm_lin_l1(): | |||
| assert all_almost_equal(x, data_1, places=2) | |||
|
|
|||
|
|
|||
| def test_admm_lin_vs_simple(): | |||
There was a problem hiding this comment.
I guess these are old, or do we intend for them to stay? If so, I'd move the "simple" guy into this file
kohr-h
force-pushed
the
issue-1201__nonlinear_admm
branch
from
62a6ec4
to
4bf95b1
Compare
There was a problem hiding this comment.
Submitting a review I started long ago.
Overall my main issue with this PR is that, as we discussed with @mehrhardt before (I think he even made an issue), is that the example is very complicated and I'm not at all sure if this fits in the standard examples folder. Do we have any plan for this?
| where ``||.||_{nuc, 1, 2}`` is the nuclear L1-L2 norm and ``A`` a 2 x 4 | ||
| operator matrix | ||
|
|
||
| (c_1 * exp(-R) c_2 * exp(-R) ) |
| A = ( c_3 * exp(-R) c_4 * exp(-R)) | ||
|
|
||
| with constants ``c_i`` and the parallel beam ray transform ``R``. In other | ||
| words, ``A`` acts on an input ``x`` with 4 components as follows: |
|
|
||
| # Create the forward operator | ||
| ray_trafo = odl.tomo.RayTransform(single_energy_space, geometry, | ||
| impl='astra_cpu') |
There was a problem hiding this comment.
Don't hardcode backend. With that said we must fix the astra issue: astra-toolbox/astra-toolbox#109, perhaps we need to do so locally.
There was a problem hiding this comment.
This issue is actually fixed now
| The linear operator that is composed with ``g`` in the problem | ||
| definition. It must fulfill ``L.domain == f.domain`` and | ||
| ``L.range == g.domain``. | ||
| Linear operator composed with ``g`` in the problem definition. |
TODO: