ENH: add decorator to fix adjoint weightings by kohr-h · Pull Request #1177 · odlgroup/odl

kohr-h · 2017-10-05T12:11:36Z

Update: The general logic is settled, details are in discussion.

~~This is a discussion PR.~~ I think we could go somewhere along these lines to make adjoints "automatically correct" with respect to weighted spaces, in particular when the weights differ.

This solution should always be correct where it applies, however it wouldn't yet help in special cases like RayTransform where the weights need to be computed per axis.
The larger question is what to do instead of the monkey-patching here. I can imagine a bunch of scenarios:

Do this automatically in Operator. This means that users implement the unweighted adjoint, and the weighted one is implemented with a helper function similar to this.
This would require us to communicate this "magic" clearly and prominently.
As in 1., but depending on some flag or whatever that is provided by the user. It could be a class attribute like _auto_weighting or so. Don't do anything by default.
This option would still need good communication but leads to less unexpected behavior on the user side.
Just offer a tool (e.g. a decorator) for adjoint, like @auto_weighting that users can choose to use.
This option has the least surprise and the simplest migration path, we would just have to use it a lot ourselves so users pick it up.

Opinions?

For the records: What I don't want is to let adjoint return an operator composition. That would be ugly and a usability nightmare, since operator properties will be buried deep down.

TODO:

What to do if one of the spaces (or both) don't define weighting at all?
We could treat no weighting as 1.0, but need to be careful with implementation since the current variant uses in-place arithmetic. It should be fine, though, since the weighting factors that are 1.0 are just skipped.
Make decorator for the (exact) Operator.norm method (should be way simpler though)
Fix existing operators.
- ~~Add domain or range where missing~~ not in this PR
- Fix operator impls for changes with differing domain/range
Add unit tests with domain and range differing in weighting
Document the decorator in the online doc.
Think about streamlining certain types, e.g. OperatorComp to avoid dumb "multiplications by 1"
There seems to be an issue when adjoint uses another property, e.g. inverse, needs investigation
Remove the repr changes that are now in DONTMERGE: improve repr by a lot #1276

adler-j · 2017-10-05T12:20:07Z

In my opinion, adding a utility is the only option that is not outrageously surprising to users. This could be, as you suggest, in the form of a decorator, but I'd simply make it a free function that is either called on the adjoint:

@property
def adjoint(self):
    ...
    return auto_weighting(adjoint)

or in the operator itself

class AdjointyOperator:
    def __call__(self, x):
        return auto_weighting(adjoint)

kohr-h · 2017-10-05T12:58:52Z

In my opinion, adding a utility is the only option that is not outrageously surprising to users.

I think the same. Only for the sake of not introducing bias in the discussion I left it kind of open.

Note that the option of decorating the adjoint, which is equivalent to your first code box, would still require some monkey-patching. However, in this case I would find it more acceptable since it happens "almost" during object creation.

The second option is cleaner in that respect, but has the disadvantage that it would make it harder to use existing (unweighted) operators off the shelf. Things like return ThisOtherOperator() in adjoint will not be "auto-weightable" since that would require a _call method to decorate.

So I'm leaning towards the operator decoration option.

And BTW, this has interesting implications for OperatorComp where you could start pushing around weighting constants all the way to the smallest space, for instance.

I'll think a bit more about this.

kohr-h · 2017-10-05T13:14:35Z

, which is equivalent to your first code box,

roughly :-)

kohr-h · 2017-10-05T14:22:03Z

Here's a more serious suggestion. Obviously needs more tests, but does the job.

kohr-h · 2017-10-05T14:33:10Z

Also obviously, the implementation does not make use of the optimize flag yet and just does it unconditionally. Further, it does it badly by computing those 1 / something arrays. That needs to change.

pep8speaks · 2017-10-05T15:56:33Z

Checking updated PR...

In the file odl/discr/discr_ops.py, following are the PEP8 issues :

Line 147:1: E303 too many blank lines (3)

Complete extra results for this file :

class ResizingOperator(Operator):^---

In the file odl/operator/tensor_ops.py, following are the PEP8 issues :

Line 517:80: E501 line too long (80 > 79 characters)

