Adaptive Primal-Dual#

This tutorial compares the traditional Chambolle-Pock Primal-dual algorithm with the Adaptive Primal-Dual Hybrid Gradient of Goldstein and co-authors.

By adaptively changing the step size in the primal and the dual directions, this algorithm shows faster convergence, which is of great importance for some of the problems that the Primal-Dual algorithm can solve - especially those with an expensive proximal operator.

For this example, we consider a simple denoising problem.

import numpy as np
import matplotlib.pyplot as plt
import pylops
from import camera

import pyproximal


def callback(x, f, g, K, cost, xtrue, err):
    cost.append(f(x) + g(K.matvec(x)))
    err.append(np.linalg.norm(x - xtrue))

Let’s start by loading a sample image and adding some noise

# Load image
img = camera()
ny, nx = img.shape

# Add noise
sigman = 20
n = np.random.normal(0, sigman, img.shape)
noise_img = img + n

We can now define a pylops.Gradient operator as well as the different proximal operators to be passed to our solvers

# Gradient operator
sampling = 1.
Gop = pylops.Gradient(dims=(ny, nx), sampling=sampling, edge=False,
                      kind='forward', dtype='float64')
L = 8. / sampling ** 2 # maxeig(Gop^H Gop)

# L2 data term
lamda = .04
l2 = pyproximal.L2(b=noise_img.ravel(), sigma=lamda)

# L1 regularization (isotropic TV)
l1iso = pyproximal.L21(ndim=2)

To start, we solve our denoising problem with the original Primal-Dual algorithm

# Primal-dual
tau = 0.95 / np.sqrt(L)
mu = 0.95 / np.sqrt(L)

