DPar2: Fast and Scalable PARAFAC2 Decomposition for Irregular Dense Tensors

We use boldface lowercases (e.g. $\mathbf{x}$) and boldface capitals (e.g. $\mathbf{X}$) for vectors and matrices, respectively. In this paper, indices start at $1$. An irregular tensor is a 3-order tensor $\boldsymbol{\mathscr{X}}$ whose $k$-frontal slice $\boldsymbol{\mathscr{X}}(:,:,k)$ is $\mathbf{X}_{k}\in\mathbb{R}^{I_{k}\times J}$. We denote irregular tensors by $\{\mathbf{X}_{k}\}^{K}_{k=1}$ instead of $\boldsymbol{\mathscr{X}}$ where $K$ is the number of $k$-frontal slices of the tensor. An example is described in Fig. 2. We refer the reader to [19] for the definitions of tensor operations including Frobenius norm, matricization, Kronecker product, and Khatri-Rao product.

Singular Value Decomposition (SVD) decomposes $\mathbf{A}\in\mathbb{R}^{I\times J}$ to $\mathbf{X}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\text{T}}$. $\mathbf{U}\in\mathbb{R}^{I\times R}$ is the left singular vector matrix of $\mathbf{A}$; $\mathbf{U}=\begin{bmatrix}\mathbf{u}_{1}\cdots\mathbf{u}_{r}\end{bmatrix}$ is a column orthogonal matrix where $R$ is the rank of $\mathbf{A}$ and $\mathbf{u}_{1}$, $\cdots$, $\mathbf{u}_{R}$ are the eigenvectors of $\mathbf{A}\mathbf{A}^{\text{T}}$. $\mathbf{\Sigma}$ is an $R\times R$ diagonal matrix whose diagonal entries are singular values. The $i$-th singular value $\sigma_{i}$ is in $\mathbf{\Sigma}_{i,i}$ where $\sigma_{1}$ $\geq$ $\sigma_{2}$ $\geq$ $\cdots$ $\geq$ $\sigma_{R}$ $\geq$ $0$. $\mathbf{V}\in\mathbb{R}^{J\times R}$ is the right singular vector matrix of $\mathbf{A}$; $\mathbf{V}=\begin{bmatrix}\mathbf{v}_{1}\cdots\mathbf{v}_{R}\end{bmatrix}$ is a column orthogonal matrix where $\mathbf{v}_{1}$, $\cdots$, $\mathbf{v}_{R}$ are the eigenvectors of $\mathbf{A}^{\text{T}}\mathbf{A}$.

Randomized SVD. Many works [21, 20, 22] have introduced efficient SVD methods to decompose a matrix $\mathbf{A}\in\mathbb{R}^{I\times J}$ by applying randomized algorithms. We introduce a popular randomized SVD in Algorithm 1. Randomized SVD finds a column orthogonal matrix $\mathbf{Q}\in\mathbb{R}^{I\times(R+s)}$ of $(\mathbf{A}\mathbf{A}^{T})^{q}\mathbf{A}\mathbf{\Omega}$ using random matrix $\mathbf{\Omega}$, constructs a smaller matrix $\mathbf{B}=\mathbf{Q}^{T}\mathbf{A}$ $(\in\mathbb{R}^{(R+s)\times J})$, and finally obtains the SVD result $\mathbf{U}$ $(=\mathbf{Q}\tilde{\mathbf{U}})$, $\mathbf{\Sigma}$, $\mathbf{V}$ of $\mathbf{A}$ by computing SVD for $\mathbf{B}$, i.e., $\mathbf{B}\approx\tilde{\mathbf{U}}\mathbf{\Sigma}\mathbf{V}^{T}$. Given a matrix $\mathbf{A}$, the time complexity of randomized SVD is $\boldsymbol{\mathscr{O}}(IJR)$ where $R$ is the target rank.

PARAFAC2 decomposition proposed by Harshman [23] successfully deals with irregular tensors. The definition of PARAFAC2 decomposition is as follows:

The objective function of PARAFAC2 decomposition [23] is given as follows.

For uniqueness, Harshman [23] imposed the constraint (i.e., $\mathbf{U}_{k}^{T}\mathbf{U}=\Phi$ for all $k$), and replace $\mathbf{U}_{k}^{T}$ with $\mathbf{Q}_{k}\mathbf{H}$ where $\mathbf{Q}_{k}$ is a column orthogonal matrix and $\mathbf{H}$ is a common matrix for all the slices. Then, Equation (1) is reformulated with $\mathbf{Q}_{k}\mathbf{H}$:

