Learning-based design of random measurement matrix for compressed sensing with inter-column correlation using copula function

: In this work, a novel learning-based approach for the design of a compressed sensing measurement matrix is proposed. In contrast with the state-of-the-art methods, the suggested approach takes into account the correlation within entries of each column of the measurement matrix, namely, the inter-column correlation (ICC). The new method makes use of a rather small number of training sparse signal vectors in a recursive scheme to obtain their corresponding measurement vectors. The latter is exploited to estimate the copula function of measurements which, in turn, is used to generate an arbitrary number of measurement vector ensembles. By using the latter, the autocorrelation of the measurement vectors is estimated precisely and then, the ICC of measurement matrix under design is obtained from the autocorrelation. Given the resulting ICC, the measurement matrix columns are to be generated independently, e.g. by employing a proper random Gaussian vector generator. Performance evaluations using both synthetic and real-world data confirm the superiority of the proposed approach to the less recent methods.


Introduction
Over the past two decades, an extensive amount of research has been devoted to the problem of capturing a compact (compressed) set of measurements from a potentially sparse signal, i.e. one with many zero coefficients in some domain. Properly designed, such a measurement set of a sparse signal preserves its essential content, whereas, at the same time, it reduces the costs in storage, processing, and transmission of the signal. This problem, often termed as compressed sensing (CS), can be represented as follows: where y ∈ ℝ n denotes the real vector of observations, Φ ∈ ℝ n × p is a known (designed) measurement matrix with n < p, x ∈ ℝ p is a so-called m-sparse signal (i.e. with m non-zero entries at most) to be estimated, and v ∈ ℝ n represents the additive measurement noise [1]. Such a problem is raised in a wide range of applications including medical resonance imaging (MRI), spectroscopy, radar, Fourier optics, shortwave-infrared cameras, facial recognition, and network tomography [2,3]. Moreover, sampling and image superresolution reconstruction in the transform domain are important categories of information processing, which are associated with the topic of CS, e.g. [4,5]. Furthermore, rather recently, electrical impedance tomography, which calculates the internal conductivity distribution within a body using electrical contact measurements, has benefited from reconstruction methods based on CS [6].
Considering the problem statement in (1), the most primary question is that how many measurements suffice to perfectly recover some sparse signal, x, with high probability. The following corollary is an answer to this remarkable question [3].
Corollary 1: Let us draw each column of the measurement matrix, Φ, independently from the standard Gaussian distribution on ℝ p . If n ≥ Cmln(p/m), it is possible to reconstruct any m-sparse signal, x ∈ ℝ p , from the observations y in (1) with the probability exceeding 1 − e −Cn .
It is conventionally known that the optimal solution to the problem in (1) can be cast as the following minimisation problem min ∥ x ∥ 0 subject to y = Φx, (2) with ∥ . ∥ 0 being the ℓ 0 -norm, i.e. the number of non-zero entries in its argument. Yet, in practice, solving this ℓ 0 -minimisation problem is NP-hard and computationally infeasible for rather large signal dimensions [1]. As an alternative, the basis pursuit (BP) solution can be used [7], which replaces the ℓ 0 -norm by an ℓ 1norm, converting the CS reconstruction problem into a convex optimisation. There are also numerous variations of the latter method in the recent literature, e.g. [8,9]. Moreover, in [10], the authors aim at directly minimising the ℓ 0 -norm by a computationally fast method called smoothed ℓ 0 -norm. On the other hand, greedy algorithms with the orthogonal matching pursuit (OMP) being the most popular, have been developed, [11]. Moreover, in the statistical framework, Bayesian formulations for the estimation of sparse signals were suggested, wherein confidence measures for the estimated signal can be calculated [12]. A few major developments of the latter work have also been proposed, e.g. in [13], for the sparse component analysis in the presence of noise. To date, the most fundamental directions in the development of the CS theory include the following: the creation of efficient yet practical sparse signal reconstruction algorithms [14][15][16][17], learningbased adoption of suitable structures/models within the sparse signal [18][19][20], taking into account the suppression of the additive noise [21][22][23], and design of the CS measurement matrix under different practical circumstances. In this work, we basically focus on the latter topic, i.e. the design of the measurement matrix for CS in a generic sense. To the best of our knowledge, compared to the tremendous amount of the existing research on CS topic, fewer methods exist that focus on the design of suitable measurement matrices under general assumptions. Traditionally, most suggested CS algorithms tend to employ the Gaussian and Bernoulli measurement matrices with independent and identically distributed (i.i.d.) entries [1]. Although this choice of the sensing matrix provides robust and desirable performance in many applications of CS to reduce the reconstruction error and the computational load involved, a number of methods to design the matrix Φ have been suggested in the past. From a general point of view, the methods of CS matrix design can be categorised into deterministic and random measurement matrices. We present a brief extract of this literature in the following.
Whereas random measurement matrices ensure sparse recovery via ℓ 1 -minimisation with an optimal bound on the number of required measurements, deterministic matrices can be more suitable for applications subject to physical or other constraints and are less computationally expensive [1]. In [19], a framework to construct deterministic measurement matrices and to optimise them along with a sparsifying dictionary was proposed. Moreover, in [24], the concept of block p-ary codes is exploited to design deterministic matrices, which are able to perform almost as well as random matrices. More recently in [25], a subsampled structure of the form P Ω Ψ was considered for Φ with Ψ and P Ω , respectively, being the basis matrix and the subsampling operator. Therein, datadriven optimisation procedures were sought to obtain the index set, Ω, of the projector P Ω .
Measurement matrices with random entries have the advantage of being robust in terms of their performance, e.g. great reconstruction power, in most CS applications. Moreover, they universally satisfy the fundamental theoretical requirements of a generic CS matrix such as the restricted isometry property and the mutual incoherence [26]. Nevertheless, regarding the design of fully random sensing matrices, limited research studies exist within the literature. In chapter 9 of [1], a thorough review of the topic of sparse signal recovery with random matrices has been presented. Therein, generic random matrices whose entries are drawn independently according to a sub-Gaussian distribution have been studied with Gaussian and Rademacher distributions as special cases. Even though the incoherence property imposes that a universally suitable measurement matrix should have uncorrelated columns, there are still a few fundamental remarks that can delineate the employment of traditional random matrices (see Section III.A of [11] for a discussion of this). One such remark is the fact that each independent column of a random Φ may have some correlation within its own entries. This correlation, namely, the inter-column correlation (ICC), has been taken into account in [27] to design the independent columns of Φ. In the latter work, based on a variant of the structure for Φ suggested in [25], namely, P(ν)Φ with ν as an unknown parameter controlling the ICC, a number of training sparse signals are used to estimate ν. Next, using the resulting P(ν)Φ, an estimation for the ICC is performed which is, in turn, exploited in a Gaussian random vector generator to produce the columns of the ultimate random Φ. Looking into the prospective directions in CS, there are more futuristic concepts to come into use in measurement matrix design, such as the learningbased method aiming at maximising the entropy of measurement vectors in [28].
In this study, towards developing a more generic structure for a fully random measurement matrix, we take into account the aforementioned notion of ICC. In this regard, we make use of the concept of copula functions which has been successfully used in [29] to develop a CS framework with a multivariate Gaussian copula. The notion of ICC, apart from providing a more robust structure for the measurement matrix, performs as a degree-offreedom to adapt the measurement matrix to input observations, i.e. the sparse signal of interest. To this end, in contrast with [27], rather than limiting the proposed solution to a particular structure in the estimation of the ICC as well as only using the limited number of training signals, by making use of the concept of copula functions, we first generate an arbitrary number of measurement vectors, y, and then use these pseudo-measurements to estimate the present correlation across the elements of the vector y. This correlation is subsequently exploited to calculate the ICC by means of a straightforward yet accurate relationship. The i.i.d. columns of the measurement matrix under design are then generated using a random vector generator with the estimated ICC.
The remainder of this paper is organised as follows. In Section 2, the proposed algorithm for the design of Φ is outlined, which is based on the concept of copula from the statistics literature. In Section 3, the proposed approach to design the random measurement matrix is investigated and compared against the robust choices for the random Φ. Therein, the advantage of the proposed method is validated by means of objective error measures as well as subjective assessment. Section 4 is devoted to concluding remarks.

