\section{A time-fractional  PDE} 
In this paper, we consider the following time-fractional parabolic equation for $u(\bm x, t)$ :
\begin{align}\label{Fractional_PDEs}
\begin{split}
\partial_{t}^{\alpha}  u(\bm {x}, t)-\Delta u(\bm {x}, t)&=f(\bm{x}, t), \qquad  (\bm{x}, t) \in \Omega \times(0, T], \\
u(\bm{x}, t)&=0,  \qquad \qquad~(\bm{x}, t) \in \partial \Omega \times(0, T], \\
u(\bm{x}, 0)&=u_0(\bm{x}), \qquad ~\bm{x} \in \Omega, 
\end{split}
\end{align}
where $\Omega$ is a bounded convex polygonal domain in $\mathbb{R}^{d}~ (d=2,3)$ with a boundary $\partial \Omega$ and $u_0$ is a given function defined on the domain $\Omega$ and $T>0$ is a fixed value. Here $\partial_{t}^{\alpha} u~(0<\alpha<1)$ denotes the left-sided Caputo fractional derivative of order $\alpha$ with respect to $t$ and it is defined by 
$$
\partial_{t}^{\alpha} u(\bm x, t)=\frac{1}{\Gamma(1-\alpha)} \int_{0}^{t}(t-s)^{-\alpha} \frac{\mathrm{d}}{\mathrm{d} s} u(\bm x, s) \mathrm{d} s,
$$
where $\Gamma(\cdot)$ is Euler's Gamma function defined by $\Gamma(x)=\int_{0}^{\infty} s^{x-1} \mathrm{e}^{-s} \mathrm{~d} s$ for $x>0 .$


Let $\mathcal{T}_{h}$ be a shape regular and quasiuniform triangulation of domain $\Omega$ and define a continuous piecewise linear finite element space $V_{h}$ by
$$
V_{h}=\left\{v_{h} \in H_{0}^{1}(\Omega):\left.v_{h}\right|_{K} \text { is a linear function, } \forall K \in \mathcal{T}_{h}\right\}.
$$
Then the semidiscrete Galerkin scheme for problem \eqref{Fractional_PDEs} reads: find $u_{h}(t) \in V_{h}$ such that
$$
\left(\partial_{t}^{\alpha} u_{h}, \chi\right)+\left(\nabla u_{h}, \nabla \chi\right)=(f, \chi) \quad \forall \chi \in V_{h}, t>0.
$$

For the time discretization, we divide the interval $[0, T]$ into $N$ equally spaced subintervals with a time step size $\Delta t=T / N$, and $t_{n}=n \Delta t, n=0, \ldots, N$. Then the $L_1$ scheme \cite{Lin_Xu_JCP2007,Li_SISC2011} approximates the Caputo fractional derivative $\partial_{t}^{\alpha} u\left(\bm x, t_{n}\right)$ by
$$
\begin{aligned}
\partial_{t}^{\alpha} u\left(\bm x, t_{n}\right) &=\frac{1}{\Gamma(1-\alpha)} \sum_{j=0}^{n-1} \int_{t_{j}}^{t_{j+1}} \frac{\partial u(\bm x, s)}{\partial s}\left(t_{n}-s\right)^{-\alpha} \mathrm{d} s \\
& \approx \frac{1}{\Gamma(1-\alpha)} \sum_{j=0}^{n-1} \frac{u\left(\bm x, t_{j+1}\right)-u\left(\bm x, t_{j}\right)}{\Delta t} \int_{t_{j}}^{t_{j+1}}\left(t_{n}-s\right)^{-\alpha} \mathrm{d} s \\
&=\sum_{j=0}^{n-1} \beta_{j} \frac{u\left(\bm x, t_{n-j}\right)-u\left(\bm x, t_{n-j-1}\right)}{(\Delta t)^{\alpha} \Gamma(2-\alpha)}=: \bar{\partial}_{\Delta t}^{\alpha} u\left(t_{n}\right),
\end{aligned}
$$
where the weights $\left\{\beta_{j}\right\}$ are given by
$$
\beta_{j}=(j+1)^{1-\alpha}-j^{1-\alpha}, \quad j=0,1, \ldots, n-1.
$$

