Rémi Gribonval
email: remi.gribonval@inria.fr

Rodolphe Jenatton
email: jenat-ton@amazon.com

Francis Bach
email: francis.bach@ens.fr

Martin Kleinsteuber

Matthias Seibert

Sample Complexity of Dictionary Learning and other Matrix Factorizations

come   L'archive ouverte pluridisciplinaire

Sample Complexity of Dictionary Learning and other Matrix Factorizations

Rémi Gribonval, IEEE Fellow, Rodolphe Jenatton, Francis Bach, Martin Kleinsteuber, Matthias Seibert

Abstract-Many modern tools in machine learning and signal processing, such as sparse dictionary learning, principal component analysis (PCA), non-negative matrix factorization (NMF), K-means clustering, etc., rely on the factorization of a matrix obtained by concatenating high-dimensional vectors from a training collection. While the idealized task would be to optimize the expected quality of the factors over the underlying distribution of training vectors, it is achieved in practice by minimizing an empirical average over the considered collection. The focus of this paper is to provide sample complexity estimates to uniformly control how much the empirical average deviates from the expected cost function. Standard arguments imply that the performance of the empirical predictor also exhibit such guarantees. The level of genericity of the approach encompasses several possible constraints on the factors (tensor product structure, shift-invariance, sparsity . . . ), thus providing a unified perspective on the sample complexity of several widely used matrix factorization schemes. The derived generalization bounds behave proportional to log(n)/n w.r.t. the number of samples n for the considered matrix factorization techniques.

Index Terms-Dictionary learning, sparse coding, principal component analysis, K-means clustering, non-negative matrix factorization, structured learning, sample complexity.

I. INTRODUCTION

T HE fact that a signal x ∈ R m which belongs to a certain class has a representation over some class dependent dictionary D ∈ R m×d is the backbone of many successful signal reconstruction and data analysis algorithms [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]- [START_REF] Elad | Sparse and Redundant Representations, ser. From Theory to Applications in Signal and Image Processing[END_REF].

That is, x is the linear combination of columns of D, referred to as atoms. Formally, this reads as

x = Dα, (1)
where the coefficient vector α ∈ R d as well as the dictionary D are subject to some constraints. Such a setting covers prominent examples like Principal Component Analysis (PCA), where D has orthogonal columns, thus representing the subspace where the signal in the given class is contained.

Another example is the sparse synthesis model, also known as sparse coding, where typically D consists of normalized columns that form an overcomplete basis of the signal space, and α ∈ R d is assumed to be sparse. The task of learning such dictionaries from a given set of training data is related to matrix factorization problems. Important examples include Higher-Order SVD (also known as multilinear SVD) [START_REF] Tucker | Some mathematical notes on three-mode factor analysis[END_REF], sparse coding also called dictionary learning [START_REF] Olshausen | Sparse coding with an overcomplete basis set: A strategy employed by VI?[END_REF]- [START_REF] Rubinstein | Dictionaries for sparse representation modeling[END_REF], its variants with separable [START_REF] Hawe | Separable dictionary learning[END_REF] or sparse [START_REF] Rubinstein | Double sparsity: Learning sparse dictionaries for sparse signal approximation[END_REF] dictionaries, Non-negative Matrix Factorization (NMF) [START_REF] Seung | Algorithms for non-negative matrix factorization[END_REF], K-means clustering [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF], sparse PCA [START_REF] Aspremont | A direct formulation for sparse PCA using semidefinite programming[END_REF]- [START_REF] Zhang | Large-scale sparse principal component analysis with application to text data[END_REF], and more. The learning task is expressed formally as follows.

Let X = [x 1 , . . . , x n ] ∈ R m×n be the matrix containing the n training samples arranged as its columns, and let A = [α 1 , . . . , α n ] ∈ R d×n contain the corresponding n coefficient vectors, a common approach to the dictionary learning process is the optimization problem Therein, g : R d → R is a function that promotes the constraints for the coefficient vectors, e.g., sparsity or positivity, and D is some predefined admissible set of solutions for the dictionary. Note that for PCA there are no constraints on the coefficient vectors, which implies the penalty function g ≡ 0 in this case. A fundamental question in such a learning process is the sample complexity issue. Assuming that the training samples are drawn according to some distribution P representing the class of signals of interest, one would ideally like to select the dictionary D ⋆ yielding the minimum expected value of (2). However, having only access to n training samples, one can at best select an empirical minimizer D. Is this empirical minimizer useful beyond the training set from which it has been selected? This depends on how much the empirical cost function deviates from its expectation.

State of the art sample complexity estimates [START_REF] Maurer | K-dimensional coding schemes in hilbert spaces[END_REF], [START_REF] Vainsencher | The sample complexity of dictionary learning[END_REF] primarily consider the case where g is the indicator function of a set, such as an ℓ 1 or an ℓ 0 ball, D is the set of all unit norm dictionaries or a subset with a restricted isometry property, and the distribution P is in the class P of distributions on the unit sphere of R m . We generalize these results to:

• General classes of penalty functions. Examples covered by our results include: the ℓ 1 penalty and its powers; any mixed norm [START_REF] Jenatton | Structured variable selection with sparsity-inducing norms[END_REF] or quasi-norm [START_REF] Chartrand | Exact reconstruction of sparse signals via nonconvex minimization[END_REF]; the characteristic function of the ℓ 1 -ball, of the set of k-sparse vectors [START_REF] Vainsencher | The sample complexity of dictionary learning[END_REF], or of non-negative vectors [START_REF] Seung | Algorithms for non-negative matrix factorization[END_REF].

• Various classes of dictionaries D that can incorporate structures. Examples covered include: dictionaries with unit norm atoms which are used in many dictionary learning schemes, e.g., K-SVD [START_REF] Aharon | K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF]; sparse dictionaries [START_REF] Rubinstein | Double sparsity: Learning sparse dictionaries for sparse signal approximation[END_REF]; shift-invariant dictionaries [START_REF] Mailhé | Shift-invariant dictionary learning for sparse representations: extending K-SVD[END_REF]; tensor product dictionaries [START_REF] Hawe | Separable dictionary learning[END_REF]; orthogonal dictionaries; non-negative dictionaries [START_REF] Seung | Algorithms for non-negative matrix factorization[END_REF]; topic models [START_REF] Jenatton | Proximal methods for hierarchical sparse coding[END_REF]; and tensor products of Stiefel matrices used for Higher-order SVDs [START_REF] Tucker | Some mathematical notes on three-mode factor analysis[END_REF], [START_REF] Lathauwer | A multilinear singular value decomposition[END_REF].

• Various classes P of probability distributions P. Examples include distributions on the unit sphere which are tackled in [START_REF] Vainsencher | The sample complexity of dictionary learning[END_REF], but also certain distributions built using sub-Gaussian random variables [START_REF] Jenatton | Local stability and robustness of sparse dictionary learning in the presence of noise[END_REF]. For more information on sub-Gaussian random variables, see [START_REF] Buldygin | Metric characterization of random variables and random processes[END_REF]. After formalizing the problem and setting our notations in Section II, we state our main results in Section III and illustrate them with families of examples in Section IV. The obtained sample complexity estimates depend on the Lipschitz constant of the optimized cost function [START_REF] Elad | On the role of sparse and redundant representations in image processing[END_REF], seen as a function of the unknown dictionary D. This Lipschitz constant is essentially driven by the penalty function g, as discussed in Section V. Our results rely on the assumption that the distributions P in the class P satisfy a concentration of measure property, and the structured constraint set D plays a role in the sample complexity estimates through its covering number, as discussed in Section VI where the main theoretical results are established. Covering numbers for a number of structure sets D are considered in Section VII. Section VIII gathers worked examples in relation to previous work. Finally, Section IX discusses sharper sample complexity estimates in the high-dimensional setting m ≫ d.

II. NOTATIONS

Sets are denoted with gothic face as D, P, A. Matrices are written as boldface capital letters like X, D, column vectors are denoted by boldfaced small letters, e.g., α, d whereas scalars are either capital or small letters like n, N . By v i we denote the i th element of the vector v, v ij denotes the i th element in the j th column of a matrix V, while the i th column is referred to by v i . The Frobenius norm for matrices is denoted • F and the corresponding inner product between matrices is denoted •, • F . The operator norm

• 1→2 of an m×d matrix ∆ = [δ 1 , . . . , δ d ] is ∆ 1→2 sup α 1 ≤1 ∆α 2 = max 1≤i≤d δ i 2 .
Finally, log is the natural logarithm so that e log t = t for t > 0.

Given a dictionary D ∈ R m×d that fulfills certain structural properties and a signal x ∈ R m , a representation vector α ∈ R d is typically obtained by solving the minimization problem α ∈ arg min

α∈R d 1 2 x -Dα 2 2 + g(α),
where g : R d → R + ∪ {+∞} is a penalty function promoting constraints for the coefficient vector. For our purposes, the question of whether a minimizer α actually exists is irrelevant, and we define the quality of how well a signal x can be coded by a dictionary D by

f x (D) inf α∈R d L x (D, α), with L x (D, α) 1 2 x -Dα 2 2 + g(α).
Given n training samples X = [x 1 , . . . , x n ] ∈ R m×n , the average quality of how well a dictionary allows a representation for X while considering the constraints imposed on the coefficients is

F X (D) 1 n n i=1 f xi (D).
The cost function F X can be written in short form as

F X (D) = inf A∈R d×n L X (D, A) with L X (D, A) 1 2n X -DA 2 F + G(A)
where A [α 1 , . . . , α n ] and G(A)

1 n n i=1 g(α i ).

III. MAIN RESULTS & OUTLINE OF OUR APPROACH

The main contribution of this paper is a general framework to establish the uniform convergence of the function F X (•) to its expectation when the samples x i are all drawn independently according to an unknown probability distribution P. We show in this work that

sup D∈D |F X (D) -E x∼P f x (D)| ≤ η n (g, D, P).
holds with "overwhelming" probability (that will be controlled explicitly in due time as well as the precision bound η n (g, D, P) ∝ log n/n).

A particular consequence is in terms of generalization bound. An ideal objective in dictionary learning and related matrix factorization problems would be to select the optimal dictionary for the underlying distribution P of training samples,

D ⋆ ∈ arg min D∈D E x∼P f x (D).
In practice, one can at best access an empirically optimal dictionary Dn ∈ D, which minimizes F X (D) given n training samples gathered in X. How close is its performance to that of D ⋆ ? A consequence of our main result is a generalization bound for the empirical optimum Dn : with controlled high probability,

E x∼P f x ( Dn ) ≤ E x∼P f x (D ⋆ ) + 2η n .
Note that the uniform convergence result holds for all dictionaries D ∈ D, and not only at the global optimum of the learning problem. Moreover, even though finding the empirical minimizer Dn often means facing a difficult nonconvex optimization problem, recent work on sparse dictionary learning establishes that certain polynomial time algorithms can provably find it with high probability in certain scenarii [START_REF] Spielman | Exact recovery of sparselyused dictionaries[END_REF]- [START_REF] Agarwal | Learning Sparsely Used Overcomplete Dictionaries via Alternating Minimization[END_REF]. Uniform convergence results as established here should enable further investigation of the finite sample analysis of such algorithms.

The main result, formally stated in Theorem 1 in Section III-E, combines a standard argument based on covering numbers with an estimation of the expected Lipschitz constant of the function D → F X (D). Our main technical contribution is the estimation of this Lipschitz constant, which is summarized in Theorems 2 and 3 in Section V-B. To prove the main result, we first show that D → F X (D) is Lipschitz with a constant expressed in terms of X and g. To be more precise, we prove the property

