TITLE: Injection from computable numbers into natural numbers
QUESTION [2 upvotes]: Each Turing machine which writes an infinite sequence of 1 and 0 can be regarded as
representing a (computable) real number (and of course each Turing machine represents a natural number by its machine table, or program).
The question is, how many Turing jumps do we need to construct an injection from such
computable numbers into natural numbers.
Since there are an infinite number of Turing machines that compute one and the same computable real number, it seems we need at least one Turing jump. Is only once is enough?
If not, how many? Even transfinite times?

REPLY [0 votes]: I believe one Turing jump suffices. Each computable real could be represented by a {0,1}-Turing machine therefore they are no more than the total number of computer programs.
Consider the map $f: \omega-K \mapsto \omega$ where $K$ is the halting set, that is set of the indices of all programs that halt. Then $f(x)=least \ e\in \omega$ such that $\exists z \varphi_e(z)\neq \varphi_i(z) \forall i<e$. This is a $\Sigma_1^0$ question thus can be answer in $\emptyset'$. Then it is easy to see the construction of f is recursive in $K$ and it is injective also.


TITLE: Calculating Lambert W Function
QUESTION [1 upvotes]: I'm trying to evaluate Lambert W Function , I used the formula
$$ W(z)e^{W(z)} = z \Rightarrow W(z) = \frac{z}{W(z)} $$
$$ W(z) \approx ln(z)-ln(ln(z)-ln(...)) $$
But the result is very bad If I used ln(z)-ln(ln(z)) when it used in solving many cases like :
$$ x^x = 100 \Rightarrow x = 3.015 $$ which is very far from 100 (27.68) . It works good in big values as shown in this graph.
I've searched and found this equation :
$$ L_1 - L_2 + \sum_{i=0}^\infty \sum_{j=0}^\infty (-1)^i \begin{bmatrix} i+j \\ i+1 \end{bmatrix}\frac{L_2\strut^{j}}{L_1\strut^{(i+j)}(i)!} \\ $$
$$ \text{Where} \hspace{10pt}L_1:ln(x), \hspace{10pt} L_2:ln(ln(x)) $$
Tried it but worse than the previous one and wolfram can't solve more than 15 terms but it can solve W(x) with prefect accuracy although i found the last formula on wolfram reference, so my question is " Is there a better way or formula to calculate Lambert W Function with high accuracy ? "
another small question : I've used Taylor Series Expansion for Lambert W Function
$$ W(z) = \sum_{k=0}^\infty (-k)^{(k-1)}\frac{z^k}{k!} $$
but the radius of convergence is 1/e ,so it's not useful for real computations and wolfram formula works for x > e is there a way to get the values in between 1/e and e

REPLY [1 votes]: $\require{begingroup} \begingroup$
$\def\e{\mathrm{e}}\def\W{\operatorname{W}}\def\Wp{\operatorname{W_0}}\def\Wm{\operatorname{W_{-1}}}$
Is there any need to reinvent the wheel?
Just use
Halley's method,
which provides successive approximations to
$w = \W(z)$ (so $z = w\exp(w)$) being
\begin{align}
w_{j+1}
&=
w_j-{\frac {w_j \exp(w_j)-z}
{\exp(w_j)(w_j+1)-
\displaystyle\frac{(w_j+2)(w_j\exp(w_j)-z)}{2w_j+2}}}
\end{align}
as an established very efficient way
to calculate
$\W(x)$.
You can even check out open-source code, like for example, specfunc/lambert.c,
which is part of GSL - GNU Scientific Library, for the details of realization.
$\endgroup$


