# Revisiting Co-Occurring Directions: Sharper Analysis and Efficient Algorithm
for Sparse Matrices

Luo Luo Cheng Chen Guangzeng Xie Haishan Ye  Department of Mathematics, The
Hong Kong University of Science and Technology<EMAIL_ADDRESS>Department of
Computer Science and Engineering, Shanghai Jiao Tong University;
<EMAIL_ADDRESS>Academy for Advanced Interdisciplinary Studies,
Peking University<EMAIL_ADDRESS>Shenzhen Research Institute of Big Data,
The Chinese University of Hong Kong, Shenzhen<EMAIL_ADDRESS>
###### Abstract

We study the streaming model for approximate matrix multiplication (AMM). We
are interested in the scenario that the algorithm can only take one pass over
the data with limited memory. The state-of-the-art deterministic sketching
algorithm for streaming AMM is the co-occurring directions (COD), which has
much smaller approximation errors than randomized algorithms and outperforms
other deterministic sketching methods empirically. In this paper, we provide a
tighter error bound for COD whose leading term considers the potential
approximate low-rank structure and the correlation of input matrices. We prove
COD is space optimal with respect to our improved error bound. We also propose
a variant of COD for sparse matrices with theoretical guarantees. The
experiments on real-world sparse datasets show that the proposed algorithm is
more efficient than baseline methods.

##  1 Introduction

A large scale machine learning system usually receives data sequentially and
it is often impossible to exactly store the entire data set. Thus, the
approximate matrix multiplication (AMM) in the streaming fashion is an
important and fundamental task for scientific computation and big data
analysis. For example, the product of matrices from multi-modal datasets
captures the correlation between different modalities. In addition, many
classical algorithms including canonical correlation analysis [14],
generalized eigenvector decomposition [9], partial least squares [34],
spectral co-clustering [6] require to perform approximate matrix
multiplication when the data set is very large. On the other hand, data
matrices from real-world are usually low-rank and sparse, which motivated us
to design efficient and effective sparse algorithms.

This paper considers streaming AMM problem as follows. Give two large matrices
${\bf X}\in{\mathbb{R}}^{n\times d_{x}}$ and ${\bf Y}\in{\mathbb{R}}^{n\times
d_{y}}$, we are interested in finding a low-rank estimator ${\bf A}^{\top}{\bf
B}$ to approximate ${\bf X}^{\top}{\bf Y}$, where ${\bf
A}\in{\mathbb{R}}^{m\times d_{x}}$, ${\bf B}\in{\mathbb{R}}^{m\times d_{y}}$
and $m$ is much smaller than $n$, $d_{x}$ and $d_{y}$. We focus on the row
update model, that is, the algorithm receives rows of ${\bf X}$ and ${\bf Y}$
sequentially and it only takes one pass over input matrices with limited
memory. The key challenge for this problem is to reduce the space/time
complexity while maintaining the approximation error.

Inspired by the idea of finding frequent items [24], Liberty [19] proposed
frequent directions algorithm (FD), which considers the symmetric case of AMM
such that ${\bf X}={\bf Y}$ (a.k.a., the covariance sketching). FD achieves
optimal tradeoffs between space cost and approximation error [36, 12, 10].
Moreover, we can combine FD with subspace power iteration [35, 26] to design
an algorithm which is efficient for sparse matrix multiplication [11], called
sparse frequent directions (SFD). Recently, Huang [15] integrated random
sampling [7] into FD to reduce its time complexity. Luo et al. [22] introduced
a regularization term for FD, which makes the estimator is more friendly to
inverse operation. FD technique can also be used to accelerate many popular
machine learning models, such as convex online optimization [20, 22],
factorization machine [21], linear contextual bandits [18, 3] and ridge
regression [32].

Mroueh et al. [25] proposed a variant of FD called co-occurring directions
(COD) for streaming AMM. COD shrinks the singular values of input matrices
${\bf X}$ and ${\bf Y}$ simultaneously at each iteration. It is shown that COD
has significantly better performance than other sketching algorithms [37, 7,
4, 31] on AMM problem empirically. However, the existing spectral error bound
of COD can not completely explain its high performance. It depends on the
Frobenius norm of ${\bf X}$ and ${\bf Y}$, which ignores the potential low-
rank structure of the data matrix. Specifically, in the case of ${\bf X}={\bf
Y}$, the procedure of COD degrades to FD, but its error bound is worse than
that of FD. Another deterministic sketching method for AMM, which we call FD-
AMM [37], directly adopts FD to sketch the concatenated matrix ${\bf Z}=[{\bf
X},{\bf Y}]$. The output of the algorithm is an approximation of ${\bf
Z}^{\top}{\bf Z}$, whose sub-matrix corresponds to an estimator of ${\bf
X}^{\top}{\bf Y}$.

In this paper, we provide a sharper analysis for co-occurring directions
(COD). We give a new spectral norm error bound which considers the potential
low-rank structure of the target matrix. Our bound could be much tighter than
Mroueh et al.’s (2016) results when the spectrum of the exact matrix product
is dominated by its top singular values. In addition, we prove that the space
complexity of COD is optimal to attain our improved error bound. Furthermore,
in the case of ${\bf X}={\bf Y}$, our result matches the error bound of FD.

We further propose sparse co-occurring directions (SCOD) algorithm and prove
its error bound is similar to standard COD while its running time mainly
depends on the non-zero entries of input matrices. We conduct numerical
experiments on cross-language datasets to show that SCOD has better
performance than state-of-the-art algorithms empirically.

The rest of the paper is organized as follows. In Section 2, define the
notation used in this paper and introduce the background of related algorithms
for streaming AMM. In Section 3, we provide our new error bound for COD
algorithm and show the corresponding space lower bound. In Section 5, we
propose SCOD and give its theoretical guarantees. In Section 6, we conduct the
numerical experiments to show the superiority of SCOD. We defer detailed proof
of some lemmas and theorems into appendix. We conclude our work in Section 7.

##  2 Notations and Preliminaries

In this section, we first introduce the notation will be used in this paper.
Then we give the backgrounds of frequent directions and related algorithms for
AMM.

###  2.1 Notations

We let ${\bf I}_{p}$ be the $p\times p$ identity matrix and ${\bf 0}_{p\times
q}$ be the $p\times q$ matrix of all zeros. For an $p\times q$ matrix ${\bf
A}=[A_{ij}]$, we denote $({\bf{a}}^{(i)})^{\top}$ be its $i$-th row, ${\rm
nnz}({\bf A})$ be the number of non-zero entries of ${\bf A}$. The condensed
singular value decomposition (SVD) of ${\bf A}$ is defined as ${\bf
U}{\bf{\Sigma}}{\bf V}^{\top}$ where ${\bf U}\in{\mathbb{R}}^{m\times r}$ and
${\bf V}\in{\mathbb{R}}^{m\times r}$ are column orthogonal,
${\bf{\Sigma}}=\operatorname{diag}(\sigma_{1}({\bf A}),\sigma_{2}({\bf
A}),\dots,\sigma_{r}({\bf A}))$ with $\sigma_{1}({\bf A})\geq\sigma_{2}({\bf
A})\geq\dots\geq\sigma_{r}({\bf A})>0$ places the nonzero singular values on
its diagonal entries and $r$ is the rank of ${\bf A}$. We have
$\sigma_{i}({\bf A})=0$ for any $i>r$. Additionally, we let $\left\|{\bf
A}\right\|_{F}=\sqrt{\sum_{i,j}A_{ij}^{2}}=\sqrt{\sum_{i=1}^{r}\sigma_{i}^{2}({\bf
A})}$ be the Frobenius norm, $\left\|{\bf A}\right\|_{2}=\sigma_{1}({\bf A})$
be the spectral norm, $\left\|{\bf
A}\right\|_{*}=\sum_{i=1}^{r}\sigma_{i}({\bf A})$ be the nuclear norm and
$\left\|{\bf A}\right\|_{k}=\sum_{i=1}^{k}\sigma_{i}({\bf A})$ be the Ky Fan
$k$-norm. We also denote ${\bf A}_{k}$ as the best rank-$k$ approximation to
${\bf A}$ for any unitary invariant norms, that is, ${\bf
A}_{k}=\sum_{i=1}^{k}\sigma_{i}({\bf A}){\bf{u}}_{i}{\bf{v}}_{i}^{\top}$,
where ${\bf{u}}_{i}$ and ${\bf{v}}_{i}$ are the $i$-th column of ${\bf U}$ and
${\bf V}$ respectively.

###  2.2 Frequent Directions

Frequent directions [19, 12] is a deterministic algorithm for covariance
sketching. Given any matrix ${\bf X}\in{\mathbb{R}}^{n\times d}$ and sketch
size $m$ that is much smaller than $n$ and $d$, FD processes the rows of ${\bf
X}$ one by one and produces a sketch matrix ${\bf A}\in{\mathbb{R}}^{2m\times
d}$ to approximate ${\bf X}^{\top}{\bf X}$ by ${\bf A}^{\top}{\bf A}$. We
present the details of FD in Algorithm 1, which requires ${\mathcal{O}}(md)$
space and ${\mathcal{O}}(mnd)$ time complexity. The algorithm has the
following theoretical guarantees.

######  Lemma 1 ([10, 12]).

The output ${\bf A}$ of Algorithm 1 satisfies

$\displaystyle\left\|{\bf X}^{\top}{\bf X}-{\bf A}^{\top}{\bf
A}\right\|_{2}\leq\frac{1}{m-k}\left(\left\|{\bf
X}\right\|_{F}^{2}-\left\|{\bf X}_{k}\right\|_{F}^{2}\right)$ (1)

for any $k<m$.

Ghashami et al. [12] also prove FD is space optimal with respect to the
guaranteed accuracy in Lemma 1. Note that the shrinking step in line 8 of the
algorithm is necessary because the output could be extremely worse without
this operation [5, 22].

Algorithm 1 Frequent Directions (FD)

1: Input: ${\bf X}\in{\mathbb{R}}^{n\times d}$ and sketch size $m$

2: ${\bf A}\leftarrow{\bf{0}}_{2m\times d}$

3: for $t=1,2,\dots,n$

4: insert $({\bf{x}}^{(t)})^{\top}$ into a zero valued row of ${\bf A}$

5: if ${\bf A}$ has no zero valued rows then

6: $[{\bf U},{\bf{\Sigma}},{\bf V}]\leftarrow{\rm SVD}({\bf A})$

7: $\delta\leftarrow\sigma^{2}_{m}({\bf A})$

8:
$\widehat{{\bf{\Sigma}}}\leftarrow\sqrt{\max\left({\bf{\Sigma}}^{2}-\delta{\bf
I}_{2m},{\bf{0}}_{2m\times 2m}\right)}$

9: ${\bf A}\leftarrow\widehat{{\bf{\Sigma}}}{\bf V}^{\top}$

10: end if

11: end for

12: Output: ${\bf A}$

###  2.3 Sketching Algorithms for AMM

It is natural to exploit the idea of FD to solve general AMM problem [37]. We
can concatenates the input matrix ${\bf X}\in{\mathbb{R}}^{n\times d_{x}}$ and
${\bf Y}\in{\mathbb{R}}^{n\times d_{y}}$ to construct a larger matrix ${\bf
Z}=[{\bf X},{\bf Y}]\in{\mathbb{R}}^{n\times(d_{x}+d_{y})}$, and then apply FD
on ${\bf Z}$ to approximate ${\bf Z}^{\top}{\bf Z}$ by ${\bf C}^{\top}{\bf
C}$, where ${\bf C}=[{\bf A},{\bf B}]$, ${\bf A}\in{\mathbb{R}}^{n\times
d_{x}}$ and ${\bf B}\in{\mathbb{R}}^{n\times d_{y}}$. The top right sub-matrix
of the ${\bf C}^{\top}{\bf C}$, i.e.,the matrix ${\bf A}^{\top}{\bf B}$ is an
approximation of ${\bf X}^{\top}{\bf Y}$. Intuitively, this algorithm wastes a
large proportion of cost to approximate ${\bf X}^{\top}{\bf X}$ and ${\bf
Y}^{\top}{\bf Y}$ (the other sub-matrices of ${\bf Z}^{\top}{\bf Z}$), which
is unnecessary for the AMM task.

