Issue1101 pdhg examples by mehrhardt · Pull Request #1191 · odlgroup/odl

mehrhardt · 2017-10-11T17:15:34Z

Here are the two examples I promised a few weeks ago. I suppose the documentation can be improved but first I would like to have some more general feedback. For these, I also needed the Huber L1 norm for which I can compute the prox of the convex conjugate and its function evaluation.

Please let me know what you think about these.

pep8speaks · 2017-10-11T17:15:37Z

Hello @mehrhardt! Thanks for updating the PR.

No PEP8 issues.

Comment last updated on October 26, 2017 at 09:26 Hours UTC

adler-j

Very good addition.

Some minor comments on style etc, and I also need a discussion on how we should handle these rather long examples, we need to somehow indicate that they are "advanced" in some sense.

adler-j · 2017-10-11T19:18:15Z

+
+        min_{x >= 0} ||x - d||_1 + lam TV_gamma(x)
+
+where ``grad`` the spatial gradient and ``d`` is given noisy data.


there is no "grad" above, but i personally like the wrong version which is ||grad(x)||_1

True, modified it similar to your suggestion.

adler-j · 2017-10-11T19:18:42Z

+# Rescale max to 1
+image /= image.max()
+
+# Discretized spaces


adler-j · 2017-10-11T19:19:58Z

+# Create space element of ground truth
+orig = space.element(image.copy())
+
+# Add noise and convert to space element


there is no conversion in the line below

adler-j · 2017-10-11T19:20:16Z

+noisy = odl.phantom.salt_pepper_noise(orig)
+
+# Gradient operator
+gradient = odl.Gradient(space, method='forward')


method='forward' is the default, so no need to add it

adler-j · 2017-10-11T19:20:36Z

+reg_param = 1
+
+# l1 data matching
+l1_norm = 1 / reg_param * odl.solvers.L1Norm(space).translated(noisy)


why 1 / reg_param

The regularization parameter changes the strong convexity of f*. Thus, it is sometimes better to have it as part of g. Here it is 1, so it doesn't matter.

adler-j · 2017-10-11T19:30:17Z

                (1 / self.sigma) * self.functional.proximal(self.sigma))


+class HuberL1L2(Functional):


Adding this partially solves #597 which should be mentioned i guess

Also not 100% sure about the naming, is there no way to create something like Huber and control the L1L2 part via parameters?

Good point, maybe this is possible. So far, I only came across the L1L2 version. Maybe we can add a TODO note and implement this when there is need for this (which may be never?!?)?

I still feel that we need to improve the naming here before we go ahead. Is there any "reasonable" way to make this simply Huber?

adler-j · 2017-10-11T19:30:37Z

+
+        Parameters
+        ----------
+        space : `DiscreteLp` or `FnBase`


DiscreteLp is a subclass of FnBase so no need to mention

adler-j · 2017-10-11T19:31:18Z

+
+        Examples
+        --------
+


Remove extra blankline

adler-j · 2017-10-11T19:32:00Z

