A Fast Frequent Directions Algorithm for Low Rank Approximation

TL;DR

This paper introduces a fast, randomized frequent directions algorithm that improves low rank approximation efficiency by integrating sparse subspace embedding, demonstrating strong theoretical guarantees and practical performance.

Contribution

It presents a novel fast algorithm combining frequent directions with sparse subspace embedding for improved low rank approximation efficiency.

Findings

Produces accurate low rank approximations with linear sketch size

Demonstrates efficiency on synthetic and real datasets

Effective in network analysis applications

Abstract

Recently a deterministic method, frequent directions (FD) is proposed to solve the high dimensional low rank approximation problem. It works well in practice, but experiences high computational cost. In this paper, we establish a fast frequent directions algorithm for the low rank approximation problem, which implants a randomized algorithm, sparse subspace embedding (SpEmb) in FD. This new algorithm makes use of FD's natural block structure and sends more information through SpEmb to each block in FD. We prove that our new algorithm produces a good low rank approximation with a sketch of size linear on the rank approximated. Its effectiveness and efficiency are demonstrated by the experimental results on both synthetic and real world datasets, as well as applications in network analysis.

Figures

Tables

Table 1. TABLE I: Flop counts with sketch size ℓ ℓ \ell

	SpEmb (Algorithm 1)	FD (Algorithm 2)	SpFD (Algorithm 3) ($q>1$)
1. Form the sketch $SA$ or $B$	$2\text{nnz}\left(A\right)$	$\left[6d{\left(2\ell \right)}^2+20{\left(2\ell \right)}^3\right]\left(\frac{n}{\ell }-1\right)$	$2\text{nnz}\left(A\right)$
2. Construct $V$	$4d{\ell }^2-\frac{4}{3}{\ell }^3$		$+\left[6d{\left(2\ell \right)}^2+20{\left(2\ell \right)}^3\right]\left(q-1\right)$
(through QR for SpEmb)
3. Compute ${\left[AV\right]}_k$ (via SVD)	$2\text{nnz}\left(A\right)\ell +6n{\ell }^2+20{\ell }^3$	$2\text{nnz}\left(A\right)\ell +6n{\ell }^2+20{\ell }^3$	$2\text{nnz}\left(A\right)\ell +6n{\ell }^2+20{\ell }^3$
and obtain ${\left[AV\right]}_k{V}^T$	$+2n\ell k+2nd\ell$	$+2n\ell k+2nd\ell$	$+2n\ell k+2nd\ell$
Total	$2nd\ell +2\text{nnz}\left(A\right)\left(\ell +1\right)+6n{\ell }^2$	$26nd\ell +2\text{nnz}\left(A\right)\ell +166n{\ell }^2$	$2nd\ell +2\text{nnz}\left(A\right)\left(\ell +1\right)+6n{\ell }^2$
	$+2n\ell k+\left(4d+18⁤\frac{2}{3}\ell \right){\ell }^2$	$+2n\ell k-\left(24d+140\ell \right){\ell }^2$	$+2n\ell k+24d\left(q-1\right){\ell }^2$
			$+\left(160q-140\right){\ell }^3$

Table 2. TABLE II: Data statistics.

	No. of	Data	nnz	Approximating	Sketch
	data points ($n$)	dimension ($d$)	%	rank ($k$)	size ($\ell$)
w8a [37]	64700	300	3.88	20	20:10:150
Birds [43]	11788	312	10.09	50	50:10:150
Protein [44]	24387	357	28.20	50	50:10:170
MNIST-all [28]	70000	784	19.14	100	100:10:200
amazon7-small [4, 14]	262144	1500	0.017	100	100:20:300
rcv1-small [30]	47236	3000	0.14	100	100:20:300

Table 3. TABLE III: The top 10 hubs, top 10 authorities and running time used to derive the hub and authority scores for networks Computational Complexity and Death Penalty .

	HITS	$e^{𝒜}$	NormSamp	DCT	SpEmb	FD	SpFD5	SpFD10	SpFD50
	57	57	57	57	57	57	57	57	57
	634	17	17	17	17	644	644	634	644
	644	644	634	21	21	643	17	644	119
	721	643	644	644	529	634	643	643	643
Top Ten Hubs	643	634	643	643	54	106	634	17	106
(Computational Complexity)	554	106	106	54	106	529	106	106	634
	632	119	119	255	119	119	21	119	17
	801	529	721	106	634	86	162	255	86
	640	86	51	119	51	162	529	86	162
	629	162	86	529	668	520	86	529	21
	no. nodes shared w. $e^{𝒜}$		8	7	6	9	9	9	9
	719	1	717	1	1	1	1	1	1
	717	315	2	2	673	315	315	717	315
	727	673	716	45	2	673	719	315	148
	723	148	718	50	49	148	673	673	2
Top Ten Authorities	808	719	719	315	341	719	148	148	717
(Computational Complexity)	735	717	720	433	45	717	717	2	719
	737	2	721	341	729	45	363	719	673
	1	45	722	742	721	98	727	374	727
	722	727	723	673	664	727	737	808	723
	770	534	724	98	148	2	534	727	735
	no. nodes shared w. $e^{𝒜}$		3	5	5	9	8	8	8
Running Time (s)	0.0388	0.5248	0.0090	0.0300	0.0100	0.0814	0.0079	0.0069	0.0062
	210	210	210	210	210	210	210	210	210
	637	637	637	637	637	637	637	637	637
	413	413	413	413	413	413	413	413	413
	1586	1586	1586	1586	552	1586	1586	1586	1586
Top Ten Hubs	552	552	552	542	1586	552	552	552	552
(Death Penalty)	462	462	930	930	545	930	462	462	618
	930	930	542	1275	618	542	618	542	930
	542	542	1275	545	462	462	542	618	542
	618	618	462	80	578	618	930	930	462
	1275	1275	80	1258	1275	1275	1275	1275	1275
	no. nodes shared w. $e^{𝒜}$		9	7	8	10	10	10	10
	4	4	16	4	1	4	4	4	1
	1	1	4	6	4	1	1	1	4
	6	6	14	1	3	6	6	6	6
	7	7	27	27	5	7	16	16	7
Top Ten Authorities	10	10	10	16	27	10	2	7	10
(Death Penalty)	16	16	1	130	21	16	10	2	16
	2	2	13	12	14	3	14	3	3
	3	3	9	513	1632	2	44	10	2
	44	44	3	2	394	44	7	27	27
	27	27	2	394	12	27	384	9	44
	no. nodes shared w. $e^{𝒜}$		7	6	4	10	8	9	10
Running Time (s)	0.0723	5.1974	0.0282	0.1312	0.0431	0.2868	0.0213	0.0197	0.0188

Equations

$${\Vert A-[AV{]}_k{V}^T\Vert }_F^2$$

$$S=\left[\begin{array}{ccc}{S}_1& & \\ & \ddots & \\ & & S_q\end{array}\right]\in {\mathbb{R}}^{q\ell \times n},$$

$$q\ell \ge \frac{k^2+k}{{\epsilon }^2\delta }.$$

$${\Vert U^T(SP{)}^{T}SPU-I\Vert }_2\le \epsilon .$$

$$S=\left[\begin{array}{ccc}{S}_1& & \\ & \ddots & \\ & & S_q\end{array}\right]=\left[\begin{array}{ccc}{\varphi }_1& & \\ & \ddots & \\ & & {\varphi }_q\end{array}\right]\left[\begin{array}{ccc}{D}_1& & \\ & \ddots & \\ & & D_q\end{array}\right]=:\varphi D$$