|F X (D)-F X (D ′ )| ≤ L • D-D ′ 1→2
for all admissible dictionaries D, D ′ ∈ D, where L is dependent on X and g. In the course of our work we provide more tangible bounds for the Lipschitz constant. However, this requires certain assumptions on the penalty function g. These assumptions will be motivated and introduced in the remainder of this section. Therefore, we refrain from presenting the concrete bounds at this point. They can be found in Corollary 3 and Lemma 6. The discussion is concluded using an argument based on covering numbers and concentration of measure. While the latter technique is fairly standard, a major contribution of this paper lies in the identification of two large classes of penalty functions for which the desired Lipschitz property holds and is nicely controlled. In the following, we will give a detailed insight into the required components to provide the generalization bound. Furthermore, we motivate the assumptions that we enforce on the penalty function g (these are labeled with A), or jointly on the penalty function g and the set of admissible dictionaries D (labeled with B), and on the probability distribution (labeled with C).

A. Role of the constraint set D

The structure of the constraint set D is incorporated in the analysis by employing a standard ǫ-net argument. A compact constraint set D ⊂ R m×d can be covered with a finite number of balls of small radius, i.e., it has a finite covering number

N (D, ǫ) min{♯Q : Q ⊂ R m×d , D ⊂ q∈Q B ǫ (q)}
where ♯ denotes the cardinality of a set. For further background about covering numbers, we refer the reader to [START_REF] Cucker | On the mathematical foundations of learning[END_REF] and references therein. In our setting the ǫ-balls B ǫ (•) are defined with respect to the metric • 1→2 since the unit ball in this metric is closely connected with the common constraint that a dictionary has unit norm columns.

We postpone explicit estimates of covering numbers to Section VII, but notice that covering numbers for all the considered classes satisfy the generic upper bound

N (D, ǫ) ≤ (C/ǫ) h , ∀ 0 < ǫ ≤ 1 (3)
where h ≥ 1 is a measure of the "dimension" of D, and C ≥ 1. In particular, for the set of all dictionaries in R m×d with unit-norm columns, which is one of the most common constraint sets, these constants take the values h = md and C = 3. The sample complexity will be essentially driven by the "dimension" h, while the constant C will only play a logarithmic role.

B. Lipschitz property for "norm-like" penalty functions g

In short, we will prove in Section V that under assumptions • A1: g is non-negative;

• A2: g is lower semi-continuous;

• A3: g is coercive: g(α) → ∞ as α → ∞; the function D → F X (D) is Lipschitz over R m×d with controlled constant. With the additional assumption

• A4: g(0) = 0; the Lipschitz constant is bounded by

L X (ḡ) 1 n n i=1 x i 2 • ḡ( x i 2 2 /2) (4)
where we define the auxiliary function ḡ as follows.

Definition 1. For penalty functions g satisfying A1-A3,

ḡ(t) sup α∈R d g(α)≤t α 1 , t ≥ 0.
The presence of the ℓ 1 norm in this definition results from the choice of the metric • 1→2 to measure covering numbers. It will become clearer in the proof of Lemma 5 where we need to control α 1 where α is a (near) minimizer of L x (D, •).

The following properties of this function will be useful.

Lemma 1. The function ḡ as defined in Definition 1 is nonincreasing and takes its values in R + .

Proof of Lemma 1: The sublevel sets of g are nested, hence ḡ is non-decreasing. Since g is coercive, its sublevel sets are bounded, hence ḡ takes its values in R + .

Note that Assumption A4 is a convenience that should not be taken too literally, as penalty functions that fulfill A1-A3 and have a global minimum at 0 can be manipulated to fulfill A4 by subtracting the value of g at 0.

C. Lipschitz property under joint assumptions on g and D

While Assumptions A1-A4 cover a wide range of penalty functions, they do not cover popular penalties related to the ℓ 0 quasi-norm, such as the indicator function (also called characteristic function) of k-sparse vectors:

χ k-sparse (α) 0, if α 0 ≤ k, +∞, otherwise;
Vainsencher et al. [START_REF] Vainsencher | The sample complexity of dictionary learning[END_REF] deal with the latter case under an incoherence assumption on D, i.e., by restricting D to be the class of dictionaries with small (cumulative) coherence.

A careful study of their technique shows that the results are actually valid under an assumption related to the well-known restricted isometry property (RIP) [START_REF] Elad | Sparse and Redundant Representations, ser. From Theory to Applications in Signal and Image Processing[END_REF], [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF]. In fact, while the RIP is usually expressed in its symmetric form, only its lower bound actually plays a role in the considered context. This justifies the following definition.

Definition 2 (Lower RIP). For k ≤ min(m, d), we denote δ k (D) the smallest 0 < δ ≤ 1 such that for all k-sparse vectors α

(1 -δ) α 2 2 ≤ Dα 2 2 .
For k ≤ min(m, d) and 0 < δ < 1 we define the compact set

D LRIP (k, δ) D ∈ R m×d : δ k (D) ≤ δ . (5)
Similarly, Assumptions A1-A4 do not cover the indicator function of non-negative coefficients,

χ NMF (α) 0, if α i ≥ 0, ∀1 ≤ i ≤ d, +∞, otherwise;
which is typically used in conjunction with a non-negativity constraint on D as well, namely D ∈ D NMF , where

D NMF D ∈ R m×d : d ij ≥ 0, ∀i, j . (6)
To unify the treatment of such penalty functions not covered by Assumptions A1-A4, we develop complementary results based on joint assumptions on g and the constraint set D:

• B1: g = χ K is the indicator function of a set K; • B2: there is κ > 0 such that for any α ∈ K and D ∈ D, κ α 2 1 ≤ Dα 2 2
• B3: K contains the origin: 0 ∈ K. Note that B2 is a type of restricted eigenvalue condition, see e.g., [START_REF] Negahban | A unified framework for high-dimensional analysis of m-estimators with decomposable regularizers[END_REF]Equation (30)].

These assumptions hold for g = χ k-sparse and D = D LRIP (k, δ) with κ = (1δ)/d, and they will be shown to also hold for g = χ NMF and D = D NMF (κ) ⊂ D NMF where

D NMF (κ) D ∈ R m×d : d ij ≥ 0, d j 2 2 ≥ κd, ∀i, j . (7)
Under Assumptions B1-B3 we show that the function D → F X (D) is Lipschitz over D. Its Lipschitz constant can be expressed as L X (ḡ) (see Eq. ( 4)) where ḡ is defined in the context of assumptions B1-B3 as: Definition 3. For penalty functions g and D satisfying B1-B3, we define ḡ(t) 2 2t/κ

Occasionally, an additional assumption will be helpful:

• B4: D is convex.

D. Role of the class P of probability distributions

Finally, the results rely on two assumptions on the probability distribution P from which the training samples x i are assumed drawn i.i.d.:

• first, we need to control the Lipschitz constant L X (ḡ) when the sample size n is large. • second, given D, we need to control the concentration of the empirical average F X (D) around its expectation; By the law of large numbers, the first condition holds under assumption

• C1: bounded moment L P (ḡ) E x∼P x i 2 ḡ xi 2 2 2 < +∞.
We will see on many examples (Table I) that tḡ(t 2 /2) ∝ t 2 or tḡ(t 2 /2) ∝ t 3 , hence this is a relatively mild condition. From a quantitative perspective our results will be expressed using

Λ n (L, ḡ) P 1 n n i=1 x i 2 ḡ xi 2 2 2 > L . (8)
By abuse of notation we will simply write Λ n (L) and will exploit the fact that lim n→∞ Λ n (L) = 0 for L > L P (ḡ).

The second condition is measured through

Γ n (γ) sup D∈D P 1 n n i=1 f xi (D) -Ef x (D) > γ . ( 9
)
As discussed in Section VI, our main results exploit C1 and

• C2: there are c > 0 and T ∈ (0, +∞] such that

Γ n (cτ ) ≤ 2 exp(-nτ 2 ), ∀0 ≤ τ ≤ T, ∀n.
As shown in Section VI, this covers the case of probability distributions on the unit sphere in R m (see, e.g., [START_REF] Vainsencher | The sample complexity of dictionary learning[END_REF]),

S m-1 {x ∈ R m : x 2 = 1},
and more generally on Euclidean balls B R of given radius R, Definition 4 (Probability distributions on a ball). The set of probability distribution within a ball of radius R is given by

P BR {P : P( x 2 ≤ R) = 1} . For P ∈ P BR , C2 holds with c = R 2 / √ 8, T = +∞, and C1 holds, with Λ n (L) = 0 for L = Rḡ(R 2 /2) (Lemma 11).
Assumption C2 also covers the following classes which contain the sub-Gaussian sparse signal model proposed in [START_REF] Jenatton | Local stability and robustness of sparse dictionary learning in the presence of noise[END_REF]. More details will be given in Section IV.

Definition 5. A distribution is in P

A , A > 0 if P( x 2 2 ≥ At) ≤ exp(-t), ∀t ≥ 1. (10)
For P ∈ P A , C2 holds with c = 12A, T = 1 (Lemma 12), and C1 holds as soon as ḡ has at most polynomial growth.

E. Main result

Our main result is obtained using a standard union bound argument. The details are in Section VI. In short, under assumptions A1-A4 or B1-B3, together with C1-C2, we show:

Theorem 1. Consider L > L P (ḡ) and define β h • max log 2LC c , 1 , (11)
η n (g, D, P) 3c β log n n + c β+x n . ( 12
)
Then, given 0 ≤ x ≤ nT 2β log n we have: except with probability at most

Λ n (L) + 2e -x , sup D∈D |F X (D) -E x f x (D)| ≤ η n (g, D, P). ( 13
)
Note that Λ n (L) is primarily characterized by the penalty function g and the class of probability distributions P (see [START_REF] Aharon | K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF]), while the constants C, h ≥ 1 depend on the class of dictionaries D, see (3), and c > 0, 0 < T ≤ ∞ depend on the class of probability distributions P, see C2.

The constant 3c in the first term of the right hand side of ( 12) is improved to 2c for penalties that satisfy A1-A4 or B1-B4, see Lemma 10 in Section VI where Theorem 1 is proved.

IV. EXAMPLES

We now discuss concrete settings covered by our results. More detailed worked examples will be given in Section VIII.

A. Penalty functions satisfying A1-A4

Many classical penalty functions are covered by Assumptions A1-A4: norms, quasi-norms, their powers, indicator functions of compact sets containing the origin, and more.

As a side note, if the penalty function g is invariant by sign flips (this is the case, e.g., for standard norms, mixed norms, etc.), one may notice that ḡ(t) is related to the Fenchel conjugate [START_REF] Boyd | Convex Optimization[END_REF] of the function h(α) = χ g(α)≤t , which reads as

h ⋆ (q) sup α ( q, α -h(α)) since ḡ(t) = h ⋆ (1)
where 1 is the vectors with all ones. 1) Norms and quasi-norms: All standard ℓ p norms (1 ≤ p ≤ ∞), quasi-norms (0 < p < 1), and mixed ℓ p,q norms used in structured sparse regularization satisfy A1-A4, and we have ḡ(t) = C g t with

C g = sup α =0 α 1 g(α) < ∞
For example, for g(α) = α p with α ∈ R d , we have

C g = d (1-1/p)+
where (•) + max(•, 0). We let the reader check that for mixed ℓ1,2 norms, α 1,2 = i α Ji 2 , where the sets J i partition the index set 1; d , we have C g = j 1/2 where j = max i ♯J i .

2) Indicator functions of compact sets: Assumptions A1-A3 also hold for the indicator function of a compact set K,

g(α) = χ α∈K 0, if α ∈ K, +∞, otherwise.
For such penalties we have ḡ(t) = sup α∈K α 1 for all t ≥ 0. Assumption A4 further requires that the compact set contains the origin, 0 ∈ K.

In particular, assumptions A1-A4 hold for the indicator function of a ball defined by a (quasi-)norm, i.e., with K = {α : α ≤ λ}. For an ℓ p (quasi-)norm, we have

