Added Huber loss function class

[msaca-okse]

Description

The Huber loss function behaves like the L2 loss close to 0, and like L1 loss for large values. The loss function is differentiable and can thus be used in GD or FISTA. It's an alternative to L2 loss that is more robust to extreme values. The implementation mirrors that of "cil.optimisation.functions.LeastSquares". The Huber loss has a parameter that selects when the cutoff between L2 and L1 loss happens.

Example Usage

G = IndicatorBox(lower=0.0)
huber_delta = 1.0
F = HuberLoss(A,b,huber_delta)
myFISTA = FISTA(f=F, 
                  g=G, 
                  initial=x0)

Contribution Notes

• The content of this Pull Request (the Contribution) is intentionally submitted for inclusion in CIL (the Work) under the terms and conditions of the Apache-2.0 License
• I confirm that the contribution does not violate any intellectual property rights of third parties

[msaca-okse]

Fixes #2281

[lauramurgatroyd]

Hi @msaca-okse I just added the license header and made some edits to ensure the documentation renders well:

I wonder whether c should be included in the above equation?

We will need to add unit tests before merging into CIL. Would you be interested in working on this or are you happy for the CIL developer team to add these?

[lauramurgatroyd]

Hi @msaca-okse, thanks for your contribution. I am now working on getting this in CIL. I am about to push some unit tests which I have developed.

[lauramurgatroyd]

I've updated the equation, rendered it looks like this:

[MargaretDuff]

A really picky review - apologies.

However, I am happy with the code, the gradient and Lipscitz calculations!

[MargaretDuff]

I know that this is copied from Least squares but for a linear operator, calling A.norm() either accesses a stashed value or calls the power method itself. I don't understand in what cases the if statement on 157 will be triggered.

[MargaretDuff]

Similarly, not sure what this else will pick up. If the operator is not linear and calculate_norm is not defined, the user will get a NotImplementedError. Perhaps we catch and return that with a bit more of an explanation?

[MargaretDuff]

This would do the same thing without building a new operator - I see arguments both ways

Suggested change

	D = DiagonalOperator(self.weight)
	self._weight_norm = D.norm()
	self._weight_norm = self.weight.abs().max()