Proposed algorithm for measurement matrix design
In this section, we present in detail our approach for designing the random measurement matrix with ICC by leveraging the concept of copula. We present our methodology in the real field, ℝ, although it can be generalised to the complex field. We make no restricting assumptions on the structure/nature of the underlying sparse signal nor the measurement matrix in general.
Denoting the kth column of the measurement matrix , ϕ 2k , …, ϕ nk ] T for 1 ≤ k ≤ p, with T denoting the matrix transpose, we have for the elements of the ICC matrix of interest: r i j (k) ≜ E{ϕ ik ϕ jk } for 1 ≤ i, j ≤ n, or in the matrix form: R k = E{φ k φ k T } with E{ . } denoting the statistical expectation. Now, in accordance with [1,3], assuming i.i.d. and zero-mean columns for the measurement matrix under study, the subscript k can be dropped and we have that r i j (k) ≡ r i j or equivalently, R k ≡ R. Our goal within the rest of this section is to provide an appropriate approach for the estimation of the ICC to generate the random Φ. To this end, considering the problem statement in (1), by denoting Using (3) to calculate the autocorrelation of y gives where we used the fact that the additive noise, v i , is zero-mean and independent from the summation in (3) and also that Φ and x are independent. Assuming the independence of columns of Φ and having zero-mean elements at each column, it follows that E{ϕ ik ϕ jk′ } = E{ϕ ik }E{ϕ jk′ } = 0 unless k = k′. Applying the latter identity to (4) results in the following: Considering that E{ϕ ik ϕ jk } is independent of k and equal to the ICC term r i j , it is inferred that To obtain the ICC, i.e. r i j , using the above, the correlation of the additive noise has to be known or estimated in advance. Also, it is observed that the ICC of the measurement matrix scales with the autocorrelation of the measurements, E{y i y j }, whereas the latter is unknown. Directly estimating this term would require an enough number of training samples of the measurement y i , which is not available prior to designing the measurement matrix. To overcome this complication, we employ the concept of copula functions to properly generate a sufficient number of samples of the measurement vector, y, and then, use the generated samples in an averaging scheme to estimate the correlation E{y i y j } in (6). Algorithm 1 in the sequel outlines the suggested algorithm to determine the a priori measurements whereas Algorithm 2 details the overall framework proposed for the design of the measurement matrix. In the following subsections, we will elaborate on the different steps of these algorithms. Algorithm 1: Suggested approach for obtaining the a priori measurements Input: the training sparse signal samples, {x λ (t) } λ = 1 Λ , with λ and Λ denoting the sample index and the number of training signal vectors, respectively, the initial measurement matrix with i.i.d. Gaussian entries, Φ (0) , and the number of desired iterations, J. Output: the desired a priori measurements, {y λ (6) to obtain the ICC as • Employ a Gaussian random vector generator to generate the updated measurement matrix, Φ ( j) , with the ICC matrix as in Λ is the set of desired a priori measurements.

Determination of the a priori measurements
According to (1), without an assumption for the CS measurement matrix, there will be no prior knowledge of training samples for y obtained by employing the available training samples of the sparse signal of interest, x. To overcome this limitation, in an iterative fashion, we start by assuming an initial measurement matrix, Φ (0) , e.g. the Gaussian matrix with i.i.d. entries herein, and use it to sense the training sparse signal, {x λ (t) } λ = 1 Λ , to obtain the corresponding set of training, or namely, a priori measurements, at the first iteration. Next, we use the latter to provide an estimate of the autocorrelation of measurements as This correlation is then plugged into (6) to obtain a preliminary estimate of the ICC of the measurement matrix, Φ, which is exploited by a proper Gaussian random vector generator to generate the i.i.d. columns of the new matrix Φ. Next, again to obtain the updated measurement vectors, which are then used in the averaging scheme to estimate the updated autocorrelation, E{yy T }. The latter is again used in (6), resulting in the updated ICC that is used for generating the next measurement matrix, Φ, to sense {x λ (t) } λ = 1 Λ . At the end of a few iterations, the ultimate set of the measurement vectors, {y λ ( †) } λ = 1 Λ , is to be used in the estimation of the copula function of y, as will be discussed in the following subsection. • Use Algorithm 1 to obtain the a priori measurements, y λ † , from x λ (t) and Φ (0) , for λ = 1, 2, …, Λ. • Estimate the marginal cumulative distribution function's (CDF's) and probability distribution function's (PDF's) of the entries of y as well as the copula density function of y by using {y λ † } λ = 1 Λ , and adopt them to obtain the empirical joint PDF of y, as explained in Section 2.2.
• Exploit the empirical joint PDF of y to generate the pseudomeasurement vectors, {y ℓ ( ⋆ ) } ℓ = 1 L , with ℓ and L denoting the sample index and the number of generated pseudo-measurement samples, respectively, as discussed in Section 2.2.
in a smoothing scheme to estimate the ICC matrix, R = [r i, j ], as explained in Section 2.3.
• Employ a Gaussian random vector generator to independently generate the p columns of the measurement matrix, Φ, with the correlation matrix equal to R in the above.

