{'fraction_non_alphanumeric': 0.04060568734470236, 'fraction_numerical': 0.008558625522968122, 'mean_word_length': 4.3351461590754585, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Blind Compressed Sensing Over a Structured Union of Subspaces

This paper addresses the problem of simultaneous signal recovery and dictionary learning based on compressive measurements. Multiple signals are analyzed jointly, with multiple sensing matrices, under the assumption that the unknown signals come from a union of a small number of disjoint subspaces. This problem is important, for instance, in image inpainting applications, in which the multiple signals are constituted by (incomplete) image patches taken from the overall image. This work extends standard dictionary learning and block-sparse dictionary optimization, by considering compressive measurements, e.g., incomplete data). Previous work on blind compressed sensing is also generalized by using multiple sensing matrices and relaxing some of the restrictions on the learned dictionary. Drawing on results developed in the context of matrix completion, it is proven that both the dictionary and signals can be recovered with high probability from compressed measurements. The solution is unique up to block permutations and invertible linear transformations of the dictionary atoms. The recovery is contingent on the number of measurements per signal and the number of signals being sufficiently large; bounds are derived for these quantities. In addition, this paper presents a computationally practical algorithm that performs dictionary learning and signal recovery, and establishes conditions for its convergence to a local optimum. Experimental results for image inpainting demonstrate the capabilities of the method.


