TITLE: Why is $\int x^2e^{x^2}\,d(x^2)$ not proper notation?
QUESTION [4 upvotes]: So, as the title states. My question is about mathematical notation. Whenever my lecturer in calculus one computes an integral they do some in my opinion very redundant steps.
When they attempt to compute for example $\int x^3e^{x^2}\,dx$, they will write down $u = x^2; du=2xdx$ and proceed to turn the integral into $\frac{1}{2}\int ue^u\,du$.
I don't see why it is required to hide this simple mathematical truth that $d(x^2)=2xdx$ behind a substitution, yet they tells me that not doing so would technically be incorrect.
By the way, I am not asking this to proof my lecturer wrong, since they are an epic person and when I told them that I found it easier without the substitution I was given permision to do so. I am merely asking why it would be "technically incorrect" to do so.

REPLY [3 votes]: The notation is proper.  The Riemann-Sieltjes for a real-valued function $f(x)$ for $x\in [a,b]$ with respect to a real-valued function $g(x)$ is written
$$I=\int_a^b f(x) dg(x)$$
If $g$ is differentiable on $[a,b]$, then $I$ can be written as
$$I=\int_a^b f(x)\frac{dg(x)}{dx}\,dx$$
However, $g(x)$ need not be differentiable and, in fact, can even be discontinuous.
In the case of interest, $g(x)=x^2$ and $f(x)=x^2e^{x^2}$.  We have, therefore
$$\int_a^b x^2e^{x^2}d(x^2)=2\int_a^b x^3 e^{x^2}\,dx$$


text:
\begin{document}
\nolinenumbers

\maketitle

\begin{abstract}
  We propose using greedy and randomized Kaczmarz inner-iterations as preconditioners for the right preconditioned flexible GMRES method to solve consistent linear systems, with a parameter tuning strategy for adjusting the number of inner iterations and the relaxation parameter.
  We also present theoretical justifications of the right-preconditioned flexible GMRES for solving consistent linear systems. Numerical experiments on overdetermined and underdetermined linear systems show that the proposed method is superior to the GMRES method preconditioned by NE-SOR inner iterations in terms of total CPU time. \blfootnote{\textbf{Funding}: This work was funded by JSPS KAKENHI Grant (No.~15K04768 and No.~20K14356), the National Natural Science Foundation of China (No.~11971354) and the China Scholarship Council (No.~201906260146).}
\end{abstract}

\begin{keywords}
  Kaczmarz method, randomized algorithm, linear system, overdetermined system, underdetermined system, least squares problem, iterative method, inner-outer iteration, preconditioner, GMRES, flexible GMRES
\end{keywords}

\begin{AMS}
  65F08, 65F10, 65F50, 15A06
\end{AMS}

\section{Introduction}
Consider solving consistent linear systems
\begin{equation}\label{l1}
	 A x = b, \quad b \in {\mathcal{R}({A})},
\end{equation}
where $A \in {\mathbb{R}^{m \times n}}$ is not necessarily of full rank and ${\mathcal{R}({A})}$ is the range space of $A$.
In particular, consider the minimum Euclidean-norm solution
\begin{equation}\label{l2}
	\operatorname*{min}_{x \in {\mathbb{R}^{n}}} \| x \| _ { 2 }
  \quad \text { s.t. } \quad A x = b, \quad b \in {\mathcal{R}({A})}.
\end{equation}
The problem~\eqref{l2} is equivalent to the normal equations of the second kind
\begin{equation}\label{n2}
	A{A^\mathsf{T}}u = b, \quad  {x = }{A^\mathsf{T}}u, \quad b \in {\mathcal{R}({A})},
\end{equation}
where $(\cdot)^\mathsf{T}$ denotes the transpose.

Direct methods for solving problem~\eqref{l2} or~\eqref{n2} are generally expensive when the coefficient matrix is large and sparse. A well-established iterative method for solving~\eqref{l2} is the (preconditioned) CGNE method~\cite{C55,S03}, which is mathematically equivalent to the (preconditioned) Conjugate Gradient (CG) method~\cite{HS52} applied to~\eqref{n2}.
Another method is the (preconditioned) MRNE method~\cite{MH15, CMTH19}, which applies the (preconditioned) MINRES method~\cite{PS75} to~\eqref{n2}.
Note that iterative methods may be slow to converge for ill-conditioned problems since the condition number of $AA^\mathsf{T}$ is the square of that of $A$, and preconditioning becomes necessary.
In~\cite{HYI10}, Hayami, Yin and Ito proposed a
right-preconditioned generalized minimal residual (GMRES) method called the AB-GMRES method by applying GMRES to $\operatorname*{min}_{u \in {\mathbb{R}^{m}}} \| b-ABu \| _ { 2 }$, where $B \in {\mathbb{R}^{n \times m}}$.

In order to accelerate the convergence of iterative methods and save the storage requirement, inner iterations can be applied as a preconditioner inside the Krylov subspace methods instead of applying preconditioning matrices explicitly. Such techniques are often called inner-outer iteration methods~\cite{S93}. Morikuni and Hayami~\cite{MH13,MH15} proposed a class of inner-iteration Krylov subspace methods by applying stationary inner iterations as implicit preconditioners, and showed their efficiency particularly for ill-conditioned and rank-deficient problems. (See also~\cite{CMTH19}.)

In AB-GMRES, common choices for stationary inner iterations are the
normal error Gauss-Seidel (NE-GS) and normal error successive overrelaxation (NE-SOR) methods~\cite{BE79,S03}, which are also commonly referred to as the Kaczmarz~\cite{K37} or row action methods~\cite{A84,BL13,HY12,C81,G03,S40}. Since it was proposed in the 1930s, the Kaczmarz method has gained great theoretical development and plentiful practical applications~\cite{B03,C88,E40,E10,PCB12,T71}. Research on the Kaczmarz method was reignited in 2006 and 2009 when Strohmer and Vershymin~\cite{SV06,SV09} proposed the randomized Kaczmarz method with expected exponential rate of convergence. In~\cite{BW18}, Bai and Wu constructed a greedy randomized Kaczmarz method by proposing a more effective probability criterion. In~\cite{P18}, Popa summarized convergence rates for Kaczmarz-type methods, including greedy Kaczmarz~\cite{A84}, randomized Kaczmarz methods~\cite{SV09} and so on.
For more literature on Kaczmarz-type methods, we refer the reader to~\cite{BWoc18,BWA18,BWO18,BW19}.
Numerical results show that these randomized or greedy Kaczmarz-type methods accelerate the original Kaczmarz method and reduce the required number of iterations and CPU time effectively. Inspired by this randomized framework, we replace the NE-SOR method by the greedy and randomized Kaczmarz methods in AB-GMRES preconditioned by stationary inner iterations.

A motivation for developing such a greedy/randomized inner-iteration preconditioning arises in applications where the operation on a row of a matrix is relatively expensive, such as in basis pursuit problems~\cite{CDM01,CMTH19}. We mention related previous work~\cite{AMT10, MS14, RT08} on randomized preconditioners for least squares problems.

When the randomized or greedy Kaczmarz method is used as the inner iteration, the rows of $A$ are selected randomly or greedily in each iteration and thus the preconditioner is not fixed during the outer iteration. Therefore, we use the flexible GMRES method~\cite{S93} as the outer iteration and propose a new algorithm called the flexible AB-GMRES method with Kaczmarz-type inner iterations.

Theoretically, an optimality property of minimizing residual norm can be given under the framework of flexible GMRES. We also propose a parameter tuning procedure for adjusting the number of inner iterations and the relaxation parameter for the new method. Numerical results show that flexible AB-GMRES preconditioned by Kaczmarz-type methods outperform the AB-GMRES method preconditioned by NE-SOR iterations~\cite{MH13,MH15} in terms of total CPU time.

The organization of this paper is as follows. In section 2, we review the AB-GMRES method. In section 3, we present the flexible AB-GMRES method for consistent linear systems, and give an optimality property of the proposed method. In section 4, we propose a parameter tuning procedure for the new method and present numerical experiment results. In section 5, we conclude the paper.

\section{AB-GMRES method}
In this section, the inner-iteration preconditioned AB-GMRES method is briefly introduced.
Consider solving equation~\eqref{l2} using AB-GMRES. AB-GMRES corresponds to GMRES applied to $\mathop {\min }_{u \in {\mathbb{R}^m}} {\left\| {b - ABu} \right\|_2}$ with $x = Bu$ and works in an $m$-dimensional space~\cite{HYI10}.
In order to achieve fast convergence of AB-GMRES and to avoid storing the preconditioner $B$, stationary inner iterations in combination with AB-GMRES were proposed in~\cite{MH13,MH15}. This algorithm can be described as follows.

\begin{algorithm}[H]
	\caption{AB-GMRES method preconditioned by inner iterations~\cite{MH15}}
	\label{ABPI}
	\begin{algorithmic}[1]
		\STATE Let $x_0$ be the initial approximate solution and $r_0=b-Ax_0$.
		\STATE $\beta  = {\left\| {{{r}_0}} \right\|_2}$, ${v_1} = {{r}_0}/\beta$
		\FOR{$k = 1,2, \ldots$ until convergence}
		\STATE Apply $\ell$ iterations of a stationary iterative method to $AA^\mathsf{T}y=v_k$, $z=A^\mathsf{T}y$ to obtain $z_k=B^{(\ell)}v_k$.
		\STATE ${w_k} = A{z_k}$
		\FOR{$i = 1,2, \ldots ,k,$}
		\STATE ${h_{i,k}} = {w_k}^{\mathsf{T}}{v_i}$, ${w_k} = {w_k} - {h_{i,k}}{v_i}$
		\ENDFOR
		\STATE ${h_{k + 1,k}} = {\left\| {{w_k}} \right\|_2}$, ${v_{k + 1}} = {w_k}/{h_{k + 1,k}}$
		\ENDFOR
		\STATE ${y_k} \!=\! \arg \mathop {\min }_{y \in {\mathbb{R}^k}} {\left\| {\beta {e_1} \!-\! {{\bar H}_k}y} \right\|_2}$, ${u_k} \!=\! \left[ {{v_1},{v_2}, \ldots ,{v_k}} \right]{y_k}$, where ${{\bar H}_k} \!=\! {\left\{ {{h_{ij}}} \right\}_{1 \leq i \leq k + 1;1 \leq j \leq k}}$
		\STATE Apply $\ell$ iterations of a stationary iterative method to $AA^\mathsf{T}y=u_k$, $z=A^\mathsf{T}y$ to obtain $z_k=B^{(\ell)}u_k$. \\
		\STATE $x_k=x_0+z_k$
	\end{algorithmic}
\end{algorithm}

In the AB-GMRES preconditioned by inner iterations, one common choice for stationary inner iterations is the NE-SOR method, which is mathematically equivalent to the SOR method applied to the normal equations of the second kind~\cite{BE79,S03}.
More specifically, if we use ${\alpha _i^\mathsf{T}}$ to represent the $i$th row of the matrix $A$, and $v_i$ to represent the $i$th entry of the vector $v$, then the NE-SOR method for $AA^\mathsf{T}y=v, z=A^\mathsf{T}y$ can be described as follows.

\begin{algorithm}[H]
	\caption{NE-SOR method~\cite{S03}, Kaczmarz method with relaxation}
	\label{NE-SOR}
	\begin{algorithmic}[1]
		\STATE Let $z^{(0)}$ be the initial approximate solution and $\omega \in {\mathbb{R}}$ be the relaxation parameter.  \\
		\FOR{$p = 0,1,2, \ldots, \ell-1$}
		\FOR{$i = 1,2, \ldots, m$}
		\STATE ${z^{(p)}} = {z^{(p)}} + \omega \frac{{{v_i} - \alpha _i^\mathsf{T}{z^{(p)}}}}{{\left\| {{\alpha _i}} \right\|_2^2}}{\alpha _i}$
		\ENDFOR
		\STATE ${z^{(p + 1)}} = {z^{(p)}}$
		\ENDFOR
	\end{algorithmic}
\end{algorithm}

The Kaczmarz method~\cite{K37} is equivalent to Algorithm~\ref{NE-SOR} with $\omega=1$.
In fact, the iteration scheme of NE-GS (NE-SOR) is exactly the same as that of the Kaczmarz method (relaxed Kaczmarz method)~\cite{K37}.
The relaxed Kaczmarz (NE-SOR) method is one of the most efficient row action methods. The method cycles through the rows of the linear system and form each iterate by orthogonally projecting the current point onto the hyperplane $\alpha _{i}^\mathsf{T}z^{(p+1)}=v_{i}$ formed by the active row, and all the $m$ equations in the linear system are swept through consecutively in $m$ iterations.

The convergence theorem of AB-GMRES preconditioned by NE-SOR
is precisely restated below.

\begin{theorem}~\cite[Theorem~5.6]{MH15}.
	AB-GMRES preconditioned by NE-SOR inner iterations with $0<\omega<2$, determines the minimum-norm solution of $Ax=b$ without breakdown for all $b \in \mathcal{R}(A)$ and for all ${x_0} \in \mathcal{R}(A^\mathsf{T})$.
