TITLE: What does it mean to normalize a combination of Spherical Harmonics?
QUESTION [0 upvotes]: Using the following as an example:

Show that the combinations $$-\frac{1}{\sqrt{2}}\left(Y_{11}-Y_{1-1}\right)\quad\text{&}\quad\frac{i}{\sqrt{2}}\left(Y_{11}+Y_{1-1}\right)$$ are real and normalized.
Where $$Y_{1\pm 1}=\mp\sqrt{\frac{3}{8\pi}}\sin\theta e^{\pm i\phi}$$

So for the first combination we have $$\begin{align}-\frac{1}{\sqrt{2}}\left(Y_{11}-Y_{1-1}\right)&=-\frac{1}{\sqrt{2}}\left[-\sqrt{\frac{3}{8\pi}}\sin\theta \,e^{i\phi}-\sqrt{\frac{3}{8\pi}}\sin\theta \,e^{-i\phi}\right]\\&=\frac{1}{\sqrt{2}}\left[\sqrt{\frac{3}{8\pi}}\sin\theta \left(e^{i\phi}+e^{-i\phi}\right)\right]\\&=\sqrt{\frac{3}{4\pi}}\sin\theta\cos\phi\\&=\sqrt{\frac{3}{4\pi}}\frac{x}{r}\end{align}$$
where in the last step I used the fact that $x=r\sin\theta\cos\phi$ in spherical polar coordinates.

For the second combination
$$\begin{align}\frac{i}{\sqrt{2}}\left(Y_{11}+Y_{1-1}\right)&=\frac{i}{\sqrt{2}}\left[-\sqrt{\frac{3}{8\pi}}\sin\theta \,e^{i\phi}+\sqrt{\frac{3}{8\pi}}\sin\theta \,e^{-i\phi}\right]\\&=-\frac{i}{\sqrt{2}}\left[\sqrt{\frac{3}{8\pi}}\sin\theta \left(e^{i\phi}-e^{-i\phi}\right)\right]\\&=\sqrt{\frac{3}{4\pi}}\sin\theta\sin\phi\\&=\sqrt{\frac{3}{4\pi}}\frac{y}{r}\end{align}$$
where in the last step I used the fact that $y=r\sin\theta\sin\phi$ in spherical polar coordinates.

The calculation is obviously not the problem here. It is the meaning of the result that concerns me. The author of this question mentions in the solution that:

They are obviously real. The spherical harmonics are orthogonal and normalized, $\color{blue}{\text{so the square integral of the two new functions will just give }}\color{blue}{\frac12(1+1)=1}$.

I know that the word 'normalize' is quite ubiquitous as it seems to have different meanings depending on the context.
For example:

The wavefunction for a particle confined to $0\le x \le a$ in the ground state is $$\Psi(x)=A\sin\left(\frac{\pi\, x}{a}\right)$$ where $A$ is the normalization constant. Find $A$ and hence the normalized wavefunction $\psi(x)$.

Omitting all calculations using the 'Normalization Condition':
$$\int_{x=-\infty}^{\infty}{\big|\Psi(x)\big|}^2\,dx=1\tag{1}$$
The Normalized wavefunction is $$\psi(x)=\sqrt{\frac{2}{a}}\sin\left(\frac{\pi\,x}{a}\right)$$

The words "Normalized wavefunction" make perfect sense to me as it has been normalized such that the condition $(1)$ is satisfied.

But what does it mean to say that $$-\frac{1}{\sqrt{2}}\left(Y_{11}-Y_{1-1}\right)=\sqrt{\frac{3}{4\pi}}\frac{x}{r}$$ is normalized?
Or put in another way; Where is the 'Normalization Condition' for this?
At a guess, I think the answer lies somewhere in what the Professor wrote in the quote above which I have marked blue. But I don't fully understand the logic behind what the professor is saying.

REPLY [3 votes]: There are several choices for the normalization of spherical harmonics, but in your case it seems to be that $Y_{l,m} (\theta,\phi)$ is normalized if
$$
\int_0^\pi \mathrm{d}\theta \int_0^{2\pi} \mathrm{d}\phi \sin(\theta)  \left|Y_{l,m} (\theta,\phi) \right|^2 = 1
$$
In problems like these it is often much easier to use the orthonormality of the spherical harmonic functions,
$$
\int_0^\pi \mathrm{d}\theta \int_0^{2\pi} \mathrm{d}\phi \sin(\theta)  Y_{l,m} (\theta,\phi) Y_{l',m'}^{*} (\theta,\phi) = \delta_{l,l'}\delta_{m,m'}
$$
to test whether linear combinations of spherical harmonic functions are normalized.


text:
\begin{document}
\maketitle
\slugger{mms}{xxxx}{xx}{x}{x--x}

\begin{abstract}
Krylov subspace recycling is a process for accelerating the convergence of sequences of linear systems. Based on this technique, the recycling BiCG algorithm has been developed recently. Here, we now generalize and extend this recycling theory to BiCGSTAB. Recycling BiCG focuses on efficiently solving sequences of dual linear systems, while the focus here is on efficiently solving sequences of single linear systems (assuming non-symmetric matrices for both recycling BiCG and recycling BiCGSTAB).

As compared with other methods for solving sequences of single linear systems with non-symmetric matrices (e.g., recycling variants of GMRES), BiCG based recycling algorithms, like recycling B\-i\-C\-G\-S\-T\-A\-B, have the advantage that they involve a short-term recurrence, and hence, do not suffer from storage issues and are also cheaper with respect to the orthogonalizations.

We modify the BiCGSTAB algorithm to use a recycle space, which is built from left and right approximate invariant subspaces. Using our algorithm for a parametric model order reduction example gives good results. We show about 40\% savings in the number of matrix-vector products and about 35\% savings in runtime.
\end{abstract}

\begin{keywords}
Krylov subspace recycling, deflation, BiCGSTAB, model reduction, rational Krylov.
\end{keywords}

\begin{AMS}
65F10, 65N22, 93A15, 93C05.
\end{AMS}

\pagestyle{myheadings}
\thispagestyle{plain}
\markboth{AHUJA, BENNER, DE STURLER, FENG}{RECYCLING BICGSTAB}

\section{Introduction}
We focus on efficiently solving sequences of linear systems of the following type:
\begin{equation}\label{eq:initial3}

    {A}^{(\iota)}\ x^{(\iota, \kappa)} = {b}^{(\iota, \kappa)},
\end{equation}
where ${A}^{(\iota)} \in \mathbb{R}^{n \times n}$ varies with $\iota$; ${b}^{(\iota, \kappa)} \in \mathbb{R}^n$ varies with both $\iota$ and $\kappa$; the matrices ${A}^{(\iota)}$ are large, sparse, and non-symmetric; and the change from one system to the next is small.

Krylov subspace methods are usually used for solving such large and sparse linear systems. For linear systems with non-symmetric matrices, GMRES~\cite{saad86gmres} is one of the first choices, but it is generally not optimal with respect to the runtime. BiCGSTAB~\cite{vorst92bicgstab}

is competitive with GMRES, and in many cases performs better than GMRES in time. Also, it does not suffer from storage issues, which is a problem in GMRES.

Krylov subspace recycling is a technique for efficient solution of sequences of linear systems. Here, while solving one system in the sequence, approximate invariant subspaces of the matrix are selected and used to accelerate the convergence of the next system in the sequence. Since the matrices in the sequence do not change much, this provides substantial reduction in both the number of matrix-vector products and time. See~\cite{parks06gcrodr} and~\cite{ahuja12rbicg} for more about Krylov subspace recycling.

Here, we have a sequence of linear systems with non-symmetric matrices, and hence, GCRO-DR~\cite{parks06gcrodr} and GCROT~\cite{parks06gcrodr} algorithms, which are recycling variants of GMRES, are more suited.

However, since there is no optimal method in time for solving linear systems with non-symmetric matrices, there is no optimal method in time for solving sequences of such linear systems. Like GMRES, its recycling variants may also suffer from storage issues. Hence, we develop a recycling variant of BiCGSTAB based on the work in~\cite{ahuja11phd, ahuja09ms}.

We explore the usage of recycling B\-i\-C\-G\-S\-T\-A\-B for parametric model order reduction (PMOR)~\cite{baur2011PMOR,feng:proceedings} that requires solution of systems of the form (\ref{eq:initial3}). We show about 40\% reduction in the number of matrix-vector products when using recycling as compared with not using recycling in B\-i\-C\-G\-S\-T\-A\-B. In terms of time, this translates to about 35\% savings in runtime.

In related work in this area (specific to BiCGSTAB),~\cite{gutnket2014bicg} discusses a variant of recycling BiCGSTAB  (using the terminology of deflation and augmentation instead of Krylov subspace recycling). There are three main differences between that approach and ours. First, we use a different inner product in our derivation as compared with~\cite{gutnket2014bicg}. Second, the focus in~\cite{gutnket2014bicg} is on only using a recycle space, while here we discuss both using and generating a recycle space. Third, we also give numerical experiments demonstrating the usefulness of our approach, while~\cite{gutnket2014bicg} discusses only a theoretical framework. Also,~\cite{abdel2007lefteig, morgan2010nlandr, abdel2014rhs} focus on variants of deflated BiCGSTAB for multiple right hand sides and do not discuss changing matrices.

Since we propose recycling BiCGSTAB as an alternative for GMRES-based recycling solvers, we also compare with GCRO-DR. In the context of PMOR, simplified versions of GCRO-DR have also been proposed~\cite{feng:proceedings, feng2012Mems}. We show that our recycling BiCGSTAB is 10\% more efficient in time than GCRO-DR for our test problem.

To simplify notation, we drop the superscripts $\iota$ and $\kappa$ in (\ref{eq:initial3}).

Throughout the paper, $||\cdot||$ refers to the two-norm, $(\cdot, \cdot)$ refers to the standard inner product, $*$ indicates the conjugate transpose operation, $\bar{\cdot}$ indicates complex conjugation, and
$\underline{\cdot}$ is used to signify a rectangular matrix.

The rest of the paper is divided into six more sections. The bi-Lanczos algorithm~\cite{lanczos1952base} and recycling BiCG~\cite{ahuja12rbicg} form the basis of our recycling BiCGSTAB. Hence, we revisit these in Sections \ref{sec:bilanczos} and \ref{sec:rbicg}, respectively. In Section \ref{sec:rbicg}, we also give a new result related to recycling BiCG. Next, we derive recycling BiCGSTAB in Section \ref{sec:rbicgstab}. In Section \ref{sec:analysis}, we analyze the subspaces that can be used in recycling BiCGSTAB. Finally, we discuss the application of recycling BiCGSTAB to PMOR in Section \ref{sec:appl}, and give concluding remarks in Section \ref{sec:conclusion}.

\section{The Bi-Lanczos Algorithm}\label{sec:bilanczos}
Consider a primary system $Ax = b$, with $x_0$ the initial guess and $r_0 = b - A x_0$ the residual. Also, consider an auxiliary dual system $A^*\tilde{x} = \tilde{b}$, with $\tilde{b}$ a random vector, $\tilde{x}_0$ the initial guess, and $\tilde{r}_0 = \tilde{b} - A^* \tilde{x}_0$ the residual. This dual system is termed auxiliary because for this work we are not interested in its solution (although the system is real). The bi-Lanczos algorithm remains the same even when the dual system is of interest.

Let the columns of $V_{i} = [v_{1}\ v_{2}\ \ldots\ v_i]$ define the basis of the primary system Krylov space $\mathcal{K}^{i}(A,r_{0}) \equiv span\{r_{0},\ Ar_{0},\ A^2 r_{0},\ \cdots,\ A^{i-1} r_{0}\}$. Also, let the columns of ${\tilde V}_{i} = [{\tilde v}_{1}\ {\tilde v}_{2}\ \ldots\ \tilde{v}_i]$ define the basis of the dual system Krylov space $\mathcal{\tilde{K}}^{i}(A^*,\tilde{r}_{0}) \equiv span\{\tilde{r}_{0},\ A^*\tilde{r}_{0},\ A^{2*} \tilde{r}_{0},\ \cdots,\ A^{(i-1)*} \tilde{r}_{0}\}$.

The bi-Lanczos algorithm computes the columns of $V_i$ and $\tilde{V}_i$ such that, in exact arithmetic, $V_{i} \perp_{b} \tilde{V}_{i}$, where $\perp_{b}$ is referred to as bi-orthogonality; this implies that $\tilde{V}_i^*V_i$ is a diagonal matrix. The columns of $V_i$ and $\tilde{V}_i$   are called Lanczos vectors. There is a degree of freedom in choosing the scaling of the Lanczos vectors~\cite{greenbaum1997book, gutnket1997survey, saad03book}. Using the scaling
\begin{equation}\label{scalingfactor}
    ||v_{i}|| = 1, \quad \quad (v_i, \tilde{v}_i) = 1,
\end{equation}
we initialize the Lanczos vectors as follows:
\begin{align*}
    \begin{array}
        [l]{ll}
        v_1 = \frac{r_0}{|| r_0 ||}, & \tilde{v}_{1} = \frac{\tilde{r}_{0}}{(v_{1}, \tilde{r}_{0})}.
    \end{array}
\end{align*}
The ($i+1$)-at Lanczos vectors are given by
\begin{align}\label{eq:biLanczosCompute}
    \begin{array}
        [l]{lll}
        \gamma v_{i+1} & = & A v_i - V_i \tau \perp \tilde{V}_i, \\
        \tilde{\gamma}\tilde{v}_{i+1} & = & A^* \tilde{v}_i - \tilde{V_i}\tilde{\tau}  \perp {V}_i,
    \end{array}
\end{align}
where $\gamma$, $\tilde{\gamma}$ and $\tau$, $\tilde{\tau}$ are are determined by the biorthogonality condition (\ref{eq:biLanczosCompute}) and the normalization condition (\ref{scalingfactor}). The computation of the ($i+1$)-st Lanczos vectors requires only the $i$-th and the ($i-1$)-st Lanczos vectors (see~\cite{saad03book}). These $3$-term recurrences are called the bi-Lanczos relations, and are defined as follows:
\begin{align}\label{eq:biLanczos}
    \begin{array}
        [l]{lllll}
            AV_{i} & = & V_{i+1}\underline{T}_{i} & =  & V_{i}{T}_{i} + t_{i+1,i}v_{i+1}e_i^T, \\
            A^*\tilde{V}_{i} & = & \tilde{V}_{i+1}\tilde{\underline{T}}_{i} & = & \tilde{V}_{i}{\tilde{T}}_{i} + \tilde{t}_{i+1,i}\tilde{v}_{i+1}e_i^T, \\
    \end{array}
\end{align}
where $T_i$, $\tilde{T}_i$ are $i \times i$ tridiagonal matrices, $t_{i+1,i}$ is the last element of the last row of $\underline{T}_i \in \mathbb{C}^{(i+1) \times i}$, and $\tilde{t}_{i+1,i}$ is the last element of the last row of $\tilde{\underline{T}}_i \in \mathbb{C}^{(i+1) \times i}$.

