---
abstract: 'We present a matrix factorization algorithm that scales to input matrices that are large in both dimensions (i.e., that contains more than 1TB of data). The algorithm streams the matrix columns while subsampling them, resulting in low complexity per iteration and reasonable memory footprint. In contrast to previous online matrix factorization methods, our approach relies on low-dimensional statistics from past iterates to control the extra variance introduced by subsampling. We present a convergence analysis that guarantees us to reach a stationary point of the problem. Large speed-ups can be obtained compared to previous online algorithms that do not perform subsampling, thanks to the feature redundancy that often exists in high-dimensional settings.'
author:
- |
    Arthur Mensch\
    Inria Parietal\
    Saclay, France\
    Julien Mairal\
    Inria Thoth\
    Grenoble, France\
    Gaël Varoquaux\
    Inria Parietal\
    Saclay, France\
    Bertrand Thirion\
    Inria Parietal\
    Saclay, France
bibliography:
- 'opt\_bib.bib'
title: |
    Subsampled online matrix factorization\
    with convergence guarantees
---

`firstname.lastname@inria.fr`

#### Setup.

The goal of matrix factorization is to decompose a matrix $\X \in \RR^{p \times n}$ – typically $n$ signals of dimension $p$ – as a product of two smaller matrices: $$\X \approx \D \A
 \quad \text{with}\quad\D \in \RR^{p \times k}, \;\A \in \RR^{k \times n},$$ with potential sparsity or structure requirements on $\D$ and $\A$. We consider a sample stream $(\x_t)_{t \geq 0}$ that cycles into the columns ${\{\x^{(i)}\}}_i$ of $\X$. Matrix factorization can be formulated as a non-convex optimization problem, where the factor $\D$ (the *dictionary*) minimizes the following empirical risk: $$\label{eq:empirical-risk}
    \D = \argmin_{\D \in \mathcal{C}}\;
    \bar f \equaldef \frac{1}{n} \sum_{i=1}^n
    f^{(i)}(\D),\quad\text{where}\quad
    f^{(i)}(\D) =
      \min_{\balpha \in \RR^k} \frac{1}{2}
     \bigl\|
     \x^{(i)}
     - \D \balpha
     \bigr\|_2^2 + \lambda \, \Omega(\balpha).$$ $\mathcal{C}$ is a column-wise separable convex set of $\RR^{p \times
k}$, and $\Omega : \RR^p \rightarrow \RR$ is a penalty over the code. The problem of *dictionary learning* [@olshausen_sparse_1997; @agarwal_learning_2014] sets $\mathcal{C} = \mathcal{B}_2^k$ and $\Omega = \Vert \cdot \Vert_1$. Due to the sparsifying effect of $\ell_1$ penalty [@tibshirani_regression_1996], the dictionary forms a basis in which the data admit a *sparse* representation. Setting $\mathcal{C} = \mathcal{B}_ 1^k$ and $\Omega = \Vert \cdot \Vert_2^2$ yields a data-adapted *sparse basis*, akin to sparse PCA[@zou_sparse_2006]. The algorithm presented here accommodates *elastic-net* penalties $\Omega(\balpha) \triangleq (1- \nu) \Vert \balpha \Vert_2^2 + \nu \Vert \balpha \Vert_1$, and elastic-net ball-constraints $\mathcal{C} \triangleq \{ \D \in \RR^{p \times k}/\:
 \Vert \d^{(j)} \Vert \triangleq \Vert \d^{(j)} \Vert_1 + (1 - \mu) \Vert \d^{(j)} \Vert_2^2 \leq 1 \}$.

#### Problem.

For many applications of matrix factorization, datasets are growing in both sample number and sample dimension. The online matrix factorization algorithm [@mairal_online_2010] can handle large numbers of samples but was designed to work in relatively small dimension. Recent work [@mensch_dictionary_2016] has adapted this algorithm to handle very high dimensional dataset. Though it demonstrates good empirical performance, the proposed algorithm yields sequences of iterates with non vanishing variance and is not asymptotically convergent.