\end{theorem}

\subsection{Flexible AB-GMRES method}
We adopt the flexible preconditioning framework proposed in FGMRES~\cite{S93} to AB-GMRES, and consider using Kaczmarz-type inner iterations in it.

\subsubsection{Outer-iteration algorithm}
It is well known that the preconditioning matrix needs to be fixed in the preconditioned GMRES method for all the outer iterations.
In fact, in order to keep the preconditioner ${B^{(\ell)}}$ in Algorithm~\ref{ABPI} fixed, the number of inner iterations in AB-GMRES should not be changed for each outer iteration.
However, if we were to adopt the randomized or greedy algorithm as the preconditioner in AB-GMRES, the preconditioning matrix for each outer iteration may change even though the number of inner iterations for each outer iteration is fixed.
In~\cite{S93}, Saad presented a variant of the GMRES algorithm called flexible GMRES (FGMRES) for solving square linear systems $Ax=b$, which allows the preconditioning matrix to change for each outer iteration. This algorithm can be described as follows.
	\begin{algorithm}[H]
		\caption{Flexible GMRES (FGMRES)~\cite{S93}}
		\label{Flexible GMRES(FGMRES)}
		\begin{algorithmic}[1]
			\STATE Let $x_0$ be the initial approximate solution and $r_0=b-Ax_0$.
			\STATE $\beta  = {\left\| {{r_0}} \right\|_2}$, ${v_1} = {r_0}/\beta$
			\FOR{$k = 1,2, \ldots$ until convergence}
			\STATE $z_k=B^{(k)}{v_k}$, ${w_k} = Az_k$
			\FOR{$i = 1,2, \ldots ,k,$}
			\STATE ${h_{i,k}} = ({w_k },{v_i})$, ${w_k } = {w_k } - {h_{i,k}}{v_i}$
			\ENDFOR
			\STATE ${h_{k + 1,k}} = {\left\| {{w_k }} \right\|_2}$, ${v_{k + 1}} = {w_k }/{h_{k + 1,k}}$
			\STATE Define ${Z_k} = \left[ {{z_1}, {z_2}, \ldots ,{z_k}} \right]$.
			\ENDFOR
			\STATE ${y_k} = \arg \mathop {\min }_{y \in {\mathbb{R}^k}} {\left\| {\beta {e_1} - {{\bar H}_k}y} \right\|_2}$, where ${{\bar H}_k} = {\left\{ {{h_{ij}}} \right\}_{1 \leq i \leq k + 1;1 \leq j \leq k}}$
			\STATE ${x_k} = {x_0} + Z_k{y_k}$
		\end{algorithmic}
	\end{algorithm}

\subsubsection{Randomized inner-iteration algorithms}

Instead of using the rows of the coefficient matrix $A$ in the consecutive order, Ansorge~\cite{A84} proposed selecting the $i_k$th row corresponding to the residual component with maximum absolute value. We will call this method the greedy Kaczmarz method, or GK method for short.
See also~\cite{S40} for earlier work.
Another approach is to choose the $i_k$th row randomly~\cite{SV09}, which is called the randomized Kaczmarz (RK) method.
Finally, Bai and Wu~\cite{BW18} combined these two ideas to propose a more effective probability criterion, which aims to
diminish entries of the residual vector with relatively large absolute value at each iteration. We will call the corresponding algorithm the greedy randomized Kaczmarz (GRK) method. These algorithms are called Kaczmarz-type methods. It can be proved that GK, RK and GRK methods converge to the minimum-norm solution whether the system is overdetermined or underdetermined and the coefficient matrix has full rank or is rank deficient~\cite{P18}. The algorithms of GK, RK and GRK are given below.

\begin{algorithm}[H]
	\caption{Greedy Kaczmarz (GK) method~\cite{A84}}
	\label{The Greedy Kaczmarz Method}
	\begin{algorithmic}[1]
		\STATE Let ${z^{(0)}}$ be the initial approximate solution.
		\FOR{$p = 0,1,2, \ldots, \ell-1$}
		\STATE Select ${i_p} \in \left\{ {1,2, \ldots m} \right\}$ according to ${i_p} = \arg \mathop {\max }\limits_i  {{{| {v_i} - \alpha _i^\mathsf{T}{z^{(p)}} |}}}$.
		\STATE ${z^{(p + 1)}} = {z^{(p)}} + \omega \frac{{{v_{{i_p}}} - \alpha _{{i_p}}^\mathsf{T}{z^{(p)}}}}{{\left\| {{\alpha _{{i_p}}}} \right\|_2^2}}{\alpha _{{i_p}}}$
		\ENDFOR
	\end{algorithmic}
\end{algorithm}

\begin{algorithm}[H]
	\caption{Randomized Kaczmarz (RK) method~\cite{SV09}} \label{RKMM}
	\begin{algorithmic}[1]
		\STATE Let ${z^{(0)}}$ be the initial approximate solution.
		\FOR{$p = 0,1,2, \ldots, \ell-1$}
		\STATE Select ${i_p} \in \left\{ {1,2, \ldots m} \right\}$ with probability $\Pr (\text{row} = {i_p}) = \frac{{\left\| {{\alpha _{{i_p}}}} \right\|_2^{2}}}{{\left\| A \right\|_\text{F}^{2}}}$.
		\STATE ${z^{(p + 1)}} = {z^{(p)}} + \omega \frac{{{v_{{i_p}}} - \alpha _{{i_p}}^\mathsf{T}{z^{(p)}}}}{{\left\| {{\alpha _{{i_p}}}} \right\|_2^2}}{\alpha _{{i_p}}}$
		\ENDFOR
	\end{algorithmic}
\end{algorithm}

\begin{algorithm}[H]
	\caption{Greedy Randomized Kaczmarz (GRK) method~\cite{BW18}} \label{GRKMM}
	\begin{algorithmic}[1]
		\STATE Let ${z^{(0)}}$ be the initial approximate solution.
		\FOR{$p = 0,1,2, \ldots, \ell-1$}
		\STATE ${\varepsilon _p} = \frac{1}{2}\left( {\frac{1}{{\left\| {v - A{z^{(p)}}} \right\|_2^2}}\mathop {\max }\limits_{1 \le {i_p} \le m} \left\{ {\frac{{{{\left| {{v_{{i_p}}} - \alpha _{{i_p}}^\mathsf{T}{z^{(p)}}} \right|}^2}}}{{\left\| {{\alpha _{{i_p}}}} \right\|_2^{\rm{2}}}}} \right\} + \frac{1}{{\left\| A \right\|_\mathsf{F}^{\rm{2}}}}} \right)$
		\STATE Determine the index set of positive integers\\
		$${U_p} = \left\{ { {{i_p}} : {{{| {{v_{{i_p}}} - \alpha _{{i_p}}^\mathsf{T}{z^{(p)}}} |}^2}} \ge {\varepsilon _p}\| {v - A{z^{(p)}}} \|_2^2}{{\| {{\alpha _{{i_p}}}} \|_2^{2}}} \right\}.$$
		\STATE Compute the $i$th entry $\tilde s_i^{(p)}$ of the vector ${{\tilde s}^{(p)}}$ according to\\
		\begin{equation*}
		\tilde s_i^{(p)}=
		\begin{cases}
		{v_i} - \alpha _i^\mathsf{T}{z^{(p)}},& \text{if} \quad i \in {U_p},\\
		0,& \text{otherwise.}
		\end{cases}
		\end{equation*}
		\STATE Select ${i_p} \in {U_p}$ with probability $\Pr (\text{row} = {i_p}) = \frac{{{{\left| {\tilde s_{i_p}^{(p)}} \right|}^2}}}{{\left\| {{{\tilde s}^{(p)}}} \right\|_2^{\text{2}}}}$.
		\STATE ${z^{(p + 1)}} = {z^{(p)}} + \omega \frac{{{v_{{i_p}}} - \alpha _{{i_p}}^\mathsf{T}{z^{(p)}}}}{{\left\| {{\alpha _{{i_p}}}} \right\|_2^2}}{\alpha _{{i_p}}}$
		\ENDFOR
	\end{algorithmic}
\end{algorithm}

In Algorithms~\ref{RKMM} and~\ref{GRKMM}, Pr($\text{row}=i$) represents the probability of selecting the $i$th row of the matrix $A$ as the working row of this iteration.

We remark that $s=v-Az$ needs to be calculated at each inner iteration for GK and GRK methods. This additional computational work cannot be ignored. On the other hand, we may update the residual vector $s$ recursively as follows~\cite{BW18}:
\begin{align}\label{r}
{s^{(p+1)}} &= v - A{z^{(p+1)}}\nonumber\\
&=v - A\left( {{z^{(p)}} + \frac{{{v_{{i_p}}} - \alpha _{{i_p}}^\mathsf{T}{z^{(p)}}}}{{\left\| {{\alpha _{{i_p}}}} \right\|_2^2}}\alpha _{{i_p}}} \right)\nonumber\\
&= v - A{z^{(p)}} - \frac{{s_{{i_p}}^{\left( p \right)}}}{{\left\| {{\alpha _{{i_p}}}} \right\|_2^{\rm{2}}}}A\alpha _{{i_p}}\nonumber\\
&={s^{\left( p \right)}}-\frac{{s_{{i_p}}^{\left( p \right)}}}{{\left\| {{\alpha _{{i_p}}}} \right\|_2^{\rm{2}}}}{C_{({i_p})}}.
\end{align}
Here, ${C_{({i_p})}}$ is the ${i_p}$th column of $C=AA^\mathsf{T}$.
Hence, if the matrix product $AA^\mathsf{T}$ is computed and stored once in the beginning, the computational work can be reduced, assuming that the total number of inner iterations is more than the number of rows of $A$. (See Appendix.) This condition was satisfied in all our numerical experiments.

\section{Flexible AB-GMRES preconditioned by Kaczmarz-type methods as inner iterations}
In FGMRES,  we can change the preconditioner for each outer iteration.
Hence, consider using $\ell_{k}$ iterations of a Kaczmarz-type method as the preconditioner for each outer iteration of the flexible AB-GMRES (F-AB-GMRES). We denote the preconditioning matrix given by the $\ell_{k}$ inner iterations by ${B^{(\ell_{k})}}$. The algorithm is given as follows.
\begin{algorithm}[H]
		\caption{Flexible AB-GMRES preconditioned by Kaczmarz-type methods}
		\label{AGPK}
		\begin{algorithmic}[1]
			\STATE Let $x_0$ be the initial approximate solution and $r_0=b-Ax_0$.
			\STATE $\beta  = {\left\| {{{r}_0}} \right\|_2}$, ${v_1} = {{r}_0}/\beta$
			\FOR{$k = 1,2, \ldots$ until convergence}
			\STATE Apply $\ell_{k}$ iterations of a Kaczmarz-type method to $Az\!=\!v_k$ to obtain $z_k\!=\!{B^{(\ell_{k})}}v_k$, where $\ell_{k}$ is the smaller of the maximum number of inner iterations $\ell_\text{max}$ and the smallest $\ell$ such that\\
			$${\left\| {v_k - A{B^{(\ell)}}v_k} \right\|_2} \le \max \left\{ {{{0.9}^k},{{10}^{ - 1}}} \right\} {\left\| v_k \right\|_2}.$$
			\STATE ${w_k} = A{z_k}$
			\FOR{$i = 1,2, \ldots ,k,$}
			\STATE ${h_{i,k}} = {w_k}^\mathsf{T}{v_i}$, ${w_k} = {w_k} - {h_{i,k}}{v_i}$
			\ENDFOR
			\STATE ${h_{k + 1,k}} = {\left\| {{w_k}} \right\|_2}$, ${v_{k + 1}} = {w_k}/{h_{k + 1,k}}$
			\ENDFOR
			\STATE ${y_k} \!=\! \arg \mathop {\min }_{y \in {\mathbb{R}^k}} {\left\| {\beta {e_1} - {{\bar H}_k}y} \right\|_2}$, ${u_k} \!=\! \left[ {{z_1},{z_2}, \ldots ,{z_k}} \right]{y_k}$, where ${{\bar H}_k} \!=\! {\left\{ {{h_{ij}}} \right\}_{1 \leq i \leq k + 1;1 \leq j \leq k}}$
			\STATE $x_k=x_0+u_k$
		\end{algorithmic}
\end{algorithm}
Here, ${{\bar H}_k} = \left\{ {{h_{i,j}}} \right\} \in {\mathbb{R}^{(k + 1) \times k}}$. Since the number of inner iterations in each outer iteration does not have to be fixed, we proposed a new inner iterations stopping criterion that varies with the outer iteration to accelerate the convergence, which is given in line 4 of Algorithm~\ref{AGPK}.
An optimality property similar to GMRES is given as in~\cite[Proposition~2.1]{S93} for FGMRES.

\begin{theorem}\label{the1}
	The approximate solution $x_k$ obtained at the $k$th iteration of F-AB-GMRES minimizes the residual norm ${\left\| {b - A{x_k}} \right\|_2}$ over ${x_0} + \mathrm{Span}\left\{ {z_1},{z_2}, \ldots ,{z_k} \right\}$.
