CGLS is similar to CG. Mathematically, is it equivalent to solve a
problem Ax = b, via the normal equation A^TAx = A^Tb. However, it avoid
explicitly creating the operator A^TA, and hence is both
numerically more stable, and more efficient.

The previous version of ISTA, only solved least squares problem with L1
regularization. However, the algorithm is capable of handling more
general problems. This is implemented now, and to reflect that change,
the name is changed from ISTA to PGD (Proximal gradient descent).

Previously, the names reflected more the implementation. Now they show
the correspondance with their functional. I.e. the proximal of the L1
norm, is now `ProximalL1` instead of `SoftThresholding`.

Type erasure is a method to achieve runtime polymorphism, similar to
classical C++ inheritance. However, it has a couple of nice aspects.

The main difference is, that any class can "bind" to the base class, as
long as it adheres to the interface. This mostly means, that as long as
a class doesn't bind to the base class, it's an independend class.
Without the cost of virtual functions, complicated cloning and such and
such. Hence, it improves local reasoning of types. You only need to know
about the interface, when you want to "bind" against it, e.g. as a
function parameter.

The one major drawback of type erasure is the rather complicated
implementation of the base class. However, this only needs to be done
for interfaces, and so it's rather rare.

This certainly is an experiment, however, I'm quite confident, that it
will improve our codebase overall (due to an increased local reasoning,
most types are just regular and such).

Proximal operators are currently still rather small, and hence it is a
good playground without ripping apart the whole framework. We can try it
here, and see if we want to further use it.

NumPy's minimum and maximum functions take an array and a scalar,
and return a new vector for which each element is either the minimum or
maximum of the array element and the provided scalar.

This is usefull for e.g. non-negativity constraints or upper bounds.

David Alkier authored Dec 27, 2022 and David Frank committed Jan 09, 2023

The rotation matrix was set up the wrong way.
It was transposed and did not work as expected.
The old shepp logan approach deals with the
wrong rotation matrix by flipping on the y axis.
Therefore `example3d,py` creates a mirror
image when comparing old and new.