Tucker Tensor Decomposition on FPGA

TL;DR

This paper introduces a FPGA-based hardware accelerator for Tucker tensor decomposition, optimizing key modules and demonstrating significant speedups over CPU and GPU implementations on real and synthetic data.

Contribution

It presents a novel FPGA implementation of Tucker decomposition modules, including a warm-start SVD algorithm, for resource-constrained hardware acceleration.

Findings

Achieved 2.16 to 30.2x speedup over CPU and GPU implementations.

Validated the design using synthetic and real MRI data sets.

Demonstrated effective hardware optimization for tensor computations on FPGA.

Abstract

Tensor computation has emerged as a powerful mathematical tool for solving high-dimensional and/or extreme-scale problems in science and engineering. The last decade has witnessed tremendous advancement of tensor computation and its applications in machine learning and big data. However, its hardware optimization on resource-constrained devices remains an (almost) unexplored field. This paper presents an hardware accelerator for a classical tensor computation framework, Tucker decomposition. We study three modules of this architecture: tensor-times-matrix (TTM), matrix singular value decomposition (SVD), and tensor permutation, and implemented them on Xilinx FPGA for prototyping. In order to further reduce the computing time, a warm-start algorithm for the Jacobi iterations in SVD is proposed. A fixed-point simulator is used to evaluate the performance of our design. Some synthetic data sets and a real MRI data set are used to validate the design and evaluate its performance. We compare our work with state-of-the-art software toolboxes running on both CPU and GPU, and our work shows 2.16 - 30.2x speedup on the cardiac MRI data set.

Tables

Table 1. TABLE I : Performance of the TTM unit.

$q$	$r$	LUTs	Registers	DSPs	Clock rate	Power
16	16	46,056	24,556	256	212MHz	2.008W
16	32	99,384	48,357	512	200MHz	2.395W
32	16	99,505	48,189	512	203MHz	2.503W
32	32	198,269	95,743	1,024	187MHz	3.141W

Table 2. TABLE II : Performance of the SVD unit (fixed point).

$p$	LUTs	Registers	DSPs	Clock rate	Power
16	8,711	12,284	128	209MHz	0.477W
32	11,134	13,965	256	207MHz	0.683W
64	16,127	17,532	512	208MHz	1.095W
128	25,360	24,631	1,024	203MHz	1.871W

Table 3. TABLE III : Performance of the tensor permute/reshape unit.

$q$	$p$	LUTs	Registers	Clock rate	Power
16	64	29,929	27,718	241MHz	1.342W
32	64	59,308	55,366	209MHz	1.961W
32	128	115,749	110,662	205MHz	2.981W

Equations

$$j=1+\sum _{m=1,m\ne k}^d\left((i_m-1)\prod _{n=1,n\ne k}^{m-1}{I}_n\right).$$

$$\text{mathbfcal}Y=\text{mathbfcal}X{\times }_k\mathbf{A}\in {\mathbb{R}}^{I_1\times \cdots I_{k-1}\times J\times I_{k+1}\cdots \times I_d}.$$

$${\mathbf{Y}}^{(k)}=\mathbf{A}{\mathbf{X}}^{(k)}$$

$$||\text{mathbfcal}X|{|}_{\mathrm{F}}=\sqrt{\sum _{i_1,i_2,\dots i_d}{x}_{i_1,i_2,\dots i_d}^2}.$$

$$\text{mathbfcal}X\approx \text{mathbfcal}G{\times }_1{\mathbf{A}}_1{\times }_2{\mathbf{A}}_2\times \cdots {\times }_d{\mathbf{A}}_d,$$

$$\text{mathbfcal}G=\text{mathbfcal}X{\times }_1{\mathbf{A}}_1^T{\times }_2{\mathbf{A}}_2^T\times \cdots {\times }_d{\mathbf{A}}_d^T.$$

$${\mathbf{X}}^{(k)}={\mathbf{U}}_k{\mathbf{S}}_k{\mathbf{V}}_k^T.$$

$$\underset{\{{\mathbf{A}}_k{\}}_{k=1}^d}{\min }||\text{mathbfcal}X-\text{mathbfcal}G{\times }_1{\mathbf{A}}_1{\times }_2{\mathbf{A}}_2\times \cdots {\times }_d{\mathbf{A}}_d|{|}_F$$

$$y_{\hat{i},r_j,i_{j+1}\dots i_d}=\sum _j^{I_j}{x}_{\hat{i},i_j,i_{j+1}\dots i_d}^{\prime }{\mathbf{A}}_j(i_j,r_j).$$

$${x^{\prime }}_{\hat{i},i_j,\dots ,i_d}\times a_{i_j,r_j},\mathrm{with}\hat{i}=l+nq\le {\hat{I}}^{\prime }=I_1{I}_2\cdots I_{j-1}.$$

$$y_{r_1,i_2,\dots ,i_d}=\sum _{i_1=1}^{I_1}{x}_{i_1,i_2,\dots ,i_d}^{\prime }{\mathbf{A}}_1(i_1,r_1).$$

$${x^{\prime }}_{i_1,i_2,\dots ,i_d}\times a_{i_1,r_1},\mathrm{with}{i}_1=l+nq\le I_1.$$