ḡ(t) = d (1-1/p)+ λ.
3) More examples: Assumptions A1-A4 hold when g is a power of a (quasi-)norm, g(α) = α/λ r p , with r > 0, leading to ḡ(t) = d (1-1/p)+ λt 1/r . Note that the indicator function of K = {α : α ≤ λ} reads as the limit when r → ∞.

There are of course measures that do not fit in our framework, such as the arctangent sparsity measure, which mimics the ℓ 0 one [START_REF] Hage | Robust PCA and subspace tracking from incomplete observations using ℓ 0 -surrogates[END_REF] and is defined as

g atan (α) d j=1 arc tan 2 (α j ) .
It is not coercive and does thereby not meet assumption A3.

B. Penalty functions and constraint sets that satisfy B1-B3

Consider g(α) = χ k-sparse (α) = χ K with the set K = {α : α 0 ≤ k} and D = D LRIP (k, δ). Assumptions B1 and B3 obviously hold since K is a finite union of subspaces. Moreover, by definition of D LRIP (k, δ), for any α ∈ K and D ∈ D we have

Dα 2 2 ≥ (1 -δ) α 2 2 ≥ 1-δ k α 2 1
where the rightmost inequality follows from the fact that α is k-sparse. Hence assumption B2 holds with κ = (1δ)/k. Consider now g(α) = χ NMF (α) = χ K with K = {α : α i ≥ 0, 1 ≤ i ≤ d}, and D = D NMF (κ). Assumptions B1 and B3 obviously hold since K is an orthant. Moreover, since D ∈ D has non-negative entries 1 , for any α ∈ K we have

Dα 2 2 = d i,j=1 α i α j d i , d j ≥ d i=1 α 2 i d i 2 2 ≥ α 2 2 • min 1≤j≤d d j 2 2 ≥ α 2 1 • min 1≤j≤d d j 2 2 /d.
Hence assumption B2 holds. Table I summarizes, for standard penalty functions satisfying either A1-A4 or B1-B3, the expression of ḡ and that of tḡ(t 2 /2) which appears in the expression of L X (ḡ) (see ( 4)).

C. Covering numbers

Table II summarizes the covering numbers of the structured classes of dictionaries considered in Section VII. The provided covering numbers all depend on the signal dimension m, which may lead to sub-optimal sample complexity estimates for high-dimensional problems, i.e., when m ≫ d. This issue is discussed in section Section IX.

D. Probability distributions

While previous work [START_REF] Maurer | K-dimensional coding schemes in hilbert spaces[END_REF], [START_REF] Vainsencher | The sample complexity of dictionary learning[END_REF] only covers distributions on the unit sphere or in the volume of the unit ball, our results cover more complex models, such as the model of sparse signals with sub-Gaussian non-zero coefficients introduced in [START_REF] Jenatton | Local stability and robustness of sparse dictionary learning in the presence of noise[END_REF].

Definition 6 (Sub-Gaussian model of sparse signals). Given a fixed reference dictionary D 0 ∈ D, each noisy sparse signal x ∈ R m is built from the following steps:

1) Support generation: Draw uniformly without replacement k atoms out of the d available in D 0 . This procedure thus defines a support J {j ∈ 1; d : δ(j) = 1} whose size is |J| = k, and where δ(j) denotes the indicator function equal to one if the j-th atom is selected, zero otherwise. Note that E[δ(j)] = k d and for i = j we further have

E[δ(j)δ(i)] = k(k-1)
d(d-1) . 2) Coefficient generation: Define a sparse vector α 0 ∈ R d supported on J whose entries in J are generated i.i.d. according to a sub-Gaussian distribution: for j not in J, [α 0 ] j is set to zero; otherwise, we assume there exists some σ > 0 such that for j ∈ J we have, for all t

∈ TABLE I PENALTY FUNCTIONS g : R d → R AND ASSOCIATED ḡ (FROM EITHER DEFINITION 1 OR DEFINITION 3) g(α) remark ḡ(t) tḡ(t 2 /2) (quasi)norm g(α) Cg sup α =0 α 1 /g(α) Cg t Cg 2 • t 3 α/λ r p , 0 < p ≤ ∞, 0 < r < ∞ d (1-1/p) + • λ • t 1/r d (1-1/p) + • λ • t • (t 2 /2) 1/r χ α p≤λ , 0 < p ≤ ∞ d (1-1/p) + • λ d (1-1/p) + • λ • t χ k-sparse (α) on D LRIP (k, δ), see (5) 2 k/(1 -δ) √ 2t 2t 2 k/(1 -δ) χ NMF (α) on D NMF (κ), see (7) 2 2t/κ 2t 2 / √ κ TABLE II CLASSES OF STRUCTURED DICTIONARIES D WITH COVERING NUMBERS N (D, ǫ) ≤ (C/ǫ) h , 0 < ǫ ≤ 1. SEE SECTION VII D exponent h constant C Unit norm D(m, d) md 3 Separable Dsep i m i d i 3 Sparse Dsparse(m, d) sd 3 m s 1/s Orthogonal D orth (m) m(m -1)/2 3πe π Stiefel D St (m, d) md -d(d + 1)/2 3πe π Stiefel Tensor D ⊗St i m i d i -d i (d i + 1)/2 3πe π R, E{exp(t[α 0 ] j )} ≤ exp(σ 2 t 2 /2
) . We denote σ α the smallest value of σ such that this property holds. For background on sub-Gaussianity, see [START_REF] Buldygin | Metric characterization of random variables and random processes[END_REF], [START_REF] Vershynin | Introduction to the non-asymptotic analysis of random matrices[END_REF]. 3) Noise: Eventually generate the signal x = D 0 α 0 + ǫ, where the entries of the additive noise ǫ ∈ R m are assumed i.i.d. sub-Gaussian with parameter σ ε .

Remark 1. The model in [START_REF] Jenatton | Local stability and robustness of sparse dictionary learning in the presence of noise[END_REF] is more restricted: it assumes that, for j ∈ J, |[α 0 ] j | ≥ α almost surely, where α > 0.

The distribution P of this sub-Gaussian sparse signal model belongs to P A with A = kσ 2 α + mσ 2 ε , as shown below (the argument can be originally found in [START_REF] Jenatton | Local stability and robustness of sparse dictionary learning in the presence of noise[END_REF]).

Lemma 2 (From [START_REF] Hsu | A tail inequality for quadratic forms of subgaussian random vectors[END_REF]). Let us consider z ∈ R m a random vector of independent sub-Gaussian variables with parameters upper bounded by σ > 0. Let A ∈ R m×p be a fixed matrix. For any t ≥ 1, we have

P Az 2 2 > 5σ 2 A 2 F t ≤ exp(-t).
Corollary 1. Let x be a signal following the model of Definition 6. For any t ≥ 1 we have

P x 2 2 > 5(kσ 2 α + mσ 2 ε )t ≤ exp(-t).
Proof: The considered norm can be expressed as follows

x 2 2 = σ α [D 0 ] J σ ε Id 1 σα [α 0 ] J 1 σε ǫ 2 2 .
The result is a direct application of Lemma 2 conditioned to the draw of J, and the observation that

σ α [D 0 ] J σ ε Id 2 F = [D 0 ] J 2 F •σ 2 α +mσ 2 ε = kσ 2 α +mσ 2 ε .
The bound being independent of J, the result is also true without conditioning. Section VIII will detail worked examples with results for specific combinations of dictionary class D, penalty g, and probability distribution P, relating our results to prior work.

V. LIPSCHITZ CONTINUITY OF F X

Under appropriate assumptions on the penalty function g, we prove below that the function F X (D) is Lipschitz continuous with a controlled constant L X (ḡ). We begin by a one-sided Lipschitz property with an additional quadratic term that we will soon get rid of. Lemma 3. Let • be some norm for m × d matrices and

• ⋆ its dual norm 2 . For any ǫ > 0, the set

A ǫ (X, D) = {[α 1 , . . . , α n ] : α i ∈ R d , L xi (D, α i ) ≤ f xi (D) + ǫ}.
is not empty, and for any D ′ we have

F X (D ′ ) ≤ F X (D)+L X (D)• D ′ -D +C X (D)• D ′ -D 2 .
with

L X (D) inf ǫ>0 sup A∈Aǫ 1 n • (X -DA)A ⊤ ⋆ , (14)
C X (D) inf ǫ>0 sup A∈Aǫ C 2n n i=1 α i 2 1 . ( 15
)
The constant C in [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF] depends only on the norm • and the dimensions m, d.

When • = • 1→2 we have C = 1. When A 0 = ∅, L X (D) and C X (D) are suprema over A 0 .
The choice of the ℓ 1 norm in the bound ( 15) is an arbitrary convenience with no impact on the nature of the results to follow. More important will be the choice of the metric • which will be related to how we measure covering numbers for D. This will be discussed later.

Proof: Fix ǫ > 0 and A ∈ A ǫ . From the identity

L X (D ′ , A) = L X (D, A) + 1 n X -DA, (D ′ -D)A F + 1 2n (D ′ -D)A 2 F ( 16
) 2 U ⋆ sup D, D ≤1 U, D F with •, • F the Frobenius inner product.
and the definition of the dual norm we have

| X -DA, (D ′ -D)A F | D ′ -D ≤ (X -DA)A ⊤ ⋆ .
Moreover, by the equivalence of all norms on the finite dimensional space of m × d matrices, there is a constant C (equal to one when

• = • 1→2 ) such that 1 2n (D ′ -D)A 2 F ≤ D ′ -D 2 1→2 • 1 2n n i=1 α i 2 1 ≤ D ′ -D 2 • C 2n n i=1 α i 2 1 .
Therefore, we deduce from ( 16) that

F X (D ′ ) = inf A L X (D ′ , A) ≤ sup A∈Aǫ L X (D ′ , A) ≤ sup A∈Aǫ L X (D, A) + sup A∈Aǫ 1 n (X -DA)A ⊤ ⋆ • D ′ -D + sup A∈Aǫ C 2n n i=1 α i 2 1 • D ′ -D 2 (17)
By definition, sup A∈Aǫ L X (D, A) ≤ F X (D)+ǫ. We conclude by taking the infimum of ( 17) over all ǫ > 0.

Corollary 2. Let • be some norm for m×d matrices and • ⋆ its dual norm, and D a class of dictionaries. If

sup D∈D L X (D) ≤ L X (D); (18)
sup D∈D C X (D) ≤ C X (D), (19)
then for any D = D ′ ∈ D we have

|F X (D ′ ) -F X (D)| D ′ -D ≤ L X (D) • 1 + CX(D) LX(D) D ′ -D .
In particular, if we can establish the bounds ( 18) and ( 19) then F X is uniformly locally Lipschitz (with respect to • ) over the class D for any constant L > L X (D).

A. Uniform local Lipschitz constants

Given Corollary 2 we now bound L X (D) and C X (D) when the norm • is the operator norm • 1→2 . This is motivated by the fact that standard sparse coding is often performed with dictionaries constrained to have unit ℓ 2 norm, which are closely connected with the unit sphere of this norm. With this choice, the shape of ḡ in Definition 1 should appear natural in the context of the proof of Lemma 5 below.

Extensions to other definitions of ḡ can be envisioned if we change the metric. In particular, when g itself is a norm, one could consider the operator norm • g→2 , at the price of possibly more tricky estimates of covering numbers.

1) Penalty functions satisfying A1-A3, with arbitrary D: Lemma 4. If g satisfies Assumptions A1-A3, then the set A 0 is not empty, and it is bounded.

Proof: By the non-negativity (A1) and coercivity (A3) of g(•), L X (D, A) is non-negative and lim k→∞ L X (D, A k ) = ∞ whenever lim k→∞ A k = ∞. Therefore, the function A → L X (D, A) has bounded sublevel sets. Moreover, since g is lower semi-continuous (A2), then so is A → L X (D, A), therefore it attains its infimum value.