The bi-Lanczos algorithm breaks down when at any step $i$, $\tilde{v}_i^*v_i$ = 0. There exist so-called look-ahead strategies~\cite{freund:lookahead, gutnket1997survey}, that take multiple Lanczos vectors together in-succession and make them block bi-orthogonal, to avoid this breakdown.

\section{Recycling BiCG Revisited}\label{sec:rbicg}
We first introduce a generalization of the bi-Lanczos algorithm~\cite{ahuja11phd}. We show that even for a pair of matrices that are not conjugate transposes of each other,

one can build bi-orthogonal bases (for the associated two Krylov subspaces) using a short-term recurrence.

Expanding the search space to include a recycle space leads to an augmented bi-orthogonality condition. The augmented bi-Lanczos algorithm, as derived for recycling BiCG~\cite{ahuja12rbicg}, computes bi-orthogonal bases for the two Krylov subspaces such that this augmented bi-orthogonality condition is satisfied. Next, we revisit augmented bi-Lanczos~\cite{ahuja09ms} and show that it is a special case of generalized bi-Lanczos.  Finally, we list the recycling BiCG algorithm from~\cite{ahuja12rbicg}.

There are numerous ways of computing good bases for Krylov subspaces $\mathcal{K}^{m}(B,v_{1})$ and $\mathcal{K}^{m}(\tilde{B},\tilde{v}_{1})$, where $B$ and $\tilde{B}$ are $n \times n$ general matrices, and $v_1$ and $\tilde{v}_1$ are any two $n$ dimensional vectors. Let the columns of $V_{m} = [v_{1}\ v_{2}\ \ldots\ v_m]$ and $\tilde{V}_{m} = [\tilde{v}_{1}\ \tilde{v}_{2}\ \ldots\ \tilde{v}_m]$ define one such pair of good bases for $\mathcal{K}^{m}(B,v_{1})$ and $\mathcal{K}^{m}(\tilde{B},\tilde{v}_{1})$, respectively. We compute these bases using the following, in principle, full recurrences:
\begin{align}\label{eq:basis1}
    \begin{array}
        [l]{lll}
	\beta_{i+1,i}{v}_{i+1} & = & {B}{v}_i - {\beta}_{ii}{v}_i -  {\beta}_{i-1,i}{v}_{i-1} - \ldots - {\beta}_{1i}{v}_1,
	\end{array}
\end{align}
\begin{align}\label{eq:basis2}
    \begin{array}
        [l]{lll}
	\tilde{\beta}_{i+1,i}\tilde{v}_{i+1} & = & \tilde{B}\tilde{v}_i - \tilde{\beta}_{ii}\tilde{v}_i -  \tilde{\beta}_{i-1,i}\tilde{v}_{i-1} - \ldots - \tilde{\beta}_{1i}\tilde{v}_1,
	\end{array}
\end{align}
where $i \in \{1, 2, 3, \ldots,m-1\}$ and $\{\beta_{ij}\}$, $\{\tilde{\beta}_{ij}\}$ are scalars to be determined. We assume that for $i < m$, $\mathcal{K}^{i}(B,v_{1})$ is not an invariant subspace of $B$ (similarly, $\mathcal{K}^{i}(\tilde{B},\tilde{v}_{1})$ is not an invariant subspace of $\tilde{B}$ for $i < m$). We can rewrite (\ref{eq:basis1}) as follows:
\begin{equation*}
	{B}{v}_i  = {\beta}_{1i}{v}_1 + {\beta}_{2i}{v}_2 + \ldots + {\beta}_{i-1,i}{v}_{i-1} + {\beta}_{ii}{v}_i  +  {\beta}_{i+1,i}{v}_{i+1}.
\end{equation*}
Combining these equations, for $i \in \{1, 2, 3, \ldots, m-1\}$, into matrix form we get
\begin{align*}
   \begin{array}
	[l]{l}
	B[v_1 \ v_{2} \ldots v_{m-1}] =

   	\left[v_1 \ v_{2} \ldots v_{m-2} \ v_{m-1} \ v_{m}\right]
	\left[
       \begin{array}
        	[c]{ccccc}
        	\beta_{11} & \beta_{12}  & \ldots      & \beta_{1,m-1} \\
        	\beta_{21} & \beta_{22}  & \ldots      & \beta_{2,m-1} \\
        	0         	   & \beta_{32}   & \ldots     & \beta_{3,m-1} \\
        	0         	   & 0                 & \ldots     & \beta_{4,m-1} \\
        	\vdots       & \ddots          & \ddots    & \vdots \\
        	0               &  0                &  \ddots 	 & \beta_{m-1,m-1} \\
        	0               &  0                &  \ldots     & \beta_{m,m-1} \\
      \end{array}
      \right],
   \end{array}
\end{align*}
or
$$BV_{m-1} = V_{m}\underline{H}_{m-1},$$
where $\underline{H}_{m-1}$ is an $m \times (m-1)$ upper Hessenberg matrix. This result also holds for each $i \in \{1, 2, 3, \ldots, m-1\}$, i.e.,
\begin{align}\label{eq:generalRelations1}
    \begin{array}
       [l]{lll}
	BV_i & = & V_{i+1}\underline{H}_i.
    \end{array}
\end{align}
Similarly, using (\ref{eq:basis2}) and following the steps above, we get the following relation for the dual system:
\begin{align}\label{eq:generalRelations2}
    \begin{array}
       [l]{lll}
	\tilde{B}\tilde{V}_i & = & \tilde{V}_{i+1}\tilde{\underline{H}}_i.
    \end{array}
\end{align}
The scalars $\{\beta_{ij}\}$ and $\{\tilde{\beta}_{ij}\}$ are determined by a choice of constraints. One option is to enforce that the columns of $V_i$ (and $\tilde{V}_i$) are orthonormal vectors (as in the Arnoldi algorithm). Another option, as in the bi-Lanczos algorithm, is to enforce\footnote{In this paper, for ease of exposition, we assume breakdowns do not happen. Hence, $(\tilde{v}_i, v_i) \neq  0$.}

\begin{equation*}
	V_i \perp_b \tilde{V}_i, \quad ||v_i|| = 1, \quad \text{and} \quad  (v_i, \tilde{v}_i) = 1,
\end{equation*}
or
\begin{equation}\label{eq:findScalars}
	\tilde{V}_{i}^*V_{i} = I \quad \text{and} \quad  ||v_i|| = 1.
\end{equation}
If $\tilde{B} = B^*$, then (\ref{eq:generalRelations1}), (\ref{eq:generalRelations2}), and (\ref{eq:findScalars}) lead to the bi-Lanczos relations (\ref{eq:biLanczos}), which consist of three-term recurrences. Our goal here is to relax the condition $\tilde{B} = B^*$ and still obtain short-term recurrences.

\begin{theorem}\label{thm:genBi}
	Let $B, \tilde{B} \in \mathbb{C}^{n \times n}$, and let the following conditions

	hold:
	\begin{enumerate}[(a)]
		\item $B  - \tilde{B}^* = \tilde{F}_k \tilde{C}_k^* - C_k F_k^*$, where $C_k, \tilde{C}_k, F_k, \tilde{F}_k  \in \mathbb{C}^{n \times k}$,
		\item $\forall{x}: B x \perp \tilde{C}_k$, $\forall\tilde{x}: \tilde{B} \tilde{x} \perp C_k$,
	 	\item $v_1 \perp \tilde{C}_k$, and $\tilde{v}_1 \perp C_k$.
	\end{enumerate}
	Also, let (\ref{eq:findScalars}) be used as the set of constraints for (\ref{eq:generalRelations1}) and (\ref{eq:generalRelations2}). Then, $\beta_{ij} = 0$ and $\tilde{\beta}_{ij} = 0$  for $j > {i+1}$, which leads to the following three-term recurrences:
	\begin{align*}
	    \begin{array}
	    [l]{lll}
	    \beta_{i+1,i}{v}_{i+1} & = & {B}{v}_i - {\beta}_{ii}{v}_i -  {\beta}_{i-1,i}{v}_{i-1}, \\
	    \tilde{\beta}_{i+1,i}\tilde{v}_{i+1} & = & \tilde{B}\tilde{v}_i - \tilde{\beta}_{ii}\tilde{v}_i -  \tilde{\beta}_{i-1,i}\tilde{v}_{i-1},
	    \end{array}
	\end{align*}
	for $i \in \{1, 2, 3, \ldots, m-1\}$.
\end{theorem}
\begin{proof}
	Using (b) and (c) we can show that
	\begin{equation}\label{eq:CVBiOrth}
		C_k^* \tilde{V}_i = 0 \quad \text{and} \quad \tilde{C}_k^* V_i = 0.
	\end{equation}
	We show $C_k^* \tilde{V}_i =0$ by induction. One can similarly show that $\tilde{C}_k^* V_i = 0$. $C_k^* \tilde{v}_1 = 0$ by (c). Let $C_k^*\tilde{v}_l = 0$ for $l = \{1, 2,\ldots,i\}$, and consider the case $l = i+1$. From (\ref{eq:basis2}) we know that
	\begin{equation*}
		\tilde{\beta}_{i+1,i}\tilde{v}_{i+1} = \tilde{B}\tilde{v}_i - \tilde{\beta}_{ii}\tilde{v}_i -  \tilde{\beta}_{i-1,i}\tilde{v}_{i-1} - \ldots - \tilde{\beta}_{1i}\tilde{v}_1.
	\end{equation*}
	Then, $C_k^* \tilde{v}_{i+1}  = 0$ since $C_k^* \tilde{B}\tilde{v}_i  = 0$ using (b) and $\tilde{\beta}_{li}C_k^*\tilde{v}_l = 0$ for $l \in \{1, 2,\ldots,i\}$ by the induction hypothesis\footnote{Note that our earlier assumption, $\mathcal{K}^{i}(\tilde{B},\tilde{v}_{1})$ is not an invariant subspace of $\tilde{B}$ for $i < m$, shows that $\tilde{\beta}_{i+1,i} \ne 0$.}. This proves (\ref{eq:CVBiOrth}). Multiplying both sides in (\ref{eq:generalRelations1}) by $\tilde{V}_i^*$ and using (\ref{eq:findScalars}) we get
	\begin{equation*}
		\tilde{V}_i^*BV_i = H_i.
	\end{equation*}
	Substituting (a) in the above equation leads to
	\begin{align*}
		\begin{array}
		[l]{llll}
		\tilde{V}_i^*\left(\tilde{B}^* + \tilde{F}_k \tilde{C}_k^* - C_k F_k^*\right)V_i & = & H_i & \Longleftrightarrow \\
		\tilde{V}_i^*\tilde{B}^*V_i + \tilde{V}_i^*\tilde{F}_k \tilde{C}_k^*V_i - \tilde{V}_i^*C_k F_k^*V_i & = & H_i. &
		\end{array}
	\end{align*}
	Using (\ref{eq:CVBiOrth}) we get
	\begin{align*}
		\begin{array}
		[l]{llll}
		\tilde{V}_i^*\tilde{B}^*V_i & = & H_i & \Longleftrightarrow \\
		(\tilde{B}\tilde{V}_i)^*V_i & = & H_i. &
		\end{array}
	\end{align*}
	Finally, using (\ref{eq:generalRelations2}) and (\ref{eq:findScalars}) in the above equation gives
	\begin{equation*}
		\tilde{H}_i^* = H_i.
	\end{equation*}
	This implies both $H_i$ and $\tilde{H}_i$ are tridiagonal matrices, and hence $\beta_{ij} = 0$ and $\tilde{\beta}_{ij} = 0$  for $j > {i+1}$.
\end{proof}

We now revisit augmented bi-Lanczos~\cite{ahuja09ms} and show that it is a special case of generalized bi-Lanczos. The BiCG algorithm is primarily used where the dual system is not auxiliary. That is, one needs to solve both a primary system and a dual system. The recycling BiCG algorithm (also termed RBiCG) was developed to accelerate the convergence of sequences of such systems.

In RBiCG, we use the matrix $U$ to define the primary system recycle space, and compute $C = A^{(\iota+1)}U$, where $U$ is derived from an approximate right invariant subspace of $A^{(\iota)}$ and $\iota$ denotes the index of the linear system in the sequence of linear systems; see (\ref{eq:initial3}). Similarly, we use the matrix $\tilde{U}$ to define the dual system recycle space, and compute $\tilde{C} = A^{(\iota+1)*} \tilde{U}$, where $\tilde{U}$ is derived from an approximate left invariant subspace of $A^{(\iota)}$. $U$ and $\tilde{U}$ are computed such that $C$ and $\tilde{C}$ are bi-orthogonal (see page 35 of~\cite{ahuja11phd}). The number of vectors selected for recycling is denoted by $k$, and hence, $U$, $\tilde{U}$, $C$, and $\tilde{C} \in \mathbb{C}^{n \times k}$.

The bi-Lanczos algorithm was modified to compute the columns of $V_i$ and $\tilde{V}_i$ such that
\begin{equation*}\label{biorth}
    \left[C\ V_{i}\right] \perp_{b} \left[\tilde{C}\ \tilde{V}_{i}\right].
\end{equation*}

Using the scaling (\ref{scalingfactor}), we initialize the Lanczos vectors as
\begin{align*}
    \begin{array}
        [c]{cc}
        v_1 = \frac{\textstyle \left(I - C \mathcal{D}_c^{-1} \tilde{C}^*\right)r_0}{\textstyle \left|\left|\left(I - C \mathcal{D}_c^{-1} \tilde{C}^*\right)r_0\right|\right|}, & \tilde{v}_{1} = \frac{\displaystyle \left(I - \tilde{C} \mathcal{D}_c^{-1} C^*\right)\tilde{r_0}}{\displaystyle \left(v_{1}, \left(I - \tilde{C} \mathcal{D}_c^{-1} C^*\right)\tilde{r_0}\right)}.
    \end{array}
\end{align*}
Here $\mathcal{D}_c = \tilde{C}^* C$ is a diagonal matrix (implied by $C \perp_b \tilde{C}$; we also enforce $\mathcal{D}_c$ to have positive, real coefficients). As for the bi-Lanczos algorithm in (\ref{eq:biLanczosCompute}), the ($i+1$)-st Lanczos vectors here are given by
\begin{align*}
    \begin{array}
        [l]{lll}
	\gamma {v}_{i+1} & = & A v_i - V_i \tau - C \rho \perp \left[\tilde{C}\ \tilde{V}_i\right], \\
	\tilde{\gamma} \tilde{v}_{i+1} & = & A^* \tilde{v}_i - \tilde{V}_i \tilde{\tau} - \tilde{C} \tilde{\rho} \perp \left[{C}\ {V}_i\right], \\
    \end{array}
\end{align*}
where $\gamma$, $\tilde{\gamma}$, $\tau$, $\tilde{\tau}$, $\rho$, and $\tilde{\rho}$ are to be determined. The computation of the ($i+1$)-st Lanczos vector for the primary system now requires the $i$-th and ($i-1$)-st Lanczos vectors and $C$ (see~\cite{ahuja09ms}). This gives a ($3+k$)-term recurrence, where $k$ is the number of columns of $C$. Similarly, we get a ($3+k$)-term recurrence for computing the Lanczos vectors for the dual system. We refer to this pair of ($3+k$)-term recurrences as the augmented bi-Lanczos relations, and they are given by
\begin{align*}
	\begin{array}[l]{lll}
	(I-C \hat{C}^*)AV_i & = & V_{i+1}\underline{T}_i, \\[5pt]
	(I-\tilde{C}\check{C}^*)A^* \tilde{V}_i & = & \tilde{V}_{i+1}\tilde{\underline{T}}_i,
	\end{array}
