Huber.prox uses the wrong branch threshold (|x| > alpha instead of |x| > alpha + tau) #275

pyproximal.Huber.prox decides which branch to apply from |x| > alpha, but for the proximal operator the switch happens at |x| = alpha + tau. Inputs with alpha < |x| <= alpha + tau are soft-thresholded when they should be scaled. The result is not the minimizer, and the prox becomes discontinuous and non-monotone in x. This is in pyproximal/proximal/Huber.py (line 63 on main, also in 0.12.0). Related: #189, which split Huber and HuberCircular but left this threshold as it was.

Reproducer

import numpy as np
import pyproximal

alpha, tau = 1.0, 1.0
h = pyproximal.Huber(alpha=alpha)
x = np.array([-2.5, -1.5, -0.5, 0.5, 1.0, 1.2, 1.5, 1.9, 2.0, 2.5, 4.0])

# brute-force argmin_y H(y) + (y - x)^2 / (2 tau)
ys = np.linspace(-5, 5, 1000001)
H = lambda t: np.where(np.abs(t) <= alpha, t**2 / (2 * alpha), np.abs(t) - alpha / 2)
brute = np.array([ys[np.argmin(H(ys) + (ys - xi) ** 2 / (2 * tau))] for xi in x])

print(np.round(h.prox(x, tau), 4))  # [-1.5 -0.5  -0.25 0.25 0.5 0.2 0.5  0.9  1. 1.5 3.]
print(np.round(brute, 4))           # [-1.5 -0.75 -0.25 0.25 0.5 0.6 0.75 0.95 1. 1.5 3.]

For example, x = 1.2 returns 0.2, which is smaller than the 0.5 returned for x = 1.0.

Fix

The prox of tau * H_alpha is

prox(x) = x * alpha / (alpha + tau)    if |x| <= alpha + tau
          x - tau * sign(x)            otherwise

which can also be written in closed form as (1 - tau / max(|x|, alpha + tau)) * x. That is the alternative already in the commented-out code; the comment says it "does not provide the same result", and that is because the current code is the one that is wrong. A minimal change to prox:

mask = xabs > self.alpha + tau

Equivalently, replace the body with y = (1.0 - tau / np.maximum(np.abs(x), tau + self.alpha)) * x. With either change the output matches the brute-force minimizer at all the points above. A test against a brute-force or scipy.optimize minimizer on a grid of x values that includes alpha < |x| < alpha + tau would catch a regression.

Related suggestion: HuberCircular applies the Huber function to the norm of the whole vector. For isotropic Huber-TV you need it per group, the way L21 works, e.g. a HuberL21(ndim, sigma, alpha) whose prox is v_i * (1 - tau*sigma / max(||v_i||, tau*sigma + alpha)) for each pixel.