$$x_{st}=$$

$$=$$

$$\mathbb{E}[x_{st}]=$$

$$=$$

$$\mathbb{E}[x_{st}^2]=$$

$$\mathbb{E}[d_{j_1}{d}_{j_2}{d}_{j_1^{\prime }}{d}_{j_2^{\prime }}]\mathbb{E}[(PU{)}_{j_{1}s}(PU{)}_{j_{2}t}(PU{)}_{j_1^{\prime }s}(PU{)}_{j_2^{\prime }t}].$$

$$\mathbb{E}\left[d_{j_1}{d}_{j_2}{d}_{j_1^{\prime }}{d}_{j_2^{\prime }}\right]$$

$$\sum _{i,i^{\prime }=1}^{q\ell }\sum _{j=1}^n\mathbb{E}[{\varphi }_{ij}^2{\varphi }_{i^{\prime }j}^2]\mathbb{E}[(PU{)}_{js}^2(PU{)}_{jt}^2]$$

$$=$$

$$\mathbb{E}[(PU{)}_{js}(PU{)}_{jt}(PU{)}_{j^{\prime }s}(PU{)}_{j^{\prime }t}]$$

$$=$$

$$\sum _{i,i^{\prime }=1}^{q\ell }\sum _{j\ne j^{\prime }}^n\mathbb{E}[{\varphi }_{ij}^2{\varphi }_{i^{\prime }{j}^{\prime }}^2]$$

$$=$$

$$+\sum _{i\ne i^{\prime },{\mathcal{L}}_i={\mathcal{L}}_{i^{\prime }}}\sum _{j,j^{\prime }\in {\mathcal{L}}_i,j\ne j^{\prime }}\mathbb{E}[{\varphi }_{ij}^2{\varphi }_{i^{\prime }{j}^{\prime }}^2])$$

$$=$$

$$\sum _{i,i^{\prime }=1}^{q\ell }\sum _{j\ne j^{\prime }}^n\mathbb{E}[{\varphi }_{ij}^2{\varphi }_{i^{\prime }{j}^{\prime }}^2]\mathbb{E}[(PU{)}_{js}(PU{)}_{jt}(PU{)}_{j^{\prime }s}(PU{)}_{j^{\prime }t}]$$

$$=$$

$$\mathbb{E}[(PU{)}_{j_{1}s}^2(PU{)}_{j_{2}t}^2]=\frac{1}{n(n-1)}\sum _{j_1\ne j_2}{U}_{j_{1}s}^2{U}_{j_{2}t}^2;$$

$$\sum _{i,i^{\prime }=1}^{q\ell }\sum _{j_1\ne j_2}^n\mathbb{E}[{\varphi }_{i{j}_1}{\varphi }_{i{j}_2}{\varphi }_{i^{\prime }{j}_1}{\varphi }_{i^{\prime }{j}_2}]$$

$$=$$

$$\sum _{i,i^{\prime }=1}^{q\ell }\sum _{j_1\ne j_2}^n\mathbb{E}[{\varphi }_{i{j}_1}{\varphi }_{i{j}_2}{\varphi }_{i^{\prime }{j}_1}{\varphi }_{i^{\prime }{j}_2}]\mathbb{E}[(PU{)}_{j_{1}s}^2(PU{)}_{j_{2}t}^2]$$

$$=$$

$$\mathbb{E}[(PU{)}_{js}(PU{)}_{j^{\prime }t}(PU{)}_{j^{\prime }s}(PU{)}_{jt}]$$

$$=$$

$$\sum _{i,i^{\prime }=1}^{q\ell }\sum _{j\ne j^{\prime }}^n\mathbb{E}[{\varphi }_{ij}{\varphi }_{i{j}^{\prime }}{\varphi }_{i^{\prime }{j}^{\prime }}{\varphi }_{i^{\prime }j}]$$

$$=$$

Full text

1 Introduction

Given an input matrix $A\in\mathbb{R}^{n\times d}$ ($n\geq d$) and a rank parameter $k\leq$ rank($A$), the low rank approximation problem is to find a matrix $\tilde{A}_{k}\in\mathbb{R}^{n\times d}$ of rank-$k$ to approximate the original matrix $A$. Low-rank approximation is a major topic in scientific computing and it is an essential tool in many applications including principal component analysis, data compression, spectral clustering and etc [38, 25, 13].

The best low rank approximation can be obtained through singular value decomposition (SVD) [23], rank revealing QR factorization [8] or two-sided orthogonal factorizations [20]. However, with rapid growth on data dimension, these traditional and deterministic methods are expensive and inefficient [33]. For a matrix $A\in\mathbb{R}^{n\times d}$ $(n\geq d)$, the SVD computation requires $O(nd^{2})$ time which is unfit for large $n$ and $d$. Many randomized algorithms have been developed to target at this large dimension issue. The idea of a randomized algorithm is to trade accuracy with efficiency. A good randomized algorithm tends to have high accuracy with low running time. Random Sampling and random projection are two commonly used randomization techniques which can be easily employed to solve the low rank approximation problem. Random sampling finds a small subset of rows or columns based on a pre-defined probability distribution to form a sketch and it is also known as the column subset selection problem [5]. Possible techniques include subspace sampling [17], norm sampling [15, 26, 16] and adaptive sampling [45, 12]. Random projection combines rows or columns together to produce a small sketch [31, 39]. Possible techniques include subsampled randomized Hadamard transform [40, 46, 6], discrete Fourier transforms [2, 41], sparse subspace embedding (SpEmb) [9, 34, 36], Gaussian projection [10, 24], subspace iterations [25, 24, 35] and the fast Johnson-Lindenstrass transform [1]. Among them, SpEmb is easy to implement, respects the sparsity of the input matrix and achieves similar error bound with the lowest running time.

Apart from randomized algorithms, another approach to produce a sketch of a matrix is through frequent directions (FD) [32, 21]. It is motivated by the idea of frequent items in item frequency approximation problem. This algorithm requires the sketch size to be linear on the rank approximated to reach the relative-error bound for the low rank approximation. It works well in practice, but with high computational cost. An additional point about FD is that besides producing the sketch, it also returns a matrix consisting of its right singular vectors during the process. This is useful for low rank approximation which we will discuss in later sections.

Is there any method that achieves high accuracy comparable to FD and low computational cost similar as SpEmb for the low rank approximation problem? In this paper, we establish a fast frequent directions algorithm, which provides an answer to this question. This new algorithm performs FD by sending more information to each iteration which makes the SVD computation in each iteration more valuable and accurate. In the meanwhile, it can be seen as performing SpEmb on the input matrix to obtain an intermediate sketch, then applying FD on the sketch which reduces the number of iterations required. Our new algorithm preserves the high accuracy of FD with the sketch size required to reach a good error bound linear on the rank approximated; it also reduces the computational cost substantially by cutting down the number of iterations in FD and exploiting the sparsity of the input matrix.

The rest of this paper is organized as follows: We review some related work on low rank approximation, SpEmb and FD in Section 2, they provide the foundation and motivation in developing the new algorithm. Then we introduce our new algorithm with theoretical analysis in Section 3. In Section 4, we evaluate the performance of our new algorithm in comparison with FD and three well-known randomized algorithms including SpEmb on both synthetic and real world datasets. We explore the application on network problems in Section 5. Finally, concluding remarks are given in Section 6.

2 Related Work