$$\left\{\begin{array}[]{ll}(p+1,q-1)&\text{if }q-p>2,\\ (1,p+q)&\text{if }q-p\leq 2\;\text{and}\;p+q\leq n,\\ (p+q+1-n,n)&\text{if }q-p\leq 2\;\text{and}\;n<p+q<2n-1,\\ (1,2)&\text{if }q=n\;\text{and}\;p=n-1.\end{array}\right.$$

$$T_{\text{total}}=\sum _{k=1}^d(\sum _{j=1}^k{T}_{\text{TTM},k,j}+T_{\text{SVD},k}+\sum _{i=1}^2{T}_{\text{permute},i,k}).$$

$$T_{\text{TTM},k,j}=\left\{\begin{array}[]{ll}I_{j}\prod\limits_{j^{\prime}=j+1}^{d}R_{j^{\prime}}\times\lceil{\frac{\prod\limits_{j^{\prime}=1}^{j-1}I_{j^{\prime}}}{q}}\rceil\times\lceil\frac{R_{j}}{r}\rceil&j>k\\ I_{j}I_{k}\prod\limits_{j^{\prime}=j+1,j^{\prime}\neq k}^{d}R_{j^{\prime}}\times\lceil{\frac{\prod\limits_{j^{\prime}=1}^{j-1}I_{j^{\prime}}}{q}}\rceil\times\lceil\frac{R_{j}}{r}\rceil&1<j<k\\ I_{k}\prod\limits_{j^{\prime}=2,j^{\prime}\neq k}^{d}R_{j^{\prime}}\times\lceil\frac{I_{1}}{q}\rceil\times\lceil\frac{R_{1}}{r}\rceil&j=1.\end{array}\right.$$

$$T_{\text{TTM},k,k}=\left\{\begin{array}[]{ll}I_{k}\prod\limits_{j^{\prime}=k+1}^{d}R_{j^{\prime}}\times\lceil{\frac{\prod\limits_{j^{\prime}=1}^{k-1}R_{j^{\prime}}}{q}}\rceil\times\lceil\frac{I_{k}}{r}\rceil&k\neq 1\\ \prod\limits_{j^{\prime}=2}^{d}R_{j^{\prime}}\times\lceil\frac{I_{1}}{q}\rceil\times\lceil\frac{I_{1}}{r}\rceil&k=1.\end{array}\right.$$

$$T_{\text{SVD},k}=I_k(I_k-1)\lceil \frac{I_k+R_{/k}}{p}\rceil$$

$$T_{\text{permute},1,k}\approx 5(I_k\lceil \frac{R_{/k}}{q^{\prime }}\rceil +R_{/k}\lceil \frac{I_k}{p^{\prime }}\rceil )$$

$$T_{\text{permute},2,k}\approx 5I_k(\lceil \frac{I_k}{q^{\prime }}\rceil +\lceil \frac{I_k}{p^{\prime }}\rceil )$$

$$\frac{\Vert \text{mathbfcal}X-\text{mathbfcal}G{\times }_1{\mathbf{A}}_1{\times }_2{\mathbf{A}}_2\times \cdots {\times }_d{\mathbf{A}}_d{\Vert }_{\mathrm{F}}}{\Vert \text{mathbfcal}A{\Vert }_{\mathrm{F}}}\times 100\%.$$

Full text

I Introduction

As a multi-dimension extension of matrices, tensors are a natural tool to represent and process multi-way data arrays [1]. For instance, an MRI sequence with three spatial dimensions and a time dimension can be naturally represented as a 4-way tensor. The convolution layer in a neural network is also a tensor. Leveraging various tensor decompositions [2, 3, 4], many high-dimensional data mining [5], machine learning [6, 7, 8, 9, 10, 11] and EDA [12, 13] problems have been solved efficiently without suffering from the curse of dimensionality.

Due to their superior performance in processing high-volume data, tensors have emerged as a promising tool to enable real-time machine learning and data analysis. Recently, tensor algorithms have achieved great success in training and compressing deep neural networks [6, 7, 8, 9, 10, 11]. The resulting tensorized models consume tremendously less memory, run-time and power than the original deep neural network models. Tensor algorithms have also been successfully employed to accelerate medical image analysis [14], anomaly detection [15] and speech recognition [16].

Despite the rapid progress of tensor algorithms, the hardware/algorithm co-design of tensor computation on resource-constrained platforms remains a new and (almost) unexplored field. Some tensor libraries [17, 18, 19] have been developed for high-performance platforms like clusters and super computers. However, little algorithm-architecture co-design targeting on power and cost-limited devices has been done, which has limited the application of tensor-based data analysis and machine learning on IoT and edge devices. Although many hardware accelerators are available for matrix and vector computations [20, 21] and have been applied to machine learning [22, 23, 24] and signal processing [25], they cannot handle tensor data, because the underlying theory and numerical procedures are fundamentally different. It is inefficient and error-prone to process tensor data by matrix- or vector-computation accelerators. Therefore, resource-constrained hardware accelerators for tensor computation are highly desired.

This paper presents, for the first time, a hardware accelerator for one of the most important tensor algorithms: Tucker decomposition [3]. Tucker decomposition is a high-order generalization of singular value decomposition (SVD) and principal component analysis (PCA), and it often achieves orders-of-magnitude higher data compression ratio than matrix compression algorithms on multi-way data. This method has been widely used in facial recognition [26], signal processing [27], deep learning [28] and data mining [5]. Tucker decomposition is often implemented via the high-order orthogonal iteration (HOOI) [29]. This algorithm involves some computation-intensive operations such as the tensor-times-matrix (TTM) and matrix SVD. Meanwhile, handling the huge amount of tensor data on FPGA or ASIC is a challenging task.

The contributions of our paper are summarized below:

• •

On the hardware side, we present an hardware architecture for Tucker decomposition. We describe the design and data communication of three units: TTM, SVD via Jacobi iterations, and tensor permutation/reshaping.

• •

On the algorithm side, we propose a warm-start algorithm to reduce the cost of the Jacobi iterations.

• •

We analyze the performance of our accelerator, implement it on a Xilinx FPGA, and show the implementation results.

• •

We compare our FPGA accelerator with some state-of-the-art algorithms on both CPU and GPU, and demonstrate its application on an MRI compression task. Our accelerator shows up to $30.2\times$ speedup on the MRI data set.

II Algorithm Background

In this section, we introduce the necessary background of tensors and Tucker decomposition.

II-A Tensors and Basic Tensor Operations

Notations: We use boldface lower-case letters (e.g., $\mathbf{x}$) to denote vectors, boldface capital letters (e.g. $\bf{X}$) to denote matrices, and boldface Euler script letters (e.g. $\mathbfcal{X}$) to denote tensors.

A tensor $\mathbfcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\cdots\times I_{d}}$ is a multidimensional data array. Here $d$ is called the order or way of $\mathbfcal{X}$. An integer $k\in[1,d]$ can be used as the index of a specific dimension or mode with size $I_{k}$. An entry of a tensor can be specified by an index vector. For instance, the $(i_{1},i_{2},\dots i_{d})$-th entry in tensor $\mathbfcal{X}$ is denoted by $x_{i_{1},i_{2},\dots i_{d}}$. Clearly, a vector and a matrix are order-1 and order-2 tensors, respectively.

Tensor Fiber and Slice: A fiber is a one-dimensional fragment of a tensor, obtained by fixing all indices but one. Tensor fibers are higher-order extension of matrix rows and columns. A third-order tensor has fibers that can be denoted by $\mathbf{x}_{:jk}$, $\mathbf{x}_{i:k}$ or $\mathbf{x}_{ij:}$ correspondingly. A tensor slice is a two-dimensional fragment of a tensor, obtained by fixing all indices but two. For instance, $\mathbf{X}_{i::}$, $\mathbf{X}_{:j:}$ and $\mathbf{X}_{::k}$ denote the horizontal, lateral, and frontal slices of a $3$-way tensor, respectively. Fig. 1 shows the slices and fibers of a tensor.

Tensor permutation: Tensor permutation changes the mode order of a tensor. It is a high-order extension of the matrix transpose. For instance, given $\mathbfcal{X}\in\mathbb{R}^{5\times 10\times 3}$, ${\bf permute}(\mathbfcal{X},[2,3,1])$ generates a new tensor $\mathbfcal{Y}\in\mathbb{R}^{10\times 3\times 5}$ with $y_{i_{2},i_{3},i_{1}}=x_{i_{1},i_{2},i_{3}}$.

Unfolding: Unfolding (or matricization) as shown in Fig. 2, is the process of transforming a tensor into a matrix. Unfolding reorders the elements of an $d$-way data array into a matrix. For instance, the mode-$k$ unfolding of a tensor $\mathbfcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\cdots\times I_{d}}$ is denoted as $\mathbf{X}^{(k)}\in\mathbb{R}^{I_{k}\times\prod\limits_{m\neq k}I_{m}}$. The mapping from the $(i_{1},\dots,i_{d})$-th element $\mathbfcal{X}$ to the $(i_{k},j)$-th element of $\mathbf{X}^{(k)}$ is given as follows

[TABLE]

TTM: Tensor-times-matrix (TTM) in short, is a high-dimensional expansion of matrix multiplication. Given a tensor $\mathbfcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\cdots\times I_{d}}$ and a matrix $\mathbf{A}\in\mathbb{R}^{J\times I_{k}}$, their $k$-mode product is a new tensor