\end{theorem}
\begin{proof}
It is easy to show that

\begin{align} \label{z}
A{Z_k} &= {V_{k + 1}}{{\bar H}_k} \nonumber \\
&= {V_k}{H_k} + {h_{k + 1,k}}{v_{k + 1}}e_k^\mathsf{T}.
\end{align}

Here, $H_k$ denotes the $k \times k$ matrix obtained from ${{\bar H}_k}$ by deleting its last row, $V_{k}$ denotes the $n \times k$ matrix with column vectors $v_1$, $v_2$, $\ldots$, $v_k$, and $e_k$ denotes the last column of the identify matrix $\mathrm{I_k}$ of size $k$.

Consider the residual vector for an arbitrary vector $x = x_0 +{Z_k}y$ in the affine space $x_0 + \mathrm{Span}\left\{ {{Z_k}} \right\}$, where ${Z_k} = [{z_1},{z_2}, \ldots ,{z_k}]$. The optimality property is based on the relations

\begin{align}\label{b1}
b-A x &=b-A\left(x_{0}+Z_{k} y\right) \nonumber \\
&=r_{0}-A Z_{k} y \nonumber \\
&=\beta v_{1}-V_{k+1} \bar{H}_{k} y \nonumber \\
&=V_{k+1}\left(\beta e_1-\bar{H}_{k} y\right),
\end{align}
where $e_1$ is the first column of $\mathrm{I_{k+1}}$.

If $J_k(y)$ denotes the function
\begin{equation}
{J_k}(y) = {\left\| {b - A\left( {{x_0} + {Z_k}y} \right)} \right\|_2},
\end{equation}
observe that by~\eqref{b1} and the fact that the column vectors of
$V_{k+1}$ are orthogonal,
\begin{equation}
{J_k}(y) =  {\left\| {\beta {e_1} - {{\bar H}_k}y} \right\|_2}.
\end{equation}
Since the algorithm minimizes this norm over all vectors in ${\mathbb{R}^{k}}$ to yield $y_k$, it is clear that the approximate solution $x_k = x_0 +Z_ky_k$ has the smallest residual norm in $x_0 + \mathrm{Span}\left\{ {{Z_k}} \right\}$. \end{proof}

Next, we consider the possibility of breakdown in F-AB-GMRES preconditioned by Kaczmarz-type methods as inner iterations.
A breakdown occurs when the vector $v_{k+1}$ cannot be computed in line 9 of Algorithm~\ref{AGPK} because $h_{k+1, k} = 0$. For AB-GMRES with $B$ satisfying the convergence conditions~\cite[Theorem~5.2]{MH15},~\cite[Corollary~3.8]{HYI10}, this is not a problem because when this happens, the approximate solution $x_k$ satisfies $Ax_k=b$.
The situation for F-AB-GMRES is slightly different, as in~\cite[Proposition~2.2]{S93} for FGMRES.

\begin{theorem}\label{the2}
	Assume that $\beta=\left\|r_{0}\right\|_{2} \neq 0$ and that $k-1$ steps of F-AB-GMRES have been successfully performed, i.e., that $h_{i+1, i} \neq 0$ for $i < k$. In addition, assume that the matrix $H_k$ is nonsingular. Then $x_k$ is a solution of $Ax=b$ if and only if $h_{k+1, k} = 0$.
\end{theorem}
\begin{proof}
If $h_{k+1, k} = 0$, then $A Z_{k}=V_{k} H_{k}$, and we have
$$J_{k}(y)=\left\|\beta v_{1}-A Z_{k} y \right\|_{2}=\left\|\beta v_{1}-V_{k} H_{k} y \right\|_{2}=\left\|\beta e_{1}-H_{k} y \right\|_{2},$$
where $e_1$ is the first column of $\mathrm{I_k}$.

If $H_k$ is nonsingular, then the above function is minimized for $y=y_k=H_{k}^{-1}\left(\beta e_{1}\right)$ and
the corresponding minimum norm is zero, that is to say, $x_k$ is a solution of $Ax=b$.

On the other hand, if $x_k$ satisfies $Ax_k=b$, then from~\eqref{z} and~\eqref{b1}, we have
\begin{equation}\label{r2}
	0=b-A x_{k}=V_{k}\left(\beta e_{1}-H_{k} y_{k}\right)- w_k e_{k}^{\mathsf{T}} y_{k}.
\end{equation}
We show that $h_{k + 1,k}=0$ by contradiction. Suppose that $h_{k + 1,k} \neq 0$ holds. If the last component of $y_k$ is equal to zero $e_{k}^{\mathsf{T}} y_{k}=0$, then a simple backward substitution for the system $H_{k} y_{k}=\beta e_{1}$, starting from the last equation, will show that all components of $y_k$ are zero. This would imply that $\beta=0$, which contradicts the assumption. Hence, $h_{k + 1,k}=0$ holds.
\end{proof}

The only difference of this result from that of AB-GMRES is that the additional assumption that $H_k$ is nonsingular must be made since it is no longer implied by the algorithm.
In fact, the following holds

\begin{theorem}
Assume $\rank A = n$, $\rank Z_k =k$ and $h_{k+1,k} = 0.$ Then, $H_k$ is nonsingular.
\end{theorem}

\begin{proof}
Let $c_1 A z_1 + \cdots + c_k A z_k = A(c_1 z_1 + \cdots + c_k z_k ) = 0$. Then, $ \rank A = n$ implies $c_1 z_1 + \cdots + c_k z_k = 0$, and $ \rank Z_k = k $ implies $c_1 = \cdots =c_k = 0$. Hence, $ \rank(AZ_k)=k$. Since, $h_{k+1,k} = 0$, $AZ_k = V_k H_k$. Thus, $k=\rank(AZ_k)=\rank(V_k H_k) \le \min(\rank V_k, \rank H_k ) = \min (k, \rank H_k)$. Hence, $\rank H_k = k$, and $H_k$ is nonsingular.
\end{proof}

See also~\cite[Theorem~3]{M17} for convergence conditions of FGMRES preconditioned by multistep matrix splitting iterations, which ensure the nonsingularity of $H_k$.

The additional cost of the flexible variant over AB-GMRES is only the extra memory required to save the set of vectors $\{z_j\}_{j=1,2,\ldots,m}$.
However, the added advantage of flexibility may be worth this extra cost.

\section{Numerical experiments}
We compare the proposed F-AB-GMRES preconditioned by Kaczmarz-type methods as inner iterations with AB-GMRES preconditioned by NE-SOR inner iterations~\cite{MH13, MH15} in terms of the central processing unit (CPU) time by numerical experiments on underdetermined and overdetermined problems.

Table~\ref{tab:1} gives the number of rows $m$, the number of columns $n$, the density of the nonzero elements, the rank and the condition number $\kappa(A)$ on the overdetermined test matrices.
The matrices RANDL$i$, $i=1,2,\ldots,6$ were randomly generated using the MATLAB function \texttt{sprandn}, as in~\cite{HYI10,MH13,MH15}. The Maragal$\_j$, $j= 3,4,5$ are rank-deficient matrices from~\cite{DH11}.
These matrices were transposed to form underdetermined problems. Table~\ref{tab:1} shows the effective size of the matrices after removing all zero rows. (If the matrix $A$ has a zero row, then the Kaczmarz-type methods can not work.) The rank and condition number were computed using the MATLAB functions \texttt{sprank}~\cite{FL} and \texttt{svd}.

In our implementations, a solution vector $x_{\star} \in {\mathbb{R}^{n}} $ is randomly generated by using the MATLAB function \texttt{randn}, and the right-hand side $b \in {\mathbb{R}^{m}}$ is taken to be $Ax_{\star}$. All computations were started from the initial vector $x_0=0$, and the stopping criterion is the relative residual
\begin{equation*}
\frac{{{{\left\| {{{r}_k}} \right\|}_2}}}{{{{\left\| b \right\|}_2}}} \le {10^{ - 6}},
\end{equation*}
where $r_k=b-Ax_k$ is the $k$th residual vector, or the maximum number of outer iterations 2000.
No restarts were used for GMRES.

All experiments were carried out using MATLAB (version R2017b) on a personal computer
with 2.80 GHz CPU (Intel(R) Core(TM) i7-7700HQ Q9550), 16
GB memory, and Windows 10 Home version 1909.

\subsection{Automatic parameter tuning for Kaczmarz-type methods}
The proposed method requires two preconditioning parameters: the maximum number of inner iterations $\ell_\text{max}$ in line 4 of Algorithm~\ref{AGPK} and the relaxation parameter $\omega$ used for the Kaczmarz-type inner iterations. Since the CPU time for the proposed method varies with the values of these parameters, it is desirable to determine the values automatically for any problem. Inspired by the idea in~\cite{MH13,MH15}, we perform the following procedure given as Algorithm~\ref{PTPK} using Kaczmarz-type methods alone for $Az=b$ before starting the outer iterations to determine the values of these parameters $\ell_\text{max}$ and $\omega_\text{opt}$.
	\begin{algorithm}[H]
			\caption{Parameter tuning procedure of Kaczmarz-type methods}
			\label{PTPK}
			\begin{algorithmic}[1]
				\STATE Set $\omega=1$ and $z^{(0)}=0$.
				\FOR{$\ell = 1, 2, \ldots $}
				\STATE Apply $\ell$ iterations of a Kaczmarz-type method to $Az=b$ to obtain $z^{(\ell)}$.
				\IF{${\left\| {b - A{z^{(\ell)}}} \right\|_2} \leq \eta {\left\| b \right\|_2}$}
				\STATE break
				\ENDIF
				\ENDFOR
				\STATE $\ell_\text{max} = \ell$
				\FOR{$\omega = 0.1, 0.2, \ldots, 1.9$}
				\STATE Apply $\ell_\text{max}$ iterations of a Kaczmarz-type method to $Az=b$.
				\ENDFOR
				\STATE ${\omega _\text{opt}} = \arg \mathop {\min }_{\omega=0.1,0.2,\ldots,1.9} {\left\| {{r^{(\ell_\text{max})}}} \right\|_2}$
			\end{algorithmic}
	\end{algorithm}
We remark that Morikuni and Hayami~\cite{MH13} evaluated the performance of AB-GMRES preconditioned by NE-SOR inner iterations for
different values of $\eta$ and finally chose $\eta=10^{-1}$.

Here, we also choose $\eta=10^{-1}$.

Since the number of inner iterations in F-AB-GMRES does not have to be fixed for each outer iteration, we set $\ell_\text{max}$ determined by Algorithm~\ref{PTPK} to be the maximum number of inner iterations for each outer iteration.

\begin{table}[!t]
	\caption{Information of the matrices.}\label{tab:1}
	\begin{center}
		\begin{tabular}{rrrrrr}
			\hline
			\multicolumn{1}{c}{Name} & \multicolumn{1}{c}{$m$} & \multicolumn{1}{c}{$n$} &\multicolumn{1}{c}{ Density[\%]} & \multicolumn{1}{c}{Rank} & \multicolumn{1}{c}{$\kappa(A)$} \\
			\hline
			RANDL1 & 5000 & 500 & 20 & 500 & 1.00$\times 10^{1}$ \\
			RANDL2 & 5000 & 500 & 20 & 500 & 1.00$\times 10^{2}$ \\
			RANDL3 & 5000 & 500 & 20 & 500 & 1.00$\times 10^{3}$ \\
			RANDL4 & 5000 & 500 & 20 & 500 & 1.00$\times 10^{4}$ \\
			RANDL5 & 5000 & 500 & 20 & 500 & 1.00$\times 10^{5}$ \\
			RANDL6 & 5000 & 500 & 20 & 500 & 1.00$\times 10^{6}$  \\
			Maragal$\_$3 & 1682 & 858 & 1.27 & 613 & 1.10$\times 10^{3}$ \\
			Maragal$\_$4 & 1964 & 1027 & 1.32 & 801 & 9.33$\times 10^{6}$ \\
			Maragal$\_$5 & 4654& 3296 & 0.61 & 2147 & 1.19$\times 10^{5}$ \\
			\hline
		\end{tabular}
		\begin{tablenotes}
			Name: name of the matrix, $m$: number of rows of the matrix, $n$: number of columns of the matrix, Density: density of the nonzero components of the matrix, Rank: maximum number of linearly independent columns of the matrix, $\kappa(A)$: condition number of the matrix ${\sigma _{\max }}/{\sigma _{\min }}$, where ${\sigma _{\max }}$ and ${\sigma _{\min }}$ are the largest and smallest nonzero singular values of the matrix, respectively.
		\end{tablenotes}
	\end{center}
\end{table}

\subsection{Underdetermined problems}
We first present numerical experiment results on underdetermined problems ($m < n$).
Tables~\ref{tab:2} and~\ref{tab:3} give the numerical experiment results on full-rank matrices and rank-deficient matrices, respectively.
The letter T at the end of the name of a matrix denotes the transposition of the matrix.
The first row in each cell in Tables~\ref{tab:2} and~\ref{tab:3} gives the number of outer iterations outside parentheses and gives the total number of inner iterations (the sum of the number of inner iterations in each outer iteration)

