# On a special class of Schwartz functions

Milutin Obradović  Department of Mathematics, Faculty of Civil Engineering,
University of Belgrade, Bulevar Kralja Aleksandra 73, 11000, Belgrade, Serbia.
<EMAIL_ADDRESS>and  Nikola Tuneski  Department of Mathematics and
Informatics, Faculty of Mechanical Engineering, Ss. Cyril and Methodius
University in Skopje, Karpoš II b.b., 1000 Skopje, Republic of North
Macedonia<EMAIL_ADDRESS>
###### Abstract.

In this paper we study functions
$\omega(z)=c_{1}z+c_{2}z^{2}+c_{3}z^{3}+\cdots$ analytic in the open unit disk
${\mathbb{D}}$ and such that $|\omega^{\prime}(z)|\leq 1$ for all
$z\in{\mathbb{D}}$. For these functions we give estimates (sometimes sharp)
for the following moduli: $|c_{3}-c_{1}c_{2}|$, $|c_{1}c_{3}-c_{2}^{2}|$, and
$|c_{4}-c_{2}^{2}|$.

###### Key words and phrases:

Schwartz functions, coefficient, estimate

###### 2020 Mathematics Subject Classification:

30C45, 30C50

##  1\. Introduction and definitions

For a function $\omega$, analytic in the open unit disk
${\mathbb{D}}=\\{z:|z|<1\\}$ and of the form

(1)
$\omega(z)=c_{1}z+c_{2}z^{2}+c_{3}z^{3}+\cdots,\qquad(c_{1},c_{2},\ldots\in{\mathcal{C}})$

we say that is Schwartz function if $|\omega(z)|<1$, $z\in{\mathbb{D}}$. We
denote by $\mathcal{B}_{0}$ the class of all such functions.

In his paper [4], Zaprawa gave many different inequalities for the
coefficients $c_{1},c_{2},\ldots$ for the functions of the class
$\mathcal{B}_{0}$.

In this paper we study the class of functions $\mathcal{B}^{\prime}_{0}$ of
type (1) such that $|\omega^{\prime}(z)|\leq 1$ for all $z\in{\mathbb{D}}$.
Since

(2) $z\omega^{\prime}(z)=c_{1}z+2c_{2}z^{2}+3c_{3}z^{3}\cdots,$

and $|z\omega^{\prime}(z)|=|z|\cdot|\omega^{\prime}(z)|\leq|z|<1$,
$z\in{\mathbb{D}}$, it means that $z\omega^{\prime}(z)$ belongs to
$\mathcal{B}_{0}$. Also, since $\omega(z)=\int_{0}^{z}\omega^{\prime}(z)dz$,
then $|\omega(z)|\leq\int_{0}^{|z|}|\omega^{\prime}(z)|\,d|z|\leq|z|<1$ for
all $z\in{\mathbb{D}}$, i.e., $\omega\in\mathcal{B}_{0}$. So,
$|\omega^{\prime}(z)|\leq 1$, $z\in{\mathbb{D}}$, is a sufficient condition
for $\omega\in\mathcal{B}_{0}$, i.e., $\mathcal{B}^{\prime}_{0}$ is subclass
of the class $\mathcal{B}_{0}$.

For the functions from $\mathcal{B}^{\prime}_{0}$ we try to find properties
for the coefficients $c_{1},c_{2},c_{3},\ldots$ that correspond to the
properties for the functions from $\mathcal{B}_{0}$.

For our considerations we will need the next lemma originating from [1].

######  Lemma 1.

Let $\omega\in\mathcal{B}_{0}$ is given by (1). Then

(3) $\begin{array}[]{c}|c_{1}|\leq 1,\qquad|c_{2}|\leq
1-|c_{1}|^{2},\\\\[11.38109pt] |c_{3}|\leq
1-|c_{1}|^{2}-\frac{|c_{2}|^{2}}{1+|c_{1}|},\\\\[11.38109pt] |c_{4}|\leq
1-|c_{1}|^{2}-|c_{2}|^{2}.\end{array}$

We showed that when $\omega$ given by (1) is in $\mathcal{B}_{0}$, then
$z\omega^{\prime}(z)$ is in $\mathcal{B}^{\prime}_{0}$. Thus, Lemma 1,
together with (2), directly brings

######  Lemma 2.

Let $\omega\in\mathcal{B}^{\prime}_{0}$ is given by (1). Then

(4) $\begin{split}|c_{1}|&\leq
1,\qquad|c_{2}|\leq\frac{1}{2}\left(1-|c_{1}|^{2}\right),\\\\[11.38109pt]
|c_{3}|&\leq\frac{1}{3}\left(1-|c_{1}|^{2}-\frac{4|c_{2}|^{2}}{1+|c_{1}|}\right),\\\\[11.38109pt]
|c_{4}|&\leq\frac{1}{4}\left(1-|c_{1}|^{2}-4|c_{2}|^{2}\right).\end{split}$

##  2\. Main results

We begin with partly sharp estimate of the modulus $|c_{3}-c_{1}c_{2}|$ for
functions from $\mathcal{B}^{\prime}_{0}$ with expansion (1).

######  Theorem 1.

If $\omega\in\mathcal{B}^{\prime}_{0}$ is of form (1), then

