Iterative Cauchy Thresholding: Regularisation with a heavy-tailed prior
Perla Mayo Robin Holmes Alin Achim
eess.IV



Evolution of photospheric flows under an erupting filament in the quiet-Sun region
    Jiří Wollmann1
          Michal Švanda1,2Corresponding author, svanda@sirrah.troja.mff.cuni.cz
          David Korda1
          Thierry Roudier3
          
    Received –; accepted –
===============================================================================================================================================================



In the machine learning era, sparsity continues to attract significant interest due to the benefits it provides to learning models. 
Algorithms aiming to optimise the ℓ_0- and ℓ_1-norm are the common choices to achieve sparsity. In this work, an alternative algorithm is proposed, which is derived based on the assumption of a Cauchy distribution characterising the coefficients in sparse domains. The Cauchy distribution is known to be able to capture heavy-tails in the data, which are linked to sparse processes. We begin by deriving the Cauchy proximal operator and subsequently propose an algorithm for optimising a cost function which includes a Cauchy penalty term. We have coined our contribution as Iterative Cauchy Thresholding (ICT). Results indicate that sparser solutions can be achieved using ICT in conjunction with a fixed over-complete discrete cosine transform dictionary under a sparse coding methodology.



iterative Cauchy thresholding, ISTA, IHT, sparsity, proximal operator




§ INTRODUCTION


Representation learning (RL) involves the acquisition of features that enable the reconstruction or recovery of a signal of interest. Such features can be learnt in several ways. In dictionary learning it is assumed that an observed signal y∈IR^M can be approximately represented or reconstructed by a linear combination of features in a matrix A∈IR^MxN with a vector of coefficients x∈IR^N. This is expressed as:


    ŷ = Ax


where ŷ∈IR^M is the estimation of the observed signal y. The error between this approximation and the true signal is thus given by ℒ(A, x) = || y - ŷ||_2^2. In Eq. (<ref>), each coefficient x_i provides the contribution of a column vector (feature) a_i in the reconstruction task for i=1,..., N. Additionally, along with learning of the vectors in A, it is of interest to understand their contribution. The goal is then to recover the coefficients within x, which is done through the optimisation of 


    x^*    = xargmin      ℒ(A, x)  
    		             = xargmin       ||y - Ax||_2^2


When the matrix A is over-complete (N >> M), the problem in Eq. (<ref>) becomes ill-posed with an infinite number of solutions. To alleviate this, a regularisation term is added in Eq. (<ref>) to enforce sparsity among the elements of the coefficients x_i:


    x^*    = xargmin      ℒ(A, x) + λφ (x)


where λ is a trade-off parameter between data fidelity and the penalty term.

In Eq. (<ref>), the ℓ_0-norm is the desired choice for φ(·). This norm is ideal since it effectively counts the number of non-zero entries in x, therefore, the problem is to find the optimal set of coefficients x_i that will allow the reconstruction of y by using as few as possible column vectors from A. Since this is a combinatorial (NP-hard) problem, its optimisation becomes challenging. An alternative to this is to replace the ℓ_0- by its relaxed version, the ℓ_1-norm. Furthermore, different penalty functions can be obtained if a different approach is followed. When the coefficients are assumed to follow a heavy-tailed distribution centred at 0 it is possible to derive alternative penalties. 
If a Laplace distribution is assumed for the coefficients, then the penalty function corresponds to the ℓ_1-norm <cit.>. 
The use of heavy tailed distributions has also been considered, e.g. in <cit.>, where it was assumed that coefficients follow the Symmetric α-Stable distribution. 

Our aim is to prove the power of the Cauchy distribution in modelling data with heavy tails, and consequently design novel sparse coding approaches.

The rest of the paper is structured as follows: in section <ref> a brief summary of related work is presented along with the motivation behind it, section <ref> describes the proposed algorithm in detail, section <ref> presents the experiments carried out along with the obtained results. Lastly, a conclusion is presented in section <ref> along with a discussion of future work.



§ RELATED WORK


The optimisation of either the ℓ_0- or ℓ_1-norms can be done through iterative thresholding approaches such as, e.g., the Iterative Hard Thresholding (IHT) <cit.> and the Iterative Soft Thresholding (IST) <cit.> algorithms. They are expressed as shown in equations (<ref>) and (<ref>) respectively.


    x = 
        
             x,    if  |x| > x_0 
    
             0,    if  |x| ≤ x_0



    x = 
        
             x - x_0,    if  x > x_0 
    
             x + x_0,    if  x < -x_0 
    
             0,    otherwise


