ENH/DOC/TST: add FunctionalRightVectorMult and in-dept doc for functionals.

aringh · aringh · commit f22e76288dfc · 2016-07-25T18:21:02.000+02:00
FunctionalRightVectorMult was implemeted before via OperatorRightVectorMult.
However, this is an operator and not a functional, why
FunctionalRightVectorMult has been explicitly added.

The commit also add tests for FunctionalRightVectorMult.

The commit also contains a (start for) an in-dept doc for functionals.
diff --git a/doc/source/guide/in_depth/functional_guide.rst b/doc/source/guide/in_depth/functional_guide.rst
@@ -0,0 +1,149 @@
+.. _functional_in_depth:
+
+#####################
+Functional
+#####################
+
+A *functional* is an operator ``f`` that maps from some vector space ``X`` to the field of scalars ``F`` associated with the vector space:
+
+.. math::
+   
+   f : X \to F.
+
+In the ODL solver package, functionals are implemented in the `Functional` class, a subclass to `Operator`.
+
+From a mathematical presepective, the above is a valide definition of a functional. However, since the purpose of these functionals are primarely to be used for solving optimization problems, the following assumptions are made:
+
+ * the vector space ``X`` is a Hilbert space.
+ * the field of scalars ``F`` are the real numbers.
+
+It is possible to use the class without fulfilling these assumtptions, however in this case some of the mathematical results might not hold.
+
+
+Implementation of functionals
+=============================
+
+To define your own functional, you start by writing::
+
+    class MyFunctional(odl.solvers.Functional):
+        ...
+
+Since `Functional` is a subclass of `Operator`, it has the *abstract method* ``domain`` which need to be explicitly overridden by any subclass.
+
+``domain``: `Set`
+    The domain of this functional, i.e., the set of elements to which this functional can be applied.
+
+Moreover, there are several optional parameters associated with a functional. These are
+
+``linear`` : `bool`, optional
+    If `True`, the functional is considered as linear. In this case, ``domain`` and ``range`` have to be instances of `LinearSpace`, or `Field`.
+    Default: ``False``
+smooth : `bool`, optional
+    If `True`, assume that the functional is continuously differentiable.
+    Default: ``False``
+concave : `bool`, optional
+    If `True`, assume that the functional is concave.
+    Default: ``False``
+convex : `bool`, optional
+    If `True`, assume that the functional is convex.
+    Default: ``False``
+grad_lipschitz : 'float', optional
+    The Lipschitz constant of the gradient.
+    Default: infinity
+
+A functional also has a number of optional methods and properties associated with it. The default value of these are all to raise a `NotImplemetedError`. The properties are:
+
+ * ``gradient``. This returns the gradient operator of the functional, i.e., the operator that corresponds to the mapping
+
+.. math::
+
+    x \\to \\nabla f(x)
+
+where :math:`\\nabla f(x)` is the element used to evaluate derivatives in a direction :math:`d` by :math:`\\langle \\nabla f(x), d \\rangle`.
+ * ``conjugate_functional``. This is the convex conjugate functional, also known as the Legendre transform or Fenchel conjugate. It is defined as 
+
+.. math::
+
+    f^*(x^*) = \sup_{x \in X} \{ \langle x^*,x \rangle - f(x)  \}.
+
+For general linear spaces :math:`X`, :math:`x^* \in X^*` the (continuous/normed) dual space of :math:`X`, i.e., the space of all continuous linear functionals defined on :math:`X`. Then :math:`\langle x^*,x \rangle` is the "dual pairing", i.e., the evaluation of the linear functional :math:`x^*` in the point :math:`x`. However, Hilbert spaces are self-dual, meaning :math:`X^* = X`, and :math:`\langle x^*,x \rangle` is the inner product.
+
+The optional method is
+
+ * ``proximal(sigma)``. This returns the proximal operator of the functional, where ``sigma`` is a nonnegative step-size like parameter. The proximal operator is defined as
+
+.. math::
+
+    \text{prox}_{\sigma f}(x) = \text{arg min}_{y \in X} \{ f(y) + \frac{1}{2\sigma} \|y - x\|_2^2 \}
+
+The `Functional` class also contains two default implementations of two help functions:
+
+* ``derivative(point)``. Given an implementation of the gradient, this method return the (directional) derivative operator in ``point``. This is the linear operator 
+
+.. math::
+    x \to \langle x, \nabla f(point) \rangle,
+
+where :math:`\nabla f(point)` is the gradient of the functional in the point :math:`point`.
+* ``translate(shift)``. Give a functional :math:`f(.)`, this method creates the functional :math:`f(. - shift)`
+
+
+Functional arithmetics
+======================
+It is common in applications to perform arithmetic operations with functionals, for example adding two functionals :math:`f` and :math:`g`
+
+.. math::
+   [f+g](x) = f(x) + g(x),
+
+or multiplication of a functional by a scalar
+
+.. math::
+   [\alpha f](x) = \alpha f (x).
+
+Another example is translating a functional with a vecotr :math:`y`
+
+.. math::
+   f(x - y),
+
+or given an `Operator` :math:`A` whose range is the same as the domain as the functional we also have composition
+
+.. math::
+    [f * A](x) = f(A(x)). 
+
+In some of these cases, properties and methods such as ``gradient``, ``convex_conjugate`` and ``proximal`` can be calculated automatically given a default implementation of the corresponding property in :math:`f`.
+
+All available functional arithmetic, including which properties and methods that automatically can be calculated, is shown below. ``f``, ``g`` represent `Functional`'s, and ``A`` an `Operator` whose range is the same as the domain as the functional. `` a`` is a scalar in the field of the domain of ``f`` and ``g``, and ``y`` is a vector in the domain of ``f`` and ``g``.
+
++------------------+-----------------+-------------------------------------------------+
+| Code             | Meaning         | Class                                           |
++==================+=================+=================================================+
+| ``(f + g)(x)``   | ``f(x) + g(x)`` | `FunctionalSum`                                 |
+|                  |                 | - Retains `gradient`.                           |
++------------------+-----------------+-------------------------------------------------+
+| ``(f + a)(x)     | ``f(x) + a``    | `FunctionalScalarSum`                           | 
+|                  |                 | - Retains all properties.                       |
++------------------+-----------------+-------------------------------------------------+
+| ``(f * A)(x)``   | ``f(A(x))``     | `FunctionalComp`                                |
+|                  |                 | - Retains gradient                              |
++------------------+-----------------+-------------------------------------------------+
+| ``(a * f)(x)``   | ``a * f(x)``    | `FunctionalLeftScalarMult`                      |
+|                  |                 | - Retains all properties, if ``a`` is positive. |
++------------------+-----------------+-------------------------------------------------+
+| ``(f * a)(x)``   | ``f(a * x)``    | `FunctionalRightScalarMult`                     |
+|                  |                 | - Retains all properties                        |
++------------------+-----------------+-------------------------------------------------+
+| ``(v * f)(x)``   | ``v * f(x)``    | `FunctionalLeftVectorMult`                      |
+|                  |                 | - Note that this is not a functional anymore.   |
++------------------+-----------------+-------------------------------------------------+
+| ``(f * v)(x)``   | ``f(v * x)``    | `FunctionalRightVectorMult`                     |
+|                  |                 | - Retains gradient and convex conjugate.        |
++------------------+-----------------+-------------------------------------------------+
+| ``f.translate(y) | ``f(. - y)``    | `TranslatedFunctional`                          |
+|                  |                 | - Retains all properties.                       |
++------------------+-----------------+-------------------------------------------------+
+
+
+
+
+
+
+
diff --git a/doc/source/guide/in_depth/index.rst b/doc/source/guide/in_depth/index.rst
@@ -12,3 +12,4 @@ This is a more in depth guide to the different parts of ODL.
     vectorization_guide
     numpy_guide
     chambolle_pock_guide
