Proximal operator of Huber broken?

[mehrhardt]

I experienced some issue with the Huber functional. I guess this has something to do with the introduction of tensor spaces and indexing.

The code that I used worked fine for L1, L2Squared etc but broke for Huber. If needed, I can provide a minimal example but guess that this might not be necessary.

Do you know what is going wrong? Also, do you know why this has not been picked up by the tests?

Traceback (most recent call last):

  File "<ipython-input-62-2a507361e5bc>", line 2, in <module>
    pdhg_primal.run(200)

  File "/Users/me549/Documents/Nextcloud_personal/projects/201902_nonconvex_pdhg/python/misc.py", line 958, in run
    self.update()

  File "/Users/me549/Documents/Nextcloud_personal/projects/201902_nonconvex_pdhg/python/misc.py", line 987, in update
    self.prox_g(dual_tmp, out=b)

  File "/Users/me549/repositories/github_myodl/odl/operator/operator.py", line 678, in __call__
    result = self._call_in_place(x, out=out, **kwargs)

  File "/Users/me549/repositories/github_myodl/odl/operator/pspace_ops.py", line 302, in _call
    op(x[j], out=out[i])

  File "/Users/me549/repositories/github_myodl/odl/operator/operator.py", line 678, in __call__
    result = self._call_in_place(x, out=out, **kwargs)

  File "/Users/me549/repositories/github_myodl/odl/solvers/nonsmooth/proximal_operators.py", line 1921, in _call
    out[mask] = gamma / (gamma + self.sigma) * x[mask]

  File "/Users/me549/repositories/github_myodl/odl/space/pspace.py", line 947, in __getitem__
    raise TypeError('bad index type {}'.format(type(indices)))

TypeError: bad index type <class 'odl.discr.lp_discr.DiscreteLpElement'>

[adler-j]

Could you give us a minimal working example of this so we can reproduce the error?

[mehrhardt]

Here you go:

import odl

n = 20
X = odl.uniform_discr([0, 0], [1, 1], [n, n])

grad = odl.Gradient(X)
f = odl.solvers.Huber(grad.range, gamma=1)
f.proximal(1)(f.domain.one())

Traceback (most recent call last):

File "", line 8, in f.proximal(1)(f.domain.one())

File "/Users/me549/repositories/github_myodl/odl/operator/operator.py", line 686, in call out = self._call_out_of_place(x, **kwargs)

File "/Users/me549/repositories/github_myodl/odl/operator/operator.py", line 50, in _default_call_out_of_place result = op._call_in_place(x, out, **kwargs)

File "/Users/me549/repositories/github_myodl/odl/solvers/nonsmooth/proximal_operators.py", line 1921, in _call out[mask] = gamma / (gamma + self.sigma) * x[mask]

File "/Users/me549/repositories/github_myodl/odl/space/pspace.py", line 947, in getitem raise TypeError('bad index type {}'.format(type(indices)))

TypeError: bad index type <class 'odl.discr.lp_discr.DiscreteLpElement'>

[adler-j]

Well that was depressingly simple. I wonder why the tests have not caught this.

Edit: Does it work in 1d?

[mehrhardt]

No, same error for X = odl.uniform_discr([0, 0], [1, 1], [n, n]). For grad = Id there is no error.

EDIT: Meant to paste: X = odl.uniform_discr(0, 1, n).

[mehrhardt]

Well that was depressingly simple.

Does that mean you know what is going wrong and know how to fix it?

[aringh]

As the output indicates, the problem occurs on line 1921 in proximal_operators.py. In particular, it is x[mask] that causes the crash. In the 1D instance, x lives in ProductSpace(uniform_discr(0.0, 1.0, 20), 1) while mask lives in uniform_discr(0.0, 1.0, 20, dtype=bool). In the 1D case, one could argue that this should not give an error.

However, in the 2D case in the above example, x lives in ProductSpace(uniform_discr([ 0., 0.], [ 1., 1.], (20, 20)), 2) and mask lives in uniform_discr([ 0., 0.], [ 1., 1.], (20, 20), dtype=bool). In this case, I am not sure how x[mask] is supposed to behave...?

@mehrhardt when you say that

For grad = Id there is no error.

on what space is Id defined?

[aringh]

@mehrhardt, correct me if I am wrong, but what you want is a Huber-functional that works similar to the GroupL1Norm functional. Or? The proximal operator also seems to be implemented in a similar manner - otherwise I do not see why the PointwiseNorm is taken here.

[aringh]

I though the following would fix the problem

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

            mask = norm.ufuncs.less_equal(gamma + self.sigma)
            inv_mask = mask.ufuncs.logical_not()
            sign_x = x.ufuncs.sign()

            if isinstance(self.domain, ProductSpace):
                for (out_i, x_i, sign_xi) in zip(out, x, sign_x):
                    out_i[mask] = gamma / (gamma + self.sigma) * x_i[mask]
                    out_i[inv_mask] = x_i[inv_mask] - self.sigma * sign_xi[inv_mask]
            else:
                 out[mask] = gamma / (gamma + self.sigma) * x[mask]
                 out[inv_mask] = x[inv_mask] - self.sigma * sign_x[inv_mask]

            return out

However, when adding the product space to the slow test, it systemically fails the last assertion of test_proximal_defintion for the following two setups

[ functional='huber' - stepsize=10.0 - linear_offset=True - quadratic_offset=False - dual=False - product_space=True ]
[ functional='huber' - stepsize=10.0 - linear_offset=False - quadratic_offset=False - dual=False - product_space=True ]

Conclusion seems to be: something is wrong for stepsizes > 1 on product spaces.