Skip to content

Commit

Permalink
DOC: improve doc of proximal interface
Browse files Browse the repository at this point in the history
  • Loading branch information
adler-j committed Aug 16, 2016
1 parent ccef55e commit e8858a3
Show file tree
Hide file tree
Showing 4 changed files with 58 additions and 5 deletions.
1 change: 1 addition & 0 deletions doc/source/guide/in_depth/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -11,4 +11,5 @@ This is a more in depth guide to the different parts of ODL.
linearspace_guide
vectorization_guide
numpy_guide
proximal_lang_guide
chambolle_pock_guide
52 changes: 52 additions & 0 deletions doc/source/guide/in_depth/proximal_lang_guide.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,52 @@
.. _numpy_in_depth:

#######################
Using ODL with ProxImaL
#######################

`Proximal
<http://www.proximal-lang.org/en/latest/>`_ is a Python-embedded modeling language for image optimization problems and can be used with ODL to solve typical inverse problems phrased as optimization problems. The package is especially suited for non-differentiable problems such as total variance denoising.

Here is a minimal example of solving Poisson's equation equation on an interval with a TV type regularizer (:math:`\min_x \ 10||-\Delta x - rhs||_2^2 + ||\nabla x||_1`)::

>>> space = odl.uniform_discr(0, 1, 5)
>>> op = -odl.Laplacian(space)
>>> proximal_lang_op = odl.as_proximal_lang_operator(op)
>>> rhs = space.element(lambda x: (x>0.4) & (x<0.6)) # indicator function on [0.4, 0.6]
>>> x = proximal.Variable(space.shape)
>>> prob = proximal.Problem([10 * proximal.sum_squares(x - rhs.asarray()),
>>> proximal.norm1(proximal.grad(x))])
>>> prob.solve()
>>> x.value
array([ 0.02352054, 0.02647946, 0.9 , 0.02647946, 0.02352054])

Notable differences between ODL and ProxImaL
============================================

It may be tempting to try to convert an arbitrary problem from ODL into ProxImaL, but some differences exist.

Norms
-----
Norms in ODL are scaled according to the underlying function space. Hence a sequence of statements converging discretizations give rise to a converging norm::
>>> for n in range(2, 10000):
... X = odl.uniform_discr(0, 1, n)
... print(X.element(lambda x: x).norm())
0.559016994375
0.576628129734
0.577343052266
0.577350268468
>>> 1 / np.sqrt(3) # exact result
0.577350269189

this is not the case in proximal, where the norm depends on the number of discretization points. Hence a scaling that is correct for a problem in ODL needs not be correct in proximal. This also changes the definition of things like the operator norm.

This also has the added effect of changing the definition of derived features, like the spectral norm of operators.

Spaces
------
ODL can represent some complicated spaces, like :math:`\mathbb{R}^3 \times \mathbb{C}^2` through the `ProductSpace` class::

>>> space = odl.ProductSpace(odl.rn(3), odl.cn(2))

This can then be used in solvers and other structures.
6 changes: 3 additions & 3 deletions examples/solvers/proximal_lang_poisson.py
Original file line number Diff line number Diff line change
Expand Up @@ -47,15 +47,15 @@
rhs += odl.phantom.white_noise(space) * np.std(rhs) * 0.1
rhs.show('rhs')

# Convert laplacian to cvx operator
cvx_laplacian = odl.as_proximal_lang_operator(laplacian)
# Convert laplacian to ProxImaL operator
proximal_lang_laplacian = odl.as_proximal_lang_operator(laplacian)

# Convert to array
rhs_arr = rhs.asarray()

# Set up optimization problem
x = proximal.Variable(space.shape)
funcs = [10 * proximal.sum_squares(cvx_laplacian(x) - rhs_arr),
funcs = [10 * proximal.sum_squares(proximal_lang_laplacian(x) - rhs_arr),
proximal.norm1(proximal.grad(x))]

# Solve the problem
Expand Down
4 changes: 2 additions & 2 deletions examples/solvers/proximal_lang_tomography.py
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@
ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl=impl)

# Convert ray transform to proximal language operator
cvx_ray_trafo = odl.as_proximal_lang_operator(ray_trafo)
proximal_lang_ray_trafo = odl.as_proximal_lang_operator(ray_trafo)

# Create sinogram of forward projected phantom with noise
phantom = odl.phantom.shepp_logan(reco_space, modified=True)
Expand All @@ -79,7 +79,7 @@
# Note that proximal is not aware of the problem scaling in ODL, hence the norm
# does not match the norm in the space.
x = proximal.Variable(reco_space.shape)
funcs = [proximal.sum_squares(cvx_ray_trafo(x) - rhs_arr),
funcs = [proximal.sum_squares(proximal_lang_ray_trafo(x) - rhs_arr),
0.2 * proximal.norm1(proximal.grad(x)),
proximal.nonneg(x),
proximal.nonneg(1 - x)]
Expand Down

0 comments on commit e8858a3

Please sign in to comment.