# Multi-antenna Coded Caching at Finite-SNR:
Breaking Down the Gain Structure

MohammadJavad Salehi and Antti Tölli
This work is supported by the Academy of Finland under grants no. 346208 (6G
Flagship), 319059 (CCCWEE), and 343586 (CAMAIDE), and by the Finnish Research
Impact Foundation under the project 3D-WIDE.  Centre for Wireless
Communications, University of Oulu, 90570 Oulu, Finland
E-mail<EMAIL_ADDRESS>
###### Abstract

Multi-antenna coded caching (CC) techniques are considered viable options for
achieving higher data rates in future networks, especially for the prominent
use case of multimedia-driven applications. However, despite their
information-theoretic analyses, which are thoroughly studied in the
literature, the research on the finite-SNR performance of multi-antenna CC
techniques is not yet mature. In this paper, we try bridging this gap by
breaking down, categorizing, and studying the effect of six crucial parameters
affecting the finite-SNR performance of multi-antenna CC schemes. We also
investigate the interaction of different parameters and clarify how they could
affect the implementation complexity in terms of the necessary computation and
subpacketization. Theoretical discussions are followed and verified by
numerical analysis.

###### Index Terms:

coded caching; multi-antenna communications; multiplexing; multicasting;
beamforming

##  I Introduction

Wireless communication networks are under mounting pressure to support higher
data rates, and this trend will continue by emerging new services such as
extended reality (XR) applications [1]. The coded caching (CC) technique [2]
is considered a viable candidate to address this issue. The critical property
of CC is that it enables the data storage memory in network devices to be
efficiently used as a communication resource, especially for the prominent use
case of multimedia-driven applications [3, 4]. Interestingly, the caching gain
of CC scales with the cumulative cache size of all users in the network [2],
and can be combined with the spatial multiplexing gain of incorporating
multiple antennas at the transmitter [5, 6] or receivers [7]. However, this
promising CC gain is followed by crucial practical bottlenecks hindering its
applicability. Two important such bottlenecks are the subpacketization issue
and the optimized beamformer design (for improved performance at the finite-
SNR regime). Former stems from the fact that the original multi-antenna CC
scheme in [6] required splitting files into exponentially-growing numbers of
smaller parts, and the latter is due to the fact that designing optimized
beamformers for this scheme requires solving complex non-convex optimization
problems [8].

While the subpacketization issue seems to be fundamental in single-antenna
systems [9], multi-antenna setups provide flexibility for reducing the
required subpacketization without altering the theoretical degrees-of-freedom
(DoF) gains [10, 11, 12]. The key to achieving this reduction is in using a
new CC approach where the cache-aided interference cancellation is done
_before_ decoding the data and in the signal domain [10]. Following [4], we
use the term _signal-level_ while referring to CC schemes with this new
approach and denote the traditional CC schemes (with cache-aided interference
cancellation after decoding the received signal) as _bit-level_.
Interestingly, even though signal-level CC schemes were initially developed
for subpacketization reduction, the work in [11] showed that they also allow
simplifying the optimized beamformer design by altering the underlying
multicast structure of bit-level schemes.

However, achieving the same DoF and applicability of optimized beamformers
does _not_ reveal the whole story about the finite-SNR performance of signal-
level CC schemes. For example, it was shown in [13] that with the same DoF,
moving eventually from a pure signal-level approach (smallest
subpacketization) to a pure bit-level one (largest subpacketization), the
performance is improved at the finite-SNR regime. Also, the work in [8]
revealed the fact that being DoF-optimal is not necessarily preferred in
finite-SNR as the DoF can be traded-off with the more effective beamformer
directivity parameter. Such results clarify a need to better understand all
the parameters affecting the finite-SNR performance of CC schemes and to
investigate the possibility of designing new CC schemes and beamforming
techniques to better tune the performance for any given complexity constraint.

This paper breaks down and classifies all the parameters affecting the finite-
SNR performance of multi-antenna CC schemes. Specifically, we identify six
different parameters, as shown in Figure 1, and analyze how they affect the
achievable rate of CC schemes in different SNR values. For each parameter, we
provide the definition and a brief discussion on how it can be adjusted,
followed by theoretical insights and numerical simulation results. As
discussed, this paper is not the first to identify all these parameters and
study their effects. However, it is the first to gather, organize, and study
the interaction of all such performance-affecting parameters. We also
carefully discuss how each parameter could affect the implementation
complexity and provide guidelines for selecting the appropriate scheme
following the available transmission power and acceptable complexity.

Figure 1: Parameters affecting the finite-SNR performance of MISO CC schemes

In this paper, we have used the following notations. Bold-face lower- and
upper-case letters denote vectors and matrices, respectively. Calligraphic
letters are used to represent sets. For an integer $K$, $[K]$ is the set of
numbers $\\{1,2,\cdots,K\\}$. The symbol $\oplus$ denotes the bit-wise XOR
operation. For two sets ${\mathcal{A}}$ and ${\mathcal{B}}$,
${\mathcal{A}}\backslash{\mathcal{B}}$ is the set of elements of
${\mathcal{A}}$ that are not in ${\mathcal{B}}$. $|{\mathcal{A}}|$ is the
number of elements in ${\mathcal{A}}$. Set brackets and element separators are
sometimes dropped for notational simplicity.

##  II System Model

We consider a multi-input single-output (MISO) cache-aided communication
setup. A single, $L$-antenna transmitter with spatial multiplexing gain of
$\alpha\leq L$ communicates with $K$ single-antenna users over a shared
wireless link. Every user has a cache memory of size $MF$ bits and requests
files from a library ${\mathcal{F}}$ of $N$ files, each with the size of $F$
bits. The CC gain is defined as $t\equiv\frac{KM}{N}$. For notational
simplicity, we use a normalized data unit and drop $F$ in our subsequent
notations. We also use CC-$(K,t,L,\alpha)$ to refer to this setup.

The system operation consists of two phases, content placement and delivery.
In the content placement phase, which is done at the low network traffic time,
cache memories of the users are filled up by data chunks of files in
${\mathcal{F}}$. Then, in the delivery phase, after each user $k\in[K]$
reveals its requested file $W(k)\in{\mathcal{F}}$, the server creates and
transmits a set of codewords such that all the users can decode their
requested files using their cache contents and the received signals. The goal
is to design content placement and delivery phases such that the symmetric
rate, defined as the amount of data delivered to all the users over a defined
time frame, is maximized. In this paper, we consider only CC schemes with
uncoded data placement and single-shot data delivery, for which the MISO-CC
scheme in [6] is shown to be DoF-optimal [14].

Data delivery is done using a set of transmission vectors, which are sent, for
example, in consecutive time slots. A transmission vector
${\mathbf{x}}_{{\mathcal{S}}}$ delivers (parts of) the requested data to every
user in the subset ${\mathcal{S}}\subseteq[K]$ of users, where
$|{\mathcal{S}}|=t+\alpha$. This is done through the parallel transmission of
codewords $X_{{\mathcal{T}}}$, where ${\mathcal{T}}\subseteq{\mathcal{S}}$.
Every codeword $X_{{\mathcal{T}}}$ includes useful data for every user
$k\in{\mathcal{T}}$ and is precoded with a beamforming vector
${\mathbf{w}}_{{\mathcal{R}}({\mathcal{T}})}$ that nulls-out or suppresses the
interference of $X_{{\mathcal{T}}}$ at every user
$k\in{\mathcal{R}}({\mathcal{T}})$. Using $\bar{{\mathcal{T}}}({\mathcal{S}})$
to denote the set of subsets ${\mathcal{T}}\subseteq{\mathcal{S}}$ for which a
codeword $X_{{\mathcal{T}}}$ is built, we can write

${\mathbf{x}}_{{\mathcal{S}}}=\sum_{{\mathcal{T}}\in\bar{{\mathcal{T}}}({\mathcal{S}})}X_{{\mathcal{T}}}{\mathbf{w}}_{{\mathcal{R}}({\mathcal{T}})}.$
(1)

Then, after the transmission of ${\mathbf{x}}_{{\mathcal{S}}}$, the received
signal at user $k\in{\mathcal{S}}$ can be modeled as

$y_{{\mathcal{S}}}(k)=\sum_{{\mathcal{T}}\in\bar{{\mathcal{T}}}({\mathcal{S}})}X_{{\mathcal{T}}}{{\mathbf{h}}_{{\mathcal{S}}}(k)}^{H}{\mathbf{w}}_{{\mathcal{R}}({\mathcal{T}})}+z_{{\mathcal{S}}}(k),$
(2)

where ${\mathbf{h}}_{{\mathcal{S}}}(k)\in\mathbb{C}^{L\times 1}$ and
$z_{{\mathcal{S}}}(k)$ represent the channel vector and noise at user $k$
during the transmission of ${\mathbf{x}}_{{\mathcal{S}}}$. In this paper, we
use assume ${\mathbf{h}}_{{\mathcal{S}}}(k)$ does not change during the
transmission of ${\mathbf{x}}_{{\mathcal{S}}}$, and full channel state
information (CSI) is available at the transmitter. We also use
${\mathbf{h}}_{k}\equiv{\mathbf{h}}_{{\mathcal{S}}}(k)$ and $z_{k}\equiv
z_{{\mathcal{S}}}(k)$ for notational simplicity.

Let us denote the maximum time required for every user $k$ in ${\mathcal{S}}$
to decode $y_{{\mathcal{S}}}(k)$ as $T_{{\mathcal{S}}}$. Then, the total
delivery time is $T=\sum_{{\mathcal{S}}}T_{{\mathcal{S}}}$, and the symmetric
channel rate, used as the comparison metric here, can be defined as

$R_{\mathrm{sym}}=\frac{K-t}{T}.$ (3)

Note that this metric reflects how well the wireless communication channel is
used and excludes the effect of the local caching gain. This is unlike the
total symmetric rate $\frac{K}{T}$ widely used in the literature [6]. Of
course, the two metrics can be converted into each other with a simple
coefficient.

##  III MISO-stemmed parameters

As shown in Figure 1, we study three MISO-stemmed parameters affecting the
finite-SNR performance of coded caching schemes. For this study, we choose the
MISO scheme in [6] as the baseline. This is because this scheme enables the
maximum performance for all CC-stemmed parameters and hence, provides a fair
structure for capturing the effect of the MISO-stemmed parameters considered
here.

With the MISO scheme in [6] as the baseline, during the placement phase, we
split every file $W\in{\mathcal{F}}$ into $\binom{K}{t}$ subfiles
$W_{{\mathcal{U}}}$, where ${\mathcal{U}}\subseteq[K]$ and
$|{\mathcal{U}}|=t$. Then, at the cache memory of user $k\in[K]$, we store
$W_{{\mathcal{U}}}$ for every $W\in{\mathcal{F}}$ and ${\mathcal{U}}\ni k$.
Subsequently, during the delivery phase, we create a separate transmission
vector ${\mathbf{x}}_{{\mathcal{S}}}$ for _every_ subset of users
${\mathcal{S}}\subseteq[K]$ with size $|{\mathcal{S}}|=t+\alpha$. Given one
such set ${\mathcal{S}}$, the transmission vector
${\mathbf{x}}_{{\mathcal{S}}}$ includes a codeword for _every_ user subset of
${\mathcal{S}}$ with $t+1$ users, i.e.,
$\bar{{\mathcal{T}}}({\mathcal{S}})=\\{{\mathcal{T}}\subseteq{\mathcal{S}}\;\mid\;|{\mathcal{T}}|=t+1\\}$.111It
should be noted that during the delivery phase, we should further split every
subfile $W_{{\mathcal{U}}}$ into $Q=\binom{K-t-1}{\alpha-1}$ smaller
subpackets $W_{{\mathcal{U}}}^{q}$, $q\in[Q]$, to ensure new data chunks are
delivered during every transmission. In this paper, we have ignored noting
this second-level subpacketization as it does not affect the discussions. The
effects are considered in numerical simulations, though.

For better clarification, let us consider an example network of
CC-$(6,2,4,\alpha)$. During the placement phase, we split each file in
${\mathcal{F}}$ into $\binom{6}{2}=15$ subfiles and store five of them in the
cache memory of each user. For example, a file $A\in{\mathcal{F}}$ is split
into subfiles $A_{12}$, $A_{13}$, $\cdots$, $A_{56}$, from which $A_{12}$,
$A_{13}$, $A_{14}$, $A_{15}$, and $A_{16}$ (i.e., all subfiles
$A_{{\mathcal{U}}}$ for which $1\in{\mathcal{U}}$) are stored in the cache
memory of user 1. Now, we review the delivery phase, assuming $\alpha=1,4$.
For notational simplicity, let us denote the files requested by users 1-6 as
$A,B,\cdots,F$, respectively.

$\bullet$ $\alpha=1$. In this case, we have to create
$\binom{K}{t+\alpha}=\binom{6}{3}=20$ transmission vectors
${\mathbf{x}}_{{\mathcal{S}}}$. We also have
$\bar{{\mathcal{T}}}({\mathcal{S}})=\\{{\mathcal{S}}\\}$, i.e., every vector
${\mathbf{x}}_{{\mathcal{S}}}$ includes only a single codeword
$X_{{\mathcal{T}}}$. As a result, there is no need to null-out (or suppress)
the interference from $X_{{\mathcal{T}}}$ at other users, and
${\mathcal{R}}({\mathcal{T}})=\varnothing$. For example, considering
${\mathcal{S}}=\\{1,2,3\\}$, the transmission vector ${\mathbf{x}}_{123}$ is
built as

$\displaystyle{\mathbf{x}}_{123}=X_{123}{\mathbf{w}}_{\varnothing},$ (4)

where the codeword $X_{123}=A_{23}\oplus B_{13}\oplus C_{12}$ includes data
for all users 1, 2, and 3.

Now, let us consider the decoding process at user 1 after transmitting
${\mathbf{x}}_{123}$. Using (2), this user receives

$y_{123}(1)=X_{123}{\mathbf{h}}_{1}^{T}{\mathbf{w}}_{\varnothing}+z_{1},$ (5)

and can extract $X_{123}$ with a maximum rate of

$r_{123}(1)\leq\log\big{(}1+\frac{|{\mathbf{h}}_{1}^{T}{\mathbf{w}}_{\varnothing}|}{N_{0}}\big{)}.$
(6)

Then, user 1 has to use its cached contents to remove unwanted data terms
$B_{23}$ and $C_{12}$ from $X_{123}$. Similarly, users 2 and 3 can also get
their requested data terms.