Mroueh et al. [25] proposed the co-occurring directions (COD) for AMM. We
present its detailed procedure in Algorithm 2. Each iteration of COD
constructs the column basis of ${\bf A}$ and ${\bf B}$ by QR factorization
independently and executes the shrinkage step on the small interaction matrix
${\bf R}_{x}{\bf R}_{y}^{\top}$. We point out that both COD and FD-AMM
requires ${\mathcal{O}}(m(d_{x}+d_{y}))$ space and
${\mathcal{O}}(mn(d_{x}+d_{y}))$ time complexity. However, COD looks more
reasonable than FD-AMM since all of its operations surrounds approximating
${\bf X}^{\top}{\bf Y}$. The numerical experiments [25] show that COD performs
significantly better than FD-AMM [37] and other AMM algorithms [7, 4, 31] when
input matrices is dense. We can prove that COD holds the guaranteed accuracy
as follows.

######  Lemma 2 ([25]).

The output ${\bf A}$ and ${\bf B}$ of Algorithm 2 satisfies

$\displaystyle\left\|{\bf X}^{\top}{\bf Y}-{\bf A}^{\top}{\bf
B}\right\|_{2}\leq\frac{\left\|{\bf X}\right\|_{F}\left\|{\bf
Y}\right\|_{F}}{m}.$ (2)

Unfortunately, the result of Lemma 2 does not reveal the advantage of COD
entirely. Consider that case of ${\bf X}={\bf Y}$, the procedure of COD will
reduce to FD, but the error bound of (2) becomes a special case of (1) in
Lemma 1 with $k=0$. The real-world dataset typically enjoys some approximate
low-rank structure, which leads to the right-hand side of bound (1) could be
much smaller than the one of (2). Hence although COD has better empirical
performance, the existing error bounds are not tight enough.

Algorithm 2 Co-Occurring Directions (COD)

1: Input: ${\bf X}\in{\mathbb{R}}^{n\times d_{x}}$, ${\bf
Y}\in{\mathbb{R}}^{n\times d_{y}}$ and sketch size $m$

2: ${\bf A}\leftarrow{\bf{0}}_{2m\times d_{x}}$

3: ${\bf B}\leftarrow{\bf{0}}_{2m\times d_{y}}$

4: for $t=1,2,\dots,n$

5: insert $({\bf{x}}^{(t)})^{\top}$ into a zero valued row of ${\bf A}$

6: insert $({\bf{y}}^{(t)})^{\top}$ into a zero valued row of ${\bf B}$

7: if ${\bf A}$ or ${\bf B}$ has no zero valued rows then

8: $({\bf Q}_{x},{\bf R}_{x})\leftarrow{\rm QR}\left({\bf A}^{\top}\right)$

9: $({\bf Q}_{y},{\bf R}_{y})\leftarrow{\rm QR}\left({\bf B}^{\top}\right)$

10: $[{\bf U},{\bf{\Sigma}},{\bf V}]\leftarrow{\rm SVD}({\bf R}_{x}{\bf
R}_{y}^{\top})$

11: $\delta\leftarrow\sigma_{m}({\bf R}_{x}{\bf R}_{y}^{\top})$

12: $\widehat{{\bf{\Sigma}}}\leftarrow\max\left({\bf{\Sigma}}-\delta{\bf
I}_{2m},{\bf{0}}_{2m\times 2m}\right)$

13: ${\bf A}\leftarrow\widehat{{\bf{\Sigma}}}^{1/2}{\bf U}^{\top}{\bf
Q}_{x}^{\top}$

14: ${\bf B}\leftarrow\widehat{{\bf{\Sigma}}}^{1/2}{\bf V}^{\top}{\bf
Q}_{y}^{\top}$

15: end if

16: end for

17: Output: ${\bf A}$ and ${\bf B}$

##  3 Sharper Analysis for COD

In this section, we provide a tighter error bound for COD. We let
$\delta^{(t)}$ be the value of $\delta$ at time step $t$. If the algorithm
does not enter the “if” section in the $t$-th step, then we have
$\delta^{(t)}=0$. Similarly, let ${\bf A}^{(t)}$, ${\bf B}^{(t)}$, ${\bf
Q}_{x}^{(t)}$, ${\bf Q}_{y}^{(t)}$, ${\bf U}^{(t)}$, ${\bf{\Sigma}}^{(t)}$,
${\bf V}^{(t)}$ and $\widehat{{\bf{\Sigma}}}^{(t)}$ be the corresponding
variables after the main loop has been executed for $t$ times. Additionally,
we use ${\widehat{\bf A}}^{(t)}$ and ${\widehat{\bf B}}^{(t)}$ to represent
the matrices after insert operations (line 5-6) have been executed at the
$t$-th iteration. We need the following two lemmas for proving our main
results.

######  Lemma 3 ([25]).

The output matrices ${\bf A}$ and ${\bf B}$ of Algorithm 2 satisfy

$\displaystyle\left\|{\bf X}^{\top}{\bf Y}-{\bf A}^{\top}{\bf
B}\right\|_{2}\leq\sum_{t=1}^{n}\delta^{(t)}$ (3)

and

$\displaystyle\left\|{\bf A}^{\top}{\bf B}\right\|_{*}\leq\left\|{\bf
X}\right\|_{F}\left\|{\bf Y}\right\|_{F}-m\sum_{t=1}^{n}\delta^{(t)}.$ (4)

######  Lemma 4.

The output of Algorithm 2 holds that

$\displaystyle\left\|{\bf X}^{\top}{\bf Y}\right\|_{*}-\left\|{\bf
A}^{\top}{\bf B}\right\|_{*}\leq\sum_{i=k+1}^{d}\sigma_{i}({\bf X}^{\top}{\bf
Y})+k\sum_{t=1}^{n}\delta^{(t)}.$ (5)

Lemma 4 is the key lemma of our proof. It improves the result in analysis of
COD [25]. The term $\sum_{i=k+1}^{d}\sigma_{i}({\bf X}^{\top}{\bf Y})$ on the
right-hand side of (5) considers the potential approximate low-rank structure
of ${\bf X}^{\top}{\bf Y}$, which leads to a tighter error bound of COD as
follows.

######  Theorem 1.

The output of Algorithm 2 holds that

$\displaystyle\left\|{\bf X}^{\top}{\bf Y}-{\bf A}^{\top}{\bf
B}\right\|_{2}\leq\frac{1}{m-k}\Big{(}\left\|{\bf X}\right\|_{F}\left\|{\bf
Y}\right\|_{F}-\big{\|}{\bf X}^{\top}{\bf Y}\big{\|}_{k}\Big{)}.$

###### Proof.

Let $\Delta=\sum_{t=1}^{n}\delta^{(t)}$. Connecting inequality (4) in Lemma 3
and the result of Lemma 4, we have

$\displaystyle m\Delta+\left\|{\bf X}^{\top}{\bf Y}\right\|_{*}-\left\|{\bf
X}\right\|_{F}\left\|{\bf Y}\right\|_{F}$ $\displaystyle\leq$
$\displaystyle\left\|{\bf X}^{\top}{\bf Y}\right\|_{*}-\left\|{\bf
A}^{\top}{\bf B}\right\|_{*}$ $\displaystyle\leq$
$\displaystyle\sum_{i=k+1}^{d}\sigma_{i}({\bf X}^{\top}{\bf Y})+k\Delta,$

that is

$\displaystyle\Delta\leq\dfrac{1}{m-k}\left(\left\|{\bf
X}\right\|_{F}\left\|{\bf Y}\right\|_{F}-\left\|{\bf X}^{\top}{\bf
Y}\right\|_{k}\right).$

Substituting above bound of $\Delta$ into inequality (4) of Lemma 3, we finish
the proof of this theorem. ∎

To achieve the accuracy that $\left\|{\bf X}^{\top}{\bf Y}-{\bf A}^{\top}{\bf
B}\right\|_{2}\leq{\varepsilon}$, the previous error bound (Lemma 2) requires
the sketch size to be at least $m_{1}=\frac{1}{{\varepsilon}}\left\|{\bf
X}\right\|_{F}\left\|{\bf Y}\right\|_{F}$ while Theorem 1 only requires the
sketch size $m_{2}=k+\frac{1}{{\varepsilon}}\left(\left\|{\bf
X}\right\|_{F}\left\|{\bf Y}\right\|_{F}-\left\|{\bf X}^{\top}{\bf
Y}\right\|_{k}\right)$. When input matrices ${\bf X}$ and ${\bf Y}$ have
strongly correlation and approximate low-rank structure, $m_{2}$ could be much
smaller than $m_{1}$.

In addition, the error bound of Theorem 1 matches that of FD (Lemma 1) when
${\bf X}={\bf Y}$:

$\displaystyle\left\|{\bf X}^{\top}{\bf X}-{\bf A}^{\top}{\bf A}\right\|_{2}$
$\displaystyle\leq$ $\displaystyle\frac{1}{m-k}\Big{(}\left\|{\bf
X}\right\|_{F}\left\|{\bf X}\right\|_{F}-\sum_{i=1}^{k}\sigma_{i}({\bf
X}^{\top}{\bf X})\Big{)}$ $\displaystyle=$
$\displaystyle\frac{1}{m-k}\Big{(}\left\|{\bf
X}\right\|_{F}^{2}-\sum_{i=1}^{k}\sigma_{i}^{2}({\bf X})\Big{)}$
$\displaystyle=$ $\displaystyle\frac{1}{m-k}\left(\left\|{\bf
X}\right\|_{F}^{2}-\left\|{\bf X}_{k}\right\|_{F}^{2}\right).$

On the other hand, the previous error bound (Lemma 2) is worse than that of FD
(Lemma 1) in the symmetric case of ${\bf X}={\bf Y}$.

##  4 Space Lower Bounds Analysis

In this section, we show that COD is space optimal with respect to our new
error bound in Theorem 1. We first introduce the following lemma for low-rank
matrices.

######  Lemma 5 ([16]).

For each $\delta>0$ there exits a set of matrices ${\mathcal{Q}}=\\{{\bf
Q}_{1},\cdots,{\bf Q}_{N}\\}$ and $N=2^{\Omega(\ell(d-\ell)\log(1/\delta))}$,
where ${\bf Q}_{i}\in{\mathbb{R}}^{\ell\times d}$ with ${\bf Q}_{i}{\bf
Q}_{i}^{\top}={\bf I}_{\ell}$, such that $\left\|{\bf Q}_{i}{\bf
Q}_{j}^{\top}\right\|_{2}<1-\delta$.

By using Lemma 5, we can construct a sets contains exponential number of
matrices that each pair of them are not “too close”. The formalized result is
shown in Lemma 6.

######  Lemma 6.

For each $\delta>0$ and $d_{x}\leq d_{y}$ there exits a set of matrices

$\displaystyle\widehat{\mathcal{Z}}_{\ell}=\left\\{({\widehat{\bf
X}}^{(1)},{\widehat{\bf Y}}^{(1)}),\cdots,({\widehat{\bf
X}}^{(N)},{\widehat{\bf Y}}^{(N)})\right\\},$

where $N=2^{\Omega(\ell(d_{y}-\ell)\log(1/\delta))}$ and ${\widehat{\bf
X}}^{(i)}\in{\mathbb{R}}^{\ell\times d_{x}},{\widehat{\bf
Y}}^{(i)}\in{\mathbb{R}}^{\ell\times d_{y}}$ satisfy ${\widehat{\bf
X}}^{(i)}{\widehat{\bf X}}^{{(i)}\top}={\bf
I}_{\ell}~{}\text{and}~{}{\widehat{\bf Y}}^{(i)}{\widehat{\bf
Y}}^{{(i)}\top}={\bf I}_{\ell}$ for any $i=1,\dots,n$ and