These approaches achieve sparsity by snapping to zero any value below a given threshold. In the case of IST the values are shrunk by x_0. In addition, the function in (<ref>) can be derived by computing the proximal operator of the penalty function, obtained by solving the expression in Eq. (<ref>) below.


    prox_λφ(x)    = zargmin       (z-x)^2 + λφ (z)


 For IST, using φ(·) = |·| results in the expression in (<ref>) with x_0 = λ/2. For more details on proximal algorithms, please refer to <cit.>.



 §.§ Cauchy distribution as prior for sparsity


As previously said, the penalty function φ(·) can be defined from a maximum a posteriori (MAP) approach by assuming the coefficients in x follow a heavy tailed distribution centred at 0. The Laplacian distribution is implicitly chosen when the regularisation term is the ℓ_1-norm. The Cauchy distribution has been found to be an interesting alternative <cit.>. It belongs to the family of the α-Stable distributions and it is one of the few that has a closed form expression for its probability density function (p.d.f.), which is defined in Eq. (<ref>).


    p(x) = γ/π(γ^2 + (x-δ)^2)


The Cauchy distribution is described by two parameters for location (δ) and scale (γ). The behaviour of these parameters is illustrated in Fig. <ref>. It can be seen that a smaller value for γ corresponds to sparser densities and a larger concentration around the location  δ. When δ=0, most values will concentrate around the origin. In addition, the smaller γ, the more aggressive the sparsity enforcing property of the corresponding algorithm will be.





§ A NEW SHRINKAGE ALGORITHM


The shrinkage algorithm proposed in this work relies on the proximal splitting approach for a Cauchy penalty function. This results in the optimisation of (<ref>) one function at a time. Furthermore, φ(·) is solved using the shrinkage operator obtained from the computing of the proximal operator for the Cauchy penalty function. This is further reviewed in section <ref> whilst the full algorithm is described in detail in section <ref>.



 §.§ The Cauchy Proximal Operator


The Cauchy proximal operator is obtained when Eq. (<ref>) is solved using φ(x) = -log p(x) (from Eq. (<ref>)) and assuming δ = 0. Therefore, the proximal operator for this penalty function is computed from:


    prox_λφ(x)    = argmin    (z - x)^2 + λφ(z) 
       = argmin    (z - x)^2 - λlog(γ/π(γ^2 + z^2 ) )


up to a constant, which results in:


    z^3 - x z^2 + (γ^2+ λ) z - γ^2 x = 0


Using Cardano's formula to solve this cubic expression, we obtain:


    z = x/3 + t


where

    t    = √(-q/2+√(Δ))
    Δ   = q^2/4+p^3/27
    
            p    = λ + γ^2-x^2/3
    
            q    = -2/27x^3 + 1/3(λ-2γ^2)x


This will provide only one real root when Δ > 0, and three real roots otherwise. In the latter case the root with the largest absolute value is chosen. 

The shrinkage behaviour offered by the Cauchy proximal operator does not require one to explicitly specify a threshold value, as must be done for IST and IHT. However, the shape of the Cauchy proximal operator still depends on the parameters selection. Such shape is an indicator of how aggressive the cut-off will be. Figure <ref>(a) shows the impact of the parameter λ in the Cauchy proximal operator, whilst Figure <ref>(b) does the same for γ. It can be seen that λ determines the location of the threshold and γ its shape. A more aggressive cut-off will be exhibited when γ→ 0.

In fact, when γ→ 0, its value will be so small that it will barely contribute anything to the function to minimise. If γ = 0, the function to minimise becomes:




    prox_λφ(x) 
                = argmin    (z - x)^2 + λlog(z^2)


which is now a non-smooth non-convex function. Its solutions, however, can be computed following the proximal operator approach once again:


    z = 
        x/2 + √(x^2-4λ)/2   , x ≥ 2λ
    x/2 - √(x^2-4λ)/2   , x ≤ -2λ
    
             0,    otherwise