$\bullet$ $\alpha=4$. In this case, we create only a single transmission
vector ${\mathbf{x}}_{1\cdots 6}$. However, this transmission vector is
comprised of $\binom{6}{3}=20$ codewords, where each codeword delivers data to
three users and its interference is suppressed by beamforming at the other
three users. In other words,

$\displaystyle{\mathbf{x}}_{1\cdots
6}=X_{123}{\mathbf{w}}_{456}+X_{124}{\mathbf{w}}_{356}+\cdots+X_{456}{\mathbf{w}}_{123}.$
(7)

The codewords are similar to the $\alpha=1$ case. For example, we have
$X_{123}=A_{23}\oplus B_{13}\oplus C_{12}$, $X_{124}=A_{24}\oplus B_{14}\oplus
D_{12}$, and $X_{456}=D_{56}\oplus E_{46}\oplus F_{45}$.

Now, let us consider the decoding process at user 1. The received signal at
this user is

$\displaystyle y_{1\cdots
6}(1)=X_{123}{\mathbf{h}}_{1}^{T}{\mathbf{w}}_{456}+\cdots+X_{456}{\mathbf{h}}_{1}^{T}{\mathbf{w}}_{123}+z_{1}.$
(8)

Unlike the $\alpha=1$ case, user 1 now has to decode its requested data from a
multiple access channel (MAC) of size 10, using a more complex successive
interference cancellation (SIC) receiver, as there exist $\binom{5}{2}=10$
codewords $X_{{\mathcal{T}}}$ for which $1\in{\mathcal{T}}$. The maximum
decoding rate in this case is determined by a rate region. Assume
$\hat{{\mathcal{V}}}(1)$ and $\bar{{\mathcal{V}}}(1)$ include the user subsets
${\mathcal{T}}\subseteq{\mathcal{S}}$ of size $t+1=3$, that include and not
include user 1, respectively. In other words

$\displaystyle\hat{{\mathcal{V}}}(1)$
$\displaystyle=\\{123,124,\cdots,156\\},$ (9)
$\displaystyle\bar{{\mathcal{V}}}(1)$
$\displaystyle=\\{234,235,\cdots,456\\}.$ (10)

Then, the SINR term $\gamma_{\hat{{\mathcal{V}}}}$ for a desired term
$A_{\hat{{\mathcal{V}}}}$, $\hat{{\mathcal{V}}}\in\hat{{\mathcal{V}}}(1)$, is
calculated as

$\gamma_{\hat{{\mathcal{V}}}}=\frac{|{\mathbf{h}}_{1}^{T}{\mathbf{w}}_{[6]\backslash\hat{{\mathcal{V}}}}|^{2}}{\sum_{\bar{{\mathcal{V}}}\in\bar{{\mathcal{V}}}(1)}|{\mathbf{h}}_{1}^{T}{\mathbf{w}}_{[6]\backslash\bar{{\mathcal{V}}}}|^{2}+N_{0}},$
(11)

and the decoding rate $r_{1\cdots 6}(1)$ satisfies

$r_{1\cdots
6}(1)\leq\\!\frac{1}{|{\mathcal{W}}|}\log\\!\big{(}1\\!+\\!\\!\\!\sum_{\hat{{\mathcal{V}}}\in{\mathcal{W}}}\gamma_{\hat{{\mathcal{V}}}}\big{)},\;\forall{\mathcal{W}}\subseteq\hat{{\mathcal{V}}}(1),|{\mathcal{W}}|>0.$
(12)

Finding optimal beamformers maximizing the rate in this case is difficult as
it requires solving a non-convex optimization problem with the number of
constraints growing exponentially with the MAC size [8]. We will discuss the
effect of different beamforming strategies in Section III-C.

###  III-A Spatial multiplexing effect

Spatial multiplexing gain $\alpha\leq L$ determines the number of users at
which we can null-out or suppress the interference caused by each term. As a
larger $\alpha$ means serving more users in parallel (and hence, a larger
DoF), most works in the literature have assumed $\alpha=L$. However, as
discussed in [8], the DoF is not the best metric at the finite-SNR regime;
setting $\alpha<L$ results in better performance due to an increased
_beamformer directivity_ gain.

In this section, to have a fair analysis of the effect of the spatial
multiplexing gain, we set $\alpha=L$ to remove the effect of beamformer
directivity. We perform numerical simulations for CC-$(6,2,L,L)$ setup, where
$L\in\\{1,2,3,4\\}$. Simulation results are shown in Figure 2 and Table I. In
all simulations, we use optimized beamformers for maximizing performance.
Highlights of the results are:

  1. 1.

If $\alpha=L$, a larger spatial multiplexing gain always improves the
performance. This is because the DoF is increased, and beamformer directivity
has no effect;

  2. 2.

Performance improvement (over the $\alpha=L=1$ case) is more prominent in the
finite-SNR regime. This is because the $\alpha=L=1$ case represents omni-
casting data in all directions, which is very inefficient in finite-SNR;

  3. 3.

Using multi-antenna transmission techniques is important in finite-SNR, even
when they are used only for spatial multiplexing and not to improve beamformer
directivity.

$0$$5$$10$$15$$20$$25$$30$$10$$20$$30$SNR [dB]Symmetric Rate
[nats/s]$L=4$$L=3$$L=2$$L=1$

Figure 2: Spatial multiplexing effect, $K=6$, $t=2$, $\alpha=L$ SNR ($\mathrm{dB}$) | 0 | 5 | 10 | 15 | 20 | 25 | 30
---|---|---|---|---|---|---|---
$L=4$ | 405 | 297 | 212 | 164 | 140 | 130 | 124
$L=3$ | 317 | 231 | 164 | 126 | 107 | 98 | 92
$L=2$ | 138 | 109 | 77 | 55 | 45 | 41 | 40
TABLE I: Rate benefit (%) over $L=1$ case, $K=6$, $t=2$, $\alpha=L$

###  III-B Beamformer directivity effect

As discussed, although choosing $\alpha<L$ decreases the DoF, it may improve
the finite-SNR performance by enhancing beamformer directivity. In fact, as we
decrease $\alpha$ (compared with $L$), the ratio of the variables to
constraints in the beamformer optimization problem grows larger, enabling
designing narrower beams that better direct data signals to end users [8],
thus improving the performance, especially in finite-SNR.

To investigate the effect of beamforming directivity, we have provided
simulation results for a CC-$(6,2,4,\alpha)$ setup, where
$\alpha\in\\{1,2,3,4\\}$, in Figure 3 and Table II. Again, optimized
beamformers are used to maximize the performance. Highlights of the results
are:

  1. 1.

Beamformer directivity has a very strong effect in finite-SNR. It can even
fully compensate for the performance loss due to decreased DoF. These results
comply with [8];

  2. 2.

The best choice seems to be choosing $\alpha$ to be slightly smaller than $L$.
This choice also slightly simplifies the beamformer design problem. Of course,
in smaller SNR values (the negative range, which is not considered here), this
recommendation may change.

$0$$5$$10$$15$$20$$25$$30$$10$$20$$30$SNR [dB]Symmetric Rate
[nats/s]$\alpha=4$$\alpha=3$$\alpha=2$$\alpha=1$

Figure 3: Beamformer directivity effect, $K=6$, $t=2$, $L=4$ SNR ($\mathrm{dB}$) | 0 | 5 | 10 | 15 | 20 | 25 | 30
---|---|---|---|---|---|---|---
$\alpha=4$ | -1 | 7 | 15 | 25 | 36 | 47 | 56
$\alpha=3$ | 7 | 13 | 20 | 29 | 37 | 43 | 47
$\alpha=2$ | 5 | 10 | 15 | 20 | 23 | 25 | 26
TABLE II: Rate benefit (%) over $\alpha=1$ case, $K=6$, $t=2$, $L=4$

###  III-C Beamformer structure effect

The beamforming vectors ${\mathbf{w}}_{{\mathcal{R}}({\mathcal{T}})}$ in (1)
can be designed in different ways. The simplest strategy is zero-forcing,
i.e., to design ${\mathbf{w}}_{{\mathcal{R}}({\mathcal{T}})}$ such that
${\mathbf{h}}_{k}^{T}{\mathbf{w}}_{{\mathcal{R}}({\mathcal{T}})}=0$ for every
user $k\in{\mathcal{R}}({\mathcal{T}})$. In other words,
${\mathbf{w}}_{{\mathcal{R}}({\mathcal{T}})}$ should lie in the null-space of
the matrix
${\mathbf{H}}=[{\mathbf{h}}_{k_{1}},...,{\mathbf{h}}_{k_{\alpha-1}}]$ formed
by concatenating channel vectors of users in set
${\mathcal{R}}({\mathcal{T}})=\\{k_{1},...,k_{\alpha-1}\\}$. This is
straightforward if $\alpha=L$, as in this case, the null-space is of dimension
one. However, if $\alpha<L$, the null-space has higher dimensions and the
beamformer vector can be any vector in that space. In this paper, we assume
the best vector (for maximizing the symmetric rate) is found in the null-space
by solving an optimization problem. The details are removed due to the lack of
space.

Of course, zero-forcing is not the optimum strategy in the finite-SNR regime
[8]. Instead, one needs to use optimized beamformers by maximizing the
symmetric rate given the SNR constraints. Indeed, this can result in non-
convex optimization problems, necessitating solutions such as successive
convex approximation (SCA) [8]. Of course, the underlying scheme can be
tweaked to reduce the optimized beamformer design complexity (e.g., using the
signal-level scheme in [15]). However, we still use the baseline MISO scheme
in [6] here, as it maximizes the positive effects of all CC-stemmed
parameters.

To investigate the effect of the beamformer structure, we have provided
numerical simulation results for zero-force and optimized beamforming
strategies for a CC-$(6,2,4,\alpha)$ setup, $\alpha\in\\{2,4\\}$, in Figure 4
and Table III. Highlights of the results are:

  1. 1.

Performance gap of the two strategies is very big in low-SNR but narrows down
as the SNR grows;

  2. 2.

Performance gap at low-SNR gets smaller if $\alpha<L$. This is because with
$\alpha<L$, zero-force beamformers are selected from a larger null-space [8];

  3. 3.

If computation capability is a bottleneck, choosing $\alpha<L$ and zero-force
beamforming is a wise choice for finite-SNR communications, as we also benefit
from improved beamformer directivity.

$0$$5$$10$$15$$20$$25$$30$$10$$20$$30$SNR [dB]Symmetric Rate
[nats/s]Optimized, $L=\alpha=4$Zeroforce, $L=\alpha=4$Optimized,
$L=4,\alpha=2$Zeroforce, $L=4,\alpha=2$

Figure 4: Beamformer structure effect, $K=6$, $t=2$, $L=4$ SNR ($\mathrm{dB}$) | 0 | 5 | 10 | 15 | 20 | 25 | 30
---|---|---|---|---|---|---|---
$\alpha=4$ | 582 | 347 | 167 | 69 | 27 | 9 | 2
$\alpha=2$ | 81 | 47 | 26 | 16 | 11 | 8 | 7
TABLE III: Rate benefit (%) over ZF beamforming, $K=6$, $t=2$, $L=4$

##  IV CC-stemmed parameters

In this section, we review the three CC-stemmed parameters in Figure 1 that
affect the finite-SNR performance of cache-aided MISO communications.

###  IV-A CC-aided multiplexing effect

CC-aided multiplexing is determined by the CC gain $t$ and indicates how much
interference can be removed by the cache contents of target users. As the
value of $t$ directly affects DoF, we expect the performance to improve as $t$
grows larger.

To investigate the effect of the CC-aided multiplexing gain, we use numerical
simulations for a CC-$(6,t,4,4)$ setup, where $t\in\\{0,1,2\\}$. We assume the
MISO-CC scheme of [6] is used as the baseline, and optimized beamformers are
applied to maximize the performance. The results are provided in Figure 5 and
Table IV. Highlights of the results are:

  1. 1.

As expected, better performance is attained with larger $t$;

  2. 2.

Performance improvement by cC-aided multiplexing is independent of the SNR
value. This is because cache-aided interference removal is also independent of
SNR;

  3. 3.

Cache-aided communications is a strong tool at finite-SNR as it does not
impose similar complexities of beamforming in that regime.

$0$$5$$10$$15$$20$$25$$30$$10$$20$$30$SNR [dB]Symmetric Rate
[nats/s]$t=2$$t=1$$t=0$

Figure 5: Cache-aided multiplexing effect, $K=6$, $\alpha=L=4$ SNR ($\mathrm{dB}$) | 0 | 5 | 10 | 15 | 20 | 25 | 30
---|---|---|---|---|---|---|---
$t=2$ | 75 | 90 | 101 | 105 | 105 | 101 | 94
$t=1$ | 48 | 53 | 57 | 61 | 62 | 60 | 56
TABLE IV: Rate benefit (%) over $t=0$ case, $K=6$, $\alpha=L=4$

###  IV-B CC-aided multicasting effect

So far, we have only considered bit-level MISO-CC schemes where the cache-
aided interference cancellation is performed after the received signal is
decoded. As mentioned earlier, another class of MISO-CC schemes with cache-
aided interference cancellation before decoding the received signal has
recently gained popularity due to nice properties such as reduced
subpacketization [10, 15], simpler beamformer design [15], and applicability
to dynamic setups [11]. However, it is also discussed in [13, 4] that these
nice properties are accompanied by penalties; most notably in reduced finite-
SNR performance and more complex signaling requirements in the control plane.
In this paper, we investigate the former issue and show that it is related to
a reduction in the coded multicasting gain. In [13], this gain is referred to
as efficiency index.

Let’s consider a CC-$(6,2,4,4)$ setup and review how the subpacketization can
be reduced from 15 to 3 using a signal-level CC scheme.
Placement phase. We split each file $W\in{\mathcal{F}}$ into $P=3$ subpackets
$W_{1}$, $W_{2}$, $W_{3}$, and store $W_{1}$ in users 1 and 2, $W_{2}$ in
users 3 and 4, and $W_{3}$ in users 4 and 5.
Delivery phase. Let us denote the files requested by users 1-6 as
$A,B,\cdots,F$, respectively. We only need the following transmission vector:

$\displaystyle{\mathbf{x}}_{1\cdots 6}=$ $\displaystyle
A_{2}{\mathbf{w}}_{256}+A_{3}{\mathbf{w}}_{234}+B_{2}{\mathbf{w}}_{156}+B_{3}{\mathbf{w}}_{134}$
(13)
$\displaystyle+C_{1}{\mathbf{w}}_{456}+C_{3}{\mathbf{w}}_{412}+D_{1}{\mathbf{w}}_{356}+D_{3}{\mathbf{w}}_{312}$
$\displaystyle+E_{1}{\mathbf{w}}_{634}+E_{2}{\mathbf{w}}_{612}+F_{1}{\mathbf{w}}_{534}+F_{2}{\mathbf{w}}_{512}.$

