Singular Value Decomposition for general matrix (GESVJ)

Overview

One-sided Jacobi alogirhm is a classic and robust method to calculate Singular Value Decompositon (SVD), which can be applied for any dense matrix with the size of \(M \times N\). In this library, it is denoted as GESVJ, same as the API name used in LAPACK.

\[A = U \Sigma V^T\]

where \(A\) is a dense symmetric matrix of size \(m \times n\), \(U\) is \(m \times m\) matrix with orthonormal columns and and \(V\) is \(n \times n\) matrix with orthonormal columns, and \(\Sigma\) is diagonal matrix. The maximum matrix size supported in FPGA is templated by NRMAX, NCMAX.

Algorithm

The calculation process of one-sided Jacobi SVD method is as follows:

  1. Initialize matrix V = I
  2. Select two columns (i, j), i < j, of matrix A, namely \(A_i\) and \(A_j\). Accumulate two columns of data by fomular
\[\begin{split}&b_{ii} = A_i^TA_i = ||A_i||^2 \\ &b_{jj} = A_j^TA_j = ||A_j||^2 \\ &b_{ij} = A_i^TA_j\end{split}\]

A \(2 \times 2\) matrix can be obtained and noted as:

\[\begin{split}\begin{bmatrix} b_{ii}\ b_{ij}\\ b_{ji}\ b_{jj}\\ \end{bmatrix}\end{split}\]

where \(b_{ij}\) equals \(b_{ji}\).

  1. Solve the \(2 \times 2\) symmetric SVD with Jacobi rotation:
\[\begin{split}&\tau = (b_{ii} - b_{jj}) / (2 * b_{ij})) \\ &t = sign(\tau)/(|\tau| + \sqrt{(1 + \tau^2)}) \\ &c = 1 / \sqrt{(1+t^2)} \\ &s = c * t \\\end{split}\]

if we put s and c in a matrix as J, then J equals

\[\begin{split}\begin{bmatrix} 1\ \ \ \ \ 0\ ...\ 0\ \ 0\ \\ 0\ \ \ \ \ c\ ...\ s\ \ 0\ \\ 0\ \ \ \ \ 0\ ...\ 0\ \ 0\ \\ 0\ -s\ ...\ c\ \ 0\ \\ 0\ \ \ \ \ 0\ ...\ 0\ \ 1\ \\ \end{bmatrix}\end{split}\]
  1. Update the i and j columns of matrix A and V by
\[\begin{split}A = A J\\ V = V J\end{split}\]
  1. Calculate converage of \(2 \times 2\) matrix by
\[conv = |b_{ij}| / \sqrt{(b_{ii}b_{jj})}\]

and select the max converage of all pairs (i, j).

  1. Repeat steps 2-4 until all pairs of (i, j) are calculated and updated.
  2. If the max converage among all paris of (i, j) is bigger than 1.e-8, repeat steps 2-6 again, which is also called one sweep.
  3. When the converge is small enough, calculate the matrix U and S for each \(i = 1, ..., n\) by
\[\begin{split}&s_i =||a_i||_2 \\ &u_i = a_i / s_i\end{split}\]

Architecture

From the algorithm, we know that the core part of the computation is two columns of data read-and-accumulate, and then updating corresponding two columns of data in Matrix A and V. In this library, we implement this core module in the following architecture.

../../../_images/gesvj_structure.png

It can be seen from the architecture, steps 2-5 (each pair of i and j) of the algorithm is divided into three stages:

Stage 1:
  1. read two columns of data of A to BRAM and accumulate to \(b_{ii}\), \(b_{jj}\) and \(b_{ij}\).
  2. preload two columns of data of matrix V to BRAM.
Stage 2:
Calculate SVD for \(2 \times 2\) matrix
Stage 3:
  1. Update two columns of data in matrix A
  2. Update two columns of data in matrix V.
  3. Meanwhile, calculate converage for current pair (i, j).

Since operating data of matrix A and V are independent, two modules of stage 1 are running in parallel. Meanwhile, thress modules of stage 3 run in parallel. The last module of stage 3 calculates converage using \(2 \times 2\) matrix data. This converage computing process is in read-and-accu module of stage 1 according to the algorithm. However, it requires ~60 cycles, which is also a lot after partitioning matrix A by row. Therefore, this calculation process is extracted as a submodule in stage 3.

Note

Why updating matrix V is divided into two modules?

From the figure, we can see that there are two modules related to matrix V, preload two columns of data of V to BRAM, and updating V. In our design, matrix A and V are all saved in URAM. And for each URAM, only 2 ports of read/write are supported. Since matrix A cummulated \(2 \times 2\) data need 100+ cycles to do SVD. We may preload two columns of V into BRAMs via 2-reading ports of URAM. And using two writting ports when updating data in V.

Besides, in order to speed up the data reading and updating of matrix V data, the matrix V is partitioned by NCU through its row. For each CU, matrix V is read/written using 2 URAM ports.

Note

Supported data size:

The supported maximum size of matrix A that templated by NRMAX and NCMAX is 512. The partitioning number MCU and NCU can support up to 16, respectively.