\end{align*}
where
\begin{align}\label{eq:CcheckChat}
	\begin{array}[l]{lll}
	\hat{C}=\left[\frac{1}{c_{1}^*\tilde{c}_{1}}\tilde{c}_{1}\ \frac{1}{c_{2}^*\tilde{c}_{2}}\tilde{c}_{2}\ \cdots\ \frac{1}{c_{k}^*\tilde{c}_{k}}\tilde {c}_{k}\right] & = & \tilde{C}\mathcal{D}_c^{-1}, \\[5pt]
	\check{C}=\left[\frac{1}{\tilde{c}_{1}^*c_{1}}c_{1}\ \frac{1}{\tilde{c}_{2}^*c_{2}}c_{2}\ \cdots\ \frac{1}{\tilde{c}_{k}^*c_{k}}c_{k}\right] & = & {C}\mathcal{D}_c^{-1}.
	\end{array}
\end{align}

\begin{theorem}\label{thm:genToAug}
 	Let $v_{1} = \eta (I - C \mathcal{D}_c^{-1} \tilde{C}^*)r_0$, $\tilde{v}_1 =  \tilde{\eta}(I - \tilde{C} \mathcal{D}_c^{-1} C^*)\tilde{r}_0$, $B = (I - C \mathcal{D}_c^{-1} \tilde{C}^*)A$, and $\tilde{B} = (I - \tilde{C} \mathcal{D}_c^{-1} C^*)A^*$, where $\eta$, $\tilde{\eta}$ are scalars and $C$, $\tilde{C} \in \mathbb{C}^{n \times k}$ s.t. $\mathcal{D}_c = \tilde{C}^* C$ is a diagonal matrix with positive, real coefficients. 	Also, let (\ref{eq:findScalars}) be used as the set of constraints for (\ref{eq:generalRelations1}) and (\ref{eq:generalRelations2}). Then, $\beta_{ij} = 0$ and $\tilde{\beta}_{ij} = 0$  for $j > {i+1}$, which leads to the following short-term recurrences:
	\begin{align*}
	    \begin{array}
	    [l]{lll}
	    \beta_{i+1,i}{v}_{i+1} & = & {B}{v}_i - {\beta}_{ii}{v}_i -  {\beta}_{i-1,i}{v}_{i-1}, \\
	    \tilde{\beta}_{i+1,i}\tilde{v}_{i+1} & = & \tilde{B}\tilde{v}_i - \tilde{\beta}_{ii}\tilde{v}_i -  \tilde{\beta}_{i-1,i}\tilde{v}_{i-1},
	    \end{array}
	\end{align*}
	for $i \in \{1, 2, 3, \ldots, m-1\}$.
\end{theorem}
\begin{proof}
	We show that conditions (a) -- (c) of Theorem \ref{thm:genBi} are satisfied. This demonstrates that augmented bi-Lanczos is a special case of generalized bi-Lanczos. We have $B, \tilde{B} \in \mathbb{C}^{n \times n}$ such that
	\begin{align*}
		\begin{array}
		[l]{lll}
		B - \tilde{B}^* & = & A -  C \mathcal{D}_c^{-1} \tilde{C}^*A - A + A C\mathcal{D}_c^{-1}\tilde{C}^* \\
					& = & \left(A C\mathcal{D}_c^{-1}\right)\tilde{C}^* -  C \left(A^*\tilde{C}\mathcal{D}_c^{-1}\right)^*.
		\end{array}
	\end{align*}
	Defining $F = A^*\tilde{C}\mathcal{D}_c^{-1}$ and $\tilde{F} = A C\mathcal{D}_c^{-1}$, we get
	$$B - \tilde{B}^* = \tilde{F} \tilde{C}^* - C F^* \quad \text{where} \quad C, \tilde{C}, F, \tilde{F}  \in \mathbb{C}^{n \times k}.$$
	Hence (a) is satisfied. For any $\tilde{x}$ consider the following:
	\begin{align*}
		\begin{array}
		[l]{lll}
		C^*\tilde{B}\tilde{x} & = & C^*(I - \tilde{C} \mathcal{D}_c^{-1} C^*)A^*\tilde{x}\\
						  & = & (C^* - \mathcal{D}_c \mathcal{D}_c^{-1} C^*)A^*\tilde{x} = 0.
		\end{array}
	\end{align*}
	Similarly, for any $x$ consider the following:
	\begin{align*}
		\begin{array}
		[l]{lll}
		\tilde{C}^*Bx & = & \tilde{C}^*(I - C \mathcal{D}_c^{-1} \tilde{C}^*)Ax\\
					& = & (\tilde{C}^* - \mathcal{D}_c \mathcal{D}_c^{-1} \tilde{C}^*)Ax = 0.
		\end{array}
	\end{align*}
Hence (b) is satisfied. Similarly, for $v_1$ and $\tilde{v}_1$ chosen in the theorem, $\tilde{C}^*v_1 = 0$ and $C^*\tilde{v}_1 = 0$. Hence, (c) is satisfied.
 \end{proof}

 For ease of future derivations, we introduce a slight change of notation. Let $x_{-1}$ and $\tilde{x}_{-1}$ be the initial guesses and $r_{-1} = b - A x_{-1}$ and $\tilde{r}_{-1} = \tilde{b} - A^* \tilde{x}_{-1}$ the corresponding initial residuals. We define
\begin{align}\label{minusXR}
    \begin{array}{lcr}
	    x_0 = x_{-1} + U \hat{C}^*r_{-1}, & \qquad & r_0 = (I - C\hat{C}^*)r_{-1}, \\
	    \tilde{x}_0 = \tilde{x}_{-1} + \tilde{U} \check{C}^*\tilde{r}_{-1}, & \qquad & \tilde{r}_0 = (I - \tilde{C}\check{C}^*)\tilde{r}_{-1},
    \end{array}
\end{align}
and follow this convention for $x_0$, $\tilde{x}_0$, $r_0$, and $\tilde{r}_0$ for the rest of the paper. Algorithm 1 gives the RBiCG algorithm from~\cite{ahuja12rbicg}. Here, we have not given details on how the recycle space is computed in RBiCG. For that we refer the reader to~\cite{ahuja12rbicg}.

Like BiCG, breakdowns can happen in RBiCG as well. Besides the breakdown in the underlying augmented bi-Lanczos algorithm ($\tilde{v}_i^*v_i$ = 0 at step $i$; called {\it serious breakdown}), a breakdown can happen when pivotless LDU decomposition of the  tridiagonal matrix (as in the augmented bi-Lanczos relations discussed earlier) does not exist. This is referred to as a breakdown of the {\it second kind}.

The breakdown in the augmented bi-Lanczos algorithm can be avoided by using look-ahead strategies~\cite{freund:lookahead, gutnket1997survey} (as applied for the bi-Lanczos algorithm). The second breakdown can also be avoided in the same way as in BiCG. That is,  by performing the LDU decomposition with $2 \times 2$ block diagonal elements~\cite{bank1993csb}.

\begin{figure}
{\bf Algorithm 1.} {\it RBiCG~\cite{ahuja12rbicg}} \\
    1. Given $U$ (also $C = AU$) and $\tilde{U}$ (also $\tilde{C} = A^* \tilde{U}$) s.t. $C \perp_b \tilde{C}$, compute $\check{C}$ and $\hat{C}$ using (\ref{eq:CcheckChat}).  If $U$ and $\tilde{U}$ are not available, then initialize $U$, $\tilde{U}$, $\check{C}$, and $\hat{C}$ to empty matrices. \\
    2. Choose $x_{-1}$, $\tilde{x}_{-1}$ and compute $x_0$, $\tilde{x}_0$, $r_0$, and $\tilde{r}_{0}$ using (\ref{minusXR}). \\
    3. {\bf if} $(r_0, \tilde{r}_0) = 0$ {\bf then} initialize $\tilde{x}_{-1}$ to a random vector. \\
    4. Set $p_0 = 0$, $\tilde{p}_0 = 0$, $\zeta_c = 0$, $\tilde{\zeta}_c = 0$, and $\beta_0 = 0$.
    Choose {\tt tol} and {\tt max\_itn}. \\
    5. {\bf for} $i = 1 \ldots$ {\tt max\_itn} {\bf do} \\
    \begin{tabular}[b]{ll}
    $\diamond$ \;  \;  $p_i = r_{i-1} + \beta_{i-1} p_{i-1}$; &
                       $\tilde{p}_i = \tilde{r}_{i-1} + \bar{\beta}_{i-1} \tilde{p}_{i-1}$ \\
    $\diamond$ \;  \;  $z_i = Ap_i$; &
                       $\tilde{z}_i = A^*\tilde{p}_i$ \\
    $\diamond$ \;  \;  $\zeta_i = \hat{C}^*z_i$; &
                       $\tilde{\zeta}_i = \check{C}^*\tilde{z}_i$ \\
    $\diamond$ \;  \;  $q_i = z_i - C\zeta_i$; &
                       $\tilde{q}_i = \tilde{z}_i - \tilde{C}\tilde{\zeta}_i$ \\
    $\diamond$ \;  \;  $\alpha_i = (\tilde{r}_{i-1}, r_{i-1}) / (\tilde{p}_i, q_i)$; &
                       $\tilde{\alpha}_i = \bar{\alpha}_i$ \\
    $\diamond$ \;  \;  $\zeta_c = \zeta_c + \alpha_i \zeta_i$; &
                       $\tilde{\zeta}_c = \tilde{\zeta}_c + \tilde{\alpha}_i \tilde{\zeta}_i$ \\
    $\diamond$ \;  \;  $x_i = x_{i-1} + \alpha_i p_i$; &
                       $\tilde{x}_i = \tilde{x}_{i-1} + \tilde{\alpha}_i \tilde{p}_i$ \\
    $\diamond$ \;  \;  $r_i = r_{i-1} - \alpha_i q_i$; &
                       $\tilde{r}_i = \tilde{r}_{i-1} - \tilde{\alpha}_i \tilde{q}_i$ \\
    \multicolumn{2}{l}{$\diamond$ \;  \;  {\bf if} $||r_i|| \leq $ {\tt tol} and
       $||\tilde{r}_i|| \leq $ {\tt tol} {\bf then break}} \\
    $\diamond$ \;  \;  $\beta_i = (\tilde{r}_i, r_i) / (\tilde{r}_{i-1}, r_{i-1})$ & \\
    \end{tabular} \\
    6. {\bf end for} \\
    7. $x_i = x_i - U \zeta_c$; \quad
       $\tilde{x}_i = \tilde{x}_i - \tilde{U}\tilde{\zeta}_c$
\end{figure}

\section{Recycling BiCGSTAB}\label{sec:rbicgstab}
In RBiCG~\cite{ahuja11phd, ahuja12rbicg}, the iteration vectors $p$, $\tilde{p}$, $r$, and $\tilde{r}$ are updated using the following recurrences:
\begin{align*}
	\begin{array}
	[l]{llllll}
    p_i & = & r_{i-1} + \beta_{i-1} p_{i-1}, & \tilde{p}_i & = & \tilde{r}_{i-1} + \tilde{\beta}_{i-1} \tilde{p}_{i-1}, \\
    r_i & = & r_{i-1} - \alpha_i B p_i, & \tilde{r}_i & = & \tilde{r}_{i-1} - \tilde{\alpha}_i \tilde{B}\tilde{p}_i,
    \end{array}
\end{align*}

where $B = (I - C\hat{C}^*)A$ and $\tilde{B} = (I - \tilde{C}\check{C}^*)A^*$. We first give the polynomial representations of these iteration vectors.
\begin{theorem}\label{thm:polyExp}
    Let $r_i$, $p_i$, $\tilde{r}_i$, and $\tilde{p}_i$ be defined as above. Then, for the primary system
    \begin{equation*}\label{eq:polynomialRep}
        r_i = \Theta_i(B)r_0, \quad p_i = \Pi_{i-1}(B)r_0,
    \end{equation*}
    where $\Theta_i(K)$ and $\Pi_{i-1}(K)$ are $i$-th and $(i-1)$-st degree polynomials, for an arbitrary square matrix $K$, that satisfy the following recurrences:
    \begin{align*}
        \Theta_i(K) = \Theta_{i-1}(K) - \alpha_i K \Pi_{i-1}(K), \\
        \Pi_{i-1}(K) = \Theta_{i-1}(K) + \beta_{i-1} \Pi_{i-2}(K).
    \end{align*}
    Similarly, for the dual system
\begin{equation*}
 	\tilde{r}_i = \bar{\Theta}_i(\tilde{B})\tilde{r}_0, \quad \tilde{p}_i = \bar{\Pi}_{i-1}(\tilde{B})\tilde{r}_0,
\end{equation*}
where $\bar{\Theta}_i(K)$ and $\bar{\Pi}_i(K)$ satisfy the following recurrences:
    \begin{align*}
        \bar{\Theta}_i(K) = \bar{\Theta}_{i-1}(K) - \bar{\alpha}_i K\bar{\Pi}_{i-1}(K), \\
        \bar{\Pi}_{i-1}(K) = \bar{\Theta}_{i-1}(K) + \bar{\beta}_{i-1} \bar{\Pi}_{i-2}(K).
    \end{align*}
\end{theorem}
\begin{proof}
	This can be proved by induction, following the derivation in~\cite{sonneveld89cgs} (Section 2; pages 37--40), but use $B$ instead of $A$ and $\tilde{B}$ instead of $A^*$.
\end{proof}

From RBiCG we know $r_i \perp \tilde{r}_j$ for $j<i$. Using Theorem \ref{thm:polyExp} we get that
\begin{equation*}
	(\bar{\Theta}_j(\tilde{B})\tilde{r}_0, \Theta_i(B)r_0) = 0 \text{ for } j < i.
\end{equation*}
This implies $\Theta_{i}(B)r_0$ $\perp$ $\mathcal{K}^i(\tilde{B}, \tilde{r}_{0})$, where $\tilde{r}_0$, $\tilde{B}\tilde{r}_0$, $\ldots$, $\tilde{B}^{i-1}\tilde{r}_0$ span the subspace $\mathcal{K}^i(\tilde{B}, \tilde{r}_{0})$. As observed in~\cite{vorst92bicgstab}, the above orthogonality conditions must be satisfied by other bases of $\mathcal{K}^i(\tilde{B}, \tilde{r}_{0})$, too. So, other polynomials can be used as well~\cite{zhang97gbicgstab}. That is,
\begin{align}\label{newPolyEquation}
	(\bar{\Omega}_j(\tilde{B})\tilde{r}_0, \Theta_i(B)r_0) = 0 \text{ for } j < i.