#### Contribution.

We address this issue and correct some aspects of the algorithm in [@mensch_dictionary_2016] to establish convergence and correctness. We thus introduce a new method to efficient solve for numerous high-dimensional data — large $p$, large $n$ — with theoretical guarantees.

As in [@mensch_dictionary_2016], we perform *subsampling* at each iteration of the online matrix factorization algorithm. The sample stream is downsampled in a *stochastic* manner —we observe different subsets of successive colunms. We thus each step of the online algorithm in a space of reduced dimension $q < p$. Unlike [@mensch_dictionary_2016], we control the variance introduced by the subsampling to obtain convergence. For this, we rely on low-dimensional statistics kept from the past, as many recent stochastic algorithms [@schmidt_minimizing_2013; @johnson_accelerating_2013; @defazio_saga:_2014].

Algorithm
=========

#### Original online algorithm.

The problem can be solved online, following [@mairal_online_2010]. At each iteration $t$, a sample $\x_t$ is drawn from one of the columns $\{ \x^{(i)} \}_{1 \leq i \leq n}$ of $\X$. Its code $\balpha_t$ is computed from the previous dictionary $\D_{t-1}$: $\balpha_s \triangleq \argmin_{\balpha \in \RR^p} \frac{1}{2}
 \bigl\|
 \x_s
 - \D_{s-1} \balpha
 \bigr\|_2^2 + \lambda \, \Omega(\balpha)$. Then, $\D_t$ is updated as $$\label{eq:minimization}
    %  \begin{split}
     \D_t \in \argmin_{\D \in \mathcal{C}} \bar g_t(\D)
     \triangleq \Big(
     \frac{1}{t}\sum_{s=1}^t \frac{1}{2}
     \bigl\|
     \x_s
     - \D \balpha_s
     \bigr\|_2^2 + \lambda \Omega(\balpha_s) \Big),
    %  \balpha_s &\triangleq \argmin_{\balpha \in \RR^p} \frac{1}{2}
    %  \bigl\|
    %  \x_s
    %  - \D_{s-1} \balpha
    %  \bigr\|_2^2 + \lambda \, \Omega(\balpha).
    %  \end{split}$$ In other words, $\D_t$ is chosen to be the best dictionary that relates past codes $(\balpha_s)_{s \leq t}$ to past samples $(\x_t)_{s \leq t}$. Codes are not recomputed from the current iterate $\D$, which would be necessary to compute the true past loss fonction $\bar f_t(\D)$, of which $\bar g_t$ is a strongly-convex upper-bound: $$\begin{aligned}
    \label{eq:bar_ft}
    \bar f_t(\D) &\triangleq \frac{1}{t} \sum_{s=1}^t \min_{\balpha \in \RR^p} \frac{1}{2}
        \bigl\|
        \x_s
        - \D \balpha
        \bigr\|_2^2\!+ \lambda \Omega(\balpha) \leq \bar g_t(\D)\end{aligned}$$ It can be shown [see @mairal_stochastic_2013 for theoretical grounding] that minimizing $(\bar g_t)_t$ yields a sequence of iterates that is asymptotically a critical point of $\bar f$ defined in . $\bar g_t$ can be minimized efficiently by projected block coordinate descent, which makes it useful in practice. Indeed, minimizing $\bar g_t$ is equivalent to minimizing the quadratic function $$\label{eq:full_quadratic}
    \D \to\frac{1}{2} \trace (\D^\top \D \bar \C_t)
      - \trace (\D^\top \B_t),\:\text{where}\quad
      \bar \B_t = \frac{1}{t} \sum_{s=1}^t \x_s \balpha_s^\top, \quad
      \bar \C_t = \frac{1}{t} \sum_{s=1}^t \balpha_s \balpha_s^\top.$$ Its gradient $\nabla \bar g_t: \D \to \D \bar \C_t - \bar \B_t$ can be tracked online by updating $\bar \C_t$ and $\bar \B_t$ at each iteration: $$\label{eq:parameter-aggregation}
    \bar \C_t = (1 - \frac{1}{t}) \bar \C_{t-1}
    + \frac{1}{t} \balpha_t \balpha_t^\top \qquad
    \bar \B_t = (1 - \frac{1}{t}) \bar \B_{t-1}
    + \frac{1}{t} \x_t \balpha_t^\top.$$ Those two statistics are thus sufficient to yield the sequence $(\D_t)_t$. The weight $\frac{1}{t}$ used above can be replaced by a more general $w_t$. In addition, the online algorithm has a minibatch extension.