Complete extra results for this file :

    ProductSpace(uniform_discr([-1., -1.], [ 1.,  1.], (1, 2)), 2).element([                                                                               ^---

Comment last updated on February 15, 2018 at 22:06 Hours UTC

kohr-h · 2017-10-05T16:02:36Z

So this is now a serious PR. I've added tests and an example, and I think this concept would actually work. Comments welcome.

kohr-h · 2017-10-05T16:05:26Z

Once we settle on a solution, I think this PR should also go through the "standard library" and fix all instances of operators that need fixing. I'll make a TODO list at the top.

kohr-h · 2017-10-05T16:10:50Z

Open question: what to do with spaces that have no weighting attribute? (See TODOs)

kohr-h · 2017-10-05T16:31:10Z

Ping @matthiasje this may be interesting to you, too. You also had trouble with unscaled adjoints, right?

mehrhardt · 2017-10-06T09:02:56Z

Yes, in the most elegant definition of TGV I would like to use weighted spaces for the second derivative. It then needs proper adjoints and prox operators. This should be used in
examples/solvers/pdhg_denoising_tgv.py. What would be an easy way for me to test this? Can I checkout your branch?

kohr-h · 2017-10-06T09:06:27Z

Sure, you can do

git remote add kohr-h https://github.com/odlgroup/odl.git
git fetch kohr-h
git checkout issue-1068__ajoint_weighting  # sorry about the typo :-)

mehrhardt · 2017-10-12T12:43:02Z

@kohr-h, I tried the TGV example again with weighted spaces. It does not work as weighting is not implemented for product spaces?
NotImplementedError: weighted spaces not supported

kohr-h · 2017-10-12T12:48:06Z

What happens if you remove that check and add the @auto_weighting decorator to ProductSpaceOperator.adjoint?

mehrhardt · 2017-10-18T15:33:13Z

I tried the following but it did not help significantly. It does then result in other errors about things not being well-defined (see below).

+++ b/odl/operator/pspace_ops.py
@@ -18,6 +18,7 @@ from odl.operator.operator import Operator
 from odl.operator.default_ops import ZeroOperator
 from odl.space import ProductSpace
 
+from odl.space.space_utils import auto_weighting
 
 __all__ = ('ProductSpaceOperator',
            'ComponentProjection', 'ComponentProjectionAdjoint',
@@ -148,15 +149,11 @@ class ProductSpaceOperator(Operator):
             if not isinstance(domain, ProductSpace):
                 raise TypeError('`domain` {!r} not a ProductSpace instance'
                                 ''.format(domain))
-            if domain.is_weighted:
-                raise NotImplementedError('weighted spaces not supported')
 
         if range is not None:
             if not isinstance(range, ProductSpace):
                 raise TypeError('`range` {!r} not a ProductSpace instance'
                                 ''.format(range))
-            if range.is_weighted:
-                raise NotImplementedError('weighted spaces not supported')
 
         # Convert ops to sparse representation
         self.ops = scipy.sparse.coo_matrix(operators)
@@ -312,6 +309,7 @@ class ProductSpaceOperator(Operator):
         return ProductSpaceOperator(deriv_matrix, self.domain, self.range)
 
     @property
+    @auto_weighting
     def adjoint(self):
         """Adjoint of this operator.

normE = E.norm(estimate=True)
normE2 = E2.norm(estimate=True)


Traceback (most recent call last):

  File "<ipython-input-2-b66a9577355f>", line 1, in <module>
    normE = E.norm(estimate=True)

  File "/home/ehrhardt/repositories/git_myODL/odl/operator/operator.py", line 751, in norm
    return power_method_opnorm(self, **kwargs)

  File "/home/ehrhardt/repositories/git_myODL/odl/operator/oputils.py", line 229, in power_method_opnorm
    op.adjoint(tmp, out=x)

  File "/home/ehrhardt/repositories/git_myODL/odl/operator/operator.py", line 686, in __call__
    result = self._call_in_place(x, out=out, **kwargs)

  File "/home/ehrhardt/repositories/git_myODL/odl/space/space_utils.py", line 421, in _call
    x = new_ran_w * x

TypeError: unsupported operand type(s) for *: 'numpy.ndarray' and 'ProductSpaceElement'

mehrhardt · 2017-10-18T15:34:58Z

Maybe it would be the best if you have a look yourself. I tried to run pdhg_denoising_tgv.py with the following modification after the definition of derivates:

# Create symmetrized operator and weighted space.
# TODO: As the weighted space is currently not supported in ODL we find a
# workaround.
W = odl.ProductSpace(U, 3, weighting=[1, 1, 2])
E = odl.operator.ProductSpaceOperator(
    [[Dx, 0], [0, Dy], [0.5*Dy, 0.5*Dx]], range=W)

E2 = odl.operator.ProductSpaceOperator(
    [[Dx, 0], [0, Dy], [0.5 * Dy, 0.5 * Dx], [0.5 * Dy, 0.5 * Dx]])
W2 = E.range

normE = E.norm(estimate=True)
normE2 = E2.norm(estimate=True)

kohr-h · 2017-10-18T21:27:00Z

Without actually looking, I believe this has to do with #1062 in some sense. What happens here is that the code tries to multiply the result, a product space element, with a Numpy array that has the weights, one for each product space component. That doesn't work, unfortunately, which is part of a bigger issue.

~~I can't find the discussion we had before, but let me make a suggestion how to solve this and #1062, and probably some other issues, too. I would suggest this:~~

>>> pspace = odl.rn(2) * odl.rn(3)
>>> z = pspace.one()
>>> z * [1, 2]  # should work, doesn't currently
ProductSpace(rn(2), rn(3)).element([
    [1.0, 1.0],
    [2.0, 2.0, 2.0]
])
>>> z * np.array([1, 2])  # same as above
ProductSpace(rn(2), rn(3)).element([
    [1.0, 1.0],
    [2.0, 2.0, 2.0]
])
>>> pspace = odl.rn(2) ** 2
>>> z = pspace.one()
>>> z * np.array([1, 2])[None, :]  # broadcast along pspace axis, now produces Numpy array
ProductSpace(rn(2), rn(3)).element([
    [1.0, 2.0],
    [1.0, 2.0]
])

~~The principle would be as follows:~~

~~For non-power product spaces, only 1D arrays of the same length as the space make sense. So the only sensible thing is to do the equivalent of z[i] * arr[i] for i in range(length).~~
~~Since power spaces are just a special case, the same thing must work for them.~~
~~To broadcast in a different manner, one has to explicitly add empty axes into the array, and the same broadcasting rules as in Numpy apply.~~

~~That would make it unambiguous as to how to interpret~~

>>> pspace = odl.rn(2) ** 2
>>> z = pspace.one()
>>> z * [1, 2]  # along axis 0 or axis 1??

It would always be axis 0, and for axis 1 one would have to do z * [[1, 2]] or with Numpy, z * np.array([1, 2])[None, :]. And it would solve the current issue that the latter variant produces a Numpy array, which I think is a bug.

kohr-h · 2017-10-23T08:20:28Z

In view of the current issue with product spaces, I think I'll implement a workaround for that case to move things forward.

kohr-h · 2017-11-11T15:20:39Z

In view of the current issue with product spaces, I think I'll implement a workaround for that case to move things forward.

Will do that now.

kohr-h · 2017-11-11T21:50:40Z

Now the decorator works also for operators on product spaces. Needs a proper review now. @mehrhardt if you like, you can test your code again.

Anyway, I think this is a rather important issue since many adjoints are wrong as of now, we should give it some priority.

mehrhardt · 2017-11-13T11:12:52Z

I have almost good news. After removing a few lines in pspace_ops.py

        # Validate input data
        if domain is not None:
            if not isinstance(domain, ProductSpace):
                raise TypeError('`domain` {!r} not a ProductSpace instance'
                                ''.format(domain))
#            if domain.is_weighted:
#                raise NotImplementedError('weighted spaces not supported')

        if range is not None:
            if not isinstance(range, ProductSpace):
                raise TypeError('`range` {!r} not a ProductSpace instance'
                                ''.format(range))
#            if range.is_weighted:
#                raise NotImplementedError('weighted spaces not supported')

the following example

# two models of the 'same' operator

import odl

U = odl.uniform_discr(min_pt=[-20, -20], max_pt=[20, 20], shape=[300, 300],
                      dtype='float32')
G = odl.Gradient(U, method='forward', pad_mode='symmetric')
V = G.range

Dx = odl.PartialDerivative(U, 0, method='backward', pad_mode='symmetric')
Dy = odl.PartialDerivative(U, 1, method='backward', pad_mode='symmetric')

E = odl.operator.ProductSpaceOperator(
    [[Dx, 0], [0, Dy], [0.5 * Dy, 0.5 * Dx], [0.5 * Dy, 0.5 * Dx]])
W = E.range

W2 = odl.ProductSpace(U, 3, weighting=[1, 1, 2])
E2 = odl.operator.ProductSpaceOperator(
        [[Dx, 0], [0, Dy], [0.5*Dy, 0.5*Dx]], range=W2)

l1 = odl.solvers.L1Norm(W)
l12 = odl.solvers.L1Norm(W2)


w = W.one()
w2 = W2.one()

v = odl.phantom.white_noise(V)

print(w2.inner(E2(v)))
print(w.inner(E(v)))

print(l1(w))
print(l12(w2))

print(E.norm(estimate=True, xstart=v))
print(E2.norm(estimate=True, xstart=v))

gives the output:

10.128983855247498
10.128984451293945
6400.0
6400.0
21.1059752145
18.2776606142

Only the norm estimate is wrong. Maybe you have an idea what is wrong here? Also it might be the easiest if you commit the small changes in pspace_ops.py yourself rather me making a PR.

kohr-h · 2017-12-01T16:37:17Z

OK I've addressed all the comments I think. The only TODO (if the design is ok) is to properly document it in the online doc.

kohr-h · 2017-12-16T11:15:44Z

Bump for review.

See also the weight handling in #1238 with per-axis weighting. The weight merging there is done in the adjoint_weightings function.

adler-j

Some overall questions:

Does this commute with @property?

How does this handle adjoints that are defined in terms of other operators?

Shouldn't this be a operator util instead of a space util?

adler-j · 2017-12-16T23:33:49Z

@@ -0,0 +1,58 @@
+"""Example demonstrating the usage of the ``auto_weighting`` decorator."""


This needs to be somewhat fleshed out, especially given that new users come here

adler-j · 2017-12-16T23:39:21Z

    return rn


+class auto_weighting(OptionalArgDecorator):


Suggest rename to auto_adjoint_weighting to make this more clear

adler-j · 2017-12-16T23:40:36Z

+    ``M: R^n -> R^m``, and its unweighted adjoint is defined by the
+    transposed matrix ``M^T: R^m -> R^n``.
+
+    Similarly, when an integration operator ``A: L^2(I) -> R`` on an


Wouldnt always using the R^n -> R style be better?

What do you mean by "using"? In this example?

adler-j · 2017-12-16T23:42:39Z

+
+    Parameters
+    ----------
+    unweighted_adjoint : `Operator`


Shouldn't the type of this be "method" or something like that?

Technically yes. I'll change it.

adler-j · 2017-12-16T23:50:40Z

+        elif dom_w_type == 'const' and ran_w_type == 'const':
+            if unweighted_adjoint.domain.size < unweighted_adjoint.range.size:
+                new_dom_w = 1.0
+                new_ran_w = ran_w / dom_w


Are we sure this always makes sense? What about operators C -> R etc?

Yes, because weightings can never be complex. (Note to self: make sure that this is checked.)

adler-j

Some further reviewing and a question.

Could this somehow become buggy if adjoint is an attribute instead of a method?

adler-j · 2017-12-17T00:01:25Z

+            return ScalingOp(self.range, self.domain, self.c)
+
+    op1 = ScalingOp(rn, rn, 1.5)
+    op2 = ScalingOp(rn_w, rn_w, 1.5)


Prefer not using 1.5 everywhere

adler-j · 2017-12-17T00:01:36Z

+def test_auto_weighting_cached_adjoint():
+    """Check if auto_weighting plays well with adjoint caching."""
+    rn = odl.rn(2)
+    rn_w = odl.rn(2, weighting=2)


Use a "random" value here

adler-j · 2017-12-17T00:02:10Z

+        def __init__(self, dom, ran, c):
+            super(InvalidScalingOp, self).__init__(dom, ran, linear=True)
+            self.c = c
+            self._adjoint = None


What does this line do?

That must be copy-paste from the caching test, should be removed.

kohr-h · 2017-12-17T00:07:15Z

Could this somehow become buggy if adjoint is an attribute instead of a method?

That cannot happen since adjoint is a property on Operator and cannot be overwritten by an attribute:

>>> class A(odl.Operator):
...     def __init__(self):
...         super(A, self).__init__(odl.rn(2), odl.rn(2))
...         self.adjoint = None
>>>
>>> a = A()
Traceback (most recent call last):
    [...]

    self.adjoint = None

AttributeError: can't set attribute

kohr-h · 2017-12-17T00:54:11Z

Does this commute with @Property?

No. I tried to hack it in, but the wrapping does not seem to work as expected. I'll leave it as it is, but I can improve the error message to give users a hint.

kohr-h · 2017-12-17T00:54:59Z

How does this handle adjoints that are defined in terms of other operators?

Not sure what you mean. Like compositions and such?

kohr-h · 2017-12-17T00:55:17Z

Shouldn't this be a operator util instead of a space util?

Yes I guess so.

kohr-h · 2018-08-27T08:47:59Z

Current status here is that there's some mysterious issue when adjoint defers to inverse.

I'll think a bit about better ways to do this, without doing brain surgery on the operator returned by adjoint. Maybe a metaclass after all.

kohr-h force-pushed the issue-1068__ajoint_weighting branch from 9dbb725 to 2a92199 Compare October 5, 2017 14:21

kohr-h force-pushed the issue-1068__ajoint_weighting branch 2 times, most recently from e25fa7d to c211b66 Compare October 5, 2017 15:56

kohr-h force-pushed the issue-1068__ajoint_weighting branch 2 times, most recently from 64e16b2 to 85c64fb Compare October 5, 2017 16:01

kohr-h changed the title ~~WIP: fix adjoint weightings~~ ENH: add decorator to fix adjoint weightings Oct 5, 2017

kohr-h added the status: review needed label Oct 5, 2017

kohr-h mentioned this pull request Oct 19, 2017

Issue1101 pdhg examples #1191

Merged

kohr-h force-pushed the issue-1068__ajoint_weighting branch from 85c64fb to e7c7da4 Compare November 11, 2017 15:20

adler-j approved these changes Dec 16, 2017

View reviewed changes

adler-j reviewed Dec 17, 2017

View reviewed changes

kohr-h force-pushed the issue-1068__ajoint_weighting branch 2 times, most recently from 3fe7cec to 962a8bc Compare December 24, 2017 11:14

kohr-h removed the status: review needed label Feb 14, 2018

Holger Kohr and others added 15 commits February 15, 2018 22:10

ENH: add auto_weighting decorator for adjoints

174e547

ENH: add auto_weighting example

969a7d6

TST: add tests for auto_weighting decorator

515c3c5

MAINT: minor changes and comments in pspace

692bcfa

ENH: implement auto-weighting on product spaces

92c96c4

TST: add tests for product space auto weighting

4016839

MAINT: simplify auto-weighting example

c161b66

MAINT: move auto-weighting example to operator examples

2013ec4

DOC: improve doc of auto_weighting decorator

a706037

ENH: make auto_weighting work with nested product spaces

651f16c

DOC: add example to auto-weighting

7c9c252

MAINT: rename auto[_adjoint]_weighting decorator and move to oputils

af58e3d

MAINT: move adjoint weighting tests to oputils_test

6d57cad

TST: add test for adjoint deco with adjoint from other class

8ac484d

WIP: port existing operators to use auto adjoint weighting

9a0a073

kohr-h force-pushed the issue-1068__ajoint_weighting branch from f67f8cf to 9a0a073 Compare February 15, 2018 22:06

kohr-h mentioned this pull request Jun 30, 2018

ENH: Tensor-valued DiscreteLp #1238

Open

21 tasks

kohr-h mentioned this pull request Aug 27, 2018

Op repr updates #1394

Open

3 tasks

		@@ -0,0 +1,58 @@
		"""Example demonstrating the usage of the ``auto_weighting`` decorator."""

Uh oh!

Conversation

kohr-h commented Oct 5, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

adler-j commented Oct 5, 2017

Uh oh!

kohr-h commented Oct 5, 2017

Uh oh!

kohr-h commented Oct 5, 2017

Uh oh!

kohr-h commented Oct 5, 2017

Uh oh!

kohr-h commented Oct 5, 2017

Uh oh!

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

Uh oh!

Comment last updated on February 15, 2018 at 22:06 Hours UTC

Uh oh!

kohr-h commented Oct 5, 2017

Uh oh!

kohr-h commented Oct 5, 2017

Uh oh!

kohr-h commented Oct 5, 2017

Uh oh!

kohr-h commented Oct 5, 2017

Uh oh!

mehrhardt commented Oct 6, 2017

Uh oh!

kohr-h commented Oct 6, 2017

Uh oh!

mehrhardt commented Oct 12, 2017

Uh oh!

kohr-h commented Oct 12, 2017

Uh oh!

mehrhardt commented Oct 18, 2017

Uh oh!

mehrhardt commented Oct 18, 2017

Uh oh!

kohr-h commented Oct 18, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

kohr-h commented Oct 23, 2017

Uh oh!

kohr-h commented Nov 11, 2017

Uh oh!

kohr-h commented Nov 11, 2017

Uh oh!

mehrhardt commented Nov 13, 2017

Uh oh!

kohr-h commented Dec 1, 2017

Uh oh!

kohr-h commented Dec 16, 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!

adler-j left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

kohr-h commented Oct 5, 2017 •

edited

Loading

pep8speaks commented Oct 5, 2017 •

edited

Loading

kohr-h commented Oct 18, 2017 •

edited

Loading