I. INTRODUCTION
The problem of learning a dictionary for a set of signals has received considerable attention in recent years. This problem is known to be ill-posed in general, unless constraints are imposed on the dictionary and signals. One such constraint is the assumption that the signals x i ∈ R n , i = 1, . . . , N , can be sparsely represented under an unknown dictionary, i.e., each vector x i can be written as where D ∈ R n×r is the dictionary, and s i ∈ R r are sparse coefficient vectors satisfying s i 0 r. The residual energy 2 is assumed to be small. For example, [1] derives conditions to ensure the uniqueness of D and the representations s i , for the case in which all the coordinates of the signals x i are observed. In this setting, the problem is known as dictionary learning (DL) or, sometimes, collaborative DL to emphasize the fact that multiple signals are considered. Note that, in DL and in most related literature, "uniqueness" is defined up to an equivalence class involving permutations and rotations of the dictionary atoms; we follow the same convention throughout this paper. In the DL setting, several algorithms have been proposed for estimating D and s i . Examples includes K-SVD [1] and the method of optimal directions (MOD) [18], both of which enjoy local convergence properties. More recently, [16] has derived local convergence guarantees for 1 minimization applied to DL when the signals are sparse.
A more structured type of sparsity is considered in blocksparse dictionary optimization (BSDO) [26], in which it is assumed that the nonzero coordinates of s i (active atoms) occur in blocks rather than in arbitrary locations. This property is called block sparsity [12], and is important for the analysis of signals that belong to a union of a small number of subspaces, as described in [13]. The standard BSDO framework, like DL, assumes that all coordinates of the signals x i are observed. The block structure, i.e., the number of atoms and composition of each block of the dictionary, is in general not known a priori and should be estimated as part of the DL process. The BSDO algorithm reduces to K-SVD when the block size is one (standard sparsity).
While, technically, the set of all k-sparse signals in R r is itself a union of r k subspaces, it greatly simplifies the problem if one considers the more structured setting of block sparsity. This reduces the number of subspaces to L K , where L is the total number of blocks in the dictionary and K is the number of blocks that are active (block sparsity reduces to standard sparsity when all blocks are singletons). Block sparsity is closely related to the group LASSO in statistics literature [20], and also to the mixture of factor analyzers (MFA) statistical model, as noted in [8]. The particular case for which only one block is active is called one-block sparsity and corresponds to a union of L 1 = L possible subspaces, which is a dramatic simplification compared to the unstructured sparsity case. While this might at first sound limiting, the authors in [30] have obtained state-of-the-art image restoration results with a one-block sparsity model, and a variety of alternative methods have been proposed for this special case [14], [29], [30].
In applications we often do not have access to data x i , but only to compressive measurements where A ∈ R m×n has m < n rows, so that y i ∈ R m . In order to enable recovery of s i from y i even when D is known, the sensing matrix A must satisfy incoherence properties with respect to D, as prescribed by compressed sensing (CS) theory [6], [9]. Learning D from compressive measurements is called blind compressed sensing (blind CS), and does not admit a unique solution unless additional structure is imposed on D [17]. See also the approach in [10] for simultaneous dictionary learning and sensing matrix design in the CS scenario.
In this paper, we address simultaneous estimation of oneblock-sparse signals and the corresponding dictionary, given only compressive measurements. This unifies blind CS and BSDO (for the one-block-sparse case in particular). It is known that the standard blind CS problem, where a fixed sensing matrix A is used, does not admit a unique solution in general, although a number of special cases of dictionaries for which such a solution exists have been identified in [17]. These special cases are: (i) finite sets of bases (i.e., it is known that the dictionary is one member of a finite set of known bases for R n ); (ii) sparse dictionaries (i.e., the atoms of the dictionary themselves admit sparse representations) and (iii) block-diagonal dictionaries, with each block composed of orthogonal columns.
We are motivated in part by the problem of inpainting and interpolating an image [2], [31], where one observes an incomplete image, i.e., we only know the intensity values at a subset of pixel locations (or a subset of their linear combinations). Additionally, the image is processed in (often overlapping) patches, which we convert to n-dimensional vectors (our signals of interest). We observe m i < n pixels in each patch, indexed by i, with these m i pixels selected at random. Therefore, the locations of the missing pixels are in general (at least partially) different for each patch. This means that, unlike the classical blind CS setting, the sensing matrix is not the same for all signals; we denote the sensing matrix for patch i as A i ∈ R mi×n . The assumption of multiple and distinct A i is crucial for solving the problem of interest, as it enables the use of existing results from matrix completion [24]. Moreover, it has been demonstrated in [29] that using multiple A i improves inpainting and interpolation performance.
In conventional image inpainting, each A i consists of a randomly chosen subset of m i rows of the n × n identity matrix (thereby selecting m i pixels). Successful estimation of the missing pixel intensity values is contingent on each patch being representable in terms of a small subset of the columns in the dictionary D [31]. In other words, there needs to exist some D such that the s i are (at least approximately) sparse. However, our analysis is not restricted to A i being defined as in conventional inpainting. Other constructions typically used in CS, such as matrices with i.i.d. random-subgaussian [27] entries, or with rows drawn uniformly at random from an orthobasis [6], can be employed. Moreover, in many settings it is appropriate to assume that the patches belong to a union of disjoint subspaces [23], [30], with the number of subspaces being small. This motivates the one-block-sparsity assumption underlying our theoretical and computational results. This assumption, coupled with the use of multiple A i , allows us to ensure recovery of D and the s i under milder conditions on D and on the number of vectors y i , as compared to standard blind CS.
We show that unique recovery of the dictionary and the oneblock-sparse signals is guaranteed with high probability, albeit with high computational effort, if the number of measurements per signal and the number of signals are sufficiently large. We derive algorithm-independent bounds for these quantities, thereby extending DL by considering compressive measurements and establishing a connection with matrix-completion theory. Our results reduce to those known for DL when the signals are fully observed.
Additionally, we present a computationally feasible algorithm that performs DL and recovery of block-sparse signals based on compressive measurements with multiple sensing matrices. We automatically learn the block structure, in the same way as the BSDO algorithm [26]; the size (number of atoms) and composition of each block is not known a priori and we only need to specify a maximum block size. Our algorithm is closely related to BSDO, the main differences being that one-block sparsity and compressive measurements are considered. The estimates of D, s i and the block structure are found by alternating least-squares minimization.
It is shown that the algorithm converges to a local optimum under mild conditions, which we derive. This convergence analysis is not available for most other one-block-sparsity promoting methods, such as [8], or for most methods that rely on standard sparsity, such as [31]. An exception is the analysis in [16] pertaining to 1 minimization, where local convergence is proven when the number of measurements is at least O(n 3 k) for signals with sparsity level k. However, compressive measurements are not considered. Besides our global uniqueness result, our method provably attains local convergence for, at most, O(nk log n) measurements, although our one-block sparsity assumption is stronger than in [16]. The approach proposed in [30] also converges, but involves an even stronger assumption akin to intra-block sparsity, and does not include an algorithm-independent uniqueness analysis. Compelling experimental results are presented below, demonstrating the ability of our algorithm to perform inpainting of real images with a significant fraction of missing pixels.
The remainder of the paper is organized as follows. Section II presents preliminary definitions and a formulation of the optimization problem. Our uniqueness result is presented in Section III. Sections IV and V respectively describe the proposed algorithm and the corresponding proof of convergence to a local optimum. Experimental results are described in Section VI and concluding remarks are given in Section VII.