\end{align}
Similar to the derivation in~\cite{vorst92bicgstab}, we define
\begin{align*}
    	\begin{array}
        [l]{lll}
	\bar{\Omega}_i(\tilde{B}) & = & (I-\bar{\omega}_1 \tilde{B})(I-\bar{\omega}_2 \tilde{B}) \cdots (I-\bar{\omega}_i \tilde{B}),
    	\end{array}
\end{align*}
where $\omega_i$ is selected to minimize the residual $r_i$ w.r.t. $\omega_i$. Then, as first proposed in~\cite{sonneveld89cgs}, instead of (\ref{newPolyEquation}), we use the following form of inner product:
\begin{equation*}
	(\tilde{r}_0, {\Omega}_j(B)\Theta_i(B)r_0) = 0 \text{ for }  j < i,
\end{equation*}
with
\begin{align*}
    	\begin{array}
        [l]{lll}
	{\Omega}_i({B}) & = & (I-{\omega}_1 {B})(I-{\omega}_2 {B}) \cdots (I-{\omega}_i {B}).
    	\end{array}
\end{align*}
This inner product does not require the transpose of $B$, and hence, is appropriate when there is no dual system to solve. Computing the inner product in this fashion, we obtain the recycling BiCGSTAB algorithm (similar to the way BiCGSTAB is obtained from BiCG in~\cite{vorst92bicgstab}). We term our recycling BiCGSTAB as RBiCGSTAB. The algorithm is given in Algorithm 2. Some algorithmic improvements to make the code faster (similar to those discussed in section 6.2 of~\cite{mello2010opt}) are not given here.

Breakdowns in RBiCG (as discussed in the end of Section 3), lead to breakdowns in the RBiCGSTAB algorithm as well. This is similar to how breakdowns in BiCG lead to breakdowns in BiCGSTAB. A breakdown free BiCGSTAB algorithm is proposed in~\cite{cao1993breakdownbicgstab}, which uses the theory of formal orthogonal polynomials. The same theory can be applied to the RBiCGSTAB algorithm.

The BiCGSTAB algorithm also breaks down when the minimization with respect to $\omega_i$ fails. This problem can be avoided by minimizing in two or more dimensions. This led to the development of BiCGSTAB2~\cite{gutnket93bicgstab2} and BiCGSTAB({\em l})~\cite{sleijpen93bicgstabl}. In~\cite{ahuja11phd}, a Recycling BiCGSTAB2 that can be extended to a Recycling BiCGSTAB({\em l}) is proposed.

\begin{figure}
{\bf Algorithm 2.} {\it RBiCGSTAB} \\
    1. Given $U$ (also $C = AU$) and $\tilde{U}$ (also $\tilde{C} = A^* \tilde{U}$) s.t. $C \perp_b \tilde{C}$, compute $D_c = \tilde{C}^*C$ and $\hat{C} = \tilde{C}\mathcal{D}_c^{-1}$. \\
    2. Choose $x_{-1}$; initialize $\tilde{r}_{-1}$ to a random vector; and compute $x_0$, $r_0$, $\tilde{r}_{0}$ using (\ref{minusXR}). \\
    3. {\bf if} $(r_0, \tilde{r}_0) = 0$ {\bf then} rechoose $x_{-1}$ or reinitialize $\tilde{r}_{-1}$ to avoid this condition. \\
    4. Set scalars $\beta_0$ and $\omega_0$ as well as vectors $p_0$, $q_0$, and $x_c$ to zero. \\
    5. Choose {\tt tol} and {\tt max\_itn}. \\
    6. {\bf for} $i = 1 \ldots$ {\tt max\_itn} {\bf do} \\
    \begin{tabular}[b]{l}
    	   $\diamond$\; \; $ p_i = r_{i-1} + \beta_{i-1}p_{i-1} - \beta_{i-1}\omega_{i-1}q_{i-1} $ \\
    	   $\diamond$\; \; $ q_i = Ap_i$ \\
    	   $\diamond$\; \; $ \zeta_i = \hat{C}^*q_i$ \\
    	   $\diamond$\; \; $ q_i = q_i - C \zeta_i$ \\

	   $\diamond$\; \; $ \alpha_{i} = \frac{\left(\tilde{r}_0,\ r_{i-1}\right)}{\left(\tilde{r}_0,\ q_i\right)} $ \\
	   $\diamond$\; \; $ s_i = r_{i-1} - \alpha_i q_i $ \\
	   $\diamond$\; \; $ t_i = A s_i $ \\
	   $\diamond$\; \; $ \gamma_i  = \hat{C}^*t_i $ \\
	   $\diamond$\; \; $ t_i = t_i - C \gamma_i$ \\

	   $\diamond$\; \; $ \omega_i = \frac{(s_i,t_i)}{(t_i,t_i)} $ \\
	   $\diamond$\; \; $ x_i = x_{i-1} + \alpha_i p_i + \omega_i s_i $ \\
	   $\diamond$\; \; $ x_c = x_c + \alpha_i \zeta_i + \omega_i \gamma_i $ \\
	   $\diamond$\; \; $ r_i = r_{i-1} - \alpha_i q_i - \omega_i t_i $ \\
	   $\diamond$ \;  \;  {\bf if} $||r_i|| \leq $ {\tt tol} {\bf then break} \\
	   $\diamond$\; \; $ \beta_{i} = \frac{(\tilde{r}_0,\ r_{i})}{(\tilde{r}_0,\ r_{i-1})} \cdot \frac{\alpha_{i}}{\omega_{i}} $ \\
    \end{tabular} \\
    7. {\bf end for} \\
    8.  $x_i = x_i - U x_c$ \\
    \small{Note: When A does not change for multiple systems, several changes should be made to make this algorithm substantially more efficient. This is done in our implementation for the model reduction problem of Section \ref{sec:experiments}.}
\end{figure}

\section{Numerical Experiments}\label{sec:analysis}

For BiCG, it has been shown that including a left eigenvector into the search space leads to the removal of the corresponding right eigenvector from the right residual (and vice versa)~\cite{desturler99householder}. In our experiments we demonstrate that recycling a left invariant subspace may improve the convergence rate in the RBiCGSTAB algorithm. We consider two examples. For the first example, we perform vertex centered finite volume discretization of the PDE
\begin{align*}
	\mathsf{-(u_x)_x - (u_y)_y + 10u_x -10u_y = 0},
\end{align*}
on a $42 \times 42$ grid unit square resulting in a $1600 \times 1600$ linear system. We use the following boundary conditions: $\mathsf{u}_\text{south} = 1, \mathsf{u}_\text{west} = 1, \mathsf{u}_\text{north} = 0, \mathsf{u}_\text{east} = 0$. We do not use a preconditioner in this example, the initial guess is a vector of all ones, and the relative convergence tolerance is $10^{-10}$.

For the second example, we perform finite difference discretization of the partial differential equation~\cite{vorst92bicgstab}
\begin{align*}
    \mathsf{-({A}v_x)_x - ({A}v_y)_y + {B}(x,y)v_x = {F}},

\end{align*}
with $\mathsf{A}$ as shown in Figure \ref{figure:secondpde}, $\mathsf{B(x,y) = 2e^{2(x^2 + y^2)}}$, and $\mathsf{F} = 0$ everywhere except in a small square in the center where $\mathsf{F} = 100$ (see Figure \ref{figure:secondpde}). The $(0,1) \times (0,1)$ domain is discretized on a $129 \times 129$ grid resulting in a  $16129 \times 16129$ linear system. We use the following boundary conditions:
\begin{align*}
    \begin{array}
    [l]{l}
    \mathsf{v(0, y) = v(1, y) = v(x, 0) = 1}, \\
    \mathsf{v(x, 1) = 0}.
    \end{array}
\end{align*}

We use an ILUTP~\cite{saad03book} preconditioner with a drop tolerance of $0.1$ (split-preconditioned). The initial guess is $0.5$ times a vector of all ones, and the relative convergence tolerance is $10^{-8}$.

\begin{figure}
    \centering
    \includegraphics[angle=90,scale=0.5]{secondpde.eps}
    \caption{Coefficients for the PDE.}\label{figure:secondpde}
\end{figure}

\begin{figure}
    \subfigure[Example 1: Left eigenvectors not needed, but recycling effective.]{
        	\scalebox{0.6}{\includegraphics{both.eps}}
    }
    \subfigure[Example 2: Left eigenvectors needed for recycling to be effective.]{
        	\scalebox{0.6}{\includegraphics{right.eps}}
    }
    \caption{Convergence curves for two examples using RBiCGSTAB. The 2nd example demonstrates that recycling a left invariant subspace may improve the convergence rate in the RBiCGSTAB algorithm.}\label{figure:needLeftEigVecs}
\end{figure}

For each example we do three experiments. First, we solve the system without recycling. Second, we use the right invariant subspace (corresponding to the smallest magnitude eigenvalues) as the recycle space. This is implemented by setting $\tilde{U} = U$. Finally, we use both the left and right invariant subspaces (again, corresponding to the smallest magnitude eigenvalues) as the recycle space.

For the first example's second set of experiments, we use five exact right eigenvectors computed using the MATLAB function {\tt eigs}. For the first example's third set of experiments, we use five exact left eigenvectors and five exact right eigenvectors (for a total of ten), again computed using the MATLAB function {\tt eigs}.

For the second example's second set of experiments, we use twenty approximate right eigenvectors obtained by solving the problem twice with RBiCG. For the second example's third set of experiments, we use twenty approximate left eigenvectors and twenty approximate right eigenvectors (for a total of forty), again obtained by solving the problem twice with RBiCG.

The results are shown in Figures \ref{figure:needLeftEigVecs} (a) and (b). For the first example, using the right invariant subspace or using both the left and right invariant subspaces works equally well. However, for the second example, we see that using only the right invariant subspace leads to convergence that is worse than BiCGSTAB without recycling, and much worse than RBiCGSTAB using both the left and right invariant subspaces. This shows that recycling a left invariant subspace may improve the convergence rate in the RBiCGSTAB algorithm.

Next, we analyze why the first example does not need left invariant subspace and the second example does, by considering the cosines of the principal angles between the left and right invariant subspaces associated with the ten smallest magnitude eigenvalues. Table \ref{table:analysis} lists these angles.

\begin{table}[tbh]
   \footnotesize
    \begin{center}
        \setlength{\extrarowheight}{1.5pt}
        \begin{tabular}{|m{1in}|m{1in}|}
            \hline
            {Example 1} & {Example 2} \\
            \hline
	   0.9998	&	0.1039    \\
	   0.9837	&	0.0302    \\
	   0.9207	&	0.0195    \\
	   0.7617	&	0.0106	  \\
	   0.4352	&	0.0089	  \\
	   0.3987	&	0.0049	  \\
	   0.2273	&	0.0043	  \\
	   0.0963	&	0.0027	  \\
	   0.0235	&	0.0018	  \\
	   0.0064	&	0.0012	  \\
          \hline
        \end{tabular}
    \end{center}
         \quad \\
    \caption{\label{table:analysis} For each example, we give the cosines of the principal angles between the  exact left and right invariant subspaces of dimension $10$, associated with the smallest magnitude eigenvalues. For the first example, we compute invariant subspaces of the matrix obtained after discretization. While for the second example, we compute invariant subspaces of the preconditioned matrix since we split-precondition the linear system obtained after discretization.}
\end{table}

From the table, we see that the principal angles between the left and right invariant subspaces for the second example are substantially larger than those for the first example (since the cosines of the principal angles is lesser for second example as compared with the first example). Since for a normal matrix left and right invariant subspaces are identical, we conclude that the second example is ``more'' non-normal than the first.

\section{Applications}\label{sec:appl}
We first discuss several techniques for parametric model order reduction and the one we are using (in Section \ref{sec:pmor}) followed by the description of how recycling BiCGSTAB is applied for sequences of linear systems arising in parametric model order reduction (in Section \ref{sec:experiments}).

\subsection{Parametric Model Order Reduction}\label{sec:pmor}
Numerical simulation is an essential tool for solving science and engineering problems. However, simulating large-scale models leads to overwhelming demands on computational resources. This is the main motivation for model reduction. The goal is to produce a surrogate model of much smaller dimension that provides a high-fidelity approximation of the input-output behavior of the original model. Often the models have design parameters associated with them, e.g., boundary conditions, geometry, material properties etc. Changes in these design parameters require generation of new reduced models, which makes the model reduction process very cumbersome. One practical application where such a challenge arises is micro-electro-mechanical systems (MEMS) design~\cite{feng2012Mems, feng2013book}. The goal of parametric model order reduction (PMOR)~\cite{baur2011PMOR,

feng:proceedings} is to generate a reduced model such that parametric dependence, as in the original model, is preserved (or retained).

We focus on physical processes that are modeled as parameterized partial differential equations (PDEs). For PMOR, the PDE is first semi-discretized using classical techniques (e.g., finite differences, finite elements, etc.), and then model reduction is applied to the resulting parameterized state-space model~\cite{baur2011PMOR}:
\begin{align}\label{eq:dynamicalSys1}
    \begin{array}
        [l]{ll}
        \mathsf{G}:
        \begin{cases}
           \mathsf{E(p)\ {\dot x}(t) = A(p)\ x(t) + B(p)\ u(t)} \\
            \mathsf{y(t) = C(p)\ x(t)}
        \end{cases}
    \end{array}
\end{align}
or
\begin{align}\label{eq:dynamicalSys2}
    \mathsf{G(s, p) = C(p)\left(sE\left(p\right) - A\left(p\right)\right)^{-1}B(p)}.
\end{align}
Sometimes, the physical process is directly available in the form of a parametrized state-space model. Equations (\ref{eq:dynamicalSys1}) -- (\ref{eq:dynamicalSys2}) represent a multiple input multiple output (MIMO) linear dynamical system, where $\mathsf{p}$ is the parameter vector; $\mathsf{u(t) = [u_1(t), \ldots}$, $\mathsf{u_m(t)]^T}$: $\mathsf{\mathbb{R} \rightarrow \mathbb{R}^m}$ is the input; $\mathsf{y(t)}$: $\mathsf{\mathbb{R} \rightarrow \mathbb{R}^l}$ is the output; $\mathsf{x(t)}$: $\mathsf{\mathbb{R} \rightarrow \mathbb{R}^n}$ is the state vector; $\mathsf{E(p)}$, $\mathsf{A(p) \in \mathbb{R}^{n \times n}}$, $\mathsf{B(p) \in \mathbb{R}^{n \times m}}$, and $\mathsf{C(p) \in \mathbb{R}^{l \times n}}$ are the system matrices; and $\mathsf{s}$ is the frequency domain variable corresponding to $\mathsf{t}$ in the time domain.

Above, (\ref{eq:dynamicalSys1}) denotes the dynamical system, and (\ref{eq:dynamicalSys2}) gives the transfer function of the system obtained after Laplace transformation. By a common abuse of notation, we denote both with $\mathsf{G}$. The dimension of the underlying state-space, $\mathsf{n}$, is called the dimension or order of $\mathsf{G}$. Unless explicitly stated, for the rest of this paper all dynamical systems are assumed to be of the above type.

