1 Introduction↩︎

Sparse dictionary-based coding has been employed for signal and image processing [1], [2], machine learning [3]–[5], compressed sensing [6] and data analytics [7]–[10]. In the sparse-coding framework, observed data is represented as a linear combination of vectors called dictionary atoms. Dictionaries can be either derived analytically or learned from data. Commonly adopted analytical dictionaries include the discrete Fourier transform (DFT), wavelets, and the Ramanujan periodic basis [11]. While they provide structured priors such as signal smoothness over a graph structure via the GFT  [12] or periodicity via the Ramanujan dictionary [11], they may fall short in enabling sparse and accurate representations for data with patterns that do not align well with the predefined atoms. An alternative approach is to learn the dictionaries directly from data which has been shown to enable higher compression rates and better representation quality [13]. The input in the dictionary learning problem is a set of (training) signals, and the goal is to learn both a dictionary and corresponding coding coefficients for the input  [14]–[17]. While many existing techniques focus on one-way (vector) input signals, multi-mode samples require learning dictionaries for each mode. For example, spatiotemporal signals could be sparsely coded by employing jointly a spatial dictionary with atoms corresponding to spatial localities and temporal dictionary with atoms corresponding to trends in time [7]. Multi-dictionary sparse coding with fixed (analytical) dictionaries has been demonstrated beneficial for a range of downstream tasks like compression, missing value imputation, and community detection [7], [8], [18]. Learning dictionaries from data in the multi-way setting requires i) learning multiple dictionaries and ii) estimating coding coefficients which correspond to combinations of atoms. An early two-way (2D) dictionary learning method SeDiL [15] promotes full rank dictionaries and atom coherence via a geometric conjugate gradient optimizer for the dictionary learning step. Follow-up works [19], [20] improve on the original method by adopting 2D-OMP [21] and FISTA [22] for the sparse coding step. The latest methods in this category MOD and CMOD [16] employ gradient projection for dictionary learning. A common limitation of all existing methods is that they do not impose any structure (beyond the usual sparsity) on the coding matrices. In the two-way case the number of coefficients grows quadratically with the sizes (number of atoms) of the left and right dictionaries. This rate of growth of the coding matrices makes the iterative learning of dictionaries and coding challenging.

To address the above challenge we propose to learn two-way dictionaries for low-rank coding matrices. The data model is schematically depicted in Fig. 1 (a) where the input signals \(X\) correspond to different snapshots of user-transportation preference matrices from some organization. We model the inputs as a product of shared left (user) \(L\) and right (transportation type) \(R\) dictionaries and a set of low-rank coding matrices represented as the product \(YW\). The inner dimension (columns of \(Y\) and rows of \(W\)) determines the rank of the coding matrix. Beyond limiting the number of coefficients generating the data, this representation allows us to capture patterns of shared behavior in the learned dictionary atoms, e.g, persistent groups of users who prefer specific types of transportation for the example from Fig. 1 (a). In addition, the low-rank coding model enables more succinct representations of spatio-temporal data as demonstrated in Fig. 1 (b) which compares the error and model sizes of the low rank model (\(LYWR^T\)) and an unconstrained rank model (\(LZR^T\)) with fixed dictionaries (only sparse coding).

Beyond introducing the low-rank two-way dictionary learning problem, we derive a data-distribution-independent sample complexity upper bound for it and back it up empirically. We generalize the bound for loss functions that enforce spread of \(L_1\) mass across the columns/rows of coding matrices to further probe the role of the rank constraint in our model. We also derive a sample complexity lower bound establishing that the rank constraint affects the sample complexity by at most a polylogarithmic factor in the number of dictionary atoms. We show that the general sparsity constrained problem can be solved as an \(L_1\) regularization objective. Informed by the above theoretical results, we propose an efficient optimization procedure \(\text{AODL}\) which outperforms baselines in terms of data reconstruction quality and missing value imputation on multiple datasets.

Our contributions in this work are as follows:
To the best of our knowledge, we are the first to pose the problem of multi-dictionary learning for the setting where coding matrices are constrained to be low-rank and sparse.
We derive both upper and lower bounds on the data complexity of our low-rank dictionary learning problem, elucidating the role of the rank constraint in the sample complexity of the problem, and back our bounds experimentally.
We prove that the general sparsity constrained problem can be solved as an \(L_1\) regularization objective for which we propose an efficient optimization procedure \(\text{AODL}\).
Our proposed method \(\text{AODL}\) consistently produces encodings with lower representation error compared to the closest baselines for a fixed number of coding coefficients (model size). It also imputes missing values more accurately than alternative techniques.
To achieve the same level of reconstruction quality, our approach saves up to \(90 \%\) of the coefficients compared to the closest baselines in large real-world datasets.

2 Related Work↩︎

Sparse coding is widely employed in signal processing [13], [23], [24], image analysis [25] and computer vision [26]. Existing methods can be grouped into three main categories: convex optimization solutions, non-convex techniques, and greedy algorithms [27]. Relaxation techniques impose sparsity on the coding coefficients via L1 regularizers [7], [8], while greedy algorithms select one atom at a time [28]–[30]. Most existing methods focus on 1D signals while we focus on 2D signals, in the sense that most existing methods do not enforce a Kronecker structure for learned dictionaries. methods generalize the one dimensional setting by employing separate dictionaries for each dimension of the data [7], [8], [21], [31], [32]. Some methods in this group place no assumptions about the rank of the encoding matrix [21], [31]–[33], while others impose a low rank on the learned encoding [7], [8]. Our coding coefficient model is inspired by the low-rank models above; however, while the methods above adopt (fixed) analytical dictionaries we learn the dictionaries from data which as we show through experimental comparisons enables more accurate and compact data encoding.

algorithms aim to find one (for vector data) or multiple (for multi-way data) shared dictionaries directly from training samples. The majority of existing work tackles vector signals [13], [34] and iterates between the sparse coding and dictionary learning stages similar to the seminal K-SVD method [14]. Dictionary learning methods for two-way (matrix) data follow a similar alternating process to learn two dictionaries [15], [16], [20], [35]. They solve the same problem as ours except that they do not impose a low-rank structure on the encoding matrices of samples. Some methods impose orthogonality [36], [37] or low-rank [38] constraints on the learned dictionaries motivated by the specific domain applications (e.g., video processing). Another separable dictionaries model [39] learns a mixture of Kronecker-structured dictionaries for tensor data but can be viewed as a general solution for [16]. We compare experimentally to the state-of-the-art methods from this group that learn general unconstrained dictionaries and demonstrate that our low-rank model enables more accurate and compact sparse coding matrices and scales better than dictionary learning alternatives.

3 Preliminaries↩︎

The goal of 1D sparse coding is to represent a vector signal \(x \in \mathbb{R}^{N}\) as a sparse linear combination \(y \in \mathbb{R}^{P}\) of the atom columns of a dictionary \(L \in \mathbb{R}^{N \times P}\) by solving the following optimization problem: \(\min_y f(y)~~\text{s.t.}~~ x=L y,\) where \(f(y)\) is a sparsity-promoting function (e.g., the \(L_1\) norm). In the 2D setting the input is a real-valued matrix \(X \in \mathbb{R}^{N \times M}\) which can be represented as a sparse encoding matrix \(Z\in \mathbb{R}^{P \times Q}\) via a left (column atoms) dictionary \(L \in \mathbb{R}^{N \times P}\) and a right (row atoms) dictionary \(R^T \in \mathbb{R}^{Q \times M}\) by optimizing the following: \[\min_Z f(Z)~~\text{s.t.}~~ X=L Z R^T, \label{eq:2dform}\tag{1}\] where \(f(Z)\) is a sparsity promoting function. A recent alternative 2D model, called TGSD [7], introduced a low-rank structure for the encoding matrix \(Z\), i.e. \(X \approx L Y W R^T\), where \(Y \in \mathbb{R}^{P \times k}\) and \(W \in \mathbb{R}^{k \times Q}\) are sparse dictionary-specific encoding matrices and \(k\) controls the rank of the encoding. All above methods do not learn dictionaries, but instead estimate the sparse coding coefficients.

the goal is to jointly learn the dictionaries from multiple data samples and estimate the per-sample sparse coding coefficients. Existing 2D approaches generalize Eq. (1 ): \[\begin{align} \underset{L,R,Z} {\mathrm{argmin}} || X - L Z R^T||_F^2 + \lambda || Z ||_1, \label{eq:1D-admm} \end{align}\tag{2}\] to learn \(L,R\) in addition to \(Z\). Solvers alternate between i) sparse coding with fixed dictionaries by employing 2D-OMP [21], FISTA [22] and others; and ii) dictionary updates for fixed coding matrices by employing conjugate gradient updates [15] or direct solutions [16].

4 Problem formulation and sample complexity↩︎

4.1 Problem formulation↩︎

The input to our problem is a set of \(S\) samples of \(2\)-way (matrix) data \(\mathcal{X} \in \mathbb{R}^{N \times M \times S}\). Each data sample \(X_s\) is a matrix in \(\mathbb{R}^{N \times M}, s \in [1, \cdots, S]\). The left \(L \in \mathbb{R}^{N \times P}\) and right \(R \in \mathbb{R}^{M \times Q}\) dictionaries have \(P\) and \(Q\) atoms respectively and the two encoding matrices for each sample are denoted as \(Y_s \in \mathbb{R}^{P \times k}\) and \(W_s \in \mathbb{R}^{k \times Q}\), where \(k\) is an encoding rank parameter. Our objective is to learn a set of two-way dictionary atoms \((L, R)\) for which there exists low-rank, sparse encodings of the data samples \(\mathcal{X}\). We measure the quality of the learned dictionaries on a data sample \(X_s \in \mathbb{R}^{N\times M}\) via the following loss function \(\ell:(\mathbb{R}^{N\times P} \times \mathbb{R}^{Q\times M}) \times \mathbb{R}^{N\times M} \to \mathbb{R}\): \[\begin{align} \ell((L, R), X_s) := & \underset{Y_s, W_s} \min \| X_s - LY_sW_sR^T \|_{F}^2, \\ \text{~~s.t.~~} &\mathrm{sparse}(Y_s, W_s) \leq \kappa, \label{eq:obj} \end{align}\tag{3}\] where \(\mathrm{sparse}(\cdot, \cdot)\) is a sparsity-promoting function. In this work, we choose \(\mathrm{sparse}(Y, W) := \max\{ \|Y\|_1, \|W\|_1 \}\) (noting that a bound on \(\mathrm{sparse}(Y, W)\) implies a bound on the sparsity parameter of \(Y\cdot W\)). The task, then, is to use the training set \(\mathcal{X}\) to learn a dictionary pair that has a small expected loss on a new sample \(X\). We denote by \(\mathrm{CoeffMat}(k, \kappa)\) the set of matrices \(Z \in \mathbb{R}^{P\times Q}\) that can be decomposed as \(Z = YW\) satisfying the above dimension and sparsity constraints with constants \(k\) and \(\kappa\).