Given an input matrix $A\in\mathbb{R}^{n\times d}$, let $A=U_{A}\Sigma_{A}V_{A}^{T}$ be the singular value decomposition of $A$, where the singular values of $A$ are arranged in a non-increasing order in $\Sigma_{A}$. By the Eckart-Young-Mirsky theorem, the best rank-$k$ ($k\leq$ rank($A$)) approximation for both Frobenius norm ($F$-norm) and spectral norm ($2$-norm) is $A_{k}=(U_{A})_{k}(\Sigma_{A})_{k}(V_{A})_{k}^{T}$, where $(\Sigma_{A})_{k}\in\mathbb{R}^{k\times k}$ consists of the top $k$ singular values in its diagonal entries, $(U_{A})_{k}\in\mathbb{R}^{n\times k}$ and $(V_{A})_{k}\in\mathbb{R}^{d\times k}$ consist of the $k$ left and $k$ right singular vectors of $A$ which correspond to the top $k$ singular values, respectively. For randomized algorithms and FD method, they aim to find a low rank approximation $\tilde{A}_{k}$ to minimize $\mathinner{\!\left\lVert A-\tilde{A}_{k}\right\rVert}^{2}_{F}$.

To construct $\tilde{A}_{k}$, the first step is to produce a sketch of the input matrix and it could be formulated as premultiplying a sketching matrix $S\in\mathbb{R}^{\ell\times n}$ to $A$ with sketch size $\ell\ll n$. In particular, for Gaussian projection, $S_{ij}\sim\mathcal{N}(0,1)$ i.i.d. with $i\in[\ell],j\in[n]$; for uniform sampling, $S$ consists of $\ell$ rows uniformly selected from $I_{n}$. $S$ is designed to form $SA$ efficiently and also to preserve the important information in $A$. It can be found in [33, 47, 16] on how such an $S$ is constructed. For the deterministic method FD, the construction of a sketch is discussed in Subsection 2.2.

After forming the sketch $B=SA\in\mathbb{R}^{\ell\times d}$, the next step is to construct the low rank approximation $\tilde{A}_{k}$. Let $\Pi_{B,k}^{F}(A)\in\mathbb{R}^{n\times d}$ denote the best rank-$k$ approximation to $A$ in the row space of $B$, with respect to $F$-norm. It has been proved [5, 6, 47] that

[TABLE]

where $V\in\mathbb{R}^{d\times\ell}$ satisfies $V^{T}V=I_{\ell}$ and its columns span the row space of $B$ and $[AV]_{k}$ denotes the best rank-$k$ approximation of $AV$. Thus, $\tilde{A}_{k}$ could be set as $[AV]_{k}V^{T}$.

2.1 Sparse Subspace Embedding (SpEmb)

Since its introduction by Clarkson and Woodruff [9], SpEmb has received great attention due to its efficient time complexity. It has been extensively used as a randomization technique to solve numerical linear algebraic problems, here we look at its application on low rank approximation.

Define a sketch matrix $S=\phi D\in\mathbb{R}^{\ell\times n}$, such that

• •

$\phi\in\{0,1\}^{\ell\times n}$ is an $\ell\times n$ binary matrix with $\phi_{h(i),i}=1$ and remaining entries [math], where $i\in[n]$ and $h\mathrel{\mathop{\ordinarycolon}}[n]\rightarrow[\ell]$ is a random map so that for each $i\in[n]$, $h(i)=\gamma$ for $\gamma\in[\ell]$ with probability $\frac{1}{\ell}$;

• •

$D$ is an $n\times n$ diagonal matrix with each diagonal entry independently chosen to be $\pm 1$ with equal probability.

$S$ is called a SpEmb matrix, each column of $S$ contains exactly one nonzero entry $-1$ or $+1$, with equal probability. Due to the structure of $S$, the sketch $SA$ can be formed efficiently by exploiting the sparsity of $A$ [9, 47, 36, 34]. SpEmb is described in Algorithm 1.

We can use the QR factorization $B^{T}=VR$ to construct $V$, while $[AV]_{k}$ could be obtained by the SVD on $AV$. The reader may refer to [9, 34, 47] for more details on SpEmb.

2.2 Frequent Directions (FD)

Compared to most randomized algorithms, FD approaches the matrix sketching problem from a different perspective. It is deterministic, space efficient and produces more accurate estimates given a fixed sketch size. It was introduced by Liberty [32], later extended by Ghashami et al. [21] with a comprehensive analysis, and a few variations [11, 22] have been developed since then. The algorithm extends the idea of frequent items for item frequency approximation problem to a general matrix. Given an input matrix $A\in\mathbb{R}^{n\times d}$ and space parameter $\ell$, at first, the algorithm considers the first $2\ell$ rows in $A$ and shrinks its top $\ell$ orthogonal vectors by the same amount to obtain an $\ell\times d$ matrix; then combines them with the next $\ell$ rows in $A$ for the next iteration, and repeat the procedure. It is illustrated in Fig. 1. The detailed algorithm is given in Algorithm 2.

Algorithm 2 is slightly different from the original one in [32, 21] in the sense that it performs one more round of SVD and $B$ update at the end when $B$ contains zero rows. This is to maintain a sketch $B$ of size $\ell\times d$ and to acquire a matrix $V$ whose columns form an orthonormal basis for the row space of $B$.

3 Fast Frequent Directions via Sparse Subspace Embedding (SpFD)

Notice that each iteration of FD requires $O(d\ell^{2})$ running time to compute the SVD and there are $\lceil n/\ell\rceil-1$ iterations. The main cost comes from the large number of iterations when $n$ is extremely large. Here note at the $i$th iteration ($i>1$), each $B_{i}^{0}$ is composed of the approximate $B_{i-1}^{1}$ from the previous step and a new set of data in $A$, from which FD extracts and shrinks the top $\ell$ important right singular vectors to form $B_{i}^{1}$ (refer to Fig. 1).

The goal is to ensure that the top $k$ right singular vectors of $A$ are preserved through every iteration. If the number of iterations is large, there is a higher risk of discarding useful information during the extractions via SVD with shrinking. Since each iteration is considered as an approximation step, to make this approximation more efficient, instead of bringing just $\ell$ rows from $A$ to form $B_{i}^{0}$ in the $i$th iteration, we bring more information of $A$ to $B_{i}^{0}$ in a compact way (squeeze more rows of $A$ in the lower $\ell$ rows of $B_{i}^{0}$). This could make each SVD computation more valuable and reduce the number of iterations which in a way lessens the approximation error. It obviously introduces a new type of approximation error as the information being analyzed in each iteration is no longer of $A$, but rather a sketch of $A$. Therefore, it requires balancing of the sketch error from compressing more information in $\ell$ rows during each iteration and the shrinking error from performing more iterations of SVD computations, i.e., balancing the error incurred in the sketching and the error incurred in FD.

Here SpEmb is used to perform the intermediate sketching. It is efficient and respects the possible sparsity of the input matrix.

Now let’s illustrate the method in detail. Define

[TABLE]

where each $S_{i}$ ($i=1,\ldots,q$) is a SpEmb matrix of dimension $\ell\times n/q$. Here we assume $n$ is a multiple of $q$. If not, we can append zero rows to the end of $A$ until $n$ is. Let $P\in\mathbb{R}^{n\times n}$ denote a random permutation matrix constructed via Fisher-Yates-Knuth shuffle, for each $e_{j}$ (the $j$th column of $I_{n}$) with $j\in[n]$, $\mathbb{P}[Pe_{j}=e_{i}]=1/n$ for $i\in[n]$. $P$ is generated using randperm in MATLAB. Then we apply FD to the intermediate sketch $SPA\in\mathbb{R}^{q\ell\times d}$ to obtain $B\in\mathbb{R}^{\ell\times d}$ and the corresponding $V\in\mathbb{R}^{d\times\ell}$. Fig. 2 illustrates the process.

