Fast Updating Truncated SVD for Representation Learning with Sparse Matrices

Projection Viewpoint. Recent perspectives (Vecharynski & Saad, 2014; Kalantzis et al., 2021) frame the prevailing methodologies for updating the truncated SVD as instances of the Rayleigh-Ritz projection process, which can be characterized by following steps.

1. 1.

   Augment the singular vector ${\bm{U}}_{k}$ and ${\bm{V}}_{k}$ by adding extra columns (or rows), resulting in $\widehat{{\bm{U}}}$ and $\widehat{{\bm{V}}}$, respectively. This augmentation is performed so that $range(\widehat{{\bm{U}}})$ and $range(\widehat{{\bm{V}}})$ can effectively approximate and capture the updated singular vectors’ column space.

2. 2.

   Compute ${\bm{F}}_{k},\mathbf{\Theta}_{k},{\bm{G}}_{k}=\text{SVD}_{k}(\widehat{{\bm{U}}}^{\top}{\bm{A}}\widehat{{\bm{V}}})$.

3. 3.

   Approximate the updated truncated SVD by $\widehat{{\bm{U}}}{\bm{F}}_{k},\mathbf{\Theta}_{k},\widehat{{\bm{V}}}{\bm{G}}_{k}$.

Zha-Simon’s Scheme. For $\widetilde{{\bm{A}}}_{k}={\bm{U}}_{k}\mathbf{\Sigma}_{k}{\bm{V}}_{k}^{\top}$, let $({\bm{I}}-{\bm{U}}_{k}{\bm{U}}_{k}^{\top}){{{\bm{E}}_{c}}}={\bm{Q}}{\bm{R}}$, where ${\bm{Q}}$’s columns are orthonormal and ${\bm{R}}$ is upper trapezoidal. Zha & Simon (1999) method updates the truncate SVD after appending new columns ${\bm{E}}_{c}$ $\in\mathbb{R}^{m\times s}$ to ${\bm{A}}\in\mathbb{R}^{m\times n}$ by

where ${\bm{F}}_{k},\mathbf{\Theta}_{k},{\bm{G}}_{k}=\text{SVD}_{k}(\begin{bmatrix}\mathbf{\Sigma}_{k}&{\bm{U}}_{k}^{\top}{{{\bm{E}}_{c}}}\\ &{\bm{R}}\end{bmatrix})$. The updated approximate left singular vectors, singular values and right singular vectors are $\begin{bmatrix}{\bm{U}}_{k}~{}~{}{\bm{Q}}\end{bmatrix}{\bm{F}}_{k},\mathbf{\Theta}_{k},\begin{bmatrix}{\bm{V}}_{k}&\\ &{\bm{I}}\end{bmatrix}{\bm{G}}_{k}$, respectively.

When it comes to weight update, let the QR-decomposition of the augmented matrices be $({\bm{I}}-{\bm{U}}_{k}{\bm{U}}_{k}^{\top}){\bm{D}}={\bm{Q}}_{\bm{D}}{\bm{R}}_{\bm{D}}$ and $({\bm{I}}-{\bm{V}}_{k}{\bm{V}}_{k}^{\top}){{{\bm{E}}_{m}}}={\bm{Q}}_{\bm{E}}{\bm{R}}_{\bm{E}}$, then the update procedure is

where ${\bm{F}}_{k},\mathbf{\Theta}_{k},{\bm{G}}_{k}=\text{SVD}_{k}(\begin{bmatrix}\mathbf{\Sigma}_{k}&\mathbf{0}\\ \mathbf{0}&\mathbf{0}\\ \end{bmatrix}+\begin{bmatrix}{\bm{U}}_{k}^{\top}{\bm{D}}\\ {\bm{R}}_{\bm{D}}\end{bmatrix}\begin{bmatrix}{\bm{V}}_{k}^{\top}{{{\bm{E}}_{m}}}\\ {\bm{R}}_{\bm{E}}\end{bmatrix}^{\top})$. The updated approximate truncated SVD is $\begin{bmatrix}{\bm{U}}_{k}~{}~{}{\bm{Q}}_{\bm{D}}\end{bmatrix}{\bm{F}}_{k},\mathbf{\Theta}_{k},\begin{bmatrix}{\bm{V}}_{k}~{}~{}{\bm{Q}}_{\bm{E}}\end{bmatrix}{\bm{G}}_{k}$.

Orthogonalization of Augmented Matrix. The above QR decomposition of the augmented matrix and the consequent compact SVD is of high time complexity and occupies the majority of the time in the whole process when the granularity of the update is large (i.e., $s$ is large). To reduce the time complexity, a series of subsequent methods have been developed based on a faster approximation of the orthogonal basis of the augmented matrix. Vecharynski & Saad (2014) used Lanczos vectors conducted by a Golub-Kahan-Lanczos (GKL) bidiagonalization procedure to approximate the augmented matrix. Yamazaki et al. (2015) and Ubaru & Saad (2019) replaced the above GKL procedure with Randomized Power Iteration (RPI) and graph coarsening, respectively. Kalantzis et al. (2021) suggested determining only the left singular projection subspaces with the augmented matrix set as the identity matrix, and obtaining the right singular vectors by projection.