In a typical statistical learning scenario, it is assumed that the training samples \(\mathcal{X}\) and the new (test) sample \(X\) are drawn independently and identically distributed from an unknown data-generating distribution \(\mathcal{D}\). In this setting, one measures the performance of a dictionary \((L, R)\) via its expected loss on \(X\), known as its risk: \[\begin{align} R_{\mathcal{D}}((L, R)):= R((L, R)) := \mathbb{E}_{X\sim \mathcal{D}}[ \ell((L, R), X) ]. \end{align}\] A learning rule that produces a dictionary \((L, R)\) given \(\mathcal{X}\) is called \((\epsilon, \delta)\)-probably approximately correct (PAC) if, for every data-generating distribution \(\mathcal{D}\), there is a sample size \(m_0 := m_0(\epsilon, \delta)\) such that with probability \(\geq 1-\delta\) over the choice of a size \(m \geq m_0\) training set sampled iid from \(\mathcal{D}\), the learning rule outputs a dictionary \((L, R)\) satisfying: \[\begin{align} R_{\mathcal{D}}((L, R)) \leq \epsilon + \inf_{(L_*, R_*)} R_{\mathcal{D}}((L_*, R_*)). \end{align}\] The minimum \(m_0\) for which this bound holds is called the sample complexity of dictionary learning. We say that the learning rule is PAC if it is \((\epsilon, \delta)\)-PAC for all \((\epsilon, \delta)\) arbitrarily close to \(0\). We denote by \(\mathcal{H}\) (for “hypothesis class”) the set of dictionaries over which the infimum is taken. The most fundamental learning rule is called empirical risk minimization (ERM), which we define in our context next. The empirical risk \(\hat{R}(h, \mathcal{X})\) of a dictionary matrix pair \(h\) on a dataset \(\mathcal{X}:= (X_1, ..., X_S)\) is: \[\begin{align} \hat{R}(h, \mathcal{X}) := S^{-1} \sum_{s \in [1..S]} \ell(h, X_s). \end{align}\] The ERM learning rule chooses a hypothesis that minimizes the empirical risk given an input dataset. ERM is foundational to statistical learning because it is a universal learning rule – whenever a hypothesis class is learnable with finitely many samples, it is learnable via ERM.

4.2 Sample complexity upper bounds↩︎

We next show two upper bounds for the sample complexity of dictionary learning. Each takes a form of a generalization or uniform convergence bound [40]: a high-probability upper bound on \(\sup_{h} |R(h) - \hat{R}(h, \mathcal{X})|\) that holds for all data-generating distributions of interest, where \(R(h)\) is the risk of hypothesis \(h\). Our results provide an accuracy guarantee on the ERM learning rule in terms of the number of samples \(S\).

For large enough sample sizes \(S\), the bound in Theorem 1 decreases monotonically with \(S\) (when all other parameters are fixed). Thus, if one wishes to ensure a risk bound of \(\epsilon\) with probability at least \(1-\delta\) for every possible data-generating distribution, one can set \(e^{-x} = \delta\) in the above expression and the entire expression equal to \(\epsilon\), and solve for \(S\) in terms of \(\epsilon\) and \(\delta\) (possibly using a numerical method). Thus, as is well known in statistical learning theory, a generalization bound translates to a sample complexity bound.

We contrast with similar results in [41]: our bound is in terms of hard sparsity and rank constraints on the coding matrix, whereas the bound in [41] is for the \(L_1\) regularized version of the problem, with no rank constraint. Thus, it is difficult to interpret their bound in terms of exact constraints.

We note that in Theorem 1, the sparsity constraint \(\kappa\) plays a logarithmic role, but the rank constraint does not appear to. To probe the rank constraint further, we note that in the proof, the possibility of involvement of \(k\) disappears when we appeal to the submultiplicativity of the \(L_1\) norm to conclude that for a pair of coding matrices \(Y, W\), \(\|YW\|_1 \leq \|Y\|_1 \|W\|_1\). This inequality is only tight when the \(L_1\) mass of \(Y\) and \(W\) is concentrated in a single column and row, respectively. There exist data-generating distributions for which this occurs, and so we cannot tighten this inequality in the above bound. However, data-generating distributions for which the submultiplicativity bound is not tight arise in reality. To understand this phenomenon further and determine to what extent tightening this inequality may involve the rank constraint, we next formulate a modification of the loss function to include an additional entropy constraint on the coding matrices, along with a corresponding generalization bound.

For matrices \(Y \in \mathbb{R}^{P\times k}, W \in \mathbb{R}^{k \times Q}\), we denote by \(c(Y), r(W) \in \mathbb{R}^{k}\) the following vectors: \(c(Y)_\ell := \sum_{i=1}^P |Y_{i,\ell}|, r(W)_\ell := \sum_{j=1}^Q |W_{\ell,j}|.\) We furthermore define \(\hat{c}(Y) := \frac{c(Y)}{\|Y\|_1}, \hat{r}(W) := \frac{r(W)}{\|W\|_1}\). We note that these are stochastic vectors.

For a stochastic vector \(p \in \mathbb{R}^{k}\), we denote by \(H_2(p)\) the Rényi entropy of order \(2\) of the probability distribution corresponding to \(p\).

With this definition, our original loss function satisfies \(\ell := \ell_0\). We have the following generalization of Theorem 1:

The proof is in Section 8.

In Theorem 2, when \(\rho \geq c\log k\) for some \(c < 2\) (a mild constraint that enforces a very approximately uniform spread of \(L_1\) mass across the columns of \(Y\) and the rows of \(W\)), we have \(\exp(-\rho/2) \leq k^{-c/2}\). When, furthermore, \(\kappa \leq \hat{\kappa} k\) for some fixed \(\hat{\kappa} > 0\), we have \(\kappa^2 \exp(-\rho/2) \leq \hat{\kappa}^2 k^{2-c/2},\) so that we see that the sample complexity increases with \(\hat{\kappa}\) (which one may think of as a per-column/row sparsity constraint) and \(k\). Empirically, we find that for appropriately chosen \(k\), the coding matrices derived via our algorithm applied to real data have nearly maximal \(\rho\), and \(\hat{\kappa} := 10/3\) yields competitive performance with the rank-unconstrained model. Thus, it is empirically reasonable to assume the parameter regime previously discussed, in which decreasing \(k\) decreases the sample complexity.

4.3 Sample complexity lower bound.↩︎

We note that the rank constraint does not appear in Theorem 1 and only appears in the logarithm of Theorem 2. To determine whether or not this is inherent or, alternatively, the result of a loose bound, we next present a sample complexity lower bound for the rank and sparsity-constrained dictionary learning problem.

The proof is in Section 9.

An implication of Theorem 3 is that the rank constraint on the coding matrix cannot improve the sample complexity uniformly over all data-generating distributions except by at most a polylogarithmic factor in \(P, Q\). In fact, the “hard” data-generating distribution that we use in the proof is quite natural: a single “ground-truth” two-way dictionary is fixed, and each data sample is a randomly chosen atom, perturbed by Gaussian noise.

The empirical test error improvements that we see when imposing the rank constraint are thus likely not due to sample complexity improvements, but instead due to improvements in the algorithmic tractability of the optimization problem that we solve, compared to those of competing methods. We also emphasize that in addition to test error improvement, the rank constraint has the advantage of reducing the number of coding coefficients required to describe a dictionary representation of a data matrix.

5 \(\text{AODL}\): Dictionary learning for low-rank sparse coding↩︎

In this section, we describe our algorithm for two-way dictionary learning, establish the convergence of the objective function values of its iterates and discuss potential limitations. Since the direct empirical risk minimization to solve the constrained learning problem is difficult due to the \(L_1\) constraints, we first reformulate the objective as an \(L_1\) regularization and show that its solution is also a solution to the original problem. The regularized objective is as follows: \[\begin{align} \underset{L, R, Y, W} {\mathrm{argmin}} \sum_{s \in [1..S]} ( &|| X_s - L Y_s W_s R^T ||_F^2 + \lambda_1 ||Y_s||_1 + \lambda_2 ||W_s||_1). \end{align} \label{eq:obj95aodl}\tag{4}\]

Our next theorem shows that exactly solving this regularized version of the problem provides an exact solution to the original constrained problem (proof available in Appendix 13). Similar, but not identical, statements are well known to hold for a single optimization variable [42]. We include the statement for our case and its proof for completeness.

Although the objective in (4 ) is not jointly convex in \((L, R, Y, W)\), each sub-problem is convex, leading to a natural alternating optimization solver we call Alternating Optimization (low rank) Dictionary learning (\(\text{AODL}\)). Our algorithm alternates between sparse coding and dictionary learning:

For fixed \(L, R\) and for each \(s \in \{1, ..., S\}\), estimate \(Y_s, W_s\): \[\begin{align} \underset{Y_s, W_s} {\mathrm{argmin}} ( &|| X_s - L Y_s W_s R^T ||_F^2 + \lambda_1 ||Y_s||_1 + \lambda_2 ||W_s||_1). \end{align} \label{eq:obj95l1}\tag{5}\]

For fixed sparse coding matrices \(Y_{s}, W_s\), in each iteration of our algorithm, we then solve the following optimization problem: \[\begin{align} \underset{L, R} {\mathrm{argmin}} \sum_{s \in [1..S]}\| X_s - L Y_s W_s R^T \|_F^2. \label{eq:obj95dicts} \end{align}\tag{6}\]

The overall algorithm is provided in Alg. 2. After initialization (Step 3), we first perform low-rank sparse coding by solving Eq. 5 (Steps 5 - 7) for individual samples \(X_s\) via ADMM (detailed steps of LRSC in Appendix 15). In the second stage (Steps 8-9), we update the two dictionaries \(L, R\) using the gradient projection method to solve Eq. 6 . The sparse coding stage is dominated by an eigendecomposition (see Appendix 15) with a complexity of \(O(P^3 + Q^3 + k^3)\) per sample sparse coding update in the worst case. The dictionary learning stage is dominated by the matrix inversions with complexity \(O(T + P^3 + Q^3)\), where \(T\) is the product of the maximum \(3\) values among \(\{N, M, P, Q, k\}\). Details of all derivations and dictionary initialization strategies are provided in Appendix 15.

We can show that the objective function values of iterates of AODL converge (proof available in Appendix 14):

We note that it is generally difficult to obtain much stronger convergence results, such as results about convergence of iterates themselves, for alternating minimization procedures. In particular, the lack of strict convexity of each coordinate block of the objective function in our setting makes it impossible to apply existing alternating minimization iterate convergence results. Despite this, such procedures are commonly used heuristics with good empirical performance.

To handle samples with missing values and perform imputation, we also introduce a version of our problem with a sample-specific \(0-1\) mask \(\Omega_s\), where the data fit term from Eq. 4 is replaced by \(|| \Omega_s \odot (X_s - L Y_s W_s R^T) ||_F^2\). We derive an ADMM solution for the missing value objective and detail it in Appendix 17.

\(\text{AODL}\) learns good dictionaries and succinct codes when samples allow encoding with low rank \(k\). If the data input is not low-rank, \(\text{AODL}\) will require \(k \approx min(P, Q)\) in the worst case, resulting in coding matrices \(Y\) and \(W\) that exceed the size of the single matrix \(Z\) in the CMOD model. In such sub-optimal settings our scalability and model size advantages would diminish compared to the CMOD baseline. Additionally, while we establish sample complexity and objective function convergence for our problem and method, establishing iterate convergence in this setting, even to a local optimum, remains a challenging open problem.

6 Experimental evaluation↩︎

We characterize \(\text{AODL}\)’s strengths and weaknesses in comparison to state-of-the-art baselines on a range of datasets. We quantify reconstruction quality (as RMSE) and compactness (as number of non-zero coefficients NNZ) of competing models as well as their ability to impute missing values. We also investigate the patterns in the learned dictionaries and empirically test our theoretical results. All tests are conducted on an Intel(R) Xeon(R) Gold 6138 CPU @ 2.00GHz and 251 GB memory server using MATLAB’s R2023a 64-bit version. An implementation of our method is available at: https://tinyurl.com/AODL-demo.

6.1 Datasets and Experiment Setup↩︎

We employ synthetic and real-world datasets summarized in Tbl. 1. The real-world datasets span multiple domains: content exchange (Twitch [44]), web traffic (Wiki [45]), sensor network readings (Road [43], Air [47]) and social interactions (MIT [46]).