SpFD is a combination of FD and SpEmb. Each $S_{i}$ maps $n/q$ rows of $PA$ to $\ell$ rows in the intermediate sketch through SpEmb, then perform FD on the intermediate sketch $SPA$. The magnitude of $q$ indicates the weight distributed between FD and SpEmb. When $q=1$, $S$ has only one block, then SpFD becomes SpEmb. When $q=n/\ell$, $S$ is of dimension $n\times n$, SpFD will become FD if $SPA$ contains the same rows as $A\in\mathbb{R}^{n\times d}$. Here note as a SpEmb matrix, each $S_{i}\in\mathbb{R}^{\ell\times\ell}$ may contain many zero rows and so will $S\in\mathbb{R}^{n\times n}$ and the sketch $SPA\in\mathbb{R}^{n\times d}$ which implies SpFD is unlikely to become FD. This will be discussed more while analyzing SpFD’s running time in Section 4. The new method SpFD is described in Algorithm 3.

The random permutation matrix $P$ is constructed to mix the blocks of $A$. Denote $A=[a_{1}~{}a_{2}~{}\ldots~{}a_{n}]^{T}$ where $a_{i}\in\mathbb{R}^{d}$ and let $[i_{1}~{}i_{2}~{}\ldots~{}i_{n}]=$ randperm ($n$), then $PA=[a_{i_{1}}~{}a_{i_{2}}~{}\ldots~{}a_{i_{n}}]^{T}$, it is only required to load $(PA)_{j}=[a_{i_{(j-1)n/q+1}}~{}a_{i_{(j-1)n/q+2}}~{}\ldots~{}a_{i_{jn/q}}]^{T}$ sequentially for each $j=1,\ldots,q$. Therefore, although SpFD requires a permutation step, it is not necessary to load the whole data in memory.

3.1 Main Theorem

Next we are going to prove SpFD’s accuracy. We first look at $SP\in\mathbb{R}^{q\ell\times n}$; show that if $q\ell$ satisfies certain lower bound, $SP$ is a good sketching matrix. $SP$ is in fact a variation of a SpEmb matrix. The difference is that $SP$ does not allow more than $n/q$ rows mapping to the same row in the sketch.

Lemma 1**.**

Given a matrix $U\in\mathbb{R}^{n\times k}$ with orthonormal columns $(k\leq n)$. For any $\varepsilon,\delta\in(0,1)$, let $S\in\mathbb{R}^{q\ell\times n}$ be defined as (2) with

[TABLE]

Let $P$ be an $n\times n$ random permutation matrix (i.e., for each $e_{j}$ with $j\in[n]$, $\mathbb{P}[Pe_{j}=e_{i}]=1/n$ for $i\in[n]$), then with probability at least $1-\delta$,

[TABLE]

Proof.

To show (3), we bound $\mathbb{E}\left[\mathinner{\!\left\lVert U^{T}(SP)^{T}SPU-I\right\rVert}_{F}^{2}\right]$ and then apply Markov’s inequality.

Let $X=U^{T}(SP)^{T}SPU\in\mathbb{R}^{k\times k}$ and $x_{st}$ denote the $(s,t)$th-entry of $X$. Let $\mathcal{L}_{i}$ denote the set of integers from domain $\left[(\lceil\frac{i}{\ell}\rceil-1)\frac{n}{q}+1,\lceil\frac{i}{\ell}\rceil\frac{n}{q}\right]$, the cardinality of $\mathcal{L}_{i}$ is $\frac{n}{q}$. Note that for $i\in\{C\ell+1,\ldots,(C+1)\ell\}$ with an integer $C\geq 0$, the corresponding $\mathcal{L}_{i}$’s are the same. Here note $S\in\mathbb{R}^{q\ell\times n}$ can be written as

[TABLE]

where each $S_{\alpha}=\phi_{\alpha}D_{\alpha}\in\mathbb{R}^{\ell\times n/q}$ ($\alpha\in[q]$) is a SpEmb matrix defined as in Section 2.1. Note that $\phi$ is a block diagonal matrix and $\phi_{ij}=0$ when $j\notin\mathcal{L}_{i}$. For $1\leq s,t\leq k$,

[TABLE]

and its expectation:

[TABLE]

The second equality holds since $\mathbb{E}[d_{j_{1}}d_{j_{2}}]=1$ if $j_{1}=j_{2}$ and [math] otherwise; the third equality is because $\phi_{ij}=0$ for $j\notin\mathcal{L}_{i}$ and $\mathbb{E}[(PU)_{js}(PU)_{jt}]=1/n\sum_{j=1}^{n}U_{js}U_{jt}$. We also need

[TABLE]

Note that

[TABLE]

thus it is required to discuss $\mathbb{E}[\phi_{ij_{1}}\phi_{ij_{2}}\phi_{i^{\prime}j_{1}^{\prime}}\phi_{i^{\prime}j_{2}^{\prime}}]$ and $\mathbb{E}[(PU)_{j_{1}s}(PU)_{j_{2}t}(PU)_{j_{1}^{\prime}s}(PU)_{j_{2}^{\prime}t}]$ for each case:

1. $j_{1}=j_{2}=j_{1}^{\prime}=j_{2}^{\prime}$: Note that each column of $\phi$ has exactly one $1$, so $\phi_{ij}^{2}\phi_{i^{\prime}j}^{2}=0$ if $i\neq i^{\prime}$, then

[TABLE]

2. $j_{1}=j_{2}\neq j_{1}^{\prime}=j_{2}^{\prime}$:

[TABLE]

where the second equality holds since $\phi_{ij}^{2}\phi_{i^{\prime}j^{\prime}}^{2}=0$ if $j\notin\mathcal{L}_{i}$ or $j^{\prime}\notin\mathcal{L}_{i^{\prime}}$; the third equality is by partitioning $i\neq i^{\prime}$ into $\mathcal{L}_{i}\neq\mathcal{L}_{i^{\prime}}$ and $\mathcal{L}_{i}=\mathcal{L}_{i^{\prime}}$. Note there are $q$ different $\mathcal{L}_{i}$’s and $\ell$ choices of $i$ fall in the same $\mathcal{L}_{i}$; for each $i$, as long as $j\in\mathcal{L}_{i}$, then $\mathbb{E}[\phi_{ij}^{2}]=\dfrac{1}{\ell}$, and for $j\neq j^{\prime}$, $\phi_{ij}$ and $\phi_{i^{\prime}j^{\prime}}$ are independent, then the fourth equality follows. Therefore,

[TABLE]

3. $j_{1}=j_{1}^{\prime}\neq j_{2}=j_{2}^{\prime}$:

[TABLE]

and since $\phi$ has only one $1$ in each column, then for $i\neq i^{\prime}$, at least one of $\phi_{ij_{1}},\phi_{i^{\prime}j_{1}}$ is [math], and $\phi_{ij_{1}}\phi_{ij_{2}}\phi_{i^{\prime}j_{1}}\phi_{i^{\prime}j_{2}}=0$, thus

[TABLE]

Therefore,

[TABLE]

4. $j_{1}=j_{2}^{\prime}\neq j_{1}^{\prime}=j_{2}$:

[TABLE]

Therefore,

[TABLE]

Now combining (3.1), (3.1), (3.1) and (7), we obtain

[TABLE]

The first inequality holds since $\dfrac{n/q-1}{(n-1)\ell}\geq 0$; the last inequality follows since $\dfrac{n/q-1}{(n-1)\ell}=\dfrac{1}{q\ell}\dfrac{n-q}{n-1}\leq\dfrac{1}{q\ell}$ as $q\geq 1$. Therefore, we have

[TABLE]

Applying Markov’s inequality, we obtain

[TABLE]

By taking $q\ell\geq\dfrac{k^{2}+k}{\varepsilon^{2}\delta}$, we conclude that

[TABLE]

∎

Remark 1**.**

Lemma 1 is to find a condition for a SpFD matrix to qualify as a $(1\pm\varepsilon)$ $\ell_{2}$-subspace embedding for $U$. By setting $q=1$, $SPA\in\mathbb{R}^{\ell\times d}$ and no FD is required, SpFD becomes SpEmb. Then the bound derived here for SpFD agrees with that for a general SpEmb matrix (Theorem 2.5 in [47]).

Next, we move on to discuss the FD part in our algorithm and to determine the bound on the final sketch size $\ell$. After obtaining $SPA$, the next step is applying FD on $SPA$ to generate the sketch $B$, the orthonormal basis $V$ and eventually the approximation $\tilde{A}_{k}$.

As mentioned in Section 2.2, the FD process in Algorithm 2 and Algorithm 3 performs one last round of SVD and $B$ update at the end. Nevertheless, the same properties can be easily derived following the same procedure in [32, 21]. Then for the FD process embedded in Algorithm 3 of SpFD, we have the following properties: let

[TABLE]

Property 1: For any vector $x\in\mathbb{R}^{d}$, $\mathinner{\!\left\lVert SPAx\right\rVert}_{2}^{2}-\mathinner{\!\left\lVert Bx\right\rVert}_{2}^{2}\geq 0$.

Property 2: For any unit vector $x\in\mathbb{R}^{d}$, $\mathinner{\!\left\lVert SPAx\right\rVert}_{2}^{2}-\mathinner{\!\left\lVert Bx\right\rVert}_{2}^{2}$ $\leq\Delta$.

Property 3: $\mathinner{\!\left\lVert SPA\right\rVert}_{F}^{2}-\mathinner{\!\left\lVert B\right\rVert}_{F}^{2}\geq\ell\Delta$.

Besides Lemma 1 and the properties listed above, we need one more lemma to derive the main theorem.

Lemma 2**.**

For any matrix $C\in\mathbb{R}^{n\times d}$, let $S\in\mathbb{R}^{q\ell\times n}$ be defined as (2) where $q$ and $\ell$ are positive integers and $P$ be an $n\times n$ random permutation matrix defined earlier, then $\mathbb{E}\left[\mathinner{\!\left\lVert SPC\right\rVert}_{F}^{2}\right]=\mathinner{\!\left\lVert C\right\rVert}_{F}^{2}$. Furthermore, by applying Markov’s inequality, with probability at least $1-\delta$ ($0<\delta<1$),

[TABLE]

This lemma states that without any requirement on the size of the sketching matrix, the expected value of the $F$-norm of the sketched matrix is the $F$-norm of the original matrix. The proof is provided in Section 1 of the supplemental material. Now we are ready to derive our main theorem.

Theorem 1**.**

Given a matrix $A\in\mathbb{R}^{n\times d}$ and let $A_{k}$ denote the best rank-$k$ approximation of $A$. For any $\varepsilon,\delta\in(0,1)$, let $S\in\mathbb{R}^{q\ell\times n}$ be defined as (2) with

[TABLE]

and $P$ be an $n\times n$ random permutation matrix (i.e., for each $e_{j}$ with $j\in[n]$, $\mathbb{P}[Pe_{j}=e_{i}]=1/n$ for $i\in[n]$). Let $\tilde{A}_{k}$ and $B\in\mathbb{R}^{\ell\times d}$ be the output of Algorithm 3, then with probability at least $1-\delta$,

[TABLE]

Proof.

Let

[TABLE]

be the singular value decompositions of $A\in\mathbb{R}^{n\times d}$ and $B\in\mathbb{R}^{\ell\times d}$, $\rho$ = rank($A$) and $\beta$ = rank($B$) and the singular values are listed in a non-increasing order in $\Lambda$ and $\Sigma$. Thus,

[TABLE]

The inequality holds by taking $X=AV_{k}\Sigma_{k}^{-1}U_{k}^{T}$; the last equality holds since $(A-AV_{k}V_{k}^{T})V_{k}V_{k}^{T}A^{T}=0$ which implies $\mathinner{\!\left\lVert A-AV_{k}V_{k}^{T}\right\rVert}_{F}^{2}=\mathinner{\!\left\lVert A\right\rVert}_{F}^{2}-\mathinner{\!\left\lVert AV_{k}\right\rVert}_{F}^{2}$ (by Lemma 2.1 in [5]).

Let $U_{AV_{k}}\Sigma_{AV_{k}}V_{AV_{k}}$ be the SVD of $AV_{k}$ where $U_{AV_{k}}\in\mathbb{R}^{n\times k}$ is a matrix whose columns composed of the left singular vectors of $AV_{k}$; assume

[TABLE]

for some $\varepsilon_{0},\delta_{0}\in(0,1)$, then by Lemma 1, with probability at least $1-\delta_{0}$, all singular values of $SPU_{AV_{k}}$ are in $[\sqrt{1-\varepsilon_{0}},\sqrt{1+\varepsilon_{0}}]$. Then with probability at least $1-\delta_{0}$,

[TABLE]

which implies

[TABLE]

Conditioned on (10), then (8) becomes

[TABLE]

The second inequality is based on $\mathinner{\!\left\lVert A_{k}\right\rVert}_{F}^{2}=\mathinner{\!\left\lVert AQ_{k}\right\rVert}_{F}^{2}$ and (Property 1). The third inequality follows since $V_{k}$ contains the top $k$ right singular vectors of $B$ and thus $\mathinner{\!\left\lVert BV_{k}\right\rVert}_{F}^{2}\geq\mathinner{\!\left\lVert BQ_{k}\right\rVert}_{F}^{2}$.

To bound $\mathinner{\!\left\lVert AQ_{k}\right\rVert}_{F}^{2}-\mathinner{\!\left\lVert BQ_{k}\right\rVert}_{F}^{2}$, we use a bridge term $\mathinner{\!\left\lVert SPAQ_{k}\right\rVert}_{F}^{2}$ and bound the singular values of $SPW_{k}$ using Lemma 1 again. Still assume (9), thus with probability at least $1-\delta_{0}$, for all $1\leq i\leq k$,

[TABLE]

This implies that $SPW_{k}\in\mathbb{R}^{q\ell\times k}$ is of full column rank and $(SPW_{k})^{\dagger}SPW_{k}=I_{k}$, then

[TABLE]

Together they imply with probability at least $1-\delta_{0}$,

[TABLE]

Conditioned on (10) and (12), then (11) continues

[TABLE]

The last inequality follows from (Property 2). To bound $\Delta$, we consider

[TABLE]

where the last inequality follows from (Property 2). On the other hand, recall that $A=W_{k}\Lambda_{k}Q_{k}^{T}+W_{\rho-k}\Lambda_{\rho-k}Q_{\rho-k}^{T}$, then

[TABLE]

The second equality follows from