Let us review the decoding process at user 1, which is interested in the first
and second data terms in (13). Using (2), this user receives

$y_{1\cdots
6}(1)=A_{2}{\mathbf{h}}_{1}^{T}{\mathbf{w}}_{256}+A_{3}{\mathbf{h}}_{1}^{T}{\mathbf{w}}_{234}+{\mathbf{h}}_{1}^{T}I_{BF}+{\mathbf{h}}_{1}^{T}I_{CC}+z_{1},$

where the interference terms are

$\displaystyle I_{BF}=$ $\displaystyle
B_{2}{\mathbf{w}}_{156}+B_{3}{\mathbf{w}}_{134}+C_{3}{\mathbf{w}}_{412}$
$\displaystyle+D_{3}{\mathbf{w}}_{312}+E_{2}{\mathbf{w}}_{612}+F_{2}{\mathbf{w}}_{512},$
$\displaystyle I_{CC}=$ $\displaystyle
C_{1}{\mathbf{w}}_{456}+D_{1}{\mathbf{w}}_{356}+E_{1}{\mathbf{w}}_{634}+F_{1}{\mathbf{w}}_{534}.$

However, the interference terms in $I_{BF}$ are nulled-out or suppressed by
beamforming, and the ones in $I_{CC}$ could be regenerated and removed from
$y_{1\cdots 6}(1)$ using the cached contents of user 1. So, user 1 can decode
$A_{2}$ and $A_{3}$ using a SIC receiver after cache-aided interference
cancellation is done in the signal domain. Comparing ${\mathbf{x}}_{1\cdots
6}$ in (13) with its bit-level counterpart in (7), it is seen that in the
signal-level scheme, we have fewer codewords and hence, need less
multicasting.

As shown in [13], for our example CC-$(6,2,4,4)$ setup, different
subpacketization levels of $P\in\\{3,6,9,12,15\\}$ are possible. Moreover, as
we move from a pure signal-level scheme (smallest $P$) to a pure bit-level one
(largest $P$), we can use more multicasting in the transmission. Here, we have
used numerical simulations to compare the system performance for
$P\in\\{3,6,12,15\\}$. The results are provided in Figure 6 and Table V.
Highlights of the results are:

  1. 1.

Performance is improved as subpacketization is increased and more multicasting
is supported. This is because multicasting is more power-efficient than
unicasting;

  2. 2.

Performance improvement is more prominent in finite-SNR, as the rate is power-
limited in this regime.

$0$$5$$10$$15$$20$$25$$30$$10$$20$$30$SNR [dB]Symmetric Rate
[nats/s]$P=15$$P=12$$P=6$$P=3$

Figure 6: Cache-aided multicasting effect, $K=6$, $t=2$, $\alpha=L=4$ SNR ($\mathrm{dB}$) | 0 | 5 | 10 | 15 | 20 | 25 | 30
---|---|---|---|---|---|---|---
$P=15$ | 48 | 48 | 46 | 42 | 36 | 30 | 25
$P=12$ | 44 | 40 | 36 | 32 | 28 | 24 | 20
$P=6$ | 10 | 10 | 10 | 11 | 11 | 10 | 9
TABLE V: Rate benefit (%) over $P=3$ case, $K=6$, $t=2$ ,$\alpha=L=4$

$0$$5$$10$$15$$20$$25$$30$$10$$20$$30$SNR [dB]Symmetric Rate
[nats/s]$\beta=10$$\beta=5$$\beta=2$$\beta=1$

Figure 7: Superposition (MAC) effect, $K=6$, $t=2$, $\alpha=L=4$

###  IV-C Superposition (MAC) effect

In the transmission vector model in (1), each codeword $X_{{\mathcal{T}}}$,
${\mathcal{T}}\in\bar{{\mathcal{T}}}({\mathcal{S}})$, includes data for every
user $k\in{\mathcal{T}}$. Let us consider a specific user
$\bar{k}\in{\mathcal{S}}$. After the transmission of
${\mathbf{x}}_{{\mathcal{S}}}$, if there exist multiple sets
${\mathcal{T}}\in\bar{{\mathcal{T}}}({\mathcal{S}})$ that include $\bar{k}$,
this user has to decode its requested data from a multiple access channel
using a SIC receiver. However, SIC receivers are complex to implement, and
hence, more effort has been recently put into designing CC schemes without a
SIC requirement [12, 15]. Interestingly, it is also shown in [16] that
removing the SIC requirement could greatly simplify optimized beamformer
design through iterative optimization methods.

Removing the SIC requirement (or reducing the MAC size) is possible in both
bit- and signal-level schemes. However, signal-level schemes generally provide
more flexibility as they are less constrained by multicasting [12, 15]. Here,
we only consider controlling the MAC size in bit-level schemes, as they enable
CC-aided multicasting gain to be achieved at its full capacity. For example,
let us consider a CC-$(6,2,4,4)$ setup and its transmission vector
${\mathbf{x}}_{1\cdots 6}$ in 7. As discussed in Section III, every user needs
to decode its requested data from a MAC of size 10 after
${\mathbf{x}}_{1\cdots 6}$ is transmitted. However, for this network, we may
control the MAC size (and even avoid it altogether) by scheduling the
codewords sent by ${\mathbf{x}}_{1\cdots 6}$ into more transmissions. For
example, to avoid the MAC, we can use the following ten transmission vectors
instead:

$\displaystyle{\mathbf{x}}_{1\cdots 6}^{1}$
$\displaystyle=X_{123}{\mathbf{w}}_{456}+X_{456}{\mathbf{w}}_{123},$ (14)
$\displaystyle{\mathbf{x}}_{1\cdots 6}^{2}$
$\displaystyle=X_{124}{\mathbf{w}}_{356}+X_{356}{\mathbf{w}}_{124},$
$\displaystyle\cdots$ $\displaystyle{\mathbf{x}}_{1\cdots 6}^{10}$
$\displaystyle=X_{156}{\mathbf{w}}_{234}+X_{234}{\mathbf{w}}_{156}.$

Controlling the MAC size is first introduced in [8] using a $\beta$ parameter
for reducing beamformer design complexity. Accordingly, we also use $\beta$ to
denote the MAC size. Simulation results for the example CC-$(6,2,4,4)$ setup
are provided in Figure 7 and Table VI. It is assumed that
$\beta\in\\{1,2,5,10\\}$, and optimized beamformers are used to maximize the
performance. Highlights of the results are:

  1. 1.

Increasing the MAC size generally improves performance (we suspect the
deviations are due to numerical errors). This is due to improved superposition
coding [8]. However, compared with other parameters, the effect is small;

  2. 2.

Given the complexity of SIC receivers and the small gain of superposition
coding, expanding the literature on MAC-avoiding bit-level CC schemes is a
promising direction.

SNR ($\mathrm{dB}$) | 0 | 5 | 10 | 15 | 20 | 25 | 30
---|---|---|---|---|---|---|---
$\beta=10$ | 3 | 8 | 11 | 13 | 16 | 19 | 22
$\beta=5$ | 11 | 12 | 12 | 12 | 14 | 17 | 18
$\beta=2$ | 4 | 6 | 6 | 6 | 7 | 8 | 10
TABLE VI: Rate benefit (%) over $\beta=1$ case, $K=6$, $t=2$ ,$\alpha=L=4$

##  V Conclusion and Future Work

In this paper, we studied six parameters affecting the finite-SNR performance
of coded caching schemes in MISO setups. Out of the six parameters, three
stemmed from the multi-antenna transmission part and the rest relied on cache-
aided communications. For each parameter, we provided a brief explanation,
followed by simulation results and insights.

As a quick summary, both multi-antenna and cache-aided communication
techniques play important roles in finite-SNR. However, being DoF-optimal is
not always preferred; one can slightly reduce the spatial multiplexing gain
and still get better results due to improved beamformer directivity. Deviation
from the optimal point also narrows the performance gap of simple zero-force
beamformers with more complex optimized ones. On the other hand, CC gain is
independent of the SNR regime, bit-level schemes have better performance at
finite-SNR due to improved multicasting gain, and designing bit-level schemes
without SIC requirement is of interest.

Future research directions include expanding the results to multi-input multi-
output (MIMO) communication setups where the users are also equipped with
antenna arrays, a more thorough mathematical analysis of the effect of the
parameters, and expanding the results considering a more diverse set of
baseline CC schemes.

## References

  * [1] N. Rajatheva and et al., “White Paper on Broadband Connectivity in 6G,” _arXiv preprint arXiv:2004.14247_ , 2020. [Online]. Available: http://arxiv.org/abs/2004.14247
  * [2] M. A. Maddah-Ali and U. Niesen, “Fundamental limits of caching,” _IEEE Transactions on Information Theory_ , vol. 60, no. 5, pp. 2856–2867, 2014.
  * [3] H. B. Mahmoodi, M. J. Salehi, and A. Tolli, “Non-Symmetric Coded Caching for Location-Dependent Content Delivery,” _IEEE International Symposium on Information Theory - Proceedings_ , vol. 2021-July, pp. 712–717, 2021.
  * [4] M. Salehi, K. Hooli, J. Hulkkonen, and A. Tolli, “Enhancing Next-Generation Extended Reality Applications with Coded Caching,” _arXiv preprint arXiv:2202.06814_ , 2022. [Online]. Available: http://arxiv.org/abs/2202.06814
  * [5] S. P. Shariatpanahi, S. A. Motahari, and B. H. Khalaj, “Multi-server coded caching,” _IEEE Transactions on Information Theory_ , vol. 62, no. 12, pp. 7253–7271, 2016.
  * [6] S. P. Shariatpanahi, G. Caire, and B. Hossein Khalaj, “Physical-Layer Schemes for Wireless Coded Caching,” _IEEE Transactions on Information Theory_ , vol. 65, no. 5, pp. 2792–2807, 2019.
  * [7] M. J. Salehi, H. B. Mahmoodi, and A. Tölli, “A Low-Subpacketization High-Performance MIMO Coded Caching Scheme,” in _WSA 2021 - 25th International ITG Workshop on Smart Antennas_ , 2021, pp. 427–432.
  * [8] A. Tölli, S. P. Shariatpanahi, J. Kaleva, and B. H. Khalaj, “Multi-antenna interference management for coded caching,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 3, pp. 2091–2106, 2020\.
  * [9] Q. Yan, X. Tang, Q. Chen, and M. Cheng, “Placement Delivery Array Design Through Strong Edge Coloring of Bipartite Graphs,” _IEEE Communications Letters_ , vol. 22, no. 2, pp. 236–239, 2018.
  * [10] E. Lampiris and P. Elia, “Adding transmitters dramatically boosts coded-caching gains for finite file sizes,” _IEEE Journal on Selected Areas in Communications_ , vol. 36, no. 6, pp. 1176–1188, 2018.
  * [11] M. J. Salehi, E. Parrinello, H. B. Mahmoodi, and A. Tolli, “Low-Subpacketization Multi-Antenna Coded Caching for Dynamic Networks,” _2022 Joint European Conference on Networks and Communications and 6G Summit, EuCNC/6G Summit 2022_ , pp. 112–117, 2022.
  * [12] S. Mohajer and I. Bergel, “MISO Cache-Aided Communication with Reduced Subpacketization,” in _IEEE International Conference on Communications_ , vol. 2020-June. IEEE, 2020, pp. 1–6.
  * [13] M. Salehi, A. Tolli, S. P. Shariatpanahi, and J. Kaleva, “Subpacketization-rate trade-off in multi-antenna coded caching,” in _2019 IEEE Global Communications Conference, GLOBECOM 2019 - Proceedings_. IEEE, 2019, pp. 1–6.
  * [14] E. Lampiris and P. Elia, “Resolving a Feedback Bottleneck of Multi-Antenna Coded Caching,” _arXiv_ , 2018. [Online]. Available: http://arxiv.org/abs/1811.03935
  * [15] M. J. Salehi, E. Parrinello, S. P. Shariatpanahi, P. Elia, and A. Tolli, “Low-Complexity High-Performance Cyclic Caching for Large MISO Systems,” _IEEE Transactions on Wireless Communications_ , vol. 21, no. 5, pp. 3263–3278, 2022.
  * [16] H. B. Mahmoodi, B. Gouda, M. Salehi, and A. Tolli, “Low-complexity Multicast Beamforming for Multi-stream Multi-group Communications,” in _2021 IEEE Global Communications Conference, GLOBECOM 2021 - Proceedings_ , 2022, pp. 01–06.


Thus, we use the results from Lemma <ref> and Lemma <ref>, which yields that Equation (<ref>) is bounded by

\begin{align*}
\|\mathbf{G}_1(\mathcal{T}^l;\mathcal{S}) - \mathbf{G}_1(\mathcal{T}^l;\mathcal{S}')\|_F \leq& M_8 \|\mathbf{A} - \mathbf{A}'\|_F + (M_9+M_{10})\|\mathbf{X}-\mathbf{X}',\mathbf{B}-\mathbf{B}',\mathbf{C}-\mathbf{C}'\|_F \\
&+(M_9+M_{10}) \|\tilde{\mathbf{X}}-\tilde{\mathbf{X}}',\mathbf{B}-\mathbf{B}',\mathbf{C}-\mathbf{C}'\|_F
\end{align*}

There, with Lemma <ref>, we have with some constant,

\begin{align*}
C_9 = M_8 + (M_9+M_{10})\left(\sqrt{C_7^2+1}+\sqrt{C_8^2+1}\right),
\end{align*}

such that

\begin{align}
\|\mathbf{G}_1(\mathcal{T}^l;\mathcal{S}) - \mathbf{G}_1(\mathcal{T}^l;\mathcal{S}')\|_F \leq  C_9 \|\mathbf{A}-\mathbf{A}',\mathbf{B}-\mathbf{B}',\mathbf{C}-\mathbf{C}'\|_F. \label{eq:result_of_first_order_derivative}
\end{align}

Proof for the Newton Step.
The proof of the Newton step is built on the bound of gradient (we add and subtract terms and use triangular inequality),

\begin{align}
&\|\mathbf{G}_1(\mathcal{T}^l;\mathcal{S})\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S})-\mathbf{G}_1(\mathcal{T}^l;\mathcal{S}')\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S}')\|_F \notag\\
=&\|\mathbf{G}_1(\mathcal{T}^l;\mathcal{S})\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S})-\mathbf{G}_1(\mathcal{T}^l;\mathcal{S}')\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S}) + \mathbf{G}_1(\mathcal{T}^l;\mathcal{S}')\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S}) - \mathbf{G}_1(\mathcal{T}^l;\mathcal{S}')\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S}')\|_F \notag\\
\leq&\|\mathbf{G}_1(\mathcal{T}^l;\mathcal{S})-\mathbf{G}_1(\mathcal{T}^l;\mathcal{S}')\|_F\|\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S})\|_F + \|\mathbf{G}_1(\mathcal{T}^l;\mathcal{S}')\|_F\|\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S}) -\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S}')\|_F. \label{eq:result_of_second_order_derivative}
\end{align}