Synthetic data. We generate synthetic data according to the model \(X_s = L Y_s W_s R^T + \epsilon\), where \(\epsilon\) is Gaussian noise. Dictionaries \(L \in \mathbb{R}^{20 \times 20}\) and \(R \in \mathbb{R}^{30 \times 30}\) contain random unit-norm atoms and encodings are of rank \(k = 3\), i.e., \(Y \in \mathbb{R}^{20 \times 3}, W \in \mathbb{R}^{3 \times 30}\). For each sample \(Y_s\) and \(W_s\) each contains \(15\) randomly selected coefficients with normally distributed in \(\mathcal{N}(0,1)\). There are a total of \(100\) training samples which are used to learn the dictionaries.

Real-world datasets. We employ \(5\) real-world datasets with temporal and spatial dimensions to be able to evaluate against competing techniques like TGSD [7] employing analytical temporal and graph dictionaries. To prepare the datasets we follow the same protocols as prior work. The datasets span multiple domains: content exchange (Twitch [44]), web traffic (Wiki [45]), sensor network readings (Road [43] and Air [47]) and a social interaction (MIT [46]). In order to create multiple samples for dictionary learning, we slice the data on the larger of its two dimensions (time or spatial/graph extent). We also consider alternative slicing (see Fig. 9 (b)) that confirms our comparative analysis findings. Twitch [44] contains viewer-streamer temporal interactions. We create a graph among viewers and add an edge between a pair of viewers if they viewed the same stream at least 3 times. We use the largest connected component of the co-viewing graph. Values in data samples \(X\in\mathbb{R}^{9000 \times 512}\) represent the number of minutes in any given hour that a viewer spent viewing any streams (i.e., their level of activity). We slice the data randomly into \(30\) samples along the graph dimension.

The Wiki dataset [45] records hourly number of views of Wikipedia articles over \(792\) hours. A co-click graph among articles is constructed by placing edges between articles with at least \(10\) pairwise click events (clicked by the same IPs) within a day. A breadth-first-search (snowball) subgraph of \(11400\) around the China article is selected and then sliced into \(38\) samples along the graph dimension.

The Road [43] dataset consists of \(2780\) highway speed sensors in the LA area. We use the average speed for \(30\) days at \(5\)-minute interval (\(8640\) timesteps) as our signal matrix. The graph is based on connected road segments. We slice the data into 30 samples on the time dimension.

MIT [46] is a communication dataset of timestamped messages between users and a weighted social graph. We split the data along its time dimension.

The Air [47] dataset contains the the number of flights between connected airports (nodes), while the edges connect flight origin and destination. Temporal snapshots represent the number of incoming flights over a 6 hour window. We remove small airports with less than \(8\) flights a day and split the data along its graph dimension obtaining 25 samples.

A note on the selection of real-world datasets: We employ the above spatio(graph)-temporal (ST) data as they have been shown to align to low-rank encoding models due to clustered (shared) behavior present in the temporal and spatial mode observed in the TGSD [7] baseline. Different from the baseline, however, our method AODL does not use graph or temporal information associated with the data but learns temporal and spatial dictionaries from scratch. The learned dictionaries perform better than analytical ones employed by TGSD [7] as demonstrated in in our comparative analysis in the main paper. The ST datasets also allow us to perform qualitative analysis (case studies) and also create controlled synthetic data from GT dictionaries akin to the setup in TGSD.

We compare against one analytical (fixed) dictionary baseline TGSD [7] and three multi-dictionary learning approaches CMOD [16], SeDiL [15], and OSubDiL [39] .

TGSD employs analytical dictionaries but has a low-rank model for the encoding similar to \(\text{AODL}\). Within this baseline we employ the authors’ implementation and the GFT dictionary based on data graphs for graph dimensions and the Ramanujan periodic dictionary for the time dimensions.

Shahriari-Mehr et Al. [16] proposed two methods 2D-CMOD and 2D-MOD for 2D dictionary learning among which 2D-CMOD converged faster to a better solution according to the authors’ experiments. Hence we adopt it as a baseline, and we call it CMOD for brevity in all experiments. We experiment with two versions of CMOD: CMOD-OMP which is the originally proposed method that uses 2D-OMP as a sparse coding solver; and a variant CMOD-ADMM employing an ADMM solver for the sparse coding step. While they produce similar quality results (see columns 11 and 12 in Tbl. 1), the OMP version is about 3 orders of magnitude slower and required over \(47\) hours for a single run on some datasets when the target number of coding coefficients is large. As a result, in all experiments (apart from Tbl. 1) we employ the CMOD-ADMM version.

SeDiL [15] is an older baseline which learns dictionaries employing a conjugate gradient approach. In our experiments, we found that it is sensitive to its hyperparameters and even when tuned extensively, it produces similar or worse results than the newer baseline CMOD [16] while requiring orders of magnitude more time to complete on some datasets (see Tbl. 1). Furthermore, our observations of SeDiL’s performance are consistent with those reported by the authors of CMOD [16]. As a result, we report SeDiL results only in Tbl. 1 and omit it from the comparisons in the rest of the experiments.

OSubDiL [39] learns a Kronecker product of the two sub-dictionaries online, and it requires less memory since it processes one data sample per iteration. However, because it maintains a Kronecker product of sub-dictionaries instead of the sub-dictionaries themselves, it still requires a large amount of memory. The other two models in the same family (STARK and TeFDiL) are omitted since they use much more resources than OSubDiL. In contrast, CMOD and \(\text{AODL}\) also process each sample separately and store the sub-dictionaries directly. As a result, we only report OSubDiL results for the small datasets MIT and AIR in Tbl. 1.

We tune the hyperparameters of all competing techniques by an extensive grid search. Details of the grid search and best parameter values for all baselines are presented in Appendix 19.

We measure the reconstruction quality as the element-wise root mean squared error between a sample \(X_s\) and its reconstruction \(X_s'\): \(\text{RMSE} = \frac{1}{S}\sum_{s}^{S} \sqrt{\frac{\sum_{(i,j)} ({X_s}_{(i,j)} - {X_s'}_{(i,j)})^2}{|X_s|}},\) where \(|X_s|\) denotes the number of elements in \(X_s\). We employ the average number of non-zero coefficients (NNZ) across samples to quantify the size of the encodings produced by competing techniques. We also measure actual running time for competing techniques to compare their scalability.

While the data dimensions \(N\) and \(M\) of the learned dictionaries \(L \in \mathbb{R}^{N \times P}\) and \(R \in \mathbb{R}^{M \times Q}\) are predetermined by the size of the input signals \(X_s\), one has a choice when it comes to the number of atoms in each dictionary (\(P\) and \(Q\)). We employ square dictionaries in all experiments (i.e., \(P=N\) and \(Q=M\)) as this is the minimum number of atoms to form a basis for each of the data dimensions. We also keep the sizes of the analytical dictionaries employed by TGSD the same as those learned by the rest of the competing techniques (GFT is square and for Ramanujan we employ the first \(Q\) atoms when ordered from low to high periods).

6.2 Experiments↩︎

Sparse coding can be viewed as a compressive lossy reconstruction of the input data. We first compare the quality of reconstruction and running time of competing techniques on all datasets for a fixed model size. We report this comparison in Tbl. 1 where the maximal model size (Max. NNZ) is listed in the \(8\)-th column. TGSD is the fastest among competing techniques as it only performs sparse coding and no dictionary learning. However, its quality is dominated by \(\text{AODL}\) in all but the Synthetic dataset since the learned dictionaries allow for a more accurate representation of the data. TGSD has better RMSE than \(\text{AODL}\) on the synthetic dataset since we equip it with the ground truth dictionaries used for generating the data, while these dictionaries are not provided to any of the dictionary learning techniques. \(\text{AODL}\) achieves the smallest error at fixed model size on all real-world datasets and its running time is similar to the ADMM version of CMOD. While CMOD-OMP is slightly better than CMOD-ADMM regarding RMSE, it requires orders of magnitude more time (10s of hours on some datasets) to converge and this gap grows with the NNZ. SeDiL’s running time is inconsistent since the required number of iterations to converge to a low-RMSE solution varies widely across datasets while its best RMSE is similar to that of CMOD and dominated by \(\text{AODL}\).

We next evaluate the trade-off between reconstruction error (RMSE) and coding coefficient size (NNZ) for all techniques on the real-world datasets in Figs. 3 (a)-3 (e). Recall that all competing techniques employ (dense and square) dictionaries of the same size, and hence, the dictionary size is not reported as part of the NNZs. To obtain different points in the RMSE-NNZ space we vary the sparsity regularization hyper-parameters for all competing techniques and we also consider variants of \(\text{AODL}\) with different rank \(k\). \(\text{AODL}\) outperforms baselines on all datasets by consistently producing more accurate models (lower RMSE) at the same level of NNZ. Note that a larger \(k\) and hence larger coding coefficient matrices enables lower RMSE, but the RMSE reduction diminishes with \(k\). In the Road (Fig. 3 (a)) and Wiki’s (Fig. ) datasets \(\text{AODL}\) has the largest relative advantage. For example, to match the best achieved RMSE of CMOD (at \(105k\) NNZ), \(\text{AODL}\) requires an order of magnitude fewer (around 10k) coefficients, while in Wiki \(\text{AODL}\) can save close to \(80\%\) of the coefficients to match CMOD’s quality (\(7k\) vs \(37k\) NNZ). The advantages of \(\text{AODL}\) stem from i) the low-rank model which aligns well with conserved temporal behaviors for subgraphs (spatial sub-regions) and ii) its learned dictionaries specifically tailored to low-rank encoding matrices. Note that TGSD also employs a low-rank encoding model and it tends to work on par or even better than CMOD in low-NNZ regimes. However, data-driven dictionaries (even for a non-low-rank CMOD model) offer an advantage for higher NNZ.

We also evaluate the quality of learned dictionaries for missing value imputation. We employ \(\text{AODL}\) handling missing values and develop a similar missing-value-aware version of CMOD-ADMM (details in Appendix 17). TGSD supports missing value imputation by design. We manually remove a fraction of values in random locations from each data sample and report the RMSE between imputed and actual data for increasing fraction of missing values (ranging from \(10\%\) to \(60\%\)). \(\text{AODL}\) achieves the best RMSE among all competing techniques and across all missing value levels. CMOD is the second best technique demonstrating the benefits of dictionary learning compared to analytical dictionaries employed by TGSD. An exception to these trends is Twitch (Fig. 3 (g)) at high (0.5-0.6) rates of missing values, where TGSD works on par and even better than the CMOD and \(\text{AODL}\). This maybe due to Twitch data being sparser than other datasets limiting the benefits of data-driven dictionaries (CMOD and \(\text{AODL}\)) due to insufficient observed data.

We also study the performance of \(\text{AODL}\) with different number of samples as an empirical validation of Thm. 1. Specifically, we treat \(\text{AODL}\) as a proxy for the empirical risk minimization rule that leads to the sample complexity bound. Given a number of samples from a natural distribution for which the minimum possible risk is \(0\), we demonstrate that \(\text{AODL}\) achieves an error that is less than or equal to the bound. The experimental setup is as follows: (i) sample data generation: \(N\), \(M\), \(P\), \(Q\) are set to \(20\), \(30\), \(20\), \(30\), respectively. Both \(L\) and \(R\) are almost orthogonal with length-1 atoms. For each sample, the coefficient matrices \(Y\) and \(W\) contain 20 and 30 non-zero values (rank \(k=3\)), and we restrict \(||Y||_1=||W||_1=10\). (ii) Theory parameters: max Frobenius norm of the samples \(C=10\), \(1 - e^{-2} = 0.8647\) and \(\kappa = 10\). The test is conducted in a training-testing manner where we (i) generate the training and testing samples using the same ground truth dictionaries and constraints on coefficients; (ii) next we learn the dictionaries from increasing number of training samples; and (iii) use the learned dictionaries to sparse code a fixed testing dataset to calculate the RMSE error.