TITLE: Would two astronauts one in a satellite, one on top of a tower have the same experience?
QUESTION [4 upvotes]: Imagine two astronauts.
The first is in a satellite in geostationary orbit around the earth.
The second is on top of a 35,786 km tall tower that has been built from the equator to space. This means that that the top of the tower is at the same height as the geostationary satellite.
Naive logic would seem to indicate that the astronaut in the satellite would experience free-fall (as do all satellites in orbit), whereas the astronaut on top of the tower would feel the effects of gravity (albeit at a reduced level due to being 35,786 km above the surface of the earth.
However, I can't believe that this would be the case. What's the flaw in my logic? Is the acceleration due to the earth's gravity zero at 35,786 km above the earth?

REPLY [6 votes]: They would both be free-falling.  YOur mountain dweller could stand on top of the mountain, and would feel no upward force from the ground.  The magic of geostationary orbit is that the gravitational force $G\frac{mM}{r^{2}}$ on an object is exactly equal to the centripetal force $\frac{mv^{2}}{r}$ required to keep something moving at the speed $v$ equivalent to one lap around the Earth per day.  Since gravity provides no less or more force than exactly the centripetal force required to keep your person above the same spot on the ground, your person would not need the ground to hold herself up.  She would feel weightless while standing on the mountaintop.


TITLE: Value of $\sin(\pi/18)$ by nested radicals
QUESTION [3 upvotes]: I recently asked a question about how to find value of $\sin(\pi/18)$ and I understand that there is no expression for $\sin(\pi/18)$ that uses the ordinary arithmetical operations. but I know that we have exact value by nested radicals. I want to know how we get this value ?
http://mathworld.wolfram.com/TrigonometryAnglesPi18.html

REPLY [1 votes]: Deriving $\sin\frac{\pi}{18}$ will be simpler when we implement Half angle cosine formula
$$2\sin\frac{\theta}{2} = \sqrt{2-2\cos\theta}$$
$$2\cos\frac{\theta}{2} = \sqrt{2+2\cos\theta}$$
and cosine angle identity
$$2\cos(\pi-\theta) = -2\cos\theta$$
We can derive as follows
$2\sin(\frac{\pi}{18})$  = $\sqrt{2-2\cos(\frac{\pi}{9})}$ = $\sqrt{2-\sqrt{2+2\cos(\frac{2\pi}{9}})}$ = $\sqrt{2-\sqrt{2+\sqrt{2+2\cos(\frac{4\pi}{9}}})}$ = $\sqrt{2-\sqrt{2+\sqrt{2+\sqrt{2+2\cos(\frac{8\pi}{9}}}})}$ = $\sqrt{2-\sqrt{2+\sqrt{2+\sqrt{2-2\cos(\frac{\pi}{9}}}})}$
Now you can observe $2\cos(\frac{\pi}{9})$ repeats in cycle infinitely
$$\therefore 2\sin\frac{\pi}{18} = \sqrt{2-\sqrt{2+\sqrt{2+...}}}$$ $[-++]$ is the repeating pattern in infinite nested radical. Refer here. Welcome to the world of infinite nested square roots of 2 as solution to cosine angles!


TITLE: Power of approximate exponential of operator sum
QUESTION [1 upvotes]: Say $A$ and $B$ are two non-commuting operators in an Hilbert space (think matrices, for instance). Then, one way we can approximate $C=e^{A+B}$ is as $\widetilde{C}=e^{A/2}e^{B/2}e^{A/2}e^{B/2}\simeq C$.
Question: Does it hold that $\widetilde{C}^{2}=e^{A}e^{B}e^{A}e^{B}$? It doesn't, right?

REPLY [0 votes]: Your question is a bit open ended and dependent on the matrix norms involved, etc... Purely formally, if you introduce a parameter t, something tells me you are interested in the Zassenhaus formula,
$$
e^{t(A+B)}= e^{tB}~ e^{tB} ~e^{-\frac{t^2}{2} [A,B]} \\ \times ~
e^{\frac{t^3}{6}(2[B,[A,B]]+ [A,[A,B]] )}   ~
e^{\frac{-t^4}{24}([[[A,B],A],A] + 3[[[A,B],A],B] + 3[[[A,B],B],B]) } \cdots
$$
Anyway, the expression you write is, by straightforward application of the CBH formula
$$
\widetilde{C}=e^{tA/2}e^{tB/2}e^{tA/2}e^{tB/2}=  (e^{tA/2}e^{tB/2})^2\\
= ( e^{\frac{t}{2}(A+B)+\frac{t^2}{8}[A,B]+ O(t^3)  } )^2 \\
 =e^{t(A+B)+ \frac{t^2}{4}[A,B]+ O(t^3) },
$$
which does not look like a good "approximation". The last step follows from the fact that the commutator of the exponent with itself upon squaring vanishes, so the squaring act has a degenerate, terminating CBH expansion!!
You could easily eliminate the $O(t^2)$ terms in the exponent by considering
$$
e^{tA/2}e^{tB}e^{tA/2}  =e^{t(A+B)+  O(t^3) },
$$
instead.


text:
\begin{document}

\begin{frontmatter}

\title{Tensor-Krylov method for computing eigenvalues of parameter-dependent matrices}

\author[KUL]{Koen Ruymbeek\corref{cor1}}
\cortext[cor1]{Corresponding author}
\ead{koen.ruymbeek@kuleuven.be}

\author[KUL]{Karl Meerbergen}
\author[KUL]{Wim Michiels}
\address[KUL]{Department of Computer Science, KU Leuven, Celestijnenlaan 200A, 3001 Leuven, Belgium}

\begin{abstract}
In this paper we extend the Residual Arnoldi method for calculating an extreme eigenvalue (e.g. largest real part, dominant,...) to the case where the matrices depend on parameters. The difference between this Arnoldi method and the classical Arnoldi algorithm is that in the former the residual is added to the subspace. We develop a Tensor-Krylov method that applies the Residual Arnoldi method (RA) for a grid of parameter points at the same time. The subspace contains an approximate Krylov space for all these points. Instead of adding the residuals for all parameter values to the subspace we create a low-rank approximation of the matrix consisting of these residuals and add only the column space to the subspace. In order to keep the computations efficient, it is needed to limit the dimension of the subspace and to restart once the subspace has reached the prescribed maximal dimension. The novelty of this approach is twofold. Firstly, we observed that a large error in the low-rank approximations is allowed without slowing down the convergence, which implies that we can do more iterations before restarting. Secondly, we pay particular attention to the way the subspace is restarted, since classical restarting techniques give a too large subspace in our case. We motivate why it is good enough to just keep the approximation of the searched eigenvector. At the end of the paper we extend this algorithm to shift-and-invert Residual Arnoldi method to calculate the eigenvalue close to a shift $\sigma$ for a specific parameter dependency. We provide theoretical results and report numerical experiments. The Matlab code is publicly available.
\end{abstract}

\begin{keyword}
Eigenvalues \sep Tensor-Krylov \sep Parameter-dependent matrices
\MSC[2010] 65F15
\end{keyword}

\end{frontmatter}

\linenumbers

\section{Introduction}

In this paper, we consider parameter-dependent eigenvalue problems
\begin{equation} \label{eqn:VW_A}
\VM A(\VM \omega) \VM x(\VM \omega) = \lambda(\VM \omega) \VM x(\VM \omega), \quad \VM \omega \in \Omega, \quad \VM A: \Omega \rightarrow \R^{n \times n}: \VM \omega \mapsto \sum_{i=1}^m f_i(\VM \omega) \VM A_i
\end{equation}
on a compact parameter set $\Omega \subset \R^d$ and functions $f_i, i=1,2,\hdots,m$ where $m$ is small relative to $n$. We look for an eigenpair $(\lambda, \VM x)$ over $\Omega$ for which $\lambda$ is extreme following some ordering, e.g. we sort by modulus if $\lambda$ is the dominant eigenvalue. We assume that the wanted eigenvalue converges for any $\VM \omega$ in a few iterations of the Arnoldi algorithm. In practice, this means that $\lambda$ should be a well separated extreme eigenvalue of $\VM A(\VM \omega)$. Here we want to develop a numerical method to efficiently calculate $\lambda$ for a large set of points in $\Omega$ as calculating this eigenvalue for each point individually using standard eigenvalue solvers becomes soon infeasible. At the end of the paper, we extend our method to a case for which the eigenvalue closest to a shift $\sigma$ is wanted.

If $\VM A(\VM \omega)$ is symmetric for all $\VM \omega \in \Omega$, a sampling-based subspace method \cite{ruymbeek2019calculating} and the Successive Constraint method (SCM) \cite{sirkovic2016} were recently proposed. In the latter, the desired eigenvalue is approximated by iteratively computing lower-and upper bounds for each parameter value. Another class of solvers are based on polynomial approximation, see for example \cite{Andreev2012}. This technique is also used for the related problem of calculating statistical moments if the parameter is stochastic, see \cite{Benner2019}. The disadvantage of using polynomials is that they suffer from possible lack of smoothness of the desired eigenvalue function \cite{Fenzi2019}.

The method that we propose is an example of a Tensor-Krylov method. In these methods, the iteration vectors in the Krylov method are represented by a low rank tensor. Tensor-Krylov techniques are already discussed for various problems, for example to solve a large system of non-linear equations \cite{Bouaricha1996, Bouaricha1998} and parametrized linear systems \cite{Kressner2011}.

There are two important aspects in our method. First the iteration vector is represented by a tensor which is truncated to a low-rank representation, for reasons of storage and computational efficiency. Second, the Ritz pairs are obtained from projection of $\VM A(\VM \omega)$ on the sum of the Krylov spaces for all sample points. A practical issue is the growth of the dimension of this subspace as we add multiple vectors in each iteration. Therefore we often need to restart to keep the algorithm feasible. We gain in efficiency if we need to add less vectors to the subspace (and so introducing a larger error) and have a good way to restart the subspace. We demonstrate that adding fewer vectors in a careful way does not lead to slower convergence and we give an explanation. We give practical bounds for the introduced errors.

The plan of the paper is as follows. In Section 2 we discuss the Arnoldi method and analyse how errors are propagated in the algorithm and give convergence results. We explain the proposed algorithm in the one parameter case and discuss why it is an extension of the Residual Arnoldi algorithm in Section 3. Subsequently, we illustrate the convergence and error propagation with some numerical experiments. In Section 5, we revisit the method for multiple parameters and report a numerical experiment. We extend our method to shift-and-invert Residual Arnoldi when all matrices except one is of rank one and give some numerical examples. We finish this paper with some concluding remarks. The Matlab code can be downloaded following \href{http://twr.cs.kuleuven.be/research/software/delay-control/TensorKrylov}{this link}.

Throughout the paper, we denote vectors by small bold letters and matrices by capital bold letters. The subspace of $\R^n$ spanned by the columns of an orthonormal matrix $\VM V$ is denoted by the calligraphic letter $\subs{V}$. An index $k$ may be added to refer to the iteration where this subspace is made.
With calligraphic letters other than $\subs{V}$, we denote tensors of order strictly larger than two. Further, we denote by $\VM{I}$ the identity matrix of appropriate dimensions and by the function $\text{diag}(\VM{v})$ a diagonal matrix with the elements of the vector $\VM{v}$ on its diagonal. With $\norm{.}$ we denote the Euclidean norm and we use  $\norm{.}_F$ to refer to the Frobenius norm.

We end this section by introducing the concept of reduced eigenvalue problem. Let $\subs{V} \subset \R^n$ be a subspace and $\VM V$ an orthonormal matrix where the columns form a basis for $\subs{V}$, then we call $\VM{V}^T \VM{A} \VM{V}$ the \emph{reduced eigenvalue problem}. We denote by $\left( \lambda^\subs{V}, \VM{x}^\subs{V} \right)$ an eigenpair of this reduced eigenvalue problem and the corresponding Ritz vector by $\VM y = \VM V \VM x^\subs{V}$. A subindex may be used to indicate the iteration where the subspace is made. We define the corresponding residual as $\VM r = \VM A \VM y - \lambda^\subs{V} \VM y$. We further assume that an eigenvector has norm $1$. All proposed algorithms in this paper are implemented in Matlab version R2018a. All experiments are performed on an Intel i5-6300U with $2.5$ GHz and $8$ GB RAM.

\section{Arnoldi method}
In this section we introduce the Arnoldi algorithm in a general form, we explain how convergence is measured and give some results regarding its performance when the iterations are contaminated by error. The Arnoldi algorithm itself can be found in pseudo-code in Algorithm \ref{algo:Arn}. In the Arnoldi algorithm, a subspace $\subs{V}$ is iteratively built and we look for an approximation of the searched eigenvector in this subspace. In the algorithm, we need to decide about the choice of the continuation vector $\VM t$, see line \ref{line:tk} in Algorithm \ref{algo:Arn}, which, after multiplying by $\VM A \VM V$, is added to the subspace.

\begin{algorithm}[h]
\hspace*{\algorithmicindent} \textbf{Input:} Matrix $\VM{A}$, normalized vector $\VM{u} \in \R^n$ \\
\hspace*{\algorithmicindent} \textbf{Output:} Approximate eigenpair $(\hat \lambda, \VM{y})$ of $\VM A$
\begin{algorithmic}[1]
\State $\VM w^\bot = \VM u$
\For{$k=1$ to $n-1$}
\State Add $\VM w^\bot$ to the subspace $\subs V$
\State $\VM w = \VM A \VM V \VM t$  \label{line:tk}
\State $\VM w^\bot =  \left( \VM I - \VM V \VM V^T \right) \VM w$
\State Calculate eigenpair  $\left( \lambda^{\subs{V}}, \VM x^{\subs{V}} \right)$ of $\VM V^T \VM A \VM V$
\State Break down if $\left( \lambda^{\subs{V}}, \VM x^{\subs{V}} \right)$ meets a certain convergence criterium
\EndFor
\State $\hat \lambda = \lambda^{\subs{V}}$ and Ritz vector $\VM y = \VM V \VM x^{\subs{V}}$
\end{algorithmic}
\caption{ Arnoldi algorithm}  \label{algo:Arn}
\end{algorithm}

It is not difficult to prove that, in exact arithmetic, the eigenvalue approximations do not depend on the used continuation vector as long as the last element in $\VM t$ is non-zero. This changes when an error (an intended error or a numerical error) is introduced in one of the iterations. In these cases the choice of the continuation vector may be decisive for convergence.

In the original version of the Arnoldi algorithm \cite{Arnoldi1951}, the $k$th standard basis vector is used as continuation vector in the $k$th iteration, but this choice performs badly if there is a significant error in one of the iterations, see \cite{Lee2007}. In the next section, another continuation vector is introduced which appears to be more robust against introduced errors. We motivate our choice for the continuation vector and give some convergence results. The convergence behaviour under errors is closely related to inexact Krylov methods, see \cite{Simoncini2003, Simoncini2005}.

In this section we want to answer the question which vector, after multiplying by $\VM A$, is the best to add to the subspace when this is a Krylov space contaminated with errors. This was the research topic in the papers \cite{Ye2008, Wu2014}. In both papers, the authors try to search for the continuation vector $\VM{t}$ that maximizes
\begin{equation} \label{eqn:VW_tk}
 \cos \left( \theta( \VM x, \text{Range} \{ \VM V, \VM A \VM V \VM t \} ) \right) := \max_{ \VM z \in \text{Range} \{ \VM V, \VM A \VM V \VM t \} } \cos( \theta( \VM x, \VM z))
\end{equation}
The following result is proved in \cite{Ye2008}.
\begin{theorem} \label{the:optim_tk}
 Let $\VM R = \VM A \VM V - \VM V \left( \VM V^T \VM A \VM V \right)$, then the continuation vector $\VM t$ that maximizes \eqref{eqn:VW_tk} is
\begin{equation} \label{eqn:optim_tk}
\VM t = \VM R^+ \VM x + \VM b
\end{equation}
with $\VM b$ in the nullspace of $\VM R$, $\VM R^+$ the pseudo-inverse of $\VM R$ and $\VM x$ the eigenvector associated with $\lambda$.
\end{theorem}

Although this result gives us insight in the choice of the optimal continuation vector, it is not computable as the formula contains the exact eigenvector. The author of \cite{Ye2008} argued that in practice, the eigenvector \begin{equation} \label{eqn:t_k}
\VM t = \VM x^{\subs{V}}.
\end{equation}
associated with the desired eigenvalue of the reduced problem $\VM V^T \VM A \VM V$ is a good practical choice for the continuation vector. In the remainder of this paper, this is used as continuation vector, despite the fact that in the case of large non-Hermitian matrices, it may be better to use the \emph{refined Ritz vector} \cite{Wu2014}.

Before the suitability as continuation vector was discussed in \cite{Ye2008}, it was already observed that this choice leads to a more robust version of the Arnoldi algorithm and was called \emph{Residual Arnoldi}  \cite{Lee2007, LeeStewart2007}. Throughout this paper, we use the same name. The name comes from the property that the vector with which we expand the subspace is the residual vector as follows from Lemma \ref{lem:resid}. The following result is well-known in the literature \cite[Proposition 6.8]{Saad2011} and is frequently used in the analysis of Krylov methods although it holds for all subspace methods.

\begin{lemma} \label{lem:resid}
Let $\left( \lambda^\subs{V}, \VM y\right)$ be a Ritz vector of $\VM A$ with respect to the subspace $\subs{V}$, then the residual satisfies
\begin{equation} \label{eqn:resid}
\VM r = \left( \VM I - \VM V \VM V^T \right) \VM A \VM y
\end{equation}
so the residual is the result of applying matrix $\VM A$ on $\VM y$ orthogonalized with respect to $\VM V$.
\begin{eproof}
We can express
\begin{align}
\VM r & = \VM A \VM y - \lambda^{\subs{V}} \VM y \nonumber \\
& = \VM A \VM V \VM x^{\subs{V}} - \VM V \lambda^{\subs{V}} \VM x^{\subs{V}} \nonumber \\
& = \VM A \VM V \VM x^{\subs{V}} - \VM V \left( \VM V^T \VM A \VM V \right) \VM x^{\subs{V}} \nonumber \\
& = \left( \VM I - \VM V \VM V^T \right) \VM A \VM V \VM x^{\subs{V}}. \nonumber
\end{align}
\end{eproof}
\end{lemma}

Lemma \ref{lem:resid} implies that choosing \eqref{eqn:t_k} as continuation vector can simplify the Arnoldi algorithm (see Algorithm \ref{algo:Arn}) provided the norm of the residual is used as convergence criterium.

The Residual Arnoldi algorithm is shown in Algorithm \ref{algo:resid_Arn}.

\begin{algorithm}[h]
\hspace*{\algorithmicindent} \textbf{Input:} Matrix $\VM{A}$, initial vector $\VM{u} \in \R^n$ \\
\hspace*{\algorithmicindent} \textbf{Output:} Approximate eigenpair $(\hat \lambda, \VM{y})$ of $\VM A$
\begin{algorithmic}[1]
\State $\VM r= \VM A \VM u - \lambda^{\VM u} \VM u$, with $\lambda^{\VM u}$ the Ritz value associated to the subspace .
\For{$k = 1$ to $n$}
\State add $\VM r$ to the subspace $\subs{V}$ \label{line:V_k_resarn}
\State Calculate eigenpair $\left( \lambda^{\subs{V}}, \VM x^{\subs{V}} \right)$ of $\VM V^T \VM A \VM V$
\State $\VM y = \VM V \VM x^\subs{V}$ \label{line:w_resarn}
\State Calculate the residual $\VM r$ of $\left( \lambda^\subs{V}, \VM y\right)$ via Lemma \ref{lem:resid} \label{line:resid}
\State Stop if $\norm{ \VM r}$ is small enough
\EndFor
\State $\hat \lambda = \lambda^{\subs{V}}$ and Ritz vector $\VM y = \VM V \VM x^{\subs{V}}$
\end{algorithmic}
\caption{ Residual Arnoldi (RA) Algorithm}  \label{algo:resid_Arn}
\end{algorithm}

\subsection{Convergence of the Residual Arnoldi algorithm} \label{sect:conv_res_arn}

In this section, we model the error in the RA algorithm, we discuss some convergence results and give tolerances for the different errors such that convergence is assured.
We model the possible error $\VM e$ as an error on $\VM r$, i.e on the vector that we want to add to the subspace. The error $\VM e$ consists of two sources: on the one hand, the error $\VM e_1 $ on $\VM x^\subs{V}$ (before the orthogonalization on line \ref{line:w_resarn} in the RA algorithm) and, on the other hand, the error $\VM e_2$ on $\VM r$ (after the orthogonalization on line \ref{line:resid}). Note that $\VM e_2 $ is already orthogonal to $ \subs{V}$, because otherwise the added vector is no longer orthogonal to $\subs{V}$. To resume, the vector that we add to the subspace can be written as
\begin{align}
\hat{ \VM r} & = \left( \VM I - \VM V \VM V^T \right) \VM A \VM V \left( \VM x^\subs{V} + \VM e_1\right) + \VM e_2  \nonumber \\
& = \VM r + \underbrace{ \underbrace{\left( \VM I - \VM V \VM V^T \right) \VM A \VM V \VM e_1}_{\VM e_1^\bot} + \VM e_2}_{\VM e}. \label{eqn:e_w_bot}
\end{align}

The convergence of this method was analysed in \cite{Lee2007, LeeStewart2007}. The authors assume that in each iteration $k$ the error $\xkb e$  relative to the residual, is bounded by some $\eta$, in the sense
\begin{equation}
\dfrac{ \norm{ \xkb e}}{ \norm{ \xkb{r}}} \leq \eta. \label{eqn:rel_error}
\end{equation}
The main result on which the analysis relies, is that in iteration $k$ there always exists $\VM B_k$ such that $k$ iterations of the original Arnoldi method on the matrix
$$ \hat{ \VM A}_k := \VM{A} + \VM B_k,$$
in the absence of errors results in the same subspace as performing $k$ iterations of RA on matrix $\VM A$ where in each iteration an error $\xkb{e}$ is made. In these papers, the authors make the assumption that the size of $\norm{ \VM B_k}$ is proportional to $\eta$, i.e. there exist $C_\eta$ such that
$$\norm{ \VM B_k} \leq \eta C_\eta.
$$
The authors indicate that they cannot prove this assumption, although they have some empirical evidence that it is satisfied. Other assumptions were made on the separation of the spectrum of both $\VM A$ and $ \hat{ \VM A}$.

\begin{theorem}
Assume that the above conditions are satisfied, let $C_X$ be a constant depending on $C_\eta$ and the spectrum of $\VM A$.  If $2 \eta C_X < 1$, then
$$\norm{ \xkb y - \VM x} \leq  \dfrac{\kappa_k}{1-2 \eta C_X}$$
where $\kappa_k$ is a bound for the difference between the approximation and the searched eigenvector when no error is made. This value goes to zero as $k$ increases.
\end{theorem}

We refer to \cite{LeeStewart2007} for a proof and more details. In practice this bound is not (or hardly) computable, so a workable value of $\eta$ needs to be used in practice.

If the method converges, then the question arises how it depends on $\eta$. Empirically, the authors of \cite{Lee2007, LeeStewart2007} noticed that the speed of convergence for $\eta$ lower than or equal to $10^{-3}$ stays unchanged but that for $1 \gg \eta > 10^{-3}$ the convergence can slow down a bit. For $\eta$ around $1$, we do not have convergence, which can be explained by the fact that the error is as large as the residual.

We now bound the norm of the errors $\VM e_1$ and $\VM e_2$ such that convergence is still preserved, i.e. such that condition \eqref{eqn:rel_error} is fulfilled with $\eta < 10^{-3}$. Let
\begin{equation} \label{eqn:e_X-e_R}
\dfrac{ \norm{ \VM e_1}}{\norm{ \VM x^\subs{V}}} = \norm{ \VM e_1} \leq \eta_X \text{ and } \dfrac{ \norm{ \VM e_2}}{ \norm{ \VM r + \VM e_1^\bot}} \leq \eta_R
\end{equation}
be tolerances for the relative error, then we obtain starting from \eqref{eqn:e_w_bot}

\begin{align*}
    \dfrac{ \norm{\VM e} }{\norm{ \VM r}} & \leq \dfrac{\norm{ \VM e_1^\bot} }{\norm{ \VM r}} + \dfrac{\norm{ \VM e_2} }{\norm{ \VM r}} \\
    & = \dfrac{ \norm{ (\VM I - \VM V \VM V^T) \VM A \VM V \VM e_1}}{\norm{ \VM r}} + \dfrac{\norm{ \VM e_2} }{\norm{ \VM r}} \\
    & \leq \underbrace{ \dfrac{ \norm{\VM A} \norm{ \VM e_1}}{\norm{ \VM r}}}_{\leq \eta/2} + \underbrace{ \dfrac{ \norm{ \VM e_2}}{\norm{ \VM r}}}_{\leq \eta/2}.
\end{align*}

The aim is now to bound both terms by $\eta/2$. Remark that the first term means that the error propagation due to the matrix-vector product needs to be limited. From \eqref{eqn:e_X-e_R} it follows that
\begin{equation} \label{eqn:bound_eX}
\eta_X  \leq \dfrac{ \eta \norm{ \VM r}  }{2 \norm{\VM A}}
\end{equation}

which means that $\eta_X$ scales with the residual norm and the norm of the matrix $\VM A$. The lower the residual is, the lower $\eta_X$ needs to be. In the example in Section \ref{sect:Error_an}, we demonstrate that if we keep $\eta_X$ constant, then convergence stagnates if the residual becomes small.

Using $\norm{ \VM r + \VM e^\bot_1} \leq \norm{ \VM r} + \norm{ \VM e^\bot_1}$, we obtain
\begin{align*}
\dfrac{ \norm{ \VM e_2}}{\norm{ \VM r}} & = \dfrac{ \norm{ \VM e_2}}{ \norm{ \VM r + \VM e_1^\bot}}  \dfrac{ \norm{ \VM r + \VM e_1^\bot} }{\norm{ \VM r}}  \\
    & \leq \dfrac{ \norm{ \VM e_2}}{ \norm{ \VM r + \VM e_1^\bot}} \dfrac{ \norm{ \VM r} + \norm{e^\bot_1} }{\norm{ \VM r}} \\
    & \leq \dfrac{ \norm{ \VM e_2}}{ \norm{ \VM r + \VM e_1^\bot}} (1 + \eta/2).
\end{align*}
As this needs to be smaller than $\eta/2$, we must impose that
\begin{equation} \label{eqn:bound_eR}
    \eta_R \leq  \dfrac{\eta }{2 + \eta}.
\end{equation}

\section{Tensor-Arnoldi for one parameter} \label{sec:Tensor-Arn_1}
We are ready to propose our algorithm for calculating the wanted eigenvalue over the parameter space. First, a general description of the algorithm is given, then we describe in more detail some specific steps of the algorithm. We explain the algorithm for the case where we have one parameter. For the extension to multiple parameters we refer to Section \ref{sect:Tensors}. First, we discretize $\Omega$ into $n_1$ points $\omega_1, \omega_2, \hdots, \omega_{n_1}$. We want to compute the desired eigenvalue $\lambda$ at these points. For this, we consider the linear operator $\Ten A$ such that

\begin{equation} \label{eqn:diagAX}
\Ten A \VM x = \begin{bmatrix}
\VM A(\omega_1) & & \\
& \ddots & \\
& & \VM A(\omega_{n_1})
\end{bmatrix} \begin{bmatrix}
\VM x(\omega_1) \\
\VM x(\omega_2) \\
\vdots \\
\VM x(\omega_{n_1})
\end{bmatrix} = \begin{bmatrix}
\VM A(\omega_1) \VM x(\omega_1) \\
\vdots \\
\VM A(\omega_{n_1}) \VM x(\omega_{n_1})
\end{bmatrix} \in \R^{n n_1}.
\end{equation}
If $\VM x$ is reformulated as matrix $\VM X = [\VM x(\omega_1), \VM x(\omega_2), \hdots, \VM x(\omega_{n_1})] \in \C^{n \times n_1}$, then \eqref{eqn:diagAX} is equivalent to an operator $\Ten{A}( \VM X)$ such that $\text{vec}( \Ten{A}( \VM X)) = \Ten A \VM x$. Conceptually, our method boils down to applying RA to \eqref{eqn:diagAX} which is nothing more than applying RA for each $\omega_i$ individually.

In order to efficiently compute $\Ten A(\VM X)$ we need the expansion in \eqref{eqn:VW_A}. If we have a decomposition $\VM X = \VM U \VM Z$ with $\VM U \in \R^{n \times n_U}$ and $\VM Z \in \C^{n_U \times n_1}$, then in the same way as in \cite{Kressner2011} we obtain

\begin{equation} \label{eqn:A(X)}
\Ten{A}(\VM X) = \underbrace{[\VM A_1 \VM U, \VM A_2 \VM U, \hdots, \VM A_m \VM U]}_{\VM W_U \in \R^{n \times n_U m}} \underbrace{\begin{bmatrix}
\VM Z \cdot \text{diag}( [f_1(\omega_1), f_1(\omega_2), \hdots, f_1(\omega_{n_1})]) \\
\VM Z \cdot \text{diag}( [f_2(\omega_1), f_2(\omega_2), \hdots, f_2(\omega_{n_1})])\\
\vdots \\
\VM Z \cdot  \text{diag}( [f_m(\omega_1), f_m(\omega_2), \hdots, f_m(\omega_{n_1})])\\
\end{bmatrix}}_{\VM W_Z \in \C^{n_U m \times n_1}}.
\end{equation}

Remark that the term $$\VM Z \cdot \text{diag}( [f_i(\omega_1), f_i(\omega_2), \hdots, f_i(\omega_{n_1})]) = [ f_i(\omega_1) \VM Z(:,1), f_i(\omega_2) \VM Z(:,2), \hdots, f_i(\omega_{n_1}) \VM Z(:,n_1)],$$ so column $j$ of $\VM Z$ is multiplied with the function $f_i$ evaluated in $\omega_j$.
Let now $\VM X^\subs{V}$ be the matrix storing in each column $i$ the eigenvector $\VM x^\subs{V}(\omega_i)$ of the reduced eigenvalue problem $\VM V^T \VM A(\omega_i) \VM V$ with respect to the subspace $\subs{V}$. We first make a low-rank decomposition $\VM U \VM Z$ of $\VM X^\subs{V}$ and subsequently we compute $\VM W_U \VM W_Z= \Ten{A}( (\VM V \VM U) \VM Z)$ via \eqref{eqn:A(X)} and we orthogonalise the result with respect to $\subs{V}$ to obtain an approximate residual for each parameter value, see Lemma \ref{lem:resid}. This is analogeous to lines \ref{line:w_resarn} and \ref{line:resid} in the RA algorithm (see Algorithm \ref{algo:resid_Arn}). Instead of adding all column vectors of the residual, we make a low-rank approximation and add in the next iteration the resulting basis $\VM U$ to $\subs{V}$. The algorithm for the low-rank decomposition is the same as in \cite{Kressner2011}.

After calculating new eigenvalues and eigenvectors we can repeat the same procedure untill all residuals are below a certain tolerance $\epsilon$. As we are adding multiple vectors to the subspace in each iteration and the maximal dimension of the subspace is assumed to be prescribed, the question arises how we restart the subspace when it attains its maximal dimension. In the literature, several techniques are developed for the Arnoldi algorithm like polynomial restart methods \cite{Saad1980, Saad1984} or the implicitly restarted method \cite{Sorensen1992}. All of them try to delete the part of the subspace that corresponds to an unwanted part of the spectrum, but still keeps a significant part of the subspace. The reason is that not only the Ritz vector associated with the wanted eigenvalue, but also some other Ritz vectors whose eigenvalues are close to $\lambda^{\subs{V}}$ help improving the convergence to the wanted eigenvalue. This is known as \emph{thick restarting}. In our case of parametric eigenvalue problems, it is slightly more complicated as eigenvectors vary over the parameter. This means that for just representing the wanted eigenvalue we already need quite a large subspace. If also eigenvectors associated with other eigenvalues need to be present in the restarted subspace, there is less room for iterations within the same subspace. Therefore, we choose to just use a basis for the first eigenvectors as restarted subspace. This is a good strategy since it is known from \cite{ruymbeek2019calculating, seyranian2003multiparameter} that the derivative of an eigenvector, which is a measure for the change of the eigenvector, depends the most on the eigenvectors whose eigenvalues are the closest to $\lambda(\VM \omega)$. As we estimate the eigenvector over the parameter space, we therefore also estimate the eigenvector whose eigenvalues are the closest to $\lambda(\VM \omega)$. The algorithm can be found in Algorithm \ref{algo:resid_Arn_param}.

\begin{algorithm}[h]
\hspace*{\algorithmicindent} \textbf{Input:} $\VM{A}(\omega)$, orthonormal matrix $\VM U \in \R^{n \times p}$, tolerance $\epsilon$, tolerances $\eta_X$ and $\eta_R$ for the two decompositions, maximal dimension of subspace $n^\subs{V}$ \\
\hspace*{\algorithmicindent} \textbf{Output:} Approximate eigenpair $(\hat \lambda(\omega_i), \VM{y}(\omega_i)), i = 1,2, \hdots, n_1$
\begin{algorithmic}[1]
\State $\hat{ \VM R} = [1, 1, \hdots, 1]$, \quad $k = 1$
\While{$\norm{ \hat{ \VM R}}_F/\sqrt{n_1} > \epsilon$}
\State Add $\VM U$ to the subspace $\subs{V}$
\State Calculate eigenpair $\left( \lambda^{\subs{V}}(\omega_i), \VM x^{\subs{V}}(\omega_i) \right)$ of $\VM V^T \VM A (\omega_i) \VM V, i = 1,2, \hdots, n_1$ \label{line:reduc_eigv_probl}
\State $\VM X^{\subs{V}} = [\VM x^{\subs{V}}(\omega_1), \VM x^{\subs{V}}(\omega_2), \hdots, \VM x^{\subs{V}}(\omega_{n_1})]$  \label{line:XV}
\If{ dimension $\subs{V} > n^\subs{V}$}
\State $[\VM U, \VM Z] = \text{lowrank}( \VM X^{\subs{V}}, 10^{-12})$ \label{line:lowrank1a}
\State Make $\subs{V}$ empty
\Else
\State $[\VM U, \VM Z] = \text{lowrank}( \VM{ X}^{\subs{V}}, \eta_X)$ \label{line:lowrank1b}

\State $[\VM W_U, \VM W_Z] = \Ten{A}( \left( \VM V \VM U \right) \VM Z)$ using \eqref{eqn:A(X)} \label{line:A(X)}
\State $\VM W^\bot_U = (\VM I - \VM V \VM V^T) \VM W_U$ ( Orthogonalize $\VM W_U$ w.r.t $\VM V$)
\State Residual $ \hat{ \VM R} := [\hat{\VM{r}}(\omega_1), \hat{\VM{r}}(\omega_2), \hdots, \hat{\VM{r}}(\omega_{n_1})] = \VM W^\bot_U \VM W_Z$ \label{line:residual}
\State $[\VM U, \sim] = \text{lowrank}( \VM W^\bot_U, \VM W_Z, \eta_R)$; \label{line:lowrank2}
\EndIf
\State $k = k+1$
\EndWhile
\State $\hat \lambda(\omega_i) = \lambda^{\subs{V}}(\omega_i)$ and Ritz vector $\VM y(\omega_i) = \VM V \VM X^{\subs{V}}(\omega_i), i= 1,2,\hdots,n_1$
\end{algorithmic}
\caption{ Parametric Residual Arnoldi (PRA) Algorithm}  \label{algo:resid_Arn_param}
\end{algorithm}

For the calculation of $\VM V^T \VM A(\omega_i) \VM V, i = 1, \hdots, n_1$, we first calculate $\VM V^T \VM A_i \VM V, i = 1, \hdots, m$ and sum these smaller matrices appropriately. The operation $\VM A_i \VM V, i = 1,\hdots m$ is not expensive as these $\VM A_i$ is usually sparse, but the inner product $\VM V^T \left( \VM A_i \VM V\right)$ need $2 n^2 n_v^2$ operations with $n_v$ the dimension of $\VM V$. As $\VM V$ is built iteratively, results from the previous iteration can be used if the subspace was not restarted. The lower $\eta_X$ and $\eta_R$, the higher the number of these inner products. In the example in the Section \ref{sect:Error_an}, we count the number of these inner products and compare them for different values of $\eta_X$ and $\eta_R$.

\section{Error analysis} \label{sect:Error_an}

We illustrate the performance of the proposed algorithm through an example. We show here what the impact is of the choice of the tolerance levels $\eta_R$ resp. $\eta_X$ in the PRA algorithm and of the maximal size of the subspace $\subs{V}$. In particular, we explain why choosing $\eta_R$ relatively high gives the best results, i.e we need less iterations and less expensive matrix-vector products without affecting the convergence.

\begin{eexample} \label{ex:conv_diff}
Consider the convection-diffusion differential equation
\begin{equation} \label{eqn:conv_diff}
f = c_1 u_x + c_2 u_y + d_{11} u_{xx} + 2d_{12} u_{xy} + d_{22} u_{yy}
\end{equation}
for $(x, y) \in \Xi= (0, 1) \times (0, 1)$. The real function $u$ has the argument $(x, y)$.
The constants $c_1$ and $c_2$ are associated with the convection part and the constants $d_{11}$, $d_{12}$ and $d_{22}$ are associated with the diffusion part of the equation. It is always assumed that matrix
$$\VM D = \begin{bmatrix}
d_{11} & d_{12} \\
d_{12} & d_{22} \end{bmatrix}
$$
is positive semi-definite. For the discretisation of $\Xi$ we use second-order central finite difference formulas and we assume Dirichlet boundary conditions. In this way, we obtain a matrix $\VM A$ depending on the values $c_1$,$c_2$,$d_{11}$,$d_{12}$, $d_{22}$ such that
\begin{equation} \label{eqn:conv_diff_A}
\VM A(c_1, c_2, d_{11}, d_{22}, d_{12}) = c_1 \VM{ \hat{ A}}_1 + c_2 \VM{ \hat{ A}}_2 + d_{11} \VM{ \hat{ A}}_3 + d_{22} \VM{ \hat{ A}}_4 + d_{12} \VM{ \hat{ A}}_5.
\end{equation}
For this first example we take $c_2 = 1, d_{11} = 1.1, d_{22} = 1, d_{12} = 1$ and we let $c_1$ vary in the interval $\Omega = [-2.5, 2.5]$. We discretize this interval in $100$ points. This means we can rewrite \eqref{eqn:conv_diff} as
$$\VM A(c_1) = \VM A_1 + c_1 \VM A_2.$$
 The desired eigenvalue is the eigenvalue with largest real part, in this case the same as the eigenvalue closest to zero. The region $\Xi$ is discretised into $100$ points in each dimension, so, the resulting matrices $\VM A_i$ are of size $10^4 \times 10^4$. The results can be found in Figure \ref{fig:conv_diff_results} and Table \ref{tbl:conv_diff}. Remark that we need many iterations but tests with individual parameters values reveal that also the Arnoldi algorithm needs many iterations.
\end{eexample}
In Section \ref{sect:conv_res_arn} we observed that, in the RA algorithm, we add the residual in each iteration to the subspace and conclude that this is more robust than classical Arnoldi. In the PRA algorithm, however, the subspace is enriched with an approximate residual for all sample points by making low-rank approximations. For this, we rely on the error analysis made in Section \ref{sect:conv_res_arn}. At line~\ref{line:lowrank1b}, on the one hand, and at line \ref{line:lowrank2} of PRA algorithm on the other, we make low-rank approximations of the matrix consisting of the eigenvector resp. the approximate residual and so introducing errors $\VM E_1 := [\VM e_1(\omega_1), \VM e_1(\omega_2), \hdots, \VM e_1(\omega_{n_1})]$ and $\VM E_2 := [\VM e_2(\omega_1), \VM e_2(\omega_2), \hdots, \VM e_2(\omega_{n_1})]$.
We first discuss the error $\VM E_1$. The low-rank approximation is such that the relative error is bounded by the tolerance $\eta_X$, i.e
\begin{equation} \label{eqn:eps_X}
\dfrac{ \norm{ \VM E_1}_F}{\norm{ \VM X^\subs{V}}_F} \leq \eta_X.
\end{equation}
Since $\norm{ \VM X^\subs{V}}_F = \sqrt{n_1}$, it follows that $\norm{ \VM E_1}_F \leq \eta_X \sqrt{n_1}$.
Technically, we can only prove from this that
$$\norm{ \VM e_1(\omega_i)}^2 \leq \eta_X^2 n_1, i = 1, 2, \hdots, n_1,$$
but assuming the error is uniformly distributed over the sample points, we can assume that
\begin{equation} \label{eqn:e_X}
\norm{ \VM e_1(\omega_i)}  \lessapprox \eta_X.
\end{equation}

\begin{figure}[h]
\begin{center}
\setlength{\figW}{7cm}
  \hspace{-0.5\linewidth}
 \begin{subfigure}[h]{0.3\textwidth}

  \hspace{1.6cm}
  \input{Afbeeldingen/rel_residu_conv_diff.tex}
  \caption{}
\end{subfigure}
\hspace{0.15\linewidth}
 \begin{subfigure}[h]{0.3\textwidth}

  \hspace{1.6cm}
  \input{Afbeeldingen/rel_E1_conv_diff.tex}
  \caption{}
\end{subfigure} \\
\begin{tikzpicture}
\hspace{0.05\linewidth}
\begin{customlegend}[legend columns=2,legend style={align=left,draw=none,column sep=2ex},legend entries={$\eta_X= 10^{-3} \quad \eta_R = \dfrac{ 10^{-3} }{2 + 10^{-3} }$ , $\eta_X= 10^{-12} \quad \eta_R = 5 \times 10^{-4}$, $\eta_X= 10^{-3} \quad \eta_R = \dfrac{ 10^{-3} }{2 + 10^{-3} }$, $\eta_X= 10^{-12} \quad \eta_R = 5 \times 10^{-13}$}]
\addlegendimage{color=mycolor1, line width=1.0pt,line legend}
\addlegendimage{color=mycolor2,solid,line width=1.0pt}
\addlegendimage{color=mycolor3,solid,line width=1.0pt}
\addlegendimage{color=mycolor4,solid,line width=1.0pt}
\end{customlegend}
\end{tikzpicture}
\end{center}
\caption{Results of the PRA algorithm for different tolerances $\eta_X$ and $\eta_R$ with maximal subspace dimension $150$ for Example \ref{ex:conv_diff}.  In (a) we plot the Frobenius norm of the residual, scaled with the number of sample points. In (b) the average first error $\VM e_1$ relative to the average residual, scaled with the norm of $\VM A_0$ is shown. As long as the value is below or around $10^{-3}$, the algorithm converges as if there were no error, see the criterium in \eqref{eqn:bound_eX}. If the value is large $(> 1)$ convergence is lost. } \label{fig:conv_diff_results}
\end{figure}

Figure \ref{fig:conv_diff_results} (a) shows the average residual over the parameter and we observe that for $\eta_X = 10^{-3}$ the convergence slows down after $400$ iterations. This is related to the average relative error $\VM e_1$ plotted in Figure \ref{fig:conv_diff_results} (b), as this becomes too high such that condition \eqref{eqn:bound_eX} is not fulfilled.
Analoguously to \eqref{eqn:e_X} we can prove that $$\dfrac{ \norm{ \VM e_2(\omega_i)}}{\norm{ \VM r(\omega_i) + \VM e_1^\bot }} \lessapprox \eta_R,  i = 1, 2, \hdots, n_1.$$
As seen in \eqref{eqn:bound_eR}, we only need to make sure that $\eta_R \leq \dfrac{\eta}{2 + \eta}$, where $\eta$ can be chosen as high as $10^{-3}$. This explains why in Figure \ref{fig:conv_diff_results} (a) choosing $\eta_R = \dfrac{10^{-3}}{2+10^{-3}}$ does not harm convergence. Furthermore, it seems even that convergence is faster for higher $\eta_R$.
This phenomenon is intriguing and cannot be explained only by a relative error plot. It is generally known in Arnoldi-type algorithms that, with each restart, convergence can slow down. The higher the tolerance $\eta_R$ the less vectors are added to the subspace, so we do not need to restart so often, hence we can do more iterations within the same subspace and convergence is faster.

In Table \ref{tbl:conv_diff}, we see statistics about the computational cost. As expected, we observe that the number of restarts increases as the tolerances decrease and that the average size of the subspace is independent of the chosen tolerances. As indication of the increased computational cost when taking the tolerances smaller, we count the number of inner products when creating the reduced eigenvalue problem (line \ref{line:reduc_eigv_probl} in the PRA algorithm). A more detailed explanation on this was given at the end of Section \ref{sec:Tensor-Arn_1}.
\begin{table}[h]
\input{Tables/table_conv_diff} \label{tbl:conv_diff}
\end{table}

In the experiments in Figure \ref{fig:conv_diff_results}, we observed that if $\eta_X$ is chosen high, this harms the convergence when the norm of the error $\VM e_1$ starts to be high relative to the residual. The solution for this is to choose the tolerance $\eta_X$ in iteration $k$ relative to the residual in the previous iteration. In Figure \ref{fig:conv_diff_results_vartol}, we compare results when holding $\eta_X$ constant and when we adapt its value relatively with respect to the residual. We see that we obtain similar results. This is explained by Figure \ref{fig:conv_diff_results_vartol}(b) as we see that the error $\VM e_1$ relative to the residual remains constant.

\begin{figure}[h]

\begin{center}
\setlength{\figW}{7cm}
  \hspace{-0.5\linewidth}
 \begin{subfigure}[h]{0.3\textwidth}

  \hspace{1.6cm}
  \input{Afbeeldingen/rel_residu_conv_diff_vartol_X.tex}
  \caption{}
\end{subfigure}
\hspace{0.15\linewidth}
 \begin{subfigure}[h]{0.3\textwidth}

  \hspace{1.6cm}
  \input{Afbeeldingen/rel_E1_conv_diff_vartol_X.tex}
  \caption{}
\end{subfigure} \\
\begin{tikzpicture}
\hspace{0.05\linewidth}
\begin{customlegend}[legend columns=2,legend style={align=left,draw=none,column sep=2ex},legend entries={$\eta_X= 10^{-12} \quad \eta_R = \dfrac{10^{-3}}{2+10^{-3}}$ , $\eta_X= 10^{-3} \frac{ \norm{\hat{ \VM R}} }{\sqrt{n1} \norm{ \VM A_0}}
 \quad \eta_R = \dfrac{10^{-3}}{2+10^{-3}}$, $\eta_X= 10^{-12} \quad \eta_R = 5 \times 10^{-13}$, $\eta_X= 10^{-3} \frac{ \norm{\hat{ \VM R}} }{\sqrt{n1} \norm{ \VM A_0}}
 \quad \eta_R = 5 \times 10^{-13}$}]
\addlegendimage{color=mycolor1, line width=1.0pt,line legend}
\addlegendimage{color=mycolor2,solid,line width=1.0pt}
\addlegendimage{color=mycolor3,solid,line width=1.0pt}
\addlegendimage{color=mycolor4,solid,line width=1.0pt}
\end{customlegend}
\end{tikzpicture}
\end{center}
\caption{Results of the PRA algorithm for different tolerances $\eta_X$ and $\eta_R$ with maximal subspace dimension $150$ for Example \ref{ex:conv_diff}. We compare the results when taking a constant low tolerance $\eta_X$ with the case that $\eta_X$ is relative to the residual. In (a) we plot the Frobenius norm of the residual, scaled with the number of sample. We see that in all four cases the algorithm converges. In (b), the average error $\VM e_1$ relative to the average residual, scaled with the norm of $\VM A_1$ is shown. As expected, we see that this error relative to the residual stays constant. } \label{fig:conv_diff_results_vartol}
\end{figure}

The conclusion is that we obtain the best results when taking $\eta_R$ high and taking $\eta_X$ relative to the residual. In this way the norm of both errors $\VM e_1$ and $\VM e_2$ stay constant with respect to the residual, which preserves convergence if this constant is low enough ($< 10^{-3}$), see Section \ref{sect:Error_an}.

\section{More than one parameter} \label{sect:Tensors}

In this section, we extend the PRA algorithm to the case of more than one parameter. Not only all computations need to have an equivalent tensor counterpart, we also need a suitable representation of the tensors as the storage requirement grows exponentially with the dimension. In this section, we first introduce some notation, then we discuss the chosen tensor-format and explain how all operations are adapted.

We assume that $\Omega$ is a Cartesian product of intervals in $\R$, so $$\Omega = \Omega_1 \times \Omega_2 \times \hdots \times \Omega_d = [a^1, b^1] \times [a^2, b^2] \times \hdots \times [a^d, b^d].$$ Each dimension $i, i=1, 2, \hdots, d$ is discretised into $n_i$ points, making a total of $n_0 = n_1 n_2 \hdots n_d$ sample points where we want to approximate the desired eigenvalue. A point $\VM \omega$ on this grid is denoted as $\left( \omega_{i_1}, \omega_{i_2}, \hdots, \omega_{i_d} \right)$ for $i_j = 1,2, \hdots, n_j, j = 1,2, \hdots, d$.
All tensor equivalents are the calligraphic version of the letters used in the one parameter case. For example, the eigenvectors are stored in tensor $\Ten X \in \C^{n \times n_1 \times \hdots \times n_d}$ where
$$\Ten X(:,i_1, i_2, \hdots, i_d) = \VM x(\omega_{i_1}, \omega_{i_2}, \hdots, \omega_{i_d}).$$
We define tensor $\Ten{F} \in \R^{m \times n_1 \times n_2 \times \hdots \times n_d}$ as the tensor with all function values, so
$$\Ten{F}(i_0,i_1, \hdots, i_d) = f_{i_0} \left( \omega_{i_1}, \omega_{i_2}, \hdots, \omega_{i_d} \right)$$
for a point $\VM \omega =  \left( \omega_{i_1}, \omega_{i_2}, \hdots, \omega_{i_d} \right).$ If $d = 1$, this is the matrix
$$\VM F = \begin{bmatrix}
f_1(\omega_1) & f_1(\omega_2) & \hdots & f_1(\omega_{n_1}) \\
f_2(\omega_1) & f_2(\omega_2) & \hdots & f_2(\omega_{n_1}) \\
\vdots \\
f_m(\omega_1) & f_m(\omega_2) & \hdots & f_m(\omega_{n_1}) \\
\end{bmatrix}.$$

\subsection{Tensors}
We first state some basic operations that can be performed on tensors, then we describe how to adapt the PRA Algorithm. A detailed overview on tensors can be found in \cite{Kolda2009}.

\begin{definition}
\begin{enumerate}
\item The \textit{mode-$i$ matriciziation} or \textit{unfolding} of a tensor $\Ten{Z} \in \C^{n \times n_1 \times \hdots \times n_d}$ is a matrix $\VM{Z}_{(i)} \in \C^{n_i \times (n n_1 \hdots n_{i-1} n_{i+1}\hdots n_d)}$ where the columns are the mode-$i$ vectors. The order in which they are placed in  $\VM{Z}_{(i)}$ is chosen to be canonical.
\item The \textit{mode-$i$ multiplication} of a matrix $\VM{M} \in \R^{k \times n_i}$ with tensor $\Ten{Z}$ is defined as
$$\Ten{Y} = \left( \VM{I},  \hdots, \VM{I},\VM{M},\VM{I}, \hdots, \VM{I} \right) . \Ten{Z} \Leftrightarrow \VM{Y}_{(i)} = \VM{M} \VM{Z}_{(i)}$$
where $\VM{M}$ is at the $i+1$th position in the tuple.
\item The \emph{Hadamard product} between two tensors $\Ten U$ and $\Ten W$ of the same dimension, is obtained by the elementwise multiplication and is denoted by $\Ten U ._H \Ten W$.
\end{enumerate}
\end{definition}

If $\VM U$ is the column basis for the tensor $\Ten X$, then we can write
\begin{equation} \label{eqn:decomp}
\Ten{X} := \left( \VM{U}, \VM{I}, \VM{I}, \hdots, \VM{I} \right) . \Ten{Z}
\end{equation}
with $\VM{U} \in \R^{n \times k}$ and $\Ten{Z} = \R^{k \times n_1 \times n_2 \times \hdots \times n_d}$, where $\VM{U}$ is a basis for $\VM{X}_{(0)}$.

We look now for a suitable tensor representation such that all the computations in the iteration can be done in an efficient way. A strong requirement for the decomposition of $\Ten{X}^\subs{V}$ besides the storage requirement is that the following operations in PRA Algorithm can be computed efficiently within this decomposition:
\begin{enumerate}
\item Calculating a real basis for the column space, even if the tensor contains complex numbers (line \ref{line:lowrank1a}, \ref{line:lowrank1b} and \ref{line:lowrank2} in the PRA Algorithm). \label{enum:2}
\item Approximating a tensor up to a given tolerance (line \ref{line:lowrank1a}, \ref{line:lowrank1b} and \ref{line:lowrank2}), in the same sense as in \eqref{eqn:eps_X}.
\item Computing a decomposition of $\Ten{W}_Z$ given the decomposition of $\Ten{Z}$ (line \ref{line:A(X)}). \label{enum:1}
\end{enumerate}

The Tensor-Train format (TT-format) satisfies all the above conditions. Before we illustrate this we briefly recall the TT-format. A more detailed overview can be found in \cite{Oseledets2011}.
A tensor $\Ten{Y} \in \C^{n_0 \times n_1 \times \hdots \times n_d}$ is in TT-format if
\begin{equation}
\Ten{Y}(i_0,i_1,\hdots,i_d) = \VM G^Y_0(i_0) \VM G^Y_1(i_1) \hdots \VM G^Y_d(i_d)
\end{equation}
with $\VM G^Y_j(i_j)$ a $r_{j-1}^Y \times r_j^Y$ matrix.
The length $r^Y_{j-1}$ and width $r^Y_{j}$ of these matrices $r^Y_i, i= -1,0,\hdots, d$ are called the \emph{TT-ranks}. It can immediately be observed that $r^Y_{-1} = r^Y_{d} = 1$. We call $r^Y_i, i = 0,1, \hdots, d-1$ the \emph{interior} ranks.
The storage of the TT-representation is lower than $(d+1)n r^2$ with $r^Y_i \leq r$ and $n_i \leq n$, which is for $r$ small significantly smaller than the storage of the full tensor. The matrices $\VM G^Y_i$ can be stored in a tensor $\Ten{G}^Y_j \in \R^{r^Y_{j-1} \times n_j \times r^Y_j}, j=0,1,\hdots,d$ as $\Ten{G}^Y_j(:,i_j,:) = \VM G^Y(i_j), j = 1,2, \hdots, n_j$.

For computing a TT-decomposition of a tensor, we perform for each mode consecutively a matricization, an SVD and a reshape, see \cite[Algorithm 1]{Oseledets2011}. Approximating a tensor for a given tolerance in the same sense as \eqref{eqn:eps_X} can also be done in the TT-decomposition, see \cite[Theorem 2.2]{Oseledets2011}.

The last requirement was that we need to be able to calculate a TT-decomposition of $\Ten{W}_Z$ (the equivalent of $\VM W_Z$ in \eqref{eqn:A(X)}), directly from the decomposition of $\Ten{Z}$. To see this, we recall the case of $d = 1$. Let $\VM 1_q$ be a vector with all ones of length $q$.

We first note that $$\VM Z \cdot \text{diag}([f_i(\omega_1), f_i(\omega_2), \hdots, f_i(\omega_{n_1})]) = \VM Z ._H ( \VM 1_{n_U} [f_i(\omega_1), f_i(\omega_2), \hdots, f_i(\omega_{n_1})])$$ from which it follows that
$$\VM W_Z = \left( (\VM 1_{n_U} \otimes \VM I_m) \VM F \right) ._H \left( (\VM I_{n_U} \otimes \VM 1_m) \VM Z \right).$$
Note that we use here the left-kronecker product $\otimes$. This means that in case of multiple parameters, $\Ten{W}_Z$ can be generalised to
$$ \Ten{W}_Z = \left( (\VM 1_{n_U} \otimes \VM I_m), \VM I, \hdots, \VM I \right).\Ten{F} ._H \left( \VM I_{n_U} \otimes \VM 1_m, \VM I, \hdots, \VM I \right).\Ten{Z}.$$
From \cite{Oseledets2011}, we obtain that
\begin{equation} \label{eqn:TT_W_Z}
\left\{\begin{matrix*}[c]
\VM G^{W_Z}_j(i_j) & = & \VM G^Z_j(i_j) \otimes \VM G^F_j(i_j) & & j = 1,2, \hdots, d, i_j = 1,2, \hdots, n_j, \\
\VM G^{W_Z}_0 & = & \VM G^Z_0 \otimes \VM G^F_0 &  & j = 0.
\end{matrix*}\right.
\end{equation}

Note that the TT-decomposition is not unique. Moreover, it is possible that the ranks of the decomposition in \eqref{eqn:TT_W_Z} of the tensor $\Ten{W}_Z$ are overestimated.

If all functions $f_i, i = 1, 2, \hdots, m$ are separable, then there is no need to first calculate a TT-decomposition of $\Ten{F}$. We can get a decomposition directly from the separability of the functions $f_i$. As the functions $f_{i}$ are separable there are functions $g_{i_j}^{i}(\omega_{i_j})$ such that
$f_{i} \left( \omega_{i_1}, \omega_{i_2}, \hdots, \omega_{i_d} \right) = g_{i_1}^{i}(\omega_{i_1}) g_{i_2}^{i}(\omega_{i_2}) \hdots  g_{i_d}^{i}(\omega_{i_d})$. It can be seen that a possible TT-decomposition is of the form
\begin{equation}
\left\{\begin{matrix*}[c]
\VM G^F_0 & = & \VM I_m \\
\VM G^F_j(i_j) & = & \text{diag}( [g_1^j(\omega_{i_j}), g_2^j(\omega_{i_j}), \hdots, g_m^j(\omega_{i_j})]) & & j = 1, 2, \hdots, d-1, i_j = 1,2, \hdots, n_j, \\
\VM G^F_d(i_d) & = & [g_1^d(\omega_{i_d}), g_2^d(\omega_{i_d}), \hdots, g_m^d(\omega_{i_d})]^T & & i_d = 1,2, \hdots, n_d.
\end{matrix*}\right.
\end{equation}
In this case, all interior ranks are $m$. This means that the interior ranks of the tensor $\Ten W_Z$ are maximal $m r^Z_i, i=0, 1, \hdots, d$. Note that also here the ranks of the decomposition are possibly overestimated. The lower are the ranks of the tensor $\Ten Z$, the more efficient the other operations are.
To have this, it is important how the eigenvectors are stored in the tensor $\Ten{X}^\subs{V}$. Eigenvectors are only unique up to a complex constant, but for the ranks of the tensor $\Ten X^{\subs V}$ it makes a difference which constant is used. The best we can do is to compare the eigenvector $\VM x(\VM \omega)$ for a point $\VM \omega$ with a certain vector $\VM a$, and store instead of $\VM x(\VM \omega)$ the vector $\dfrac{(\VM a^H \VM x(\VM \omega))}{\norm{ \VM a^H \VM x(\VM \omega)}} \VM x(\VM \omega)$. Remark that if $\VM a$ and the eigenvector have no complex elements then this boils down to multiplying with $1$ or $-1$. A possible choice for $\VM a$ is to take the eigenvector corresponding to a sample point.

If there are not many points in each dimension, it might be better not to use tensors but to reshape the tensor $\Ten{X}^\subs{V}$ as a matrix and to apply the PRA algorithm as described in Algorithm~ \ref{algo:resid_Arn_param}. The advantage of using a decomposition of the tensor is to reduce the storage and to reduce the number of operations, but if the number of parameter values is low in each dimension, then we cannot gain from the TT-representation as the cost of decomposing is too high in comparison with the possible gain in the remaining steps of the algorithm.

\subsection{Numerical Experiments} \label{sect:Tensors_NumExp}
In this section we extend Example \ref{ex:conv_diff} to multiple parameters and give the convergence results. The code can be downloaded following \href{http://twr.cs.kuleuven.be/research/software/delay-control/TensorKrylov}{this link}.

\begin{eexample} \label{ex:conv_diff_3D}
Consider the convection-diffusion differential equation \eqref{eqn:conv_diff} in Example \ref{ex:conv_diff}. We consider the case where $d_{12} = 0$, we vary the variables $c_1$ and $c_2$ both in the interval $[4,6]$ and we let the variable $d := d_{11} = d_{22}$ vary in the interval $[-1.1, -0.9]$. In this way we can rewrite \eqref{eqn:conv_diff_A} into the form
$$\VM A(c_1, c_2, d) = c_1 \VM A_1 + c_2 \VM A_2 + d \VM A_3.$$
The desired eigenvalue is the eigenvalue with smallest real part, in this case the same as the eigenvalue closest to zero. The matrices are of size $10^3 \times 10^3$ and we choose $n_1 = n_2 = n_3 = 30$ resulting in a total of $n_0 = 27 \times 10^{3}$ parameter values. As tolerances we choose $\eta_R = \dfrac{10^{-3}}{2+10^{-3}}$ and $\eta^{(k)}_X = 10^{-3} \norm{ \hat{\VM R}^{(k)}}_F/\sqrt{n_0}$
\end{eexample}
We first solve the problem discretised in each dimension by $10$ interior discretisation points. We see the result in Figure \ref{fig:conv_diff_3D}~(a). We observe that we need many iterations for convergence to the solution which is not surprising as using Arnoldi's method for all parameter points separately needs approximately the same number of iterations.  The obtained subspace $\subs{V}$ is then used as initial subspace when discretising each dimension in $30$. We see in Figure \ref{fig:conv_diff_3D} (b) that we only need a few iterations.

\begin{figure}
\setlength{\figW}{7cm}
 \begin{subfigure}[b]{0.5\textwidth}

  \input{Afbeeldingen/conv_diff_10_meanresidu.tex}
  \caption{}
\end{subfigure}
 \begin{subfigure}[b]{0.5\textwidth}

  \input{Afbeeldingen/conv_diff_30_meanresidu.tex}
  \caption{}
\end{subfigure}
\caption{Results of the PRA algorithm with maximal subspace dimension $100$ for Example \ref{ex:conv_diff_3D}. In (a), we see the convergence when taking $n_1 = n_2 = n_3 = 10$. We see that we need quite a large number of iterations but tests showed us that  using Arnoldi for all parameter points separately needs approximately the same number of iterations. If we take the resulting subspace $\subs{V}$ and we use this as the initial subspace for $n_1 = n_2 = n_3 = 30$, we see that we just need a few iterations} \label{fig:conv_diff_3D}
\end{figure}

\section{Parametric shift-inverted Residual Arnoldi method}
in each iteration of the residual Arnoldi method, we applied the operator $\Ten A$ (see \eqref{eqn:diagAX}) or its tensor counterpart to Ritz vectors contaminated with error. In this way we are able to compute eigenvalues which are extreme in some sense. Users are often interested in the eigenvalue close to a certain shift $\sigma$. For this, the shift-inverted Arnoldi method can be used. For the ease of notation we assume that $\sigma = 0$.

In this case we need to be able to efficiently compute linear systems with $\VM A(\VM \omega)$ for all parameter samples. The inverse of $\VM A(\VM \omega)$ can in general not be written in the form \eqref{eqn:VW_A}.  However if all matrices except one is of rank one, i.e $\VM A(\VM \omega) = f_0(\VM \omega) \VM A_0 + \sum_{i=1}^m f_i(\VM \omega) \VM A_i$ with $\VM A_i = \VM u_i \VM q_i^T, i = 1, \hdots, m$ with $m$ a small number, it can be proven via the Sherman-Morrison-Woodbury formula \cite{Hager1989, Woodbury1950} that for a vector $\VM q$, it holds that

\begin{equation} \label{eqn:Ainv_subs}
    \VM A(\VM \omega)^{-1} \VM q \in \text{span}( \VM A_0^{-1} \VM q, \VM A_0^{-1} \VM u_1, \VM A_0^{-1} \VM u_2, \hdots, \VM A_0^{-1} \VM u_m).
\end{equation}
We see that there is only one basis vector of the subspace that depends on $\VM q$, while other vectors are the same for any $\VM q$ and for all parameter samples. Hence, one could keep all other vectors during all iterations in subspace $\subs{V}$.

It is possible to directly apply the PRA algorithm on $\VM A(\VM \omega)^{-1}$, however, in \cite{LeeStewart2007, Lee2007}, the authors proposed an adaptation of RA. In exact arithmetic, if we apply the RA algorithm on $\VM A^{-1}$, we add the residual vector $\VM r_{A^{-1}} = \VM A^{-1} \VM y - \dfrac{1}{\lambda} \VM y$ with $(\lambda, \VM y)$ a Ritz-pair of $\VM A$. The authors proposed to first calculate the residual $\VM r = \VM A \VM y - \lambda \VM y,$ to solve the system
\begin{equation} \label{eqn:system_A}
    \VM A \VM u = \VM r.
\end{equation}
and to add the resulting $\VM u$ to the subspace. The authors in \cite{LeeStewart2007, Lee2007} proved that $\VM u$ is parallel to $\VM r_{A^{-1}}$. This has the advantage that solving \eqref{eqn:system_A} inexactly, produces an error proportional to $\VM r$.

As the authors of \cite{LeeStewart2007, Lee2007} claim that also here the error proportional to the residual can be as high as $10^{-3}$ without harming the convergence, it is more appropriate to calculate the vector via \eqref{eqn:system_A}. Their method was called the Shift-and-invert Residual Arnoldi algorithm and similarly to the PRA algorithm, we extend this method to the parametric case. Instead of \eqref{eqn:system_A} we only need to calculate
$$\VM A_0 \VM u(\VM \omega) = \VM r(\VM \omega)$$
upto a relative tolerance $\eta_R$, as we keep the other basis vectors in \eqref{eqn:Ainv_subs} in the subspace $\subs{V}$ in all iterations.
The tolerance for the relative error when solving \eqref{eqn:system_A} is denoted with $\eta_R$ and the tolerance for the low-rank approximation of the eigenvector over the parameter, is denoted with $\eta_X$.

\subsection{Numerical experiments} \label{sec:PSIRA_numexp}

The following examples arise from spectrum-based stabilization methods for time-delay systems \cite{Fenzi2017,Michiels2011}. Time-delay stability can be inferred from the spectral abscissa (real part of the rightmost eigenvalue) of an eigenvalue problem, discretizing the infinitesimal generator operator by collocation approach \cite{Breda2015}. Spectrum-based stabilization is based on minimizing the spectral abscissa as a function of controller parameter and favors situations where the spectral abscissa is not smooth, since minima of the spectral abscissa are often characterized by defective rightmost eigenvalues \cite{Fenzi2019}. In this context, our novel approach may lead to a significant contribution in the probabilistic spectrum-based stabilization \cite{Fenzi2017} since it permits to tackle larger scale problem with higher stochastic dimension. The code can be downloaded following \href{http://twr.cs.kuleuven.be/research/software/delay-control/TensorKrylov}{this link}.

\begin{eexample} \label{ex:FbackCtrl1}
The spectral abscissa function, associated to the oscillator with feedback delay system \cite[Example~2]{Fenzi2019},  is approximated by the Infinitesimal Generator approach  \cite{Breda2015} as the spectral abscissa of the following parametric matrix:
\begin{equation}\label{eq:Aex2}
\VM A(\omega_1, \omega_2;K_1, K_2) = \VM A_0(K_1, K_2) -\omega_1^2  \VM A_1 +-2\omega_1\omega_2\VM A_2, \text{ with } \omega_1 \in [0.9, 1.1],\; \omega_2 \in [0.1, 0.2].
\end{equation}
The non-singular matrix $\VM A_0$ and the rank-$1$ matrices $\VM A_1$, $\VM A_2$ have dimension $1002 \times 1002$ and are given by:
\[
\VM A_0(K_1, K_2)=\begin{pmatrix}
0 & 1 & 0 & \cdots & 0 & 0 & 0\\
0 & 0 & 0 & \cdots & 0 & K_1 & K_2\\
  &   &   & \VM I_2 \otimes \VM D_M  &   &     &
\end{pmatrix},\;

\VM A_1=\begin{pmatrix}
0 \\ 1 \\ 0 \\ \vdots \\ 0
\end{pmatrix}
\begin{pmatrix}
1 \\ 0 \\ 0 \\ \vdots \\ 0
\end{pmatrix}^T,\;

\VM A_2=\begin{pmatrix}
0 \\ 1 \\ 0 \\ \vdots \\ 0
\end{pmatrix}
\begin{pmatrix}
0 \\ 1 \\ 0 \\ \vdots \\ 0
\end{pmatrix}^T
\]
where $\VM D_M\in\mathbb{R}^{M\times M}$ is the differentiation matrix \cite[Chapter~6]{Trefethen2000} obtained with $M=500$ Chebyshev nodes as collocation points. The controller parameters $K_1$ and $K_2$ are selected so that the spectral abscissa function presents different behaviors on the parametric domain  $(\omega_1,\omega_2)\in[0.9, 1.1]\times [0.1, 0.2]$ as illustrated in Figures \ref{fig:FbackCtrl1_eigw_conv} (a) and (b): for $(K_1, K_2) = [0.2, 0.2]$  the spectral abscissa smoothly behaves, while for $(K_1, K_2) = [0.6179, -0.0072]$  the spectral abscissa presents non-Lipschitz continuous points, due to multiple non-semi-simple rightmost eigenvalues. We refer to \cite{Fenzi2019} for further information on the example.

We use $\eta_R = \dfrac{ 10^{-3}}{2+10^{-3}}$ fixed and $\eta^{(k)}_X = 10^{-3} \norm{ \hat{\VM R}^{(k)}}_F/\sqrt{n_0}$ relative to the residual as explained in Section \ref{sect:Error_an}. The parameter intervals are both discretised into $30$ points, so we have a total of $n_0 = 900$ points for which the spectral abscissa needs to be found. The spectral abscissa is plotted in Figure \ref{fig:FbackCtrl1_eigw_conv} (a)-(b). The convergence results are shown in Figure \ref{fig:FbackCtrl1_eigw_conv} (c)-(d).
\end{eexample}

\begin{figure}

\setlength{\figW}{7cm}
 \begin{subfigure}[b]{0.5\textwidth}

  \input{Afbeeldingen/FbackCtrl_1_eigenw.tex}
  \caption{}
\end{subfigure}
 \begin{subfigure}[b]{0.5\textwidth}

  \input{Afbeeldingen/FbackCtrl_3_eigenw.tex}
  \caption{}
\end{subfigure} \\
\setlength{\figW}{7cm}
 \begin{subfigure}[b]{0.5\textwidth}

  \input{Afbeeldingen/FbackCtrl_1_meanresidu.tex}
  \caption{}
\end{subfigure}
 \begin{subfigure}[b]{0.5\textwidth}

  \input{Afbeeldingen/FbackCtrl_3_meanresidu.tex}
  \caption{}
\end{subfigure}
\caption{The spectral abscissa over the parameters for the two different cases in Example \ref{ex:FbackCtrl1}. We see in (a) that the spectral abscissa varies smoothly over the parameters, in (b) this is not the case anymore. The convergence results for the two different cases in Example \ref{ex:FbackCtrl1} are plotted in (c) and (d). We see that in the two cases the algorithm converges within eleven iterations.} \label{fig:FbackCtrl1_eigw_conv}
\end{figure}

\begin{eexample}  \label{ex:Heating}

	We analyze the spectral abscissa function associated to the heat transfer set up with two static feedback controllers $\VM K$, designed by probabilistic stabilization method considering the temperature on the left cooler and both cooler affected by uncertainty \cite[Section~5.3]{Fenzi2017}.
	The spectral abscissa function is determined by the eigenvalues of the following $1030 \times 1030$ parametric matrix
	\[\VM A(\omega_1, \omega_2; \VM K) = \VM A_0(\VM K) +  \dfrac{1-\omega_1}{\omega_1} \VM A_1 + \dfrac{1-\omega_2}{\omega_2} \VM A_2, \text{ where }\omega_1\in[13.5, 16.5],\, \omega_2\in[15.3, 18.7].\]

	The two static controllers $\VM K$, designed in order to minimize the mean of the spectral abscissa, modelling the uncertainty on the temperatures by independent uniform random variables \cite[Section~5.3]{Fenzi2017},  lead to different behaviors of the spectral abscissa.  In Figure \ref{fig:Heating_eigw_conv} (a) resp. (b), the spectral abscissa for $ \VM K = [0.5105, -9.1810 \times 10^{-2}]$ resp. $\VM K=[0.45043,-0.1631]$ is used. The parameter intervals are both discretised into $30$ points, so we have a total of $n_0 = 900$ parameter values for which the spectral abscissa needs to be found. We use $\eta_R = \dfrac{ 10^{-3}}{2+10^{-3}}$ fixed and $\eta^{(k)}_X = 10^{-3} \norm{ \hat{\VM R}^{(k)}}_F/\sqrt{n_0}$ relative to the residual as explained in Section \ref{sect:Error_an}.  The convergence results are shown in Figure \ref{fig:Heating_eigw_conv} (c)-(d).
\end{eexample}

\begin{figure}

\setlength{\figW}{7cm}
 \begin{subfigure}[b]{0.5\textwidth}

  \input{Afbeeldingen/Heating_2_eigenw.tex}
  \caption{}
\end{subfigure}
 \begin{subfigure}[b]{0.5\textwidth}

  \input{Afbeeldingen/Heating_4_eigenw.tex}
  \caption{}
\end{subfigure} \\
\setlength{\figW}{7cm}
 \begin{subfigure}[b]{0.5\textwidth}

  \input{Afbeeldingen/Heating_2_meanresidu_fix.tex}
  \caption{}
\end{subfigure}
\begin{subfigure}[b]{0.5\textwidth}

  \input{Afbeeldingen/Heating_4_meanresidu_fix.tex}
  \caption{}
\end{subfigure}
\caption{In (a)-(b), The spectral abscissa over the parameters for the two different cases in Example \ref{ex:Heating}. The convergence results for the two different cases in Example \ref{ex:Heating} can be seen in (c)-(d).} \label{fig:Heating_eigw_conv}
\end{figure}

We see in Figures \ref{fig:FbackCtrl1_eigw_conv} (c)-(d) and \ref{fig:Heating_eigw_conv} (c)-(d) that in both examples only a few iterations were needed. In Figure \ref{fig:FbackCtrl1_eigw_conv} (b) and \ref{fig:Heating_eigw_conv} (a)-(b) it can be observed that the spectral abscissa over the parameter is not differentiable as in some points (we only count eigenvalues with imaginary part equal or greater than zero) there are multiple eigenvalues for which the real part collides. We need to be careful when some eigenvalues of the large eigenvalues are close to each other. In our experience, subspace methods for calculating eigenvalues of non-normal matrices, can fail to converge, if the searched eigenvalue is close to other eigenvalues. Therefore, it is possible that due to the projection an eigenvalue is selected which is not the one we look for. The technique that we apply, is to save for every parameter sample in iteration $k$ five eigenvectors of the projected problem that best satisfies our criterium. In the folowing iteration $k + 1$ only an eigenvalue with associated eigenvector is selected if the angle between this eigenvector and one of the five selected eigenvectors in iteration $k$ is below a prescribed tolerance.

\section{Conclusion}

In this paper, we proposed a method to calculate eigenvalues of a parametrized matrix $\VM A(\omega)$ for a number of parameter values at the same time. We started with introducing the Residual Arnoldi algorithm, which is a version of the Arnoldi algorithm that is robust against introduced errors. We applied the Residual Arnoldi Algorithm on a large number of sample points at the same time. To keep the computations feasible, we made low-rank approximations. On the one hand, we made a low-rank approximation of the searched eigenvector of the reduced eigenvalue problem over the parameters with a relative tolerance of $\eta_X$, and on the other, we made a low-rank approximation on the residual over the parameters with a relative tolerance of $\eta_R$. We saw that the total error relative to the residual needs to be smaller than $10^{-3}$, from which we deduce that $\eta_X$ needs to be chosen proportional to the residual and that $\eta_R$ can be chosen as high as $\dfrac{ 10^{-3}}{2+10^{-3}}$, see Section \ref{sect:conv_res_arn}.

We illustrated with numerical experiments that a higher $\eta_R$ even improves performance, since a higher tolerance leads to a slower increase of the subspace dimension and consecutely less restarts. At the end of the paper, we made the extension to mutltiple parameters and to shift-inverted Arnoldi for a specific range of problems. A possible extension is to make a version of the parametric shift-inverted Arnoldi method when more than one matrix has rank higher than one. The difficulty in this problem, is that we need to be able to solve parametrized linear systems efficiently. Another possibility for further research is the calculation of statistical moments of certain eigenvalues over the parameter.

\section*{Acknowledgements}
The authors thank the referees for their helpful remarks. We also thank Nick Vannieuwenhoven for some fruitful discussions about the used tensor format and we thank Luca Fenzi for the examples used in Section \ref{sec:PSIRA_numexp}. This work was supported by the project KU Leuven Research Council grant C14/17/072 and by the project G0A5317N of the Research Foundation-Flanders (FWO - Vlaanderen). This work does not have any conflicts of interest.

\bibliographystyle{elsarticle-harv}
\bibliography{ms}

\end{document}

config: arxiv


TITLE: Find the Closest Vector to a Given Vector from a Given Set Defined by Set of Multiplications
QUESTION [3 upvotes]: Let $x=\left(x_t\right)_{t=1}^n$ be a vector such that
$$
x_t = \prod_{i=1}^t u_i, \tag{1}
$$
where each parameters $u_i$ can take any of two value
$$
u_i \in \left\{a,b \right\} = \left\{ 1.3, 0.8 \right\}.
$$
Each vector $x$ represents a path in a binomial tree.
Given an arbitrary vector $y=\left(y_t\right)_{t=1}^n$,
where each $y_t$ is in $\mathbb R^+$, how to solve the problem
$$
\min_{x \in X} \lVert y - x \rVert^2 = \sum_{t=1}^n \left(y_t - x_t \right)^2
$$
where $X$ is the set of all possible vector of the form (1).
So, given a realization of the stock, find the closest binomial path.
The only way I see to solve this problem is to do Gram-Schmidt on $X$,
and then project $y$ on the orthonormal basis found.
But, my $n$ is large, say $100$.
This induces
$$\mathrm{cardinality} X = 2^{100},$$
so I don't think Gram-Schmidt is viable.
I would use least-square or Newton-Raphson, but the variable are not continuous.

REPLY [0 votes]: I would got with one of the approaches (They are not closed form solutions but might be good enough):

Try estimating $ {u}_{i} $ by looking at the vector $ \left[ {y}_{1}, {y}_{2} / {y}_{1}, {y}_{3} / {y}_{2}, \ldots, {y}_{n} / {y}_{n - 1} \right] $ then taking the nearest value out of $ \left\{ 0.8, 1.3 \right\} $.
Apply the $ \ln \left( \cdot \right) $ on each element of $ \boldsymbol{y} $. Then each element becomes sum of $ {u}_{i} $ which means you can solve a problem of the form $ A u = \ln \left( y \right) $ (Pay attention to $ A $ as it is only 1 and 0 in a convolution like form). From the estimated $ \boldsymbol{u} $ you may build $ \boldsymbol{x} $.


TITLE: Proving some identities making use of properties of complex numbers
QUESTION [2 upvotes]: Problem: Show that for real $\phi$ we have the following identities: \begin{align*} & a) \ \sin(3 \phi) = 3 \cos^2 (\phi) \sin(\phi) - \sin^3 (\phi), \\ & b) \  1 + e^{i \phi} + e^{2 i \phi} + \ldots + e^{i n \phi} = e^{\frac{i}{2} n \phi} \frac{\sin (\frac{1}{2}(n+ 1) \phi)}{\sin \big(\frac{1}{2} \phi\big)}. \end{align*}
Attempt: For $(a)$, I did the following. We have that $$ e^{3 i \phi} = \cos (3 \phi) + i \sin (3 \phi) $$ and so \begin{align*} (e^{i \phi})^3 &= ( \ cos \phi + i \sin \phi)^3 \\ &= \cos^3 \phi + 3i \cos^2 \phi \sin \phi - 3 i \cos \phi \sin^2 \phi - i \sin^3 \phi \\ &= \underbrace{(\cos^3 \phi - 3 \cos \phi \sin^2 \phi)}_{\cos 3 \phi} + i \underbrace{(3 \cos^2 \phi \sin \phi - \sin^3 \phi)}_{\sin 3 \phi}
\end{align*} With this argument I proved (a). But I don't know how to prove (b), help would be appreciated!

REPLY [1 votes]: First apply the formula for the sum of a geometric progression(assume $\phi\ne 0 )$:
$$1 + e^{i \phi} + e^{2 i \phi} + \ldots + e^{i n \phi}=\frac{e^{i(n+1)\phi}-1}{e^{i\phi}-1}$$
Break the numerator and denominator like this:$$\frac{e^{i(n+1)\phi}-1}{e^{i\phi}-1}=\frac{e^{i\phi(n+1)/2}}{e^{i\phi /2}}\cdot\frac{e^{i\phi(n+1)/2}-e^{-i\phi(n+1)/2}}{e^{i\phi /2}-e^{-i\phi /2}}$$
Now I think you can simplify the thing to get the form you want.


TITLE: Sturm-Liouville system variable bounded in interval
QUESTION [0 upvotes]: I am solving the following Sturm-Liouville system problem
$1.$ $X''(x) + \lambda X(x)=0, 0<x<2 $
$ 2.$ $X'(0)=0 , X(2) =0$
So I have considered the cases when $ \lambda <0  , \lambda=0$ .Now I am at the 3rd case when $\lambda >0$ then
the solutions is $X(x)= c_1 \cos(\lambda x)+c_2 \sin(\lambda x )$ .
From $X'(0) =0$ we get that $c_2=0$ and from $X(2)=0 \implies c_1=0 $ is one solution,the other is $\cos(\lambda x)=0$.
Now when $\cos(\lambda x)=0$,how should I proceed with the solution considering $0<x<2$ ?

REPLY [1 votes]: For $\lambda \gt 0$, general solution of the given Sturm-Liouville problem is $$X(x)= c_1 \cos(\sqrt \lambda x)+c_2 \sin(\sqrt \lambda x )$$
Now $X'(0)=0 , X(2) =0 \quad$ gives $c_2=0$ and $c_1\cos(2\sqrt \lambda )=0$
Now if $c_1=0$, then we get the trivial solution, so we take $c_1 \neq 0 $ and
$\cos(2\sqrt \lambda )=0\implies 2\sqrt \lambda=(2n+1)\frac{\pi}{2}\quad$ where $n=0,\pm1,\pm2,. . . $
So $\sqrt \lambda = (2n+1)\frac{\pi}{4} \implies \lambda= (2n+1)^2 \frac{\pi^2}{16}\quad$  where $n=0,\pm1,\pm2,. . . $
Hence the general solution is $$X(x)= c_1 \cos((2n+1)\frac{\pi}{4} x),\qquad \text{where} \quad n=0,\pm1,\pm2,. . . $$


TITLE: How to explain the following definition of $\lim\sup s_n$ intuitively?
QUESTION [1 upvotes]: How to explain the following definition of $\lim\sup s_n$ intuitively?
let $\lim\sup s_n$=$l$
The definition is ($\forall t<l,\forall N,\exists n>N, s_n>t$) and ($\forall t>l,\exists N,\forall n>N, s_n<t$)
For the first one, if we have a increasing sequence from very negative to $0$, how is it still true. i.e. $\{-1, -0.5, ... , -0.00001,...\}$, here $\lim\sup s_n=0$. let $t=-1,$ so $t<0$, but not all terms are greater than $t$.
For the second one, I am also confused. I think It would be great if someone could show some valid examples.
Also when we prove the such limit exists, is it necessary to show both characteristics?

REPLY [4 votes]: Informally speaking the $\limsup s_n$ is the limit of the sequence defined from the supremum of the tails of the sequence $(s_n)$.
Here the intuitive key is the concept of tail of a sequence. A tail, as it descriptive name want to define, is a sequence resulting of discard the first $m$ elements of the sequence (where $m$ is finite, consequently a tail ever have infinite elements).
We will represent the sequence visualized as
$$\overbrace{s_0,s_1,\underbrace{s_2,\ldots,s_k,\ldots}_{\text{a tail of }(s_n)}}^{\text{the sequence }(s_n)}$$
Hence we can define the sequence of tails from the original sequence $(s_n)$ as
$$t_0=s_0,s_1,\ldots$$
$$t_1=s_1,s_2,\ldots$$
$$t_k=s_k,s_{k+1},\ldots$$
Then $\limsup s_n=\lim(\sup t_n)$.

REPLY [0 votes]: Your first criterion means that for all $t<l$, there exists a subsequence $s_{n_k}$ of $s_n$ such that all members of the subsequence are greater than $t$.
Your second criterion means that for every $t>l$ there exists a point in the sequence $s_n$ from which all sequence members are smaller than $t$.


TITLE: Preparing for exam and don't understand how to calculate this integral via integrating by parts twice
QUESTION [2 upvotes]: It seems to be a pretty simple integral but I am blanking out on how to find this integral:
$$\int_0^1 e^{-x}\cos(n\pi x) \;dx$$
I tried Mathematica, but it seems to use some formula instead of integrating by parts.
Thanks for your help!

REPLY [1 votes]: Here is another method that applies if you know some complex numbers.
Recall
$$e^{i \theta} = \cos \theta + i \sin \theta$$
In this case we have a $\cos (n \pi x)$, which is the real part of $e^{i n \pi x}$
$$\mathscr{R}(e^{in \pi x}) = \cos(n \pi x)$$
Then we have $$\int_{0}^1 e^{-x} \cos (n \pi x)dx = \mathscr{R} \int_{0}^1 e^{-x}\cdot e^{i n \pi x}$$


TITLE: Theorem 1.17 in Rudin's RCA
QUESTION [2 upvotes]: Here is theorem 1.17 in Rudin's Real and Complex Analysis

Let $f : X \to [0,\infty]$ be measurable. There exist simple measurable functions $s_n$ on $X$ such that (a) $0 \le s_1 \le s_2 \le ... \le f$ and (b) $s_n(x) \to f(x)$ as $n \to \infty$ for every $x \in X$.

I am working through both Rudin's proof and this proof of the above theorem. I just want to go over a few things to make sure I understand everything.
First, I think I can see why $\varphi_n$ is a Borel function: both $[k2^{-n},(k+1)2^{-n})$ and $[n,\infty]$ are Borel sets, so the characteristic function on each is a Borel function; finally, the sum of Borel function is a Borel function, making $\varphi_n(t)$ a Borel function. Second, I think I can see that $\varphi_n(t) \to t$. Given $t \ge 0$ and $n \in \Bbb{N}$, there exists a nonnegative integer $k_n$ such that $k_n < t 2^{n} < k_n+1$ or $k_n 2^{-n} < t < k_n 2^{-n} + 2^{-n}$ and therefore $0 < t - k_n 2^{-n} < 2^{-n}$. But $\varphi_n(t) = k_n 2^{-n}$, so that taking the limit $n \to \infty$, we get $\varphi_n(t) \to t$. Doing the the same for $f(x) \ge 0$, we can show that $\varphi_n(f(x)) \to f(x)$ for every $x \in X$. Hopefully someone having greater familiarity with the subject and Rudin's book RCA will correct me if I am wrong.
However, I am at the moment having trouble seeing why $\{\varphi_n\}$ is a monotonically increasing sequence.

REPLY [2 votes]: For an interval $I:=\left[k2^{-n},(k+1)2^{-n}\right)$, we have
\begin{align*}
I&=\left[2k\cdot 2^{-(n+1)},(2k+2)\cdot 2^{-(n+1)}\right)\\
&=\left[2k\cdot 2^{-(n+1)},(2k+1)\cdot 2^{-(n+1)}\right)\bigcup\left[(2k+1)\cdot 2^{-(n+1)},(2k+2)\cdot 2^{-(n+1)}\right),
\end{align*}
you may now see that for a typical interval of $\phi_{n+1}$, it belongs to either the left part or the right part of a typical interval of $\phi_{n}$, and if you compare the values that $\phi_{n+1}$ and $\phi_{n}$ taken over there, they are $l2^{-(n+1)}$ and $k2^{-n}$ respectively, if it were belonging to the left part, then $l=2k$, if it were right part, then $l=2k+1$, in either case, $\phi_{n+1}(x)\geq\phi_{n}(x)$ over there.


TITLE: Natural deduction technique/question
QUESTION [2 upvotes]: Current Knowledge:

Know off by heart the standard rules + variant rules for the inference system (natural deduction).
Know about predicate logic and proportional logic.

Does anyone have a technique in how to approach natural deduction questions, because honestly I can't rack my mind around it. I know practice will help me get better an all, and I have been doing tons of practice. But ever time a new question comes up I just get completely stuck....
If someone can give me some tips/heuristics or something to help me, that would be awesome. Something like what to watch out for if the conclusion is supposed to be this, then approach it in this manner...etc . Basic mathematical language type explanations would be preferable.
I have been having trouble with the following example, as of late. If you could base your answers on this,that would be awesome.
Example:
¬(P∧ ¬Q) , R ⊢ P → (Q ∧ R)
Side question:

Could someone tell me the distinct difference between ⊢ and ⊨ ?
What is meant by soundness and completeness ? if you could link this too both predicate logic and proportional logic, that would be awesome.

Again basic mathematical language type explanations would be preferable.
Thanks in advance.

REPLY [1 votes]: In order to apply Natural Deduction to the example, we have to start from the sought conclusion : $P \to (Q \land R)$.
It is a conditional, so we have to assume $P$ and try to derive $Q \land R$.
We have $R$ but not $Q$; thus, we need a way to derive it: if we assume $\lnot Q$ and are able to derive a contradiction, we can use double negation elimination to get $Q$.
Thus, we need a further assumption : $P$.

1) $¬(P∧ ¬Q)$ --- 1st premise
2) $R$ --- 2nd premise

3) $P$ --- "temporary" assumption, starting sub-proof [a]

4) $\lnot Q$ --- "temporary" assumption, starting sub-proof [b]
5) $P \land \lnot Q$ --- $\land$-intro, from 3) and 4)
6) $\bot$ --- from 1) and 5), by $\lnot$-elim
7) $Q$ --- from 4) and 6) by double negation elimination, closing sub-proof [b]