Fig. 2 shows an example of PARAFAC2 decomposition for a given irregular tensor. A common approach to solve the above problem is ALS (Alternating Least Square) which iteratively updates a target factor matrix while fixing all factor matrices except for the target. Algorithm 2 describes PARAFAC2-ALS. First, we update each $\mathbf{Q}_{k}$ while fixing $\mathbf{H}$, $\mathbf{V}$, $\mathbf{S}_{k}$ for $k=1,...,K$ (lines 4 and 5). By computing SVD of $\mathbf{X}_{k}\mathbf{V}\mathbf{S}_{k}\mathbf{H}^{T}$ as $\mathbf{Z}^{\prime}_{k}\mathbf{\Sigma^{\prime}}_{k}\mathbf{P}_{k}^{{}^{\prime}T}$, we update $\mathbf{Q}_{k}$ as $\mathbf{Z}^{\prime}_{k}\mathbf{P}_{k}^{{}^{\prime}T}$, which minimizes Equation (3) over $\mathbf{Q}_{k}$ [11, 24, 25]. After updating $\mathbf{Q}_{k}$, the remaining factor matrices $\mathbf{H}$, $\mathbf{V}$, $\mathbf{S}_{k}$ is updated by minimizing the following objective function:

Minimizing this function is to update $\mathbf{H}$, $\mathbf{V}$, $\mathbf{S}_{k}$ using CP decomposition of a tensor $\boldsymbol{\mathscr{Y}}\in\mathbb{R}^{R\times J\times K}$ whose $k$-th frontal slice is $\mathbf{Q}_{k}^{T}\mathbf{X}_{k}$ (lines 8 and 10). We run a single iteration of CP decomposition for updating them [24] (lines 11 to 16). $\mathbf{Q}_{k}$, $\mathbf{H}$, $\mathbf{S}_{k}$, and $\mathbf{V}$ are alternatively updated until convergence.

Iterative computations with an irregular dense tensor require high computational costs and large intermediate data. RD-ALS [18] reduces the costs by preprocessing a given tensor and performing PARAFAC2 decomposition using the preprocessed result, but the improvement of RD-ALS is limited. Also, recent works successfully have dealt with sparse irregular tensors by exploiting sparsity. However, the efficiency of their models depends on the sparsity patterns of a given irregular tensor, and thus there is little improvement on irregular dense tensors. Specifically, computations with large dense slices $\mathbf{X}_{k}$ for each iteration are burdensome as the number of iterations increases. We focus on improving the efficiency and scalability in irregular dense tensors.

Before describing main ideas of our method, we present main challenges that need to be tackled.

• C1.

  Dealing with large irregular tensors. PARAFAC2 decomposition (Algorithm 2) iteratively updates factor matrices (i.e., $\mathbf{U}_{k}$, $\mathbf{S}_{k}$, and $\mathbf{V}$) using an input tensor. Dealing with a large input tensor is burdensome to update the factor matrices as the number of iterations increases.

• C2.

  Minimizing numerical computations and intermediate data. How can we minimize the intermediate data and overall computations?

• C3.

  Maximizing multi-core parallelism. How can we parallelize the computations for PARAFAC2 decomposition?

The main ideas that address the challenges mentioned above are as follows:

• I1.

  Compressing an input tensor using randomized SVD considerably reduces the computational costs to update factor matrices (Section III-B).

• I2.

  Careful reordering of computations with the compression results minimizes the intermediate data and the number of operations (Sections III-C to III-E).

• I3.

  Careful distribution of work between threads enables DPar2 to achieve high efficiency by considering various lengths $I_{k}$ for $k=1,...,K$ (Section III-F).

As shown in Fig. 3, DPar2 first compresses each slice of an irregular tensor using randomized SVD (Section III-B). The compression is performed once before iterations, and only the compression results are used at iterations. It significantly reduces the time and space costs in updating factor matrices. After compression, DPar2 updates factor matrices at each iteration, by exploiting the compression results (Sections III-C to III-E). Careful reordering of computations is required to achieve high efficiency. Also, by carefully allocating input slices to threads, DPar2 accelerates the overall process (Section III-F).