[TABLE]

It can be expressed as matrix-matrix multiplication

[TABLE]

However, implementing it in this way directly can be inefficient on hardware because tensor permutation will be needed. This needs to move almost all the data in a tensor.

Tensor Norms: The norm of a tensor is a function which maps any tensor to a non-negative scalar. A widely used norm of tensors is the Frobenius norm, defined as

[TABLE]

II-B Tenser Tucker Decomposition

Given an $d$-way tensor $\mathbfcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\dots\times I_{d}}$, we may compress it by the Tucker decomposition [30]. As shown in Fig. 3, the Tucker decomposition approximates a large-size tensor with a small-size core tensor $\mathbfcal{G}\in\mathbb{R}^{R_{1}\times R_{2}\cdots\times R_{d}}$ and $d$ factor matrices $\{\mathbf{A}_{k}\in\mathbb{R}^{I_{k}\times R_{k}}\}_{k=1}^{d}$ as follows:

[TABLE]

where all columns are orthonormal inside each factor matrix $\mathbf{A}_{k}$. We call $\mathbf{R}=[R_{1},R_{2},\cdots,R_{d}]$ as the multilinear rank of $\mathbfcal{X}$. When $R_{k}\ll I_{k}$, $\mathbfcal{X}$ can be represented with the above Tucker format at a much lower cost. Once all factor matrices are obtained, the core tensor can be computed as