Then the fully discrete scheme reads: given $u_{h}^{0}= \mathcal I_h u_0 \in V_{h}$ with $c_{\alpha}=\Gamma(2-\alpha)$, find $u_{h}^{n} \in V_{h}$ for $n=1,2, \ldots, N$ such that
\begin{align*}
&\beta_{0}\left(u_{h}^{n}, v_h\right)+c_{\alpha} (\Delta t)^\alpha\left(\nabla u_h^{n}, \nabla v_h\right)\\
&=\beta_{n-1}\left(u_h^{0}, v_h\right)+\sum_{j=1}^{n-1}\left(\beta_{j-1}-\beta_{j}\right)\left(u_h^{n-j}, v_h\right)+c_{\alpha} (\Delta t)^{\alpha}\left(f\left(t_{n}\right), v_h\right).
\end{align*}
Here $\mathcal I_h$ is the standard $L^2$ projection from $L^2(\Omega)$ into $V_h$. The implementationof the scheme is given in \Cref{fractionalPDEs}, and we use $M$ and $A$  to denote  the standard mass and stiffness matrices. 
\begin{algorithm}[h]
	\caption{ (Finite element and $L_1$ scheme for solving time fractional heat equation)}
	\label{fractionalPDEs}
	{\bf Input:} $c_\alpha$, $\alpha$, $\Delta t$, $M\in \mathbb R^{m\times m}$, $A\in \mathbb R^{m\times m}$, $u_0$
	\begin{algorithmic}[1]
		\State  Set $\widetilde u_0  = Mu_0$; $\widetilde A = M + c_{\alpha} (\Delta t)^{\alpha} A$; $\bm \beta = \texttt{zeros}(m,1)$; $\widetilde U = \texttt{zeros}(m,n)$
	
	
	
		\For {$i=1$ to $m$} 
		\State $\bm \beta (i) = i^{1-\alpha} - (i-1)^{1-\alpha}$; $\bm \beta(1)=1$;
		\State Get the load vector at $t_i$, denoted by $b_i$;
		\State $b_i =\bm  \beta(i)\widetilde u_0 + c_{\alpha} (\Delta t)^{\alpha}  b_i + \sum_{j=1}^{i-1} (\bm\beta(j) - \bm \beta(j+1))\widetilde U(:,i-j)$;
		\State Solve $\widetilde Au_i = b_i$;
		\State $\widetilde U(:,i) = Mu_i$;
		\EndFor
	\end{algorithmic}
	{\bf Output:} $u_m$
\end{algorithm}

The computational cost of \Cref{fractionalPDEs} is obvious: To compute the numerical solution at $t_{n}$, the solutions  at all preceding time instances are required, which poses a significant challenge especially for high-dimensional problems and multi-query applications. 





\section{Incremental SVD}
In 2002, Brand \cite{Brand} proposed an incremental  algorithm to find the SVD of a low rank matrix.  
However, the algorithm needs to compute  thousands or millions of orthogonal matrices and their
multiplications together. Those multiplications may corrode the orthogonality and  many reorthogonalizations are needed in practice.  

Very recently, Zhang \cite{Zhang2022isvd} improved the incremental SVD and could avoid computing those large amount  orthogonal matrices. Hence, the expensive  reorthogonalizations are not needed. Furtheremore, it is proved that the modification does not change the accuarcy of the algorithm. 

Let us start by introducing some notation needed throughout this work. For convenience, we adopt 
the Matlab notation herein. Given a matrix $U$, we use $U(:, 1:r)$ and $U_r$ to denote the first $r$ columns of $U$, and $U(1:r,1:r)$ be the leading principal minor of $U$ of order $r$. 


Next, we follow \cite{Zhang2022isvd} and  briefly discuss the incremental  SVD algorithm.

\subsection{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  u_1)^{1/2}, \quad Q=u_1 \Sigma^{-1}, \quad R=1.
\end{align*}
The algorithm is shown in \Cref{algorithm1}.
\begin{algorithm}[H]
	\caption{ (Initialize  ISVD)}
	\label{algorithm1}
	{\bf Input:} $ u_1\in \mathbb R^m$ 
	\begin{algorithmic}[1]
		\State $\Sigma = \left(u_1^\top  u_1\right)^{1/2}$;  $Q = u_1\Sigma^{-1} $;  $R=1$;
	\end{algorithmic}
	{\bf Output:} $Q,\Sigma,R$.
\end{algorithm}