8) $Q \land R$ --- $\land$-intro
9) $P \to (Q \land R)$ --- $\to$-intro, from 3) and 8), closing sub-proof [a].

The above derivation amounts to showing that:

$¬(P∧ ¬Q), R \vdash P \to (Q \land R)$

i.e. that $P \to (Q \land R)$ is derivable (in the calculus of Natural Deduction) from the set of premises: $\{ ¬(P∧ ¬Q), R \}$.
In general, to say that:

$\Gamma \vDash \varphi$,

where $\Gamma$ is a set of formulas and $\varphi$ a formula, means that the conclusion $\varphi$ is a logical consequence of the premises.
The calculus is (stongly) complete when :

if $\Gamma \vDash \varphi$, then $\Gamma \vdash \varphi$.

The calculus is sound when :

if $\Gamma \vdash \varphi$, then $\Gamma \vDash \varphi$.


TITLE: Conditional Distribution of a Brownian Motion
QUESTION [0 upvotes]: For a fixed $t>0$, how to compute the conditional distribution of a Brownian motion at time $t$, $B(t)$, given that $B(2t)=x$?
My reasoning is that this is just equivalent to the distribution of $B(t)$ given $B(0)=x$. So I am guessing that $\mathcal{N}(x,t)$ is the answer, but I am not sure. What looks really odd in my reasoning is the reversed time.
Any thoughts on this? Any help is appreciated.