Lemma 5 (Penalty functions that satisfy A1-A3). Under Assumptions A1-A3, for any training set X and dictionary D, the constants defined in [START_REF] Seung | Algorithms for non-negative matrix factorization[END_REF] and [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF] with the norm

• = • 1→2 satisfy the bounds L X (D) ≤ 1 n n i=1 2f xi (D) • ḡ(f xi (D)) (20)
C X (D) ≤ 1 2n n i=1 [ḡ(f xi (D))] 2 . ( 21
)
Proof: By Lemma 4 the set A 0 is not empty, and by definition, for any

A [α 1 , . . . , α n ] ∈ A 0 , 1 2 x i -Dα i 2 2 + g(α i ) ≤ f xi (D). (22)
Therefore, A1 (non-negativity of g) implies

0 ≤ g(α i ) ≤ f xi (D) (23)
x i -Dα i 2 ≤ 2f xi (D), (24)
for i = 1, . . . , n. Combined with Definition 1, (23) implies

α i 1 ≤ ḡ(f xi (D)), (25)
C X (D) ≤ 1 2n n i=1 [ḡ(f xi (D))] 2 , (26)
which proves the bound in [START_REF] Jenatton | Structured variable selection with sparsity-inducing norms[END_REF]. We now prove the inequality [START_REF] Vainsencher | The sample complexity of dictionary learning[END_REF]. For any m×d matrix ∆ we have

X -DA, ∆A F = n i=1 x i -Dα i , ∆α i ≤ n i=1 x i -Dα i 2 • ∆α i 2 ≤ n i=1 x i -Dα i 2 • ∆ 1→2 • α i 1 . (27)
This allows us to provide an upper bound for the dual norm by exploiting Equations ( 24) and ( 25)

1 n (X -DA)A ⊤ ⋆ ≤ 1 n n i=1 2f xi (D) • ḡ(f xi (D)) (28)
which shows [START_REF] Vainsencher | The sample complexity of dictionary learning[END_REF].

Corollary 3 (Penalty functions that satisfy A1-A4). Under Assumptions A1-A4, for any training set X and dictionary D, the constants defined in [START_REF] Seung | Algorithms for non-negative matrix factorization[END_REF] and [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF] with the norm

• = • 1→2 satisfy the bounds L X (D) ≤ 1 n n i=1 x i 2 • ḡ( 1 2 x i 2 2 ) L X (ḡ) (29)
C X (D) ≤ 1 2n n i=1 [ḡ( x i 2 2 /2)] 2 C X (ḡ). ( 30
)
Proof: If the condition A4 (g(0) = 0) is fulfilled, we can extend Equations [START_REF] Mailhé | Shift-invariant dictionary learning for sparse representations: extending K-SVD[END_REF][START_REF] Jenatton | Proximal methods for hierarchical sparse coding[END_REF] to

0 ≤ g(α i ) ≤ f xi (D) ≤ L xi (D, 0) = 1 2 x i 2 2
(31)

x i -Dα i 2 ≤ 2f xi (D) ≤ x i 2 , (32)
for i = 1, . . . , n. Due to the fact that ḡ is non-decreasing, the inequality in [START_REF] Agarwal | Learning Sparsely Used Overcomplete Dictionaries via Alternating Minimization[END_REF] follows directly from Equation [START_REF] Jenatton | Local stability and robustness of sparse dictionary learning in the presence of noise[END_REF].

It remains to show the inequality in [START_REF] Arora | New Algorithms for Learning Incoherent and Overcomplete Dictionaries[END_REF], which is verified by plugging Equations [START_REF] Cucker | On the mathematical foundations of learning[END_REF][START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF] in [START_REF] Spielman | Exact recovery of sparselyused dictionaries[END_REF]. This yields the upper bound

1 n n i=1 2f xi (D) • ḡ(f xi (D)) ≤ 1 n n i=1 x i 2 • ḡ( 1 2 x i 2 2 )
which is independent of the constraint set D.

2) Penalty functions and constraint sets that satisfy B1-B3:

Lemma 6 (Penalty functions g and constraint sets D that satisfy B1-B3). Under assumptions B1-B3, for any training set X and dictionary D ∈ D, the constants defined in [START_REF] Seung | Algorithms for non-negative matrix factorization[END_REF] and [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF] with the norm • = • 1→2 satisfy the bounds

L X (D) ≤ 2 n √ κ X 2 F L X (ḡ) (33)
C X (D) ≤ 2 nκ X 2 F C X (ḡ). ( 34
)
Proof: Fix ǫ > 0 and consider α i , i = 1, . . . , n such that A

[α 1 , . . . , α n ] ∈ A ǫ . One can check that by B1, Equations ( 22)-( 24) hold up to an additive O(ǫ) term 3 on the right hand side, as well as [START_REF] Buldygin | Metric characterization of random variables and random processes[END_REF]. By B3 we also have ( 31)-( 32) with a similar additive term, yielding

Dα i 2 ≤ x i 2 + x i -Dα i ≤ 2 x i + O(ǫ). (35)
Moreover, by B2, we have

α i 2 1 ≤ Dαi 2 2 κ ≤ 4 x i 2 2 κ + O(ǫ). (36)
Taking the infimum over ǫ > 0 shows [START_REF] Boyd | Convex Optimization[END_REF].

To prove [START_REF] Negahban | A unified framework for high-dimensional analysis of m-estimators with decomposable regularizers[END_REF] we combine [START_REF] Buldygin | Metric characterization of random variables and random processes[END_REF] with [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF] and [START_REF] Vershynin | Introduction to the non-asymptotic analysis of random matrices[END_REF], yielding an upper bound to 1 n (X -DA)A ⊤ ⋆ :

1 n n i=1 x i -Dα i 2 • α i 1 ≤ 2 n √ κ n i=1 x i 2 2 + O(ǫ). ( 37
)
Observe the similarity between the expressions (33) and ( 4). It justifies Definition 3 which yields the identity L X (D) = L X (ḡ). With this notation we will be able to jointly cover the different penalty functions dealt with so far.

3 f (ǫ) = O(ǫ) if there exists C < ∞ such that |f (ǫ)| ≤ Cǫ.

B. Main results on Lipschitz properties

We are ready to state the main results on the Lipschitz property of F X (D). First we go for a global Lipschitz property. Note that Theorem 2 below is expressed over a convex class of dictionaries: under assumptions A1-A4 the result is valid uniformly over D ∈ R m×d (in particular, it is not restricted to dictionaries with unit norm columns) 4 ; under assumptions B1-B3, we add the explicit assumption B4 that D is convex.

Theorem 2. Assume either A1-A4 or B1-B4. Then, for any X, and any D, D ′ ∈ D,

|F X (D) -F X (D ′ )| ≤ L X (ḡ) • D -D ′ 1→2 . ( 38
)
Proof: Fix ǫ > 0. By Corollaries 2 and 3 (resp. Lemma 6):

|F X (D ′ ) -F X (D)| ≤ (1 + ǫ)L X (ḡ) • D ′ -D 1→2 (39)
whenever δ = D ′ -D 1→2 ≤ ǫLX(ḡ) CX(ḡ) . When δ exceeds this bound, we choose an integer k ≥ 1 such that δ/k ≤ ǫLX(ḡ) CX(ḡ) and define

D i = D + i(D ′ -D)/k, 0 ≤ i ≤ k.
Note that this sequence is a priori constructed in the surrounding space of D. Under assumptions B1-B4, the convexity of D (B4) ensures D i ∈ D. Under assumptions A1-A4, the local Lipschitz constant L X (ḡ) defined in [START_REF] Arora | New Algorithms for Learning Incoherent and Overcomplete Dictionaries[END_REF] actually holds independently of a particular dictionary structure D. Hence, in both cases, since D i+1 -D i 1→2 ≤ ǫLX(ḡ) CX(ḡ) for i = 0, . . . , k -1, the bound we just obtained yields:

|F X (D i+1 ) -F X (D i )| ≤ (1 + ǫ) L X (ḡ) D i+1 -D i 1→2 ≤ (1 + ǫ) L X (ḡ) D ′ -D 1→2 /k |F X (D ′ ) -F X (D)| ≤ k-1 i=0 |F X (D i+1 ) -F X (D i )| ≤ (1 + ǫ) L X (ḡ) D ′ -D 1→2 .
Thus, the bound (39) can be extended to any pair D, D ′ ∈ D.

Since the choice of ǫ > 0 is arbitrary, Equation (38) follows.

As an example, consider NMF expressed with g = χ NMF and D NMF&ℓ 1 the class of non-negative dictionaries with ℓ 1 normalized columns, used in topic models [START_REF] Jenatton | Proximal methods for hierarchical sparse coding[END_REF],

D NMF&ℓ 1 {D : d i ∈ R m ; d i 1 = 1; d ij ≥ 0},
which is convex (B4) as it is the Cartesian product of d copies of the simplex. Since D ⊂ D NMF (κ) with κ = 1/md, assumptions B1-B3 hold.

A slightly weaker result expressed in terms of uniform local Lipschitz property holds under B1-B3 for non-convex D.

|F X (D ′ ) -F X (D)| D ′ -D 1→2 ≤ L X (ḡ) • 1 + 1 κ D ′ -D 1→2 .
Proof: This is the direct result of Lemma 3. The upper bounds for C X (D), L X (D) provided in Lemma 6 yield the factor CX(ḡ) LX(ḡ) = 1 κ .

Remark 2. The reader may wonder why we state a global Lipschitz property in Theorem 2 but only a uniform local one in Theorem 3. First, the triangle inequality argument does not seem to extend: the line joining D to D ′ cannot be cut in small segments so that each intermediate pair remains in D.

In fact, there may even be several connected components in the class. Moreover, even in a given connected component, it is not clear what the length of the geodesics that would mimic such a triangle inequality argument is. Answering such questions would require a deeper investigation of the topology and geometry of D. Since Theorem 3 is sufficient for our needs, we leave them to future work.

C. Lipschitz-continuity for high-dimensional settings

When the signal dimension m is larger than the number d of dictionary elements, then the Lipschitz-continuity studied above is not appropriate as it leads to sample complexities that grow with m. Instead, we may consider Lipschitz-continuity with respect to D ⊤ D and D ⊤ X, which we show below and use in Section IX to sketch improved results in such highdimensional settings.

Lemma 7. Under the same assumptions as in Lemma 3, for any D ′ , we have

F X (D ′ ) ≤ F X (D)+C (1) X (D)• (D ′ ) ⊤ D ′ -D ⊤ D F +C (2) X (D)• X ⊤ D ′ -X ⊤ D F . with C (1) X (D) inf ǫ>0 sup A∈Aǫ 1 2n • AA ⊤ F , C (2)
X (D) inf ǫ>0 sup A∈Aǫ 1 n • A F .
Proof: We follow the same principle as in the proof of Lemma 3, but we use the following equality

L X (D ′ , A) = L X (D, A) + 1 2n (D ′ ) ⊤ D ′ -D ⊤ D, AA ⊤ F -1 n (D ′ ) ⊤ X -D ⊤ X, A F .
Note that we use the Frobenius norm for this lemma, but that we could use any norm on matrices.

D. Other losses

One can envision extensions of Lemmata 3 and 7 and their consequences to matrix factorization problems where the ℓ 2 data-fidelity term, 1 2 x-Dα 2 2 , is replaced by a more general loss ℓ(x, Dα). In particular, Lemma 3 can be extended to losses ℓ(•, •) satisfying for any x, y ℓ(x, y + h) ≤ ℓ(x, y) + ∇ y ℓ(x, y), h + C h 2 2 , i.e., with a local quadratic behavior. This would be reminiscent of the work of Negahban et al [START_REF] Negahban | A unified framework for high-dimensional analysis of m-estimators with decomposable regularizers[END_REF] which covers M-estimation problems where the loss, which is convex in the unknown parameter is regularized with a convex decomposable penalty.