Making use of copula to generate pseudo measurement vectors
As the next step of the proposed algorithm, it is required to generate a sufficient number of samples for the pseudomeasurement vectors, {y ℓ ( ⋆ ) } ℓ = 1 L , having the set of a priori measurements, {y λ ( †) } λ = 1 Λ . The most straightforward technique to fulfil this purpose is to estimate the parametric joint distribution of the measurement vector, y, using {y ℓ ( ⋆ ) } ℓ = 1 L . However, due to the inherent complexity in doing so, given the dimension of y in practice, we leverage the concept of copula and make use of the empirical copula to estimate the joint distribution. The latter is used to generate an arbitrary number of samples for y.
In the probability theory and statistics, a copula is defined as a multivariate probability distribution, c(Y 1 , Y 2 , …, Y n ), for which the marginal distribution of each random variable, Y i , is uniform [30]. According to Sklar's theorem, any multivariate joint distribution can be expressed in terms of univariate marginal distributions and a copula function, as the following [31]: where f (y) is the joint distribution of y, F i (y i ) and f i (y i ) are, respectively, the marginal CDF and PDF of y i , and c( . ) is the socalled copula density, which is the joint distribution of F i (y i )'s.
There are a few famous parametric copula families in the literature, which include parameters controlling the strength of dependence, e.g. the Archimedean copulas [30]. Although the estimation of copula parameters has been thoroughly studied in the literature, e.g. in [32], in our application, due to the high dimensionality of y, direct estimation of parametric copulas was found to be a tedious task. Thus, to avoid further complexity, the empirical (nonparametric) copula is preferred for the purpose of this work.
To empirically estimate the copula-based joint PDF in (8) } λ = 1 Λ , we shall have for the empirical marginal CDF's the following: IET Signal Process., 2020, Vol. 14 Iss. 6, pp. 385-395 © The Institution of Engineering and Technology 2020 for λ = 1, 2, …, Λ − 1, where y λ, i (S) is the ith entry of y λ (S) for i = 1, 2, …, n. Subsequently, for the empirical marginal PDF's, it follows that where we used the fact that f i (y i ) is the derivative (forward difference) of F i (y i ). To derive an estimate for the copula density function, c( . ), in (8), we resort to the empirical copula proposed in [33]. In this sense, by taking u i = F i (y i ) for i = 1, 2, …, n, and also defining the positive integer K as a precision parameter, with ζ = 1/K, the following expression can be deduced for the empirical copula density function [33] c^(u) = where u = [u 1 , u 2 , …, u n ] T , ). By this definition, N I K is in fact equal to the number of u λ vectors belonging to the set of intervals I K . In (11), it is observed that the copula density estimate at u located in a specific subcube, i.e. I K , is proportional to the number of measurements, u λ , located in I K divided over the entire number of measurements, Λ. This actually justifies the expression for c^(u) as the joint density of the vector u, which is in fact a piecewise constant function in each subcube of the unit n-cube.
At this stage, having estimated the marginal CDF's and PDF's as well as the copula density through (9)-(11), the joint PDF of measurements in (8) can be estimated. Next, the estimate of f (y) can be utilised in a convenient random vector generator to produce an arbitrary number of pseudo-measurement vectors, {y ℓ ( ⋆ ) } ℓ = 1 L . In this sense, we exploit the MATLAB tool provided in [34], namely, pwlCopula, wherein the generation algorithm for the pseudomeasurement vectors consists of two steps. First, the random vector u = [u 1 , u 2 , …, u n ] T is produced using the aforementioned estimated quantities, and then, each of its elements is transformed by the inverse of the corresponding marginal CDF, F i −1 (u i ). This MATLAB tool is moderately fast for implementation and it allows for the generation of random vectors with high dimensions as large as a few hundreds. A sufficient number, L, of the generated vector samples at the output of this tool will be used in the following subsection to estimate the ICC of the measurement matrix under design.