cost_fixed = []
err_fixed = []
iml12_fixed = \
    pyproximal.optimization.primaldual.PrimalDual(l2, l1iso, Gop,
                                                  tau=tau, mu=mu, theta=1.,
                                                  gfirst=False, niter=300, show=True,
                                                  callback=lambda x: callback(x, l2, l1iso,
                                                                              Gop, cost_fixed,
iml12_fixed = iml12_fixed.reshape(img.shape)
Primal-dual: min_x f(Ax) + x^T z + g(x)
Proximal operator (f): <class 'pyproximal.proximal.L2.L2'>
Proximal operator (g): <class 'pyproximal.proximal.L21.L21'>
Linear operator (A): <class 'pylops.basicoperators.gradient.Gradient'>
Additional vector (z): None
tau = 0.33587572106361          mu = 0.33587572106361
theta = 1.00            niter = 300

   Itn       x[0]          f           g          z^x       J = f + g + z^x
     1   2.26495e+00   1.148e+08   1.329e+05   0.000e+00       1.150e+08
     2   4.70719e+00   1.119e+08   1.386e+05   0.000e+00       1.120e+08
     3   7.26347e+00   1.090e+08   1.220e+05   0.000e+00       1.091e+08
     4   9.83797e+00   1.062e+08   1.118e+05   0.000e+00       1.063e+08
     5   1.23754e+01   1.035e+08   1.110e+05   0.000e+00       1.036e+08
     6   1.48911e+01   1.009e+08   1.145e+05   0.000e+00       1.010e+08
     7   1.73926e+01   9.828e+07   1.189e+05   0.000e+00       9.840e+07
     8   1.98782e+01   9.576e+07   1.243e+05   0.000e+00       9.589e+07
     9   2.23416e+01   9.331e+07   1.306e+05   0.000e+00       9.344e+07
    10   2.47768e+01   9.092e+07   1.376e+05   0.000e+00       9.106e+07
    31   6.82254e+01   5.304e+07   2.882e+05   0.000e+00       5.332e+07
    61   1.11156e+02   2.521e+07   4.538e+05   0.000e+00       2.567e+07
    91   1.40820e+02   1.268e+07   5.666e+05   0.000e+00       1.325e+07
   121   1.60553e+02   7.029e+06   6.424e+05   0.000e+00       7.672e+06
   151   1.73850e+02   4.475e+06   6.933e+05   0.000e+00       5.168e+06
   181   1.82651e+02   3.317e+06   7.274e+05   0.000e+00       4.044e+06
   211   1.88592e+02   2.789e+06   7.503e+05   0.000e+00       3.540e+06
   241   1.92588e+02   2.547e+06   7.656e+05   0.000e+00       3.313e+06
   271   1.95246e+02   2.435e+06   7.759e+05   0.000e+00       3.211e+06
   292   1.96549e+02   2.395e+06   7.810e+05   0.000e+00       3.176e+06
   293   1.96603e+02   2.393e+06   7.812e+05   0.000e+00       3.174e+06
   294   1.96655e+02   2.392e+06   7.814e+05   0.000e+00       3.173e+06
   295   1.96707e+02   2.390e+06   7.816e+05   0.000e+00       3.172e+06
   296   1.96759e+02   2.389e+06   7.818e+05   0.000e+00       3.171e+06
   297   1.96809e+02   2.388e+06   7.820e+05   0.000e+00       3.170e+06
   298   1.96859e+02   2.386e+06   7.822e+05   0.000e+00       3.169e+06
   299   1.96909e+02   2.385e+06   7.824e+05   0.000e+00       3.167e+06
   300   1.96957e+02   2.384e+06   7.826e+05   0.000e+00       3.166e+06

Total time (s) = 13.50

We do the same with the adaptive algorithm

cost_ada = []
err_ada = []
iml12_ada, steps = \
    pyproximal.optimization.primaldual.AdaptivePrimalDual(l2, l1iso, Gop,
                                                          tau=tau, mu=mu,
                                                          niter=45, show=True, tol=0.05,
                                                          callback=lambda x: callback(x, l2, l1iso,
                                                                                      Gop, cost_ada,
iml12_ada = iml12_ada.reshape(img.shape)
Adaptive Primal-dual: min_x f(Ax) + x^T z + g(x)
Proximal operator (f): <class 'pyproximal.proximal.L2.L2'>
Proximal operator (g): <class 'pyproximal.proximal.L21.L21'>
Linear operator (A): <class 'pylops.basicoperators.gradient.Gradient'>
Additional vector (z): None
tau0 = 3.358757e-01     mu0 = 3.358757e-01
alpha0 = 5.000000e-01   eta = 9.500000e-01
s = 1.000000e+00        delta = 1.500000e+00
niter = 45              tol = 5.000000e-02

   Itn       x[0]          f           g          z^x       J = f + g + z^x
     2   2.26495e+00   1.148e+08   1.329e+05   0.000e+00       1.150e+08
     3   7.08552e+00   1.090e+08   1.628e+05   0.000e+00       1.092e+08
     4   1.61658e+01   9.892e+07   2.032e+05   0.000e+00       9.913e+07
     5   3.15279e+01   8.319e+07   2.861e+05   0.000e+00       8.348e+07
     6   5.47617e+01   6.215e+07   4.080e+05   0.000e+00       6.255e+07
     7   8.55772e+01   3.917e+07   5.569e+05   0.000e+00       3.972e+07
     8   1.09818e+02   2.502e+07   6.657e+05   0.000e+00       2.569e+07
     9   1.28922e+02   1.631e+07   7.410e+05   0.000e+00       1.706e+07
    10   1.44003e+02   1.095e+07   7.922e+05   0.000e+00       1.174e+07
    13   1.70313e+02   4.735e+06   8.578e+05   0.000e+00       5.593e+06
    17   1.77332e+02   3.762e+06   8.466e+05   0.000e+00       4.609e+06
    21   1.81822e+02   3.311e+06   8.340e+05   0.000e+00       4.145e+06
    25   1.86603e+02   2.935e+06   8.434e+05   0.000e+00       3.779e+06
    29   1.90464e+02   2.700e+06   8.342e+05   0.000e+00       3.534e+06
    33   1.93277e+02   2.560e+06   8.168e+05   0.000e+00       3.377e+06
    37   1.94751e+02   2.474e+06   8.061e+05   0.000e+00       3.281e+06
    38   1.94975e+02   2.459e+06   8.045e+05   0.000e+00       3.263e+06
    39   1.95168e+02   2.445e+06   8.032e+05   0.000e+00       3.248e+06
    40   1.95341e+02   2.432e+06   8.021e+05   0.000e+00       3.234e+06
    41   1.95503e+02   2.421e+06   8.013e+05   0.000e+00       3.222e+06
    42   1.95660e+02   2.411e+06   8.006e+05   0.000e+00       3.212e+06
    43   1.95814e+02   2.403e+06   8.000e+05   0.000e+00       3.203e+06
    44   1.95968e+02   2.395e+06   7.996e+05   0.000e+00       3.195e+06
    45   1.96120e+02   2.389e+06   7.993e+05   0.000e+00       3.188e+06
    46   1.96269e+02   2.383e+06   7.990e+05   0.000e+00       3.182e+06

Total time (s) = 2.50

Let’s now compare the final results as well as the convergence curves of the two algorithms. We can see how the adaptive Primal-Dual produces a better estimate of the clean image in a much smaller number of iterations

fig, axs = plt.subplots(1, 4, figsize=(16, 4))
axs[0].imshow(img, cmap='gray', vmin=0, vmax=255)
axs[1].imshow(noise_img, cmap='gray', vmin=0, vmax=255)
axs[2].imshow(iml12_fixed, cmap='gray', vmin=0, vmax=255)
axs[3].imshow(iml12_ada, cmap='gray', vmin=0, vmax=255)
axs[3].set_title('Adaptive PD')

fig, axs = plt.subplots(2, 1, figsize=(12, 7))
axs[0].plot(cost_fixed, 'k', label='Fixed step')
axs[0].plot(cost_ada, 'r', label='Adaptive step')
axs[1].plot(err_fixed, 'k', label='Fixed step')
axs[1].plot(err_ada, 'r', label='Adaptive step')

fig, axs = plt.subplots(3, 1, figsize=(12, 7))
axs[0].plot(steps[0], 'k')
axs[1].plot(steps[1], 'k')
axs[2].plot(steps[2], 'k')
  • Original, Noisy, PD, Adaptive PD
  • Functional, MSE
  • $\tau^k$, $\mu^k$, $\alpha^k$

Total running time of the script: (0 minutes 17.188 seconds)

Gallery generated by Sphinx-Gallery