Here the considered problem is intrinsically non convex in D.

A generalization of Corollary 3 would further require assuming that ℓ(x, y) ≥ 0 for any x, y, to replace [START_REF] Cucker | On the mathematical foundations of learning[END_REF] with 0 ≤ g(α i ) ≤ ℓ(x, 0), and identifying conditions on the loss ensuring that for some function B(.), an analogue of the bound (32) holds:

∇ y ℓ(x i , Dα i ) 2 ≤ B(x i ). The resulting Lipschitz constant would read L X = 1 n n i=1 B(x i )ḡ(ℓ(x i , 0)).
The full characterization of such extensions and of the families of losses that can be considered is, however, beyond the scope of this paper.

VI. SAMPLE COMPLEXITY (PROOF OF THEOREM 1)

Given the global (resp. uniform local) Lipschitz property of F X , a standard route (see, e.g., [START_REF] Maurer | K-dimensional coding schemes in hilbert spaces[END_REF], [START_REF] Vainsencher | The sample complexity of dictionary learning[END_REF]) to control the sample complexity via a uniform convergence result is to rely on concentration of measure and covering numbers for the considered class of dictionaries (see Section VII for details on covering numbers).

A. Lipschitz property of the expected cost function

Mild assumptions on the quantities Λ n (L) and Γ n (γ) as defined in ( 8) and ( 9) are sufficient to control the Lipschitz constants of F X and of its expectation: Lemma 8. Under assumptions A1-A4 or B1-B4, and C1:

1) the function D → F X (D) is Lipschitz with constant L , with probability at least 1 -Λ n (L). 2) the expected cost function D → Ef x (D) is Lipschitz with constant L as soon as L > L P (ḡ) and: C2' there exists a sequence γ n such that

lim n→∞ γ n = 0; lim n→∞ Γ n (γ n ) = 0
Proof: The first result trivially follows from Section V and the definition of Λ n (L). For the second one, given D, D ′ , consider an i.i.d. draw X of n samples from P. We have

|Ef x (D ′ ) -Ef x (D)| ≤ |Ef x (D ′ ) -F X (D ′ )| + |F X (D ′ ) -F X (D)| + |F X (D) -Ef x (D)| ≤ L D ′ -D 1→2 + 2γ n
except with probability at most Λ n (L) + 2Γ n (γ n ). The limit for large n yields the desired bound with probability one. Remark 3. Under B1-B3 only (without the convexity assumption on D), an analogon of Lemma 8 holds where the conclusion that D → F X (D) (resp. its expectation) is "Lipschitz with constant L" is replaced with:

|F X (D ′ ) -F X (D)| D -D ′ 1→2 ≤ L • 1 + 1 κ D -D ′ 1→2 . (40)
(resp. [START_REF] Blanchard | Statistical properties of kernel principal component analysis[END_REF] holds with F X (•) replaced by Ef x (•)).

We are now ready to state a first uniform convergence result.

B. Abstract uniform convergence result

Lemma 9. Assume C1-C2. For any ǫ, γ > 0 we have, except with probability at most

Λ n (L) + N (D, ǫ) • Γ n (γ) (41)
• under A1-A4 or B1-B4:

sup D∈D |F X (D) -Ef x (D)| ≤ 2Lǫ + γ (42)
• under B1-B3 only:

sup D∈D |F X (D) -Ef x (D)| ≤ 2Lǫ • 1 + 1 κ ǫ + γ (43)
Proof: We only give the proof under A1-A4 or B1-B4. A straightforward adaptation yields the proof under B1-B3 only.

Fix ǫ > 0 and consider an ǫ-cover of D with respect to the • 1→2 metric with at most N = N (D, ǫ) elements {D j }. Fixing D ∈ D there is an index j such that D j -D 1→2 ≤ ǫ, hence

|F X (D) -Ef x (D)| ≤ |F X (D) -F X (D j )| + |F X (D j ) -Ef x (D j )| + |Ef x (D j ) -Ef x (D)| ≤ |F X (D) -F X (D j )| + sup 1≤j≤N |F X (D j ) -Ef x (D j )| + Lǫ
where we used the Lipschitz property of the expected cost function (Lemma 8). By a union bound, for any γ > 0 we obtain the bound [START_REF] Linder | On the training distortion of vector quantizers[END_REF] except with probability at most [START_REF] Pollard | A central limit theorem for k-means clustering[END_REF].

C. Main uniform convergence result

The above sample complexity estimate is not quite explicit since it combines covering numbers and concentration bounds. In the following, we give more explicit sample complexity estimates under more specific assumptions on N (D, ǫ) and Γ n (γ). The assumption on N (D, ǫ) is justified by the fact that all examples of classes D developed in Section VII have covering numbers with bounds of the type expressed in [START_REF] Elad | Sparse and Redundant Representations, ser. From Theory to Applications in Signal and Image Processing[END_REF].

The assumption C2 on Γ n (γ) is further discussed at the end of this section.

Lemma 10. Assume C1-C2 holds and that N (D, ǫ) satisfy the bound [START_REF] Linder | Rates of convergence in the source coding theorem, in empirical quantizer design, and in universal lossy source coding[END_REF], i.e., for c > 0, T ∈ (0, ∞], and C, h ≥ 1:

Γ n (cτ ) ≤ 2 exp(-nτ 2 ), ∀0 ≤ τ ≤ T ; (44)
N (D, ǫ) ≤ C ǫ h , ∀0 < ǫ ≤ 1. ( 45
) Define β h • max(log 2LC c , 1) and D = 1 (under A1-A4 or B1-B4) or D = max(1/κ, 1) (under B1-B3). Assume that the sample size n satisfies n log n ≥ max 8, β T 2 , D • c 2L 2 β . ( 46
)
Then, for any

0 ≤ x ≤ nT 2 -β log n ( 47
)
we have, except with probability at most

Λ n (L) + 2e -x sup D∈D |F X (D) -Ef x (D)| ≤ η n (L, C, h, c).
In the case of A1-A4 or B1-B4, we have

η n 2c β log n n + c β+x n ,
whereas for B1-B3, we obtain η n 3c β log n n + c β+x n . Note that Theorem 1 follows from Lemma 10.

Proof: First we observe that the condition C2 on Γ n implies C2' (defined in Lemma 8), and since we assume C1 we can apply Lemma 9. Notice that [START_REF] Biau | On the performance of clustering in hilbert spaces[END_REF] implies that one can indeed find x that satisfies assumption [START_REF] Sridharan | Fast Rates for Regularized Objectives[END_REF]. We set

ǫ = c √ β 2L log n n , τ = h[log C/ǫ]+x n = h[log(2LC/c √ β)]+ h 2 log n log n +x n . Since (46) implies n ≥ 8, hence log n 2 ≥ 1, we have 0 ≤ x ≤ nT 2 -β 2 log n -β 2 log n ≤ nT 2 -h 2 log n -β.

By definition of β and h we have

β ≥ h ≥ 1 hence nτ 2 -h 2 log n log n = h log 2LC c √ β + x ≤ h log 2LC c + x ≤ β + x ≤ nT 2 -h 2 log n log n
. This shows that 0 ≤ τ ≤ T . Moreover, by [START_REF] Biau | On the performance of clustering in hilbert spaces[END_REF] we further have 0 < ǫ ≤ 1, hence we can apply ( 44) and ( 45) to obtain

N (D, ǫ) • Γ n (cτ ) ≤ 2 • (C/ǫ) h • exp(-nτ 2 ) = 2 • e -x .
Using A1-A4 or B1-B4 we conclude as follows: since β ≥ 1 and log n ≥ 1 we have

2Lǫ + cτ = c β log n n + c • h log 2LC c √ β + h 2 log n log n + x • 1 √ n ≤ c β log n n + c • β 2 log n + β + x • 1 √ n ≤ c β log n n • (1 + 1 √ 2 ) + c • β + x • 1 n ≤ 2c β log n n + c β+x n .
Under B1-B3, the definition of D and assumption (46) imply 0 < 1 κ ǫ ≤ 1 hence we get similarly:

2Lǫ • 1 + 1 κ ǫ + cτ ≤ 3c β log n n + c β+x n .

D. On assumptions C1 and C2

Assumptions C1-C2 are actually satisfied under rather standard hypotheses: Lemma 11. Assume that P ∈ P RS m-1 . Then

Γ n (R 2 τ / √ 8) ≤ 2 exp(-nτ 2 ), ∀n, ∀τ ≥ 0 (48) Λ n (R ḡ(R 2 /2)) = 0, ∀n.
In other words, C2 holds with c = R 2 / √ 8 and T = +∞, and C1 holds with L P (ḡ) ≤ Rḡ(R 2 /2)).

Proof: For any D, the random variables

y i = f xi (D) satisfy 0 ≤ y i ≤ 1 2 x 2 2 ≤ R 2 2
almost surely. Applying Hoeffding's inequality yields [START_REF] Levrard | Fast rates for empirical vector quantization[END_REF]. We conclude by observing that since x i 2 ≤ R almost surely, we have L X (ḡ) ≤ R ḡ(R 2 /2) almost surely. Lemma 12. Assume there is a constant A > 0 such that P ∈ P A and that the penalty function g satisfies A1 (nonnegativity) and A4 (g(0) = 0). Then

Γ n (12Aτ ) ≤ 2 exp(-nτ 2 ), ∀0 ≤ τ ≤ 1.
In other words, C2 holds with c = 12A and T = 1.

Moreover, in all examples of penalties g considered in Table I, the growth of tḡ(t 2 /2) for large t is at most polynomial, so that L P (ḡ) < +∞, and C1 holds.

Proof: We follow the argument in [START_REF] Jenatton | Local stability and robustness of sparse dictionary learning in the presence of noise[END_REF]. This primarily results from Bernstein's inequality, see, e.g., [START_REF] Bennett | Probability inequalities for the sum of independent random variables[END_REF].

Lemma 13 (Bernstein's Inequality). Let {z j } j∈ 1;n be independent, zero-mean random variables. If there exist v, M > 0 such that for any integer q ≥ 2 and any j ∈ 1; n , it holds

E[|z j | q ] ≤ q! 2 M q-2 v 2 ,
then we have for any γ ≥ 0,

P 1 n n j=1 z j > γ ≤ exp - nγ 2 2(v 2 + M γ) .
In particular, for any τ ≤ v 2M , we have

P 1 n n j=1 z j > 2vτ ≤ exp -nτ 2 .
We can now state a simplification of [START_REF] Jenatton | Local stability and robustness of sparse dictionary learning in the presence of noise[END_REF]Lemma 24].

Lemma 14. Let y be a random variable and assume there is a constant B > 0 such that, for any t ≥ 1,

P (|y| > Bt) ≤ exp(-t). ( 49
)
Then, for any u ≥ 1, any integer q ≥ 1, and 0 < p ≤ 1, we have

E |y| pq ≤ q! B p u q [1 + exp(3 -u)] (50)
E |y| p -E |y| p q ≤ q! 2B p u q [1 + exp(3 -u)]. ( 51
)
To keep the flow of the paper we postpone the proof to the appendix. We now have the tools to state a simplified version of [START_REF] Jenatton | Local stability and robustness of sparse dictionary learning in the presence of noise[END_REF]Corollary 6] which suits our needs. Corollary 4. Consider n independent draws {y i } i∈ 1;n satisfying the hypothesis [START_REF]Non Asymptotic Bounds for Vector Quantization in Hilbert Spaces[END_REF] for some B > 0. Then, for any 0 < p ≤ 1 and 0 ≤ τ ≤ 1, we have