Smoothing scheme for the estimation of ICC
In this part, we use a straightforward smoothing scheme on the set of pseudo-random measurement vectors, {y ℓ ( ⋆ ) } ℓ = 1 L , to obtain the autocorrelation of measurements, E{y i y j }, and thus the ICC, r i j , given by (6). For doing so, given a sufficient number of samples, L, we approximate the expectation, E{ . }, by an averaging, resulting in the estimation, E{y i y j } = (1/L)∑ ℓ = 1 L y ℓ, i ( ⋆ ) y ℓ, j , with y ℓ, i ( ⋆ ) denoting the ith element of y ℓ ( ⋆ ) , 1 ≤ i, j ≤ n. In the case of complex field, however, to further reduce the estimation error, one can take advantage of the fact that the correlation matrix, E{yy H }, with H denoting the matrix Hermitian, is a Hermitian symmetric matrix, and therefore, it follows where y¯j denotes the complex conjugate of y j and E{y i y j } denotes the estimation of E{y i y j } by the aforementioned averaging scheme for the upper-diagonal elements, i.e. for j > i. Obviously, the above strategy is useful only for the off-diagonal elements of the correlation matrix. Finally, exploiting the estimated autocorrelation from (12) in (6) provides the desired ICC to be used as the correlation matrix in a proper Gaussian random vector generator to generate the columns of the measurement matrix, i.e. φ k for 1 ≤ k ≤ p.
It is worth mentioning that since each column of the proposed measurement matrix is generated independently, there will be no mutual coherence between each pair of columns. Therefore, the proposed ICC-based measurement matrix satisfies the mutual incoherence condition within the CS theory.

