# coding=utf-8
# Copyright 2014-2019 The ODL contributors
#
# This file is part of ODL.
#
# This Source Code Form is subject to the terms of the Mozilla Public License,
# v. 2.0. If a copy of the MPL was not distributed with this file, You can
# obtain one at https://mozilla.org/MPL/2.0/.
from __future__ import print_function, division, absolute_import
import numpy as np
from odl.operator.operator import (
Operator, OperatorComp, OperatorLeftScalarMult, OperatorRightScalarMult,
OperatorRightVectorMult, OperatorSum, OperatorPointwiseProduct)
from odl.operator.default_ops import (IdentityOperator, ConstantOperator)
from odl.solvers.nonsmooth import (proximal_arg_scaling, proximal_translation,
proximal_quadratic_perturbation,
proximal_const_func, proximal_convex_conj)
from odl.util import signature_string, indent
__all__ = ('Functional', 'FunctionalLeftScalarMult',
'FunctionalRightScalarMult', 'FunctionalComp',
'FunctionalRightVectorMult', 'FunctionalSum', 'FunctionalScalarSum',
'FunctionalTranslation', 'InfimalConvolution',
'FunctionalQuadraticPerturb', 'FunctionalProduct',
'FunctionalQuotient', 'BregmanDistance', 'simple_functional')
[docs]class Functional(Operator):
"""Implementation of a functional class.
A functional is an operator ``f`` that maps from some domain ``X`` to the
field of scalars ``F`` associated with the domain:
``f : X -> F``.
Notes
-----
The implementation of the functional class assumes that the domain
:math:`X` is a Hilbert space and that the field of scalars :math:`F` is a
is the real numbers. It is possible to create functions that do not fulfil
these assumptions, however some mathematical results might not be valid in
this case. For more information, see `the ODL functional guide
<http://odlgroup.github.io/odl/guide/in_depth/functional_guide.html>`_.
"""
[docs] def __init__(self, space, linear=False, grad_lipschitz=np.nan):
"""Initialize a new instance.
Parameters
----------
space : `LinearSpace`
The domain of this functional, i.e., the set of elements to
which this functional can be applied.
linear : bool, optional
If `True`, the functional is considered as linear.
grad_lipschitz : float, optional
The Lipschitz constant of the gradient. Default: ``nan``
"""
# Cannot use `super(Functional, self)` here since that breaks
# subclasses with multiple inheritance (at least those where both
# parents implement `__init__`, e.g., in `ScalingFunctional`)
Operator.__init__(self, domain=space, range=space.field, linear=linear)
self.__grad_lipschitz = float(grad_lipschitz)
@property
def grad_lipschitz(self):
"""Lipschitz constant for the gradient of the functional."""
return self.__grad_lipschitz
@grad_lipschitz.setter
def grad_lipschitz(self, value):
"""Setter for the Lipschitz constant for the gradient."""
self.__grad_lipschitz = float(value)
@property
def gradient(self):
r"""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`.
"""
raise NotImplementedError(
'no gradient implemented for functional {!r}'
''.format(self))
@property
def proximal(self):
r"""Proximal factory of the functional.
Notes
-----
The proximal operator of a function :math:`f` is an operator defined as
.. math::
prox_{\sigma f}(x) = \sup_{y} \left\{ f(y) -
\frac{1}{2\sigma} \| y-x \|_2^2 \right\}.
Proximal operators are often used in different optimization algorithms,
especially when designed to handle nonsmooth functionals.
A `proximal factory` is a function that, when called with a step
length :math:`\sigma`, returns the corresponding proximal operator.
The nonsmooth solvers that make use of proximal operators to solve a
given optimization problem take a `proximal factory` as input,
i.e., a function returning a proximal operator. See for example
`forward_backward_pd`.
In general, the step length :math:`\sigma` is expected to be a
positive float, but certain functionals might accept more types of
objects as a stepsize:
- If a functional is a `SeparableSum`, then, instead of a positive
float, one may call the `proximal factory` with a list of positive
floats, and the stepsize are applied to each component individually.
- For certain special functionals like `L1Norm` and `L2NormSquared`,
which are not implemented as a `SeparableSum`, the proximal factory
will accept an argument which is `element-like` regarding the domain
of the functional. Its components must be strictly positive floats.
A stepsize like :math:`(\sigma_1, \ldots, \sigma_n)` coincides
with a matrix-valued distance according to Section XV.4 of [HL1993]
and the rule
.. math::
M = \mathrm{diag}(\sigma_1^{-1}, \ldots, \sigma_n^{-1})
or the Bregman-proximal according to [E1993] and the rule
.. math::
h(x) = \langle x, M x \rangle.
References
----------
[HL1993] Hiriart-Urruty J-B, and Lemaréchal C. *Convex analysis and
minimization algorithms II. Advanced theory and bundle methods.*
Springer, 1993.
[E1993] Eckstein J. *Nonlinear proximal point algorithms using Bregman
functions, with applications to convex programming.* Mathematics of
Operations Research, 18.1 (1993), pp 202--226.
"""
raise NotImplementedError(
'no proximal operator implemented for functional {!r}'
''.format(self))
@property
def convex_conj(self):
r"""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.
For literature references see, e.g., [Lue1969], [Roc1970], the
wikipedia article on `Convex conjugate
<https://en.wikipedia.org/wiki/Convex_conjugate>`_ or the wikipedia
article on the `Legendre transformation
<https://en.wikipedia.org/wiki/Legendre_transformation>`_.
References
----------
[Lue1969] Luenberger, D G. *Optimization by vector space methods*.
Wiley, 1969.
[Roc1970] Rockafellar, R. T. *Convex analysis*. Princeton
University Press, 1970.
"""
return FunctionalDefaultConvexConjugate(self)
[docs] def derivative(self, point):
"""Return the derivative operator in the given point.
This function returns the linear operator given by::
self.derivative(point)(x) == self.gradient(point).inner(x)
Parameters
----------
point : `domain` element
The point in which the gradient is evaluated.
Returns
-------
derivative : `Operator`
"""
return self.gradient(point).T
[docs] def translated(self, shift):
"""Return a translation of the functional.
For a given functional ``f`` and an element ``translation`` in the
domain of ``f``, this operation creates the functional
``f(. - translation)``.
Parameters
----------
translation : `domain` element
Element in the domain of the functional
Returns
-------
out : `FunctionalTranslation`
The functional ``f(. - translation)``
"""
return FunctionalTranslation(self, shift)
[docs] def bregman(self, point, subgrad):
r"""Return the Bregman distance functional.
Parameters
----------
point : element of ``functional.domain``
Point from which to define the Bregman distance.
subgrad : element of ``functional.domain``
A subgradient of ``functional`` in ``point``. If it exists, a
valid option is ``functional.gradient(point)``.
Returns
-------
out : `BregmanDistance`
The Bregman distance functional.
Notes
-----
Given a functional :math:`f`, a point :math:`y`, and a (sub)gradient
:math:`p \in \partial f(y)`, the Bregman distance functional
:math:`D_f^p(\cdot, y)` in a point :math:`x` is given by
.. math::
D_f^p(x, y) = f(x) - f(y) - \langle p, x - y \rangle.
For mathematical details, see
`[Bur2016] <https://arxiv.org/abs/1505.05191>`_. See also the Wikipedia
article: https://en.wikipedia.org/wiki/Bregman_divergence
References
----------
[Bur2016] Burger, M. *Bregman Distances in Inverse Problems and Partial
Differential Equation*. In: Advances in Mathematical Modeling,
Optimization and Optimal Control, 2016. p. 3-33.
"""
return BregmanDistance(self, point, subgrad)
def __mul__(self, other):
"""Return ``self * other``.
If ``other`` is an `Operator`, this corresponds to composition with the
operator:
``(func * op)(x) == func(op(x))``
If ``other`` is a scalar, this corresponds to right multiplication of
scalars with functionals:
``(func * scalar)(x) == func(scalar * x)``
If ``other`` is a vector, this corresponds to right multiplication of
vectors with functionals:
``(func * vector) == func(vector * x)``
Note that left and right multiplications are generally different.
Parameters
----------
other : `Operator`, `domain` element or scalar
`Operator`:
The `Operator.range` of ``other`` must match this functional's
`domain`.
`domain` element:
``other`` must be an element of this functionals's
`Functional.domain`.
scalar:
The `domain` of this functional must be a
`LinearSpace` and ``other`` must be an element of the `field`
of this functional's `domain`. Note that this `field` is also this
functional's `range`.
Returns
-------
mul : `Functional`
Multiplication result.
If ``other`` is an `Operator`, ``mul`` is a
`FunctionalComp`.
If ``other`` is a scalar, ``mul`` is a
`FunctionalRightScalarMult`.
If ``other`` is a vector, ``mul`` is a
`FunctionalRightVectorMult`.
"""
if isinstance(other, Operator):
return FunctionalComp(self, other)
elif other in self.range:
# Left multiplication is more efficient, so we can use this in the
# case of linear functional.
if other == 0:
from odl.solvers.functional.default_functionals import (
ConstantFunctional)
return ConstantFunctional(self.domain,
self(self.domain.zero()))
elif self.is_linear:
return FunctionalLeftScalarMult(self, other)
else:
return FunctionalRightScalarMult(self, other)
elif other in self.domain:
return FunctionalRightVectorMult(self, other)
else:
return super(Functional, self).__mul__(other)
def __rmul__(self, other):
"""Return ``other * self``.
If ``other`` is an `Operator`, since a functional is also an operator
this corresponds to operator composition:
``(op * func)(x) == op(func(x))``
If ``other`` is a scalar, this corresponds to left multiplication of
scalars with functionals:
``(scalar * func)(x) == scalar * func(x)``
If ``other`` is a vector, since a functional is also an operator this
corresponds to left multiplication of vectors with operators:
``(vector * func)(x) == vector * func(x)``
Note that left and right multiplications are generally different.
Parameters
----------
other : `Operator`, `domain` element or scalar
`Operator`:
The `Operator.domain` of ``other`` must match this functional's
`Functional.range`.
`LinearSpaceElement`:
``other`` must be an element of this functionals's
`Functional.range`.
scalar:
The `Operator.domain` of this operator must be a
`LinearSpace` and ``other`` must be an
element of the ``field`` of this operator's
`Operator.domain`.
Returns
-------
rmul : `Functional` or `Operator`
Multiplication result.
If ``other`` is an `Operator`, ``rmul`` is an `OperatorComp`.
If ``other`` is a scalar, ``rmul`` is a
`FunctionalLeftScalarMult`.
If ``other`` is a vector, ``rmul`` is a
`OperatorLeftVectorMult`.
"""
if other in self.range:
if other == 0:
from odl.solvers.functional.default_functionals import (
ZeroFunctional)
return ZeroFunctional(self.domain)
else:
return FunctionalLeftScalarMult(self, other)
else:
return super(Functional, self).__rmul__(other)
def __add__(self, other):
"""Return ``self + other``.
If ``other`` is a `Functional`, this corresponds to
``(func1 + func2)(x) == func1(x) + func2(x)``
If ``other`` is a scalar, this corresponds to adding a scalar to the
value of the functional:
``(func + scalar)(x) == func(x) + scalar``
Parameters
----------
other : `Functional` or scalar
`Functional`:
The `Functional.domain` and `Functional.range` of ``other``
must match this functional's `Functional.domain` and
`Functional.range`.
scalar:
The scalar needs to be in this functional's `Functional.range`.
Returns
-------
sum : `Functional`
Addition result.
If ``other`` is in ``Functional.range``, ``sum`` is a
`FunctionalScalarSum`.
If ``other`` is a `Functional`, ``sum`` is a `FunctionalSum`.
"""
if other in self.domain.field:
return FunctionalScalarSum(self, other)
elif isinstance(other, Functional):
return FunctionalSum(self, other)
else:
return super(Functional, self).__add__(other)
# Since addition is commutative, right and left addition is the same
__radd__ = __add__
def __sub__(self, other):
"""Return ``self - other``."""
return self + (-1) * other
[docs]class FunctionalLeftScalarMult(Functional, OperatorLeftScalarMult):
"""Scalar multiplication of functional from the left.
Given a functional ``f`` and a scalar ``scalar``, this represents the
functional
``(scalar * f)(x) == scalar * f(x)``.
``Functional.__rmul__`` takes care of the case scalar = 0.
"""
[docs] def __init__(self, func, scalar):
"""Initialize a new instance.
Parameters
----------
func : `Functional`
Functional to be scaled.
scalar : float, nonzero
Number with which to scale the functional.
"""
if not isinstance(func, Functional):
raise TypeError('`func` {!r} is not a `Functional` instance'
''.format(func))
Functional.__init__(
self, space=func.domain, linear=func.is_linear,
grad_lipschitz=np.abs(scalar) * func.grad_lipschitz)
OperatorLeftScalarMult.__init__(self, operator=func, scalar=scalar)
@property
def functional(self):
"""The original functional."""
return self.operator
@property
def gradient(self):
"""Gradient operator of the functional."""
return self.scalar * self.functional.gradient
@property
def convex_conj(self):
"""Convex conjugate functional of the scaled functional.
``Functional.__rmul__`` takes care of the case scalar = 0.
"""
if self.scalar <= 0:
raise ValueError('scaling with nonpositive values have no convex '
'conjugate. Current value: {}.'
''.format(self.scalar))
return self.scalar * self.functional.convex_conj * (1.0 / self.scalar)
@property
def proximal(self):
"""Proximal factory of the scaled functional.
``Functional.__rmul__`` takes care of the case scalar = 0
See Also
--------
odl.solvers.nonsmooth.proximal_operators.proximal_const_func
"""
if self.scalar < 0:
raise ValueError('proximal operator of functional scaled with a '
'negative value {} is not well-defined'
''.format(self.scalar))
elif self.scalar == 0:
# Should not get here. `Functional.__rmul__` takes care of the case
# scalar = 0
return proximal_const_func(self.domain)
else:
def proximal_left_scalar_mult(sigma=1.0):
"""Proximal operator for left scalar multiplication.
Parameters
----------
sigma : positive float, optional
Step size parameter. Default: 1.0
"""
return self.functional.proximal(sigma * self.scalar)
return proximal_left_scalar_mult
[docs]class FunctionalRightScalarMult(Functional, OperatorRightScalarMult):
"""Scalar multiplication of the argument of functional.
Given a functional ``f`` and a scalar ``scalar``, this represents the
functional
``(f * scalar)(x) == f(scalar * x)``.
`Functional.__mul__` takes care of the case scalar = 0.
"""
[docs] def __init__(self, func, scalar):
"""Initialize a new instance.
Parameters
----------
func : `Functional`
The functional which will have its argument scaled.
scalar : float, nonzero
The scaling parameter with which the argument is scaled.
"""
if not isinstance(func, Functional):
raise TypeError('`func` {!r} is not a `Functional` instance'
''.format(func))
scalar = func.domain.field.element(scalar)
Functional.__init__(
self, space=func.domain, linear=func.is_linear,
grad_lipschitz=np.abs(scalar) * func.grad_lipschitz)
OperatorRightScalarMult.__init__(self, operator=func, scalar=scalar)
@property
def functional(self):
"""The original functional."""
return self.operator
@property
def gradient(self):
"""Gradient operator of the functional."""
return self.scalar * self.functional.gradient * self.scalar
@property
def convex_conj(self):
"""Convex conjugate functional of functional with scaled argument.
`Functional.__mul__` takes care of the case scalar = 0.
"""
return self.functional.convex_conj * (1 / self.scalar)
@property
def proximal(self):
"""Proximal factory of the functional.
See Also
--------
odl.solvers.nonsmooth.proximal_operators.proximal_arg_scaling
"""
return proximal_arg_scaling(self.functional.proximal, self.scalar)
[docs]class FunctionalComp(Functional, OperatorComp):
"""Composition of a functional with an operator.
Given a functional ``func`` and an operator ``op``, such that the range of
the operator is equal to the domain of the functional, this corresponds to
the functional
``(func * op)(x) == func(op(x))``.
"""
[docs] def __init__(self, func, op):
"""Initialize a new instance.
Parameters
----------
func : `Functional`
The left ("outer") operator
op : `Operator`
The right ("inner") operator. Its range must coincide with the
domain of ``func``.
"""
if not isinstance(func, Functional):
raise TypeError('`fun` {!r} is not a `Functional` instance'
''.format(func))
OperatorComp.__init__(self, left=func, right=op)
Functional.__init__(self, space=op.domain,
linear=(func.is_linear and op.is_linear),
grad_lipschitz=np.nan)
@property
def gradient(self):
"""Gradient of the compositon according to the chain rule."""
func = self.left
op = self.right
class FunctionalCompositionGradient(Operator):
"""Gradient of the compositon according to the chain rule."""
def __init__(self):
"""Initialize a new instance."""
super(FunctionalCompositionGradient, self).__init__(
op.domain, op.domain, linear=False)
def _call(self, x):
"""Apply the gradient operator to the given point."""
return op.derivative(x).adjoint(func.gradient(op(x)))
def derivative(self, x):
"""The derivative in point ``x``.
This is only defined
"""
if not op.is_linear:
raise NotImplementedError('derivative only implemented '
'for linear opertors.')
else:
return (op.adjoint * func.gradient * op).derivative(x)
return FunctionalCompositionGradient()
[docs]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)``.
"""
[docs] def __init__(self, func, vector):
"""Initialize a new instance.
Parameters
----------
func : `Functional`
The domain of ``func`` must be a ``vector.space``.
vector : `domain` element
The vector to multiply by.
"""
if not isinstance(func, Functional):
raise TypeError('`fun` {!r} is not a `Functional` instance'
''.format(func))
OperatorRightVectorMult.__init__(self, operator=func, vector=vector)
Functional.__init__(self, space=func.domain)
@property
def functional(self):
return self.operator
@property
def gradient(self):
"""Gradient operator of the functional."""
return self.vector * self.operator.gradient * self.vector
@property
def convex_conj(self):
"""Convex conjugate functional of the functional.
This is only defined for vectors with no zero-elements.
"""
return self.functional.convex_conj * (1.0 / self.vector)
[docs]class FunctionalSum(Functional, OperatorSum):
"""Expression type for the sum of functionals.
``FunctionalSum(func1, func2) == (x --> func1(x) + func2(x))``.
"""
[docs] def __init__(self, left, right):
"""Initialize a new instance.
Parameters
----------
left, right : `Functional`
The summands of the functional sum. Their `Functional.domain`
and `Functional.range` must coincide.
"""
if not isinstance(left, Functional):
raise TypeError('`left` {!r} is not a `Functional` instance'
''.format(left))
if not isinstance(right, Functional):
raise TypeError('`right` {!r} is not a `Functional` instance'
''.format(right))
Functional.__init__(
self, space=left.domain,
linear=(left.is_linear and right.is_linear),
grad_lipschitz=left.grad_lipschitz + right.grad_lipschitz)
OperatorSum.__init__(self, left, right)
@property
def gradient(self):
"""Gradient operator of functional sum."""
return self.left.gradient + self.right.gradient
[docs]class FunctionalScalarSum(FunctionalSum):
"""Expression type for the sum of a functional and a scalar.
``FunctionalScalarSum(func, scalar) == (x --> func(x) + scalar)``
"""
[docs] def __init__(self, func, scalar):
"""Initialize a new instance.
Parameters
----------
func : `Functional`
Functional to which the scalar is added.
scalar : `element` in the `field` of the ``domain``
The scalar to be added to the functional. The `field` of the
``domain`` is the range of the functional.
"""
from odl.solvers.functional.default_functionals import (
ConstantFunctional)
if not isinstance(func, Functional):
raise TypeError('`fun` {!r} is not a `Functional` instance'
''.format(func))
if scalar not in func.range:
raise TypeError('`scalar` {} is not in the range of '
'`func` {!r}'.format(scalar, func))
super(FunctionalScalarSum, self).__init__(
left=func,
right=ConstantFunctional(space=func.domain, constant=scalar))
@property
def scalar(self):
"""The scalar that is added to the functional"""
return self.right.constant
@property
def proximal(self):
"""Proximal factory of the FunctionalScalarSum."""
return self.left.proximal
@property
def convex_conj(self):
"""Convex conjugate functional of FunctionalScalarSum."""
return self.left.convex_conj - self.scalar
[docs]class FunctionalTranslation(Functional):
"""Implementation of the translated functional.
Given a functional ``f`` and an element ``translation`` in the domain of
``f``, this corresponds to the functional ``f(. - translation)``.
"""
[docs] def __init__(self, func, translation):
"""Initialize a new instance.
Given a functional ``f(.)`` and a vector ``translation`` in the domain
of ``f``, this corresponds to the functional ``f(. - translation)``.
Parameters
----------
func : `Functional`
Functional which is to be translated.
translation : `domain` element
The translation.
"""
if not isinstance(func, Functional):
raise TypeError('`func` {!r} not a `Functional` instance'
''.format(func))
translation = func.domain.element(translation)
super(FunctionalTranslation, self).__init__(
space=func.domain, linear=False,
grad_lipschitz=func.grad_lipschitz)
# TODO: Add case if we have translation -> scaling -> translation?
if isinstance(func, FunctionalTranslation):
self.__functional = func.functional
self.__translation = func.translation + translation
else:
self.__functional = func
self.__translation = translation
@property
def functional(self):
"""The original functional that has been translated."""
return self.__functional
@property
def translation(self):
"""The translation."""
return self.__translation
[docs] def _call(self, x):
"""Evaluate the functional in a point ``x``."""
return self.functional(x - self.translation)
@property
def gradient(self):
"""Gradient operator of the functional."""
return (self.functional.gradient *
(IdentityOperator(self.domain) - self.translation))
@property
def proximal(self):
"""Proximal factory of the translated functional.
See Also
--------
odl.solvers.nonsmooth.proximal_operators.proximal_translation
"""
return proximal_translation(self.functional.proximal,
self.translation)
@property
def convex_conj(self):
r"""Convex conjugate functional of the translated functional.
Notes
-----
Given a functional :math:`f`, the convex conjugate of a translated
version :math:`f(\cdot - y)` is given by a linear pertubation of the
convex conjugate of :math:`f`:
.. math::
(f( . - y))^* (x) = f^*(x) + <y, x>.
For reference on the identity used, see [KP2015].
References
----------
[KP2015] Komodakis, N, and Pesquet, J-C. *Playing with Duality: An
overview of recent primal-dual approaches for solving large-scale
optimization problems*. IEEE Signal Processing Magazine, 32.6 (2015),
pp 31--54.
"""
return FunctionalQuadraticPerturb(
self.functional.convex_conj,
linear_term=self.translation)
def __repr__(self):
"""Return ``repr(self)``."""
return '{!r}.translated({!r})'.format(self.functional,
self.translation)
def __str__(self):
"""Return ``str(self)``."""
return '{}.translated({})'.format(self.functional,
self.translation)
[docs]class InfimalConvolution(Functional):
"""Functional representing ``h(x) = inf_y f(x-y) + g(y)``."""
[docs] def __init__(self, left, right):
"""Initialize a new instance.
Parameters
----------
left : `Functional`
Function corresponding to ``f``.
right : `Functional`
Function corresponding to ``g``.
Examples
--------
>>> space = odl.rn(3)
>>> l1 = odl.solvers.L1Norm(space)
>>> l2 = odl.solvers.L2Norm(space)
>>> f = odl.solvers.InfimalConvolution(l1.convex_conj, l2.convex_conj)
>>> x = f.domain.one()
>>> f.convex_conj(x) - (l1(x) + l2(x))
0.0
"""
if not isinstance(left, Functional):
raise TypeError('`func` {} is not a `Functional` instance'
''.format(left))
if not isinstance(right, Functional):
raise TypeError('`func` {} is not a `Functional` instance'
''.format(right))
super(InfimalConvolution, self).__init__(
space=left.domain, linear=False, grad_lipschitz=np.nan)
self.__left = left
self.__right = right
@property
def left(self):
"""Left functional."""
return self.__left
@property
def right(self):
"""Right functional."""
return self.__right
@property
def convex_conj(self):
"""Convex conjugate functional of the functional.
Notes
-----
The convex conjugate of the infimal convolution
.. math::
h(x) = inf_y f(x-y) + g(y)
is the sum of it:
.. math::
h^*(x) = f^*(x) + g^*(x)
"""
return self.left.convex_conj + self.right.convex_conj
def __repr__(self):
"""Return ``repr(self)``."""
posargs = [self.left, self.right]
inner_str = signature_string(posargs, [], sep=',\n')
return '{}(\n{}\n)'.format(self.__class__.__name__, indent(inner_str))
def __str__(self):
"""Return ``str(self)``."""
return repr(self)
[docs]class FunctionalQuadraticPerturb(Functional):
"""The functional representing ``F(.) + a * <., .> + <., u> + c``."""
[docs] def __init__(self, func, quadratic_coeff=0, linear_term=None,
constant=0):
"""Initialize a new instance.
Parameters
----------
func : `Functional`
Function corresponding to ``f``.
quadratic_coeff : ``domain.field`` element, optional
Coefficient of the quadratic term. Default: 0.
linear_term : `domain` element, optional
Element in domain of ``func``, corresponding to the translation.
Default: Zero element.
constant : ``domain.field`` element, optional
The constant coefficient. Default: 0.
"""
if not isinstance(func, Functional):
raise TypeError('`func` {} is not a `Functional` instance'
''.format(func))
self.__functional = func
quadratic_coeff = func.domain.field.element(quadratic_coeff)
if quadratic_coeff.imag != 0:
raise ValueError(
"Complex-valued quadratic coefficient is not supported.")
self.__quadratic_coeff = quadratic_coeff.real
if linear_term is not None:
self.__linear_term = func.domain.element(linear_term)
else:
self.__linear_term = func.domain.zero()
if linear_term is None:
grad_lipschitz = func.grad_lipschitz
else:
grad_lipschitz = (func.grad_lipschitz + self.linear_term.norm())
constant = func.domain.field.element(constant)
if constant.imag != 0:
raise ValueError(
"Complex-valued `constant` coefficient is not supported.")
self.__constant = constant.real
super(FunctionalQuadraticPerturb, self).__init__(
space=func.domain,
linear=func.is_linear and (quadratic_coeff == 0),
grad_lipschitz=grad_lipschitz)
@property
def functional(self):
"""Original functional."""
return self.__functional
@property
def quadratic_coeff(self):
"""Cofficient of the quadratic term."""
return self.__quadratic_coeff
@property
def linear_term(self):
"""Linear term."""
return self.__linear_term
@property
def constant(self):
"""The constant coefficient."""
return self.__constant
[docs] def _call(self, x):
"""Apply the functional to the given point."""
return (self.functional(x) +
self.quadratic_coeff * x.inner(x) +
x.inner(self.linear_term) + self.constant)
@property
def gradient(self):
"""Gradient operator of the functional."""
return (self.functional.gradient +
(2 * self.quadratic_coeff) * IdentityOperator(self.domain) +
ConstantOperator(self.linear_term))
@property
def proximal(self):
"""Proximal factory of the quadratically perturbed functional."""
if self.quadratic_coeff < 0:
raise TypeError('`quadratic_coeff` {} must be non-negative'
''.format(self.quadratic_coeff))
return proximal_quadratic_perturbation(
self.functional.proximal,
a=self.quadratic_coeff, u=self.linear_term)
@property
def convex_conj(self):
r"""Convex conjugate functional of the functional.
Notes
-----
Given a functional :math:`f`, the convex conjugate of a linearly
perturbed version :math:`f(x) + <y, x>` is given by a translation of
the convex conjugate of :math:`f`:
.. math::
(f + \langle y, \cdot \rangle)^* (x^*) = f^*(x^* - y).
For reference on the identity used, see `[KP2015]`_. Moreover, the
convex conjugate of :math:`f + c` is by definition
.. math::
(f + c)^* (x^*) = f^*(x^*) - c.
References
----------
[KP2015] Komodakis, N, and Pesquet, J-C. *Playing with Duality: An
overview of recent primal-dual approaches for solving large-scale
optimization problems*. IEEE Signal Processing Magazine, 32.6 (2015),
pp 31--54.
.. _[KP2015]: https://arxiv.org/abs/1406.5429
"""
if self.quadratic_coeff == 0:
cconj = self.functional.convex_conj.translated(self.linear_term)
if self.constant != 0:
cconj = cconj - self.constant
return cconj
else:
return super(FunctionalQuadraticPerturb, self).convex_conj
def __repr__(self):
"""Return ``repr(self)``."""
return '{}({!r}, {!r}, {!r}, {!r})'.format(self.__class__.__name__,
self.functional,
self.quadratic_coeff,
self.linear_term,
self.constant)
def __str__(self):
"""Return ``str(self)``."""
return '{}({}, {}, {}, {})'.format(self.__class__.__name__,
self.functional,
self.quadratic_coeff,
self.linear_term,
self.constant)
[docs]class FunctionalProduct(Functional, OperatorPointwiseProduct):
"""Product ``p(x) = f(x) * g(x)`` of two functionals ``f`` and ``g``."""
[docs] def __init__(self, left, right):
"""Initialize a new instance.
Parameters
----------
left, right : `Functional`
Functionals in the product. Need to have matching domains.
Examples
--------
Construct the functional || . ||_2^2 * 3
>>> space = odl.rn(2)
>>> func1 = odl.solvers.L2NormSquared(space)
>>> func2 = odl.solvers.ConstantFunctional(space, 3)
>>> prod = odl.solvers.FunctionalProduct(func1, func2)
>>> prod([2, 3]) # expect (2**2 + 3**2) * 3 = 39
39.0
"""
if not isinstance(left, Functional):
raise TypeError('`left` {} is not a `Functional` instance'
''.format(left))
if not isinstance(right, Functional):
raise TypeError('`right` {} is not a `Functional` instance'
''.format(right))
OperatorPointwiseProduct.__init__(self, left, right)
Functional.__init__(self, left.domain, linear=False,
grad_lipschitz=np.nan)
@property
def gradient(self):
r"""Gradient operator of the functional.
Notes
-----
The derivative is computed using Leibniz's rule:
.. math::
[\nabla (f g)](p) = g(p) [\nabla f](p) + f(p) [\nabla g](p)
"""
func = self
class FunctionalProductGradient(Operator):
"""Functional representing the gradient of ``f(.) * g(.)``."""
def _call(self, x):
return (func.right(x) * func.left.gradient(x) +
func.left(x) * func.right.gradient(x))
return FunctionalProductGradient(self.domain, self.domain,
linear=False)
[docs]class FunctionalQuotient(Functional):
"""Quotient ``p(x) = f(x) / g(x)`` of two functionals ``f`` and ``g``."""
[docs] def __init__(self, dividend, divisor):
"""Initialize a new instance.
Parameters
----------
dividend, divisor : `Functional`
Functionals in the quotient. Need to have matching domains.
Examples
--------
Construct the functional || . ||_2 / 5
>>> space = odl.rn(2)
>>> func1 = odl.solvers.L2Norm(space)
>>> func2 = odl.solvers.ConstantFunctional(space, 5)
>>> prod = odl.solvers.FunctionalQuotient(func1, func2)
>>> prod([3, 4]) # expect sqrt(3**2 + 4**2) / 5 = 1
1.0
"""
if not isinstance(dividend, Functional):
raise TypeError('`dividend` {} is not a `Functional` instance'
''.format(dividend))
if not isinstance(divisor, Functional):
raise TypeError('`divisor` {} is not a `Functional` instance'
''.format(divisor))
if dividend.domain != divisor.domain:
raise ValueError('domains of the operators do not match')
self.__dividend = dividend
self.__divisor = divisor
super(FunctionalQuotient, self).__init__(
dividend.domain, linear=False, grad_lipschitz=np.nan)
@property
def dividend(self):
"""The dividend of the quotient."""
return self.__dividend
@property
def divisor(self):
"""The divisor of the quotient."""
return self.__divisor
[docs] def _call(self, x):
"""Apply the functional to the given point."""
return self.dividend(x) / self.divisor(x)
@property
def gradient(self):
r"""Gradient operator of the functional.
Notes
-----
The derivative is computed using the quotient rule:
.. math::
[\nabla (f / g)](p) = (g(p) [\nabla f](p) -
f(p) [\nabla g](p)) / g(p)^2
"""
func = self
class FunctionalQuotientGradient(Operator):
"""Functional representing the gradient of ``f(.) / g(.)``."""
def _call(self, x):
"""Apply the functional to the given point."""
dividendx = func.dividend(x)
divisorx = func.divisor(x)
return ((1 / divisorx) * func.dividend.gradient(x) +
(- dividendx / divisorx**2) * func.divisor.gradient(x))
return FunctionalQuotientGradient(self.domain, self.domain,
linear=False)
def __repr__(self):
"""Return ``repr(self)``."""
return '{}({!r}, {!r})'.format(self.__class__.__name__,
self.dividend, self.divisor)
def __str__(self):
"""Return ``str(self)``."""
return '{}({}, {})'.format(self.__class__.__name__,
self.dividend, self.divisor)
[docs]class FunctionalDefaultConvexConjugate(Functional):
r"""The `Functional` representing ``F^*``, the convex conjugate of ``F``.
This class does not provide a way to evaluate the functional, it is rather
intended to be used for its `proximal`.
Notes
-----
The proximal is found by using the Moreau identity
.. math::
\text{prox}_{\sigma F^*}(y) = y -
\sigma \text{prox}_{F / \sigma}(y / \sigma)
which allows the proximal of the convex conjugate to be calculated without
explicit knowledge about the convex conjugate itself.
"""
[docs] def __init__(self, func):
"""Initialize a new instance.
Parameters
----------
func : `Functional`
Functional corresponding to F.
"""
if not isinstance(func, Functional):
raise TypeError('`func` {} is not a `Functional` instance'
''.format(func))
super(FunctionalDefaultConvexConjugate, self).__init__(
space=func.domain, linear=func.is_linear)
self.__convex_conj = func
@property
def convex_conj(self):
"""The original functional."""
return self.__convex_conj
@property
def proximal(self):
"""Proximal factory using the Moreu identity.
Returns
-------
proximal : proximal_convex_conj
Proximal computed using the Moreu identity
"""
return proximal_convex_conj(self.convex_conj.proximal)
def __repr__(self):
"""Return ``repr(self)``."""
return '{!r}.convex_conj'.format(self.convex_conj)
def __str__(self):
"""Return ``str(self)``."""
return '{}.convex_conj'.format(self.convex_conj)
[docs]class BregmanDistance(Functional):
r"""The Bregman distance functional.
The Bregman distance, also refered to as the Bregman divergence, is similar
to a metric but satisfies neither the triangle inequality nor symmetry.
Nevertheless, the Bregman distance is used in variational regularization of
inverse problems, see, e.g., `[Bur2016]`_.
Notes
-----
Given a functional :math:`f`, a point :math:`y`, and a (sub)gradient
:math:`p \in \partial f(y)`, the Bregman distance functional
:math:`D_f^p(\cdot, y)` in a point :math:`x` is given by
.. math::
D_f^p(x, y) = f(x) - f(y) - \langle p, x - y \rangle.
For mathematical details, see
`[Bur2016] <https://arxiv.org/abs/1505.05191>`_. See also the Wikipedia
article: https://en.wikipedia.org/wiki/Bregman_divergence
References
----------
[Bur2016] Burger, M. *Bregman Distances in Inverse Problems and Partial
Differential Equation*. In: Advances in Mathematical Modeling, Optimization
and Optimal Control, 2016. p. 3-33.
"""
[docs] def __init__(self, functional, point, subgrad):
"""Initialize a new instance.
Parameters
----------
functional : `Functional`
Functional on which to base the Bregman distance.
point : element of ``functional.domain``
Point from which to define the Bregman distance.
subgrad : element of ``functional.domain``
A subgradient of ``functional`` in ``point``. If it exists, a
valid option is ``functional.gradient(point)``.
Examples
--------
Example of initializing the Bregman distance functional:
>>> space = odl.uniform_discr(0, 1, 10)
>>> l2_squared = odl.solvers.L2NormSquared(space)
>>> point = space.one()
>>> subgrad = l2_squared.gradient(point)
>>> bregman_dist = odl.solvers.BregmanDistance(
... l2_squared, point, subgrad)
This is gives squared L2 distance to the given point, ||x - 1||^2:
>>> expected_functional = l2_squared.translated(point)
>>> bregman_dist(space.zero()) == expected_functional(space.zero())
True
"""
if not isinstance(functional, Functional):
raise TypeError('`functional` {} not an instance of ``Functional``'
''.format(functional))
self.__functional = functional
if point not in functional.domain:
raise ValueError('`point` {} is not in `functional.domain` {}'
''.format(point, functional.domain))
self.__point = point
if subgrad not in functional.domain:
raise TypeError(
'`subgrad` must be an element in `functional.domain`, got '
'{}'.format(subgrad))
self.__subgrad = subgrad
self.__constant = -functional(point) + subgrad.inner(point)
self.__bregman_dist = FunctionalQuadraticPerturb(
functional, linear_term=-subgrad, constant=self.__constant)
grad_lipschitz = functional.grad_lipschitz + subgrad.norm()
super(BregmanDistance, self).__init__(
space=functional.domain, linear=False,
grad_lipschitz=grad_lipschitz)
@property
def functional(self):
"""The functional used to define the Bregman distance."""
return self.__functional
@property
def point(self):
"""The point used to define the Bregman distance."""
return self.__point
@property
def subgrad(self):
"""The subgradient used to define the Bregman distance."""
return self.__subgrad
[docs] def _call(self, x):
"""Return ``self(x)``."""
return self.__bregman_dist(x)
@property
def convex_conj(self):
"""The convex conjugate"""
return self.__bregman_dist.convex_conj
@property
def proximal(self):
"""Return the ``proximal factory`` of the functional."""
return self.__bregman_dist.proximal
@property
def gradient(self):
"""Gradient operator of the functional."""
try:
op_to_return = self.functional.gradient
except NotImplementedError:
raise NotImplementedError(
'`self.functional.gradient` is not implemented for '
'`self.functional` {}'.format(self.functional))
op_to_return = op_to_return - ConstantOperator(self.subgrad)
return op_to_return
def __repr__(self):
'''Return ``repr(self)``.'''
posargs = [self.functional, self.point, self.subgrad]
optargs = []
inner_str = signature_string(posargs, optargs, sep=',\n')
return '{}(\n{}\n)'.format(self.__class__.__name__, indent(inner_str))
[docs]def simple_functional(space, fcall=None, grad=None, prox=None, grad_lip=np.nan,
convex_conj_fcall=None, convex_conj_grad=None,
convex_conj_prox=None, convex_conj_grad_lip=np.nan,
linear=False):
"""Simplified interface to create a functional with specific properties.
Users may specify as many properties as-is needed by the application.
Parameters
----------
space : `LinearSpace`
Space that the functional should act on.
fcall : callable, optional
Function to evaluate when calling the functional.
grad : callable or `Operator`, optional
Gradient operator of the functional.
prox : `proximal factory`, optional
Proximal factory for the functional.
grad_lip : float, optional
lipschitz constant of the functional.
convex_conj_fcall : callable, optional
Function to evaluate when calling the convex conjugate functional.
convex_conj_grad : callable or `Operator`, optional
Gradient operator of the convex conjugate functional
convex_conj_prox : `proximal factory`, optional
Proximal factory for the convex conjugate functional.
convex_conj_grad_lip : float, optional
lipschitz constant of the convex conjugate functional.
linear : bool, optional
True if the operator is linear.
Examples
--------
Create squared sum functional on rn:
>>> def f(x):
... return sum(xi**2 for xi in x)
>>> def dfdx(x):
... return 2 * x
>>> space = odl.rn(3)
>>> func = simple_functional(space, f, grad=dfdx)
>>> func.domain
rn(3)
>>> func.range
RealNumbers()
>>> func([1, 2, 3])
14.0
>>> func.gradient([1, 2, 3])
rn(3).element([ 2., 4., 6.])
"""
if grad is not None and not isinstance(grad, Operator):
grad_in = grad
class SimpleFunctionalGradient(Operator):
"""Gradient of a `SimpleFunctional`."""
def _call(self, x):
"""Return ``self(x)``."""
return grad_in(x)
grad = SimpleFunctionalGradient(space, space, linear=False)
if (convex_conj_grad is not None and
not isinstance(convex_conj_grad, Operator)):
convex_conj_grad_in = convex_conj_grad
class SimpleFunctionalConvexConjGradient(Operator):
"""Gradient of the convex conj of a `SimpleFunctional`."""
def _call(self, x):
"""Return ``self(x)``."""
return convex_conj_grad_in(x)
convex_conj_grad = SimpleFunctionalConvexConjGradient(
space, space, linear=False)
class SimpleFunctional(Functional):
"""A simplified functional for examples."""
def __init__(self):
"""Initialize an instance."""
super(SimpleFunctional, self).__init__(
space, linear=linear, grad_lipschitz=grad_lip)
def _call(self, x):
"""Return ``self(x)``."""
if fcall is None:
raise NotImplementedError('call not implemented')
else:
return fcall(x)
@property
def proximal(self):
"""Return the proximal of the operator."""
if prox is None:
raise NotImplementedError('proximal not implemented')
else:
return prox
@property
def gradient(self):
"""Return the gradient of the operator."""
if grad is None:
raise NotImplementedError('gradient not implemented')
else:
return grad
@property
def convex_conj(self):
return simple_functional(space, fcall=convex_conj_fcall,
grad=convex_conj_grad,
prox=convex_conj_prox,
grad_lip=convex_conj_grad_lip,
convex_conj_fcall=fcall,
convex_conj_grad=grad,
convex_conj_prox=prox,
convex_conj_grad_lip=grad_lip,
linear=linear)
return SimpleFunctional()
if __name__ == '__main__':
from odl.util.testutils import run_doctests
run_doctests()