+        >>> H = alpha * odl.solvers.HuberL1L2(X, gamma)
+        >>> L1 = alpha * odl.solvers.GroupL1Norm(X, 2)
+        >>> abs(H(x) - L1(x)) < 1e-10
+        """


It would be tremendous with a Notes section here where you give the definition of the functional

Done, please have a look.

adler-j · 2017-10-11T19:37:29Z

+    @property
+    def convex_conj(self):
+        '''The convex conjugate'''
+        return FunctionalQuadraticPerturb(GroupL1Norm(self.domain,


split this on two lines for readability

mehrhardt · 2017-10-18T10:23:41Z

Anything left to do here before merging?

adler-j

Some minor comments, overall the only major thing is fixing the HuberL1L2 naming

adler-j · 2017-10-12T10:11:20Z

+        --------
+        Compare HuberL1L2 and L1 for vanishing smoothing ``\\gamma=0``
+
+        >>> import odl


odl is auto-imported in doctests, no need to import it

I know. I thought we discussed this a few months ago and you preferred this so that one can copy-paste the examples into an ipython shell?

adler-j · 2017-10-12T10:11:30Z

+        Compare HuberL1L2 and L1 for vanishing smoothing ``\\gamma=0``
+
+        >>> import odl
+        >>> X = odl.uniform_discr([0, 0], [1, 1], [5, 5])


we usually call it "space"

adler-j · 2017-10-12T10:11:55Z

+        >>> import odl
+        >>> X = odl.uniform_discr([0, 0], [1, 1], [5, 5])
+        >>> x = odl.phantom.white_noise(X)
+        >>> alpha = 2


alpha here is not really needed to convey the example

adler-j · 2017-10-12T10:12:25Z

+        >>> X = odl.uniform_discr([0, 0], [1, 1], [5, 5])
+        >>> x = odl.phantom.white_noise(X)
+        >>> alpha = 2
+        >>> gamma = 0


Might as well write gamma=0 on the line below. Keeping the examples compact is good for readability

adler-j · 2017-10-12T10:13:07Z

+
+    def _call(self, x):
+        '''Return the HuberL1L2-norm of ``x``.'''
+        if isinstance(self.domain, ProductSpace):


What about product-space of product-space, etc?

Good point. I don't know how to handle this. You are the expert here I guess.

adler-j · 2017-10-18T14:01:53Z

                (1 / self.sigma) * self.functional.proximal(self.sigma))


+class HuberL1L2(Functional):


I still feel that we need to improve the naming here before we go ahead. Is there any "reasonable" way to make this simply Huber?

mehrhardt · 2017-10-18T14:32:05Z

I made the changes you requested. The product space issue is hard to handle I suppose. It needs to be clear what "norm" is locally defined and not for all norms does this construction make sense. Probably it is best to check for these things and to throw an error if one tries something strange?

Regarding the naming, I changed it to "Huber" as this is the only HuberNorm I know of that is being used. One could then generalize this later if need be without breaking backward compatibility by ensuring that the default parameter will refer to this version.

It would also be good to merge this version with the one of Axel #1195 to get the best of two worlds.

kohr-h

A somewhat larger PR, therefore more comments. All in all a very nice addition, this is for the final polish.

kohr-h · 2017-10-19T13:33:58Z

@@ -0,0 +1,135 @@
+"""Total variation denoising using PDHG.
+
+This exhaustive example solve the L1-HuberTV problem


solve -> solves

kohr-h · 2017-10-19T13:39:33Z

+obj_fun = l1_norm + huber * gradient
+
+# Strong convexity of "f*"
+strong_convexity = 1 / huber.grad_lipschitz


Maybe our initial examples encouraged this, but I feel that this style is a bit too much "comment-code-staccato". I would prefer if there were some logical blocks like "space and data creation", "setup of the functionals and operators", etc.
Also, some of the comments only repeat what is already obvious from the code, like "Gradient operator", they can just be left out.

Agreed. In particular as this example is not very introductory.

kohr-h · 2017-10-19T13:41:35Z

+        self.obj_function_values = []
+
+
+callback = (odl.solvers.CallbackPrintIteration() & CallbackStore())


Using a step of 5 or 10 makes these kinds of examples run much faster. In 2D, plotting usually takes most of the runtime with step=1.

kohr-h · 2017-10-19T14:32:33Z

+# Assign operator and functionals
+op = gradient
+f = huber
+g = l1_norm


I would prefer using the names f and g immediately farther up instead of aliasing here.

kohr-h · 2017-10-19T14:33:07Z

+tau = 1.0 / gradient.norm  # Step size for primal variable
+sigma = 1.0 / gradient.norm  # Step size for dual variable
+
+# Run algorithms 2 and 3


Huh? I only see 1 algorithm being run.

kohr-h · 2017-10-19T16:25:46Z

+            i.ufuncs.logical_not(out=i)
+            out += i * (n - self.gamma / 2)
+        else:
+            out = n


This is a matter of taste, but a good principle for local variable names is "the larger the scope, the more verbose the name".
In practice that means that names like i and n are fine for short loops or list comprehensions. As soon as the scope spans more than, say, 5 lines of code, though, one-letter names make code less readable simply because they are harder to distinguish from numbers, e.g. i or l from 1.

Here, n could easily be norm without making lines annoyingly long. Also, the i multiplication can be replaced by boolean array indexing by going through Numpy:

norm = self.local_norm(x).asarray() with writable_array(out) as out_arr: if self.gamma > 0: # Quadratic part i = (norm < gamma) out_arr[i] = norm[i] ** 2 / (2 * self.gamma) # Absolute value part np.logical_not(i, out=i) out_arr[i] = norm[i] - self.gamma / 2

This is slightly slower due to boolean indexing, but more memory-friendly. You can scrape off even more copies by changing norm in-place before assigning to out, at the expense of an additional (slow) boolean indexing each time.

I agree with the variable names. When you write this code for yourself I guess one is just too lazy :)

About the computations, I am not so sure. I would like to keep the code numpy independent so that you can also use it on the gpu which I believe should be the case at the moment.

Fair enough. Just be aware that most of these ufuncs are not implemented in the CUDA backend and fall back to "copy to CPU, then use Numpy, then copy back to GPU". So you will probably end up being slower AND use more memory.