REPLY [2 votes]: The joint distribution of $\begin{bmatrix} B(s) \\ B(t)\end{bmatrix} $ is a joint Gaussian distribution with mean $\begin{bmatrix} 0 \\ 0 \end{bmatrix}$, and covariance matrix $\begin{bmatrix} s & s \wedge t \\ t \wedge s & t \end{bmatrix}$.
Note that the conditional distribution is still Gaussian with some well known mean and covariance (which you can calculate via Proposition 3.9 in these notes; Eq. 3.10 has a typo it should read $\hat{E}[X|Y=y] = E[X] + \ldots$)


TITLE: How does this integration pass?
QUESTION [2 upvotes]: Can someone explain the passage in the red-box for me?

I am getting
\begin{align}
& -\frac{1}{2\pi}\int_{0}^{2\pi} \log r  \;d\theta  - \frac{1}{2\pi} \int \log \| e^{i\theta_1} - e^{i\theta}\|  \;d\theta + \log r \\
= {} & -\frac{1}{2\pi} \int \log \| e^{i\theta_1} - e^{i\theta}\|  \\
\stackrel{?}{=} {} & - \frac{1}{2\pi} \int \log \| 1 - e^{i\theta}\| \, d\theta
\end{align}
It briefly mentions "periodicity" of $e^{i\theta}$.
But I don't see how that fits in here...