E{|y i | p } ≤ 6B p (52) P 1 n n i=1 (|y i | p -E |y i | p ) ≥ 24B p • τ ≤ 2 exp(-nτ 2 ) ( 53
)
Proof: We apply Lemma 14 with u = 3. For q = 1, Equation ( 50) yields [START_REF] Schnass | On the identifiability of overcomplete dictionaries via the minimisation principle underlying K-SVD[END_REF]. For q ≥ 2, Equation ( 14) shows that we can apply Bernstein's inequality (Lemma 13) with

z i = |y i | p -E |y i | p , M = 2B p u = 6B p and v = √ 2M √
1 + e 3-u = 2M = 12B p . This shows that for 0 ≤ τ ≤ v/2M = 1 we have [START_REF] Gautschi | The incomplete Gamma functions since Tricomi[END_REF].

We are now ready to prove Lemma 12. First, we notice that for any D, 0 ≤ y i = f x (D) ≤ 1 2 x 2 2 , hence we have

sup D∈D P f x (D) ≥ A 2 t ≤ exp(-t), ∀t ≥ 1.
Applying Corollary 4 with B = A/2 and p = 1 yields Γ n (12Aτ ) ≤ 2 exp(-nγ 2 ) for 0 ≤ τ ≤ 1.

VII. CONSTRAINT SET STRUCTURES AND COVERING

NUMBERS

Dictionaries learned with standard methods are unstructured matrices D ∈ R m×d that allow factored representations of the signals of interest. However, in practice the dimension of the signals which are being represented and consequently the possible dictionaries' dimensions are inherently restricted by limited memory and limited computational resources. Furthermore, when used within signal reconstruction algorithms where many matrix vector multiplications have to be performed, those dictionaries are computationally expensive to apply. In the following we first describe various classes of dictionaries and their motivation in signal processing. Then we turn to the estimation of their covering numbers.

A. Unstructured Dictionaries with Unit Norm Atoms

Imposing unit norm on the atoms is a common technique to avoid trivial solutions for the learning problem. An exemplary learning algorithm that employs this kind of restriction is the famous K-SVD algorithm proposed in [START_REF] Aharon | K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF]. All following dictionary structures base upon this restriction. The class of dictionaries of this form is a product of spheres defined as

D(m, d) {D ∈ R m×d : ddiag(D ⊤ D) = I d }.
The operator ddiag maps matrices R d×d to R d×d by leaving the diagonal entries unaltered while setting all others to zero.

In NMF, this dictionary structure is further limited by requiring non-negative entries, yielding the class

D NMF ∩ D(m, d), or D NMF&ℓ 1 if the ℓ 1 normalization is preferred.
Lemma 15 (Covering number bound for D(m, d)). For the Euclidean metric, S m-1 the unit sphere in R m and for any ǫ > 0, we have

N (S m-1 , ǫ) ≤ 1 + 2/ǫ m .
Moreover, for the metric induced on R m×d by • 1→2 , and for any ǫ > 0, we have

N (D(m, d), ǫ) ≤ 1 + 2/ǫ md .
Proof: Lemma 2 in [START_REF] Vershynin | Introduction to the non-asymptotic analysis of random matrices[END_REF] gives the first conclusion for the sphere in R m . As for the second result, remember that the set D(m, d) is a Cartesian product of d (m-1)-spheres.

B. Sparse Dictionaries

The authors of [START_REF] Rubinstein | Double sparsity: Learning sparse dictionaries for sparse signal approximation[END_REF] propose a technique to learn so-called sparse dictionaries. A sparse dictionary is the representation of a learned dictionary D and a given base dictionary Θ (e.g., DCT, ODCT, Wavelet) with the sparse matrix D via D = Θ D. This learning algorithm offers a reduction in training complexity while aiming at maintaining the data modeling abilities of unconstrained dictionary learning. We show in the appendix that:

N (D sparse (m, d, s, Θ), ǫ) ≤ m s (1 + 2/ǫ) s d .

C. Orthogonal Dictionaries

Orthogonal dictionaries are square matrices belonging to

D orth (d) {D ∈ R d×d : D ⊤ D = I d , det(D) = 1}. ( 54
)
They represent the simplest form of dictionaries that allow a unique representation of the signal. Computing the inner product of the signal and the atoms of the dictionary provides the representation coefficients. We show in the appendix that:

N (D orth (d), ǫ) ≤ (3πe π /ǫ) d(d-1)/2 .

D. Stiefel manifold

In PCA (see Section VIII), the dictionary is a matrix consisting of d orthonormal columns where d ≤ m, i.e., D is a tall matrix. This constraint set is the so called Stiefel manifold

D St (m, d) {D ∈ R m×d : D ⊤ D = I d }. (55)
We show in the appendix that:

N (D St (d, m), ǫ) ≤ (3πe π /ǫ) md-d(d+1)/2 .

E. Separable Dictionaries

In [START_REF] Hawe | Separable dictionary learning[END_REF] a dictionary learning scheme is developed that a separable dictionary for sparse image representations: instead of sparsely representing x = vec(Y) as x = Dα, where Y is a matrix containing an image, one directly represents Y sparsely as B via two dictionaries D 1 , D 2 by the equation

Y = D 1 BD ⊤ 2 . (56)
Note that in order to relate this method to standard dictionary learning techniques Equation (56) can be rewritten as

vec(Y) = (D 2 ⊗ D 1 ) • vec(B).
We show in the appendix that:

N (D sep , ǫ) ≤ (3/ǫ) i midi .

F. Other Tensor Dictionaries

Similar to the structure of Separable dictionaries we can define other constraint sets that are combinations of dictionary structures. Examples for this are tensor products of Stiefel matrices defined as

D ⊗St D St (m 1 , d 1 ) ⊗ . . . ⊗ D St (m z , d z ) (57)
which are used for Higher-order SVD. For the sake of readability the dimensions are omitted in this definition. We show in the appendix that:

N (D ⊗St , ǫ) ≤ (3πe π /ǫ) i midi-di(di+1)/2 .

VIII. WORKED EXAMPLES

In this section we propose some settings that fit to our framework and provide the corresponding sample complexity estimates, i.e., the value of β.

Example 1 (PCA). PCA can be expressed as the optimization problem min

A;D∈DSt(m,d) 1 2 X -DA 2 F + 1 n g(α i )
with the penalty function g = χ d-sparse . Since α i ∈ R d , the penalty function is equivalent to 0. Its sole purpose is to shape the problem in a way that fits into the general framework [START_REF] Elad | On the role of sparse and redundant representations in image processing[END_REF]. The constraint set is D St (m, d), the set of m×d-matrices with orthonormal columns. a) Training data in the unit ball: Initially, we consider a distribution P in the unit ball

(R = 1, therefore, c = 1/ √ 8, cf. Lemma 11). Since D St (m, d) ⊂ D LRIP (k, δ) with k = d and δ = 0, the Lipschitz constant is L = 2 √ d, cf. Table I. Table II yields h = md -d(d + 1
)/2 and C = 3πe π . Hence, the constant driving the sample complexity is

β PCA = md -d(d+1) 2 • log(12πe π √ 8d).
Compared to existing sample complexity results for PCA [START_REF] Shawe-Taylor | On the eigenspectrum of the gram matrix and the generalization error of kernel-PCA[END_REF], [START_REF] Blanchard | Statistical properties of kernel principal component analysis[END_REF], this seems rather pessimistic: in [START_REF] Blanchard | Statistical properties of kernel principal component analysis[END_REF] the sample complexity is driven by β ∝ d, and we lose a large dimensiondependent m factor 5 as well as an additional mild log n factor. In fact, the assumption that the data is distributed in the unit ball is crucial in the results of [START_REF] Shawe-Taylor | On the eigenspectrum of the gram matrix and the generalization error of kernel-PCA[END_REF], [START_REF] Blanchard | Statistical properties of kernel principal component analysis[END_REF] as it allows the use of the McDiarmid concentration inequality, while, in our general context, we exploit a standard but somewhat crude argument using union bounds and covering numbers. This however means we can handle data distributions not restricted to live on the unit ball, as we will see next.

b) Sub-Gaussian training data:

The sub-optimality of our results for PCA with respect to state of the art is probably the price to pay for the generality of the approach. For example, in contrast to previous work, our approach provides results for other distributions. For example, given the distribution P ∈ P A , we obtain the constant c = 12A (Lemma 12), while all other constants remain the same. Thus, we get

β PCA = md -d(d+1) 2 • max log πe π √ d A , 1 , (58)
Example 2 (Sparse coding, training data in the unit ball). Consider P a distribution in the unit ball (R = 1). We have c = 1/ √ 8 (Lemma 11), and Table I immediately provides values of the Lipschitz constant L = ḡ(1/2) such that Λ n (L) = 0. For example, for 0 < p < ∞

• for g χp (α) = χ α p ≤λ , we have L p = d (1-1/p)+ λ; II provides h = md and C = 3, and Theorem 1 yields

• for g χ0 (α) = χ α 0≤k , we have L 0 = 2 k 1-δ For D = D(m, d) (resp. D = D LRIP (k, δ)), Table
β χp = md • max log 6 √ 8λd (1-1/p)+ , 1 , β χ0 = md • log 12 8k/(1 -δ).
5 See Section IX for techniques to handle this high-dimensional scaling.

The cases g χ1 and g χ0 are dealt with by Vainsencher et al. [START_REF] Vainsencher | The sample complexity of dictionary learning[END_REF]Theorem 7] who provide a comparable result, however restricted to distributions on the unit sphere whereas our framework allows for more complex distributions. Moreover while they make assumptions on the cumulative coherence of D for their results with the ℓ 0 constraint, we only rely on the restricted isometry property.

The bound of Maurer and Pontil [START_REF] Maurer | K-dimensional coding schemes in hilbert spaces[END_REF] for g χp , p ≥ 1, λ = 1, is of the order β ′ p ∝ d 4-2/p . We extend it to p < 1 as well, and improve the case p = 1 as soon as m ≤ d.

In addition to the above penalty functions associated to the indicator functions of the ℓ 1 (resp. ℓ 0 ) ball, our results allow dealing with the ℓ 1 norm itself -a very popular proxy for the ℓ 0 "norm"-or its powers as penalty functions. For D = D(m, d) (resp. D = D LRIP (k, δ)), the sparsity measures

• g ℓ 1 (α) = α 1 /λ, yielding L = λ 2 ; • g (ℓ 1 ) 2 (α) = α 2 1 /λ 2 , yielding L = λ √ 2
; are associated to the sample complexities

β ℓ 1 = md • max log 3 √ 8λ, 1 β (ℓ 1 ) 2 = md • max (log 12λ, 1)
for probability distributions in the unit ball.

Example 3 (Sparse coding, sub-Gaussian training data). For a distribution P satisfying assumption [START_REF] Tošić | Dictionary learning[END_REF] we have c = 12A. for the product of spheres. This yields the constants C = 3 and h = md, and together with the constant c = 1/ √ 8, which we get for distributions on the unit sphere, we obtain

β NMF = md • log 12 √ 8md. (60)
In a similar setting Maurer and Pontil [START_REF] Maurer | K-dimensional coding schemes in hilbert spaces[END_REF] obtain a sample complexity estimate β ∝ d 3 which has the advantage of being independent of the signal dimension m. Yet, our result seems sharper unless the signal dimension m exceeds d 2 , in which case the approach sketched in Section IX could be followed to recover improved estimates in this high-dimensional scaling.

Example 6 (K-means clustering, training data in the unit ball).

Given a set of n data points K-means clustering algorithms learn a dictionary with K columns and assign each data point to one of these columns while minimizing the sum of squared distances. These atoms represent the K distinct classes. In our notation so far d = K. Then the problem of K-means clustering can be expressed as the minimization problem