The first term is bounded by Equation (<ref>), and as for the second term, we use the result from Equation (<ref>),

\begin{align*}
&\|\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S}) -\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S}')\|_F \\
=&\|\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S})\left[\mathbf{H}_{1}(\mathcal{T}^l;\mathcal{S}') - \mathbf{H}_{1}(\mathcal{T}^l;\mathcal{S})\right]\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S}')\|_F \\
=&\left\|\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S})\left[\Pi(\mathbf{X}',\mathbf{B}',\mathbf{C}')+\Pi(\tilde{\mathbf{X}}',\mathbf{B}',\mathbf{C}') - \Pi(\mathbf{X},\mathbf{B},\mathbf{C})-\Pi(\tilde{\mathbf{X}},\mathbf{B},\mathbf{C})\right]\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S}')\right\|_F \\
\leq&\|\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S})\|_F\left(\left\|\Pi(\mathbf{X}',\mathbf{B}',\mathbf{C}')- \Pi(\mathbf{X},\mathbf{B},\mathbf{C})\right\|_F+\left\|\Pi(\tilde{\mathbf{X}}',\mathbf{B}',\mathbf{C}')-\Pi(\tilde{\mathbf{X}},\mathbf{B},\mathbf{C})\right\|_F\right)\|\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S}')\|_F.
\end{align*}

By using the results from Lemma <ref> and Lemma <ref>, we have

\begin{align}
&\|\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S}) -\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S}')\|_F \notag\\
\leq&C_4\max_{\mathcal{S}}\left\{\|\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S})\|_F^2\right\}\left(\|\mathbf{X}-\mathbf{X}',\mathbf{B}-\mathbf{B}',\mathbf{C}-\mathbf{C}'\|_F+\|\tilde{\mathbf{X}}-\tilde{\mathbf{X}}',\mathbf{B}-\mathbf{B}',\mathbf{C}-\mathbf{C}'\|_F\right) \notag\\
\leq&C_4\left(\sqrt{C_7^2+1}+\sqrt{C_8^2+1}\right)\max_{\mathcal{S}}\left\{\|\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S})\|_F^2\right\}\|\mathbf{A}-\mathbf{A}',\mathbf{B}-\mathbf{B}',\mathbf{C}-\mathbf{C}'\|_F. \label{eq:genius_step_also}
\end{align}

Therefore, by inserting the results from Equation (<ref>) and Equation (<ref>) into Equation (<ref>), the proof is completed for some constants

\begin{align*}
C_{10} = C_9\max_{\mathcal{S}}\{\|\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S})\|_F\} + C_4\left(\sqrt{C_7^2+1}+\sqrt{C_8^2+1}\right)\max_{\mathcal{S}'}\{\|\mathbf{G}_1(\mathcal{T}^l;\mathcal{S}')\|_F\}\max_{\mathcal{S}}\left\{\|\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S})\|_F^2\right\},
\end{align*}

such that

\begin{align*}
\|\mathbf{G}_1(\mathcal{T}^l;\mathcal{S})\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S})-\mathbf{G}_1(\mathcal{T}^l;\mathcal{S}')\mathbf{H}_{1}^{-1}(\mathcal{T}^l;\mathcal{S}')\|_F \leq C_{10} \|\mathbf{A}-\mathbf{A}',\mathbf{B}-\mathbf{B}',\mathbf{C}-\mathbf{C}'\|_F.
\end{align*}

For a up-to-date factor set $\mathcal{S}^{l,0}=\{\mathbf{A}^l,\mathbf{B}^l,\mathbf{C}^l\}$, for a new tensor and the augmented tensor batch $\mathcal{T}^j,\tilde{\mathcal{T}}^j$, we first conduct the Cold Start and Auxiliary Step to obtain $\{\mathbf{X}^j,\tilde{\mathbf{X}}^j\}$. Then, there exist constants, $0\leq C_{11}, C_{12}< \infty$, such that the gradient and the Newton step from the batch (e.g., a small scale of the objective in main text Equation (8)) are bounded,

\begin{align*}
\|\mathbf{G}_1(\mathcal{T}^j;\mathcal{S}^{l,0})-\mathbf{G}_1(\mathcal{T}^k;\mathcal{S}^{l,0})\|_F &\leq C_{11} \eta^l, \\
\|\mathbf{G}_1(\mathcal{T}^j;\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})-\mathbf{G}_1(\mathcal{T}^k;\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^k;\mathcal{S}^{l,0})\|_F &\leq C_{12} \eta^l.
\end{align*}

This bound states that given different batch data at the $l$-th batch, the gradient and the Newton step will be bounded by the learning rate. We only prove for $\mathbf{A}$. Similar derivations could be applied for $\mathbf{B}$ (given $\mathcal{S}^{l,1}$) and $\mathbf{C}$ (given $\mathcal{S}^{l,2}$).

Proof for the gradient. Following the moving average setting and Definition <ref>, the gradient is,

\begin{align*}
&\mathbf{G}_1(\mathcal{T}^j;\mathcal{S}^{l,0}) = \mathbf{A}^l\left[\Pi(\mathbf{X}^j,\mathbf{B}^l,\mathbf{C}^l)+\Pi(\tilde{\mathbf{X}}^j,\mathbf{B}^l,\mathbf{C}^l)+\alpha\mathbf{I}\right]- \left[\left(\eta^l\mathbf{T}^j_2 +(1-\eta^l)\mathbf{T}^l_2\right) \kr(\mathbf{X}^j,\mathbf{B}^l,\mathbf{C}^l)\right. \\
&\quad\quad\quad\quad\quad\quad\quad\quad\quad \left.+ \left(\eta^l\tilde{\mathbf{T}}^j_2 +(1-\eta^l)\tilde{\mathbf{T}}^l_2\right) \kr(\tilde{\mathbf{X}}^j,\mathbf{B}^l,\mathbf{C}^l)\right].
\end{align*}

where $\{\mathbf{X}^j,\tilde{\mathbf{X}}^j\}$ are the outputs of the Auxiliary Step from $\mathcal{T}^j,\tilde{\mathcal{T}}^j,\mathcal{S}^{l,0}$.

We can bound the difference of the first order derivatives by (using triangular inequality)

\begin{align}
&\|\mathbf{G}_1(\mathcal{T}^j;\mathcal{S}^{l,0}) - \mathbf{G}_1(\mathcal{T}^k;\mathcal{S}^{l,0})\|_F \notag\\
= & \left\|\mathbf{A}^l\left(\Pi(\mathbf{X}^j,\mathbf{B}^l,\mathbf{C}^l) -\Pi(\mathbf{X}^k,\mathbf{B}^l,\mathbf{C}^l) + \Pi(\tilde{\mathbf{X}}^j,\mathbf{B}^l,\mathbf{C}^l) - \Pi(\tilde{\mathbf{X}}^k,\mathbf{B}^l,\mathbf{C}^l)\right) \right.\notag\\
&+\eta^l\left(\mathbf{T}^k_2-\mathbf{T}^j_2\right)\kr(\mathbf{X}^j,\mathbf{B}^l,\mathbf{C}^l) + \left(\eta^l\mathbf{T}^k_2 +(1-\eta^l)\mathbf{T}^l_2\right)\kr(\mathbf{X}^k-\mathbf{X}^j,\mathbf{B}^l,\mathbf{C}^l) \notag \\
&\left.+\eta^l\left(\tilde{\mathbf{T}}_2^k-\tilde{\mathbf{T}}_2^j\right)\kr(\tilde{\mathbf{X}}^j,\mathbf{B}^l,\mathbf{C}^l) + \left(\eta^l\tilde{\mathbf{T}}^k_2 +(1-\eta^l)\tilde{\mathbf{T}}^l_2\right)\kr(\tilde{\mathbf{X}}^k-\tilde{\mathbf{X}}^j,\mathbf{B}^l,\mathbf{C}^l)\right\|_F \notag \\
\leq & \|\mathbf{A}^l\|_F\left(\left\|\Pi(\mathbf{X}^j,\mathbf{B}^l,\mathbf{C}^l) -\Pi(\mathbf{X}^k,\mathbf{B}^l,\mathbf{C}^l)\right\|_F + \left\| \Pi(\tilde{\mathbf{X}}^j,\mathbf{B}^l,\mathbf{C}^l) - \Pi(\tilde{\mathbf{X}}^k,\mathbf{B}^l,\mathbf{C}^l)\right\|_F\right)\notag\\
&+\eta^l\left\|\mathcal{T}^k-\mathcal{T}^j\right\|_F\left\|\kr(\mathbf{X}^j,\mathbf{B}^l,\mathbf{C}^l)\right\|_F + \left\|\eta^l\mathcal{T}^k +(1-\eta^l)\mathcal{T}^l\right\|_F\left\|\kr(\mathbf{X}^k-\mathbf{X}^j,\mathbf{B}^l,\mathbf{C}^l)\right\|_F \notag \\
&+\eta^l\left\|\tilde{\mathcal{T}}^k-\tilde{\mathcal{T}}^j\right\|_F\left\|\kr(\tilde{\mathbf{X}}^j,\mathbf{B}^l,\mathbf{C}^l)\right\|_F + \left\|\eta^l\tilde{\mathcal{T}}^k +(1-\eta^l)\tilde{\mathcal{T}}^l\right\|_F\left\|\kr(\tilde{\mathbf{X}}^k-\tilde{\mathbf{X}}^j,\mathbf{B}^l,\mathbf{C}^l)\right\|_F. \notag
\end{align}

For the first term, we use the results from Lemma <ref> and Lemma <ref>, which yields

\begin{align*}
\left\|\Pi(\mathbf{X}^j,\mathbf{B}^l,\mathbf{C}^l) -\Pi(\mathbf{X}^k,\mathbf{B}^l,\mathbf{C}^l)\right\|_F\leq C_4\|\mathbf{X}^j-\mathbf{X}^k\|_F\leq C_4C_5\eta^l, \\
\left\|\Pi(\tilde{\mathbf{X}}^j,\mathbf{B}^l,\mathbf{C}^l) -\Pi(\tilde{\mathbf{X}}^k,\mathbf{B}^l,\mathbf{C}^l)\right\|_F\leq C_4\|\tilde{\mathbf{X}}^j-\tilde{\mathbf{X}}^k\|_F\leq C_4C_6\eta^l.
\end{align*}

For the third and fifth term, we use the results from Lemma <ref> and Lemma <ref>, which yields

\begin{align*}
\left\|\kr(\mathbf{X}^k-\mathbf{X}^j,\mathbf{B}^l,\mathbf{C}^l)\right\|_F \leq C_3\|\mathbf{X}^k-\mathbf{X}^j\|_F \leq C_3C_5\eta^l, \\
\left\|\kr(\tilde{\mathbf{X}}^k-\tilde{\mathbf{X}}^j,\mathbf{B}^l,\mathbf{C}^l)\right\|_F \leq C_3\|\tilde{\mathbf{X}}^k-\tilde{\mathbf{X}}^j\|_F \leq C_3C_6\eta^l.
\end{align*}

We combine the results and get the following for some constant,

\begin{align*}
&C{11} = C_4(C_5+C_6)\max_{l}\{\|\mathbf{A}^l\|_F\}  \\
\end{align*}

such that

\begin{align}
\|\mathbf{G}_1(\mathcal{T}^j;\mathcal{S}^{l,0})-\mathbf{G}_1(\mathcal{T}^k;\mathcal{S}^{l,0})\|_F &\leq C_{11} \eta^l. \label{eq:bound_on_first_order}
\end{align}

Proof for the Newton Step.
The proof of the Newton step is built on the bound of gradient (we add and subtract terms and use triangular inequality),

\begin{align*}
&\|\mathbf{G}_1(\mathcal{T}^j;\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})-\mathbf{G}_1(\mathcal{T}^k;\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^k;\mathcal{S}^{l,0})\|_F \\
=&\|\mathbf{G}_1(\mathcal{T}^j;\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})-\mathbf{G}_1(\mathcal{T}^k;\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0}) \\
&\quad\quad\quad\quad\quad+ \mathbf{G}_1(\mathcal{T}^k;\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0}) - \mathbf{G}_1(\mathcal{T}^k;\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^k;\mathcal{S}^{l,0})\|_F \\
\leq&\|\mathbf{G}_1(\mathcal{T}^j;\mathcal{S}^{l,0})-\mathbf{G}_1(\mathcal{T}^k;\mathcal{S}^{l,0})\|_F\|\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})\|_F \\
&\quad\quad\quad\quad\quad + \|\mathbf{G}_1(\mathcal{T}^k;\mathcal{S}^{l,0})\|_F\|\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0}) -\mathbf{H}_1^{-1}(\mathcal{T}^k;\mathcal{S}^{l,0})\|_F.
\end{align*}

The first term is bounded by the gradient, and as for the second term, we have

\begin{align}
&\|\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0}) -\mathbf{H}_1^{-1}(\mathcal{T}^k;\mathcal{S}^{l,0})\|_F \notag\\
=&\|\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})\left[\mathbf{H}(\mathcal{T}^k;\mathcal{S}^{l,0}) - \mathbf{H}_1(\mathcal{T}^j;\mathcal{S}^{l,0})\right]\mathbf{H}_1^{-1}(\mathcal{T}^k;\mathcal{S}^{l,0})\|_F \notag\\
=&\left\|\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})\left[\Pi(\mathbf{X}^k,\mathbf{B}^l,\mathbf{C}^l)+\Pi(\tilde{\mathbf{X}}^k,\mathbf{B}^l,\mathbf{C}^l)\right.\right. \notag\\
&\quad\quad\quad\quad\left.\left.- \Pi(\mathbf{X}^j,\mathbf{B}^l,\mathbf{C}^l)-\Pi(\tilde{\mathbf{X}}^j,\mathbf{B}^l,\mathbf{C}^l)\right]\mathbf{H}_1^{-1}(\mathcal{T}^k;\mathcal{S}^{l,0})\right\|_F \notag\\
\leq&\|\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})\|_F\left(\left\|\Pi(\mathbf{X}^k,\mathbf{B}^l,\mathbf{C}^l)- \Pi(\mathbf{X}^j,\mathbf{B}^l,\mathbf{C}^l)\right\|_F\right. \notag\\
&\quad\quad\quad\quad \left.+\left\|\Pi(\tilde{\mathbf{X}}^k,\mathbf{B}^l,\mathbf{C}^l)-\Pi(\tilde{\mathbf{X}}^j,\mathbf{B}^l,\mathbf{C}^l)\right\|_F\right)\|\mathbf{H}_1^{-1}(\mathcal{T}^k;\mathcal{S}^{l,0})\|_F. \notag
\end{align}

