Constained MOD

This example demonstrates the use of admm.cmod.CnstrMOD for computing a dictionary update via a constrained variant of the method of optimal directions [20]. This problem is mainly useful as a component within dictionary learning, but its use is demonstrated here since a user may wish to construct such objects as part of a custom dictionary learning algorithm, using dictlrn.DictLearn.

from __future__ import print_function
from builtins import input

import numpy as np

from sporco.admm import bpdn
from sporco.admm import cmod
from sporco import util
from sporco import array
from sporco import plot
plot.config_notebook_plotting()

Load training images.

exim = util.ExampleImages(scaled=True, zoom=0.25, gray=True)
S1 = exim.image('barbara.png', idxexp=np.s_[10:522, 100:612])
S2 = exim.image('kodim23.png', idxexp=np.s_[:, 60:572])
S3 = exim.image('monarch.png', idxexp=np.s_[:, 160:672])
S4 = exim.image('sail.png', idxexp=np.s_[:, 210:722])
S5 = exim.image('tulips.png', idxexp=np.s_[:, 30:542])

Extract all 8x8 image blocks, reshape, and subtract block means.

S = array.extract_blocks((S1, S2, S3, S4, S5), (8, 8))
S = np.reshape(S, (np.prod(S.shape[0:2]), S.shape[2]))
S -= np.mean(S, axis=0)

Load initial dictionary.

D0 = util.convdicts()['G:8x8x64']
D0 = np.reshape(D0, (np.prod(D0.shape[0:2]), D0.shape[2]))

Compute sparse representation on current dictionary.

lmbda = 0.1
opt = bpdn.BPDN.Options({'Verbose': True, 'MaxMainIter': 200,
                         'RelStopTol': 1e-3})
b = bpdn.BPDN(D0, S, lmbda, opt)
X = b.solve()
Itn   Fnc       DFid      Regℓ1     r         s         ρ
----------------------------------------------------------------
   0  2.68e+04  2.07e+04  6.15e+04  3.57e-01  2.33e+00  6.00e+00
   1  2.48e+04  1.74e+04  7.39e+04  1.49e-01  1.02e+00  6.00e+00
   2  2.40e+04  1.62e+04  7.80e+04  8.52e-02  6.64e-01  6.00e+00
   3  2.35e+04  1.54e+04  8.06e+04  6.04e-02  5.31e-01  6.00e+00
   4  2.32e+04  1.50e+04  8.21e+04  4.53e-02  4.28e-01  6.00e+00
   5  2.30e+04  1.46e+04  8.32e+04  3.59e-02  3.68e-01  6.00e+00
   6  2.28e+04  1.44e+04  8.40e+04  2.92e-02  3.17e-01  6.00e+00
   7  2.27e+04  1.42e+04  8.46e+04  2.45e-02  2.81e-01  6.00e+00
   8  2.26e+04  1.41e+04  8.51e+04  2.08e-02  2.50e-01  6.00e+00
   9  2.25e+04  1.40e+04  8.55e+04  1.80e-02  2.26e-01  6.00e+00
  10  2.23e+04  1.35e+04  8.74e+04  6.77e-02  1.82e-01  1.20e+00
  11  2.22e+04  1.34e+04  8.81e+04  5.10e-02  1.30e-01  1.20e+00
  12  2.21e+04  1.33e+04  8.81e+04  3.72e-02  9.44e-02  1.20e+00
  13  2.21e+04  1.33e+04  8.81e+04  2.77e-02  7.19e-02  1.20e+00
  14  2.21e+04  1.32e+04  8.83e+04  2.11e-02  5.61e-02  1.20e+00
  15  2.21e+04  1.32e+04  8.84e+04  1.65e-02  4.48e-02  1.20e+00
  16  2.21e+04  1.32e+04  8.85e+04  1.30e-02  3.64e-02  1.20e+00
  17  2.20e+04  1.32e+04  8.85e+04  1.05e-02  3.00e-02  1.20e+00
  18  2.20e+04  1.32e+04  8.86e+04  8.55e-03  2.51e-02  1.20e+00
  19  2.20e+04  1.32e+04  8.86e+04  7.07e-03  2.13e-02  1.20e+00
  20  2.20e+04  1.32e+04  8.87e+04  1.36e-02  1.74e-02  4.89e-01
  21  2.20e+04  1.32e+04  8.88e+04  1.14e-02  1.34e-02  4.89e-01
  22  2.20e+04  1.31e+04  8.88e+04  9.16e-03  1.00e-02  4.89e-01
  23  2.20e+04  1.31e+04  8.89e+04  7.23e-03  7.61e-03  4.89e-01
  24  2.20e+04  1.31e+04  8.89e+04  5.65e-03  5.88e-03  4.89e-01
  25  2.20e+04  1.31e+04  8.89e+04  4.42e-03  4.61e-03  4.89e-01
  26  2.20e+04  1.31e+04  8.89e+04  3.48e-03  3.66e-03  4.89e-01
  27  2.20e+04  1.31e+04  8.89e+04  2.75e-03  2.94e-03  4.89e-01
  28  2.20e+04  1.31e+04  8.89e+04  2.20e-03  2.38e-03  4.89e-01
  29  2.20e+04  1.31e+04  8.89e+04  1.76e-03  1.95e-03  4.89e-01
  30  2.20e+04  1.31e+04  8.89e+04  1.95e-03  1.59e-03  3.29e-01
  31  2.20e+04  1.31e+04  8.89e+04  1.64e-03  1.26e-03  3.29e-01
  32  2.20e+04  1.31e+04  8.89e+04  1.36e-03  9.99e-04  3.29e-01
  33  2.20e+04  1.31e+04  8.89e+04  1.11e-03  7.99e-04  3.29e-01
  34  2.20e+04  1.31e+04  8.89e+04  9.05e-04  6.45e-04  3.29e-01
