Compressed Sensing and Sparse Recovery

6 Compressed Sensing and Sparse Recovery Most of us have noticed how saving an image in JPEG dramatically reduces the space it occupies in our hard drives (as oppose to file types that save the pixel value of each pixel in the image). The idea behind these compression methods is to exploit known structure in the images; although our cameras will record the pixel value (even three values in RGB) for each pixel, it is clear that most collections of pixel values will not correspond to pictures that we would expect to see. This special structure tends to exploited via sparsity. Indeed, natural images are known to be sparse in certain bases (such as the wavelet bases) and this is the core idea behind JPEG (actually, JPEG2000; JPEG uses a different basis). Let us think of x ∈ RN as the signal corresponding to the image already in the basis for which it is sparse. Let’s say that x is s-sparse, or kxk0 ≤ s, meaning that x has, at most, s non-zero components and, usually, s  N . The `0 norm26 kxk0 of a vector x is the number of non-zero entries of x. This means that, when we take a picture, our camera makes N measurements (each corresponding to a pixel) but then, after an appropriate change of basis, it keeps only s  N non-zero coefficients and drops the others. This motivates the question: “If only a few degrees of freedom are kept after compression, why not measure in a more efficient way and take considerably less than N measurements?”. This question is in the heart of Compressed Sensing [CRT06a, CRT06b, CT05, CT06, Don06, FR13]. It is particularly important in MRI imaging [?] as less measurements potentially means less measurement time. The following book is a great reference on Compressed Sensing [FR13]. More precisely, given a s-sparse vector x, we take s < M  N linear measurements yi = aTi x and the goal is to recover x from the underdetermined system:      y = A           x .         Last lecture we
used Gordon’s theorem to show that, if we took random measurements, on the order of s log Ns measurements suffice to have all considerably different s-sparse signals correspond
to considerably different sets of measurements. This suggests that ≈ s log Ns may be enough to recover x, we’ll see (later) in this lecture that this intuition is indeed correct. Since the system is underdetermined and we know x is sparse, the natural thing to try, in order to recover x, is to solve min kzk0 (55) s.t. Az = y, and hope that the optimal solution z corresponds to the signal in question x. Unfortunately, (55) is known to be a computationally hard problem in general. Instead, the approach usually P taken in sparse recovery is to consider a convex relaxation of the `0 norm, the `1 norm: kzk1 = N i=1 |zi |. Figure 19 26 The `0 norm is not actually a norm though. depicts how the `1 norm can be seen as a convex relaxation of the `0 norm and how it promotes sparsity. Figure 19: A two-dimensional depiction of `0 and `1 minimization. In `1 minimization (the picture of the right), one inflates the `1 ball (the diamond) until it hits the affine subspace of interest, this image conveys how this norm promotes sparsity, due to the pointy corners on sparse vectors. This motivates one to consider the following optimization problem (surrogate to (55)): min kzk1 s.t. Az = y, (56) In order for (56) we need two things, for the solution of it to be meaningful (hopefully to coincide with x) and for (56) to be efficiently solved. We will formulate (56) as a Linear Program (and thus show that it is efficiently solvable). Let us think of ω + as the positive part of x and ω − as the symmetric of the negative part of it, meaning that x = ω + − ω − and, for each i, either ωi− or ωi+ is zero. Note that, in that case, kxk1 = N X
ωi+ + ωi− = 1T ω + + ω − . i=1 Motivated by this we consider: min 1T (ω + + ω − ) s.t. A (ω + − ω − ) = y ω+ ≥ 0 ω − ≥ 0, (57) which is a linear program. It is not difficult to see that the optimal solution of (57) will indeed satisfy that, for each i, either ωi− or ωi+ is zero and it is indeed equivalent to (56). Since linear programs are efficiently solvable [VB04], this means that (56) efficiently. 6.1 Duality and exact recovery The goal now is to show that, under certain conditions, the solution of (56) indeed coincides with x. We will do this via duality (this approach is essentially the same as the one followed in [Fuc04] for the real case, and in [Tro05] for the complex case.) Let us start by presenting duality in Linear Programming with a game theoretic view point. The idea will be to reformulate (57) without constraints, by adding a dual player that wants to maximize the objective and would exploit a deviation from the original constraints (by, for example, giving the dual player a variable u and adding to to the objective uT (y − A (ω + − ω − ))). More precisely consider the following