+    functional_guide
diff --git a/odl/solvers/functional/functional.py b/odl/solvers/functional/functional.py
@@ -35,10 +35,12 @@
 
 
 # TODO: Add missing functionals here
-__all__ = ('Functional', 'ConvexConjugateArgScaling',
-           'FunctionalLeftScalarMult',
-           'ConvexConjugateFuncScaling', 'ConvexConjugateLinearPerturb',
-           'ConvexConjugateTranslation', 'TranslatedFunctional')
+__all__ = ('Functional', 'FunctionalLeftScalarMult',
+           'FunctionalRightScalarMult', 'FunctionalComp',
+           'FunctionalRightVectorMult', 'FunctionalSum', 'FunctionalScalarSum',
+           'TranslatedFunctional', 'ConvexConjugateTranslation',
+           'ConvexConjugateFuncScaling', 'ConvexConjugateArgScaling',
+           'ConvexConjugateLinearPerturb')
 
 
 class Functional(Operator):
@@ -278,6 +280,8 @@ def __mul__(self, other):
                 return FunctionalLeftScalarMult(self, other)
             else:
                 return FunctionalRightScalarMult(self, other)
+        elif isinstance(other, LinearSpaceVector) and other in self.domain:
+            return FunctionalRightVectorMult(self, other)
         else:
             return super().__mul__(other)
 
@@ -728,6 +732,104 @@ def _call(self, x):
         return CompositGradient()
 
 