min A;D∈R m×K 1 2 X -DA 2 F + n i=1 g(α i ) (61)
with the penalty function

g(α) = χ α 0 =1 + χ k α k =1 . (62)
This penalty function ensures that all columns of A have exactly one entry with the value 1. The centers of the clusters are represented in the columns of D. This penalty function fulfills A1-A3, and its auxiliary function is ḡ(t) ≡ 1 for all t ∈ R. However, g does not fulfill A4. This has the consequence that we cannot merely apply Corollary 3. Instead, we have to rely on the preceding Lemma 5 which leaves us with the upper bound x -Dα 2 2 + g(α) to be finite, the coefficient vector α has to fulfill g(α) = 0. Due to the construction of g, there is only a finite number of possible choices for α (namely α is a standard basis vector in R d ), and the Lipschitz constant can be upper bounded by

L X (D) ≤ 1 n n i=1 max j=1,...,K x i -d j 2 ≤ max i,j x i -d j 2 .
Restricting our analysis to the case of data samples x i lying within the unit ball ( x i 2 ≤ 1), the optimum dictionary coincides with the optimum where the matrix D is constrained to belong to D ball (m, K). It follows that the Lipschitz constant is simply bounded by L = 2. Moreover, D ball (m, K) has the same covering number as D(m, K). From the structural constraints to the dictionary and the signal distribution we obtain the constants L = 2, C = 3, h = mK, c = 1/ √ 8. The above implies the constant

β K-means = mK • log(12 √ 8).
The literature on the sample complexity of K-means is abundant. Following the early work of Pollard [START_REF] Pollard | A central limit theorem for k-means clustering[END_REF], a series of authors have established worst-case lower bounds on the excess risk decaying at best as 1/

√ n [42]- [START_REF] Antos | Improved minimax bounds on the test and training distortion of empirically designed vector quantizers[END_REF], as well as upper bounds of the order of O(1/ √ n) (for a given dimension m) [START_REF] Linder | Rates of convergence in the source coding theorem, in empirical quantizer design, and in universal lossy source coding[END_REF]. More concrete estimates, where the role of K and m is explicated, are needed to compare with our sample complexity estimate. Bartlett et al. [START_REF] Bartlett | The minimax distortion redundancy in empirical quantizer design[END_REF] bound the excess risk by a constant times min Km/n, K 1-2/m m log n/n which, for large m and arbitrary K, is essentially matched by our estimate. Estimates independent of the dimension m have been obtained by Biau et al. [START_REF] Biau | On the performance of clustering in hilbert spaces[END_REF] (see also [START_REF] Maurer | K-dimensional coding schemes in hilbert spaces[END_REF]) of the order K/ √ n. These correspond to β ∝ K 2 , hence our result seems sharper when m K, i.e., when the number of clusters exceeds the dimension.

Example 7 (Higher-order SVD). A natural extension of PCA is the so called Higher-order SVD (HOSVD) or multilinear SVD introduced in [START_REF] Lathauwer | A multilinear singular value decomposition[END_REF]. Let the training signals x i be real ztensors, i.e., X is an element of R m1×m2×...×mz ×n with X = [X 1 . . . X n ] where X i ∈ R m1×m2×...×mz . The HOSVD X is obtained by solving the minimization problem min A,D∈D⊗St

1 2 n i=1 X i -D • A i 2 F .
In this equation D is an element of the constraint set D ⊗St as defined in (57) and therefore a concatenation of Stiefel matrices D 1 , . . . , D z with D i ∈ D St (m i , d i ), and the coefficient tensor A ∈ R d1×...dz×n has the form A = [A 1 . . . A n ] with A i ∈ R d1×...dz . The operator • denotes the operation

D • A i A i × 1 D 1 × 2 D 2 . . . × z D z ,
where × N denotes the so-called N -mode product introduced in [START_REF] Lathauwer | A multilinear singular value decomposition[END_REF]. Just as for PCA the trivial penalty G(A) ≡ 0 can be written using g(α) = χ α 0 ≤ z i=1 di where α i is the vector resulting from the unfolding of the i-th component of A, 1 ≤ i ≤ n, along its rightmost dimension.

Table II yields the constants C = 3πe π and h = z i=1 m i d id i (d i + 1)/2. To obtain the Lipschitz constant L we unfold the tensor A to the matrix R ( z i=1 di)×n and then use a penalty definition similar to the one used for PCA. This yields the constant L = 2 • z i=1 d i . For distributions P in the ball with R = 1, Theorem 1 provides

β HOSVD = z i=1 m i d i -di(di+1) 2 • log   12πe π 8 z i=1 d i   .
In light of state of the art results [START_REF] Blanchard | Statistical properties of kernel principal component analysis[END_REF] on PCA, this apparently new sample complexity estimate is probably quite pessimistic.

IX. HIGH-DIMENSIONAL SETTINGS

When m gets larger, the reasoning in Section VI based on covering numbers of the set of admissible dictionaries D leads to unfavorable scaling with respect to m. This factor is introduced through the ǫ-net argument following Lemma 9. In order to avoid the issues that arise in the high-dimensional setting, we could envision using an alternative to Lemma 3 which served as a precursor to the discussion of the Lipschitz property. Such an alternative is proposed in Lemma 7, where we derive a Lipschitz property of F X w.r.t. the Frobenius norms of D ⊤ D and D ⊤ X which are in R d×d and R d×n , respectively. Covering net arguments revolving around these quantities clearly become independent of the signal dimension m, but would yield difficulties due to their dependence on the draw of X. Providing concrete results for this cases is beyond the scope of this paper, instead, we sketch below an alternative approach. Another option to handle the high-dimensional setting is to extend the approach proposed in [START_REF] Maurer | K-dimensional coding schemes in hilbert spaces[END_REF]. We provide a rough outline below. First, recall that the initial problem is to provide an upper bound for the expression

sup D∈D 1 n n i=1 f xi (D) -E x f x (D) .
We can use McDiarmid's inequality in combination with a symmetrization argument and Rademacher averages to provide the upper bound

sup D∈D 1 n n i=1 f xi (D) -E x f x (D) =C 1 • E η sup D∈D A D (η) + C 2 • 2t/n
with probability at least 1e -t and the Gaussian process

A D (η) = 1 n n i=1 η i f xi (D)
where the variables η i are i.i.d. standard Gaussian. Therefore, it remains to find an upper bound for E η sup D∈D A D (η). Slepian's Lemma states that for any Gaussian process B D that fulfills the condition

E η |A D (η) -A D ′ (η)| 2 ≤ E ξ |B D (ξ) -B D ′ (ξ)| 2 the bound E η [sup D∈D A D (η)] ≤ E ξ [sup D∈D B D (ξ)] holds. In the case at hand we can observe that E η |A D (η) -A D ′ (η)| 2 = 1 n 2 n i=1 |f xi (D) -f xi (D ′ )| 2 . Lemma 7 yields the upper bound |f x (D) -f x (D ′ )| ≤C (1) X D ⊤ D -(D ′ ) ⊤ D ′ F + C (2) X (D -D ′ ) ⊤ x 2 , (63)
which implies that the Gaussian process

B D (ξ) = C ′ D ⊤ D, ξ (1) F + C ′′ D ⊤ x, ξ (2) 2 (64)
fulfills the condition of Slepian's Lemma for appropriate choices of the scalars C ′ , C ′′ . Here, the first scalar product operates in R d×d , whereas the second one is defined on R d . Thus, computing the supremum over all D would result in an expression which is independent of the signal dimension m. Especially in the case of PCA, the final sample complexity result would benefit since the first term in Equation ( 63), and thereby the first term in (64), vanishes, yielding a tighter bound. However, it is not within the scope of this work to provide concrete results for this alternative approach.

X. CONCLUSION AND DISCUSSION

We proposed a general framework to determine the sample complexity of dictionary learning and related matrix factorization problems. The generality of the framework makes it applicable for a variety of structure constraints, penalty functions, and signal distributions beyond previous work. In particular, it covers formulations such as principal component analysis, sparse dictionary learning, non-negative matrix factorization, or K-means clustering, for which we provide sample complexity bounds in the worked examples section.

To keep the exposition of our results simple, we focused on matrix factorization where the data fidelity term is expressed with squared ℓ 2 norm. A straightforward adaptation of the computation of the Lipschitz constant of Section V can easily be made to achieve similar results with the plain ℓ 2 norm [START_REF] Vainsencher | The sample complexity of dictionary learning[END_REF], and easy adaptations can be envisioned when the data fidelity term is ℓ(x-Dα) where the gradient ∇ℓ(•) is Lipschitz. More general loss functions are briefly discussed in Section V-D.

The obtained sample complexity results applied to sparse coding extend those of Maurer and Pontil [START_REF] Maurer | K-dimensional coding schemes in hilbert spaces[END_REF] and Vainsencher et al. [START_REF] Vainsencher | The sample complexity of dictionary learning[END_REF] in primarily two ways.

First, we relax the assumption that the training data lives in the unit ball [START_REF] Maurer | K-dimensional coding schemes in hilbert spaces[END_REF] or even the unit Euclidean sphere [START_REF] Vainsencher | The sample complexity of dictionary learning[END_REF] by showing that it is sufficient to have sufficient decay of the probability of drawing training samples with "large" norm. This is essentially achieved by replacing Hoeffding's inequality with a more refined Bernstein inequality argument.

Second, and more importantly, we handle penalty functions g beyond indicator functions of compact sets [19, Theorem 1], or of sets K such that sup D∈D,α∈K Dα 2 < ∞ [19, Theorem 2], or of ℓ 1 or ℓ 0 balls [START_REF] Vainsencher | The sample complexity of dictionary learning[END_REF]. Indeed, the first generic case dealt with in this paper involves penalty functions that only need to be non-negative, lower semi-continuous, coercive (and for convenience be minimum at zero). The second generic case also covers as particular cases the indicator function of k-sparse vectors and that of non-negative vectors, with restrictions on D.

Beyond sparse coding, our results provide apparently new sample complexity estimates for K-means. Compared to state of the art results [START_REF] Bartlett | The minimax distortion redundancy in empirical quantizer design[END_REF] and [START_REF] Biau | On the performance of clustering in hilbert spaces[END_REF] we pay an additional log n factor in our excess risk analysis, but improve the sample complexity estimate when the number of clusters K exceeds the ambient dimension m. A brief discussion on possible strategies of how to obtain results independent of m is provided in Section IX.

Similarly, our results provide new sample complexity estimates for NMF improving over state of the art results [START_REF] Maurer | K-dimensional coding schemes in hilbert spaces[END_REF] unless the ambient dimension m exceeds the order d 2 . Again, we refer to Section IX for a discussion of the high-dimensional setting.

Despite its successes, the main limitation of the proposed approach seems to be in the rather crude use of union bounds and covering numbers. One can envision sharper sample complexity estimates using more refined concentration tools. For example, sample complexity results for PCA for signal distributions in the unit ball were studied in [START_REF] Blanchard | Statistical properties of kernel principal component analysis[END_REF]. Although our approach achieves results for more general distributions such as the class P A , the sample complexity results we obtain for this setting seem disappointingly pessimistic, both paying an additional log n factor and overestimating the sample complexity as md-d(d+1)/2 rather than d. We expect this may be due to the use by the authors in [START_REF] Blanchard | Statistical properties of kernel principal component analysis[END_REF] of more refined tools such as McDiarmid's bounded difference concentration inequality for the supremum of the deviation, and more importantly of a clever dimension-independent bound (i.e., independent of m) on the expectation of the supremum of the deviation. We dedicated Section IX to a brief discussion of how we would envision to apply these tools to the more general context considered here, to obtain dimension-independent results.

