Added Huber loss function class

Wrappers/Python/cil/optimisation/functions/HuberLoss.py

@@ -0,0 +1,186 @@
+#  Copyright 2026 United Kingdom Research and Innovation
+#  Copyright 2026 The University of Manchester
+#
+#  Licensed under the Apache License, Version 2.0 (the "License");
+#  you may not use this file except in compliance with the License.
+#  You may obtain a copy of the License at
+#
+#      http://www.apache.org/licenses/LICENSE-2.0
+#
+#  Unless required by applicable law or agreed to in writing, software
+#  distributed under the License is distributed on an "AS IS" BASIS,
+#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#  See the License for the specific language governing permissions and
+#  limitations under the License.
+#
+# Authors:
+# CIL Developers, listed at: https://github.com/TomographicImaging/CIL/blob/master/NOTICE.txt
+# Martin Sæbye Carøe (Technical University of Denmark, DTU Compute)
+
+import numpy as np
+import warnings
+from numbers import Number
+
+from cil.optimisation.functions import Function
+from cil.optimisation.operators import DiagonalOperator, LinearOperator
+from cil.framework import DataContainer
+
+class HuberLoss(Function):
+    r"""
+    (Weighted) Huber loss
+
+    For residual :math:`r = Ax - b`:
+
+    .. math::
+        \phi_\delta(r) =
+        \begin{cases}
+            0.5 * r^2 & \text{if } |r| \leq \delta \\
+            \delta * (|r| - 0.5*\delta) & \text{otherwise}
+        \end{cases}
+
+    Parameters
+    ----------
+    A : LinearOperator
+    b : Data, DataContainer
+    huber_delta : float
+        Transition point between L2 and L1 behaviour
+    c : float, default 1.0
+        Scaling constant
+    weight : DataContainer, optional
+        Positive diagonal weights
+    """
+
+    def __init__(self, A, b, huber_delta, c=1.0, weight=None):
+        super(HuberLoss, self).__init__()
+
+        if huber_delta <= 0:
+            raise ValueError("huber_delta must be positive")
+
+        self.A = A
+        self.b = b
+        self.c = c
+        self.huber_delta = huber_delta
+
+        self.weight = weight
+        self._weight_norm = None
+
+        if weight is not None:
+            if (self.weight < 0).any():
+                raise ValueError("Weight contains negative values")
+
+    def __call__(self, x):
+
+        r = self.A.direct(x)
+        r.subtract(self.b, out=r)
+
+        abs_r = r.abs()
+
+        # m = min(|r|, delta)
+        m = abs_r.copy()
+        m.minimum(self.huber_delta, out=m)
+
+        # 0.5 * m^2
+        val = m.power(2)
+        val.multiply(0.5, out=val)
+
+        # delta * (|r| - m)
+        lin = abs_r.copy()
+        lin.subtract(m, out=lin)
+        lin.multiply(self.huber_delta, out=lin)
+
+        val.add(lin, out=val)
+
+        if self.weight is not None:
+            val.multiply(self.weight, out=val)
+
+        return self.c * val.sum()
+
+
+
+    def gradient(self, x, out=None):
+
+        if out is None:
+            out = x * 0.0
+
+        r = self.A.direct(x)
+        r.subtract(self.b, out=r)
+
+        abs_r = r.abs()
+
+        # m = min(|r|, delta)
+        m = abs_r.copy()
+        m.minimum(self.huber_delta, out=m)
+
+        # grad wrt residual: sign(r) * m
+        grad_r = r.sign()
+        grad_r.multiply(m, out=grad_r)
+
+        if self.weight is not None:
+            grad_r.multiply(self.weight, out=grad_r)
+
+        self.A.adjoint(grad_r, out=out)
+        out.multiply(self.c, out=out)
+
+        return out
+
+
+
+    @property
+    def L(self):
+        if self._L is None:
+            self.calculate_Lipschitz()
+        return self._L
+
+    @L.setter
+    def L(self, value):
+        warnings.warn("You should set the Lipschitz constant with calculate_Lipschitz().")
+        if isinstance(value, Number) and value >= 0:
+            self._L = value
+        else:
+            raise TypeError("The Lipschitz constant must be non-negative")
+
+    def calculate_Lipschitz(self):
+        """
+        Lipschitz constant of gradient.
+
+        For Huber:
+            .. math:: \max \phi'' = 1
+        so:
+            .. math:: L = c * ||A||^2
+        (weighted: multiplied by :math:`||W||`)
+        """
+        try:
+            self._L = np.abs(self.c) * (self.A.norm() ** 2)
+        except AttributeError:
+            if self.A.is_linear():
+                Anorm = LinearOperator.PowerMethod(self.A, 10)[0]
+                self._L = np.abs(self.c) * (Anorm * Anorm)
+            else:
+                warnings.warn(
+                    f"{self.__class__.__name__} could not calculate Lipschitz Constant."
+                )
+
+        if self.weight is not None:
+            self._L *= self.weight_norm
+
+    @property
+    def weight_norm(self):
+        if self.weight is not None:
+            if self._weight_norm is None:
+                D = DiagonalOperator(self.weight)
+                self._weight_norm = D.norm()
+        else:
+            self._weight_norm = 1.0
+        return self._weight_norm
+
+    def __rmul__(self, other):
+        if not isinstance(other, Number):
+            raise NotImplemented
+
+        return HuberLoss(
+            A=self.A,
+            b=self.b,
+            huber_delta=self.huber_delta,
+            c=self.c * other,
+            weight=self.weight
+        )

Wrappers/Python/cil/optimisation/functions/__init__.py

@@ -40,4 +40,5 @@
 from .SVRGFunction import SVRGFunction, LSVRGFunction
 from .SAGFunction import SAGFunction, SAGAFunction
 from .AbsFunction import FunctionOfAbs
+from .HuberLoss import HuberLoss
 

docs/source/optimisation.rst

@@ -535,6 +535,12 @@ Mixed L11 norm
    :members:
    :inherited-members:
 
+Huber Loss
+---------------
+.. autoclass:: cil.optimisation.functions.HuberLoss
+   :members:
+   :inherited-members:
+
 Total variation
 ---------------