DPar2 (see Algorithm 3) is a fast and scalable PARAFAC2 decomposition method based on ALS described in Algorithm 2. The main challenge that needs to be tackled is to minimize the number of heavy computations involved with a given irregular tensor $\{\mathbf{X}_{k}\}^{K}_{k=1}$ consisting of slices $\mathbf{X}_{k}$ for $k=1,...,K$ (in lines 4 and 8 of Algorithm 2). As the number of iterations increases (lines 2 to 17 in Algorithm 2), the heavy computations make PARAFAC2-ALS slow. For efficiency, we preprocess a given irregular tensor into small matrices, and then update factor matrices by carefully using the small ones.

Our approach to address the above challenges is to compress a given irregular tensor $\{\mathbf{X}_{k}\}^{K}_{k=1}$ before starting iterations. As shown in Fig. 4, our main idea is two-stage lossy compression with randomized SVD for the given tensor: 1) DPar2 performs randomized SVD for each slice $\mathbf{X}_{k}$ for $k=1,...,K$ at target rank $R$, and 2) DPar2 performs randomized SVD for a matrix, the horizontal concatenation of singular value matrices and right singular vector matrices of slices $\mathbf{X}_{k}$. Randomized SVD allows us to compress slice matrices with low computational costs and low errors.

First Stage. In the first stage, DPar2 compresses a given irregular tensor by performing randomized SVD for each slice $\mathbf{X}_{k}$ at target rank $R$ (line 3 in Algorithm 3).

where $\mathbf{A}_{k}\in\mathbb{R}^{I_{k}\times R}$ is a matrix consisting of left singular vectors, $\mathbf{B}_{k}\in\mathbb{R}^{R\times R}$ is a diagonal matrix whose elements are singular values, and $\mathbf{C}_{k}\in\mathbb{R}^{J\times R}$ is a matrix consisting of right singular vectors.

Second Stage. Although small compressed data are generated in the first step, there is a room to further compress the intermediate data from the first stage. In the second stage, we compress a matrix $\mathbf{M}=\mathbin{\|}_{k=1}^{K}{\left(\mathbf{C}_{k}\mathbf{B}_{k}\right)}$ which is the horizontal concatenation of $\mathbf{C}_{k}\mathbf{B}_{k}$ for $k=1,...,K$. Compressing the matrix $\mathbf{M}$ maximizes the efficiency of updating factor matrices $\mathbf{H}$, $\mathbf{V}$, and $\mathbf{W}$ (see Equation (3)) at later iterations. We construct a matrix $\mathbf{M}\in\mathbb{R}^{J\times KR}$ by horizontally concatenating $\mathbf{C}_{k}\mathbf{B}_{k}$ for $k=1,...,K$ (line 5 in Algorithm 3). Then, DPar2 performs randomized SVD for $\mathbf{M}$ (line 6 in Algorithm 3):

where $\mathbf{D}\in\mathbb{R}^{J\times R}$ is a matrix consisting of left singular vectors, $\mathbf{E}\in\mathbb{R}^{R\times R}$ is a diagonal matrix whose elements are singular values, and $\mathbf{F}\in\mathbb{R}^{KR\times R}$ is a matrix consisting of right singular vectors.

With the two stages, we obtain the compressed results $\mathbf{D}$, $\mathbf{E}$, $\mathbf{F}$, and $\mathbf{A}_{k}$ for $k=1,...,K$. Before describing how to update factor matrices, we re-express the $k$-th slice $\mathbf{X}_{k}$ by using the compressed results:

where $\mathbf{F}^{(k)}\in\mathbb{R}^{R\times R}$ is the $k$th vertical block matrix of $\mathbf{F}$:

Since $\mathbf{C}_{k}\mathbf{B}_{k}$ is the $k$th horizontal block of $\mathbf{M}$ and $\mathbf{D}\mathbf{E}\mathbf{F}^{(k)T}$ is the $k$th horizontal block of $\mathbf{D}\mathbf{E}\mathbf{F}^{T}$, $\mathbf{B}_{k}\mathbf{C}_{k}^{T}$ corresponds to $\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}$. Therefore, we obtain Equation (6) by replacing $\mathbf{B}_{k}\mathbf{C}_{k}^{T}$ with $\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}$ from Equation (4).

