Proximal of LInfty wrong?

[mehrhardt]

odl/odl/solvers/nonsmooth/proximal_operators.py

Line 1497 in a269486

radius = 1

As it is implemented, the proximal of the Linfty-norm does not depend on its scaling factor in front of it (see line above). Not sure what the correct formulation is but the one implemented has to be wrong. This is related to the missing tests as in #1385.

[mehrhardt]

My current guess is that it should be

        def _call(self, x, out):
            """Return ``self(x)``."""

            if x is out:
                x = x.copy()

            proj_l1(x, self.sigma, out)

but I didn't have time to verify it.

Anyone knows the answer?

[mehrhardt]

As it is implemented, the proximal of the Linfty-norm does not depend on its scaling factor in front of it (see line above). Not sure what the correct formulation is but the one implemented has to be wrong.

One way to see this is for sigma = 0 where one expects to get the identity as the prox.

[ozanoktem]

@JevgenijaAksjonova or @sbanert, do you have any input regarding the mathematical question that @mehrhardt posted?

[mehrhardt]

I think I figured what the prox should be:

        def _call(self, x, out):
            """Return ``self(x)``."""

            if x is out:
                x = x.copy()

            proj_l1(x, self.sigma, out)

            out.lincomb(-1, out, 1, x)

but is there a good way to test whether this is correct?

[mehrhardt]

It turns out that the issue is even deeper than I thought. The following example reveals that the projection onto L1 balls is incorrect.

import odl

X = odl.rn(2)

in0 = X.one().copy()

print(f(in0))
for sigma in [0, 1, 10]:

    f = odl.solvers.L1Norm(X)
    out = odl.solvers.nonsmooth.proj_l1(in0, sigma)
        
    print(sigma, f(out))

2.0
0 0.0
1 1.0
10 10.0

[mehrhardt]

OK, I believe I figured it out :D

The issue is the projection onto the L1-ball which should have been

    if out is None:
        out = x.space.element()

    u = x.ufuncs.absolute()
    
    if u.ufuncs.sum() <= radius:
        out[:] = x
    else:
        v = x.ufuncs.sign()
        proj_simplex(u, radius, out)
        out *= v

    return out

where in the original ODL code the "if" was missing. I'll make a PR.

[mehrhardt]

PR is now open #1614