[TABLE]

which implies

[TABLE]

again by Lemma 2.1 in [5]; and $\mathinner{\!\left\lVert SPW_{k}\Lambda_{k}Q_{k}^{T}\right\rVert}_{F}^{2}=\mathinner{\!\left\lVert SPAQ_{k}Q_{k}^{T}\right\rVert}_{F}^{2}=\mathinner{\!\left\lVert SPAQ_{k}\right\rVert}_{F}^{2}$. The first inequality is based on (Property 3). And the last inequality

[TABLE]

holds with probability at least $1-\delta_{1}$ where $\delta_{1}\in(0,1)$, followed from Lemma 2 by setting $\delta=\delta_{1}$ and $C=W_{\rho-k}\Lambda_{\rho-k}$. Conditioned on (16), we combine (14) and (15), then

[TABLE]

Now we substitute (17) to (13) and apply union bound on (10), (12) and (16), thus we conclude that with probability at least $1-2\delta_{0}-\delta_{1}$,

[TABLE]

By (1),

[TABLE]

Taking $\varepsilon_{0}=\dfrac{1}{2}\varepsilon$, $\delta_{0}=\delta_{1}=\dfrac{1}{3}\delta$ and $\ell=\left\lceil k+\dfrac{3k}{\varepsilon\delta}\right\rceil$, then

[TABLE]

Hence by taking

[TABLE]

we conclude that with probability at least $1-\delta$,

[TABLE]

∎

Next, we analyze the computational complexity of Algorithm 3 for SpFD, comparing with Algorithm 1 for SpEmb and Algorithm 2 for FD. The analysis is in most detail as it will be used for explanation of the numerical results in Section 4 and 5. With the same sketch size $\ell$, the number of flops for each algorithm are summarized in TABLE I.

Here note that QR factorization on a $d\times\ell$ $(d\geq\ell)$ matrix requires $4d\ell^{2}-\dfrac{4}{3}\ell^{3}$ flops; SVD on an $n\times\ell$ $(n\geq\ell)$ matrix requires $6n\ell^{2}+20\ell^{3}$ flops and the flop counts follow from [23]. For SpEmb, as the sketching matrix is generated on a probability distribution, we count the flops in worst case of constructing $SA$ 111Note $SA$ is formed by combining the rows of $A$ directly and no explicit $S$ is created.. Then we can conclude:

• •

To achieve the best error bounds, SpFD requires $\ell$ to be linear on $k$ which is the same as for FD [32, 21], and lower compared to SpEmb [47, 9].

• •

Since $n$ is large, $n\geq d\geq\ell\geq k$, $n\gg\ell$ and $q\ell$ is only required to be $O(k^{2})$ by Theorem 1, the number of iterations have been significantly reduced from $n/\ell$ to $q=q\ell/\ell$. Thus SpFD runs much faster than FD for the same sketch size. Nevertheless, when $q$ approaches $n/\ell$, $S\in\mathbb{R}^{q\ell\times n}$ becomes more square-like and contains many zero rows with a high probability, so will $SPA$. Hence performing FD on $SPA$ would still be much more efficient than on $A$. A detailed discussion on the number of zero rows in $S$ is provided in Section 2 of the supplemental material.

4 Experimental Results

In this section, we evaluate the performance of the proposed SpFD algorithm by comparing it with four other commonly used algorithms on both synthetic and real datasets. All algorithms are implemented in MATLAB R2016a, and all experiments are conducted on a PC Dell OptiPlex 7040. The competing algorithms are:

1. Norm Sampling (NormSamp) [15, 26, 16] : It is a sampling method. The sketch is composed of $\ell$ rows chosen i.i.d. from the $n$ rows of $A$ and scaled properly. Each row in the sketch takes $A_{(i)}/\sqrt{\ell p_{i}}$ with probability $p_{i}=\mathinner{\!\left\lVert A_{(i)}\right\rVert}_{2}^{2}/\mathinner{\!\left\lVert A\right\rVert}_{F}^{2}$. Norm sampling requires $O(nd)$ time to form the probability distribution $p_{i}$’s and construct the sketch [16].

2. Discrete Cosine Transform (DCT) [2, 41]: It is a random projection method, a type of fast Fourier transform which has more practical advantages compared to SRHT [2, 48]. The sketch is formed as $\sqrt{n/\ell}RFDA\in\mathbb{R}^{\ell\times d}$ with

• •

$R\in\mathbb{R}^{\ell\times n}$ is a subset of $\ell$ rows from $I_{n}$, where the rows are chosen uniformly and without replacement;

• •

$F\in\mathbb{R}^{n\times n}$ is a normalized discrete cosine transform matrix;

• •

$D\in\mathbb{R}^{n\times n}$ is a random diagonal matrix whose entries are independently chosen to be $\pm 1$ with equal probability.

3. SpEmb: Follow Algorithm 1.

4. FD: Follow Algorithm 2.

5. SpFD: Follow Algorithm 3. Here we employ three different $q$ for comparison, refer to Fig. 2 for an easy understanding (in later discussions, the following three variations together are referred to as SpFDq):

• •

$q=5$ (SpFD5): For $\alpha=1,\ldots,5$, each $S_{\alpha}\in\mathbb{R}^{\ell\times n/q}$ combines $n/5$ rows in $PA$ to $\ell$ rows in each iteration for SVD computation. In other words, $S\in\mathbb{R}^{q\ell\times n}$ consists of $5$ diagonal blocks, sketching $PA\in\mathbb{R}^{n\times d}$ into $5$ $\ell\times d$ blocks, then perform FD.

• •

$q=10$ (SpFD10): $S\in\mathbb{R}^{q\ell\times n}$ consists of $10$ diagonal blocks, sketching $PA\in\mathbb{R}^{n\times d}$ into $10$ $\ell\times d$ blocks, then perform FD.

• •

$q=50$ (SpFD50): $S\in\mathbb{R}^{q\ell\times n}$ consists of $50$ diagonal blocks, sketching $PA\in\mathbb{R}^{n\times d}$ into $50$ $\ell\times d$ blocks, then perform FD.

For randomized algorithms, after forming the sketch matrix $B$, we follow the same procedure as SpEmb, use QR factorization to construct $V$ and SVD to construct $[AV]_{k}$, and then obtain $\tilde{A}_{k}=[AV]_{k}V^{T}$.

The randomized algorithms with dense embedding matrices used in practice use sketch sizes of $\ell=k+p$, where $p$ is a small nonnegative integer [25]. It has also been demonstrated in [41] that $\ell=4k$ is empirically sufficient for sparse embeddings. These results have been taken into considerations while setting the sketch sizes in our numerical experiments below.

Synthetic Data