In this section, we introduce an inner product space that is isomorphic to the column space of the augmented matrix, where each element is a tuple of a sparse vector and a low-dimensional vector. The orthogonalization process in this space can exploit sparsity and low dimension, and the resulting orthonormal basis is bijective to the orthonormal basis of the column space of the augmented matrix.

Previous methods perform QR decomposition of the augmented matrix with ${\bm{Q}}{\bm{R}}=({\bm{I}}-{\bm{U}}_{k}{\bm{U}}_{k}^{\top}){\bm{B}}$, to obtain the updated orthogonal matrix. The matrix ${\bm{Q}}$ derived from the aforementioned procedure is not only orthonormal, but its column space is also orthogonal to the column space of ${\bm{U}}_{k}$, implying that matrix $\begin{bmatrix}{\bm{U}}_{k}~{}~{}{\bm{Q}}\end{bmatrix}$ is orthonormal.

Let ${\bm{Z}}=({\bm{I}}-{\bm{U}}_{k}{\bm{U}}_{k}^{\top}){\bm{B}}={\bm{B}}-{\bm{U}}_{k}({\bm{U}}_{k}^{\top}{\bm{B}})={\bm{B}}-{\bm{U}}_{k}{\bm{C}}$ be the augmented matrix, then each column of ${\bm{Z}}$ can be written as a sparse $m$-dimensional matrix of column vectors minus a linear combination of the column vectors of the $m$-by-$k$ orthonormal matrix ${\bm{U}}_{k}$ i.e., ${\bm{z}}_{i}={\bm{b}}_{i}-{\bm{U}}_{k}{\bm{C}}[i]$. We refer to the form of a sparse vector minus a linear combination of orthonormal vectors as SV-LCOV, and its definition is as follows:

With the scalar multiplication, addition and inner product operations of SV-LCOV, we can further study the inner product space of SV-LCOV.

The isometry of an inner product space here is a transformation that preserves the inner product, thereby preserving the angles and lengths in the space. From Lemma 3, we can see that in SV-LCOV, since the dot product remains unchanged, the orthogonalization process of an augmented matrix can be transformed into an orthogonalization process in the inner product space.

Approximating the augmented matrix with SV-LCOV. The Gram-Schdmit process is less efficient when $s$ is large. To this end, Vecharynski & Saad (2014) and Yamazaki et al. (2015) approximated the orthogonal basis of the augmented matrix with Golub-Kahan-Lanczos bidiagonalization(GKL) procedure (Golub & Kahan, 1965) and Randomized Power Iteration (RPI) process. We find they consists of three operations in Lemma 2 and can be transformed into SV-LCOV for efficiency improvement. Limited by space, we elaborate the proposed algorithm of appoximating the augmented matrix with SV-LCOV in Appendix D.1 and Appendix D.2, respectively.

Low-dimensional matrix mutliplication and sparse matrix addition. Suppose an orthonormal basis $({\bm{b}}_{1},{\bm{c}}_{1})_{{\bm{U}}_{k}},...,({\bm{b}}_{s},{\bm{c}}_{s})_{{\bm{U}}_{k}}$ of the augmented matrix in the SV-LCOV is obtained, according to Definition 1, this orthonormal basis corresponds to the matrix ${\bm{B}}-{\bm{U}}_{k}{\bm{C}}$ where the $i$-th column of ${\bm{B}}$ is ${\bm{b}}_{s}$. Regarding the final step of the Rayleigh-Ritz process for updating the truncated SVD by adding columns, we have the update procedure for the left singular vectors:

where ${\bm{F}}_{k}[:k]$ denotes the submatrix consisting of the first $k$ rows of ${\bm{F}}_{k}$, and ${\bm{F}}_{k}[k:]$ denotes the submatrix consisting of rows starting from the $(k+1)$-th rows of ${\bm{F}}_{k}$.

Equation (3) shows that each update of the singular vectors ${\bm{U}}_{k}$ consists of the following operations:

1. 1.

   A low-dimensional matrix multiplication to ${\bm{U}}_{k}$ by a $k$-by-$k$ matrix $({\bm{F}}_{k}[:k]-{\bm{C}}{\bm{F}}_{k}[k:])$.