#### Subsampled algorithm.

We adapt the algorithm from [@mairal_online_2010] to handle large sample dimension $p$. The complexity of this algorithm linearly depends on the dimension $p$ in three aspects:

-   $\x_t \in \RR^p$ is used to compute the code $\balpha_t$,

-   it is used to update the surrogate parameters $\bar \C_t \in \RR^{p\times k}$,

-   $\D_t \in \RR^{p\times k}$ is fully updated at each iteration.

Our new *subsampling online matrix factorization* algorithm () reduces the dimensionality of each of these steps, so that the single-iteration complexity in $p$ depends on $q = \frac{p}{r}$ rather than $p$. $r > 1$ is a *reduction factor* that is close to the computational speed-up per iteration in the large dimensional regime $p \gg k$. Formally, we randomly draw, at iteration $t$, a mask $\M_t$ that “selects” a random subset of $\x_t$. $\M_t$ is a $\RR^{p\times p}$ random diagonal matrix, such that each coefficient is a Bernouilli variable with parameter $\frac{1}{r}$, normalized to be $1$ in expectation. With this definition at hand, $\M_t \x_t$ constitutes a non-biased, low-dimensional estimator of $\x_t$: $
    \EE[\Vert \M_t \x_t \Vert_0] = \frac{p}{r} = q$, and $
    \EE[\M_t \x_t] = \x_t
$, with $\Vert \cdot \Vert_0$ counting the number of non-zero coefficients. Thus, $r$ is the average proportion of observed features at each iteration. We further define the pair of orthogonal projectors $\P_t \in \RR^{q \times p}$ and $\P_t^\perp \in \RR^{(p - q)\times p}$ that projects $\RR^p$ onto $\mathrm{Im}(\M_t)$ and $\mathrm{Ker}(\M_t)$, which we will use for the dictionary update step.

In brief, , defined in Alg. \[alg:somf\], follows the outer loop of online matrix factorization, with the following major modifications at iteration $t$:

-   it uses $\M_t \x_t$ and low-size statistics instead of $\x_t$ to estimate the code $\balpha_t$ and the surrogate $g_t$,

-   it updates a subset of the dictionary $\P_t \D_{t-1}$ to reduce the surrogate value $\bar g_t(\D)$. Relevant parameters of $\bar g_t$ are computed using $\P_t \x_t$ and $\balpha_t$ only.

We describe in detail the new code computation and dictionary update steps. We then state convergence guarantees for . Those are non trivial to obtain as is not an exact (stochastic) majorization-minimization algorithm.