In updating factor matrices, we use $\mathbf{A}_{k}\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}$ instead of $\mathbf{X}_{k}$. The two-stage compression lays the groundwork for efficient updates.

Our goal is to efficiently update factor matrices, $\mathbf{H}$, $\mathbf{V}$, and $\mathbf{S}_{k}$ and $\mathbf{Q}_{k}$ for $k=1,...,K$, using the compressed results $\mathbf{A}_{k}\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}$. The main challenge of updating factor matrices is to minimize numerical computations and intermediate data by exploiting the compressed results obtained in Section III-B. A naive approach would reconstruct $\tilde{\mathbf{X}}_{k}=\mathbf{A}_{k}\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}$ from the compressed results, and then update the factor matrices. However, this approach fails to improve the efficiency of updating factor matrices. We propose an efficient update rule using the compressed results to 1) find $\mathbf{Q}_{k}$ and $\mathbf{Y}_{k}$ (lines 5 and 8 in Algorithm 2), and 2) compute a single iteration of CP-ALS (lines 11 to 13 in Algorithm 2).

There are two differences between our update rule and PARAFAC2-ALS (Algorithm 2). First, we avoid explicit computations of $\mathbf{Q}_{k}$ and $\mathbf{Y}_{k}$. Instead, we find small factorized matrices of $\mathbf{Q}_{k}$ and $\mathbf{Y}_{k}$, respectively, and then exploit the small ones to update $\mathbf{H}$, $\mathbf{V}$, and $\mathbf{W}$. The small matrices are computed efficiently by exploiting the compressed results $\mathbf{A}_{k}\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}$ instead of $\mathbf{X}_{k}$. The second difference is that DPar2 obtains $\mathbf{H}$, $\mathbf{V}$, and $\mathbf{W}$ using the small factorized matrices of $\mathbf{Y}_{k}$. Careful ordering of computations with them considerably reduces time and space costs at each iteration. We describe how to find the factorized matrices of $\mathbf{Q}_{k}$ and $\mathbf{Y}_{k}$ in Section III-D, and how to update factor matrices in Section III-E.

The first goal of updating factor matrices is to find the factorized matrices of $\mathbf{Q}_{k}$ and $\mathbf{Y}_{k}$ for $k=1,...,K$, respectively. In Algorithm 2, finding $\mathbf{Q}_{k}$ and $\mathbf{Y}_{k}$ is expensive due to the computations involved with $\mathbf{X}_{k}$ (lines 4 and 8 in Algorithm 2). To reduce the costs for $\mathbf{Q}_{k}$ and $\mathbf{Y}_{k}$, our main idea is to exploit the compressed results $\mathbf{A}_{k}$, $\mathbf{D}$, $\mathbf{E}$, and $\mathbf{F}^{(k)}$, instead of $\mathbf{X}_{k}$. Additionally, we exploit the column orthogonal property of $\mathbf{A}_{k}$, i.e., $\mathbf{A}_{k}^{T}\mathbf{A}_{k}=\mathbf{I}$, where $\mathbf{I}$ is the identity matrix.

We first re-express $\mathbf{Q}_{k}$ using the compressed results obtained in Section III-B. DPar2 reduces the time and space costs for $\mathbf{Q}_{k}$ by exploiting the column orthogonal property of $\mathbf{A}_{k}$. First, we express $\mathbf{X}_{k}\mathbf{V}\mathbf{S}_{k}\mathbf{H}^{T}$ as $\mathbf{A}_{k}\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}\mathbf{V}\mathbf{S}_{k}\mathbf{H}^{T}$ by replacing $\mathbf{X}_{k}$ with $\mathbf{A}_{k}\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}$. Next, we need to obtain left and right singular vectors of $\mathbf{A}_{k}\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}$ $\mathbf{V}\mathbf{S}_{k}\mathbf{H}^{T}$. A naive approach is to compute SVD of $\mathbf{A}_{k}\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}\mathbf{V}\mathbf{S}_{k}\mathbf{H}^{T}$, but there is a more efficient way than this approach. Thanks to the column orthogonal property of $\mathbf{A}_{k}$, DPar2 performs SVD of $\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}\mathbf{V}\mathbf{S}_{k}$ $\mathbf{H}^{T}$ $\in\mathbb{R}^{R\times R}$, not $\mathbf{A}_{k}\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}\mathbf{V}\mathbf{S}_{k}\mathbf{H}^{T}\in\mathbb{R}^{I_{k}\times R}$, at target rank $R$ (line 9 in Algorithm 3):