We report the average test reconstruction error as a function of the number of training samples as well as the distribution of sample errors (inset) in Fig. 4 (a). As expected, the test error decreases with the number of training samples and the error remains relatively stable beyond \(1000\) samples. To support the theoretical bound, we need to show that the test error is less than or equal to the error predicted by the theory. Substituting \(S=1000\) and the remaining theory parameters mentioned above, we obtain an error (RMSE) bound of roughly \(0.64\), which is slightly larger than the error in the figure. Small perturbations of the confidence probability (governed by the parameter \(x\)) changes this number very little. We simplified the data distribution for the purposes of this test in order to be able to reach a test error close to zero. The bar chart in Fig. 4 (a) shows the test error distribution when \(S=1000\) and empirical error of less than \(1\) for all samples, which indicates that the learned dictionaries fit the data well. The gap between the test error achieved by our algorithm and the bound motivates further theoretical investigation discussed at the end of Appendix 7.

We next conduct empirical convergence analysis of \(\text{AODL}\) on synthetic data. The synthetic parameters are set based on the default values outlined in Appendix 19. The ground truth (GT) dictionaries are almost orthogonal. We are interested in quantifying the rate at which the learned dictionaries converge to the GT and quantify this in terms of F1 score. Specifically, we compute inner products between all pairs of learned and GT atoms in a dictionary and pick, without replacement, the top pairs in turn keeping the inner product value as a fractional success (TP) count. The F1 score is then computed based on precision and recall. We measure the F1 score as a function of the number of iterations for \(L\), \(R\) and the 2D atoms computed as the Kronecker products of the two dictionaries (Fig. 4 (b) (green, yellow and red curves). \(\text{AODL}\) recovers the GT dictionaries in about 140 iterations for this setting (\(F1=1\)). We also plot the corresponding RMSE as a function of the iterations in Fig. 4 (b) (blue curve). The RMSE drops significantly in the first few iterations, then becomes relatively stable and converges at around iteration \(150\) when the GT dictionaries are recovered (note that the final RMSE>0 since the data contains noise). Next, we show empirically that we can effectively control the sparsity (i.e., \(\kappa_*\)) of learned coding matrices via regularization. We first predefine some target \(\kappa_*\) values for a given dataset. When learning the dictionaries, we (i) grid search the regularizers (\(\lambda_1, \lambda_2\)) in \(\text{AODL}\); (ii) pick the regularizers that produce a \(\kappa_* = max\{\|Y_s\|_1, \|W_s\|_1\}\) that is close to the target \(\kappa_*\) values; and (iii) report the RMSE and NNZ v.s. \(\kappa_*\) in Figs. 4 (c),. The RMSE decreases for increasing \(\kappa_*\) (blue curves) while the NNZ increases (red curves). The model is less sparsity-constrained at larger \(\kappa_*\) allowing denser coding matrices. The fact that we succeed in finding \((\lambda_1, \lambda_2)\) that yield each target \(\kappa_*\) (even those below the maximum \(L_1\) norm of the coding matrices that generated the data) is an empirical observation that aligns with Thm. 4 (although Thm. 4 cannot be completely confirmed empirically, since we cannot efficiently solve the constrained problem). In a practical scenario with a desired budget (\(\kappa\)) of coefficients, one can perform bisection search on the regularizers to satisfy the desired sparsity level.

Our expectation is that the learned dictionaries by \(\text{AODL}\) capture patterns that are representative of typical sample behavior. To qualitatively investigate this hypothesis, we visualize the learned atoms by \(\text{AODL}\) that are most used for sparse coding. We learn the dictionaries \(L\) (atoms are spatial patterns) and \(R\) (atoms are temporal daily patterns) from the 30 daily samples of the Road dataset. In addition to the atoms, \(\text{AODL}\) has also estimated the encoding matrices \(Y_s \in \mathbb{R}^{P \times k}, W_s \in \mathbb{R}^{k \times Q}\). For a given sample, if a row in \(Y_s\) has at least one non-zero coefficient, we count the corresponding atom as being used for this sample and find the top atom used by the most number of samples. We similarly find the most frequently used right atom based on the columns in \(W_s\) for each sample. We also investigate where and when the top temporal and spatial atoms are used respectively. Assuming that the top graph atom is \(l\), we calculate its alignment \(a(l)\) with each data timestep in \(X_s \in \mathbb{R}^{N \times M}\) as follows: \(a(l) = l^T X_s\). We retain the top 24 best aligned timesteps (spanning 2 hours) for each sample, and sum the total occurrence of these timesteps in all samples. This cumulative score is proportional to the frequency with which the top atom is used at different times of the day. We similarly quantify the locations whose daily time series best align with the top right (temporal) atom \(r\). We then retain the top \(20\) locations by alignment from each sample and quantify their frequency in space.

We visualize the top atoms based on the above frequency definitions in Fig. 5. The top temporal atom is plotted in Fig. 5 (a) where the vertical axis is a proxy for (or the magnitude of) the speed as a function of time of the day (horizontal axis). This atom captures an expected daily commute pattern. High (average) speeds (or low congestion) occur at night (between 8pm and 6am) while the lowest speed coincide with morning (7am-8am) and afternoon (5pm-6pm) rush hours. The atom also has relatively low value for mid-day hours. We also visualize the typical locations where this atom is used in Fig. 5 (b) in which circles designate the locations of road sensors and the circle color designates how often the top temporal atom is used at a given location (more purple colors designates locations that employ the atom more frequently). Among the top locations are intersection of highways 10 and 57 and several in downtown LA (dark blues).

We similarly visualize the top spatial atom in Fig. 5 (c) where the atom values are color-coded. Redder colors correspond to low atom values (i.e., lower speeds) and higher values are color-coded by less red colors like blue and green. Heavy traffic ares areas such as downtown LA and segments of Highway 405 are as expected redder unlike areas far from downtown. We also plot the typical times of the day when the top spatial atom is being used to encode samples in Fig. 5 (d). The top spatial atom is predominantly used to encode daily traffic with a peak in the afternoon hours.

7 Proof of Theorem 1↩︎

In this section, we prove Theorem 1, which gives a generalization bound for two-way dictionary learning. Omitted proofs of auxiliary lemmas and propositions are in Section 11.

Our sample complexity analysis generalizes the one in [48] (see also [41]), which is for the case of a single dictionary, with real-valued vector signals. Our generalization proceeds by viewing our signal matrices as vectors in an appropriate space and dictionary pairs as linear operators on that space.

To derive a sample complexity bound, we will derive what is known as a uniform convergence bound. This is a high-probability bound on the generalization gap function, defined as follows: \(\Psi(\mathcal{X}) := \sup_{h \in \mathcal{H}} |R(h) - \hat{R}(h, \mathcal{X})|.\) The bound must hold regardless of the data-generating distribution \(\mathcal{D}\). Such a bound constitutes a generalization bound for empirical risk minimization, which yields a sample complexity bound.

We first collect standard notions required for our proof.

For a metric space \(M\), we denote its \(\gamma\)-covering number by \(\mathcal{N}(M, \gamma)\). For the purpose of deriving a generalization bound, we must bound the covering numbers of the following loss class, which is a class of functions \(\mathcal{F}:= \{f_h:\Omega \to \mathbb{R}~|~ h \in \mathcal{H}\}\) induced by composing hypotheses in \(\mathcal{H}\) with the loss function.

The relevant metric on \(\mathcal{F}\) is the one induced by the \(L_{\infty}\) norm: for \(f \in \mathcal{F}\), \(\|f\|_{\infty} := \sup_{x \in \Omega} |f(x)|.\) We note that every element of the loss class satisfies \(f_h(x) \in [0, B],\) where \(B := C^2\), the constraint on \(\|x\|_F^2\) imposed by Theorem 1. This is because we can upper bound the minimization in the loss function by setting \(Y = 0, W=0\) in their respective spaces.

We next give a standard generalization bound [48] based on the \(L_{\infty}\) covering numbers of a function class \(\mathcal{F}\).

The \(L_{\infty}\) covering number of \(\mathcal{F}\) can be upper bounded using the covering numbers of \(\mathcal{H}\) along with a bound on the Lipschitz constant of \(f_h\) as a function of \(h\). In order to make this precise, we must specify a norm on \(\mathcal{H}\). Recall the normalization conditions on the outer products of columns of \(L\) with rows of \(R^T\). These translate to an operator norm condition on the mappings induced by pairs \((L, R)\). We have the following standard definition.

It is well known that the \(L_{1}\to L_{2}\) operator norm of a matrix is equal to the maximum \(L_2\) norm of any of its columns. We thus will use the following norm for \(\mathcal{H}\): for any \((L, R) \in \mathcal{H}\), \[\begin{align} \label{expr:hyp-norm} \|(L, R)\|_{op,1,2} := \max_{Z} \frac{\|LZR^T\|_F}{\|Z\|_{L_1}}. \end{align}\tag{7}\] The advantage of this norm is that it turns \(\mathcal{H}\) into a unit ball in a finite-dimensional Banach space, for which covering number bounds are known. Specifically, we have the following lemma.

We next turn to translating the covering number bounds for \(\mathcal{H}\) to bounds for \(\mathcal{F}\). The following well known lemma is the vehicle for this translation.

In the next lemma, we bound the Lipschitz constant for the map \(G:\mathcal{H}\to \mathcal{F}\) defined by \(G(h) := f_h\). We recall that the norm on \(\mathcal{F}\) is the \(L_{\infty}\) norm.

Applying Lemma 2, Lemma 3, and Lemma 4 we get the following bound on the covering numbers of \(\mathcal{F}\).

Applying the bound in Proposition 1 to Lemma 1, we finally get the generalization bound given in Theorem 1.

8 Proof of Theorem 2↩︎

The proof of Theorem 2 is the same as that of Theorem 1 except for two parts: upper bounding \(\|YW\|_1\) after (9 ) and computing a probability-\(1\) upper bound \(B\) on the loss function.

In particular, steps after (9 ) require modification because we require a tighter bound on the Lipschitz constant of the loss function, and we derive this from a tighter bound than submultiplicativity of the \(L_1\) norm for products of matrices affords us.

To that end, we have \[\begin{align} \|YW\|_1 &= \sum_{i=1}^P\sum_{j=1}^Q \left| \sum_{\ell=1}^k Y_{i,\ell}W_{\ell,j} \right| \\ &\leq \sum_{i=1}^P\sum_{j=1}^Q \sum_{\ell=1}^k \left| Y_{i,\ell}W_{\ell,j} \right| \\ &= \sum_{\ell=1}^k \sum_{i=1}^P |Y_{i,\ell}| \sum_{j=1}^Q |W_{\ell,j}| \\ &= \langle {c(Y)}, {r(W)} \rangle \\ &\leq \|c(Y)\|_2 \|r(W)\|_2, \end{align}\] where the first inequality is the triangle inequality, and the second equality is by interchanging summations. The second inequality is by Cauchy-Schwarz.

Next, by multiplying and dividing by \(\|c(Y)\|_1 \|r(Y)\|_1 = \|Y\|_1 \|W\|_1\), we get \[\begin{align} \|c(Y)\|_2 \|r(W)\|_2 &= \|Y\|_1 \|W\|_1 \|\hat{c}(Y)\|_2 \|\hat{r}(W)\|_2 \\ &\leq \kappa^2 \|\hat{c}(Y)\|_2 \|\hat{r}(W)\|_2 \\ &= \kappa^2 \exp(\log( \|\hat{c}(Y)\|_2 \|\hat{r}(W)\|_2)) \\ &= \kappa^2 \exp( -\frac{ H_2(\hat{c}(Y)) + H_2(\hat{r}(W))}{2}). \end{align}\] Finally, applying the constraint \(\rho \leq H_2(\hat{c}(Y)) + H_2(\hat{r}(W))\) completes the modification of the upper bound on the Lipschitz constant of the loss function.

Turning to calculation of the value \(B\), we cannot upper bound the loss function by plugging in \(Y = 0, W=0\), because the stochastic vectors \(\hat{c}(Y)\) and \(\hat{r}(W)\) are not well-defined. We instead plug in \(Y = \alpha \mathbf{1}_{P\times k}, W = \alpha \mathbf{1}_{k \times Q}\), where \(\mathbf{1}_{n\times m}\) denotes the matrix of all \(1\)s in \(\mathbb{R}^{n\times m}\) and \(0 < \alpha < \frac{\kappa}{\max\{P, Q\}k}\) so that the sparsity constraint is satisfied. Note that the Rényi entropies of \(\hat{c}(Y)\) and \(\hat{r}(W)\) are maximized, so that the Rényi entropy constraint is also satisfied. Then \(YW = \alpha^2 k^2 \mathbf{1}_{P\times Q}\). This yields an upper bound on the loss function of \[\begin{align} \|X - LYWR^{T}\|_{F}^2 &\leq (\|X\|_F + \|LYWR^T\|_F)^2 \\ &\leq (C + \alpha^2 k^2 PQ)^2 \\ &= C^2 + \alpha^2\cdot(2C k^2 PQ + \alpha^2 k^4 (PQ)^2). \end{align}\] The first inequality is the triangle inequality. The second is using the upper bound on the Frobenius norm of \(X\) and the operator norm of the dictionary. Then, choosing \(\alpha\) sufficiently small implies that we can set \(B = C^2 + \hat{\alpha}\) for any \(\hat{\alpha} > 0\), so that we can set \(B = C^2\).

This completes the proof of Theorem 2.

9 Proof of Theorem 3↩︎

Here we prove the sample complexity lower bound, Theorem 3. Omitted proofs of lemmas and propositions are in Section 12.

In order to derive a sample complexity lower bound, we must exhibit a “hard” data-generating distribution. We will choose one such that, given \(s\) samples, a successful dictionary learner can be used to implement a successful test for a multiple hypothesis testing problem in which hypotheses are a collection of dictionary atoms, and observations are their noise-perturbed versions. We prove an error probability lower bound for this testing problem using Fano’s inequality.

We define the following set of dictionaries. \[\begin{align} \mathbb{D}:= \{ &(L, R) ~|~ \forall i,j, \|L_i\|_F = \|R_j\|_F = 1, \\ &\mu(L~\otimes~R) \leq \mu = \Theta(1) \}. \end{align}\] Here, \(\mu(M)\) is the pairwise coherence of a matrix \(M\); that is, it is the maximum absolute value of the dot product between any two columns of \(M\).

We next define the sense in which we approximate dictionaries. We consider two dictionaries to be equivalent if one can be obtained from the other by permutation of atoms. This motivates the following dictionary equivalence class metric.

To define the hard data-generating distribution, we construct a packing \(\mathcal{P}\) of the space of dictionaries \(\mathbb{D}\). We guarantee the existence of a large packing in the next lemma.

9.0.0.1 The hard data-generating distribution

Throughout, we fix the rank constraint \(k = 1\) and the sparsity constraint \(\kappa = 1\).

Consider the following data-generating distribution for a sample matrix \(X \in \mathbb{R}^{N\times M}\): a ground-truth dictionary \((L_*, R_*)\) is fixed from \(\mathcal{P}\). We generate a coefficient matrix \(Z_X \in \mathbb{R}^{P\times Q}\) by selecting an element \((i, j) \in [P]\times [Q]\) uniformly at random and placing a \(1\) there. The remaining entries are set to \(0\). We then denote by \(\vec{\epsilon}_X\in \mathbb{R}^{N\times M}\) a matrix whose entries are iid \(\mathcal{N}(0, \sigma^2)\) random variables, for some \(\sigma > 0\). We will choose \(\sigma = \frac{1}{\log(PQ)}\) with foresight. Then \(X = LZ_XR^{T} + \vec{\epsilon}_X\). We also denote by \(X^* := X - \vec{\epsilon}_X\) the noiseless version of \(X\). In other words, \(X\) is a random dictionary atom perturbed entrywise by mean-\(0\) Gaussian noise, and \(X^*\) is its noiseless version.

9.0.0.2 A successful dictionary learner can be used as a successful multiple hypothesis tester from noisy samples

We first show that if we have \(s\) samples \(X := (X_1, ..., X_s)\) from the hard distribution, for some unknown \(D_*\), then an \((\epsilon, \delta)\)-PAC dictionary learner \(\mathcal{A}\) can be used to estimate \(D_*\) to arbitrary precision in the sense of Definition 4. This is the content of Proposition 2 below. We show this by showing that the excess risk \(\Delta R(h) := R(h) - \inf_{h_* \in \mathcal{H}} R(h_*)\) of a hypothesis \(h\) is lower bounded by the squared distance between the hypothesis and the ground truth dictionary.

9.0.0.3 A sample complexity lower bound for the multiple hypothesis testing problem

Our next goal is to prove an error probability lower bound for the multiple hypothesis testing problem associated with \(\mathcal{P}\).

9.0.0.4 Putting everything together for the learning lower bound

We can now combine Propositions 2 and Proposition 3 to prove Theorem 3, the lower bound for dictionary learning in the rank-constrained setting.

Let \(\mathcal{A}\) be a learning rule as in the statement of the theorem. We will use it to construct a tester for the multiple hypothesis testing problem associated with the packing \(\mathcal{P}\) with minimum distance \(\sqrt{PQ}\gamma\). In particular, given a dataset \(X\) generated by a ground-truth dictionary \(D_* \in \mathcal{P}\), \(\mathcal{A}\) outputs a dictionary \(D \in \mathcal{P}\). The test then outputs a dictionary \(\hat{D} \in \mathcal{P}\) in the set \(\mathop{\mathrm{arg\,min}}_{d \in \mathcal{P}} d_{dict}(d, D)\). From Proposition 2, we have that with probability at least \(1-\delta\), \(d_{dict}(D, D_*) \leq C\cdot \sqrt{PQ\epsilon}.\) This implies by the triangle inequality that \[\begin{align} d_{dict}(\hat{D}, D_*) &\leq d_{dict}(\hat{D}, D) + d_{dict}(D, D_*) \\ &\leq \sqrt{PQ}\gamma/2 + C \sqrt{PQ\epsilon}. \end{align}\] Provided that \(\sqrt{\epsilon} < \gamma/2\), we have that \(d_{dict}(\hat{D}, D_*) < \sqrt{PQ}\gamma,\) which implies by the minimum distance property of the packing that \(\hat{D} = D_*\), so that the test constructed from \(\mathcal{A}\) is correct with probability at least \(1-\delta\).

From Proposition 3, we then have that \(\delta \geq 1 - \frac{S\log(PQ)}{\Omega((NP+MQ)\log(1/\epsilon)}.\) This completes the proof of Theorem 3.

10 Conclusion↩︎

In this paper we introduced \(\text{AODL}\), a 2D dictionary learning model with low rank sparse coding for 2D data samples. We established theoretical sample complexity bounds (both upper and lower) for our proposed problem that elucidate the effect on sample complexity of rank and sparsity constraints on coding matrices in the dictionary learning problem. We also proposed an optimization approach, called \(\text{AODL}\), which learns both the low-rank coding matrices and the dictionaries by alternating optimization. We showed that the objective function values of the alternating optimization converge and that a solution to the \(L_1\)-regularized problem solves the original dictionary learning problem with \(L_1\) constraints. We demonstrated the quality of \(\text{AODL}\) on five real-world datasets in comparison to analytical dictionary baselines as well as dictionary learning methods. \(\text{AODL}\) outperformed all state-of-the-art baselines in data reconstruction and missing value imputation tasks. Compared to the best dictionary learning baselines, \(\text{AODL}\) obtained up to \(10\times\) model size reduction for the same representation quality in real-world datasets. It also outperformed baselines for missing value imputation. The atoms learned by our model represented intuitive road traffic patterns.

Appendix↩︎

In this Appendix we provide supplemental content including omitted proofs, optimization steps for \(\text{AODL}\) and baselines, as well as \(\text{AODL}\) in the presence of missing values, additional evaluation results, and hyperparameter tuning protocol to aid reproducibility.

11 Omitted proofs of lemmas used in the proof of Theorem 1↩︎

11.1 Proof of Lemma 4↩︎

We expand out the definition of \(G\) and the \(L_{\infty}\) norm: \[\begin{align} &\|G(L, R) - G(L', R')\|_{\infty} \\ &= \max_{X} | \min_{Y,W} \|X - LYWR^{T}\|_F^2 - \min_{Y',W'} \|X - L'Y'W'R'^T\|_F^2 | \label{expr:starting-lipschitz} \end{align}\tag{8}\] Now, note that \(\min_{Y,W} \|X - L'Y'W'R'^T\|_F^2 \leq \|X - L'YWR'^T\|_{F}^2\), where \((Y, W)\) minimizes the \((L, R)\) term in the above expression. Thus, if we can upper bound \[\begin{align} |\|X - LYWR^{T}\|_F^2 - \|X - L'YWR'^T\|_F^2| \end{align}\] for arbitrary \(Y, W\), this gives us an upper bound on (8 ). We next have \[\begin{align} &|\|X - LYWR^{T}\|_F^2 - \|X - L'YWR'^T\|_F^2| \\ &\leq 2 |\|X - LYWR^{T}\|_F - \|X - L'YWR'^T\|_F| \\ &\leq 2\| LYWR^T - L'YWR'^T\|_{F} \\ &\leq 2 \|(L, R) - (L', R')\|_{op,1,2} \cdot \|YW\|_{1} \label{expr:using-submultiplicativity} \\ &\leq 2 \cdot \|(L, R) - (L', R')\|_{op,1,2}\cdot \|Y\|_{1} \|W\|_1 \\ &\leq 2 \cdot \kappa^2 \|(L, R) - (L', R')\|_{op,1,2}. \end{align}\tag{9}\] This implies the claimed Lipschitz bound.

12 Omitted proofs of lemmas used in the proof of Theorem 3↩︎

12.1 Proof of Lemma 5↩︎

Step 1: Reduce the problem of packing of the quotient space to packing of the ground space. We first claim the following relationship between the packing number of \(\mathbb{D}\) and that of the space \(\mathbb{D}_{ord}\) of ordered dictionaries (of which \(\mathbb{D}\) is a quotient space with respect to equivalence under permutation of atoms): \[\begin{align} \mathcal{M}(\mathbb{D}, \sqrt{PQ}\gamma) \geq \mathcal{M}(\mathbb{D}_{ord}, \sqrt{PQ}\gamma)/(PQ)!. \label{expr:reduction-to-dictsord} \end{align}\tag{10}\]

The proof is as follows: in any packing of \(\mathbb{D}_{ord}\), there are at most \((PQ)!\) members of each dictionary equivalence class. So taking a maximal packing of \(\mathbb{D}_{ord}\) and deleting redundant dictionaries leaves a packing of \(\mathbb{D}\) with cardinality at least the packing number of \(\mathbb{D}_{ord}\) divided by \((PQ)!\).

Thus, it is sufficient for us to lower bound the packing number of \(\mathbb{D}_{ord}\). In what follows, it will also be useful to consider the set \(\mathbb{D}_{ord,unconst}\) of ordered dictionaries without the mutual coherence constraint: \[\begin{align} \mathbb{D}_{ord,unconst} := \{ L ~\otimes~R ~|~ \|L\|_F = \|R\|_F = 1\}. \label{expr:dicts-ord-unconstrained} \end{align}\tag{11}\] We also denote by \(\mathbb{D}_{ord,unconst,L}\) and \(\mathbb{D}_{ord,unconst,R}\) the spaces of left and right ordered dictionaries unconstrained by mutual coherence, respectively.

Step 2: Lower bound the packing number of the Kronecker product space by (approximately) the product of the packing numbers of the left and right dictionary spaces. We next show that if there exist Frobenius norm packings \(\mathcal{P}_L, \mathcal{P}_R\) of the left and right dictionary spaces with distance \(\sqrt{P}\gamma\) and \(\sqrt{Q}\gamma\), respectively, then there exists a packing \(\mathcal{P}\) of \(\mathbb{D}_{ord}\) with cardinality \(|\mathcal{P}_L|\cdot |\mathcal{P}_R|\) and minimum distance \(\sqrt{PQ}\gamma \cdot (1 - o(1))\) as \(\gamma\to 0\).

The proof is as follows: we consider the following packing: \[\begin{align} \mathcal{P}:= \{ L ~\otimes~R ~|~ L \in \mathcal{P}_L, R \in \mathcal{P}_R \}. \end{align}\]

The cardinality of \(\mathcal{P}\) is \(|\mathcal{P}| = |\mathcal{P}_L| \cdot |\mathcal{P}_R|,\) by the unique factorization property of the tensor product.

We next need to verify the distance property of \(\mathcal{P}\). To this end, let \(L_1 ~\otimes~R_1\) and \(L_2 ~\otimes~R_2 \in \mathcal{P}\). Then we have, by direct calculation, \[\begin{align} &\| L_1 ~\otimes~R_1 - L_2 ~\otimes~R_2\|_F^2 \\ &= \|L_1\|_F^2 \cdot \|R_1\|_F^2 - 2 \mathrm{Tr}(L_1^T L_2) \mathrm{Tr}(R_1^T R_2). \label{expr:packing-dist-1} \end{align}\tag{12}\] Now, for each \(b \in \{1, 2\}\), \(\|L_b\|_F^2 = P, \|R_b\|_F^2 = Q,\) by the column normalization condition on the dictionaries. We thus have \[\begin{align} (\ref{expr:packing-dist-1}) = 2PQ - 2 \mathrm{Tr}(L_1^T L_2) \mathrm{Tr}(R_1^T R_2). \label{expr:packing-dist-2} \end{align}\tag{13}\]

We next show that the two trace expressions in (12 ) have simpler expressions in terms of the distance between the marginal dictionaries. In particular, we claim that \[\begin{align} \mathrm{Tr}(L_1^T L_2) = P - \frac{1}{2} \|L_1 - L_2\|_F^2, \\ \mathrm{Tr}(R_1^T R_2) = Q - \frac{1}{2} \|R_1 - R_2\|_F^2. \label{expr:trace-to-distance} \end{align}\tag{14}\] We show this only for \(\mathrm{Tr}(L_1^T L_2)\), as the proof for \(R_1^TR_2\) is entirely analogous. We have \[\begin{align} &\mathrm{Tr}(L_1^T L_2) \\ &= \mathrm{Tr}((L_1 - L_2+ L_2)^T (L_2 - L_1 + L_1)) \\ &= \mathrm{Tr}((L_1 - L_2)^T (L_2 - L_1)) \\ &~~~+ \mathrm{Tr}(L_1^TL_1) + \mathrm{Tr}(L_2^TL_2) - \mathrm{Tr}(L_1^TL_2) \\ &= - \|L_1 - L_2\|_F^2 + 2P - \mathrm{Tr}(L_1^TL_2). \end{align}\] Rearranging and dividing both sides by \(2\) completes the derivation of (14 ).

We thus have \[\begin{align} &(\ref{expr:packing-dist-2}) = 2PQ - 2 \cdot (P - \frac{1}{2} \|L_1 - L_2\|_F^2)(Q - \frac{1}{2} \|R_1 - R_2\|_F^2) \\ &= P \|R_1 - R_2\|_F^2 + Q \|L_1 - L_2\|_F^2 \\ &~~~ - \frac{1}{2}\|L_1 - L_2\|_F^2 \|R_1 - R_2\|_F^2. \label{expr:packing-dist-3} \end{align}\tag{15}\]

Next, we appeal to the distance property of the left and right packings. \((\ref{expr:packing-dist-3}) \geq PQ \gamma^2 - \frac{1}{2}PQ\gamma^4 = PQ \gamma^2 (1 - \frac{1}{2} \gamma^2).\) This completes the proof of the claim reducing the problem of constructing \(\mathcal{P}\) from left and right dictionary packings. In particular, we have shown, for any \(\hat{\gamma}\) small enough and \(> 0\), \[\begin{align} &\mathcal{M}(\mathbb{D}_{ord}, \sqrt{PQ}\hat{\gamma}\sqrt{1 - \hat{\gamma}/2}) \\ &\geq \mathcal{M}(\mathbb{D}_{ord,L}, \sqrt{P}\hat{\gamma})\mathcal{M}(\mathbb{D}_{ord,R}, \sqrt{Q}\hat{\gamma}). \label{expr:reduction-to-marginal-packings} \end{align}\tag{16}\]

Step 3: Lower bound the packing number of the left dictionary space. Next, we claim that the \(\sqrt{P}\gamma\)-packing number of the space \(\mathbb{D}_{ord,L}\) of left dictionaries is at least \(\left( \frac{C}{\gamma} \right)^{\Omega(NP)}\), and the \(\sqrt{Q}\gamma\)-packing number of the space of right dictionaries is at least \(\left(\frac{C}{\gamma} \right)^{\Omega(MQ)}\). The proof for both cases is the same, so we focus only on the one for the left dictionaries.

Letting \(B(r)\) denote the Frobenius ball with radius \(r\) in \(\mathbb{R}^{N\times P}\), we have that \[\begin{align} &\mathcal{M}(\mathbb{D}_{ord,L}, \sqrt{P}\gamma) \geq \frac{\mathrm{Vol}(\mathbb{D}_{ord,L})}{\mathrm{Vol}(B(\sqrt{P}\gamma)} \\ &= \frac{\mathrm{Vol}(\mathbb{D}_{ord,L})}{\mathrm{Vol}(\mathbb{D}_{ord,unconst,L})}\cdot \frac{\mathrm{Vol}(\mathbb{D}_{ord,unconst,L})}{B(\sqrt{P}\gamma)}. \end{align}\] By a standard concentration of measure argument, the ratio \(\frac{\mathrm{Vol}(\mathbb{D}_{ord,L})}{\mathrm{Vol}(\mathbb{D}_{ord,unconst,L})}\) is \(1-o(1)\), and the second ratio is \((C/\gamma)^{(N-1)P}\) for some positive constant \(C\), by known asymptotics for volumes of balls. The same derivation implies analogous asymptotics for the packing numbers of \(\mathbb{D}_{ord,R}\). This implies \[\begin{align} \mathcal{M}(\mathbb{D}_{ord,L}, \sqrt{P}\gamma) &= (1 - o(1))\cdot (C/\gamma)^{(N-1)P} \\ \mathcal{M}(\mathbb{D}_{ord,R}, \sqrt{Q}\gamma) &= (1 - o(1))\cdot (C/\gamma)^{(N-1)Q}. \label{expr:marginal-ordered-packing-nums} \end{align}\tag{17}\]

Step 5: Put everything together. We put everything together as follows: from (16 ) and (17 ), we have \[\begin{align} \mathcal{M}(\mathbb{D}_{ord}, \sqrt{PQ}\hat{\gamma}\sqrt{1-\hat{\gamma}^2/2}) \geq (C_1/\hat{\gamma})^{(N-1)P + (M-1)Q}, \end{align}\] for some constant \(C_1 > 0\). Setting \(\gamma := \hat{\gamma}\sqrt{1-\hat{\gamma}/2}\) and noting that \(\gamma \geq const \cdot \hat{\gamma}\) provided that \(\hat{\gamma}\) is small enough, this implies that \[\begin{align} \mathcal{M}(\mathbb{D}_{ord}, \sqrt{PQ}\gamma) \geq (C/\gamma)^{(N-1)P + (M-1)Q}, \end{align}\] for some \(C>0\).

Finally, applying (10 ) completes the proof of Lemma 5.

12.2 Proof of Proposition 2↩︎

To prove Proposition 2, we need some lemmas. The first lemma gives a characterization of the excess risk of a hypothesis in terms of the expected squared Frobenius distance between a random ground truth dictionary atom and the set of matrices in the range of the hypothesis. To state it, we define the range of a dictionary. \[\begin{align} \mathrm{range}((L, R)) := \{ X ~|~ X = LZR^{T}, Z \in \mathrm{CoeffMat}(k, \kappa) \}. \end{align}\] This is simply the set of matrices \(X\) that can be represented exactly by the dictionary \((L, R)\). We also define \(\mathrm{FullRange}((L, R))\) to be the range of the given dictionary, except that we drop the rank constraint. We note that \(\mathrm{range}((L, R)) \subset \mathrm{FullRange}((L, R))\), and \(\mathrm{FullRange}((L, R))\) is a polytope whose vertices correspond to the atoms of the dictionary.

The second lemma gives the promised bound on \(d_{dict}^2\) in terms of the excess risk.

Lemma 7 has an easy consequence: the ground truth dictionary \(D_*\) is the unique global minimizer of the risk. Moreover, empirical risk minimization recovers it with sufficient data.

We finally have enough to prove Proposition 2 reducing the multiple hypothesis testing problem to PAC learning.

12.3 Proof of Proposition 3↩︎

To prove Proposition 3, we will need a standard corollary of Fano’s inequality [50].

We will use \(\Omega := \mathcal{P}\). To apply Fano’s inequality, we need an upper bound on the KL divergence between the observation distributions induced by pairs of dictionaries in \(\mathcal{P}\). This is the content of the next lemma.

We can now present the proof of Proposition 3, the lower bound for the multiple hypothesis testing problem.

13 Proof of Theorem [THM:REGULARIZATION-SOLVES-CONSTRAINED]↩︎

Here we describe in detail the relationship between the reduction of the constrained optimization problem for sparse coding to the regularized one.

We start with a lemma.

Now, it is easy to check the following monotonicity property of \(F(\gamma)\): if \(\gamma_j \leq \hat{\gamma}_j\) for each \(j\), and if we define \(\kappa := \min\{ k ~|~ F(\gamma) \subseteq G(k)\}\) and \(\hat{\kappa} := \min\{ k ~|~ F(\hat{\gamma}) \subseteq G(k)\}\), then \(\kappa \geq \hat{\kappa}\) and \(F(\hat{\gamma}) \subseteq F(\gamma).\) Moreover, setting \(\gamma = \boldsymbol{0}\) eliminates the constraint, meaning that \(F(\boldsymbol{0}) \subseteq G(\infty)\). This allows for a bisection search to select each regularization parameter to enforce desired constraints. It also has the following consequence: for every \(\kappa > 0\), there exists a \(\gamma\) for which \(F(\gamma) \subseteq G(\kappa)\).

We can apply Lemma 10 to our problem formulation as follows to show that it can be solved via the regularization approach. The constrained problem takes the form \[\begin{align} &\mathop{\mathrm{arg\,min}}_{L,R} \sum_{j=1}^S \ell((L, R), X_j) \\ &= \mathop{\mathrm{arg\,min}}_{L, R} \sum_{j=1}^S \min_{Y_j, W_j ~:~ \|Y_j\|_1, \|W_j\|_1 \leq \kappa} \| X_j - L Y_j W_j R^T \|_{F}^2. \label{expr:constrained} \end{align}\tag{22}\] Using Lemma 10, we can solve the inner minimization by solving the regularized version: defining \[\begin{align} &Y_{*,j}, W_{*,j} \\ &:= \mathop{\mathrm{arg\,min}}_{Y_j, W_j} \|X_j - LYWR^T\|_{F}^2 + \lambda_1 \|Y_j\|_1 + \lambda_2 \|W_j\|_1, \end{align}\] for appropriately chosen \(\lambda_1(L, R, X_j)\) and \(\lambda_2(L, R, X_j)\), we get that \[\begin{align} (\ref{expr:constrained}) = \mathop{\mathrm{arg\,min}}_{L, R} \sum_{j=1}^S \| X_j - LY_{*,j}W_{*,j}R^T \|_{F}^2. \end{align}\]

This shows that solving the sparse coding problem via the regularized objective function provides an exact solution to the original problem, which completes the proof.

14 Proof of Theorem [THM:CONVERGENCE-AODL]↩︎

Here we prove that the sequence of objective function values in our alternating minimization scheme converges. This follows from the fact that any monotone decreasing sequence of real numbers that is bounded below converges to its infimum. To apply this to our case, we simply establish that the sequence of iterates of the dictionaries \((L^{(k)}, R^{(k)})\) and the sparse coding matrices \((Y^{(k)}, W^{(k)}\) yields a monotone decreasing sequence of objective function values. The objective function is bounded below by \(0\), and so this establishes convergence.

For convenience, we introduce notation for the objective function: we denote the objective of the inner minimization by \(J_{coding,s}(X, L, R, Y, W)\), and then the entire objective function becomes \[\begin{align} J(X, L, R, Y, W) := \sum_{s=1}^S J_{coding,s}(X, L, R, Y, W). \end{align}\] For fixed \(L, R, X\), replacing \(Y_s, W_s\) by
\(\mathop{\mathrm{arg\,min}}_{Y_s, W_s}J_{coding,s}(X, L, R, Y, W)\) decreases the value of the objective, and the same holds when we replace \(L, R\) by \(\mathop{\mathrm{arg\,min}}_{L, R} J(X, L, R, Y, W)\). Thus, the value of the objective function is non-increasing in each iteration of the algorithm. This implies convergence of the sequence of objective function values, which completes the proof.

15 \(\text{AODL}\) algorithm and derivation details↩︎

We experiment with two different approaches to initialize the dictionaries in Step 3 of Alg. 2: normally random and on tensor decomposition. In the former approach we sample each element interdependently from a normal distribution \(\mathcal{N}(0, 1)\), and then normalize atoms to unit length. In the latter, we stack all samples in a 3-way tensor and employ Tucker decomposition, which decomposes a tensor into a product of 3 factor matrices and a core tensor with user-defined inner dimensions. We use the mode-1 factor with selected second dimension \(P\) to initialize \(L\), and mode-2 factor with selected second dimension \(Q\) to initialize \(R\). Note that the Tucker initialization of the dictionaries works only when we are aiming to learn a complete or under-complete dictionaries since Tucker requires \(P \leq N, Q \leq M\). We used random initialization for all experiments (and all baselines that learn dictionaries) with synthetic data; and tensor decomposition initialization for all real-world experiments (and all baselines that learn dictionaries). These decisions were based on empirically faster convergence (fewer iterations) for all competing methods with the corresponding initialization schemes.

Given fixed dictionaries, we estimate the sparse coding coefficients one sample at a time. Since for any \(X_s\), the problem is convex in \(Y_s\) and \(W_s\) individually, we can employ the Alternating Direction Method of Multipliers (ADMM) [51] to solve for \(Y_s, W_s\) similar to [7]. We first introduce intermediate variables \(U_s = Y_s, V_s = W_s\), and rewrite the objective for sample \(s\) as:

\[\begin{align} \underset{Y_s, W_s, U_s, V_s} {\mathrm{argmin}}& || X_s - L Y_s W_s R^T ||_F^2 + \lambda_1 ||U_s||_1 + \lambda_2 ||V_s||_1 \\ s.t.&Y_s = U_s, W_s = V_s \end{align} \label{eq:obj95l195admm}\tag{23}\]

We form the corresponding Lagrangian function:

\[\begin{align} &\mathcal{L}(Y_s, W_s, U_s, V_s) \\ = & || X_s - L Y_s W_s R^T ||_F^2 + \lambda_1 ||U_s||_1 + \lambda_2 ||V_s||_1 \\ &+ \frac{\rho_1}{2} || U_s - Y_s + \frac{\Gamma_1}{\rho_1} ||_F^2 + \frac{\rho_2}{2} || V_s - W_s + \frac{\Gamma_2}{\rho_2} ||_F^2. \end{align} \label{eq:obj95lag}\tag{24}\]

We alternate between direct updates of \(Y_s, W_s, U_s\) and \(V_s\) obtained by setting gradients w.r.t. each variable to zero. To update \(Y_s\), we have the following optimization problem: \[\begin{align} \underset{Ys} {\mathrm{argmin}} || X_s - L Y_s W_s R^T ||_F^2 + \frac{\rho_1}{2} || U_s - Y_s + \frac{\Gamma_1}{\rho_1} ||_F^2 \end{align} \label{eq:update95Ys}\tag{25}\]

Setting the gradient with respect to \(Y_s\) to zero, we obtain: \[\begin{align} 2 L^T L Y_s B B^T + \rho_1 Y_s = 2 L^T X_s B^T + \rho_1 U_s + \Gamma_1, \end{align} \label{eq:update95Ys95grad}\tag{26}\] where \(B = W_s R^T\). To simplify we use the following eigendecompositions: \(L^T L = Q_1 \Lambda_1 Q_1^T, B B^T = Q_2 \Lambda_2 Q_2^T\), and set \(\Pi_1 = 2 L^T X_s B^T + \rho_1 U_s + \Gamma_1\). We can then rewrite the above equation as: \[\begin{align} 2 L^T L Y_s B B^T + \rho_1 Y_s =& \Pi_1 \\ 2 Q_1 \Lambda_1 Q_1^T Y_s Q_2 \Lambda_2 Q_2^T + \rho_1 Y_s =& \Pi_1 \\ 2 Q_1^T Q_1 \Lambda_1 Q_1^T Y_s Q_2 \Lambda_2 Q_2^T Q_2 + \rho_1 Q_1^T Y_s Q_2 =& Q_1^T \Pi_1 Q_2 \\ 2 \Lambda_1 Q_1^T Y_s Q_2 \Lambda_2 + \rho_1 Q_1^T Y_s Q_2 =& Q_1^T \Pi_1 Q_2\\ \end{align}\] Allowing \(E_1 = Q_1^T Y_s Q_2\), we obtain \(2 \Lambda_1 E_1 \Lambda_2 + \rho_1 E_1 = Q_1^T \Pi_1 Q_2\), and an element-wise solution for \(E_1\) as follows: \[\begin{align} {[E_1]}_{i,j} = \frac{{[Q_1^T \Pi_1 Q_2]}_{i,j}}{2 {[\Lambda_1]}_{i,i} {[\Lambda_2]}_{j,j} + \rho_1}, \end{align}\] and \(Y_s\) can then be recovered as \(Y_s = Q_1 E_1 Q_2^T\). We follow a similar procedure to derive an analogous update for \(W_s\).

The optimization sub-problems for \(U_s\) and \(V_s\) are: \[\begin{align} \underset{U_s} {\mathrm{argmin}} & \lambda_1 ||U_s||_1 + \frac{\rho_1}{2} || U_s - Y_s + \frac{\Gamma_1}{\rho_1} ||_F^2 \\ \underset{V_s} {\mathrm{argmin}} & \lambda_2 ||V_s||_1 + \frac{\rho_2}{2} || V_s - W_s + \frac{\Gamma_2}{\rho_2} ||_F^2, \end{align}\] with existing closed-form solutions due to [52]: \[\begin{align} {[U_s]}_{i,j} =& \text{sign}({[H_1]}_{i,j}) \times \max{(|{[H_1]}_{i,j}| - \frac{\lambda_1}{\rho_1}, 0)}\\ {[V_s]}_{i,j} =& \text{sign}({[H_2]}_{i,j}) \times \max{(|{[H_2]}_{i,j}| - \frac{\lambda_2}{\rho_2}, 0)}, \end{align}\] where \(H_1 = Y_s - \frac{\Gamma_1}{\rho_1}, H_2 = W_s - \frac{\Gamma_2}{\rho_2}\).

The overall LRSC algorithm is listed in Alg. 6. We first initialize all variables by sampling from a normal distribution \(\mathcal{N}(0, 1)\) (Step 3) and then iterate over the derived 0-gradient updates for each variable in turn: update \(Ys\) (line 5-9); update \(Ws\) (line 10-14); update \(U_s, V_s\) (line 15-17); and update \(\Gamma_1, \Gamma_1\) in the end. The eigendecomposition steps (line 6 and line 11) are the most expensive steps (cubic in their input) since eigendecomposition compared to matrix multiplications requires iterations (depending on the solver). As a result, eigendecomposition is dominating the running time of each iteration of LRSC with a complexity of \(O(P^3 + Q^3 + k^3).\)

We employ gradient projection for dictionary updates given fixed \(Y_s, W_s, \forall s \leq S\) with objective \[\begin{align} \underset{L, R} {\mathrm{argmin}} \sum_{s=1}^{S} (&|| X_s - L Y_s W_s R^T ||_F^2, \end{align} \label{eq:obj95LR}\tag{27}\] and by setting gradients w.r.t. \(L\) to zero, we obtain: \[\begin{align} \sum_{s=1}^{S} -2 (X_s - L Y_s W_s R^T) R W_s^T Y_s^T = 0, \end{align}\] with a closed-form solution for \(L\): \[\begin{align} L = (\sum X_s R W_s^T Y_s^T)(\sum Y_s W_s R^T R W_s^T Y_s^T)^{-1}. \end{align}\] \(R\)’s update is derived in a similar manner and is also listed in Alg. 2 Atoms of both dictionaries are normalized by normalize(\(\cdot\)) in Alg. 2 so the magnitude in representing samples is fully represented in the coding matrices as opposed to the dictionaries. The dictionary updates involve matrix multiplications and matrix inversions, which run in cubic time in the input size. However, matrix inversion is in general much slower than regular multiplication. Thus, the matrix inversion term is dominating the running time with a complexity of \(O(T + P^3 + Q^3)\), where \(T\) is the product of the maximum \(3\) values among \(\{N, M, P, Q, k\}\). \(T\) represents the run time of matrix multiplication inside the inversion term. \(P^3\) is the inversion in the solution of \(L\), and \(Q^3\) is the inversion in the solution of \(R\).

16 CMOD-ADMM derivation↩︎

While the original paper introducing the CMOD method employs 2D-OMP for the sparse coding step, our experimental analysis (Tbl. 1) demonstrated that the reliance on OMP limits the method’s scalability as NNZ grows. For this reason we derive and employ a version of CMOD with ADMM sparse coding in the sparse coding subproblem (with fixed dictionaries \(L\) and \(R\)) is as follows: \[\begin{align} \underset{Z} {\mathrm{argmin}} \sum_{s=1}^{S} (&|| X_s - L Z_s R^T ||_F^2 + \lambda_1 ||Z_s||_1 \end{align} \label{eq:obj95missing95appendix}\tag{28}\]

To obtain an ADMM solution we introduce a proxy variable \(U_s=Z_s\) and obtain the following Lagrangian for sample \(s\):

\[\begin{align} &\mathcal{L}(Z_s, U_s) = || X_s - L Z_s R^T ||_F^2 + \lambda_1 ||U_s||_1 \\ &+ \frac{\rho}{2} || U_s - Z_s + \frac{\Gamma}{\rho} ||_F^2 \end{align} \label{eq:obj95lag95Z}\tag{29}\]

The above equation is similar to Eq. 25 . By simply replacing \(B\) with \(R\), \(Y_s\) with \(Z_s\), we can obtain closed-form updates for \(Z_s\) and \(U_S\) using the same steps.

17 \(\text{AODL}\)(and CMOD) with missing values↩︎

The missing values objective for our problems is:

\[\begin{align} \underset{L, R, Y, W} {\mathrm{argmin}} \sum_{s=1}^{S} (&|| \Omega_s \odot (X_s - L Y_s W_s R^T) ||_F^2 + \lambda_1 ||Y_s||_1 + \lambda_2 ||W_s||_1), \end{align} \label{eq:obj95missing}\tag{30}\] where \(\Omega_s\) is a sample-specific missing value \(1-0\) mask and \(\odot\) denotes the element-wise product. To optimize the missing values objective from Eq. 30 we introduce additional proxy variables \(D_s = X_s, U_s = Y_s, V_s = W_s\), arriving at the following objective:\[\begin{align} \underset{L, R, Y, W} {\mathrm{argmin}} \sum_{s=1}^{S} (&|| D_s - L Y_s W_s R^T ||_F^2 + \lambda_1 ||U_s||_1 \\ &+ \lambda_2 ||V_s||_1 + \lambda_3 || \Omega_s \odot (D_s - X_s)||_F^2) \\ s.t. ~~D_s = X_s, &Y_s = U_s, W_s = V_s, \end{align} \label{meftugid}\tag{31}\]

We form the Lagrangian function for sample \(s\):

\[\begin{align} &\mathcal{L}(D_s, Y_s, W_s, U_s, V_s) \\ = & || D_s - L Y_s W_s R^T ||_F^2 + \lambda_1 ||U_s||_1 \\ &+ \lambda_2 ||V_s||_1 + \lambda_3 || \Omega_s \odot (D_s - X_s)||_F^2 \\ & + \frac{\rho_1}{2} || U_s - Y_s + \frac{\Gamma_1}{\rho_1} ||_F^2 + \frac{\rho_2}{2} || V_s - W_s + \frac{\Gamma_2}{\rho_2} ||_F^2. \end{align} \label{eq:obj95lag95missing}\tag{32}\]

To update \(D_s\), we have the following optimization problem: \[\begin{align} \underset{D_s} {\mathrm{argmin}} || D_s - L Y_s W_s R^T ||_F^2 + \lambda_3 || \Omega_s \odot (D_s - X_s)||_F^2 \end{align} \label{eq:obj95D}\tag{33}\]

Taking the gradient and setting the equation to zero, we get

\[\begin{align} D_s = (L Y_s W_s R^T + \lambda_3 \Omega_s \odot X_s ) \oslash (I + \lambda_3 \Omega_s), \end{align} \label{eq:solve95D}\tag{34}\] where \(\oslash\) is element-wise division.

The update for other variables \(Y_s, W_s, U_s, V_s\) are exactly the same as in Alg. 6, we will omit here.

To update the dictionaries, we have the following objective function which differs slightly from the one in the main paper: \[\begin{align} \underset{L, R} {\mathrm{argmin}} \sum_{s=1}^{S} (&|| D_s - L Y_s W_s R^T ||_F^2. \end{align} \label{eq:obj95LR95missing}\tag{35}\] By setting \(\partial f / \partial L=0\), we get a closed-form solution for \(L\): \[\begin{align} L = (\sum D_s R W_s^T Y_s^T)(\sum Y_s W_s R^T R W_s^T Y_s^T)^{-1}, \end{align}\] and \(R\) can be solved for similarly. Both dictionaries are normalized to unit atom length.

Note that the sparse coding stage of CMOD-ADMM is similar to our objective with the key difference of a single sample-wise encoding matrix \(Z_s\). To derive a missing-value-aware version of CMOD-ADMM we employ the same ADMM approach as the one outlined above for \(\text{AODL}\), with proxy variables for \(Z_s\) and \(X_s\) only: \[\begin{align} \underset{Z} {\mathrm{argmin}} \sum_{s=1}^{S} (&|| D_s - L Z_s R^T ||_F^2 \\ &+ \lambda_1 ||U_s||_1 + \lambda_2 || \Omega_s \odot (D_s - X_s)||_F^2) \\ s.t. ~~D_s = X_s, &Z_s = U_s. \end{align} \label{eq:obj95missing95Z95appendix}\tag{36}\] Updates for the above objective when the dictionaries are fixed are obtained in the same manner as those for \(\text{AODL}\). For the dictionary update stage, again, we just replace \(X_s\) with \(D_s\).

18 Additional experiments↩︎

Next we evaluate the ability of \(\text{AODL}\)(and baselines) to encode data in the presence of noise and also their ability to recover dictionaries similar to the ground truth (GT) dictionaries used to generate the data. In all synthetic tests, TGSD has the advantage of encoding with the GT dictionaries and as a result could be viewed as a GT baseline that CMOD and \(\text{AODL}\) can approach, but not necessarily beat.

In Fig. 7 (a) we vary the noise level quantified as SNR (signal to noise ratio) in the range \([30, 10, 5, 2]\)(dB). Each sample of the ground truth signal is produced via the low-rank encoding model with a total of \(30\) NNZ coefficients and by employing random GT dictionaries. We report the NNZ of each method when for hyperparameter settings with which both CMOD and \(\text{AODL}\) can perfectly recover the GT dictionaries (best matching pairs of GT and learned atoms have a cosine similarity exceeding \(0.99\)). While both CMOD and \(\text{AODL}\) can both recover the GT dictionaries, \(\text{AODL}\) produces this result with less NNZs in its encoding (with largest advantage gap in the noisiest setting SNR = 2). Both \(\text{AODL}\) and CMOD use more coefficients than TGSD to represent the data (the latter uses GT dictionaries), however, thanks to \(\text{AODL}\)’s low rank model it requires fewer coefficients. We further visualize the alignment of learned and GT dictionaries from this experiment in Fig. 8. The NNZ are fixed (up to 80), we plot the alignment of the learned dictionaries with the GT dictionaries. It is clear that \(\text{AODL}\) can almost perfectly recover the GT dictionaries. Meanwhile, the dictionaries learned by CMOD is much noisier. For more complex scenarios, GT dictionaries may not be identifiable.

In Fig. 7 (b), we plot the reconstruction error (RMSE) of all methods at different sparsity levels for SNR=30db. \(\text{AODL}\) is closely aligned with TGSD which uses the GT dictionaries while CMOD requires more coefficients to achieve the same RMSE levels. Figs. 7 (c) and 7 (d) present the RMSE v.s. NNZ trade-off in a much noisier setting (SNR = 2db). While the methods are executed on noisy data (sparse coding from GT dictionary + noise) we seek to quantify the quality of fit to the “clean” component of the samples as well as the noisy samples. Thus, we report the RMSE computed with respect to noisy data (i.e. clean + noise) in Fig. 7 (c); and with respect to only the clean component of the data in Fig. 7 (d). CMOD requires significantly more coefficients than \(\text{AODL}\) to achieve the same RMSE levels. This is likely due to the low-rank encoding model in \(\text{AODL}\) acting as a noise filter. Fig. 7 (d) also suggests that CMOD likely uses coefficients to represent noise since its quality in clean data is worse than that in noisy data at high NNZs. \(\text{AODL}\) is able to capture the clean component in the data much better, and its curve is closer to that of TGSD which employs the GT dictionaries.

To demonstrate the learned dictionary does help in producing better reconstruction error since low rank model by nature provide better NNZ. Instead of comparing \(\text{AODL}\) with the ground truth dictionaries, we compare it with random generated dictionaries. In this test, the basic settings are the same with synthetic test settings. We generated two random dictionaries that serve as the left and the right dictionaries, and we use them directly with our low rank sparse coding model, which is named RAND. We can clearly see in Fig 9, with the same level of NNZ, RAND’s reconstruction error is much higher than \(\text{AODL}\).

As mentioned in the real-world test, we slice the data into multiple samples on the longer dimension to reduce the calculation cost and also prevent learning large dictionaries. However, our goal is just dividing the data into multiple samples to fit our model and the baselines, we are free to slice on any dimension. In the previous test, we slice Twitch data who has \(9000\) nodes and \(512\) timesteps into \(30\) samples on the node direction. Here, we slice it on the time direction into \(16\) samples, each sample has a size of \(\mathbb{R}^{9000 \times 32}\). From Fig 9 (b), we can see our method \(\text{AODL}\) still has some advantage regarding NNZ vs RMSE, however, the difference is vary limited comparing with CMOD. This is because the dictionary \(R \in \mathbb{R}^{32 \times 32}\) has a very small size (atom number), and it is close to the rank parameter \(k\) we choose. We know that the number of coefficients of \(\text{AODL}\) is NNZ(Y) + NNZ(W), while for CMOD, the number of coefficients is NNZ(Z). The size of Y is \(P \times k\); and size of W is \(k \times Q\). So, the total size of \(\text{AODL}\) would be \(P k + k Q = k (P+Q)\). While the size of \(Z\) in CMOD is just \(P \times Q\). The advantage of our method \(\text{AODL}\) will be large when \(k (P+Q) << PQ\). Since we always prefer smaller \(k\), if either \(P\) or \(Q\) is small as well, our advantage will be vanished. As a result, our model will always prefer the two learned dictionaries to have relatively large size, and a small rank \(k\).

19 Hyperparameter tuning and selection↩︎

The parameter settings for all competing techniques unless otherwise specified are as follows. For low rank models like TGSD and \(\text{AODL}\), we set \(k=3\) in synthetic test, and \(k\) takes values in \([15, 30, 45]\) in real-world reconstruction test and missing value imputations tests. In the case study, the dictionaries are calculated at \(k = 15\). In addition to grid search for \(k\), as mentioned in [7], one could estimate the rank \(k\) using SVD. The \(\lambda\)s are sparsity parameter for ADMM method. In CMOD, we only have \(\lambda_1\) since only one encoding matrix is calculated, and in TGSD, we have \(\lambda_1\) and \(\lambda_2\). We grid search these parameter in range \([10^{-3}, 10^{-2}, \cdots, 10^{3}]\) to produce curve with different NNZ and RMSE values. In CMOD-OMP, \(T_0\) is the only parameter, which indicates the target number of coefficients in the encoding matrix. OMP models works much slower as \(T_0\) increases, as a result, we choose some relatively small targets in tests and only report them in Tbl. 1.

SeDiL is another baseline model that requires intensive parameter search. \(q\) means the weight of mixed sparsity measure, which indicates how the sparsity term is being adjusted using power of \(q\). \(\mu\) is the multiplier of the sparse matrix. \(\lambda\) is the Lagrange multiplier of the sparsity term. \(\beta\) is the step size when updating the sparse matrix. \(iter\) is the maximum iteration of the model. We can see that all the above parameters are affecting the performance of the sparsity of the model. We picked the values that can help us to reach the target NNZ, which is defined in Tbl. 1, for a fair comparison of all models.

All hyperparameter ranges are listed in Tbl. 2.

References↩︎

---

1. Contact emails: {bma, amagner, mmcneil2, pbogdanov}albany.edu?↩︎

2. University at Albany – SUNY, Computer Science, Albany, NY, USA↩︎

3. AM funded by NSF CAREER CCF 2338855 and CCF 2212327.↩︎