Convolutional Sparse Coding (CSC)

CSC is defined by the optimization problem

\[ \begin{eqnarray*} \min_{ \{ \mathbf{u}_{l}\} } & & \frac{1}{2} \left\| \sum_{l} h^{(l)} * \mathbf{u}^{(l)} - \mathbf{b} \right\|_2^2 + \lambda \sum_{l} \| \mathbf{u}^{(l)} \|_1, \label{eqn:csc}\tag{1} \end{eqnarray*} \]

where

\(\{ h^{(l)} \}\) is the set of convolutional dictionaries,
\(\mathbf{b}\) is the input image, and
\(\{ \mathbf{u}^{(l)} \}\) is the corresponding set of coefficient maps that encodes image \(\mathbf{b}\) via \(\{ H_l \}\).

Problem (1) is a particular case of composite unconstrained convex programming (see Section 5.2 [BCG11]), which in general can be written as

\[ \begin{eqnarray*} \min_{ \{ \mathbf{u}\} } & & f( \mbox{Op}(\mathbf{u}) ) + \lambda g( \mathbf{u} ), \label{eqn:genProb}\tag{2} \end{eqnarray*} \]

and be solved via several numerical algorithms, e.g. Douglas-Rachford splitting [EB92], forward-backward splitting [CW05], ADMM [BPC+11], FISTA [BT09], etc.

Solving (1) using the F2O's framework in the frequency domain

Summary

\[\begin{split} \begin{array}{l} \hline \; \\ \qquad f( \mbox{Op}(\mathbf{u}) ) + \lambda g( \mathbf{u} ) \; \longleftrightarrow \; 0.5\cdot\| \mbox{Op}(\mathbf{u}) - \mathbf{b}\|_2^2 + \lambda \cdot \| \mathbf{u} \|_1 \\ \; \\ % --- \begin{array}{clll} \mbox{Op}(\mathbf{u}) & = & \displaystyle \sum_{l=1}^L h^{(l)} * \mathbf{u}^{(l)} \\ f(\cdot) & = & 0.5 \| \cdot \|_2^2 & \color{darkgray} \mbox{quadratic cost function} \\ g(\cdot) & = & \| \cdot \|_1 & \color{darkgray} \mbox{sparsity-inducing penalty} \\ \mathop{\mathrm{prox}}_{g,\lambda}(\mathbf{z}) & = & \displaystyle \mathop{\mathrm{argmin}}_{\mathbf{x}} \; \frac{1}{2} \| \mathbf{x} - \mathbf{z} \| + g( \mathbf{x} ) \\ & = & \mbox{soft}(\mathbf{z}, \lambda) \end{array} \\ \\ \hline \mathbf{\mbox{Inputs}} \color{white}{ \Large A} \\ % --- \begin{array}{llll} \mathbf{v} & = & \mathbf{v}_{LP} + \mathbf{b} & \color{darkgray} \mbox{Input image} \\ \mathbf{b} & : & \mbox{high-pass version of }\mathbf{v} & \color{darkgray}\left(\mbox{Input to } (\ref{eqn:csc}), B_F = {\cal{F}}\{b\} \right) \\ \{ h^{(l)} \} & : & \mbox{2D filter bank} & \color{darkgray} \left( H^{(l)} = {\cal{F}}\{h^{(l)}\} \right) \\ {\cal{F}}\{\cdot\} & : & \mbox{Fourier transform} & \color{darkgray} \mbox{(includes a zero-pad argument)} \\ & \end{array} \\ \hline % % --- \end{array} \end{split}\]

About \(\{ \mbox{u}^{(l)} \}\) and \(\{ h^{(l)} \}\)