There are various ways of performing PMOR~\cite{benner2013survey, baur2011PMOR, patera06reducedbasis, benner2014moment}. This includes moment matching, local $\Htwo$-optimality, and reduced basis approaches. For this work, we focus on moment matching based PMOR because of its flexibility (few limits on the system properties) and low computational cost in many industrial applications.

Moment matching based PMOR algorithms~\cite{benner2014moment, feng2012Mems} require solution of sequences of linear systems of the type (\ref{eq:initial3}), which is a key bottleneck when using these algorithms for reducing larger models. Specifically, the systems have the form as follows:
\begin{align*}
    {A}^{(1)}\ x^{(1, 1)} & = {b}^{(1, 1)} \\
    {A}^{(1)}\ x^{(1, 2)} & = {b}^{(1, 2)} \\
    {A}^{(1)}\ x^{(1, 3)} & = {b}^{(1, 3)} \\
    & \vdotswithin{=}   \\
    {A}^{(2)}\ x^{(2, 1)} & = {b}^{(2, 1)} \\
    {A}^{(2)}\ x^{(2, 2)} & = {b}^{(2, 2)} \\
    {A}^{(2)}\ x^{(2, 3)} & = {b}^{(2, 3)} \\
    & \vdotswithin{=}
\end{align*}
In the example used here, each matrix has the form
\begin{align*}
{A}^{(i)} = A_0 + \sigma_iA_1 + \delta_i A_2 \quad \text{with} \quad i \in \mathbb{N}.
\end{align*}
Please note that $\sigma_iA_1$ and $\delta_i A_2$ are small perturbations to $A_0$. As discussed in the introduction, our goal here is to use Krylov subspace recycling (recycling BiCGSTAB specifically) to efficiently solve such systems.

\subsection{Application to PMOR}\label{sec:experiments}
Our test dynamical system comes from a silicon nitride membrane model~\cite{bechtold2010SiN}. Such a membrane can be part of many devices, e.g., a gas sensor chip, a microthruster, an optical filter etc.

We use the moment matching based PMOR algorithm, described in~\cite{benner2014moment, feng2012Mems}, to compute a reduced model. This leads to a sequence of linear systems of the form (\ref{eq:initial3}) and size $60,020$.

Whenever the matrix changes in the sequence, we call RBiCG to perform the linear solve. This helps to approximate both left and right invariant subspaces, which are not easily available from the RBiCGSTAB iterations (sometimes a left invariant subspace is available from a right invariant subspace~\cite{abdel2007lefteig, morgan2010nlandr}). The primary system right-hand side comes from the PDE. We take a vector of all ones as the dual system right-hand side. We call RBiCGSTAB for all remaining systems with the same matrix. This corresponds to linear systems where only the right-hand sides change. This is an effective strategy because it has been shown that the recycle space can be useful for multiple consecutive systems~\cite{parks06gcrodr, kilmer2006dot, morgan2010nlandr,  abdel2014rhs}. Moreover, using RBiCG for all systems will be expensive since an unnecessary dual system would be solved at each step in the sequence. It needs to be emphasized here that these is a need to implement this efficiently (as discussed in Algorithm 2).

While solving a linear system with RBiCG, Lanczos vectors are generated at each iterative step. These Lanczos vectors are used to build the recycle space. We have the flexibility in deciding when to build the recycle space. One option is to wait for the RBiCG to converge, save all the Lanczos vectors, and then build the recycle space. The problem with this approach is that this requires large amounts of memory and may be computationally expensive as well. Instead, we divide the RBiCG iteration in cycles of a certain number of iterations (to be chosen, e.g., $50$ iterations). At the end of each cycle, we use the stored Lanczos vectors from that cycle to build or improve the recycle space, and then discard these Lanczos vectors, except for the last few that are needed to continue the Lanczos iteration. At the end of the first cycle, we may build a new recycle space or update a recycle space constructed for a previous linear system. In the following discussion, we denote the length of a cycle by $s$.

For this experiment, we take $s = 25$ and $k = 20$ (the number of vectors selected for recycling as defined earlier in Section 3). These values are chosen based on experience with other recycling algorithms~\cite{parks06gcrodr}. The linear systems are split-preconditioned with an incomplete LU preconditioner with threshold and pivoting (ILUTP)~\cite{saad03book}. The drop tolerance is taken as $10^{-4}$. For RBiCG, we take a vector of all zeros as the initial guess for both the primary system and the dual system. For RBiCGSTAB we take a vector of all zeros as the initial guess as well (we are only solving the primary system in RBiCGSTAB). Using the solution from the previous linear system an initial guess leads to poor convergence for BiCGSTAB with and without recycling. In general, a better initial guess may be based on knowledge of the system. The relative convergence tolerance for the iterative solves is taken as $10^{-8}$.

The number of matrix-vector products required to solve systems $1$ through $63$ are given in Figure \ref{figure:matvecs}, and the corresponding timing data is given in Figure \ref{figure:time}. In both the figures, the peaks in the recycling BiCGSTAB plot correspond to when the matrix changes and RBiCG is called (three times; at the 1st, 22nd, and 43rd linear system). For all other steps, when only the right-hand side changes, RBiCGSTAB is called. When recomputing the recycle space for the 22nd system we use the recycle space generated while solving the 1st system, and make it better. This is evident in Figure 3 and 4 where both the number of matrix-vector products and time for solving the 22nd linear system are less as compared with the 1st system. The same process happens when recomputing the recycle space for the 43rd system. Here, the recycle space from the 22nd linear system is improved. Again, the Figures 3 and 4 show the decrease in the matrix-vector product count and time for the 43rd system as compared with the 22nd system.

First, we compare our results with BiCGSTAB. The RBiCG and RBiCGSTAB combination requires about 40\% fewer matrix-vector products in total. Also, computing the reduced model with recycling takes about 35\% less time than without recycling. This demonstrates the effectiveness of recycling Krylov subspaces for PMOR.

Second, we compare our results with GCRO-DR. Although the RBiCG and RBiCGSTAB combination does not beat GCRO-DR in the number of matrix vector products, it is 10\% more efficient than GCRO-DR in time and reduces storage requirements. This is similarly seen in BiCGSTAB and GMRES comparison for some examples. That is, although BiCGSTAB is expensive than GMRES in the number of matrix vector products, it is cheaper than GMRES in time. The reason being that BiCGSTAB is cheaper than GMRES with respect to the orthogonalizations.

\begin{figure}
    \centering
    \includegraphics[scale=0.7]{matvec.eps}
    \caption{Comparison of matrix-vector product count when using BiCGSTAB, RBiCGSTAB, and GCDR-DR as the linear solvers for PMOR.}\label{figure:matvecs}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[scale=0.7]{time.eps}
    \caption{Comparison of time (in seconds) when using BiCGSTAB, RBiCGSTAB, and GCRO-DR as the linear solvers for PMOR.}\label{figure:time}
\end{figure}

\section{Conclusions}\label{sec:conclusion}

For solving linear systems with non-symmetric matrices, B\-i\-C\-G\-S\-T\-A\-B is one of the best available algorithms. As compared with GMRES, which is the most commonly used algorithm for such linear systems, it has the advantage of a short term recurrence, and hence, does not suffer from storage issues.

For solving sequences of linear systems with non-symmetric matrices, it is advantageous to use Krylov subspace recycling for the BiCGSTAB algorithm, and hence we propose the RBiCGSTAB algorithm. We have demonstrated the usefulness of RBiCGSTAB for a parametric model order reduction example.

Here, we have used RBiCG to solve a linear system and generate the recycle space. We have then used RBiCGSTAB, that uses this generated space, to solve the subsequent systems until the matrix changes. RBiCG in this approach can be replaced with other solvers like GCRO-DR~\cite{parks06gcrodr} as well.

In the future, we plan to test RBiCGSTAB for other application areas (e.g., acoustics problems). We also plan to extend the recycling framework of RBiCGSTAB to BiCGSTAB({\em l})~\cite{sleijpen93bicgstabl} and IDR~\cite{wesseling80idr, sonneveld2008idr}.

In section \ref{sec:rbicgstab}, we saw that BiCGSTAB (and RBiCGSTAB) performs one-dimensional minimization of the residual. This minimization can be done in higher dimensions as well (say {\em l}), leading to BiCGSTAB({\em l})\footnote{One-dimensional minimization was first extended to two dimensions in~\cite{gutnket93bicgstab2}}. Like BiCGSTAB, the induced dimension reduction (IDR) method involves a short term recurrence and has been shown to perform better than BiCGSTAB in many cases~\cite{sonneveld2008idr}.

\vspace{6pt}
\noindent
{\bf Acknowledgments.} We thank the anonymous reviewers for their careful and helpful suggestions, which greatly helped us to improve this paper.

\bibliographystyle{abbrv}
\bibliography{rbicgstab}

\end{document}

config: arxiv


text:
\begin{definition}[Definition:Gentzen Proof System/Instance 1/Beta-Rule]
$(\beta)$: For any [[Definition:Beta-Formula|$\beta$-formula]] $\mathbf B$ and associated $\mathbf B_1, \mathbf B_2$ from the [[Definition:Table of Beta-Formulas|table of $\beta$-formulas]]:
:: Given $U_1 \cup \left\{{\mathbf B_1, \mathbf B_2}\right\}$, one may infer $U_1 \cup \left\{{\mathbf B}\right\}$.
\end{definition}

config: wiki


text:
\begin{document}
\IEEEoverridecommandlockouts
\maketitle
\thispagestyle{empty}

\begin{abstract}
Consider a scenario where a source continuously monitors an object and sends time-stamped status updates to a destination through a rate-limited link. In order to measure the ``freshness" of the status information available at the destination, we adopt the metric called Age of Information (AoI). We assume all updates are of the same size, and arrive randomly at the source according to a Bernoulli process. Due to the link capacity constraint, it takes $d$ ($d\geq 2$) time slots for the source to complete the transmission of an update. Therefore, when a new update arrives at the source during the transmission of another update, the source needs to decide whether to skip the new arrival or to switch to it, in order to minimize the expected average AoI at the destination. We prove that within a broadly defined class of online policies, the optimal policy should be a renewal policy, and has a sequential switching property. We then show that the optimal decision of the source in any time slot has a multiple-threshold structure, and only depends on the age of the update being transmitted and the AoI in the system. The thresholds are then numerically identified by formulating the problem as a Markov Decision Process (MDP).
\end{abstract}
\begin{IEEEkeywords}
 Age of information, online scheduling.
  \end{IEEEkeywords}

\section{Introduction}
Enabled by the proliferation of ubiquitous sensing devices and the pervasive wireless data connectivity, real-time monitoring has become a reality in large-scale cyber-physical systems, such as power grids, manufacturing facilities, and smart transportation systems. However, the unprecedented high-dimensionality and generation rate of the sensing data also impose critical challenges on its timely delivery. In order to measure and ensure the freshness of information available to the central controller, a metric called Age of Information (AoI) has been introduced and analyzed in various networks~\cite{infocom/KaulYG12}. Specifically, at time $t$, the AoI in the system is defined as $t-u(t)$, where $u(t)$ is the time stamp of the latest received update at the destination. Since AoI depends on data generation as well as queueing and transmission, it exhibits fundamental differences between traditional network performance metrics, such as throughput and delay.

Modeling the status updating process as a queueing process, time average AoI has been analyzed in systems with a single server~\cite{infocom/KaulYG12,ciss/KaulYG12,isit/YatesK12,YatesK16,Pappas:2015:ICC,isit/NajmN16,isit/KamKNWE16,isit/ChenH16}, and multiple servers \cite{isit/KamKE13, isit/KamKE14,tit/KamKNE16}. Peak Age of Information (PAoI) has been introduced and studied in \cite{isit/CostaCE14,tit/CostaCE16,isit/HuangM15}. The optimality properties of a preemptive Last Generated First Served service discipline are identified in \cite{isit/BedewySS16}.

AoI minimization has also been investigated, either by controlling the generation process of the updates~\cite{infocom/SunUYKS16,SunPU17,isit/Yates15,ita/BacinogluCU15,Yang:2017:AoI,BacinogluU17,Ahmed:2018:ICC,Ahmed:2018:ITA,Uysal:2018:EH,Feng:2018:INFOCOM,Feng:2018:ISIT}, or by scheduling the transmission of updates that have already been generated~\cite{Modiano:2018:BC,Hsu:2017:ISIT,Hsu:2018:ISIT,Kadota:2018:INFOCOM,He:TIT:2017}. Optimal status updating policy with knowledge of the server state has been studied in~\cite{infocom/SunUYKS16}. AoI-optimal sampling of a Wiener process is investigated in \cite{SunPU17}. Under an energy harvesting setting, optimal status updating have been studied in~\cite{isit/Yates15,ita/BacinogluCU15,Yang:2017:AoI,BacinogluU17,Ahmed:2018:ICC,Ahmed:2018:ITA,Uysal:2018:EH,Feng:2018:INFOCOM,Feng:2018:ISIT}. Transmission scheduling in a broadcast channel has been studied in \cite{Modiano:2018:BC,Hsu:2017:ISIT,Hsu:2018:ISIT}. Reference \cite{Modiano:2018:BC} shows that a greedy policy which always tries to update the most outdated client is optimal in a symmetric setting. Reference \cite{Hsu:2017:ISIT} formulates the problem as a Markov Decision Process (MDP), and show that the optimal policy is a switch-type. It also proposes a sequence of finite-state approximations for the infinite-state MDP and proves its convergence.
A restless bandits based formulation and a Whittle's index based scheduling have been studied in \cite{Hsu:2018:ISIT}. Different transmission scheduling policies for AoI minimization in a multiple access channel under throughput constraints on individual nodes have been analyzed in \cite{Kadota:2018:INFOCOM}.
Age-optimal link scheduling in a multiple-source system with conflicting links is studied in~\cite{He:TIT:2017}, and the problem is shown to be NP-complete in general.
Head-of-line age-based scheduling algorithms have been shown to be throughput optimal in wireless networks in~\cite{Srikant:2015:ageSchedule}.

In this paper, we investigate the optimal online transmission scheduling for a single link under the assumption that the link capacity is limited and each update takes multiple time slots to transmit. During the transmission of an update, new updates may arrive. Therefore, the source has to decide whether to switch to the new arrival, or to continue its current transmission and drop the new update. What makes the problem challenging is that the impact of a decision on the AoI evolution won't become clear immediately. This is because the instantaneous AoI at the destination will be reset only when a transmission is completed. Even if the source decides to transmit an update at an earlier time, it may drop the update later before the transmission is complete, leading to uncertain AoI evolution in the system. To overcome this challenge, we first prove that within a broadly defined class of online policies, the optimal policy should be a renewal policy, and the decision-making over each renewal interval only depends on the arrival time of the updates in that interval. Then, we show that the optimal renewal policy has a multiple-threshold structure, which enables us to formulate the problem as an MDP, and identify the thresholds numerically through structured value iteration.