----------------------------------------------------------------

Update dictionary for training image set.

opt = cmod.CnstrMOD.Options({'Verbose': True, 'MaxMainIter': 100,
                             'RelStopTol': 1e-3, 'rho': 4e2})
c = cmod.CnstrMOD(X, S, None, opt)
D1 = c.solve()
print("CMOD solve time: %.2fs" % c.timer.elapsed('solve'))
Itn   DFid      Cnstr     r         s         ρ
------------------------------------------------------
   0  1.13e+04  6.86e-07  4.56e-01  2.38e+00  4.00e+02
   1  1.12e+04  7.37e-07  2.08e-01  2.09e-01  4.00e+02
   2  1.12e+04  7.79e-07  1.01e-01  7.24e-02  4.00e+02
   3  1.12e+04  7.73e-07  4.97e-02  3.14e-02  4.00e+02
   4  1.12e+04  7.00e-07  2.72e-02  1.44e-02  4.00e+02
   5  1.12e+04  6.73e-07  1.53e-02  8.43e-03  4.00e+02
   6  1.12e+04  7.52e-07  8.87e-03  4.74e-03  4.00e+02
   7  1.12e+04  7.31e-07  5.22e-03  3.04e-03  4.00e+02
   8  1.12e+04  8.13e-07  3.12e-03  1.91e-03  4.00e+02
   9  1.12e+04  7.16e-07  1.88e-03  1.29e-03  4.00e+02
  10  1.12e+04  7.56e-07  1.15e-03  8.55e-04  4.00e+02
  11  1.12e+04  7.63e-07  7.08e-04  5.95e-04  4.00e+02
------------------------------------------------------
CMOD solve time: 0.11s

Display initial and final dictionaries.

D0 = D0.reshape((8, 8, D0.shape[-1]))
D1 = D1.reshape((8, 8, D1.shape[-1]))
fig = plot.figure(figsize=(14, 7))
plot.subplot(1, 2, 1)
plot.imview(util.tiledict(D0), title='D0', fig=fig)
plot.subplot(1, 2, 2)
plot.imview(util.tiledict(D1), title='D1', fig=fig)
fig.show()
../../_images/cmod_13_0.png

Get iterations statistics from CMOD solver object and plot functional value, ADMM primary and dual residuals, and automatically adjusted ADMM penalty parameter against the iteration number.

its = c.getitstat()
fig = plot.figure(figsize=(20, 5))
plot.subplot(1, 3, 1)
plot.plot(its.DFid, xlbl='Iterations', ylbl='Functional', fig=fig)
plot.subplot(1, 3, 2)
plot.plot(np.vstack((its.PrimalRsdl, its.DualRsdl)).T,
          ptyp='semilogy', xlbl='Iterations', ylbl='Residual',
          lgnd=['Primal', 'Dual'], fig=fig)
plot.subplot(1, 3, 3)
plot.plot(its.Rho, xlbl='Iterations', ylbl='Penalty Parameter', fig=fig)
fig.show()
../../_images/cmod_15_0.png