and relaxation parameter in parentheses.
Here, the number of inner iterations means $\ell m$ in Algorithm~\ref{NE-SOR}, and $\ell_{k}$ in line 4 of Algorithm~\ref{AGPK}.
The second row in each cell gives the total CPU time in seconds including the parameter tuning time and the formation of $C$ outside parentheses, and the parameter tuning time in seconds in parentheses. Here, the CPU time means the median of the elapsed CPU time for 10 times of repeated runs for the same $b$

for the corresponding method, respectively. The $\star$ indicates the fastest method regarding the CPU time.
In order to improve computing efficiency for the GK and GRK methods, we used the recursive formula~\eqref{r} to update the residual vector $s^{(p)}$ for each inner iteration.
To do so, we need to compute $C=AA^\mathsf{T}$ in advance.
The CPU time in seconds for computing matrix $C$ is given below the name of the matrix.
The total CPU time for F-AB-GMRES preconditioned by GK and GRK inner iterations includes the time for computing $C$. We also remark that the matrix $A^\mathsf{T}$ instead of $A$ was used throughout the programs for efficient data access with MATLAB. (The CPU time required to transpose $A$ is negligible.)

\begin{table}[!t]
    {\footnotesize
	\setlength{\tabcolsep}{4.8pt}
	\caption{Results for full-rank matrices (underdetermined
		problems)}\label{tab:2}
	\begin{center}
		\begin{tabular}{ll|rrrrrr}
			\hline
			Outer & Inner & \multicolumn{2}{c}{RANDL1T} &
			\multicolumn{2}{c}{RANDL2T} & \multicolumn{2}{c}{RANDL3T} \\
			iteration & iteration & \multicolumn{2}{c}{0.16} &
			\multicolumn{2}{c}{0.15} & \multicolumn{2}{c}{0.15} \\
			\hline
			AB-GMRES & NE-SOR &  11 & (16500, 1.1) & 29 & (101500,
			1.0) & 107 & (267500, 1.0) \\
			& & 1.36 & (0.76) & 4.38 & (1.70) & 7.87 & (1.21) \\
			\hline
			F-AB-GMRES & RK & 20.5 & (34450, 1.2) & 47.5 & (265240, 1.3) & 232.5 & (1074371,
			0.9) \\
			& & 4.07 & (2.41) & 17.21 & (5.67) & 49.41 & (3.99) \\
			\cline{2-8}
			& GRK & 20 & (5371.5, 1.1) & 38.5 & (29305, 1.3) & 122.5 & (83258, 1.3) \\
			& & 1.01 & (0.54) & 2.82 & (1.04) & 5.57 & (0.76) \\
			\cline{2-8}
			& GK & 29 & (6553, 1.1) & 38 & (31080, 1.2) & 119 & (82532, 1.1) \\
			& & \color{red}$\star$0.65 & (0.29) & \color{red}$\star$1.60 &
			(0.57) & \color{red}$\star$2.91 & (0.41) \\
			\hline
		\end{tabular}
		\end{center}}
		\medskip
		{\footnotesize
	    \setlength{\tabcolsep}{4.6pt}
	    \begin{center}
		\begin{tabular}{ll|rrrrrr}
			\hline
			Outer & Inner & \multicolumn{2}{c}{RANDL4T} &
			\multicolumn{2}{c}{RANDL5T} & \multicolumn{2}{c}{RANDL6T} \\
			iteration & iteration & \multicolumn{2}{c}{0.17} &
			\multicolumn{2}{c}{0.16} & \multicolumn{2}{c}{0.16} \\
			\hline
			AB-GMRES & NE-SOR & 72 & (180000, 0.9)& 199 & (497500,
			0.9) & 132 & (330000, 1.1) \\
			& & 5.85 & (1.22) & 13.67 & (1.28) & 9.47 & (1.22) \\
			\hline
			F-AB-GMRES & RK & 163 & (672860, 1.3) & 345.5 & (1560800, 1.0) & 243 &
			(987859, 1.0) \\
			& & 33.14 & (3.68) & 72.7 & (4.03) & 46.69 & (3.55) \\
			\cline{2-8}
			& GRK & 122 & (79897, 1.2) & 327.5 & (233590, 1.1) & 157 & (116580, 1.2) \\
			& & 5.29 & (0.73) & 13.23 & (0.76) & 7.35 & (0.83) \\
			\cline{2-8}
			& GK & 86 & (53544, 1.2) & 264 & (189135, 1.2) & 149 & (118184, 1.3) \\
			& & \color{red}$\star$2.07 & (0.36) & \color{red}$\star$5.78 &
			(0.40) & \color{red}$\star$3.89 & (0.45) \\
			\hline
		\end{tabular}
		\end{center}}
		{\footnotesize
	    \begin{center}
		\begin{tablenotes}
			\item[1] First row: Number of (outer) iterations (the total number of inner iterations, $\omega$).
			\item[2] Second row: Total CPU time, which includes parameter tuning
			time in parentheses, in seconds.
		\end{tablenotes}
	    \end{center}}
\end{table}

Table~\ref{tab:2} shows that F-AB-GMRES preconditioned by the GK inner iterations is the fastest among all the methods for the full-rank matrices RANDL$i$, $i=1,2,\ldots,6$.
We remark that the total number of inner iterations of F-AB-GMRES with GRK and GK is also smaller than that of AB-GMRES with NE-SOR. This may imply that F-AB-GMRES with GRK and GK has a smaller workload than AB-GMRES with NE-SOR.

In Figure~\ref{fig:1}, we plot the relative residual norm ${{{{\left\| {{{r}_k}} \right\|}_2}}}/{{{{\left\| b \right\|}_2}}}$ versus total number of inner iterations and CPU time in seconds for the matrix RANDL6T.
The figure shows that F-AB-GMRES preconditioned by the GRK inner iterations is best among all the methods regarding the total number of inner iterations, and the GK inner iterations is the fastest among all the methods regarding CPU time. These results are in accordance with Table~\ref{tab:2}.

\begin{figure}[!t]
    \centering
	\begin{minipage}[t]{0.46\linewidth}
		\centering
		\includegraphics[width=2.3in]{R6T_IT.eps}
	\end{minipage}
	\begin{minipage}[t]{0.46\linewidth}
		\centering
		\includegraphics[width=2.3in]{R6T_CPU.eps}
	\end{minipage}
	\caption{Relative residual norm vs.\ total number of inner iterations (left) and relative residual norm vs.\ CPU time (right) for RANDL6T.}
	\label{fig:1}
\end{figure}

\begin{table}[!t]
    {\footnotesize
	\setlength{\tabcolsep}{3.8pt}
	\caption{Results for rank-deficient matrices (underdetermined problems)}\label{tab:3}
	\begin{center}
		\begin{tabular}{ll|rrrrrr}
			\hline
			Outer & Inner & \multicolumn{2}{c}{Maragal$\_$3T } &
			\multicolumn{2}{c}{Maragal$\_$4T } & \multicolumn{2}{c}{Maragal$\_$5T } \\
			iteration & iteration & \multicolumn{2}{c}{0.01} &
			\multicolumn{2}{c}{0.02} & \multicolumn{2}{c}{0.21} \\
			\hline
			AB-GMRES & NE-SOR &  57 & (195624, 1.2) & 51 & (209508,
			1.1) & 144 & (1898496, 1.2) \\
			& & 3.05 & (0.79) & 3.69 & (1.08) & 44.7 & (5.05) \\
			\hline
			F-AB-GMRES & RK & 305 & (1352502, 1.0) & 179 & (1161812, 1.1) & 514 & (8856928,
			0.9) \\
			& & 48.59 & (3.11) & 46.71 & (4.79) & 730.07 & (27.05) \\
			\cline{2-8}
			& GRK & 104.5 & (64322, 1.2) & 132 & (124250, 1.2) & 250.5 & (578360, 1.1) \\
			& & 3.50 & (0.62) & 7.05 & (1.01) & 73.57 & (6.05) \\
			\cline{2-8}
			& GK & 121 & (73165, 1.0) & 99 & (84543, 1.1) & 228 & (555285, 1.1) \\
			& & \color{red}$\star$1.32 & (0.20) & \color{red}$\star$1.67 &
			(0.31) & \color{red}$\star$19.13 & (1.44) \\
			\hline
		\end{tabular}
		\end{center}}
		{\footnotesize
	    \begin{center}
		\begin{tablenotes}
		\item[1] First row: Number of (outer) iterations (the total number of inner iterations, $\omega$).
		\item[2] Second row: Total CPU time, which includes parameter tuning time in parentheses, in seconds.
		\end{tablenotes}
	\end{center}}
\end{table}

Table~\ref{tab:3} shows that F-AB-GMRES preconditioned by the GK inner iterations is also the fastest among all the methods for the rank-deficient matrices Maragal$\_$3T, Maragal$\_$4T and Maragal$\_$5T.

\begin{figure}[!t]
    \centering
	\begin{minipage}[t]{0.46\linewidth}
		\centering
		\includegraphics[width=2.3in]{M5T_IT.eps}
	\end{minipage}
	\begin{minipage}[t]{0.46\linewidth}
		\centering
		\includegraphics[width=2.3in]{M5T_CPU.eps}
	\end{minipage}
	\caption{Relative residual norm vs.\ total number of inner iterations (left) and relative residual norm vs.\ CPU time (right) for Maragal$\_$5T.}
	\label{fig:2}
\end{figure}

In Figure~\ref{fig:2}, we plot the relative residual norm versus total number of inner iterations and CPU time in seconds for the matrix Maragal$\_$5T.
Figure~\ref{fig:2} shows that F-AB-GMRES preconditioned by the GRK inner iterations is best among all the methods regarding the total number of inner iterations, and the GK inner iterations is the fastest among all the methods regarding CPU time for the matrix Maragal$\_$5T. These results are in accordance with Table~\ref{tab:3}.

\subsection{Overdetermined problems}
Next, we present numerical experiment results on overdetermined problems ($m>n$). Table~\ref{tab:4} and~\ref{tab:5} give the numerical experiment results for full-rank matrices and rank-deficient matrices, respectively similarly to Tables~\ref{tab:2} and~\ref{tab:3}.

\begin{table}[!t]
    {\footnotesize
	\setlength{\tabcolsep}{4.0pt}
	\caption{Results for full-rank matrices (overdetermined problems)}\label{tab:4}
	\begin{center}
			\begin{tabular}{ll|rrrrrr}
				\hline
				Outer & Inner & \multicolumn{2}{c}{RANDL1} &
				\multicolumn{2}{c}{RANDL2} & \multicolumn{2}{c}{RANDL3} \\
				iteration & iteration & \multicolumn{2}{c}{2.26} &
				\multicolumn{2}{c}{2.32} & \multicolumn{2}{c}{2.24} \\
				\hline
				AB-GMRES & NE-SOR &  2 & (10000, 1.0) & 2 & (10000, 0.9) & 2 & (10000, 0.7) \\
				& & \color{red}$\star$6.10 & (3.32) & \color{red}$\star$6.27 & (3.43) & 6.10 & (3.34) \\
				\hline
				F-AB-GMRES & RK & 20 & (42352, 0.9) & 146 & (703920, 0.9) & 328.5 & (1298700,
				0.7) \\
				& & 16.25 & (8.87) & 103.41 & (13.23) & 174.79 & (9.59) \\
				\cline{2-8}
				& GRK & 60 & (14281, 1.0) & 79.5 & (25718, 1.0) & 53 & (11622, 1.0) \\
				& & 6.24 & (1.26) & 8.83 & (1.54) & \color{red}$\star$5.81 & (1.24) \\
				\cline{2-8}
				& GK & 220 & (71784, 1.0) & 383 & (178359, 1.0) & 469 & (140146, 1.0) \\
				& & 6.81 & (0.35) & 13.97 &
				(0.95) & 12.56 & (0.28) \\
				\hline
		\end{tabular}
		\end{center}}
		\medskip
		{\footnotesize
	    \setlength{\tabcolsep}{3.2pt}
    	\begin{center}
			\begin{tabular}{ll|rrrrrr}
				\hline
				Outer & Inner & \multicolumn{2}{c}{RANDL4} &
				\multicolumn{2}{c}{RANDL5} & \multicolumn{2}{c}{RANDL6} \\
				iteration & iteration & \multicolumn{2}{c}{2.30} &
				\multicolumn{2}{c}{2.35} & \multicolumn{2}{c}{2.34} \\
				\hline
				AB-GMRES & NE-SOR & 2 & (10000, 1.0)& 2 & (10000,
				1.0) & 2 & (10000, 1.0) \\
				& & \color{red}$\star$6.27 & (3.44) & 6.24 & (3.37) & 6.26 & (3.40) \\
				\hline
				F-AB-GMRES & RK & 278 & (1295432, 0.7) & 330 & (1544889, 0.8) & 259.5 &
				(1524497, 1.2) \\
				& & 176.5 & (11.43) & 207.13 & (11.27) & 204.36 & (14.46) \\
				\cline{2-8}
				& GRK & 74.5 & (17678, 1.0) & 80.5 & (21406, 1.0) & 53.5 & (11422, 1.0) \\
				& & 7.04 & (1.20) & 7.9 & (1.29) & 5.94 & (1.27) \\
				\cline{2-8}
				& GK & 241 & (114389, 1.0) & 129 & (43046, 1.0) & 136 & (35842, 1.0) \\
				& & 9.21 & (0.32) & \color{red}$\star$5.23 &
				(0.33) & \color{red}$\star$4.81 & (0.29) \\
				\hline
		\end{tabular}
		\end{center}}
		{\footnotesize
	    \begin{center}
		\begin{tablenotes}
			\item[1] First row: Number of (outer) iterations (the total number of inner iterations, $\omega$).
			\item[2] Second row: Total CPU time, which includes parameter tuning time in parentheses, in seconds.
		\end{tablenotes}
	\end{center}}