Initial iterate $\D_0$, weight sequences ${(w_t)}_{t>0}$, ${(\gamma_c)}_{c>0}$, sample set ${\{\x^{(i)}\}}_{i> 0}$, \# iterations $T$. Draw $\x_t = \x^{(i)}$ at random and $\M_t$ (see text) Update the regression parameters for sample $i$: $c^{(i)} \gets c^{(i)} + 1$, $\gamma \gets \gamma_{c^{(i)}}$ $$(\bbeta_t^{(i)}, \G_t^{(i)}) \gets (1 - \gamma) (\bbeta_{t-1}^{(i)}, \G_{t-1}^{(i)})
             + \gamma (\D_{t-1}^\top \M_t \x^{(i)}, \D_{t-1}^\top \M_t \D_{t-1}),\,
             ( \bbeta_t, \G_t) \gets (\bar \bbeta_t^{(i)}, \bar \G_t^{(i)})$$ Compute the approximate code for $\x_t$: $\balpha_t \gets \argmin_{\balpha \in \RR^k}
            \frac{1}{2} \balpha^\top  \G_t \balpha -
            \balpha^\top  \bbeta_t + \lambda \, \Omega(\balpha).$ Update the parameters of the aggregated surrogate $\bar g_t$: $$\label{eq:somf_partial}
             \bar \C_t \gets (1 - w_t) \bar \C_t + w_t  \balpha_t  \balpha_t^\top. \qquad
             \P_t \bar \B_t \gets (1 - w_t ) \P_t \bar \B_t + w_t \P_t \x_t  \balpha_t^\top.$$ Compute simultaneously (using (using [@mensch_dictionary_2016 Alg 2]) for expression): $$\label{eq:somf_minimization}
            \P_t \D_t \gets \argmin_{\D^r \in \mathcal{C}^r}
            \frac{1}{2} \trace ({\D^r}^\top (\D^r \bar \C_t - \bar \B_t)),\:
             \P_t^\perp \bar \B_t \gets (1 - w_t ) \P_t^\perp \bar \B_{t-1} + w_t \P_t^\perp \x_t  \balpha_t^\top.$$ Final iterate $\D_T$.

#### Code computation.

In the online algorithm, $\balpha_t$ is obtained solving the linear regression problem $$\label{eq:regression}
    \balpha_t = \argmin_{\balpha \in \RR^k} \frac{1}{2} \balpha^\top \G_t^\star \balpha - \balpha^\top
\bbeta_t^\star + \lambda \Omega(\balpha),\quad\text{where}\quad
    \G_t^\star = \D_{t-1}^\top \D_{t-1}\:\text{and}\:\bbeta_t^\star = \D_{t-1}^\top \x_t$$ For large $p$, computing $\G_t^\star$ and $\bbeta_t^\star$ dominates the complexity of the code computation step. To reduce this complexity, we introduce *estimators* for $\G_t$ and $\bbeta_t$, computable at a cost proportional to $q$, whose use does not break convergence. Recall that the sample $\x_t$ is drawn from a finite set of samples ${\{\x^{(i)}\}}_i$. We estimate $\G_t^\star$ and $\bbeta_t^\star$ from $\M_t \x_t$ and data from previous iterations $s < t$ for which $\x^{(i)}$ was drawn. Namely, we keep in memory $2n$ estimators, written ${( \G_t^{(i)}, \bbeta^{(i)}_t)}_{1\leq i \leq
n}$, observe the sample $i = i_t$ at iteration $t$ and use it to update the $i$-th estimators $\bar \G_t^{(i)}$, $\bar \bbeta^{(i)}_t$ following $$\bbeta_t^{(i)} = (1 - \gamma)  \G_{t-1}^{(i)} + \gamma \D_{t-1}^\top \M_t \x^{(i)},\qquad
        \G_t^{(i)} = (1 - \gamma)  \G_{t-1}^{(i)} + \gamma \D_{t-1}^\top \M_t \D_t^{(i)},$$ where $\gamma$ is a weight factor determined by the number of time sample $i$ has been previously observed at time $t$. Precisely, given ${(\gamma_c)}_c$ a decreasing sequence of weights, $\gamma \triangleq \gamma_{c^{(i)}_t}$, where $
    c^{(i)}_t =
         | \lbrace s \leq t, \x_s = \x^{(i)} \rbrace |