[TABLE]

Two popular methods can be used to compute a Tucker decomposition. The first well-known method is the high-order SVD (HOSVD) [3]. The idea of HOSVD is simple:

• •

One first unfolds $\mathbfcal{X}$ along mode $k$ to get $\mathbf{X}^{(k)}$;

• •

Then, perform a singular value decomposition (SVD)

[TABLE]

• •

Finally, pick $\mathbf{A}_{k}$ as the first $R_{k}$ columns of $\mathbf{U}_{k}$.

This method is easy to implement, however it is not optimal in fitting the data. Alternatively, an alternative least-square method called HOOI is widely used to get a better solution.

HOOI: The High Order Orthogonal Iteration (HOOI) [29, 31, 32] method aims to minimize the approximation error

[TABLE]

through the iterative process as shown in Alg. 1. Each iteration of the HOOI consists of two steps for every mode index $k$: (1) obtain a tensor $\mathbfcal{B}$ via a power iteration (TTM along all modes except $k$), which can be done from mode $d$, then $d-1$, … until mode 1. (2) an SVD of the mode-$k$ unfolded matrix of $\mathbfcal{B}$ to extract a mode-$k$ factor matrix $\mathbf{A}_{k}$.

In practice, the initialization process via HOSVD can be very time-consuming, because it needs $d$ SVD operations, and each of it works on a matrix whose size is equal to the original tensor. Therefore, some random orthonormal matrices are often used as the initial factor matrices. In this case, the total number of iterations needed may increase slightly, but the total runtime can decrease significantly. Even though, when the size of the tensor $\mathbfcal{X}$ is large, the time and energy comsumed to compute the Tucker decomposition can be very high.

III Proposed Architecture

In this section, we propose a new hardware architecture to perform Tucker decomposition via HOOI.

Fig. 4 shows the system-level diagram of our architecture. In order to load and accommodate the huge amount of tensor data elements, our design stores the tensor in an external DRAM. The HOOI engine consists of three parts: a tensor-times-matrix (TTM) unit, a singular value decomposition (SVD) unit, and a tensor permuting/reshaping unit. The data elements of the three units are stored in different memories. Because the size of a tensor can be very large, all tensors (including the intermediate and final results of TTM) are stored in an external DRAM. The matrix for each SVD is stored in an on-chip memory to reduce latency and to achieve a maximum throughput. Both parts are parallel and pipelined to achieve the maximum performance. The tensor permuting unit moves the data between the on-chip memory used by the SVD unit and the external DRAM used by the TTM unit, and it permutes the tensor accordingly.

III-A Tensor-Times-Matrix (TTM) Unit

The TTM unit can be implemented as a matrix-matrix multiplication that is available in some common computational libraries. However, we need to permute tensors before it can be implemented in this way. This is time- and memory-consuming because almost all $\prod\limits_{k=1}^{d}I_{k}$ data elements of $\mathbfcal{X}$ have to be moved. In this work, we develop a TTM unit without tensor permutations. For simplicity, the tensors are always stored by incrementing the mode-1 index $i_{1}$, then the second index $i_{2}$, and so forth. Since the size of tensor $\mathbfcal{X}$ is often beyond the capacity of an on-chip memory, all data elements (including the input and output data, and the intermediate result of a TTM operation) are stored in an external DRAM.

TTM with a 2-D PE Array. Our TTM unit is shown in Fig. 5. In order to maximize the throughput of the TTM module, we design a 2-D array of processing elements (PE) with $q$ columns and $r$ rows. Each PE computes the product of one scalar element of the tensor (black arrow) and one scalar from the matrix (either from the blue line or stored in the PE) at each clock cycle. The PEs in a single column are always connected to the same bus so they handle the same element in the tensor $\mathbfcal{X^{\prime}}$ each time. Here $\mathbfcal{X}^{\prime}$ can be either the original tensor $\mathbfcal{X}$ or an intermediate tensor after the TTM of $\mathbfcal{X}$ with some factor matrices. With $q$ columns, this module can handle up to $q$ neighbouring elements of a tensor fiber in total at the same time. Due to our method of storing $\mathbfcal{X^{\prime}}$, the fibers are always obtained by only changing the first mode index. Each row of this array handles one column of the factor matrix $\mathbf{A}_{j}$.

We compute the mode-$j$ TTMs in a decreasing order from $j=d$ to $j=1$. This PE array works in different ways dependent on the value of $j$, which is explained below.

• •

Assume that $j\neq 1$ and that we have done TTMs for all modes $>j$ except mode $k$. Let us ignore mode $k$ for simplicity, and the size of $\mathbfcal{X^{\prime}}$ becomes $I_{1}\times I_{2}\times\dots I_{j}\times R_{j+1}\times\dots\times R_{d}$. To simplify the expression, we fold the modes $1,2,\dots j-1$ into a single mode and use $\hat{i}$ as its index. The element-wise expression of this operation is