\end{table}

Table~\ref{tab:4} shows that AB-GMRES preconditioned by the NE-SOR inner iterations is the fastest among all the methods for matrices RANDL1, RANDL2 and RANDL4. F-AB-GMRES preconditioned by the GRK inner iterations is the fastest among all the methods for matrices RANDL3, and the GK inner iterations is the fastest among all the methods for matrices RANDL5 and RANDL6.

\begin{table}[!t]
    {\footnotesize
	\setlength{\tabcolsep}{3.5pt}
	\caption{Results for rank-deficient matrices (overdetermined problems)}\label{tab:5}
	\begin{center}
		\begin{tabular}{ll|rrrrrr}
			\hline
			Outer & Inner & \multicolumn{2}{c}{Maragal$\_$3 } &
			\multicolumn{2}{c}{Maragal$\_$4 } & \multicolumn{2}{c}{Maragal$\_$5 } \\
			iteration & iteration & \multicolumn{2}{c}{0.02} &
			\multicolumn{2}{c}{0.02} & \multicolumn{2}{c}{0.07} \\
			\hline
			AB-GMRES & NE-SOR &  149 & (751854, 1.1) & 97 & (571524,
			1.1) & 330 & (4607460, 1.1) \\
			& & 9.35 & (1.16) & 7.67 & (1.33) & 78.40 & (4.21) \\
			\hline
			F-AB-GMRES & RK & 238 & (1303000, 0.9) & 165.5 & (896606, 1.0) & 468 & (8781000,
			0.9) \\
			& & 56.60 & (4.71) & 42.60 & (5.08) & 906.15 & (37.92) \\
			\cline{2-8}
			& GRK & 262.5 & (146302, 1.0) & 173.5 & (124610, 1.1) & 457 & (933609, 1.1) \\
			& & 8.95 & (0.71) & 8.42 & (0.92) & 152.79 & (6.38) \\
			\cline{2-8}
			& GK & 253 & (156155, 1.0) & 200 & (161556, 1.1) & 520 & (1192959, 1.1) \\
			& & \color{red}$\star$2.98 & (0.24) & \color{red}$\star$3.08 &
			(0.29) & \color{red}$\star$45.46 & (1.38) \\
			\hline
		\end{tabular}
		\end{center}}
		{\footnotesize
	    \begin{center}
		\begin{tablenotes}
			\item[1] First row: Number of (outer) iterations (the total number of inner iterations, $\omega$).
			\item[2] Second row: Total CPU time, which includes parameter tuning time in parentheses, in seconds.
		\end{tablenotes}
	    \end{center}}
\end{table}

Table~\ref{tab:5} shows that F-AB-GMRES preconditioned by
the GK method is the fastest regarding the CPU time among all the methods for matrices Maragal$\_$3, Maragal$\_$4 and Maragal$\_$5.

In Figure~\ref{fig:4}, we plot the relative residual norm versus the total number of inner iterations and CPU time in seconds for the matrix Maragal$\_$5. Figure~\ref{fig:4} shows that F-AB-GMRES preconditioned by the GRK inner iterations is best among all the methods when comparing the total number of inner iterations, and the GK inner iterations is the fastest among all the methods regarding the CPU time for the matrix Maragal$\_$5. These results are in accordance with Table~\ref{tab:5}.

\begin{figure}[!t]
    \centering
	\begin{minipage}[t]{0.46\linewidth}
		\centering
		\includegraphics[width=2.3in]{M5_IT.eps}
	\end{minipage}
	\begin{minipage}[t]{0.46\linewidth}
		\centering
		\includegraphics[width=2.3in]{M5_CPU.eps}
	\end{minipage}
	\caption{Relative residual norm vs.\ total number of inner iterations (left) and relative residual norm vs.\ CPU time (right) for Maragal$\_$5.}
	\label{fig:4}
\end{figure}

We have tried to further speed up the methods based on GK and GRK by computing approximations of $AA^\mathsf{T}$ using the method in~\cite{HI15} for over- and underdetermined systems, but so far we have not been successful, and this is left for future research.

\section{Conclusion}
In this paper, we proposed replacing the NE-SOR method by Kaczmarz-type methods in the previous AB-GMRES method preconditioned by stationary inner iterations for solving consistent systems of linear equations.
To do so, we developed a new algorithm called flexible AB-GMRES method preconditioned by Kaczmarz-type methods as inner iterations.
An optimality property of minimizing residuals was given for the proposed method.
We also proposed a tuning procedure for adjusting the maximum number of inner iterations and value of the relaxation parameter in the method.
Numerical experiment results showed that flexible AB-GMRES preconditioned by Kaczmarz-type methods converge faster than the previous method in terms of total CPU time.

\section*{Appendix} \textbf{Comparison of computational work of GK and modified GK inner-iteration preconditioning.}

We compare the total computational work (the number of floating-point operations) for the greedy Kaczmarz inner-iteration preconditioning using Algorithm~\ref{The Greedy Kaczmarz Method} (GK), with GK modified by precomputing and storing $C=AA^\mathsf{T}$ once beforehand and updating the residual vector $s$ using equation~\eqref{r} (MGK).

Let $A$ be an $m \times n$ matrix. Assume that the number of outer iterations of the F-AB-GMRES is $k$, and that the number of inner Kaczmarz iterations is fixed at $\ell$ for each outer iteration.

First, consider the case when $A$ is dense. Then, the total work for GK is given by
\begin{equation*}
w_{\text{GK}}^\text{d} = k \ell (mn+n).
\end{equation*}
The first term corresponds to step 3, and the second corresponds to step 4 of Algorithm~\ref{The Greedy Kaczmarz Method}, respectively.
The total work for MGK is
\begin{equation*}
w_{\text{MGK}}^\text{d} = m^2 n + k \ell (m+n).
\end{equation*}
The first term corresponds to computing  $C=AA^\mathsf{T}$ once beforehand, the second is for step 4 of Algorithm~\ref{The Greedy Kaczmarz Method}, and the third corresponds to the update in equation~\eqref{r}. Hence,
\begin{equation*}
w_{\text{GK}}^\text{d} - w_{\text{MGK}}^\text{d} = m \left[ k \ell (n-1) - mn \right].
\end{equation*}
Therefore,
\begin{equation}\label{wd}
w_{\text{MGK}}^\text{d}  < w_{\text{GK}}^\text{d} \Longleftrightarrow k \ell > m \left( 1 + \frac{1}{n-1} \right).
\end{equation}

Next, consider the case when $A$ is sparse and the position of the nonzero elements are random. Let $\mathrm{nz}$ be the number of nonzero elements of $A$. Define ${\displaystyle q= \mathrm{nz}/m}$ as the average number of nonzero elements per row of $A$. Thus, the density of $A$ is $d={\displaystyle q/n}$. Assume that the computational work to compute $C = (c_{ij}) = AA^\mathsf{T}$ is $qm^2$.
Let the density of $C$ be $p$.

Then, the total work for GK is
\begin{equation*}
w_{\text{GK}}^\text{s} = k \ell (\mathrm{nz} + q) = k \ell q (m+1).
\end{equation*}
The first term corresponds to step 3, and the second corresponds to step 4 of Algorithm~\ref{The Greedy Kaczmarz Method}, respectively.
The total work for MGK is
\begin{equation*}
w_{\text{MGK}}^\text{s} = qm^2+ k \ell \left(q+mp\right).
\end{equation*}
The first term is for computing $AA^\mathsf{T}$, the second is for step 4 of Algorithm~\ref{The Greedy Kaczmarz Method}, and the third is for the update in equation~\eqref{r}, respectively.
Hence,
\begin{equation*}
 w_{\text{GK}}^\text{s} - w_{\text{MGK}}^\text{s} = m\left[k \ell \left(q-p\right)  -mq\right].
\end{equation*}
Therefore,
\begin{equation}\label{ws}
w_{\text{MGK}}^\text{s}  <  w_{\text{GK}}^\text{s} \Longleftrightarrow k \ell > m\left(1 + \frac{p}{q-p}\right).
\end{equation}

The density $p$ of $C=AA^\mathsf{T}$ can be estimated as follows. The probability that $c_{ij} \ne 0$ for $i \ne j$ is ${p_{\text{nd}}} = 1 - {\left(1 - {d^2}\right)^n}$, and the probability that $c_{ii} \ne 0$ is ${p_\text{d}} = 1 - {(1 - d)^n}$. Therefore the probability that $c_{ij} \ne 0$ (or the density of $C$) is
\begin{align} \label{p}
    p &= \frac{{\left( {{m^2} - m} \right){p_{\text{nd}}} + m{p_\text{d}}}}{{{m^2}}} \nonumber\\
    &= 1 - \left( {1 + \frac{1}{m}} \right){\left( {1 - {d^2}} \right)^n} + \frac{1}{m}{\left( {1 - d} \right)^n}.
\end{align}
If $d=1$ ($A$ is dense),~\eqref{p} implies $p=1$. Then, also $q=n$, so that~\eqref{ws} agrees with~\eqref{wd}.
If $d \ll 1$, $p$ can be approximated as
\begin{equation*}
    p \sim 1 - \left( {1 + \frac{1}{m}} \right){e^{ - n{d^2}}} + \frac{1}{m}{e^{ - q}}.
\end{equation*}

Table~\ref{tab:6} gives estimated (using~\eqref{p}) and actual values of the density of $C=AA^\mathsf{T}$ for the matrices used in our experiments.

\begin{table}[!t]
	\caption{The estimated and actual densities of the matrix $C=AA^\mathsf{T}$.}\label{tab:6}
	\begin{center}
		\begin{tabular}{rrrrrr}
			\hline
			\multicolumn{1}{c}{matrix} & \multicolumn{1}{c}{$m$} & \multicolumn{1}{c}{$n$} &\multicolumn{1}{c}{$d$} & \multicolumn{1}{c}{$p(\text{estimated})$} & \multicolumn{1}{c}{$p(\text{actual})$} \\
			\hline
			RANDL1 & 5000 & 500 & 0.2 & 1.00 & 0.745 \\
			RANDL2 & 5000 & 500 & 0.2 & 1.00 & 0.724\\
			RANDL3 & 5000 & 500 & 0.2 & 1.00 & 0.731\\
			RANDL4 & 5000 & 500 & 0.2 & 1.00 & 0.714\\
			RANDL5 & 5000 & 500 & 0.2 & 1.00 & 0.719\\
			RANDL6 & 5000 & 500 & 0.2 & 1.00 & 0.710\\
			Maragal$\_$3 & 1682 & 858 & 1.27$\times 10^{-2}$ & 0.129 & 0.160\\
			Maragal$\_$4 & 1964 & 1027 & 1.32$\times 10^{-2}$ & 0.163 & 0.126\\
			Maragal$\_$5 & 4654& 3296 & 6.10$\times 10^{-3}$ & 0.115 &  0.0728\\
			RANDL1T  & 500 & 5000 & 0.2 & 1.00 & 0.940 \\
			RANDL2T  & 500 & 5000 & 0.2 & 1.00 & 0.940\\
			RANDL3T  & 500 & 5000 & 0.2 & 1.00 & 0.926\\
			RANDL4T  & 500 & 5000 & 0.2 & 1.00 & 0.927\\
			RANDL5T  & 500 & 5000 & 0.2 & 1.00 & 0.932\\
			RANDL6T  & 500 & 5000 & 0.2 & 1.00 & 0.940\\
			Maragal$\_$3T  & 858 & 1682 & 1.27$\times 10^{-2}$ & 0.237 & 0.562\\
			Maragal$\_$4T  & 1027 & 1964 & 1.32$\times 10^{-2}$ & 0.289 & 0.669\\
			Maragal$\_$5T & 3296 & 4654 & 6.10$\times 10^{-3}$ & 0.159 & 0.461\\
			\hline
		\end{tabular}
		\begin{tablenotes}
		    $m$: number of rows of $A$, $n$: number of columns of $A$,
		    $d$: density of nonzero elements of $A$, $p$: density of nonzero elements of $C=AA^\mathsf{T}$.
		\end{tablenotes}
	\end{center}