REPLY [1 votes]: Note that
\begin{align}
\log |e^{i\theta_1} - e^{i\theta}|
& = \log |e^{i\theta_1}||1 - e^{i(\theta - \theta_1)}|\\
& = \log |e^{i\theta_1}| + \log |1 - e^{i(\theta - \theta_1)}|\\
& = \log 1 + \log |1 - e^{i(\theta - \theta_1)}|\\
& = 0 + \log |1 - e^{i(\theta - \theta_1)}|\\
& = \log |1 - e^{i(\theta - \theta_1)}|
\end{align}
Hence
$$\int_0^{2\pi} \log |e^{i\theta_1} - e^{i\theta}| \, \text{d}\theta
 = \int_0^{2\pi} \log |1 - e^{i(\theta - \theta_1)}| \, \text{d}\theta$$
Substituting $\phi = \theta - \theta_1$ and by periodicity and another substitutionin the third equal sign, we have
\begin{align}
\int_0^{2\pi} \log |1 - e^{i(\theta - \theta_1)}| \, \text{d}\theta
& = \int_{-\theta_1}^{2\pi - \theta_1} \log |1 - e^{i\phi}| \, \text{d}\phi\\
& = \int_{-\theta_1}^{0} + \int_{0}^{2\pi-\theta_1} \log |1 - e^{i\phi}| \, \text{d}\phi\\
& = \int_{2\pi-\theta_1}^{2\pi} + \int_{0}^{2\pi-\theta_1} \log |1 - e^{i\phi}| \, \text{d}\phi\\
& = \int_{0}^{2\pi} \log |1 - e^{i\phi}| \, \text{d}\phi\\
\end{align}


