An answer to an open question in the incremental SVD

Yangwen Zhang
Department of Mathematics Science, Carnegie Mellon University, Pittsburgh, PA, USA ().

Incremental singular value decomposition (SVD) was proposed by Brand to efficiently compute the SVD of a matrix. The algorithm needs
to compute  thousands or millions of orthogonal matrices and to multiply them together. However, the multiplications may corrode the orthogonality. Hence many reorthogonalizations are needed in practice.  In [Linear Algebra and its Applications 415 (2006) 20–30], Brand asked “It is an open question how often this is necessary to guarantee a certain overall level of numerical precision; it does not change the overall complexity.” In this paper, we answer this question and the answer is we can avoid computing the large amount of those  orthogonal matrices and hence the  reorthogonalizations are not necessary by modifying his algorithm. We prove that the modification does not change the outcomes of the algorithm. Numerical experiments are presented to illustrate the performance of our modification.

§ INTRODUCTION

The singular value decomposition (SVD) of a matrix has many applications, such as proper orthogonal decomposition model order reduction  and  principal component analysis.

One drawback of the SVD is the cost of its computation.  Let $U$ be a $m$ by $n$ dense matrix of low rank $r$,  the  computational complexity of traditional methodsis $O(mn^2+m^2 n + n^3)$ time, which is unfeasible for a large size matrix. Lanczos methods [2] yield thin svds and the complexity is $O(mnr^2)$ time, but we need to know the rank in advance.
These  methods are referred to as batch methods because they need to store the matrix.

In some scenarios the data sets will be produced incrementally, such as the snapshots of a time dependent partial differential equations (PDEs). It may be advantageous to perform the SVD as the columns of a matrix  become available, instead of waiting until all data sets are available before doing any computation. These characteristics on the availability of data sets have given rise to a class of incremental methods.

In 2002, Brand [3] proposed a new  algorithm to find the SVD of  a matrix incrementally.  Given the SVD of a matrix $U = Q\Sigma R^\top$, the goal is to update  the SVD of the related matrix $\left[
U ~\big|~	c\right]$, by using the SVD of $U$ and the new adding vector $c$.  Then the SVD of $\left[
U ~\big|~	c\right]$ can be constructed by

1. letting  $e = c - QQ^\top c$ and $p = \|e\|$,

2. finding the full SVD of $\left[\begin{array}{cc}
\Sigma & Q^\top c \\
0 & p
\end{array}\right] = \widetilde  Q \widetilde \Sigma \widetilde R^\top$, and then

3. updating the SVD of $\left[
U ~\big|~	c\right]$ by

\left[
U ~\big|~	c\right]=\left(\left[
Q ~\big|~  e/p \right] \widetilde Q\right) \widetilde \Sigma \left(\left[\begin{array}{cc}
R & 0 \\
0 & 1
\end{array}\right] \widetilde R \right)^{\top}.

For many of the motivating applications, only the dominant singular vectors and values of $U$ are needed. Hence, in practice, we perform truncation when $p$ or the last singular value  of $\widetilde \Sigma$ is small.  The incremental SVD  has been applied to many different areas; see [8, 10, 11, 12] for more details.  There are many other incremental methods, see [1] for a survey.

There are two ways to update the left and right singular subspace. The first way is very straightforward and the update is given by

\begin{align}\label{rankaddone3}
Q\longleftarrow \left[
Q~ \big|~   e/p
\right]\widetilde Q,  \qquad   \Sigma \longleftarrow  \widetilde \Sigma, \qquad   R \longleftarrow  \left[\begin{array}{cc}
R & 0 \\
0 & 1
\end{array}\right] \widetilde R.
\end{align}

However, this update has two defects:

(1) The computational cost of (<ref>) is expensive since the updates involves a very  tall thin
matrix $Q$ and $R$.

(2) Reorthogonalizations are needed for $Q$ almost after each update.	Theoretically, the matrix $\left[
Q  ~\big|~ e/p\right] \widetilde Q$ is  orthogonal. However, in practice, small numerical errors cause a loss of orthogonality. Without a reorthogonalization, the decomposition of the incremental SVD algorithm is not  orthogonal. Therefore,  the singular values are not true and  truncation is not reliable. Hence a reorthogonalization is needed and the computational cost can be high if the rank of the matrix is not extremely low. Specifically, Fareed et al.  recommended a Gram-Schmidt orthogonalization [7] while Oxberry et al.  chose thin QR  in [13].

Numerical experiments in <Ref> show that the above two items take the majority  of the CPU time.

To avoid the huge computational cost, Brand [4] suggested  to update the matrix $Q$ and $R$ indirectly.  His idea is to leave  the SVD decomposed into the five matrices

\begin{align*}
U = Q_{\texttt{old}}  \widetilde Q_{\texttt{old}}  \Sigma_{\texttt{old}}  \widetilde R_{\texttt{old}}^\top  R_{\texttt{old}}^\top
\end{align*}

with orthonormal $Q_{\texttt{old}}  \widetilde Q_{\texttt{old}}  $, $R_{\texttt{old}} \widetilde R_{\texttt{old}} $, $Q_{\texttt{old}}$ and $ \widetilde Q_{\texttt{old}}  $ (but not $ \widetilde R_{\texttt{old}}$ or $ R_{\texttt{old}}^\top$). The large outer matrices only record the span of the left and right subspaces and are built by appending columns to $Q_{\texttt{old}} $ and rows to $R_{\texttt{old}} $.  The updates only involve with the small matrices $\widetilde Q_{\texttt{old}} $ and $\widetilde R_{\texttt{old}}  $. This makes the update much faster.

However, over thousands or millions of updates, the multiplications may erode the orthogonality of  through numerical error; see the right figure in <Ref>. Hence reorthogonalizations are needed for these small orthonormal matrices. In [4], Brand asked: “It is an open question how often this is necessary to guarantee a certain overall level of numerical precision; it does not change the overall complexity.”  We also note that only reorthogonalize the small matrix $\widetilde Q$ is not enough for some data sets. In <Ref>, we show that without reorthogonalizations for the larger outside matrix $Q_{\texttt{old}}$, the algorithm provides untruthful SVD.

In this paper,  we first answer the open question,  and the answer is we can avoid computing the large amount of those  orthogonal matrices and hence the  reorthogonalizations for those small orthonormal matrices are not necessary by modifying Brand's algorithm.  Second, we proposed an efficient reorthogonalization for the larger out matrix $Q_{\texttt{old}} $, numerical experiments in <Ref> show that the reorthogonalization works well in practice. Third, we prove that our modification does not change the outputs of Brand's algorithm.

§ INCREMENTAL SVD (I)

We begin by introducing notation needed throughout the paper.  Let $I_k$ be a $k\times k$ identity matrix and the $W$-weighted inner product:

\begin{align*}
(a, b)_W = a^\top W b,\quad a, b\in \mathbb R^m.
\end{align*}

Let $U $ be a $m$ by $n$ dense matrix  and $u_j$ denote the $j$-th column of $U$, i.e.,

\begin{align*}
U = \left[u_1 ~ \big|~ u_2  ~ \big|~\cdots~ \big|~u_n  \right].
\end{align*}

For convenience, we adopt Matlab notation herein. We use both $U(:, 1:\ell)$ and $U_\ell$ to denote the first $\ell$ columns of $U$, and $U(1:\ell,1:\ell)$ be the $\ell-$th leading principal minor of $U$. The function : $\mathbb{R}^{k} \rightarrow \mathbb{R}^{k \times k}$ takes as input a vector and outputs a square, diagonal matrix with that vector's entries on its main diagonal as follows: if $\alpha =\left(a_{1}, a_{2}, \ldots, a_{k}\right)^\top$, then

\texttt{diag}(\alpha)=\left[\begin{array}{llll}
a_{1} & & & \\
& a_{2} & & \\
& & \ddots & \\
& & & a_{k}
\end{array}\right] .

If the data arising from a Galerkin-type simulation of a PDE, one has to deal with Cholesky factorization of  a weighted matrix if one directly apply Brand's algorithm.  To avoid this,
Fareed et al. [7] extended Brand's algorithm to accommodate data of this type without computing Cholesky factorization.  More specifically,  if the data lies in a finite element space and is expressed using a collection of basis functions. Then  the SVD of the data is equivalent to find a so-called core SVD of their coefficient matrix.

A core SVD of a matrix $U\in \mathbb R^{m\times n}$ is a decomposition $U=Q \Sigma R^{\top}$, where $Q \in \mathbb{R}^{m \times d}, \Sigma \in \mathbb{R}^{d \times d}$, and $R \in \mathbb{R}^{n \times d}$ satisfy

Q^{\top} W Q=I, \quad R^{\top} R=I, \quad \Sigma=\textup{\texttt{diag}}\left(\sigma_{1}, \ldots, \sigma_{d}\right),

where $\sigma_{1} \geq \sigma_{2} \geq \cdots \geq \sigma_{d}>0$. The values $\left\{\sigma_{i}\right\}$ are called the (positive) singular values of $U$ and the columns of $Q$ and $R$ are called the corresponding singular vectors of $U$.

We note that if $W = I$, the above core SVD is reduced to the standard SVD of the matrix $U$. In [7], $W$ is the standard finite element mass or stiffness matrix.  For more information about the incremental SVD in a weighted norm setting, see [6, 5].

Next, we introduce Brand's incremental  SVD algorithm and follow the discussions in [7].

Step 1: Initialization. Assume that the first column of $U$  is nonzero, i.e., $u_1\neq 0$, we initialize the core SVD of $u_1$ by setting

\begin{align*}
\Sigma=(u_1^\top W u_1)^{1/2}, \quad Q=u_1 \Sigma^{-1}, \quad R=1.
\end{align*}

The algorithm is shown in <Ref>.
(Initialize incremental SVD)

Input: $ u_1\in \mathbb R^m$, $W\in \mathbb R^{m\times m}$
[1]
$\Sigma = \left(u_1^\top W u_1\right)^{1/2}$;  $Q = u_1\Sigma^{-1} $;  $R=1$;
return $Q,\Sigma,R$

Step 2: Core SVD of $U_\ell$. Suppose we already have a rank- $k$ truncated core SVD of $U_\ell$:

\begin{align}\label{SVD01}
U_\ell = Q \Sigma R^{\top},\quad \textup{with}\;\; 	Q^{\top} W Q=I_k, \quad R^{\top} R=I_k, \quad \Sigma=\texttt{diag}\left(\sigma_{1}, \ldots, \sigma_{k}\right),
\end{align}

where $\Sigma \in \mathbb{R}^{k \times k}$ is a diagonal matrix with the $k$ (ordered) singular values of $U_\ell$ on the diagonal, $Q \in \mathbb{R}^{m \times k}$ is the matrix containing the corresponding $k$ left singular vectors of $U_\ell$, and $R \in \mathbb{R}^{\ell \times k}$ is the matrix of the corresponding $k$ right singular vectors of $U_\ell$.

Step 3: Updates the core SVD of $U_{\ell+1}$. Our goal next is to update the above SVD, i.e., we want to find the core SVD of $U_{\ell+1}$ by using the core SVD of $U_\ell$ and the new adding vector  $u_{\ell+1}$.

To do this, let $e \in \mathbb{R}^{m}$ be the residual  of $u_{\ell+1}$ onto the  space spanned by the columns of $Q$ with the weighted inner product $(\cdot, \cdot)_W$. Therefore,

\begin{align*}
e =    u_{\ell+1}-QQ^{\top}W u_{\ell+1}.
\end{align*}

Let $p$ be the magnitude of $e$, i.e., $p=(e^\top W e)^{1/2}$. If $p>0$, we let $\widetilde e$ be the unit vector in the direction of $e$, i.e., $\widetilde e=e / p$. If $p=0$, we set $\widetilde e=0$. Then we have the fundamental identity:

\begin{align}\label{fundamental_increment_SVD}
\begin{split}
U_{\ell+1} =
\left[
U_\ell ~ \big|~ u_{\ell+1}\right]=\left[
Q\Sigma R^{\top} ~ \big|~ u_{\ell+1}\right]=\left[
Q ~ \big|~  \widetilde e \right]\underbrace{\left[\begin{array}{cc}
\Sigma & Q^{\top}W u_{\ell+1} \\
0 & p
\end{array}\right]}_{Y}\left[\begin{array}{cc}
R & 0 \\
0 & 1
\end{array}\right]^{\top}.
\end{split}
\end{align}

We can find the SVD of the updated matrix $U_{\ell+1}$ by finding the full SVD of the  matrix $Y$ in the right hand side of the above identity.

(1) If $p<\texttt{tol}$, we approximate and set $p=0$ and $\widetilde e = 0$.  Let $Q_Y \Sigma_Y R_Y^\top$ be the full SVD of $
Y = \left[\begin{array}{cc}
\Sigma & Q^{\top}W u_{\ell+1} \\
0 & 0
\end{array}\right]$.
Then the  core SVD of $U_{\ell+1}\in \mathbb R^{m\times (\ell+1)}$ is given by

\left[
U_\ell ~ \big|~	u_{\ell+1}\right]=\left(\left[
Q ~ \big|~0\right] Q_Y\right) \Sigma_{Y}\left(\left[\begin{array}{cc}
R & 0 \\
0 & 1
\end{array}\right] R_{Y}\right)^{\top}.

This suggests the following update

\begin{align}\label{ranksame}
Q\longleftarrow  QQ_Y(1:k,1:k), \quad  \Sigma \longleftarrow  \Sigma_Y(1:k,1:k), \quad  R \longleftarrow \left[\begin{array}{cc}
R & 0 \\
0 & 1
\end{array}\right] R_Y(:,1:k).
\end{align}

(2) If $p\ge \texttt{tol}$, we let $Q_Y \Sigma_Y R_Y^\top$ be the full SVD of the middle matrix $Y$ in (<ref>), and  then update the SVD of $U_{\ell+1}$ by

\begin{align}\label{rankaddone}
Q\longleftarrow \left[
Q~ \big|~  \widetilde e
\right]Q_Y\quad   \Sigma \longleftarrow  \Sigma_Y\quad   R \longleftarrow  \left[\begin{array}{cc}
R & 0 \\
0 & 1
\end{array}\right] R_Y.
\end{align}

In [7], the authors considered a  data from the simulation of 2D laminar flow around a cylinder with circular cross-section.  The  incremental SVD algorithms kept computing very small singular values (near machine precision) - these required huge computational cost and these values were not needed in their application.   Hence, they  suggested another truncation in <cit.> if the last few singular values are less than a  tolerance. That is, find  the minimal $r$ such that $\Sigma_Y(r,r)\ge \texttt{tol}$ and $\Sigma_Y(r+1,r+1)\ge \texttt{tol}$  and set the following truncation:

\begin{align}\label{Updatesnohope1}
Q\longleftarrow Q(:,1:r),\quad   \Sigma \longleftarrow  \Sigma(1:r,1:r), \quad   R \longleftarrow  R(:,1:r),
\end{align}

where $Q, \Sigma$ and $R$ were defined in (<ref>). Now we summarize the above discussions  in <Ref>.

(Update Incremental SVD (I))