\section{System Model and Problem Formulation} \label{sec:model}
We consider a single-link status monitoring system where the source keeps sending time-stamped status updates to a destination through a rate-limited link. We assume the time axis is discretized into time slots, which are labeled as $t=1, 2, 3,\cdots$. At the beginning of time slot $t$, an update packet is generated and arrives at the source according to an independent and identically distributed (i.i.d.) Bernoulli process $A(t)$ with parameter $p$. We assume each update is of the same size, and it takes exactly $d$ time slots, $d\geq 2$, to transmit one update to the destination. Similar to \cite{Modiano:2018:BC,Hsu:2017:ISIT,Hsu:2018:ISIT}, we assume that at most one update can be transmitted during each time slot, and there is {\it no buffer} at the source to store the updates that are not being transmitted. Therefore, once an update arrives at the source, it needs to decide whether to transmit it and drop the one currently under transmission if there is any, or to drop the new arrival.

A status update policy is denoted as $\pi$, which consists of a sequences of transmission decisions $\{D(t)\}$. We let $D(t)\in\{0,1\}$. Specifically, when $A(t)=1$, $D(t)$ can take both values 1 and 0: If $D(t)=1$, the source will start transmitting the new arrival in time slot $t$ and drop the unfinished update if necessary. We term this as {\it switch}; Otherwise, if $D(t)=0$, the source will drop the new arrival, and continue the unfinished transmission. We term this as {\it skip}. When $A(t)=0$, we can show that dropping the update being transmitted is sub-optimal. Thus, we restrict to the policies under which $D(t)$ can only take value 0, i.e., to continue transmitting the unfinished update if there is one, or to idle.

Let $S_n$ be the the time slot when an update is completely transmitted to the destination. Then,
the inter-update delays can be denoted as $X_n:=S_n-S_{n-1}$, for $n=1,2,\ldots$. Without loss of generality, we assume $S_0=0$. Note that under the bufferless assumption, the AoI after a completed transmission is always equal to $d$. An example sample path of the AoI evolution under a given status update policy is shown in Fig.~\ref{fig:AoI}.
As illustrated, some updates are skipped when they arrive, while others are transmitted partially or completely.

We use $N(T)$ to denote the total number of successfully delivered status updates over $(0,T]$.
Define $R(T)$ as the total age of information experienced by the system over $[0,T]$. Denote $R_n:=(2d+X_n)X_n/2$, i.e., the total AoI experienced by the receiver over the $n$th epoch $X_n$. Then,
\begin{align*}
R(T)&=\sum_{n=1}^{N(T)}R_n+ \frac{1}{2}(d+T-S_{N(T)})(T-S_{N(T)}).
\end{align*}
We focus on a set of {\it online} policies $\Pi$, in which the information available for determining $D(t)$ includes the decision history $\{D(i)\}_{i=1}^{t-1}$, the update arrival profile $\{A(i)\}_{i=1}^{t}$, as well as the update statistics (i.e., $p$ in this scenario).
The optimization problem can be formulated as
\begin{eqnarray}\label{eqn:opt}
\underset{\pi\in\Pi}{\min} & &\limsup_{T\rightarrow \infty} \Eb\left[\frac{R(T)}{T}\right]\end{eqnarray}
where the expectation in the objective function is taken over all possible update arrival sample paths.

\begin{figure}[t]
	\centering
\includegraphics[width=6.5cm]{AoI.pdf}
\vspace{-0.1in}
	\caption{\small{AoI evolution with $d=3$. Circles represent transmitted updates, and crosses represent skipped ones. Red dashed curve indicates the transmitted portion of the corresponding update.} }
	\label{fig:AoI}
	\vspace{-0.1in}
\end{figure}

\section{Structure of the Optimal Policy}
Consider the $n$th epoch, i.e., the duration between time slots $S_{n-1}+1$ and $S_n$ under any online policy in $\Pi$. Let $a_{n,k}$ be the time slot when the $k$th update after $S_{n-1}$ arrives, and let $x_{n,k}:=a_{n,k}-S_{n-1}$. Denote the update arrival profile in epoch $n$ as $\xv_n:=(x_{n,1},x_{n,2},\ldots)$. Then, we introduce the following definition.

\begin{Definition}[Uniformly Bounded Policy]
Under an online policy $\pi\in\Pi$, if there exists a function $g(\xv)$, such that for any $\xv_n=\xv$, the length of the corresponding epoch $X_n$ is upper bounded by $g(\xv)$, and $\Eb[g^2(\xv)]<\infty$, then this policy is a uniformly bounded policy.
\end{Definition}

Denote the subset of uniformly bounded policies as $\Pi'$. Then, using techniques similar to the proof of Theorem 1 in \cite{Ahmed:2018:ICC}, we can show the following theorem.
\begin{Theorem}\label{thm:renewal}
Any uniformly bounded policy $\pi\in\Pi'$ is sub-optimal to a renewal policy. That is, $\{S_n\}_{n=1}^{\infty}$ form a renewal process. Besides, the decision $D(t)$ over the $n$th renewal epoch only depends on $\xv_n$ causally.
\end{Theorem}

Due to space limitation, the proof of Theorem~\ref{thm:renewal}, as well as the proofs of Lemma~\ref{lemma:idle}, Lemma~\ref{lemma:threshold} and Theorem~\ref{thm:threshold} are omitted.

\if{0}
\begin{Proof}
Assume $\{X_i\}$ satisfies uniformly bounded condition.
For a fixed history $\Hc^{i-1}$, let us group all the status updating sample paths that have the same $\tauv_i$ and perform a statistical averaging over all of them to get the following average age in the $i$th epoch
\begin{align}
\hat{R}_i (\gammav, \Hc^{i-1}) := \Eb [R_i |\tauv_i=\gammav, \Hc^{i-1}]
\end{align}
Then, we have
\begin{align}
\Eb[\hat{R}_i \lv_{i\leq N(T)+1}]=\Eb_{\Hc^{i-1}}[\Eb_{\gammav}[\hat{R}_i (\gammav, \Hc^{i-1})] \lv_{i\leq N(T)+1}|\Hc^{i-1}]
\end{align}
where the equality follows from the fact that $\lv_{i\leq N(T)+1}$ is independent of $\tauv_i$ given $\Hc^{i-1}$.

Similarly, define the average $i$th epoch length as
\begin{align}
\hat{X}_i(\gammav, \Hc^{i-1}):=\Eb [X_i |\tauv_i=\gammav, \Hc^{i-1}]
\end{align}

Note that for any uniformly bounded policy, we have $\lim_{T\rightarrow\infty}\frac{X^2_{N(T)+1}}{T}=0$ \cite{Feng1806:Minimizing}.
Then,
\begin{align}
&\lim_{T\rightarrow\infty}\frac{\Eb[R(T)]}{T}=\lim_{T\rightarrow\infty}\frac{\Eb\left[\sum_{i=1}^{N(T)+1}R_i\right]}{T}  \label{eqn:thm-1}\\
&\geq \frac{\Eb\left[\sum_{i=1}^{N(T)+1}R_i\right]}{2\Eb[S_{N(T)+1}]}=\frac{ \sum_{i=1}^{\infty} \Eb\left[R_i\lv_{i\leq N(T)+1}\right]}{ \sum_{i=1}^{\infty}\Eb[X_i\lv_{i\leq N(T)+1}]}\\
&=\frac{ \sum_{i=1}^{\infty} \Eb_{\Hc^{i-1}}\left[\Eb_{\gammav}\left[\hat{R}_i(\gammav,\Hc^{i-1})\right]\lv_{i\leq N(T)+1}|\Hc^{i-1} \right]}{ \sum_{i=1}^{\infty}\Eb[X_i\lv_{i\leq N(T)+1}]}\\
&=\frac{ \sum_{i=1}^{\infty} \Eb_{\Hc^{i-1}}\left[\Eb_{\gammav} [\hat{X}_i(\gammav, \Hc^{i-1})]\frac{\Eb_{\gammav}\left[\hat{R}_i(\gammav,\Hc^{i-1})\right]\lv_{i\leq N(T)+1}}{\Eb_{\gammav} [\hat{X}_i(\gammav, \Hc^{i-1})]} |\Hc^{i-1}\right]}{\sum_{i=1}^{\infty}\Eb[X_i\lv_{i\leq N(T)+1}]}\\
&\geq\frac{ \sum_{i=1}^{\infty} \Eb_{\Hc^{i-1}}\left[\Eb_{\gammav} [\hat{X}_i(\gammav, \Hc^{i-1})] R^*\lv_{i\leq N(T)+1} |\Hc^{i-1}\right]}{\sum_{i=1}^{\infty}\Eb[X_i\lv_{i\leq N(T)+1}]}\\
&= R^*
\end{align}
where $R^*:=\min_{i,\Hc^{i-1}}\frac{\Eb_{\gammav}\left[\hat{R}_i(\gammav,\Hc^{i-1})\right]}{\Eb_{\gammav} [\hat{X}_i(\gammav, \Hc^{i-1})]}$, which corresponds to the expected average AoI achieved by a renewal policy under which $X_i$ only depends on $\gammav$.
\end{Proof}
\fi

Based on Theorem~\ref{thm:renewal}, in the following, we will focus on renewal policies that depend on $\xv_n$ only.

\begin{Lemma}\label{lemma:idle}
If the source is idle when an update arrives, it should start transmitting the update immediately.
\end{Lemma}

\if{0}
\begin{Proof}
 If $t+d=S_n$, the file arrives at $t$ should transmit immediately to reach the receiver at $S_n$. If $t+d<S_n$, transmit the file immediately will not change $S_n$ due to the fact that a new arrival file can replace the transmitting one. Besides, the age of the source keeps growing no matter the source is idle or not. Therefore, if a files arrives at the transmitter when it is idle, the transmitter should start transmitting it immediately.
\end{Proof}

We note the the optimal policy is to decide whether the transmitter should switch to a new update when it arrives.
\fi
\begin{Definition}[Sequential Switching Policy]
A sequential switching (SS) policy is a renewal policy under which the source switches to an update arriving at time slot $t$ \textbf{only if} it switches to all update arrivals prior to $t$ in the same epoch.
\end{Definition}

\textbf{Remark:} The definition of SS policy implies that once a source skips a new update arrival at $t$, it will skip all of the upcoming update arrivals until it finishes the one being transmitted at $t$. We point out that an SS policy is in general different from threshold type of policies, as it does not impose any threshold structure on when the source should skip or switch to a new update arrival.

\begin{Lemma}\label{lemma:ssp}
The optimal renewal policy in $\Pi'$ is an SS policy.
\end{Lemma}
\begin{Proof}
We prove this lemma through contradiction. Now assume the optimal policy $\pi_0$ is not an SS policy. Without loss of generality, we consider the first renewal epoch starting at time 0 (the beginning of time slot $1$). We assume under $\pi_0$ there exists a sample path under which the source transmits the new update arrival at time slot $i$ and does not switch to the next arrival at time slot $j$ in the same epoch, i.e., $i<j<i+d$. Depending on the upcoming random arrivals, the sample path may evolve into different sample paths. Denote the set of such sample paths as $\Fc_j$, as they share the same history up to time slot $j$. We can partition $\Fc_{j}$ into two subsets:
\begin{itemize}
\item $\Fc_{j,1}$: The source skips all the upcoming arrivals and finishes transmitting the update arrives at $i$.
\item $\Fc_{j,2}$: The source switches to some later arrival.
\end{itemize}
Let $X^{\pi_0}$ be the corresponding length of the renewal epoch under policy $\pi_0$. Then, $X^{\pi_0}=i+d-1$ for sample paths in $\Fc_{j,1}$, and $X^{\pi_0}>j+d-1$ for sample paths in $\Fc_{j,2}$.

We now construct two policies $\pi_1$ and $\pi_2$ as follows. Under both $\pi_1$ and $\pi_2$, the source will behave exactly the same as under $\pi_0$ for all sample paths not in $\Fc_j$. However, for the sample paths in $\Fc_j$, the actions the source will take after $j$ will be different. Specifically, under $\pi_1$, the source will finish the update that arrives at time slot $i$ irrespective of other factors. Therefore, for all sample paths in $\Fc_j$ under $\pi_0$, the corresponding length of the renewal epoch under $\pi_1$ will be $X^{\pi_1}=i+d-1$ under $\pi_1$. For $\pi_2$, we will let the source first switch to the arrival at time slot $j$, and then switch to a later arrival whenever the source switches under $\pi_0$. Then, for the sample paths in $\Fc_{j,1}$ under $\pi_0$, the corresponding length of renewal epoch will be changed to $X^{\pi_2}=j+d-1$ under $\pi_2$; while for those in $\Fc_{j,2}$,  $X^{\pi_2}=X^{\pi_0}$.

Therefore, considering all possible sample paths under those policies, we have

$\Eb[X^{\pi_1}]<\Eb[X^{\pi_0}]<\Eb[X^{\pi_2}]$,

which implies that there must exist a $\rho$, $0<\rho<1$, such that
\begin{align}\label{eqn:E[X]}
\rho \Eb[X^{\pi_1}]+(1-\rho)\Eb[X^{\pi_2}]=\Eb[X^{\pi_0}].
\end{align}