where $\mathbf{\Sigma}_{k}$ is a diagonal matrix whose entries are the singular values of $\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}\mathbf{V}\mathbf{S}_{k}\mathbf{H}^{T}$, the column vectors of $\mathbf{Z}_{k}$ and $\mathbf{P}_{k}$ are the left singular vectors and the right singular vectors of $\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}\mathbf{V}\mathbf{S}_{k}\mathbf{H}^{T}$, respectively. Then, we obtain the factorized matrices of $\mathbf{Q}_{k}$ as follows:

where $\mathbf{A}_{k}\mathbf{Z}_{k}$ and $\mathbf{P}_{k}$ are the left and the right singular vectors of $\mathbf{A}_{k}\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}\mathbf{V}\mathbf{S}_{k}\mathbf{H}^{T}$, respectively. We avoid the explicit construction of $\mathbf{Q}_{k}$, and use $\mathbf{A}_{k}\mathbf{Z}_{k}\mathbf{P}_{k}^{T}$ instead of $\mathbf{Q}_{k}$. Since $\mathbf{A}_{k}$ is already column-orthogonal, we avoid performing SVD of $\mathbf{A}_{k}\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}\mathbf{V}\mathbf{S}_{k}\mathbf{H}^{T}$, which are much larger than $\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}\mathbf{V}\mathbf{S}_{k}\mathbf{H}^{T}$.

Next, we find the factorized matrices of $\mathbf{Y}_{k}$. DPar2 re-expresses $\mathbf{Q}_{k}^{T}\mathbf{X}_{k}$ (line 8 in Algorithm 2) as $\mathbf{Q}_{k}^{T}\mathbf{A}_{k}\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}$ using Equation (6). Instead of directly computing $\mathbf{Q}_{k}^{T}\mathbf{A}_{k}\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}$, we replace $\mathbf{Q}_{k}^{T}$ with $\mathbf{P}_{k}\mathbf{Z}_{k}^{T}\mathbf{A}_{k}^{T}$. Then, we represent $\mathbf{Y}_{k}$ as the following expression (line 12 in Algorithm 3):

Note that we use the property $\mathbf{A}_{k}^{T}\mathbf{A}_{k}=\mathbf{I}_{R\times R}$, where $\mathbf{I}_{R\times R}$ is the identity matrix of size $R\times R$, for the last equality. By exploiting the factorized matrices of $\mathbf{Q}_{k}$, we compute $\mathbf{Y}_{k}$ without involving $\mathbf{A}_{k}$ in the process.

The next goal is to efficiently update the matrices $\mathbf{H}$, $\mathbf{V}$, and $\mathbf{W}$ using the small factorized matrices of $\mathbf{Y}_{k}$. Naively, we would compute $\boldsymbol{\mathscr{Y}}$ and run a single iteration of CP-ALS with $\boldsymbol{\mathscr{Y}}$ to update $\mathbf{H}$, $\mathbf{V}$, and $\mathbf{W}$ (lines 11 to 13 in Algorithm 2). However, multiplying a matricized tensor and a Khatri-Rao product (e.g., $\mathbf{Y}_{(1)}(\mathbf{W}\odot\mathbf{V})$) is burdensome, and thus we exploit the structure of the decomposed results $\mathbf{P}_{k}\mathbf{Z}_{k}^{T}\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}$ of $\mathbf{Y}_{k}$ to reduce memory requirements and computational costs. In other word, we do not compute $\mathbf{Y}_{k}$, and use only $\mathbf{P}_{k}\mathbf{Z}_{k}^{T}\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}$ in updating $\mathbf{H}$, $\mathbf{V}$, and $\mathbf{W}$. Note that the $k$-th frontal slice of $\boldsymbol{\mathscr{Y}}$, $\boldsymbol{\mathscr{Y}}(:,:,k)$, is $\mathbf{P}_{k}\mathbf{Z}_{k}^{T}\mathbf{F}^{(k)}\mathbf{E}\mathbf{D}^{T}$.