Input: $Q\in \mathbb R^{m\times k}$, $\Sigma\in \mathbb R^{k\times k}$, $R\in \mathbb R^{\ell\times k}, u_{\ell+1}\in \mathbb R^m$, $W\in \mathbb R^{m\times m}$, $\texttt{tol}$
[1]
Set $d = Q^\top (W u_{\ell+1})$;   $e = u_{\ell+1} - Qd$;   $p = (e^\top W e)^{1/2}$;
Set $p=0$;
Set  $  e = e/p$;
Set $Y = \left[\begin{array}{cc}
\Sigma & d\\
0 & p
\end{array}\right]$;

$[Q_Y, \Sigma_Y, R_Y] = \texttt{svd}(Y)$;
Set $Q = QQ_Y(1:k, 1:k)$,   $\Sigma = \Sigma_Y(1:k,1:k)$,    $R =  \left[\begin{array}{cc}
R & 0 \\
0 & 1
\end{array}\right] R_Y(:,1:k)$;
Set $
Q= \left[
Q~ \big|~   e
\right]Q_Y, \quad   \Sigma = \Sigma_Y, \quad  R =  \left[\begin{array}{cc}
R & 0 \\
0 & 1
\end{array}\right] R_Y$.
$\Sigma(r,r)\ge \texttt{tol}$&& $\Sigma(r+1,r+1)\ge \texttt{tol}$
Set $Q = Q(:,1:r)$,     $\Sigma = \Sigma(1:r,1:r)$,    $R = R(:,1:r)$;

return $Q, \Sigma, R$

Next, we discuss the complexity of Brand's algorithm (<Ref>) with $W = I$ and  note that $U\in \mathbb R^{m\times n}$ with low rank $r$. The total time complexity is $\mathcal O((m+n)r^2)$ and the details are shown below.

1. Line 1: this needs to be executed $n$ times and the computational cost of each time is about $\mathcal O(mr)$, the whole complexity is $\mathcal O(mnr)$.

2. Lines 8:  since $Y$ is a one-column bordered diagonal matrix, then it can be bidiagonalized in $\mathcal O(r^2)$ time and  the whole time is $\mathcal O(nr^2)$.

3. Lines 10-12:  the computational cost of each time is about $\mathcal O(mr^2 + \ell r^2)$, hence the total complexity is about $\mathcal O\left(\sum_{\ell=1}^n (mr^2 + \ell r^2)\right) = \mathcal O((m+n)nr^2)$.

Step 4: Reorthogonalization. Theoretically, the above  SVD update yields orthonormal left and right singular
vectors. However, in practice, small numerical errors cause a loss of orthogonality. Then the incremental SVD algorithm is a non-orthogonal decomposition, and the singular values are not true hence truncation is not reliable. Therefore a reorthogonalization is needed and the computational cost can be high if the rank of the matrix is not extremely low. Specifically, Brand recommended a Gram-Schmidt orthogonalization [3] while Oxberry et al. chose thin QR  in [13]. Our numerical test show that thin QR is  faster than the Gram-Schmidt when the weighted matrix is an identity, see <Ref>.  However, since there is no thin QR with a weighted norm setting at hand, Fareed et al. in [7] use the $W$-weighted Gram–Schmidt for the  reorthogonalization.

Input: $Q\in \mathbb R^{m\times k}$, $W\in \mathbb R^{m\times m}$,  $\texttt{tol}$
[1]
$| (Q(:, \texttt{end}))^\top W Q(:,1))|> \texttt{tol}$
$i=1$ to $k$
Set $\alpha = Q(:, i)$;
$j=1$ to $i-1$
$Q(:,i) = Q(:,i) - (\alpha^\top W Q(:,j))Q(:,j)$;
Set  $= ({(Q(:,i))^\top W Q(:,i)})^{1/2}$;
Set $Q(:,i) =Q(:,i)/\texttt{norm}$;
return $Q$

We do not discuss the time complexity of the reorthogonalization since there is no clue that how many times we need to execute the <Ref>. However, numerical experiments in <Ref> show that the orthogonalization step described above is a large part of the computational cost of the incremental SVD algorithm.

Finally we conclude the full implementation of the incremental SVD with the $W$-weighted Gram–Schmidt for the  reorthogonalization in <Ref>.
(Fully incremental SVD (I))

Input:  $W\in \mathbb R^{m\times m}$, $\texttt{tol}$
[1]
Get $u_1$;
$[Q, \Sigma, R]=$  ($u_1$, $W$);   % <Ref>
$\ell=2,\ldots, n$
Get $u_\ell$;
$[Q, \Sigma, R]=$  ($Q, \Sigma, R, u_{\ell}$, $W$, $\texttt{tol}$);  % <Ref>
$Q =$  ($Q$, $W$, $\texttt{tol}$);  % <Ref>
return $Q, \Sigma, R$

In this example, we test the <Ref>. We let $\Omega \subset \mathbb R^2$ be an unit square and we partition  it into  $524288$ uniform triangles.  Let $\{\varphi_i\}_{i=1}^m$ be the standard linear finite element basis functions. Next we  let $\{(x_i, y_j)\}_{i,j=1}^m$ be the grids of $\Omega$,  $\{t_i\}_{i=1}^{n}$ be an equally space grid in $[0, 10]$ and  time step $\Delta t = 1/10^3$.  Define

\begin{gather*}
f(t,x, y) = \cos(t(x+y)),\quad b_i = [(f(t_i),\varphi_j)]_{j=1}^m, \quad B = [b_1\mid b_2\mid \ldots\mid b_n],\\
M = [(\varphi_j, \varphi_i)]_{i,j=1}^m,\quad u_k = [f(t_k, x_i, y_j)]_{i,j=1}^m,\quad  U= [u_1\mid u_2\mid \ldots\mid u_n].
\end{gather*}

Here $M$ is the finite element mass matrix and
$u_k$ is the coefficient of the linear Lagrange interpolation of $f(t,x,y)$ at $t = t_k$.  We  compute  the SVD of the matrix $B$ with $W = I$ and of the matrix $U$ with $W = M$ by setting $\texttt{tol}=10^{-12}$. Let  $\|A\|_F$ be  the  Frobenius norm of matrix $A$, and we use

\begin{align*}
\mathcal E_W = \|I - Q^\top W Q\|
\end{align*}

to measure the error of the  orthogonality of the orthonormal matrix $Q$ under the weight inner product $(\cdot,\cdot)_W$.

In <Ref> we show the orthogonality of the left singular matrix $Q$ with different weighted norms. For the case  $W = M$, the orthogonality of $Q$ is eroded (even with reorthogonalizations).

<Ref>: Left:  $W = I$. Right: $W = M$.

The whole simulations take about 154 and  4000 seconds for  $W = I$ and $W = M$, respectively. 	Next, we show the CPU[All the code for all examples in the paper has been made by the author using MATLAB R2020b and has been run on a laptop with MacBook Pro, 2.3 Ghz8-Core Intel Core i9 with 64GB 2667 Mhz DDR4.] time (seconds)  of the main part of the <Ref> in <Ref>. We see that the  orthogonalization step  is a large part of the computational cost of the incremental SVD algorithm, especially for the  case $W = M$.

Line 1 in <Ref>      Lines 9-13 in <Ref>     <Ref>

$W = I$
32.034    	60.213       60.129

$W = M$
110.58    	115.96       3761.6

The CPU time  (seconds)  for different weighted inner products.

§ INCREMENTAL SVD (II)

<Ref> is easy to implement while the computational cost is high, especially for the weighted norm setting. For a $m\times n$ dense matrix with low rank $r$, the time  complexity of <Ref>  is $\mathcal O(m(m+n)r^2)$.  In [4], Brand reduced the complexity to $\mathcal O(mnr)$. However, to update the left singular space, one has to  compute  thousands or millions of orthogonal matrices and multiply them together. It is more complicate to update the right singular space  which need to compute the pseudo-inverse of a matrix  for thousands or millions times.

Next, we follow [4] to extend  Brand's improvements to compute the core SVD of $U\in \mathbb R^{m\times n}$ with respect to
the weighted inner product $(\cdot, \cdot)_W$.  Instead of updating the large singular vector matrices as prescribed in (<ref>) or (<ref>), we
suppose the core SVD  of $U_\ell\in \mathbb R^{m\times \ell}$  can be written into the  following five matrices

\begin{align}\label{Newformal}
U_{\ell} = Q_{\texttt{old}}  \widetilde Q_{\texttt{old}}  \Sigma_{\texttt{old}}  \widetilde R_{\texttt{old}}^\top  R_{\texttt{old}}^\top.
\end{align}

Here $Q_{\texttt{old}}\in \mathbb R^{m\times k}$ is a
orthonormal matrix under the weighted inner product,  $ \widetilde Q_{\texttt{old}}$ and $\widetilde R_{\texttt{old}}^\top  R_{\texttt{old}}^\top  $ are standard orthonormal matrices, but $\widetilde R_{\texttt{old}}^\top $ and $ R_{\texttt{old}}^\top  $ are not. The large outer matrices only record the span of the left and right subspaces and are built by appending columns to $Q_{\texttt{old}}$ and rows to $R_{\texttt{old}}^\top$. And the update only relates to the small matrices  $ \widetilde Q_{\texttt{old}}$ and $\widetilde R_{\texttt{old}}^\top $, this makes the update much faster.

Next, we consider the core SVD of $U_{\ell+1} = \left[U_\ell ~\big|~ u_{\ell+1}\right]\in \mathbb R^{m\times(\ell+1)}$. Similar to  <Ref>, we let $e \in \mathbb{R}^{m}$ be the residual  of $u_{\ell+1}$ onto the  space spanned by the columns of $Q_{\texttt{old}}\widetilde Q_{\texttt{old}}$ with the weighted inner product $(\cdot, \cdot)_W$, i.e.,

\begin{align*}
e =    u_{\ell+1}-( Q_{\texttt{old}}   \widetilde Q_{\texttt{old}})( Q_{\texttt{old}}   \widetilde Q_{\texttt{old}})^\top W u_{\ell+1} = u_{\ell+1}- Q_{\texttt{old}}  Q_{\texttt{old}}^\top W u_{\ell+1}.
\end{align*}

Let $p=(e^\top W e)^{1/2}$. If $p>0$, we let  $\widetilde e=e / p$, otherwise we set $\widetilde e=0$. Then by (<ref>) and  the fundamental identity (<ref>) we have

\begin{align*}
U_{\ell+1} =\left[
Q_{\texttt{old}}  \widetilde Q_{\texttt{old}} ~ \big|~  \widetilde e \right]\underbrace{\left[\begin{array}{cc}
\Sigma_{\texttt{old}} & \widetilde Q_{\texttt{old}} Q_{\texttt{old}}^{\top}W u_{\ell+1} \\
0 & p
\end{array}\right]}_{Y}\left[\begin{array}{cc}
R_{\texttt{old}}  \widetilde R_{\texttt{old}} & 0 \\
0 & 1
\end{array}\right]^{\top}.
\end{align*}

Assume that  $Q_Y \Sigma_Y R_Y^\top$ be the full SVD of $Y$, then the core SVD of  $U_{\ell+1}$ is given by

\begin{align}\label{NeedR}
Q_{\texttt{old}}  \widetilde Q_{\texttt{old}} ~ \big|~\widetilde e\right] Q_Y\right) \Sigma_{Y}\left(\left[\begin{array}{cc}
R_{\texttt{old}}  \widetilde R_{\texttt{old}} & 0 \\
0 & 1
\end{array}\right]R_{Y}\right)^{\top}.
\end{align}

(1) If $p<\texttt{tol}$, we set $p=0$ and $\widetilde e = 0$,  then the update is not rank-increasing. We update the small orthonormal matrix $\widetilde Q_{\texttt{new}}\in\mathbb R^{m\times k}$ and the diagonal matrix $\Sigma_{\texttt{new}}\in \mathbb R^{k\times k}$  by

\begin{align}\label{ranksame1}
\widetilde Q_{\texttt{new}} \longleftarrow  	\widetilde Q_{\texttt{old}} Q_Y(1:k,1:k), \qquad  \Sigma_{\texttt{new}} \longleftarrow  \Sigma_Y(1:k,1:k).
\end{align}

(2) If $p\ge \texttt{tol}$, we  then update the large outer matrix $Q_{\texttt{new}} $ by adding the vector $\widetilde e$ to $Q_{\texttt{old}} $, and update the small orthonormal matrix $\widetilde Q_{\texttt{new}}\in\mathbb R^{m\times (k+1)}$ and the diagonal matrix $\Sigma_{\texttt{new}}\in \mathbb R^{(k+1)\times (k+1)}$  by

\begin{align}\label{rankaddone1}
Q_{\texttt{new}} \longleftarrow \left[Q_{\texttt{old}}~\big|~ \widetilde e\right],\qquad  \widetilde Q_{\texttt{new}} \longleftarrow  	\left[\begin{array}{cc}
\widetilde Q_{\texttt{old}} & 0 \\
0 & 1
\end{array}\right]  Q_Y, \qquad  \Sigma_{\texttt{new}} \longleftarrow  \Sigma_Y.
\end{align}

The time complexity of (<ref>) or (<ref>) is about $\mathcal O(r^2)$ since the update of the left singular vector matrix is completed by the multiplication of two small orthnormal matrices. This is a major difference   comparing to   (<ref>)  or (<ref>).

However, we have to update the small orthonormal matrices over thousands or millions of times, the multiplications may erode the orthogonality of $\widetilde Q$ through numerical error. Hence reorthogonalizations are  needed when the  inner product between the first and last left singular vectors of $\widetilde Q$ is greater than some tolerance.

In [4], Brand asked “It is an open question how often this is necessary to guarantee a certain overall level of numerical precision; it does not change the overall complexity”. In the next section, we answer this question and the answer is we can avoid computing the large amount of those  orthogonal matrices and those  multiplications are not necessary. Hence, the reorthogonalizations are not needed.

Theoretically, the new adding vector $\widetilde e$ is orthogonal to each column of $Q_{\texttt{old}}$. However, in practice,  the matrix  $\left[Q_{\texttt{old}}~\big|~ \widetilde e\right]$ loss its orthogonality very  quickly for some data sets. Hence reorthogonalizations are needed for this matrix. We note that in [4], Brand did not include this. In <Ref>, we show how this is important and in the next section we propose a very efficient way to remedy this issue.

The updates of right singular vector matrix $R_{\texttt{new}}$ and $\widetilde R_{\texttt{new}}$  are  more complicated because of  adding rows to $R_{\texttt{old}}$ while guaranteeing that the columns of the product $R_{\texttt{new}}\widetilde R_{\texttt{new}}$ are orthonormal. From the fundamental identity (<ref>) and (<ref>), the right-side update must satisfy

\begin{align*}
R_{\texttt{new}}\widetilde R_{\texttt{new}} = 	\left[\begin{array}{cc}
R_{\texttt{old}}\widetilde R_{\texttt{old}}& 0 \\
0 & 1
\end{array}\right] R_Y.
\end{align*}

In the following updates, we need to compute the pseudo-inverse of $\widetilde R_{\texttt{old}}$ and a sub-matrix  of $R_Y$. We use $X^+$ to denote the  pseudo-inverse of $X$.

(3) If $p< \texttt{tol}$, the rank does not increase, the last column of $R_Y$ represents an unused subspace dimension and should be suppressed. Let $\mathcal R_1 =R_Y(1:k,1:k)$, $\mathcal R_2 =R_Y(k+1,1:k)$,