TITLE: Maximal induced cycles on $n$-clique graphs
QUESTION [0 upvotes]: For any set $X$ we set $[X]^2 = \big\{\{a,b\}: a, b\in X\text{ and } a\neq b\big\}$.
We say a simple undirected graph $G=(V,E)$ is an $n$-clique graph if there are $S_1,\ldots,S_n\subseteq V$ such that

$|S_k| = n$ for all $k=1,\ldots, n$;
$V = \bigcup_{k=1}^n S_k$;
$i\neq k \in \{1,\ldots,n\}$ implies $|S_i \cap S_k| = 1$.
$E = \bigcup_{k=1}^n [S_k]^2$ (that is, all the $S_k$ are complete, and there are no edges between different $S_k$.)

Let $c(n)$ be the maximum length of an induced cycle that any $n$-clique graph $G$ can have. Is there an explicit formula for $c(n)$, and if not, what is $\lim_{n\to\infty}\frac{c(n)}{n}$?

REPLY [2 votes]: I am not sure that I understand all conditions correctly, since the question looks too straightforward, but if yes, then there are no cycles for $n\leq 2$ and for $n>2$ I claim that $c(n)=n$.
At first, $c(n)\leq n$. Indeed, no two edges of our cycle $x_1\dots x_{c(n)}x_1$ may belong to the same clique (if $c(n)\geq n+1\geq 4$), else the cycle is not induced.
Example of induced cycle $1\dots n$ of length $n$: $S_i$ contains vertices $i,i+1$ (modulo $n$, of course), and $n-2$ vertices $v_{i,j}$, $j\in \{1,\dots,n\}\setminus \{i+1,i-1\}$.