Experimental results
In this section, we present a summary of the experiments demonstrating the advantage of the proposed approach. To implement the suggested approach in Section 2, we took the number of training sparse signals, Λ, in the range of [100,200], the precision K = 20, and we subsequently generated L = 10 4 ensembles of the pseudo-measurement vectors, y ℓ ( ⋆ ) . Also, it was experimentally observed that having only J = 2 iterations suffices to obtain the set of a priori measurements in Algorithm 1, with higher J values providing no sensible improvement. In all conducted experiments, the set of training and test signals were chosen to be different but with similar signal parameters including d and m for the simulated signals. Finally, to independently generate the columns of the measurement matrix, Φ, we used the mvnrnd function of MATLAB with the parameters MU and SIGMAm respectively, set to [0] n × 1 and the ICC matrix, R, as obtained in (6).
To thoroughly investigate the advantage of the proposed measurement matrix, Φ, in Section 2 w.r.t. other measurement matrices, we evaluate their performance by using them in several major CS methods, namely, the OMP [11], Compressive sampling matching pursuit (CoSaMP) [35], Bayesian compressed sensing (BCS) [12], and BP [7] methods. We consider the traditional Φ with i.i.d. Gaussian elements, the deterministic measurement matrix suggested in [25], a Gaussian Φ with fixed ICC, as well as the suggested ICC-based Φ in [27]. The measurement matrix with deterministic elements in [25] is a subsampled structure of the form P Ω Ψ with Ψ and P Ω , respectively,being the basis matrix and subsampling operator. Therein, a data-driven optimisation procedure is sought to obtain the index set, Ω, of projector P Ω . In our implementation of the aforementioned method, we take the corresponding objective function to be f avg , which suggests an averaging over training signals prior to optimising for the index set Ω used in P Ω . As for the Φ with a fixed ICC, we use an exponentially decaying correlation structure, namely, e −ατ , with α and τ, respectively, being a fixed adjustable constant and the lag between each two entries of a column. Also, as with the suggested Φ in [27] with random elements, some structures are assumed for the ICC within the columns of Φ, which is controllable through a parameter that is optimised in a training-based framework.
For the evaluation of the methods under study, we consider both synthetic and real-world signals. In case of synthetic signals, we distinguish two different scenarios, one with a randomly generated binary signal and the other with a random zero-mean Gaussian signal. In the former case, we generate a binary signal of length p, consisting of m ones (sparsity degree) randomly distributed along with the signal array with the rest of its elements set to zero. This signal is to be recovered using n linear measurements as given by (1). Since perfect reconstruction is possible in this case, the evaluations can be performed in terms of the perfect reconstruction probability of the sparse signal.
We first investigate the performance of the proposed measurement matrix in Section 2 with different numbers of the training sparse signal vectors, Λ. Shown in Fig. 1 is the reconstruction error obtained by exploiting the proposed Φ in the OMP method versus Λ.
In this experiment, we considered reconstructing randomly generated Gaussian sparse signals contaminated by white noise in an signal-to-noise ratio (SNR) value of 10 dB; and we computed the normalised mean squared error (MSE) between the true and reconstructed signal. It is seen that by increasing the number of training signals from 20 to 100, the normalised MSE reduces dramatically but it remains almost fixed with higher Λ values for the parameter setting under study. This shows the number of required training samples for precisely estimating the underlying copula function, and therefore, the joint distribution of the measurement vector, which is considered in performing the foregoing experiments.
Next, we study the effect of the number of iterations used in Algorithm 1, namely, J, on the overall performance of the proposed measurement matrix in Section 2. To this end, using the previous setup, we obtained the reconstruction errors versus J, as shown in Fig. 2.
It is inferred that using only J = 2 iterations suffices for the aforementioned algorithm to provide the desired performance and that increasing J to higher values only results in trivial improvements. This fact led us to use no more than two iterations for all the other experiments and enabled us to avoid extra computational burden. Unless otherwise stated, we used the parameter settings J = 2 and Λ = 120.
Next, in Fig. 3, the reconstruction probability of a randomly generated binary signal is illustrated versus the sparsity level for the different measurement matrices. As viewed, independent of the four underlying CS approaches, the proposed Φ is capable of providing larger reconstruction probabilities for a reasonable range of the sparsity level, m. Furthermore, with increasing m, it is observed that the reconstruction probability obtained by the suggested Φ decays less sharply as compared to the other methods.
These demonstrate the robust advantage in taking into account the ICC in designing the measurement matrix with different reconstruction methods.
Moreover, in the case of non-binary (continuous) signals, we reconstruct randomly generated Gaussian samples in white noise at different SNR values. Next, we calculated the normalised MSE between the true and reconstructed signals versus SNR, as illustrated in Fig. 4. It is observed that using the proposed ICCbased measurement matrix helps further in reducing the reconstruction error w.r.t. the other methods. This advantage is more pronounced at higher SNR values, since none of the measurement matrices are basically designed to deal with adverse noisy conditions. Such a concept is to be considered as a topic of future research. It should be noted that, however, to fairly compare the proposed approach with the other methods in Fig. 3, we did not take into account the noise correlation matrix, E{vv T }, in implementing Algorithm 1; whereas integrating a proper estimate of E{vv T } into the latter, as performed before and shown in Fig. 4 (the plot named as proposed Φ with known noise), leads to superior performance.
Often in the relevant literature, the performance of CS methods against the number of observations is of high interest. For this reason, in Fig. 5, we have demonstrated the error performance of the proposed measurement matrix along with others with the number of measurements, n, changing from 60 to 140. As viewed in this experiment, the performance advantage of the proposed measurement matrix is apparent w.r.t. the less recent methods in consistence with n.
The majority of the superiority of the proposed Φ in Section 2 is due to making use of the learning-based process in Algorithm 2 to update the present ICC in accordance with the given training sparse signals, as against the Φ containing a fixed ICC. Moreover, as opposed to the deterministic learning-based Φ suggested in [25], the randomness of the proposed Φ allows for a better IET Signal Process., 2020, Vol. 14 Iss. 6, pp. 385-395 © The Institution of Engineering and Technology 2020 389 reconstruction accuracy. Pursuing the performance evaluation for the case of binary signals, we here investigate the number of required measurements, n, to identify a sparse signal with a fixed reconstruction probability. Fig. 6 displays the relationship between n and the sparsity level, m, to achieve a recovery probability of 95% with p = 250 using the OMP method with different measurement matrices presuming the noise-free condition. A trend of almost linearly growing n with increasing m is observed in this experiment. The advantage of using the proposed measurement matrix is evident by the smaller number of measurements it requires.
Subsequently, to investigate the relative performance of the proposed measurement matrix using a real-world data set, we implement the four reconstruction methods, namely, OMP, CoSaMP, BCS, and BP, using different measurement matrices on a dataset consisting of MRI images of the knee available in 'mridata' as discussed in [36]. For this purpose, we represent the MRI images in the discrete cosine transform (DCT) domain where the corresponding coefficient vectors of the images are sparse, and therefore, CS measurements can be realised efficiently. We also downsampled the original image files from 512 × 512 to smaller sizes to examine the results for different image resolutions. To avoid excessive complexity, we exploited the block-based CS, as

