Distributed Convolutional Dictionary Learning (DiCoDiLe): Pattern Discovery in Large Images and Signals

TL;DR

This paper introduces DiCoDiLe, a distributed and asynchronous convolutional dictionary learning algorithm that efficiently discovers patterns in large-scale images and signals, demonstrated on Hubble Space Telescope images.

Contribution

It presents a novel distributed, asynchronous optimization algorithm for convolutional dictionary learning that scales linearly with data size, enabling pattern discovery in very large images.

Findings

Algorithm scales linearly with data size

Effective pattern discovery on large images

Validated on Hubble Space Telescope images

Abstract

Convolutional dictionary learning (CDL) estimates shift invariant basis adapted to multidimensional data. CDL has proven useful for image denoising or inpainting, as well as for pattern discovery on multivariate signals. As estimated patterns can be positioned anywhere in signals or images, optimization techniques face the difficulty of working in extremely high dimensions with millions of pixels or time samples, contrarily to standard patch-based dictionary learning. To address this optimization problem, this work proposes a distributed and asynchronous algorithm, employing locally greedy coordinate descent and an asynchronous locking mechanism that does not require a central server. This algorithm can be used to distribute the computation on a number of workers which scales linearly with the encoded signal's size. Experiments confirm the scaling properties which allows us to learn patterns on large scales images from the Hubble Space Telescope.

Figures

Equations

$$(X\ast \mathbf{Y})[\omega ]=\sum _{\tau \in \Omega }X[\omega -\tau ]\cdot \mathbf{Y}[\tau ],\forall \omega \in \Omega .$$

$$\text{ST}(u,\lambda )=\text{sign}(u)\max (|u|-\lambda ,0).$$

$$X=Z\ast \mathbf{D}+\xi ,$$

$$\underset{\begin{array}{c}Z,\mathbf{D}\\ \Vert {\mathbf{D}}_k\Vert \le 1\end{array}}{\min }\underset{F(Z,\mathbf{D})}{\underbrace{\frac{1}{2}{\Vert X-Z\ast \mathbf{D}\Vert }_2^2}}+\underset{G(Z)}{\underbrace{\lambda \Vert Z{\Vert }_1}}.$$

$$Z^{\star }=\underset{Z}{\mathrm{argmin}}\frac{1}{2}{\Vert X-Z\ast \mathbf{D}\Vert }_2^2+\lambda {\Vert Z\Vert }_1.$$

$${\nabla }_{Z}F(0,\mathbf{D})=X\ast {\mathbf{D}}^{\Lsh }\in \partial G(0)=[-\lambda ,\lambda {]}^{K|\Omega |}.$$

$${\mathbf{D}}^{\ast }=\underset{\begin{array}{c}\mathbf{D}\\ \Vert {\mathbf{D}}_k{\Vert }_2\le 1\end{array}}{\mathrm{argmin}}\frac{1}{2}{\Vert X-Z\ast \mathbf{D}\Vert }_2^2.$$

$$Z_{k_0}^{\prime }[{\omega }_0]=\frac{1}{\Vert {\mathbf{D}}_{k_0}{\Vert }_2^2}\text{ST}({\beta }_{k_0}^{(q)}[{\omega }_0],\lambda ),$$

$${\beta }_k^{(q)}[\omega ]=\left(\left(X-Z^{(q)}\ast \mathbf{D}+Z_k^{(q)}[\omega ]e_{\omega }\ast {\mathbf{D}}_k\right)\ast {\mathbf{D}}_k^{\Lsh }\right)[\omega ],$$

$${\beta }_k^{(q+1)}[\omega ]={\beta }_k^{(q)}[\omega ]-({\mathbf{D}}_{k_0}\ast {\mathbf{D}}_k^{\Lsh })[\omega -{\omega }_0]\Delta Z_{k_0}^{(q)}[{\omega }_0],$$

