Adaptive Algorithm in Image Reconstruction Based on Information Geometry

Abstract

Compressed sensing in image reconstruction has attracted attention and many studies are proposed. As we know, adding prior knowledge about the distribution of the support on the original signal to CS can improve the quality of reconstruction. However, it is still difficult for a recovery framework adjusts its strategy for exploiting the prior knowledge efficiently according to the current estimated signals in serial iterations. With the theory of information geometry, we propose an adaptive strategy based on the current estimated signal in each iteration of the recovery. We also improve the performance of existing algorithms through the adaptive strategy for exploiting the prior knowledge according to the current estimated signal. Simulations are presented to validate the results. In the end, we also show the application of the model in the image.

1. Introduction

Compressed Sensing (CS) is a signal processing technique to recover sparse signals of interest from far fewer samples and has been widely used in image reconstruction due to the booming of media [1]. Chai et al. proposed an image encryption algorithm based on the memristive chaotic system, elementary cellular automata and CS [2]. Zha et al. aproposed a group-based sparse representation method with non-convex regularization for image CS reconstruction [3]. Furthermore, CS was applied to reconstruction of under-sampled atomic images [4]. In 2018, the scholars proposed an image encryption method integrating CS and detour cylindrical diffraction [5].

The CS is to acquire and reconstruct a signal x∈R^N by finding a solution to underdetermined linear system [6-8]

\(y=\Phi x+e\)       (1)

where y∈R^M is the compressed representation of x, Φ=(Φ_ij)_M×N is the random projection matrix, e∈R^M represents the noise, and M<<N. x is inherently sparse or sparse in some domain. That is to say, most of its entries are zero, or there exists an orthogonal transformational matrix ψ satisfying most of the entries of ψ^-1x are zero. One may recover an estimated x as the solution to a l0-minimization problem P0 [9], which is discontinuous and NP-hard. Considering that under the restricted isometry property [10], the P0 is usually relaxed to :

\(\text { P1 } \quad \hat{x}=\underset{x \in R^{N}}{\arg \min }\|x\|_{l_{1}} \quad \text { s.t. } y=\Phi x\)       (2)

To solve the two problems (P0 and P1), many achievements have been obtained. Recovery frameworks fall into three categories: greedy algorithms [11], convex optimization algorithms [12-14], and statistical algorithms [15-17].

Because of lower computational complexity, greedy algorithms are used widely in CS. A typical greedy algorithm is Orthogonal Matching Pursuit (OMP) [18]. OMP can accommodate various scenarios. And it can exactly recover a signal with S nonzero entries in dimension N by O(S ln N) random linear measurements. So OMP has been studied by many scholars and a family of variants are supported, such as Stagewise Orthogonal Matching Pursuit (StOMP) [19], Regularized Orthogonal Matching Pursuit (ROMP) [20], Compressed Sampling Matching Pursuit (CoSaMP) [21], Iterative Hard Threshold (IHT) [22], Gradient Descent with Sparsification (GraDeS) [23] and Generalized Orthogonal Matching Pursuit (GOMP) [24], and so on.

A well-known convex optimization algorithm is the Least Absolute Shrinkage and Selection Operator (LASSO) [25]

\(\begin{aligned} &\tilde{x}=\underset{x \in R^{N}}{\arg \min } \frac{1}{2}\|y-\Phi x\|_{l_{2}}^{2}\\ &\text { s.t. } \quad\|x\|_{l_{1}} \leq t \end{aligned}\)       (3)

where t≥0 is a tuning parameter. A stage-wise fast LASSO algorithm for the image reconstruction from CS optimizes an insensitive Huber objective function to achieve a decision function [26]. Lian et al. [27] exploit a multilevel prior support information model and incorporate it into the LASSO using a weighted l1-norm penalty function. Least Angle Regression (LARS) provides an efficient algorithm for computing the solution paths of LASSO [28]. LARS is an algorithm for fitting linear regression models to high-dimensional data [29]. In 2007, Keerthi et al. [30] applied a fast tracking algorithm for LARS to sparse kernel logistic regression.