\begin{align}\label{rankaddone2}
\widetilde R_{\texttt{new}} \longleftarrow  	\widetilde R_{\texttt{old}}\mathcal R_1 , \quad 	\widetilde R_{\texttt{new}}^+ \longleftarrow
\mathcal R_1^+ 	\widetilde R_{\texttt{old}}^+, \quad R_{\texttt{new}} \longleftarrow  	\left[\begin{array}{c}
R_{\texttt{old}} \\
\mathcal R_2\widetilde R_{\texttt{new}}^+
\end{array}\right].
\end{align}

(4) If $p\ge \texttt{tol}$,   the update is  rank-increasing. We update the small orthonormal matrix $\widetilde Q_{\texttt{new}}\in\mathbb R^{m\times k}$ and the diagonal matrix $\Sigma_{\texttt{new}}\in \mathbb R^{k\times k}$  by

\begin{align}\label{ranksame2}
\widetilde R_{\texttt{new}} \longleftarrow  	\left[\begin{array}{cc}
\widetilde R_{\texttt{old}}& 0 \\
0 & 1
\end{array}\right] R_Y, \quad 	\widetilde R_{\texttt{new}}^{+} \longleftarrow  	R_Y^\top \left[\begin{array}{cc}
\widetilde R_{\texttt{old}}^+& 0 \\
0 & 1
\end{array}\right],
\quad  R_{\texttt{new}} \longleftarrow  	\left[\begin{array}{cc}
R_{\texttt{old}}& 0 \\
0 & 1
\end{array}\right].
\end{align}

The implementation of this incremental SVD algorithm is shown in <Ref>.
(Incremental SVD (II))

Input: $Q\in \mathbb R^{m\times k}$, $\Sigma\in \mathbb R^{k\times k}$, $R\in \mathbb R^{\ell\times k}, u_{\ell+1}\in \mathbb R^m$, $W\in \mathbb R^{m\times m}$, $\texttt{tol}$, $\widetilde Q$, $\widetilde R$, $\widetilde R^+$
[1]
Set $d = Q^\top (W u_{\ell+1})$; $e = u_{\ell+1} - Qd$;   $p = (e^\top W e)^{1/2}$;
Set $Y = \left[\begin{array}{cc}
\Sigma & d\\
0 & 0
\end{array}\right]$;
$[Q_Y, \Sigma_Y, R_Y] =  \texttt{svd}(Y)$;
Set $\widetilde Q = \widetilde Q Q_Y$,   $\Sigma = \Sigma_Y$,    $\mathcal  R_1= R_Y(1:k,1:k)$,    $\mathcal  R_2= R_Y(k+1,1:k)$,

$\widetilde R= \widetilde R \mathcal  R_1$,   $\widetilde R^+ =  \mathcal  R_1^+ \widetilde R^+$,   $R=  \left[\begin{array}{c}
R \\
\mathcal R_2 \widetilde R^+
\end{array}\right]$;
Set  $  e = e/p$;

Set $Y = \left[\begin{array}{cc}
\Sigma & d\\
0 & p
\end{array}\right]$;

$[Q_Y, \Sigma_Y, R_Y] = $  $(Y)$;
Set $
\widetilde Q=  \left[\begin{array}{cc}
\widetilde Q & 0 \\
0 & 1
\end{array}\right]Q_Y, \;  Q = [Q ~\big|~ e],\;  \Sigma = \Sigma_Y, \;  R =  \left[\begin{array}{cc}
R & 0 \\
0 & 1
\end{array}\right]$,

$\widetilde R =  \left[\begin{array}{cc}
\widetilde R & 0 \\
0 & 1
\end{array}\right]R_Y$,     $\widetilde R^+ = R_Y^\top \left[\begin{array}{cc}
\widetilde R^+ & 0 \\
0 & 1
\end{array}\right]$;
$|(\widetilde Q(:,1))^\top \cdot \widetilde Q(:,\texttt{end})|>\texttt{tol}$
Apply Gram-Schmidt to Reorthogonalize the matrix $\widetilde Q$;
return $Q$, $\Sigma$, $R$, $\widetilde Q$, $\widetilde R$,  $\widetilde R^+$

(Fully incremental SVD (II))

Input:  $W\in \mathbb R^{m\times m}$, $\texttt{tol}$
[1]
Get $u_1$;
$[Q, \Sigma, R]$=($u_1$, $W$);
% <Ref>
Set $\widetilde Q = \widetilde R = 1$;
$\ell=2,\ldots, n$
Get $u_\ell$
$[ Q, \Sigma, R,\widetilde Q, \widetilde R, \widetilde R^+]$=($Q, \Sigma, R, \widetilde Q, \widetilde R, \widetilde R^+, u_{\ell},W, \texttt{tol}$);    % <Ref>
return $Q$, $\Sigma$, $R$, $\widetilde Q$, $\widetilde R$

In this example, we first  show that the reorthogonalizations are needed for the matrix $Q$ in <Ref>.   Let $\Omega \subset \mathbb R^2$ be an unit square and  partition  it into  $512$ uniform triangles.   Let $\{\varphi_i\}_{i=1}^m$ be the standard linear finite element basis functions, $\{t_i\}_{i=1}^{n}$ be an equally space grid in $[0, 10]$ and  time step $\Delta t = 1/10^2$.  Define

\begin{gather*}
f(t,x, y) = \cos(t(x+y)),\quad b_i = [(f(t_i),\varphi_j)]_{j=1}^m, \quad B = [b_1\mid b_2\mid \ldots\mid b_n].
\end{gather*}

We then use the <Ref> to compute the SVD of $B$ with the weighted matrix $W = I$. To make a comparison, we use the MATLAB built-in function $\texttt{svd()}$ to find the “exact” singular values of $B$. From <Ref> we see that Brand's algorithm  provides  unrelated  singular values.

Singular values of $B$: Left: <Ref>.  Right:   built-in function $\textup{\texttt{svd()}}$.

Next, we add  reorthogonalizations (details will be given in the next section) to <Ref> and plot the errors (left in <Ref>) between the modified <Ref> and the MATLAB built-in function $\texttt{svd()}$. We see that the new adding  reorthogonalizations works well in practice.

Finally, we test the orthogonality of the small orthonormal matrices $\widetilde Q$ at each update,  we use $\mathcal E_I(\widetilde Q) = \|I - \widetilde Q^\top \widetilde Q\|$ to measure the orthogonality of $\widetilde Q$, which is shown  in the right of <Ref>,  we see that the matrix $\widetilde Q$ loss its orthogonality quickly after $n=600$ even performed the reorthogonalizations.

<Ref>: Left: the errors of the first 34 singular values  of $B$ by adding the reorthogonalizations. Right:   the error of the orthogonality of $\widetilde Q$ at each update.

§ OUR IMPROVEMENTS

As it is pointed  in <Ref>, the <Ref> needs to update the small orthonormal matrices $\widetilde Q$ over thousands or millions of times, the multiplications may erode the orthogonality of $\widetilde Q$ through numerical error. To overcome this issue,  Brand proposed reorthogonalizations when the  inner product between the first and last left singular vectors of $\widetilde Q$ is greater than some tolerance.  However, as we have seen in <Ref>, even we performed the reorthogonalizations, the orthonormal matrix $\widetilde Q$ could still lose its orthogonality.  Furthermore,  to update the right singular space, the <Ref> needs to compute the pseudo-inverse of a matrix for thousands or millions of times.

First, we summarize our contributions in this section.

(1) We  answer the question that Brand asked in [4]: “It is an open question how often this is necessary to guarantee a certain overall level of numerical precision; it does not change the overall complexity”.

(2) We prove that the output of our new  algorithm is the same with  Brand's.

(3) We show that it is not necessary to compute  thousands or millions of the pseudo-inverse of a matrix to update the right singular subspace.

(4) We propose a  reorthogonalization for the matrix $\left[Q_{\texttt{old}}~\big|~ \widetilde e\right]$, we note that this step is not in the original algorithm in [4]. In <Ref> we have showed that the reorthogonalizations are necessary.

To begin,  we introduce the following lemma:

Let $A\in \mathbb R^{k\times k}$, $\alpha\in \mathbb R^k$, and let $Q_1\Sigma_1R_1^\top$ be a thin SVD of $\left[A ~\big|~ \alpha\right]$, then

\begin{align}\label{SVDnew}
\left[\begin{array}{cc}
Q_1 & 0 \\
0 & 1
\end{array}\right] \left[\begin{array}{cc}
\Sigma_1 & 0 \\
0 & 0
\end{array}\right]\left[\begin{array}{c}
R_1^\top \\
\end{array}\right]
\end{align}

is a SVD of $	\left[\begin{array}{cc}
A & \alpha\\
0 & 0
\end{array}\right]$, where $r_{k+1}\in \mathbb R^{k+1}$ is an unit vector and orthogonal to each column of $R_1$.

Next, we use <Ref> to  modify the updates of the left singular subspace if the rank does not increase. Assume that  $Q \Sigma R^\top$ is a  rank $k$ truncated core SVD of $U_\ell$,   the norm  of $u_{\ell+1}- QQ^\top W u_{\ell+1}$ is very small and we approximate it by $0$.  Therefore,

\begin{align}\label{fundamental_increment_SVD10}
\begin{split}
U_{\ell+1} = \left[U_\ell ~ \big | ~ u_{\ell+1} \right] =
Q  \left[	\Sigma  ~ \big | ~  Q^{\top}W u_{\ell+1}\right] \left[\begin{array}{cc}
R & 0 \\
0 & 1
\end{array}\right]^{\top}.
\end{split}
\end{align}

Let $Q_1\Sigma_1R_1^\top$ be a thin SVD of $ \left[	\Sigma  ~ \big | ~  Q^{\top}W u_{\ell+1}\right]$, then the core SVD of $U_{\ell+1}$ is given by

\begin{align}
U_{\ell+1}  = QQ_1 \Sigma_1 R_1^\top  \left[\begin{array}{cc}
R & 0 \\
0 & 1
\end{array}\right]^{\top}.
\end{align}

Hence the following update

\begin{align}\label{SVD}
Q \longleftarrow  Q Q_1,\qquad \Sigma  =  \Sigma_1, \qquad  R =\left[\begin{array}{cc}
R& 0 \\
0 & 1
\end{array}\right]  R_1,
\end{align}

is exactly the same with (<ref>) due to <Ref>.

Next, we consider a more general situation.   Assume that the matrix $C$ has $s$ columns, and the norm of $C - QQ^\top W C$ is zero or very small, then we approximate it by $0$ and  we then have the following fundamental identity:

\begin{align}\label{fundamental_increment_SVD2}
\begin{split}
U_{\ell+1} = \left[U_\ell ~ \big | ~ C \right] =
Q \underbrace{\left[	\Sigma  ~ \big | ~  Q^{\top}W C\right]}_{Y} \left[\begin{array}{cc}
R & 0 \\
0 & I_s
\end{array}\right]^{\top}.
\end{split}
\end{align}