Suppose we already have the rank- $k$ truncated SVD of the first $\ell$ columns of $U$: 
\begin{align}\label{SVD01}
U_\ell = Q \Sigma R^{\top},\quad \textup{with}\;\; 	Q^{\top}   Q=I, \quad R^{\top} R=I, \quad \Sigma=\operatorname{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$.

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  $s$ vectors, $u_{\ell+1}$, $\ldots$, $u_{\ell+s}$,  their residuals are less than the tolerance  when project them onto the  space spanned by the columns of $Q$. While  the residual of $u_{\ell+s+1}$ is greater than the tolerance, i.e.,
\begin{align*}
\|u_{i}- QQ^\top u_{i}\| &<\texttt{tol},\qquad i=\ell+1, \ldots, \ell+s,\\
\|u_{i}- QQ^\top u_{i}\| &\ge \texttt{tol}, \qquad i = \ell+s+1.
\end{align*}

\subsection{Update the SVD of $U_{\ell+s}$}
By \eqref{SVD01} and the above assumption, we have 
\begin{align*}
U_{\ell+s} &=[U_\ell ~\big | ~ u_{\ell+1} ~\big |~ \ldots~\big | ~u_{\ell+s}]\\
& = [Q \Sigma R^{\top}~\big | ~ u_{\ell+1} ~\big |~ \ldots~\big | ~u_{\ell+s}] \\
& \approx [Q \Sigma R^{\top}~\big | ~ QQ^\top u_{\ell+1} ~\big |~ \ldots~\big | ~QQ^\top u_{\ell+s}]\\
& = Q \underbrace{[\Sigma  ~\big | ~ Q^\top u_{\ell+1} ~\big |~ \ldots~\big | ~Q^\top u_{\ell+s}]}_Y  \left[\begin{array}{cc}
R& 0 \\
0 & I_s
\end{array}\right]^\top
\end{align*}

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. We let $Y= Q_Y\Sigma_Y R_Y^\top$ and split $ R_Y$ into $\left[\begin{array}{cc}
R_Y^{(1)}\\
R_Y^{(2)}
\end{array}\right] $,  we then  update the SVD of $U_{\ell+s}$ by 
\begin{align*}
Q\longleftarrow QQ_Y, \quad   \Sigma \longleftarrow  \Sigma_Y,\quad   R \longleftarrow  \left[\begin{array}{cc}
RR_Y^{(1)}\\
R_Y^{(2)}
\end{array}\right]  \in \mathbb R^{k\times (k+s)}.
\end{align*}




\subsection{Update the SVD of $U_{\ell+s+1}$} The SVD of $U_{\ell+s+1}  $ can be constructed by
\begin{itemize}
	\item[(1)] letting   $e =u_{\ell+s+1} -    Q   Q^\top  u_{\ell+s+1}$ and let $p = \|e\| $ and  $\widetilde e = e/p$,
	\begin{align*}
	U_{\ell+s+1} &=
	\left[ 
	U_{\ell+s} ~ \big|~ u_{\ell+1}\right]=\left[ 
	Q\Sigma R^{\top} ~ \big|~ p\widetilde e +QQ^\top u_{\ell+s+1}\right]\\
	&=\left[ 
	Q ~ \big|~  \widetilde e \right]\underbrace{\left[\begin{array}{cc}
		\Sigma & Q^{\top}u_{\ell+s+1} \\
		0 & p
		\end{array}\right]}_{\bar Y}\left[\begin{array}{cc}
	R & 0 \\
	0 & 1
	\end{array}\right]^{\top},	
	\end{align*}
	
	
	
	\item [(2)] finding the full SVD of $\bar Y = \left[\begin{array}{cc}
	\Sigma &   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
	\item[(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[ 
	Q  ~\big |~ \widetilde  e
	\right] \bar Q\right) \bar \Sigma \left(\left[\begin{array}{cc}
	R & 0 \\
	0 & 1
	\end{array}\right] \bar R \right)^{\top}.
	\end{align}
\end{itemize}

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. Hence, we  need to set another truncation if the last few singular values are less than a  tolerance. 
\begin{lemma}\label{com_Sin}\textup{(\cite[Lemma 5.1]{Zhang2022isvd})}
	Assume that $\Sigma=\textup{\texttt{diag}}\left(\sigma_{1}, \sigma_2, \ldots, \sigma_{k}\right)$ with $\sigma_{1}\ge \sigma_2\ge \ldots \ge \sigma_{k}$, and  $\bar \Sigma=\textup{\texttt{diag}}\left(\mu_{1}, \mu_2, \ldots, \mu_{k+1}\right)$ with $\mu_{1}\ge \mu_2\ge \ldots \ge \mu_{k+1}$. Then we have 
	\begin{subequations}
		\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}
	\end{subequations}
\end{lemma} 

Inequality \eqref{Singular1} implies that the last singular value of $\bar 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 \eqref{Singular2} guarantees us  that \emph{only} the last singular value of $\bar Y$ can possibly less than the tolerance. Hence, we only need to check the last one.

\begin{itemize}
	\item[(i)] If $\Sigma_Y(k+1,k+1)\ge\texttt{tol}\sum_{j=1}^{k+1} \Sigma_Y(j,j)$, 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}
	
	\item[(ii)] If $\Sigma_Y(k+1,k+1)<\texttt{tol}\sum_{j=1}^{k+1} \Sigma_Y(j,j)$, 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} 	