+class FunctionalRightVectorMult(Functional, OperatorRightVectorMult):
+
+    """Expression type for the functional right vector multiplication.
+
+    Given a functional ``func`` and a vector ``y`` in the domain of ``func``,
+    this corresponds to the functional
+
+        ``(func * y)(x) = func(y*x)``.
+    """
+
+    def __init__(self, func, vector):
+        """Initialize a new instance.
+
+        Parameters
+        ----------
+        func : `Functional`
+            The domain of ``func`` must be a ``vector.space``.
+        vector : `LinearSpaceVector` in ``func.domain``
+            The vector to multiply by.
+        """
+        if not isinstance(func, Functional):
+            raise TypeError('functional {!r} is not a Functional instance.'
+                            ''.format(func))
+
+        OperatorRightVectorMult.__init__(self, operator=func, vector=vector)
+
+        # TODO: can some of the parameters convex, etc. be decided?
+        Functional.__init__(self, domain=func.domain)
+
+    @property
+    def gradient(self):
+        """Gradient operator of the functional.
+
+        Notes
+        -----
+        The operator that corresponds to the mapping
+
+        .. math::
+
+            x \\to \\nabla f(x)
+
+        where :math:`\\nabla f(x)` is the element used to evaluate
+        derivatives in a direction :math:`d` by
+        :math:`\\langle \\nabla f(x), d \\rangle`.
+        """
+        return self.vector * self.operator.gradient * self.vector
+
+    # TODO: can this be computed?
+    def proximal(self, sigma=1.0):
+        """Return the proximal operator of the functional.
+
+        Parameters
+        ----------
+        sigma : positive float, optional
+            Regularization parameter of the proximal operator.
+
+        Returns
+        -------
+        out : `Operator`
+            Domain and range equal to domain of functional.
+
+        Notes
+        -----
+        The nonsmooth solvers that make use of proximal operators in order to
+        solve a given optimization problem, see for example
+        `forward_backward_pd`, take a `proximal factory` as input. Note that
+        ``Functional.proximal`` is in fact a `proximal factory`.
+        """
+        raise NotImplementedError('there is no known expression for this')
+
+    @property
+    def conjugate_functional(self):
+        """Convex conjugate functional of the functional.
+
+        Notes
+        -----
+        The convex conjugate functional of a convex functional :math:`f(x)`,
+        defined on a Hilber space, is defined as the functional
+
+        .. math::
+
+            f^*(x^*) = \\sup_{x} \{ \\langle x^*,x \\rangle - f(x)  \}.
+
+        The concept is also known as the Legendre transformation.
+
+        References
+        ----------
+        Wikipedia article on `Convex conjugate
+        <https://en.wikipedia.org/wiki/Convex_conjugate>`_.
+
+        Wikipedia article on `Legendre transformation
+        https://en.wikipedia.org/wiki/Legendre_transformation
+
+        For literature references see, e.g., [Lue1969]_, [Roc1970]_.
+        """
+        return self.operator.conjugate_functional * (1.0 / self.vector)
+
+
 class FunctionalSum(Functional, OperatorSum):
 
     """Expression type for the sum of functionals.
diff --git a/test/solvers/functional/functional_test.py b/test/solvers/functional/functional_test.py
@@ -330,17 +330,30 @@ def test_multiplication_with_vector():
 
     x = example_element(space)
     y = example_element(space)
-    func = odl.solvers.L1Norm(space)
+    func = odl.solvers.L2NormSquare(space)
 
     wrong_space = odl.uniform_discr(1, 2, 10)
     y_other_space = example_element(wrong_space)
 
-    # Multiplication from the right. Make sure it is a OperatorRightVectorMult
+    # Multiplication from the right. Make sure it is a
+    # FunctionalRightVectorMult
     func_times_y = func * y
-    assert isinstance(func_times_y, odl.OperatorRightVectorMult)
+    assert isinstance(func_times_y, odl.solvers.FunctionalRightVectorMult)
 
     expected_result = func(y * x)
-    assert almost_equal((func * y)(x), expected_result, places=PLACES)
+    assert almost_equal(func_times_y(x), expected_result, places=PLACES)
+
+    # Test for the gradient.
+    # Explicit calculations: 2*y*y*x
+    expected_result = 2.0 * y * y * x
+    assert all_almost_equal(func_times_y.gradient(x), expected_result,
+                            places=PLACES)
+
+    # Test for conjugate_functional
+    cc_func_times_y = func_times_y.conjugate_functional
+    # Explicit calculations: 1/4 * ||x/y||_2^2
+    expected_result = 1.0 / 4.0 * (x / y).norm()**2
+    assert almost_equal(cc_func_times_y(x), expected_result, places=PLACES)
 
     # Make sure that right muliplication is not allowed with vector from
     # another space