$ All others estimators $\{\G^{(j)}_t, \bbeta^{(j)}_t\}_{j \neq i}$ are left unchanged from iteration $t-1$. The set ${\{ \G_t^{(i)}, \bbeta^{(i)}_t\}}_{1\leq i\leq n}$ is used to define the *averaged* estimators at iteration $t$, related to sample $i$: $$\label{eq:agg-estimates}
     \G_t \triangleq  \G_t^{(i)} = \sum_{s \leq t, \x_s = \x^{(i)}} \gamma_{s,t}^{(i)} \D_{s-1}^\top \M_s \D_{s-1},\quad
     \bbeta_t \triangleq  \bbeta_t^{(i)} = \sum_{\substack{s \leq t, \x_s = \x^{(i)}}} \gamma_{s,t}^{(i)} \D_{s-1}^\top \M_s \x^{(i)},$$ where $\gamma_{s,t}^{(i)} = \gamma_{c^{(i)}_t} \prod_{s < t, \x_s = \x^{(i)}} (1 -
\gamma_{c^{(i)}_s})$. Replacing $(\G_t^\star, \bbeta_t^\star)$ by $(\G_t, \bbeta_t)$ in , $\balpha_t$ minimizes the masked loss averaged over the previous iterations where sample $i$ appeared: $$\label{eq:approx-regression}
     \min_{\balpha \in \RR^k} \sum_{s\leq t\\
     \x_s = \x^{(i)}} \frac{\gamma_{s,t}^{(i)}}{2}
     \Vert \M_s(\x^{(i)} - \D_{s-1}^\top \balpha) \Vert_2^2
      + \lambda \Omega( \balpha ).$$ The sequences ${(\G_t)}_t$ and ${(\bbeta_t)}_t$ are *consistent* estimations of ${(\G_t^\star)}_t$ and ${(\bbeta_t^\star)}_t$ — consistency arises from the fact that a single sample $\x^{(i)}$ is observed with different masks along iterations. This was not the case in the algorithm proposed in [@mensch_dictionary_2016], which use estimators that does involve averaging from past data. Solving  is made closer and closer to solving , in a manner that ensures the correctness of the algorithm. Yet, computing the estimators  is $r$ times as costly as computing $\G_t^\star$ and $\bbeta_t^\star$ from  and permits to speed up the code computation steop close to $r$ times. The weight sequences $(w_t)_t$ and $(\gamma_c)_c$ are selected appropriately to ensure convergence. For instance, we can set $w_t = \frac{1}{t^v}, \gamma_c = \frac{1}{c^{2.5 - 2v}}$, with $v \in (\frac{3}{4}, 1)$.

#### Dictionary update.

In the original online algorithm, the whole dictionnary $\D_{t-1}$ is updated at iteration $t$. To reduce the time complexity of this step, we add a “freezing” constraint to the minimization of the quadratic function . Every row $r$ of $\D$ that corresponds to an *unseen* row $r$ at iteration $r$ (such that $\M_t[r, r] = 0$) remains unchanged. $\D_t$ is obtained thus solving $$\label{eq:dict-update-cons}
    \D_t \in \argmin_{\substack{\D \in \mathcal{C}\\
    \P_t^\perp \D = \P_t^\perp \D_{t-1}}}
     \frac{1}{2} \trace (\D^\top \D \bar \C_t)
      - \trace (\D^\top \bar \B_t),\,\text{with } \P_t\text{ orth. projector on }\mathrm{Im}(\M_t)$$ With elastic-net ball constraints, solving reduces to performing the partial dictionary update  (Alg. \[alg:somf\]), with $\mathcal{C}^r = \{\D^r \in \RR^{k \times q}, \Vert {(\d^r)}^{(j)} \Vert \leq 1 -
\Vert \d^{(j)}_{t-1} \Vert + \Vert \P_t \d_{t-1}^{(j)} \Vert\}$. We perform this update using a single pass of projected block coordinate descent with blocks in the reduced space $\RR^{q}$. The dictionary update step is thus performed $r$ times faster than the original algorithm.

#### Surrogate computation