[TABLE]

In each clock cycle, each vertical bus can carry some neighbouring elements in a tensor fiber $\mathbf{x}^{\prime}_{:i_{j}\dots i_{d}}$ and each horizontal bus can carry an element in the factor matrix. Specifically, In the $n$-th cycle of the $i_{j}$-th round, the $l$-th vertical bus and the $r_{j}$-th horizontal bus carry scalars ${x^{\prime}}_{\hat{i},i_{j},\dots,i_{d}}$ and $\mathbf{A}_{j}(i_{j},r_{j})$, respectively. Consequently the $(r_{j},l)$-th PE calculates

[TABLE]

Finally $y_{\hat{i},r_{j},i_{j+1},\dots,d}$ is obtained by summing the above product terms over $i_{j}$ (in each round). Note that the index of matrix elements in $\mathbf{A}_{j}$ used at each PE does not depend on $l$. Therefore, we can multiply the whole fiber $\mathbf{x}^{\prime}_{:i_{j}\dots i_{d}}$ with the same matrix element, and all PEs in the same row share the same data elements of $\mathbf{A}_{j}$ by connecting them to the same horizontal bus. Because the size the new dimension, $\hat{I}=I_{1}I_{2}\dots I_{k-1}$, is very large, we divide the partly-folded tensor $\mathbfcal{X}^{\prime}$ into some small sub-tensors of size $m\times I_{j}\times R_{j+1}\dots\times R_{d}$, and the resulting tensor $\mathbfcal{Y}$ into some sub-tensors of size $m\times R_{j}\times R_{j+1}\dots\times R_{d}$. We can compute the sub-tensors of $\mathbfcal{Y}$ one by one to reduce the buffer size.

• •

When computing $\mathbfcal{X^{\prime}}\times_{1}\mathbf{A}_{1}^{T}$, the element-wise result is

[TABLE]

The $r_{1}$-th row of the PE array computes $q$ product terms of (10). Because $q$ is usually smaller than the fiber length, we need to compute all product terms by several cycles. In the $n$-th cycle, the $(r_{1},l)$-th PE calculates

[TABLE]

Note that in mode-1 TTM, the RAM inside the $(r_{1},l)$-th PE stores all $a_{i_{1},r_{1}}$’s for all $i_{1}=l+nq\leq I_{1}$. The TTM element $y_{r_{1},i_{2},\dots i_{d}}$ is obtained by accumulating all products terms in the $r_{1}$-th row during all cycles.

Processing Element (PE). As shown in Fig. 6, each PE consists of a multiplier, a small RAM, and another memory used as an output buffer storing the result before writing it to a DRAM. The RAM (marked in red) is only used to store $\mathbf{A}_{1}$ when computing the mode-1 TTM. Otherwise, the bus at each row imports data elements of $\mathbf{A}_{j}$ when $j\neq 1$. Therefore, a mutiplexer (MUX) is used to select the correct data to perform product operations. After computing one batch of the data (a tensor fiber if $j=1$ and a tensor slice if $j\neq 1$), the result is written to the DRAM and then reset to zero. The buffer temporarily holds the intermediate results, and keeps updating during the multiplication and sum operations. The buffer stops updating when its data is written into an external memory. In order to avoid timing conflicts and to increase throughput, another buffer is used (not shown in the figure) for transferring data to an external memory. These two buffers switch their roles after processing each batch of data.

In-place Adder Tree. As mentioned above, we need an adder tree for the mode-1 TTM. If the adder tree is implemented as a pipeline, a lot of adders and registers will be used. Given that the product terms need to be summed up only once per batch of data instead of per clock cycle, we only need an in-place adder tree. We split a adder tree into multiple stages. After each stage, each two elements are summed up so the total number of data elements is reduced by $50\%$. The registers and adders are shared among different stages. The data elements are read from and write back to the same group of registers after each clock cycle. This is why we call it an “in-place” adder tree. This treatment can reduce the number of adders and registers by $50\%$.

III-B Singular Value Decomposition (SVD) Unit

Our SVD unit employs the Jacobi iterations [33, 34, 35, 36, 37]. Both single-side and double-side Jacobi iterations are widely used. We use the single-side version because of its higher accurate, ease to parallelize and less data dependency. The whole framework is summarized in Alg. 2. Given a matrix $\mathbf{B}$, this algorithms computes its left singular vectors by orthogonizing every two rows (i.e., $\mathbf{b}_{i,:}$ and $\mathbf{b}_{j,:}$) of the matrix iteratively. The iterative process involves the norms and inner products of the row vectors and performing some rotations.

We can select the order of $(i,j)$ in different ways. A natural choice is to increment $j$ in the inner loop and increment $i$ in the outer loop, or vice versa. However, because of the data dependency between two adjacent operations, this choice makes it impossible to implement parallel or pipe-lined design. In order to overcome this issue, we employ the round-robin ordering suggested in [38], which eliminates the data dependency between adjacent iterations. This method starts by dividing all indices into $n/2$ pairs $(1,2),(3,4),\dots(n-1,n)$ where $n$ is the total number of rows. After orthogonalizing all the pair of rows specified above, we generate new index pairs in this way: suppose the pair in the previous round is $(p,q)$, this pair index is updated in the next round as