2. 2.

   A sparse matrix addition to ${\bm{U}}_{k}$ by of a sparse matrix ${\bm{B}}{\bm{F}}_{k}[k:]$ that has at most $nnz({\bm{B}})\cdot k$ non-zero entries. While ${\bm{B}}{\bm{F}}_{k}[k:]$ may appear relatively dense in the context, considering it as a sparse matrix during the process of sparse matrix addition ensures asymptotic complexity.

An extended decomposition for reducing complexity. To reduce the complexity, we write the truncated SVD as a product of the following five matrices:

with orthonormal ${\bm{U}}={\bm{U}}^{\prime}\cdot{\bm{U}}^{\prime\prime}$ and ${\bm{V}}^{\prime}\cdot{\bm{V}}^{\prime\prime}$ (but not ${\bm{V}}^{\prime}$ or ${\bm{V}}^{\prime\prime}$), that is, the singular vectors will be obtained by the product of the two matrices. Initially, ${\bm{U}}^{\prime\prime}$ and ${\bm{V}}^{\prime\prime}$ are set by the $k$-by-$k$ identity matrix, and ${\bm{U}}^{\prime}$, ${\bm{V}}^{\prime}$ are set as the left and right singular vectors. Similar additional decomposition was used in Brand (2006) for updating of SVD without any truncation and with much higher complexity. With the additional decomposition presented above, the two operations can be performed to update the singular vector:

which is a low-dimensional matrix multiplication and a sparse matrix addition. And the update of the right singular vectors is

where ${\bm{G}}_{k}[:k]$ denotes the submatrix consisting of the first $k$ rows of ${\bm{G}}_{k}$ and ${\bm{G}}_{k}[k:]$ denotes the submatrix cosisting of rows starting from the $(k+1)$-th rows of ${\bm{G}}_{k}$. Now this can be performed as

Above we have only shown the scenario of adding columns, but a similar approach can be used to add rows and modify weights. Such an extended decomposition reduces the time complexity of updating the left and right singular vectors so that they can be deployed to large-scale matrix with large $m$ and $n$. In practical application, the aforementioned extended decomposition might pose potential numerical issues when the condition number of the matrix is large, even though in most cases these matrices have relatively low condition numbers. One solution to this issue is to reset the $k$-by-$k$ matrix to the identity matrix.

Algorithm 3 and Algorithm 4 are the pseudocodes of the proposed algorithm for adding columns and modifying weights, respectively.

The proposed method can theoretically produce the same output as Zha & Simon (1999) method with a much lower time complexity.

The proposed variants with the approximate augmented space. To address updates with coarser granularity (i.e., larger $s$), we propose two variants of the proposed algorithm based on approximate augmented spaces with GKL and RPI (see section 3.1) in SV-LCOV, denoted with the suffixes -GKL and -RPI, respectively. The proposed variants procude theoretically the same output as Vecharynski & Saad (2014) and Yamazaki et al. (2015), repectively. We elaborate the proposed variants in Appendix D.3.

Table. 1 presents the comparison of the time complexity of our proposed algorithm with the previous algorithms when updating rows ${{{\bm{E}}_{r}}}\in\mathbb{R}^{s\times n}$. To simplify the results, we have assumed $s<n$ and omitted the big-$O$ notation. $l$ denotes the rank of approximation in GKL and RPI. $t$ denotes the number of iteration RPI performed.

Our method achieves a time complexity of $O(nnz({\bm{E}}))$ for updating, without any $O(n)$ or $O(m)$ terms when $s$ and $k$ are constants (i.e., the proposed algorithm is at the update sparsity time complexity). This is because SV-LCOV is used to obtain the orthogonal basis, eliminating the $O(n)$ or $O(m)$ terms in processing the augmented matrix. Additionally, our extended decomposition avoids the $O(n)$ or $O(m)$ terms when restoring the SV-LCOV and eliminates the projection in all rows of singular vectors. Despite the increase in time complexity for querying from $O(k)$ to $O(k^{2})$, this trade-off remains acceptable considering the optimizations mentioned above.

For updates with one coarse granularity (i.e. larger $s$), the proposed method of approximating the augmented space with GKL or RPI in the SV-LCOV space eliminates the squared term of $s$ in the time complexity, making the proposed algorithm also applicable to coarse granularity update.

Baselines. We evaluate the proposed algorithm and its variants against existing methods, including Zha & Simon (1999), Vecharynski & Saad (2014) and Yamazaki et al. (2015). Throughout the experiments, we set $l$, the spatial dimension of the approximation, to $10$ based on previous settings (Vecharynski & Saad, 2014). The required number of iterations for the RPI, denoted by $t$, is set to $3$. In the methods of Vecharynski & Saad (2014) and Yamazaki et al. (2015), there may be differences in running time between 1) directly constructing the augmented matrix, and 2) accessing the matrix-vector multiplication as needed without directly constructing the augmented matrix. We conducted tests under both implementations and considered the minimum value as the running time.