The shape of ICT under these conditions is shown in figure <ref>b): as γ→ 0, the threshold → 2 (λ = 1). On the other hand, when λ = 0, the solutions are given by γ i, -γ i, and x, which implies that there is no shrinkage. This comes from the fact that λ=0 removes the penalty on the coefficients from x, and so the learning relies purely on the reconstruction error.








 §.§ Iterative Cauchy Thresholding


Once the proximal operator has been obtained, it is our goal to solve Eq. (<ref>) for the coefficients using the Cauchy penalty function. Explicitly, the cost function is now:


    f(x) = ||y-Ax||_2^2 - ∑_ilog(γ/π(γ^2 + x_i^2))





This is solved by means of a proximal splitting approach. Under this framework, the function f(x)=ℒ(x) + λφ(x) can be optimised by separating ℒ(x) and λφ(x), therefore, the solutions for each function are projected into the next one. For this, only the gradient of ℒ(x) is computed to solve ℒ(x). Such solution is then passed to the solver for the Cauchy penalty function, which is the shrinkage function given by Eq. (<ref>). The shrinkage is applied in an entry-wise fashion. The process is iterated until a stopping criterion is met, that is, after a fixed number of iterations or when the error has reached a certain predefined value. This is the Iterative Cauchy Thresholding (ICT) algorithm, and its pseudo-code is shown in Algorithm <ref>.







§ EXPERIMENTAL RESULTS


The performance of ICT has been assessed for a reconstruction task via Sparse Coding (SC). The problem to solve for SC is the same as the one posed in Eq. (<ref>). The performance of ICT is compared against the ones for IST and IHT, where only the penalty function changes during the optimisation accordingly. The dictionary used is the fixed over-complete discrete cosine transform (DCT) using a patch size of 8x8 and a total of 144 atoms. Due to the nature of SC, the image to reconstruct is thus split in T overlapping patches of the same size as the atoms (8x8). There is one vector of sparse coefficients x for every one of these patches. The reconstructed image corresponds to the average of the reconstructed overlapping patches. Further details on sparse coding are beyond the scope of this work and can be found elsewhere <cit.>.

The images to reconstruct correspond to the AT&T faces dataset [www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.html], the MNIST[http://yann.lecun.com/exdb/mnist/] dataset as well as the Lena and the Shepp-Logan phantom images. The thresholding algorithms were evaluated with varying values for their parameters. To obtain different levels of sparsity the parameter λ was modified for the three algorithms. 
The learning rate η was fixed for the three algorithms to 0.005. We randomly chose 500 observations for the MNIST dataset and 50 for the AT&T faces. 
Tables <ref> and <ref> report the average results on this reconstruction task for the three thresholding algorithms, PSNR values along with the % of non-zero elements (ℓ_0-norm) are reported, whilst Figure <ref> shows the performance as function of the iterations for the four datasets.



The performance of ICT is very similar to the one produced by IHT, both achieving better PSNRs than the  reconstruction results obtained with IST by a considerable margin, for most levels of sparsity considered. The best PSNR values obtained by the three algorithms correspond to the lowest degree of sparsity for each and vice-versa. It is the latter case the one of interest. 

The similarity in behaviour between IHT and ICT could be attributed to the shape of their functions. In Fig. <ref> it can be seen that the ICT shrinkage function lies between those of both IHT and IST and it affects all coefficients regardless of their amplitude.
It is important to note as well that the reconstruction itself corresponds to the average of all estimations obtained for every overlapping patch extracted from the original observations. Hence there will always be a subtle blurring effect affecting the final image. This can be alleviated by employing a different framework, such as convolutional sparse coding <cit.>.



§ CONCLUSIONS AND FUTURE WORK


 In this paper we introduced a new algorithm to recover sparse coefficients in regularised inverse problems. The proposed approach is based on the assumption that sparse coefficients can be accurately  modelled by Cauchy distributions. Consequently, we introduce the Cauchy proximal operator and describe a corresponding proximal splitting approach, which we prefer to refer to as ICT. As an illustrative example, we chose to address the problem of image reconstruction from learnt sparse representations, and we achieve a performance which compares very favourably with the IHT algorithm. We conjecture that the proposed ICT algorithm has the potential to perform better in other inverse problems, including, e.g. image denoising or deconvolution, due to the optimal penalty applied to large coefficients. This is the focus of our current research endeavours.
 






IEEEbib