The construction of synthetic datasets follows from [32, 21] which were used to test FD. Each row of the generated input matrix $A\in\mathbb{R}^{n\times d}$ consists of $k$ dimensional signal and $d$ dimensional noise ($k\ll d$). To put it in an accurate form, $A=S^{c}DU+N/\zeta$. $S^{c}\in\mathbb{R}^{n\times k}$ is the signal coefficient matrix and $S^{c}_{ij}\sim\mathcal{N}(0,1)$ i.i.d. for $i\in[n]$ and $j\in[k]$. $D\in\mathbb{R}^{k\times k}$ is a diagonal matrix with $D_{ii}=1-(i-1)/m$ which gives linearly diminishing signal singular values. $U\in\mathbb{R}^{k\times d}$ is the signal row space matrix which contains a random $k$ dimensional subspace in $\mathbb{R}^{d}$ and $UU^{T}=I_{k}$. The matrix $S^{c}DU$ is of exactly rank-$k$ and constitutes the signal we wish to recover. $N\in\mathbb{R}^{n\times d}$ contributes additive Gaussian noise $N_{ij}\sim\mathcal{N}(0,1)$. According to [42] and as mentioned in [32, 21], to ensure the signal is recoverable, it is required that $c_{1}\leq\zeta\leq c_{2}\sqrt{d/k}$ for constants $c_{1}$ and $c_{2}$ close to $1$. In the experiments, we consider $n=10000$ and $d=1000$, $k\in\{10,20,50\}$ ($k=10$ as default value) and $\zeta\in\{5,10,15,20\}$ ($\zeta=10$ as default value). We set the sketch size $\ell\in\{k\mathrel{\mathop{\ordinarycolon}}10\mathrel{\mathop{\ordinarycolon}}200\}$. For each dataset setting, we generate $3$ matrices and for each of these matrices, we run each method $5$ times and the results are taken as the median of all $15$ repetitions for each dataset setting.

Real Data

We consider the six real world datasets in TABLE II. For each dataset, we run $10$ times for each method and the results presented are the median of the $10$ repetitions.

Results

To measure the accuracy, we set the best rank-$k$ approximation $A_{k}$ obtained through SVD of $A$ as our benchmark. For each method, we examine the following three measures:

• •

$F$-norm Error:  $\mathinner{\!\left\lVert A-\tilde{A}_{k}\right\rVert}_{F}/\mathinner{\!\left\lVert A-A_{k}\right\rVert}_{F}$;

• •

$2$-norm Error:  $\mathinner{\!\left\lVert A-\tilde{A}_{k}\right\rVert}_{2}/\mathinner{\!\left\lVert A-A_{k}\right\rVert}_{2}$;

• •

Running time to construct $\tilde{A}_{k}$ in seconds.

For synthetic datasets of different signal dimension $k$ with $k=10,20,50$ while $\zeta=10$, and different noise ratio $\zeta$ with $\zeta=5,15,20$ while $k=10$, the results are displayed in Fig. 3 (a) - (f). The results for real datasets are displayed in Fig. 3 (g) - (h) and Fig. 4 (i) - (l).

Regarding on the accuracy of the competing algorithms for both $F$-norm and $2$-norm, the relative errors for all algorithms decline as the sketch size $\ell$ increases. Among them, NormSamp, DCT and SpEmb achieve similar relative errors; our algorithm SpFD5, SpFD10 and SpFD50 obtain much more accurate results than the other randomized algorithms; and in general FD is the most accurate. The accuracy level of SpFD5, SpFD10 and SpFD50 are close to that of FD.

Among SpFD5, SpFD10 and SpFD50, the ordering of their accuracy level cannot be determined precisely. SpFD5 contains larger error from the SpEmb it incorporates due to the smaller intermediate sketch size but less error from FD as the number of SVD computations is reduced. SpFD50 on the other hand, with larger error from FD and less error from the intermediate sketching. The accuracy level of each SpFDq depends on which error dominates. Nevertheless, since FD is considered more accurate than SpEmb, we can find from our results that the accuracy level of SpFDq $(q=5,10,50$) improves as $q$ increases for most datasets.

For some datasets, SpFDq achieve higher accuracy level than FD. This suggests it is possible that integrating certain level of SpEmb into FD generates a better result than the original FD.

From the running time perspective, all algorithms’ running time shows a trend of inclining as the sketch size $\ell$ increases. Among the algorithms, FD is the slowest; then followed by DCT; NormSamp, SpEmb and SpFDq are gradually increasing with the sketch size and they are far more efficient compared to FD.

A close-up view of the running time for SpFDq ($q=5,10,50$), NormSamp and SpEmb on all datasets is displayed in Fig. 1 in Section 3 of the supplemental material. NormSamp computes the probability distribution by calculating the norm square of each row that respects the sparsity of $A$. It is generally faster than projection methods like SpEmb since it does not require any row combinations to form the sketch. SpFDq ($q=5,10,50$) has similar running time as NormSamp and faster than SpEmb; in addition, the running time of SpFDq ($q=5,10,50$) reduces as $q$ increases for most datasets although the difference among them is not obvious. But flop counts in TABLE I show the opposite. To understand this, possible reasons are:

• •

Demonstrated by the two tables in Section 4 of the supplemental material, the running time for computing $[AV]_{k}$ and $\tilde{A}_{k}=[AV]_{k}V^{T}$ are similar for SpEmb and SpFDq as predicted, however the running time for constructing $V$ for SpEmb, SpFD5, SpFD10 and SpFD50 declines according to this order. Note that this running time consists of two parts – forming $SA$ or $SPA$ and performing $1$ round of QR or $q-1$ rounds of SVD. One possible reason is the inclining number of zero rows in $S$, for which the reader may refer to Section 2 of the supplemental material. This will cut down the cost of the $q-1$ rounds of SVD for SpFDq with large $q$. The other possible reason is that as a SpEmb (-like) matrix ($S\in\mathbb{R}^{\ell\times n}$ for SpEmb, $S\in\mathbb{R}^{q\ell\times n}$ for SpFDq), the smaller the sketch size, the larger number of row combinations is required while forming $SA$ or $SPA$.

• •

Note that the earlier versions of MATLAB counted the flops. With the incorporation of LAPACK in MATLAB 6 (and newer versions), this is no longer practical since many optimization techniques are adopted in LAPACK, and furthermore the memory latency to fetch anything not in cache is much greater than the cost of a flop, so the flop counts in TABLE I for each algorithm and the running time in Fig. 1 of Section 3 in the supplemental material may not be consistent.

Based on the observation and discussions above, as the efficiency is not compromised much for large $q$, then from the accuracy perspective, we may choose a relatively large $q$ for SpFD in practice. In conclusion, our new algorithm SpFD can be seen as an integration of FD and SpEmb in favor of FD’s accuracy and SpEmb’s efficiency, as a result, it is the best among the algorithms we compared in terms of accuracy and efficiency.

5 Applications in Network Analysis

Large-scale networks arise in many applications, one basic question often asked in network analysis is to identify the most ‘important’ nodes in a network [7, 19]. In a directed network, two types of nodes are of interest: hubs and authorities [27, 3]. Hubs are nodes which point to many nodes that are considered important while authorities are these important nodes. Thus, important hubs are those which point to many important authorities and important authorities are those pointed to by many important hubs. To measure the importance, each node has a ‘hub score’ and an ‘authority score’. It is necessary to have algorithms computing the hub and authority scores associated with each node and thus determine the important nodes.

Given a directed graph $G$, $G=(V,E)$ is formed by a set of vertices $V$ and a set of edges $E$ which is formed by ordered pairs of vertices. $(i,j)\in E$ implies there is a link from node $i$ to $j$. The adjacency matrix of $G$ is a matrix $A\in\mathbb{R}^{\mathinner{\!\left\lvert V\right\rvert}\times\mathinner{\!\left\lvert V\right\rvert}}$ defined in the following way:

[TABLE]

Method I (HITS) [27, 3]: A well-known algorithm in ranking the hubs and authorities is Hypertext Induced Topics Search (HITS) algorithm. HITS is in fact an iterative power method to compute the dominant eigenvectors for $A^{T}A$ and $AA^{T}$ which correspond to the authority and hub scores respectively. It depends on the choice of the initial vectors and only the information obtained from the dominant eigenvectors of $A^{T}A$ or $AA^{T}$ is used. In our experiments, the starting vectors we use are with all entries i.i.d. from $\mathcal{N}(0,1)$ and the stopping criterion is $\mathinner{\!\left\lVert x^{(k)}-x^{(k-1)}\right\rVert}_{2}\leq 10^{-3}$ where $x^{(k)}$ represents the resulting vector at the $k$th iteration.