TITLE: Probability that a Multivariate Normal RV lies within a Spherical Region of Radius R
QUESTION [1 upvotes]: I am currently using different procedures to estimate the probability that a $D$-dimensional Gaussian random variable with mean $\mu$ and covariance $\Sigma$ lies within a sphere of radius $R$ that centered about the origin.  That is, I am estimating $P(|| X ||_2  < R)$  where  $X \sim N(\mu, \Sigma)$  and $X \in \mathbb{R}^D$.
I am wondering whether there is a way to obtain the exact value of this probability analytically (i.e. without using numerical integration or Monte Carlo)? I currently have two basic approaches to follow:
Approach 1
Find a way to analytically evaluate the integral:
$\int_{x \in S} (2\pi)^{-\frac{D}{2}}|\Sigma|^{-\frac{1}{2}} \exp(-\frac{1}{2} (x-\mu)^T \Sigma^{-1} (x-\mu) dx $
over the spherical region:
$S = \{||x|| < R \} = \{x^Tx < R^2\}$
Approach 2
Exploit the fact that given $X \sim N(\mu,\Sigma)$, the random variable
$(x-\mu)^T \Sigma^{-1} (x-\mu) \sim \chi^2(D) $
This implies that
$P( (x-\mu)^T \Sigma^{-1} (x-\mu) < R^2 ) = P(\chi^2(D) < R^2)$
which is very simple to calculate... though what I will actually need to evaluate is:
$P( x^T x < R^2) $

REPLY [1 votes]: There is not a closed form solution.  So there is no better way to calculate it than numerical integration.  The Monte Carlo method as I proposed it will approximate the integral and the approximation will be close to exact if you generate lots of samples like 100,000 or 1 million.


TITLE: For which x the matrix above is the covariance matrix?
QUESTION [0 upvotes]: For which x the matrix above is the covariance matrix?
\begin{pmatrix} 1 & x & x \\ x & 1 & x \\ x & x & 1 \\ \end{pmatrix}
I only understand it must be positively-definite, but how to prove that all these matrices are realizable with some vector $(\xi_1, \xi_2, \xi_3)$?

REPLY [0 votes]: If it is positive definite and symmetric then it can be diagonalised by a rotation in 3-space so it is sufficient to prove that all diagonal matrices with positive eigenvalues are realizable. But that can easily be achieved using 3 independent variables with the appropriate variances.
The solution to the original problem is then obtained by rotating the 3 variables back, which gives linear combinations of them.
Thus it is sufficient to find all $x$ for which the matrix is positive definite.


TITLE: Understanding the Concept of Monodromy; case of Lefschetz Fibrations.
QUESTION [1 upvotes]: My question is on the concept of monodromy around  critical points in a Lefschetz fibration $p: M^4 \rightarrow S^2$  (and monodromy in general), where $M^4$ is a 4-manifold and $S^2$ is the 2-sphere.
It is  known that the monodromy about any  critical point in a Lefschetz fibration (there are finitely-many critical points) is described by a Dehn twist $D$. Is the following an accurate description of the concept of monodromy by D about a critical point x in this type of fibration:
Let x be a critical  point in the base $S^2$. Let  $x':= p^{-1}(x)$ be the lift of x by the fibration map p ,x' lives in the 4-manifold M. Then, if we draw a small loop $B(x,r)$ (small-enough so that $B(x,r)$ does not intersect any other critical point in the base; critical points for a Lefschetz fibration are isolated ) in the base, centered at x, and we wind around the loop once, say from $a=0$ to $b=2\pi$, then $a':= p^{-1}(a)$ is related to $b':=p^{-1}(b)$ by the Dehn twist $D$, i.e., is it the case that $a'=D(b')$?
Thanks for any help, references, etc