By using the results from Lemma <ref> and Lemma <ref>, we have

\begin{align}
&\|\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0}) -\mathbf{H}_1^{-1}(\mathcal{T}^k;\mathcal{S}^{l,0})\|_F \notag\\
\leq&C_4\max_{j,l}\left\{\|\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})\|^2_F\right\}\left(\|\mathbf{X}^j-\mathbf{X}^k\|_F+\|\tilde{\mathbf{X}}^j-\tilde{\mathbf{X}}^k\|_F\right) \notag\\
\leq&C_4(C_5+C_6)\max_{j,l}\left\{\|\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})\|^2_F\right\}\eta^l. \label{eq:bound_on_H_inverse}
\end{align}

Therefore, by inserting the results of Equation (<ref>) and Equation (<ref>), we can find the constant,

\begin{align*}
C_{12} = C_{11}\max_{j,l}\left\{\|\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})\|_F\right\} + C_4(C_5+C_6)\max_{k,l}\{\|\mathbf{G}_1(\mathcal{T}^k;\mathcal{S}^{l,0})\|_F\}\max_{j,l}\left\{\|\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})\|^2_F\right\},
\end{align*}

such that

\begin{align*}
\|\mathbf{G}_1(\mathcal{T}^j;\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})-\mathbf{G}_1(\mathcal{T}^k;\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^k;\mathcal{S}^{l,0})\|_F \leq C_{12}\eta^l.
\end{align*}

(Adopted from Lemma 3.12 in <cit.>) Assume $\mathcal{S}^{l,0}=\{\mathbf{A}^{l},\mathbf{B}^l,\mathbf{C}^l\}$ is up-to-date base set, given $\mathcal{H}^{l}=\{\mathcal{T}^1,\dots,\mathcal{T}^{l-1}\}\bigcup \{\tilde{\mathcal{T}}^1,\dots,\tilde{\mathcal{T}}^{l-1}\}$ is the notation to represent historical information up to the $(l-1)$-th batch. $\{\mathbf{X}^j,\tilde{\mathbf{X}}^j\}$ and $\{\mathbf{X},\tilde{\mathbf{X}}\}$ are the outputs of the Auxiliary Step given a new batch $\mathcal{T}^j,\tilde{\mathcal{T}}^j$ or the whole tensor and the augmented tensor $\mathcal{T},\tilde{\mathcal{T}}$, separately, and the batch size is $b$, the size of the dataset is $N$. Then, there exists a constant $0\leq C_{13}<\infty$, given $\mathcal{H}^l$,

\begin{equation} \label{lemma9}
\mathbb{E}\left[\left\|\left(\mathbf{G}_1(\mathcal{T}^j;\mathcal{S}^{l,0}) - \frac{b}{N}\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\right)\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})\right\|_F \right] \leq C_{13}\eta^l.
\end{equation}

We show the proof for $\mathbf{A}$, while similar derivations can be applied to $\mathbf{B}$ (given $\mathcal{S}^{l,1}$) and $\mathbf{C}$ (given $\mathcal{S}^{l,2}$). We have

\begin{align}
&\mathbb{E}\left[\left\|\left(\mathbf{G}_1(\mathcal{T}^j;\mathcal{S}^{l,0}) - \frac{b}{N}\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\right)\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})\right\|^2_F \right] \label{eq:first_grad_replace_before}\\
=& \mathbb{E}_j\left[\left\|\left(\mathbf{G}_1(\mathcal{T}^j;\mathcal{S}^{l,0}) - \mathbb{E}_k\left[\mathbf{G}_1(\mathcal{T}^k;\mathcal{S}^{l,0})\right]\right)\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})\right\|^2_F \right] \label{eq:first_grad_replace_after}\\
\leq& \mathbb{E}_j\left[\left\|\mathbf{G}_1(\mathcal{T}^j;\mathcal{S}^{l,0}) - \mathbb{E}_k\left[\mathbf{G}_1(\mathcal{T}^k;\mathcal{S}^{l,0})\right]\right\|^2_F\left\|\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})\right\|^2_F \right] \\
\leq& \mathbb{E}_j\left[\left\|\mathbf{G}_1(\mathcal{T}^j;\mathcal{S}^{l,0}) - \mathbb{E}_k\left[\mathbf{G}_1(\mathcal{T}^k;\mathcal{S}^{l,0})\right]\right\|^2_F\right]\max_{j}\left\{\left\|\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})\right\|^2_F\right\} \label{eq:Jesen_before}\\
\leq& \mathbb{E}_j\left[\mathbb{E}_k\left[\left\|\mathbf{G}_1(\mathcal{T}^j;\mathcal{S}^{l,0}) - \mathbf{G}_1(\mathcal{T}^k;\mathcal{S}^{l,0})\right\|^2_F\right]\right]\max_{j}\left\{\left\|\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})\right\|^2_F\right\} \label{eq:use_of_gradient_bound_before}\\
\leq& \mathbb{E}_j\left[\mathbb{E}_k\left[(C_{11}\eta^l)^2\right]\right]\max_{j}\left\{\left\|\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})\right\|^2_F\right\} = C^2_{11}\max_{j}\left\{\left\|\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})\right\|^2_F\right\}(\eta^l)^2, \label{eq:use_of_gradient_bound_after}
\end{align}

where from Equation (<ref>) to Equation (<ref>), we use Lemma <ref>, from Equation (<ref>) to Equation (<ref>), we replace $\mathbf{G}_1(\mathcal{T}^j;\mathcal{S}^{l,0})$ with the expected gradient over the batch, from Equation (<ref>) to Equation (<ref>), we use Jensen's inequality (given that the square function being a convex function).

To this end, we achieve Equation (<ref>) using the above bound,

\begin{align*}
&\mathbb{E}\left[\left\|\left(\mathbf{G}_1(\mathcal{T}^j;\mathcal{S}^{l,0}) - \frac{b}{N}\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\right)\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})\right\|_F \right] \\
=&\mathbb{E}\left[\sqrt{\left\|\left(\mathbf{G}_1(\mathcal{T}^j;\mathcal{S}^{l,0}) - \frac{b}{N}\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\right)\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})\right\|^2_F} \right] \\
\leq&\sqrt{\mathbb{E}\left[\left\|\left(\mathbf{G}_1(\mathcal{T}^j;\mathcal{S}^{l,0}) - \frac{b}{N}\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\right)\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})\right\|^2_F \right]} \\
\leq& C_{11}\sqrt{\max_{j}\left\{\left\|\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})\right\|^2_F\right\}}\eta^l.
\end{align*}

The above inequality holds due to the Jensen's inequality (square root being a concave function).
The lemma is completed by defining

\begin{equation*}
C_{13} = C_{11}\sqrt{\max_{j}\left\{\left\|\mathbf{H}_1^{-1}(\mathcal{T}^j;\mathcal{S}^{l,0})\right\|^2_F\right\}}.
\end{equation*}

[Lemma A.5 in <cit.>]
Let $(a^l)_{l\geq 1}$ and $(b^l)_{l\geq 1}$ be two non-negative real sequences such that the series $\sum_{l=1}^\infty a^l$ diverges, the series $\sum_{l=1}^\infty a^l b^l$ converges, and there exists $K>0$ suth that $|b^{l+1}-b^l|\leq Ka^l$. Then, the sequence $(b^l)_{l\geq 1}$ converges to 0.

§.§ Proof of the Main Theorem

We now prove the convergence of $\mathbf{A}$ in the main theorem, and similar derivations can be applied to $\mathbf{B},\mathbf{C}$. We let $a^l =\eta^l$ and $b^l=\mathbb{E}\left[\|\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\|_F^2\right]$, which is the short-hand quantity of gradient defined in Definition <ref>.

Proof of the first condition. First, we know that the harmonic sequence $\{\frac1l\}_{l\geq 1}$ diverges, which satisfies the first condition,

\begin{equation}
\sum_{l=1}^\infty a^l \geq C_{min} \sum_{l=1}^\infty \frac1l\rightarrow\infty.
\end{equation}

Proof of the second condition. Second, we prove that $\sum_{l=1}^\infty a^l b^l$ converges. According to the main paper, the update rule of $\mathbf{A}$ is as follows,

\begin{equation*}
\mathbf{A}^{l+1} = \eta^l\mathbf{A}^l+(1-\eta^l)\mathbf{A}^*,
\end{equation*}

where $\mathbf{A}^*$ is the closed-form solution of the least square problem. In fact, our update is one Newton step with learning rate $\eta^l$, where the original form is