The reader may have noticed that the rate of convergence of our estimates is in 1/ √ n, which is typical when using techniques based on empirical processes. In certain settings it is possible with substantially more work to achieve fast rates in 1/n, see, e.g., [START_REF] Sridharan | Fast Rates for Regularized Objectives[END_REF] in a much simpler setting, or the discussions of fast rates in [START_REF] Maurer | K-dimensional coding schemes in hilbert spaces[END_REF], [START_REF] Vainsencher | The sample complexity of dictionary learning[END_REF], [START_REF] Blanchard | Statistical properties of kernel principal component analysis[END_REF], [START_REF] Biau | On the performance of clustering in hilbert spaces[END_REF]. It is not clear at this stage whether it is a realistic objective to achieve fast rates in the investigated general setting. In any case this is expected to require well chosen "margin conditions" to hold, see, e.g., [START_REF] Levrard | Fast rates for empirical vector quantization[END_REF], [START_REF]Non Asymptotic Bounds for Vector Quantization in Hilbert Spaces[END_REF] for the case of K-means.

Improving the estimate of Lipschitz constants is another avenue for improvement, although its role in our sample complexity estimates is already only logarithmic. In fact, the generality of our approach comes at a price. The bound on L X (D) in Lemma 5 (resp. Lemma 6) is a worst case estimate given the bounds ( 24) and (25) (resp. ( 36)), and therefore rather crude. In a probabilistic setting, where the training samples x i are drawn i.i.d. according to some distribution P, one can envision much sharper bounds for 1 n (X-DA)A ⊤ ⋆ using matrix concentration inequalities. A key difficulty will come from the control of A which is dependent on X and D. Moreover, while we concentrated on the exploitation of Lemma 3 and its corollary for the special metric • 1→2 on dictionaries, better choices may be possible, including metrics defined in terms of the penalty g when applicable (potentially leading to alternate definitions of ḡ), or metrics dependent on the data distribution P.

The reader may object that the role of the Lipschitz constant L in our final results (see, e.g., Theorem 1) is only logarithmic so the added value of such technicalities might be limited. Here, we notice that for the sake of simplicity we expressed Corollary 3 with a uniform Lipschitz constant independent of D. Clearly, some intermediate steps yield finer estimates L X (D) that depend on the considered dictionary D, i.e., local Lipschitz constants. In certain scenarii these may provide more precise estimates that may turn out to be useful especially in the analysis of local properties of F X (D). This is best illustrated with an example. This is likely to have an impact when performing a local stability analysis of sparse coding with the penalty λ α 2 1 rather than λ α 1 , in the spirit of [START_REF] Jenatton | Local stability and robustness of sparse dictionary learning in the presence of noise[END_REF], see also [START_REF] Gribonval | Dictionary identification -Sparse matrixfactorisation via ℓ 1 minimisation[END_REF]- [START_REF] Schnass | On the identifiability of overcomplete dictionaries via the minimisation principle underlying K-SVD[END_REF].

Finally, observing that this paper provides sample complexity estimates of matrix factorization with penalties g(α) much beyond indicator functions of (often compact) sets, it is natural to wonder if we could also replace the constraint D ∈ D, i.e., the penalty χ D (D), by a more general penalty on D promoting certain dictionaries. This would somehow lead to a model selection process where one could envision, e.g., adapting the effective number of nonzero columns of D, but would raise hard questions regarding the usability of covering numbers.

2) Orthogonal Dictionaries: The orthogonal dictionaries as defined in Equation ( 54) are the group of special orthogonal matrices. The special orthogonal group is a Lie group, and the corresponding Lie algebra is the vector space of skew symmetric matrices so(d), i.e., exp : so(d) → D orth (d) is a surjective map. Furthermore, it can be shown that it is sufficient to restrict so(d) to a ball around 0 with radius π in order to obtain a surjective map. The exponential mapping is Lipschitz since it holds that ρ(exp(S 1 ), exp(S 2 )) ≤ e a S 1 -S 2 , ( 3) Stiefel matrices: The Stiefel manifold, defined in Equation ( 55), can also be defined as the quotient space D orth (m)/ D orth (md), cf. [START_REF] Edelman | The geometry of algorithms with orthogonality constraints[END_REF]. In the previous section, we already mentioned that the Lie algebra of D orth (m) is given by so(m). Now, we define the function Φ : so(m) → D St (d, m), Φ q • exp, where q : D orth (m) → D orth (m)/ D orth (md) denotes the quotient mapping. The Lipschitz constant for the exponential mapping has been previously established in Equation (66). The function q is a quotient map and therefore a contraction, and as the mapping Φ is a combination of two Lipschitz mappings it is Lipschitz itself. Thus, conditions [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] and [START_REF] Elad | On the role of sparse and redundant representations in image processing[END_REF] As this is just a special case of the mapping in the previous constraint set, this mapping is surjective as well as a contraction. Hence, according to Lemmata 17 and16 the covering number bound is given by N (D ⊗St , ǫ) ≤ (3πe π /ǫ) i midi-di(di+1)/2 .

  i ) s. t. D ∈ D. (2)

Theorem 3 .

 3 Under assumptions B1-B3, for any training set X and any dictionaries D ′ = D ∈ D we have

2 β = md • max log λ 2 √2A , 1 ..Example 5 (A;D∈D NMF&ℓ 1 1 2 X -DA 2 F + n i=1 χ

 2215122i=1 For D = D(m, d) we again have the constants C = 3, h = md and Theorem 1 then yields sample complexity for g (ℓ 1 ) Example 4 ("Doubly sparse" coding, training data in the unit ball). Our results also allow to deal with other classes of structured dictionaries, such as the sparse dictionary D sparse (m, d) defined in Equation (65), yielding h = sd and C = 3 m s 1/d . With g χ1 (respectively g χ0 with a restriction to D sparse (m, d) ∩ D LRIP (k, δ)), using Stirling's approximation, Theorem 1 yieldsβ s χ1 ≤ sd • max log 6 √ 8λ + log me s , 1 β s χ0 ≤ sd • log 12 8k/(1δ) + log mes Non-negative Matrix Factorization, training data in the unit ball). In Non-negative Matrix Factorization (NMF) (cf.[START_REF] Seung | Algorithms for non-negative matrix factorization[END_REF]) a matrix X ∈ R m×n is decomposed into a product of two matrices with non-negative entries. It can be expressed as the optimization problem min NMF (α i ).(59)Since the constraint set D NMF&ℓ 1 is a subset of D NMF (1/md), the penalty g = χ NMF and the set D NMF&ℓ 1 satisfy B1-B3 with κ = 1/md, and Table I gives the Lipschitz constant L = 2 √ md. Moreover, since D NMF&ℓ 1 ⊂ D ball (m, d), where D ball (m, d) {D ∈ R m×d : d i 2 ≤ 1} has the same covering number as D(m, d), we get N (D, ǫ) ≤ 3 ǫ md as a direct consequence of the covering number bound

  xi (D) • ḡ(f xi (D)) = 1 n n i=1 2f xi (D) for the Lipschitz constant. To proceed, we recall the definition f x (D) = inf α 1 2 x -Dα 2 2 + g(α). For the function 1 2

Example 8 .

 8 When g(α) = λ α 2 1 , we have ḡ(t) = t/λ hence L X (D) = 1 n n i=1 2f xi (D) • f xi (D)/λ = F X (D) 2/λ. As a result, Corollary 2 implies that we locally have F X (D ′ )/F X (D) ≤ 1 + 2/λ D ′ -D 1→2. This implies that log F X (D) is uniformly Lipschitz with constant L = 2/λ. In other words, F X (D) is more regular where it takes small values.

with S 1 ,

 1 S 2 ∈ B π and a = max ( S 1 , S 2 ) ≤ π and we obtain the Lipschitz constant L = e π . Finally, the dimension of D orth (d) is d(d -1)/2 which yields the covering number bound N (D orth (d), ǫ) ≤ (3πe π /ǫ) d(d-1)/2 .

  in Lemma 17 are fulfilled and taking into account that dim(D St (d, m)) = mdd(d + 1)/2 we obtain the covering number boundN (D St (d, m), ǫ) ≤ (3πe π /ǫ) md-d(d+1)/2 .4) Separable Dictionaries: The mappingΦ : D(m 1 , d 1 ) × . . . × D(m z , d z ) → D sep , (D 1 , D 2 , . . . , D 2 ) → D 1 ⊗ D 2 ⊗ . . . ⊗ D zis surjective and Lipschitz with a Lipschitz constant smaller than one. Thus, Lemmata 17 & 16 yieldN (D sep , ǫ) ≤ (3/ǫ) i midi .5) Tensor Product ofStiefel matrices: Analogous to above we introduce the mapping Φ : D St (m 1 , d 1 ) × . . . × D St (m z , d z ) → D ⊗St , (D 1 , . . . , D z ) → D 1 ⊗ . . . ⊗ D z .

In fact, as in[START_REF] Maurer | K-dimensional coding schemes in hilbert spaces[END_REF] we observe it is sufficient to have d i , d j ≥ 0, i = j

This requirement will only arise from the need to have a finite covering number for the sample complexity estimate

ACKNOWLEDGMENTS

R. Gribonval would like to thank Daniel Vainsencher and Karin Schnass for memorable discussions motivating this work and Gilles Blanchard for invaluable comments.

This work was partially supported by the EU FP7, SMALL project, FET-Open grant number 225913, and by the European Research Council, PLEASE project (ERC-StG-2011-277906). R. Jenatton contributed during his Ph.D. with the Inria -SIERRA Project-Team, This work has been supported by the German DFG funded Cluster of Excellence CoTeSys -Cognition for Technical Systems, and by the DFG grant KL 2189/8-1. {kleinsteuber, m.seibert}@tum.de This paper was presented in part at SSP 2014 -the IEEE Workshop on Statistical Signal

APPENDIX A PROOF OF

 LEMMA 14To begin with, let us notice that by invoking twice the triangle inequality, we have

+ (E {(E {|y| p }) q }) 1/q , so that by using Jensen's inequality, we obtain

thus proving [START_REF] Geng | On the local correctness of ℓ 1 minimization for dictionary learning[END_REF] provided that (50) holds. We now focus on these raw moments. Let us fix some u ≥ 1. We introduce the event

and define l u as the largest integer such that u ∈ [l u , l u + 1). We can then "discretize" the event K c as

We have

where in the last line we used u p ≤ u since u ≥ 1 and p ≤ 1. Using the hypothesis (49), we continue

Upper bounding the discrete sum by a continuous integral, we recognize here the incomplete Gamma function [START_REF] Gautschi | The incomplete Gamma functions since Tricomi[END_REF],

where again we used t pq ≤ t q for t ≥ 1. A standard formula [53, equation (1.3)] leads to, for u ≥ 1, Γ(q + 1, u) = q! exp(-u) q j=0 u j j! ≤ e q! exp(-u)u q . Putting all the pieces together we reach the conclusion.

APPENDIX B COVERING NUMBERS

To estimate covering numbers in Lemma 15 we used implicitly the following Lemma, which follows directly from the definition of covering numbers and will serve to extend this result to other constraint sets.

In order to further extend these results to more elaborate constraint sets we resort to a result from [START_REF] Szarek | Metric entropy of homogeneous spaces and finsler geometry of classical lie groups[END_REF].

Lemma 17. Given a constraint set D with distance measure ρ 1 , a normed space M with distance measure ρ, and a mapping Φ : M → D that fulfill the conditions 1) There exists an r ∈ R, r > 0 such that Φ(B r ) ⊃ D where B r is a ball around 0 ∈ M with radius r with respect to the distance metric ρ (Surjectivity). 2) There exists an

Then the covering number of D has the upper bound

where the exponent h is the dimension of the constraint set D in the sense of its manifold structure.

These lemmata imply covering number bounds for all the constraint sets introduced in the previous section.