In many applications, it may be possible to catch some prior knowledge about the vector representation of a natural image in the domain of some linear transform [31,32]. For example, it can be known that some regions are more likely to be highlight than the others in magnetic resonance images [33]. In statistical algorithms, adding the prior knowledge to CS can improve the quality of reconstruction. There are many studies on sparse recovery with prior knowledge, which enhance the signal recovery by statistic in various degrees. Next, we introduce previous typical statistical algorithms. Babacan et al. [34] proposed an excellent greedy algorithm using Laplace priors based on Bayesian framework. Besides that, the main idea of model-CS is to introduce an extra structure to exploit the distribution of the original signal [35]. In [36], a method based on approximate message passing and Markov-tree prior was supported. Furthermore, neural networks were trained to learn the prior knowledge for recover signal [37]. However, the computational complexity is high because of the extra structure as a model. For the tree based algorithms, the reconstruction involves the travel across the whole tree in every iteration. As far as we know, the average complexity per iteration for the fastest tree-based algorithm for CS is O(MN²)[38]. In addition, the algorithm based on neural networks involves the forwarding from the bottom layer to the top layer of the networks. And the lower bound on average complexity per iteration for the algorithm is O(2N) for a single forwarding computation with a simplest architecture of a hidden layer of 2[(α-2)/M+1] hidden neurons, where  is the size of the training set and   N [39]. Instead of the models, a modified method works by exerting punishment on the selection for a wrong support [40]. In [41], relations between entries of the original signal were built by Kalman filter based on the assumption that support changes slowly. But the Kalman filter does not always use the actual co-variances of the changes of the estimated signals. Furthermore, weighted minimization provided better upper bounds on the reconstruction error where the weights were determined by prior knowledge about a specific probability of nonzero entries [42,43]. In 2018, the weights were chosen to minimize the expected statistical dimension of the descent cones of a weighted cross-polytope [44]. Moreover, such optimal weights have been used to improve greedy algorithms and convex optimization algorithms [45], such as the Weighted Orthogonal Matching Pursuit (WOMP) and Weighted LARS (W-LARS). The algorithms combine either greedy algorithms or convex optimization algorithms with statistic by weights. In the algorithms, a new estimated signal is generated in an iteration according to the result from its previous iterations. However, the optimal weights are fixed in all iterations in the whole procedure of the reconstruction without the consideration of the changing of the estimated signal. Thus, it is more reasonable to adjust the strategy for exploiting the prior knowledge according to the current estimated signal.

Information geometry [46] is a mathematical branch that applies the techniques of differential geometry to study probability theory and statistic. It has been widely used to many applications, such as asymptotic theory of statistical inference [47], the expectation maximisation (EM) algorithm [48], the learning of neural networks [49-51], and many others [52-56]. Wang et al. proposed a model-CS method with neural networks based on information geometry [57]. However, the computational complexity is high because of the neural networks as an extra structure.

To address the problems of the previous algorithms mentioned above, we introduce information geometry to improve the performance of the greedy algorithms and convex optimization algorithms with no extra structure. It is possible to reflect the image features with a specific manifold structure [58]. Based on the fact, we evaluate the candidates of the estimated signals with the Fisher distance [59] between the distribution over estimated support set and the known distribution over the true support on the original signal. Then, we choose the recovery signal which has a distribution close to the true one according to the evaluations. The strategy can be seen as adaptive to the developed support sets generated in different iterations, because we calculate the estimated distribution according to the current estimated support set in each iteration. Experiments are presented to validate the results.

2. Methods

2.1 Construction of Geometry Model

As discussed above, the true distribution of the indexes of the non-zero entries of x is able to be known by statistic. The index distribution can be described by a specific Probability Density Function (PDF) p(k|θ) where θ=(θ₁,θ₂,…,θ_N)'.