Tasks and Settings. For the adjacency matrix, we initialize the SVD with the first 50% of the rows and columns, and for the user-item matrix, we initialize the SVD with the first 50% of the columns.

The experiment involves $\phi$ batch updates, adding $n/\phi$ rows and columns each time for a $2n$-sized adjacency matrix. For a user-item matrix with $2n$ columns, $n/\phi$ columns are added per update.

• •

  Link Prediction aims at predicting whether there is a link between a given node pair. Specifically, each node’s representation is obtained by a truncated SVD of the adjacency matrix. During the infer stage, we first query the representation obtained by the truncated SVD of the node pair $(i,j)$ to predict. Then a score, the inner product between the representation of the node pair

  $${\bm{U}}[i]^{\top}\Sigma{\bm{V}}[j]$$

  is used to make predictions. For an undirected graph, we take the maximum value in both directions. We sort the scores in the test set and label the edges between node pairs with high scores as positive ones.

  Following prior research, we create the training set $\mathcal{G^{\prime}}$ by randomly removing 30% of the edges from the original graph $\mathcal{G}$. The node pairs connected by the eliminated edges are then chosen, together with an equal number of unconnected node pairs from $\mathcal{G}$, to create the testing set $\mathcal{E}_{test}$.

• •

  Collaborative Filtering in recommender systems is a technique using a small sample of user preferences to predict likes and dislikes for a wide range of products. In this paper, we focus on predicting the values of in the normalized user-item matrix.

  In this task, a truncated SVD is used to learn the representation for each user and item. And the value is predicted by the inner product between the representation of $i$-th user and $j$-th item with

  $${\bm{U}}[i]^{\top}\Sigma{\bm{V}}[j]$$

  The matrix is normalized by subtracting the average value of each item. Values in the matrix are initially divided into the training and testing set with a ratio of $8:2$.

Datasets. The link prediction experiments are conducted on three publicly available graph datasets, namely Slashdot (Leskovec et al., 2009) with $82,168$ nodes and $870,161$ edges, Flickr (Tang & Liu, 2009) with $80,513$ nodes and $11,799,764$ edges, and Epinions (Richardson et al., 2003) with $75,879$ nodes and $1,017,674$ edges. The social network consists of nodes representing users, and edges indicating social relationships between them. In our setup, all graphs are undirected.

For the collaborative filtering task, we use data from MovieLens (Harper & Konstan, 2015). The MovieLens25M dataset contains more than $2,500,000$ ratings across $62,423$ movies. According to the dataset’s selection mechanism, all selected users rated at least $20$ movies, ensuring the dataset’s validity and a moderate level of density. All ratings in the dataset are integers ranging from $1$ to $5$.

To study the efficiency of the proposed algorithm, we evaluate the proposed method and our optimized GKL and RPI methods in the context of link prediction and collaborative filtering tasks.

Experiments are conducted on the undirected graphs of Slashdot (Leskovec et al., 2009), Flickr (Tang & Liu, 2009), and Epinions (Richardson et al., 2003) to investigate link prediction. The obtained results are presented in Table 2, Fig. 1 and Fig. 2. Our examination metrics for the Efficiency section include runtime and Average Precision (AP), which is the percentage of correctly predicted edges in the predicted categories to the total number of predicted edges. Besides, we report the Frobenius norm between the ${\bm{U}}_{k}\mathbf{\Sigma}_{k}{\bm{V}}_{k}^{\top}$ and the train matrix. The results demonstrate that across multiple datasets, the proposed method and its variants have demonstrated significant increases in efficiency without compromising AP, reducing time consumption by more than $85\%$.

Table 3 demonstrates the results of four experiments conducted for the collaborative filtering task: batch update($\phi=2000$) and streaming update($\phi=\#entries$) with $k=16$ and $k=64$. Our evaluation metrics include runtime and mean squared error (MSE). The results show that our method significantly reduces runtime while maintaining accuracy. Existing methods struggle with updating large and sparse datasets, especially for streaming updates, making real-time updates challenging. The methodology employed in our study effectively decreases the runtime by a factor of $10$ or more.

We conduct link prediction experiments on three datasets for $k$ and $\phi$, aiming to explore the selection of variants and approaches in various scenarios. The results in Fig. 1 show that our optimized methods are all significantly faster than the original ones for different choices of $k$. All methods become somewhat slower as $k$ increases, consistent with the asymptotic complexity metrics in Table 1.

We experimentally evaluate the performance of each method for different batch sizes, $\phi$. As shown in Fig. 2, our methods (ours-GKL and ours-RPI) outperform others when a large number of entries are updated simultaneously. For $\phi=10^{4}$, the three methods exhibit similar runtime. These experiments provide guidance for selecting the most suitable model based on specific context. For tasks with varying requirements for single updates, it is recommended to choose the most appropriate method.