That said, Jonas is doing the (I hope) final round of review of #1088 which will make all this work much nicer.

kohr-h · 2017-10-19T16:30:59Z

+        f = FunctionalQuadraticPerturb(n.convex_conj,
+                                       quadratic_coeff=self.gamma / 2)
+
+        f.strong_convexity = 1 / self.grad_lipschitz


See above. Assigning attributes to foreign objects is not pretty (although I'm guilty of this, see #1177, but that's an exception 😇 ). We should aim for a solution that takes this as an optional initalization argument.

kohr-h · 2017-10-19T16:32:06Z

+    Parameters
+    ----------
+    space : `FnBase`
+        Space X which is the domain of the functional F


X and F don't refer to anything. I'd go for "Domain of the functional."

Too much copy-pasting :)

Yeah we sure have some cleanup to do in the older parts.

kohr-h · 2017-10-19T16:33:50Z

+    -----
+    The proximal operator is given by given by the proximal operator of
+    1/(2*gamma) * L2 norm in points that are <= gamma, and by the
+    proximal operator of the l1 norm in points that are > gamma.


Can you use backticks for things like gamma and formula-like stuff?

kohr-h · 2017-10-19T16:34:24Z

+            if isinstance(self.domain, ProductSpace):
+                return PointwiseNorm(self.domain, 2)(x)
+            else:
+                return x.ufuncs.absolute()


As above. Used only once, so just inline the code.

mehrhardt · 2017-10-20T14:17:15Z

All implemented. The only thing that is left is the actual computation of the function call, gradient, prox etc where I would like not to involve numpy. In the end @adler-j, @kohr-h you are the experts how to implement this most efficiently.

mehrhardt · 2017-10-20T14:23:10Z

OK, this looks weird. In the log of the files, these are not "new" commits.

What I did was to rebase the current branch onto upstream/master as

git fetch upstream
git rebase upstream/master

It then complained about the branches having diverged. I force-pushed my local version into the repo which "looks" good but github can't handle it (see fake commits above). What is the proper way to handle this situation?

kohr-h · 2017-10-20T14:29:53Z

It then complained about the branches having diverged. I force-pushed my local version into the repo which "looks" good but github can't handle it (see fake commits above). What is the proper way to handle this situation?

Which ones do you mean? Isn't 23fe237 the first commit in this PR?

mehrhardt · 2017-10-20T14:31:55Z

It is, but my github reminds me again that I "added some commits 9 days ago". Is that different for you?

kohr-h · 2017-10-20T14:33:58Z

It is, but my github reminds me again that I "added some commits 9 days ago". Is that different for you?

Ah, I think I know what you mean. When you rebase, GitHub adds the commits as new ones below the last comment. If you reload the page, the old ones go away.

kohr-h

I just have a few further comments and one hiccup with the norm thing, and I would like to settle this discussion about the two examples is a good way. Otherwise (without looking at all details) it looks good to me.

kohr-h · 2017-10-24T07:46:45Z

+        self.obj_function_values = []
+
+
+callback = (odl.solvers.CallbackPrintIteration(step=10) & CallbackStore())


Remove outer parens

I am OK with removing these but the two of you should settle on this. A while ago @adler-j told me to put parentheses for readability somewhere.

For stuff on one line? I'd be surprised. It makes sense when you have stuff that spans 2 or more lines.

kohr-h · 2017-10-24T07:48:50Z

+        min_{x >= 0}  1/2 ||x - d||_2^2
+            + lam * sum_i eta_gamma(||grad(x)_i||_2)
+
+where grad the spatial gradient and d is given noisy data. Here eta_gamma


Can you put the math stuff in double backticks, like ``grad``? Not that it's rendered anywhere, but it sticks a bit more out in the plain text. And who knows, maybe one day we can actually also render this text somewhere :-)

Also "is" missing after "grad"

kohr-h · 2017-10-24T08:00:05Z