As prior knowledge, we know that the true parameters for the original signals are θ¹. In addition, θ² are estimated parameters in the model. In the iterative algorithm, we estimate the distribution by computing θ² over the set of optional indexes in each iteration. We choose the recovery signal in the sense that its support has a distribution closer to the known distribution. The details are described as follows.

A distribution from the exponential family [60] is defined as a function in the following form

\(p(k \mid \theta)=\exp \{\theta \cdot k-\varphi(\theta)\}\)       (4)

where ψ is a function of θ.

As a discrete distribution, the position distribution is from the exponential family naturally. Suppose that the frequency of the occurrence of index i∈{0,1,…,N} is p(k=i|θ)=p_i, where θ=(θ₁,θ₂,…,θ_N)'. And we denote the i th element of θ as θ₁. Then, it can be rewritten as

\(p(\xi \mid \theta)=\exp \{\theta \cdot \xi-\varphi(\theta)\}\)       (5)

where

\(\xi_{i}=\delta_{i}(k)=\left\{\begin{array}{l} 1, k=i \\ 0, k \neq i \end{array}\right.\)       (6)

\(\theta_{i}=\log \frac{p_{i}}{p_{0}}\)       (7)

then,

\(p_{0}=1-\sum_{i=1}^{N} p_{i}=1-p_{0} \sum_{i=1}^{N} e^{\theta_{i}}\)       (8)

From (8), we obtain

\(\varphi(\theta)=\log \left(1+\sum_{i=1}^{N} e^{\theta_{i}}\right)\)       (9)

A classic parametric space for this family of PDF’s is

\(H=\left\{\theta \in R^{N} \mid \theta \geq 0\right\}\)       (10)

The Fisher Information Matrix (FIM) is defined as

\(G(\theta)=\left[g_{i j}(\theta)\right]\)       (11)

\(g_{i j}(\theta)=E\left[\frac{\partial \log p(k \mid \theta)}{\partial \theta_{i}} \cdot \frac{\partial \log p(k \mid \theta)}{\partial \theta_{j}}\right]\)       (12)

where E means the expectation. The FIM is a way of measuring the amount of information of the parameter. And a distance called Fisher distance arises from the FIM [61]. The distance between two points θ¹=θ(t1) and θ²=θ(t2) in the half-plane H measures the dissimilarity between the associated distributions p(x|θ¹) and p(x|θ²). Furthermore, the Fisher distance is given by the value which is the minimum of the lengths of all the piece-wise smooth paths θ(t) that joins θ¹ and θ², t₁≤t≤t₂. That is 

\(D_{F}\left(\theta^{1}, \theta^{2}\right):=\min _{\left\{\theta(t): \theta\left(t_{1}\right)=\theta^{1}, \theta\left(t_{2}\right)=\theta^{2}\right\}} \int_{t_{1}}^{t_{2}} \sqrt{\left(\frac{d \theta}{d t}\right)^{\prime} G(\theta) \frac{d \theta}{d t}} d t\)       (13)

where t is the parameter of the curve θ(t). The above equations can be transformed to the Euler-Lagrange equations as

\(\frac{d^{2} \theta_{h}}{d t^{2}}+\sum_{q=1}^{n} \sum_{z=1}^{n}\left[\frac{1}{2} \sum_{l=1}^{n} g^{h l}\left(\frac{\partial g_{q l}}{\partial \theta_{z}}+\frac{\partial g_{z l}}{\partial \theta_{q}}-\frac{\partial g_{q z}}{\partial \theta_{l}}\right)\right] \frac{d \theta_{q}}{d t} \frac{d \theta_{z}}{d t}=0, \forall h, l, q, z \in\{1,2, \ldots, N\}\)       (14)

where the subscripts have the same meanings as above, while the superscripts index the elements in the inverse matrix of G(θ) in (11).

2.2 Construction of Geometry Model

How to choose support is important to the recovery algorithms. The common method is described as follows. Firstly, the correlation between each column of measurement matrix and residual is calculated. Secondly, the index (or indexes) which has (have) the highest correlation is (are) added to the support. So the correlation is the indicator of determining support.

In this paper, we improve the traditional indicator by a reward term. If the addition of an index into the support leads to a decrease of D_F between the estimated distribution and the known distribution, the reward will be set higher. Define the reward term by:

\(\operatorname{reward}(\tilde{A})=\frac{1}{\lambda\left(\theta^{1}, S\right)+D_{F}\left(\theta^{1}, \theta^{2}\right)}\)       (15)

where λ(θ¹, S) is a constant which offsets D_F’s influence on the indicator. It depends on the parameters and the sparsity. The details will be discussed in section 3.

θ² is determined with the support \(\widetilde{A}\) generated in each iteration. In an iteration, the support set \(\tilde{A}=\left\{k^{1}, k^{2}, \ldots, k^{j}\right\}\) has j elements. We obtain θ² by the Maximum Likelihood Estimate (MLE) [62].

\(\theta^{2}=\underset{\theta}{\arg \max } \sum_{i=1}^{j} p\left(k^{i} \mid \theta\right)\)       (16)

The reward term is used to modify the evaluation of the support in each iteration. An indicator is multiplied by its reward term. It is adaptive according to the different estimated support set in each iteration. In another word, it is an adaptive strategy based on the current estimated signal in each iteration of the recovery. The improved algorithm for greedy algorithm is shown in Algorithm 1. The improved algorithm for convex optimization algorithm is shown in Algorithm 2.

Algorithm 1

Algorithm 2

fun() in (*) is defined by

\(\operatorname{fun}\left(r^{j-1}, \phi_{i}\right)=\left\{\begin{array}{lr} \left\langle W r^{j-1}, \phi_{i}>,\right. & j \leq \text { thres }_{\mathrm{min}} \\ <W r^{j-1}, \phi_{i}>\cdot \operatorname{reward}\left(A^{j-1} \cup k^{i}\right), j>\text { thres }_{\mathrm{min}} \end{array}\right.\)       (17)

where <> represents inner product in Euclidean space, W is a positive definite matrix and thres_min is a threshold. For non-weighted algorithms, W is a unit matrix and for algorithms with weights, W is a diagonal.

In the algorithms discussed above, the computing power mainly focuses on the computation of (15) to (17). In the j th iteration, the computation described above will be implemented for N-j+1 candidates. However, the computation doesn’t involve any travel or propagation across a model with a complex architecture. Thus, we infer that the computational complexity of the proposed algorithm is lower than that of the model-based algorithms introduced in section 1.

3. Experimental Results and Discussion

Generally, the analysis of various applications allows the extraction of prior knowledge about the specific distribution over the supports of the signal’s sparse representation. In this section, we take a series of experiments with the simulated signals and imaging signals.

3.1 Simulation

We choose two PDF randomly in Table 1 as the known index distributions to generate 1000 simulated signals following each of them, respectively. We suppose that the latitude of each nonzero entry follows a standard normal distribution. And N is set to 300. We set iter_max=8×M. The res_min is 10^-5, while the thres_min is 30 (it is big enough for static). λ is determined based on different sparsity for each PDF respectively, as shown in Table 1.

Table 1. Settings for PDF

^aFor simplicity, λ is set to the expectation of DF over the samples subtracted by 1.

The signal is acquired from noisy measurements with SNR_mes=30dB. Accurate recovery is declared when

\(\|x-\tilde{x}\|_{l_{2}} /\|x\|_{l_{2}} \leq 10^{-1}\)       (18)

During the experiment, M/N is increasing gradually. Meanwhile, we compute the accurate recovery rate. Firstly, the proposed method is used to improve OMP and WOMP, respectively. The improved algorithms are called IG-OMP and IG-WOMP, respectively. We illustrate the results in Fig. 1 and Fig. 2. Secondly, proposed method is used to improve LARS and WLARS, respectively. The improved algorithms are called IG-LARS and IG-WLARS, respectively. We illustrate the results in Fig. 3 and Fig. 4.

Fig. 1. Comparation among OMP, IG-OMP, WOMP and IG-WOMP with the distribution described by f₁

Fig. 2. Comparation among OMP, IG-OMP, WOMP and IG-WOMP with the distribution described by f₂

Fig. 3. Comparation among LARS, IG-LARS, WLARS and IG-WLARS with the distribution described by f₁

Fig. 4. Comparation among LARS, IG-LARS, WLARS and IG-WLARS with the distribution described by f₂

To validate the superiority of the proposed algorithm in the computational complexity, we compare the performance of an algorithm based on information geometry with a typical model-based algorithm for CS (i.e., neural networks for CS introduced in section 1) in both of the average running time per iteration and the accurate recovery rate. The experimental settings are as the same as above. For each algorithm, we take the average over all conditions. Naturally, the improved algorithms are more complex than the original algorithms. In addition, OMP has a lower time complexity than LARS. Similarly, IG-WOMP has a lower time complexity than IG-WLARS. Hence, IG-WLARS is the most complex one among the algorithms (OMP, IG- OMP, WOMP, IG-WOMP, LARS, IG-LARS, WLARS and IG-WLARS). So the running time of the lower complex algorithms are meaningless in the validation. Then, we choose the result of the most complex one (i.e., IG-WLARS) for the comparison in Table 2. For simplicity, Neural Networks for CS is abbreviated as NNCS.

Table 2. Results

3.2 Application Examples

Imaging is a typical signal whose support follows a specific distribution. Here, x is an imaging signal in DCT domain. The image is blocked, so that N is 32 ×32. And the M is set to 200.

In Table 3 and Table 4, the results are shown numerically. Moreover, the results are shown intuitively in Fig. 5 and Fig. 6.

Table 3. Results

^aThe definition of the Peak-Signal-to-Noise Ratio (PSNR) is given in [63].

Table 4. Results

Fig. 5. (a) Original image (b) IG-WOMP (c) WOMP (d) IG-OMP (e) OMP

Fig. 6. (a) Original image (b) IG-WLARS (c) WLARS (d) IG-LARS (e) LARS

To validate the superiority of the proposed algorithm in the computational complexity, we compare the performance of the most complex algorithm based on information geometry in the paper (i.e., IG-WLARS) with a typical model based algorithm for CS (i.e., NNCS) in both of the average running time and the recovery qualities of the images (i.e., the average PSNR over all images). And we list the results in Table 5.

Table 5. Results

3.3 Discussion

Firstly, we discuss Fig. 1, Fig. 2, Fig. 3, and Fig. 4. S-sparse means x is with S nonzero entries. Under different sparsities, algorithms improved by information geometry outperform original algorithms constantly. The results of two types of decoders are shown: original decoders (OMP, WOMP, LARS, and WLARS) and improved decoders (IG-OMP, IG-WOMP, IG-LARS, and IG-WLARS). Original decoders have worse performance.

Secondly, we discuss the results in Fig. 5, Fig. 6, Table 3, and Table 4. As it is shown in Fig. 5, Fig. 6, Table 3, and Table 4, the recovery is applied to different standard test images (fruits, boat, camera man). In Table 3 and Table 4, a higher PSNR indicates that the reconstruction is of higher quality. In Fig. 5 and Fig. 6, the reconstructions of images show that algorithms improved by the proposed method outperform the original algorithms in real data.

Thirdly, we discuss the results in Table 2 and Table 5. In Table 2, the running time of IGWLARS is lower than NNCS, while the values of accurate recovery rate are close. In Table 5, the running time of IG-WLARS is lower than NNCS, while the recovery quality of the IGWLARS is better. The results of running time of these algorithms contrast sharply with each other, while the values of the recovery quality are very close. The results support the conclusion that the proposed algorithm has a lower computational complexity, while keeping a commensurate recovery accuracy.

4. Conclusion

With the theory of information geometry, we propose an adaptive strategy based on the current estimated signal in each iteration of the recovery. And we improve the performance of recovery algorithms through the adaptive strategy for exploiting the prior knowledge about the index distribution. Simulations are presented to validate the results. In the end, we also show the application of the model in the image.