$\displaystyle\left\|{\widehat{\bf X}}^{{(i)}\top}{\widehat{\bf
Y}}^{{(i)}\top}-{\widehat{\bf X}}^{{(j)}\top}{\widehat{\bf
Y}}^{{(j)}\top}\right\|_{2}>\sqrt{2\delta}$

for any $j\neq i$.

###### Proof.

Based on Lemma 5, there exist a set of matrices
${\mathcal{Y}}=\\{{\widehat{\bf Y}}^{(1)},\cdots,{\widehat{\bf Y}}^{(N)}\\}$,
$N=2^{\Omega(\ell(d-\ell)\log(1/\delta))}$, where ${\widehat{\bf
Y}}^{(i)}\in{\mathbb{R}}^{\ell\times d}$ with ${\widehat{\bf
Y}}^{(i)}{\widehat{\bf Y}}^{{(i)}\top}={\bf I}_{\ell}$, such that
$\big{\|}{\widehat{\bf Y}}^{(i)}{\widehat{\bf
Y}}^{{(j)}\top}\big{\|}<1-\delta$. We further set ${\widehat{\bf
X}}^{(i)}=[{\bf I}_{\ell},{\bf 0}_{\ell\times(d_{x}-\ell)}]$. We have

$\displaystyle\left\|{\widehat{\bf X}}^{{(i)}\top}{\widehat{\bf
Y}}^{(i)}-{\widehat{\bf X}}^{{(j)}\top}{\widehat{\bf
Y}}^{(j)}\right\|_{2}^{2}$ $\displaystyle=$ $\displaystyle\left\|{\widehat{\bf
Y}}^{(i)}-{\widehat{\bf Y}}^{(j)}\right\|_{2}^{2}$ $\displaystyle=$
$\displaystyle\left\|({\widehat{\bf Y}}^{(i)}-{\widehat{\bf
Y}}^{(j)})({\widehat{\bf Y}}^{{(i)}\top}-{\widehat{\bf
Y}}^{{(j)}\top})\right\|_{2}$ $\displaystyle\geq$ $\displaystyle
2-\left\|{\widehat{\bf Y}}^{(j)}{\widehat{\bf Y}}^{{(i)}\top}+{\widehat{\bf
Y}}^{(i)}{\widehat{\bf Y}}^{{(j)}\top}\right\|_{2}$ $\displaystyle\geq$
$\displaystyle 2\delta,$

where we use the definition of ${\widehat{\bf X}}^{(i)}$, ${\widehat{\bf
Y}}^{(i)}$ and the fact $\left\|{\bf A}^{\top}{\bf A}\right\|_{2}=\left\|{\bf
A}\right\|_{2}^{2}$. ∎

Algorithm 3 Subspace Power Method (SPM)

1: Input: ${\bf M}\in{\mathbb{R}}^{d_{1}\times d_{2}}$, target rank $m$ and
integer $q>0$

2: ${\bf G}=[G_{ij}]\in{\mathbb{R}}^{d_{2}\times m}$, where
$G_{ij}\sim{\mathcal{N}}(0,1)$ i.i.d

3: ${\bf K}=\left({\bf M}{\bf M}^{\top}\right)^{q}{\bf M}{\bf
G}\in{\mathbb{R}}^{d_{1}\times m}$

4: ${\bf Z}\leftarrow$ orthonormal column basis of ${\bf K}$

5: Output: ${\bf Z}$

Then we present a lower bound of space complexity for approximate matrix
multiplication, which matches the memory cost of COD. Hence, we can conclude
that COD is space optimal with respect to the guaranteed accuracy in Theorem
1.

######  Theorem 2.

We consider any matrix sketching algorithm with inputs as ${\bf
X}\in{\mathbb{R}}^{n\times d_{x}}$ and ${\bf Y}\in{\mathbb{R}}^{n\times
d_{y}}$ and outputs ${\bf A}\in{\mathbb{R}}^{m\times d_{x}}$ and ${\bf
B}\in{\mathbb{R}}^{m\times d_{y}}$ with guarantee

$\displaystyle\left\|{\bf X}^{\top}{\bf Y}-{\bf A}^{\top}{\bf
B}\right\|_{2}\leq\frac{1}{m-k}\Big{(}\left\|{\bf X}\right\|_{F}\left\|{\bf
Y}\right\|_{F}-\big{\|}{\bf X}^{\top}{\bf Y}\big{\|}_{k}\Big{)}$

for any $k<\ell$. Assuming that a constant number of bits is required to
describe a word (i.e., a unit of memory), then the algorithm requires at least
$\Omega(m(d_{x}+d_{y}))$ bits of space.

###### Proof.

Without loss of generality, we suppose $d_{y}\geq d_{x}$. Let
$\widehat{\mathcal{Z}}_{\ell}=\\{({\widehat{\bf X}}^{(1)},{\widehat{\bf
Y}}^{(1)}),\cdots,({\widehat{\bf X}}^{(N)},{\widehat{\bf Y}}^{(N)})\\}$ be the
set of matrices defined in Lemma 6 with $\ell=m/4$, $\delta=1/8$ and
$N=2^{\Omega(\frac{m}{4}\cdot(d_{y}-m/4)\log(8))}$.

We construct the matrices ${\bf X}^{(i)}=[{\widehat{\bf X}}^{(i)};{\bf
0}_{(n-m/4)\times d_{x}}]\in{\mathbb{R}}^{n\times d_{x}}$ and ${\bf
Y}^{(i)}=[{\widehat{\bf Y}}^{(i)};{\bf 0}_{(n-m/4)\times
d_{y}}]\in{\mathbb{R}}^{n\times d_{y}}$ for any $i=1,\dots,N$. Then we have
${\mathcal{Z}}_{\ell}=\\{({\bf X}^{(i)},{\bf Y}^{(i)})\\}_{i=1}^{N}$ which
satisfies

$\displaystyle\left\|{\bf X}^{{(i)}\top}{\bf Y}^{(i)}-{\bf X}^{{(j)}\top}{\bf
Y}^{(j)}\right\|_{2}>1/2$

for each $i\neq j$. Let $({\bf A},{\bf B})$ be the output of the matrix
sketching algorithm with input $({\bf X}^{(i)},{\bf Y}^{(i)})$. The guarantee
of the algorithm indicates

$\displaystyle\left\|{\bf X}^{{(i)}\top}{\bf Y}^{(i)}-{\bf A}^{\top}{\bf
B}\right\|_{2}\leq\frac{1}{m-k}\left(\frac{m}{4}-k\right)\leq\frac{1}{4}.$

Hence, each $({\bf A},{\bf B})$ only encodes one matrix pencil in
${\mathcal{Z}}_{\ell}$ (the product of the matrices), which means that the
lower bound of space to attach the desired accuracy is
$\log_{2}N=\Omega(md_{y})=\Omega(m(d_{x}+d_{y}))$ bits. ∎

##  5 Sparse Co-Occurring Directions

In this section, we proposed a variant of COD for sparse AMM. We also prove
its error bound is similar to our improved result of COD.

###  5.1 The Algorithm

We describe details of our sparse co-occurring directions (SCOD) in Algorithm
4. The procedure of SCOD maintains the sparse data in two buffer matrices
${\bf X}^{\prime}$ and ${\bf Y}^{\prime}$. The algorithm restricts the non-
zero entries in buffers to be less than $m(d_{x}+d_{y})$ and the number of row
of each buffer is at most $d_{x}+d_{y}$. When the buffers are full, we perform
subspace power method (SPM) [35, 26] to approximate the data in the buffers by
low-rank matrices ${\widetilde{\bf X}}\in{\mathbb{R}}^{m\times d_{x}}$ and
${\widetilde{\bf Y}}\in{\mathbb{R}}^{m\times d_{y}}$ such that ${\bf
X}^{\prime\top}{\bf Y}^{\prime}\approx{\widetilde{\bf
X}}^{\top}{\widetilde{\bf Y}}$. We present the procedure of SPM in Algorithm
3.

Let ${\widetilde{\bf X}}^{(i)}$ and ${\widetilde{\bf Y}}^{(i)}$ be the results
of ${\widetilde{\bf X}}$ and ${\widetilde{\bf Y}}$ after Algorithm 4 has
executed “then” section for $i$-times. Define ${\bf C}=[{\widetilde{\bf
X}}^{(1)};\cdots;{\widetilde{\bf X}}^{(T)}]$ and ${\bf D}=[{\widetilde{\bf
Y}}^{(1)};\cdots;{\widetilde{\bf Y}}^{(T)}]$ where $T$ is the number of total
times we enter “then” section of the algorithm. Then ${\bf C}$ and ${\bf D}$
are the estimators of ${\bf X}$ and ${\bf Y}$ respectively and the procedure
of SCOD can be regarded as as running standard COD on input matrices ${\bf C}$
and ${\bf D}$ in streaming fashion. Since the row numbers of buffers ${\bf
X}^{\prime}$ and ${\bf Y}^{\prime}$ could be much larger than $m$, the
operations on dense matrices (line 14-20) will not be executed frequently.
Hence, SCOD is much more efficient than COD for sparse inputs.

###  5.2 Analysis of Error Bound

The analysis of SCOD is more challenging than sparse frequent directions (SFD)
[11] which only addresses the case of ${\bf X}={\bf Y}$. The reason is the
“mergeability property” of FD [12, 5, 11] only works for Frobenius norm and it
is not applicable to COD.

The approximation error of SCOD comes from two parts: the compressing error
from sub-routine SPM and the merge error from estimating ${\bf C}^{\top}{\bf
D}$ by ${\bf A}^{\top}{\bf B}$. We first consider a single call of SPM, which
approximation error can be bounded as follows.

######  Lemma 7.

Let $q=\tilde{\Theta}(\log(md_{1}/p)/{\varepsilon})$ for Algorithm 3, then the
output ${\bf Z}$ satisfies

$\displaystyle\left\|{\bf M}-{\bf Z}{\bf Z}^{\top}{\bf
M}\right\|_{2}\leq(1+{\varepsilon})\sigma_{m+1}({\bf M})$

with probability at least $1-p$.

Based on Lemma 7, we can bound the total compressing error of SOCD by the
following lemma.

######  Lemma 8.

Setting $q_{i}=\tilde{\Theta}(\log(md_{1}/p_{i})/{\varepsilon})$ and
$p_{i}=\delta/2i^{2}$, then we have then Algorithm 4 holds that

$\displaystyle\left\|{\bf X}^{\top}{\bf Y}-{\bf C}^{\top}{\bf
D}\right\|_{2}\leq\dfrac{1+{\varepsilon}}{m-k}\left(\left\|{\bf
X}\right\|_{F}\left\|{\bf Y}\right\|_{F}-\big{\|}{\bf X}^{\top}{\bf
Y}\big{\|}_{k}\right),$

for any $k<m$ and ${\varepsilon}>0$ with probability $1-\delta$.

###### Proof.

Let ${\bf X}^{\prime{(i)}}$ and ${\bf Y}^{\prime{(i)}}$ be the value of ${\bf
X}^{\prime}$ and ${\bf Y}^{\prime}$ when we execute subspace power methods at
$i$-th time in line 8 of Algorithm 4, then we have ${\bf X}=[{\bf
X}^{\prime(1)};\dots;{\bf X}^{\prime(T)}]$ and ${\bf Y}=[{\bf
Y}^{\prime(1)};\dots;{\bf Y}^{\prime(T)}]$. Using Lemma 7 with ${\bf
M}={\widetilde{\bf X}}^{{(i)}\top}{\widetilde{\bf Y}}^{(i)}$ and $q=q_{i}$,
then with probability $1-p_{i}$, we have