\end{itemize}






\subsection{Reorthogonalization}	
Theoretically, the matrix $\left[ 
Q  ~\big|~ \widetilde e\right] Q_Y$ 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. 

For the reorthogonalization to $\left[ 
Q  ~\big|~ \widetilde e\right] Q_Y$, Oxberry et al. \cite{Oxberry} performed QR and Fareed et al. \cite{Fareed} used modified Gram–Schmidt  if the  inner product between $\widetilde e$ and the first column of $ Q$ is larger than some tolerance. The computational cost can be high if the rank of the matrix is not extremely low. In \cite{Zhang2022isvd}, Zhang suggests  reorthogonalize it recursively: only apply the Gram-Schmidt to  the new  adding vector $\widetilde e$ since the matrix $ Q$ has already been  reorthogonalized in previous steps.


\begin{algorithm}[H]
	\caption{ (Update ISVD)}
	\label{algorithm2}
	{\bf 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$, $\texttt{tol}$,  $V$, $Q_0$, $q$
	\begin{algorithmic}[1]
		\State  Set $d = Q^\top u_{\ell+1}$; $e = u_{\ell+1} - Qd$;   $p = (e^\top  e)^{1/2}$;
		\If  {$p<\texttt{tol}$}
		\State $q = q+1$;
		\State Set $V\{q\} = d$;
		\Else 
		\If{$q>0$}
		\State  Set $Y = \left[
		\Sigma ~\big |~ \texttt{cell2mat}(V)\right]$;
		\State 	 $[Q_Y, \Sigma_Y, R_Y] =  \texttt{svd(Y, {\color{magenta} `econ'})}$;
		\State Set  $Q_0 = Q_0Q_Y$,   $\Sigma = \Sigma_Y$,   $\mathcal R_1 = R_Y(1:k,:)$,  $\mathcal R_2 = R_Y(k+1:{\texttt{end}},:)$, $R =  \left[\begin{array}{c}
		R \mathcal R_1\\
		\mathcal R_2
		\end{array}\right]$;
		\State Set $d = Q_Y^\top d$
		\EndIf 
		
		\State Set  $  e = e/p$;
		\If {$|  e^\top W Q(:,1) |> \texttt{tol}$}
		\State $e = e - Q(Q^\top e)$;  $p_1= (e^\top e)^{1/2}$;  $ e = e/p_1$;
		\EndIf
		
		
		\State  Set $Y = \left[\begin{array}{cc}
		\Sigma & d\\
		0 & p
		\end{array}\right]$;
		
		
		\State $[Q_Y, \Sigma_Y, R_Y] = $ \texttt{svd} $(Y)$;
		\State Set $
		Q_0=  \left[\begin{array}{cc}
		Q_0 & 0 \\
		0 & 1
		\end{array}\right]Q_Y$;
		\If {$\Sigma_Y(k+1,k+1)>\texttt{tol} \sum_{j=1}^{k+1} \Sigma_Y(j,j)$}
		\State  $Q = [Q ~\big|~ e]Q_0,\quad  \Sigma = \Sigma_Y$, \quad $R =  \left[\begin{array}{cc}
		R & 0 \\
		0 & 1
		\end{array}\right] R_Y $, \quad $Q_0 = I_{k+1}$;
		\Else
		\State Set $Q = [Q ~\big|~ e]Q_0(:,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), \; Q_0 = I_{k}$
		\EndIf
		\State $V =$ []; $q = 1$;
		\EndIf
	\end{algorithmic}
	{\bf Output:} $Q$, $\Sigma$, $R$, $V$, $Q_0$, $q$
\end{algorithm}
\begin{remark}
	Some care should be taken to assemble the matrix $Y$ in line 7 of \Cref{algorithm2}. In Matlab, it is extremely expensive if we update the matrix $V$ by 
	\begin{align*}
	V \longleftarrow \left[
	V~ \big|~  d
	\right].
	\end{align*}
	We save the vector $ d$ in a cell array  and then use  $\texttt{cell2mat}$ to convert it into a matrix.
\end{remark}


\subsection{Final check and full implementation}

We note that the output of \Cref{algorithm2}, $V$,  may be  not empty. This implies that the output of \Cref{algorithm2} 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 \Cref{algorithm3}.
\begin{algorithm}[H]
	\caption{ ( ISVD final check)}
	\label{algorithm3}
	{\bf Input:} $Q\in \mathbb R^{m\times k}$, $\Sigma\in \mathbb R^{k\times k}$, $R\in \mathbb R^{\ell\times k}$, $V$, $q$
	\begin{algorithmic}[1]
		\If{$q>0$}
		\State  Set $Y = \left[
		\Sigma ~\big |~ \texttt{cell2mat}(V)\right]$;
		\State 	 $[Q_Y, \Sigma_Y, R_Y] =  \texttt{svd(Y, {\color{magenta} `econ'})}$;
		\State  	Set $Q = QQ_Y$,   $\Sigma = \Sigma_Y$,    $\mathcal R_1 = R_Y(1:k,:)$,   $\mathcal R_2 = R_Y(k+1:{\texttt{end}},:)$, $R =  \left[\begin{array}{c}
		R \mathcal R_1\\
		\mathcal R_2
		\end{array}\right]$;
		\EndIf 
	\end{algorithmic}
	{\bf Output:} $Q$, $\Sigma$, $R$.
\end{algorithm}

Next,  we complete the full implementation in \Cref{algorithm4}.

\begin{algorithm}[H]
	\caption{ (Full ISVD)}
	\label{algorithm4}
	{\bf Input:}  $\texttt{tol}$
	\begin{algorithmic}[1]
		\State  Get $u_1$;
		\State $[Q, \Sigma, R]$=\texttt{InitializeISVD}($u_1$);  \hspace{6.5cm} \% \Cref{algorithm1}
		\State Set $V = $ []; $Q_0=1$; $q=0$;
		\For{$\ell=2,\ldots, n$}
		\State Get $u_\ell$
		\State $[q,V, Q_0, Q, \Sigma, R]$=\texttt{UpdateISVD}($q,V, Q_0, Q, \Sigma, R, u_{\ell}, \texttt{tol}$);  \hspace{1.75cm}  \% \Cref{algorithm2}
		\EndFor
		\State $[Q, \Sigma, R]$=\texttt{ISVDfinalcheck}($q,V,Q_0, Q, \Sigma, R$);   \hspace{4.4cm}\% \Cref{algorithm3}
	\end{algorithmic}
	{\bf Output:} $Q$, $\Sigma$, $R$.
\end{algorithm}















\section{Incremental SVD for the time fractional PDEs}
In this section, we shall apply the incremental SVD for the time fractional heat equation \eqref{Fractional_PDEs}. 

Suppose we have already computed $u_1$, $\ldots$, $u_{\ell}$, and  have the rank- $k$ truncated SVD of $U_\ell = [u_1\mid u_2\mid \cdots u_{\ell}]$ by using the incremental SVD algorithm : 
\begin{align}\label{SVD001}
U_\ell = Q \Sigma R^{\top},\quad \textup{with}\;\; 	Q^{\top}   Q=I, \quad R^{\top} R=I, \quad \Sigma=\operatorname{diag}\left(\sigma_{1}, \ldots, \sigma_{k}\right).	
\end{align}

We do not store the data $\{u_i\}_{i=1}^\ell$, instead, we keep the matrices $Q\in \mathbb R^{n\times k}$, $\Sigma\in \mathbb R^{k\times k}$ and $R\in \mathbb R^{\ell\times k}$. Since $k\ll \ell$, then we remedy the memory issue of the standard method. 

To compute $u_{\ell+1}$, we need to extract the data from the matrices $Q$, $\Sigma$ and $R$, i.e.,
\begin{align*}
u_j = Q\Sigma R_j^\top,\quad j=1,2,\ldots, \ell,
\end{align*}
and then compute
\begin{align*}
Mu_j = MQ\Sigma R_j^\top,\quad j=1,2,\ldots, \ell.
\end{align*}
The matrix  $MQ$ stays unchanged unless the residual of the new coming data $u_{\ell+1}$ is greater than the tolerance,  when we project it onto the  space spanned by the columns of $Q$. 


\begin{algorithm}[H]
	\caption{ (Solve a fractional PDE with ISVD)}
	\label{fpdeandisvd}
	{\bf Input:} $c_\alpha$, $\alpha$, $\Delta t$, $m$, $M\in \mathbb R^{m\times m}$, $A\in \mathbb R^{m\times m}$, $u_0$, $u_1$, $\texttt{tol}$
	\begin{algorithmic}[1]
		\State  Set $V =$ []; $ Q_0=1$; $q=0$; $\widetilde u_0  = Mu_0$; $\widetilde A = M + c_{\alpha} (\Delta t)^{\alpha} A$; $\bm \beta = \texttt{zeros}(m,1)$;
		\State $[Q, \Sigma, R]$=\texttt{InitializeISVD}($u_1$);  \hspace{6.5cm} \% \Cref{algorithm1}
		\For {$i=2$ to $m$} 
		\State $\bm \beta (i) = i^{1-\alpha} - (i-1)^{1-\alpha}$;
		\State Get the load vector at $t_i$, denoted by $b_i$;
		\If {$q==0$}
		\State $\texttt{MQ} = MQ$;  $C = \Sigma R^\top$
		\State $b_i =\bm  \beta(i)\widetilde u_0 + c_{\alpha} (\Delta t)^{\alpha}  b_i + \texttt{MQ} \left(\sum_{j=1}^{i-1} (\bm\beta(j) - \bm \beta(j+1)) C(:,i-j)\right)$;
		\State Solve $\widetilde Au_i = b_i$;
		\State 	 $[q,V, Q_0, Q, \Sigma, R]$=\texttt{UpdateISVD}($q,V, Q_0, Q, \Sigma, R, u_{i},\texttt{tol}$);  \hspace{1.2cm}  \% \Cref{algorithm2}
		\If {$q>0$}
		\State $C =  \Sigma R^\top$
		\EndIf
		\Else
		\State $D = [C \quad   \texttt{cell2mat}(V)]$;
		\State $b_i =\bm  \beta(i)\widetilde u_0 + c_{\alpha} (\Delta t)^{\alpha}  b_i + \texttt{MQ} \left(\sum_{j=1}^{i-1}  D(:,i-j) \bm\beta(j)\right)$;
		\State Solve $\widetilde Au_i = b_i$;
		\State 	 $[q,V, Q_0, Q, \Sigma, R]$=\texttt{UpdateISVD}($q,V, Q_0, Q, \Sigma, R, u_{i},\texttt{tol}$);  \hspace{1.2cm}  \% \Cref{algorithm2}
		\EndIf
		\EndFor
		\State $[Q, \Sigma, R]$=\texttt{UpdateISVDcheck}($q,V,Q_0, Q, \Sigma, R$);   \hspace{4.2cm}\% \Cref{algorithm3}
	\end{algorithmic}
	{\bf Input:} $Q$, $ \Sigma$, $R$
\end{algorithm}


\section{Numerical experiments}
In this section, we test the efficiency of our \Cref{fpdeandisvd}. 
\begin{example}\label{Example1}
	Let  $\Omega = (0,1)\times (0,1)$, we consider the equation \eqref{Fractional_PDEs} with
	\begin{align*}
	\alpha = 0.1, \quad T=1,\quad  u_0 = xy(1-x)(1-y), \quad f = 100\sin(2\pi t(x+y))x(1-x)y(1-y).	
	\end{align*}
	We use the linear finite element for the spatial discretization with different  time steps $\Delta t$, here $h$ is	the mesh size (max diameter of the triangles in the mesh). 					
	For the incremental SVD, we take \texttt{tol} $= 10^{-12}$ in \Cref{fpdeandisvd}.  For solving linear systems, we apply the Matlab built-in solver backslash ($\backslash$). We report the wall time (seconds) \footnote{All the code for all examples in the paper has been created by the authors 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. We use the Matlab built-in function \texttt{tic}-\texttt{toc} to denote  the real simulation time.} of \Cref{fractionalPDEs} and \Cref{fpdeandisvd}  in \Cref{table_0}. We also compute the $L^2$-norm error between the solutions of \Cref{fractionalPDEs} and \Cref{fpdeandisvd}  at the final time $T=1$, we see that the error is close to the machine error. Furthermore,  \Cref{fpdeandisvd} is more efficient than \Cref{fractionalPDEs}.  In \Cref{table_1}, we see that the standard solver has an issue to save the data. However, our new  \Cref{fpdeandisvd} avoids storing the huge data.  
	
	\begin{table}[H]
		\centering
		{
			\begin{tabular}{c|c|c|c|c|c|c|c}
				\Xhline{1pt}
				
				\cline{3-8}
				$h$	& $1/2^3$ &$1/2^4$ 
				&$1/2^5$ 
				&$1/2^6$  &$1/2^7$ &$1/2^8$  &$1/2^{9}$ 
				\\
				\cline{1-8}
				{Wall time of \Cref{fractionalPDEs} }
				& 	 0.3&   0.6&   2&   10&   36&   162&   680\\ 
				\cline{1-8}
				{Wall time of \Cref{fpdeandisvd}}
				& 	  0.4&   0.8&   2&  8&  36&   146&   670\\ 
				\cline{1-8}
				{Error}
				&8E-15& 2E-14&	  3E-14&   3E-14&   3E-14&  3E-14&   3E-14\\ 
				\Xhline{1pt}
				
				
			\end{tabular}
		}
		\caption{\Cref{Example1} with $\Delta t = 10^{-3}$: The wall time  and the  $L^2$-norm error between the solutions of  \Cref{fractionalPDEs}  and \Cref{fpdeandisvd}.}\label{table_0}
	\end{table}		
	
	
	
	
	
	\begin{table}[H]
		\centering
		{
			\begin{tabular}{c|c|c|c|c|c|c|c|c}
				\Xhline{1pt}
				
				\cline{3-9}
				$h$	& $1/2^3$ &$1/2^4$ 
				&$1/2^5$ 
				&$1/2^6$  &$1/2^7$ &$1/2^8$  &$1/2^{9}$&$1/2^{10}$
				\\
				\cline{1-9}
				{Wall time of \Cref{fractionalPDEs} }
				& 	 6&   17&   52&  257&   700&   2410&   9519&---\\ 
				\cline{1-9}
				{Wall time of \Cref{fpdeandisvd}}
				& 	  11&   15&   26&  83&  330&   1405&   6424&25835\\ 
				\cline{1-9}
				{Error}
				&1E-12& 2E-14&	  3E-14&   3E-14&   3E-14&  3E-14&   3E-14&---\\ 
				\Xhline{1pt}
				
				
			\end{tabular}
		}
		\caption{ \Cref{Example1} with $\Delta t = 10^{-4}$: The wall time   and the  $L^2$-norm error between the solutions of  \Cref{fractionalPDEs}  and \Cref{fpdeandisvd},  --- denotes we have the memory issue.}\label{table_1}
	\end{table}		
\end{example}

\begin{example}\label{Example2}
	In this example, we use the incremental SVD algorithm to study the  time fractional diffusion-wave equation:
	$$
	\begin{aligned}
	\partial_{t}^{\alpha}  u(\bm {x}, t)-\Delta u(\bm {x}, t)&=f(\bm{x}, t), &(\bm{x}, t) \in \Omega \times(0, T], \\
	u(\bm{x}, t)&=0, &(\bm{x}, t) \in \partial \Omega \times(0, T], \\
	u(\bm{x}, 0)&=u_0(\bm{x}), & \bm{x} \in \Omega, \\
	u_t(\bm{x}, 0)&=v_0(\bm{x}),  & \bm{x} \in \Omega,
	\end{aligned}
	$$
	where $\Omega \subset \mathbb R^2$ is a bounded convex polygonal region with boundary $\partial \Omega$, $\bm {x}=(x, y)$, $u_0(\mathbf{x})$, $v_0(\bm {x})$ and $f(\bm{x}, t)$ are given functions, and $\partial_{t}^{\alpha} u~(1<\alpha<2)$ denotes the left-sided Caputo fractional derivative of order $\alpha$ with respect to $t$ and it is defined by 
	$$
	\partial_{t}^{\alpha}  u(\mathbf{x}, t)=\frac{1}{\Gamma(2-\alpha)} \int_0^t \frac{\partial^2 u(\bm {x}, s)}{\partial s^2} \frac{\mathrm{d} s}{(t-s)^{\alpha-1}}, \quad 1<\alpha<2.
	$$
	
	Let $0=t_0<t_1<\cdots<t_N=T$ be a given partition of the time interval. Then we have the time step $\Delta t =T / N$, $t_n=n \Delta t(n=0,1, \ldots, N)$, and define  $\bar{\partial}_t^\alpha u^{n+1} $ by 
	$$
	\begin{aligned} \frac{\partial^\alpha u\left(\bm{x}, t_{n+1}\right)}{\partial t^\alpha} &=\frac{1}{\Gamma(2-\alpha)} \int_0^{t_{n+1}} \frac{\partial^2 u(\bm{x}, \xi)}{\partial s^2}\left(t_{n+1}-s\right)^{1-\alpha} \mathrm{d} s \\ &=\frac{1}{\Gamma(2-\alpha)} \sum_{k=0}^n \int_{t_k}^{t_{k+1}} \frac{\partial^2 u(\bm {x}, s)}{\partial s^2}\left(t_{n+1}-s\right)^{1-\alpha} \mathrm{d} s\\
	& \approx  \frac{1}{\Gamma(2-\alpha)} \sum_{k=0}^n \frac{u_{k+1}-2 u_{k}+u_{k-1}}{\Delta t^2} \int_{t_k}^{t_{k+1}}\left(t_{n+1}-s\right)^{1-\alpha} \mathrm{d} s,
	\end{aligned}
	$$
	where   $u_k=u\left(\bm{x}, t_k\right), k=0,1,2, \ldots, N$. In the above equation, the integral is easily obtained as
	$$
	\int_{t_k}^{t_{k+1}}\left(t_{n+1}-s\right)^{1-\alpha} \mathrm{d} s=\frac{1}{(2-\alpha)} \Delta t^{2-\alpha}\left[(n-k+1)^{2-\alpha}-(n-k)^{2-\alpha}\right] .
	$$
	Let  $\beta_k=(k+1)^{2-\alpha}- k^{2-\alpha}$  and define $ \bar{\partial}_t^\alpha  u_h^{n+1}$ by 
	$$
	\begin{aligned}
	\bar{\partial}_t^\alpha  u_h^{n+1}  =\frac{\Delta t^{-\alpha}}{\Gamma(3-\alpha)} \sum_{k=0}^n \beta_k\left[u_{n-k+1}-2 u_{n-k}+u_{n-k-1}\right],
	\end{aligned}
	$$
	where $u_{h}^{0}= \mathcal I_h u_0$, $u_h^{-1} = u_h^0-\Delta t  \mathcal I_h v_0$. Then the fully discrete scheme reads: given $u_h^{-1}, u_{h}^{0} \in V_{h}$, find $u_{h}^{n+1} \in V_{h}$ for $n=0,1,2, \ldots$ such that
	\begin{align*}
	\left(\bar{\partial}_t^\alpha u_h^{n+1}, v_h\right)+\left(\nabla u_h^{n+1}, \nabla v_h\right)=\left(f^{n+1}, v_h\right), \quad \forall v_h \in V_h.
	\end{align*}
	
	Next,  we consider the above wave equation  with 
	\begin{align*}
	\Omega = (0,1)^2, \;\; 	\alpha = 1.5, \;\; u_0 = xy(1-x)(1-y),\;\; v_0 = 0,\;\; f=1.
	\end{align*}
	The mesh and the finite element space are the same with \Cref{Example1}.  We report the wall time (seconds)   and  errors at the final time $T=1$ in \Cref{table_3}.  Again, our incremental SVD not only memory free, but also more efficient than the standard solver. 
	
	
	\begin{table}[H]
		\centering
		{
			\begin{tabular}{c|c|c|c|c|c|c|c|c}
				\Xhline{1pt}
				
				\cline{3-9}
				$h$	& $1/2^3$ &$1/2^4$ 
				&$1/2^5$ 
				&$1/2^6$  &$1/2^7$ &$1/2^8$  &$1/2^{9}$&$1/2^{10}$
				\\
				\cline{1-9}
				{Wall time of \Cref{fractionalPDEs} }
				& 	 6&   19&   77&  330&   1002&   3183&   13175&---\\ 
				\cline{1-9}
				{Wall time of \Cref{fpdeandisvd}}
				& 	  24&   63&   111&  293&  862&   2645&   12748&60234\\ 
				\cline{1-9}
				{Error}
				&1E-09& 4E-10&	  4E-10&   3E-10&   3E-10&  3E-10&   6E-10&---\\ 
				\Xhline{1pt}
				
				
			\end{tabular}
		}
		\caption{\Cref{Example2} with $\Delta t = 10^{-4}$: The wall time   and the  $L^2$-norm error between the solutions of  \Cref{fractionalPDEs}  and \Cref{fpdeandisvd},  --- denotes we have the memory issue.}\label{table_3}
	\end{table}		
\end{example}



\section{Conclusion} 
In this paper, we propose an efficient and memory free algorithm for solving the time-fractional PDEs by using the incremental SVD algorithm. Although we only test the fractional heat equation, we beleive this algorithm could be easily to extend to other time-fractional PDEs such as the Cole-Cole Maxwell's equations \cite{Li_SISC2011}. Furthermore, we will investigate time-dependent PDEs with memory effects, since they share the same challenge that a whole past history of the solution is needed to be stored and used. 