T + − + T + − T − T + − min max 1 ω + ω v ω − − v ω + u y − A ω − ω . u ω+ + ω − v − ≥0 v ≥0 (58) Indeed, if the primal player (picking ω + and ω − and attempting to minimize the objective) picks variables that do not satisfy the original constraints, then the dual player (picking u, v + , and v − and trying to maximize the objective) will be able to make the objective value as large as possible. It is then clear that (57) = (58). Now image that we switch the order at which the players choose variable values, this can only benefit the primal player, that now gets to see the value of the dual variables before picking the primal variables, meaning that (58) ≥ (59), where (59) is given by:

T
T

max min 1T ω + + ω − − v + ω + − v − ω − + uT y − A ω + − ω − . u (59)
T +
T − + T − T max min 1 − v − A u ω + 1 − v + A u ω + uT y u (60) ω+ v + ≥0 ω − v − ≥0 Rewriting ω+ v + ≥0 ω − − v ≥0 With this formulation, it becomes clear that the dual players needs to set 1 − v + − AT u = 0, 1 − v − + AT u = 0 and thus (60) is equivalent to max u uT y v + ≥0 v − ≥0 1−v + −AT u=0 1−v − +AT u=0 or equivalently, maxu uT y s.t. −1 ≤ AT u ≤ 1. (61) The linear program (61) is known as the dual program to (57). The discussion above shows that (61) ≤ (57) which is known as weak duality. More remarkably, strong duality guarantees that the optimal values of the two programs match. There is a considerably easier way to show weak duality (although not as enlightening as the one above). If ω − and ω + are primal feasible and u is dual feasible, then




0 ≤ 1T − uT A ω + + 1T + uT A ω − = 1T ω + + ω − − uT A ω + − ω − = 1T ω + + ω − − uT y, (62) showing that (61) ≤ (57). 6.2 Finding a dual certificate In order to show that ω + − ω − = x is an optimal solution27 to (57), we will find a dual feasible point u for which the dual matches the value of ω + − ω − = x in the primal, u is known as a dual certificate or dual witness.