\begin{equation} \label{eq:newton_step}
\mathbf{A}^{l+1} = \mathbf{A}^l-\eta^l\mathbf{G}_1(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{H}^{-1}_1(\mathcal{T}^l;\mathcal{S}^{l,0}),
\end{equation}

and $\mathbf{A}^*$ is,

\begin{equation*}
\mathbf{A}^* = \mathbf{A}^l - \mathbf{G}_1(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{H}^{-1}_1(\mathcal{T}^l;\mathcal{S}^{l,0}).
\end{equation*}

Taylor Expansion. We have that for any quadratic function $f$,

\begin{align*}
&f(\mathbf{x_1},\ldots,\mathbf{x_k} +\mathbf{\Delta x_k}, \ldots, \mathbf{x_n})  = f(\mathbf{x_1},\ldots,\mathbf{x_k}, \ldots, \mathbf{x_n}) + \mathbf{\Delta x_k}^T\mathbf{\nabla_{\mathbf{x_k}} }f+ \frac{1}{2}\mathbf{\Delta x_k}^T\mathbf{H}_f(\mathbf{x_k})\mathbf{\Delta x_k},
\end{align*}

where $\nabla_{\mathbf{x_k}}f$ is the gradient vector and $\mathbf{H}_f(\mathbf{x_k})$ is the Hessian matrix.

Using the above Taylor expansion and using the  matrix form we get

\begin{align}
&\mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,1}) - \mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,0}) \notag\\
% =& \left(\left\|\mathcal{T}-\llbracket\mathbf{X}^{l+1},\mathbf{A}^{l+1},\mathbf{B}^l,\mathbf{C}^l\rrbracket\right\|_F^2 + \left\|\tilde{\mathcal{T}}-\llbracket\tilde{\mathbf{X}}^{l+1},\mathbf{A}^{l+1},\mathbf{B}^l,\mathbf{C}^l\rrbracket\right\|_F^2+\alpha \|\mathbf{A}^{l+1}\|_F\right) \\
% &- \left(\left\|\mathcal{T}-\llbracket\mathbf{X}^{l+1},\mathbf{A}^l,\mathbf{B}^l,\mathbf{C}^l\rrbracket\right\|_F^2 + \left\|\tilde{\mathcal{T}}-\llbracket\tilde{\mathbf{X}}^{l+1},\mathbf{A}^l,\mathbf{B}^l,\mathbf{C}^l\rrbracket\right\|_F^2+\alpha \|\mathbf{A}^{l}\|_F\right)\notag \\
% =& \tr\left(\left(\mathbf{T}_2-\mathbf{A}^{l+1}\kr^\top(\mathbf{X}^{l+1},\mathbf{B}^l,\mathbf{C}^l)\right)^\top\left(\mathbf{T}_2-\mathbf{A}^{l+1}\kr^\top(\mathbf{X}^{l+1},\mathbf{B}^l,\mathbf{C}^l)\right)\right) \\
% &-\tr\left(\left(\mathbf{T}_2-\mathbf{A}^{l}\kr^\top(\mathbf{X}^{l+1},\mathbf{B}^l,\mathbf{C}^l)\right)^\top\left(\mathbf{T}_2-\mathbf{A}^{l}\kr^\top(\mathbf{X}^{l+1},\mathbf{B}^l,\mathbf{C}^l)\right)\right) \\
% &+ \tr\left(\left(\tilde{\mathbf{T}}_2-\mathbf{A}^{l+1}\kr^\top(\tilde{\mathbf{X}}^{l+1},\mathbf{B}^l,\mathbf{C}^l)\right)^\top\left(\tilde{\mathbf{T}}_2-\mathbf{A}^{l+1}\kr^\top(\tilde{\mathbf{X}}^{l+1},\mathbf{B}^l,\mathbf{C}^l)\right)\right) \\
% &-\tr\left(\left(\tilde{\mathbf{T}}_2-\mathbf{A}^{l}\kr^\top(\tilde{\mathbf{X}}^{l+1},\mathbf{B}^l,\mathbf{C}^l)\right)^\top\left(\tilde{\mathbf{T}}_2-\mathbf{A}^{l}\kr^\top(\tilde{\mathbf{X}}^{l+1},\mathbf{B}^l,\mathbf{C}^l)\right)\right) \\
% &+ \alpha \tr\left((\mathbf{A}^{l+1}-\mathbf{A}^{l})^\top(\mathbf{A}^{l+1}+\mathbf{A}^{l})\right).
% =& \tr\left((\eta^l\mathbf{G}_1(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{H}^{-1}_1(\mathcal{T}^l;\mathcal{S}^{l,0})\kr^\top(\mathbf{X}^{l+1},\mathbf{B}^l,\mathbf{C}^l))^\top\right) \\
=&-\eta^l \tr\left(\mathbf{G}_1(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{G}_1^\top(\mathcal{T};\mathcal{S}^{l,0})\right) \notag\\
&+ \frac{(\eta^l)^2}{2} \tr\left(\mathbf{G}_1(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{H}_1(\mathcal{T};\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{G}^\top_1(\mathcal{T}^l;\mathcal{S}^{l,0})\right). \label{eq:taylor_expansion}
\end{align}

Following the Definition <ref>, we sum these three results and get Equation (<ref>).

By rearranging the terms in Equation (<ref>), we have

\begin{align}
&\eta^l \tr\left(\frac{b}{N}\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{G}_1^\top(\mathcal{T};\mathcal{S}^{l,0})\right) \notag\\
=& -\eta^l Tr\left(\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^l;\mathcal{S}^{l,0})\left(\mathbf{G}^\top_1(\mathcal{T}^l;\mathcal{S}^{l,0})-\frac{b}{N}\mathbf{G}^\top_1(\mathcal{T};\mathcal{S}^{l,0})\right)\right) \notag\\
&+\mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,0}) -\mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,1}) \notag\\
&+\frac{(\eta^l)^2}{2} \tr\left(\mathbf{G}_1(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{H}_1(\mathcal{T};\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{G}^\top_1(\mathcal{T}^l;\mathcal{S}^{l,0})\right). \label{eq:taylor}
\end{align}

F-Norm Inequality. For now, let us consider the following problem: given a symmetric positive-definite matrix, $\Sigma\in\mathbb{R}^{d\times d}$, its inverse is also symmetric positive-definite. Suppose $\Sigma^{-1}\succ a>0$, given any non-zero matrix $\mathbf{Y}\in\mathbb{R}^{s\times d}$, we have the following trace inequality,

\begin{align}
\tr\left(\mathbf{Y}(\Sigma^{-1}-a)\mathbf{Y}^\top\right) >0~~~\Rightarrow~~~\tr\left(\mathbf{Y}\Sigma^{-1}\mathbf{Y}^\top\right) >a\|\mathbf{Y}\|_F^2. \label{eq:f_norm_inequality}
\end{align}

We use the result of this F-norm Inequality. By Definition <ref>, $\mathbf{H}_1(\mathcal{T}^l;\mathcal{S}^{l,0})$ is a positive-definite matrix. We let

\begin{align*}
\Sigma &= \mathbf{H}_1(\mathcal{T}^l;\mathcal{S}^{l,0}) \\
\end{align*}

Then, we have, with some constant $C_{14}>0$,

\begin{align*}
\eta^l \tr\left(\frac{b}{N}\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{G}_1^\top(\mathcal{T};\mathcal{S}^{l,0})\right) > \frac{1}{C_{14}}\eta^l \left\|\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\right\|_F^2.
\end{align*}

The LHS can be replaced by using the results in Equation (<ref>), and the expectation of the RHS is $a^lb^l$ scaled by a constant. We thus take the expectation on both side for the above inequality and let

\begin{align*}
Q^l_1 &= \mathbb{E}\left[-\eta^l Tr\left(\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^l;\mathcal{S}^{l,0})\left(\mathbf{G}^\top_1(\mathcal{T}^l;\mathcal{S}^{l,0})-\frac{b}{N}\mathbf{G}^\top_1(\mathcal{T};\mathcal{S}^{l,0})\right)\right)\right], \\
Q^l_2 &= \mathbb{E}\left[\mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,0}) - \mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,1})\right], \\
Q^l_3 &=  \mathbb{E}\left[\frac{(\eta^l)^2}{2} \tr\left(\mathbf{G}_1(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{H}_1(\mathcal{T};\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{G}^\top_1(\mathcal{T}^l;\mathcal{S}^{l,0})\right)\right].
\end{align*}

Since $a^lb^l=\mathbb{E}\left[\eta^l \left\|\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\right\|_F^2\right]$, if we can bound $\sum_{l=1}^\infty Q_1^l$, $\sum_{l=1}^\infty Q_2^l$, $\sum_{l=1}^\infty Q_3^l$, then, the second condition will be satisfied: $\sum_{l=1}^\infty a^l b^l$ converges. Now, let us show how to bound three summation terms.

For $\sum_{l=1}^\infty Q_1^l$, we have

\begin{align}
\sum_{l=1}^\infty Q_1^l =&  \sum_{l=1}^\infty \mathbb{E}\left[-\eta^l Tr\left(\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^l;\mathcal{S}^{l,0})\left(\mathbf{G}^\top_1(\mathcal{T}^l;\mathcal{S}^{l,0})-\frac{b}{N}\mathbf{G}^\top_1(\mathcal{T};\mathcal{S}^{l,0})\right)\right)\right] \notag\\
\leq& \sum_{l=1}^\infty \eta^l\mathbb{E}\left[ \left\|\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\right\|_F\left\|\mathbf{H}_1^{-1}(\mathcal{T}^l;\mathcal{S}^{l,0})\left(\mathbf{G}^\top_1(\mathcal{T}^l;\mathcal{S}^{l,0})-\frac{b}{N}\mathbf{G}^\top_1(\mathcal{T};\mathcal{S}^{l,0})\right)\right\|_F\right] \notag\\
\leq& \sum_{l=1}^\infty \eta^l\max_{l}\{ \left\|\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\right\|_F\}\mathbb{E}\left[ \left\|\mathbf{H}_1^{-1}(\mathcal{T}^l;\mathcal{S}^{l,0})\left(\mathbf{G}^\top_1(\mathcal{T}^l;\mathcal{S}^{l,0})-\frac{b}{N}\mathbf{G}^\top_1(\mathcal{T};\mathcal{S}^{l,0})\right)\right\|_F\right]\label{eq:Q1_before} \\
\leq& \sum_{l=1}^\infty C_{13}(\eta^l)^2\max_{l}\{ \left\|\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\right\|_F\} <\infty \label{eq:Q1_after},
\end{align}

where from Equation (<ref>) to Equation (<ref>), we use the results from Lemma <ref>.

For $\sum_{l=1}^\infty Q_2^l$, we have,

\begin{align}
\sum_{l=1}^\infty Q_2^l =& \sum_{l=1}^\infty\mathbb{E}\left[\mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,0}) - \mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,1})\right] \notag\\
=& \mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,0}) + \sum_{l=1}^\infty\mathbb{E}\left[\mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l+1,0}) - \mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,1})\right] \notag\\
=& \mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,0}) + \sum_{l=1}^\infty\mathbb{E}\left[\mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l+1,0}) - \mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,2})+\mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,2}) - \mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,1})\right] \label{eq:change_S_before}\\
\leq& \mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,0}) + \sum_{l=1}^\infty\mathbb{E}\left[\mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,3}) - \mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,2})+\mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,2}) - \mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,1})\right] \label{eq:change_S_after}\\
=& \mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,0}) + \sum_{l=1}^\infty\mathbb{E}\left[\mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,3}) - \mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,2})\right]+\sum_{l=1}^\infty\mathbb{E}\left[\mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,2}) - \mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,1})\right] \label{eq:loss_decrease}\\
\leq & \mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,0}) + 2M_9 \label{eq:loss_final_bound}.
\end{align}

From Equation (<ref>) to Equation (<ref>),

$\mathcal{S}^{l+1,0} = \mathcal{S}^{l,3}= \{\mathbf{A}^{l+1},\mathbf{B}^{l+1},\mathbf{C}^{l+1}\}$, however, the loss $\mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l+1,0})$ is calculated based on the updated coefficient matrices $\{\mathbf{X}^{l+2},\tilde{\mathbf{X}}^{l+2}\}$, which is optimized w.r.t. $\{\mathbf{A}^{l+1},\mathbf{B}^{l+1},\mathbf{C}^{l+1}\}$, while $\mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,3})$ uses old coefficient set $\{\mathbf{X}^{l+1},\tilde{\mathbf{X}}^{l+1}\}$. Thus,

\begin{equation*}
\mathbb{E}\left[\mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l+1,0})\right] < \mathbb{E}\left[\mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,3})\right].
\end{equation*}

From Equation (<ref>) to Equation (<ref>), we show that

\begin{align}
&\mathbb{E}\left[\mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,3}) - \mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,2})\right] \notag\\
=&-\eta^l \mathbb{E}\left[\tr\left(\mathbf{G}_3(\mathcal{T}^l;\mathcal{S}^{l,2})\mathbf{H}_3^{-1}(\mathcal{T}^l;\mathcal{S}^{l,2})\mathbf{G}_3^\top(\mathcal{T};\mathcal{S}^{l,2})\right)\right] \label{eq:split_before}\\
&+ \frac{(\eta^l)^2}{2} \mathbb{E}\left[\tr\left(\mathbf{G}_3(\mathcal{T}^l;\mathcal{S}^{l,2})\mathbf{H}_3^{-1}(\mathcal{T}^l;\mathcal{S}^{l,2})\mathbf{H}_3(\mathcal{T};\mathcal{S}^{l,2})\mathbf{H}_3^{-1}(\mathcal{T}^l;\mathcal{S}^{l,2})\mathbf{G}^\top_3(\mathcal{T}^l;\mathcal{S}^{l,2})\right)\right]\\
=&-\eta^l \mathbb{E}\left[\tr\left(\left(\mathbf{G}_3(\mathcal{T}^l;\mathcal{S}^{l,2})-\frac{b}{N}\mathbf{G}_3(\mathcal{T};\mathcal{S}^{l,2})\right)\mathbf{H}_3^{-1}(\mathcal{T}^l;\mathcal{S}^{l,2})\mathbf{G}_3^\top(\mathcal{T};\mathcal{S}^{l,2})\right)\right] \label{eq:split_after1}\\
&-\eta^l \frac{b}{N}\mathbb{E}\left[\tr\left(\mathbf{G}_3(\mathcal{T};\mathcal{S}^{l,2})\mathbf{H}_3^{-1}(\mathcal{T}^l;\mathcal{S}^{l,2})\mathbf{G}_3^\top(\mathcal{T};\mathcal{S}^{l,2})\right)\right] \label{eq:split_after2}\\
&+ \frac{(\eta^l)^2}{2} \mathbb{E}\left[\tr\left(\mathbf{G}_3(\mathcal{T}^l;\mathcal{S}^{l,2})\mathbf{H}_3^{-1}(\mathcal{T}^l;\mathcal{S}^{l,2})\mathbf{H}_3(\mathcal{T};\mathcal{S}^{l,2})\mathbf{H}_3^{-1}(\mathcal{T}^l;\mathcal{S}^{l,2})\mathbf{G}^\top_3(\mathcal{T}^l;\mathcal{S}^{l,2})\right)\right] \notag\\
\leq&~ C_{13}(\eta^l)^2\max_l\left\{\|\mathbf{G}_3(\mathcal{T};\mathcal{S}^{l,2})\|_F\right\} \label{eq:after_lemma_8}\\
&+ \frac{(\eta^l)^2}{2} \max_{l}\left\{\tr\left(\mathbf{G}_3(\mathcal{T}^l;\mathcal{S}^{l,2})\mathbf{H}_3^{-1}(\mathcal{T}^l;\mathcal{S}^{l,2})\mathbf{H}_3(\mathcal{T};\mathcal{S}^{l,2})\mathbf{H}_3^{-1}(\mathcal{T}^l;\mathcal{S}^{l,2})\mathbf{G}^\top_3(\mathcal{T}^l;\mathcal{S}^{l,2})\right)\right\} \notag.
\end{align}

From Equation (<ref>) to Equation (<ref>)(<ref>), we add and subtract one term. From Equation (<ref>) to Equation (<ref>), we apply Lemma (<ref>). For Equation (<ref>), since the term is always negative (refer to Equation (<ref>)), we directly throw it away. Similar analysis can be applied for $\sum_{l=1}^\infty\mathbb{E}\left[\mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,2}) - \mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,1})\right]$, so we have with some constant, $M_9>0$,

\begin{align*}
&\sum_{l=1}^\infty\mathbb{E}\left[\mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,3}) - \mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,2})\right] \leq M_9, \\
&\sum_{l=1}^\infty\mathbb{E}\left[\mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,2}) - \mathcal{L}(\mathcal{T},\tilde{\mathcal{T}};\mathcal{S}^{l,1})\right] \leq M_9.
\end{align*}

For $\sum_{l=1}^\infty Q_3^l$, we have

\begin{align*}
\sum_{l=1}^\infty Q_3^l =&\sum_{l=1}^\infty\mathbb{E}\left[\frac{(\eta^l)^2}{2} \tr\left(\mathbf{G}_1(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{H}_1(\mathcal{T};\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{G}^\top_1(\mathcal{T}^l;\mathcal{S}^{l,0})\right)\right] \\
\leq &\sum_{l=1}^\infty\frac{(\eta^l)^2}{2}\max\left\{ \left\|\mathbf{G}_1(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{H}_1(\mathcal{T};\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T}^l;\mathcal{S}^{l,0})\mathbf{G}^\top_1(\mathcal{T}^l;\mathcal{S}^{l,0})\right\|_F\right\} < \infty.
\end{align*}

Proof of the third condition. Third, we prove that $|b^{l+1}-b^l|$ is bounded by $a^l$,

\begin{align}
&\left|\mathbb{E}\left[\|\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l+1,0})\|_F^2\right]-\mathbb{E}\left[\|\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\|_F^2\right]\right| \notag\\
=& \left|\mathbb{E}\left[\|\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l+1,0})\|_F^2-\|\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\|_F^2\right]\right| \notag\\
=& \mathbb{E}\left[\tr\left(\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l+1,0})^\top \mathbf{G}_1(\mathcal{T};\mathcal{S}^{l+1,0})-\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})^\top \mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\right)\right] \notag\\
=& \mathbb{E}\left[\tr\left(\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l+1,0})+\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\right)^\top\left( \mathbf{G}_1(\mathcal{T};\mathcal{S}^{l+1,0})-\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\right)\right] \notag\\
\leq & 2\max_l\{\|\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\|_F\} \mathbb{E}\left[\|\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l+1,0})-\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\|_F\right] \label{eq:third_condition_before}\\
\leq & 2\max_l\{\|\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\|_F\}C_7 \mathbb{E}\left[\|\mathbf{A}^{l+1}-\mathbf{A}^{l},\mathbf{B}^{l+1}-\mathbf{B}^{l},\mathbf{C}^{l+1}-\mathbf{C}^{l}\|_F\right] \label{eq:third_condition_after},
\end{align}

where from Equation (<ref>) to Equation (<ref>), we use Lemma <ref>.

Finally, by the (Newton step based) update rule of the factors in Equation (<ref>), we have that

\begin{align*}
&\mathbf{A}^{l+1}-\mathbf{A}^{l} = -\eta^l\mathbf{G}_1\mathcal{L}(\mathcal{T};\mathcal{S}^{l,0})\mathbf{H}_1^{-1}(\mathcal{T};\mathcal{S}^{l,0}),\\
&\mathbf{B}^{l+1}-\mathbf{B}^{l} = -\eta^l\mathbf{G}_2\mathcal{L}(\mathcal{T};\mathcal{S}^{l,1})\mathbf{H}_2^{-1}(\mathcal{T};\mathcal{S}^{l,1}),\\
&\mathbf{C}^{l+1}-\mathbf{C}^{l} = -\eta^l\mathbf{G}_3\mathcal{L}(\mathcal{T};\mathcal{S}^{l,2})\mathbf{H}_3^{-1}(\mathcal{T};\mathcal{S}^{l,2}).
\end{align*}