\end{table}

As for the CPU time, the computation of $AA^\mathsf{T}$ should perform relatively more efficiently than the flops count suggests, especially for the dense case, due to fast memory access.

\section*{Acknowledgement}
We would like to thank Professor Zhong-Zhi Bai for stimulating discussions and valuable remarks.

\bibliographystyle{siamplain}
\bibliography{references}

\end{document}

config: arxiv


TITLE: Seating Arrangement and Independent Events?
QUESTION [0 upvotes]: PROBLEM: Alice, Bob, Catherine, Doug, and Edna are randomly assigned seats at a circular table in a perfectly circular room.  Assume that rotations of the table do not matter, so there are exactly 24 possible outcomes in the sample space.
Bob and Catherine are married.  Doug and Edna are married.  When people are married they love to sit beside each other.
Let T denote the event that Bob and Catherine are sitting next to each other.  Let U be the event that Alice and Bob are sitting next to each other. Are events T and U independent? Why?  Justify your answer.
My QUESTIONS:
1. Problem statement says couples love to sit beside each other. Why did author include this statement? Does he mean by "love to sit beside each other" that "they are always sitting beside each other" or he is just trying to confuse you ?
2. How to start solving the problem. I mean I can't even think of with what I should starts ;(
If by "love to" author meant they are "always sitting together" then nothing matters and $P(T) = 1$  and Alice has choice to sit next to rest 4 people: $P(U) = \frac14$
Since they always sit together and Alice got its own wish to sit anywhere, no one is forcing her, they are independent events. This is the reasoning I could come up with.

REPLY [1 votes]: $1)$ The author means nothing relevant when he says that. It's just an an inane piece of information meant to provide some idea of a context for the actual problem statement. It has nothing to do with solving the problem, i.e. even if the author had not put this line, the solution to the problem would not have changed one iota.
$2)$ I am going to give a mathematical proof here. There might be a logical way to do this too. To find whether the events are independent or not, we are going to use the fact that if two events are independent, then
$$P(A\cap B) = P(A)P(B)$$
Now to calculate $P(T)$ and $P(U)$.
Starting with $T$, we see that Bob and Catherine are sitting together. So we combine them in one unit and calculate number of ways of seating them around the table (BC),E,D,A (Here letters represent the initials of the names of people). Since there are 4 entities, there are $3!$ ways to do it. Also, we can shift Bob and Catherine around (as in (BC) and (CB)) in two ways giving us
$$|T| = 3!\cdot2!$$
$$P(T) = \frac{3!\cdot2!}{4!}=\frac12$$
Applying a similar logic for $U$, we see that Bob and Alice are sitting together. So we combine them in one unit and calculate number of ways of seating them around the table (BA),E,D,C (Here letters represent the initials of the names of people). Since there are 4 entities, there are $3!$ ways to do it. Also, we can shift Bob and Alice around (as in (BA) and (AB)) in two ways giving us
$$|U| = 3!\cdot2!$$
$$P(U) = \frac{3!\cdot2!}{4!}=\frac12$$
Now we calculate $P(T\cap U)$. For this we put $T$ and $U$ together, i.e. both Alice and Catherine are sitting together with Bob. So we combine the three of them in one unit and try to seat 3 entities around a table (CBA),D,E. This can be done in $2!$ ways. Also Catherine and Alice can switch places inside the entity giving us a multiplicative factor of $2$. Now the total number of ways is
$$|T\cap U| = 2!\cdot 2!$$
$$P(T\cap U) = \frac{2!\cdot 2!}{4!} = \frac16$$
We can see that $P(T\cap U) \ne P(T)P(U)$, meaning they are not independent.
I am attaching a few links regarding independent events in probability and circular permutations for your understanding
https://brilliant.org/wiki/probability-independent-events/
https://www.toppr.com/guides/business-mathematics-and-statistics/permutations-and-combinations/permutations-and-circular-permutation/


TITLE: Double partial derivative at maximum or minimum points
QUESTION [0 upvotes]: Is this true?
Let $U\subset \mathbb R^n$ be open set, and $u:U \to \mathbb R$ be differentiable every order . If $x_0\in U$ is a maximum point of $u$, then
$$u_{x_i}(x)=0$$ and $$u_{x_ix_i}(x)<0.$$
Thanks

REPLY [0 votes]: The fact that $U$ is an open sets makes the "gradient equal to zero at the maximum" a necessary condition (if U was a closed set, then it would not be a necessary condition). The second derivative (Hessian) being negative (negative definite) is just a sufficient condition, so the condition might not hold, as the examples in the other answers have showed.


TITLE: Find matrix $\mathscr{M}(T)$ of the following linear map.
QUESTION [1 upvotes]: Let $u_1,u_2,u_3,u_4$ be the basis of $V$ an let $T \in \mathscr{L}(v)$ be given by
$Tu_1 = u_1$
$Tu_i = u_i - u_{i-1} , i > 1$
Find $\mathscr{M}(T)$.
My response to this incorrect. I was presumed that this would be a diagonal matrix and I ended up writing
\begin{bmatrix}1&0&0&0\\0&0&0&0\\0&0&1&0\\0&0&0&0\end{bmatrix}
My reasoning was that 1 minus 0 was 1 and 1 minus 1 was 0...
The actual answer is
 \begin{bmatrix}1&-1&0&0\\0&1&-1&0\\0&0&1&-1\\0&0&0&1\end{bmatrix}
My understanding of this solution is that there is some computation that I can't seem recall. My guess is that row1column1 is just 1 given by the defined map. then the ith entry row1column2 1-1 = 0. Not sure how to get -1.

REPLY [2 votes]: In general, say $T$ is a linear transformation from $V$ to $V$ and $\mathcal B=\{v_1,\ldots,v_n\}$ is a basis of $V$. Then there are constants $c_{i,j}$ with $i,j=1,\ldots,n$ such that
$$Tv_1 = c_{1,1}v_1+c_{1,2}v_2+\ldots+c_{1,n}v_n$$ $$\vdots$$
$$Tv_n = c_{n,1}v_1+c_{n,2}v_2+\ldots+c_{n,n}v_n$$
And we consider $[T]_{\mathcal B}$, the matrix representation of $T$ under basis $\mathcal B$ as
$$[T]_{\mathcal B}=\begin{pmatrix}
c_{1,1}& c_{2,1}&\ldots&c_{n,1}\\
c_{1,2}&c_{2,2}&\ldots&c_{n,2}\\
\vdots &\vdots&\ddots&\vdots\\
c_{1,n}&c_{2,n}&\ldots&c_{n,n}
\end{pmatrix}$$

REPLY [2 votes]: We have $u_1=\begin{bmatrix}1\\0\\0\\0\end{bmatrix}$, $u_2=\begin{bmatrix}0\\1\\0\\0\end{bmatrix}$,
$u_3=\begin{bmatrix}0\\0\\1\\0\end{bmatrix}$,
$u_4=\begin{bmatrix}0\\0\\0\\1\end{bmatrix}$.  To apply our linear map to a vector, we multiply $\mathscr{M}(T)$ by that vector, written as a column on the right.
$\begin{bmatrix}1&0&0&0\\0&0&0&0\\0&0&1&0\\0&0&0&0\end{bmatrix} u_2=\begin{bmatrix}1&0&0&0\\0&0&0&0\\0&0&1&0\\0&0&0&0\end{bmatrix}\begin{bmatrix}0\\1\\0\\0\end{bmatrix}=\begin{bmatrix}0\\0\\0\\0\end{bmatrix}$
which is not what we want.
However $\begin{bmatrix}1&-1&0&0\\0&1&-1&0\\0&0&1&-1\\0&0&0&1\end{bmatrix} u_2=\begin{bmatrix}1&-1&0&0\\0&1&-1&0\\0&0&1&-1\\0&0&0&1\end{bmatrix} \begin{bmatrix}0\\1\\0\\0\end{bmatrix}=\begin{bmatrix}-1\\1\\0\\0\end{bmatrix}=-u_1+u_2$


TITLE: Measuring angle between E and B fields
QUESTION [2 upvotes]: How can we actually measure the angle between the $E$ and $B$ fields in a light beam? I know this angle is 90 degrees but I have hard time accepting it without experimental proof. Can one prove for example this is not 91 degrees?

REPLY [5 votes]: One can receive the electromagnetic plane wave with an antenna that is responsive to either the electric field or to the magnetic field but not to both within a reasonable engineering accuracy. For example, the standard electric dipole (thin wire linear antenna) or magnetic dipole (thin wire loop antenna) is essentially unresponsive to the magnetic or the electric field, respectively, when are properly aligned in space. It just turns out experimentally for a linearly polarized EM wave the maximum energy is received when the electric dipole is parallel with the electric field and the magnetic dipole is perpendicular to the electric field of the plane wave. Granted the accuracy of such orthogonality for practical antennas is not within 1 degree for simple dipoles but if you instead consider really high gain antennas whose beamwidth is fraction of a degree you do get to angular precision on the order of millidegrees.


TITLE: Calculation of Wilson Coefficients
QUESTION [2 upvotes]: In Flavour physics, amplitude of decay processes is generally expressed in terms of effective operators (reference).
In the framework of effective field theory amplitude $\propto C_iO_i$, where $O_i$ are the operators (currents) and $C_i$ are the wilson coefficients (numbers).
Now, effective field theory and full theory are matched at a particular scale to find values of Wilson coefficients. In this particular paper, where they are studying a $b\rightarrow s~ l^+l^-$ transition, they choose the matching scale to be $m_W$ (Mass of W boson).Why is that? Whay are they not choosing matching scale to be 246 GeV which is the electroweak scale?

REPLY [0 votes]: Because W boson is the lightest particle which is being integrated out. Below this scale Standard Model picture and effective picture vary but the physics remain same. So, the two pictures should match at the boundary.


TITLE: Proper usage of $\lim_{n\to\infty} a_n^{b_n} = \lim_{n\to\infty}a_n^{\lim_{n\to\infty}b_{n}}$
QUESTION [2 upvotes]: Let's consider the following limit
$$\lim_{n\to\infty} a_n^{b_n}. \tag{1}$$
Of course, in general, the following statement is false
$$\lim_{n\to\infty} a_n^{b_n} = \lim_{n\to\infty}a_n^{\lim_{n\to\infty}b_{n}},$$
we can give an obvious counterexample, i.e.
$$e = \lim_{n\to\infty} \bigg(1 + \frac{1}{n} \bigg)^{n} \neq \lim_{n\to\infty} \bigg(1 + \frac{1}{n} \bigg)^{\lim_{n\to\infty} n} = 1.$$
However there are plenty other examples which can be solved using $(1)$ and the result is correct. Is it enough to assume that $(1)$ holds if and only if $\lim_{n\to\infty} a_n = a$ and $\lim_{n\to\infty} b_n = b$, where $a, b \neq 0$ and $a, b$ are finite?
Let's consider one more example
$$\lim_{n\to\infty} \bigg(\frac{n+4}{n+3}\bigg)^{5-2n}.$$
After several simple transformations we get
$$\lim_{n\to\infty} \bigg(\frac{n+4}{n+3}\bigg)^{5-2n} = \lim_{n\to\infty} \Bigg(\bigg(1 + \frac{1}{n+3}\bigg)^{n+3}\Bigg)^{\frac{5-2n}{n+3}}.$$
Let
$$a_n = \bigg(1 + \frac{1}{n+3}\bigg)^{n+3}, \qquad b_n = \frac{5-2n}{n+3}.$$
Of course
$$a_n \to e, \qquad b_n \to -2,$$
due to the fact that both limits differ from zero and are finite we conclude that
$$\lim_{n\to\infty} \bigg(\frac{n+4}{n+3}\bigg)^{5-2n} = e^{-2}.$$
Is this reasoning correct?

REPLY [2 votes]: To elaborate on what Fishbane pointed out, we have that since $\log$ is a continuous function on $(0, \infty)$, then $\lim_{n \rightarrow \infty} \log(x_n) = \log(\lim_{n \rightarrow \infty} x_n)$ if $\lim_{n \rightarrow \infty} x_n > 0$. Then if $a_n \rightarrow a$ where $a \neq 0$, then
$$
\begin{align*}
\lim_{n \rightarrow \infty} \log(a_n^{b_n}) &= \lim_{n \rightarrow \infty} b_n \log(a_n) \\
&= \lim_{n \rightarrow \infty} b_n \ \cdot \lim_{n \rightarrow \infty} \log(a_n) \\
&= b \ \cdot \log(a).
\end{align*}
$$
Then since $e^x$ is continuous on all of $\mathbb{R}$, then
$$
\begin{align*}
\lim_{n \rightarrow \infty} a_n^{b_n} &= e^{b \log a} \\
&= a^b.
\end{align*}
$$
Note that throughout we never needed that $b \neq 0$; from the argument here we see that if $b = 0$, then $\lim_{n \rightarrow \infty} a_n^{b_n} = 1$. Hope this helps! Let me know if anything is unclear.