Method II ($e^{\mathcal{A}}$) [18, 3]: The hub and authority scores of node $i$ are given by $(e^{\mathcal{A}})_{ii}$ and $(e^{\mathcal{A}})_{n+i,n+i}$ respectively for $i\in[n]$, where $\mathcal{A}=\begin{bmatrix}0&A\\ A^{T}&0\end{bmatrix}\in\mathbb{R}^{2n\times 2n}$. Compared to HITS, this method takes spectral information from all eigenvectors of $A^{T}A$ and $AA^{T}$ (with declining weights) into consideration which may lead to improved results.

Method III: Low rank approximation techniques could be applied to compute the hub and authority scores. For a randomized algorithm with sketching matrix $S$, the orthogonal basis $V\in\mathbb{R}^{n\times\ell}$ for the row space of the sketch $SA\in\mathbb{R}^{\ell\times n}$ is obtained similar as in Algorithm 1. The orthogonal bases for FD and SpFD are obtained following Algorithm 2 and Algorithm 3 respectively. The sketch size $\ell$ is chosen based on $k$, the number of top hubs and authorities required in a question. Then following Algorithm 4, we could obtain the approximated SVD of $A$, $\tilde{U}\tilde{\Sigma}\tilde{V}^{T}$ where $\tilde{U},\tilde{V}\in\mathbb{R}^{n\times\ell}$ and $\tilde{\Sigma}\in\mathbb{R}^{\ell\times\ell}$.

Note that since $e^{\mathcal{A}}=~{}\begin{bmatrix}U\cosh(\Sigma)U^{T}&U\sinh(V^{T})\\ V\sinh(\Sigma)U^{T}&V\cosh(\Sigma)V^{T}\end{bmatrix}$ where $A=U\Sigma V^{T}$ is the SVD of $A$. We could use $\left(\tilde{U}\cosh(\tilde{\Sigma})\tilde{U}^{T}\right)_{ii}$ and $\left(\tilde{V}\cosh(\tilde{\Sigma})\tilde{V}^{T}\right)_{ii}$ to approximate the hub and authority scores of the $i$th node, respectively. In the experiments, we use $\ell=k+p$ where $p$ is set to be $5$.

We again compare the five algorithms NormSamp, DCT, SpEmb, FD and SpFDq ($q=5,10,50$) against the HITS and $e^{\mathcal{A}}$ methods. All methods are applied to the following two real world networks.

• •

Computational Complexity: The dataset contains $884$ nodes and $1616$ directed edges.

• •

Death Penalty: The dataset contains $1850$ nodes and $7363$ directed edges.

The rank of the decomposition used in the standard algorithms for the original hubs and authorities is $k=1$ [27]. In our experiments, we consider a more general case, that is, for each network we obtain the top $10$ hubs and authorities ($k=10$). The results are summarized in TABLE III.

For network Computational Complexity, HITS and $e^{\mathcal{A}}$ produce different rankings, they share 4 out of the top 10 hubs and 4 out of the top 10 authorities. This discrepency is expected with reasons as discussed earlier. As the low rank approximation methods follow the computation of $e^{\mathcal{A}}$, the results are compared against $e^{\mathcal{A}}$. For the hub ranking, SpFD5, SpFD10, SpFD50 and FD generate 9, 9, 9 and 9 out of 10 top hubs, respectively. In terms of authority ranking, SpFD5, SpFD10, SpFD50 and FD have 8, 8, 8, 9 out of 10 top authorities respectively projected by $e^{\mathcal{A}}$.

For network Death Penalty, HITS and $e^{\mathcal{A}}$ generate the same top 10 hubs and 10 authorities with the same orderings. FD, SpFD5, SpFD10 and SpFD50 all have the same top 10 hubs as HITS and $e^{\mathcal{A}}$ with slight variations on the ordering. Similarly for authority ranking, FD and SpFD50 identify all top 10 authorities, followed by SpFD10 and SpFD5 with 9 and 8 out of 10 top authorities identified respectively.

In comparison, SpFDq and FD are more accurate than NormSamp, DCT and SpEmb and this outcome aligns with the results obtained in Section 4. In terms of running time, SpFDq are more efficient than DCT and FD and similar as NormSamp and SpEmb which is also in line with the results in Section 4. From this exercise, it is obvious that the low rank approximation could be applied to the network problem of identifying hubs and authorities. According to the accuracy level and efficiency, we have seen that SpFDq ($q=5,10,50$) perform the best among the competing algorithms.

6 Concluding Remarks

In this paper, we have developed a fast frequent directions algorithm SpFD for the low rank approximation problem. It makes use of the natural block structure of FD and incorporates the idea of SpEmb in it and because of that, SpFD is able to accelerate the relatively inefficient algorithm FD by a great amount. On the other hand, SpFD obtains a good low rank approximation with sketch size linear on the rank approximated which implies it is more accurate than SpEmb. These two points are well supported by our experimental results and the application on network problems.

Many interesting issues are worthy of further study. These issues include:

• •

Based on the performance of SpFD in our experiments, the proved error bound in Theorem 1 does not seem to be tight. The techniques we used to derive SpFD’s error bound still largely depend on the properties derived for FD, it is important to seek new ways to derive a better error bound by further exploiting the well-defined structure of SpFD, this is still under investigation.

• •

An alternative to using FD would be to form the strong rank-revealing QR decomposition for a sparse sketch of $A$ (which could be terminated once the diagonal entries of the triangular/trapezoidal factor are small enough). Another thread would be to replace FD with other randomized algorithms such as in [25]. An advantage of FD and our SpFD is that the matrix can be processed sequentially. To improve the performance, other randomized techniques [33] or sparse sketching methods [41] could be used to form the sketch. They would have similar theoretical guarantees proved in this paper. This is a subject for our future work.

• •

All algorithms for SVD are necessarily iterative (the number of iterations required may be small and the flop counts depend on the machine precision). The data needs to be clean in order for SpFD and FD to work well for big datasets. In practice, noisy big datasets require extra power iterations to attain acceptable accuracy [25]. This has been demonstrated, for example, by a method named SFD in [22]. SFD in [22] is a variation of FD which introduces a lower bound on the number of nonzero entries being processed in each iteration and it adopts power iteration to boost the accuracy of the shrinking. Numerical results in [22] has shown that SFD is very efficient for highly sparse matrices, so we only compare the accuracy of SFD, FD and SpFDq for three sparse datasets w8a, Birds and rcv1 in Fig. 5 (the number of power iterations in SFD is set to be $1$). Clearly SFD achieves a higher accuracy level for both $2$-norm and $F$-norm compared to FD and SpFDq.

It is critical to measure accuracy for big data using the spectral norm, as discussed in the appendix of [29]. To the best of our knowledge, there is still a lack of good error bounds measured by spectral norm for FD and many randomized methods. According to our numerical results, SpFD and FD work well for both Frobenius norm and spectral norm approximation. We are now working on incorporating power iterations into randomized sketch to derive good spectral norm error bounds.

Acknowledgments

The authors would like to thank the associate editor and anonymous reviewers for their valuable comments and suggestions on earlier versions of this paper.