We insert the results into Equation (<ref>) and define

\begin{align*}
K =2\sqrt{3}C_7\max_l\left\{\|\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\|_F\right\}\max_{i\in[0,1,2], l}\left\{\|\mathbf{G}_1\mathcal{L}(\mathcal{T};\mathcal{S}^{l,i})\mathbf{H}_1^{-1}(\mathcal{T};\mathcal{S}^{l,i})\|_F\right\},
\end{align*}

so that,

\begin{equation}
|b^{l+1}-b^l|=\left|\mathbb{E}\left[\|\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l+1,0})\|_F^2\right]-\mathbb{E}\left[\|\mathbf{G}_1(\mathcal{T};\mathcal{S}^{l,0})\|_F^2\right]\right| \leq K\eta^l=Ka^l.
\end{equation}

Consequently, we have proven that three conditions of Lemma <ref> are satisfied, which completes the proof of the main convergence theorem.


# Infinite derivative gravity resolves nonscalar curvature singularities

Ivan Kolář<EMAIL_ADDRESS>Institute of Theoretical Physics, Faculty of
Mathematics and Physics, Charles University, V Holešovičkách 2, Prague 180 00,
Czech Republic  Van Swinderen Institute, University of Groningen, 9747 AG,
Groningen, Netherlands  Tomáš Málek<EMAIL_ADDRESS>Institute of
Mathematics of the Czech Academy of Sciences, Žitná 25, 115 67 Prague 1, Czech
Republic

###### Abstract

We explicitly demonstrate that the nonlocal ghost-free ultraviolet
modification of general relativity (GR) known as the infinite derivative
gravity (IDG) resolves nonscalar curvature singularities in exact solutions of
the full theory. We analyze exact pp-wave solutions of GR and IDG describing
gravitational waves generated by null radiation. Curvature of GR and IDG
solutions with the same energy-momentum tensor is compared in parallel-
propagated frames along timelike and null geodesics at finite values of the
affine parameter. While the GR pp-wave solution contains a physically
problematic nonscalar curvature singularity at the location of the source, the
curvature of its IDG counterpart is finite.

##  I Introduction

Terms with higher derivatives appear in many effective descriptions of quantum
gravity. They mitigate issues related to ultraviolet (UV) incompleteness of
general relativity (GR) such as the existence of curvature singularities and
quantum non-renormalizability [1, 2]. Although higher-derivative theories of
finite order already regularize the linearized solutions (and their linearized
curvature) to a certain order [3], they suffer from the Ostrogradsky
instability and ghosts. Infinite derivative gravity (IDG) [4, 5, 6, 7, 8, 9]
is a promising UV modification of GR that resolves these issues by introducing
nonlocal analytic operators $\mathcal{F}(\Box)$ in the action. At the
linearized level, this theory is ghost-free, singularity-free, renormalizable,
and asymptotes to linearized GR in the infrared (IR). Naturally, the major
open question is whether the nonlocality of IDG can resolve curvature
singularities at the level of exact solutions of the full theory.

Due to the immense complexity of the field equations [10], the known exact
solutions of IDG are very scarce. Putting aside the universal spacetimes [11]
that solve practically every gravitational theory in vacuum (up to an
appropriate cosmological constant), the known exact solutions are either i)
FLRW spacetimes satisfying recursive ansatz such as $\Box R=\alpha R+\beta$
[7, 12, 13] or ii) pp-wave and Kundt geometries with trace-free (TF) Ricci and
Weyl of types III/N [14, 15, 16, 17]. Since the above recursive ansatz
effectively ‘localizes’ the field equations, the former are actually quite
insensitive to nonlocality.111They are affected by $\mathcal{F}(\Box)$ only
through a finite number of values, $\mathcal{F}(0)$, $\mathcal{F}(\alpha)$,
and $\mathcal{F}^{\prime}(\alpha)$, meaning that one can always find an
equivalent local theory (with polynomial $\mathcal{F}$) that admits the same
spacetimes with the same matter content. On the other hand, the latter still
retain strong dependence on the nonlocal operators $\mathcal{F}(\Box)$, so
they provide better testbeds of the singularity resolution due to nonlocality.

In this Letter, we study exact pp-wave solutions of GR and IDG generated by
the same energy-momentum tensor — null radiation with the point-like
distribution in the transverse 2-space and an arbitrary wave profile. In the
context of IDG, these Aichelburg–Sexl-type solutions were first studied at the
linearized level in [18] and later promoted to the exact solutions of the full
theory in [14]. Note that we actually consider an improved version of IDG. [In
particular, our theory (1) contains the Weyl term instead of the Ricci term.
We also demand suppression of the propagator at both ends of the spectrum of
$\bar{\Box}$ in (4), which should improve the behavior of the theory not only
for short spacelike but also timelike distances.] Nevertheless, our
conclusions remain valid also for the solutions from [18, 14].

The goal of this Letter is to analyze the curvature of these pp-wave solutions
and explicitly show the presence of curvature singularities in the GR solution
and their absence in the IDG solution. Recall that curvature singularities can
be either scalar or nonscalar. Both of them are equally unwanted and
physically problematic since they are associated with unboundedly large tidal
forces experienced by an observer in a finite time. The pp-waves of type N are
free from scalar curvature singularites (since all scalar invariant vanish)
but they still may contain a nonscalar curvature singularity. This is defined
as a point at the boundary of the manifold with diverging components of the
Riemann tensor when measured in a parallel-propagated (PP) frame along a
timelike or null geodesic that reaches the point with a finite value of affine
parameter. In what follows, we explicitly demonstrate that the curvature in
such PP frames blow up for the GR solution but remain finite for the IDG
solution.

##  II Infinite derivative gravity

Consider a general class of gravitational theories described by the Lagrangian

$\displaystyle L$
$\displaystyle=\tfrac{1}{2}\big{[}\varkappa^{-1}R-\tfrac{1}{3}R\mathcal{F}_{0}(\square)R+C^{abcd}\mathcal{F}_{2}(\square)C_{abcd}\big{]}+L_{\textrm{m}}\;,$
(1)

where $\mathcal{F}_{i}$, $i=0,2$, are analytic functions of the wave operator
$\Box$, $R$ is the Ricci scalar and $\boldsymbol{C}$ is the Weyl tensor.222We
use the bold font for tensors and their abstract indices. The regular font is
used for scalars, e.g., coordinates and tensor components. GR is recovered by
setting $\mathcal{F}_{i}(\Box)=0$, while IDG corresponds to non-polynomial
functions $\mathcal{F}_{i}$, which are typically chosen by demanding specific
properties of the linearized theory.

To fix a particular physically interesting IDG, let us consider the metric
perturbations $\boldsymbol{\gamma}$,
$\boldsymbol{g}=\bar{\boldsymbol{g}}+\boldsymbol{\gamma}$, around Minkowski
spacetime $\bar{\boldsymbol{g}}$. The linearized field equations in the Landau
gauge, ${\bar{\nabla}_{a}\gamma^{a}_{b}}=\frac{1}{4}\bar{\nabla}_{b}\gamma$,
read

$-\bar{\Box}\mathcal{E}_{2}(\bar{\Box}){\gamma}_{ab}^{\perp}=2\varkappa{T_{ab}^{\perp}}\;,\quad-\bar{\Box}\mathcal{E}_{0}(\bar{\Box})\gamma=-\tfrac{4}{3}\varkappa
T\;,$ (2)

where the symbol $\perp$ stands for the covariant transverse-traceless
projection of the symmetric rank-2 tensor and the analytic operators
$\mathcal{E}_{i}(\bar{\Box})$ are defined by

$\mathcal{E}_{i}(\bar{\Box})=1+2\varkappa\bar{\Box}\mathcal{F}_{i}(\bar{\Box})\;.$
(3)

Interesting choices of $\mathcal{F}_{i}(\Box)$ can be characterized by the
following conditions on $\mathcal{E}_{i}(\bar{\Box})$:

$\mathcal{E}_{i}(\bar{\Box})=e^{\mathcal{A}_{i}(\bar{\Box})}\;,\quad\lim_{\mathclap{\bar{\Box}\to\pm\infty}}\mathcal{A}_{i}(\bar{\Box})=+\infty\;,\quad\mathcal{A}_{i}\left(0\right)=0\;,$
(4)

where $\mathcal{A}_{i}$ are entire functions. The first condition guarantees
the absence of ghosts and other extra degrees of freedom (because the
exponential function is non-zero in $\mathbb{C}$); the second condition
represents the desired suppression of the propagator at short spacetime
distances (i.e., in the UV); and the last condition is necessary to recover GR
solutions at long spacetime distances (i.e., in the IR).

The simplest choice satisfying (4) is

$\mathcal{E}_{i}(\bar{\Box})=e^{\ell^{4}\bar{\Box}^{2}}\;.$ (5)

It has been demonstrated on various examples that (5) (and similar choices)
regularize singular GR solutions at this linearized level. The choice (5) also
fixes the non-local operators $\mathcal{F}_{i}(\Box)$ of IDG in the full
theory (1),

$\mathcal{F}_{i}(\Box)=\frac{e^{\ell^{4}\Box^{2}}-1}{2\varkappa\Box}\;.$ (6)

In what follows, we would like to demonstrate that the exact pp-wave solution
of such a theory are also regular in contrast to their singular GR
counterparts.

##  III Exact PP-wave solutions

The field equations of (1) reduce immensely for various subclasses of Kundt
spacetimes [16, 17]. One of the simplest cases is the pp-wave metric of type
N,

$\displaystyle\boldsymbol{g}=-\boldsymbol{\mathrm{d}}u\vee\big{(}\boldsymbol{\mathrm{d}}r+H\,\boldsymbol{\mathrm{d}}u\big{)}+\boldsymbol{\mathrm{d}}\rho^{2}+\rho^{2}\boldsymbol{\mathrm{d}}\varphi^{2}$
(7)

where $H=H(u,\rho,\varphi)$ is the only unknown function. Choosing this metric
ansatz and the pure-radiation energy-momentum tensor,
$\boldsymbol{T}={E}\boldsymbol{\mathrm{d}}u^{2}$, ${E}={E}(u,\rho,\varphi)$,
the field equations of the full theory reduce to

$\triangle\mathcal{E}_{2}(\triangle){H}=\varkappa{E}\;,$ (8)

where $\triangle$ is the Laplace operator on transverse 2-space
${\boldsymbol{\mathrm{d}}\rho^{2}+\rho^{2}\boldsymbol{\mathrm{d}}\varphi^{2}}$.
The similarity of (8) with the first equation of (2) reflects the fact that
the exact solutions of the full theory within this class of metrics are also
solutions of the linearized theory.333Note that this is not true for the
gyraton pp-wave metrics of more general type III [17] corresponding to the
spinning null matter. This feature is purely due to the simplicity of the
metric ansatz (7) — no approximation was assumed in the derivation of (8).

Consider a source located at ${\rho=0}$ in the transverse 2-space with an
arbitrary wave profile $w(u)$, ${\varkappa{E}=w(u)\delta(\rho)}$. Due to the
symmetry of the source, the field equation contain only derivatives with
respect to $\rho$, so the metric is characterized by one single-variable
function, $h=h(\rho)$,

$H=w(u)h(\rho)\;.$ (9)

The equation (8) can be solved for an arbitrary analytic choice of
$\mathcal{F}_{i}(\Box)$ using the Fourier transform in polar coordinates that
takes the form of the Hankel transform,444This integral usually contains a
diverging part that is constant in $\rho$ and of no physical significance. An
easy way to extract the convergent part is to replace
$J_{0}\to\partial_{\rho}J_{0}=-sJ_{1}$, evaluate the $s$-integral, and
calculate the primitive function in $\rho$.

$h=\int_{0}^{\infty}\\!\\!\\!ds\frac{J_{0}(s\rho)}{-s\mathcal{E}_{2}\left(-s^{2}\right)}\;,$
(10)

where $J_{0}$ is the Bessel function of the first kind.

The solutions for GR and IDG from the previous section read

$\displaystyle h_{\textrm{GR}}$ $\displaystyle=\log\rho+C\;,$ (11)
$\displaystyle h_{\textrm{IDG}}$
$\displaystyle=\tfrac{\sqrt{\pi}\rho^{2}}{2^{4}\ell^{2}}\,_{1}F_{3}\left(\tfrac{1}{2};1,\tfrac{3}{2},\tfrac{3}{2};\tfrac{\rho^{4}}{2^{8}\ell^{4}}\right)$
$\displaystyle\mathrel{\phantom{=}}-\tfrac{\rho^{4}}{2^{8}\ell^{4}}\,_{2}F_{4}\left(1,1;\tfrac{3}{2},\tfrac{3}{2},2,2;\tfrac{\rho^{4}}{2^{8}\ell^{4}}\right){+}C\;,$

where $C$ is a constant of no physical significance [19]. The metric
components of the GR pp-wave solution diverge at the location of source
($\rho=0$) because ${h_{\textrm{GR}}\propto\log\rho}$. On the other hand, the
metric components of the IDG solution are smooth at ${\rho=0}$ because
$h_{\textrm{IDG}}$ is analytic as a function of $\rho^{2}$, which follows from
analyticity of generalized hypergeometric functions ${}_{p}F_{q}$ and
$h_{\textrm{IDG}}^{(2n+1)}(0)=0$. As we will see, these properties of $h$ are
intimately related to the presence/absence of curvature singularities.

##  IV Geodesics and PP frames

Since all scalar curvature invariants vanish, we have to analyze the
presence/absence of nonscalar curvature singularities. Hence, we have to
investigate the geodesic motion of massive and massless test particles near
$\rho=0$. For an arbitrary geodesic
$\mathrm{x}(\lambda)=(u(\lambda),r(\lambda),\rho(\lambda),\varphi(\lambda))$,
the corresponding tangent vector $\boldsymbol{v}$ is given by

$\boldsymbol{v}=\tfrac{\boldsymbol{\mathrm{D}}\mathrm{x}}{\mathrm{d}\lambda}=u^{\prime}\boldsymbol{\partial}_{u}+r^{\prime}\boldsymbol{\partial}_{r}+\rho^{\prime}\boldsymbol{\partial}_{\rho}+\varphi^{\prime}\boldsymbol{\partial}_{\varphi}\;,$
(12)

where $\lambda$ is the affine parameter and ${}^{\prime}=d/d\lambda$.