390
IET Signal Process., 2020, Vol. 14 Iss. 6, pp. 385-395 © The Institution of Engineering and Technology 2020 explained in Section 3.4 of [37]. In this strategy, an image is divided into B × B non-overlapping blocks and then acquired using the appropriately-sized measurement matrix. In this sense, we implemented the CS method on 16 × 16 image blocks with overlapping of 50% to reduce artefacts, resulting in the entire length of p = 256 for the sparse signal under reconstruction. In Fig. 7, the normalised MSE in the reconstruction of the MRI images versus the image resolution is shown for a fixed compression ratio, i.e. p/n = 4. It can be viewed that the proposed approach consistently provides less reconstruction errors for different image sizes. In a similar fashion, the normalised MSE versus the compression ratio for 256 × 256 images is demonstrated in Fig. 8. Under the constraint of p = 256, the compression ratio, in fact, controls the number of measurements, n. Therefore, this experiment also investigates the impact of the number of measurements on the performance. It is observed in Fig. 8 that by increasing the compression ratio, or equivalently, reducing the number of measurements taken, the reconstruction errors monotonically increase. Comparatively, the smaller MSE values observed consistently confirm the advantage of using the proposed measurement matrix design.
To conduct further analysis of the real-world data, we here take advantage of an image quality model, the peak SNR (PSNR) which is defined as [38] where b is the number of bits of each sample value of the images under test and MSE is that calculated as the mean of the squared errors between the corresponding pixels of the original and reconstructed image files. In Fig. 9, the PSNR values are shown versus the sampling rate, n/p, used in the underlying OMP method to recover the images in the DCT domain. We performed this experiment for three different 512 × 512 image samples. Numerical values of PSNR shown for the range of 0.1-0.8 of the sampling rate validate the superior quality of CS reconstruction by using the proposed measurement matrix. Next, to provide an estimate for the processing costs involved in implementing the measurement matrices, we employ them in the OMP method to reconstruct binary signals with specified parameters and measure the execution times in MATLAB. The corresponding results for the case of 1000 trials of signal reconstruction are shown in Fig. 10. Please note that this chart displays the execution time to complete 1000 trials, which includes the time to generate 1000 sparse signals in addition to the time required by 1000 invocations of the OMP algorithm. For the most computationally intensive method, each trial takes an average of 0.33 s. It should be noted that the smaller computational effort required by the method in [25] is mostly due to the deterministic nature of the suggested measurement matrix therein as opposed to the random structure in others. We believe that the relatively higher complexity imposed by the correlation-based measurement matrices is an acceptable compromise, considering their superior reconstruction performance.
Finally, to visually demonstrate the effectiveness of the proposed measurement matrix as compared with other methods, we show the recovered 512 × 512 MRI images of the knee from 'mridata' using the OMP algorithm with different measurement matrices under study in Figs. 11-16. For this experiment, we used the parameter setting n = 100 and m = 50 within the block-based implementation of the CS for the images. The advantage in the reconstruction by using the proposed measurement matrix compared with the other CS matrices is observable in this figure by means of the higher quality of the recovered images. It should be