@@ -0,0 +1,135 @@
+"""Total variation denoising using PDHG.


Another possibility is to make a big

if data_term == 'L1': # Do the L1 case, including plotting and whatnot elif data_term == 'L2': # Other case

Since it's a rather elaborate example anyway, making one large example from two slightly-less-than-large would still be good IMO.
What's your opinion @adler-j?

kohr-h · 2017-10-24T08:00:59Z

+
+# Define objective functional
+op = odl.Gradient(space)  # operator
+op.norm = np.sqrt(8) + 1e-4  # norm with forward differences is well-known


We shouldn't be setting attributes on instances like this. In fact, this overrides the method Operator.norm with a constant.
Until Gradient.norm is implemented please make this a free variable instead. Do we have an open issue on that somewhere? I can't find one.

OK, changed.

kohr-h · 2017-10-24T08:14:15Z

+        space : `FnBase`
+            Domain of the functional.
+        gamma : float
+            Smoothing parameter of Huberization. If ``gamma = 0``, then


kohr-h · 2017-10-24T08:14:24Z

+            Domain of the functional.
+        gamma : float
+            Smoothing parameter of Huberization. If ``gamma = 0``, then
+            functional is non-smooth corresponds to the usual L1 norm. For


kohr-h · 2017-10-24T08:15:05Z

+            Smoothing parameter of Huberization. If ``gamma = 0``, then
+            functional is non-smooth corresponds to the usual L1 norm. For
+            ``gamma > 0``, it has a ``1/gamma``-Lipschitz gradient so that
+            its convex conjugate is ``gamma``-strongly convex.


You surely prepare for this addition 😛

kohr-h · 2017-10-24T08:15:44Z

+        >>> abs(huber_norm(x) - l1_norm(x)) < tol
+        True
+
+        Redo previous example for a product space in two dimensions.


Colons at the end of those intermediate texts.

kohr-h · 2017-10-25T21:14:15Z

I made a small PR due to some things I noticed when doing the final skim. When that's in (you can also pick only the stuff that you like) I'll hit the merge button.

MAINT: small fixes to Huber doc and examples

mehrhardt · 2017-10-26T09:35:27Z

I liked all your changes @kohr-h so I merged it without cherry picking.

kohr-h · 2017-10-26T09:46:43Z

Thanks for the work @mehrhardt!

adler-j reviewed Oct 11, 2017

View reviewed changes

adler-j reviewed Oct 18, 2017

View reviewed changes

mehrhardt mentioned this pull request Oct 18, 2017

Issue 597 add the huber norm #1195

Closed

mehrhardt force-pushed the issue1101_pdhg_examples branch from 857b736 to 5b65eaf Compare October 19, 2017 09:38

kohr-h requested changes Oct 19, 2017

View reviewed changes

mehrhardt force-pushed the issue1101_pdhg_examples branch from 9b029d1 to 76be4e5 Compare October 20, 2017 14:19

kohr-h approved these changes Oct 24, 2017

View reviewed changes

mehrhardt added 14 commits October 24, 2017 09:41

BUG: remove bug in example

ceacc54

new examples for pdhg and huber functional

af5b002

STY: change linebreak for style

0f9c102

BUG: fix super issue flagged by travis

53fe427

STY: improve style of examples

6db63e4

STY: improve style of Huber functional

a606b87

BUG: Extend HuberL1L2 to non-product spaces

03a815a

STY: add solution to examples

fead912

Small changes based on Jonas' comments

7b7cfd3

manually merge this version of huber with Axels

2019b6f

STY: add blank line

d04365c

BUG: fix doc tests and make local_norm visible

ee952dc

Changes as requested by Holger

b723953

STY: improve style based on @kohr-h's comments

d7ba623

mehrhardt force-pushed the issue1101_pdhg_examples branch from 76be4e5 to d7ba623 Compare October 24, 2017 08:42

mehrhardt and others added 2 commits October 24, 2017 10:00

Change L1-Huber example

f27ae33

MAINT: small fixes to Huber doc and examples

744f175

Merge pull request #2 from kohr-h/issue1101_pdhg_examples

674cd96

MAINT: small fixes to Huber doc and examples

kohr-h merged commit 3934c20 into odlgroup:master Oct 26, 2017

mehrhardt deleted the issue1101_pdhg_examples branch October 26, 2017 09:47


		min_{x >= 0} \|\|x - d\|\|_1 + lam TV_gamma(x)

		where ``grad`` the spatial gradient and ``d`` is given noisy data.

		(1 / self.sigma) * self.functional.proximal(self.sigma))


		class HuberL1L2(Functional):

		@@ -0,0 +1,135 @@
		"""Total variation denoising using PDHG.

		This exhaustive example solve the L1-HuberTV problem

		self.obj_function_values = []


		callback = (odl.solvers.CallbackPrintIteration() & CallbackStore())

Uh oh!

Conversation

mehrhardt commented Oct 11, 2017

Uh oh!

pep8speaks commented Oct 11, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Comment last updated on October 26, 2017 at 09:26 Hours UTC

Uh oh!

adler-j left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mehrhardt commented Oct 18, 2017

Uh oh!

adler-j left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mehrhardt commented Oct 18, 2017

Uh oh!

kohr-h left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

pep8speaks commented Oct 11, 2017 •

edited

Loading