Updating $\mathbf{H}$. In $\mathbf{Y}_{(1)}(\mathbf{W}\odot\mathbf{V})(\mathbf{W}^{T}\mathbf{W}*\mathbf{V}^{T}\mathbf{V})^{\dagger}$, we focus on efficiently computing $\mathbf{Y}_{(1)}(\mathbf{W}\odot\mathbf{V})$ based on Lemma 1. A naive computation for $\mathbf{Y}_{(1)}(\mathbf{W}\odot\mathbf{V})$ requires a high computational cost $\boldsymbol{\mathscr{O}}(JKR^{2})$ due to the explicit reconstruction of $\mathbf{Y}_{(1)}$. Therefore, we compute that term without the reconstruction by carefully determining the order of computations and exploiting the factorized matrices of $\mathbf{Y}_{(1)}$, $\mathbf{D}$, $\mathbf{E}$, $\mathbf{P}_{k}$, $\mathbf{Z}_{k}$, and $\mathbf{F}^{(k)}$ for $k=1,...,K$. With Lemma 1, we reduce the computational cost of $\mathbf{Y}_{(1)}(\mathbf{W}\odot\mathbf{V})$ to $\boldsymbol{\mathscr{O}}(JR^{2}+KR^{3})$.

As shown in Fig. 5, we compute $\mathbf{Y}_{(1)}(\mathbf{W}\odot\mathbf{V})$ column by column. In computing $\mathbf{G}^{(1)}(:,r)$, we compute $\mathbf{E}\mathbf{D}^{T}\mathbf{V}(:,r)$, sum up $\mathbf{W}(k,r)\left(\mathbf{P}_{k}\mathbf{Z}_{k}^{T}\mathbf{F}^{(k)}\right)$ for all $k$, and then perform matrix multiplication between the two preceding results (line 14 in Algorithm 3). After computing $\mathbf{G}^{(1)}\leftarrow\mathbf{Y}_{(1)}(\mathbf{W}\odot\mathbf{V})$, we update $\mathbf{H}$ by computing $\mathbf{G}^{(1)}(\mathbf{W}^{T}\mathbf{W}*\mathbf{V}^{T}\mathbf{V})^{\dagger}$ where ${\dagger}$ denotes the Moore-Penrose pseudoinverse (line 15 in Algorithm 3). Note that the pseudoinverse operation requires a lower computational cost compared to computing $\mathbf{G}^{(1)}$ since the size of $(\mathbf{W}^{T}\mathbf{W}*\mathbf{V}^{T}\mathbf{V})\in\mathbb{R}^{R\times R}$ is small.

Updating $\mathbf{V}$. In computing $\mathbf{Y}_{(2)}(\mathbf{W}\odot\mathbf{U})(\mathbf{W}^{T}\mathbf{W}*\mathbf{U}^{T}\mathbf{U})^{\dagger}$, we need to efficiently compute $\mathbf{Y}_{(2)}(\mathbf{W}\odot\mathbf{U})$ based on Lemma 2. As in updating $\mathbf{H}$, a naive computation for $\mathbf{Y}_{(2)}(\mathbf{W}\odot\mathbf{U})$ requires a high computational cost $\boldsymbol{\mathscr{O}}(JKR^{2})$. We efficiently compute $\mathbf{Y}_{(2)}(\mathbf{W}\odot\mathbf{U})$ with the cost $\boldsymbol{\mathscr{O}}(JR^{2}+KR^{3})$, by carefully determining the order of computations and exploiting the factorized matrices of $\mathbf{Y}_{(2)}$.

Updating $\mathbf{W}$. In computing $\mathbf{Y}_{(3)}(\mathbf{V}\odot\mathbf{H})(\mathbf{V}^{T}\mathbf{V}*\mathbf{H}^{T}\mathbf{H})^{\dagger}$, we efficiently compute $\mathbf{Y}_{(3)}(\mathbf{V}\odot\mathbf{H})$ based on Lemma 3. As in updating $\mathbf{H}$ and $\mathbf{V}$, a naive computation for $\mathbf{Y}_{(3)}(\mathbf{V}\odot\mathbf{H})$ requires a high computational cost $\boldsymbol{\mathscr{O}}(JKR^{2})$. We compute $\mathbf{Y}_{(3)}(\mathbf{V}\odot\mathbf{H})$ with the cost $\boldsymbol{\mathscr{O}}(JR^{2}+KR^{3})$ based on Lemma 3. Exploiting the factorized matrices of $\mathbf{Y}_{(3)}$ and carefully determining the order of computations improves the efficiency.