$\displaystyle\begin{split}&\left\|{\widetilde{\bf
X}}^{{(i)}\top}{\widetilde{\bf Y}}^{(i)}-{\bf X}^{\prime{(i)}\top}{\bf
Y}^{\prime{(i)}}\right\|_{2}\\\ \leq&(1+{\varepsilon})\sigma_{m+1}\left({\bf
X}^{\prime{(i)}\top}{\bf Y}^{\prime{(i)}}\right)\\\
\leq&\frac{1+{\varepsilon}}{m-k}\left(\big{\|}{\bf
X}^{\prime{(i)}}\big{\|}_{F}\big{\|}{\bf
Y}^{\prime{(i)}}\big{\|}_{F}-\big{\|}{\bf X}^{\prime{(i)}\top}{\bf
Y}^{\prime{(i)}}\big{\|}_{k}\right),\end{split}$ (6)

where the last step use Srebro et al.’s (2005) Lemma 1 such that $\|{\bf
X}^{\prime{(i)}\top}{\bf Y}^{\prime{(i)}}\|_{*}\leq\|{\bf
X}^{\prime{(i)}}\|_{F}\|{\bf Y}^{\prime{(i)}}\|_{F}$.

Summing over inequality (6) with $i=1,\dots,T$, we have

$\displaystyle\left\|{\bf X}^{\top}{\bf Y}-{\bf C}^{\top}{\bf
D}\right\|_{2}\leq\sum_{t=1}^{T}\left\|{\bf X}^{\prime{(i)}}{\bf
Y}^{\prime{(i)}}-{\widetilde{\bf X}}^{{(i)}\top}{\widetilde{\bf
Y}}^{(i)}\right\|_{2}$ $\displaystyle\leq$
$\displaystyle\frac{1+{\varepsilon}}{m-k}\sum_{t=1}^{T}\left(\big{\|}{\bf
X}^{\prime{(i)}}\big{\|}_{F}\big{\|}{\bf
Y}^{\prime{(i)}}\big{\|}_{F}-\big{\|}{\bf X}^{\prime{(i)}\top}{\bf
Y}^{\prime{(i)}}\big{\|}_{k}\right)\leq\frac{1+{\varepsilon}}{m-k}\left(\big{\|}{\bf
X}\big{\|}_{F}\big{\|}{\bf Y}\big{\|}_{F}-\big{\|}{\bf X}^{\top}{\bf
Y}\big{\|}_{k}\right)$

with probability $1-\delta$. The last inequality is based on the
Cauchy–Schwarz inequality and the triangle inequality of Ky Fan $k$-norm. Note
that the failure probability is no more than
$p_{1}+\dots+p_{T}=\frac{\delta}{2}\sum_{i=1}^{T}1/i^{2}\leq\delta$. ∎

Unlike SFD [12] which introduces a verifying step to boost the success
probability, our method instead requires $q_{i}$ to be increased
logarithmically to ensure the error bound of SCOD holds with probability at
least $1-\delta$ for given $\delta\in(0,1)$.

Another important property of SCOD is that the compression step shrink the
magnitude of the product of input matrices. The steps in line 10-11 of
Algorithm 4 balance the singular values of ${\widetilde{\bf X}}$ and
${\widetilde{\bf Y}}$, which leads to the following lemma:

######  Lemma 9.

Algorithm 4 holds that

$\displaystyle\big{\|}{\widetilde{\bf
X}}^{(i)}\big{\|}_{F}\big{\|}{\widetilde{\bf
Y}}^{(i)}\big{\|}_{F}\leq\big{\|}{\bf X}^{\prime{(i)}}\big{\|}_{F}\big{\|}{\bf
Y}^{\prime{(i)}}\big{\|}_{F}.$

Since the analysis of merging error is similar to standard COD, we can
establish the error bound of SCOD by using above lemmas.

######  Theorem 3.

Setting $q_{i}=\tilde{\Theta}(\log(md_{1}/p_{i})/{\varepsilon})$ with constant
${\varepsilon}>0$ and $p_{i}=\delta/2i^{2}$, with probability $1-\delta$, the
outputs ${\bf A}$ and ${\bf B}$ of Algorithm 4 hold that

$\displaystyle\left\|{\bf X}^{\top}{\bf Y}-{\bf A}^{\top}{\bf
B}\right\|_{2}\leq\left(\frac{2+{\varepsilon}}{m-k}+\frac{(1+{\varepsilon})k}{(m-k)^{2}}\right)\left(\left\|{\bf
X}\right\|_{F}\left\|{\bf Y}\right\|_{F}-\big{\|}{\bf X}^{\top}{\bf
Y}\big{\|}_{k}\right)$

for all $k<m$.

###### Proof.

Consider that ${\bf A}$ and ${\bf B}$ can be viewed as the output of running
Algorithm 2 with input matrices

$\displaystyle{\bf C}=[{\widetilde{\bf X}}^{(1)};\cdots;{\widetilde{\bf
X}}^{(T)}]\text{~{}~{}and~{}~{}}{\bf D}=[{\widetilde{\bf
Y}}^{(1)};\cdots;{\widetilde{\bf Y}}^{(T)}].$

Following the proof of Theorem 1, we have

$\displaystyle\begin{split}&\left\|{\bf C}^{\top}{\bf D}-{\bf A}^{\top}\bf{\bf
B}\right\|_{2}\\\
\leq&\frac{1}{m-k}\left(\sum_{i=1}^{T}\big{\|}{\widetilde{\bf
X}}^{(i)}\big{\|}_{F}\big{\|}{\widetilde{\bf
Y}}^{(i)}\big{\|}_{F}-\big{\|}{\bf C}^{\top}{\bf D}\big{\|}_{k}\right)\\\
\leq&\frac{\sum_{i=1}^{T}\big{\|}{\widetilde{\bf
X}}^{(i)}\big{\|}_{F}\big{\|}{\widetilde{\bf
Y}}^{(i)}\big{\|}_{F}-\big{\|}{\bf X}^{\top}{\bf Y}\big{\|}_{k}+\big{\|}{\bf
X}^{\top}{\bf Y}-{\bf C}^{\top}{\bf D}\big{\|}_{k}}{m-k}\\\
\leq&\frac{\sum_{i=1}^{T}\big{\|}{\widetilde{\bf
X}}^{(i)}\big{\|}_{F}\big{\|}{\widetilde{\bf
Y}}^{(i)}\big{\|}_{F}-\big{\|}{\bf X}^{\top}{\bf Y}\big{\|}_{k}+k\big{\|}{\bf
X}^{\top}{\bf Y}-{\bf C}^{\top}{\bf D}\big{\|}_{2}}{m-k}\\\
\leq&\frac{1}{m-k}\left(\sum_{i=1}^{T}\big{\|}{\bf
X}^{\prime{(i)}}\big{\|}_{F}\big{\|}{\bf
Y}^{\prime{(i)}}\big{\|}_{F}-\big{\|}{\bf X}^{\top}{\bf
Y}\big{\|}_{k}\right)\\\
&\quad+\frac{(1+{\varepsilon})k}{(m-k)^{2}}\left(\left\|{\bf
X}\right\|_{F}\left\|{\bf Y}\right\|_{F}-\big{\|}{\bf X}^{\top}{\bf
Y}\big{\|}_{k}\right)\\\
\leq&\left(\frac{1}{m-k}+\frac{(1+{\varepsilon})k}{(m-k)^{2}}\right)\left(\left\|{\bf
X}\right\|_{F}\left\|{\bf Y}\right\|_{F}-\big{\|}{\bf X}^{\top}{\bf
Y}\big{\|}_{k}\right)\end{split}$

where we use Lemma 8, 9 and triangle inequality.

Combing above results and Lemma 8, we have

$\displaystyle\left\|{\bf X}^{\top}{\bf Y}-{\bf A}^{\top}{\bf B}\right\|_{2}$
$\displaystyle\leq$ $\displaystyle\left\|{\bf X}^{\top}{\bf Y}-{\bf
C}^{\top}{\bf D}\right\|_{2}+\left\|{\bf C}^{\top}{\bf D}-{\bf A}^{\top}{\bf
B}\right\|_{2}$ $\displaystyle\leq$
$\displaystyle\left(\frac{2+{\varepsilon}}{m-k}+\frac{(1+{\varepsilon})k}{(m-k)^{2}}\right)\left(\left\|{\bf
X}\right\|_{F}\left\|{\bf Y}\right\|_{F}-\left\|{\bf X}^{\top}{\bf
Y}\right\|_{k}\right),$

with probability at least $1-\delta$. ∎

Algorithm 4 Sparse Co-Occurring Directions (SCOD)

1: Input: ${\bf X}\in{\mathbb{R}}^{n\times d_{x}}$, ${\bf
Y}\in{\mathbb{R}}^{n\times d_{y}}$, sketch size $m$, failure probability
$\delta$ and sequence $\\{q_{i}\\}_{i=1,2\dots}$

2: $i=0$

3: ${\bf A}\leftarrow{\bf{0}}_{m\times d_{x}}$, ${\bf
B}\leftarrow{\bf{0}}_{m\times d_{y}}$

4: ${\bf X}^{\prime}\leftarrow{\rm empty}$, ${\bf Y}^{\prime}\leftarrow{\rm
empty}$

5: for $t=1,2,\dots,n$

6: ${\bf X}^{\prime}\leftarrow[{\bf X}^{\prime};({\bf{x}}^{(t)})^{\top}]$,
${\bf Y}^{\prime}\leftarrow[{\bf Y}^{\prime};({\bf{y}}^{(t)})^{\top}]$

7: if  ${\rm nnz}({\bf X}^{\prime})+{\rm nnz}({\bf
Y}^{\prime})>m(d_{x}+d_{y})$ or $t=n$ or ${\rm rows}({\bf
X}^{\prime})=d_{x}+d_{y}$ or ${\rm rows}({\bf Y}^{\prime})=d_{x}+d_{y}$ then

8: ${\bf Z}={\rm{SubspacePowerMethod}}({\bf X}^{\prime\top}{\bf
Y}^{\prime},m,q_{i})$

9: $[\widetilde{{\bf U}},\widetilde{{\bf{\Sigma}}},\widetilde{{\bf V}}]={\rm
SVD}\left({\bf Z}^{\top}{\bf X}^{\prime\top}{\bf Y}^{\prime}\right)$

10: ${\widetilde{\bf
X}}\leftarrow\widetilde{{\bf{\Sigma}}}^{1/2}\widetilde{{\bf U}}^{\top}{\bf
Z}^{\top}$

11: ${\widetilde{\bf
Y}}\leftarrow\widetilde{{\bf{\Sigma}}}^{1/2}\widetilde{{\bf V}}^{\top}$

12: ${\bf A}\leftarrow[{\bf A};{\widetilde{\bf X}}]$

13: ${\bf B}\leftarrow[{\bf B};{\widetilde{\bf Y}}]$

14: $({\bf Q}_{x},{\bf R}_{x})\leftarrow{\rm QR}\left({\bf A}^{\top}\right)$

15: $({\bf Q}_{y},{\bf R}_{y})\leftarrow{\rm QR}\left({\bf B}^{\top}\right)$

16: $[{\bf U},{\bf{\Sigma}},{\bf V}]\leftarrow{\rm SVD}({\bf R}_{x}{\bf
R}_{y}^{\top})$

17: $\delta\leftarrow\sigma_{m}({\bf R}_{x}{\bf R}_{y}^{\top})$

18: $\widehat{{\bf{\Sigma}}}\leftarrow\max\left({\bf{\Sigma}}-\delta{\bf
I}_{m},{\bf{0}}_{m\times m}\right)$

19: ${\bf A}\leftarrow\widehat{{\bf{\Sigma}}}^{1/2}{\bf U}^{\top}{\bf
Q}_{x}^{\top}$

20: ${\bf B}\leftarrow\widehat{{\bf{\Sigma}}}^{1/2}{\bf V}^{\top}{\bf
Q}_{y}^{\top}$

21: ${\bf X}^{\prime}\leftarrow{\rm empty}$, ${\bf Y}^{\prime}\leftarrow{\rm
empty}$