\[\begin{split} \begin{array}{l} \hline \mathbf{\mbox{Coefficient map}} \color{white}{ \Large A} \\ % --- \begin{array}{llll} \{\mathbf{u}^{(l)}_0 \} & : & \mbox{Initial guess} & \quad \\ U_{F,0} & = & {\cal{F}}\left\{ \{\mathbf{u}^{(l)}_0 \} \right\} = \left[ \begin{array}{c} {\cal{F}}\{ \mathbf{u}_0^{(0)} \} \\ {\cal{F}}\{ \mathbf{u}_0^{(1)} \} \\ \vdots \\ {\cal{F}}\{ \mathbf{u}_0^{(L)} \} \end{array} \right] & \color{darkgray}\begin{array}{l} (\mbox{may include boundary} \\ \mbox{conditions; see this} ) \end{array}\\ \end{array} \\ \hline \mathbf{\mbox{Filterbank}} \color{white}{ \Large A} \\ \begin{array}{llll} H^{(l)} & = & {\cal{F}}\{ h^{(l)}, \mbox{shape}=\cdot \} & \color{darkgray} \begin{array}{l} ( \mbox{zero-pad Fourier transform}; \\ \mbox{it matches }\mathbf{u}^{(l)}\mbox{'s shape.} )\end{array}\\ H_F & = & \left[ \begin{array}{c} H^{(1)} \\ H^{(2)} \\ \vdots \\ H^{(L)} \end{array} \right] \end{array} \\ \hline % % --- \end{array} \end{split}\]

Op\((\cdot)\)’s adjoint

\[\begin{split} \begin{array}{l} \hline \mathbf{\mbox{Spatial domain}} \color{white}{ \Large A} \\ % %--- \begin{array}{llll} \mbox{Op}(\mathbf{u}) & = & \displaystyle \sum_{l=1}^L h^{(l)} * \mathbf{u}^{(l)} & \\ % \color{darkgray}\begin{array}{l} ( \mbox{recall that }\mathbf{u} \mbox{ represents} \\ % \mbox{the coefficient maps} ) \end{array} \\ % {\cal{A}}\{ \mbox{Op}(\mathbf{z}) \} & = & \left[ \begin{array}{c} g^{(1)} * \mathbf{z} \\ g^{(2)} * \mathbf{z} \\ \vdots \\ g^{(L)} * \mathbf{z} \end{array} \right] & \color{darkgray}\begin{array}{l} \bullet\; \mathbf{z}\mbox{'s shape is equal to }\mathbf{b}\mbox{'s} \\ \bullet\; g^{(l)} = {\cal{A}}\{ h^{(l)} \} =\, \scriptstyle\texttt{np.flip(np.flip(h,1),0)} \\ \phantom{\bullet}\; \mbox{(adjoint of filter } h^{(l)} ) \end{array} \end{array} \\ \; \\ \hline %--- \mathbf{\mbox{Frequency domain}} \color{white}{ \Large A} \\ % \begin{array}{llll} \mbox{Op}(U_F) & = & H_F \odot U_F = \displaystyle \sum_{l=1}^L H^{(l)}_F \odot U_F^{(l)} & \quad \\ % {\cal{A}}\{ \mbox{Op}(Z_F) \} & = & H_F^* \odot Z_F = \left[ \begin{array}{c} \left\{ H^{(1)}_F \right\}^* \odot Z_F \\ \left\{ H^{(2)}_F \right\}^* \odot Z_F\\ \vdots \\ \left\{ H^{(L)}_F \right\}^* \odot Z_F \end{array} \right] & \color{darkgray} \{ \cdot \}^* : \mbox{complex conjugate}\\ \end{array} \\ % --- \; \\ \hline \end{array} \end{split}\]

Gradient \(\nabla f\)

\[\begin{split} \begin{array}{l} \hline \mathbf{\mbox{Spatial domain}} \color{white}{ \Large A} \\ %--- \begin{array}{llll} \nabla f(\mathbf{u}) & = & \left[ \begin{array}{c} g^{(1)} * \mathbf{z} \\ g^{(2)} * \mathbf{z} \\ \vdots \\ g^{(L)} * \mathbf{z} \end{array}\right] \quad & \color{darkgray} \begin{array}{l} \bullet\; \mathbf{z} = \mbox{Op}(\mathbf{u}) - \mathbf{b} \\ \bullet\; g^{(l)} = {\cal{A}}\{ h^{(1)} \} \end{array} \end{array} \\ \\ \hline %--- \mathbf{\mbox{Frequency domain}} \color{white}{ \Large A} \\ %--- \begin{array}{llll} \nabla f( U_F ) & = & H_F^* \odot (H_F \odot U_F - B_F) & \color{darkgray} \{ \cdot \}^* : \mbox{complex conjugate} \end{array} \\ \\ \hline %--- \end{array} \end{split}\]

FISTA

\[\begin{split} \begin{array}{l} \hline \\ % --- \mathbf{\mbox{FISTA}} \qquad \color{darkgray} \mbox{Init :} \left\{ \begin{array}{rcl} U_{F.0} & : & \mbox{Initial guess} \\ Y_F & = & U_{F,0} \\ \alpha_k & : & \mbox{step-size policy (includes cte. case)} \\ \gamma_k & : & \mbox{inertial sequence} \end{array} \right. \\ \\ \hline % \\ % \mathbf{\mbox{for }} k=0,1,\ldots \\ \begin{array}{llll} & \nabla f & = & H_F^* \odot ( H_F \odot Y_F - B_F ) \\ % & \alpha_k & = & \mbox{step-size policy (includes cte. case)} \\ & \mathbf{x}^{(l)} & = & {\cal{F}^{-1}} \{ Y_F - \alpha_k \nabla f \} \\ & U_{F,k+1} & = & {\cal{F}} \left\{ \mbox{soft}\left(\mathbf{x}^{(l)}, \alpha_k \cdot \lambda \right) \right\} \\ & Y_F & = & U_{F,k+1} + \gamma_k \cdot \left( U_{F,k+1} - U_{F,k} \right) \end{array} % % --- \end{array}\end{split}\]

ADMM

\[\begin{split} \begin{array}{l} \hline \\ % --- \mathbf{\mbox{ADMM}} \qquad \color{darkgray} \mbox{Init :} \left\{ \begin{array}{rcl} U_{F.0} & : & \mbox{Initial guess} \\ Z_F, C_F & : & \mbox{ADMM variables} \\ \rho_k & : & \mbox{ADMM's parameter} \\ \beta & : & \mbox{relaxation parameter} \end{array} \right. \\ \\ \hline % \\ % \mathbf{\mbox{for }} k=0,1,\ldots \\ \begin{array}{llll} % & U_{F,k+1} & = & \displaystyle \frac{ H_F^* \odot B_F + \rho_k \cdot (C_{F,k} - Z_{F,k}) }{ H_F^* \odot H_F + \rho_k } & \\ & \mathbf{x}^{(l)} & = & {\cal{F}^{-1}} \{ \beta\cdot U_{F,k+1} - (1-\beta)\cdot Z_{F,k}\} \\ & Z_{F,k+1} & = & {\cal{F}} \left\{ \mbox{soft}\left(\mathbf{x}^{(l)}, \lambda/ \rho_k \right) \right\} \\ & C_{F,k+1} & = & C_{F,k} + \beta\cdot U_{F,k+1} - (1-\beta)\cdot Z_{F,k} - Z_{F,k+1} & \\ \end{array} % % --- \end{array}\end{split}\]

Code snippet

The corresponding code snippets for FISTA and ADMM are access via the corresponding tabs.

FISTA

  for k in range(nIter):
        
    grad = computeGrad(y, args.Bf, x_sptl)      # compute the gradient
                                                # Depends on the definition of f()
                                                
    alpha, grad = computeSS(k, grad, y, x_sptl) # compute step-size

    x0 = x1.copy()

    x_sptl = prox(irfft2( y - alpha*grad, axes=(0,1) ), alpha*lmbda)     
    x1     = rfft2( x_sptl, axes=(0,1) ).ravel()
        
    t1 = iseq(k, t0)                            # compute inertial sequence
    gamma = (t0-1)/t1
    
    y = x1 + gamma*(x1-x0)                      # extrapolation step

    t0 = t1

ADMM

ToDo

Simulation

The collapsible items below sets up a simulation to test diferent algorithms that target the solution of (1).

Simulation setup

1. Load F2O

# Generic F2O imports
import F2O.F2O_utils as F2O
import F2O.fwOp.fwOperator as fwOp

# F2O 'utils' for reading / ploting images
from F2O.imgUtils.image_utils import ImgPlot, ImgRead

# F2O 'utils' for (i) adding noise, (ii) collecting image metrics and (iii) applying the forward model 
from F2O.imgUtils.image_utils import ImgMetrics, ImgApplyFwNoise
from F2O.noise.noiseModels import noiseModels

# Other packages
import matplotlib.pylab as PLT
import numpy as np

# If you get an error while loading the F2O, then
# 
# * Exit Jupyter
# * Go to the F2O root dir, and execute
#      export PYTHONPATH=$PYTHONPATH:`pwd`
# * Relaunch Jupyter

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
/tmp/ipykernel_251/1396103323.py in <module>
      1 # Generic F2O imports
----> 2 import F2O.F2O_utils as F2O
      3 import F2O.fwOp.fwOperator as fwOp
      4 
      5 # F2O 'utils' for reading / ploting images

ModuleNotFoundError: No module named 'F2O'

2. Load test images (SIPI Image database)

# Imports needed to read images from an URL address
import requests
from io import BytesIO

# Define an 'F2O-image' object 
testImgs = ImgRead()
        
# Test images from the SIPI image database
fname = {0: requests.get('http://sipi.usc.edu/database/misc/5.2.10.tiff'),   # bridge (grayscale)
         1: requests.get('http://sipi.usc.edu/database/misc/boat.512.tiff'), # boat (grayscale)
         2: requests.get('http://sipi.usc.edu/database/misc/4.2.03.tiff'),   # mandrill (color)
        }

testImgs.list.append([BytesIO(fname.get(1).content),'g'])
testImgs.list.append([BytesIO(fname.get(2).content),'c'])

u = testImgs.readListImgs()   # read list of images, normalize them between 0 and 1

pltImg  = ImgPlot()
pltImg.plotNImgs(u, len(u), None, 5)

3. Compute low / high pass filtered image; if any, add noise.

# Noise model

noise = noiseModels(model=F2O.f2oDef.noise_Gaussian, sigma=0.0)
addNoise = noise.sel_NoiseModel()

# Select one image from list
uObs = addNoise(u[0])

uLP, uHP = testImgs.computeLPHP( uObs , padFlag=True)

pltImg  = ImgPlot()
pltImg.plotNImgs(uHP, len(uHP), None, 5)

4. Load convolutional dictionary.

# Load library to read .mat files
import scipy.io as sio

# Load pre-computed dictionary
Dict = sio.loadmat('../data/nsp_36x36x12.mat')  
H = Dict['FB']

FISTA-based simulation

1. Define optimization variable; set properties for FISTA

# Load F2O specific modules
from F2O.fwOp.fwOperator import fwOp_f
from F2O.PG_freq import apg as fAPG

# Create an 'optimization variable' in freq.; set generic properties
argsF = F2O.argsF2O()

argsF.verbose         = False
argsF.fCostClass      = argsF.f2oDef.cost_L2L1_lin   # F(x) = 0.5|| Op(x) - b ||_2^2 + \lambda||x||_1, 
                                                     # where Op(.) is linear
argsF.padFlag         = True
argsF.padMode         = 'symmetric'


# Set Op(.) -- in this particular case, the convolutional dictionary

freqOp         = fwOp_f(linOp=argsF.f2oDef.fAx_FB2D, A=H)

freqOp.label     = "nsp_36x36x12"
freqOp.vecFlag   = True            # input / output data are asumed / force to be vectorized

1. Execute FISTA

# Set particular FISTA properties
argsF.ssPolicy   = argsF.f2oDef.ss_CauchyLagged     # step-size policy
argsF.ISeqPolicy = argsF.f2oDef.iseq_ntrv           # inertial sequence

lmbda = 0.1
nIter = 40

x, Stats = fAPG(freqOp, uHP, lmbda, nIter, argsF)


fig = PLT.figure(figsize=(24, 16))
ax1 = fig.add_subplot(2, 1, 1)
PLT.plot(Stats[:,2], Stats[:,0], label=r'$\alpha_k$ : {0}'.format(argsF.f2oDef.ss_list[argsF.ssPolicy]) )

PLT.legend(loc='upper right',fontsize=20)
PLT.ylabel(r'$f(x) = \frac{1}{2} \|\|$Op$(\mathbf{u}) - \mathbf{b} \|\|_2^2$',fontsize=20)
PLT.xlabel('Time',fontsize=20);

imgs = []
uRec = freqOp.fRec(x, uLP, getROI=argsF.getROI )
imgs.append(uRec)
imgs.append(uRec[:,:,0]-uObs)

txt = []
txt.append('Reconstructed')
txt.append('Difference {:1.2f}'.format(np.linalg.norm(uRec.ravel() - uObs.ravel()) ) )

pltImg  = ImgPlot()
pltImg.plotNImgs(imgs, len(imgs), txt, 5)