A. Preliminaries
Assume vectors y i ∈ R mi , i = 1, . . . , N are observed, such that where x i ∈ R n is an unknown signal, A i ∈ R mi×n is a known sensing matrix, D ∈ R n×r is an unknown dictionary, s i ∈ R r is an unknown sparse vector of coefficients such that x i = Ds i + i , with small residual i 2 . We focus on the noiseless case, although our analysis can be extended straightforwardly to include observation noise, following [4]. Given y i we would like to estimate x i by finding D and s i . Thus, we are interested in solving s.t. s i one-block sparse with our estimates denoted D, s 1 , . . . , s N . It is assumed that each x i lives in a subspace specified by a subset of the columns of D, so that there exists a permutation of the columns such that x i = Ds i and the coefficient vector s i is one-block sparse, as defined below. It is also assumed that D is composed of blocks (subsets of columns) corresponding to the blocks of s i , with the cardinality of the blocks summing to r. The cardinality and composition of each block are unknown, although we fix a maximum cardinality as in BSDO. The atoms in each block are assumed orthonormal, in order to avoid scaling indeterminacies and also to obviate concerns with sub-coherence (see [12] for a discussion). The expression in (4), particularly the 2 norm, is equivalent to a maximum-likelihood solution for the unknown model parameters, assuming that the components of i are i.i.d. Gaussian and the model deviation is negligible.
In general, the block size k is not known a priori, although it is common to define a maximum size k max .
We are specifically interested in one-block sparsity, where K = 1. An interpretation of one-block sparsity is that it corresponds to signals that live in a union of L linear subspaces, with each subspace spanned by the columns of one block D[ ]. If block D[ ] has k columns, then its subspace is of dimension k . Figure 1 illustrates the concept of oneblock sparsity. The index set of all signals which use block is ω = {i : s i [ ] = 0}. Note that for one-block-sparse signals there can be no overlap between ω 1 and ω 2 for 1 = 2 .
The optimization problem (4) can be decomposed as  { s 1 s 2 s 3 Illustration of the concept of one-block sparsity. The figure represents a dictionary with r = 16 columns, arranged consecutively in blocks D [1], . . . , D [4], all of size four. Also represented is a coefficient matrix with columns s 1 , s 2 , s 3 and s 4 . Non-zero coefficient values are colored, while zero-valued coefficient values are blank. All signals are one-block sparse, as they only use one block each. Note that there exists no overlap between index sets ω 1 , ω 2 , ω 3 and ω 4 . Our first goal is to develop conditions under which (4) admits a unique solution. To this end we note that the index sets ω are unknown, and therefore estimates ω must be found. This is due to the fact that, a priori, we do not know the correct assignments of columns to blocks. In our uniqueness analysis, we follow the same strategy as [1] and assume that all possible index sets can be tested. This is admittedly impractical in general, and it is done only for the purpose of constructing the uniqueness proof. The algorithm proposed in Section IV provides estimates of ω without exhaustive search.

class. Definition 2 (Equivalence class for
For a given ω , (5) is an instance of the weighted orthonormal Procrustes problem (WOPP), which in general requires iterative optimization without guarantee of convergence to a global optimum [21]. In the following, we establish conditions under which uniqueness can be guaranteed. To this end we show that (5) can be written as a matrix-completion problem, and then rely on known results developed in that context.

B. Formulation as a matrix completion problem
To reformulate (5), define X ω ∈ R n×|ω | as the matrix whose columns are the signals x i that use block , i.e., with i ∈ ω , and similarly define S ω [ ] ∈ R k ×|ω | as the matrix whose columns are the coefficient blocks Since D[ ] has only k columns and S ω [ ] has k rows, it is clear that X ω = D[ ] S ω [ ] has rank at most k . In fact, unless the null space of D[ ] intersects the range of S ω , the rank is exactly k . This should not happen in general, barring degeneracies that are precluded by the conditions assumed in our results. This suggests the strategy of treating each subproblem (5), assuming ω is known and hence the subspace selection has been determined (in our method we address estimation of ω ), as a low-rank matrix completion problem, as we illustrate next.
DefineÃ ∈ R M ×n with M ≥ n, as the matrix constructed by taking the union of the unique rows of all sensing matrices A i , for i ∈ ω (in the interest of notation brevity, we omit the dependence ofÃ on ). For example, in the aforementioned inpainting problem, for which each A i is defined by selecting rows at random from the n × n identity matrix, denoted I n×n , we haveÃ equal to a subset of I n×n and in the limit, given enough A i ,Ã = I n×n (up to row permutation). However, the discussion below considers more-general A i , for example composed in terms of draws from a subgaussian distribution. Let Y ω ∈ R M ×|ω | be defined as Y ω =Ã X ω . When performing measurements, we do not observe all elements of Y ω ; for each column i ∈ ω , we only have access to the entries selected by A i , these corresponding to the matching observed vector y i . This is illustrated in Figure 2, which also shows a pictorial comparison with the cases of fully measured signals (DL) and a single measurement matrix A (standard blind CS).
Denoting the locations of the observed entries of Y ω as the observation model can be written as where P Ω (·) is the operator that extracts the values of its argument at locations indexed by Ω. Each subproblem (5) can then be reformulated as where · F is the Frobenius norm. If X ω has rank k (the block size), assumingÃ has rank greater or equal to k and does not reduce the rank of X ω (since k is generally small andÃ is assumed to have rank n, these conditions forÃ are anticipated to be typical), we have that Y ω is also of rank k . This establishes the fact that, when X ω has low rank, so does Y ω (matrix completion results are most useful when the rank is low, although they hold for any value of the rank). Each subproblem (8) can then be solved by first completing the matrix Y ω and then obtaining X ω and, subsequently, D[ ] and S ω [ ]. Note that, although this is a two-step process, if A has rank n then it can be inverted to find a unique mapping from any estimated complete matrix Y ω to a corresponding estimated X ω . As long as X ω can be correctly estimated, the solution to (8) can then be found by a singular value decomposition (SVD) of X ω . This is due to the fact that, by the Eckart-Young theorem [11]. Following matrix completion theory [5], [7], [24], [25], if Y ω is truly of low rank (and additional technical conditions are met), it can be correctly completed with high probability by solving the convex program where · * is the nuclear norm, which is defined for a generic matrix Z of rank k as the sum of its singular values, i.e., with γ i (Z) indicating the i-th singular value of Z. Importantly, problem (9) is not equivalent to (8), since a solution of (8) may not be a solution of (9). This is a reflection of the fact that there may exist a high-rank Y ω which solves (8) but has high nuclear norm. However, the opposite holds: a solution of (9) is always a solution of (8), and therefore of (5). Moreover, the solution of (9) will have low rank and therefore produce D[ ] with the smallest possible block size k , which is clearly desirable. Thus, the nuclear norm formulation yields those solutions of the original problem (8) that are useful.
In the following section, we state conditions on the number of observed entries and onÃ for successful recovery of Y ω , X ω and, subsequently, of D[ ] and s 1 [ ], . . . , s N [ ], assuming exhaustive search over ω , by exploiting the connection between (5), (8) and (9). In Section IV we propose an algorithm that does not require exhaustive search and is guaranteed to converge to a local minimum under mild conditions here established. The notation used in this paper is summarized in Table 1 below, for ease of reference.
Dictionary block , i.e. the -th subset of the columns of D

III. UNIQUENESS RESULT
Our results rely on the concepts of spark and coherence, which we now define. The spark σ(·) of a matrix is the smallest number of columns that are linearly dependent; the number of linearly independent columns is the rank. Additionally, let U be a subspace of R n with dimension k and spanned by vectors {z u } u=1,...,k . The coherence of U vis-à-vis the standard basis {e v } v=1,...,n is defined as [24] This quantity is very similar to the coherence of a matrix, as defined in [28]. Our results also build on the following DL uniqueness conditions. Definition 3 (Uniqueness conditions for DL [1]).
• Support: There exist at least k + 1 signals for every possible combination of k atoms from D. With regular sparsity level k , this amounts to at least (k + 1) r k signals. For one-block sparsity with L blocks of size k , = 1, . . . , L, however, we need L =1 (k +1) signals, which is typically a far smaller number. • Non-degeneracy: Given k + 1 signals from the same combination of k atoms, their rank is exactly k (general position within the subspace). Similarly, any k +1 signals from different combinations must have rank k + 1. These conditions apply to the case of fully measured signals and constitute a limiting case, since our problem reduces to DL when all A i are equal to the identity matrix. Even in this limiting case, however, the one-block sparsity condition has the advantage of greatly reducing the number of possible subspaces and therefore the required number of signals, as stated in the richness property definition above. Note that we do not require prior knowledge of the block size k , as explained below. We now present our main result.
Theorem 1 (Uniqueness conditions for blind CS with multiple measurement matrices). Let y i ∈ R mi , with i = 1, . . . , N , be a set of observed vectors, obtained through projections of unknown signals x i ∈ R n so that each y i = A i x i and A i ∈ R mi×n is a known sensing matrix. Furthermore, let x i = Ds i where s i ∈ R r is an unknown vector of coefficients obeying one-block sparsity according to Definition 1 and D ∈ R n×r is an unknown dictionary having a block structure such that D = The vectors x i , i ∈ ω , are non-degenerate, i.e., any subset of k ≤ k vectors span a subspace of rank at least k. (iv) 32µ k (M 1 + M 2 )β log(2M 2 ) entries of the matrix Y ω are observed uniformly at random. The total number of observed entries is, thus, Proof of Theorem 1: Condition (i)-(iii) are analogous to the DL uniqueness conditions in Definition 3, which are proven in [1], and are always required even if all y i are fully measured. Condition (i) is adapted to our setting by imposing the spark restriction onÃD; this condition ensures that no two dictionary blocks are linearly dependent, and that the sensing matrices are sufficiently incoherent with respect to the dictionary. This is the case, with high probabilty, wheñ A obeys a random subgaussian or orthobasis construction [6], [9].
Condition (ii) stems from the fact that any k -dimensional hyperplane is uniquely determined by k +1 vectors, and since a subspace will always contain the origin, only k additional vectors are required. We write |ω | > k instead of |ω | ≥ k because, in the setting of one-block sparsity, there exist N = L =1 |ω | vectors from a union of L different subspaces, and we need to be able to assign vectors to their respective subspaces.
If condition (iii) holds, which precludes spurious colinearities or coplanarities, then the assignment can be done by the (admittedly impractical) procedure of testing the rank of all possible N k +1 matrices constructed by concatenating subsets of k + 1 column vectors, as assumed in [1]. If the rank of such a matrix is k , then all its column vectors belong to the same subspace. Otherwise, the rank will be k + 1. It is not necessary to know k a priori, provided that we test all possible values k = 1, . . . , n − 1 in the following way. Start by setting k = 1 and testing all possible subsets containing k + 1 = 2 signals. If any such subset has rank one, then we have found a (singleton) block. Any other possibility is precluded by the non-degeneracy condition (iii). If multiple subsets have rank one, then test the subspace angles in order to determine how many distinct singleton blocks exist. Record the corresponding blocks and signal assignments and remove the signals involved from further consideration (due to oneblock sparsity, each signal belongs to only one block). Iterate this procedure until either: (a) we have exhausted the signals; (b) all r atoms have been clustered into blocks, or (c) we have reached k = n.
We now address Condition (iv). We have |ω | > k vectors from subspace and we assume thatÃ has rank n and is, therefore, invertible. As discussed in Section II, Y ω is an incomplete low-rank matrix. Under conditions studied in [24], the convex program (9) will recover the complete Y ω with high probability. Therefore, by showing uniqueness of the complete low-rank matrix Y ω , we can show the uniqueness of the corresponding subspace spanned by D[ ], which can be found by first solving the optimization (9) and then taking the estimate Y ω and performing the SVD ofÃ + Y ω .
We now invoke Theorem 1.1 in [24] which states that, under the conditions and with the probability specified in our theorem, the minimizer of problem (9) is unique and equal to the true Y ω . This concludes the proof.
Note that the aforementioned theorem from [24], like preceding results [7], is based on Ω having at least one entry for each row and for each column; this is already accounted for in the derivation of the bound on |Ω| and associated probability of successful recovery under a uniformly-at-random pattern for the missing data. If there is a fixed A i = A (with less than n rows), or if there is any row or column entirely missing from Ω, then there is no recovery guarantee. This is the case of standard blind CS, where there is a fixed A which is the same for all signals, and thus there are entire rows without observations, as illustrated in Figure 2. Hence, blind CS requires additional constraints on D, as explained in [17]. Using multiple A i allows us to avoid those constraints. Also, even when M ≥ n, we are still subsampling because we do not observe all entries of Y ω .
The limiting case when all blocks are singletons, i.e., k = 1 for all , coupled with the one-block sparsity assumption, means that each patch belongs to a straight line in R n , which is a very strong restriction. On the other hand, if all k = 1 but the one-block sparsity assumption is removed and there are K active blocks (singleton atoms when k = 1), then we revert to standard sparsity; Theorem 1 still applies, but we would need to be test r K possible combinations of atoms, which may be an extremely large number. In fact, the conditions in Theorem 1, with standard sparsity and additionally with all A = I n×n , reduce to those in DL. One-block sparsity makes Theorem 1 more appealing, since then there are only L possible combinations instead of r K . Under Theorem 1, recovery is contingent on having the correct clustering given by ω . In the above-presented analysis, it is assumed that we have the computational ability to exhaustively search over all index sets ω . Under the conditions of Theorem 1, this search will achieve the correct assignment of signals to subspaces (clustering), and we need only concern ourselves with vectors that all come from the same subspace. This is done only for the purpose of constructing a proof, following the same strategy as in [1]. Also, like the DL uniqueness result in [1], Theorem 1 is overly pessimistic in the required number of measurements. Specifically, for image inpainting applications, where Y ω typically has far more columns then rows, the prescribed number of measurements can be excessive (often orders of magnitude larger than the total number of pixels). This reflects a limitation of the current state-of-the-art in matrix completion theory, and we stress that other uniqueness results, such as [1] for the simpler DL setting, suffer from the same problem. Nevertheless, Theorem 1 provides peace of mind by guaranteeing the existence of a unique optimal solution given enough measurements. It remains necessary to address the practical issue of finding a solution with reasonable computational effort, from realistic amounts of unclustered data.
Toward that end, in the next section we propose an algorithm that can estimate the clustering and each of the subspaces without combinatorial search, although it is not guaranteed to reach the globally optimum solution. The algorithm uses alternating least squares, similarly to BSDO, K-SVD and Fig. 3. Block diagram of the overall dictionary learning procedure. The analysis in this paper focuses on the shaded block, where we plug our Algorithm 1 instead of the BK-SVD step in [26].
MOD. It is computationally efficient and enjoyes local convergence guarantees given much fewer measurements than those prescribed by Theorem 1. In Section V we show convergence of the algorithm to a local optimum, as long as specified algorithm-specific conditions hold.

IV. ALGORITHM
Our algorithm expands on the iterative procedure described in [26]. This procedure alternates between two major steps: (i) inferring the block structure and coefficients using sparse agglomerative clustering (SAC) [26] and block-orthogonal matching pursuit (BOMP) [12], and (ii) updating the dictionary blocks using block-K-SVD (BK-SVD) [26]. The latter step assumes fully measured data. Briefly, SAC progressively merges pairs of blocks 1 , 2 for which the index sets ω 1 and ω 2 are most similar, up to an upper limit (k max ) on the block size. BOMP sequentially selects the blocks that best match the observed signals y i and can be viewed as a generalization of OMP [19], [22] for the block-sparse case. Similarly, BK-SVD is an extension of K-SVD for block-sparsifying dictionaries. A more detailed explanation of these methods may be found in [26] and the references therein.
We follow the procedure in [26], but replace the BK-SVD step by our Algorithm 1 in order to take compressive measurements into account, as the block diagram in Figure 3 illustrates. We employ alternating least-squares minimization to find D and S in a manner similar to existing methods such as K-SVD or MOD but adapted to the CS setting.
Recall that we are solving the optimization problem (5). Assuming the current estimate of ω and the block structure are fixed (as estimated in the previous SAC step), then we need to solve, for each ,

The optimization alternates between D[ ] and S ω [ ].
For now we will set the constraint D[ ] T D[ ] = I aside and focus on the unconstrained problem; we will return to the constraint later. A description of the two main algorithm steps is presented below, followed by the derivation of conditions for convergence.

A. Dictionary update
Holding s i fixed, we can analytically find D[ ] using properties of the Kronecker product and the vectorization operator, denoted as vec(·). This is similar to MOD [18] in that the entire dictionary block is updated at once. For a generic matrix Z, vec(Z) is the column vector constructed by vertically stacking the columns of Z. For simplicity of presentation, and without loss of generality, assume that all m i = m. LetỸ ω denote the matrix whose columns are the vectors y i , so that we have |vec(Ỹ ω )| = m|ω |. We rewrite the unconstrained problem over D[ ] as where B is defined as with ⊗ denoting the Kronecker product. This is a least-squares optimization problem for a standard system of linear equations. The solution is well-defined and unique (recall that we are holding the s i , and therefore B, fixed) as long as B has rank equal to k n, the number of unknowns in vec (D[ ]).

B. Coefficient update
After obtaining D[ ] for the current iteration, the coefficient vectors s i [ ] are found analytically by least squares as The factor This is equivalent to Gram-Schmidt orthogonalization, and will not change the locations of the nonzero elements of the s i , although the blocks D[ ] will not have orthogonal columns during the iterative procedure. This procedure has worked well in practice, and we therefore do not pursue a more direct constrained minimization.
In addition to the proposed Algorithm 1, we have also implemented and tested an alternative method where, instead of alternating least squares, we perform convex optimization of D and S jointly (given the estimate of ω provided by the current SAC step) by completing matrix Y ω , which is closer in spirit to the ideas underlying the uniqueness conditions in Theorem 1. However, existing general-purpose convex optimization packages do not readily scale up to matrices containing millions of entries as in our inpainting experiments, as noted in [3]. Using mixed nuclear/Frobenius norm minimization via the singular value thresholding (SVT) approach proposed in [3], we have attained reconstruction performance comparable to our Algorithm 1, but at much higher computational cost. We have also found that successful convergence with this alternative SVT-based approach depends on careful tuning of step-size and regularization parameters. In contrast, our alternating least squares method does not have any tuning parameters (other than the maximum block size k max ) and is computationally far less demanding.
The following result establishes convergence conditions for our proposed alternating least-squares algorithm.
Proposition 1 (Convergence of Algorithm 1). Let y i ∈ R mi , with i = 1, . . . , N be a set of vectors that satisfy the conditions stated in Theorem 1. Let A i ∈ R mi×N be sensing matrices and x i ∈ R n signals such that y i = A i x i . Assume each signal is of the form x i = Ds i , with D a dictionary comprised of blocks D [1], . . . , D[L], where each block has k atoms, and s i ∈ R r a one-block-sparse vector. Let ω be the index set of vectors x i from dictionary block , and let Y ω ∈ R M ×|ω | be an incomplete data matrix such that P Ω (Y ω ) = P Ω (ÃD[ ]S ω [ ]), with Ω the locations of observed entries, andÃ ∈ R M ×n the union of the rows of the A i . Then, the alternating minimization procedure described in Section IV will converge to a local minimum of (4), with probability specified below, if the following conditions hold: (i) For all and all i ∈ ω , s i 0 = k < σ(ÃD) if the rows of A i are selected uniformly at random from A and rank(Ã) = n. (iv) Each vector y i has m i ≥ k measured values. The probability of convergence, for each block , is equal to one with i.i.d. Gaussian A i , and equal to 1 − n −β+1 for random selection of rows fromÃ, where β is the constant such that |Ω| = βn log n.
Proof of Proposition 1: Since condition (i) is the same as in Theorem 1, we focus on conditions (ii)-(iv). As mentioned above, the analytic solution (15) is well-defined and unique as long as B has rank equal to k n. Examining (13) and (14) we can see that, if the s i [ ] are non-degenerate (i.e., rank(S ω [ ])=k , with S ω [ ] the subset of nonzero rows of S ω ), then the rank condition is equivalent to having the number of observed elements |vec(Ỹ ω )| = |Ω| ≥ k n and also that is, the vertical concatenation of all A i , having rank n. This is due to the fact that rank(s i ⊗A i ) = rank(s i )×rank(A i ), ∀i.
The number of observed values needed to ensure this rank condition depends on the probabilistic mechanism that generates the A i . We will analyze two such mechanisms: random Gaussian and random subset of projections. Without loss of generality, assume all m i = m (we can always treat m as a lower limit on the m i ). Random Gaussian: We assume that all A i have elements independently drawn from a Gaussian distribution. Then, with probability one, any subset of size n of the m|ω | rows of Γ are linearly independent. Therefore, the condition m|ω | ≥ n ensures rank(Γ) = n; since we need m|ω | ≥ k n (and therefore |ω | > n) in order to have enough observed entries anyway, the latter condition ensures that the rank of Γ is adequate. Note that this reasoning assumes that the A i do not share rows; if they do share rows, e.g., in the case when there is a finite pool of random-Gaussian projections and the rows of A i are subsets from the pool, then we must use the following result instead. Random subset of projections: We analyze the case for which the rows of each sensing matrix A i are randomly drawn from a pool of linearly independent rows, the union of which constitutesÃ. The typical example of this situation is when the rows ofÃ form a random orthobasis, although here orhonormality is not required (linear independence suffices). This includes the case whenÃ is the identity matrix, which is most relevant to the inpainting/interpolation problem. Note that the problem can be treated as an instance of the classic Coupon Collector problem (see, e.g., [15]). There are m|ω | rows (the number of rounds our "collector" draws a randomly chosen row, or "coupon") that have to span a space of n dimensions (the different types of "coupons" we need to collect). Denote as Z m|ω | i the event that no coupons of type i have been drawn after m|ω | rounds. It is easy to see that where the inequality comes from our ignoring, for simplicity, the fact that the rows are drawn without replacement within each sensing matrix (stated otherwise, there are no repeated rows within each A i ). Then, by the union bound, The left-hand side is precisely the probability that Γ is not rank n. To look at how this bound scales relative to n, we can apply the inequality Then, consider m|ω | = βn log n, with β constant, and define T to be the minimum value of m|ω | before we achieve full rank, so that P Therefore, we need O(n log n) observed values in order to achieve rank n with arbitrarily high probability. Hence, we need O(nk log n) for solving (13) for each block. Although we do not impose the orthogonality contraint directly in the minimization, any solution of (13) can be orthonormalized and thus become a solution of (12) without changing the subspace defined by D[ ].
Concerning the coefficient update expression (16), the issue of invertibility comes up due to the term Note that the matrix being inverted has size k × k . The rank of A T i A i is equal to m i , the number of rows of A i . For invertibility of the above expression, we need to have k ≤ m i , so that there is no loss of rank. This means that our method requires a lower bound on m i . We conclude the proof by invoking, as in [1], the fact that all of the steps of this alternating minimization method are optimal under the above assumptions, and hence cannot increase the value of the objective function. Therefore, since the objective function is bounded from below by zero, the algorithm will converge.
We have proven convergence in the objective function, which is a weak form of convergence in the sense that the objective function will attain a minimum but the estimates might not reach a stopping point. However, convergence in the objective function is still a useful result, and it is the same guarantee as K-SVD and MOD. In practice, the above result ensures that the required number of signals scales linearly with k and also linearly (up to a log factor in the case of measurements using an orthobasis ensemble) with signal dimensionality n. These are more favorable bounds than those in Theorem 1 for uniqueness. However, this is not surprising since Theorem 1 is based on matrix completion theory, which establishes slightly pessimistic sufficient conditions for strong unique global recovery guarantees, while Proposition 1 only establishes convergence of our algorithm to a local minimum. It is also interesting that, in the case of random subsets of projections, the O(k n log n) requirement matches the information-theoretic limit established in [7] for matrix completion. This similarity is due to the fact that, as in [7] and other work, we make use of the Coupon Collector model. Regarding the lower bound on m i , since we typically expect the block size to be small, the condition m i ≥ k is relatively mild and should only be of concern for extremely small measurement fractions.

VI. EXPERIMENTAL RESULTS
The algorithm proposed in Section IV is validated by inpainting the well-known "Barbara" (512 × 512 pixels) and "house" (256 × 256 pixels) images, for varying percentages of observed pixels (these two examples are representative of many others we have considered). The images are processed in 8 × 8 overlapping patches, treated as vectors of dimension n = 64. In Figures 4 and 5 we present the original images and the test versions, with 25%, 50% and 75% of the pixel values observed (selected uniformly at random). In all experiments, the total number of dictionary elements is set to r = 256. Figures 6 and 7 show the inpainting results achieved by our algorithm from 50% observed pixels, using maximum block sizes k max = 4 and k max = 8. The peak signal-to-noise ratios (PSNR) for this case are shown in Table II. The PSNR values for varying percentages of observed pixels are plotted in Figure 8, for k max equal to four and eight.
Results are averaged over ten runs with different random locations for the missing pixels, with error bars showing the one standard deviation interval. These experiments are intended to validate our method, rather than claiming outperformance of the state-of-the-art in image inpainting and interpolation. Our performance is comparable to that of the algorithm described in [31], which we have used on the same images. However, our model and algorithm are significantly simpler and, unlike [31], our algorithm has convergence guarantees. Also, while the PSNR achieved in our experiments is lower than in the state-of-the-art approach in [30], the results are not directly comparable due to the fact that additional structure is assumed in [30], including fast decay of the singular values within each block (which is analogous to internal sparsity). Like other authors (e.g., [16]), we have observed a phase transition phenomenon where the reconstruction performance falls sharply for measurement percentages under a rankdependent threshold (near 40% observed pixels in our case). This is shown in Figure 9, which was obtained using as input subsets of pixels from a union-of-subspaces approximation of the "Barbara" image for block sizes four and eight, thus matching our model exactly. We plot the empirical frequency of successful reconstructions over ten realizations of the missing pixel locations, for varying measurement percentages. The reconstruction is deemed successful when the PSNR exceeds 40 dB.
In addition, the learned dictionaries for k max = 4 and k max = 8 are depicted in Figures 10 and 11. It is possible to observe that the atoms within each dictionary block are more similar to each other than those in different blocks. This is expected and desired clustering behavior, as one would like the blocks to reflect different image properties. This behavior is further illustrated in Figure 12, where the index of the block used by each patch is represented by a corresponding color (to avoid clutter, this is shown for the most frequently used dictionary blocks). It is apparent that patches with similar texture tend to share blocks.
In all cases tested, convergence occurs after a maximum of ten iterations. A non-optimized MATLAB implementation of Algorithm 1, including the SAC/BOMP steps, requires approximately 12 hours for inpainting the full 512 × 512 "Barbara" image (worst case) on a computer with CPU clock frequency of 2.52 GHz. This computation time is similar to the BSDO algorithm described in [26]. We employ overlapping patches in order to avoid block artifacts, leading to a total of N = 255, 025 vectors being processed for the "Barbara" image, with N = 62, 001 for the "house" image. Note that the overlapping procedure follows the same "non-overlapping sensing with overlapping reconstruction" strategy as [29], [31]; this procedure employs a distinct sensing matrix per nonoverlapping patch, with the image reconstruction performed by averaging the overlapping estimated patches.  Fig. 9. Illustration of the phase transition phenomenon, using as input subsets of pixels from an approximation of the "Barbara" image that exactly obeys the union-of-subspaces model, with block size (rank) four and eight. The plot shows the normalized frequency of successful reconstructions over ten realizations of missing pixel locations, for varying percentages of observed pixels. We declare a successful reconstruction when the PSNR of the estimate is greater than 40 dB.

VII. CONCLUSION
We have proposed a framework for simultaneous estimation of one-block-sparse signals and the corresponding dictionary, based on compressive measurements. Multiple sensing matrices are employed, which allows use of low-rank matrix completion results to guarantee unique recovery (up to a specified equivalence class) for a sufficiently large number of signals and measurements per signal; bounds are derived for the number of measurements. The assumption of oneblock sparsity is related to the union-of-subspaces signal model and to the well-known mixture of factor analyzers (MFA) statistical framework. Existing results for the related problems of dictionary learning (DL) and block-sparse dictionary optimization (BSDO) are extended, due to the analysis of compressive measurements, and blind compressed sensing (blind CS) is also extended by considering a broader class of dictionaries.
Additionally, a practical algorithm based on alternating least-squares minimization has been proposed. The algorithm is related to BSDO, and has been shown to converge to a local optimum under mild conditions. We observe encouraging performance on image inpainting tasks without the need for parameter tuning or careful initialization. An avenue for further research is the treatment of observation noise, which constitutes a natural extension of our work, as well as extension of the theory beyone one-block sparsity.    6. Inpainted 512 × 512 "Barbara" image with 50% observed pixels. The peak signal-to-noise ratio (PSNR) of this estimate is 27.93 dB. Maximum block size (kmax) is eight and the total of dictionary elements is set to r = 256, divided in L = 32 blocks. Image patches have size 8 × 8 and are treated as vectors of dimension n = 64. As we employ overlapping patches, the total number of vectors is N = 255, 025.   Fig. 8. Peak signal-to-noise ratio (PSNR) for inpainting the "Barbara" (left) and "house" (right) images with 25%, 50%, 75% and 100% observed pixel values, using maximum block size (kmax) of four and eight. Results are averaged over ten realizations of the randomly missing pixel locations. The error bars depict one standard deviation. Using higher kmax yields better PSNR (but also a less parsimonious model).    Fig. 12. Dictionary block usage for the "Barbara" image with kmax = 8. The colors represent the index of the block used by each patch, for a subset containing the 22 most frequently used blocks, out of a total of 32. The subset is shown to the right, aligned with the indices of the colorbar, with each block corresponding of a row of subimages (atoms). It is apparent that patches with similar textures (e.g., the pants, portions of the scarf and table cloth) tend to share blocks (best viewed in color).