The gradient we use to solve  requires to know only $\bar \C_t$ and $\P_t \bar \B_t$. We thus parallelize the partial update of the dictionary and the update of $\P_t^\perp \bar \B_t$, using a second thread. The update of $\P_t \bar \B_t$ is performed in the main thread at a cost proportional to $q$. As the parallel computation is dominated by dictionary update, this is enough to effectively reduce the computation time of $\bar \B_t$ computation.

#### Convergence guarantees

All in all, the three steps whose complexity depends on $p$ in the original algorithm now depends on $q$, which speeds-up a single iteration by a factor close to $r$. Yet retain convergence guarantees. We assume that there exists $\nu$ such that for all $t > 0$, $\D_t^\top \D_t \succ \nu \I$ (met in practice or by adding a small ridge regularization to ), and make a technical data-independent hypothesis on $(w_t)_t$ and $(\gamma_c)_c$ decay. Then iterates converge in the same sense as in the original algorithm [@mairal_online_2010]: $\D_t \to \D_\infty \in \RR^{p\times k}$ and $\nabla \bar f(\D_\infty, \D - \D_\infty) \geq 0$ for any $\D \in \mathcal{C}$, where $\bar f$ is the empirical risk defined in .

We use the aggregated nature of $\bar g_t$ and the fact that we can obtain a *geometric* rate of convergence for a single pass of projected block coordinate descent [*e.g* see @richtarik_iteration_2014] to control the terms $\D_t - \D_{t-1}$ and $\bar g_t(\D_t) - \bar g_t(\D_t^*)$, where $\D_t^* =
\argmin_{\D \in \mathcal{C}} \bar g_t$. We obtain $\bar g_t(\D_t) - \bar
g_t(\D_{t-1}) = \mathcal{O}(w_t)$, a crucial result on estimate stability. Simultaneously, we show that the partial minimization yields the same result as the full minimization for $t \to \infty$, as $\theta_t - \theta_t^\star \to 0$. Using this result, with appropriate selection of $(w_t)_t$ and $(\gamma_t)_t$, the noise induced by the use of estimators in can be bounded in the derivations of [@mairal_online_2010]. We then write the difference - as the sum of a *lag* term and an empirical mean over ${\{\M_s\}}_{\x_s =\x_t}$. Both can be bounded with appropriate selection of weights. Formally, ${(\bar g_t)}_t$ are no longer upper-bounds of ${(\bar f_t)}_t$, but become so for $t \to \infty$, at a sufficient rate to guarantee convergence.

Experiments
===========

#### Hyperspectral images.

We benchmark our algorithm by performing dictionary learning on a large hyperspectral image. Dictionary learning is indeed used on patches of hyperspectral images, as in [@maggioni_nonlocal_2013]. Extracting $16{\times}16$ patches from an 1GB hyperspectral image from the AVIRIS project with $224$ channels yields samples of dimension $p = 57,000$. Figure 1 demonstrates that the newly proposed algorithm is faster than the original, non-subsampled algorithm from [@mairal_online_2010] by a factor close to $r = 4$. This speed-up is helped by the redundancy in the different channels of the hyperspectral patches. In the first epochs, the proposed method also outperforms the recent subsampled algorithm from [@mensch_dictionary_2016], thanks to the introduction of consistent estimators in the code computation step. All algorithms are implemented in Cython and benched on two cores. They cycle over $100,\!000$ normalized samples with minibatches of size 50. The code for reproduction is available at [github.com/arthurmensch/modl](github.com/arthurmensch/modl)..

#### Figure 1

Performing dictionary learning on hyperspectral data (224 channels, $16\times16$ patches) is faster with stochastic subsampling, and even faster with the newly proposed variance control. \[fig:opt\]

![image](opt_bench.pdf){width="\textwidth"}

#### Conclusion

The new algorithm can efficiently factorize large and tall matrices. It preserves the speed gains of [@mensch_dictionary_2016] but has convergence guarantees that are beneficial in practice.