392
IET Signal Process., 2020, Vol. 14 Iss. 6, pp. 385-395 © The Institution of Engineering and Technology 2020 mentioned that the set of training images was taken from the same database with the test and training images both being the figure of a knee MRI image captured at different shots.

Concluding remarks
A novel approach for the design of random CS measurement matrices, Φ, was proposed. The proposed Φ considers correlations within each column of Φ, namely, the ICC, in a generalised framework. The ICC is estimated using training sparse signals by means of a statistical tool, i.e. the copula function. The copula, employed in this work as an empirical estimate of the joint distribution of the measurement vector, is exploited to estimate the autocorrelation of the measurement vector, and subsequently, the ICC. We considered a random Gaussian measurement matrix as the designed Φ with the obtained ICC along each of its independent columns and compared its reconstruction performance with some of the state-of-the-art measurement matrices using different CS methods. Our results consistently confirmed the performance advantage of the suggested Φ. Moreover, by integrating a proper measurement noise estimator into the proposed approach, even further improvements can be achieved in the case of highly noisy measurements. Prospective research topics in this direction can include • Investigating application-oriented structures for the ICC in Φ as opposed to blindly estimating it. • Making use of the specific features of the sparse signal of interest in designing the measurement matrix. • Looking into the portion of correlation between the entries of the measurement vector due to the correlation present within the non-zero samples of the sparse signal. • A few practical extensions including the development and integration of a blind estimator for the measurement noise, using imperfectly-known training signals, and application-specific constraints on the measurement matrix.