[TABLE]

Fig. 7 shows the block diagram of our SVD unit. In each step, two vectors are orthogonalized. The on-chip memory provides two ports to operate independently. In this part, one port takes the two vectors from the memory, and another port writes the orthogonalized vectors back into the memory. Because we use only one port to input data and the other for output, the two vectors have to be fetched through the same port alternatively. Three sum-of-products are needed to calculate the rotation angle $\theta$. Given that the two vectors are fetched alliteratively, the multiplier and adder to get $\|\mathbf{b}_{i:}\|_{2}$ and $\|\mathbf{b}_{j:}\|_{2}$ can be shared. Therefore, only two sets of multipliers and adders are used. In order to get $\theta$, $\sin\theta$ and $\cos\theta$, we employ the CORDIC algorithm [39], which uses the rotations of some fixed angles to approximate the rotation of any angle. This algorithm is efficient to calculate the trigonometric functions on hardware, and an FPGA IP core is available. After these two vectors are fetched from the memory, they are stored in a local buffer until $\sin\theta$ and $\cos\theta$ are calculated, then they are rotated accordingly. In this way, it is guaranteed each vector is read from the memory only once when orthogonalizing each pair of rows.

Besides the input matrix $\mathbf{B}$, the orthogonal matrix $\mathbf{U}$ also needs to be rotated in the same way. We store $\mathbf{U}$ and $\mathbf{B}$ in the same memory and handle them in the same way, except that $\mathbf{U}$ is not used for calculating $\alpha,\beta$ and $\gamma$. Since $\mathbf{U}$ has a much smaller size than $\mathbf{B}$, such a design only causes negligible run-time overhead but saves half of the area and power.

III-C Tensor Permutation and Reshaping

Once $\mathbfcal{B}=\mathbfcal{X}\times_{1}\mathbf{A}_{1}^{T}\dots\times_{k-1}\mathbf{A}_{k-1}^{T}\times_{k+1}\mathbf{A}_{k+1}^{T}\dots\times_{d}\mathbf{A}_{d}^{T}$ is computed, we need to permute and reshape the tensor $\mathbfcal{B}$ into $\mathbf{B}^{(k)}$ before performing an SVD. Since $\mathbfcal{B}$ is stored in an external DRAM and the matrix $\mathbf{B}^{(k)}$ is stored in an on-chip memory, we need an extra module to move the data between the DRAM and the on-chip memory and re-organize the data elements to match $\mathbf{B}^{(k)}$. After SVD, the factor matrix need to be transposed and moved from the on-chip memory to DRAM, which is done by this module as well.

When moving the data from on-chip memory and external DRAM, the data is first read from the its original memory, written to a local buffer with size $p^{\prime}\times q^{\prime}$, then written to the destination memory. The size of the buffer determines the size of data set that can be moved in every batch.

IV Implementation Details

IV-A Fixed Point Number

Floating-point numbers usually cause higher hardware cost and longer latency. Therefore, we use a fixed-point number system based on the trade-off between accuracy and hardware complexity.

We decide the fixed-point precision based on some hardware constrains. Because the data width at the interface of a DRAM controller is fixed as 512 bits, the memory is most efficient if the data width is a factor of 512 (i.e., $2^{n}$ with integer $n\leq 9$). Meanwhile, each DSP our FPGA can perform an multiplication with data sizes up to 27 bits $\times$ 18 bits. Considering these constraints, we use 16-bit numbers to represent all tensor data elements, and store them in an external DRAM. On the other hand, we use 27-bit numbers to represent the matrix data in both TTM and SVD in order to achieve higher accuracy and to avoid excessive use of multipliers. Note that we use a smaller data width for tensor data in order to save some memory space when processing the huge amount of data in a tensor. In this case, each multiplier in the TTM unit requires one DSP block, and each multiplier in the SVD unit requires two DSP blocks.

There are many sum-of-product operations in both TTM and SVD. The small error in the product terms will accumulate when calculating the sum. In order to address this issue, we use 48-bit numbers to represent the product terms. We use 27-bit numbers to represent most of the intermediate results, except for the 32-bit $\alpha,\beta,\gamma,\theta$ in the SVD unit.

IV-B HOOI with A Warm-start Algorithm

We observe that the SVD $\mathbf{B}^{(k)}=\mathbf{U}_{k}\mathbf{S}_{k}\mathbf{V}_{k}$ via the Jacoboi iterations is the most time-consuming part of HOOI. Therefore, we employ a warm-start strategy to reduce the number of Jacobi iterations. In the standard Jacobi iterations, the initial guess is chosen as an identity matrix. In our implementation, we use the orthonormal matrix $\mathbf{U}_{k}$ obtained from the previous iteration of HOOI as the initial guess of the Jacobi iteration. Thanks to the good initial guess, we only need to perform one or two rounds of Jacobi iterations inside each SVD, and the whole HOOI still converges after an enough number of power iterations. The optimized algorithm is shown in Alg. 3.