We will then construct a randomized policy $\pi'$, under which it follows $\pi_1$ with probability $\rho$ and follows $\pi_2$ with probability $1-\rho$. Apparently, the expected length of the renewal epoch under $\pi'$, denoted as $X^{\pi'}$, will be the same as that under $\pi_0$.

Next, we will show that $\Eb[(X^{\pi'})^2]\leq \Eb[(X^{\pi_0})^2]$. Denote $P_1:=\frac{\Pb_{\pi_0}[\Fc_{j,1}]}{\Pb_{\pi_0}[\Fc_j]},P_{2}:=\frac{\Pb_{\pi_0}[\Fc_{j,2}]}{\Pb_{\pi_0}[\Fc_j]}$. Then, (\ref{eqn:E[X]}) can be expressed as
\begin{align}\label{eqn:E[X]_2}
&\rho(i+d-1)+(1-\rho)\left[(j+d-1)P_1+ \Eb[X^{\pi_0}|\Fc_{j,2}]P_{2}\right]\nonumber\\
&=(i+d-1)P_1+\Eb[X^{\pi_0}|\Fc_{j,2}]P_{2},
\end{align}
which can be reduced to
\begin{align}\label{eqn:E[X]_3}
&(1-\rho)P_{1}(j+d-1)\nonumber\\
&=(P_1-\rho) (i+d-1)+\rho P_{2}\Eb[X^{\pi_0}|\Fc_{j,2}].
\end{align}
Since $\Eb[X^{\pi_0}|\Fc_{j,2}]>j+d-1$, $(1-\rho)P_{1}=(P_1-\rho)+\rho P_{2}$, (\ref{eqn:E[X]_3}) implies that

$P_1-\rho>0$.

Dividing both sides of (\ref{eqn:E[X]_3}) by $(1-\rho)P_{1}$, we have
\begin{align*}
j+d-1&=\frac{P_1-\rho}{(1-\rho)P_{1}}(i+d-1)+\frac{\rho P_{2}}{(1-\rho)P_{1}} \Eb[X^{\pi_0}|\Fc_{j,2}].
\end{align*}
Note that $\frac{P_1-\rho}{(1-\rho)P_{1}}$ and $\frac{\rho P_{2}}{(1-\rho)P_{1}}$ form a valid distribution. Therefore, based on Jensen's inequality, we have
\begin{align}
&(j+d-1)^2\nonumber\\
& < \frac{P_1-\rho}{(1-\rho)P_{1}} (i+d-1)^2+ \frac{\rho P_{2}}{(1-\rho)P_{1}} \left(\Eb[X^{\pi_0}|\Fc_{j,2}]\right)^2\\
&\leq\frac{P_1-\rho}{(1-\rho)P_{1}} (i+d-1)^2+ \frac{\rho P_{2}}{(1-\rho)P_{1}} \Eb[\left(X^{\pi_0}\right)^2|\Fc_{j,2}],
\end{align}
which is equivalently to
\begin{align}\label{eqn:E[X^2]}
&\rho(i+d-1)^2+(1-\rho)\left[(j+d-1)^2P_1+ \Eb[\left(X^{\pi_0}\right)^2|\Fc_{j,2}]P_{2}\right]\nonumber\\
&<(i+d-1)^2P_1+\Eb[(X^{\pi_0})^2|\Fc_{j,2}]P_{2}.
\end{align}
I.e.,
\begin{align}\label{eqn:E[X^2]_2}
\rho \Eb[(X^{\pi_1})^2]+(1-\rho)\Eb[(X^{\pi_2})^2]<\Eb[(X^{\pi_0})^2].
\end{align}
Combining (\ref{eqn:E[X]}) and (\ref{eqn:E[X^2]_2}), we have
\begin{align}
\frac{1}{2}\frac{\rho \Eb[(X^{\pi_1})^2]+(1-\rho)\Eb[(X^{\pi_2})^2]}{\rho \Eb[X^{\pi_1}]+(1-\rho)\Eb[X^{\pi_2}]}< \frac{1}{2}\frac{\Eb[(X^{\pi_0})^2]}{\Eb[X^{\pi_0}]},
\end{align}
i.e., the new policy $\pi'$ achieves a lower expected average AoI than $\pi_0$, which contradicts with the assumption that $\pi_0$ is optimal.
\end{Proof}

\begin{Lemma}\label{lemma:threshold}
Under the optimal SS policy in $\Pi'$, if the source is transmitting an update that arrives at the $i$th time slot in an epoch when the new update arrives, then, there exists a threshold $\tau_i$, $\tau_i\geq i$, which depends on $i$ only, such that if the new update arrives before or at the $\tau_i$th time slot in that epoch, the source will switch to the new arrival; otherwise, it will skip the new arrival and complete the current transmission.
\end{Lemma}

\if{0}
\begin{Proof}
We prove this lemma through contradiction. Again, we focus on the first renewal epoch that starts at time 0. Assume under the optimal SS policy $\pi_0$, the transmitter transmits the update that arrives in time slot $i$. It will then skip the next arrival if it arrives at time slot $j$, and switches to it if it arrives at time slot $k$, with $j<k<i+d$. We aim to perturb $\pi_0$ a bit to obtain a new policy $\pi'$ and show that it achieves a lower expected average AoI.

Similarly to the proof of Lemma~\ref{lemma:ssp}, we consider all sample paths under $\pi_0$ that choose to transmit the update that arrives at time slot $i$, while the next arrival time is either $j$ or $k$. Denote the subset of such sample paths as $\Fc_i$. Again, we can divide $\Fc_i$ into two subsets.
\begin{itemize}
\item $\Fc_{i,1}$: The next arrival time is $j$. Then, the sample paths will skip $j$ and finish transmitting $i$. The corresponding length of renewal interval $X^{\pi_0}=i+d-1$.
\item $\Fc_{i,2}$: The next arrival time is $k$. The sample paths will switch to the new arrival, and the corresponding length of renewal interval $X^{\pi_0}\geq k+d-1$.
\end{itemize}

Then, we construct the new policy $\pi'$ in this way: For all sample paths in $\Fc_{i,1}$, the new policy will let the transmitter randomly switch to the next arrival at $j$ and finish transmitting it with probability $P_j$, and let the rest sample paths follow the original policy $\pi_0$, i.e., skip the new arrival and finish transmitting the arrival at $i$. For all sample paths in $\Fc_{a,2}$, the new policy will let the transmitter randomly skip the next arrival at $k$ with probability $P_k$, and let the rest sample paths follow the original policy $\pi_0$. Then, under $\pi'$, we have
\begin{align}
\Eb[X^{\pi'}|\Fc_{i,1}]&= P_j (j+d-1) +(1-P_j)(i+d-1)\\
\Eb[X^{\pi'}|\Fc_{i,2}]&= P_k (i+d-1) +(1-P_k)\Eb[X^{\pi_0}|\Fc_{i,2}]
\end{align}
We pick $P_j$ and $P_k$ in a way such that
\begin{align}
\Eb[X^{\pi'}|\Fc_i]&=\Eb[X^{\pi_0}|\Fc_i]  \label{eqn:lemma3-11}
\end{align}
i.e.,
\begin{align*}
P_j (j-i) \Pb[\Fc_{i,1}] &=P_k [\Eb[X^{\pi_0}|\Fc_{i,2}]- (i+d-1)] \Pb[\Fc_{i,2}]
\end{align*}
Following similar argument as in the proof of Lemma~\ref{lemma:ssp}, we can show that
\begin{align}
\Eb[(X^{\pi'})^2|\Fc_i]&\leq\Eb[(X^{\pi_0})^2|\Fc_i]   \label{eqn:lemma3-0}
\end{align}

\if{0}
Separating term $(b+d)$ from equation (\ref{eqn:lemma3-1}), we have
\begin{align}
b+d=\frac{P_b \Pb[\Fc^I_a] \hspace{-0.03in} - \hspace{-0.03in} P_c\Pb[\Fc^{II}_a]}{P_b \Pb[\Fc^I_a]}(a+d)
+ \frac{P_c\Pb[\Fc^{II}_a]}{P_b \Pb[\Fc^I_a]} \Eb[X^{\pi_0}|\Fc^{II}_a]
\end{align}
Since $\Eb[X^{\pi_0}|\Fc^{II}_a]\geq c+d$ and $a<b<c$, from equation (\ref{eqn:lemma3-1}) we must have
\begin{align}
\frac{P_c\Pb[\Fc^{II}_a]}{P_b \Pb[\Fc^I_a]} = \frac{b-a}{\Eb[X^{\pi_0}|\Fc^{II}_a]-(a+d)} \leq \frac{b-a}{c-a}
\leq 1
\end{align}
Thus $\frac{P_b \Pb[\Fc^I_a]- P_c\Pb[\Fc^{II}_a]}{P_b \Pb[\Fc^I_a]}$ and $\frac{P_c\Pb[\Fc^{II}_a]}{P_b \Pb[\Fc^I_a]}$ form a valid distribution.

From Jensen's inequality, we have
\begin{align}
(b+d)^2 &\leq \frac{P_b \Pb[\Fc^I_a]\hspace{-0.03in} -\hspace{-0.03in} P_c\Pb[\Fc^{II}_a]}{P_b \Pb[\Fc^I_a]}(a+d)^2 \nonumber \\
&\qquad \qquad + \frac{P_c\Pb[\Fc^{II}_a]}{P_b \Pb[\Fc^I_a]} \left( \Eb[X^{\pi_0}|\Fc^{II}_a] \right)^2 \\
&\leq \frac{P_b \Pb[\Fc^I_a]\hspace{-0.03in} -\hspace{-0.03in} P_c\Pb[\Fc^{II}_a]}{P_b \Pb[\Fc^I_a]}(a+d)^2 \nonumber \\
&\qquad \qquad + \frac{P_c\Pb[\Fc^{II}_a]}{P_b \Pb[\Fc^I_a]}  \Eb[(X^{\pi_0})^2|\Fc^{II}_a] \label{eqn:lemma3-2}
\end{align}
Plug (\ref{eqn:lemma3-2}) in the left hand side of (\ref{eqn:lemma3-0}), we have
\begin{align}
LHS&=\Pb[\Fc^I_a]\left( P_b(b+d)^2 + (1-P_b) (a+d)^2 \right)  \nonumber \\
&\quad + \Pb[\Fc^{II}_a]\left( P_c (a+d)^2 + (1-P_c) \Eb[(X^{\pi_0})^2|\Fc^{II}_a] \right)  \\
&\leq \Pb[\Fc^I_a] (a+d)^2 + \Pb[\Fc^{II}_a] \Eb[(X^{\pi_0})^2|\Fc^{II}_a] \\
&=RHS
\end{align}
which implies (\ref{eqn:lemma3-0}) holds.
\fi

Thus policy $\pi'$ achieves lower expected average AoI which follows from (\ref{eqn:lemma3-11}), which contradicts with the assumption that $\pi_0$ is optimal. Therefore, if the source skips the next update if it arrives at time slot $j$, it must skip it if it arrives later than $j$. Picking the minimum $j$ ($j>i$) in which a skip would occur, denote it as $j^*$. Then, the threshold $\tau_i$ equals $j^*-1$. Since the proof does not depend other information except $i$, $\tau_i$ depends on $i$ only.
\end{Proof}
\fi

\begin{Theorem}\label{thm:threshold}
Under the optimal policy in $\Pi'$, there exists a sequence of thresholds $\tau_1\geq\tau_2\geq \cdots \geq \tau_K$, such that if the source is transmitting an update that arrives in the $i$th ($i\leq K$) time slot in an renewal epoch when a new update arrives, and the arrival time of the new update is before or at the $\tau_i$th time slot in the epoch, the source will switch to the new arrival; Otherwise, if the next update arrives after $\tau_i$, or the update being transmitted arrives after $K$, the source will skip all upcoming arrivals until it finishes the current transmission.
\end{Theorem}

Theorem~\ref{thm:threshold} indicates that the optimal decision of the source only depends on two parameters: the arrival time of the update being transmitted, and the arrival time of the new update, both relative to the beginning of the renewal epoch. Therefore, the problem is essentially an MDP. In Sec.~\ref{sec:mdp}, we will cast the problem as an MDP, and numerically search for the optimal thresholds $\tau_1,\tau_2, \cdots \tau_K$ and $K$.

\if{0}
\begin{Proof}
	For successful transmitting time $d=1$, we should always switch according to Lemma \ref{lemma:idle}. Thus we consider $d>1$ in the following proof.
	We consider two scenarios.

	1) $\tau_a>a,\tau_{a+1}>a+1$.

	We prove in this case the switch thresholds has monotonicity property, i.e., $\tau_{a+1}\leq \tau_a$, $\forall a\geq 0$.

	Note that when $\tau_{a}=a+1$, the result is obvious. We consider the case when $\tau_{a}>a+1$ which implies $\tau_{a+1}>a+2$.

	Assume under optimal policy $\pi_0$ there exists two adjacent arrival times $a$ and $a+1$ such that $\tau_{a+1}>\tau_a>a+1>a$. Consider two subset of sample paths as follows.

	\begin{itemize}
		\item[I)] $\Fc_{a,\tau_a+1}$: The transmitter transmits the arrival at $a$ and the next arrival is at $\tau_a+1$. Then, under $\pi_0$, it will skip the arrival and finish transmitting the arrival at $a$. The corresponding $X^{\pi_0}=a+d$.
		\item[II)] $\Fc_{a+1,\tau_a+1}$: The transmitter transmits the arrival at $a+1$ and the next arrival is at $\tau_a+1$. Since  $\tau_a<\tau_{a+1}$, under $\pi_0$, it will switch to the new arrival. The corresponding $X^{\pi_0}\geq \tau_{a}+1+d$.
	\end{itemize}
	Then, we will modify $\pi_0$ to $\pi'$ as follows:
	For sample paths in $\Fc_{a,\tau_a+1}$, the transmitter will randomly switch to the new arrival at $\tau_a$ with probability $P_a$ and then follow policy $\pi_0$ afterwards. Note that the sample paths after switching at $\tau_a+1$ should be the same as those in $\Fc_{a+1,\tau_a+1}$.  Then, we have
	\begin{align}
	\Eb[X^{\pi'}|\Fc_{a,\tau_a+1}]&=P_a \Eb[X^{\pi_0}|\Fc_{a+1,\tau_a+1}]+(1-P_a) (a+d)
	\end{align}
	For sample paths in $\Fc_{a+1,\tau_a+1}$, the transmitter will randomly skip the new arrival at $\tau_a+1$ with probability $P_{a+1}$ and finish transmitting the arrival at $a+1$; It will stick with policy $\pi_0$ for the rest sample paths. Then
	\begin{align}
	&\Eb[X^{\pi'}|\Fc_{a+1,\tau_a+1}]=P_{a+1} (a+1+d)  \nonumber \\
	&\qquad \qquad\qquad \qquad\quad+(1-P_{a+1}) \Eb[X^{\pi_0}|\Fc_{a+1,\tau_a+1}]
	\end{align}
	We pick $P_a$ and $P_{a+1}$ in a way such that
	\begin{align}\label{eqn:mean}
	&\Eb[X^{\pi'}|\Fc_{a,\tau_a+1}]P_I+\Eb[X^{\pi'}|\Fc_{a+1,\tau_a+1}]P_{II}\\
	&=\Eb[X^{\pi_0}|\Fc_{a,\tau_a+1}]P_I+\Eb[X^{\pi_0}|\Fc_{a+1,\tau_a+1}]P_{II}
	\end{align}
	where we define $P_I:=\frac{\Pb[\Fc_{a,\tau_a+1}]}{\Pb[\Fc_{a,\tau_a+1}]+\Pb[\Fc_{a+1,\tau_a+1}]}$, $P_{II}:=\frac{\Pb[\Fc_{a+1,\tau_a+1}]}{\Pb[\Fc_{a,\tau_a+1}]+\Pb[\Fc_{a+1,\tau_a+1}]}$.
	Equation (\ref{eqn:mean}) implies that
	\begin{align}
	&P_{a+1}P_{II}(a+1+d) =(P_{a+1}P_{II}-P_aP_I)\Eb[X^{\pi_0}|\Fc_{a+1,\tau_a+1}] \nonumber \\
	&\qquad \qquad\qquad \qquad\qquad \qquad+ P_a P_I (a+d) \label{eqn:lemma4-1}
	\end{align}
	Since $a+d<b+d<\Eb[X^{\pi_0}|\Fc_{a+1,\tau_a+1}]$, from equation (\ref{eqn:lemma4-1}) we must have
	\begin{align}
	\frac{P_aP_I}{P_{a+1}P_{II}}=\frac{\Eb[X^{\pi_0}|\Fc_{a+1,\tau_a+1}]-(a+1-d)}{\Eb[X^{\pi_0}|\Fc_{a+1,\tau_a+1}]-(a+d)}<1
	\end{align}
	Thus $\frac{P_{a+1}P_{II}-P_aP_I}{P_{a+1}P_{II}}$ and $\frac{P_aP_I}{P_{a+1}P_{II}}$ form a valid distribution.
	Following similar steps as in the proof of Lemma~\ref{lemma:ssp}, we have
	\begin{align}
	&P_{a+1}P_{II}(b+d)^2 \leq(P_{a+1}P_{II}-P_aP_I)\Eb[(X^{\pi_0})^2|\Fc_{a+1,\tau_a+1}] \nonumber \\
	& \qquad\qquad\qquad\qquad\qquad+ P_a P_I (a+d)^2
	\end{align}
	which implies that
	\begin{align}\label{eqn:second_moment}
	\Eb[(X^{\pi'})^2\leq\Eb[(X^{\pi_0})^2]
	\end{align}
	Thus, $\pi_0$ cannot be optimal.

	2) $\tau_{a}=a$.

	We prove in this case $\tau_{a+1}=a+1$.

	Assume under optimal policy $\pi_0$ there exists $\tau_{a+1}>a+1$. Consider two subset of sample paths as follows.
	\begin{itemize}
		\item[I)] $\Fc_{a,a+2}$: The transmitter transmits the arrival at $a$ and the next arrival is at $a+2$. Then, under $\pi_0$, it will skip the arrival and finish transmitting the arrival at $a$. The corresponding $X^{\pi_0}=a+d$.
		\item[II)] $\Fc_{a+1,a+2}$: The transmitter transmits the arrival at $a+1$ and the next arrival is at $a+2$. Since $\tau_{a+1}>a+1$, under $\pi_0$, it will switch to the new arrival. The corresponding $X^{\pi_0}\geq a+2+d$.
	\end{itemize}
	Then, we modify $\pi_0$ to $\pi'$ as follows. For sample paths in $\Fc_{a,a+2}$, the transmitter will randomly switch to the new arrival at $a+2$ with probability $P_a$ and then follow policy $\pi_0$ afterwards. Note that the sample paths after switching at $a+2$ should be the same as those in $\Fc_{a+1,a+2}$.

	For sample paths in $\Fc_{a+1,a+2}$, the transmitter will randomly skip the new arrival at $a+2$ with probability $P_{a+1}$ and finish transmitting the arrival at $a+1$; Otherwise it will follow policy $\pi_0$ for the rest sample paths.

	Apply the same technique in 1). Pick $P_a$ and $P_{a+1}$ such that $E[X]$ does not change but $\Eb[(X^{\pi'})^2\leq\Eb[(X^{\pi_0})^2]$. Thus $\pi_0$ is not optimal.

	From 1) and 2) we conclude lemma \ref{lemma:threshold} holds.
\end{Proof}
\fi

\section{MDP based Scheduling}\label{sec:mdp}
\subsection{MDP formulation}
Motivated by the Markovian structure of the optimal policy in Theorem~\ref{thm:threshold}, we formulate the problem as an MDP as follows.

\textbf{States:} We define the state $S(t):=(\Delta(t), L(t), A(t))$, where $\Delta(t)$ and $L(t)$ are the AoI in the system, and the age of the unfinished update, at the beginning of time slot $t$, respectively. $A(t)$ is the update arrival status. \yj{Then, $\Delta(t)\geq d$, $0\leq L(t)\leq d-1$, and the state space $\Sc$ can be determined accordingly.}

\textbf{Actions:} $D(t)\in\{0,1\}$, as defined in Sec.~\ref{sec:model}.

\textbf{Transition probabilities:} The transition probability from a state $\sv:=(\delta, l, \lambda)$ to another state $\sv'$ under action $a$, denoted as $P_{\sv\sv'}(a)$, is shown in Table \ref{Tab:tranP}.

\textbf{Cost:} Let $C(S(t), D(t))$ be the immediate cost after the action $D(t)$ is taken at $t$ under state $S(t)$. We consider the instantaneous AoI after the action as the immediate cost, i.e.,
\begin{align*}
 C(S(t), D(t))&=\left\{ \begin{array}{ll}
 d& \mbox{if $L(t)=d-1, D(t)=0$};\\
 \Delta(t)+1& \mbox{otherwise}.\\
 \end{array}
 \right.
\end{align*}

In order to reduce the computational complexity, we define an approximate MDP as follows: We define $\delta_m$ as the boundary AoI, and truncate the state space of the original MDP as $\Sc_m=\{\sv\in \Sc: \delta\leq \delta_m\}$. In the transition probabilities, we bound $\delta+1$ by $\delta_m$, i.e., $[\delta+1]^+_m=\min{(\delta+1,\delta_m)}$.

Then, the optimal policy can be determined through relative value iteration as follows:
 \begin{align}\label{eqn:value iteration}
V_{n+1}(\sv)\hspace{-0.03in}=\hspace{-0.03in}\min_{a\in\{0, 1\}} \hspace{-0.03in} C(\sv, a)\hspace{-0.03in}+\hspace{-0.03in}\sum_{s'} P_{\sv\sv'}(a)V_n(\sv')\yj{-V_n(\sv_0)},
 \end{align}
where $\sv_0$ is a reference state and we set it as $\sv_0:=(d,0,0)$. For each iteration $n$, we need to update the optimal cost function for all states $\sv\in \Sc_m$ by minimizing the right hand side of (\ref{eqn:value iteration}), which causes a high computational complexity as the number of states increases. Motivated by \cite{Hsu:2017:ISIT}, we then leverage the multi-threshold structure of the optimal policy to reduce the computational complexity, as detailed in the structured value iteration algorithm in Algorithm \ref{algorithm:ssa}.

With the multiple-threshold structure, Algorithm \ref{algorithm:ssa} does not need to seek the optimal action by equation (\ref{eqn:value iteration}) for all states in each iteration as the traditional value iteration algorithm does. Specifically, if the optimal action for a state $(\delta',l,1)$ is to {\it skip} the new arrival, the optimal action for state $(\delta, l ,1)$, $\delta>\delta'$ must be to skip as well. Similarly, if the optimal action for a state $(\delta,l',1)$ is to switch to the new arrival, the optimal action for state $(\delta, l ,1)$, $l<l'$, must be to switch.

 \begin{table}[t]
 \vspace{0.05 in}
	\centering
	\resizebox{9cm}{!}{\begin{tabular}{c|c|c}
		\hline\hline
		$P_{\sv\sv'}(a)$& $a=0$ & $a=1$\\
		\hline
		$l=0$& $\begin{aligned}
			P[(\delta+1, 0, 1)|(\delta, l, \lambda)]&= p \\
			P[(\delta+1, 0, 0)|(\delta, l, \lambda)]&= 1-p \end{aligned}$
		& $\begin{aligned}
			P[(\delta+1, 0, 1)|(\delta, l, 0)] &= p\\
			P[(\delta+1, 0, 0)|(\delta, l, 0)] &= 1-p\\
			P[(\delta+1, 1, 1)|(\delta, l, 1)] &= p\\
			P[(\delta+1, 1, 0)|(\delta, l, 1)] &= 1-p\end{aligned}$\\
			\hline
				$0<l<d-1$& $\begin{aligned}
					P[(\delta+1, 0, 1)|(\delta, l, \lambda)]&= p \\
					P[(\delta+1, 0, 0)|(\delta, l, \lambda)]&= 1-p \end{aligned}$
				& $\begin{aligned}
					P[(\delta+1, 1, 1)|(\delta, l, 1)] &= p\\
					P[(\delta+1, 1, 0)|(\delta, l, 1)] &= 1-p\\
						P[(\delta+1, l+1, 1)|(\delta, l, 0)] &= p\\
						P[(\delta+1, l+1, 0)|(\delta, l, 0)] &= 1-p\end{aligned}$\\
				\hline
					$l=d-1$& $\begin{aligned}
						P[(d, 0, 1)|(\delta, l, \lambda)] &= p \\
						P[(d, 0, 0)|(\delta, l, \lambda)]&= 1-p \end{aligned}$
					& $\begin{aligned}
						P[(d, 0, 1)|(\delta, l, 0)] &= p\\
						P[(d, 0, 0)|(\delta, l, 0)] &= 1-p\\
						P[(\delta+1, 1, 1)|(\delta, l, 1)] &= p\\
						P[(\delta+1, 1, 0)|(\delta, l, 1)] &= 1-p\end{aligned}$\\
		\hline\hline
	\end{tabular}}
	\caption{Transition probabilities.}\label{Tab:tranP}

\end{table}

 \begin{algorithm}[t]
 	\caption{Structured Value Iteration.}\label{algorithm:ssa}
 	\begin{algorithmic}[1]
 		\State Initialize: $V_0(\sv)=0,\forall \sv\in \Sc_m$.
 		\For{$i=0 : n$}
 		\For{$\forall \sv\in \Sc_m$}
 		\If{$\lambda=0$}
 		\State $a^*(\sv)=0$;
 		\ElsIf{$\exists \delta'<\delta, a^*(\delta', l, 1)=0$}
 		\State{$a^*(\sv)=0$};
 		\ElsIf{$\exists l'\yj{>}l, a^*(\delta, l', 1)=1$}
 		\State{$a^*(\sv)=1$};
 		\Else
 		\State $a^*(\sv)\hspace{-0.05in}=\hspace{-0.05in}\arg\min_{a\in\{0, 1\}}\hspace{-0.04in}C(\sv,a)\hspace{-0.04in}+\hspace{-0.04in}\sum_{\sv'}\hspace{-0.04in} P_{\sv\sv'}\hspace{-0.02in}(\hspace{-0.02in}a\hspace{-0.02in})\hspace{-0.02in}V_i(\sv')$
 		\EndIf
 		\State $V_{i+1}(\sv)\hspace{-0.03in}=\hspace{-0.03in}C(\sv, a^*(\sv))+\sum_{\sv'} P_{\sv\sv'}(a)V_i(s')-V_i(\sv_0)$
        \EndFor
 		\EndFor
 		\State \textbf{return} $a^*(\sv), V(\sv).$

 	\end{algorithmic}
 \end{algorithm}

\subsection{Numerical results}
We then search for the optimal policy numerically using Algorithm~\ref{algorithm:ssa}.
We set $d=10$, $p=0.07$, and \yj{the number of iterations} to be $10,000$.\yjc{?} We set $\delta_m=50$ \yjc{too small?} for the approximate MDP. Fig.~\ref{fig:switch}(a) shows the optimal action for each state $(\Delta(t),L(t),1)$. We then plot the optimal action for each pair of arrival time of the update being transmitted and that of the new arrival in a renewal epoch in Fig.~\ref{fig:switch}(b). We note the thresholds {$\tau_1=9$, $\tau_2=8$, $\tau_3=7$, $\tau_4=6$}. They are monotonically decreasing, as predicted by Theorem~\ref{thm:threshold}. When the update being transmitted arrives after the fourth time slot in that epoch, all upcoming updates will be skipped.

\begin{figure}[t]
	\centering
	\begin{minipage}[t]{4.3cm}
		\centering
		\centerline{\includegraphics[width=4.2cm]{switch.pdf}}
		\vspace{-0.05in}
		\centerline{\small{(a)}}
		\label{fig:mdp}

	\end{minipage}
	\begin{minipage}[t]{4.2cm}
		\centering
		\centerline{\includegraphics[width=4.3cm]{switch1.pdf}}
				\vspace{-0.05in}
				\centerline{\small{(b)}}
		\label{fig:amdp}
	\end{minipage}
	\caption{The optimal policy when $p=0.07$, $d=10$. Circles represent {\it switch}, while crosses represent {\it skip}.}
	\label{fig:switch}
			\vspace{-0.1in}
\end{figure}

Then, we compare the average AoI under the optimal policy identified by Algorithm~\ref{algorithm:ssa} and a myopic policy over $10,000$ time slots. Under the myopic policy, the source will never switch to a new update arrival until it finishes the one being transmitted. The performance gap is plotted in Fig. \ref{fig:ageResults2}. As we observe, the optimal policy always outperforms the myopic policy. Although the greedy policy minimizes the length of the each epoch greedily, it does not render the minimum average AoI. This is because $X_i$ has a larger second moment in this case, leading to higher AoI. We note that when $p$ gets sufficiently small or large, the performance gap between both policies becomes close to zero. This is because for such extreme cases, the multiple-threshold policy and the myopic policy become identical to each other.

\if{0}
\begin{figure}[t]
	\centering
	\centerline{\includegraphics[width=8cm]{./graph/result1.eps}}
	\caption{Average age under different transmitting time.}
	\label{fig:ageResults1}
\end{figure}
\fi

\begin{figure}[t]
	\centering
	\includegraphics[width=7cm]{result3.pdf}
		\vspace{-0.1in}
	\caption{Performance gap between the optimal multiple-threshold policy and a myopic policy.}
	\label{fig:ageResults2}
	\vspace{-0.1in}
\end{figure}

\if{0}
Tab. \ref{Tab:runningtime} shows the average running times of the algorithms. We note that the SWS has significant improvement over the traditional value iteration, as the approximate MDP reduces the state spaces and the reduced iteration according to the sequential switching feature in the SWS.

\begin{table}[h]
	\begin{center}
		\caption{Average running times in seconds.}\label{Tab:runningtime}
		{\begin{tabular}{r||c| c| c}
				\hline
				      & T=100  &T=1000 & T=10000\\
				[0.5ex]
				\hline\hline
				VI               & 0.1258 & 1.41      & 12.28    \\
				SWS               & 0.045& 0.056     & 0.079        \\
			 \hline
			\end{tabular}}
		\end{center}
		\vspace{-0.2in}
	\end{table}
\fi

\bibliographystyle{IEEEtran}
\bibliography{IEEEabrv,AgeInfo,ener_harv}

\end{document}


TITLE: Analytic family of functions and Lipschitz continuity
QUESTION [2 upvotes]: Let $f_t:U_t\rightarrow V_t$ be an analytic family of homeomorphisms, where $0\leq t\leq 1$, and $U_t, V_t\subset \mathbb{R}^n$ open and bounded domains. Moreover, there exist $c,C>0$ such that for all x,y,t:
$$c\cdot \vert x-y \vert\leq \vert f_t(x)-f_t(y) \vert\leq C\cdot \vert x-y \vert.$$
I want to show that all derivatives $\frac{d^n}{dt^n} f _{\vert t=t_0}:U_{t_0}\rightarrow \mathbb{R}^n$ are Lipschitz continuous for all $n\in\mathbb{N}$.
Unfortunately, I do not know how to prove this claim. Since $f$ is analytic, there is an interval $I_{\epsilon}:=(t_0-\epsilon, t_0+\epsilon)$ such that for all $t\in I_{\epsilon}$:
$$f_t(x)=\sum_{n=0}^{\infty} \frac{d^n}{dt^n} f _{\vert t=t_0}(x)\cdot (t-t_0)^n.$$
I think this might help.But I do not know exactly what to do next. If you know how to prove the result, please help me.
Best regards

REPLY [1 votes]: This is false in general. Let $U_t=V_t=\mathbb{R}^2$ for all $t$.  For $0\le t\le 1$  define $f_t$ by
$$
f_t(x_1, x_2) = (x_1, x_2 + \sin (tx_1))
$$
This is an analytic function of $x_1, x_2, t$. It is Lipschitz continuous with respect to $(x_1, x_2)$ with Lipschitz constant at most $3$. The same holds for its inverse, which is
$$
f_t^{-1}(x_1, x_2) = (x_1, x_2 - \sin (tx_1))
$$
However,
$$
\frac{d}{dt}f_t(x_1, x_2) = (x_1, x_2 + x_1\cos (tx_1))
$$
is not Lipschitz continuous, as its derivative with respect to $x_1$ is unbounded.