(5)
$|c_{3}-c_{1}c_{2}|\leq\left\\{\begin{array}[]{cc}\frac{1}{48}(1+|c_{1}|)\left[9|c_{1}|^{2}-16|c_{1}|+16\right],&0\leq|c_{1}|\leq\frac{4}{7}\\\\[5.69054pt]
\frac{5}{6}|c_{1}|(1-|c_{1}|^{2}),&\frac{4}{7}\leq|c_{1}|\leq
1\end{array}\right..$

The estimate is sharp for $|c_{1}|=0$ and for $\frac{4}{7}\leq|c_{1}|\leq 1$.

###### Proof.

For $\omega\in\mathcal{B}^{\prime}_{0}$ and $\omega$ given by (1) we apply the
inequalities (4):

$\begin{split}|c_{3}-c_{1}c_{2}|&\leq|c_{3}|+|c_{1}||c_{2}|\leq\frac{1}{3}\left(1-|c_{1}|^{2}-\frac{4|c_{2}|^{2}}{1+|c_{1}|}\right)+|c_{2}||c_{2}|\\\
&=-\frac{4}{3(1+|c_{1}|)}|c_{2}|^{2}+|c_{1}||c_{2}|+\frac{1}{3}(1-|c_{1}|^{2}).\end{split}$

If we consider the last expression as a function of $|c_{2}|$,
$0\leq|c_{2}|\leq\frac{1}{2}(1-|c_{1}|^{2})$, then we easily obtain the
estimate given by (5), depending on its maximum which in the case
$0\leq|c_{1}|\leq\frac{4}{7}$ is attained for
$|c_{2}|=\frac{3}{8}|c_{1}|(1+|c_{1}|)$ lying in the interval
$\left(0,\frac{1}{2}(1-|c_{1}|^{2})\right)$, and in the case
$\frac{4}{7}\leq|c_{1}|\leq 1$ is attained for
$|c_{2}|=\frac{1}{2}(1-|c_{1}|^{2})$.

For $|c_{1}|=0$ and for $\frac{4}{7}\leq|c_{1}|\leq 1$ the result is sharp
with extremal functions $\omega_{1}(z)=\frac{1}{3}z^{3}$ and

$\omega_{2}(z)=\int_{0}^{z}\frac{|c_{1}|+z}{1+|c_{1}|z}\,dz=|c_{1}|z+\frac{1}{2}(1-|c_{1}|^{2})z^{2}-\frac{1}{3}|c_{1}|(1-|c_{1}|^{2})z^{3}+\cdots,$

respectively. ∎

######  Remark 1.

Theorem 1 brings:

$\omega\in\mathcal{B}^{\prime}_{0}\qquad\Rightarrow\qquad|c_{3}-c_{1}c_{2}|\leq\frac{1}{3},$

while

$\omega\in\mathcal{B}_{0}\qquad\Rightarrow\qquad|c_{3}-c_{1}c_{2}|\leq 1$

follows from [4].

Similarly as Theorem 1 we can prove the next theorem.

######  Theorem 2.

If $\omega\in\mathcal{B}^{\prime}_{0}$ is of form (1) and
$\mu\in{\mathcal{C}}$, then

(6) $|c_{3}-\mu
c_{1}c_{2}|\leq\left\\{\begin{array}[]{cc}\frac{1}{48}(1+|c_{1}|)\left[9|\mu|^{2}|c_{1}|^{2}-16|c_{1}|+16\right],&0\leq|c_{1}|\leq\frac{1}{1+3/4|\mu|}\\\\[5.69054pt]
\left(\frac{1}{3}+\frac{1}{2}|\mu|\right)|c_{1}|(1-|c_{1}|^{2}),&\frac{1}{1+3/4|\mu|}\leq|c_{1}|\leq
1\end{array}\right..$

The estimate is sharp for $|c_{1}|=0$, and for
$\frac{1}{1+3/4|\mu|}\leq|c_{1}|\leq 1$ when $\mu$ is nonnegative real number.
The extremal functions are $\omega_{1}$ and $\omega_{2}$, respectively
($\omega_{1}$ and $\omega_{2}$ as defined in the proof of Theorem 1).

For $\mu=2$ in Theorem 2 we receive

######  Corollary 1.

If $\omega\in\mathcal{B}^{\prime}_{0}$ is of form (1). Then

$|c_{3}-2c_{1}c_{2}|\leq\left\\{\begin{array}[]{cc}\frac{1}{12}(1+|c_{1}|)\left[9|c_{1}|^{2}-4|c_{1}|+4\right],&0\leq|c_{1}|\leq\frac{2}{5}\\\\[5.69054pt]
\frac{4}{3}|c_{1}|(1-|c_{1}|^{2}),&\frac{2}{5}\leq|c_{1}|\leq
1\end{array}\right..$

The estimate is sharp for $|c_{1}|=0$ and for $\frac{2}{5}\leq|c_{1}|\leq 1$,
with extremal functions $\omega_{1}$ and $\omega_{2}$, respectively
($\omega_{1}$ and $\omega_{2}$ as defined in the proof of Theorem 1).

Next, for the modulus $|c_{1}c_{3}-c_{2}^{2}|$ we have the following sharp
estimate.

######  Theorem 3.

If $\omega\in\mathcal{B}^{\prime}_{0}$ is of form (1), then the following
estimate is sharp

(7)
$|c_{1}c_{3}-c_{2}^{2}|\leq\frac{1}{42}(1-|c_{1}|^{2})(3+|c_{1}|^{2}),\quad
0\leq|c_{1}|\leq 1.$

###### Proof.

Using Lemma 2 we have

$\begin{split}|c_{1}c_{3}-c_{2}^{2}|&\leq|c_{1}||c_{3}|+|c_{2}|^{2}\\\
&\leq|c_{1}|\cdot\frac{1}{3}\left(1-|c_{1}|^{2}-\frac{4|c_{2}|^{2}}{1+|c_{1}|}\right)+|c_{2}|^{2}\\\
&=\frac{1}{3}|c_{1}|(1-|c_{1}|^{2})+|c_{2}|^{2}\cdot\frac{3-|c_{1}|}{3(1+|c_{1}|)}\\\
&\leq\frac{1}{3}|c_{1}|(1-|c_{1}|^{2})+\frac{3-|c_{1}|}{3(1+|c_{1}|)}\cdot\frac{1}{4}(1-|c_{1}|^{2})^{2}\\\
&=\frac{1}{12}(1-|c_{1}|^{2})(2+|c_{1}|^{2}).\end{split}$

The equality in (7) is obtained for the function $\omega_{2}(z)$ given in
Theorem 1,

$\omega_{2}(z)=\int_{0}^{z}\frac{|c_{1}|+z}{1+|c_{1}|z}\,dz=|c_{1}|z+\frac{1}{2}(1-|c_{1}|^{2})z^{2}-\frac{1}{3}|c_{1}|(1-|c_{1}|^{2})z^{3}+\cdots.$

∎

######  Remark 2.

From (7) we have taht for every $0\leq|c_{1}|\leq 1$,

$|c_{1}c_{3}-c_{2}^{2}|\leq\frac{1}{12}\left(3-2|c_{1}|^{2}-|c_{1}|^{4}\right)\leq\frac{1}{4}.$

######  Remark 3.

As it is shown in [2], for the class ${\mathcal{U}}$ of functions
$f(z)=z+a_{2}z^{2}+a_{3}z^{3}+\cdots$ defined by the condition

$\left|\left(\frac{z}{f(z)}\right)^{2}f^{\prime}(z)-1\right|<1,\qquad
z\in{\mathbb{D}},$

we have

(8) $\frac{z}{f(z)}=1-a_{2}z-z\omega(z),$

where $\omega\in\mathcal{B}^{\prime}_{0}$ and
$\omega(z)=c_{1}z+c_{2}z^{2}+\cdots$. From (8) we can express the coefficients
$a_{3}$, $a_{4}$, and $a_{5}$, of the function $f$, depending on $a_{2}$,
$c_{1}$, $c_{2},$, $c_{3}$,…. After some calculations we receive

$|H_{3}(1)(f)|=\left|c_{1}c_{3}-c_{2}^{2}\right|\leq\frac{1}{4},$

where $H_{3}(1)(f)$ is the Hankel determinant of third order (see [3]) and
that result is the best possible. This property was the inspiration to study
the class $\mathcal{B}^{\prime}_{0}$ as a continuation of the study of the
class $\mathcal{B}_{0}$ in [4].

Similarly as in Theorem 2 we get

######  Theorem 4.

If $\omega\in\mathcal{B}^{\prime}_{0}$ is of form (1) and
$\mu\in{\mathcal{C}}$, then

(9) $|c_{1}c_{3}-\mu
c_{2}^{2}|\leq\left\\{\begin{array}[]{cc}\frac{1}{3}|c_{1}|(1-|c_{1}|^{2}),&|\mu|\leq\frac{4}{3}\frac{|c_{1}|}{1+|c_{1}|}\\\\[5.69054pt]
\frac{1}{12}\left[3|\mu|+2(2-3|\mu|)|c_{1}|^{2}-(4-3|\mu|)|c_{1}|^{4}\right],&|\mu|\geq\frac{4}{3}\frac{|c_{1}|}{1+|c_{1}|}\end{array}\right..$

Finally, for the modulus $|c_{4}-c_{2}^{2}|$ we have

######  Theorem 5.

If $\omega\in\mathcal{B}^{\prime}_{0}$ is of form (1), then

(10) $|c_{4}-c_{2}^{2}|\leq\frac{1}{4}(1-|c_{1}|^{2})$

and the estimate is sharp as the function

$\omega(z)=\int_{0}^{z}\frac{|c_{1}|+z^{3}}{1+|c_{1}|z^{3}}\,dz=|c_{1}|z+\frac{1}{4}(1-|c_{1}|^{2})z^{4}-\frac{1}{6}|c_{1}|(1-|c_{1}|^{2})z^{6}+\cdots.$

shows.

###### Proof.

Using Lemma 2, we easily get

$|c_{4}-c_{2}^{2}|\leq|c_{4}|+|c_{2}|^{2}\leq\frac{1}{4}(1-|c_{1}|^{2}-4|c_{2}|^{2})+|c_{2}|^{2}=\frac{1}{4}(1-|c_{1}|^{2}).$

∎

## References

  * [1] F. Carlson, Sur les coefficients d’une fonction born´ee dans le cercle unit´e, Ark. Mat. Astr. Fys., 27A, 1, (1940), 8 pp.
  * [2] M. Obradović and N. Tuneski, Some properties of the class U, Ann. Univ. Mariae Curie-Skłodowska Sect. A 73 (2019), 49–56.
  * [3] D.K. Thomas, N. Tuneski, and A. Vasudevarao, Univalent Functions: A Primer, De Gruyter Studies in Mathematics 69, De Gruyter, Berlin, Boston, 2018.
  * [4] P. Zaprawa, Inequalities for the Coefficients of Schwartz Functions, Bull. Malays. Math. Sci. Soc. 46 (2023), 144.


# Block Discrete Empirical Interpolation Methods

Perfect Y. Gidisu<EMAIL_ADDRESS>Michiel E. Hochstenbach
<EMAIL_ADDRESS>
Department of Mathematics and Computer Science
TU Eindhoven
The Netherlands

###### Abstract

We present two block variants of the discrete empirical interpolation method
(DEIM); as a particular application, we will consider a CUR factorization. The
block DEIM algorithms are based on the rank-revealing QR factorization and the
concept of the maximum volume of submatrices. We also present a version of the
block DEIM procedures, which allows for adaptive choice of block size.
Experiments demonstrate that the block DEIM algorithms may provide a better
low-rank approximation, and may also be computationally more efficient than
the standard DEIM procedure.

Keywords: Block DEIM, MaxVol, CUR decomposition, rank-revealing QR
factorization, low-rank approximation

##  1 Introduction

A CUR decomposition approximates a data matrix using a subset of its rows and
columns. Such factorization preserves in the reduced matrix, properties such
as interpretability, sparsity, and nonnegativity of the original data matrix.
In machine learning, one can use a CUR decomposition as an unsupervised
feature or sample selection technique. Given an $m\times n$ matrix $A$ and a
target rank $k$, a CUR factorization takes the form

$\begin{array}[]{ccccc}A&\approx&C&M&R~{},\\\ m\times n&&m\times k&k\times
k&k\times n\end{array}$

where $C$ and $R$ are small subsets of the columns and rows of $A$,
respectively. The selected columns and rows tend to capture the important
information of $A$. One standard way of constructing the factors $C$ and $R$
of a CUR factorization is to apply the discrete empirical interpolation index
selection method (DEIM) to the $k$ dominant right and left singular vectors of
$A$, respectively (Sorensen and Embree, 2016). Given the factors $C$ and $R$,
following Sorensen and Embree (2016) and others (Mahoney and Drineas, 2009;
Stewart, 1998), the middle matrix $M$ can be computed as
$(C^{T}C)^{-1}C^{T}AR^{T}(RR^{T})^{-1}$.

In this work, we propose three different block DEIM procedures. There are
three motivating factors for a block DEIM scheme:

  * •

The DEIM algorithm has a greedy nature: it locally selects the index
corresponding to the largest magnitude element of a vector. A block method
shares the same principle, but it may be less greedy since the optimization is
done over more indices instead of just one.

  * •

The block DEIM procedures that we propose will generally have higher flop
counts than the standard DEIM algorithm. However, in practice, we expect that
the block DEIM schemes will have higher flop performance than the classical
DEIM algorithm as they mainly involve BLAS-3 operations.

  * •

A block DEIM scheme may be a good solution where the DEIM procedure faces a
difficult choice when the local maximizer is (nearly) nonunique. Indeed, in
(Sorensen and Embree, 2016, Footnote 3) it is already hinted that a
potentially problematic case is if multiple entries in the vector have nearly
the same magnitude, e.g.,
$|(\mathbf{v}_{1})_{i}|\approx|(\mathbf{v}_{1})_{j}|$ for $i\neq j$, then the
DEIM scheme will pick the smaller index. However, if the other large-magnitude
entries in $\mathbf{v}_{1}$ are not significant in the subsequent $\mathbf{v}$
vectors, the corresponding indices will never be selected although they are
equally important (see Example 1.1).

Our selection criterion for the block-$b$ case will consist of maximizing the
modulus of the determinant or volume over many possible $b\times b$
submatrices. Although $b=5$ is the default, we also discuss options for
approximating the maximum volume using larger values of $b$ and taking
adaptive values of $b$.

To motivate the discussion, let us first consider a small illustrative
example.

###### Example 1.1

Consider the matrix of left singular vectors (where, e.g.,
$\varepsilon=10^{-15}$)

$U={\mbox{\scriptsize$\left[\\!\\!\begin{array}[]{cc}\frac{1}{3}\sqrt{3}+{\varepsilon}&\phantom{-}0\\\\[2.84526pt]
\frac{1}{3}\sqrt{3}&\phantom{-}\frac{1}{2}\sqrt{2}+\varepsilon\\\\[2.84526pt]
\frac{1}{3}\sqrt{3}&-\frac{1}{2}\sqrt{2}\end{array}\\!\\!\right]$}},\quad\text{corresponding
to singular values}\ \sigma_{1},\ \sigma_{2}.$

It may easily be checked that, independent of the $\sigma_{i}$, standard DEIM
will pick the first index, followed by the second. However, it is clear that
the second and third indices are a better choice, especially if $\sigma_{2}$
is close to $\sigma_{1}$. For instance, if
$A=U\cdot{\mbox{\scriptsize$\left[\\!\\!\begin{array}[]{cc}1&0\\\\[2.84526pt]
0&0.99\end{array}\\!\\!\right]$}}$, working with a determinant of $2\times 2$
submatrices gives a better result. Since

$\Big{|}\det{\mbox{\scriptsize$\left[\\!\\!\begin{array}[]{rr}\frac{1}{3}\sqrt{3}&\phantom{-}\frac{1}{2}\sqrt{2}\\\\[2.84526pt]
\frac{1}{3}\sqrt{3}&\phantom{-}\frac{1}{2}\sqrt{2}+\varepsilon\end{array}\\!\\!\right]$}}\Big{|}>\Big{|}\det{\mbox{\scriptsize$\left[\\!\\!\begin{array}[]{cc}\frac{1}{3}\sqrt{3}+{\varepsilon}&\phantom{-}0\\\\[2.84526pt]
\frac{1}{3}\sqrt{3}&\phantom{-}\frac{1}{2}\sqrt{2}+\varepsilon\end{array}\\!\\!\right]$}}\Big{|},$

a block DEIM variant with block size 2 picks the more appropriate indices 2
and 3.

The outline of the paper is as follows. Section 2 provides a brief literature
review of some existing deterministic CUR approximation algorithms. Section 3
describes how the DEIM scheme can be blocked using either a column-pivoted QR
factorization or the concept of the maximum determinant or volume of
submatrices. We discuss the newly proposed block DEIM algorithms and their
computational complexities, and state a well-known error bound for the
approximation. Section 4 reports the results from numerical experiments
evaluating the computational efficiency and approximation quality of the
various block DEIM procedures. Section 5 summarizes the key points and results
of this paper. Some main properties of the methods are summarized in Table 1.

Table 1: Overview of the various variants with their properties. All methods have the SVD as a basis. The quantity $\eta$ appears in error upper bounds and is discussed in Section 3.4.  Method $\backslash$ Properties | SVD | Block | Adapt. | Tunable | Greediness | $\eta$ | Speed
---|---|---|---|---|---|---|---
Standard DEIM | $+$ | $-$ | $-$ | $-$ | High | Medium | Moderate
B-MaxVol (Section 3.1) | $+$ | $+$ | $-$ | $+$ | Low | Low | Fast
B-RRQR (Section 3.2) | $+$ | $+$ | $-$ | $+$ | Low | Low | Fast
B-Adaptive (Section 3.3) | $+$ | $+$ | $+$ | $+$ | Medium | Medium | Variable

##  2 Review of CUR factorization algorithms

In this section, we summarize some known index selection algorithms for
computing a CUR factorization. We use MATLAB notations to index vectors and
matrices, i.e., $A(:,\mathbf{p})$ denotes the $k$ columns of $A$ whose
corresponding indices are in vector $\mathbf{p}\in\mathbb{N}_{+}^{k}$. We
denote the 2-norm norm by $\lVert\cdot\rVert$ and the notation $R^{+}$ denotes
the Moore–Penrose pseudoinverse of $R$.

###  2.1 MaxVol Algorithm

The MaxVol algorithm (Goreinov et al., 1997, 2010) is one of the most common
greedy deterministic procedures for a CUR factorization. Goreinov et al.
(1997) propose a rank-$k$ CUR factorization then referred to as pseudoskeleton
approximation based on performing the MaxVol procedure on the top $k$ right
and left singular vectors. Oseledets and Tyrtyshnikov (2010) develop a cross-
approximation scheme that alternates between searching over rows and columns
using the MaxVol method.

Given a tall-thin matrix $U\in\mathbb{R}^{m\times k}$, the MaxVol algorithm
searches for $k$ row indices such that the $k\times k$ upper submatrix
$\widehat{U}$ is dominant in $U$ (Goreinov et al., 2010). Although the
dominant property does not imply that $\widehat{U}$ has the maximum volume, it
guarantees that $\widehat{U}$ is locally optimal so that replacing any row in
$\widehat{U}$ with a row in $U$ that is already not present in $\widehat{U}$,
does not increase the volume (Goreinov et al., 2010). Algorithm 1 presents the
procedure to find a dominant submatrix using the MaxVol algorithm. In
practice, a good initialization step for the MaxVol algorithm is the $k$ rows
of the input matrix that are the pivoted rows from the LU decomposition of the
input matrix (Goreinov et al., 2010).

1:$U\in\mathbb{R}^{m\times k}$ with $m>k$, convergence tolerance $\delta$
(default 0.01)

2:$\mathbf{s}\in\mathbb{N}_{+}^{k}$ indices

3:$\mathbf{s}\leftarrow k$ first indices of the pivoted rows from the LU
decomposition of $U$

4:repeat

5: Set $\widehat{U}\leftarrow U(\mathbf{s},:)$ and $B\leftarrow
U\widehat{U}^{-1}$

6: Find the element of maximum absolute value in $B$:
$(i,j)\leftarrow\operatorname*{\arg\\!\max}|b_{ij}|$

7: if $|b_{ij}|>1$, swap rows $i$ and $j$ in $B$: $\mathbf{s}(j)=i$

8:until $\forall$ $(i,j)$: $|b_{ij}|<1+\delta$

Algorithm 1 MaxVol: Approximation to dominant submatrix (Goreinov et al.,
2010)

The parameter $\delta$ is the convergence tolerance to find pivot elements;
this parameter serves as a stopping criterion and should be sufficiently small
(a good choice can be $0.01$) (Goreinov et al., 2010). Note that by swapping
the rows as done in 7, the volume of the upper submatrix in $B$ is increased,
and also in $U$ until convergence. The most expensive part of the iterations
is 5; this needs a $k\times k$ matrix inversion and $mk^{2}$ operations for
the matrix multiplication. Goreinov et al. (2010) describe a speed
optimization process that avoids the expensive matrix multiplications and
inversions. We refer the reader to (Goreinov et al., 2010) for a more detailed
explanation of the MaxVol algorithm.

###  2.2 DEIM

The DEIM point selection method is a deterministic greedy index selection
algorithm originally presented in the context of model order reduction for
nonlinear dynamical systems (Chaturantabut and Sorensen, 2010). Sorensen and
Embree (2016) show that this procedure is a viable index selection algorithm
for constructing a CUR factorization. To select the indices for a rank-$k$ CUR
factorization, first, compute a rank-$k$ singular value decomposition (SVD) of
the original matrix. Using as input, the top-$k$ right and left singular
vectors contained in $V$ and $U$, respectively, the DEIM algorithm selects $k$
column and row indices, $\mathbf{p}$ and $\mathbf{s}$, respectively, as in
Algorithm 2 111Note that the backslash operator used in the algorithms is a
Matlab type notation for solving linear systems and least-squares problems..
The DEIM procedure selects the indices by processing the singular vectors one
at a time. We describe this method using the left singular vectors, and the
procedure on the right singular vectors follows similarly.

Starting from the leading right singular vector $\mathbf{u}_{1}$, the first
row index $s_{1}$ corresponds to the entry in $\mathbf{u}_{1}$ with the
largest magnitude, i.e.,
$|\mathbf{u}_{1}(s_{1})|=\lVert\mathbf{u}_{1}\rVert_{\infty}$ (where
$\lVert\cdot\rVert_{\infty}$ denotes the infinity norm). The remaining indices
$s_{j}$ for $j=2,\dots,k$ are selected so that each index corresponds to the
largest magnitude entry in the residual
$\mathbf{r}_{j}=\mathbf{u}_{j}-\mathbb{S}_{j-1}\mathbf{u}_{j}$, where
$\mathbb{S}_{j-1}$ is an interpolatory projector computed as
$U_{j-1}(S_{j-1}^{T}U_{j-1})^{-1}S_{j-1}^{T}$. (The linear independence of the
columns of $U$ guarantees that $S_{j-1}^{T}U_{j-1}$ is nonsingular.) As
mentioned in our third motivating factor in Section 1, in the case where we
have multiple index options, e.g.,
$|(\mathbf{r}_{j})_{1}|=|(\mathbf{r}_{j})_{2}|$, the smaller index is picked.
For further details about the DEIM scheme, we refer the reader to
(Chaturantabut and Sorensen, 2010; Sorensen and Embree, 2016).

1:$U\in\mathbb{R}^{m\times k}$, $V\in\mathbb{R}^{n\times k}$, $k\leq\min(m,n)$

2:Indices $\mathbf{p},\mathbf{s}\in\mathbb{N}_{+}^{k}$ with non-repeating
entries

3:for $j=1,\dots,k$ do

4: $\mathbf{s}(j)$ = $\operatorname*{\arg\\!\max}_{1\leq i\leq
m}~{}|(U(:,j))_{i}|$ (Iteratively pick indices)

5: $\mathbf{p}(j)$ = $\operatorname*{\arg\\!\max}_{1\leq i\leq
n}~{}|(V(:,j))_{i}|$

6: if $j<k$ thenM (Update remaining columns)

7: $U(:,j+1)=U(:,j+1)-U(:,~{}1:j)\cdot(U(\mathbf{s},\,1:j)\ \backslash\
U(\mathbf{s},\,j+1))$

8: $V(:,j+1)=V(:,j+1)-V(:,~{}1:j)\cdot(V(\mathbf{p},\,1:j)\ \backslash\
V(\mathbf{p},\,j+1))$

9: end if

10:end for

Algorithm 2 DEIM Index Selection (Sorensen and Embree, 2016)

###  2.3 Rank-revealing QR factorization

The classical truncated pivoted QR factorization is another approach for
computing a CUR factorization (Voronin and Martinsson, 2017). The QR
factorization of a matrix $A\in\mathbb{R}^{m\times n}$ with $m\geq n$ is of
the form

$\begin{array}[]{ccccc}A&\Pi&=&Q&T,\\\ m\times n&n\times n&&m\times n&n\times
n\end{array}$

where $\Pi$ is a permutation matrix, $Q$ is an orthonormal matrix, and $T$ is
an upper triangular matrix that satisfies the condition

$|T_{kk}|^{2}\geq\sum_{i=k}^{j}|T_{ij}|^{2},~{}j=k+1,\dots,n,\quad
k=1,\dots,n.$

Given $\mathbf{p}$, a vector of indices, we can express $\Pi$ as
$\Pi=I(:,\mathbf{p})$ where $I$ is the $n\times n$ identity matrix. The matrix
$A$ can then be written as a sum of rank-1 matrices

$A(:,\mathbf{p})\approx\sum_{j=1}^{n}Q(:,j)~{}T(:,j).$

Let us partition $Q$ and $T$ so that

$A(:,\mathbf{p})=\begin{bmatrix}{Q_{1}}&{Q_{2}}\end{bmatrix}\begin{bmatrix}{T_{11}}&{T_{12}}\\\
0&T_{22}\end{bmatrix}=Q_{1}\begin{bmatrix}T_{11}&T_{12}\end{bmatrix}+Q_{2}\begin{bmatrix}0&T_{22}\end{bmatrix},$

where $Q_{1}\in\mathbb{R}^{m\times k}$, $Q_{2}\in\mathbb{R}^{m\times(n-k)}$,
$T_{11}\in\mathbb{R}^{k\times k}$, $T_{12}\in\mathbb{R}^{k\times(n-k)}$, and
$T_{22}\in\mathbb{R}^{(n-k)\times(n-k)}$. We have that $\lVert A\Pi-
Q_{1}\begin{bmatrix}T_{11}&T_{12}\end{bmatrix}\rVert\leq\lVert T_{22}\rVert$,
which implies that $Q_{1}$ is an approximation of the range of $A$ and as long
as $\lVert T_{22}\rVert$ is small, $A(:,\mathbf{p})$ can be approximated by
$Q_{1}\begin{bmatrix}T_{11}&T_{12}\end{bmatrix}$. For an arbitrary $k$, the
best rank-$k$ approximation of $A$ $(A_{k})$ from the SVD gives $\lVert
A-A_{k}\rVert=\sigma_{k+1}(A)$ where $\sigma_{k+1}$ is the $(k+1)$st singular
value of $A$. It is always the case that $\sigma_{k+1}(A)\leq\lVert
T_{22}\rVert$ (Voronin and Martinsson, 2017). A factorization is said to be a
rank-revealing QR (RRQR) factorization if the ratio $\lVert
T_{22}\rVert\,/\,\sigma_{k+1}(A)$ is guaranteed to be bounded (Gu and
Eisenstat, 1996). Note that a rank-revealing QR factorization is based on
selecting certain well-conditioned submatrices (Chandrasekaran and Ipsen,
1994); which is the underlying objective of the maximum volume concept.

To compute a pivoted QR factorization, one may use a column-pivoted
Gram–Schmidt algorithm (Golub and Van Loan, 2012). Although there are counter-
examples, the classical Gram–Schmidt process usually results in an RRQR
factorization (Gu and Eisenstat, 1996). Notice that from the matrix partition,
$Q_{1}T_{11}$ equals the first $k$ columns of $A(:,\mathbf{p})$. To construct
a CUR decomposition, one may apply a Gram–Schmidt-based pivoted QR algorithm
to the matrices $A$ and $A^{T}$ to obtain the matrices $C$ and $R$,
respectively (Berry et al., 2005; Stewart, 1998). Alternatively, Voronin and
Martinsson (2017) show how to compute a CUR factorization via a two-sided
interpolative decomposition (ID), which in turn can be constructed from a
pivoted QR factorization. Using the Q-DEIM approach (Drmac and Gugercin,
2016), one can apply a column-pivoted QR procedure on the transpose of the
leading $k$ right and left singular vectors to find the indices for
constructing the factors $C$ and $R$ in CUR approximation.

We will now use the tools described in this section to design our block
variants of the standard DEIM algorithm in the following section.

##  3 Block DEIM

This section introduces three new block variants of the DEIM procedure. To
design our block DEIM algorithms, we combine the index selection algorithms
discussed in Section 2. The block DEIM algorithms select a $b$-size set of
indices at each step using a block of singular vectors and then update the
subsequent block of vectors using the oblique projection technique in the
standard DEIM procedure. The main difference between the block DEIM and the
traditional DEIM scheme resides in the $b$-size indices selection. In the
standard DEIM, we select the indices by processing one singular vector at a
time, each iteration step produces an index whilst the block version processes
a block of singular vectors at a time, and each step provides an index set of
size $b$. As a result, this may lead to a different selection of indices. For
ease of presentation, we assume that the number of singular vectors $k$ is a
multiple of the block size $b$ and $2\leq b<k$.

###  3.1 Block DEIM based on Maximum Volume (B-DEIM-MaxVol)

We first introduce a block DEIM procedure based on the MaxVol scheme (Goreinov
et al., 2010) described in Section 2 for finding a submatrix with the largest
possible determinant in a thin-tall matrix. This algorithm takes as inputs the
left and right singular vectors, a block size $b$, a rank parameter $k$, and a
convergence tolerance $\delta$. For a $b$-block of singular vectors, the
procedure returns $b$ row indices at each step such that the final $b\times b$
upper submatrix in the $b$-block of singular vectors is dominant. We then use
the oblique projection technique as in the standard DEIM scheme to update the
subsequent blocks of singular vectors. At the end of the iteration, the
algorithm returns a column and row index set of size $k$.

The computational cost of this algorithm is dominated by two calls of the
MaxVol procedure and the block updates. The initialization step of Algorithm 1
requires a permutation of the input matrix, which is done via the LU
factorization. Given an $n\times b$ matrix, the LU decomposition requires
$\mathcal{O}(nb^{2})$ operations. Furthermore, the dominant cost of each
iteration in the MaxVol algorithm is the multiplication of an $n\times b$
matrix and a $b\times b$ matrix, for a cost of $\mathcal{O}(nb^{2})$
operations. Let $\kappa\in\mathbb{N}_{+}$ denote the number of iterations
performed. The computational complexity of the MaxVol procedure is
$\mathcal{O}(\kappa nb^{2})$. The cost of the two calls of MaxVol is
$\mathcal{O}(\kappa(m+n)b^{2})$ and the block updates cost
$\mathcal{O}((m+n)kb)$. Given that we have $k/b$ iterations in Algorithm 3,
the total cost of the B-DEIM-MaxVol algorithm is $\mathcal{O}((m+n)(\kappa
bk+k^{2}))$. A crude bound on the number of iterations in Algorithm 1 is
$\kappa\leq(\log|{\sf det}(\widehat{U}_{\sf dom})|-\log|{\sf
det}(\widehat{U}_{\sf ini})|)/\log(1+\delta)$, where $\widehat{U}_{\sf ini}$
is the submatrix at the initialization step and $\widehat{U}_{\sf dom}$ is the
dominant submatrix in Algorithm 1, respectively (Goreinov et al., 2010). The
procedure of block DEIM using MaxVol is summarized in Algorithm 3. In 4 and 5
of the algorithm, we apply Algorithm 1 as implemented by Kramer and Gorodetsky
(2018) to the block of singular vectors.

1:$U\in\mathbb{R}^{m\times k}$, $V\in\mathbb{R}^{n\times k}$,
$k\leq\min(m,n)$, block size $b$ (with $b\,|\,k$), convergence tolerance
$\delta$ (default 0.01)

2:Indices $\mathbf{p},\mathbf{s}\in\mathbb{N}_{+}^{k}$ with non-repeating
entries

3:for $j=1,\dots,k/b$ do

4: $\mathbf{s}((j-1)b+1:bj)={\sf MaxVol}(U(:,(j-1)b+1:bj),\delta)$

5: $\mathbf{p}((j-1)b+1:bj)={\sf MaxVol}(V(:,(j-1)b+1:bj),\delta)$

6: $U(:,bj+1:bj+b)=U(:,bj+1:bj+b)-U(:,1:bj)\cdot(U(\mathbf{s},\,1:bj)\
\backslash\ U(\mathbf{s},\,bj+1:bj+b))$

7: $V(:,bj+1:bj+b)=V(:,bj+1:bj+b)-V(:,1:bj)\cdot(V(\mathbf{p},\,1:bj)\
\backslash\ V(\mathbf{p},\,bj+1:bj+b))$

8:end for

Algorithm 3 Block DEIM index selection based on MaxVol (B-DEIM-MaxVol)

###  3.2 Block DEIM based on RRQR (B-DEIM-RRQR)

The MaxVol algorithm seeks to find a good approximation of a submatrix with
the maximum volume in a given input matrix. It is noted in (Drmac and
Gugercin, 2016, Remark 2.3) that the pivoting in the RRQR factorization can be
interpreted as a greedy volume maximizing scheme. Practical experience shows
that the index selection via a column-pivoted QR factorization may be more
computationally efficient than the MaxVol method.

Additionally, Khabou et al. (2013) propose a block LU factorization with panel
rank-revealing pivoting (LU-PRRP) procedure by employing the strong rank-
revealing QR (RRQR) panel factorization. The algorithm computes the block LU-
PRRP as follows: at each step of the block factorization, a block of columns
is factored by performing strong RRQR factorization (Gu and Eisenstat, 1996)
on its transpose. Pivoting is done by applying the permutation matrix returned
on the entire original matrix followed by an update of the trailing matrix. As
noted in (Sorensen and Embree, 2016), the DEIM index selection procedure is
equivalent to the index selection of partially pivoted LU decomposition.

With this knowledge, we propose an alternative block DEIM variant (B-DEIM-
RRQR) that merges the components of an RRQR factorization and the standard
DEIM scheme. At each step of the index selection process, given a block size
$b$ and a block of singular vectors, we perform a QR decomposition with column
pivoting on the transpose of a block of singular vectors; the first $b$
columns give the selected indices. We then update the next block of singular
vectors comparatively to how the standard DEIM algorithm updates subsequent
singular vectors by using interpolatory projection.

The cost of the B-DEIM-RRQR scheme is dominated by two QR factorizations and
block updates. Given an $n\times b$ matrix, a QR factorization requires
$\mathcal{O}(nb^{2})$ operations hence the cost of the two QR decompositions
here is $\mathcal{O}((m+n)b^{2})$ and the block updates cost
$\mathcal{O}((m+n)kb)$. Since there are $k/b$ iterations in Algorithm 4, the
total cost of the B-DEIM-RRQR scheme is therefore
$\mathcal{O}((m+n)(bk+k^{2}))$. The B-DEIM-RRQR may be viewed as a hybrid
standard DEIM and Q-DEIM scheme. The Q-DEIM selects $k$ indices by performing
one column-pivoted QR while our B-DEIM-RRQR algorithm selects the $k$ indices
by performing $k/b$ rounds of column-pivoted QR. Note that when the block size
$b=k$, this B-DEIM-RRQR algorithm is just the $Q-DEIM$ scheme. Algorithm 4
summarizes the procedure. The permutation matrices $\Pi_{\mathbf{s}}$ and
$\Pi_{\mathbf{p}}$ in 5 and 6 are from the column-pivoted QR of $U^{T}$ and
$V^{T}$, respectively.

1:$U\in\mathbb{R}^{m\times k}$, $V\in\mathbb{R}^{n\times k}$,
$k\leq\min(m,n)$, block size $b$ (with $b\,|\,k$)

2:Indices $\mathbf{p},\mathbf{s}\in\mathbb{N}_{+}^{k}$ with non-repeating
entries

3:for $j=1,\dots,k/b$ do

4: Perform a column-pivoted QR on $U(:,(j-1)b+1:bj)^{T}$ and
$V(:,(j-1)b+1:bj)^{T}$

5: $\mathbf{s}((j-1)b+1:bj)=\Pi_{\mathbf{s}}(1:b)$

6: $\mathbf{p}((j-1)b+1:bj)=\Pi_{\mathbf{p}}(1:b)$

7: $U(:,bj+1:bj+b)=U(:,bj+1:bj+b)-U(:,1:bj)\cdot(U(\mathbf{s},\,1:bj)\
\backslash\ U(\mathbf{s},\,bj+1:bj+b))$

8: $V(:,bj+1:bj+b)=V(:,bj+1:bj+b)-V(:,1:bj)\cdot(V(\mathbf{p},\,1:bj)\
\backslash\ V(\mathbf{p},\,bj+1:bj+b))$

9:end for

Algorithm 4 Block DEIM index selection based on RRQR (B-DEIM-RRQR)

###  3.3 Adaptive choice of block size (AdapBlock)

To combine the strength of standard DEIM with a block version, we also
consider adaptive choices for the block size $b$. In particular, we propose
the following variant: perform block DEIM if for the singular vector
$\mathbf{v}_{j}$ being considered the two largest elements are nearly equal,
i.e.,

$\left\\{\begin{array}[]{ll}{\rm block~{}DEIM}&\quad{\rm if}\
|(\mathbf{v}_{j})_{i}|\approx|(\mathbf{v}_{j})_{\ell}|\quad{\rm for}\
i\neq\ell,\\\\[2.84526pt] {\rm
standard~{}DEIM}&\quad\text{otherwise}.\end{array}\right.$

In Algorithm 5, we show an implementation of the adaptive block DEIM using the
block DEIM variants discussed in Section 3. The parameter $\rho$ is the
desired lower bound on the ratio
$|(\mathbf{v}_{j})_{i}|\,/\,|(\mathbf{v}_{j})_{\ell}|$ for $i\neq\ell$. We
note that although our criterion for switching from a standard DEIM scheme to
a block DEIM method is based on how close the two largest entries (magnitude)
in the vector being considered are, other criteria can be used.

1:$U\in\mathbb{R}^{m\times k}$, $V\in\mathbb{R}^{n\times k}$,
$k\leq\min(m,n)$, block size $b$, $\rho$ (default 0.95), tolerance $\delta$
(default 0.01), method (MaxVol or QR)

2:Indices $\mathbf{p},\mathbf{s}\in\mathbb{N}_{+}^{k}$ with non-repeating
entries

3:$j=1$

4:while $j\leq k$ do

5: $U(:,j)=U(:,j)-U(:,1:j-1)\cdot(U(\mathbf{s}(1:j-1),1:j-1)\ \backslash\
U(\mathbf{s}(1:j-1),j))$ ($U_{0}=0$)

6: $[\mathbf{u},\,{\sf ind}]={\sf sort}(|U(:,j)|)$ (in descending order)

7: if $(j=k)\ {\bf or}\ (u_{2}<\rho\,u_{1})$ then $\mathbf{s}(j)={\sf
ind}(1)$; $j=j+1$

8: else

9: ${\rm cols}=j+b-1$

10: $U(:,j+1:{\rm cols})=U(:,j+1:{\rm
cols})-U(:,1:j-1)\cdot(U(\mathbf{s}(1:j-1),1:j-1)\ \backslash\
U(\mathbf{s}(1:j-1),j+1:{\rm cols}))$

11: if method = MaxVol

12: $\mathbf{s}(j:{\rm cols})={\sf MaxVol}(U(:,j:{\rm cols}),\,\delta)$

13: else

14: Perform a column-pivoted QR on $U(:,j:{\rm cols})^{T}$

15: $\mathbf{s}(j:{\rm cols})=\Pi_{\mathbf{s}}(1:b)$

16: end if

17: $j=j+b$

18: end if

19:end while

20:Repeat the procedure on $V$ to get indices $\mathbf{p}$

Algorithm 5 Adaptive Block DEIM index selection (AdapBlock)

###  3.4 Error Bounds

The following proposition restates a known theoretical error bound for a CUR
approximation, which holds for the block DEIM algorithms proposed in this
paper. A detailed constructive proof is in (Sorensen and Embree, 2016), but we
provide the necessary details here for the reader’s convenience. Let
$P\in\mathbb{R}^{n\times k}$ and $S\in\mathbb{R}^{m\times k}$ be matrices with
some columns of the identity indexed by the indices selected by employing any
of the three block DEIM algorithms.

###### Proposition 3.1

(Sorensen and Embree, 2016, Thm. 4.1) Given $A\in\mathbb{R}^{m\times n}$ and a
target rank $k$, let $V\in\mathbb{R}^{n\times k}$ and $U\in\mathbb{R}^{m\times
k}$ be the leading $k$ right and left singular vectors of $A$. Suppose $C=AP$
and $R=S^{T}\\!A$ are of full rank, and $V^{T}\\!P$ and $S^{T}U$ are
nonsingular, then, with $M=C^{+}\\!AR^{+}$, a rank-$k$ CUR decomposition
constructed by either of the block DEIM schemes presented in this paper
satisfies

$\lVert
A-CMR\rVert\leq(\eta_{\mathbf{s}}+\eta_{\mathbf{p}})\,\sigma_{k+1}\quad
with\quad\eta_{\mathbf{s}}<\sqrt{\tfrac{nk}{3}}\,2^{k}~{},\quad\eta_{\mathbf{p}}<\sqrt{\tfrac{mk}{3}}\,2^{k},$

where $\eta_{\mathbf{p}}=\lVert(V^{T}\\!P)^{-1}\rVert$,
$\eta_{\mathbf{s}}=\lVert(S^{T}U)^{-1}\rVert$.

Proof  Let $\mathbb{P}=P(V^{T}P)^{-1}V^{T}$ and
$\mathbb{S}=U(S^{T}U)^{-1}S^{T}$ be oblique projectors. Note that
$V^{T}\mathbb{P}=V^{T}$ and $\mathbb{S}U=U$, implying that
$V^{T}(I-\mathbb{P})=0$ and $(I-\mathbb{S})U=0$. Using $M=C^{+}\\!AR^{+}$, we
have

$\displaystyle\lVert A-CMR\rVert$ $\displaystyle=\lVert
A-CC^{+}AR^{+}R\rVert=\lVert(I-CC^{+})A+CC^{+}A(I-R^{+}R)\rVert$
$\displaystyle\leq\lVert(I-CC^{+})A\rVert+\lVert CC^{+}\rVert\,\lVert
A(I-R^{+}R)\rVert.$

Leveraging the fact that $CC^{+}$ is an orthogonal projector, $\lVert
CC^{+}\rVert=1$ and (Sorensen and Embree, 2016, Lemma 4.1 and 4.2)

$\displaystyle\lVert(I-CC^{+})A\rVert$ $\displaystyle\leq\lVert
A(I-\mathbb{P})\rVert=\lVert A(I-VV^{T})(I-\mathbb{P})\rVert$
$\displaystyle\leq\lVert(V^{T}\\!P)^{-1}\rVert\,\lVert A(I-VV^{T})\rVert,$
$\displaystyle\lVert A(I-R^{+}\\!R)\rVert$
$\displaystyle\leq\lVert(I-\mathbb{S})A\rVert=\lVert(I-\mathbb{S})(I-UU^{T})A\rVert$
$\displaystyle\leq\lVert(S^{T}U)^{-1}\rVert\,\lVert(I-UU^{T})A\rVert,$

we have that

$\lVert A-CMR\rVert\leq\lVert(V^{T}\\!P)^{-1}\rVert\,\lVert
A(I-VV^{T})\rVert+\lVert(S^{T}U)^{-1}\rVert\,\lVert(I-UU^{T})A\rVert.$

Since $U$ and $V$ contain the leading $k$ left and right singular vectors,
respectively, $\lVert(I-UU^{T})A\rVert=\lVert A(I-VV^{T})\rVert=\sigma_{k+1}$.
Hence

$\lVert
A-CMR\rVert\leq(\lVert(V^{T}\\!P)^{-1}\rVert+\lVert(S^{T}U)^{-1}\rVert)\cdot\sigma_{k+1}.$

Proposition 3.1 suggests that the quality of an index selection method may be
assessed using the error constants $\eta_{\mathbf{s}}$ and $\eta_{\mathbf{p}}$
Sorensen. Fig. 1 illustrates the difference in the values of
$\eta_{\mathbf{s}}=\lVert(S^{T}U)^{-1}\rVert$ computed by DEIM, B-DEIM-RRQR,
and B-DEIM-MaxVol using 50 randomly generated orthonormal matrices of size
$10000\times 100$. In most cases, the block DEIM variants provide a smaller
values of $\eta_{\mathbf{s}}$; this may be an indication that the indices
picked by the block DEIM variants are better for approximating than those
selected by the DEIM scheme. We also observe that the values of
$\eta_{\mathbf{s}}$ become smaller as we increase the block sizes.

$10$$20$$30$$40$$50$$80$$100$$120$$140$$160$$\|(S^{T}U)^{-1}\|$DEIMB-DEIM-
RRQRB-DEIM-MaxVol (a) Block size 2

$10$$20$$30$$40$$50$$80$$100$$120$$140$$160$ (b) Block size 5

$10$$20$$30$$40$$50$$80$$100$$120$$140$$160$$\|(S^{T}U)^{-1}\|$ (c) Block size
10

$10$$20$$30$$40$$50$$80$$100$$120$$140$$160$ (d) Block size 20

Figure 1: Comparison of the value
$\eta_{\mathbf{s}}=\lVert(S^{T}U)^{-1}\rVert$ in DEIM, B-DEIM-RRQR, and
B-DEIM-MaxVol with different block sizes using 50 random orthonormal matrices
of size $10000\times 100$.

##  4 Numerical Experiments

In this section, we conduct several sets of illustrative experiments to show
the effectiveness of the block DEIM variants, i.e., B-DEIM-MaxVol, B-DEIM-
RRQR, and AdapBlock, proposed in this paper using synthetic and real data
sets. We evaluate the algorithms by applying them to data analysis problems in
several application domains: recommendation system analysis, text
categorization and information retrieval analysis, stocks, and image
compression. We use real-world sparse and dense data matrices with sizes
ranging from mid to large scale, an overview is presented in Table 2 . We
compare the performance of our algorithms for constructing a CUR approximation
with three state-of-the-art deterministic algorithms: DEIM (Sorensen and
Embree, 2016), Q-DEIM (Drmac and Gugercin, 2016), and MaxVol (Goreinov et al.,
2010). All these algorithms require the leading $k$ right and left singular
vectors to construct a rank-$k$ CUR factorization. We use these two evaluation
criteria: the rank-$k$ approximation relative error $\lVert
A-CMR\rVert\,/\,\lVert A\rVert$; the computational efficiency, i.e., the
runtime scaling for the rank parameter $k$. Here, the runtime measures the
time it takes each algorithm to select the desired number of column and row
indices. We do not consider the run time for computing the singular vectors
since all the methods we consider require the SVD. Our experiments are not
meant to be exhaustive; however, they provide clear evidence that the block
DEIM schemes proposed in this paper may provide a better low-rank
approximation and may be computationally more efficient.

In the implementation, we perform the column-pivoted QR factorization and the
truncated SVD using the MATLAB built-in functions qr and svds (Baglama and
Reichel, 2005), respectively (and svd for small cases). For the MaxVol
algorithm, we use a MATLAB implementation by Kramer and Gorodetsky (2018) made
available on GitHub 222github.com/bokramer/CURERA/blob/master/maxvol.m. Unless
otherwise stated, in all the experiments we use as default block size $b=5$,
the AdapBlock method parameter $\rho=0.95$, and the MaxVol scheme convergence
tolerance $\delta=0.01$.

Table 2: Various examples and dimensions considered. Exp. | Domain | Matrix | $m$ | $n$
---|---|---|---|---
1 | Recommendation system | Dense | 14116 |  12100
2 | Text categorization | Sparse |  18293 | 18933
3 | Image compression | Sparse | 10000 |  12784
4 | Stock prices | Dense |  13165 |  12440
5 | Synthetic | Dense |  12000 |  14000

###### Experiment 4.1

Our first data set is from the recommendation system analysis domain, where
one is usually interested in making service or purchase recommendations to
users. One of the most common techniques for recommendation systems is
collaborative filtering, which involves recommending to users items that
customers with similar preferences liked in the past.

The Jester data set is often used as a benchmark for recommendation system
research (Goldberg et al., 2001). This data matrix consists of 73421 users and
their ratings for 100 jokes. We only consider users who have ratings for all
100 jokes resulting in a $14116\times 100$ matrix. We center the matrix by
subtracting the mean of each column from all entries in that column.

From Fig. 2, we observe that the approximation qualities of the block DEIM
methods are more accurate compared to the state-of-art methods. In comparison
to the original DEIM scheme, both the B-DEIM-MaxVol and B-DEIM-RRQR algorithms
enjoy considerable lower runtimes while their adaptive variants have similar
runtimes to the DEIM procedure. In Fig. 2a, we notice that both the B-DEIM-
MaxVol algorithm and its adaptive variant are more efficient than the standard
MaxVol scheme. Here, by using the block DEIM variants, we gain in both
accuracy and speed compared to the standard MaxVol procedure. From Fig. 2b, we
see that the B-DEIM-RRQR method is just as efficient as the Q-DEIM procedure
while providing a more accurate approximation.

$10$$20$$30$$40$$50$$0.2$$0.3$$0.4$$0.5$$0.6$$\|A-CUR\|\ /\ \|A\|$DEIMB-DEIM-
MaxVolMaxVolAdapBlock-MaxVol

$10$$20$$30$$40$$50$$0$$5$$10$$15$$\cdot 10^{-2}$time (s)

(a) MaxVol methods

$10$$20$$30$$40$$50$$0.2$$0.3$$0.4$$0.5$$0.6$$\|A-CUR\|\ /\ \|A\|$DEIMB-DEIM-
RRQRQ-DEIMAdapBlock-RRQR

$10$$20$$30$$40$$50$$0$$2$$4$$6$$\cdot 10^{-2}$time (s)

(b) RRQR methods

Figure 2: Relative approximation errors (left) and runtimes (right) as a
function of $k$ for the block DEIM CUR approximation algorithms compared with
some standard CUR approximation algorithms using the Jester data set.

###### Experiment 4.2

The second data set is from text categorization and information retrieval
applications. In these data analysis problems, a “bag of words” is commonly
used to represent the documents. We use the data set from the Reuters-21578
text categorization collection, which is a collection of documents that
appeared on Reuters newswire in 1987. This data set has been heavily used as a
benchmark in the text classification community. The corpus contains 21578
documents in 135 categories. We perform our experiment using the preprocessed
data set containing 18933 unique terms and 8293 documents (Cai et al., 2005).
We scale the rows of this $8293\times 18933$ sparse matrix to have a unit
length. This matrix is much larger than the one from the previous section, and
is also less well-approximated by a low-rank matrix; the first 50 principal
components explain less than 50% of the overall variance in the data set and
less than 80% is captured by the top 500 principal components. (However,
spectral methods are often applied to this data.) The matrix is very sparse,
with approximately 0.25% nonzero elements.

For this data set, just as in the previous experiment we observe that the
columns and rows sampled by the block DEIM variants lead to considerably
smaller reconstruction errors than the existing methods. From Fig. 3a, it is
evident that the standard MaxVol method yields better approximation accuracy
than the block DEIM schemes when $k$ is small but as $k$ increases, the
approximation quality becomes worse than the those of the block DEIM schemes.
We also see that the error of the MaxVol and Q-DEIM approximations do not
necessarily decrease monotonically as the rank $k$ increases. The efficiency
of the algorithms is similar to previous experiments; the adaptive variants
have similar runtimes as the DEIM whilst the B-DEIM-MaxVol and B-DEIM-RRQR
schemes are more speed efficient. Once more, the MaxVol algorithm is the least
efficient method. Fig. 3b illustrates better runtimes for the AdapBlock-RRQR
algorithm as $k$ gets larger.

$20$$40$$60$$80$$100$$0.2$$0.25$$0.3$$0.35$$\|A-CUR\|\ /\ \|A\|$DEIMB-DEIM-
MaxVolMaxVolAdapBlock-MaxVol

$20$$40$$60$$80$$100$$0$$0.2$$0.4$$0.6$time (s)

(a) MaxVol methods

$20$$40$$60$$80$$100$$0.2$$0.25$$0.3$$0.35$$\|A-CUR\|\ /\ \|A\|$DEIMB-DEIM-
RRQRQ-DEIMAdapBlock-RRQR

$20$$40$$60$$80$$100$$0$$0.2$$0.4$time (s)

(b) RRQR methods

Figure 3: Relative approximation errors (left) and runtimes (right) as a
function of $k$ for the block DEIM CUR approximation algorithms compared with
some standard CUR approximation algorithms using the Reuters-21578 data set.

###### Experiment 4.3

For this experiment, we use the popular MNIST test set, which consists of
10000 images of handwritten digits from 0 to 9, each image of size $28\times
28$. The matrix is flattened and normalized, resulting in a $10000\times 784$
sparse nonnegative matrix with entries in the interval [0, 1].

From Fig. 4, we observe that the block DEIM algorithms produce slightly better
approximation accuracy than the standard DEIM and MaxVol, and enjoy
considerably lower approximation error than the Q-DEIM scheme. In terms of
efficiency, compared with the classical DEIM and MaxVol, the block DEIM
procedures have appreciable lower runtimes. The Q-DEIM demonstrates the lowest
runtimes among all the algorithms but produces the worst approximation
accuracy. The improved efficiency gained by using the block DEIM algorithms
instead of the standard DEIM scheme is apparent as $k$ increases.

$20$$40$$60$$80$$100$$10$$15$$20$$25$$\cdot 10^{-2}$$\|A-CUR\|\ /\
\|A\|$DEIMB-DEIM-MaxVolMaxVolAdapBlock-MaxVol

$20$$40$$60$$80$$100$$0$$20$$40$$\cdot 10^{-2}$time (s)

(a) MaxVol methods

$20$$40$$60$$80$$100$$10$$20$$30$$\cdot 10^{-2}$$\|A-CUR\|\ /\ \|A\|$DEIMB-
DEIM-RRQRQ-DEIMAdapBlock-RRQR

$20$$40$$60$$80$$100$$0$$5$$10$$15$$\cdot 10^{-2}$time (s)

(b) RRQR methods

Figure 4: Relative approximation errors (left) and runtimes (right) as a
function of $k$ for the block DEIM CUR approximation algorithms compared with
some standard CUR approximation algorithms using the MNIST data set.

###### Experiment 4.4

Our final data set is the historical stock prices for the S&P 500 index, one
of world’s most famous financial benchmarks. This stock market index tracks
the performance of 500 large companies listed on the stock exchanges in the
United States. We collect historical adjusted closing prices for the 500
stocks of the index from Jan 4, 2010, to Jul 29, 2022, a total of 3165 days.
The index includes multiple classes of stock of some constituent companies; as
a result, there are than 500 stocks in the gauge. We construct a $3165\times
503$ date-by-stock matrix, where the $(i,j)$th entry of this matrix represents
the price of the $j$th stock on the $i$th date. We remove stocks with missing
values resulting in a $3165\times 440$ dense matrix, which we center by
subtracting the mean of each column from its the entries. Numerically, the
resulting matrix is quite low-rank; less than ten principal components capture
99% of the total variance in the data set.

Fig. 5 reveals that the block DEIM variants provide smaller approximation
errors compared to Q-DEIM and the original DEIM schemes; with the latter
producing the worst approximation quality. We also observe that, for this
dataset, the standard MaxVol algorithm and its block variants yield similar
approximation quality. Table 3 shows the stocks that achieved the least
approximation error among all 7 methods for a rank-10 CUR approximation. Both
the B-DEIM-RRQR and AdapBlock-RRQR schemes achieved this result with the same
stocks but have slightly different rankings for the stocks. Note that four
stocks from the Consumer Cyclical sector are picked; this seems to correspond
to the fact that this sector is quite diverse, and its behavior may be
representative of the S&P 500 index as a whole.

$5$$10$$15$$20$$25$$30$$2$$4$$6$$8$$10$$\cdot 10^{-2}$$\|A-CUR\|\ /\
\|A\|$DEIMB-DEIM-MaxVolMaxVolAdapBlock-MaxVol (a) MaxVol methods

$5$$10$$15$$20$$25$$30$$2$$4$$6$$8$$10$$\cdot 10^{-2}$DEIMB-DEIM-RRQRQ-
DEIMAdapBlock-RRQR (b) RRQR methods

Figure 5: Relative approximation errors as a function of $k$ for the block DEIM CUR approximation algorithms compared with some standard CUR approximation algorithms using the S&P 500 index stocks data set. Table 3: The ten stocks that yield the least relative approximation error for $k=10$ The first two columns are the rankings of the various stocks by the B-DEIM-RRQR and the AdapBlock-RRQR methods.  B-DEIM | AdapBlock | Symbol | Company | Sector
---|---|---|---|---
11 |  11 | NVR | NVR, Inc. | Consumer Cyclical
12 |  12 | BKNG | Booking, Inc. | Consumer Cyclical
13 |  14 | AZO | AutoZone, Inc. | Consumer Cyclical
14 |  13 | CMG | Chipotle Mexican Grill, Inc. | Consumer Cyclical
15 |  18 | MKTX | MarketAxess, Inc. | Financial Services
16 |  15 | REGN | Regeneron Pharma., Inc. | Healthcare
17 |  16 | MTD | Mettler-Toledo Int. Inc. | Healthcare
18 |  17 | NFLX | Netflix, Inc. | Communication Services
19 | 10 | BIIB | Biogen Inc. | Healthcare
10 |  19 | GE | General Electric Company | Industrials

###### Experiment 4.5

Following the experiments in (Voronin and Martinsson, 2017), our test matrix
in this experiment is a full-rank data set $A\in\mathbb{R}^{2000\times 4000}$
that has the structure of the SVD, i.e., $A=U\Sigma V^{T}$. The matrices $U$
and $V$ have random orthonormal columns obtained via a QR factorization of a
random Gaussian matrix, and the diagonal matrix $\Sigma$ has entries that are
logspace ranging from 1 to $10^{-3}$.

Using two of the block DEIM algorithms proposed: the B-DEIM-MaxVol and B-DEIM-
RRQR, in Fig. 6, we investigate how varying block sizes, i.e.,
$b=(2,\,5,\,10,\,20)$ may affect their approximation quality and computational
efficiency. For each fixed block size, maintaining the properties of $A$, we
generate five different test cases and compute the averages of the evaluation
criteria for the range of $k$ values.

We observe that for increasing block sizes, both algorithms become
considerably faster. On the other hand, the approximation quality of the
varying block sizes may not degrade significantly. In this experiment, given
the various values of $k$, the errors are almost similar irrespective of the
block size.

$200$$400$$600$$800$$0$$5$$10$time (s)b=2b=5b=10b=20

$200$$400$$600$$800$$0.2$$0.4$$0.6$$0.8$$1$$\|A-CUR\|\ /\ \|A\|$

$200$$400$$600$$800$$0$$2$$4$$6$$8$$10$time (s)b=2b=5b=10b=20

$200$$400$$600$$800$$0.2$$0.4$$0.6$$0.8$$1$$\|A-CUR\|\ /\ \|A\|$

Figure 6: Runtimes and average approximation errors for the B-DEIM-RRQR (up)
and B-DEIM-MaxVol (down) CUR approximation algorithms as a function of $k$ for
varying block sizes using large matrices of size $2000\times 4000$.

##  5 Conclusions

This paper presents various block variants of the discrete empirical
interpolation method for computing CUR decompositions. We exploit the
advantages of the classical DEIM procedure, a column-pivoted QR decomposition,
and the concept of maximum determinant or volume of submatrices to develop
these block variants. We have then presented a version of the block DEIM,
which allows for adaptive choice of block size.

We perform the following procedures in the block DEIM based on RRQR; at each
iteration step, we compute a QR factorization with column pivoting on the
transpose of a block of singular vectors to obtain the indices corresponding
to the first $b$ columns. Then, we update the next block of vectors using the
interpolatory projection technique in the DEIM algorithm (repeat these two
steps until all indices are selected). A similar procedure is used in the
block DEIM based on MaxVol; the difference here is, instead of using a column-
pivoted QR decomposition, we use the MaxVol method.

Numerical experiments illustrate that the accuracy of a CUR factorization
using the newly proposed block DEIM procedures may be better than the
classical DEIM, MaxVol, and Q-DEIM schemes. The experiments also demonstrate
that the block variants, regarding computational speed, may have an advantage
over the standard DEIM and MaxVol algorithms. Relative to the Q-DEIM
algorithm, the B-DEIM-RRQR scheme yields lesser approximation errors while
maintaining comparable runtimes. Using the B-DEIM-RRQR and B-DEIM-MaxVol
methods, we have also illustrated how increasing the block size improves the
speed of the algorithms but may not necessarily degrade the approximation
quality significantly. A schematic overview of some properties of the various
methods is displayed in Table 1. A Matlab code of the proposed algorithms is
available via github.com/perfectyayra.

## References

  * Baglama and Reichel (2005) J. Baglama and L. Reichel.  Augmented implicitly restarted Lanczos bidiagonalization methods.  _SIAM J. Sci. Comput._ , 27(1):19–42, 2005.
  * Berry et al. (2005) M. W. Berry, S. A. Pulatova, and G. W. Stewart.  Algorithm 844: Computing sparse reduced-rank approximations to sparse matrices.  _ACM Trans. Math. Softw._ , 31(2):252–269, 2005\.
  * Cai et al. (2005) D. Cai, X. He, and J. Han.  Document clustering using locality preserving indexing.  _IEEE Trans. Knowl. Data Eng._ , 17(12):1624–1637, 2005.
  * Chandrasekaran and Ipsen (1994) S. Chandrasekaran and I. C. F. Ipsen.  On rank-revealing factorisations.  _SIAM J. Matrix Anal. Appl._ , 15(2):592–622, 1994.
  * Chaturantabut and Sorensen (2010) S. Chaturantabut and D. C. Sorensen.  Nonlinear model reduction via discrete empirical interpolation.  _SIAM J. Sci. Comput._ , 32:2737–2764, 2010.
  * Drmac and Gugercin (2016) Z. Drmac and S. Gugercin.  A new selection operator for the discrete empirical interpolation method—Improved a priori error bound and extensions.  _SIAM J. Sci. Comput_ , 38(2):A631–A648, 2016\.
  * Goldberg et al. (2001) K. Goldberg, T. Roeder, D. Gupta, and C. Perkins.  Eigentaste: A constant time collaborative filtering algorithm.  _Inform. Retrieval_ , 4(2):133–151, 2001.
  * Golub and Van Loan (2012) G. Golub and C. F. Van Loan.  _Matrix Computations_.  Johns Hopkins University Press, Baltimore, 4th edition, 2012.
  * Goreinov et al. (1997) S. A. Goreinov, E. E. Tyrtyshnikov, and N. L. Zamarashkin.  A theory of pseudoskeleton approximations.  _Linear Algebra Appl._ , 261(1-3):1–21, 1997\.
  * Goreinov et al. (2010) S. A. Goreinov, I. V. Oseledets, D. V. Savostyanov, E. E. Tyrtyshnikov, and N. L. Zamarashkin.  How to find a good submatrix.  In _Matrix Methods: Theory, Algorithms And Applicati_ , pages 247–256. World Scientific, Singapore, 2010.
  * Gu and Eisenstat (1996) M. Gu and S. C. Eisenstat.  Efficient algorithms for computing a strong rank-revealing QR factorization.  _SIAM J. Sci. Comput._ , 17(4):848–869, 1996\.
  * Khabou et al. (2013) A. Khabou, J. W. Demmel, L. Grigori, and M. Gu.  LU factorization with panel rank revealing pivoting and its communication avoiding version.  _SIAM J. Matrix Anal. Appl._ , 34(3):1401–1429, 2013.
  * Kramer and Gorodetsky (2018) B. Kramer and A. A. Gorodetsky.  System identification method via CUR factorization of Hankel matrix.  _SIAM J. Sci. Comput._ , 40(2):A848–A866, 2018\.
  * Mahoney and Drineas (2009) M. W. Mahoney and P. Drineas.  CUR matrix decompositions for improved data analysis.  _Proc. Natl. Acad. Sci. USA_ , 106:697–702, 2009.
  * Oseledets and Tyrtyshnikov (2010) I. Oseledets and E. Tyrtyshnikov.  TT-cross approximation for multidimensional arrays.  _Linear Algebra Appl._ , 432(1):70–88, 2010.
  * Sorensen and Embree (2016) D. C. Sorensen and M. Embree.  A DEIM induced CUR factorization.  _SIAM J. Sci. Comput._ , 33(3):A1454–A1482, 2016\.
  * Stewart (1998) G. W. Stewart.  Four algorithms for the efficient computation of truncated pivoted QR approximations to a sparse matrix.  _Numer. Math._ , 83:313–323, 1998.
  * Voronin and Martinsson (2017) S. Voronin and P. G. Martinsson.  Efficient algorithms for CUR and interpolative matrix decompositions.  _Adv. Comput. Math_ , 43:495–516, 2017.


∎

e1e-mail<EMAIL_ADDRESS>11institutetext: Departamento de Física,
Universidade Federal do Maranhão, 65080-805, São Luís, Maranhão, Brazil

# Ground state of a bosonic massive charged particle in the presence of
external fields in a Gödel-type spacetime

Edilberto O. Silvae1

(Received: date / Accepted: date)

###### Abstract

The relativistic quantum dynamics of a spinless charged particle interacting
with both Aharonov-Bohm and Coulomb-type potentials in the Gödel-type
spacetime is considered. The dynamics of the system is governed by the Klein-
Gordon equation with interactions. We verify that it is possible to establish
a quantum condition between the energy of the particle and the parameter that
characterizes the vorticist of the spacetime. We rigorously analyze the ground
state of the system and determine the corresponding wavefunctions to it.

††journal: Eur. Phys. J. C

##  1 Introduction

The relativistic quantum dynamics of spinless particles in the presence of
external fields has been an object of study for many decades. The physical
properties of the systems are accessed by the solution of the Klein-Gordon
field equation with electromagnetic interactions Book.2000.Greiner ;
Book.Landau.Vol4 . The electromagnetic interactions are introduced into the
Klein-Gordon equation through the so called minimal substitution,
$p^{\mu}\rightarrow p^{\mu}-eA^{\mu}$, where $e$ is the charge and $A^{\mu}$
is the four-potential of the electromagnetic field. It is also known that
interactions can be implemented by making a modification in the mass term as
$M\rightarrow M+S\left(\mathbf{r},t\right)$, where
$S\left(\mathbf{r},t\right)$ is a scalar potential. Contrary to the minimal
substitution, the scalar interaction $S\left(\mathbf{r},t\right)$ is
independent of the charge of the spinless particle considered and,
consequently, it has the same effect on particles and antiparticles,
respectively. In the context of recent applications, such couplings have been
used to study the motion of particles in external fields in various branches
of physics. For example, in Ref. PRD.2016.93.045033 , the minimal substitution
was used to analyze the particle scattering and vacuum instability problems in
a constant electric field for the special set of stationary solutions of the
Dirac and Klein-Gordon equations. A procedure to construct generalized ladder
operators for the Klein-Gordon equation with a scalar curvature term was
obtained in Ref. PRD.2018.97.025011 . Another interesting work addresses the
study of nonlinear quantum electrostatic waves in a pseudorelativistic quantum
semiconductor plasm PRE.2015.91.043108 . In this description, the authors
considered the substitution $p^{0}\rightarrow p^{0}-eA^{0}$ and $p\rightarrow
p$, and they showed that the system is governed by the Klein-Gordon equation
for the collective wave functions of the conducting electrons and Poisson’s
equation for the electrostatic potential. It is also common to use both the
scalar and the minimum couplings to study physical models in gravitation. In
this direction, we call attention to the relativistic quantum mechanics of
particles in external fields in the cosmic string spacetime Book.2000.Vilenkin
; EPJC.2012.72.2051 ; MPLA.2018.33.1850025 ; PRD.2014.90.125014 ;
PRD.2010.82.084025 ; EPJC.2018.78.13 ; EPJC.2016.76.61 . In general
relativity, these couplings appear in the study of the quantum dynamics of
relativistic particle in the Gödel-type spacetime EPJPlus.2015.130.36 . The
Gödel spacetime is considered a usual framework for studying physical systems
in general relativity RMP.1949.21.447 . From the topological point of view,
this spacetime is known to be geodesically complete and it is a singularity-
free solution of the Einstein field equation JHEP.2012.2012.32 ;
PRD.2013.87.087503 ; PRD.2015.92.123541 ; EPJC.2017.77.289 . The interest in
this issue led some researchers to propose a generalization of this universe
by examining the conditions for spacetime homogeneity of a Riemannian manifold
with a Gödel-type metric PRD.1983.28.1251 . The study of physical systems
involving Gödel-type solutions is a problem of current interest. The
generalization of Ref. PRD.1983.28.1251  resulted in the first exact Gödel-
type solution of the Einstein equation for a rotating universe, including the
problem of causality in (Raychaudhuri-Thakurta)-homogeneous Riemannian
manifolds PRSLA.1968.304.81 . In this context, several important contributions
have been made, as can be see in Refs. JMP.1999.40.4011 ; GRG.2008.40.2115 ;
JCAP.2004.2004.012 ; EPJC.2014.74.2935 . In this paper, we study the
relativistic quantum dynamics of a spinless charged particle interacting with
both the Aharonov-Bohm and the Coulomb-type potentials in the Gödel-type
spacetime, with line element in cylindrical coordinates written as (in natural
units, $\hbar=c=G=1$)

$ds^{2}=-\left(dt+\Omega\frac{\sinh^{2}lr}{l^{2}}d\phi\right)^{2}+\frac{\sinh^{2}2lr}{4l^{2}}d\phi^{2}+dr^{2}+dz^{2},$
(1)

with $0\leq r<\infty,0\leq\phi\leq 2\pi$ and $-\infty<(z,t)<\infty$. The
metric (1) is characterized by the parameters $(l,\Omega)$, both with the
dimensions of the inverse length, with $\Omega$ representing the vorticity of
the spacetime. The physics contained in the metric (1) is related to the
problem of a charged particle in a plane in the presence of a constant
magnetic field, with $\Omega$ playing the role of such a magnetic field
JCAP.2004.2004.012 . We address the particular case where the limit
$l\rightarrow 0$ is imposed for the metric (1), which reveals that the
resulting metric has the same geometry as that of the so called Som-
Raychaudhuri spacetime PRSLA.1968.304.81 . In our approach, the relevant
equation is the Klein-Gordon equation with interactions. We solve this
equation and determine the energy eigenvalues and the corresponding wave
functions. We evaluate the expression for the energy spectrum of the particle
and establish conditions for its validity. Finally, we consider the
fundamental state of the system and study it in detail.

##  2 The Klein-Gordon equation

After imposing the limit $l\rightarrow 0$ in (1), the metric is simplified and
it takes the form

$ds^{2}=-\left(dt+\Omega
r^{2}d\varphi\right)^{2}+r^{2}d\varphi^{2}+dr^{2}+dz^{2}.$ (2)

The relativistic quantum dynamics of a spinless charged particle of mass $M$
in an electromagnetic field is described by the Klein-Gordon equation with
minimal coupling, which is written in the covariant form as

$\left[\frac{1}{\sqrt{-g}}D_{\mu}\left(\sqrt{-g}g^{\mu\nu}D_{\nu}\right)-M^{2}\right]\psi\left(t,r,\varphi,z\right)=0,$
(3)

with $D_{\mu}=\partial_{\mu}+ieA_{\mu}$, where $e$ is the electric charge,
$A_{\mu}$ denotes the quadrivector potential associated to the electromagnetic
field and $g=\det\left(g^{\mu\nu}\right)$ is the determinant of the metric
tensor in the spacetime of the line element (2). We can exploit the
translational symmetry by considering the Klein-Gordon equation for the case
where $\partial/\partial z=0$ and $A_{z}=0$. The Aharonov-Bohm potential is
given by $eA_{\varphi}=-\Phi/r$, where $\Phi=\Phi_{B}/\left(2\pi/e\right)$ is
the magnetic flux parameter. In order to include the Coulomb potential, we put
$eA_{0}=V\left(r\right)=\xi/r$, where $\xi=\pm e^{2}$ (the choice of sign
corresponds to either repulsion or attraction, respectively). Equation (3)
assumes the form

$\displaystyle\Bigg{\\{}\frac{\partial^{2}}{\partial
r^{2}}+\frac{1}{r}\frac{\partial}{\partial r}$
$\displaystyle+\frac{1}{r^{2}}\left[\frac{\partial}{\partial\varphi}+i\Phi-\Omega
r^{2}\left(\frac{\partial}{\partial t}+\frac{i\xi}{r}\right)\right]^{2}$
$\displaystyle-\left(\frac{\partial}{\partial
t}+\frac{i\xi}{r}\right)^{2}-M^{2}\Bigg{\\}}\psi\left(t,r,\varphi\right)=0,$
(4)

It is known that the partial wave functions for the Aharonov-Bohm problem can
be obtained exactly even though a Coulomb potential is included into the
equation of motion JPA.1992.25.L183 . Then, for the wave function
$\psi\left(t,r,\varphi\right)$, we make the following ansatz:

$\psi\left(r,\varphi,t\right)=e^{-iEt+im\varphi}R\left(r\right),$ (5)

with $m=0,\pm 1,\pm 2,\ldots$. This leads to the radial equation

$\left[\frac{d^{2}}{dr^{2}}+\frac{1}{r}\frac{d}{dr}-\frac{j^{2}}{r^{2}}-\frac{a}{r}+br-\omega^{2}r^{2}+k^{2}\right]R\left(r\right)=0,$
(6)

where $j^{2}=\left(m+\Phi\right)^{2}-\xi^{2}$,
$a=2\xi\left[E-\Omega\left(m+\Phi\right)\right]$, $b=2E\xi\Omega^{2}$,
$\omega^{2}=E^{2}\Omega^{2}$ and
$k^{2}=E^{2}-M^{2}-\xi^{2}\Omega^{2}-2\Phi\Omega E-2m\Omega E$. In order to
solve the differential equation (6), it is convenient to define a new radial
function,

$R\left(r\right)=\frac{f\left(r\right)}{\sqrt{r}},$ (7)

which, together with the variable change $x=\sqrt{\omega}r$, leads to the
following Schrödinger-type equation:

$\left[\frac{d^{2}}{dx^{2}}-\frac{\left(j^{2}-\frac{1}{4}\right)}{x^{2}}-\frac{a^{\prime}}{x}+b^{\prime}x-x^{2}+\frac{k^{2}}{\omega}\right]f\left(x\right)=0,$
(8)

where $a^{\prime}={a/}\sqrt{\omega}$ and $b^{\prime}=b/\sqrt{\omega^{3}}$.
Now, let us study the asymptotic limits of Eq. (8). This analysis reveals that
it is possible to obtain normalizable eigenfunctions by writing
$f\left(x\right)$ as

$f\left(x\right)={x}^{j+\frac{1}{2}}e{{}^{-\frac{1}{2}\,x\left(x+b^{\prime}\right)}y}\left(x\right),$
(9)

where ${y}\left(x\right)$ is a new function to be determined. Substituting Eq.
(9) into the Eq. (8), we get

$\displaystyle x\frac{d^{2}y}{d{x}^{2}}$ $\displaystyle+\left(1+\alpha-\beta
x-2x^{2}\right)\frac{dy}{d{x}}$
$\displaystyle+\left\\{\left(\gamma-\alpha-2\right)x-\frac{1}{2}\left[\delta+\left(1+\alpha\right)\beta\right]\right\\}y=0,$
(10)

where $\alpha=\,{2{j}}$, $\beta=-b^{\prime}$,
$\gamma=\frac{k^{2}}{\omega}+\frac{1}{4}b^{\prime 2}$ and
$\delta=2a^{\prime}$. Equation (10) is the biconfluent Heun equation in its
canonical form ronveaux1995heun ; JCAM.37.161.1991 . This equation has a
regular singularity at $x=0$, and an irregular singularity at $\infty$ of rank
$2$. The use of the Frobenius method allows us to determine a series solution
for (10). Upon writing

$N\left(\alpha,\beta,\gamma,\delta;x\right)=\sum\limits_{\upsilon=0}^{\infty}\frac{\mathcal{A}_{\upsilon}\left(\alpha,\beta,\gamma,\delta\right)}{\left(1+\alpha\right)_{\upsilon}}\frac{x^{\upsilon}}{\upsilon!},$
(11)

where

$\left(1+\alpha\right)_{\upsilon}=\frac{\Gamma\left(\upsilon+\alpha+1\right)}{\Gamma\left(\alpha+1\right)},$
(12)

we obtain the three-term recurrence relation,

$\displaystyle\mathcal{A}_{\upsilon+2}$
$\displaystyle=\left\\{\left(\upsilon+1\right)\beta+\frac{1}{2}\left[\delta+\beta\left(1+\alpha\right)\right]\right\\}\mathcal{A}_{\upsilon+1}$
$\displaystyle-\left(\upsilon+1\right)\left(\upsilon+1+\alpha\right)\left[\gamma-2-\alpha-2\upsilon\right]\mathcal{A}_{\upsilon},$
(13)

with

$\mathcal{A}_{0}=1,$ (14)
$\mathcal{A}_{1}=\frac{1}{2}\left[\delta+\beta\left(1+\alpha\right)\right].$
(15)

It can be verified directly in the recursion relation (13) that, if
$\mathcal{A}_{\upsilon+2}=0$, the function
$N\left(\alpha,\beta,\gamma,\delta;x\right)$ becomes a polynomial of degree
$n$. This can be accomplished if and only if the following conditions are
imposed:

$\gamma-\alpha-2=2n,$ (16)

and

$\mathcal{A}_{n+1}=0,$ (17)

where $n$ is a positive integer. In this case, the $\left(n+1\right)$th
coefficient in the series expansion is a polynomial of degree $n$ in $\delta$.
When $\delta$ is a root of this polynomial, the $\left(n+1\right)$th and the
subsequent coefficients are identically null and the series truncates,
resulting in a polynomial form of degree $n$ for series solution
$N\left(\alpha,\beta,\gamma,\delta;x\right)$. The general solution of (10) is
given by

$y\left(x\right)=c_{1}\,N\left(\alpha,\beta,\gamma,\delta;x\right)+c_{2}\mathrm{\,}x^{-\alpha}N\left(-\alpha,\beta,\gamma,\delta;x\right),$
(18)

where $N\left(\alpha,\beta,\gamma,\delta;x\right)$ are the Heun functions and
$c_{1}$ and $c_{2}$ are normalization constants. The solution (9) becomes

$\displaystyle
f\left(x\right)=c_{1}\,{x}^{j+\frac{1}{2}}e{{}^{-\frac{1}{2}x\left(x+b^{\prime}\right)}}N\left({2{j}},-b^{\prime},\frac{k^{2}}{\omega}+\frac{1}{4}b^{\prime
2},2a^{\prime};x\right)$
$\displaystyle+c_{2}\,{x}^{-j+\frac{1}{2}}e{{}^{-\frac{1}{2}x\left(x+b^{\prime}\right)}}N\left(-{2{j}},-b^{\prime},\frac{k^{2}}{\omega}+\frac{1}{4}b^{\prime
2},2a^{\prime};x\right).$ (19)

By using the conditions (16) and (17), we can determine an expression for the
energy eigenvalues of the particle. From condition (16), we obtain

$E_{nm}^{2}-2\Omega\left({n}+j+m+\Phi+{1}\right)E_{nm}-M^{2}=0{.}$ (20)

By solving (20) for $E_{nm}$, we find

$\displaystyle E_{nm}$ $\displaystyle=\left({n}+j+m+\Phi+{1}\right)\Omega$
$\displaystyle\pm\sqrt{\Omega^{2}\left({n}+j+m+\Phi+{1}\right)^{2}+M^{2}}.$
(21)

As we can see, the eigenvalues in Eq. (21) are given in terms of all the
physical parameters of the problem, that is, all the physical quantities
present in Klein-Gordon equation (4). However, these eigenvalues do not
constitute the spectrum of the system. In order to make the Eq. (21) to
represent it, we need to implement the condition (17). As pointed out in Ref.
AoP.2014.347.130 , such condition (17) gives rise to a quantum condition
between the energy eigenvalues and the frequency of the model. In our case,
the frequency is characterized by the vorticity $\Omega$ of the spacetime. On
the other hand, in some particular models found in the literature, the authors
suggest that this quantum condition has been used to "restore" some physical
parameter of the problem which do not appear in the expression for the energy
derived in Eq. (16) (and which is present in the equation of motion) (see for
example the Ref. EPJC.2012.72.2051 ). In our case, the requirement (17)
appears as a necessary condition to ensure the existence of a set of energy
eigenfunctions of the system. In fact, it is well-known of quantum mechanics
that for a set of eigenvalues $E_{n}$ of a given hamiltonian operator $H_{n}$
, we have the set of energy eigenfunctions $\psi_{n}$ corresponding to each
eigenvalue. The use of the condition (17) to determine quantum conditions to
restore physical parameters, a priori, is not a general rule. Each physical
model should either require or not the establishment of such a condition.

Here, the requirement (17) also provides a quantum condition between energy
and frequency as in Ref. AoP.2014.347.130 . We would like to emphasize that,
for $n\geqslant 1$, the equations resulting from conditions (16) and (17) are
more complicated and they generally do not provide physically acceptable
energy eigenvalues. For this reason, we investigate only the ground state of
the particle. Thus, the $n=0$ state provides two quantum conditions, given by

$\displaystyle\left[\Omega\allowbreak_{0m}\right]_{1}$
$\displaystyle=\frac{2E_{0m}}{2\left(j+m+\Phi\right)+1},$ (22)
$\displaystyle\left[\Omega\allowbreak_{0m}\right]_{2}$
$\displaystyle=\frac{2E_{0m}}{2\left(-j+m+\Phi\right)-1}.$ (23)

The final expressions for the energies are obtained after substituting (22)
and (23) in (20):

$\displaystyle\left[E_{0m}\right]_{1}$ $\displaystyle=$
$\displaystyle\pm\sqrt{-\frac{2\left(j+m+\Phi\right)+1}{2\left(j+m+\Phi\right)+3}}\,M,$
(24) $\displaystyle\left[E_{0m}\right]_{2}$ $\displaystyle=$
$\displaystyle\pm{\sqrt{-\frac{2\left(-\,j+m+\Phi\right)-1}{2\left(3j+m+\,\Phi\right)+5}}\,M}.$
(25)

Thus, there exist solutions both for positive as well as for negative
energies, respectively, with the requirement that

$\displaystyle\frac{2\left(j+m+\Phi\right)+1}{2j+2m+2\Phi+3}<0,$ (26)
$\displaystyle\frac{2\left(-\,j+m+\Phi\right)-1}{2\left(3j+m+\Phi\right)+5}<0,$
(27)

to ensure that $\left[E_{0m}\right]_{1,2}$ are real. Evidently, for the ground
state, the regular Heun function, Eq. (11), becomes

$N\left({2{j}},-b^{\prime},\frac{k^{2}}{\omega}+\frac{1}{4}b^{\prime
2},2a^{\prime};x\right)=1.$ (28)

In this case, the unnormalized bound state wave function for our problem
corresponding to the regular solution at the origin is given by

$f\left(x\right)={x}^{j+\frac{1}{2}}e{{}^{-\frac{1}{2}x\left(x+b^{\prime}\right)}}.$
(29)

The conditions (26) and (27) reveal that the energies are physically
acceptable only for appropriate values of $m$, $\Phi$ and $\xi$. It is
important to note that, at least for the ground state, when $M=0$ in Eqs. (24)
and (25), the energies are identically zero. This fact that does not occur for
the spectrum of Ref. EPJC.2014.74.2935 , where $M=0$ leads immediately to the
analogue Landau levels.

Figure 1: Profile of the ground state (Eq. (24)) as a function of the magnetic
flux $\Phi$ for $M=1$ and $\xi=1$. The energy is symmetric in relation to the
quantum number $m$. The behavior for small variations of magnetic flux for
$m=-3$ is shown by the inset. Figure 2: Profile of the ground state (Eq. (25))
as a function of the magnetic flux $\Phi$ for $M=1$ and $\xi=1$. Clearly, we
have regions in which energies are not defined for some values of $m$,
existing only in small intervals.

The energy (24) is plotted as a function of the magnetic flux parameter $\Phi$
in Fig. 1 for some values from $m$. We can see that both the particle and the
antiparticle energy levels are members of the energy spectrum. Another
explicit evidence in Fig. 1 is that the energies are symmetrical about
$\left[E_{0m}\right]_{1}=0$ and, since the positive and the negative energies
never intercept each other, we can see that there is no channel for
spontaneous particle-antiparticle creation. Furthermore, we also see the
presence of regions where energies are not allowed. The appearance of such
regions is justified by the imposition of the quantum conditions (26) and
(27).

Figure 3: Ground state $\left[E_{0m}\right]_{1}$ (blue dot) and
$\left[E_{0m}\right]_{2}$ (dot red) as a function of the quantum number $m$
for $M=1$, $\xi=1$ and $\Phi=0$ (Fig. (a)) and $\Phi=4.8$ (Fig. (b)). Figure
4: Ground state $\left[E_{0m}\right]_{1,2}$ as a function of the quantum
number $m$ for $M=1$ and $\xi=1$.

This feature also reveals that the energy is bounded in the finite interval
$\left[E_{0m}\right]_{1,2}=\pm 1$. For positive values of $m$, including
$m=0$, there is a magnetic flux inversion. We can also see that for small
variations of magnetic field the energy is abruptly modified (see inset in
Fig. 1). This occurs when the magnetic flux takes on integer values. This
characteristic is manifested when the magnetic flux takes on integer values.

Figure 5: Ground state (Eq. (30)) as a function of the quantum number $m$ for
$M=1$ and $\Omega=1$.

The profile of the energy (25) is shown in Fig. 2 for the same values of $m$
considered in Fig. 1 and, it also exhibits the same characteristics. The only
difference is that in Fig. 2 appears new regions of bounded energies. We also
investigate the ground state of the system by assuming a constant magnetic
field in comparison to the case when it is identically null (Fig. 3). For this
case, the behavior of it is similar to which is observed in Figs. 1 and 2. To
show the different effects of the Coulomb-type potential and Aharonov-Bohm
flux on the ground state of the particle, we plot it in Fig. 4 for various
values of $m$ and large $\Phi$.

In particular, when $\xi=\Phi=0$, the condition (16) assumes the form
$\gamma-\alpha-2=4\bar{n}$, where $\bar{n}$ is a positive integer. The
energies in this case are given by

$\displaystyle E_{\bar{n},m}=(2\bar{n}+{\left|m\right|}$
$\displaystyle+m+1)\Omega$
$\displaystyle\pm\sqrt{\Omega^{2}\left(2\bar{n}+{\left|m\right|}+m+1\right)^{2}+M^{2}},$
(30)

which is the same energy of the particle obtained in EPJC.2014.74.2935 . It
should also be noted that when $M=0$ in Eq. (30), we obtain the Landau analog
quantization. When we compare the ground state ($\bar{n}=0$) of the result in
Eq. (30) (see Fig. (5)) with the energy levels of the Fig. 3(b), which shows a
symmetry between the energy of the particle and the anti-particle, we observe
that such symmetry is no longer present. For negative values of $m$, the
difference between energy levels is the same, whereas for positive values, the
energy levels of the particle increase while that of the antiparticle is
approximately constant.

## Conclusion

We have presented a detailed study on the ground state of a Klein-Gordon
particle in the presence of external fields in a Gödel-type spacetime. Through
minimal substitution, we implemented both the Aharonov-Bohm potential and a
Coulomb-type potential. The radial equation of motion was derived by means of
an appropriate ansatz, and it was shown that such equation is expressed in
terms of the biconfluent Heun differential equation. The expression for the
ground state energy was obtained using the relations (16) and (17), which
reveals the possibility of establishing a quantum condition between the energy
of the particle and the parameter that characterizes the vorticity $\Omega$ of
the spacetime. As we have mentioned, expressions for such energies are
obtained without major complications only for the ground state ($n=0$). We
specialized to the $n=0$ case and determined the expression for the energy and
wave function of the particle. We evaluated the expressions obtained and
established validity criteria for them (Eqs. (26) and (27)). The physical
implications imposed by it reveals that the energy is bounded in the finite
interval $\left[E_{0m}\right]_{1,2}=\pm 1$. It was also observed an inversion
of magnetic flux for certain values of quantum number $m$. Finally, it was
found that the fundamental state of both the particle and the antiparticle
belongs to the same spectrum, and there is no channel that allows the
spontaneous creation of particles.

## Acknowledgments

The author would like to thank Cleverson Filgueiras and Fabiano M. Andrade for
their critical reading of the manuscript. This work was supported by CNPq,
Grants No. 427214/2016-5 (Universal), No. 303774/2016-9 (PQ), and FAPEMA,
Grants No. 01852/14 (PRONEM) and 01202/16 (Universal).

## References

  * (1) W. Greiner, _Relativistic Quantum Mechanics. Wave Equations_ (Springer Science Business Media, 2000).  DOI 10.1007/978-3-662-04275-5
  * (2) V. Berestetskii, L. Pitaevskii, E. Lifshitz, _Quantum Electrodynamics_.  v. 4 (Elsevier Science, 2012)
  * (3) S.P. Gavrilov, D.M. Gitman, Phys. Rev. D 93, 045033 (2016).  DOI 10.1103/PhysRevD.93.045033
  * (4) W. Mück, Phys. Rev. D 97, 025011 (2018).  DOI 10.1103/PhysRevD.97.025011
  * (5) Y. Wang, X. Wang, X. Jiang, Phys. Rev. E 91, 043108 (2015).  DOI 10.1103/PhysRevE.91.043108
  * (6) A. Vilenkin, E.P.S. Shellard, _Cosmic Strings and Other Topological Defects_ (Cambridge University Pres, Canbridge, 2000)
  * (7) E.R. Figueiredo Medeiros, E.R. Bezerra de Mello, Eur. Phys. J. C 72(6), 2051 (2012).  DOI 10.1140/epjc/s10052-012-2051-9
  * (8) B.Q. Wang, Z.W. Long, C.Y. Long, S.R. Wu, Modern Physics Letters A 33(04), 1850025 (2018).  DOI 10.1142/S0217732318500256
  * (9) D. Chowdhury, B. Basu, Phys. Rev. D 90, 125014 (2014).  DOI 10.1103/PhysRevD.90.125014
  * (10) K. Bakke, C. Furtado, Phys. Rev. D 82(8), 084025 (2010).  DOI 10.1103/PhysRevD.82.084025
  * (11) L.C.N. Santos, C.C. Barros, The European Physical Journal C 78(1), 13 (2018).  DOI 10.1140/epjc/s10052-017-5476-3
  * (12) L.B. Castro, The European Physical Journal C 76(2), 61 (2016).  DOI 10.1140/epjc/s10052-016-3904-4
  * (13) Z. Wang, Z.w. Long, C.y. Long, M.l. Wu, The European Physical Journal Plus 130(3), 36 (2015).  DOI 10.1140/epjp/i2015-15036-2
  * (14) K. Gödel, Rev. Mod. Phys. 21, 447 (1949).  DOI 10.1103/RevModPhys.21.447
  * (15) C.M. Brown, O. DeWolfe, Journal of High Energy Physics 2012(1), 32 (2012).  DOI 10.1007/JHEP01(2012)032
  * (16) L. Herrera, J. Ibáñez, A. Di Prisco, Phys. Rev. D 87, 087503 (2013).  DOI 10.1103/PhysRevD.87.087503
  * (17) S. Khodabakhshi, A. Shojai, Phys. Rev. D 92, 123541 (2015).  DOI 10.1103/PhysRevD.92.123541
  * (18) S.L. Li, X.H. Feng, H. Wei, H. Lü, The European Physical Journal C 77(5), 289 (2017).  DOI 10.1140/epjc/s10052-017-4856-z
  * (19) M.J. Rebouças, J. Tiomno, Phys. Rev. D 28, 1251 (1983).  DOI 10.1103/PhysRevD.28.1251
  * (20) M.M. Som, A.K. Raychaudhuri, Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 304(1476), 81 (1968).  DOI 10.1098/rspa.1968.0073
  * (21) H.L. Carrion, M.J. Rebouças, A.F.F. Teixeira, Journal of Mathematical Physics 40(8), 4011 (1999).  DOI 10.1063/1.532939
  * (22) S. Das, J. Gegenberg, General Relativity and Gravitation 40(10), 2115 (2008).  DOI 10.1007/s10714-008-0619-3
  * (23) N. Drukker, B. Fiol, J. Simón, Journal of Cosmology and Astroparticle Physics 2004(10), 012 (2004)
  * (24) J. Carvalho, A.M. de M. Carvalho, C. Furtado, The European Physical Journal C 74(6), 2935 (2014).  DOI 10.1140/epjc/s10052-014-2935-y
  * (25) J. Law, M.K. Srivastava, R.K. Bhaduri, A. Khare, Journal of Physics A: Mathematical and General 25(4), L183 (1992)
  * (26) A. Ronveaux, F. Arscott, _Heun’s Differential Equations_.  Oxford science publications (Oxford University Press, 1995)
  * (27) E. Arriola, A. Zarzo, J. Dehesa, J. Comput. Appl. Math. 37(1-3), 161 (1991).  DOI 10.1016/0377-0427(91)90114-y
  * (28) F. Caruso, J. Martins, V. Oguri, Annals of Physics 347, 130 (2014).  DOI https://doi.org/10.1016/j.aop.2014.04.023