From (62) it is clear that u must satisfy 1T − uT A ω + = 0 and 1T + uT A ω − = 0, this is known as complementary slackness. This means that we must have the entries of AT u be +1 or −1 when x is non-zero (and be +1 when it is positive and −1 when it is negative), in other words
AT u S = sign (xS ) , where S = supp(x), and AT u ∞ ≤ 1 (in order to be dual feasible).
Remark 6.1 It is not difficult to see that if we further ask that AT u S c ∞ < 1 any optimal primal solution would have to have its support contained in the support of x. This observation gives us the following Lemma. Lemma 6.2 Consider the problem of sparse recovery discussed this lecture. Let S = supp(x), if AS is injective and there exists u ∈ RM such that
AT u S = sign (xS ) , and
AT u S c ∞ < 1, then x is the unique optimal solution to the `1 minimization program 56.
Since we know that AT u S = sign (xS ) (and that AS is injective), we’ll try to construct u by least
squares and hope that it satisfies AT u S c ∞ < 1. More precisely, we take u = ATS where ATS Corollary.
† = AS ATS AS
−1
† sign (xS ) , is the Moore Penrose pseudo-inverse of ATS . This gives the following Corollary 6.3 Consider the problem of sparse recovery discussed this lecture. Let S = supp(x), if AS is injective and  
− 1 sign (xS ) c < 1, ATS c AS ATS AS S ∞ then x is the unique optimal solution to the `1 minimization program 56. 27 For now we will focus on showing that it is an optimal solution, see Remark 6.1 for a brief discussion of how to strengthen the argument to show uniqueness Recall the definition of A ∈ RM ×N satisfying the restricted isometry property from last Lecture. Definition 6.4 (Restricted Isometry Property) A matrix A ∈ RM ×N is said to satisfy the (s, δ) restricted isometry property if (1 − δ) kxk2 ≤ kAxk2 ≤ (1 + δ) kxk2 , for all two s-sparse vectors x.
Last lecture (Lecture
5) we showed that if M  s log Ns and A ∈ RM ×N is drawn with i.i.d. 1 gaussian entries N 0, M then it will, with high probability, satisfy the (s, 1/3)-RIP. Note that, if A q
−1
−1 satisfies the (s, δ)-RIP then, for any |S| ≤ s one has kAS k ≤ 1 + 31 and l ATS AS = k ≤ 1 − 13 3 2, where k · k denotes the operator norm kBk = maxkxk=1 kBxk.
1 This means that, if we take A random with i.i.d. N 0, M entries then, for any |S ≤ s| we have that r
−1 13 √ √ T kAS AS AS sign (xS ) k ≤ 1 + = 3 s, 32 and because of the independency among the entries of A, AS c is independent of this vector and so for each j ∈ S c we have   2 
− 1 t 1 √ √ T T 3 st ≤ 2 exp − , Prob Aj AS AS AS sign (xS ) ≥ √ 2 M where Aj is the j-th column of A. Union bound gives 
−1 Prob ATS AS ATS AS sign (xS ) ∞ 1 √ √ ≥√ 3 st M  t2 ≤ 2N exp − 2   , which implies   √ 2  Prob  ATS AS
−1 ATS AS sign (xS )  ∞  ≥ 1 ≤ 2N exp − √M 3s 2    1 M  − 2 log(2N ) ,  = exp − 2 3s which means that we expect to exactly recover x via `1 minimization when M  s log(N ), similarly to what was predicted by Gordon’s Theorem last Lecture. 6.3 A different approach Given x a sparse vector, we want to show that x is the unique optimal solution to min kzk1 s.t. Az = y, (63) Let S = supp(x) and suppose that z 6= x is an optimal solution of the `1 minimization problem. Let v = z − x, it satisfies kv + xk1 ≤ kxk1 and A(v + x) = Ax, this means that Av = 0. Also, kxkS = kxk1 ≥ kv+xk1 = k (v + x)S k1 +kvS c k1 ≥ kxkS −kvS k1 +kvkS c , where the last inequality follows by triangular inequality. This means that kvS k1 ≥ kvS c k1 , but since |S|  N it is unlikely for A to have vectors in its nullspace that are this concentrated on such few entries. This motivates the following definition. Definition 6.5 (Null Space Property) A is said to satisfy the s-Null Space Property if, for all v in ker(A) (the nullspace of A) and all |S| = s we have kvS k1 < kvS c k1 . From the argument above, it is clear that if A satisfies the Null Space Property for s, then x will indeed be the unique optimal solution to (56). Also, now that recovery is formulated in terms of certain vectors not belonging to the nullspace of A, one could again resort to Gordon’s theorem. And indeed, Gordon’s Theorem can be used to understand the number of necessary measurements to do sparse recovery28 [CRPW12]. There is also an interesting
approach based on Integral Geometry [ALMT14]. As it turns out one can show that the 2s, 31 -RIP implies s-NSP [FR13]. We omit that proof as it does not appear to be as enlightening (or adaptable) as the approach that was shown here. 6.4 Partial Fourier matrices satisfying the Restricted Isometry Property While the results above are encouraging, rarely one has the capability of designing random gaussian measurements. A more realistic measurement design is to use rows of the Discrete Fourier Transform: Consider the random M × N matrix obtained by drawing rows uniformly with replacement from the N × N discrete Fourier transform matrix. It is known [CT06] that if M = Ωδ (K polylog N ), then the resulting partial Fourier matrix satisfies the restricted isometry property with high probability. A fundamental problem in compressed sensing is determining the order of the smallest number M of random rows necessary. To summarize the progress to date, Candès and Tao [CT06] first found that M = Ωδ (K log6 N ) rows suffice, then Rudelson and Vershynin [RV08] proved M = Ωδ (K log4 N ), and more recently, Bourgain [Bou14] achieved M = Ωδ (K log3 N ); Nelson, Price and Wootters [NPW14] also achieved M = Ωδ (K log3 N ), but using a slightly different measurement matrix. The latest result is due to Haviv and Regev [HR] giving an upper bound of M = Ωδ (K log2 k log N ). As far as lower bounds, in [BLM15] it was shown that M = Ωδ (K log N ) is necessary. This draws a contrast with random Gaussian matrices, where M = Ωδ (K log(N/K)) is known to suffice. Open Problem 6.1 Consider the random M × N matrix obtained by drawing rows uniformly with replacement from the N × N discrete Fourier transform matrix. How large does M need to be so that, with high probability, the result matrix satisfies the Restricted Isometry Property (for constant δ)? 6.5 Coherence and Gershgorin Circle Theorem Last lectures we discussed the problem of building deterministic RIP matrices (building deterministic RIP matrices is particularly important because checking whether a matrix is RIP is computationally 28 In these references the sets considered are slightly different than the one described here, as the goal is to ensure recovery of just one sparse vector, and not all of them simultaneously. hard [BDMS13, TP13]). Despite suboptimal, coherence based methods are still among the most popular ways of building RIP matrices, we’ll briefly describe some of the ideas involved here. Recall the definition of the Restricted Isometry Property (Definition 6.4). Essentially, it asks that any S ⊂ [N ], |S| ≤ s satisfies: (1 − δ)kxk2 ≤ kAS xk2 ≤ (1 + δ)kxk2 , for all x ∈ R|S| . This is equivalent to
xT ATS AS − I x max ≤ δ, x xT x or equivalently ATS AS − I ≤ δ. If the columns of A are unit-norm vectors (in RM ), then the diagonal of ATS AS is all-ones, this means that ATS AS − I consists only of the non-diagonal elements of ATS AS . If, moreover, for any two columns ai , aj , of A we have aTi aj ≤ µ for some µ then, Gershgorin’s circle theorem tells us that ATS AS − I ≤ δ(s − 1). More precisely, given a symmetric matrix B, the Gershgorin’s circle theorem [HJ85] tells that all of the eigenvalues in thoe so called Gershgorin discs (for each i, the Gershgorin disc n of B are contained P corresponds to λ : |λ − Bii | ≤ j=i 6 |Bij | . If B has zero diagonal, then this reads: kB k ≤ maxi |Bij |. Given a set of N vectors a1 , . . . , aN ∈ RM we define its worst-case coherence µ as µ = max aTi aj i6=j Given a set of unit-norm vectors a1 , . . . , aN RM with worst-case coherence µ, if we form a matrix
with these vectors as columns, then it will be (s, µ(s − 1)µ)-RIP, meaning that it will be s, 13 - RIP for s ≤ 13 µ1 . 6.5.1 Mutually Unbiased Bases We note that now we will consider our vectors to be complex valued, rather than real valued, but all of the results above hold for either case. Consider the following 2d vectors: the d vectors from the identity basis and the d orthonormal vectors corresponding to columns of the Discrete Fourier Transform F . Since these basis are both orthonormal the vectors in question are unit-norm and within the basis are orthogonal, it is also easy to see that the inner product between any two vectors, one from each basis, has absolute value √1d , meaning that the worst case coherence of this set of vectors is µ = √1d this corresponding matrix [I F ] √ is RIP for s ≈ d. It is easy to see that √1d coherence is the minimum possible between two orthonormal bases in Cd , such bases are called unbiased (and are important in Quantum Mechanics, see for example [BBRV01]) This motivates the question of how many orthonormal basis can be made simultaneously (or mutually) unbiased in Cd , such sets of bases are called mutually unbiased bases. Let M(d) denote the maximum number of such bases. It is known that M(d) ≤ d + 1 and that this upper bound is achievable when d is a prime power, however even determining the value of M(6) is open [BBRV01]. Open Problem 6.2 How many mutually unbiased bases are there in 6 dimensions? Is it true that M(6) < 7?29 6.5.2 Equiangular Tight Frames Another natural question is whether one can get better coherence (or more vectors) by relaxing the condition that the set of vectors have to be formed by taking orthonormal basis. A tight frame (see, for example, [CK12] for more on Frame Theory) is a set of N vectors in CM (with N ≥ M ) that spans CM “equally”. More precisely: Definition 6.6 (Tight Frame) v1 , . . . , vN ∈ CM is a tight frame if there exists a constant α such that N X |hvk , xi|2 = αkxk2 , ∀x∈CM , k=1 or equivalently N X vk vkT = αI. k=1 The smallest coherence of a set of N unit-norm vectors in M dimensions is bounded below by the Welch bound (see, for example, [BFMW13]) which reads: s N −M µ≥ . M (N − 1) One can check that, due to this limitation, deterministic constructions based on coherence cannot √ yield matrices that are RIP for s  M , known as the square-root bottleneck [BFMW13, Tao07]. There are constructions that achieve the Welch bound, known as Equiangular Tight Frames (ETFs), theseqare tight frames for which all inner products between pairs of vectors have the same modulus µ = N −M M (N −1) , meaning that they are “equiangular”. It is known that for there to exist an ETF in CM one needs N ≤ M 2 and ETF’s for which N = M 2 are important in Quantum Mechanics, and known as SIC-POVM’s. However, they are not known to exist in every dimension (see here for some recent computer experiments [SG10]). This is known as Zauner’s conjecture. Open Problem 6.3 Prove or disprove the SIC-POVM / Zauner’s conjecture: For any d, there exists an Equiangular tight frame with d2 vectors in Cd dimensions. (or, there exist d2 equiangular lines in Cd ). We note that this conjecture was recently shown [Chi15] for d = 17 and refer the reader to this interesting remark [Mix14c] on the description length of the constructions known for different dimensions. 29 The author thanks Bernat Guillen Pegueroles for suggesting this problem. 6.5.3 The Paley ETF There is a simple construction of an ETF made of 2M vectors in M dimensions (corresponding to a M × 2M matrix) known as the Paley ETF that is essentially a partial Discrete Fourier Transform matrix. While we refer the reader to [BFMW13] for the details the construction consists of picking rows of the p × p Discrete Fourier Transform (with p ∼ = 1 mod 4 a prime) with indices corresponding to quadratic residues modulo p. Just by coherence considerations this construction is known to be √ p RIP for s ≈ p but conjectured [BFMW13] to be RIP for s ≈ polylogp , which would be predicted if 30 the choice os rows was random (as discussed above) . Although partial conditional (conditioned on a number theory conjecture) progress on this conjecture has been made [BMM14] no unconditional √ result is known for s  p. This motivates the following Open Problem. Open Problem 6.4 Does the Paley Equiangular tight frame satisfy the Restricted Isometry Property pass the square root bottleneck? (even by logarithmic factors?). We note that [BMM14] shows that improving polynomially on this conjecture implies an improvement over the Paley clique number conjecture (Open Problem 8.4.) 6.6 The Kadison-Singer problem The Kadison-Singer problem (or the related Weaver’s conjecture) was one of the main questions in frame theory, it was solved (with a non-constructive proof) in the recent breakthrough of Marcus, Spielman, and Srivastava [MSS15b], using similar techniques to their earlier work [MSS15a]. Their theorem guarantees the existence of universal constants η ≥ 2 and θ > 0 s.t. for any tight frame ω1 , . . . , ωN ∈ CM satisfying kωk k ≤ 1 and N X ωk ωkT = ηI, k=1 there exists a partition of this tight frame S1 , S2 ⊂ [N ] such that each is “almost a tight frame” in the sense that, X ωk ωkT  (η − θ) I. k∈Sj However, a constructive prove is still not known and there is no known (polynomial time) method that is known to construct such partitions. Open Problem 6.5 Give a (polynomial time) construction of the tight frame partition satisfying the properties required in the Kadison-Singer problem (or the related Weaver’s conjecture). 30 We note that the quadratic residues are known to have pseudorandom properties, and indeed have been leveraged to reduce the randomness needed in certain RIP constructions [BFMM14]