22: $i\leftarrow i+1$

23: end if

24: end for

25: Output: ${\bf A}$ and ${\bf B}$

###  5.3 Complexity Analysis

We use the constant-word-size model for our analysis like that of sparse FD
[12]. We suppose floating point numbers are represented by a constant number
of bits, random access into memory requires ${\mathcal{O}}(1)$ time and
multiplying a sparse matrix ${\bf M}$ by a dense vector requires
${\mathcal{O}}({\rm nnz}({\bf M}))$ time and storing ${\bf M}$ requires
${\mathcal{O}}({\rm nnz}({\bf M}))$ space.

The procedure of SCOD (Algorithm 4) implies the buffer ${\bf X}^{\prime}$ and
${\bf Y}^{\prime}$ is sparse and contains at most $m(d_{x}+d_{y})$ non-zero
entries and it is not difficult to verify that all dense matrices in the
algorithm cost no more than ${\mathcal{O}}(m(d_{x}+d_{y}))$ space. Hence, the
space complexity of SCOD is ${\mathcal{O}}(m(d_{x}+d_{y}))$ in total which is
the same as COD (Algorithm 2).

Then we analyze the time complexity of SCOD. The constraints on buffer size
means we have

$\displaystyle T\leq\frac{{\rm nnz}({\bf X})+{\rm nnz}({\bf
Y})}{m(d_{x}+d_{y})}+\frac{n}{d_{x}+d_{y}}.$

Since each QR factorization or SVD on $m\times d$ matrix cost
${\mathcal{O}}(m^{2}d)$ time, the operation on dense matrices of Algorithm 4
from line 9-20 requires at most

$\displaystyle{\mathcal{O}}(m^{2}(d_{x}+d_{y})T)={\mathcal{O}}(m({\rm
nnz}({\bf X})+{\rm nnz}({\bf Y}))+m^{2}n).$

Note that SCOD calls SPM with input ${\bf M}={\bf X}^{\prime\top}{\bf
Y}^{\prime}$. Since both ${\bf X}^{\prime}$ and ${\bf Y}^{\prime}$ are sparse,
it is unnecessary to construct ${\bf M}$ explicitly and we can multiply ${\bf
X}^{\prime}$ and ${\bf Y}^{\prime}$ on ${\bf G}$ separately in line 3 of
Algorithm 3. Then the time complexity of executing SPM needs
${\mathcal{O}}(mq_{i}({\rm nnz}({\bf X}^{\prime{(i)}})+{\rm nnz}({\bf
Y}^{\prime{(i)}}))+m^{2}d_{x})$ when the algorithm enters “then” section at
the $i$-th time. Following the upper bound of $T$ and the setting of $q_{i}$
in Theorem 3, the calls of SPM in Algorithm 3 entirely takes at most

$\displaystyle{\mathcal{O}}\left(\sum_{i=1}^{T}\left(mq_{i}({\rm nnz}({\bf
X}^{\prime{(i)}})+{\rm nnz}({\bf Y}^{\prime{(i)}}))+m^{2}d_{x}\right)\right)$
$\displaystyle\leq$ $\displaystyle{\mathcal{O}}\left(mq_{T}({\rm nnz}({\bf
X})+{\rm nnz}({\bf Y}))+Tm^{2}d_{x}\right)$ $\displaystyle=$
$\displaystyle\tilde{\mathcal{O}}\big{(}m({\rm nnz}({\bf X})+{\rm nnz}({\bf
Y}))+m^{2}n\big{)}.$

Hence, the total time complexity of proposed SCOD is
$\tilde{\mathcal{O}}\big{(}m({\rm nnz}({\bf X})+{\rm nnz}({\bf
Y}))+m^{2}n\big{)}$.

##  6 Numerical Experiments

In this section, we empirically compare the proposed sparse co-occurring
directions (SCOD) with frequent direction based AMM (FD-AMM) [37], co-
occurring directions (COD) [25] and sparse frequent direction based AMM
algorithm (SFD-AMM)111SFD-AMM refers to the method simply replacing FD step in
FD-AMM with sparse frequent directions [12]. We provide more detailed
discussion about SFD-AMM in appendix.. Our experiments are conducted on a
desktop computer with Intel(R) Core(TM) i5-4570 CPU and 24GB memory. We use
MATLAB 2019a to run the experiments and the operating system is Windows 10.

Instead of increasing $q_{i}$ logarithmically as the analysis of Theorem 2, we
fix $q_{i}=5$ in our experiment since the empirical error arise from SPM
(Algorithm 3) is very small in practice. In the implementation, powering ${\bf
M}{\bf M}^{\top}$ in Algorithm 3 makes $({\bf M}{\bf M}^{\top})^{q}{\bf M}{\bf
G}$ could be ill-conditioned. We include an additionally orthonormalization
step after each round of multiplications to improve the stability [23, 26].
This operation does not change the column span, so it gives an equivalent
algorithm in exact arithmetic, but improves empirical performance
significantly. Since $q$ is typical a small constant in practice, the
additionally cost of orthonormalization is limited.

Dataset | $n$ | $d_{x}$ | $d_{y}$ | ${\rm density}({\bf X})$ | ${\rm density}({\bf Y})$
---|---|---|---|---|---
APR (EN-FR) | $2.32\times 10^{4}$ | $2.80\times 10^{4}$ | $4.28\times 10^{4}$ | $6.31\times 10^{-4}$ | $4.53\times 10^{-4}$
PAN (EN-FR) | $8.90\times 10^{4}$ | $5.12\times 10^{4}$ | $9.96\times 10^{4}$ | $4.38\times 10^{-4}$ | $2.43\times 10^{-4}$
JRC (EN-FR) | $1.50\times 10^{5}$ | $1.72\times 10^{5}$ | $1.87\times 10^{5}$ | $1.65\times 10^{-4}$ | $1.64\times 10^{-4}$
JRC (EN-ES) | $1.50\times 10^{5}$ | $1.72\times 10^{5}$ | $1.92\times 10^{5}$ | $1.65\times 10^{-4}$ | $1.60\times 10^{-4}$
JRC (FR-ES) | $1.50\times 10^{5}$ | $1.87\times 10^{5}$ | $1.92\times 10^{5}$ | $1.64\times 10^{-4}$ | $1.60\times 10^{-4}$
EURO (EN-FR) | $4.76\times 10^{5}$ | $7.25\times 10^{4}$ | $8.77\times 10^{4}$ | $3.46\times 10^{-4}$ | $3.65\times 10^{-4}$
EURO (EN-ES) | $4.76\times 10^{5}$ | $7.25\times 10^{4}$ | $8.80\times 10^{4}$ | $3.46\times 10^{-4}$ | $3.47\times 10^{-4}$
EURO (FR-ES) | $4.76\times 10^{5}$ | $8.77\times 10^{4}$ | $8.80\times 10^{4}$ | $3.65\times 10^{-4}$ | $3.47\times 10^{-4}$
Table 1: We present the size and density of datasets used in our experiments,
where ${\rm density}({\bf X})={\rm nnz}({\bf X})/nd_{x}$ and ${\rm
density}({\bf Y})={\rm nnz}({\bf Y})/nd_{y}$. All of these datasets are
publicly available [8], which can be downloaded from:
https://github.com/FerreroJeremy/Cross-Language-Dataset.
|  |  |
---|---|---|---
(a) APR (EN-FR) | (b) PAN (EN-FR) | (d) JRC (EN-FR) | (e) JRC (FR-ES)
|  |  |
(f) JRC (FR-ES) | (g) EURO (EN-FR) | (j) EURO (EN-ES) | (i) EURO (FR-ES)
Figure 1: The plot of sketch size against relative spectral norm error

We evaluate performance of all algorithms on cross-language datasets: Amazon
Product Reviews (APR), PAN-PC-11 (PAN), JRC Acquis (JRC) and Europarl (EURO)
which contain millions of English (EN), French (FR) and Spanish (ES) sentences
[29, 27, 28, 17]. Each of dataset has alignment information of two languages
at sentence-level and there are $n$ sentences in total. We let $t$-th row of
${\bf X}$ be the bag-of-words feature of $t$-th sentence with respect to one
language and $t$-th row of ${\bf Y}$ be the bag-of-words feature of the same
sentence respect to the other language. All of input matrices are large but
very sparse and we summary the parameters in Table 1.

We demonstrate sketch-error and time-error comparisons in Figure 1 and 2
respectively. It is apparently that SCOD always performs better than all
baseline algorithms. We present sketch-time comparison in Figure 3. The
algorithms SFD-AMM and SCOD are much more faster than FD-AMM and COD, since
FD-AMM and COD ignore the sparse structure of the input matrices. The running
time of SFD-AMM and SCOD are comparable which satisfies our complexity
analysis.

|  |  |
---|---|---|---
(a) APR (EN-FR) | (b) PAN (EN-FR) | (d) JRC (EN-FR) | (e) JRC (FR-ES)
|  |  |
(f) JRC (FR-ES) | (g) EURO (EN-FR) | (j) EURO (EN-ES) | (i) EURO (FR-ES)
Figure 2: The plot of time (s) against relative spectral norm error |  |  |
---|---|---|---
(a) APR (EN-FR) | (b) PAN (EN-FR) | (d) JRC (EN-FR) | (e) JRC (FR-ES)
|  |  |
(f) JRC (FR-ES) | (g) EURO (EN-FR) | (j) EURO (EN-ES) | (i) EURO (FR-ES)
Figure 3: The plot of sketch size against time (s)

##  7 Conclusion

In this paper, we first improved the error bound of a deterministic sketching
algorithm COD for streaming AMM problem. In symmetric case, our result matches
the error bound of classical algorithm FD. We also proved COD matches the
space lower bound complexity to achieve our error bound. In addition, we
proposed a sparse variant of COD with a reasonable error bound. The
experimental results show that the proposed algorithm has better performance
than baseline methods in practice.

It would be interesting to borrow the idea of this paper to establish better
theoretical guarantees and streaming algorithms for more classical machine
learning and statistical models such as canonical correlation analysis [14, 1,
37], generalized eigenvector decomposition [2, 9] and spectral co-clustering
[6].

