Optimization of sparse recovery based on recursive unbiased predictive risk estimate

Sparse recovery via -penalized minimization can be typically solved by iteratively reweighted least squares (IRLS). The reconstruction quality is generally sensitive to the value of regularization parameter. In this wok, we propose a novel data-driven optimization scheme based on minimization of unbiased predictive risk estimate (UPRE). For a given IRLS iterate, we decompose and linearize the update by matrix splitting strategy, which makes the computation of Jacobian matrix tractable; and the UPRE of the IRLS iterate can be recursively evaluated. The developed recursive UPRE enables us to monitor the reconstruction error during the IRLS iterations, which can be applied to optimize regularization parameter. Numerical tests demonstrate the reliability of the recursive UPRE and the optimality of sparse reconstruction.


Introduction
Consider the following standard linear inverse problem: estimate the unknown coefficients from a linear model (1,2) : (1) where is the observed data, is a deterministic design matrix, is a vector of i.i.d.zero-mean Gaussian random variable with known variance .The classical approaches, e.g.ordinary least squares and ridge regression, generally have the drawbacks in prediction accuracy and interpretation (2) .In many applications, e.g.statistical inference (1,2) , signal recovery (3,4) and compressed sensing (5) , it is preferable to promote the sparsity of in some transform domain (e.g.wavelets), which is often formulated as the -penalized linear regression (2,3) : Here, denotes a decomposition that is believed to have a sparse representation of . is a regularization parameter that balances the data fidelity and the sparsity enhancement (1,3) .
There is a vast literature dedicated to solving P (3,4) .Iteratively reweighted least squares (IRLS) is an appealing choice due to its superior convergence speed (3,6) .We use IRLS to solve P in this paper.Note that we are not going to develop a more efficient algorithm to solve P, but will focus on the reconstruction quality, i.e., how to optimize the sparse reconstruction , solved by a given IRLS iterate.It has been recognized that -the solution to P-is generally sensitive to the value of regularization parameter (7,8) .Hence, many regularized iterative reconstruction algorithms, e.g.IRLS (6) and FISTA (4) , often require selection of the appropriate value for λ.A number of criteria for this selection have been proposed, e.g.discrepancy principle (9) , L-curve method (10) and generalized cross validation (GCV) (11) .Squared-error-based criterion can be a promising alternative to others: the prediction accuracy is often quantified in terms of expected prediction error (EPE for short) (1,12) : (2 Since EPE is inaccessible in practice due to unknown μ and x, the unbiased predictive risk estimate (UPRE) provides an unbiased substitute for EPE (1,12,13) : (3) which depends solely on , thus, can be accessed in practice.Refer to Ref.(12) for the proof of the unbiasedness of UPRE.Here, is a Jacobian matrix defined as: Any estimate can be regarded as a linear or non-linear transformation of the observed data .Here, we make the dependence explicit.In particular, if is a linear estimate: , then, .For the non-linear sparse estimate considered here, the computation of becomes essential for the UPRE evaluation.
The UPRE has been widely used in the statistics and signal processing communities, as a principled and efficient way for parameter selection with a variety of linear and non-linear estimators (7,8,12,13) .In this paper, for the non-linear sparse estimate by a given IRLS iterate, we develop a recursive procedure for computing UPRE, which enables us to evaluate the (estimated) prediction error for the sparse recovery.In particular, we incorporate matrix splitting strategy (14) into the IRLS algorithms, to make the computation of Jacobian matrix tractable.The optimal is then identified by exhaustive search for minimum UPRE.
Notations: To keep consistency throughout the paper, we use lowercase letter to denote vector, and uppercase for matrix.We denote the entry index of vector or matrix by the subscripts , and .The iteration is indexed by , and .

Recursive UPRE for IRLS algorithm 2.1 Basic scheme of IRLS
To solve P with fixed , the IRLS algorithm updates the solution by solving the following data-dependent equations (3,6) : (4) at -th iteration, where is a diagonal matrix and for .
Eq.( 4) leads to the following update of : (5) By Eq.( 3), UPRE of the -th iterate is: The UPRE computation requires to evaluate , which, however, cannot be directly evaluate from Eq.( 5), since is NOT a linear transformation of .

Matrix splitting to solve IRLS
Instead of the closed-form solution (5), we apply the matrix-splitting (MS) scheme (14) to solve Eq.( 4), which, as we will see later, enables us to evaluate in a recursive manner.This method splits the matrix as: where is solved by iteration indexed on j: (7) The MS iteration ( 7) is convergent to the unique solution of Eq.( 4) for any initial , if and only if the spectral radius (14)   .To guarantee the convergence, we choose with , such that (1) is positive definite; ( 2) is symmetric, positive definite, and thus, invertible; (3)   .Notice that it is difficult to compute .We now use matrix inversion lemma (MIL) to solve it: with .Thus, where is also solved by matrix splitting method: where , with the following iterate indexed on : (9) Similar to Eq.( 7), Eq.( 9) is also guaranteed to converge.

Recursion of Jacobian matrix
In this part, we will see that the computation of Jacobian matrix becomes tractable by the above IRLS-MS algorithm.Let us begin with the inner loop of iteration (9), which is rewritten here as , where , the super-and sub-scripts are ignored for brevity.The Jacobian matrix is derived as: (10)   for the (l,m)-th entry.The first term of Eq.( 10) is: where , by the property of Jacobian matrix.Noting that , if , the second term of (10) for fixed and is: Here, is diagonal with for .Finally, Eq.( 10) becomes: (11)   which expresses the recursion of Jacobian matrix.From ( 8), we obtain that other evolutions of Jacobian matrices are: (12) Thus, we can compute the UPRE during the IRLS iterations: (13)

Short summary
We summarize the proposed IRLS-MS-UPRE as Algorithm 1, which enables us to solve P with a prescribed value of λ, and simultaneously evaluate the UPRE during the IRLS iterations.
To find the optimal value of λ, an intuitive idea is to repeatedly implement Algorithm 1 for various tentative values of λ, then, the minimum UPRE indicates the optimal λ (see Fig. 2 for example).This global search, depicted in Algorithm 2, has been frequently used (1,7,8) .

Experimental results and discussions
In this section, we are going to solve P by IRLS, and present the results by a few numerical examples.
The convergence of the algorithm does not necessarily lead to the minimum prediction error: it is possible that the best reconstruction quality can be achieved before convergence.In principle, we could optimize the number of iteration.However, in all our experiments of this paper, we implement the algorithm until the final convergence is reached, such that the reconstruction is indeed a solution to problem.
We assume that (the variance of noise in y) was known in all simulations to compute the UPRE.In practice, σ can be estimated fairly reliably using, for example, MAD method (12) .

Random example
First, we consider a random example: we randomly generate the matrix , and set as a sparse vector with very few non-zeros (in this example, 10 non-zeros).Then, we add the noise ε with noise variance to obtain the observed data y = Ax + ε.Since x itself is already sparse in the original domain, for this case, we are solving a special case of P with decomposition matrix being identity (D=I).
We apply Algorithm 1 to solve P with fixed λ = 1.Fig. 1- (1) shows that the objective value keeps decreasing to converge; Fig. 1-(2) shows the evolutions of UPRE and true EPE during the iterations.We can see that the UPRE is always a reliable substitute for EPE.

Deconvolution of 1-D signal
We now consider a typical 1-D signal x.A represents a convolution matrix, which is circulant under periodic boundary condition, constructed from a Gaussian kernel.We use decimated wavelet transform (DWT) and choose Haar basis as the analysis matrix D. Fig. 4 shows the relation between UPRE and λ.

Conclusions
UPRE has been proven a powerful tool to select regularization parameter (7,8) .In this paper, to solve -minimization problem, we proposed a recursive UPRE for IRLS algorithms, by incorporating matrix splitting scheme.It enables us to monitor the (estimated) prediction loss during the IRLS iterations, without referring to the true unknown data.
Theoretical derivations in this work related to the UPRE evaluation can be extended, in principle, to other types of regularizers and regularized iterative reconstruction algorithms.We would also like to emphasize that not limited to the simple example shown here, the developed recursive UPRE can be applied to many practical applications, e.g.compressed sensing (5) and image deconvolution (4,12) .

Figure 2
Figure 2 The relation between UPRE and λ.

Fig. 3
Fig.3shows a part of the reconstructed x, which is close to the true vector.

Figure 3 A
Figure 3 A part of the reconstruction of x.

Figure 4
Figure 4 The relation between UPRE and λ.

Fig. 5
Fig.5shows a part of the reconstructed x, which is close to the true vector.

Figure 5 A
Figure 5 A part of the reconstruction of x.