TITLE: Evaluating the integral $\int^{\infty}_{-\infty} \frac{dx}{x^4-2\cos(2\theta)x^2 +1}$
QUESTION [0 upvotes]: The first part of this question required me to find out the zeroes of the denominator, and to treat the equation as that of a complex number, which allows us to write:
$$\frac{1}{z^4-2\cos(2\theta)z^2 +1}=\frac{1}{(z-e^{i\theta})(z+e^{i\theta})(z-e^{-i\theta})(z+e^{-i\theta})}$$
As far as I can understand these zeros have the effect of giving a circle of singularities in the complex plane with radius 1.
I assume to do the integral given in the question, I would have to do some kind of contour integration over the complex plane, but I am stumped as to which contour I should use to do this. I assume a branch cut will also be needed due to the transcendental nature of the zeroes in the denominator.
As always and help is appreciated - thanks!

REPLY [0 votes]: $\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,}
 \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack}
 \newcommand{\dd}{\mathrm{d}}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\mc}[1]{\mathcal{#1}}
 \newcommand{\mrm}[1]{\mathrm{#1}}
 \newcommand{\pars}[1]{\left(\,{#1}\,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,}
 \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$
$\ds{\underline{\underline{\mbox{With}\ \theta \not\in \pi\mathbb{Z}}}}$:
\begin{align}
&\bbox[10px,#ffe]{\int_{-\infty}^{\infty}
{\dd x \over x^{4} - 2\cos\pars{2\theta}x^{2} + 1}} =
2\int_{0}^{\infty}
{x^{-2}\,\dd x \over x^{2} - 2\cos\pars{2\theta} + x^{-2}}
\\[5mm] = &\
2\int_{0}^{\infty}
{x^{-2}\,\dd x \over \pars{x - x^{-1}}^{2} +
2\bracks{1 -\cos\pars{2\theta}}}
\\[5mm] = &\
\int_{0}^{\infty}
{x^{-2}\,\dd x \over \pars{x - x^{-1}}^{2} + 4\sin^{2}\pars{\theta}} +
\int_{\infty}^{0}
{-\,\dd x \over \pars{x^{-1} - x}^{2} + 4\sin^{2}\pars{\theta}}
\\[5mm] = &\
\int_{0}^{\infty}
{\pars{1 + x^{-2}}\,\dd x \over \pars{x - x^{-1}}^{2} + 4\sin^{2}\pars{\theta}} =
\int_{-\infty}^{\infty}
{\dd x \over x^{2} + 4\sin^{2}\pars{\theta}}
\\[5mm] = &\
{1 \over 4\sin^{2}\pars{\theta}}\,
2\verts{\sin\pars{\theta}}\int_{-\infty}^{\infty}
{\dd x/\bracks{2\verts{\sin\pars{\theta}}} \over \braces{x/\bracks{2\verts{\sin\pars{\theta}}}}^{2} + 1}
\\[5mm] = &
{1 \over 2\verts{\sin\pars{\theta}}}\int_{-\infty}^{\infty}
{\dd x \over x^{2} + 1} =
\bbx{\pi \over 2\verts{\sin\pars{\theta}}}\,,\qquad
\theta \not\in \pi\mathbb{Z}
\end{align}


TITLE: Reducing an inquality with fractions
QUESTION [0 upvotes]: can you help me reduce the following inequality (i need to get a relation between x and y -- express x in terms of y)
$\frac{n}{2x} < \frac{n}{(4+\epsilon)y}+1$
I would like to show somehow that $x > (1+\epsilon) y$   OR that $x>2y$
The assumption is that $0 < \epsilon < 1$ and that $n>4x, n>4y$, x and y are positive integers.

REPLY [0 votes]: You already have what you want. For $x>0$, $y>0$, $\epsilon>0$, and $n>0$, you have: $$\dfrac{(4+\epsilon)}{n}\dfrac{n-2x}{2x}<y.$$
Hence:
$$y>\dfrac{a}{x}+b.$$
Where $a=\dfrac{4+\epsilon}{2}$ and $b=-\dfrac{4+\epsilon}{n}$


TITLE: Property of $2^n+1=xy$
QUESTION [5 upvotes]: I was wondering if the following were true. It makes sense but I'm having trouble concocting any formal reasoning.

Let $2^n+1=xy$ for some integers $x,y>1$ and $n>0$. For $a\in\mathbb{Z}^+$, does $2^a\mid (x-1)$ $\iff$ $2^a\mid (y-1)$?

Without loss of generality, one only needs to prove one direction. However, I'm not sure how to approach the problem. Here's my attempt:
 Suppose $2^a\mid (x-1)$. Then $x\equiv 1$ mod $2^a$, so
$$y\equiv xy=2^n+1\equiv 2^r+1 \hspace{5mm}(\text{mod }2^a),\hspace{5mm}\text{ where $0\leq r<a$}.$$
I'm having trouble continuing from here, since I don't know any extra information to determine $2^r\equiv 0$ mod $2^a$. It doesn't seem like that could be true for arbitrary $a\in\mathbb{Z}$, so I must have veered horribly off-track. I appreciate any help!

REPLY [0 votes]: Hint $\ $ Case $\rm\, 1\!:\ \, 2^a\mid 2^n,\, \ so\  \ mod\ 2^a\!:\,\ xy \equiv 1+2^n\equiv 1,\ $ so $\rm\,\ x\equiv 1\iff y\equiv 1.\ $
Case $\rm\,2\!:\ \ 2^n\mid 2^a,\,$ so $\rm\:xy\!-\!1 = 2^n\mid 2^a\mid x\!-\!1\:\Rightarrow\:y=1\Rightarrow\Leftarrow\,$ (by $\rm\:xy\!-\!1 > x\!-\!1\:$ for $\rm\:y>1).$


TITLE: Lagrange theorem modulo arithmetic
QUESTION [1 upvotes]: as far as I can see, Lagrange says:
IF
$p$ = prime, $p-1 = 2*q$, $q$ = prime
THEN
$g^q \mod p = 1 \implies \text{order}(g) = q$
$g^q \mod p \neq 1 \implies \text{order}(g) = p-1$
However if i try to calculate the order of $6 \mod 7$:
$p = 7 \implies p-1 = 6 = 2 \cdot 3 \implies q=3$, $6^3 \mod 7 = 6 \neq 1 \implies \text{order}(6)$ should be $6$.
But doing it by hand i get the following:
$6^1 \mod 7 = 6$
$6^2 \mod 7 = 1$
$6^3 \mod 7 = 6$
$6^4 \mod 7 = 1$
...
so order should be $2$ instead of $6$ right? What am I missing, i really don't get it :(

REPLY [0 votes]: The only divisors of $2q$ are $1$, $2$, $q$, and $2q$. Thus they are the only candidates for the order of $g$.
If $g^q\equiv 1\pmod{p}$, then $g$ has order $1$ or $q$.
If $g^q\not\equiv 1\pmod{p}$, then $g$ has order $2$ or $2q$.


text:
\section{Mathematical modelling}\label{sec:mathematical_modelling}
We consider a material body $\mathcal{B}$ composed of two constituents, tumour cells and healthy cells, which occupy a common portion of a bounded Lipschitz domain $\Omega\subset\R^d$, $d=2$, during the time $t\in[0,T]$. Nutrients such as oxygen and glucose in $\Omega$ nourish both the healthy and tumour cells. The increasing number of tumour cells by the intake of nutrients and interaction with the surrounding healthy cells increases the mechanical stress which in turn affects the mobility of the tumour.
Treatment for cancer is given by chemotherapy in which the drug diffuses through the region $\Omega$ and kills the fast growing cancerous cells.
\par The quantities of interest to us are as follows:

the mass density of the tumour cells per unit volume $\rho\phi$, where $\phi:\Omega\times[0,T]\rightarrow[0,1]$ is the volume fraction of the tumour cells in $\mathcal{B}$ and $\rho$ is the mass density of the tumour cells,
the displacement field $\vec{u}:\Omega\times[0,T]\rightarrow\R^d$, the mass density of the nutrients $\psi:\Omega\times[0,T]\rightarrow\R$  and
 the mass density of the chemotherapeutic agents $\chi:\Omega\times[0,T]\rightarrow\R$.
\subsection{Evolution of tumour}
Time evolution of the physical system must obey the laws of conservation of mass, linear and angular momentum, energy and the second law of thermodynamics. We ignore the temperature and thermal effects and proceed as done in \cite{lima2016selection}.
\begin{itemize}
  \item \textit{Conservation of mass:}
\begin{eqnarray}\label{law:phi_mass}
    \p_t(\rho \phi) + \nabla \cdot (\rho\phi\vec{v}) = \rho (S-\nabla\cdot \vec{J}),
\end{eqnarray}
where $\vec{v}$ is the velocity field, $\rho S$ is the mass density supplied by other constituents, which encompasses proliferation of tumour cells and their death due to chemotherapy treatment, and $\rho \vec{J}$ is the mass flux over the boundary of $\Omega$ which we denote by $\p\Omega$.
  \item \textit{Conservation of linear and angular momentum:}
    \begin{equation}\label{law:linear_angular_momentum}
    \begin{aligned}
    \p_t (\rho \phi\vec{v}) + \nabla \cdot (\rho \phi\vec{v}\otimes \vec{v}) = \nabla \cdot \vec{T} + \rho\phi\vec{b} + \vec{p}, \\
    \vec{T} - \vec{T}^{t} = \vec{m},
    \end{aligned}
    \end{equation}
where $\vec{T}$ is the Cauchy stress tensor for the tumour, $\vec{b}$ is the body force, $\vec{p}$ momentum supplied by other constituents, $\vec{m}$ is the intrinsic moment of momentum, and $\vec{T}^t$ denotes the transpose of $\vec{T}$.
\end{itemize}
\par The total energy of the system $\tilde{\Psi}$ consists of the Ginzburg--Landau component $\Psi(\phi,\nabla\phi)$
depending only on $\phi$ and its gradient $\nabla\phi$, and the stored energy potential $W(\phi,\symgrad(\vec{u}))$ depending on $\phi$ and the strain measure $\symgrad(u)$. Assuming small deformations, we consider the potentials

\begin{gather*}
    \tilde{\Psi} = \int_{\Omega} \Psi(\phi,\nabla\phi) + W(\phi,\symgrad(\vec{u}))\,\dd \vec{x}, \\
    \Psi(\phi,\nabla\phi) = \frac{c}{2}\phi^2,\quad W(\phi,\symgrad(\vec{u})) = \frac{1}{2} \symgrad:\vec{C}(\phi)\symgrad + \symgrad:\overline{\vec{T}}(\phi),
\end{gather*}
where $c>0$ is a constant, $\symgrad(\vec{u}) = \frac{1}{2}(\nabla\vec{u}+\nabla\vec{u}^t)$,

 $\overline{\vec{T}}(\phi)=\lambda \phi \mathbb{I}$ is the symmetric compositional stress tensor, $\lambda>0$ depending on the tumour growth rate and $\mathbb{I}$ being the identity matrix, $\vec{C}(\phi)$ is the linear elastic inhomogeneous material tensor, and the operator $:$ denotes the inner product for second-order tensors.

\par The first variations of the energy functional with respect to $\phi$ and $\symgrad$ define the chemical potential, and the stress tensor respectively
\begin{equation}\label{eqn:chempot_stress}
\begin{aligned}
\mu = \frac{\delta \Psi}{\delta \phi} + \frac{\delta W}{\delta \phi}, \quad
\vec{T} = \frac{\delta W}{\delta \symgrad}.
\end{aligned}
\end{equation}

The effects of the elastic deformation on the movement of tumour cells is prescribed by the term $\lambda\nabla\cdot\vec{u}$ in the chemical potential.
\par To incorporate subdiffusion, we introduce fractional derivatives in the mass flux and mass sources.

Modelling the subdiffusion and proliferation of cancer cells can lead to a linear fractional partial differential equation through a comb model with proliferation in one dimension \cite{iomin2015continuous}. On a microscopic level, subdiffusion-limited reaction is modelled in \cite{seki2003recombination, yuste2004reaction} by having fractional derivatives in flux and reaction terms.

\par Motivated by the previous models, the subdiffusion limited reaction for tumour mass density takes the form,
\begin{equation}\label{frac_flux_source}
\begin{aligned}
    \vec{J}=-M_\phi(\vec{x})\p_t^{1-\alpha} \nabla\mu,\quad
    S = N_\phi \p_t^{1-\alpha}f(\phi,\psi) - P_\phi \p_t^{1-\alpha} g(\phi,\chi),
\end{aligned}
\end{equation}
for $\alpha\in(0,1)$, where $M_\phi:\Omega\rightarrow\R^+$ is such that $cM_\phi$ is the mobility of tumour cells, $f,g:\Omega\times[0,T]\rightarrow\R$ model the uptake of nutrient and chemotherapic by the tumour cells, $N_\phi>0$ is the rate at which the tumour cells proliferate by using the nutrients, and $P_\phi>0$ is the rate at which the tumour cells die due to the chemotherapy treatment. The operator $\p^{1-\alpha}_t$ is the Riemann--Liouville fractional derivative and is defined for a function $\varphi:\Omega\times[0,T]\rightarrow\R$, as

\begin{equation}\label{Def:frac_Riemann}
    \begin{aligned}
    \p^{\alpha}_t \varphi(t)& = \partial_t (g_{1-\alpha}*\varphi)(t),

    \end{aligned}
\end{equation}
where the kernel is defined by
$$g_{\alpha}(t):=\begin{cases}t^{\alpha-1}/\Gamma(\alpha),&\alpha>0,\\\delta(t),&\alpha=0,\end{cases}$$ where $\delta(t)$ is the Dirac delta distribution and $*$ denotes the convolution on the positive halfline with respect to the time variable, i.e., $(g_\alpha*\varphi)(t)=\int_0^t g_\alpha(t-s)\varphi(s)\dd s$. If $\varphi$ is sufficiently smooth, then we have
\begin{equation}\label{Def:frac_Caputo}
    \begin{aligned}

    \p^{\alpha}_t(\varphi(t)-\varphi_0)=g_{1-\alpha}*\p_t\varphi(t),

    \end{aligned}
\end{equation}
where \revision{$\varphi_0$ is a given initial value}. The right-hand side is the classical Caputo fractional derivative. The formulation on the left hand side, which expresses the Caputo fractional derivative in terms of Riemann--Liouville fractional derivative, has the advantage that it requires less regularity of $\varphi$ than the classical definition.

\par We reduce the complexity of the system by using the common simplifying assumptions as in \cite{lima2016selection}: the tumour and healthy cells have constant mass density $\rho=\rho_0$, $\vec{m}=\vec{0}$, i.e., the material is monopolar, no body force, i.e., $\vec{b}=\vec{0}$,

\revision{we neglect the terms with $\vec{v}\otimes \vec{v}$ and $\vec{p}$ by not considering the inertial effects, and we further assume that the mechanical equilibrium is attained on a much faster time scale than diffusion takes place, i.e., the term $\rho_0 \vec{v} \p_t \phi$ on the left hand side in the linear momentum equation vanishes}. For ease of technical difficulties, we assume that the tumour is an isotropic and homogeneous $\vec{C}(\phi)\equiv \vec{C}$ material, and so $\vec{C}$ takes the form $\vec{C}\symgrad=2G\symgrad+\frac{2G\nu}{1-2\nu} \tr\symgrad \mathbb{I}$, where $G>0$ denotes the shear modulus, while $\nu<\frac{1}{2}$ is the Poisson's ratio. This assumption assures that the energy functional $W(\phi,\symgrad)$ is convex in both its variables, which is required in using Lemma \ref{Lem_basic_inequality_fractional_chain_rule}, which is an analogous result to the chain rule in fractional derivatives for providing estimates for $\phi$. A more general energy functional is considered in \cite{garcke2003cahn} with integer-order derivatives.

\par

\revision{Along with these assumptions, we integrate \eqref{law:phi_mass}, take $\p_t^\alpha$ on both sides of equation \eqref{law:phi_mass} and use the semigroup property of the kernel $g_{1-\alpha}*g_{\alpha}=g_1=1$ in the following way
    \begin{equation*}
        g_{1-\alpha}*\p_t(g_{\alpha}*\varphi) = \p_t(g_{1-\alpha}*g_{\alpha}*\varphi) - g_{1-\alpha}(t)(g_{\alpha}*\varphi)(0) = \p_t(1*\varphi) = \varphi,
    \end{equation*}
    assuming sufficient smoothness on the functions, to obtain from \eqref{law:phi_mass}--\eqref{frac_flux_source} the system}
\begin{subequations}\label{model:tumour_mechcoup_frac}
    \begin{align}
     \p_t^{\alpha}(\phi-\phi_0) &= \nabla\cdot\left( M_\phi(\vec{x}) \nabla\mu\right) + N_\phi f(\phi,\psi) - P_\phi g(\phi,\chi),\label{model:tumour}\\
     \mu &= c\phi + \lambda \nabla\cdot \vec{u},\label{model:chemicalpot}\\
    \vec{0}&=\nabla \cdot \left( 2G \symgrad(\vec{u}) + \frac{2G\nu}{1-2\nu}\text{tr}(\symgrad(\vec{u}))\mathbb{I} + \lambda \phi\mathbb{I}\right),\label{model:displacement}
    \end{align}
\end{subequations}
where $\phi_0$ is a given data, playing the role of initial condition.

\begin{remark}\label{rem:decoupling}
If we assume smoothness on all involved variables, we can formally take the divergence in the deformation equation \eqref{model:displacement} and conclude that
$$\left(G + \frac{G}{1-2\nu} \right)\Delta (\nabla\cdot \vec{u}) = -\lambda \Delta \phi.$$
If $M_\phi$ is a constant, then, by substitution into \eqref{model:tumour}--\eqref{model:chemicalpot}, we obtain
$$\p_t^{\alpha}(\phi-\phi_0) = M_\phi  \left( c  - \frac{\lambda^2(1-2\nu)}{2G(1-\nu)}\right)  \Delta \phi + N_\phi f(\phi,\psi) - P_\phi g(\phi,\chi).$$

We note that in this case the equation for $\phi$ is independent of $\vec{u}$. Further this also suggests that we require at least $c > \frac{\lambda^2(1-2\nu)}{2G(1-\nu)}$ to conclude the existence of a solution. We indeed see in Section \ref{sec:mathematical_analysis} that we require a slightly stronger condition on $c$.
\end{remark}

\subsection{Evolution of nutrient}
The nutrient mass density is assumed to obey a reaction-diffusion equation, as standard \cite[Ch. 5 and 10]{preziosi}
\begin{equation}\label{model:nutrient}
    \frac{\partial \psi}{\partial t}=\nabla\cdot (M_\psi(\vec{x}) \nabla \psi) + S_\psi(\vec{x},t) - N_\psi f(\phi,\psi),
\end{equation}
$M_\psi:\Omega\rightarrow\R^+$ is the mobility of the nutrients, $S_\psi:\Omega\times[0,T]\rightarrow\R$ denotes the external source of nutrients over the volume, $N_\psi>0$ denotes the rate at which nutrients are consumed by the tumour cells, $f(\phi,\psi) =  \frac{\phi (1-\phi) \psi}{K_\psi+\psi}$ is a monod equation combined with the term $(1-\phi)$ that ensures that $\phi$, which is a volume fraction, does not take values greater than 1.

The parameter $K_\psi>0$ is the monod half saturation constant, corresponding to that nutrient mass density, where the nutrient-dependent growth takes its half maximum value.

\subsection{Evolution of chemotherapy}
The mass density of chemotherapy is assumed to be governed by a reaction-diffusion equation
\begin{equation}\label{model:chemo}
\frac{\partial \chi}{\partial t}=\nabla\cdot (M_\chi(\vec{x}) \nabla \chi) - N_\chi \chi + S_\chi(\vec{x},t) - P_\chi g(\phi,\chi),
\end{equation}
where $M_\chi:\Omega\rightarrow\R^+$ is the mobility of chemotherapeutic agents, $S_\chi:\Omega\times[0,T]\rightarrow\R$ is the external supply of chemotherapeutic over the domain, $N_\chi>0$ is the rate at which the chemotherapeutic agents are degraded, $P_\chi>0$ denotes the rate at which chemotherapeutic agents act and are blocked later by killing tumour cells. The term $g(\phi,\chi) = \frac{\phi (1-\phi) \chi}{K_\chi+\chi}$ includes, analogously to the nutrient uptake, a saturation effect, including also that mainly cells in a certain growth phase are sensible to the chemotherapy. The parameter $K_\chi>0$ is the density of chemotherapeutic agents when they reach their half maximum value.

\par Finally, collecting \eqref{model:tumour_mechcoup_frac}--
\eqref{model:chemo}, the tumour evolution is governed by the system
\begin{subequations}\label{System}
    \begin{align}
     \p_t^{\alpha}(\phi-\phi_0) &= \nabla\cdot\left( M_\phi(\vec{x}) \nabla\mu\right) + N_\phi f(\phi,\psi) - P_\phi g(\phi,\chi),\label{eqn:phieq}\\
     \mu &= c\phi + \lambda \nabla\cdot \vec{u},\label{eqn:mueq}\\
    \vec{0}&=\nabla \cdot \left( 2G \symgrad(\vec{u}) + \frac{2G\nu}{1-2\nu}\text{tr}(\symgrad(\vec{u}))\mathbb{I} + \lambda \phi\mathbb{I}\right),\label{eqn:ueq}\\%
     \frac{\partial \psi}{\partial t}&=\nabla\cdot (M_\psi(\vec{x}) \nabla \psi) + S_\psi(\vec{x},t) - N_\psi f(\phi,\psi),\label{eqn:psieq}\\
     \frac{\partial \chi}{\partial t}&=\nabla\cdot (M_\chi(\vec{x}) \nabla \chi) - N_\chi \chi + S_\chi(\vec{x},t) - P_\chi g(\phi,\chi), \label{eqn:chieq}
    \end{align}
\end{subequations}

in $\Omega$. We add to this system the following initial and boundary conditions
 \begin{subequations}  \label{eq:data}
 \begin{align}
 (\phi,\psi,\chi) = (\phi_0,\psi_0,\chi_0)&\text{ on } \Omega\times\{t=0\},\label{eq:IC}\\
 \nabla \mu\cdot\vec{n}=0&\text{ on } \revision{\p\Omega}\times (0,T), \label{eq:Neumannbc_mu}\\
 \vec{u}=0&\text{ on } \Diru\times (0,T),|\Diru|>0,\label{eq:Dirichletbc_u}\\
 \left(2G \symgrad(\vec{u}) + \frac{2G\nu}{1-2\nu}\text{tr}(\symgrad(\vec{u}))\mathbb{I} + \lambda \phi\mathbb{I}\right)\cdot\vec{n}= \vec{0} &\text{ on } \p\Omega\backslash\Diru\times (0,T),\label{eq:Neumannbc_u}\\

\psi=\tilde{\psi_b},\  \chi=\tilde{\chi_b}&\text{ on } \Dirpsi\times (0,T),\label{eq:Dirichletbc_psi_chi}\\
 M_\psi\nabla\psi\cdot \vec{n}=\psi_{b},\ M_\chi\nabla\chi\cdot \vec{n}=\chi_{b} &\text{ on } \p\Omega\backslash\Dirpsi\times (0,T),\label{eq:Neumannbc_psi_chi}
 \end{align}
 \end{subequations}
 where $\vec{n}$ denotes the outer normal to $\Omega$ and $\Diru,\Dirpsi\subset\p\Omega$ are parts of the boundary $\p\Omega$ with non-zero measures.

We assume no-flux boundary conditions for the chemical potential and Dirichlet--Neumann mixed boundary condition for the displacement, density of nutrient and chemotherapy. The homogeneous Dirichlet condition on the part of the boundary $\Diru$ for displacement accounts for the presence of a rigid part of the body such as bone, which prevents the variations of the displacement. The non-homogeneous Dirichlet condition on part of the boundary $\Dirpsi$ for density of nutrient and chemotherapy accounts for the concentration supply from blood vessels. In the rest of the boundary, the more natural Neumann boundary condition is applied. \revision{The problem with non-homogeneous Dirichlet condition for density of nutrient and chemotherapy can be converted to a problem to have homogeneous Dirichlet boundary conditions by taking $\tilde{\psi}=\psi-\tilde{\psi}_b$ and similarly for $\chi$. Therefore we assume $\tilde{\psi}_b=\tilde{\chi}_b=0$.}

 \begin{remark}
We assumed that the displacement $\vec{u}$ vanishes on $\Diru \subset \partial \Omega$, that means we imposed a homogeneous Dirichlet boundary on $\Diru$, see also  \cite{garcke2019phase, faghihi2020coupled,bartkowiak2005cahn} for the same choice of boundary behaviour. This allows us to apply the well-known Korn inequality directly in the proof of existence. In the case of pure Neumann boundary conditions, one has to consider a different solution space for $\vec{u}$, see \cite{miranville2001long, miranville2003generalized,garcke2005cahn} for more details.

 \end{remark}


text:
\begin{document}
\maketitle

\begin{abstract}
\input{000abstract.tex}
\end{abstract}

\section{Introduction}
\label{sec:intro}
\input{010introduction.tex}

\section{Related Work}
\label{sec:related}
\input{020related.tex}

\section{Notation and Preliminaries}
\label{sec:not}
\input{030pre.tex}

\section{Core Idea: Reducible Gradient Perturbation Sequence}
\label{sec:rgps_ci}
\input{040core.tex}

\section{\radagrad achieves an \texorpdfstring{$\tilde{O}(T^{-1/2})$}{TEXT} Convergence Rate}
\input{050revised.tex}

\section{Experiments}
\label{sec:experiments}
\input{060expr.tex}

\input{070conclu.tex}

\bibliographystyle{unsrt}
\bibliography{references}

\newpage
\appendix
\input{0X0appendix.tex}

\end{document}