The pp-wave spacetimes possess symmetries which generate conserved quantities
and simplify the integration of the geodesic equation. If we also assume that
the wave-profile function is constant, ${w(u)=\text{const.}}$, the metric (7)
is independent of $u$, $r$ and $\varphi$, which gives rise to three Killing
vectors $\boldsymbol{\partial}_{u}$, $\boldsymbol{\partial}_{r}$,
$\boldsymbol{\partial}_{\varphi}$ and consequently three constants of motion
$c_{(u)}$, $c_{(r)}$, $c_{(\varphi)}$. In addition to that, $\boldsymbol{v}$
fulfills the normalization condition with the normalization constant
$\epsilon$, hence

$\displaystyle\boldsymbol{v}^{\flat}\cdot\boldsymbol{\partial}_{u}$
$\displaystyle=-r^{\prime}-2whu^{\prime}=c_{(u)}\;,$ (13)
$\displaystyle\boldsymbol{v}^{\flat}\cdot\boldsymbol{\partial}_{r}$
$\displaystyle=-u^{\prime}=c_{(r)}\;,$
$\displaystyle\boldsymbol{v}^{\flat}\cdot\boldsymbol{\partial}_{\varphi}$
$\displaystyle=\rho^{2}\varphi^{\prime}=c_{(\varphi)}\;,$
$\displaystyle\boldsymbol{v}^{\flat}\cdot\boldsymbol{v}$
$\displaystyle=-2r^{\prime}u^{\prime}-2whu^{\prime 2}+\rho^{\prime
2}+\rho^{2}\varphi^{\prime 2}=\epsilon\;,$

where the ${\epsilon=-1}$ for timelike geodesics and ${\epsilon=0}$ for null
geodesics.

Let us now we determine the PP frames along all causal geodesics
($\epsilon=-1,0$) reaching ${\rho=0}$. We start with the timelike geodesics
($\epsilon=-1$) and identify the first frame vector with the timelike tangent
vector $\boldsymbol{v}$. The remaining frame vectors are chosen to satisfy the
orthonormality condition, i.e.,

$\displaystyle\boldsymbol{e}_{1}$ $\displaystyle=\boldsymbol{v}\;,$
$\displaystyle\boldsymbol{e}_{3}$
$\displaystyle=\tfrac{\rho^{\prime}}{u^{\prime}}\boldsymbol{\partial}_{r}+\boldsymbol{\partial}_{\rho}\;,$
(14) $\displaystyle\boldsymbol{e}_{2}$
$\displaystyle=\boldsymbol{v}-\tfrac{1}{u^{\prime}}\boldsymbol{\partial}_{r}\;,$
$\displaystyle\boldsymbol{e}_{4}$
$\displaystyle=\tfrac{\rho\varphi^{\prime}}{u^{\prime}}\boldsymbol{\partial}_{r}+\tfrac{1}{\rho}\boldsymbol{\partial}_{\varphi}\;.$

The only non-vanishing components of the TF Ricci tensor $\boldsymbol{S}$ and
Weyl tensor $\boldsymbol{C}$ of the metric (7) in the orthonormal frame (14)
are given by

$\displaystyle{S}_{11}$ $\displaystyle={S}_{22}={S}_{12}=w\Xi_{S}u^{\prime
2}\;,$ (15) $\displaystyle{C}_{1313}$
$\displaystyle={C}_{1323}={C}_{2323}=-{C}_{1414}$
$\displaystyle=-{C}_{1424}=-{C}_{2424}=w\Xi_{C}u^{\prime 2}\;,$

where

$\Xi_{S}=h^{\prime\prime}+\tfrac{1}{\rho}h^{\prime}\;,\quad\Xi_{C}=h^{\prime\prime}-\tfrac{1}{\rho}h^{\prime}\;.$
(16)

Note that the vectors $\boldsymbol{e}_{3}$ and $\boldsymbol{e}_{4}$ are PP
along the geodesic with the tangent vector $\boldsymbol{e}_{1}$ only if
${\varphi^{\prime}=0}$; otherwise they have to be properly rotated. However,
as we show bellow, $\varphi$ is constant for all geodesics passing through the
source in the IDG solution.

In the case of null geodesics ($\epsilon=0$) with the null tangent vector
$\boldsymbol{v}$, we construct the corresponding complex null frame

$\boldsymbol{l}=\boldsymbol{v}\;,\quad\boldsymbol{n}=\tfrac{1}{u^{\prime}}\boldsymbol{\partial}_{r}\;,\quad\boldsymbol{m}=\tfrac{1}{\sqrt{2}}(\boldsymbol{e}_{3}-i\boldsymbol{e}_{4})\;,$
(17)

where the tangent vector $\boldsymbol{v}$ is identified with the null frame
vector $\boldsymbol{l}$ and $\boldsymbol{e}_{3}$, $\boldsymbol{e}_{4}$ are
defined as in (14). Similarly as for the orthonormal frame (14), the null
frame (17) is PP along $\boldsymbol{l}$ only if $\varphi^{\prime}=0$. In the
standard notation of the Newman–Penrose formalism, the only non-vanishing
components of the TF Ricci tensor $\boldsymbol{S}$ and Weyl tensor
$\boldsymbol{C}$ of the metric (7) then read

$\Phi_{22}=\tfrac{w}{2}\Xi_{S}u^{\prime
2}\;,\quad\Psi_{4}=\tfrac{w}{2}\Xi_{C}u^{\prime 2}\;.$ (18)

##  V GR pp-wave solution exhibits curvature singularity

It is obvious that ${\Xi_{\text{S}}=0}$ for $\rho>0$ while $\Xi_{\text{C}}$
diverges at the location of the source, ${\rho=0}$, for the GR solution
$h_{\textrm{GR}}$. Therefore also the non-vanishing components of the Weyl
tensor $\boldsymbol{C}$ in PP frames along any causal geodesic reaching
$\rho=0$ necessarily diverge. To prove the presence of non-scalar curvature
singularity in the GR pp-wave solution, it only remains to show that the
boundary points ${\rho=0}$ can be reached by at least one particular causal
geodesic in a finite value of affine parameter.

An explicit timelike geodesic with vanishing angular momentum per unit mass in
the transverse 2-space, ${c_{(\varphi)}=0}$, can be obtained analytically by
integration of (13). Choosing the remaining constants of motion
${c_{(u)}=-1}$, ${c_{(r)}=-1}$, and the initial conditions ${u(0)=0}$,
${r(0)=0}$, ${\rho(0)=e^{1/4}}$, one gets

$\displaystyle u$ $\displaystyle=\lambda\;,\quad
r=2\left(\lambda-\rho\operatorname{erf}^{-1}\left(\lambda/\lambda_{0}\right)\right)\;,$
(19) $\displaystyle\rho$
$\displaystyle=\exp\left[{1}/{4}-\left(\operatorname{erf}^{-1}\left(\lambda/\lambda_{0}\right)\right)^{2}\right],\quad\varphi=\text{const.}\;,$

where we denoted ${\lambda_{0}=e^{1/4}\sqrt{\pi}/2}$. This geodesic ends up at
the location of the source, ${\rho=0}$ [in particular at
$(u,r,\rho,\varphi)=(\lambda_{0},2\lambda_{0},0,\textrm{const.})$], for the
finite value of proper time, ${\lambda=\lambda_{0}}$. Graphs of the functions
(19) are depicted in Fig. 1.

Figure 1: Coordinate description of the timelike geodesic (19) in the GR pp-
wave solution that reaches the source (${\rho=0}$) for the finite value of
${\lambda=\lambda_{0}\approx 1.14}$.

##  VI IDG pp-wave solution is regular

In the previous section, we studied only one particular geodesic as it is
sufficient to show the presence of nonscalar curvature singularities. In order
to verify the absence of such singularities in the IDG pp-wave solution, we
have to discuss all possible causal geodesics passing through the source,
$\rho=0$.

First, we notice that $\varphi(\lambda)$ is constant along such geodesics in
the IDG solution. If one assumes ${\varphi^{\prime}\neq 0}$, it follows from
(13) that ${c_{(\varphi)}\varphi^{\prime}=(c_{(\varphi)}/\rho)^{2}\to+\infty}$
on the source, ${\rho=0}$. Employing the constants of motion, $r^{\prime}$ and
$\rho^{2}$ can be eliminated in the normalization condition

$2wh_{\textrm{IDG}}c_{(r)}^{2}-2c_{(u)}c_{(r)}+\rho^{\prime
2}+c_{(\varphi)}\varphi^{\prime}=\epsilon\;.$ (20)

Since $h_{\textrm{IDG}}$ is bounded and $\rho^{\prime 2}$ is positive, there
is no way how to compensate ${c_{(\varphi)}\varphi^{\prime}\to+\infty}$ on the
source and therefore necessarily ${\varphi^{\prime}=0}$ implying
${c_{(\varphi)}=0}$. On the other hand, the vanishing of the constant of
motion ${c_{(\varphi)}=0}$ then ensures that ${\varphi^{\prime}=0}$ also for
${\rho\neq 0}$ and we can thus conclude that $\varphi$ is constant along any
geodesic passing through the source.

As an immediate consequence of this statement, the orthonormal frame (14) and
the null frame (17) are PP along the corresponding geodesics. It turns out
that the components of the curvature tensors ($\boldsymbol{S}$ and
$\boldsymbol{C}$) in PP frames along any timelike and null geodesics passing
through the source are given exactly by (15) and (18), respectively. These
components are finite everywhere since $\Xi_{S}$, $\Xi_{C}$ are bounded and
$u^{\prime}=-c_{(r)}$ for IDG pp-waves solution $h_{\textrm{IDG}}$.

##  VII Conclusions

We have shown that IDG can resolve nonscalar curvature singularities in exact
pp-wave solutions. Thus, the nonlocality of IDG can prevent arbitrarily large
tidal forces that would be experienced by an observer in finite time in these
spacetimes. This is the first explicit exact demonstration of the singularity
resolution due to nonlocality where the curvature of the exact GR and IDG
solutions with the same matter content were compared side by side. These
results were achieved by analyzing components of the curvature tensors in PP-
frames along timelike and null geodesics reaching the source with the finite
value of affine parameter.

We believe that our study will stimulate further search for exact solutions of
IDG and examination of the presence/absence of (scalar or nonscalar) curvature
singularities due to nonlocality. For example, a similar calculation can be
repeated for gyraton solutions generated by spinning null matter [17] or Kundt
geometries with non-zero cosmological constant [15, 16]. The former would be
particularly interesting because the exact IDG gyratons [17] differ from the
linearized IDG gyratons [20]. A careful analysis of spacetime singularities in
nonlocal gravitational theories (and the associated matter distribution)
should also reduce the lack of clarity on this topic in the current literature
[21].

## Acknowledgements

I.K. was supported by Netherlands Organization for Scientific research (NWO)
grant no. 680-91-119 and Primus grant PRIMUS/23/SCI/005 from Charles
University. T.M. acknowledges the support of the Czech Academy of Sciences
(RVO 67985840) and the Czech Science Foundation GAČR grant no. GA19-09659S.

## References

  * Stelle [1977] K. S. Stelle, Phys. Rev. D 16, 953 (1977).
  * Stelle [1978] K. S. Stelle, Gen. Rel. Grav. 9, 353 (1978).
  * Burzillà  _et al._ [2021] N. Burzillà, B. L. Giacchini, T. d. P. Netto, and L. Modesto, Eur. Phys. J. C 81, 462 (2021), arXiv:2012.11829 [gr-qc] .
  * Krasnikov [1987] N. V. Krasnikov, Theor. Math. Phys. 73, 1184 (1987).
  * Kuzmin [1989] Y. V. Kuzmin, Soviet Journal of Nuclear Physics-Ussr 50, 1011 (1989).
  * Tomboulis [1997] E. T. Tomboulis, Superrenormalizable gauge and gravitational theories (1997), arXiv:hep-th/9702146 .
  * Biswas  _et al._ [2006] T. Biswas, A. Mazumdar, and W. Siegel, JCAP 2006 (03), 009, arXiv:hep-th/0508194 .
  * Modesto [2012] L. Modesto, Phys. Rev. D 86, 044005 (2012), arXiv:1107.2403 [hep-th] .
  * Biswas  _et al._ [2012] T. Biswas, E. Gerwick, T. Koivisto, and A. Mazumdar, Phys. Rev. Lett. 108, 031101 (2012), arXiv:1110.5249 [gr-qc] .
  * Biswas  _et al._ [2014] T. Biswas, A. Conroy, A. S. Koshelev, and A. Mazumdar, Class. Quant. Grav. 31, 015022 (2014), arXiv:1308.2319 [hep-th] .
  * Hervik  _et al._ [2014] S. Hervik, V. Pravda, and A. Pravdova, Class. Quant. Grav. 31, 215005 (2014), arXiv:1311.0234 [gr-qc] .
  * Biswas  _et al._ [2010] T. Biswas, T. Koivisto, and A. Mazumdar, JCAP 2010 (11), 008, arXiv:1005.0590 [hep-th] .
  * Koshelev  _et al._ [2018] A. S. Koshelev, K. Sravan Kumar, and A. A. Starobinsky, JHEP 2018 (3), 071, arXiv:1711.08864 [hep-th] .
  * Kilicarslan [2019] E. Kilicarslan, Phys. Rev. D 99, 124048 (2019), arXiv:1903.04283 [gr-qc] .
  * Dengiz  _et al._ [2020] S. Dengiz, E. Kilicarslan, I. Kolář, and A. Mazumdar, Phys. Rev. D 102, 044016 (2020), arXiv:2006.07650 [gr-qc] .
  * Kolář  _et al._ [2021] I. Kolář, T. Málek, and A. Mazumdar, Phys. Rev. D 103, 124067 (2021), arXiv:2103.08555 [gr-qc] .
  * Kolář  _et al._ [2022] I. Kolář, T. Málek, S. Dengiz, and E. Kilicarslan, Phys. Rev. D 105, 044018 (2022), arXiv:2107.11884 [gr-qc] .
  * Frolov and Zelnikov [2016] V. P. Frolov and A. Zelnikov, Phys. Rev. D 93, 064048 (2016), arXiv:1509.03336 [hep-th] .
  * Griffiths and Podolský [2009] J. B. Griffiths and J. Podolský,  _Exact Space-Times in Einstein’s General Relativity_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2009).
  * Boos  _et al._ [2020] J. Boos, J. Pinedo Soto, and V. P. Frolov, Phys. Rev. D 101, 124065 (2020), arXiv:2004.07420 [gr-qc] .
  * Buoninfante  _et al._ [2018] L. Buoninfante, A. S. Koshelev, G. Lambiase, J. a. Marto, and A. Mazumdar, JCAP 2018 (06), 014, arXiv:1804.08195 [gr-qc] .