For the estimation of constants in the standard restricted isometry condition for a tight frame, the techniques adapted from the closely related area of Sparse Principal Component Analysis are applied. Based on the paper of M.Journee, Yu.Nesterov, P.Richtarik, and R.Sepulchre (2010), we consider the related parametrized l1-penalty optimization problem and propose an algorithm for the solution of corresponding nonlinear eigenvalue problems arising from the optimality conditions. A generalization of the method for the case of complex-valued data is developed. The efficiency of the proposed method is verified by numerical results obtained for several important test examples. The computed bounds are compared with the exact values obtained by direct search, and with the lower bounds expressed in terms of extreme roots of the properly parametrized Jacobi polynomials.
The construction and analysis of the so-called low-coherence redundant bases is a hot topic in data retrieval and processing during past few decades. Presented in the form of rectangular matrices, such objects are used, for instance, in compressed sensing technologies, see, e.g., recent book [Foucart & Rauhut, 2013] and survey [Casazza & Kutyniok, 2013]. Considering a complex-valued rectangular matrix A ∈ Cm×n, where m < n, it is required that any m × k submatrix
AJ = [aj1 | : : : |ajk ]; 1 ≤ j1 < : : : < jk ≤ n; of A has full column rank. Here k ≤ m is as large as possible and J = {j1; : : : ; jk} is an arbitrary ordered subset of pairwise distinct indices. Moreover, certain practical considerations lead to the requirement of wellconditioning for any such submatrix AJ . The latter assumption is often characterized by the presence of the restricted isometry property (hereafter RIP) of the kth order, requiring that the inequalities 1 − 1(k; A) ≤ ∥Az∥2=∥z∥2 ≤ 1 + 2(k; A) hold with some 0 < 1(k; A) < 1 and 2(k; A) > 0 for all sparse n-vectors z ∈ Cn such that ∥z∥0 ≤ k (hereafter, ∥z∥0 denotes the number of nonzero components of the vector z). We will consider the optimum (unimprovable) (1) (2) bounds in (2), which are defined via minimum and maximum eigenvalues of k × k matrices AJH AJ , where k ≤ m, as follows: Here and further on, ∥z∥2 = zH z is the squared Euclidean norm of z, and BH is the conjugate transpose of B.
Note that the problem under consideration can be readily reformulated in terms of the Gram matrix 1 − 1(k; A) = min |J|=k as the problem of finding maximum sparse eigenpairs of a Hermitian nonnegative definite matrices G and I − G (with properly chosen positive constants and , see Section 3.2 below). Hence, as it was noticed in [d’Aspremont et al., 2008], the evaluation of the RIP constants is equivalent to the solution of a pair of sparse Principal Component Analysis (Sparse PCA) problems. A survey of sparse PCA applications and an efficient algorithm for approximate solution of the problem (in the case of real-valued matrix G) can be found in [Journee et al., 2010]. Recall that in [Tillmann & Pfetsch, 2014], the NP-hardness of the exact evaluation of RIP constants was shown.
In this paper, we will restrict the exposition by the case of “equal norm tight frames”, that is, we assume that Diag(AH A) = In and
AAH = n m
Im (such properties typically correlate with better RIP constants, and therefore are advantageous from the application viewpoint). We will present a reformulation of the Generalized Power Method of [Journee et al., 2010] applicable to the analysis of complex-valued matrices A (whereas its original version can only be applied to realvalued data). Moreover, the proposed algorithm is targeted to approximation of the functions of k = 1; 2; : : : ; m defined by (3) and (4) rather than to solving a typical Sparse PCA problem (such as extracting a single sparse dominant principal component of the data matrix A). 2
According to the above discussion, the original sparse PCA problem is set (with account of (5)) as zk = arg max ∥z∥0=k zH z
; zH Gz or zk = arg max ∥Az∥ : ∥z∥0=k ∥z∥ Next we consider an L1-relaxation of this problem and its reformulation. 2.1
Parametrized L1-relaxation of Sparse PCA Problem Following [Journee et al., 2010] and recalling that ∥z∥1 = ∑jn=1 |zj |, let us replace the original problem by z( ) = arg max (√zH Gz − ∥z∥1); ∥z∥=1 or z( ) = arg ∥mz∥a=x1 (∥Az∥ − ∥z∥1) ; where is a positive parameter which controls the sparsity degree k = ∥z∥0. As follows from the discussion in Section 2.1 of [Journee et al., 2010], if A is a column normalized matrix (that is, (G)j;j = 1 for all j by (5) and (6)), then it suffices to consider the variation of the parameter in the range 0 < < 1. Note that larger values of typically correspond to smaller values of k. 2.2
Reformulation of the Optimization Problem Using the representation ∥v∥ = max∥x∥=1 |xH v| with v = Az, one has, similar to [Journee et al., 2010] and [d’Aspremont et al., 2008],
( ∥z∥=1 ∥mx∥a=x1 |xH Az| − ∥z∥1 ∥mz∥a=x1 (∥Az∥ − ∥z∥1) = max )
∥x∥=1 ∥mz∥a=x1 (|xH Az| − ∥z∥1) = max (3) (4) (5) (6) n = max max ∑ |zj | (|(AH x)j | − ∥x∥=1 ∥z∥=1 j=1 n ∑ (|(AH x)j | − )= ∥mx∥a=x1 j=1 =|ajH x|)+ : Here aj is the jth column of A, so that (AH x)j = ajH x, and is a positive scalar.
Recall that our purpose is to find the spectral bounds (3) and (4) as functions of the discrete argument k. Therefore, one needs a reasonable strategy to reduce the number of different problems (7) to solve. 3
Below we propose a construction of rather small discrete set of values for the penalty parameter still allowing to estimate the required spectral bounds for all k = 2; : : : ; m.
In accordance with the hardness result [Tillmann & Pfetsch, 2014], the iterative approximation to the solution of equation (7) may not, in general, deliver the global optimum. Therefore, the use of many iteration restarts is necessary to obtain a satisfactory result. On the other hand, a high-performance implementation of such method on modern supercomputers is straitforward due to its inherently parallel structure. 3.1
The method is based on the variation of the penalty parameter using the formula (l) = 1 + ( 2 − 1)l=lmax; l = 0; : : : ; lmax; where 0 < 1 < 2 < 1. The variation of quasirandom initial guess for the vector x is made by the use of the so-called logistic sequence with starting term equal to l: i = 1 − 2 i2−1; i = 1; : : : ; n; 0 = l; and the sth initial guess is formed as yi = i + i i−1; x(0) = y=∥y∥ (here and elsewhere, we denote i = √−1).
For each l we perform a sufficient number of power iterations (exactly the same as in Algorithm 2 in [Journee et al., 2010] except of specialization (8) for g(·)): As the main result of iterations, one can consider the index set x(t+1) = g(x(t))=∥g(x(t))∥;
t = 0; 1; : : : ; tmax: J (t) = {j1; : : : ; jk} iff |ajHs x(t)| > ; s = 1; : : : ; k: Using J (t) corresponding to a well-converged value of (t) = ∥g(x(t))∥, one then finds the maximum eigenvalue of the matrix AJH AJ . The latter is taken as an approximation from below to the quantity max|J|=k max(AJH AJ ). (7) (8) (9) In order to reduce the problem size, we will use the standard Naimark complement techniques, see, e.g., [Casazza et al., 2013]. Recall that the Naimark complement of an m × n tight frame A is the (n − m) × n tight frame B satisfying, in particular, the identity mAH A + (n − m)BH B = nIn. Thus, we first evaluate the Gram matrix of B: where U is an upper trapezoidal (n − m) × n-matrix and P is a permutation matrix. Finally, we apply the algorithm of the above Subsection to the matrix B (of course, with replacement of m by n − m in (??)), and the required result readily follows from the obvious relation holding for any index subset J , which yields
AJH AJ = n m Ik − n − m m
BJH BJ ; min |J|=k min(AJH AJ ) = n m − n − m m max max(BJH BJ ): |J|=k In fact, this means that one should use the index subsets J arising in the generalized power algorithm applied to B in order to form matrices AJH AJ and to determine their minimum eigenvalues. 4
In absence of fast algorithms for determining the exact values of RIP constants, it is useful to have at least “inner” bounds for these quantities. Such bounds (in general, not tight for k ≥ 4) were found in [Kaporin, 2017].
Theorem. Let A be an m × n tight frame, where m < n. Then, in the RIP condition of the form ∥z∥0 ≤ k →
0 < 1 − 1(k) ≤ ∥Az∥2=∥z∥2 ≤ 1 + 2(k); the best possible RIP constants must satisfy 1(k) = 1 − |mJ|=ink 1(k) ≥ 1 − 2nm (1 − max(k; m; n)); 2(k) ≥ 2nm (1 − min(k; m; n))− 1; where min(k; m; n) and max(k; m; n) are the smallest and the largest roots of the Jacobi polynomial P (m−k;n−m−k)(x), respectively.
k
In other words, for a prescribed size parameters k; m; n the quantities 1 and 2 will be separated from zero below as indicated above for any choice of elements in A. 5
We will consider several tight frames of the sizes n ≈ m2, as well as another cases with n ≈ 2m, all of them having rather good conditioning. Recall that the coherence between the columns of normalized frame A is measured by the parameter (A) = max |akH al| ≥ k̸=l √ n − m ; (n − 1)m which is always positive if n > m. The smaller this value, the better is the performance of compressed sensing algorithms, at least within the limits predicted by related theory. The inequality presents the Welch bound [Welch, 1974]. The equality in the Welch bound is attained for certain m and n if there exists the corresponding optimum equiangular tight frame (further on ETF). (10) (11) 5.1
The partial Discrete Fourier Transform (further on DFT) matrices are both well explored theoretically and important in applications. Generally, these are determined by
1 (A)i;j = √m exp that is, A is formed by picking out of the n × n DFT matrix a subset of its m pairwise different rows. Moreover, it is known that such matrix is an optimum equiangular tight frame if n = m2 − m + 1 and the sequence {d(j)|j = 0; : : : ; m − 1} forms the so called difference set, see, e.g., [Strohmer & Heath, 2003].
We will consider the following equiangular tight frames based on difference sets. Using in (12) the Singer difference set for m = 17 and n = 273 given by see [Hall, 1956, Xia et al., 2005, LJ2017], we derive complex-valued test matrix denoted as SINGER17x273. Recall that it possess the ETF property with (A) = √m − 1=m.
Another test matrix is defined by the choice m = 31 and n = 63 in (12) with
d(i) ∈ {0; 1; 2; 4; 5; 7; 8; 9; 10; 14; 15; 16; 17; 18; 20; 27; 28; 30; 32; 34; 35; 36; 39; 40; 45; 49; 51; 54; 56; 57; 60}; and we will refer it to as SINGER31x63. 5.2
Biangular Tight Frame Based on Chirp Functions Following constructions described in [Strohmer & Heath, 2003, Applebaum et al., 2009], we consider normalized tight frame determined as
1 (A)i;j = √m exp ( 2 i m
) (ij1 + i2j2) ; 0 ≤ i; j1; j2 ≤ m − 1; j = j1 + mj2; (13) where m is a prime number not smaller than 5, and n = m2. Here we use m = 17, and the resulting normalized tight frame is referred below as CHIRP17x289. Note that the pairwise correlation between columns of such matrices takes values 0 or 1=√m, i.e. the frame is biangular. 5.3
Equivolume Tight Frame Based on Conference Matrices The following construction for n = 2m, where m = 2q, and q = 1; 2; : : :was proposed in [Kaporin, 2015]. These matrices have an unique property of coincidence of values of all principal minors of the third order in the Gram matrix AH A. (Recall that ETF property is equivalent to the coincidence of valied for all second order minors of AH A). Thus, A has optimum RIP constants 1(k) and 2(k) not only for k = 2, but also for k = 3. Such matrix A is constructed as follows: where
A = [( + (m 1) )Im + i( + )Cm j ((m 1) i )Im + ( + i )Cm]; = √ 2m 1 + p2m2 2m m ; = √ 2m 1 2m2 p2m2 m m ;
1 = p2m 1 ; and the skew-symmetric conference matrices Cm are determined recursively as
C2 = [ 0 1 Simultaneous change of sign at √2m2 − m in the expressions for and gives another family of equivolume tight frames with the same properties. For testing, we used m = 32 and the corresponding matrix is referred to as CONFSK32x64.
In Tables below, we present the following three types of data:
(i) the exact values of emxianct and emxaaxct representing the RIP bounds (3) and (4), respectively, obtained (if available) by the direct search over all possible index subsets J (shown in the 2nd and 7th columns of tables); (ii) the approximations imteinr and imtearx for the RIP bounds (3) and (4), respectively, obtained by the use of the generalized power method of Section 3 with lmax = 10000 and tmax = 200 (shown in the 3rd and 6th columns of tables);
(iii) the values of theoretical estimates emsitn and emsatx for the RIP bounds (3) expressed via the extreme roots of the Jacobi polynomials as asserts Theorem in Section 4 (shown in the 4th and 5th columns of tables).
Note that, by the definition, for every fixed k it holds exact iter min ≤ min ≤ emsitn < est iter max ≤ max ≤ emxaaxct: 6.1
Test Results for SINGER17x273 and CHIRP17x289 An interesting observation following from the data presented in Tables 1 and 2 is that the equiangularity property is not sufficient for the best conditioning of the tight frame (compare the data for k = 6). One can also notice that the theoretical bounds emsitn and emsatx are getting the more and more loose as the sparsity parameter k increases. Additionally, the four last lines in the Table clearly indicate that the RIP property cannot be adequately characterised by the single quantity (k) since max(AJH AJ ) may simply be larger than 2.
The approximations to the RIP spectral bounds obtained by the generalized power method are remarkably close to their exact values (whenever the latter are available). the current frame-related research is biased towards various generalizations of equiangular frame constructions such as fusion frames, biangilar frames, etc.
In the case n ≈ 2m, one observes a much better agreement of theoretical bounds emsitn and est max with the exact spectral bounds.
In comparison with the cases n ≈ m2, one can notice an even better quality of the approximations obtained by the generalized power method for the RIP spectral bounds, especially for the test case CONFSK32x64. First one uses a uniform sweep over (0,1) in order to roughly localize the desirable range of values for the penalty parameter (e.g. resulting in 1 < k < m). Next one uses the results of the first rough sweep to determine a much more tight range of values [ 1; 2] for the penalty parameter: 7
Based on analysis of exactly evaluated RIP constants in the standard restricted isometry condition, it was observed that the conditioning of equiangular tight frames based on partial DFT and difference sets may be inferior compared to other frame designs. This well agrees with the known fact that equiangularity alone does not guarantee good RIP constants in general.
For the larger values of k (when the direct search for the exact spectral bounds has prohibitely large costs), a version of the generalized power method [Journee et al., 2010] adjusted for the processing of complex-valued data was developed, implemented and successively tested. A future research may be directed towards the certification of optimality for the obtained solutions using SDP relaxation techniques similar to [d’Aspremont et al., 2008]. Also, certain more efficient strategies for the choice of parameter and initial guess x(0) may likely be found.
This work was supported by the Russian Foundation for Basic Research (project no. 17–07–00510) and the grant NSh–8860.2016.1.