## References

  * Avron et al. [2013] Haim Avron, Christos Boutsidis, Sivan Toledo, and Anastasios Zouzias.  Efficient dimensionality reduction for canonical correlation analysis.  In _ICML_ , 2013.
  * Bhatia et al. [2018] Kush Bhatia, Aldo Pacchiano, Nicolas Flammarion, Peter L. Bartlett, and Michael I. Jordan.  Gen-Oja: Simple & efficient algorithm for streaming generalized eigenvector computation.  In _NIPS_ , 2018.
  * Chen et al. [2020] Cheng Chen, Luo Luo, Weinan Zhang, Yong Yu, and Yijiang Lian.  Efficient and robust high-dimensional linear contextual bandits.  In _IJCAI_ , 2020.
  * Clarkson and Woodruff [2017] Kenneth L. Clarkson and David P. Woodruff.  Low-rank approximation and regression in input sparsity time.  _Journal of the ACM_ , 63(6):1–45, 2017.
  * Desai et al. [2016] Amey Desai, Mina Ghashami, and Jeff M. Phillips.  Improved practical matrix sketching with guarantees.  _IEEE Transactions on Knowledge and Data Engineering_ , 28(7):1678–1690, 2016.
  * Dhillon [2001] Inderjit S. Dhillon.  Co-clustering documents and words using bipartite spectral graph partitioning.  In _SIGKDD_ , 2001.
  * Drineas et al. [2006] Petros Drineas, Ravi Kannan, and Michael W. Mahoney.  Fast monte carlo algorithms for matrices I: Approximating matrix multiplication.  _SIAM Journal on Computing_ , 36(1):132–157, 2006\.
  * Ferrero et al. [2016] Jérémy Ferrero, Frédéric Agnes, Laurent Besacier, and Didier Schwab.  A multilingual, multi-style and multi-granularity dataset for cross-language textual similarity detection.  In _LREC_ , 2016.
  * Gene and Charles [1996] H GOLUB Gene and F Charles.  Matrix computations.  _Johns Hopkins Universtiy Press, 3rd edtion_ , 1996.
  * Ghashami and Phillips [2014] Mina Ghashami and Jeff M. Phillips.  Relative errors for deterministic low-rank matrix approximations.  In _SODA_ , 2014.
  * Ghashami et al. [2016a] Mina Ghashami, Edo Liberty, and Jeff M. Phillips.  Efficient frequent directions algorithm for sparse matrices.  In _SIGKDD_ , 2016a.
  * Ghashami et al. [2016b] Mina Ghashami, Edo Liberty, Jeff M. Phillips, and David P. Woodruff.  Frequent directions: Simple and deterministic matrix sketching.  _SIAM Journal on Computing_ , 45(5):1762–1792, 2016b.
  * Horn and Johnson [1994] Roger A. Horn and Charles R. Johnson.  _Topics in matrix analysis_.  Cambridge university press, 1994.
  * Hotelling [1992] Harold Hotelling.  Relations between two sets of variates.  In _Breakthroughs in statistics_ , pages 162–190. Springer, 1992\.
  * Huang [2019] Zengfeng Huang.  Near optimal frequent directions for sketching dense and sparse matrices.  _Journal of Machine Learning Research_ , 20(56):1–23, 2019.
  * Kapralov and Talwar [2013] Michael Kapralov and Kunal Talwar.  On differentially private low rank approximation.  In _SODA_ , 2013.
  * Koehn [2005] Philipp Koehn.  Europarl: A parallel corpus for statistical machine translation.  In _MT summit_ , volume 5, pages 79–86. Citeseer, 2005.
  * Kuzborskij et al. [2019] Ilja Kuzborskij, Leonardo Cella, and Nicolò Cesa-Bianchi.  Efficient linear bandits through matrix sketching.  In _AISTATS_ , 2019.
  * Liberty [2013] Edo Liberty.  Simple and deterministic matrix sketching.  In _SGIKDD_ , 2013.
  * Luo et al. [2016] Haipeng Luo, Alekh Agarwal, Nicolo Cesa-Bianchi, and John Langford.  Efficient second order online learning by sketching.  In _NIPS_ , 2016.
  * Luo et al. [2018] Luo Luo, Wenpeng Zhang, Zhihua Zhang, Wenwu Zhu, Tong Zhang, and Jian Pei.  Sketched follow-the-regularized-leader for online factorization machine.  In _SIGKDD_ , 2018.
  * Luo et al. [2019] Luo Luo, Cheng Chen, Zhihua Zhang, Wu-Jun Li, and Tong Zhang.  Robust frequent directions with application in online learning.  _Journal of Machine Learning Research_ , 20(45):1–41, 2019.
  * Martinsson et al. [2010] Per-Gunnar Martinsson, Arthur Szlam, Mark Tygert, et al.  Normalized power iterations for the computation of svd.  In _NIPS Workshop on Low-rank Methods for Large-scale Machine Learning_ , 2010.
  * Misra and Gries [1982] Jayadev Misra and David Gries.  Finding repeated elements.  _Science of computer programming_ , 2(2):143–152, 1982.
  * Mroueh et al. [2016] Youssef Mroueh, Etienne Marcheret, and Vaibhava Goel.  Co-Occuring directions sketching for approximate matrix multiply.  In _AISTATS_ , 2016.
  * Musco and Musco [2015] Cameron Musco and Christopher Musco.  Randomized block Krylov methods for stronger and faster approximate singular value decomposition.  In _NIPS_ , 2015.
  * Potthast et al. [2010] Martin Potthast, Benno Stein, Alberto Barrón-Cedeño, and Paolo Rosso.  An evaluation framework for plagiarism detection.  In _COLING_ , 2010.
  * Potthast et al. [2011] Martin Potthast, Alberto Barrón-Cedeño, Benno Stein, and Paolo Rosso.  Cross-language plagiarism detection.  _Language Resources and Evaluation_ , 45(1):45–62, 2011.
  * Prettenhofer and Stein [2010] Peter Prettenhofer and Benno Stein.  Cross-language text classification using structural correspondence learning.  In _ACL_ , 2010.
  * Rudelson and Vershynin [2010] Mark Rudelson and Roman Vershynin.  Non-asymptotic theory of random matrices: extreme singular values.  In _ICM_ , 2010.
  * Sarlos [2006] Tamas Sarlos.  Improved approximation algorithms for large matrices via random projections.  In _FOCS_ , 2006.
  * Shi and Phillips [2020] Benwei Shi and Jeff M. Phillips.  A deterministic streaming sketch for ridge regression.  _arXiv preprint:2002.02013_ , 2020.
  * Srebro et al. [2005] Nathan Srebro, Jason Rennie, and Tommi S Jaakkola.  Maximum-margin matrix factorization.  In _Advances in neural information processing systems_ , pages 1329–1336, 2005.
  * Wegelin [2000] Jacob A. Wegelin.  A survey of partial least squares (PLS) methods, with emphasis on the two-block case.  _University of Washington, Technical Report_ , 2000.
  * Woodruff [2014a] David P. Woodruff.  Sketching as a tool for numerical linear algebra.  _Foundations and Trends in Theoretical Computer Science_ , 10(1-2):1–157, 2014a.
  * Woodruff [2014b] David P. Woodruff.  Low rank approximation lower bounds in row-update streams.  In _NIPS_ , 2014b.
  * Ye et al. [2016] Qiaomin Ye, Luo Luo, and Zhihua Zhang.  Frequent direction algorithms for approximate matrix multiplication with applications in CCA.  In _IJCAI_ , 2016.

##  Appendix A The Proof of Lemma 3

This lemma can be proved by the analysis of Mroueh et al.’s (2016) Theorem 2.
We reformulate the details by our notations here for completeness.

###### Proof.

Let ${\widehat{\bf A}}^{(t)}=({\bf{\Sigma}}^{(t)})^{1/2}{\bf
U}^{{(t)}\top}{\bf Q}_{x}^{{(t)}\top}$ and ${\widehat{\bf
B}}^{(t)}=({\bf{\Sigma}}^{(t)})^{1/2}{\bf V}^{{(t)}\top}{\bf
Q}_{y}^{{(t)}\top}$. The first inequality can be proved as follows:

$\displaystyle\left\|{\bf X}^{\top}{\bf Y}-{\bf A}^{\top}{\bf B}\right\|_{2}$
$\displaystyle=$
$\displaystyle\left\|\sum_{t=1}^{n}{\bf{x}}^{(t)}({\bf{y}}^{(t)})^{\top}-\sum_{t=1}^{n}\left(\big{(}{\bf
A}^{(t)}\big{)}^{\top}{\bf B}^{(t)}-\big{(}{\bf A}^{(t-1)}\big{)}^{\top}{\bf
B}^{(t-1)}\right)\right\|_{2}$ $\displaystyle=$
$\displaystyle\left\|\sum_{t=1}^{n}{\bf{x}}^{(t)}({\bf{y}}^{(t)})^{\top}-\sum_{t=1}^{n}\left(\big{(}{\bf
A}^{(t)}\big{)}^{\top}{\bf B}^{(t)}-\big{(}{\bf A}^{(t-1)}\big{)}^{\top}{\bf
B}^{(t-1)}\right)\right\|_{2}$ $\displaystyle=$
$\displaystyle\left\|\sum_{t=1}^{n}\left(\big{(}{\widehat{\bf
A}}^{(t)}\big{)}^{\top}{\widehat{\bf B}}^{(t)}-\big{(}{\bf
A}^{(t)}\big{)}^{\top}{\bf B}^{(t)}\right)\right\|_{2}$ $\displaystyle\leq$
$\displaystyle\sum_{t=1}^{n}\left\|\big{(}{\widehat{\bf
A}}^{(t)}\big{)}^{\top}{\widehat{\bf B}}^{(t)}-\big{(}{\bf
A}^{(t)}\big{)}^{\top}{\bf B}^{(t)}\right\|_{2}$ $\displaystyle\leq$
$\displaystyle\sum_{t=1}^{n}\left\|{\bf Q}_{x}^{(t)}{\bf
U}^{(t)}\left({\bf\Sigma}^{(t)}-\delta^{{(t)}}{\bf
I}_{2m}-{\bf\Sigma}^{(t)}\right)({\bf V}^{(t)})^{\top}({\bf
Q}_{y}^{(t)})^{\top}\right\|_{2}$ $\displaystyle=$
$\displaystyle\sum_{t=1}^{n}\delta^{{(t)}},$

where we use triangle inequality and the definition of the notations.

Then we show the second inequality. Similar to above analysis, we have

$\displaystyle\big{\|}{\bf A}^{\top}{\bf B}\big{\|}_{*}$ $\displaystyle=$
$\displaystyle\sum_{t=1}^{n}\left(\big{\|}({\bf A}^{(t)})^{\top}{\bf
B}^{(t)}\big{\|}_{*}-\big{\|}({\bf A}^{(t-1)})^{\top}{\bf
B}^{(t-1)}\big{\|}_{*}\right)$ $\displaystyle=$
$\displaystyle\sum_{t=1}^{n}\left(\left\|({\widehat{\bf
A}}^{(t)})^{\top}{\widehat{\bf B}}^{(t)}\right\|_{*}-\left\|({\bf
A}^{(t-1)})^{\top}{\bf
B}^{(t-1)}\right\|_{*}\right)-\sum_{t=1}^{n}\left(\left\|({\widehat{\bf
A}}^{(t)})^{\top}{\widehat{\bf B}}^{(t)}\right\|_{*}-\left\|({\bf
A}^{(t)})^{\top}{\bf B}^{(t)}\right\|_{*}\right).$

The QR steps means $({\widehat{\bf A}}^{(t)})^{\top}{\widehat{\bf B}}^{(t)}$
and $({\bf A}^{(t)})^{\top}{\bf B}^{(t)}$ can be written as

$\displaystyle({\widehat{\bf A}}^{(t)})^{\top}{\widehat{\bf B}}^{(t)}={\bf
Q}_{x}^{(t)}{\bf U}^{(t)}{\bf\Sigma}^{(t)}({\bf V}^{(t)})^{\top}{\bf
Q}_{y}^{(t)}\text{~{}~{}and~{}~{}}({\bf A}^{(t)})^{\top}{\bf B}^{(t)}={\bf
Q}_{x}^{(t)}{\bf U}^{(t)}({\bf\Sigma}-\sigma_{m}^{(t)}{\bf I}_{m})^{(t)}({\bf
V}^{(t)})^{\top}{\bf Q}_{y}^{(t)}.$

which implies $\big{\|}({\widehat{\bf A}}^{(t)})^{\top}{\widehat{\bf
B}}^{(t)}\big{\|}_{*}-\left\|({\bf A}^{(t)})^{\top}{\bf
B}^{(t)}\right\|_{*}=m\sigma_{m}^{(t)}$.

Using triangle inequality, we have

$\displaystyle\left\|({\widehat{\bf A}}^{(t)})^{\top}{\widehat{\bf
B}}^{(t)}\right\|_{*}-\left\|({\bf A}^{(t-1)})^{\top}{\bf
B}^{(t-1)}\right\|_{*}$ $\displaystyle\leq$
$\displaystyle\left\|({\widehat{\bf A}}^{(t)})^{\top}{\widehat{\bf
B}}^{(t)}-({\bf A}^{(t-1)})^{\top}{\bf B}^{(t-1)}\right\|_{*}$
$\displaystyle=$
$\displaystyle\big{\|}{\bf{x}}^{(t)}({\bf{y}}^{(t)})^{\top}\big{\|}_{*}=\big{\|}{\bf{x}}^{(t)}\big{\|}_{2}\big{\|}{\bf{y}}^{(t)}\big{\|}_{2}$