Fig. 8 shows the convergence curves of the standard HOOI and our proposed warm-start HOOI, respectively. We consider a tensor of size $128\times 128\times 128$ and with a multi-linear rank $\mathbf{R}=[32,32,32]$, which is generated by a Gaussian distribution and corrupted with some Gaussian noise. The noise variance is half of that of the tensor element. The standard method converges after only two HOOI iterations, but the SVD of each mode requires about 10 rounds of Jacobi iterations. Our optimized HOOI converges after 7 iterations, but each SVD requires only one round of Jacobi iterations, leading to a significant reduction of the total cost.

V Performance Analysis

This section analyzes the hardware performance of our FPGA-accelerated Tucker decomposition.

V-A Run-time Analysis

The total runtime is the sum of each part: TTM, SVD, and tensor permuting. For a $d$-way tensor, each HOOI iteration requires $d(d-1)$ TTMs, $d$ SVDs and $2d$ tensor permuting/reshaping operations. Some intermediate results can be reused. For instance, after computing $\mathbfcal{X}\times_{3}\mathbf{A}_{3}\times_{2}\mathbf{A}_{2}$ for a 3-way $\mathbfcal{X}$, the result of $\mathbfcal{X}\times_{3}\mathbf{A}_{3}$ can be reused when computing $\mathbfcal{X}\times_{3}\mathbf{A}_{3}\times_{1}\mathbf{A}_{1}$. By considering the TTM data reuse, the actual total number of TTMs is reduced to $(d-1)(1+d/2)$. When applying the warm-start algorithm for Jacobi iteration, the unfolded matrix $\mathbf{B}^{(k)}$ need to be multiplied with $\mathbf{U}_{k}$ first, and this is done by TTM as well, causing additional $d$ TTM operations. The total run-time is given by

[TABLE]

The details of each term are provided below.

TTM Part: The run-time of TTM depends on the size of its multiplier matrix. Suppose that we have a multiplier matrix of size $q\times r$. In this case, the multiplier takes in $q$ elements of $\mathbfcal{X}^{\prime}$ per clock cycle, and each element is multiplied with $r$ elements of the factor matrix $\mathbf{A}_{j}$. At most $qr$ product and sum operations can be done per clock cycle. Therefore, the number of clock cycles is

[TABLE]

Similarly, the clock cycles required for $\mathbf{U}^{(k)}\mathbf{B}_{k}$ is

[TABLE]

Some extra clock cycles are caused by the latency of the pipeline and ping-pong buffer, but they are often less than $1\%$ of the total run-time and thus negligible. When the TTM is applied over the first mode, the data need to be copied to the local memory of each PE in advance. This causes extra $O(I_{k}\lceil\frac{J_{k}}{r}\rceil)$ clock cycles, which is negligible again and can be done in parallel with other operations.

SVD Part: When updating the $k$-th factor matrix, we need to do SVD to a matrix with size $I_{k}\times R_{/k}$, with $R_{/k}=\prod_{i\neq k}R_{i}$. In each Jacobi iteration, $I_{k}(I_{k}-1)/2$ pairs of matrices will be orthogonalized. Each orthogonalization handles $2(I_{k}+R_{/k})$ numbers, therefore it takes $2\lceil\frac{I_{k}+R_{/k}}{p}\rceil$ clock cycles, where $p$ is the degree of parallelism. As a result, each Jacobi iteration takes approximately

[TABLE]

clock cycles. Similar to the case in TTM, the extra time cased by the latency of pipeline is negligible.

Tensor Permutation: Suppose that the size of the buffer is $p^{\prime}\times q^{\prime}$. In each clock cycle, this buffer can either exchange (read or write) $p^{\prime}$ numbers with the internal memory or $q^{\prime}$ numbers with the DRAM. Note that $p^{\prime}$ and $q^{\prime}$ are not necessarily equal to $p$ and $q$, respectively, as long as the memory supports writing $p^{\prime}$ or $q^{\prime}$ elements each time. However, setting $p^{\prime}=p,q^{\prime}=q$ can simplify our design and maximize the hardware efficiency. Each tensor permutation requires $O\left(I_{k}\lceil\frac{R_{/k}}{q^{\prime}}\rceil+R_{/k}\lceil\frac{I_{k}}{p^{\prime}}\rceil\right)$ clock cycles, with $R_{/k}=\prod_{i\neq k}R_{i}$. Our simulation shows that the practical run-time is

[TABLE]

clock cycles when we move $\mathbfcal{B}$ from a DRAM to an on-chip memory and permute it. We need

[TABLE]

clock cycles to move $\mathbf{U}_{k}$ to DRAM and transpose it.

V-B Area and Power

The area and power depends on the design parameters $p$, $q$ and $r$. The higher the degree of parallelism is, the more PEs, hardware area and power will be required. In the TTM unit, the total number of multipliers, adders, accumulators and buffers are proportional to the size of multiplier matrix, $q\times r$. Therefore, the area of TTM is approximately $O(q\times r)$. Similarly, the area of the tensor permutation unit is proportional to the buffer size $p\times q$. The power also increases when the area increases.

The area of the SVD unit is independent of its input matrix size, but depends on the degree of parallelism $p$. Additionally, one arctan module and one sin/cos module are required. Therefore, the area of the SVD unit is estimated as $c_{1}p+c_{2}$, where $c_{1}$ represents the area (multipliers, adders, accumulators, etc.) required for each matrix element, and $c_{2}$ represents the area of arctan and sin/cos blocks.

