A Unified Framework for Stochastic Matrix Factorization via Variance Reduction

TL;DR

This paper introduces a unified variance reduction framework for stochastic matrix factorization, significantly improving convergence speed and accuracy across various formulations through theoretical analysis and extensive experiments.

Contribution

It presents a general framework that accelerates stochastic matrix factorization algorithms and provides non-asymptotic convergence guarantees, unifying multiple existing methods.

Findings

Faster convergence compared to state-of-the-art methods

More accurate dictionary recovery in experiments

Theoretical bounds on convergence and complexity

Abstract

We propose a unified framework to speed up the existing stochastic matrix factorization (SMF) algorithms via variance reduction. Our framework is general and it subsumes several well-known SMF formulations in the literature. We perform a non-asymptotic convergence analysis of our framework and derive computational and sample complexities for our algorithm to converge to an $\epsilon$-stationary point in expectation. In addition, extensive experiments for a wide class of SMF formulations demonstrate that our framework consistently yields faster convergence and a more accurate output dictionary vis-\`a-vis state-of-the-art frameworks.

Equations

$${\ell }_1(\mathbf{y},\mathbf{W})\triangleq \underset{\mathbf{h}\in \mathcal{H}}{\min }{\tilde{\ell }}_1(\mathbf{y},\mathbf{W},\mathbf{h}),{\tilde{\ell }}_1(\mathbf{y},\mathbf{W},\mathbf{h})\triangleq \frac{1}{2}{\Vert \mathbf{y}-\mathbf{Wh}\Vert }_2^2+\phi (\mathbf{h}),$$

$${\ell }_2(\mathbf{y},\mathbf{W})\triangleq \underset{\mathbf{h}\in \mathcal{H},\mathbf{r}\in \mathcal{R}}{\min }{\tilde{\ell }}_2(\mathbf{y},\mathbf{W},\mathbf{h},\mathbf{r}),{\tilde{\ell }}_2(\mathbf{y},\mathbf{W},\mathbf{h},\mathbf{r})\triangleq \frac{1}{2}{\Vert \mathbf{y}-\mathbf{Wh}-\mathbf{r}\Vert }_2^2+\phi (\mathbf{h})+\varphi (\mathbf{r}),$$

$$\underset{\mathbf{W}\in \mathcal{C}}{\min }\left[f(\mathbf{W})\triangleq g(\mathbf{W})+\psi (\mathbf{W})\right],g(\mathbf{W})\triangleq \frac{1}{n}\sum _{i=1}^n\ell ({\mathbf{y}}_i,\mathbf{W}),$$

$$\underset{\mathbf{W}\in \mathcal{C}}{\min }\left[\overline{f}(\mathbf{W})\triangleq {\mathbb{E}}_{\mathbf{y}\sim \nu }[\ell (\mathbf{y},\mathbf{W})]+\psi (\mathbf{W})\right],$$

$${\mathbf{h}}^{+}$$

$${\mathbf{r}}^{+}$$

$${\nabla }_{\mathbf{W}}{\ell }_2(\mathbf{y},\mathbf{W})=\left(\mathbf{W}{\mathbf{h}}^{\ast }(\mathbf{y},\mathbf{W})+{\mathbf{r}}^{\ast }(\mathbf{y},\mathbf{W})-\mathbf{y}\right){\mathbf{h}}^{\ast }(\mathbf{y},\mathbf{W}{)}^T.$$

$${\mathbb{E}}_{{\mathcal{B}}_{s,t}}[\Vert {\mathbf{V}}_{s,t}-\nabla g({\mathbf{W}}^{s,t}){\Vert }^2|{\mathcal{F}}_{s,t}]\le \alpha (n,b)L^2\Vert {\mathbf{W}}^{s,t}-{\mathbf{W}}^{s,0}{\Vert }^2,\forall s\ge 0,\forall t\in (m].$$

$$\mathbb{E}[\Vert {\Gamma }_{f^{\prime },\eta }({\mathbf{W}}_{\mathrm{final}}){\Vert }^2]\le \left(\frac{f({\mathbf{W}}^{0,0})-f^{\ast }}{\eta \left(1/2-\eta L\right)}\right)\frac{1}{mS},\text{mbox}where{f}^{\ast }\triangleq \underset{\mathbf{W}\in \mathcal{C}}{\min }f(\mathbf{W}).$$

$$f^{\prime }({\mathbf{W}}^{s,t+1})$$

$$+\left(L/2-1/(2\eta )\right)\Vert {\mathbf{W}}^{s,t+1}-{\mathbf{W}}^{s,t}{\Vert }^2+\left(L-1/(2\eta )\right)\Vert {\tilde{\mathbf{W}}}^{s,t+1}-{\mathbf{W}}^{s,t}{\Vert }^2.$$

$${\mathbb{E}}_{{\mathcal{B}}_{s,t}}[{\hat{f}}_{s,t+1}({\mathbf{W}}^{s,t+1})|{\mathcal{F}}_{s,t}]\le {\hat{f}}_{s,t}({\mathbf{W}}^{s,t})+\eta \left(\eta L-1/2\right)\Vert {\Gamma }_{f^{\prime },\eta }({\mathbf{W}}^{s,t}){\Vert }^2.$$

$${\mathbb{E}}_{{\mathcal{B}}_{s,(m]}}[f^{\prime }({\mathbf{W}}^{s+1,0})|{\mathcal{F}}_{s,0}]\le f^{\prime }({\mathbf{W}}^{s,0})+\eta \left(\eta L-1/2\right){\sum }_{t=0}^{m-1}{\mathbb{E}}_{{\mathcal{B}}_{s,(t]}}[\Vert {\Gamma }_{f^{\prime },\eta }({\mathbf{W}}^{s,t}){\Vert }^2|{\mathcal{F}}_{s,0}].$$

$$f^{\ast }\le f({\mathbf{W}}^{0,0})+\eta \left(\eta L-1/2\right){\sum }_{s=0}^{S-1}{\sum }_{t=0}^{m-1}\mathbb{E}[\Vert {\Gamma }_{f^{\prime },\eta }({\mathbf{W}}^{s,t}){\Vert }^2]+\eta {\sum }_{s=0}^{S-1}{\sum }_{t=0}^{m-1}{\Vert {\mathbf{E}}_{s,t}\Vert }^2.$$

$${\mathbf{W}}^{s,t+1}={\mathbf{prox}}_{\eta {\psi }^{\prime }}({\mathbf{W}}^{s,t}-\eta {\mathbf{V}}_{s,t}).$$

$${\tilde{\mathbf{W}}}^{s,t+1}\triangleq {\mathbf{prox}}_{\eta {\psi }^{\prime }}({\mathbf{W}}^{s,t}-\eta \nabla g({\mathbf{W}}^{s,t})),$$

$${\Gamma }_{f^{\prime },\eta }({\mathbf{W}}^{s,t})=\frac{1}{\eta }({\mathbf{W}}^{s,t}-{\tilde{\mathbf{W}}}^{s,t+1}).$$

$$f^{\prime }({\mathbf{W}}^{s,t+1})$$

$$+\left(\frac{L}{2}+\frac{1}{2\eta }\right)\Vert {\mathbf{Z}}_1-{\mathbf{W}}^{s,t}{\Vert }^2-\frac{1}{2\eta }\Vert {\mathbf{W}}^{s,t+1}-{\mathbf{Z}}_1{\Vert }^2,$$

$$f^{\prime }({\tilde{\mathbf{W}}}^{s,t+1})$$

$$-\frac{1}{2\eta }\Vert {\tilde{\mathbf{W}}}^{s,t+1}-{\mathbf{Z}}_2{\Vert }^2.$$

$$f^{\prime }({\mathbf{W}}^{s,t+1})$$

$$+\left(L-\frac{1}{2\eta }\right)\Vert {\tilde{\mathbf{W}}}^{s,t+1}-{\mathbf{W}}^{s,t}{\Vert }^2-\frac{1}{2\eta }\Vert {\mathbf{W}}^{s,t+1}-{\tilde{\mathbf{W}}}^{s,t+1}{\Vert }^2.$$

$$⟨{\mathbf{W}}^{s,t+1}-{\tilde{\mathbf{W}}}^{s,t+1},\nabla g({\mathbf{W}}^{s,t})-{\mathbf{V}}_{s,t}⟩$$

$${\mathbb{E}}_{{\mathcal{B}}_{s,t}}\left[⟨{\mathbf{W}}^{s,t+1}-{\tilde{\mathbf{W}}}^{s,t+1},\nabla g({\mathbf{W}}^{s,t})-{\mathbf{V}}_{s,t}⟩|{\mathcal{F}}_{s,t}\right]$$

$$+\frac{\eta L^2}{2}\alpha (n,b)\Vert {\mathbf{W}}^{s,t}-{\mathbf{W}}^{s,0}{\Vert }^2.$$

$${\mathbb{E}}_{{\mathcal{B}}_{s,t}}[f^{\prime }({\mathbf{W}}^{s,t+1})|{\mathcal{F}}_{s,t}]$$

$$+\left(\frac{L}{2}-\frac{1}{2\eta }\right){\mathbb{E}}_{{\mathcal{B}}_{s,t}}[\Vert {\mathbf{W}}^{s,t+1}-{\mathbf{W}}^{s,t}{\Vert }^2|{\mathcal{F}}_{s,t}]+\left(L-\frac{1}{2\eta }\right)\Vert {\tilde{\mathbf{W}}}^{s,t+1}-{\mathbf{W}}^{s,t}{\Vert }^2.$$

$${\zeta }_t=(1+1/m){\zeta }_{t+1}+\frac{\eta L^2}{2}\alpha (n,b),{\zeta }_m=0.$$

$${\hat{f}}_{s,0}({\mathbf{W}}^{s,0})=f^{\prime }({\mathbf{W}}^{s,0})\text{mbox}and{\hat{f}}_{s,m}({\mathbf{W}}^{s,m})=f^{\prime }({\mathbf{W}}^{s,m}).$$

$${\zeta }_0=\frac{\left((1+1/m{)}^m-1\right)m}{2}\alpha (n,b)\eta L^2\le \frac{\left(e-1\right)m}{2}\alpha (n,b)\eta L^2\le m\alpha (n,b)\eta L^2.$$

$${\hat{f}}_{s,t+1}({\mathbf{W}}^{s,t+1})\le f^{\prime }({\mathbf{W}}^{s,t+1})+(1+m){\zeta }_{t+1}\Vert {\mathbf{W}}^{s,t+1}-{\mathbf{W}}^{s,t}{\Vert }^2+(1+1/m){\zeta }_{t+1}\Vert {\mathbf{W}}^{s,0}-{\mathbf{W}}^{s,t}{\Vert }^2.$$

$${\mathbb{E}}_{{\mathcal{B}}_{s,t}}[{\hat{f}}_{s,t+1}({\mathbf{W}}^{s,t+1})|{\mathcal{F}}_{s,t}]$$

$$+\left(\frac{L}{2}-\frac{1}{2\eta }+(1+m){\zeta }_{t+1}\right){\mathbb{E}}_{{\mathcal{B}}_{s,t}}[\Vert {\mathbf{W}}^{s,t+1}-{\mathbf{W}}^{s,t}{\Vert }^2|{\mathcal{F}}_{s,t}]+\left(L-\frac{1}{2\eta }\right)\Vert {\tilde{\mathbf{W}}}^{s,t+1}-{\mathbf{W}}^{s,t}{\Vert }^2$$

$$={\hat{f}}_{s,t}({\mathbf{W}}^{s,t})+\left(\frac{L}{2}-\frac{1}{2\eta }+(1+m){\zeta }_{t+1}\right){\mathbb{E}}_{{\mathcal{B}}_{s,t}}[\Vert {\mathbf{W}}^{s,t+1}-{\mathbf{W}}^{s,t}{\Vert }^2|{\mathcal{F}}_{s,t}]$$

Full text

1 Introduction

Matrix factorization plays an important role in machine learning, due to its wide applications including collaborative filtering [1], parts-based learning [2] and clustering [3]. Given a data matrix $\mathbf{Y}\!\in\!\mathbb{R}^{d\times n}$ (where $n$ denotes the number of data samples and $d$ denotes the dimension of each data sample), one aims to find matrices $\mathbf{W}\!\in\!\mathbb{R}^{d\times k}$ and $\mathbf{H}\!\in\!\mathbb{R}^{k\times n}$ such that $\mathbf{Y}\!\approx\!\mathbf{W}\mathbf{H}$, where usually $k\!\ll\!\min(d,n)$. Many existing algorithms in the literature find $\mathbf{W}$ and $\mathbf{H}$ by minimizing $\|\mathbf{Y}-\mathbf{W}\mathbf{H}\|_{F}^{2}$ (possibly with regularizations on $\mathbf{W}$ or $\mathbf{H}$ to induce structural properties). However, when the number of data samples $n$ becomes large, solving this problem using batch algorithms can be highly inefficient, in terms of both computation and storage. To improve the efficiency, online (stochastic) matrix factorization algorithms [4, 5, 6, 7, 8, 9, 10, 11, 12] have been proposed to learn the dictionary (matrix) $\mathbf{W}$ from a sequence of randomly drawn data samples (possibly in mini-batches).111We use “online” and “stochastic” interchangeably. In the literature, online algorithms also cover the setting where the number of data samples $n$ is infinite. However, we only consider $n$ to be finite, but can be very large. Extensive numerical evidence has shown that the stochastic matrix factorization algorithms can learn $\mathbf{W}$ significantly faster than their batch counterparts, in terms of (empirical) convergence rate. In addition, only $O(dk)$ storage space are consumed by these algorithms, in contrast to $O(ndk)$ by batch methods. Due to these advantages, stochastic matrix factorization algorithms have gained much popularity in recent years.

From an optimization point of view, previous stochastic matrix factorization algorithms [4, 5, 6, 7, 8, 9, 10, 11, 12] aim to solve a nonconvex stochastic program involving the dictionary $\mathbf{W}$ (see (4) in Remark 1). They mainly leverage two optimization frameworks, namely stochastic majorization-minimization (SMM) [13] and stochastic (sub-)gradient descent (SGD) [14]. Based on either framework, almost sure asymptotic convergence to stationary points of the stochastic program has been shown. However, the asymptotic nature of this convergence analysis cannot provide insights into the dictionary learning process at any finite time instant. Thus, we wish to understand the non-asymptotic convergence of dictionary learning, as well as the sample complexity for learning a “reasonable” dictionary $\mathbf{W}$. Besides, we desire improved stochastic methods that yield faster convergence, at least for large but finite datasets.

Main Contributions. Inspired by the recent advances in variance reduction techniques (particularly for nonconvex problems) [15, 16, 17], we propose a unified framework to speed up the existing stochastic matrix factorization algorithms via variance reduction (VR). As shown in Section 2, our framework is general and it subsumes eight well-known SMF formulations as special instances, including those robust variants that explicitly model outliers in their objective functions. Within this framework, we derive a non-asymptotic convergence rate, measured in terms of how “fast” the algorithm converges to a stationary point in expectation (see Section 4.1 for precise statements). Accordingly, we derive sample and computational complexities for our algorithm to converge to an $\epsilon$-stationary point in expectation. To the best of our knowledge, this is the first time that a non-asymptotic convergence analysis (together with sample complexity) is established for SMF algorithms. To further improve the efficiency of our framework, we also provide an inexact convergence analysis of our framework, where in each step the coefficient (and outlier) vectors are only learned approximately. In the context of stochastic (variance-reduced) nonconvex proximal gradient algorithms, this is the first inexact analysis thus far. We show, via extensive experiments on various SMF formulations and datasets, that our variance-reduced framework consistently outperforms state-of-the-art frameworks (including SMM and SGD), in terms of both convergence rate and accuracy of learned dictionaries.

Related works. Another line of works [18, 19, 20] considers using SMF algorithms in a matrix completion setting, where the matrices $\mathbf{W}$ and $\mathbf{H}$ are jointly learned from a random sequence of entries (or columns) of a data matrix $\mathbf{Y}$. By leveraging low-rankness and the restricted isometry property (RIP) on $\mathbf{Y}$, the product of $\mathbf{W}$ and $\mathbf{H}$ can exactly recover $\mathbf{Y}$ (with high probability). However, their problem setting differs from ours in two aspects. First, we do not store or explicitly learn $\mathbf{H}$, due to space constraints. Second, our convergence results do not rely on low-rankness or RIP. Therefore, similar to [4, 5, 6, 7, 8, 9, 10, 11, 12], our results are stated in terms of stationary points, which are oftentimes appealing in real applications [4, 11]. In the context of stochastic PCA, VR techniques have been applied to yield a linear convergence [21]. Again, the problem setting therein vastly differs from ours.

Notations. We use boldface capital letters, boldface lowercase letters and lowercase letters to denote matrices, column vectors and scalars respectively. Denote $\mathbb{R}_{+}$ as $[0,\infty)$ and $\overline{\mathbb{R}}_{+}$ as $[0,\infty]$. For any $n\geq 1$, define $[n]\triangleq\{1,\ldots,n\}$ and $(n]\triangleq\{0,\ldots,n\}$. For a matrix $\mathbf{X}$, denote its $j$-th column as $\mathbf{X}_{:j}$. We use $\left\|\cdot\right\|$ to denote the $\ell_{2}$-norm of a vector and the Frobenius norm of a matrix. For a convex and closed set $\mathcal{X}$, denote its indicator function as $\delta_{\mathcal{X}}$, such that $\delta_{\mathcal{X}}(\mathbf{x})=0$ if $\mathbf{x}\in\mathcal{X}$ and $\infty$ otherwise. All the sections, lemmas and algorithms beginning with ‘S’ will appear in the supplemental material.

2 Problem Formulation

Let $\{\mathbf{y}_{i}\}_{i=1}^{n}\subseteq\mathbb{R}^{d}$ denote the set of data samples. We first define a loss function $\ell_{1}$ w.r.t. a data sample $\mathbf{y}$ and a dictionary $\mathbf{W}\in\mathbb{R}^{d\times k}$ as

[TABLE]

where $\mathcal{H}$ and $\varphi:\mathbb{R}^{k}\to\overline{\mathbb{R}}_{+}$ denote the constraint set and regularizer of the coefficient vector $\mathbf{h}$ respectively. Since large-sale datasets may contain outliers, we can define a “robust” version of $\ell_{1}$ as

[TABLE]

where $\mathcal{R}$ and $\phi:\mathbb{R}^{d}\to\overline{\mathbb{R}}_{+}$ denote the constraint set and regularizer of the outlier vector $\mathbf{r}$ respectively. For convenience, let $\ell$ denote either $\ell_{1}$ or $\ell_{2}$. Based on $\ell$, we formulate our problem as

[TABLE]

where $\mathcal{C}$ and $\psi:\mathbb{R}^{d\times k}\to\overline{\mathbb{R}}_{+}$ denote the constraint set and regularizer of the dictionary $\mathbf{W}$ respectively. Our targeted problem (3) is general and flexible, in the sense that by choosing the constraint sets $\mathcal{C}$, $\mathcal{H}$ and $\mathcal{R}$ and the regularizers $\psi$, $\varphi$ and $\phi$ in different manners, (3) encompasses many important examples in the literature of SMF. For loss function $\ell_{1}$, we have:

1. (P1)

Online DL (ODL) [4]: $\mathcal{C}\triangleq\{\mathbf{W}\in\mathbb{R}^{d\times k}:\left\|\mathbf{W}_{:i}\right\|\leq 1,\forall i\in[n]\}$, $\mathcal{H}\triangleq\mathbb{R}^{k}$, $\psi\equiv 0$ and $\varphi\triangleq\lambda\left\|\cdot\right\|_{1}$ ($\lambda>0$). Some other variants and extensions were also discussed in [4, Section 5]. 2. (P2)

Online Structured Sparse Learning (OSSL) [5]: $\mathcal{C}\triangleq\mathbb{R}^{d\times k}$, $\mathcal{H}\triangleq\mathbb{R}^{k}$, $\psi\equiv 0$ and $\varphi(\mathbf{h})\triangleq\sum_{j\in[n^{\prime}]}\lambda_{j}\|\underline{\mathbf{h}}_{j}\|$, where $\underline{\mathbf{h}}_{j}$ is a subvector of $\mathbf{h}$, $\lambda_{j}\geq 0$ and $n^{\prime}\in[n]$. 3. (P3)

Online NMF (ONMF) [6]: $\mathcal{C}\triangleq\{\mathbf{W}\in\mathbb{R}_{+}^{d\times k}\!:\!\left\|\mathbf{W}_{:i}\right\|_{1}=1,\forall i\!\in\![n]\}$, $\mathcal{H}=\mathbb{R}_{+}^{k}$, $\psi\equiv 0$ and $\varphi\triangleq\frac{\lambda}{2}\left\|\cdot\right\|^{2}$ ($\lambda>0$). See [6, Section III] for some extensions. 4. (P4)

Smooth Sparse ODL (SSODL) [10]: $\mathcal{C}\triangleq\{\mathbf{W}\in\mathbb{R}_{+}^{d\times k}\!:\!\left\|\mathbf{W}_{:i}\right\|_{1}=1,\forall i\!\in\![n]\}$, $\mathcal{H}=\mathbb{R}_{+}^{k}$, $\psi(\mathbf{W})\triangleq\frac{\lambda_{1}}{2}\operatorname{tr}(\mathbf{W}^{T}\mathbf{L}\mathbf{W})$ and $\varphi\triangleq\frac{\lambda_{2}}{2}\left\|\cdot\right\|^{2}$, where $\lambda_{1},\lambda_{2}>0$ and $\mathbf{L}$ is positive semidefinite.

For loss function $\ell_{2}$, we have:

1. (P5)

Online Robust PCA (ORPCA) [7]: $\mathcal{C}\triangleq\mathbb{R}^{d\times k}$, $\mathcal{H}\triangleq\mathbb{R}^{k}$, $\mathcal{R}\triangleq\mathbb{R}^{d}$, $\psi\triangleq\frac{\lambda_{1}}{2n}\left\|\cdot\right\|^{2}$, $\varphi\triangleq\frac{\lambda_{1}}{2}\left\|\cdot\right\|^{2}$ and $\phi\triangleq\lambda_{2}\left\|\cdot\right\|_{1}$, where $\lambda_{1},\lambda_{2}>0$. 2. (P6)

Online Max-norm Regularized Matrix Decomposition (OMRMD) [8]: $\mathcal{C}\triangleq\mathbb{R}^{d\times k}$, $\mathcal{H}\triangleq\{\mathbf{h}\in\mathbb{R}^{k}:\left\|\mathbf{h}\right\|\leq 1\}$, $\mathcal{R}\triangleq\mathbb{R}^{d}$, $\psi\triangleq\frac{\lambda_{1}}{2n}\left\|\cdot\right\|_{2,\infty}^{2}$, $\varphi\equiv 0$ and $\phi=\lambda_{2}\left\|\cdot\right\|_{1}$, where $\lambda_{1},\lambda_{2}>0$. 3. (P7)

Online Low-Rank Representation (OLRR) [9]: $\mathcal{C}\triangleq\mathbb{R}^{d\times k}$, $\mathcal{H}\triangleq\mathbb{R}^{k}$, $\mathcal{R}\triangleq\mathbb{R}^{d}$, $\varphi\triangleq\frac{\lambda_{1}}{2}\left\|\cdot\right\|^{2}$, $\phi\triangleq\lambda_{2}\left\|\cdot\right\|_{1}$ and $\psi(\mathbf{W})\triangleq\frac{1}{2n^{2}}(\|\mathbf{W}^{T}\mathbf{A}\mathbf{Y}\|^{2}+\lambda_{3}\|\mathbf{B}\mathbf{W}\|^{2})$, where $\mathbf{Y}$, $\mathbf{A}$ and $\mathbf{B}$ are given (constant) matrices, and $\lambda_{1},\lambda_{2},\lambda_{3}>0$. 4. (P8)

Online Robust NMF (ORNMF) [11]: $\mathcal{C}\triangleq\{\mathbf{W}\in\mathbb{R}_{+}^{d\times k}\!:\!\left\|\mathbf{W}_{:i}\right\|\leq 1,\forall i\!\in\![n]\}$, $\mathcal{H}\triangleq[0,M]^{k}$, $\mathcal{R}\triangleq[-M^{\prime},M^{\prime}]^{d}$, $\psi\equiv 0$, $\varphi\equiv 0$ and $\phi\triangleq\lambda\left\|\cdot\right\|_{1}$, where $M,M^{\prime},\lambda>0$.

Remark 1*.*

In almost all the works cited in this section, the optimization problem is formulated as a stochastic program that minimizes the expected risk of dictionary training. It has the form

[TABLE]

where $\nu$ is a probability distribution. The reasons that we consider (3) rather than (4) are three-fold. First, (3) can be regarded as the sample average approximation [22] of (4). Solving this approximation has been a popular and efficient approach for solving a general stochastic program like (4) [23]. Second, (3) can also be interpreted as minimizing the empirical risk incurred by any large-scale finite training dataset. Solving (3) efficiently can greatly speed up the dictionary learning process. Last but not least, since (3) is a finite-sum minimization problem, we can employ the recently developed variance reduction techniques [15, 16] to solve it in an efficient manner. See Section 3 for details.

For algorithm and convergence analysis, we will focus on the loss function $\ell_{2}$, which contains $\ell_{1}$ as a special case (by choosing $\mathcal{R}=\{\mathbf{0}\}$).

3 Algorithm

We present the pseudo-code for our stochastic matrix factorization algorithm with variance reduction techniques in Algorithm 1. To develop our algorithm, we first make two mild assumptions. We then describe our algorithm in details, with focus on learning the coefficient and outlier vectors.

Assumptions.  The following two assumptions are satisfied by all the Problems (P1) to (P8).

1. (A1)

$\psi$, $\varphi$ and $\phi$ are convex, proper and closed functions on $\mathbb{R}^{d\times k}$, $\mathbb{R}^{k}$ and $\mathbb{R}^{d}$ respectively, with proximal operators that can be evaluated efficiently.222For a function $\gamma:\mathbb{R}^{d}\to\mathbb{R}$, the “efficient” evaluation of $\mathbf{prox}_{\gamma}$ means that for any $\mathbf{x}\in\mathbb{R}^{d}$, $\mathbf{prox}_{\gamma}(\mathbf{x})$ can be computed in $O(d)$ arithemetic operations. The same definitions apply to “efficient projections” in (A2). 2. (A2)

$\mathcal{C}$, $\mathcal{H}$ and $\mathcal{R}$ are convex and admit efficient (Euclidean) projections onto them.

Algorithm Description. Variance-reduced stochastic optimization algorithms typically contain outer and inner loops [15, 24]. In Algorithm 1, we use $s$ and $t$ to denote the indices for outer and inner iterations respectively. At the beginning of each outer iteration $s$, we solve $n$ (regularized) least-square regression problems to learn the set of coefficient and outlier vectors $\{(\mathbf{h}^{s,0}_{j},\mathbf{r}^{s,0}_{j})\}_{j=1}^{n}$ based on $\mathbf{W}^{s,0}$. Based on the learned vectors, we compute $\mathbf{G}_{s,0}$, which will be shown to be $\nabla g(\mathbf{W}^{s,0})$ in Section 4. In each inner iteration $t$, we solve a (random) mini-batch of regression problems, based on $\mathbf{W}^{s,0}$ and $\mathbf{W}^{s,t}$ respectively. The regression results, together with $\mathbf{G}_{s,0}$, are used to calculate a search direction $\mathbf{V}_{s,t}$, which acts as a variance-reduced gradient estimate of the loss function $g$ at $\mathbf{W}^{s,t}$. The proximal operator of $\eta\psi+\delta_{\mathcal{C}}$ either has a closed form [25] or can be obtained using Douglas-Rachford splitting or alternating direction method of multipliers [26], both in $O(dk)$ (arithmetic) operations. We choose the final dictionary via option I or II (in lines 15 and 16 respectively).333In this work, the line numbers always refer to those in Algorithm 1. In practice, we choose the final dictionary as the last iterate $\mathbf{W}^{s-1,m}$, i.e., option II. However, to streamline the analysis, we output the dictionary using option I. This is consistent with the many works on variance-reduced SGD algorithms, e.g., [15].

Learning coefficient and outlier vectors. It can be seen that the computational complexity of Algorithm 1 crucially depends on solving the problem (2), which appears in lines 4, 8 and 9 of Algorithm 1. To solve this problem, we leverage the block successive upperbound minimization framework [27], by alternating between the following two steps, i.e.,

[TABLE]

where $\mathbf{h}^{+}$ and $\mathbf{r}^{+}$ denote the updated values of $\mathbf{h}$ and $\mathbf{r}$ respectively (in one iteration). The proximal operators of $\varphi/L^{\prime}+\delta_{\mathcal{H}}$ and $\phi+\delta_{\mathcal{R}}$ can be obtained in the same way as that of $\eta\psi+\delta_{\mathcal{C}}$. Overall, the computational complexity of solving (2) is $O(dk)$.

Remark 2*.*

In our implementation, we do not store $\{\mathbf{h}^{s,0}_{j},\mathbf{r}^{s,0}_{j}\}_{j=1}^{n}$ before computing $\mathbf{G}_{s,0}$. Instead, after obtaining $\mathbf{h}^{s,0}_{j}$ and $\mathbf{r}^{s,0}_{j}$, we add them to the running sum in line 5, and then discard them. The same procedure applies to computing $\mathbf{V}_{s,t}$. Thus the storage complexity of Algorithm 1 is $O(dk)$. This complexity is linear in the dimension of the optimization variable $\mathbf{W}$, which has $dk$ entries.

4 Convergence Analysis

Additional assumptions. For analysis purposes, we make the following three additional assumptions:

1. (A3)

The data samples $\{\mathbf{y}_{i}\}_{i=1}^{n}$ lie in a compact set $\mathcal{Y}\subseteq\mathbb{R}^{d}$. 2. (A4)

For any $\mathbf{y}\in\mathcal{Y}$ and $\mathbf{W}\in\mathcal{C}$, $\widetilde{\ell}_{2}(\mathbf{y},\mathbf{W},\cdot,\cdot)$ is jointly strongly convex on $\mathcal{H}\times\mathcal{R}$. 3. (A5)

$\mathcal{C}$ is compact; $\mathcal{H}$ is compact or $\varphi$ is coercive on $\mathbb{R}^{k}$; $\mathcal{R}$ is compact or $\varphi$ is coercive on $\mathbb{R}^{d}$.

Remark 3*.*

Assumption (A3) is natural since all real data are bounded. Assumption (A4) is common in the literature [4, 11] and can be satisfied by adding Tikhonov regularizers (see [11, Remark 4] for a detailed discussion). Assumption (A5) is satisfied by all the Problems (P1) to (P8), except in (P2), (P5), (P6) and (P7), $\mathcal{C}$ is not bounded. However, for each problem, it was proven that the sequence of basis matrices $\{\mathbf{W}_{i}\}_{i\geq 0}$ generated by the corresponding algorithm in [7], [8] or [9] belonged to a compact set almost surely. Thus the compactness of $\mathcal{C}$ is a reasonable assumption.

4.1 Key Lemmas and Main Results

We first present some definitions that are used in the statements of our results. We then state our key lemmas (Lemmas 1 and 2) and main convergence theorems (Theorems 1 and 2), together with the derived sample and computational complexities (Corollary 1). The proof sketches of all the results are shown in Section 4.2, with complete proofs deferred to Section S-1 to Section S-4.

Definitions. Define a filtration $\{\mathcal{F}_{s,t}\}_{s\in(S-1],t\in[m]}$ such that $\mathcal{F}_{s,t}$ is the $\sigma$-algebra generated by the random sets $\{\mathcal{B}_{s^{\prime},t^{\prime}}\}_{s^{\prime}\in(s-1],t^{\prime}\in[m]}\cup\{\mathcal{B}_{s,t^{\prime}}\}_{t^{\prime}\in[t-1]}$. It is easy to see that $\{\mathbf{W}^{s,t}\}_{s\in(S-1],t\in[m]}$ is adapted to to $\{\mathcal{F}_{s,t}\}_{s\in(S-1],t\in[m]}$. Consider a composite function $h\!\triangleq\!h_{1}\!+\!h_{2}$ defined on a Euclidean space $\mathcal{E}$, where $h_{1}\!:\!\mathcal{E}\!\to\!\mathbb{R}$ is differentiable and $h_{2}\!:\!\mathcal{E}\!\to\!\overline{\mathbb{R}}_{+}$ is proper, closed and convex. For any $\eta\!>\!0$ and $\mathbf{x}\!\in\!\mathcal{E}$, define the (proximal) gradient mapping [28] of $h$ at $\mathbf{x}$ with step size $\eta$, $\Gamma_{h,\eta}(\mathbf{x})\!\triangleq\!(1/\eta)(\mathbf{x}\!-\!\mathbf{prox}_{\eta h_{2}}(\mathbf{x}\!-\!\eta\nabla h_{1}(\mathbf{x})))$. A point $\mathbf{y}\in\mathcal{E}$ is $\epsilon$-stationary w.r.t. the problem $\min_{\mathbf{x}\in\mathcal{E}}\!h(\mathbf{x})$ and $\eta$ if $\|\Gamma_{h,\eta}(\mathbf{y})\|^{2}\leq\epsilon$. Finally, define $\alpha(n,b)\triangleq(n-b)/(b(n-1))$.

Lemma 1** (Smoothness of $\ell_{2}$).**

Define $(\mathbf{h}^{*}(\mathbf{y},\mathbf{W}),\mathbf{r}^{*}(\mathbf{y},\mathbf{W}))\triangleq\operatorname*{arg\,min}_{\mathbf{h}\in\mathcal{H},\mathbf{r}\in\mathcal{R}}\widetilde{\ell}_{2}(\mathbf{y},\mathbf{W},\mathbf{h},\mathbf{r})$. The loss function $\ell_{2}(\cdot,\cdot)$ in (2) is differentiable on $\mathbb{R}^{d}\times\mathbb{R}^{d\times k}$ and

[TABLE]

In addition, $(\mathbf{y},\mathbf{W})\mapsto\nabla_{\mathbf{W}}\ell_{2}(\mathbf{y},\mathbf{W})$ is continuous on $\mathcal{Y}\times\mathcal{C}$. Moreover, for any $\mathbf{y}\in\mathcal{Y}$, $\mathbf{W}\mapsto\nabla_{\mathbf{W}}\ell_{2}(\mathbf{y},\mathbf{W})$ is Lipschitz on $\mathcal{C}$ with a Lipschitz constant $L$ independent of $\mathbf{y}$.

Lemma 2** (Variance bound of $\mathbf{V}_{s,t}$).**

In Algorithm 1, $\mathbb{E}_{\mathcal{B}_{s,t}}[\mathbf{V}_{s,t}|\mathcal{F}_{s,t}]=\nabla g(\mathbf{W}^{s,t})$ and

[TABLE]

Measuring convergence rate. Since the loss function $\ell$ in (3) is nonconvex (due to bilinearity), obtaining global minima of $f$ is in general out-of-reach. Thus in almost all the previous works [4, 5, 6, 7, 8, 9, 10, 11, 12], convergence to stationary points of (3) was studied instead. Define $f^{\prime}\triangleq f+\delta_{\mathcal{C}}$. Following [29, 16], for the sequence of (random) dictionaries $\{\mathbf{W}^{s,t}\}_{s\in(S-1],t\in[m]}$ generated in Algorithm 1, we propose to use $\mathbb{E}[\|\Gamma_{f^{\prime},\eta}(\mathbf{W}^{s,t})\|^{2}]$ to measure its convergence rate (to stationary points).To see the validity of this measure, define the directional derivative of $f$ at $\mathbf{W}\in\mathcal{C}$ along $\mathbf{D}$, i.e., $\overline{\nabla}f(\mathbf{W};\mathbf{D})\triangleq\lim_{\delta\downarrow 0}(f(\mathbf{W}+\delta\mathbf{D})-f(\mathbf{W}))/{\delta}$. In previous works [4, 5, 6, 7, 8, 9, 10, 11, 12], $\mathbf{W}^{\prime}\in\mathcal{C}$ is characterized as a stationary point of (3) if and only if for any $\mathbf{W}\in\mathcal{C}$, $\overline{\nabla}f(\mathbf{W}^{\prime};\mathbf{W}-\mathbf{W}^{\prime})\geq 0$. It can be shown that this condition is equivalent to $\Gamma_{f^{\prime},\eta}(\mathbf{W}^{\prime})=\mathbf{0}$, for any $\eta>0$. In particular, if $\psi\equiv 0$, we can verify that both conditions are equivalent to the condition that $-\nabla g(\mathbf{W}^{\prime})$ lies in the normal cone of $\mathcal{C}$ at $\mathbf{W}^{\prime}$.

Equipped with the proper convergence measure, we now state our main convergence theorem.

Theorem 1**.**

Assume that 2 to (A5) hold. In Algorithm 1, if we choose $\eta=1/(\theta L)$ for some $\theta>2$ and $m$ and $b$ to satisfy $\theta(\theta-1)\geq 2m(m+1)\alpha(n,b)$, then

[TABLE]

Note that $mS$ denotes the total number of inner iterations. Thus (9) states a sublinear convergence rate of $\{\mathbf{W}^{s,t}\}_{s\geq 0,t\in[m]}$ towards the set of stationary points of (3). In addition, since $\alpha(n,b)\leq 1/b$, we can choose $\theta$, $m$ and $b$ such that $m\!\leq\!\!\sqrt{b\theta(\theta\!-\!1)/2}\!-\!1\,(\theta\!>\!2)$ to satisfy the conditions in Theorem 1.

In practice, for efficiency considerations, it is important that we learn the coefficient and outlier vectors (in lines 4, 8 and 9) approximately, by only finding inexact solutions to the subproblem (2). The errors in these inexact solutions will be propagated to the variance-reduced gradient $\mathbf{V}_{s,t}$ as a whole. In the next theorem, we show that if the (additive) error $\mathbf{E}_{s,t}$ in $\mathbf{V}_{s,t}$ is properly controlled (i.e., $\|\mathbf{E}_{s,t}\|$ decreases at a certain rate), then we can achieve the same convergence rate as in Theorem 1. Note that although inexact analyses have been done for batch [30] and incremental [31] proximal gradient (PG) algorithms, our result is the first one for a stochastic (variance-reduced) PG algorithm.

Theorem 2** (Inexact Analysis).**

Assume that 2 to (A5) hold. In Algorithm 1, if $\eta=1/(\theta L)$ $(\theta\!>\!2)$, $\theta(\theta-1)\!\geq\!4m(m+1)\alpha(n,b)$ and $\|\mathbf{E}_{s,t}\|\!=(ms+t)^{-(1/2+\tau)}$ for any $\tau>0$, then $\mathbb{E}[\|\Gamma_{f^{\prime},\eta}(\mathbf{W}_{\rm final})\|^{2}]\!=\!O(1/(mS))$.

Remark 4* (Bound $\|\mathbf{E}_{s,t}\|$).*

Denote an approximate solution of $\operatorname*{arg\,min}_{\mathbf{h}\in\mathcal{H},\mathbf{r}\in\mathcal{R}}\widetilde{\ell}_{2}(\mathbf{y}_{j},\mathbf{W}^{s,t},\mathbf{h},\mathbf{r})$ as $(\widehat{\mathbf{h}}^{s,t}_{j},\widehat{\mathbf{r}}^{s,t}_{j})$. From the definition of $\mathbf{V}_{s,t}$ and the compactness of $\mathcal{H}$ and $\mathcal{R}$, we see that $\|\mathbf{V}_{s,t}\|$ can be bounded by the learning errors in coefficient and outlier vectors, i.e., $\{(\|\mathbf{h}^{s,t}_{j}-\widehat{\mathbf{h}}^{s,t}_{j}\|,\|\mathbf{r}^{s,t}_{j}-\widehat{\mathbf{r}}^{s,t}_{j}\|)\}_{j\in\mathcal{B}_{s,t}}$. For each $j\!\in\![n]$, these learning errors can be further bounded by the (infimum of) norm of (sub-)gradient of $\widetilde{\ell}_{2}(\mathbf{y}_{j},\mathbf{W}^{s,t},\cdot,\cdot)$ at $(\widehat{\mathbf{h}}^{s,t}_{j},\widehat{\mathbf{r}}^{s,t}_{j})$, due to the first-order optimality conditions. Therefore, we are able to bound $\|\mathbf{E}_{s,t}\|$ in terms of the approximate solutions $\{(\widehat{\mathbf{h}}^{s,t}_{j},\widehat{\mathbf{r}}^{s,t}_{j})\}_{j=1}^{n}$.

Finally, based on the convergence rates in Theorems 1 and 2, we are able to derive the sample and computational complexities for Algorithm 1 to attain an $\epsilon$-stationary point (in expectation).

Corollary 1** (Sample and Computational Complexities).**

*For any $\epsilon\!>\!0$ and $\theta\!>\!2$, the sample and computational complexities for $\mathbf{W}_{\rm final}$ in Algorithm 1 to be an $\epsilon$-stationary point of (3) in expectation (i.e., $\mathbb{E}[\|\Gamma_{f^{\prime},\eta}(\mathbf{W}_{\rm final})\|^{2}]\leq\epsilon$) are $O(n^{2/3}/\epsilon)$ and $O(n^{2/3}dk/\epsilon)$ respectively. *

4.2 Proof Sketches

Proof Sketch of Lemma 1. By Assumption (A5), if $\mathcal{H}$ is not compact, then $\varphi$ is coercive. This implies the boundedness of $\mathbf{h}^{*}(\mathbf{y},\mathbf{W})$. A similar argument also applies to $\mathbf{r}^{*}(\mathbf{y},\!\mathbf{W})$. Thus it is equivalent to minimizing $\widetilde{\ell}_{2}$ over a compact set $\mathcal{H}^{\prime}\!\times\!\mathcal{R}^{\prime}\subseteq\mathcal{H}\!\times\!\mathcal{R}$. In addition, Assumptions (A3) and (A5) ensure the compactness of $\mathcal{Y}$ and $\mathcal{C}$ respectively. Since Assumptions (A2) and (A4) ensure the uniqueness of $(\mathbf{h}^{*}(\mathbf{y},\mathbf{W}),\mathbf{r}^{*}(\mathbf{y},\mathbf{W}))$, for any $(\mathbf{y},\mathbf{W})\in\mathbb{R}^{d}\times\mathbb{R}^{d\times k}$, we can invoke Danskin’s theorem (see Lemma S-1) to guarantee the differentiability of $\ell_{2}$, and compute $\nabla_{\mathbf{W}}\ell_{2}(\mathbf{y},\mathbf{W})$ as in (7). Additionally, we can invoke the Maximum theorem (see Lemma S-2) to ensure the continuity of $(\mathbf{y},\mathbf{W})\!\mapsto\!(\mathbf{h}^{*}(\mathbf{y},\mathbf{W}),\mathbf{r}^{*}(\mathbf{y},\mathbf{W}))$ on $\mathcal{Y}\!\times\!\mathcal{C}$. This implies the continuity of $(\mathbf{y},\mathbf{W})\!\mapsto\!\nabla_{\mathbf{W}}\ell_{2}(\mathbf{y},\mathbf{W})$. Based on this, we can assert the Lipschitz continuity of both $\mathbf{h}^{*}(\cdot,\cdot)$ and $\mathbf{r}^{*}(\cdot,\cdot)$ on $\mathcal{Y}\!\times\!\mathcal{C}$. The proceeding arguments, together with the compactness of $\mathcal{Y}$, $\mathcal{C}$, $\mathcal{H}^{\prime}$ and $\mathcal{R}^{\prime}$, allow us to conclude the Lipschitz continuity of $\mathbf{W}\!\mapsto\!\nabla_{\mathbf{W}}\ell_{2}(\mathbf{y},\mathbf{W})$. ∎

Proof Sketch of Lemma 2. To show (8), we first apply Lemma S-4 to $\|\mathbf{V}_{s,t}-\nabla g(\mathbf{W}^{s,t})\|^{2}$ and then make use of the $L$-Lipschitz continuity of $\mathbf{W}\!\mapsto\!\nabla_{\mathbf{W}}\ell_{2}(\mathbf{y},\mathbf{W})$ in Lemma 1. ∎

Proof Sketch of Theorem 1. From Lemma 1, we know the loss function $g$ is differentiable. Define $\psi^{\prime}\!\triangleq\!\psi+\delta_{\mathcal{C}}$ and $\widetilde{\mathbf{W}}^{s,t+1}\!\triangleq\!\mathbf{prox}_{\eta\psi^{\prime}}(\mathbf{W}^{s,t}\!-\!\eta\nabla g(\mathbf{W}^{s,t}))$, then $\Gamma_{f^{\prime},\eta}(\mathbf{W}^{s,t})\!=\!({1}/{\eta})(\mathbf{W}^{s,t}\!-\!\widetilde{\mathbf{W}}^{s,t+1})$. Using Lemma S-5 and the Lipschitz continuity of $\nabla g$ (see Lemma 1), we have the recursion

[TABLE]

Conditioning on $\mathcal{F}_{s,t}$, we then take expectations w.r.t. $\mathcal{B}_{s,t}$ on both sides of (10). Using Lemma 2, we can bound $\mathbb{E}_{\mathcal{B}_{s,t}}[\langle\mathbf{W}^{s,t+1}\!-\!\widetilde{\mathbf{W}}^{s,t+1},\!\nabla g(\mathbf{W}^{s,t})\!-\!\mathbf{V}_{s,t}\rangle|\mathcal{F}_{s,t}]$ in terms of $\|\mathbf{W}^{s,t}-\mathbf{W}^{s,0}\|^{2}$ and $\mathbb{E}_{\mathcal{B}_{s,t}}[\|\mathbf{W}^{s,t+1}\!-\!\widetilde{\mathbf{W}}^{s,t+1}\|^{2}|\mathcal{F}_{s,t}]$. Then, instead of directly telescoping (10), we construct a surrogate function $\widehat{f}_{s,t}(\mathbf{W})\triangleq f^{\prime}(\mathbf{W})+\zeta_{t}\|\mathbf{W}-\mathbf{W}^{s,0}\|^{2}$, where the sequence $\{\zeta_{t}\}_{t=0}^{m}$ is given by the recursion $\zeta_{t}=(1+1/m)\zeta_{t+1}+(\eta L^{2}/2)\alpha(n,b)$, where $\zeta_{m}=0$. Using this recursive relation, together with the conditions $\eta=1/(\theta L)$ and $\theta(\theta-1)\geq 2(m+1)\alpha(n,b)$ in Theorem 1, we manage to arrive at

[TABLE]

Define $\mathcal{B}_{s,(t]}\triangleq\{\mathcal{B}_{s,j}\}_{j=0}^{t}$. We now telescope (11) over $t=0,1,\ldots,m-1$ to obtain

[TABLE]

Finally, we telescope (12) over $s\!\in\!(S\!-\!1]$ and use option I to choose the final dictionary $\mathbf{W}_{\rm final}$.∎

Proof Sketch of Theorem 2. The proof of Theorem 2 is modified from that of Theorem 1. Due to the error $\mathbf{E}_{s,t}$,​ $\mathbf{V}_{s,t}$ in (10) is replaced by $\widehat{\mathbf{V}}_{s,t}\!\triangleq\!\mathbf{V}_{s,t}\!+\!\mathbf{E}_{s,t}$. We decouple $\mathbf{E}_{s,t}$ from $\mathbf{V}_{s,t}$ through $\langle\mathbf{W}^{s,t+1}-\widetilde{\mathbf{W}}^{s,t+1},\nabla g(\mathbf{W}^{s,t})-\widehat{\mathbf{V}}_{s,t}\rangle\!\leq\!\frac{1}{2\eta}\|\mathbf{W}^{s,t+1}-\widetilde{\mathbf{W}}^{s,t+1}\|^{2}+{\eta}\|{\nabla g(\mathbf{W}^{s,t})-\mathbf{V}_{s,t}}\|^{2}+\eta\|\mathbf{E}_{s,t}\|^{2}$. We then take expectation on both sides and carefully follow the telescoping procedure (over $t$ and $s$) used in proving Theorem 1. After some algebraic manipulations, we arrive at

[TABLE]

Since $\|\mathbf{E}_{s,t}\|=(ms+t)^{-(1/2+\tau)}\,(\tau\!>\!0)$, for any $S\geq 1$, $\sum_{s=0}^{S-1}\sum_{t=0}^{m-1}\|\mathbf{E}_{s,t}\|^{2}<\infty$. Therefore by using option I to choose $\mathbf{W}_{\rm final}$, we complete the proof. ∎

Proof of Corollary 1. By Theorem 1, we know that $O(1/\epsilon)$ inner iterations are needed for $\mathbb{E}[\|\Gamma_{\!f^{\prime}\!,\eta}(\mathbf{W}_{\rm final})\|^{2}]\!\leq\!\epsilon$. By evenly distributing the $n$ data samples drawn at the start of each outer iteration into $m$ inner iterations, we see that each inner iteration takes in $O(n/m+b)$ samples. Thus the sample complexity is $O((n/m+b)/\epsilon)$. Based on it, we can also derive the computational complexity to be $O((n/m+b)dk/\epsilon)$. This is because both subproblem (2) (in lines 4, 8 and 9) and the proximal operator of $\eta\psi+\delta_{\mathcal{C}}$ (in line 7) can be solved or evaluated in $O(dk)$ arithmetic operations (see Section 3). From the condition $\theta(\theta-1)\geq 2m(m+1)\alpha(n,b)$ in Theorem 1, we see that $m=O(\sqrt{nb/(n-b)})=O(\sqrt{b})$. Thus the sample and computational complexities become $O((n/\sqrt{b}+b)/\epsilon)$ and $O((n/\sqrt{b}+b)dk/\epsilon)$ respectively. To optimize both costs, it is necessary that $b=\Theta(n/\sqrt{b})$, which amounts to $b=\Theta(n^{2/3})$. (Note that this also implies $m=\Theta(n^{1/3})$.) ∎

5 Numerical Experiments

5.1 Tested SMF Formulations and Datasets

We tested the performance of our proposed framework (Algorithm 1) on four representative SMF formulations in Section 2, including ODL, ONMF, ORPCA and ORNMF. Among them, ODL and ONMF are non-robust SMF formulations whereas ORPCA and ORNMF are robust ones, i.e., they explicitly model outliers. Therefore, for ODL and ONMF, we tested their algorithms on the CBCL [2] and MNIST [32] datasets; while for ORPCA and ORNMF, we tested their algorithms on the Synth and YaleB [33] datasets. The CBCL and MNIST datasets are commonly used in testing (non-robust) matrix factorization algorithms, e.g., [2]. The YaleB dataset consists of face images taken under various lighting conditions. The shadow and gloss in these images caused by imperfect illumination can be regarded as outliers. The Synth dataset was generated synthetically in a similar fashion to those in [7, 8, 11]. Specifically, we first generated a $d\!\times\!k^{\prime}$ matrix $\overline{\bf W}$ and a $k^{\prime}\!\times\!n$ matrix $\overline{\bf H}$, where $d\!=\!400$, $n\!=\!1\!\times\!10^{5}$ and $k^{\prime}\!=\!10$. The entries of $\overline{\bf W}$ and $\overline{\bf H}$ were drawn i.i.d. from a normal distribution with mean 0.5 and variance $1/\!\sqrt{k^{\prime}}$. We then generated a $d\!\times\!n$ outlier matrix $\overline{\bf R}$ by first uniformly randomly selecting $\lfloor{(1\!-\!\rho_{s})dn}\rfloor$ of its entries to be zero, where $\rho_{s}$ denotes the outlier density. The remaining entries were then generated i.i.d. from a uniform distribution with support $[-M,\!M]$. We set $\rho_{s}\!=\!0.1$ and $M\!=\!1000$. Finally $\mathbf{V}\!=\!\overline{\bf W}\,\overline{\bf H}\!+\!\overline{\bf R}$.

5.2 Benchmarking Frameworks and Choice of Parameters

For each tested SMF formulation, we compared our variance-reduced SMF framework (denoted as VR) against another two optimization frameworks commonly used in previous works on SMF, namely SMM and SGD. For completeness, the pseudo-codes for these two methods are shown in Algorithms S-2 and S-3 respectively. We next describe the parameter setting in our method. From our analysis in Corollary 1, we set the mini-batch size $b=c_{1}n^{2/3}$ and the number of inner iterations $m=c_{2}n^{1/3}$, where $c_{1}=0.2$ and $c_{2}=0.5$. We set the number of outer iterations $S\!=\!10$ and the parameter $\tau$ in Theorem 2 to $1\!\times\!10^{-3}$. For each tested SMF formulation on each dataset, we chose the step size $\eta$ such that our method yielded the best (empirical) convergence rate. The plots of objective values resulting from other step sizes are shown in Figure S-3. For the latent dimension $k$, following the convention in [2, 4, 11], we set it to 49. This parameter can also be chosen from domain knowledge or a Bayesian approach, e.g., [34]. Since in practice we found the performance of all the three frameworks (VR, SMM and SGD) was not sensitive to $k$, we fixed it for simplicity. For both SMM and SGD, we used the same mini-batch size $b$ as in VR. The step sizes $\{\gamma_{t}\}_{t\geq 0}$ in SGD are chosen in the classical way [14], i.e., $\gamma_{t}=\beta/(bt+\beta^{\prime})$. Similar to VR, we chose $\beta,\beta^{\prime}>0$ such that SGD achieved the best empirical convergence rate. For simplicity, we initialized the dictionary $\mathbf{W}$ such that all its entries are equal to one. Finally, following [4, 6, 7, 11], we set the regularization weight $\lambda$ in ODL, ONMF and ORNMF to $1/\sqrt{d}$. In addition, we set $\lambda_{1}=\lambda_{2}=1/\sqrt{d}$ in ORPCA.

5.3 Plots of Objective Values

For each SMF formulation and each dataset, since we focus on convergence to a stationary point, we first run a batch (deterministic) gradient-based matrix factorization algorithm (e.g., [35]) to estimate a stationary point $\widehat{\mathbf{W}}$. This resembles the way to use batch methods to estimate a global optimum in the literature of stochastic convex optimization [15, 24]. Based on $\widehat{\mathbf{W}}$, we plot the log-suboptimality of the objective value, i.e., $\log(f(\mathbf{W})-f(\widehat{\mathbf{W}}))$ versus both running time and number of data passes, for VR, SMM and SGD respectively (see Figure 1d).444All the algorithms are implemented in Matlab® 8.6 and run on a machine with a 3.6-GHz processor. We plot objective values versus number of data passes so that the results in Figure 1d are agnostic to the actual implementation of the algorithms. From Figure 1d, in terms of both time and number of data passes, we observe that for all the SMF algorithms and datasets, our variance-reduced framework (VR) not only converges faster than SMM and SGD, but also find a more accurate approximate of the stationary point $\widehat{\mathbf{W}}$. Explanations and discussions of these observations are deferred to Section S-6.

5.4 Subspace Recovery by ORPCA

We also considered the subspace recovery task using ORPCA [7] on the Synth dataset. The ground-truth subspace $\mathcal{U}$ is given by the column space of the (ground-truth) dictionary $\overline{\bf W}$ in generating the Synth dataset (see Section 5.1). Accordingly, the estimated subspace $\mathcal{U}^{\prime}$ at any time instant is given by the column space of our learned dictionary. The similarity between $\mathcal{U}^{\prime}$ and $\mathcal{U}$ is measured by the expressed variance (EV) (see the definition in [36]). A larger value of EV indicates a higher similarity, hence a better subspace recovery result. (A unit EV indicates perfect recovery, i.e., $\mathcal{U}^{\prime}\!=\!\mathcal{U}$.) We employed VR, SMM and SGD on the ORPCA problem to recover the subspace $\mathcal{U}$. In addition to the original Synth dataset with outlier density $\rho_{s}\!=\!0.1$, we generated another one with $\rho_{s}\!=\!0.3$. From the results in Figure 2d, we observe that VR consistently outperforms the other two frameworks. Specifically, compared to SMM and SGD, VR not only recovers $\mathcal{U}$ faster, but also recovers a subspace that is closer to $\mathcal{U}$ (in terms of EV). These observations are indeed consistent with those in Section 5.3, and are explained in Section S-6.

S-1 Proof of Theorem 1

For notational convenience, define $\psi^{\prime}\triangleq\psi+\delta_{\mathcal{C}}$, then

[TABLE]

In addition, define

[TABLE]

then

[TABLE]

Applying Lemma S-5 to both (S-1) and (S-2), we have that for any $\mathbf{Z}_{1},\mathbf{Z}_{2}\in\mathbb{R}^{d\times k}$,

[TABLE]

By first setting $\mathbf{Z}_{1}=\widetilde{\mathbf{W}}^{s,t+1}$ in (S-4) and $\mathbf{Z}_{2}=\mathbf{W}^{s,t}$ in (S-5) and then summing both inequalities, we have

[TABLE]

Using Lemma S-6, we have

[TABLE]

Therefore, by Lemma 2, we have

[TABLE]

Taking expectations on both sides of (S-6) and making use of (S-8), we have

[TABLE]

Now, define a surrogate function $\widehat{f}_{s,t}(\mathbf{W})\triangleq f^{\prime}(\mathbf{W})+\zeta_{t}\|\mathbf{W}-\mathbf{W}^{s,0}\|^{2}$, where the sequence $\{\zeta_{t}\}_{t=0}^{m}$ is given by the recursion

[TABLE]

In particular, we have

[TABLE]

From the recursion (S-10), we also observe that $\{\zeta_{t}\}_{t=0}^{m}$ is deceasing and

[TABLE]

Again by using Lemma S-6, we bound

[TABLE]

Combining (S-9) and (S-13), we have

[TABLE]

Using the conditions $\eta=1/(\theta L)$ and $\theta(\theta-1)\geq 2m(m+1)\alpha(n,b)$ and (S-12), we have for any $t\in[m]$,

[TABLE]

Therefore, by (S-3), we have

[TABLE]

Telescoping (S-20) over $t=0,\ldots,m-1$ and noting (S-11), we have

[TABLE]

For any $S\geq 0$, telescope (S-21) over $s=0,1,\ldots,S-1$ and we have

[TABLE]

If we choose the final basis matrix $\mathbf{W}_{\rm final}$ uniformly randomly from $\{\mathbf{W}^{s,t}\}_{s\in(S-1],t\in(m]}$, then

[TABLE]

where (a) holds because $\mathbf{W}^{0,0}\in\mathcal{C}$.

S-2 Proof of Theorem 2

Define $\widehat{\mathbf{V}}_{s,t}\triangleq\mathbf{E}_{s,t}+\mathbf{V}_{s,t}$, then

[TABLE]

Applying Lemma S-5 to (S-24), we have that for any $\mathbf{Z}_{3}\in\mathbb{R}^{d\times k}$,

[TABLE]

By setting $\mathbf{Z}_{3}=\widetilde{\mathbf{W}}^{s,t+1}$ in (S-25) and then summing (S-25) and (S-5), we have

[TABLE]

By repeatedly applying Lemma S-6, we have

[TABLE]

Therefore, by Lemma 2, we have

[TABLE]

Taking expectations on both sides of (S-26) and making use of (S-27), we have

[TABLE]

Now, define a surrogate function $\widehat{f}_{s,t}(\mathbf{W})\triangleq f^{\prime}(\mathbf{W})+\zeta_{t}\|\mathbf{W}-\mathbf{W}^{s,0}\|^{2}$, where the sequence $\{\zeta_{t}\}_{t=0}^{m}$ is given by the recursion

[TABLE]

In particular, we have

[TABLE]

From the recursion (S-29), we also observe that $\{\zeta_{t}\}_{t=0}^{m}$ is deceasing and

[TABLE]

Again by using Lemma S-6, we bound

[TABLE]

Combining (S-28) and (S-32), we have

[TABLE]

Using the conditions $\eta=1/(\theta L)$ and $\theta(\theta-1)\geq 4m(m+1)\alpha(n,b)$ and (S-12), we have for any $t\in[m]$,

[TABLE]

Therefore, by (S-3), we have

[TABLE]

Telescoping (S-20) over $t=0,\ldots,m-1$ and noting (S-30), we have

[TABLE]

For any $S\geq 0$, telescope (S-40) over $s=0,1,\ldots,S-1$ and we have

[TABLE]

Since $\|\mathbf{E}_{s,t}\|=o(1/\!\sqrt{ms+t})$, for any $S\geq 1$, there exists a constant $E<\infty$ (independent of $S$) such that $\sum_{s=0}^{S-1}\sum_{t=0}^{m-1}\|\mathbf{E}_{s,t}\|^{2}\leq E$. Therefore by using option I to choose $\mathbf{W}_{\rm final}$, we have

[TABLE]

S-3 Proof of Lemma 1

Define $(\mathbf{h}^{*}(\mathbf{y},\mathbf{W}),\mathbf{r}^{*}(\mathbf{y},\mathbf{W}))\triangleq\min_{\mathbf{h}\in\mathcal{H},\mathbf{r}\in\mathcal{R}}\widetilde{\ell}_{2}(\mathbf{y},\mathbf{W},\mathbf{h},\mathbf{r})$. By assumption (A5), if $\mathcal{H}$ is not compact, then $\varphi$ will be coercive. This implies the boundedness of $\mathbf{h}^{*}(\mathbf{y},\mathbf{W})$. A similar argument also applies to $\mathbf{r}^{*}\!(\mathbf{y},\!\mathbf{W})$. Thus it is equivalent to consider minimizing $\widetilde{\ell}_{2}$ over some compact sets $\mathcal{H}^{\prime}\subseteq\mathcal{H}$ and $\mathcal{R}^{\prime}\subseteq\mathcal{R}$.

It is easy to verify that the following conditions hold:

$\widetilde{\ell}_{2}(\cdot,\cdot,\mathbf{h},\mathbf{r})$ is differentiable on $\mathbb{R}^{d}\times\mathbb{R}^{d\times k}$, for each $(\mathbf{h},\mathbf{r})\in\mathbb{R}^{k}\times\mathbb{R}^{d}$, 2. 2.

$\widetilde{\ell}_{2}(\cdot,\cdot,\cdot,\cdot)$ and $(\mathbf{y},\mathbf{W},\mathbf{h},\mathbf{r})\mapsto\nabla_{(\mathbf{y},\mathbf{W})}\widetilde{\ell}_{2}(\mathbf{y},\mathbf{W},\mathbf{h},\mathbf{r})$ are continuous on $\mathbb{R}^{d}\times\mathbb{R}^{d\times k}\times\mathbb{R}^{k}\times\mathbb{R}^{d}$, 3. 3.

For any $\mathbf{y}\in\mathbb{R}^{d}$ and $\mathbf{W}\!\in\!\mathbb{R}^{d\times k}$, the minimizer $(\mathbf{h}^{*}(\mathbf{y},\mathbf{W}),\mathbf{r}^{*}(\mathbf{y},\mathbf{W}))$ is unique, due to assumption (A4).

Thus, we can invoke Danskin’s theorem (see Lemma S-1) to conclude that $\ell_{2}(\cdot,\cdot)$ is differentiable on $\mathbb{R}^{d}\times\mathbb{R}^{d\times k}$, and compute

[TABLE]

Furthermore, we can show $\mathbf{h}^{*}(\mathbf{y},\mathbf{W})$ and $\mathbf{r}^{*}(\mathbf{y},\mathbf{W})$ are both continuous on $\mathcal{Y}\times\mathcal{C}$ by the maximum theorem (see Lemma S-2), since the conditions in this theorem are trivially satisfied in our case. Therefore, from (S-43), we can see that $(\mathbf{y},\mathbf{W})\mapsto\nabla_{\mathbf{W}}\ell_{2}(\mathbf{y},\mathbf{W})$ is continuous on $\mathcal{Y}\times\mathcal{C}$.

We next show that for all $\mathbf{y}\in\mathcal{Y}$, both $\mathbf{h}^{*}(\mathbf{y},\cdot)$ and $\mathbf{r}^{*}(\mathbf{y},\cdot)$ are Lipschitz on $\mathcal{C}$, with Lipschitz constants independent of $\mathbf{y}$. Fix any $\mathbf{y}_{1},\mathbf{y}_{2}\in\mathcal{Y}$ and $\mathbf{W}_{1},\mathbf{W}_{2}\in\mathcal{C}$. Define

[TABLE]

where $\mathbf{Y}_{i}=[\mathbf{W}_{i}\;\;\mathbf{I}]$, $i=1,2$ and $\mathbf{b}(\mathbf{h},\mathbf{r})=[\mathbf{h}^{T}\;\;\mathbf{r}^{T}]^{T}$. By Lemma S-3, we have for all $(\mathbf{h}_{1},\mathbf{r}_{1})$ and $(\mathbf{h}_{2},\mathbf{r}_{2})$ in $\mathcal{H}^{\prime}\times\mathcal{R}$, there exist constants $c_{1},c_{2}>0$ (independent of $\mathbf{y}_{1}$, $\mathbf{y}_{2}$, $\mathbf{Y}_{1}$ and $\mathbf{Y}_{2}$) such that

[TABLE]

In particular, we have

[TABLE]

where $\mathbf{h}^{*}_{i}=\mathbf{h}^{*}(\mathbf{y}_{i},\mathbf{W}_{i})$ and $\mathbf{r}^{*}_{i}=\mathbf{r}^{*}(\mathbf{y}_{i},\mathbf{W}_{i})$, $i=1,2$. On the other hand, by Assumption (A4), there exist a constant $\nu>0$ such that

[TABLE]

Combining (S-45) and (S-46), we have

[TABLE]

where $c^{\prime}_{i}=c_{i}/\nu$, $i=1,2$. Since

[TABLE]

we have

[TABLE]

This indeed shows both $\mathbf{h}^{*}(\cdot,\cdot)$ and $\mathbf{r}^{*}(\cdot,\cdot)$ are Lipschitz on $\mathcal{Y}\times\mathcal{C}$.

To show the Lipschitz continuity of $\mathbf{W}\mapsto\nabla_{\mathbf{W}}\ell_{2}(\mathbf{y},\mathbf{W})$ on $\mathcal{C}$, consider any $\mathbf{W}_{1},\mathbf{W}_{2}\in\mathcal{C}$ and any $\mathbf{y}\in\mathcal{Y}$. From (S-43), we have

[TABLE]

Define $c\triangleq 2\max(c^{\prime}_{1},c^{\prime}_{2})$. We bound each term on the RHS of (S-51) as follows

[TABLE]

[TABLE]

[TABLE]

Using compactness of $\mathcal{Y}$, $\mathcal{C}$, $\mathcal{H}^{\prime}$ and $\mathcal{R}^{\prime}$, we complete the proof.

S-4 Proof of Lemma 2

From Lemma 1, we know

[TABLE]

where

[TABLE]

Thus

[TABLE]

where in (a) we use (7) in Lemma 1, in (b) we use Lemma S-4 and in (c) we use the Lipschitz continuity of $\mathbf{W}\!\mapsto\!\nabla_{\mathbf{W}}\ell_{2}(\mathbf{y},\mathbf{W})$ in Lemma 1.

S-5 Algorithms for SMF via SMM and SGD

For completeness, we present the algorithms to solve SMF problems via SMM and (proximal) SGD in Algorithms S-2 and S-3 respectively.

S-6 Additional Experimental Results and Discussions

Additional experimental results are shown in Figure S-3. Next, we provide some intuitions about the observations in Section 5.3. As shown in [37], under mild conditions, the SMM method can be regarded as a spacial case of the SGD method (with diminishing step sizes). Since the step size vanishes asymptotically, both SMM and SGD will make minute progress in learning the dictionary $\mathbf{W}$ after a large number of iterations (or equivalently, data samples). Therefore, they fail to find highly accurate stationary points. In contrast, our variance-reduced method employs a constant step size, therefore it continues to make non-negligible progress asymptotically, so as to reach a much higher accuracy. Similar arguments can also explain the results in Section 5.4.

S-7 Technical Lemmas

Lemma S-1** (Danskin’s Theorem [38]).**

Let $\mathcal{X}$ be a metric space and $\mathcal{U}$ be a normed vector space. Let $f:\mathcal{X}\times\mathcal{U}\to\mathbb{R}$ have the following properties

$f(x,\cdot)$* is differentiable on $\mathcal{U}$, for any $x\in\mathcal{X}$,* 2. 2.

$f(x,u)$* and $\nabla_{u}f(x,u)$ are continuous on $\mathcal{X}\times\mathcal{U}$.*

Let $\Phi$ be a compact set in $\mathcal{X}$. Define $v(u)\triangleq\inf_{x\in\Phi}f(x,u)$ and $\mathcal{S}(u)\triangleq\operatorname*{arg\,min}_{x\in\Phi}f(x,u)$. If for some $u_{0}\in\mathcal{U}$, $\mathcal{S}(u_{0})=\{x_{0}\}$, then $v$ is (Hadamard) differentiable at $u=u_{0}$ and $\nabla v(u_{0})=\nabla_{u}f(x_{0},u_{0})$.

Lemma S-2** (The Maximum Theorem [39, Theorem 14.2.1]).**

Let $\mathcal{P}$ and $\mathcal{X}$ be two metric spaces. Consider a maximization problem

[TABLE]

where $\mathcal{B}:\mathcal{P}\rightrightarrows\mathcal{X}$ is a set-valued map and $f:\mathcal{P}\times\mathcal{X}\to\mathbb{R}$ is a function. Define the set-valued map $\mathcal{S}(p)\triangleq\operatorname*{arg\,max}_{x\in\mathcal{B}(p)}f(p,x)$. If $\mathcal{B}$ is compact-valued and continuous on $\mathcal{P}$ and $f$ is continuous on $\mathcal{P}\times\mathcal{X}$, and if for some $p_{0}\in\mathcal{P}$, $\mathcal{S}(p_{0})=\{s(p_{0})\}$, where $s:\mathcal{P}\to\mathcal{X}$ is a function, then $s$ is continuous at $p=p_{0}$. Moreover, we have the same conclusions if the maximization in (S-57) is replaced by minimization.

Lemma S-3**.**

Let $\mathbf{z}$, $\mathbf{z}^{\prime}\in\mathcal{Z}\subseteq\mathbb{R}^{m}$ and $\mathbf{A}$, $\mathbf{A}^{\prime}\in\mathcal{A}\subseteq\mathbb{R}^{m\times n}$, where both $\mathcal{Z}$ and $\mathcal{A}$ are compact sets. Let $\mathcal{B}$ be a compact set in $\mathbb{R}^{n}$, and define $g:\mathcal{B}\to\mathbb{R}$ as $g(\mathbf{b})=1/2\left\|\mathbf{z}-\mathbf{A}\mathbf{b}\right\|_{2}^{2}-1/2\left\|\mathbf{z}^{\prime}-\mathbf{A}^{\prime}\mathbf{b}\right\|^{2}_{2}$. Then $g$ is Lipschitz on $\mathcal{B}$ with Lipschitz constant $c_{1}\left\|\mathbf{z}-\mathbf{z}^{\prime}\right\|_{2}+c_{2}\left\|\mathbf{A}-\mathbf{A}^{\prime}\right\|_{2}$, where $c_{1}$ and $c_{2}$ are two positive constants. In particular, when both $\mathbf{z}^{\prime}$ and $\mathbf{A}^{\prime}$ are zero, we have that $\widetilde{g}(\mathbf{b})=1/2\left\|\mathbf{z}-\mathbf{A}\mathbf{b}\right\|_{2}^{2}$ is Lipschitz on $\mathcal{B}$ with Lipschitz constant $c$ independent of $\mathbf{z}$ and $\mathbf{A}$.

Proof.

It suffices to show $\left\|\nabla g(\mathbf{b})\right\|_{2}\leq c_{1}\left\|\mathbf{z}-\mathbf{z}^{\prime}\right\|_{2}+c_{2}\left\|\mathbf{A}-\mathbf{A}^{\prime}\right\|_{2}$ for any $\mathbf{b}\in\mathcal{B}$ and some positive constants $c_{1}$ and $c_{2}$ (independent of $\mathbf{b}$). We write $\left\|\nabla g(\mathbf{b})\right\|_{2}$ as

[TABLE]

By the compactness of $\mathcal{Z}$, $\mathcal{A}$ and $\mathcal{B}$, there exist positive constants $M_{1}$, $M_{2}$ and $M_{3}$ such that $\left\|\mathbf{z}\right\|_{2}\leq M_{1}$, $\left\|\mathbf{A}\right\|_{2}\leq M_{2}$ and $\left\|\mathbf{b}\right\|_{2}\leq M_{3}$, for any $\mathbf{z}\in\mathcal{Z}$, $\mathbf{A}\in\mathcal{A}$ and $\mathbf{b}\in\mathcal{B}$. Thus,

[TABLE]

Similarly for $\left\langle 2\right\rangle$ we have

[TABLE]

Hence $\left\langle 1\right\rangle+\left\langle 2\right\rangle\leq M_{2}\left\|\mathbf{z}-\mathbf{z}^{\prime}\right\|_{2}+(M_{1}+2M_{2}M_{3})\left\|\mathbf{A}-\mathbf{A}^{\prime}\right\|_{2}$. We now take $c_{1}=M_{2}$ and $c_{2}=M_{1}+2M_{2}M_{3}$ to complete the proof. ∎

Lemma S-4** ([40]).**

Let $\{\mathbf{z}_{i}\}_{i=1}^{n}\subseteq\mathbb{R}^{d}$ and define $\overline{\bf z}\triangleq 1/n\sum_{i=1}^{n}\mathbf{z}_{i}$. Uniformly sample a random subset $\mathcal{S}$ of $[n]$ with size $b$ without replacement. Then

[TABLE]

Lemma S-5** (Upperbound of a composite function; [29, 16]).**

Consider a composite function $h\triangleq h_{1}+h_{2}$ defined on a Euclidean space $\mathcal{E}$, where both $h_{1}:\mathcal{E}\to\mathbb{R}$ is differentiable with $L$-Lipschitz gradient and $h_{2}:\mathcal{E}\to\overline{\mathbb{R}}_{+}$ is proper, closed and convex. Fix any $\eta>0$ and $\mathbf{x}\in\mathcal{E}$. Let $y=\mathbf{prox}_{\eta h_{2}}(\mathbf{x}-\eta\mathbf{d})$, for some $\mathbf{d}\in\mathcal{E}$. Then

[TABLE]

Lemma S-6** ([41]).**

For any two vectors $\mathbf{a},\mathbf{b}$ in a Euclidean space $\mathcal{E}$ and any $\beta>0$, we have

[TABLE]

Proof.

$(\sqrt{{1}/{\beta}}\mathbf{a}-\sqrt{\beta}\mathbf{b})^{2}\geq 0$. ∎