$$\mathcal{V}({\omega }_0)=\prod _{i=1}^d[({\omega }_0{)}_i-L_i+1,({\omega }_0{)}_i+L_i[.$$

$$(k_0,{\omega }_0)=\underset{(k,\omega )\in [1,K]\times {\mathcal{C}}_m}{\mathrm{argmax}}|\Delta Z_k[\omega ]|,$$

$${\mathcal{B}}_L({\mathcal{S}}_w)=\prod _{i=1}^d[l_i,l_i+L_i[\cup [u_i-L_i,u_i[.$$

$$(k_0,{\omega }_0)=\underset{(k,\omega )\in [1,K]\times {\mathcal{C}}_m^{(w)}}{\mathrm{argmax}}|\Delta Z_k[\omega ]|,$$

$$\Delta$$

$$-\sum _{w\ne w^{\prime }}({\mathbf{D}}_{k_w}\ast {\mathbf{D}}_{k_{w^{\prime }}}^{\Lsh })[{\omega }_{w^{\prime }}-{\omega }_w]\Delta Z_{k_w}[{\omega }_w]\Delta Z_{k_{w^{\prime }}}[{\omega }_{w^{\prime }}],$$

$${\mathcal{E}}_L({\mathcal{S}}_w)=\prod _{i=1}^d[l_i-L_i,l_i[\cup [u_i,u_i+L_i[.$$

$$|\Delta Z_{k_0}[{\omega }_0]|>\underset{k,\omega \in [1,K]\times \mathcal{V}({\omega }_0)\cap {\mathcal{E}}_L({\mathcal{S}}_w)}{\max }|\Delta Z_k[\omega ]|.$$

$$P(\omega \notin \text{SL})\ge \prod _{i=1}^d(1-\frac{W_i{L}_i}{T_i})$$

$${\nabla }_{\mathbf{D}}F(Z,\mathbf{D})=Z^{\Lsh }\ast (X-Z\ast \mathbf{D})=\psi -\varphi \ast \mathbf{D},$$

$$\varphi [\tau ]=\sum _{\omega \in \Omega }Z[\omega ]Z[\tau +\omega ]=\sum _{w=1}^W\sum _{\omega \in {\mathcal{S}}_w}Z[\omega ]Z[\tau +\omega ]$$

$$\Delta E_{k_0}[{\omega }_0]$$

$${\alpha }_{k_0}[\omega ]=(X-Z\ast \mathbf{D})[\omega ]+{\mathbf{D}}_{k_0}[\omega -{\omega }_0]Z_{k_0}[{\omega }_0]=(X-Z\ast \mathbf{D})[\omega ]+Z_{k_0}[{\omega }_0](e_{{\omega }_0}\ast {\mathbf{D}}_{k_0})[{\omega }_0].$$

$$\Delta E_{k_0}[{\omega }_0]=$$

$$=$$

$$Z_{k_0}^{\prime }[{\omega }_0]=\arg \underset{y\in \mathbb{R}}{\min }{e}_{k_0,{\omega }_0}(u)\sim \arg \underset{u\in \mathbb{R}}{\min }\frac{\Vert {\mathbf{D}}_{k_0}{\Vert }_2^2}{2}{\left(u-{\beta }_{k_0}[{\omega }_0]\right)}^2+\lambda |u|.$$

$$\Delta E=\stackrel{\text{iterative steps}}{\overbrace{\sum _{i=1}^W\Delta E_w}}-\underset{\text{interference}}{\underbrace{\sum _{w\ne w^{\prime }}({\mathbf{D}}_{k_w}\ast {\mathbf{D}}_{k_{w^{\prime }}}^{\Lsh })[{\omega }_{w^{\prime }}-{\omega }_w]\Delta Z_{k_w}[{\omega }_w]\Delta Z_{k_{w^{\prime }}}[{\omega }_{w^{\prime }}]}},$$

$$\alpha [\omega ]=(X-Z\ast \mathbf{D})[t]+\sum _{w=1}^K{\mathbf{D}}_{k_w}[\omega -{\omega }_w]Z_{k_w}[{\omega }_w]=(X-Z^{(1)}\ast \mathbf{D})[t]+\sum _{w=1}^K{\mathbf{D}}_{k_w}[\omega -{\omega }_w]Z_{k_w}^{\prime }[{\omega }_w],$$

$${\alpha }_w[\omega ]=(X-Z\ast \mathbf{D})[t]+{\mathbf{D}}_{k_w}[\omega -{\omega }_w]Z_{k_w}[{\omega }_w]=\alpha [\omega ]-\sum _{w^{\prime }\ne w}{\mathbf{D}}_{k_{w^{\prime }}}[\omega -{\omega }_{w^{\prime }}]Z_{k_{w^{\prime }}}[{\omega }_{w^{\prime }}].$$

$$\Delta E=$$

Full text

1 Introduction

The sparse convolutional linear model has been used successfully for various signal and image processing applications. It was first developed by Grosse et al. (2007) for music recognition and has then been used extensively for denoising and inpainting of images (Kavukcuoglu et al., 2010; Bristow et al., 2013; Wohlberg, 2014). The benefit of this model over traditional sparse coding technique is that it is shift-invariant. Instead of local patches, it considers the full signals which allows for a sparser reconstruction, as not all shifted versions of the same patches are present in the estimated dictionary (Bergstra et al., 2011).

Beyond denoising or inpainting, recent applications used this shift-invariant model as a way to learn and localize important patterns in signals or images. Yellin et al. (2017) used convolutional sparse coding to count the number of blood-cells present in holographic lens-free images. Jas et al. (2017) and Dupré la Tour et al. (2018) learned recurring patterns in univariate and multivariate signals from neuroscience. This model has also been used with astronomical data as a way to denoise (Starck et al., 2005) and localize predefined patterns from space satellite images (del Aguila Pla and Jalden, 2018).

Convolutional dictionary learning estimates the parameters of this model. It boils down to an optimization problem, where one needs to estimate the patterns (or atoms) of the basis as well as their activations, i.e. where the patterns are present in the data. The latter step, which is commonly referred to as convolutional sparse coding (CSC), is the critical bottleneck when it comes to applying CDL to large data. Chalasani et al. (2013) used an accelerated proximal gradient method based on Fast Iterative Soft-Thresholding Algorithm (FISTA; Beck and Teboulle 2009 that was adapted for convolutional problems. The Fast Convolutional Sparse Coding of Bristow et al. (2013) is based on Alternating Direction Method of Multipliers (ADMM). Both algorithms compute a convolution on the full data at each iteration. While it is possible to leverage the Fast Fourier Transform to reduce the complexity of this step (Wohlberg, 2016), the cost can be prohibitive for large signals or images. To avoid such costly global operations, Moreau et al. (2018) proposed to use locally greedy coordinate descent (LGCD) which requires only local computations to update the activations. Based on Greedy Coordinate Descent (GCD), this method has a lower per-iteration complexity and it relies on local computation in the data, like a patch based technique would do, yet solving the non-local problem.

To further reduce the computational burden for CDL, recent studies consider different parallelization strategies. Skau and Wohlberg (2018) proposed a method based on consensus ADMM to parallelize the CDL optimization by splitting the computation across the different atoms. This limits the number of cores that can be used to the number of atoms in the dictionary, and does not significantly reduce the burden for very large signals or images. An alternative parallel algorithm called DICOD, proposed by Moreau et al. (2018) distributes GCD updates when working with 1d signals. Their results could be used with multidimensional data such as images, yet only by splitting the data along one dimension. This limits the gain of their distributed optimization strategy when working with images. Also, the choice of GCD for the local updates make the algorithm inefficient when used with low number of workers.

In the present paper, we propose a distributed solver for CDL, which can fully leverage multidimensional data. For the CSC step, we present a coordinate selection rule based on LGCD that works in the asynchronous setting without a central server. We detail a soft-lock mechanism ensuring the convergence with a grid of workers, while splitting the data in all convolutional dimensions. Then, we explain how to also leverage the distributed set of workers to further reduce the cost of the dictionary update, by precomputing in parallel some quantities used in the gradient operator. Extensive numerical experiments confirm the performance of our algorithm, which is then used to learn patterns from a large astronomical image.

In Section 2, we describe the convolutional linear model and its parameter estimation through CDL. Then, Section 3 details the LGCD algorithm for convolutional sparse coding (CSC). Our distributed algorithm DiCoDiLe is introduced in Section 4. Section 5 presents results on simulations and images.

2 Convolutional Dictionary Learning

Notations.

For a vector $U\in{\mathbb{R}}^{P}$ we denote $U_{p}\in{\mathbb{R}}$ its $p^{th}$ coordinate. For finite domain $\Omega=\prod_{i=1}^{d}[0,T_{i}[$, where $[0,T_{i}[$ are the integers from [math] to $T_{i}-1$, $|\Omega|$ is the size of this domain $\prod_{i=1}^{d}T_{i}$. We denote ${\mathcal{X}}_{\Omega}^{P}$ (resp. ${\mathcal{X}}_{\Omega}^{P\times Q}$) the space of observations defined on $\Omega$ with values in ${\mathbb{R}}^{P}$ (resp. ${\mathbb{R}}^{P\times Q}$). For instance, ${\mathcal{X}}_{\Omega}^{3}$ with $d=2$ is the space of RGB-images with height and width $T_{1},T_{2}$. The value of an observation $X\in{\mathcal{X}}_{\Omega}^{P}$ at position $\omega\in\Omega$ is given by $X[\omega]=X[\omega_{1},\dots,\omega_{d}]\in{\mathbb{R}}^{P}$, and its restriction to the $p^{th}$ coordinate at each position is denoted $X_{p}\in{\mathcal{X}}_{\Omega}^{1}$. All signals are 0-padded, i.e. for $\omega\not\in\Omega$, $X[\omega]=0$. The convolution operator is denoted $*$ and is defined for $X\in{\mathcal{X}}_{\Omega}^{K}$ and $\boldsymbol{Y}\in{\mathcal{X}}_{\Omega}^{K\times P}$ as the signal $X*\boldsymbol{Y}\in{\mathcal{X}}_{\Omega}^{P}$ with

[TABLE]

For any signal $X\in{\mathcal{X}}_{\Omega}^{P}$, the reversed signal is defined as $X^{\Lsh}[\omega]=X[T_{1}-\omega_{1},\dots,T_{d}-\omega_{d}]^{\mathsf{T}}{}$ and the $p$-norm is defined as $\|X\|_{p}=\left(\sum_{\omega\in\Omega}\|X[\omega]\|_{p}^{p}\right)^{1/p}$. We will also denote ST the soft-thresholding operator defined as

[TABLE]

Convolutional Linear Model.

For an observation $X\in{\mathcal{X}}_{\Omega}^{P}$, the convolutional linear model reads

[TABLE]

with $\boldsymbol{D}\in{\mathcal{X}}_{\Theta}^{K\times P}$ a dictionary of $K$ atoms with the same number of dimensions $P$ on a domain $\Theta=\prod_{i=1}^{d}[0,L_{i}[\subset\Omega$ and $Z\in{\mathcal{X}}_{\Omega}^{K}$ the set of $K$ activation vectors associated with this dictionary. $\xi\in{\mathcal{X}}_{\Omega}^{P}$ is an additive noise term which is typically considered Gaussian and white. The dictionary elements $\boldsymbol{D}_{k}$ are patterns present in the data and typically have a much smaller support $\Theta$ than the full domain $\Omega$. The activations $Z_{k}$ encode the localization of the pattern $k$ in the data. As patterns typically occur at few locations in each observation, the coding signal $Z$ is considered to be sparse. Figure 1 illustrates the decomposition of a temporal signal with one pattern.

The parameters of this model are estimated by solving the following optimization problem, called convolutional dictionary learning (CDL)

[TABLE]

where the first term measures how well the model fits the data, the second term enforces a sparsity constraint on the activation, and $\lambda>0$ is a regularization parameter that balances the two. The constraint $\|D_{k}\|_{2}\leq 1$ is used to avoid the scale ambiguity in our model – due to the $\ell_{1}$ constraint. Indeed, the $\ell_{1}$-norm of $Z$ can be made arbitrarily small by rescaling $D$ by $\alpha\in{\mathbb{R}}_{+}^{*}$ and $Z$ by $\frac{1}{\alpha}$. The minimization (missing) 3 is not jointly convex in $(Z,\boldsymbol{D})$ but it is convex in each of its coordinates. It is usually solved with an alternated minimization algorithm, updating at each iteration one of the block of coordinate, $Z$ or $\boldsymbol{D}$.

Convolutional Sparse Coding (CSC).

Given a dictionary of patterns $\boldsymbol{D}$, CSC aims to retrieve the sparse decomposition $Z^{*}$ associated to the signal $X$. When minimizing only over $Z$, (missing) 3 becomes

[TABLE]

The problem in (missing) 4 boils down to a LASSO problem with a Toeplitz design matrix $\boldsymbol{D}$. Therefore, classical LASSO optimization techniques can easily be applied with the same convergence guarantees. Chalasani et al. (2013) adapted FISTA (Beck and Teboulle, 2009) exploiting the convolution in the gradient computation to save time. Based on the ADMM algorithm (Gabay and Mercier, 1976), Bristow et al. (2013) solved the CSC problem with a fast computation of the convolution in the Fourier domain. For more details on these techniques, see Wohlberg (2016) and references therein. Finally, Kavukcuoglu et al. (2010) proposed an efficient adaptation of the greedy coordinate descent (Osher and Li, 2009) to efficiently solve the CSC, which has been later refined by Moreau et al. (2018) with locally greedy coordinate selection (LGCD).

It is useful to note that for $\lambda\geq\|X*\boldsymbol{D}^{\Lsh}\|_{\infty}$, [math] minimizes the equation (missing) 4. The first order condition for (missing) 4 reads

[TABLE]

This condition is always verified if $\lambda>\lambda_{\max}=\|X*\boldsymbol{D}^{\Lsh}\|_{\infty}$. In the following, the regularization parameter $\lambda$ is set as a fraction of this value.

Dictionary Update.

Given a activation signal $Z^{(q)}$, the dictionary update aims at improving how the model reconstructs the signal $X$ using the fixed codes, by solving

[TABLE]

This problem is smooth and convex and can be solved using classical algorithms. The block coordinate descent for dictionary learning (Mairal et al., 2010) can be easily extended to the convolutional cases. Recently, Yellin et al. (2017) proposed to adapt the K-SVD dictionary update method (Aharon et al., 2006) to the convolutional setting. Classical convex optimization algorithms, such as projected gradient descent (PGD) and its accelerated version (APGD), can also be used for this step (Combettes and Bauschke, 2011). These last two algorithms are often referred to as ISTA and FISTA (Chalasani et al., 2013).

3 Locally Greedy Coordinate Descent for CSC

Coordinate descent (CD) is an algorithm which updates a single coefficient at each iteration. For (missing) 4, it is possible to compute efficiently, and in closed form, the optimal update of a given coordinate when all the others are kept fixed. Denoting $q$ the iteration number, and $Z^{\prime}_{k_{0}}[\omega_{0}]$ (atom $k_{0}$ at position $\omega_{0}$) the updated version of $Z^{(q)}_{k_{0}}[\omega_{0}]$, the update reads

[TABLE]

where $\beta^{(q)}$ is an auxiliary variable in ${\mathcal{X}}^{K}_{\Omega}$ defined for $(k,\omega)\in[1,K]\times\Omega$ as

[TABLE]

and where $e_{\omega}$ is a dirac (canonical vector) with value 1 in $\omega$ and 0 elsewhere. To make this update efficient it was noticed by Kavukcuoglu et al. (2010) that if the coordinate $Z_{k_{0}}[\omega_{0}]$ is updated with an additive update $\Delta Z_{k_{0}}[\omega_{0}]=Z^{\prime}_{k_{0}}[\omega_{0}]-Z^{(q)}_{k_{0}}[\omega_{0}]$, then it is possible to cheaply obtain $\beta^{(q+1)}$ from $\beta^{(q)}$ using the relation

[TABLE]

for all $(k,\omega)\neq(k_{0},\omega_{0})$. As the support of $(\boldsymbol{D}_{k_{0}}*\boldsymbol{D}_{k}^{\Lsh})[\omega]$ is $\prod_{i=1}^{d}[-L_{i}+1,L_{i}[$, equation (missing) 8 implies that, after an update in $\omega_{0}$, the value $\beta^{(q+1)}_{k}$ only changes in the neighborhood $\mathcal{V}(\omega_{0})$ of $\omega_{0}$, where

[TABLE]

We recall that $L_{i}$ is the size of an atom in the direction $i$. This means that only $\mathcal{O}(2^{d}K|\Theta|)$ operations are needed to maintain $\beta$ up-to-date after each $Z$ update. The procedure is run until $\max_{k,t}|\Delta Z_{k}[t]|$ becomes smaller than a specified tolerance parameter $\epsilon\geq 0$.

The selection of the updated coordinate $(k_{0},\omega_{0})$ can follow different strategies. Cyclic updates (Friedman et al., 2007) and random updates (Shalev-Shwartz and Tewari, 2009) are efficient strategies which only require to access one value of (missing) 7. They have a $\mathcal{O}(1)$ computational complexity. Osher and Li (2009) propose to select greedily the coordinate which is the farther from its optimal value. In this case, the coordinate is chosen as the one with the largest additive update $\max_{(k,\omega)}|\Delta Z_{k}[\omega]|$. The quantity $\Delta Z_{k}[\omega]$ acts as a proxy for the cost reduction obtained with this update.

This strategy is computationally more expensive, with a cost of $\mathcal{O}(K|\Omega|)$. Yet, it has a better convergence rate (Nutini et al., 2015; Karimireddy et al., 2018) as it updates in priority the important coordinates. To reduce the iteration complexity of greedy coordinate selection while still selecting more relevant coordinates, Moreau et al. (2018) proposed to select the coordinate in a locally greedy fashion. The domain $\Omega$ is partitioned in $M$ disjoint sub-domains $\mathcal{C}_{m}$: $\Omega=\cup_{m}\mathcal{C}_{m}$ and $\mathcal{C}_{m}\cap\mathcal{C}_{m^{\prime}}=\emptyset$ for $m\neq m^{\prime}$. Then, at each iteration $q$, the coordinate to update is selected greedily on the $m$-th sub-domain $\mathcal{C}_{m}$

[TABLE]

with $m=q\mod M$ (cyclic rule).

This selection differs from the greedy selection as it is only locally selecting the best coordinate to update. The iteration complexity of this strategy is therefore linear with the size of the sub-segments $\mathcal{O}(K|\mathcal{C}_{m}|)$. If the size of the sub-segments is $1$, this algorithm is equivalent to a cyclic coordinate descent, while with $M=1$ and $\mathcal{C}_{1}=\Omega$ it boils down to greedy updates. The selection of the sub-segment size is therefore a trade-off between the cyclic and greedy coordinate selection. By using sub-segments of size $2^{d}|\Theta|$, the computational complexity of selecting the coordinate to update and the complexity of the update of $\beta$ are matched. The global iteration complexity of LGCD becomes $\mathcal{O}(K|\Theta|)$. Algorithm 1 summarizes the computations for LGCD.

4 Distributed Convolutional Dictionary Learning (DiCoDiLe)

In this section, we introduce DiCoDiLe, a distributed algorithm for convolutional sparse coding. Its step are summarized in Algorithm 2. The CSC part of this algorithm will be referred to as DiCoDiLe-$Z$ and extend the DICOD algorithm to handle multidimensional data such as images, which were not covered by the original algorithm. The key to ensure the convergence with multidimensional data is the use of asynchronous soft-locks between neighboring workers. We also propose to employ a locally greedy CD strategy which boosts running time, as demonstrated in the experiments. The dictionary update relies on smart pre-computation which make the computation of the gradient of (missing) 6 independent of the size of $\Omega$.

4.1 Distributed sparse coding with LGCD

For convolutional sparse coding, the coordinates that are far enough – compared to the size of the dictionary – are only weakly dependent. It is thus natural to parallelize CSC resolution by splitting the data in continuous sub-domains. It is the idea of the DICOD algorithm proposed by Moreau et al. (2018) for one-dimensional signals.

Given $W$ workers, DICOD partitions the domain $\Omega$ in $W$ disjoint contiguous time intervals $\mathcal{S}_{w}$.

Then each worker $w$ runs asynchronously a greedy coordinate descent algorithm on $\mathcal{S}_{w}$. To ensure the convergence, Moreau et al. (2018) show that it is sufficient to notify the neighboring workers if a local update changes the value of $\beta$ outside of $\mathcal{S}_{w}$.

Let us consider a sub-domain, $\mathcal{S}_{w}=\prod_{i=1}^{d}[l_{i},u_{i}[$. The $\Theta$-border is

[TABLE]

In case that a $Z$ update affects $\beta$ in the $\Theta$-border $\mathcal{B}_{L}(\mathcal{S}_{w})$ of $\mathcal{S}_{w}$, then one needs to notify other workers $w^{\prime}$ whose domain overlap with the update of $\beta$ i.e.$\mathcal{V}(\omega_{0})\cap\mathcal{S}_{w^{\prime}}\neq\emptyset$. It is done by sending the triplet $(k_{0},\omega_{0},\Delta Z_{k_{0}}[\omega_{0}])$ to these workers, which can then update $\beta$ using formula (missing) 8. Figure 2 illustrates this communication process on a 2D example. There is few inter-processes communications in this distributed algorithm as it does not rely on centralized communication, and a worker only communicates with its neighbors. Moreover, as the communication only occurs when $\omega_{0}\in\mathcal{B}_{L}(\mathcal{S}_{w})$, a small number of iterations need to send messages if $|\mathcal{B}_{L}(\mathcal{S}_{w})|\ll|\mathcal{S}_{w}|$.

As a worker can be affected by its neighboring workers, the stopping criterion of CD cannot be applied independently in each worker. To reach a consensus, the convergence is considered to be reached once no worker can offer a change on a $Z_{k}[\omega]$ that is higher than $\epsilon$. Workers that reach this state locally are paused, waiting for incoming communication or for the global convergence to be reached.

DiCoDiLe-$Z$

This algorithm, described in Algorithm 3, is distributed over $W$ workers which update asynchronously the coordinates of $Z$. The domain $\Omega$ is once also partitioned with sub-domain $\mathcal{S}_{w}$ but unlike DICOD, the partitioning is not restricted to be chosen along one direction. Each worker $w\in[1,W]$ is in charge of updating the coordinates of one sub-domain $\mathcal{S}_{w}$.

The worker $w$ computes a sub-partition of its local domain $\mathcal{S}_{w}$ in $M$ disjoint sub-domains $\mathcal{C}_{m}^{(w)}$ of size $2^{d}|\Theta|$. Then, at each iteration $q$, the worker $w$ chooses an update candidate

[TABLE]

with $m=q\mod M$. The process to accept or reject the update candidate uses a soft-lock mechanism described in the following paragraphs. If the candidate is accepted, the value of the coordinate $Z_{k_{0}}[\omega_{0}]$ is updated to $Z^{\prime}_{k_{0}}[\omega_{0}]$, beta is updated with (missing) 8 and the neighboring workers $w^{\prime}$ are notified if $\mathcal{V}(\omega_{0})\cap\mathcal{S}_{w^{\prime}}\neq\emptyset$. If the candidate is rejected, the worker moves on to the next sub-domain $\mathcal{C}_{m+1}^{w}$ without updating a coordinate.

Interferences.

When $W$ coordinates $(k_{w},\omega_{w})_{w=1}^{W}$ of $Z$ are updated simultaneously by respectively $\Delta Z_{k_{w}}[\omega_{w}]$, the updates might not be independent. The local version of $\beta$ used for the update does not account for the other updates. The cost difference resulting from these updates is denoted $\Delta E$, and the cost reduction induced by only one update $w$ is denoted $\Delta E_{k_{w}}[\omega_{w}]$. Simple computations, detailed in Appendix A, show that

[TABLE]

If for all $w$, all other updates $w^{\prime}$ are such that if $\omega_{w^{\prime}}\notin\mathcal{V}(\omega_{w})$, then $(\boldsymbol{D}_{k_{w}}*\boldsymbol{D}_{k_{w^{\prime}}}^{\Lsh})[\omega_{w^{\prime}}-\omega_{w}]=0$ and the updates can be considered to be sequential as the interference term is zero. When $\left|\left\{\omega_{w}\right\}_{w}\cap\mathcal{V}(\omega_{0})\right|=I_{0}>1$, the interference term does not vanish. Moreau et al. (2018) show that if $I_{0}<3$, then, under mild assumption, the interference term can be controlled and DICOD converges. This is sufficient for 1D partitioning as this is always verified in this case. However, their analysis cannot be extended to $I_{0}\geq 3$, limiting the partitioning of $\Omega$.

Soft Locks.

To avoid this limitation, we propose the soft-lock, a novel mechanism to avoid interfering updates of higher order. The idea is to avoid updating concurrent coordinates simultaneously. A classical tool in computer science to address such concurrency issue is to rely on a synchronization primitive, called lock. This primitive can ensure that one worker enters a given block of code. In our case, it would be possible to use this primitive to only allow one worker to make an update on $\mathcal{B}_{L}(\mathcal{S}_{w})$. This would avoid any interfering update. The drawback of this method is however that it prevents the algorithm to be asynchronous, as typical lock implementations rely on centralized communications. Instead, we propose to use an implicit locking mechanism which keeps our algorithm asynchronous. This idea is the following.

Additionally to the sub-domain $\mathcal{S}_{w}$, each worker also maintains the value of $Z$ and $\beta$ on the $\Theta$-extension $\mathcal{E}_{L}(\mathcal{S}_{w})$ of its sub-domain, defined for $\mathcal{S}_{w}=\prod_{i=1}^{d}[l_{i},u_{i}[$ as

[TABLE]

These two quantities can be maintained asynchronously by sending the same notifications on an extended border $\mathcal{B}_{2L}(\mathcal{S}_{w})$ with twice the size of $\Theta$ in each direction. Then, at each iteration, the worker $w$ gets a candidate coordinate $(k_{0},\omega_{0})$ to update from (missing) 11. If $\omega_{0}\not\in\mathcal{B}_{L}(\mathcal{S}_{w})$, the coordinate is updated like in LGCD. When $\omega_{0}\in\mathcal{B}_{L}(\mathcal{S}_{w})$, the candidate is accepted if there is no better coordinate update in $\mathcal{V}(\omega_{0})\cap\mathcal{E}_{L}(\mathcal{S}_{w})$, i.e. if

[TABLE]

In case of equality, the preferred update is the one included in the sub-domain $\mathcal{S}_{w^{\prime}}$ with the minimal index. Using this mechanism effectively prevents any interfering update. If two workers $w<w^{\prime}$ have update candidates $\omega_{0}$ and $\omega^{\prime}_{0}$ such that $\omega^{\prime}_{0}\in\mathcal{V}(\omega_{0})$, then with (missing) 14, only one of the two updates will be accepted – either the largest one if $|\Delta Z_{k_{0}}[\omega_{0}]|\neq|\Delta Z_{k_{0}^{\prime}}[\omega^{\prime}_{0}]|$ or the one from $w$ if both updates have the same magnitude – and the other will be dropped. This way the $Z$ updates can be done in an asynchronous way.

Speed-up analysis of DiCoDiLe-$Z$

Our algorithm DiCoDiLe-$Z$ has a sub-linear speed-up with the number of workers $W$. As each worker runs LGCD locally, the computational complexity of a local iteration in a worker $w$ is $\mathcal{O}(2^{d}K|\Theta|)$. This complexity is independent of the number of workers used to run the algorithm. Thus, the speed-up obtained is equal to the number of updates that are performed in parallel. If we consider $W$ update candidates $(k_{w},\omega_{w})$ chosen by the workers, the number of updates that are performed corresponds to the number of updates that are not soft locked. When the updates are uniformly distributed on $\mathcal{S}_{w}$, the probability that an update candidate located in $\omega\in\mathcal{S}_{w}$ is not soft locked can be lower bounded by

[TABLE]

Computations are detailed in Appendix B. When the ratios $\frac{T_{i}}{W_{i}L_{i}}$ are large for all $i$, i.e. when the size of the workers sub-domains $|\mathcal{S}_{w}|$ are large compared to the dictionary size $|\Theta|$, then $P(\omega\not\in\textbf{SL})\simeq 1$. In this case, the expected number of accepted update candidates is close to $W$ and DiCoDiLe-$Z$ scales almost linearly. When $W_{i}$ reaches $\frac{T_{i}}{L_{i}}(\frac{2^{1/d}}{2^{1/d}-1})$, then $P(\omega\not\in\textbf{SL})\gtrsim\frac{1}{2}$, and the acceleration is reduced to $\frac{W}{2}$.

In comparison, Moreau et al. (2018) reports a super-linear speed-up for DICOD on 1D signals. This is due to the fact that they are using GCD locally in each worker, which has a very high iteration complexity when $W$ is small, as it scales linearly with the size of the considered sub domain. Thus, when sub-dividing the signal between more workers, the iteration complexity is decreasing and more updates are performed in parallel, which explains why the speed-up is almost quadratic for low $W$. As the complexity of iterative GCD are worst than LGCD, DiCoDiLe-$Z$ has better performance than DICOD in low $W$ regime. Furthermore, in the 1D setting, DICOD becomes similar to DiCoDiLe-$Z$ once the size of the workers sub-domain becomes too small to be further partitioned with sub-segments $\mathcal{C}_{m}^{(w)}$ of size $2^{d}|\Theta|$. The local iteration of DiCoDiLe-$Z$ are the same as the one performed by DICOD in this case. Thus, DICOD does not outperform DiCoDiLe-$Z$ in the large $W$ setting either.

4.2 Distributed dictionary updates

Dictionary updates need to minimize a quadratic objective under constraint (missing) 6. DiCoDiLe uses PGD to update its dictionary, with an Armijo backtracking line-search (Wright and Nocedal, 1999). Here, the bottle neck is that the computational complexity of the gradient for this objective is $\mathcal{O}(PK^{2}|\Omega|\log(|\Omega|))$. For very large signals or images, this cost is prohibitive. We propose to also use the distributed set of workers to scale the computations. The function to minimize in (missing) 6 can be factorized as

[TABLE]

with $\phi\in{\mathcal{X}}_{\Phi}^{K\times K}$ being the restriction of the convolution $Z^{\Lsh}*Z$ to $\Phi=\prod_{i=1}^{d}[-L_{i}+1,L_{i}[$, and $\psi\in{\mathcal{X}}_{\Theta}^{K\times P}$ the convolution $Z^{\Lsh}*X$ restricted to $\Theta$. These two quantities can be computed in parallel by the workers. For $\tau\in\Phi$,

[TABLE]

and for all $\tau,\omega\in\Phi\times\mathcal{S}_{w}$, we have $\omega-\tau\in\mathcal{S}_{w}\cup\mathcal{E}_{\Theta}(\mathcal{S}_{w})$. Thus, the computations of the local $\phi^{w}$ can be mapped to each worker, before being provided to the reduce operation $\phi=\sum_{w=1}^{W}\phi^{w}$. The same remark is valid for the computation of $Z^{\Lsh}*X$. By making use of these distributed computation, $\phi$ and $\psi$ can be computed with respectively $\mathcal{O}(K^{2}|\mathcal{S}_{w}|\log(|\mathcal{S}_{w}\cup\mathcal{E}_{\Theta}(\mathcal{S}_{w})|))$ and $\mathcal{O}(KP|\mathcal{S}_{w}|\log(|\mathcal{S}_{w}\cup\mathcal{E}_{\Theta}(\mathcal{S}_{w})|))$ operations. The values of $\psi$ and $\phi$ are sufficient statistics to compute $\nabla_{\boldsymbol{D}}F$ and $F$, which can thus be computed with complexity $\mathcal{O}(K^{2}P|\Theta|\log(2^{d}|\Theta|))$, independently of the signal size. This allows to efficiently update the dictionary atoms.

5 Numerical Experiments

The numerical experiments are run on a SLURM cluster with 30 nodes. Each node has 20 physical cores, 40 logical cores and 250 GB of RAM. DiCoDiLe is implemented in Python (Python Software Fundation, 2017) using mpi4py (Dalcín et al., 2005)222Code available at github.com/tommoral/dicodile.

5.1 Performance of distributed sparse coding

The distributed sparse coding step in DiCoDiLe is impacted by the coordinate descent strategy, as well as the data partitioning and soft-lock mechanism. The following experiments show their impact on DiCoDiLe-$Z$.

Coordinate Selection Strategy

The experimental data are generated following the sparse convolutional linear model (missing) 2 with $d=1$ in ${\mathbb{R}}^{P}$ with $P=7$. The dictionary is composed of $K=25$ atoms $\boldsymbol{D}_{k}$ of length $L=250$. Each atom is sampled from a standard Gaussian distribution and then normalized. The sparse code entries are drawn from a Bernoulli-Gaussian distribution, with Bernoulli parameter $\rho=0.007$, mean [math] and standard variation $10$. The noise term $\xi$ is sampled from a standard Gaussian distribution with variance 1. The length of the signals $X$ is denoted $T$. The regularization parameter was set to $0.1\lambda_{\max}$.

Using one worker, Figure 3 compares the average running times of three CD strategies to solve (missing) 4: Greedy (GCD), Randomized (RCD) and LGCD. LGCD consistently outperforms the two other strategies for both $T=150L$ and $T=750L$. For both LGCD and RCD, the iteration complexity is constant compared to the size of the signal, whereas the complexity of each GCD iteration scales linearly with T. This makes running time of GCD explode when $T$ grows. Besides, as LGCD also prioritizes the more important coordinates, it is more efficient than uniform RCD. This explains Figure 3 and justifies the use of LGCD is the following experiments.

Figure 4 investigates the role of CD strategies in the distributed setting. The average runtimes of DICOD (with GCD on each worker), and DiCoDiLe-$Z$ (with LGCD) are reported in Figure 4. We stick to a small scale problem in 1D to better highlight the limitation of the algorithms but scaling results for a larger and 2D problems are included in Figure C.1 and Figure C.2. DiCoDiLe-$Z$ with LGCD scales sub-linearly with the number of workers used, whereas DICOD scales almost quadratically. However, DiCoDiLe-$Z$ outperforms DICOD consistently, as the performance of GCD is poor for low number of workers using large sub-domains. The two coordinate selection strategies become equivalent when $W$ reaches $T/4L$, as in this case DiCoDiLe has only one sub-domain $\mathcal{C}^{(w)}_{1}$.

Impact of 2D grid partitioning on images

The performance of DiCoDiLe-$Z$ on images are evaluated using the standard colored image Mandrill at full resolution ($512\times 512$). A dictionary $\boldsymbol{D}$ consists of $K=25$ patches of size $16\times 16$ extracted from the original image. We used $\lambda=0.1\lambda_{\max}$ as a regularization parameter.

The reconstruction of the image with no soft-locks is shown in Figure 5. To ensure the algorithm finishes, workers are stopped either if they reach convergence or if $\|Z\|_{\infty}$ becomes larger than $\min_{k}\frac{50}{\|\boldsymbol{D}_{k}\|_{\infty}}$ locally. When a coefficient reaches this value, at least one pixel in the encoded patch has a value 50 larger than the maximal value for a pixel, so we can safely assume algorithm is diverging. The reconstructed image shows artifacts around the edges of some of the sub-domains $\mathcal{S}_{w}$, as the activation are being wrongly estimated in these locations. This visually demonstrates the need for controlling the interfering updates when more than two workers might interfere. When using DiCoDiLe-$Z$, the workers do not diverge anymore due to soft-locks.

To demonstrate the impact of the domain partitioning strategy on images, Figure 6 shows the average runtime for DiCoDiLe-$Z$ using either a line of worker or a grid of workers. In this experiment $K=5$ and the atom size is $8\times 8$. Both strategies scale similarly in the regime with low number of workers. But when $W$ reaches $T_{1}/3L_{1}$, the performances of DiCoDiLe-$Z$ stops improving when using the unidirectional partitioning of $\Omega$. Indeed in this case many update candidates are selected at the border of the sub-domains $\mathcal{S}_{w}$, therefore increasing the chance to reject an update. Moreover, the scaling of the linear split is limited to $W=T_{1}/2L_{1}=32$, whereas the 2D grid split can be used with $W$ up to $1024$ workers.

5.2 Application to real signals

Learning dictionary on Hubble Telescope images

We present here results using CDL to learn patterns from an image of the Great Observatories Origins Deep Survey (GOODS) South field, acquired by the Hubble Space Telescope (Giavalisco and the GOODS Teams, 2004). We used the STScI-H-2016-39 image333The image is at www.hubblesite.org/image/3920/news/58-hubble-ultra-deep-field with resolution $6000\times 3600$, and used DiCoDiLe to learn 25 atoms of size $32\times 32$ with regularization parameter $\lambda=0.1\lambda_{\max}$. The learned patterns are presented in Figure 7. This unsupervised algorithm is able to highlight structured patterns in the image, such as small stars. Patterns 1 and 7 are very fuzzy. This is probably due to the presence of very large object in the foreground. As this algorithm is not scale invariant, the objects larger than our patterns are encoded using these low frequency atoms.

For comparison purposes, we run the parallel Consensus ADMM algorithm444Code available at https://sporco.readthedocs.io/ proposed by Skau and Wohlberg (2018) on the same data with one node of our cluster but it failed to fit the computations in 250 GB of RAM. We restricted our problem to a smaller patch of size $512\times 512$ extracted randomly in the Hubble image. Figure C.3 in supplementary material reports the evolution of the cost as a function of the time for both methods using $W=36$. The other algorithm performs poorly in comparison to DiCoDiLe.

6 Conclusion

This work introduces DiCoDiLe, an asynchronous distributed algorithm to solve the CDL for very large multidimensional data such as signals and images. The number of workers that can be used to solve the CSC scales linearly with the data size, allowing to adapt the computational resources to the problem dimensions. Using smart distributed pre-computation, the complexity of the dictionary update step is independent of the data size. Experiments show that it has good scaling properties and demonstrate that DiCoDiLe is able to learn patterns in astronomical images for large size which are not handled by other state-of-the-art parallel algorithms.

Appendix A Computation for the cost updates

When a coordinate $Z_{k}[t]$ is updated to $Z^{\prime}_{k}[t]$, the cost update is a simple function of $Z_{k}[t]$ and $Z^{\prime}_{k}[t]$.

Proposition A.1**.**

If the coordinate $(k_{0},\omega_{0})$ is updated from value $Z_{k_{0}}[\omega_{0}]$ to value $Z^{\prime}_{k_{0}}[\omega_{0}]$, then the change in the cost function $\Delta E_{k_{0}}[\omega_{0}]$ is given by:

[TABLE]

Proof.

We denote $Z^{(1)}_{k}[\omega]=\begin{cases}Z^{\prime}_{k_{0}}[\omega_{0}],&\text{ if }(k,\omega)=(k_{0},\omega_{0})\\ Z_{k}[\omega],&\text{ elsewhere }\\ \end{cases}$. Let $\alpha_{0}$ be the residual when coordinate $(k_{0},\omega_{0})$ of $Z$ is set to 0. For $\omega\in\Omega$,

[TABLE]

We also have $\alpha_{k_{0}}[\omega]=(X-Z^{(1)}*\boldsymbol{D})[\omega]+\boldsymbol{D}_{k_{0}}[\omega-\omega_{0}]Z^{\prime}_{k_{0}}[\omega_{0}]$. Then

[TABLE]

This conclude our proof. ∎

Using this result, we can derive the optimal value $Z^{\prime}_{k_{0}}[\omega_{0}]$ to update the coordinate $(k_{0},\omega_{0})$ as the solution of the following optimization problem:

[TABLE]

In the case where multiple coordinates $(k_{w},\omega_{w})$ are updated in the same iteration to values $Z^{\prime}_{k_{w}}[\omega_{w}]$, we obtain the following cost variation.

Proposition A.2**.**

The update of the $W$ coordinates $(k_{w},\omega_{w})_{w=1}^{W}$ with additive update $\Delta Z_{k_{w}}[\omega_{w}]$ changes the cost by:

[TABLE]

Proof.

We define $Z_{k}^{(1)}[t]=\begin{cases}Z_{k_{w}}[\omega_{w}],&\text{ if }(k,\omega)=(k_{w},\omega_{w})\\ Z_{k}[t],&\text{ otherwise }\end{cases}~{}.$

Let

[TABLE]

and

[TABLE]

Then

[TABLE]

∎

Appendix B Probability of a coordinate to be soft-locked

Proposition B.1**.**

The probability of an update candidate being accepted, if the updates are spread uniformly on $\Omega$, is lower bounded by

[TABLE]

Proof.

Let $w\in[1,W]$ be a worker indice and $m\in[1,M]$ be the indice of the current sub-domain being considerd by the worker $w$. We denote $\omega=\operatorname*{\arg\!\max}_{\omega\in\mathcal{C}_{m}^{(w)}}\max_{k}|\Delta Z_{k}[\omega]|$ the position of the update candidate chosen by $w$ and $j=|\{w~{}\emph{ s.t. }~{}\mathcal{S}_{w}\cap\mathcal{V}(\omega)=\emptyset\}|$ the number of workers that are impacted by this update. With the definition of the soft-lock, there exists at least one worker in $\{w~{}\emph{ s.t. }~{}\mathcal{S}_{w}\cap\mathcal{V}(\omega)=\emptyset\}$ has a coordinate which is not soft-locked. As the updates are uniformly spread across the workers, we get the probability of $\omega$ being locked is $P(\omega\not\in\textbf{SL}~{}|~{}\omega)\geq\frac{1}{j}$.

Thus, the probability of an update candidate $(k,\omega)\in[1,K]\times\mathcal{S}_{w}$ to be accepted is

[TABLE]

where (missing) B.3 results from the summation on $\Omega$. This step is clearer by looking at Figure B.1.

∎

Appendix C Extra experiments

Scaling of DiCoDiLe-$Z$ for larger 1D signals

Figure C.1 scaling properties of DICOD and DiCoDiLe-$Z$ are consistent for larger signals. Indeed, in the case where $T=750$, DiCoDiLe-$Z$ scales linearly with the number of workers. We do not see the drop in performance, as in this experiments $W$ does not reach the scaling limit $T/3L$. Also, DiCoDiLe-$Z$ still outperforms DICOD for all regimes.

Scaling of DiCoDiLe-$Z$ for 2D signals

Figure C.2 illustrates the scaling performance of DiCoDiLe-$Z$ on images. The average running times of DiCoDiLe-$Z$ are reported as a function of the number of workers $W$ for different regularization parameters $\lambda$ and for greedy and locally greedy selection. As for the 1D CSC resolution, the locally greedy selection consistently outperforms the greedy selection. This is particularly the case when using low number of workers $W$. Once the size of the sub-domains becomes smaller, the performances of both coordinate selection become closer as both algorithms become equivalent. Also, the convergence of DiCoDiLe-$Z$ is faster for large $\lambda$. This is expected, as in this case, the solution $Z^{*}$ is sparser and less coordinates need to be updated.

C.1 Comparison with Skau and Wohlberg (2018)

Figure C.3 displays the evolution of the objective function (missing) 3 as a function of time for DiCoDiLe and the Consensus ADMM algorithm proposed by Skau and Wohlberg (2018). Both algorithms were run on a single node with $36$ workers. We used the Hubble Space Telescope GOODS South image as input data. As the algorithm by Skau and Wohlberg (2018) raised a memory error with the full image, we used a random-patch of size $512\times 512$ as input data. Due to the non-convexity of the problem, the algorithms converge to different local minima, even when initialized with the same dictionary. The Figure C.3 shows 5 runs of both algorithm, as well as the median curve, obtained by interpolation. We can see that DiCoDiLe outperforms the Consensus ADMM as it converges faster and to a better local minima. To ensure a fair comparison, the objective function is computed after projecting back the atoms of the dictionary $\boldsymbol{D}_{k}$ to the $\ell_{2}$-ball by scaling it with $\|\boldsymbol{D}_{k}\|_{2}^{2}$. This implies that we also needed to scale $Z$ by the inverse of the dictionary norm i.e.$\frac{1}{\|\boldsymbol{d}_{k}\|_{2}^{2}}$. This projection and scaling step is necessary for the solver as ADMM iterations do not guarantee that all iterates are feasible (i.e., satisfy the constraints). This also explains the presence of several bumps in the evolution curve of the objective function.