VI Results

VI-A FPGA Performance Validation

We implement all parts of the optimized HOOI with different design parameters on FPGA. All the results below are based on Xilinx XCVU9P-L2FSGD2104E FPGA, which is available on a Xilinx VCU1525 acceleration board. The power is estimated with a 200-MHz clock rate. The results for different units are shown in Tables I-III. The hardware complexity, including the number of lookup tables (LUTs), registers and DSP blocks, is determined by the design parameters. The area of the TTM unit is approximately proportional to $q\times r$. The area of the SVD unit is approximately proportional to $p$. In tensor permute unit, its area is approximately proportional to $p^{\prime}\times q^{\prime}$. The power consumption increases as we increase the design parameters but the relationship is not strictly linear, since the power consumption of some parts (e.g., the clock generator) is independent of our design parameters.

VI-B Performance Comparisons

We compare our FPGA accelerator with other two commonly used toolboxes on different platforms: the Tensor toolbox [40, 41] in MATLAB on CPU, and the TensorLy toolbox [42] in Python on both CPU and GPU. The TensorLy can select NumPy, PyTorch or MXNet as its backend, and the latter two allow high-performance GPU computation. In our experiment, we use NumPy and PyTorch for comparison. The runtime on CPU is measured on a Linux desktop with Intel i5-8400 6-core CPU and 32 GB memory. The GPU experiments are conducted on a Titan V GPU. Since the accuracy and convergence of our Tucker decomposition depends on the fixed-point number system, we develop a fixed-point simulator with the Xilinx fixed-point number library to simulate the truncation error and overflow in our FPGA accelerator.

We perform the comparison by using some randomly generated low-rank tensors. For each tensor, both the core tensor and the factor matrices are generated by Gaussian distributions, and the tensor is then corrupted by some Gaussian random noise. The relative error is evaluated as

[TABLE]

Fig. 9 shows the results on a set of 3-way tensors with size $I_{1}=I_{2}=I_{3}$. Our architecture can get $1.41-9.90\times$ speedup compared with MATLAB tensor toolbox on CPU, and even more compared with the TensorLy toolbox on both CPU and GPU.The convergence behavior of our FPGA-based Tucker decomposition is shown in Fig. 9(a). It is shown that our method always converges after 6-8 HOOI iterations. Therefore, we assume 8 HOOI iterations to estimate the runtime of our FPGA architecture.

We further perform comparisons using some 4-way tensors and show the results in Fig. 10. Our PC with 32GB RAM can no longer accommodate such 4-way tensor data, therefore the results on CPU are obtained by running our experiments on Amazon AWS r4.4x large instance with 16 virtual CPUs (vCPU) and 122GB memory. Since large 4-way tensors exceed the memory of GPU, and TensorLy with the NumPy backend consumes extremely long time, their runtimes are not shown in Fig. 10. When the size along each dimension is 256, MATLAB run out of memory when computing. On these 4-way tensor data, our FPGA design can get $3.18-8.22\times$ speedup compared with the MATLAB tensor toolbox on the AWS CPU.

Remark: Our FPGA accelerator uses the Jacobi iteration to perform SVD, whereas the MATLAB tensor toolbox and TensorLy use more powerful advanced SVD algorithms. If the same SVD algorithms are used in all implementations, our FPGA design may achieve more significant speedup.

VI-C Application: MRI Compression

The accelerated Tucker decomposition can be applied to multiple application domains. Here we demonstrate its application in compressing multi-way MRI data. This dataset is a single-channel, first-pass myocardial perfusion real-time MRI sequence ($n_{x}=190,n_{y}=90,n_{t}=70$). We use the pre-processed data available in [43, 44], and set rank to $R_{x}=40,R_{y}=32,R_{t}=28$ in HOOI. The estimated runtime with our architecture is 0.0447s at a clock rate of 185MHz. In comparison, on a Linux PC with Intel i5-8500 CPU 6 core CPU, the same operation takes 0.0964s with the MATLAB tensor toolbox, 0.335s with TensorLy toolbox and NumPy backend, 1.352s with TensorLy toolbox and Pytorch backend. on a PC with NVIDIA TITAN V GPU, this operation takes 0.217s. Therefore, our method provides 2.16-30.2$\times$ speedup compared with existing frameworks on CPU, and 4.85$\times$ speedup compared with GPU. Fig. 11 shows the original MRI data and the data approximated by our low-rank Tucker decomposition on FPGA.

VII Conclusion

This paper has presented an algorithm-architecture co-design to perform tensor Tucker decomposition. We have implemented Tensor-Times-Matrix, matrix SVD with single side Jacobi iteration, and tensor permuting on FPGA. We have also proposed a warm-start algorithm for the Jacobi iterations to reduce its computation time. We have done simulations on both synthetic data sets and an MRI data set. Our method is significantly faster than existing computation frameworks on both CPU and GPU. This accelerator can be employed in a broad range of applications including data mining, scientific computing, computer vision, and deep learning.

Acknowledgment

This work was supported by NSF CCF 1817037. We would like to thank Sophia Shao (NVIDIA), Lei Deng, Zheng Qu and Bangyan Wang (UCSB) for their technical discussions.