REPLY [0 votes]: This is how I understand monodromy: consider the small ball $B(x,r)$ within a  trivialized neighborhood that you described. Then this ball, since it is  assumed to be inside of a trivialized neighborhood, is an $S^1-$ bundle $S^1 \times F$ , where $F$ is the fiber. Now, if the monodromy is trivial, you do get a trivialized $S^1 \times F$. Basically, as you wind around the boundary $S^1$ of the ball $B(x,r)$ , each point lifts to a fiber $F$ , but, as you complete a loop around $S^1$ , you end up gluing fibers $F$ to each other when $0 \rightarrow 2\pi$ . This gluing is done by an element of the mapping class group MCG(F) of the fiber. The choice of map, up to isotopy, is the monodromy.
In other words, the monodromy dictates that you glue the lift $S^1 \times F$ along a Dehn twist D , so that you end up with the lifted bundle $S^1 \times F /$~ , where '~' is the relationship $(x,0)$~$(D(x),1)$ , where $D(x)$ is the Dehn twist associated to the monodromy. So you end up gluing your fiber along the monodromy map $D$.


TITLE: What does pc(eV) = 300 Br(G cm) mean?
QUESTION [1 upvotes]: In the January 1983 issue of the American Journal of Physics, Edward Purcell published a one-page "Round-Number Handbook of Physics" (see e.g. http://lemeshko.blogspot.com/2008/06/one-page-round-number-handbook-of.html).  In the section on unit conversions, one line reads: "pc(eV) = 300 Br(G cm)".  Does anyone have any idea what this means?  I note that 1 eV/(G cm) seems to be closely related to the elementary charge, but I still can't quite figure out what this line is about.

REPLY [0 votes]: It's the cyclotron radius $r$ for a particle with momentum $p$ in a magnetic field $B$,
$$
|q|Br = p_\perp,
$$
for the common case where the particle has one quantum of charge.
The numerical conversion constant 300 has the same first significant figure as the speed of light, $c = \rm 3\times10^8\,m/s$. The "rule of thumb" multiplies both sides by $c$ so that both have energy units --- the electron-volt is a convenient unit of energy for subatomic particles. The numerical factor ends up as $3\times 10^2$ because of the difference between gauss-centimeters and the SI unit, tesla-meters.