Combing all above results, we have

$\displaystyle\left\|{\bf A}^{\top}{\bf
B}\right\|_{*}\leq\sum_{t=1}^{n}\left\|{\bf{x}}^{(t)}\right\|_{2}\left\|{\bf{y}}^{(t)}\right\|_{2}-m\sum_{t=1}^{n}\sigma_{m}^{(t)}\leq\left\|{\bf
X}\right\|_{F}\left\|{\bf Y}\right\|_{F}-m\sum_{t=1}^{n}\sigma_{m}^{(t)}.$

∎

##  Appendix B The Proof of Lemma 4

This lemma is crucial to establish the tighter bound of COD. We first
introduce the following property of Ky Fan $k$-norm.

######  Lemma 10 (Horn and Johnson 13, Theorem 3.4.1).

Given matrix ${\bf M}$, we have

$\displaystyle\left\|{\bf M}\right\|_{k}=\max\left\\{\left|{\rm tr}\left({\bf
P}^{\top}{\bf M}{\bf Q}\right)\right|:{\bf P}^{\top}{\bf P}={\bf I}_{k},{\bf
Q}^{\top}{\bf Q}={\bf I}_{k}\right\\},$

where ${\rm tr}(\cdot)$ is the trace of the matrix.

Then we prove Lemma 4 by using Lemma 10.

###### Proof.

We let $({\widehat{\bf P}},{\widehat{\bf
Q}})=\operatorname*{arg\,max}\left\\{\left|{\rm tr}\left({\bf P}^{\top}{\bf
A}{\bf Q}\right)\right|:{\bf P}^{\top}{\bf P}={\bf I},{\bf Q}^{\top}{\bf
Q}={\bf I}\right\\}$. Then we have

$\displaystyle\left\|{\bf X}^{\top}{\bf Y}\right\|_{*}-\left\|{\bf
A}^{\top}{\bf B}\right\|_{*}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{k}\sigma_{i}({\bf X}^{\top}{\bf
Y})+\sum_{i=k+1}^{d}\sigma_{i}({\bf X}^{\top}{\bf Y})-\left\|{\bf
A}^{\top}{\bf B}\right\|_{*}$ $\displaystyle=$
$\displaystyle\max\left\\{\left|{\rm tr}\left({\bf P}^{\top}({\bf
X}^{\top}{\bf Y}){\bf Q}\right)\right|:{\bf P}^{\top}{\bf P}={\bf I}_{k},{\bf
Q}^{\top}{\bf Q}={\bf I}_{k}\right\\}+\sum_{i=k+1}^{d}\sigma_{i}({\bf
X}^{\top}{\bf Y})$ $\displaystyle-\max\left\\{\left|{\rm tr}\left({\bf
P}^{\top}({\bf A}^{\top}{\bf B}){\bf Q}\right)\right|:{\bf P}^{\top}{\bf
P}={\bf I}_{k},{\bf Q}^{\top}{\bf Q}={\bf
I}_{k}\right\\}-\sum_{i=k+1}^{d}\sigma_{i}({\bf A}^{\top}{\bf B})$
$\displaystyle\leq$ $\displaystyle\sum_{i=k+1}^{d}\sigma_{i}({\bf
X}^{\top}{\bf Y})+\left|{\rm tr}\left({\widehat{\bf P}}^{\top}\left({\bf
X}^{\top}{\bf Y}\right){\widehat{\bf Q}}\right)\right|-\left|{\rm
tr}\left({\widehat{\bf P}}^{\top}\left({\bf A}^{\top}{\bf
B}\right){\widehat{\bf Q}}\right)\right|$ $\displaystyle\leq$
$\displaystyle\sum_{i=k+1}^{d}\sigma_{i}({\bf X}^{\top}{\bf Y})+\left|{\rm
tr}\left({\widehat{\bf P}}^{\top}\left({\bf X}^{\top}{\bf Y}-{\bf
A}^{\top}{\bf B}\right){\widehat{\bf Q}}\right)\right|$ $\displaystyle\leq$
$\displaystyle\sum_{i=k+1}^{d}\sigma_{i}({\bf X}^{\top}{\bf
Y})+\max\left\\{\left|{\rm tr}\left({\bf P}^{\top}\left({\bf X}^{\top}{\bf
Y}-{\bf A}^{\top}{\bf B}\right){\bf Q}\right)\right|:{\bf P}^{\top}{\bf
P}={\bf I}_{k},{\bf Q}^{\top}{\bf Q}={\bf I}_{k}\right\\}$ $\displaystyle\leq$
$\displaystyle\sum_{i=k+1}^{d}\sigma_{i}({\bf X}^{\top}{\bf
Y})+\max\left\\{\left|{\rm tr}\left({\bf
P}^{\top}\left(\sum_{t=1}^{n}\sigma_{m}^{(t)}{\bf Q}_{x}^{(t)}{\bf
U}^{(t)}{\bf V}^{{(t)}\top}{\bf Q}_{y}^{{(t)}\top}\right){\bf
Q}\right)\right|:{\bf P}^{\top}{\bf P}={\bf I}_{k},{\bf Q}^{\top}{\bf Q}={\bf
I}_{k}\right\\}$ $\displaystyle\leq$
$\displaystyle\sum_{i=k+1}^{d}\sigma_{i}({\bf X}^{\top}{\bf
Y})+\sum_{t=1}^{n}\delta^{(t)}\max\left\\{\left|{\rm tr}\left({\bf
P}^{\top}\left({\bf Q}_{x}^{(t)}{\bf U}^{(t)}{\bf V}^{{(t)}\top}{\bf
Q}_{y}^{{(t)}\top}\right){\bf Q}\right)\right|:{\bf P}^{\top}{\bf P}={\bf
I}_{k},{\bf Q}^{\top}{\bf Q}={\bf I}_{k}\right\\}$ $\displaystyle=$
$\displaystyle\sum_{i=k+1}^{d}\sigma_{i}({\bf X}^{\top}{\bf
Y})+\sum_{t=1}^{n}\delta^{(t)}\sum_{i=1}^{k}\sigma_{i}\left({\bf
Q}_{x}^{(t)}{\bf U}^{(t)}{\bf V}^{{(t)}\top}{\bf Q}_{y}^{{(t)}\top}\right)$
$\displaystyle=$ $\displaystyle\sum_{i=k+1}^{d}\sigma_{i}({\bf X}^{\top}{\bf
Y})+k\sum_{t=1}^{n}\delta^{(t)},$

where the first inequality is due to the definition of ${\widehat{\bf P}}$ and
${\widehat{\bf Q}}$; the second and the third one use triangle inequality; the
last to inequality is based on the procedure of the algorithm; all equalities
come from Lemma 10 and the procedure of COD. ∎

##  Appendix C The Proof of Lemma 7

We can prove Lemma 7 by modifying the analysis in Section 4.3 of Woodruff’s
(2014a) survey. We present the details for completeness. The proof is based on
the following lemma.

######  Lemma 11 (Rudelson and Vershynin 30, Proposition 2.4 and (3.2)).

Let ${\bf{\Omega}}\in{\mathbb{R}}^{d_{1}\times d_{2}}$ be a random matrix
whose entries are independent mean zero sub-gaussian random variables whose
subgaussian moments are bounded by 1. Then we have

  1. 1.

${\mathbb{P}}\big{(}\left\|{\bf{\Omega}}\right\|_{2}>C(\sqrt{d_{1}}+\sqrt{d_{2}}+t)\big{)}\leq
2\exp(-ct^{2})$ for any $t>0$;

  2. 2.

${\mathbb{P}}(\sigma_{\min}({\bf{\Omega}})\leq\zeta d^{-1/2})\leq\zeta$ when
$d_{1}=d_{2}=d$ for any $\zeta>0$;

where $c>0$ and $C>0$ are some constants.

Then we provide the proof of Lemma 7.

###### Proof.

Let ${\bf N}=({\bf M}{\bf M}^{\top})^{q}{\bf M}$. By Woodruff’s (2014a) Lemma
4.14, ${\bf Z}{\bf Z}^{\top}{\bf M}$ is the best rank-$m$ approximation of
${\bf M}$ in the column space of ${\bf Z}$ with respect to the spectral norm.
Hence, we have

$\displaystyle\left\|{\bf M}-{\bf Z}{\bf Z}^{\top}{\bf
M}\right\|_{2}\leq\left\|{\bf M}-({\bf Z}{\bf Z}^{\top}{\bf N})({\bf Z}{\bf
Z}^{\top}{\bf N})^{\dagger}{\bf M}\right\|_{2}\leq\left\|\left({\bf
I}_{d_{1}}-({\bf Z}{\bf Z}^{\top}{\bf N})({\bf Z}{\bf Z}^{\top}{\bf
N})^{\dagger}\right){\bf M}\right\|_{2},$

where the notation $(\cdot)^{\dagger}$ presents pseudo-inverse; the inequality
follows ${\bf Z}{\bf Z}^{\top}{\bf N}$ is of rank-$m$ and in the column space
of ${\bf Z}$.

Since ${\bf I}_{d_{1}}-({\bf Z}{\bf Z}^{\top}{\bf N})({\bf Z}{\bf
Z}^{\top}{\bf N})^{\dagger}$ is a projection matrix, we can apply Woodruff’s
(2014a) Lemma 4.15 to infer that

$\displaystyle\left\|({\bf I}_{d_{1}}-({\bf Z}{\bf Z}^{\top}{\bf N})({\bf
Z}{\bf Z}^{\top}{\bf N})^{\dagger}){\bf M}\right\|_{2}$ $\displaystyle\leq$
$\displaystyle\left\|({\bf I}_{d_{1}}-({\bf Z}{\bf Z}^{\top}{\bf N})({\bf
Z}{\bf Z}^{\top}{\bf N})^{\dagger})({\bf M}{\bf M}^{\top})^{q}{\bf
M}\right\|_{2}^{1/(2q+1)}$ $\displaystyle=$ $\displaystyle\left\|{\bf N}-({\bf
Z}{\bf Z}^{\top}{\bf N})({\bf Z}{\bf Z}^{\top}{\bf N})^{\dagger}{\bf
N}\right\|_{2}^{1/(2q+1)}$ $\displaystyle=$ $\displaystyle\left\|{\bf N}-{\bf
Z}{\bf Z}^{\top}{\bf N}\right\|_{2}^{1/(2q+1)}$

where we use that $({\bf Z}{\bf Z}^{\top}{\bf N})^{\dagger}=({\bf
Z}^{\top}{\bf N})^{\dagger}{\bf Z}^{\top}$ since ${\bf Z}$ has orthonormal
columns, and thus

$\displaystyle({\bf Z}{\bf Z}^{\top}{\bf N})({\bf Z}{\bf Z}^{\top}{\bf
N})^{\dagger}{\bf N}=({\bf Z}{\bf Z}^{\top}{\bf N})({\bf Z}^{\top}{\bf
N})^{\dagger}({\bf Z}^{\top}{\bf N})={\bf Z}{\bf Z}^{\top}{\bf N}.$

Hence, we have

$\displaystyle\left\|{\bf M}-{\bf Z}{\bf Z}^{T}{\bf
M}\right\|_{2}\leq\left\|{\bf N}-{\bf Z}{\bf Z}^{\top}{\bf
N}\right\|_{2}^{1/(2q+1)}.$ (7)

Let ${\bf U}{\bf{\Sigma}}{\bf V}^{\top}$ be the SVD of ${\bf N}$,
${\bf{\Omega}}_{U}={\bf V}_{m}^{\top}{\bf G}\in{\mathbb{R}}^{m\times m}$ and
${\bf{\Omega}}_{L}={\bf V}^{\top}_{d_{1}-m}{\bf
G}\in{\mathbb{R}}^{(d_{1}-m)\times m}$, where ${\bf V}^{\top}_{m}$ denotes the
top $m$ rows of ${\bf V}^{\top}$ and ${\bf V}^{\top}_{d_{1}-m}$ the remaining
rows. Since ${\bf V}^{\top}$ are column orthonormal, by rotational invariance
of the Gaussian distribution, both ${\bf{\Omega}}_{U}$ and ${\bf{\Omega}}_{L}$
are independent matrices of i.i.d. ${\mathcal{N}}(0,1)$ entries.