In this paper we use the MATLAB built-in function $\texttt{svd(Y, {\color{magenta} `econ'})}$ to find the SVD of $Y$ since the matrix $Y$ is always short-and-fat.

§.§ Update the left and right singular spaces

Since we  assume that the matrix $U$ is low rank, then it is reasonable to anticipate that most vectors of $\{u_{\ell+1}, u_{\ell+2},\ldots, u_n\}$ are linear dependent or almost linear dependent with  the vectors in $Q\in \mathbb R^{m\times k}$. Without loss of generality,  we assume that the next  $2$ vectors, $u_{\ell+1}$ and $u_{\ell+2}$,  their residuals are less than the $\texttt{tol}$ in <Ref> when project them onto the  space spanned by the columns of $Q$.  In other words, we assume that

\begin{align}\label{assss}
\|u_{\ell+1}- QQ^\top Wu_{\ell+1}\|_W <\texttt{tol}, \qquad \|u_{\ell+2}- QQ^\top Wu_{\ell+2}\|_W <\texttt{tol}.
\end{align}

We then follow  Brand's algorithm to find the SVD of $U_{\ell+2}$ incrementally:

(1) Find the SVD of $\left[ \Sigma   ~\big|~  Q^{\top} W u_{\ell+1} \right]$, and we assume

\begin{align}\label{SVD001}
\left[ \Sigma   ~\big|~  Q^{\top} W u_{\ell+1} \right] = Q_{(1)}\Sigma_{(1)} R_{(1)}^\top,
\end{align}

where $Q_{(1)}\in \mathbb R^{k\times k}$, $\Sigma_{(1)}\in \mathbb R^{k\times k}$, $R_{(1)}\in \mathbb R^{(k+1)\times k}$ and   $Q_{(1)}^\top Q_{(1)}  = Q_{(1)} Q_{(1)} ^\top = I_k $.

(2) Update the SVD of $U_{\ell+1}$ by

\begin{align}\label{SVD002}
Q_{\ell+1}  =   Q Q_{(1)},\qquad \Sigma_{\ell+1}  =  \Sigma_{(1)},\qquad  R_{\ell+1} =\left[\begin{array}{cc}
R& 0 \\
0 & 1
\end{array}\right]  R_{(1)}.
\end{align}

By the above discussion, the update (<ref>) is exactly the same with (<ref>). Next, according to <Ref> , we should check the magnitude of the residue of $u_{\ell+2}$ when project onto the  space spanned by the columns of $Q_{\ell+1}$. To do this, we introduce the following lemma.

Let $\mathtt{Q} \in \mathbb R^{k\times k}$ be an orthonormal matrix, then for any $x\in \mathbb R^m$, $Q\in \mathbb R^{m\times k}$ we have

\begin{align*}
QQ^\top x = (Q\mathtt{Q}) (Q\mathtt{Q})^\top x.
\end{align*}

Since $Q_{(1)}$ is an orthonormal matrix, then by <Ref> and (<ref>) we have

\begin{align*}
\|u_{\ell+2}- Q_{\ell+1}Q_{\ell+1}^\top Wu_{\ell+2}\|_W  = \|u_{\ell+2}- QQ^\top Wu_{\ell+2}\|_W <\texttt{tol}.
\end{align*}

Next we follow the same arguments as in the above steps (1) and (2) to deal with $u_{\ell+2}$.

(3) Find the SVD of $\left[ \Sigma_{\ell+1}  ~\big|~  Q_{\ell+1}^{\top} W u_{\ell+2} \right]$, and we assume

\begin{align}\label{SVD003}
\left[ \Sigma_{\ell+1}  ~\big|~  Q_{\ell+1}^{\top} W u_{\ell+2} \right] = Q_{(2)}\Sigma_{(2)} R_{(2)}^\top,
\end{align}

where $Q_{(2)}\in \mathbb R^{k\times k}$, $\Sigma_{(2)}\in \mathbb R^{k\times k}$, $R_{(2)}\in \mathbb R^{(k+1)\times k}$ and   $Q_{(2)}^\top Q_{(2)}  = Q_{(2)} Q_{(2)} ^\top = I_k $.

(4) By (<ref>) we update the SVD of $U_{\ell+2}$ by

\begin{align}\label{SVD004}
\begin{split}
Q_{\ell+2}  &=   Q_{\ell+1} Q_{(2)} = Q Q_{(1)}Q_{(2)},\\
\Sigma_{\ell+2}  &=  \Sigma_{(2)},\\
R_{\ell+2} &=\left[\begin{array}{cc}
R_{\ell+1}& 0 \\
0 & 1
\end{array}\right]  R_{(2)} = \left[\begin{array}{cc}
\left[\begin{array}{cc}
R& 0 \\
0 & 1
\end{array}\right]R_{(1)}& 0 \\
0 & 1
\end{array}\right]  R_{(2)} .
\end{split}
\end{align}

<Ref>  updates the left singular matrix $Q_{\ell+2}$ by computing $Q_{\ell+1} = Q Q_{(1)}$ first and then compute $Q_{\ell+2} = Q_{\ell+1} Q_{(2)}$. The drawback of this update is its high computational cost. To avoid this, Brand proposed to  compute $\widetilde Q = Q_{(1)}Q_{(2)}$ first and then update the matrix $Q$ by computing  $Q \widetilde Q $; see the details in <Ref>.  For the right singular space, <Ref> updates  $R_{\ell+2}$  by computing $ R_{\ell+1} =\left[\begin{array}{cc}
R& 0 \\
0 & 1
\end{array}\right]  R_{(1)}$ first and then $  R_{\ell+2} =\left[\begin{array}{cc}
R_{\ell+1}& 0 \\
0 & 1
\end{array}\right]  R_{(2)} $. To save the  CPU time, one can use Brand's algorithm to update the right singular subspace, but has to compute the pseudo-inverse of the sub-matrix of $R_{(1)}$ and $R_{(2)}$.

From <Ref>, we have seen that both algorithms have suffered in  loss of the orthogonality of the left singular vector matrix.

Our new algorithm below can avoid the above issue without changing the complexity of <Ref>. To  describe our new algorithm, we need the following  fundamental identity.

\begin{align}\label{important_identity}
\left[\begin{array}{cc}
\left[\begin{array}{cc}
A& 0 \\
0 & 1
\end{array}\right] B& 0 \\
0 & 1
\end{array}\right]=	\left[\begin{array}{cc}
A& 0 \\
0 & I_2
\end{array}\right]  \left[\begin{array}{cc}
B& 0 \\
0 & 1
\end{array}\right],  \quad A\in \mathbb R^{m\times n}, \quad B\in \mathbb R^{(n+1)\times n}.
\end{align}

By the fundamental identity (<ref>) we can rewrite $ R_{\ell+2} $ as

\begin{align*}
R_{\ell+2} = \left[\begin{array}{cc}
R& 0 \\
0 & I_2
\end{array}\right]  \left[\begin{array}{cc}
R_{(1)}& 0 \\
0 & 1
\end{array}\right] R_{(2)}.
\end{align*}

By (<ref>) we know that  the SVD of $U_{\ell+2}$ is given by

\begin{align}\label{dede}
U_{\ell+2} = Q Q_{(1)}Q_{(2)}   \Sigma_{(2)}  R_{(2)}^\top \left[\begin{array}{cc}
R_{(1)}& 0 \\
0 & 1
\end{array}\right]^\top \left[\begin{array}{cc}
R& 0 \\
0 & I_2
\end{array}\right]^\top.
\end{align}

Motivated by  (<ref>), if we can compute $Q_{(1)}Q_{(2)}  $ and $\left[\begin{array}{cc}
R_{(1)}& 0 \\
0 & 1
\end{array}\right] R_{(2)}$ directly,  then the updates of the SVD of $U_{\ell+2}$ can possibly  avoid many matrix-multiply-matrix operations.

Let $\widetilde Q = Q_{(1)}Q_{(2)}$, $\widetilde\Sigma  = \Sigma_{(2)}$ and $\widetilde R = \left[\begin{array}{cc}
R_{(1)}& 0 \\
0 & 1
\end{array}\right] R_{(2)}$, then  $\widetilde Q\widetilde\Sigma \widetilde R^\top$ is a thin SVD of

\begin{align*}
\left[ \Sigma~\big|~  Q^\top W u_{\ell+1}  ~\big |~Q^\top W u_{\ell+2}   \right].
\end{align*}

First, by using (<ref>) we have

\begin{align*}
\left[ \Sigma~\big|~  Q^\top W u_{\ell+1}  ~\big |~Q^\top W u_{\ell+2}   \right] =  	\left[ Q_{(1)}\Sigma_{(1)} R_{(1)}^\top~\big |~Q^\top W u_{\ell+2}   \right].
\end{align*}

Note that $Q_{(1)}^\top Q_{(1)}  = Q_{(1)} Q_{(1)} ^\top = I_k $ and by the fundamental identity (<ref>)  we have

\begin{align*}
\left[ \Sigma~\big|~  Q^\top W u_{\ell+1}  ~\big |~Q^\top W u_{\ell+2}   \right] &=  Q_{(1)}	\left[ \Sigma_{(1)} ~\big |~Q_{(1)}^\top Q^\top W u_{\ell+2}   \right]\left[\begin{array}{cc}
R_{(1)}& 0 \\
0 & 1
\end{array}\right]^\top\\
& = Q_{(1)}	\left[ \Sigma_{\ell+1} ~\big |~Q_{\ell+1}^\top W u_{\ell+2}   \right]\left[\begin{array}{cc}
R_{(1)}& 0 \\
0 & 1
\end{array}\right]^\top,
\end{align*}

where we used the first two identities in (<ref>). Next we use (<ref>) to obtain:

\begin{align*}
\left[ \Sigma~\big|~  Q^\top W u_{\ell+1}  ~\big |~Q^\top W u_{\ell+2}   \right] = Q_{(1)}	Q_{(2)}\Sigma_{(2)} R_{(2)}^\top\left[\begin{array}{cc}
R_{(1)}& 0 \\
0 & 1
\end{array}\right]^\top = \widetilde Q\widetilde\Sigma \widetilde R^\top.
\end{align*}

<Ref> implies that we can update the SVD of  $U_{\ell+2}$ by computing a thin  SVD  of

\begin{align*}
\left[ \Sigma~\big|~  Q^\top W u_{\ell+1}  ~\big |~Q^\top W u_{\ell+2}   \right]: = \widetilde Q\widetilde\Sigma \widetilde R^\top
\end{align*}

first, and then update them by

\begin{align}\label{update2}
Q_{\ell+2} =  Q\widetilde Q,\qquad \Sigma_{\ell+2} = \widetilde \Sigma,\qquad  R_{\ell+2} = \left[\begin{array}{cc}
R& 0 \\
0 & I_{2}
\end{array}\right] \widetilde R.
\end{align}

Now we extend the above discussions to a more general situation.

Let $Q \Sigma R^\top$ be a rank $k$ truncated core SVD of $U_\ell$,  $u_{\ell+j}$ be the $(\ell+j)$-th column of $U$, $j=1,2,\ldots,s$. If their residuals are less than $\textup{\texttt{tol}}$ in <Ref> when project them onto the  space spanned by the columns of $Q$. Let $\widetilde Q \widetilde \Sigma  \widetilde R^\top $ be the standard  thin SVD of

\begin{align}\label{middle_M}
\left[
\Sigma  ~\big |~ Q^{\top}W u_{\ell+1}  ~\big |~ Q^{\top}W u_{\ell+2}  ~\big | \cdots ~\big | ~  Q^{\top}W u_{\ell+s} \right]\in \mathbb R^{k\times (k+s)}.
\end{align}

Then we suggest the following updates

\begin{align}\label{updates2}
Q \longleftarrow  Q\widetilde Q,\qquad \Sigma \longleftarrow \widetilde \Sigma,\qquad  R \longleftarrow \left[\begin{array}{cc}
R& 0 \\
0 & I_{s}
\end{array}\right] \widetilde R.
\end{align}

Moreover, the updates in (<ref>) are exactly the same with the output of <Ref>.

In real computations, we do not compute $Q \widetilde Q$, only need to  store $Q$ and $\widetilde Q$ for future updates. Since we only need to compute $\widetilde Q$ for $r$ times, then we do not loss the orthogonality of their multiplications. Hence, our answer to Brand's question in [4] is that the reorthogonalizations to $\widetilde Q$ are not needed.  Furthermore,  to update the right singular subspace, we do not compute the  pseudo-inverse of a matrix. It is worthwhile mentioning that  in practice we do not form the matrix $\left[\begin{array}{cc}
R& 0 \\
0 & I_{s}
\end{array}\right] $  by splitting $\widetilde R$ into $\left[\begin{array}{cc}
\widetilde R_1\\
\widetilde R_2
\end{array}\right] $, and then $ R  $ is updated by $\left[\begin{array}{cc}
R\widetilde R_1\\
\widetilde R_2
\end{array}\right] $. The total computational cost to update the matrix $R$ is $\mathcal O(nr^3)$. Since $r\ll m, n$, then the cost can be neglected. However, it is worthwhile motioning that the technique in [4] can be applied to update the right singular space by computing $r$ pseudo-inverse of a small matrix, and then the computational cost can be reduced to $\mathcal O(nr^2)$.

Once  completed the SVD of $U_{\ell+s}$, we then update the SVD of $U_{\ell+s+1}$.  To better describe the algorithm. We assume that $Q\Sigma R^\top$ is the core  SVD of $U_\ell$ and $\widehat Q \widehat\Sigma \widehat R^\top$ is the core SVD of $U_{\ell+s}$, then

\begin{align}\label{updates3}
\widehat Q =  Q\widetilde Q,  \qquad  \widehat R  =  \left[\begin{array}{cc}
R& 0 \\
0 & I_{s}
\end{array}\right] \widetilde R,
\end{align}

where $\widetilde Q$ and $\widetilde R$ are the left and right singular matrix of (<ref>), respectively.

By our assumption, the  residual of $u_{\ell+s+1}$ is larger than $\texttt{tol}$ when project it onto the  space spanned by the columns of $Q$. To keep the same accuracy, we should follow   Brand's idea to  compute the residual of $u_{\ell+s+1}$ that project  it onto the  space spanned by the columns of $\widehat  Q$.

In other words, for those columns, we have to compute their residual twice. Hence the computational cost of those repetitions can be high. Next, we show that this cost can be avoided.

To this end, we follow the idea  of Brand to construct the SVD of $\left[
U_{\ell+s}~\big |~ 	u_{\ell+s+1}\right]$ by

(1) letting   $e =u_{\ell+s+1} -  \widehat Q \widehat Q^\top Wu_{\ell+s+1}$ and let $p = \|e\|_W$ and  $\widetilde e = e/p$,

(2) finding the full SVD of $\bar Y = \left[\begin{array}{cc}
\widehat \Sigma & \widehat Q^\top u_{\ell+s+1} \\
0 & p
\end{array}\right]$,  and let $\bar Q \bar \Sigma \bar R^\top$ be the SVD of $\bar Y$,  then

(3) updating the SVD of $\left[
U_{\ell+s} ~\big |~	u_{\ell+s+1}
\right]$ by

\begin{align}\label{updates1}
\left[
U_{\ell+s} ~\big |~	u_{\ell+s+1}
\right]=\left(\left[
\widehat Q  ~\big |~ \widetilde  e
\right] \bar Q\right) \bar \Sigma \left(\left[\begin{array}{cc}
\widehat 	R & 0 \\
0 & 1
\end{array}\right] \bar R \right)^{\top}.
\end{align}

Next, we discuss the above three steps.   First,  the computation in  step (1) is not needed due to <Ref>.  Second, we do need  to compute $\widehat Q^\top u_{\ell+s+1} $ in step (2). However,  since $\widehat Q^\top u_{\ell+s+1} = \widetilde Q^\top Q^\top u_{\ell+s+1}$ and we have already computed  $Q^\top u_{\ell+s+1} $, then the update of  ${\widetilde Q}^\top Q^\top u_{\ell+s+1}$ is almost nothing. For the step (3),  by using (<ref>) we have

\begin{align}\label{ded}
\begin{split}
\left[
U_{\ell+s} ~\big |~	u_{\ell+s+1}
\right]&=\left(\left[
Q\widetilde  Q  ~\big |~  \widetilde  e
\right] \bar Q\right) \bar \Sigma \left(\left[\begin{array}{cc}
\widehat 	R & 0 \\
0 & 1
\end{array}\right] \bar R \right)^{\top}\\
& = \left(\left[
Q ~\big |~  \widetilde  e
\right]\left[\begin{array}{cc}
\widetilde Q& 0 \\
0 & 1
\end{array}\right]  \bar Q\right) \bar \Sigma \left(\left[\begin{array}{cc}
\widehat 	R & 0 \\
0 & 1
\end{array}\right] \bar R \right)^{\top}.
\end{split}
\end{align}

This suggest to update the left singular matrix $Q$ by adding the vector $\widetilde e$, compute $\left[\begin{array}{cc}
\widetilde Q& 0 \\
0 & 1
\end{array}\right]  \bar Q$ and then store it, the computational cost is $\mathcal O(r^2)$. For the right singular matrix, we directly compute it and the time complexity is $\mathcal O(nr^2)$. The step (3) will be executed $r$ times, then the whole CPU time for this step is $\mathcal O(nr^3)$.  We note that if we use the technique in [4], the computational cost can be reduced to $\mathcal O(nr^2)$ by computing $r$ pseudo-inverse of a matrix.

§.§ Reorthogonalization

As we claimed in <Ref>, small numerical errors in practice  can cause a loss of orthogonality.  A reorthogonalization to the  matrix $([
\widehat   Q  ~\big |~\widetilde  e
]  \bar Q)$ in [13, 7] is needed and the computational cost is high  due to the two  loops in <Ref>.

Observing the identity (<ref>), we note that the orthogonal matrices  $\bar Q$ and $\widetilde Q$ are small and only  $r$ of them multiplied together,  they can keep the  orthogonality in a reasonable numerical  precision. Hence we  do not need to perform reorthogonalizations to these small orthonormal matrices. This answers the question which was asked by Brand in [4].

Based on our numerical experiment in  <Ref>.  The matrix $[
Q  ~\big |~ \widetilde  e
]$ could loss its orthogonality very quickly.
Therefore,  we  consider  reorthogonalizations  to   the matrix $[
Q  ~\big |~ \widetilde  e
]$. We note that this is necessary and was not included in Brand's original algorithm.

By the  construction of $[
Q  ~\big |~  \widetilde  e
]$, we can reorthogonalize it recursively: only apply the $W$-weighted Gram-Schmidt to  the new  adding vector $\widetilde e$ since the matrix $ Q$ has already been  reorthogonalized in previous steps if the weighted inner product between $\widetilde e$ and the first column of $ Q$ is larger than some tolerance. Furthermore, if $\widetilde e=0$, then the last column will be truncated. In other words,  we do not need to perform the reorthogonalization in the truncation part ($\|p\|_W<\texttt{tol}$); see <Ref> for more details.

(Incremental SVD (III))

Input: $Q\in \mathbb R^{m\times k}$, $\Sigma\in \mathbb R^{k\times k}$, $R\in \mathbb R^{\ell\times k}, u_{\ell+1}\in \mathbb R^m$, $W\in \mathbb R^{m\times m}$, $\texttt{tol}$,  $V$, $Q_0$, $q$
[1]
Set $d = Q^\top (W u_{\ell+1})$; $e = u_{\ell+1} - Qd$;   $p = (e^\top W e)^{1/2}$;
$q = q+1$;
Set $V\{q\} =   Q_0^\top d$;
Set $Y = \left[
\Sigma ~\big |~ \texttt{cell2mat}(V)\right]$;
$[Q_Y, \Sigma_Y, R_Y] =  \texttt{svd(Y, {\color{magenta} `econ'})}$;
Set $Q_0 = Q_0Q_Y$,   $\Sigma = \Sigma_Y$,   $\mathcal R_1 = R_Y(1:k,1:\texttt{end-1})$,  $\mathcal R_2 = R_Y(k+1,1:{\texttt{end-1}})$,

$R =  \left[\begin{array}{c}
R \mathcal R_1\\
\mathcal R_2
\end{array}\right]$;
Set $d = Q_Y^\top d$
Set  $  e = e/p$;
$|  e^\top W Q(:,1) |> \texttt{tol}$
$e = e - Q(Q^\top(We))$;  $p_1= (e^\top W e)^{1/2}$;  $ e = e/p_1$;
Set $Y = \left[\begin{array}{cc}
\Sigma & d\\
0 & p
\end{array}\right]$;

$[Q_Y, \Sigma_Y, R_Y] = $  $(Y)$;
Set $
Q_0=  \left[\begin{array}{cc}
Q_0 & 0 \\
0 & 1
\end{array}\right]Q_Y, \quad  Q = [Q ~\big|~ e],\quad  \Sigma = \Sigma_Y, \quad  R =  \left[\begin{array}{cc}
R & 0 \\
0 & 1
\end{array}\right] R_Y$.
Set $V =$ []; $q=0$;
return $Q$, $\Sigma$, $R$, $V$, $Q_0$, $q$

§.§ The complexity of <Ref>

We restrict ourselves the case $W = I$ to discuss the  <Ref>. Recall that $U\in \mathbb R^{m\times n}$ with low rank $r$.

1. Line 1:  this is exactly the same with <Ref>, the complexity is $\mathcal O(mnr)$.

2. Lines  3-4:  the vector $d\in \mathbb R^k$ with $k\le r$, hence this step is almost nothing for both CPU time and memory storage. Brand's method does not have this step.

3. Lines 7-11:  the main  cost is the update of $R$ and the CPU time   is less than  $\mathcal (nr^3)$. As we mentioned before, the cost can be reduced to $\mathcal O(nr^2)$. We leave it here since $r\ll m$.

4. Lines 14-16:  the CPU time of the reorthogonalization is $\mathcal O(mr)$ each time and this part will be executed $\mathcal O(r)$ times, hence the whole complexity is $\mathcal O(mr^2)$.
This part was not concluded  in <Ref> and it is necessary for some data sets.

5. Lines 17-20:  the whole complexity is less than $\mathcal O(nr^3)$. Again, this can be reduced to $\mathcal O(nr^3)$ by using the technique  in [4].

By our assumption, $r\ll m, n$, hence the CPU time of <Ref> is $\mathcal O(mnr)$.

We note that the output of <Ref>, $V$,  may be  not empty. This implies that the output of <Ref> is not the SVD of $U$. Hence we have to update the SVD for the vectors in $V$.   We give the implementation of this step  in <Ref> and the complexity is $\mathcal O(mr^2)$. Finally,  we complete the full implementation in <Ref>.
(Incremental SVD (III) final check)

Input: $Q\in \mathbb R^{m\times k}$, $\Sigma\in \mathbb R^{k\times k}$, $R\in \mathbb R^{\ell\times k}$, $V$, $Q_0$, $q$
[1]
Set $Y = \left[
\Sigma ~\big |~ \texttt{cell2mat}(V)\right]$;
$[Q_Y, \Sigma_Y, R_Y] =  \texttt{svd(Y, {\color{magenta} `econ'})}$;
Set $Q = Q(Q_0Q_Y)$,   $\Sigma = \Sigma_Y$,    $\Sigma = \Sigma_Y$,   $\mathcal R_1 = R_Y(1:k,1:\texttt{end-1})$,
$\mathcal R_2 = R_Y(k+1,1:{\texttt{end-1}})$,  $R =  \left[\begin{array}{c}
R \mathcal R_1\\
\mathcal R_2
\end{array}\right]$;
Set $Q = QQ_0$;
return $Q$, $\Sigma$, $R$.

(Fully incremental SVD (III))

Input:  $W\in \mathbb R^{m\times m}$, $\texttt{tol}$
[1]
Get $u_1$;
$[Q, \Sigma, R]$=($u_1$, $W$);   % <Ref>
Set $V = $ []; $Q_0=1$; $q = 0$;
$\ell=2,\ldots, n$
Get $u_\ell$
$[q,V, Q_0, Q, \Sigma, R]$=($q,V, Q_0, Q, \Sigma, R, u_{\ell},W, \texttt{tol}$);    % <Ref>
$[Q, \Sigma, R]$=($q,V,Q_0, Q, \Sigma, R$);   % <Ref>
return $Q, \Sigma, R$

In this example, we revisit <Ref> to test the efficiency and accuracy of  <Ref>. We show the orthogonality of the left singular vector matrix  in <Ref>, where

\begin{align*}
\mathcal E_W(Q) = \|I - Q^\top W Q\|, \qquad 	\mathcal E_I(\widetilde Q) = \|I - \widetilde Q^\top \widetilde Q\|.
\end{align*}

We see that our <Ref> keeps a good orthogonality.

<Ref>: Top: the orthogonality of $\widetilde Q$ and $Q$ with $W= I$.  Bottom: the orthogonality of $\widetilde Q$ and $Q$ with $W= M$.

Next, we show the CPU time (seconds)  of the main part of the <Ref> in <Ref>. 	 The whole simulations take about 34 and  118 seconds for  $W = I$ and $W = M$, respectively. Comparing with <Ref>, the CPU times is reduced greatly, especially with a weighted norm setting ($W\neq I$). Furthermore, the projection part (line 1 in <Ref>) takes about $99 \%$ of the whole simulation time, this implies there is litter space to improve the efficiency of <Ref>.

Line(s) in <Ref>	   	1      3-4    7-11    14-16    17-20

$W = I$
33.457    	0.2416       0.0327   0.1251   0.2882

$W = M$
116.41    0.1870      0.0650   0.6526   0.7120

<Ref>:  The CPU time  (seconds)  for two  different weighted norms by using <Ref>.

§ ANOTHER TRUNCATION

For many data sets, they may have a large number of nonzero singular values but most of them are very small. Without  truncating those small singular values, the  incremental SVD algorithm maybe  keep  computing them and hence the computational cost  is huge.  Following <cit.>, we  set another truncation if the last few singular values are less than a  tolerance.

Next, we show that it is only necessary to check the last singular value.
Let $Q\Sigma R^\top$ be a core SVD of $U_\ell$, $e = u_{\ell+1} - QQ^\top W u_{\ell+1}$ and  $p = \|e\|_W>\texttt{tol}$, then

\begin{align*}
Q ~ \big|~  \widetilde e \right]\underbrace{\left[\begin{array}{cc}
\Sigma & Q^{\top}W u_{\ell+1} \\
0 & p
\end{array}\right]}_{Y}\left[\begin{array}{cc}
R & 0 \\
0 & 1
\end{array}\right]^{\top} = \left[
Q ~ \big|~  \widetilde e \right] Q_Y \Sigma_Y R_Y^\top\left[\begin{array}{cc}
R & 0 \\
0 & 1
\end{array}\right]^{\top}.
\end{align*}

Next, we estimate the singular values of $Y$.

Assume that $\Sigma=\textup{\texttt{diag}}\left(\sigma_{1}, \ldots, \sigma_{k}\right)$ with $\sigma_{1}\ge \ldots \ge \sigma_{k}$, and  $\Sigma_Y=\textup{\texttt{diag}}\left(\mu_{1},  \ldots, \mu_{k+1}\right)$ with $\mu_{1}\ge \ldots \ge \mu_{k+1}$. Then we have

\begin{gather}
\mu_{k+1} \le p,\label{Singular1}\\
\mu_{k+1}\le \sigma_{k} \le \mu_k\le \sigma_{k-1}\le \ldots\le \sigma_1\le  \mu_1.\label{Singular2}
\end{gather}

First, we consider  $Y^\top Y $,

\begin{align*}
Y^\top Y  = \left[\begin{array}{cc}
\Sigma^2 & \Sigma Q^{\top}W u_{\ell+1} \\
u_{\ell+1}^\top WQ\Sigma & p^2+u_{\ell+1}^\top WQQ^\top Wu_{\ell+1}
\end{array}\right].
\end{align*}

By using the Cauchy interlace theorem <cit.> we have

\begin{align*}
\mu_{k+1}\le \sigma_{k} \le \mu_k\le \sigma_{k-1}\le \ldots\le \sigma_1\le  \mu_1.
\end{align*}

Next, we take $e = [0,0,\ldots, 1]$ to obtain

\begin{align*}
p^2 = e Y Y^\top e \ge \min_{\|x\|=1} x Y Y^\top x = \lambda_{\min} ( Y Y^\top)= \mu_{k+1}^2.
\end{align*}

Inequality (<ref>) implies that the last singular value of $Y$  can possibly be very small, no matter how large of  $p$. That is to say the tolerance we set for $p$ can not avoid the algorithm keeping compute very small singular values. Hence, another truncation is needed if the data has  a large number of very small  singular values.  Inequality (<ref>) guarantees us  that only the last singular value of $Y$ can possibly less than the tolerance. Hence, we only need to check the last one.

(i) If $\Sigma_Y(k+1,k+1)\ge\texttt{tol}$, then

\begin{align}\label{Updateshashope}
Q\longleftarrow \left[
Q~ \big|~  \widetilde e
\right]Q_Y\quad   \Sigma \longleftarrow  \Sigma_Y\quad   R \longleftarrow  \left[\begin{array}{cc}
R & 0 \\
0 & 1
\end{array}\right] R_Y.
\end{align}

(ii) If $\Sigma_Y(k+1,k+1)<\texttt{tol}$, then

\begin{align}\label{Updatesnohope}
Q\longleftarrow \left[
Q~ \big|~  \widetilde e
\right]Q_Y(:,1:k)\quad   \Sigma \longleftarrow  \Sigma_Y(1:k,1:k)\quad   R \longleftarrow  \left[\begin{array}{cc}
R & 0 \\
0 & 1
\end{array}\right] R_Y(:,1:k).
\end{align}

As we discussed in <Ref>, we can avoid time complexity $\mathcal O(mr^2)$ in (<ref>) by only updating the small matrix $Q_Y$ and storing $Q_Y$ and   $ \left[
Q~ \big|~  \widetilde e
\right]$ separately. Unfortunately, to truncate the very small singular values, it seems very difficult that we can find two orthonormal matrices $\widehat Q\in \mathbb R^{m\times k}$ and $\widetilde  Q\in \mathbb R^{k\times k}$ in $\mathcal O(r^2)$ time such that

\begin{align*}
\left[
Q~ \big|~  \widetilde e
\right]Q_Y(:,1:k) = \widehat Q  \widetilde Q.
\end{align*}

Hence, we directly use (<ref>) and  (<ref>) to update the matrix $Q$ and $R$.  Assume the data has $d$ nonzero singular values, then this new truncation needs $\mathcal O((m+n)r^2d)$ time, here $r$ is the number of singular values which are larger than the tolerance.

(Incremental SVD (IV))

Input: $Q\in \mathbb R^{m\times k}$, $\Sigma\in \mathbb R^{k\times k}$, $R\in \mathbb R^{\ell\times k}, u_{\ell+1}\in \mathbb R^m$, $W\in \mathbb R^{m\times m}$, $\texttt{tol}$,  $V$, $Q_0$, $q$
[1]
Set $d = Q^\top (W u_{\ell+1})$; $e = u_{\ell+1} - Qd$;   $p = (e^\top W e)^{1/2}$;
$q = q+1$;
Set $V\{q\} = d$;
Set $Y = \left[
\Sigma ~\big |~ \texttt{cell2mat}(V)\right]$;
$[Q_Y, \Sigma_Y, R_Y] =  \texttt{svd(Y, {\color{magenta} `econ'})}$;
Set $Q_0 = Q_0Q_Y$,   $\Sigma = \Sigma_Y$,    $\mathcal R_1 = R_Y(1:k,1:\texttt{end-1})$,
$\mathcal R_2 = R_Y(k+1,1:{\texttt{end-1}})$,  $R =  \left[\begin{array}{c}
R \mathcal R_1\\
\mathcal R_2
\end{array}\right]$;
Set $d = Q_Y^\top d$
Set  $  e = e/p$;
$|  e^\top W Q(:,1) |> \texttt{tol}$
$e = e - Q(Q^\top(We))$;  $p_1= (e^\top W e)^{1/2}$;  $ e = e/p_1$;
Set $Y = \left[\begin{array}{cc}
\Sigma & d\\
0 & p
\end{array}\right]$;

$[Q_Y, \Sigma_Y, R_Y] = $  $(Y)$;
Set $
Q_0=  \left[\begin{array}{cc}
Q_0 & 0 \\
0 & 1
\end{array}\right]Q_Y, \quad  Q = [Q ~\big|~ e]Q_0,\quad  \Sigma = \Sigma_Y$,   $\mathcal R_1 = R_Y(1:k,1:\texttt{end-1})$,
$\mathcal R_2 = R_Y(k+1,1:{\texttt{end-1}})$,  $R =  \left[\begin{array}{c}
R \mathcal R_1\\
\mathcal R_2
\end{array}\right]$ ,   $Q_0 = I_{k+1}$;
Set $
Q_0=  \left[\begin{array}{cc}
Q_0 & 0 \\
0 & 1
\end{array}\right]Q_Y,  \qquad  Q = [Q ~\big|~ e]Q_0(:,1:k)$,
$\Sigma = \Sigma_Y(1:k,1:k),  \qquad  R =  \left[\begin{array}{cc}
R & 0 \\
0 & 1
\end{array}\right] R_Y(:, 1:k), \quad Q_0 = I_{k}$
$V =$ []; $q = 1$;
return $Q$, $\Sigma$, $R$, $V$, $Q_0$, $q$

Next,  we complete the full implementation in <Ref>.
(Incremental SVD (IV) final check)

Input: $Q\in \mathbb R^{m\times k}$, $\Sigma\in \mathbb R^{k\times k}$, $R\in \mathbb R^{\ell\times k}$, $V$, $q$
[1]
Set $Y = \left[
\Sigma ~\big |~ \texttt{cell2mat}(V)\right]$;
$[Q_Y, \Sigma_Y, R_Y] =  \texttt{svd(Y, {\color{magenta} `econ'})}$;
Set $Q = QQ_Y$,   $\Sigma = \Sigma_Y$,    $\mathcal R_1 = R_Y(1:k,1:\texttt{end-1})$,
$\mathcal R_2 = R_Y(k+1,1:{\texttt{end-1}})$,  $R =  \left[\begin{array}{c}
R \mathcal R_1\\
\mathcal R_2
\end{array}\right]$;
return $Q$, $\Sigma$, $R$.

(Fully incremental SVD (IV))

Input:  $W\in \mathbb R^{m\times m}$, $\texttt{tol}$
[1]
Get $u_1$;
$[Q, \Sigma, R]$=($u_1$, $W$);   % <Ref>
Set $V = $ []; $Q_0=1$; $q=1$;
$\ell=2,\ldots, n$
Get $u_\ell$
$[q,V, Q_0, Q, \Sigma, R]$=($q,V, Q_0, Q, \Sigma, R, u_{\ell},W, \texttt{tol}$);    % <Ref>
$[Q, \Sigma, R]$=($q,V,Q_0, Q, \Sigma, R$);   % <Ref>
return $Q, \Sigma, R$

§ CONCLUSION

In this paper, we  answered the open question which was asked by Brand  in [4].  We prove that the output of our new  algorithm is the same with  Brand's original algorithm. Numerical experiments showed that our new algorithm not only saves the CPU, but also keeps a very good orthogonality for the left singular vector matrix.  Im the future, we will compare our new algorithm with many other incremental SVD algorithms in [1]. Besides this, we will apply this new algorithm for many different problems, such as model order reduction and fractional PDEs.

§ ACKNOWLEDGMENT
The author would like to thank Dr. Lihan Wang  for the proof of <Ref>.

[1]
C. Bach, D. Ceglia, L. Song, and F. Duddeck, Randomized low-rank
approximation methods for projection-based model order reduction of large
nonlinear dynamical problems, Internat. J. Numer. Methods Engrg., 118
(2019), pp. 209–241, <https://doi.org/10.1002/nme.6009>.

[2]
M. W. Berry, Large-scale sparse singular value computations, The
International Journal of Supercomputing Applications, 6 (1992), pp. 13–49.

[3]
M. Brand, Incremental singular value decomposition of uncertain data
with missing values, in European Conference on Computer Vision, Springer,
2002, pp. 707–720.

[4]
M. Brand, Fast low-rank modifications of the thin singular value
decomposition, Linear Algebra Appl., 415 (2006), pp. 20–30,

[5]
H. Fareed and J. R. Singler, A note on incremental POD algorithms
for continuous time data, Appl. Numer. Math., 144 (2019), pp. 223–233,

[6]
H. Fareed and J. R. Singler, Error analysis of an incremental proper
orthogonal decomposition algorithm for PDE simulation data, J. Comput.
Appl. Math., 368 (2020), pp. 112525, 14,

[7]
H. Fareed, J. R. Singler, Y. Zhang, and J. Shen, Incremental proper
orthogonal decomposition for PDE simulation data, Comput. Math. Appl., 75
(2018), pp. 1942–1960, <https://doi.org/10.1016/j.camwa.2017.09.012>.

[8]
M. Ghashami, E. Liberty, J. M. Phillips, and D. P. Woodruff, Frequent directions: simple and deterministic matrix sketching, SIAM J.
Comput., 45 (2016), pp. 1762–1792, <https://doi.org/10.1137/15M1009718>.

[9]
S.-G. Hwang, Cauchy's interlace theorem for eigenvalues of
Hermitian matrices, Amer. Math. Monthly, 111 (2004), pp. 157–159,

[10]
M. A. Iwen and B. W. Ong, A distributed and incremental SVD
algorithm for agglomerative data analysis on large networks, SIAM J. Matrix
Anal. Appl., 37 (2016), pp. 1699–1718,

[11]
L. Lin and Y. Tong, Low-rank representation of tensor network
operators with long-range pairwise interactions, SIAM J. Sci. Comput., 43
(2021), pp. A164–A192, <https://doi.org/10.1137/19M1287067>.

[12]
N. Mastronardi, E. E. Tyrtyshnikov, and P. Van Dooren, A fast
algorithm for updating and downsizing the dominant kernel principal
components, SIAM J. Matrix Anal. Appl., 31 (2010), pp. 2376–2399,

[13]
G. M. Oxberry, T. Kostova-Vassilevska, W. Arrighi, and K. Chand, Limited-memory adaptive snapshot selection for proper orthogonal
decomposition, Internat. J. Numer. Methods Engrg., 109 (2017), pp. 198–217,


11institutetext:  School of Physics, Beijing Institute of Technology, Beijing
100081, China
CAS Key Laboratory of Theoretical Physics, Institute of Theoretical Physics,
Chinese Academy of Sciences, Beijing 100190, China
Max-Planck-Institut für Festkörperforschung, Heisenbergstrasse 1, D-70569
Stuttgart, Germany

# Topological $p_{z}$-wave nodal-line superconductivity with flat surface
bands in the AHxCr3As3 (A=Na, K, Rb, Cs) superconductors

Juan-Juan Hao 111These two authors contributed equally to this work. 11 Ming
Zhang 222These two authors contributed equally to this work. 11 Xian-Xin Wu
2233 Fan Yang<EMAIL_ADDRESS>11112233

###### Abstract

We study the pairing symmetry and the topological properties of the hydrogen-
doped ACr3As3 superconductors. Based on our first-principle band structure
with spin-orbit-coupling (SOC), we construct tight-binding model including the
on-site SOC terms, equipped with the multi-orbital Hubbard interactions. Then
using the random-phase-approximation (RPA) approach, we calculate the pairing
phase diagram of this model. Our RPA results yield the triplet $p_{z}$-wave
pairing in the $\left(\uparrow\downarrow+\downarrow\uparrow\right)$ spin
channel to be the leading pairing symmetry all over the experiment relevant
hydrogen-doping regimes. This pairing state belongs to the spin-U(1)-symmetry
protected time-reversal-invariant topological nodal-line superconductivity.
Determined by the momentum-dependent topological invariant
$Z\left(k_{x},k_{y}\right)$, the whole $(001)$ surface Brillouin zone is
covered with topological flat bands with different regimes covered with
different numbers of flat surface bands, which can be detected by the scanning
tunneling microscope experiments.

##  1 Introduction

In recent years, the quasi-one-dimensional (Q1D) Cr-based superconductors have
attracted great research interests. There are two Q1D Cr-based superconductors
families, i.e. the A2Cr3As3 (abbreviated as the 233) family (A=Na, K, Rb, Cs)
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,
22, 23, 24, 25, 26, 27, 28, 29] and the ACr3As3 (abbreviated as the 133)
family [33, 34, 35, 36, 37, 39, 30, 32, 31, 38], with both sharing the same
Q1D alkaline-cations-separated [(Cr3As3)2-]∞ double-walled subnanotubes (DWSN)
structure. The 233 family were synthesized earlier, and have shown
unconventional pairing nature [2, 20, 22, 3, 4, 5, 21, 7, 40, 23, 25] with
possible triplet pairing symmetry[2, 20, 22, 21, 40, 25, 41, 18]. The 133
family were synthesized slightly later, obtained by deintercalating half of
the alkali metal atoms with ethanol[33, 34, 35, 36]. The primitive ACr3As3
have a spin glass ground state [33, 34, 39], which is consistent with the
relaxation behavior of density functional theory (DFT) calculations[42, 39,
37, 38] and the muon spin-rotation ($\mu$SR) experiments [5]. However, when
enough amount of hydrogen atoms are intercalated in the material,
superconductivity (SC) emerges [43, 44, 45, 46], with highest $T_{c}$ around 8
K[30, 31]. Therefore, the correct chemical formula for“ACr3As3” should be
AHxCr3As3. It is the difference in H content that leads to the difference in
ground state properties between superconducting and non-superconducting
samples [43, 44, 45, 46]. The proximity of the magnetic ordered state to the
SC in the phase diagram implies spin fluctuation-driven pairing mechanism,
similar to the cuprates and the iron-based superconductors.

Figure 1: (Color online). Crystal and band structures of KHCr3As3. (a) Top and
side views for the crystal structure. (b) Band structures with (black) and
without (red) SOC along the high-symmetry lines. The five bands near the Fermi
level are marked as $\alpha$ \- $\varepsilon$ respectively. (c) and (d) are
FSs without and with SOC, with the high symmetric points marked in the 3D BZ.

The density-functional theory (DFT) based calculations suggest that the low-
energy degrees of freedom for the AHxCr3As3 is the Cr-3d orbitals[43, 44, 45,
46, 47], similar with those of the 233 family[2, 20, 22, 3, 6, 24, 13, 25, 15,
26, 27, 14, 28, 16, 17, 29], implying the strong electron-electron
interactions can account for the origin of the SC in the system. Furthermore,
the presence of H induces effective electron doping, which confirms that H
doping has the effect of raising Fermi level [43, 47, 44]. A remarkable
property of the band structure of the AHxCr3As3 is the presence of the
Lifshitz transition under moderate hole doping, during which its 3D FSs are
changed to two disconnected Q1D FS sheets [44, 47]. It’s clarified in our
previous work[47] based on the random-phase approximation (RPA) [48, 49, 50,
51, 52, 54, 53, 55, 56] that the presence of the type-II VHS[57, 58, 59, 60,
61] during the Lifshitz transition doping regimes has driven the triplet
$p_{z}$-wave pairing symmetry. It’s interesting to investigate the topological
properties of this triplet pairing state. However, as the spin-orbit coupling
(SOC) has not been considered in that work, the three spin components of the
triplet pairing are exactly degenerate, which leads to uncertainty of the
superconducting ground state.

In this paper, we report our RPA study on the pairing symmetry and the
topological properties of the KHxCr3As3, based on our first-principle DFT band
structure considering SOC. Our main results are as follow. Due to the
inversion symmetry of the material, the double-folded spin degeneracy survives
the SOC. The main effect of the SOC on the band structure lies in the
elimination of the double degeneracy along the $\Gamma$-$A$ and the $K$-$H$
lines. Starting from our previous six-orbital tight-binding (TB) model[47]
without SOC, we add the dominating on-site SOC conserving the spin-U(1)
symmetry (SUS), leading to a band structure well consistent with the DFT one.
Then we add the multi-orbital Hubbard interactions into this model and perform
the RPA calculations. The leading pairing symmetry obtained from our RPA
calculations is triplet $p_{z}$-wave pairing in the $S_{z}=0$ channel, i.e.
the $\uparrow\downarrow+\downarrow\uparrow$. This pairing state belongs to the
SUS-protected time-reversal-invariant (TRI) topological nodal-line SC defined
in the Ref[62], which includes the 233 family as a good example. Determined by
the momentum-dependent topological winding number, the whole surface Brillouin
zone (BZ) is covered with topological flat bands, with different regimes
covered with different numbers of flat bands, as shown by our surface spectrum
in the $(001)$ surface BZ. Our results appeal more experimental access into
this intriguing family.

##  2 The DFT Band Structure

The crystal structure of KHCr3As3 is shown in Fig. 1(a). This compound
consists of an infinite linear chain of DWSN [(Cr3As3)2-]∞ extended along the
c-axis, forming a Q1D structure, which is interconnected by K+ cations. The
hydrogen atoms locate at the center of the Cr6 octahedrons in such a DWSN
structure[43, 44, 45, 46, 47]. Compared with the crystal structure of K2Cr3As3
[2, 20, 22, 3, 6, 24, 13, 25, 15, 26, 27, 14, 28, 16, 17, 29] and KCr3As3 [56,
33, 34, 35, 39, 43, 36, 37, 30, 32, 31, 38], KHCr3As3 is more similar to
KCr3As3, which has a centrosymmetric structure with space group $P6_{3}/m$
(No.176)) and point group $C_{6h}$.

The band structure of KHCr3As3 is calculated using DFT theory, adopting the
ultra-soft pseudopotential (USPP) of Perdew-Burke-Ernzerhof (PBE) [63] in
QUANTUM ESPRESSO (QE) [64]. We use $5\times 5\times 11$ sparse k-points in the
entire Brillouin Zone (BZ) for self-consistent calculation. We set the
truncated kinetic energy of the wave function to 60Ry and the charge density
to 400Ry. The experimental lattice parameters of the crystal structure of
KHCr3As3 were adopted: $a=9.09481$ Å and $c=4.17703$ Å [43], with only the
atomic positions within the unit cell to be determined by optimization.

The obtained band structure of KHCr3As3 without considering SOC is shown in
Fig. 1(b) by the red curves. It can be seen that there are five bands near the
Fermi level, labeled as $\alpha,\beta,\gamma,\delta$ and $\epsilon$,
respectively. The three bands across the Fermi level are the $\alpha$-,
$\beta$\- and $\gamma$\- bands. These bands form two Q1D Fermi surfaces
($\alpha$ and $\beta$) and a 3D Fermi surface ($\gamma$), as shown in Fig.
1(c). Note that the $\beta$ and $\gamma$\- Fermi surfaces (FSs) intersect at
the $k_{x}=k_{y}=0,k_{z}=k_{F}$ point on the $\Gamma$-$A$ line. Compared with
the band structure of KCr3As3 [56, 33, 34, 35, 39, 43, 36, 37, 30, 32, 31,
38], it can be found that the Fermi level rises when H is inserted, which
reflects that H is used as an electron donor in KHCr3As3, consistent with
Ref[43].

After the band structure calculation, a six-band low-energy tight-binding
model has been constructed in our previous work[47]. Briefly, the relevant
orbital degrees of freedom include the Cr-3d${}_{z^{2}}$, 3d${}_{x^{2}-y^{2}}$
and 3dxy orbitals from the two inequivalent Cr atoms within each unit cell,
and the effective hopping parameters are obtained through the calculations of
maximal localization for these orbitals using Wannier90 [65]. The obtained TB
Hamiltonian with $C_{6h}$ point-group symmetry can be expressed in the
momentum space as,

$\displaystyle H_{{\rm
TB}}=\sum_{\mathbf{k}\mu\nu\sigma}h_{\mu\nu}(\mathbf{k})c^{\dagger}_{\mathbf{k}\mu\sigma}c_{\mathbf{k}\nu\sigma},$
(1)

Here $\mu,\nu=1,\cdots,6$ indicate the orbital-sublattice indices, containing
the $d_{z^{2}}$, $d_{x^{2}-y^{2}}$ and $d_{xy}$ orbitals of A or B
sublattices. The elements of the $h(\bm{k})$ matrix is obtained through
Fourier transformation using the real-spacing hopping-integral data provided
in our previous work[47]. The band structure thus obtained is well consistent
with that of the DFT at low energy near the Fermi level, see Ref[47].

##  3 The SOC Hamiltonian

Figure 2: (Color online). Band structure (a) and FS (b) of our band model (3)
for KHCr3As3 with on-site SOC.

Using DFT calculations, we obtain the band structure with SOC. The band
structures with and without SOC are compared in Fig. 1(b). From the
comparison, we find the two band structures are overall similar except that
the degeneracies along the high-symmetric lines are lift up by the SOC.
Without SOC, the $\beta$\- and $\gamma$\- bands are degenerate at the
$\Gamma$\- and $K$\- points, or in general along the $\Gamma$-$A$ and $K$-$H$
lines. Similarly, the $\delta$\- and $\varepsilon$\- bands are also degenerate
along these high-symmetric lines. However, with SOC, the energies of the
$\beta$\- and $\gamma$\- bands at the $\Gamma$\- and $K$\- points are split by
about 0.02eV and 0.05eV, and the energies of the $\delta$\- and
$\varepsilon$\- bands at the two momenta are split by about 0.08eV and
0.006eV, respectively. Note that the spin degeneracy is maintained with SOC,
due to the combination of the inversion and time-reversal (TR) symmetries. The
FS with SOC is shown in Fig. 1(d). Now the $\beta$\- and $\gamma$\- FSs no
longer intersect at the $(0,0,k_{F})$ point on the $\Gamma$-$A$ line.

In the following studies involving electron-electron interactions, we need a
simple TB formula for the SOC part of the Hamiltonian. For simplicity, we only
consider the dominating part of the SOC, i.e. the on-site SOC term. From
combination of the rotation symmetry, the mirror-reflection symmetry about the
$xy$-plane, and the TR symmetry of the system, the only possible on-site SOC
takes the following formula[62, 27],

$\displaystyle H_{\text{SOC}}$ $\displaystyle=$ $\displaystyle
i\lambda_{\text{SOC}}\sum_{\mathbf{i}\sigma}\sigma[c^{\dagger}_{\mathbf{i}2\sigma}c_{\mathbf{i}3\sigma}-c^{\dagger}_{\mathbf{i}3\sigma}c_{\mathbf{i}2\sigma}]+i\lambda_{\text{SOC}}$
(2)
$\displaystyle\times\sum_{\mathbf{i}\sigma}\sigma[c^{\dagger}_{\mathbf{i}5\sigma}c_{\mathbf{i}6\sigma}-c^{\dagger}_{\mathbf{i}6\sigma}c_{\mathbf{i}5\sigma}].$

This SOC Hamiltonian also satisfies the inversion symmetry $\hat{R}$ which
dictates $c^{\dagger}_{\mathbf{i}\mu\sigma}\to
c^{\dagger}_{\left(\hat{R}\mathbf{i}\right)\mu\sigma}$.

Our full model band Hamiltonian with SOC is given by

$H_{\text{band}}=H_{\text{TB}}+H_{\text{SOC}}.$ (3)

Setting $\lambda_{\text{SOC}}=10$ meV, the obtained band structure and FS for
this model are shown in Fig. 2 (a) and (b) respectively. Comparing Fig. 2 (a)
with Fig. 1 (b), it’s found that this band structure is qualitatively
consistent with the DFT one with SOC. Firstly, the spin degeneracy is kept in
this band structure, i.e.
$\varepsilon_{\mathbf{k}\uparrow}^{\alpha}=\varepsilon_{\mathbf{k}\downarrow}^{\alpha}\equiv\varepsilon_{\mathbf{k}}^{\alpha}$,
due to the combined TRS and inversion symmetry, where
$\varepsilon_{\mathbf{k}\sigma}^{\alpha}$ denotes the band energy for the spin
$\sigma$ and momentum $\mathbf{k}$ in the $\alpha$-th band. Secondly, the
degeneracies between the $\beta$\- and $\gamma$\- bands and those between the
$\delta$\- and $\varepsilon$\- bands at the $\Gamma$\- and $K$\- points and
along the $\Gamma$-$A$ lines (as well as along the unshown $K$-$H$ lines) are
lift up. Furthermore, the FS for the up or down spin shown in Fig. 2 (b) for
this model is also qualitatively consistent with the DFT FS with SOC shown in
Fig. 1(d).

##  4 The RPA pairing phase diagram

According to our previous RPA result for this material without considering
SOC[47], the triplet $p_{z}$-wave pairing is the leading pairing symmetry
whose three spin components are exactly degenerate. Here we consider the SOC
to lift up such degeneracy. We adopt the following extended Hubbard model
Hamiltonian in our study,

$\displaystyle H=$ $\displaystyle H_{\text{band}}+H_{\text{int}}$
$\displaystyle H_{\text{int}}=$ $\displaystyle
U\sum_{i\mu}n_{i\mu\uparrow}n_{i\mu\downarrow}+V\sum_{i,\mu<\nu}n_{i\mu}n_{i\nu}+J_{H}\sum_{i,\mu<\nu}$
$\displaystyle\Big{[}\sum_{\sigma\sigma^{\prime}}c^{+}_{i\mu\sigma}c^{+}_{i\nu\sigma^{\prime}}c_{i\mu\sigma^{\prime}}c_{i\nu\sigma}+(c^{+}_{i\mu\uparrow}c^{+}_{i\mu\downarrow}c_{i\nu\downarrow}c_{i\nu\uparrow}+h.c.)\Big{]}$
(4)

Here, the interaction parameters $U$, $V$, and $J_{H}$ denote the intra-
orbital, inter-orbital Hubbard repulsion, and the Hund’s rule coupling (as
well as the pair hopping) respectively, which satisfy the relation
$U=V+2J_{H}$. Note that this Hamiltonian conserves the SUS, i.e. the total
$S_{z}$ is a good quantum number.

Figure 3: (Color online). Doping $x$ dependences of the largest pairing
eigenvalues $\lambda$ for the $p_{z}$-, $(d_{x^{2}-y^{2}},d_{xy})$-,
$(p_{x},p_{y})\times p_{z}$\- and $s$-wave pairings for fixed $U=0.25$
eV,$J_{H}=0.2U$ and $\lambda_{\text{SOC}}=10$ meV. Note that for the triplet
$p_{z}$-wave pairing, the $S_{z}=\pm 1$ pairing channels are degenerate, which
nevertheless are non-degenerate with the $S_{z}=0$ channel.

The pairing symmetry of this model is obtained through the standard multi-
orbital RPA approach considering SOC[27]. The procedure of this approach is
briefly introduced as follow. Firstly, we need to calculate the susceptibility
tensor as function of the momentum. As the Hamiltonian (4) conserves the SUS,
the spin indices of the nonzero susceptibility tensor elements can only be
$\uparrow\uparrow\uparrow\uparrow$, $\uparrow\uparrow\downarrow\downarrow$,
$\downarrow\downarrow\uparrow\uparrow$,
$\downarrow\downarrow\downarrow\downarrow$,
$\uparrow\downarrow\downarrow\uparrow$,
$\downarrow\uparrow\uparrow\downarrow$. The renormalized susceptibility tensor
can be obtained from the bare one via the Dyson equation in the RPA level.
Secondly, we need to calculate the effective pairing interaction mediated by
exchanging the particle-hole excitations whose propagators are represented by
these susceptibility tensor elements. Note that with SOC here, the effective
pairing interactions in the three spin channels of the triplet pairing are no
longer degenerate. Finally, we need to solve the linearized mean-field (MF)
gap equation near the SC $T_{c}$ involving the obtained effective pairing
interaction function to obtain all possible pairing eigenvalues $\lambda$ and
corresponding pairing eigenvectors. The pairing eigenvectors corresponding to
the largest pairing eigenvalue $\lambda$ determines the pairing symmetry, and
the $T_{c}$ is related to $\lambda$ via $T_{c}\sim
W_{\text{band}}e^{-1/\lambda}$, with $W_{\text{band}}$ denoting the low-energy
band width.

The doping $x$ dependences of the pairing eigenvalues $\lambda$ for the five
leading pairing symmetries are shown in Fig.3 in the doping regime
$x\in(0.4,1)$ wherein definite experimental evidences for SC have been
detected. Note that for the triplet pairings, due to the SUS, the $S_{z}$ of
the Cooper pair is a good quantum number, which suggests that we can get the
$S_{z}=1$ (i.e. the $\uparrow\uparrow$), $S_{z}=-1$ (i.e. the
$\downarrow\downarrow$) or $S_{z}=0$ (i.e. the
$\uparrow\downarrow+\downarrow\uparrow$) pairing channels in the presence of
the SOC. While the $S_{z}=1$ and $S_{z}=-1$ pairing channels are still
degenerate due to the spin degeneracy, they are not degenerate with the
$S_{z}=0$ pairing channel. As a result, Fig.3 shows that the triplet
$p_{z}$-wave pairing in the $S_{z}=0$ (i.e. the
$\uparrow\downarrow+\downarrow\uparrow$) channel is the leading pairing
symmetry. This pairing phase diagram can be compared with that without SOC
provided in Ref[47].

##  5 SUS-Protected Topological Nodal-Line SC

The BCS-mean-field (MF) Hamiltonian describing the triplet $p_{z}$-wave
pairing in the $S_{z}=0$ channel reads,

$\displaystyle H_{\text{BCS}}$
$\displaystyle=H_{\text{TB}}+H_{\text{SOC}}+H_{\Delta},$ (5) $\displaystyle
H_{\Delta}$
$\displaystyle=\sum_{\mathbf{k}\alpha}\left[\left(c^{\dagger}_{\mathbf{k}\alpha\uparrow}c^{\dagger}_{\mathbf{-k}\alpha\downarrow}+c^{\dagger}_{\mathbf{k}\alpha\downarrow}c^{\dagger}_{\mathbf{-k}\alpha\uparrow}\right)\Delta_{\mathbf{k}}^{\alpha}+h.c.\right].$

Here the real gap function $\Delta_{\mathbf{k}}^{\alpha}$ obtained from our
RPA calculations satisfies the $p_{z}$-wave pairing symmetry. This symmetry
requires the gap function to be zero on the $k_{z}=0$ plane. Therefore, the FS
cut on this plane gives the gap nodal lines. It’s interesting to investigate
the topological properties of such nodal-line SC.

It’s easily checked that the MF Hamiltonian (5) satisfies the combined TRS and
SUS. The topological classification of nodal-line superconductors hosting the
two symmetries can be attributed to the following momentum-dependent
topological invariant $Z(k_{x},k_{y})$ defined as[62]

$Z(k_{x},k_{y})=\frac{1}{2\pi}\sum_{\alpha}\int^{\pi}_{-\pi}d\theta_{\mathbf{k}\alpha}\equiv\sum_{\alpha}I_{\alpha}(k_{x},k_{y}).$
(6)

Here $\theta_{\mathbf{k}\alpha}$ denotes the complex phase angle of
$\varepsilon_{\mathbf{k}}^{\alpha}+i\Delta_{\mathbf{k}}^{\alpha}$, and the
integer
$I_{\alpha}(k_{x},k_{y})\equiv\frac{1}{2\pi}\int^{\pi}_{-\pi}d\theta_{\mathbf{k}\alpha}$
represents the winding number of $\theta_{\mathbf{k}\alpha}$ along a closed
path in the BZ which passes through $(k_{x},k_{y},0)$ and is perpendicular to
the $k_{z}=0$ plane. As topological invariant, $2Z(k_{x},k_{y})$ counts the
number of gapless surface modes at the momentum $(k_{x},k_{y})$ in the
$(001)$-surface BZ.

Figure 4: (Color online) (a) The $\gamma$-FS for one spin species considering
SOC. The paths $L_{i,i=1,2}$ are perpendicular to the $(k_{x},k_{y})$-plane.
The momenta $A_{i}\sim E_{i}$ locate on the path $L_{i}$. (b) The $k_{z}=0$
cut of the $\gamma$-FS. The path $L_{1}$ ($L_{2}$) intersects with the
$k_{z}=0$ plane at the momentum $(k_{x1},k_{y1})$ ($(k_{x2},k_{y2})$),
indicated by the center of the small blue (red) circle. The high-symmetry
lines are marked. (c) The $k_{z}$-dependence of the relative gap function
averaged on the FSs. Schematic diagrams of how the phase angle of
$\varepsilon_{\mathbf{k}\gamma}+i\Delta_{\mathbf{k}}^{\gamma}$ evolve with
$k_{z}$ along the paths $L_{1}$ for (d) and $L_{2}$ for (e). The doping level
for (a)-(e) is $x=0.85$.

As an example, we illustrate how to calculate $I_{\gamma}(k_{x},k_{y})$ for
two typical $(k_{x},k_{y})$ points at the doping level $x=0.85$. The
$\gamma$\- FS for $x=0.85$ is shown in Fig. 4(a). This FS consists of two Q1D
FS sheets nearly parallel to the $k_{z}=0$ plane, connected by six hollow
tubes (only three are visible in the figure) nearly perpendicular to the
$k_{z}=0$ plane. The six hollow tubes cut the $k_{z}=0$ plane to form six
pockets, as shown in Fig. 4(b). The $(k_{x1},k_{y1})$ locates near the
$\Gamma$\- point, passed through by the closed vertical path $L_{1}$, and the
$(k_{x2},k_{y2})$ locates within one of the six pockets, passed through by the
path $L_{2}$. The $k_{z}$-dependence of $\Delta^{\alpha}_{\mathbf{k}}$
averaged on the FSs is shown in Fig. 4(c), where the sign of
$\Delta^{\alpha}_{\mathbf{k}}$ follows that of $k_{z}$. This character is also
satisfied for the $\Delta^{\alpha}_{\mathbf{k}}$ with fixed $(k_{x},k_{y})$,
as the Q1D $p_{z}$-wave pairing gap function value doesn’t obviously depend on
$(k_{x},k_{y})$.

Let’s first evaluate $I_{\gamma}$ for $(k_{x1},k_{y1})$. Fig. 4(a) shows the
locations of five typical momenta $A_{1}\sim E_{1}$ along $L_{1}$, which
intersects with the FS twice. The energies $\varepsilon_{\mathbf{k}}^{\gamma}$
for $C_{1}$ and $D_{1}$ are positive, and those for the remaining three are
negative. The gap functions $\Delta_{\mathbf{k}}^{\gamma}$ for $D_{1}$ and
$E_{1}$ are positive, and those for $B_{1}$ and $C_{1}$ are negative. Due to
the $p_{z}$-wave symmetry, the $\Delta_{\mathbf{k}}^{\gamma}$ for $A_{1}$ is
zero. Fig. 4(d) shows the locations of the
$\varepsilon_{\mathbf{k}}^{\gamma}+i\Delta_{\mathbf{k}}^{\gamma}$ for the five
momenta $A_{1}\sim E_{1}$ in the complex plane expanded by the
$\varepsilon_{\mathbf{k}}^{\gamma}$\- axis and the
$\Delta_{\mathbf{k}}^{\gamma}$\- axis, which form a closed trajectory circling
around the original point, suggesting that the phase winding number of
$\varepsilon_{\mathbf{k}}^{\gamma}+i\Delta_{\mathbf{k}}^{\gamma}$ along
$L_{1}$ is 1. Therefore, we have $I_{\gamma}(k_{x1},k_{y1})=1$. Then let’s
investigate $I_{\gamma}(k_{x2},k_{y2})$. As $L_{2}$ locates within one of the
six hollow tubes, it does not intersect with the $\gamma$-FS. Consequently,
the $\varepsilon_{\mathbf{k}}^{\gamma}$ for the five momenta $A_{2}\sim E_{2}$
are all negative. On the other hand, the signs for the
$\Delta_{\mathbf{k}}^{\gamma}$ for the five momenta are the same as those for
$A_{1}\sim E_{1}$. As a result, the trajectory formed by the locations of the
$\varepsilon_{\mathbf{k}}^{\gamma}+i\Delta_{\mathbf{k}}^{\gamma}$ for
$A_{2}\sim E_{2}$ in the complex plane doesn’t circle around the original
point, as shown in Fig. 4(e), which leads to zero winding number of the phase
angle of $\varepsilon_{\mathbf{k}}^{\gamma}+i\Delta_{\mathbf{k}}^{\gamma}$
along $L_{2}$. Thus, we have $I_{\gamma}(k_{x2},k_{y2})=0$.

Figure 5: (Color online) (a) The energy spectrum for KH0.85Cr3As3 along the
high-symmetry lines in the $(k_{x},k_{y})$ plane, with open boundary condition
along the $z$-axis and periodic ones along the $x$\- and $y$\- axes. The
segment marked red (blue) is covered by 4 (6) flat bands. To enhance the
visibility, we have enlarged the gap amplitude by an order of magnitude.

Besides the $\gamma$\- FS, the compound with $x=0.85$ hosts two other Q1D FSs,
i.e. the $\alpha$\- and $\beta$\- FSs. These two FSs take simple shapes, with
each hosting two Q1D FS sheets almost parallel to the $k_{z}=0$ plane, similar
as those for $x=1$ shown in Fig. 1 (d) and Fig. 2 (b). For arbitrary
$(k_{x},k_{y})$, the vertical path $L$ passing through it always intersects
with the $\alpha$\- and the $\beta$\- FS twice. On the other hand, the gap
functions for these two bands are similar with that for the $\gamma$ band in
that the sign of $\Delta_{\mathbf{k}}^{\alpha/\beta}$ follows that of $k_{z}$.
Following the above illustrated procedure, we can easily obtain
$I_{\alpha/\beta}(k_{x},k_{y})=1$ for an arbitrary $(k_{x},k_{y})$. The total
topological invariant $Z(k_{x},k_{y})$ is the sum of
$I_{\alpha}(k_{x},k_{y})$, $I_{\beta}(k_{x},k_{y})$ and
$I_{\gamma}(k_{x},k_{y})$, which takes different nonzero values in different
regimes of the $(k_{x},k_{y})$ plane. For example, for the two typical momenta
shown in Fig. 4 (b), we have $Z(k_{x1},k_{y1})=3$ and $Z(k_{x2},k_{y2})=2$.

Figure 6: (Color online) (a1)-(c1) FSs for KHxCr3As3 with $x$= 1, 0.61 and
0.45. (a2)-(c2) The corresponding surface spectra for the $(001)$ surface. The
segment marked red, blue and green represents 4,6 and 8 flat bands
respectively.

The nontrivial topological invariant $Z(k_{x},k_{y})$ leads to exactly-flat
surface bands on the $(001)$ surface[62]. Concretely, for nonzero
$Z(k_{x},k_{y})$, there would be $2Z(k_{x},k_{y})$ surface bands which are
exactly flat when they pass through the momentum $(k_{x},k_{y})$ in the
surface BZ[62]. Such knowledge guides us to study the surface spectrum on the
$(001)$ surface. For this purpose, we adopt open boundary condition along the
$z$\- axis and periodic ones along the $x$\- and $y$\- axes, under which we
Fourier transform the BCS-MF Hamiltonian (5) to the
$\left(k_{x},k_{y},i_{z}\right)$ space. Then we diagonalize the Bogoliubov-de-
Genes equation for the fixed $(k_{x},k_{y})$ to get the corresponding energy
spectrum. The obtained energy spectrum as function of $k_{x}$ and $k_{y}$ is
shown in Fig. 5 along the high symmetric lines in the $(k_{x},k_{y})$ plane.
To enhance visibility, we have enlarged the gap amplitude by an order of
magnitude. As shown in Fig. 5, some regime in the $(k_{x},k_{y})$-plane is
covered by six flat bands while the remaining regime is covered by four flat
bands. We have checked that these flat bands are surface bands as the
corresponding wave functions are localized on the $(001)$ surface. Comparing
Fig. 5 with Fig. 4 (b), it’s found that the regime within the six Fermi
pockets containing the typical $(k_{x2},k_{y2})$ point wherein
$Z(k_{x},k_{y})=2$ is covered by four flat bands, and the remaining regime
containing the typical $(k_{x1},k_{y1})$ point wherein $Z(k_{x},k_{y})=3$ is
covered by six flat bands.

The FS topology of KHxCr3As3 can be drastically changed upon tuning $x$ which
is experimentally accessible[44]. On the other hand, Fig. 3 shows that the
triplet $p_{z}$-wave pairing in the $S_{z}=0$ channel is always the leading
pairing symmetry, belonging to the SUS-protected TRI topological nodal-line SC
with flat surface bands. Consequently, the distribution of the number of flat
bands in the surface BZ will be easily engineered through tuning the hydrogen
doping. In Fig. 6, the FSs and surface spectra for $x=1$, $x=0.61$ and
$x=0.45$ are shown. Clearly, the different FS topologies for the three
different hydrogen-doping levels lead to different distributions of the number
of flat bands in the surface BZ. For all the hydrogen-doping levels we
studied, the whole surface BZ is covered with flat bands. Particularly in the
case of $x=0.61$, due to the simple FS shape with three Q1D FSs, the whole
surface BZ is covered with six flat bands.

##  6 Discussion and Conclusion

The SUS-protected TRI topological nodal-line SC with flat surface bands
predicted here for the AHxCr3As3 family is similar with that predicted for the
A2Cr3As3 family, which can be detected in the surface conductance spectrum in
the scanning tunneling microscope (STM) experiment by the sharp zero-bias
peak[62]. There can be inter-site spin-flipping SOC in the two families, which
however should be much weaker than the on-site SOC adopted here. In the
presence of such weak spin-flipping SOC, the SUS would be slightly broken,
leading to slightly broadening of the zero-bias peak in the STM, which is
still well detectable[62].

In conclusion, we have performed a RPA-based study on the pairing symmetry and
topological properties of the KHxCr3As3, based on our DFT band structure
considering SOC. Our results yield triplet $p_{z}$-wave pairing in the
$(\uparrow\downarrow+\downarrow\uparrow)$ channel to be the leading pairing
symmetry. This pairing state belongs to the SUS-protected TRI topological
nodal-line SC hosting flat surface bands. Determined by the
$\mathbf{k}$-dependent topological invariant $Z(k_{x},k_{y})$, the whole
$(001)$-surface BZ is covered with topological flat bands, with different
regimes covered with different numbers of flat surface bands, as shown by our
surface spectra, which can be verified by experiments.

###### Acknowledgements.

We appreciate the help from Li-Da Zhang, Chen Lu and Yu-Bo Liu. This work is
supported by the NSFC under the Grant No. 12074031 and 11674025.

## References

  * [1] Mu Q. G., Ruan B. B., Pan B. J., Liu T., Yu J., Zhao K., Chen G. F. Ren Z. A. Phys. Rev. Materials22018034803\.
  * [2] Bao J. K, Liu J. Y., Ma C. W., Meng Z. H., Tang Z. T., Sun Y. L., Zhai H. F., Jiang H., Bai H., Feng C. M., Xu Z. A. Cao G. H. Phys. Rev. X52015011013\.
  * [3] Pang G. M., Smidman M., Jiang W. B., Bao J. K., Weng Z. F., Wang Y. F., Jiao L., Zhang J. L., Cao G. H. Yuan H. Q Phys. Rev. B912015220502(R).
  * [4] Zhi H. Z., Imai T., Ning F. L., Bao J. K. Cao G. H. Phys. Rev. Lett1142015147004\.
  * [5] Adroja D. T., Bhattacharyya A., Telling M., Feng Y., Smidman M., Pan B., Zhao J., Hillier A. D., Pratt F. L. Strydom A. M. Phys. Rev. B922015134505\.
  * [6] Kong T., Budko S. L. Canfield P. C. Phys. Rev. B912015020507(R).
  * [7] Balakirev F. F., Kong T., Jaime M., McDonald R. D., Mielke C. H., Gurevich A., Canfield P. C. Budko S. L. Phys. Rev. B912015220505\.
  * [8] Wang Z., Yi W. Wu Q. arXiv:1603.01688 (2016).
  * [9] Zhao K., Mu Q. G. Liu T. arXiv:1805.11577 (2018).
  * [10] Taddei K. M., Xing G. Sun J. Phys. Rev. L1212018(18): 187002.
  * [11] Bhattacharyya A., Adroja D. T. Feng Y. arXiv:2111.05492 (2021).
  * [12] Wang Z., Yi W. Wu Q. Scientific reports62016(1): 1-8.
  * [13] Taddei K. M., Zheng Q., Sefat A. S. Cruz C. Phys. Rev. B962017180506(R).
  * [14] Zhang L. D., Wu X. X., Fan H., Yang F. Hu J. P. Europhys. Lett113201637003\.
  * [15] Jiang H., Cao G. H. Cao C. Sci. Rep5201516054\.
  * [16] Miao J. J., Zhang F. C. Zhou Y. Phys. Rev. B942016205129\.
  * [17] Zhong H., Feng X. Y., Chen H. Dai J. Phys. Rev. Lett1152015227001\.
  * [18] Yang J., Luo J., Yi C. J., Shi Y. G., Zhou Y. Zheng G. Q. arXiv:2109.01336v1.
  * [19] Pang G., Smidman M., Jiang W., Shi Y., Bao J., Tang Z., Weng Z., Wang Y., Jiao L. Zhang J. J. Magn. Magn. Mater400201684-87.
  * [20] Tang Z. T., Bao J. K., Liu Y., Sun Y. L., Ablimit A., Zhai H. F., Jiang H., Feng C. M., Xu Z. A. Cao G. H. Phys. Rev. B912015020506(R).
  * [21] Yang J., Tang Z .T., Cao G. H. Zheng G. Q. Phys. Rev. Lett1152015147002\.
  * [22] Tang Z. T., Bao J. K., Wang Z., Bai H., Jiang H., Liu Y., Zhai H. F., Feng C. M., Xu Z. A. Cao G. H. Sci China Mater58201516-20.
  * [23] Adroja D. T., Bhattacharyya A., Smidman M., Hillier A. D., Feng Y., Pan B., Zhao J., Lees M. R., Strydom A. M., Biswas P. K. J. Phys. Soc. Jpn862017044710\.
  * [24] Wang X. F., Roncaioli c., Eckberg C., Kim H., Yong J., Nakajima Y., Saha S. R., Zavalij P. Y. Paglione J. Phys. Rev. B922015020508(R).
  * [25] Luo J., Yang J., Zhou R., Mu Q. G., Liu T., Ren Z. A., Yi C. J., Shi Y. G. Zheng G. Q. Phys. Rev.Lett1232019047001\.
  * [26] Wu X., Le C., Yuan J., Fan H. Hu J. Chin. Phys. Lett322015057401\.
  * [27] Wu X. X., Yang F., Le C. C., Fan H. Hu J. P. Phys. Rev. B922015104511\.
  * [28] Zhou Y., Cao C. Zhang F. C. Sci. Bull622017208\.
  * [29] Wu X., Yang F., Qin S., Fan H. Hu J. arXiv:1507.07451.
  * [30] Liang Q., Liu T., Xi C. Y., Han Y. Y., Mu G., Pi L., Ren Z. A. Wang Z. S. Phys. Rev. B1002019(21): 214512.
  * [31] Liu Z. X., Chen M. Y., Xiang Y., Chen X. Y., Yang H., Liu T., Ren Z. A. Wen H. H. Phys. Rev. B1002019(9): 094511
  * [32] Cuono G., Forte F. Romano A. Phys. Rev. B1032021(21): 214406.
  * [33] Bao J. K., Li L., Tang Z. T., Liu Y., Li Y. K., Bai H., Feng C.M., Xu Z. A. Cao G. H. Phys. Rev. B912015180404(R).
  * [34] Tang Z. T., Bao J. K., Liu Y., Bai H., Jiang H., Zhai H. F., Feng C. M., Xu Z. A. Cao G. H. Sci. China Mater582015543\.
  * [35] Mu Q. G., Ruan B. B., Pan B. J., Liu T., Yu J., Zhao K., Chen G. F. Ren Z. A. Phys. Rev. B962017140504(R).
  * [36] Liu T., Mu Q. G., Pan B.J., Yu J., Ruan B. B., Zhao K., Chen G. F. Ren Z. A. Europhys. Lett120201727006\.
  * [37] Cao C., Jiang H., Feng X. Y. Dai J. H. Phys. Rev. B922015235107\.
  * [38] Galluzzi A., Cuono G. Romano A. arXiv:2107.05032 (2021).
  * [39] Feng Y., Zhang X., Hao Y., Hillier A. D., Adroja D. T. Zhao J. Phys. Rev. B992019174401\.
  * [40] Cao G. H., Bao J. K., Tang Z. T., Liu Y. Jiang H. Philos. Mag972017591\.
  * [41] Chen R. Y. Wang N. L. Rep. Prog. Phys822018012503\.
  * [42] Xing G., Shang L. Fu Y. Phys. Rev. B992019(17): 174508.
  * [43] Taddei K. M., Sanjeewa L. D., Lei B. H., Fu Y. H., Zheng Q., Singh D. J., Sefat A. S. Cruz C. D. Phys. Rev. B1002019220503(R).
  * [44] Wu S. Q., Cao C. Cao G. H. Phys. Rev. B1002019155108\.
  * [45] Xiang J. J., Yu Y. L., Wu S. Q., Li B. Z., Shao Y. T., Tang Z. T., Bao J. K. Cao G. H. Phys. Rev. M32019114802\.
  * [46] Xiang J. J., Shao Y. T., Cui Y. W., Nie L. P., Wu S. Q., Li B. Z., Ren Z., Wu T. Cao G. H. Phys. Rev. M42020124802\.
  * [47] Zhang M., Hao J. J., Wu X. X. Yang F. arXiv:2111.14367 (2021).
  * [48] Yada K. Kontani H. J. Phys. Soc. Jpn7420052161\.
  * [49] Kuroki K., Onari S., Arita R., Usui H., Tanaka Y., Kontani H. Aoki H. Phys. Rev. Lett1012008087004\.
  * [50] Kubo K. Phys. Rev. B752007224509\.
  * [51] Graser S., Maier T. A., Hirschfeld P. J. Scalapino D. J. New J. Phys112009025016\.
  * [52] Maier T. A., Graser S., Hirschfeld P. J. Scalapino D. J. Phys. Rev. B832011100515(R).
  * [53] Liu F., Liu C. C., Wu K., Yang F. Yao Y. Phys. Rev. Lett1112013066804\.
  * [54] Liu C. C., Zhang L. D., Chen W. Q. Yang F. Phys. Rev. Lett1212018217001\.
  * [55] Takimoto T., Hotta T. Ueda K. Phys. Rev. B692004104504\.
  * [56] Zhang L. D., Zhang X. M., Hao J. J., Huang W. Yang F. Phys. Rev. B992019094511
  * [57] Yao H. Yang F. Phys. Rev. B922015035132\.
  * [58] Ma T. X., Yang F., Yao H. Lin H. Q. Phys. Rev. B902014245114\.
  * [59] Meng Z. Y., Yang F., Chen K. S., Yao H. Kee H. Y. Phys. Rev. B912015184509\.
  * [60] Chen X., Yao Y. G., Yao H., Yang F. Ni J. Phys. Rev. B922015174503\.
  * [61] Zhang L. D., Yang F. Yao Y. G. Phys. Rev. B922015104504\.
  * [62] Liu C. C., Lu C., Zhang L. D., Wu X. X., Fang C. Yang F. Phys. Rev. Research22020033050\.
  * [63] Perdew J. P., Burke K. Ernzerhof M. Phys. Rev. Lett7719963865\.
  * [64] Giannozzi P., Andreussi O., Brumme T., Bunau O., Nardelli M. B., Calandra M. Baroni S. Journal of Physics: Condensed Matter292017465901\.
  * [65] Mostofi A., Yates J., Pizzi G., Lee Y., Souza I., Vanderbilt D. Marzari N. Comput. Phys. Commun18520142309\.