We now apply Woodruff’s (2014a) Lemma 4.4 with the ${\bf C}$ of that lemma
equal to ${\bf Z}$ above, the ${\bf Z}$ of that lemma equal to ${\bf V}_{m}$,
and the ${\bf A}$ of that lemma equal to ${\bf N}$ above. This implies the
${\bf E}$ of that lemma is equal to ${\bf N}-{\bf N}_{m}$. Note that to apply
the lemma we need ${\bf V}^{T}_{m}{\bf G}$ to have full rank, which holds with
probability 1 since it is a $m\times m$ matrix of i.i.d. ${\mathcal{N}}(0,1)$
random variables. We thus have

$\displaystyle\begin{split}&\left\|{\bf N}-{\bf Z}{\bf Z}^{\top}{\bf
N}\right\|_{2}^{2}\\\ =&\left\|{\bf N}-{\bf
N}_{m}\right\|_{2}^{2}+\left\|({\bf N}-{\bf N}_{m}){\bf G}({\bf
V}_{m}^{\top}{\bf G})^{\dagger}\right\|_{2}^{2}\\\ =&\left\|{\bf N}-{\bf
N}_{m}\right\|_{2}^{2}+\left\|{\bf U}_{d_{1}-m}{\bf{\Sigma}}_{d_{1}-m}{\bf
V}_{d_{1}-m}^{\top}{\bf G}({\bf V}_{m}^{\top}{\bf
G})^{\dagger}\right\|_{2}^{2}\\\ =&\left\|{\bf N}-{\bf
N}_{m}\right\|_{2}^{2}+\left\|{\bf{\Sigma}}_{d_{1}-m}{\bf
V}_{d_{1}-m}^{\top}{\bf G}({\bf V}_{m}^{\top}{\bf
G})^{\dagger}\right\|_{2}^{2}\\\ \leq&\left\|{\bf N}-{\bf
N}_{m}\right\|_{2}^{2}\left(1+\left\|{\bf{\Omega}}_{L}\right\|_{2}^{2}\big{\|}{\bf{\Omega}}_{U}^{\dagger}\big{\|}^{2}\right)\end{split}$
(8)

where ${\bf{\Sigma}}_{d_{1}-m}$ denotes the $(d_{1}-m)\times(d_{1}-m)$
diagonal matrix whose entries are the bottom $d_{1}-m$ diagonal entries of
${\bf{\Sigma}}$, and ${\bf U}_{d_{1}-m}$ denotes the rightmost $d_{1}-m$
columns of ${\bf U}$. Here in the second equality we use unitary invariance of
${\bf U}_{d_{1}-m}$, while in the inequality we use sub-multiplicativity of
the spectral norm.

By using Lemma 11 with ${\bf{\Omega}}={\bf{\Omega}}_{L}$ and
$t=\sqrt{c^{-1}\log(4/p)}$, we have

$\displaystyle{\mathbb{P}}\left(\left\|{\bf{\Omega}}_{L}\right\|_{2}^{2}\leq\left(\sqrt{d_{1}-m}+\sqrt{m}+\sqrt{c^{-1}\log\left(4/p\right)}\right)^{2}\right)\geq
1-\frac{p}{2}.$ (9)

By using Lemma 11 with ${\bf{\Omega}}={\bf{\Omega}}_{U}$ and $\zeta=p/2$, we
have

$\displaystyle{\mathbb{P}}\left(\sigma_{\min}^{2}({\bf{\Omega}}_{U})\geq\frac{p^{2}}{4m}\right)\geq
1-\frac{p}{2}.$ (10)

Since ${\bf{\Omega}}_{L}$ and ${\bf{\Omega}}_{U}$ are independent, combing
inequalities (9) and (10), we have

$\displaystyle
1+\left\|{\bf{\Omega}}_{L}\right\|_{2}^{2}\|{\bf{\Omega}}_{U}^{\dagger}\|^{2}\leq$
$\displaystyle
1+\left(\sqrt{d_{1}-m}+\sqrt{m}+\sqrt{c^{-1}\log\left(4/p\right)}\right)^{2}\cdot\frac{4m}{p^{2}}\leq\frac{c_{0}(d_{1}+\log(1/p))m}{p^{2}}$
(11)

for some constant $c_{0}>0$ with probability at least $(1-p/2)^{2}>1-p$.

Combining results of (7), (LABEL:bound:MNZ) and (11), we have

$\displaystyle\left\|{\bf M}-{\bf Z}{\bf Z}^{T}{\bf
M}\right\|_{2}\leq\left\|{\bf N}-{\bf
N}_{m}\right\|_{2}^{1/(2q+1)}\cdot\left(\frac{c_{0}(d_{1}+\log(1/p))m}{p^{2}}\right)^{1/(4q+2)}.$

Noting that $\left\|{\bf N}-{\bf N}_{m}\right\|_{2}=\left\|{\bf M}-{\bf
M}_{m}\right\|_{2}^{2q+1}$ and setting

$\displaystyle
q=\frac{1}{4}\left(\frac{1}{{\varepsilon}}\log\left(\frac{c_{0}(d_{1}+\log(1/p))m}{p^{2}}\right)-2\right)=\tilde{\Theta}\left(\frac{1}{{\varepsilon}}\log\left(\frac{md_{1}}{p}\right)\right)$

we have

$\displaystyle\left\|{\bf M}-{\bf Z}{\bf Z}^{T}{\bf
M}\right\|_{2}\leq(1+{\varepsilon})\left\|{\bf M}-{\bf
M}_{m}\right\|_{2}=(1+{\varepsilon})\sigma_{m+1}({\bf M})$

with probability at least $1-p$. ∎

##  Appendix D The Proof of Lemma 9

###### Proof.

The procedure of Algorithm 4 means $\sigma^{2}_{i}({\widetilde{\bf
X}}^{(i)})=\sigma^{2}_{i}({\widetilde{\bf
Y}}^{(i)})=\sigma_{i}\big{(}{\widetilde{\bf X}}^{{(i)}\top}{\widetilde{\bf
Y}}^{(i)}\big{)}$. Consider that the output ${\bf Z}^{(i)}$ of SPM (Algorithm
3) is column orthonormal, then we have

$\displaystyle\big{\|}{\widetilde{\bf
X}}^{(i)}\big{\|}^{2}_{F}\big{\|}{\widetilde{\bf Y}}^{(i)}\big{\|}_{F}^{2}$
$\displaystyle=$ $\displaystyle\big{\|}{\widetilde{\bf
X}}^{{(i)}\top}{\widetilde{\bf Y}}^{(i)}\big{\|}_{F}^{2}$ $\displaystyle=$
$\displaystyle\big{\|}{\bf Z}^{(i)}{\bf Z}^{{(i)}\top}{\bf
X}^{\prime{(i)}\top}{\bf Y}^{\prime{(i)}}\big{\|}_{F}^{2}$ $\displaystyle=$
$\displaystyle{\rm tr}\left({\bf Y}^{\prime{(i)}\top}{\bf X}^{\prime{(i)}}{\bf
Z}^{(i)}{\bf Z}^{{(i)}\top}{\bf Z}^{(i)}{\bf Z}^{{(i)}\top}{\bf
X}^{\prime{(i)}\top}{\bf Y}^{\prime{(i)}}\right)$ $\displaystyle=$
$\displaystyle{\rm tr}\left({\bf Y}^{\prime{(i)}\top}{\bf X}^{\prime{(i)}}{\bf
Z}^{(i)}{\bf Z}^{{(i)}\top}{\bf X}^{\prime{(i)}\top}{\bf
Y}^{\prime{(i)}}\right)$ $\displaystyle\leq$ $\displaystyle{\rm tr}\left({\bf
Y}^{\prime{(i)}\top}{\bf X}^{\prime{(i)}}{\bf X}^{\prime{(i)}\top}{\bf
Y}^{\prime{(i)}}\right)$ $\displaystyle=$ $\displaystyle\big{\|}{\bf
X}^{\prime{(i)}\top}{\bf Y}^{\prime{(i)}}\big{\|}_{F}^{2}$ $\displaystyle\leq$
$\displaystyle\big{\|}{\bf X}^{\prime{(i)}}\big{\|}_{F}^{2}\big{\|}{\bf
Y}^{\prime{(i)}}\big{\|}_{F}^{2}.$

∎

##  Appendix E Additional Discussion on SFD-AMM

Ghashami et al. [12] proposed a variant of FD for sketching sparse matrices
called sparse frequent directions (SFD). Given input matrix ${\bf
Z}\in{\mathbb{R}}^{n\times d}$, the algorithm output ${\bf
C}\in{\mathbb{R}}^{m\times d}$ such that

$\displaystyle\left\|{\bf Z}^{\top}{\bf Z}-{\bf C}^{\top}{\bf
C}\right\|_{2}\leq\frac{1}{\alpha m-k}\left(\left\|{\bf
Z}\right\|_{F}^{2}-\left\|{\bf Z}_{k}\right\|_{F}^{2}\right)$

with high probability for any $k<\alpha m$, where $\alpha$ is a constant
depends on the accuracy of SPM. SFD requires
$\tilde{\mathcal{O}}\big{(}m\cdot{\rm nnz}({\bf Z})+m^{2}n\big{)}$ time
complexity and ${\mathcal{O}}(md)$ space. The procedure of SFD is similar to
SCOD in the case of ${\bf X}={\bf Y}$, but includes additional shrinking
operation on the output of SPM [11] to apply the “mergeability porperty” of FD
[12, 5] in their analysis.

For streaming AMM with sparse input, it is natural to combine the idea of SFD
with FD-AMM directly which leads to the algorithm sparse FD-AMM (SFD-AMM).
Similar to FD-AMM, SFD-AMM applies SFD on concatenated matrix ${\bf Z}=[{\bf
X},{\bf Y}]\in{\mathbb{R}}^{n\times(d_{1}+d_{2})}$ and its output ${\bf C}$
which can be written as ${\bf C}=[{\bf A},{\bf B}]$, where ${\bf
A}\in{\mathbb{R}}^{n\times d_{x}}$ and ${\bf B}\in{\mathbb{R}}^{n\times
d_{y}}$. Then we use ${\bf A}^{\top}{\bf B}$ to approximate ${\bf
X}^{\top}{\bf Y}$ that satisfies

$\displaystyle\left\|{\bf X}^{\top}{\bf Y}-{\bf A}^{\top}{\bf
B}\right\|_{2}\leq\left\|{\bf Z}^{\top}{\bf Z}-{\bf C}^{\top}{\bf
C}\right\|_{2}\leq\frac{1}{\alpha m-k}\left(\left\|{\bf
Z}\right\|_{F}^{2}-\left\|{\bf Z}_{k}\right\|_{F}^{2}\right).$ (12)

The time complexity of SFD-AMM has the same order as SCOD since ${\rm
nnz}({\bf Z})={\rm nnz}({\bf X})+{\rm nnz}({\bf Y})$. It is not easy to
compare the error bound (12) with SCOD (Theorem 3) in general because there
does not exist simple relationship between the singular values of ${\bf
Z}=[{\bf X},{\bf Y}]$ and ${\bf X}^{\top}{\bf Y}$. However, SCOD always
performs better than SFD-AMM empirically as we observed in Section 6.

In theoretical, we can improve the time complexity of SFD-AMM to achieve the
error bound of (12) by integrating random sampling [15]. However, the
implementation of this strategy requires the value of $k$ is given.
Unfortunately, it is difficult to select a suitable $k$ for streaming setting
in general. In contrast, the value of $k$ in SFD-AMM or SCOD is only for
theoretical analysis and it is no related to the implementation of algorithms.