# Matrix Multiplication Reductions

Ashish Gola<EMAIL_ADDRESS>Simon Fraser University Igor Shinkar
<EMAIL_ADDRESS>Simon Fraser University Harsimran Singh
<EMAIL_ADDRESS>Simon Fraser University

###### Abstract

In this paper we study a worst case to average case reduction for the problem
of matrix multiplication over finite fields. Suppose we have an efficient
_average case_ algorithm, that given two random matrices $A,B$ outputs a
matrix that has a non-trivial correlation with their product $A\cdot B$. Can
we transform it into a _worst case_ algorithm, that outputs the correct answer
for all inputs without incurring a significant overhead in the running time?
We present two results in this direction.

Two-sided error in the high agreement regime

We begin with a brief remark about a reduction for high agreement algorithms,
i.e., an algorithm which agrees with the correct output on a large (say
$>0.9$) fraction of entries, and show that the standard self-correction of
linearity allows us to transform such algorithms into algorithms that work in
worst case.

One-sided error in the low agreement regime

Focusing on average case algorithms with one-sided error, we show that over
$\mathbb{F}_{2}$ there is a reduction that gets an $O(T)$ time _average case_
algorithm that given a random input $A,B$ outputs a matrix that agrees with
$A\cdot B$ on at least $51\%$ of the entries (i.e., has only a slight
advantage over the trivial algorithm), and transforms it into an
$\widetilde{O}(T)$ time _worst case_ algorithm, that outputs the correct
answer for _all inputs_ with high probability.

##  1 Introduction

The problem of efficiently multiplying two matrices has been extensively
studied for decades. Improving on the straightforward $O(n^{3})$ time
algorithm, Strassen’s algorithm [Str69] computes the product of two matrices
in time $O(n^{\log_{2}7}=n^{2.807})$, and it is perhaps the most widely used
in practice. Since then, a long and exciting line of research ([Pan78],
[BCRL79], [Sch81], [Rom82], [Str86], [CW90], [Sto10], [Wil12], [LG14], [AV21])
has led to a significant improvement of the value of the optimal exponent of
the running time for matrix multiplication problem. The fastest algorithm
known today is due to Duan, Wu, and Zhou [DWZ23], and its running time is
$O(n^{2.371866})$.

Worst-case to average-case reductions serve as a means to convert algorithms
that output correct answers on a fraction of inputs into algorithms with
correct outputs on all possible inputs. These reductions can be viewed from
two different perspectives. From the hardness point of view, they can be used
to show that a problem maintains its hardness even in the average case. From
the algorithmic side, they provide a framework for developing worst-case
algorithms, by first designing weak algorithms with average case guarantees,
and then transforming them into algorithms which work on all outputs.

In this paper, we study the following variant of a worst-case to average-case
reduction for the matrix multiplication problem. Suppose we have an efficient
algorithm that given two random matrices $A,B\in\mathbb{F}^{n\times n}$,
computes a matrix $C\in\mathbb{F}^{n\times n}$ that agrees with the product
$A\cdot B$ on a large fraction of the entries of the matrix. Can we transform
such an algorithm into one that computes $A\cdot B$ correctly for all entries
of the output matrix without incurring a significant overhead in the running
time?

More formally, we define the _agreement_ between two matrices as the fraction
of entries on which the two matrices agree.

######  Definition 1.1.

Let $\mathbb{F}$ be a field, and let $A,B\in\mathbb{F}^{n\times n}$ be two
matrices. We define agreement between $A$ and $B$, denoted by ${\rm
agr}(A,B)$, as the fraction of entries $(i,j)$ on which $A_{i,j}=B_{i,j}$,
i.e.,

${\rm agr}(A,B)=\frac{|\\{(i,j):A_{i,j}=B_{i,j}\\}|}{n^{2}}\enspace.$

Then, our goal can be stated as the task of transforming an algorithm that on
a random input $A,B$ outputs a matrix $C$ such that ${\rm
agr}(C,AB)\geq\alpha$ for some parameter $\alpha\in[0,1]$ into an algorithm
that solves the matrix multiplication problem _correctly on all inputs_.

We present two results in this direction. Both results consider the matrix
multiplication problem over finite fields.

#### High agreement regime with two-sided error:

We show that any algorithm that solves the matrix multiplication problem
correctly on a high fraction of the coordinates, can be converted into a worst
case algorithm. Specifically, we prove the following theorem.

######  Theorem 1.2.

Fix a finite field $\mathbb{F}$. Let $\alpha\in(0,1/8)$. Let ${\mathrm{A}LG}$
be an algorithm that gets as input two matrices $A,B\in\mathbb{F}^{n\times
n}$, runs in time $T(n)$, and outputs a matrix $C\in\mathbb{F}^{n\times n}$.
Suppose that

$\mathbb{E}_{A,B\in\mathbb{F}^{n\times n}}[{\rm agr}(ALG(A,B),A\cdot
B)]>1-\alpha\enspace.$

Then, there is an algorithm ${\mathrm{A}LG}^{*}$ that gets as input two
matrices $A,B\in\mathbb{F}^{n\times n}$, runs in time $O(T(n)\cdot\log(n))$,
and outputs a matrix $C\in\mathbb{F}^{n\times n}$ such that _for all_ $A,B$ it
holds that

$\Pr[{\mathrm{A}LG}^{*}(A,B)=A\cdot B]>1-1/n\enspace,$

where the randomness is only over the internal coins of ${\mathrm{A}LG}^{*}$.

The proof of this result relies on rather standard ideas, and essentially uses
the self-correction of linear functions [BLR90].

#### Low agreement with one-sided error:

For this result, we restrict our discussion to the finite field
$\mathbb{F}_{2}$. Note that it is trivial to design an $O(n^{2})$ time
algorithm such that $\mathbb{E}_{A,B\in\mathbb{F}_{2}^{n\times n}}[{\rm
agr}(ALG(A,B),A\cdot B)]\geq 1/2$. Indeed, the algorithm can simply output 0
in all entries irrespective of the input. Alternatively, the algorithm can
output a random 0/1 matrix. Hence, it is natural to ask whether it is possible
to obtain a better-than-1/2 algorithm for the matrix multiplication over
$\mathbb{F}_{2}$.

Below we show that in the special case of _one-sided error_ approximation, any
better-than-1/2 approximation $O(T)$ time algorithm can be transformed into a
worst case algorithm with running time $\widetilde{O}(T)$. Formally, we prove
the following theorem.

######  Theorem 1.3.

Let ${\mathrm{A}LG}$ be an algorithm that gets input two matrices
$A,B\in\mathbb{F}_{2}^{n\times n}$, runs in time $T(n)$, and outputs a matrix
$C\in\mathbb{F}_{2}^{n\times n}$. Let $\delta>0$, and suppose that

  * •

$\mathbb{E}_{A,B\in\mathbb{F}^{n\times n}}[{\rm agr}(ALG(A,B),A\cdot B)]\geq
1/2+\delta$.

  * •

If $(AB)_{i,j}=0$, then ${\mathrm{A}LG}(A,B)_{i,j}=0$.

Then, there is an algorithm ${\mathrm{A}LG}^{*}$ that gets as input two
matrices $A,B\in\mathbb{F}^{n\times n}$, runs in time $\widetilde{O}(T(n))$,
and outputs a matrix $C\in\mathbb{F}^{n\times n}$ such that for all $A,B$ it
holds that

$\Pr[{\mathrm{A}LG}^{*}(A,B)=A\cdot B]>1-1/n,$

where the randomness is only over the internal coins of ${\mathrm{A}LG}^{*}$.

######  Remark 1.4.

Below we make several comments about Theorem 1.3.

  1. 1

Note that the conditions of the theorem can be written equivalently as
follows.

     * •

$\Pr_{\begin{subarray}{c}A,B\in\mathbb{F}_{2}^{n\times n}\\\
i,j\in[n]\end{subarray}}[ALG(A,B)_{i,j}=1]\geq\delta$.

     * •

If $(AB)_{i,j}=0$, then ${\mathrm{A}LG}(A,B)_{i,j}=0$.

  2. 2

The notion of algorithms with one-sided error is typically studied in the
context of randomized algorithms, e.g., relating to the classes $\mathcal{RP}$
(and $co\mathcal{RP}$), where the guarantee is that for _every NO input_ the
algorithm outputs the correct answer with probability 1, and for _every YES
input_ it is correct with probability at least 2/3. The error model in Theorem
1.3 is different, as we consider algorithms that are correct on _random
inputs_ on all output 0-bits, and on at least some $\alpha$-fraction of
1-bits.

  3. 3

We remark that the standard methods of self-correcting linear functions work
in the high agreement regime, but fail when the average case guarantee is low.
We apply the techniques from additive combinatorics developed in [AGGS22],
particularly a version of the probabilistic Bogolyubov-Ruzsa Lemma, to perform
a self-correction procedure which helps in this regime.

  4. 4

Our proof of Theorem 1.3 assumes that ${\mathrm{A}LG}$ is deterministic. It is
rather straightforward to extend the proof and allow it to be randomized, by
appropriately modifying the sets of good inputs ($X_{i,j}$ and $Y_{i,j}^{A}$)
to account for the randomness of the algorithm.

###  1.1 Related work

The study of average-case complexity began with Levin’s work [Lev86], followed
by subsequent works like [BDCGL92]. A substantial body of research (e.g.,
[IL90], [Imp11] and related references) identified numerous barriers in
formulating worst-case to average-case reductions for NP-complete problems.
For a comprehensive overview of this subject, see the classical surveys by
Impagliazzo [Imp95], Bogdanov and Trevisan [BT06] and Goldreich [Gol11].

Asadi et al. [AGGS22, AGG+] presented a new framework for carrying out
efficient worst-case to average case reductions for various fundamental
problems. Particularly, for the problem of matrix multiplication, they proved
that if there exists an $O(T(n))$ time algorithm $M$ for matrix multiplication
which computes the correct output on an $\epsilon$ fraction of inputs, then
there exists a randomized algorithm $M^{\prime}$ which computes the correct
output on all inputs, running in time $O(\exp(O(\log^{5}(1/\epsilon))\cdot
T(n))$. The proof relied on additive combinatorial techniques and used the
probabilistic quasi-polynomial Bogolyubov-Ruzsa Lemma.

Hirahara and Shimizu [HS23] improved the $\exp(O(\log^{5}(1/\epsilon)))$
overhead to an $\widetilde{O}(1/\epsilon)$ factor. Their idea involved
dividing the output matrix into smaller blocks and using the Direct-Product
Theorem in a black-box manner.

The aforementioned papers assume that we have access to an algorithm which
gives a fully correct output on some fraction of the inputs, i.e., for these
inputs _all entries_ in the output matrix are correct. The setting presented
in this paper, where the output of the given algorithm is not fully correct,
seems to differ significantly from the works mentioned above. In particular,
we do not see how to apply the Direct-Product theorem to our setting of the
problem.

###  1.2 Open problems

We mention the following two problems that are left open in this work.

#### Low agreement with two-sided error:

Is it possible to transform a two-sided error algorithm over $\mathbb{F}_{2}$
with a low agreement guarantee into a worst case algorithm. That is, given an
$O(T(n))$ time algorithm $ALG$ with the guarantee
$\mathbb{E}_{A,B\in\mathbb{F}^{n\times n}}[{\rm agr}(ALG(A,B),A\cdot
B)]>1/2+\delta$, can we convert it into an algorithm that correctly outputs
the correct answer on all inputs and has running time $\widetilde{O}(T(n))$?

#### Generalizing over finite fields:

Extend Theorem 1.3 in a meaningful way to work over any finite field.

##  2 Preliminaries

For a positive integer $n$ we define $[n]=\\{0,1,\dots,n-1\\}$. We index the
coordinates of our matrices starting from 0 rather than 1, which is typically
more standard. We refer to the element in the row $i$ and column $j$ of the
matrix $A$ as $A_{i,j}$.

We define a notion of _row-shift_ (or _row-rotation_) and _column-shift_ as
follows.

######  Definition 2.1.

Given a matrix $A\in\mathbb{F}^{n\times m}$, $0\leq\pi\leq n-1$, and
$0\leq\sigma\leq m-1$, define $A^{\pi,\sigma}$ to be the matrix obtained from
$A$ by cyclically rotating all its rows downwards by $\pi$ units and all its
columns rightwards by $\sigma$ units, that is,

$(A^{\pi,\sigma})_{i,j}=A_{(i-\pi)\bmod n,(j-\sigma)\bmod m}\enspace$

The following proposition is immediate from the definition above.

######  Proposition 2.2.

For any $A,B,C\in\mathbb{F}^{n\times n}$ and any $\pi,\sigma$ we have $AB=C$
if and only if $A^{\pi,0}\cdot B^{0,\sigma}=C^{\pi,\sigma}$.

###  2.1 Additive Combinatorics Tools

We now present the additive combinatorics toolkit which will be useful in the
worst-case to average-case reduction for the low agreement regime with one-
sided error.

For a set $A\subseteq\mathbb{F}_{2}^{n}$, let
$1_{A}\colon\mathbb{F}_{2}^{n}\to\\{0,1\\}$ denote the indicator function of
A. The Fourier expansion of a function
$f\colon\mathbb{F}_{2}^{n}\to\mathbb{C}$ is given by
$f(x)=\sum_{r\in\mathbb{F}_{2}^{n}}\hat{f}(r)\cdot\chi_{r}(x)$, where
$\chi_{r}(x)=(-1)^{\left\langle x,r\right\rangle}$, and the Fourier
coefficients of $f$ are defined as $\hat{f}(r)=\left\langle
f,\chi_{r}\right\rangle=\mathbb{E}_{x}[f(x)\cdot\chi_{r}(x)]$. _Parseval’s
identity_ says that
$\sum_{r\in{\mathbb{F}_{2}}^{n}}{\widehat{1_{A}}(r)}^{2}=\left\langle
1_{A},1_{A}\right\rangle=\alpha$, where $\alpha$ is the density of $A$.

Define
${\mathrm{Spec}}_{\gamma}(A)=\\{r\in\mathbb{F}_{2}^{n}:{\left|\widehat{1_{A}}(r)\right|}\geq\gamma\\}$.
Below we state Chang’s lemma, which describes a certain structure of
${\mathrm{Spec}}_{\gamma}(A)$.

######  Lemma 2.3 (Chang’s Theorem [Cha02]).

Let $A\subseteq\mathbb{F}^{n}$ be a set of size
${\left|A\right|}=\alpha\cdot{\left|\mathbb{F}\right|}^{n}$, and let
$\gamma>0$. Then

$\dim(\mathsf{span}({\mathrm{Spec}}_{\gamma\alpha}(A)))\leq
O\left(\frac{\log(1/\alpha)}{\gamma^{2}}\right)\enspace.$

Recall that the subset sum of two sets $A$ and $B$ is defined as
$A+B=\\{a+b:a\in A\;,b\in B\\}$. Analogously, we define $tA=A+A+\dots+A$ ($t$
times) as $tA=\\{a_{1}+a_{2}\dots+a_{t}:a_{1},a_{2},\dots,a_{t}\in A\\}$. The
following lemma says that for an arbitrary set $A\subseteq\mathbb{F}^{n}$, the
sumset $tA$ contains a large affine subspace.

######  Lemma 2.4 (Probabilistic quasi-polynomial Bogolyubov-Ruzsa lemma).

Let $A\subseteq\mathbb{F}_{2}^{n}$ be a set of size
${\left|A\right|}=\alpha\cdot 2^{n}$, for some $\alpha\in(0,1]$, and let
$t\geq 3$ be an integer. Then, $tA$ contains an _affine_ subspace
$V\subseteq\mathbb{F}_{2}^{n}$ of dimension $\dim(V)\geq n-O(\log(1/\alpha))$
such that for all $v\in V$ it holds that

$\Pr_{a_{1},a_{2},..,a_{t-1}\in\mathbb{F}_{2}^{n}}[a_{1},a_{2},a_{3},..,a_{t}\in
A]\geq{\alpha}^{t}\left(1+\frac{1}{2^{t-2}}\right)-\frac{{\alpha}^{t-1}}{2^{t-2}}\enspace,$

where $a_{t}=v-a_{1}-a_{2}-..-a_{t-1}$.

In particular, if $t>\log_{2}(1/\alpha)+2$, then $tA$ contains an _affine_
subspace $V\subseteq\mathbb{F}_{2}^{n}$ of dimension $\dim(V)\geq n-k$, for
$k=O(\log(1/\alpha))$, such that for all $v\in V$ it holds that

$\Pr_{\begin{subarray}{c}a_{1},a_{2},..,a_{t}\in\mathbb{F}_{2}^{n}\\\
v=\sum_{i=1}^{t}a_{i}\end{subarray}}[a_{1},a_{2},a_{3},..,a_{t}\in
A]\geq(\alpha/2)^{t}\enspace.$

Below we prove Lemma 2.4 only for odd values of $t$, which is slightly more
complicated than the case of even $t$. After the proof, we remark how to
modify the proof to work for even $t$’s.

###### Proof.

Let $A\subseteq\mathbb{F}_{2}^{n}$ be a set of size
${\left|A\right|}=\alpha\cdot{\left|\mathbb{F}\right|}^{n}$, for some
$\alpha\in(0,1]$. Consider the set

$R={\mathrm{Spec}}_{\alpha/2}\setminus\\{0\\}=\\{r\in
F_{2}^{n}\setminus\\{0\\}:{\left|\widehat{1_{A}}(r)\right|}>\frac{\alpha}{2}\\}$

Next we define an affine subspace $V=\\{v\in\mathbb{F}_{2}^{n}:\left\langle
v,r\right\rangle=s_{r}\ \forall r\in R\\}$ for some $s_{r}\in\\{0,1\\}$ to be
defined later, and claim that $V$ satisfies the conclusions of Lemma 2.4. We
will need the following two claims.

######  Claim 2.5.

For all $r\in R$ there exists $s_{r}\in\\{0,1\\}$ such that (1) $\sum_{r\in
R}\widehat{1_{A}}(r)^{t}\cdot(-1)^{s_{r}}\geq 0$ and (2) if $r^{*}\in R$ is a
linear combination $r^{*}=\sum_{r\in R}c_{r}\cdot r$ of vectors in $R$ (with
$c_{r}\in\mathbb{F}_{2}$), then $s_{r^{*}}=\sum_{r\in R}c_{r}\cdot
s_{r}\pmod{2}$.

###### Proof.

Let $R^{\prime}$ be a maximal subset of $R$ of linearly independent vectors.
Choose $s_{r^{\prime}}\in\\{0,1\\}$ independently with probability 0.5 each
for every $r^{\prime}\in R^{\prime}$. Now any $r\in R\setminus R^{\prime}$,
can be expressed as a linear combination
$r=\sum_{r^{\prime}}c_{r^{\prime}}\cdot r^{\prime}$ of vectors in $R^{\prime}$
with $c_{r^{\prime}}\in\\{0,1\\}$ define
$s_{r}=\sum_{r^{\prime}}c_{r^{\prime}}\cdot s_{r^{\prime}}$. It is immediate
to verify that condition (2) is satisfied.

In order to satisfy condition (1) note that by linearity of expectation
$\mathbb{E}[\sum_{r\in R}\widehat{1_{A}}(r)^{t}\cdot(-1)^{s_{r}}]=0$, and
hence there exists a choice of $(s_{r})_{r\in R}$ such that $\sum_{r\in
R}\widehat{1_{A}}(r)^{t}\cdot(-1)^{s_{r}}\geq 0$, as required. ∎

######  Claim 2.6.

We have $\sum_{r\not\in R,r\neq
0}{{\left|\widehat{1_{A}}(r)\right|}}^{t}\leq(\alpha/2)^{t-2}(\alpha-{\alpha}^{2})$.

###### Proof.

For $t\geq 3$, it holds that

$\displaystyle\sum_{r\not\in R,r\neq
0}{{\left|\widehat{1_{A}}(r)\right|}}^{t}\leq\max_{r\not\in R,r\neq
0}{{\left|\widehat{1_{A}}(r)\right|}}^{t-2}\sum_{r\not\in R,r\neq
0}{{\left|\widehat{1_{A}}(r)\right|}}^{2}$
$\displaystyle\leq(\alpha/2)^{t-2}\sum_{r\in{\mathbb{F}_{2}}^{n}\setminus\\{0\\}}{\widehat{1_{A}}(r)}^{2}$
$\displaystyle<(\alpha/2)^{t-2}(\alpha-{\alpha}^{2})\enspace.$

∎

Define an affine subspace $V=\\{v\in\mathbb{F}_{2}^{n}:\left\langle
v,r\right\rangle=s_{r}\ \forall r\in R\\}$, where $s_{r}\in\\{0,1\\}$ is from
2.5. Note that if the vectors in $R$ are linearly dependent, then the second
condition of 2.5 guarantees that we can define
$V=\\{v\in\mathbb{F}^{n}:\left\langle
v,r^{\prime}\right\rangle=s_{r^{\prime}}\ \forall r^{\prime}\in R^{\prime}\\}$
for a maximal set $R^{\prime}\subset R$ of linearly independent vectors in
$R$, and the remaining constraints will be satisfied by linearity. Then,
according to Lemma 2.3 we have

$\dim(V)\geq n-O\left(\log(1/\alpha)\right)\enspace.$

Using the two claims above, and noting that
$\Pr_{a_{1},a_{2},..,a_{t-1}\in\mathbb{F}^{n}}[a_{1},a_{2},a_{3},..,a_{t}\in
A]=1_{A}*1_{A}*..*1_{A}(v)$ ($t$ times), for any $v\in V$ we have

$\displaystyle\Pr_{a_{1},a_{2},..,a_{t-1}\in\mathbb{F}^{n}}[a_{1},a_{2},a_{3},..,a_{t}\in
A]$ $\displaystyle=1_{A}*1_{A}*..*1_{A}(v)$
$\displaystyle=\sum_{r\in\mathbb{F}^{n}}{\widehat{1_{A}}(r)}^{t}\chi_{r}(v)$
$\displaystyle={\widehat{1_{A}}(0)}^{t}+\sum_{r\in
R}{\widehat{1_{A}}(r)}^{t}\chi_{r}(v)+\sum_{r\not\in R,r\neq
0}{\widehat{1_{A}}(r)}^{t}\chi_{r}(v)$
$\displaystyle\geq{\alpha}^{t}+\sum_{r\in
R}{\widehat{1_{A}}(r)}^{t}\cdot(-1)^{s_{r}}-(\alpha/2)^{t-2}(\alpha-{\alpha}^{2})$
$\displaystyle\geq{\alpha}^{t}+0-(\alpha/2)^{t-2}(\alpha-{\alpha}^{2})$
$\displaystyle={\alpha}^{t}\left(1+\frac{1}{2^{t-2}}\right)-\frac{{\alpha}^{t-1}}{2^{t-2}}\enspace.$

In particular, for $t>\log_{2}(1/\alpha)+2$, we have

$\displaystyle\Pr_{a_{1},a_{2},..,a_{t-1}\in\mathbb{F}^{n}}[a_{1},a_{2},a_{3},..,a_{t}\in
A]$
$\displaystyle\geq{\alpha}^{t}\left(1+\frac{1}{2^{t-2}}\right)-\frac{{\alpha}^{t-1}}{2^{t-2}}$
$\displaystyle\geq{\alpha}^{t}\left(1+\frac{1}{2^{t-2}}\right)-{\alpha}^{t}$
$\displaystyle\geq(\alpha/2)^{t}\enspace,$

as required. ∎

######  Remark 2.7.

For even values of $t$ the lemma is slightly easier. Specifically, since
${\widehat{1_{A}}(r)}^{t}$ is always non-negative, we can take $s_{r}=0$ in
2.5, and the rest of the proof works the same.

##  3 High Agreement with Two-Sided Error

In this section, we prove Theorem 1.2. Specifically, we show that if there
exists an algorithm which, given two matrices $A,B\in\mathbb{F}^{n\times n}$,
runs in time $T(n)$ and correctly computes their product on a large fraction
of all entries of output on average, then there exists another algorithm that
runs in $\widetilde{O}(T(n))$ time and correctly computes their product on all
entries of output. The proof essentially uses the self-correction of linearity
[BLR90].

See 1.2

###### Proof.

Given the algorithm ${\mathrm{A}LG}$ as in the assumption of the theorem, we
design ${\mathrm{A}LG}^{*}$ as follows.

Input: $A,B\in\mathbb{F}^{n\times n}$, $ALG$

Output: $A\cdot B$

1

2Let $k=O(\log(n))$

3for  _$r=0$ to $k$_ do

4 Generate two random matrices $R,S\in\mathbb{F}^{n\times n}$

5 Select two random variables $\pi,\sigma\in[n]$ independently

6
$M={\mathrm{A}LG}((A+R)^{\pi,0},(B+S)^{0,\sigma})-{\mathrm{A}LG}(R^{\pi,0},(B+S)^{0,\sigma})-{\mathrm{A}LG}((A+R)^{\pi,0},S^{0,\sigma})+{\mathrm{A}LG}(R^{\pi,0},S^{0,\sigma})$

7 Let $C_{r}=M^{n-\pi,n-\sigma}$

Define the matrix $C\in\mathbb{F}^{n\times n}$ by taking the majority vote of
all $C_{r}$ in each coordinate.

Algorithm 1 Approximation for High Agreement Matrix Multiplication Algorithms

#### Correctness:

Consider an entry $(i,j)$ in the output matrix $A\cdot B$. In each iteration
we call ${\mathrm{A}LG}$ four times, and in each of the calls the input is
distributed uniformly in $\mathbb{F}^{n\times n}$. Furthermore, since
$\pi,\sigma$ are chosen uniformly, it follows that $(i+\pi,j+\sigma)$ are
distributed uniformly. Therefore, the probability that that all the four calls
of ${\mathrm{A}LG}$ produce the correct answer in this entry is at least
$1-4\alpha$. Therefore, for each repetition $r$, we have

$\Pr[(C_{r})_{i,j}=(A\cdot B)_{i,j}]\geq 1-4\alpha\enspace.$

By Chernoff bound, the probably that the majority vote of the $k$ repetition
will produce an incorrect answer is upper bounded by

$\Pr[(C_{r})_{i,j}=(A\cdot B)_{i,j}]\geq\exp\left(-\Omega((1-4\alpha-1/2)\cdot
k)\right)<1/n^{3}\enspace,$

Here we make the assumption that $\alpha$ is bounded below $1/8$.

Hence, the probability that a particular entry $(i,j)$ is incorrect after $k$
iterations is at most $n^{-3}$. By union bound over all entries, the
probability that at least one entry is incorrect in the output matrix is at
most $n^{2}\cdot n^{-3}=1/n$.

#### Running time:

The total running time is dominated by $O(\log(n))$ invocations of
${\mathrm{A}LG}$, and hence, the runtime of ${\mathrm{A}LG}^{*}$ is
$O(T(n)\cdot\log(n))$. ∎

##  4 Low Agreement with One-Sided Error

In this section we prove Theorem 1.3. We restate the theorem here for
convenience.

See 1.3

Before proving the theorem, we need some definitions. We start by defining the
notion of a good coordinate. We say a coordinate $(i,j)\in[n]\times[n]$ is
good, if $ALG$ returns $1$ at the entry $(i,j)$ for more than ${\delta}/{2}$
fraction of possible inputs.

######  Definition 4.1.

Denote by $G$ the set of _good_ coordinates, defined as

$G=\\{(i,j)\in[n]\times[n]:\Pr_{A,B\in\mathbb{F}_{2}^{n\times
n}}[ALG(A,B)_{i,j}=1]>\delta/2\\}\enspace.$

The following claim is immediate from the definition and the assumptions of
the theorem.

######  Claim 4.2.

$|G|\geq(\delta/2)\cdot n^{2}$.

###### Proof.

Let $p_{i,j}=\Pr_{\begin{subarray}{c}A,B\in\mathbb{F}_{2}^{n\times
n}\end{subarray}}[ALG(A,B)_{i,j}=1]$. Note that by the assumptions of Theorem
1.3, we have $\mathbb{E}_{i,j}[p_{i,j}]\geq\delta$. Note that

$\delta\leq\mathbb{E}_{i,j\in[n]\times[n]}[p_{i,j}]\leq\Pr_{i,j}[(i,j)\in
G]\cdot 1+\Pr_{i,j}[(i,j)\not\in G]\cdot(\delta/2)\leq\Pr_{i,j}[(i,j)\in
G]\cdot 1+1\cdot(\delta/2)\enspace,$

and hence $\Pr[(i,j)\in G]\geq\delta/2$, as required. ∎

Next, we define the notion of good input matrices with respect to a good
coordinate.

######  Definition 4.3.

For a coordinate (i,j), define $X_{i,j}$ as follows.

$X_{i,j}=\\{A:\Pr_{B}[{\mathrm{A}LG}(A,B)_{i,j}=1]\geq\delta/4\\}\enspace.$

Given a coordinate $(i,j)$ and a matrix $A$, define $Y_{i,j}^{A}$ to be the
set of matrices $B$ for which ${\mathrm{A}LG}$ returns 1 at the entry $(i,j)$.
That is,

$Y_{i,j}^{A}=\\{B:{\mathrm{A}LG}(A,B)_{i,j}=1\\}\enspace.$

We make the following claim about the densities of $X_{i,j}$ and
$Y_{i,j}^{A}$.

######  Claim 4.4.

For any $(i,j)\in G$ it holds that $\Pr_{A\in\mathbb{F}^{n\times n}}[A\in
X_{i,j}]\geq\delta/4$. Furthermore, if $A\in X_{i,j}$ then
$\Pr_{B\in\mathbb{F}^{n\times n}}[B\in Y_{i,j}^{A}]\geq\delta/4$.

###### Proof.

Fix a good coordinate $(i,j)\in G$, and for each $A\in\mathbb{F}^{n\times n}$
let $p_{A}=\Pr_{B\in\mathbb{F}^{n\times n}}[{\mathrm{A}LG}(A,B)_{(i,j)}=1]$.
From the definition of $G$ we have $\mathbb{E}_{A}[p_{A}]\geq\delta/2$.

$\delta/2\leq\mathbb{E}_{A}[p_{A}]=\Pr_{A\in\mathbb{F}^{n\times n}}[A\in
X_{i,j}]\cdot 1+\Pr_{A\in\mathbb{F}^{n\times n}}[A\not\in
X_{i,j}]\cdot\delta/4\leq\Pr_{A\in\mathbb{F}^{n\times n}}[A\in
X_{i,j}]+\delta/4\enspace,$

and hence, $\Pr[A\in X_{i,j}]\geq\delta/4$.

The furthermore part is by definition of $X_{i,j}$. ∎

######  Definition 4.5.

Denote by $L^{k}\in\mathbb{F}^{n\times n}$ a random matrix of rank at most
$k$, constructed by sampling the first $k$ columns independently uniformly at
random from $\mathbb{F}^{n}$, and then taking the remaining $n-k$ columns to
be uniformly random linear combinations of the first $k$ vectors.

The following lemma is from [AGGS22]. It shows that if $L^{2k}_{A}$ is a
random matrix of rank at most $2k$ sampled as in Definition 4.5, then
$M_{A}=A-(L^{2k}_{A})$ belongs to any subspace of matrices of co-dimension $k$
with a non-negligible probability. We provide the proof of the lemma here for
completeness.

######  Lemma 4.6 (Lemma 4.8 from [AGGS22]).

Fix a matrix $A\in\mathbb{F}_{2}^{n\times n}$, let $k$ be a parameter, and let
$\ell\geq 2k$. Let $L^{\ell}_{A}$ be a random matrix of rank at most $\ell$
sampled as in Definition 4.5, and let $M_{A}=A-(L^{\ell}_{A})$. Then, for any
subspace $V\subseteq\mathbb{F}^{n\times n}$ of $\dim(V)\geq n^{2}-k$ it holds
that

$\Pr[M_{A}\in V]\geq\frac{1}{2{\left|\mathbb{F}\right|}^{k}}\enspace.$

###### Proof.

Since $V$ has co-dimension at most $k$, a matrix in $V$ must be orthogonal to
all the basis vectors of its orthogonal complement. Since there are up to $k$
such basis vectors, the membership condition of $M_{A}$ in $V$ can be written
down in the form of $k$ linear constraints. Viewing $M_{A}$ as a vector in
$\mathbb{F}^{n^{2}}$, we can write the $k$ linear constraints on the elements
of the matrix $M^{2k}_{A}$ in the form

$\alpha_{1}\cdot(M_{A})_{i_{1},j_{1}}+\alpha_{2}\cdot(M_{A})_{i_{2},j_{2}}+\cdots+\alpha_{r}(M_{A})_{i_{t},j_{t}}=0\enspace.$

Here, $\alpha_{i}$’s are constants and $t\in[n^{2}]$ is the number of elements
upon which the constraints depend. Writing $M_{A}$ as
$m\in\mathbb{F}^{n^{2}}$, we can represent these linear constraints as a
system of equations of the form $G\cdot m=0$, where $G$ is a $k\times n^{2}$
matrix. Now, we perform Gaussian elimination on $G$, which gives us a matrix
$G^{\prime}$, where each row has a 1 entry such that all the other entries in
the column containing this 1 are 0. We refer to such 1’s as _leading 1’s_.
That is, by permuting the columns of $G^{\prime}$, we may think of it as being
of the form $G^{\prime}=[I_{k}|G^{*}]$.

Consider the set of $k$ coordinates of $m$ corresponding to the $k$ leading
$1$s in $G^{\prime}$, one from each row. These $k$ coordinates of $m$ in turn
correspond to $k$ pairs of coordinates
$\\{(i_{1},j_{1}),(i_{2},j_{2})\ldots(i_{k},j_{k})\\}$ in $M_{A}$. These $k$
pairs of coordinates in $M_{A}$ can belong to at most $k$ rows in $M_{A}$. We
now bound the probability of none of these $k$ rows in $L_{A}^{\ell}$ being a
linear combination of the other rows. Let us denote this event as $\Omega$.
Then

$\displaystyle\Pr[\Omega]$
$\displaystyle=\left(1-\frac{1}{2^{\ell}}\right)\left(1-\frac{2}{2^{\ell}}\right)\left(1-\frac{4}{2^{\ell}}\right)\cdots\left(1-\frac{2^{k-1}}{2^{\ell}}\right)$
$\displaystyle\geq\left(1-\frac{2^{k-1}}{2^{\ell}}\right)^{k}$
$\displaystyle\geq\left(1-\frac{1}{2^{k+1}}\right)^{k}$ $\displaystyle\geq
1-\frac{k}{2^{k+1}}\geq\frac{1}{2}$

If $\Omega$ happens, then we get a coordinate $(i_{r},j_{r})$ in $M_{A}$
corresponding to the $r^{th}$ linear constraint, for all $r\in[k]$, such that
no other constraint depends upon it (as it corresponds to a leading $1$) and
it is chosen uniformly at random (since the rows containing these coordinates
are linearly independent). Therefore, the probability that this random value
satisfies the $i^{th}$ constraint is $1/\lvert\mathbb{F}\rvert$. To see this,
assume that the values of all other coordinates involved in the $i^{th}$
constraint are fixed, then we are left with only one choice for the value of
the coordinate $(c_{i},c_{i}^{\prime})$ which satisfies the constraint.
Therefore, we have

$\displaystyle\Pr[M_{A}\in V]$ $\displaystyle=\Pr[\text{All $k$ linear
constraints are satisfied}]$
$\displaystyle=\Pr[\Omega]\cdot\frac{1}{\mathbb{\lvert
F\rvert}^{k}}\geq\frac{1}{2\mathbb{\lvert F\rvert}^{k}}\enspace.$

This completes the proof of Lemma 4.6. ∎

###  4.1 Computing the good coordinates

Next, we start describing the reduction guaranteed by Theorem 1.3. As a first
step we design Algorithm 2, that gets two matrices $A,B$ and outputs a matrix
$C$ with values in $\mathbb{F}_{2}\cup\\{\bot\\}$, satisfying the following
guarantees.

  1. 1

If $C_{i,j}\neq\bot$, then $C_{i,j}$ contains the correct values, i.e.,
$C_{i^{*},j^{*}}=(A\cdot B)_{i^{*},j^{*}}$.

  2. 2

For any good coordinate $(i^{*},j^{*})\in G$ we have
$\Pr[C_{i^{*},j^{*}}=(A\cdot B)_{i^{*},j^{*}}]\geq\delta_{0}$, where
$\delta_{0}$ is some constant that depends only on $\delta$. That is, with
non-negligible probability $C_{i^{*},j^{*}}$ contains the correct answer, and
not $\bot$.

Then, in Section 4.2 we use Algorithm 2 as a subroutine, in order to compute
the entire matrix $A\cdot B$ correctly.

Input: $A,B\in\mathbb{F}_{2}^{n\times n}$, ${\mathrm{A}LG}$

Output: An $n\times n$ matrix $C$ with values $\mathbb{F}_{2}\cup\\{\bot\\}$

1

2Let $t>\log(4/\delta)+2$

3Let $k=O(\log(1/{\delta}))$ from the “in particular” part of Lemma 2.4 with
$\alpha=\delta/4$ and $t$ chosen above

4Sample two random matrices $L_{A}^{2k}$ and $L_{B}^{2tk}$ of rank at most
$2k$ and $2tk$ respectively, as in Definition 4.5

5Define $M_{A}=A-L_{A}^{2k}$ and $M_{B}=B-L_{B}^{2tk}$

6Let $C$ be the $n\times n$ matrix initialized with all $\bot$

7Sample $t-1$ random matrices $R_{1},R_{2},..,R_{t-1}\in\mathbb{F}^{n\times
n}_{2}$ and set $R_{t}=M_{A}-(R_{1}+R_{2}+..+R_{t-1})$

8for  _$r=1,\dots,t$_ do

9 Sample $t-1$ random matrices
$S^{(r)}_{1},S^{(r)}_{2},..,S^{(r)}_{t-1}\in\mathbb{F}^{n\times n}_{2}$ and
set $S^{(r)}_{t}=M_{B}-(S^{(r)}_{1}+S^{(r)}_{2}+..+S^{(r)}_{t-1})$

10

11 for  _$r,s=1,\dots,t$_ do

12 Compute ${\mathrm{A}LG}(R_{r},S^{(r)}_{s})$

13

14 for  _$(i,j)\in[n]\times[n]$_ do

15 if  _${\mathrm{A}LG}(R_{r},S^{(r)}_{s})_{i,j}=1$ for all
$r,s\in\\{1,\dots,t\\}$_ then

16 Set $C_{i,j}=\sum_{r,s}{\mathrm{A}LG}(R_{r},S^{(r)}_{s})_{i,j}\pmod{2}$

17

18

return $C+M_{A}\cdot L_{B}^{2tk}+L_{A}^{2k}\cdot M_{B}+L_{A}^{2k}\cdot
L_{B}^{2tk}$ // if $C_{i,j}=\bot$, then we return $\bot$ in the coordinate
$(i,j)$

Algorithm 2 Approximating good coordinates for one-sided error algorithms

In lines 3-4 we decompose $A=L_{A}^{2k}+M_{A}$, and $B=L_{B}^{2tk}+M_{B}$ with
the intention of computing $A\cdot B$ by writing

$\displaystyle AB$ $\displaystyle=(M_{A}+L_{A}^{2k})\cdot(M_{B}+L_{B}^{2tk})$
$\displaystyle=M_{A}\cdot M_{B}+M_{A}\cdot L_{B}^{2tk}+L_{A}^{2k}\cdot
M_{B}+L_{A}^{2k}\cdot L_{B}^{2tk}\enspace.$

Lines 5-13 try to compute $C=M_{A}\cdot M_{B}$. Then, in line 14, we sum up
the 4 terms. Using the fact that multiplication of matrices of rank $k$ takes
$O(kn^{2})$ time, the last three terms can be computed in $O(tkn^{2})$ time,
and hence, it remains to compute $M_{A}\cdot M_{B}$. The remainder of this
subsection is dedicated to analyzing lines 5-13, which contain the most
involved part of the algorithm.

We would like to compute $M_{A}\cdot M_{B}$ by writing
$M_{A}=R_{1}+R_{2}+\dots+R_{t}$, and
$M_{B}=S^{(r)}_{1}+S^{(r)}_{2}+\dots+S^{(r)}_{t}$ for $r=1\dots t$, and then
computing ${\mathrm{A}LG}(R_{r},S^{(r)}_{s})$ for all $r,s$. Note that if we
could guarantee that ${\mathrm{A}LG}(R_{r},S^{(r)}_{s})$ returns $R_{r}\cdot
S^{(r)}_{s}$, then, we would have

$M_{A}\cdot M_{B}=\sum_{r,s}R_{r}\cdot
S^{(r)}_{s}=\sum_{r,s}{\mathrm{A}LG}(R_{r},S^{(r)}_{s})\enspace.$ (1)

However, ${\mathrm{A}LG}$ is not guaranteed to return the product of the
inputs correctly. Instead, we claim that (1) for _some_ good coordinates
$(i^{*},j^{*})$ it holds $C_{i^{*},j^{*}}=(M_{A}\cdot M_{B})_{i^{*},j^{*}}$,
and (2) the remaining coordinates in $C$ remain $\bot$. This is summarized
formally in the next two claims.

######  Claim 4.7.

For any $(i,j)\in[n]\times[n]$ if $C_{i,j}\in\\{0,1\\}$ (i.e.,
$C_{i,j}\neq\bot$), then $C_{i,j}=(M_{A}\cdot M_{B})_{i,j}$.

###### Proof.

Fix any coordinate $(i,j)$. Note that in line 13 we set the value of
$C_{i,j}=\sum_{r,s}{\mathrm{A}LG}(R_{r},S^{(r)}_{s})_{i,j}\pmod{2}$ only if
${\mathrm{A}LG}(R_{r},S^{(r)}_{s})_{i,j}=1$ for all $r,s$. Recall, that by the
assumption of the algorithm if ${\mathrm{A}LG}(R_{r},S^{(r)}_{s})_{i,j}=1$,
then ${\mathrm{A}LG}(R_{r},S^{(r)}_{s})_{i,j}=(R_{r}\cdot S^{(r)}_{s})_{i,j}$.
The claim follows by Eq. 1 restricted to the coordinate $(i,j)$, as

$C_{i,j}=\sum_{r,s}{\mathrm{A}LG}(R_{r},S^{(r)}_{s})_{i,j}=\sum_{r,s}(R_{r}\cdot
S^{(r)}_{s})_{i,j}=(M_{A}\cdot M_{B})_{i,j}\enspace,$

as required. ∎

######  Claim 4.8.

Fix a good coordinate $(i^{*},j^{*})\in G$. Then
$\Pr[C_{i^{*},j^{*}}=(M_{A}\cdot
M_{B})_{i^{*},j^{*}}]\geq\delta_{0}=0.5^{O(\log^{3}(1/\delta))}$.

###### Proof.

Consider the set $X_{i^{*},j^{*}}$ from Definition 4.3 for a good entry
$(i^{*},j^{*})\in G$. By 4.4, the density of $X_{i,j}$ is at least
${\delta}/4$, and hence, Lemma 2.4 guarantees the existence of an affine
subspace $V_{i^{*},j^{*}}$ of dimension $\dim(V_{i^{*},j^{*}})\geq n-k$. Then,
using Lemma 4.6 with $V_{i^{*},j^{*}}$ we have

$\Pr[M_{A}\in{V_{X_{i,j}}}]\geq\frac{1}{2\cdot 2^{k}}\enspace.$ (2)

Let us condition on the event that $M^{A}\in{V_{X_{i^{*},j^{*}}}}$. Then by
Lemma 2.4,

$\Pr_{\begin{subarray}{c}R_{1},R_{2},..,R_{t}\in\mathbb{F}_{n}^{n\times n}\\\
\sum_{r}R_{r}=M_{A}\end{subarray}}[R_{1},R_{2},..,R_{t}\in{X_{i,j}}]\geq(\delta/8)^{t}\enspace.$
(3)

For each of $R_{1},\dots,R_{t}$ define the sets
$Y^{R_{1}},Y^{R_{2}},\dots,Y^{R_{t}}$ as in Definition 4.3. (Recall $Y^{R}$ is
the set of are all matrices $S$ such that $(R\cdot S)_{i^{*},j^{*}}=1$. We
omit the subscript $(i^{*},j^{*})$ for readability.)

From 4.4, we know each of $Y^{R_{1}},\dots,Y^{R_{t}}$ has density at least
$\delta/4$. Hence, by applying Lemma 2.4 on each of them, we obtain subspaces
$V_{Y^{R_{1}}},..,V_{Y^{R_{t}}}$ of co-dimension at most $k$. Define
$V_{Y}=V_{Y^{R_{1}}}\cap V_{Y^{R_{2}}}\cap\dots\cap V_{Y^{R_{t}}}$ to be their
intersection, and note that $\dim(V_{Y})\geq n-tk$. Therefore, by applying
Lemma 4.6 on the matrix $B$ with the subspace $V_{Y}$, we get 111Note that
although the algorithm samples $M_{B}$ before $R_{1},\dots,R_{t}$, in fact
they are sampled independently of each other, and hence Lemma 4.6 is
applicable here.

$\Pr[M_{B}\in V_{Y}]\geq\frac{1}{2\cdot{2}^{tk}}\enspace.$ (4)

Conditioning further on the event that $M_{B}\in V_{Y}$, we apply Lemma 2.4,
and for each $r=1,\dots,t$ we get

$\Pr_{\begin{subarray}{c}S^{(r)}_{1},\dots,S^{(r)}_{t}\\\
\sum_{s}S^{(r)}_{s}=M_{B}\end{subarray}}[S^{(r)}_{1},\dots,S^{(r)}_{t}\in{Y^{R_{r}}}]\geq(\delta/8)^{t}\enspace.$

Since the events above are independent between different $r$’s, the
probability that the algorithm returns correct output on the entry
$(i^{*},j^{*})$ is lower bounded by the product of the probabilities in Eqs.
2, 3 and 4, and hence

$\Pr[C_{i^{*},j^{*}}=(M_{A}\cdot
M_{B})_{i^{*},j^{*}}]\geq\frac{1}{2^{k+1}}\times(\delta/8)^{t}\times\frac{1}{2\cdot{2}^{tk}}\times\left((\delta/8)^{t}\right)^{t}\geq\frac{1}{2^{O(\log^{3}(1/\delta))}}\enspace.$

This completes the proof of 4.8.

∎

###  4.2 Proof of Theorem 1.3

We are now ready to prove Theorem 1.3. Algorithm 3 uses Algorithm 2 as a
subroutine, by running it several times. We claim that Algorithm 3 correctly
computes the correct answer with high probability for any input $A,B$. Note
that Algorithm 2 is guaranteed to be correct only for good coordinates,
although it does not get the good coordinates as an input, and the guarantee
about the good coordinates only appears in the analysis.

Input: $A,B\in\mathbb{F}_{2}^{n\times n}$

Output: $A\cdot B$

1

2Let $C$ be the $n\times n$ matrix initialized with all $\bot$.

3Let $\delta_{0}$ be the constant from 4.8

4repeat  _$O\left(\frac{\log(n)}{\delta\cdot\delta_{0}}\right)$_ times

5 Sample uniformly random $\pi,\sigma\in[n]$.

6 Run Algorithm 2 with the inputs as $A^{\pi,0},B^{0,\sigma},{\mathrm{A}LG}$.

7 Let $C^{*}$ be the resulting matrix

8 for  _$(i,j)\in[n]\times[n]$_ do

9 if  _$C^{*}_{i+\pi\pmod{n},j+\sigma\pmod{n}}\neq\bot$_ then

10 Set $C_{i,j}=C^{*}_{i+\pi,j+\sigma}$

11

12

return C

Algorithm 3 Approximation for one-sided Agreement Matrix Multiplication
Algorithms

The following claim completes the proof of Theorem 1.3.

######  Claim 4.9.

Fix a coordinate $(i,j)\in[n]\times[n]$. Algorithm 3 returns the matrix $C$
such that $\Pr[C_{i,j}=(A\cdot B)_{i,j}]\geq 1-1/n^{3}$.

In particular, by taking the union bound over all coordinates $(i,j)$ it
follows that for any input $A,B$ Algorithm 3 returns their product with
probability at least $1-1/n$.

###### Proof.

Fix a coordinate $(i,j)\in[n]\times[n]$. The algorithm chooses random $\pi$
and $\sigma$, and runs Algorithm 2 on the shifted matrices $A^{\pi,0}$ and
$B^{0,\sigma}$.

Note that since $\pi$ and $\sigma$ are chosen uniformly at random, it follows
that $\Pr[(i+\pi\pmod{n},j+\sigma\pmod{n})\in
G]={\left|G\right|}/n^{2}\geq\delta/2$.

Suppose that $(i+\pi\pmod{n},j+\sigma\pmod{n})$ is indeed a good coordinate.
Then by 4.9 with probability $\delta_{0}$ we obtain the correct answer in the
coordinate $(i+\pi\pmod{n},j+\sigma\pmod{n})$, in which case we set $C_{i,j}$
to be that answer $(A\cdot B)_{i,j}$. Otherwise, $C_{i,j}$ remains $\bot$.

Therefore, with probability at least $(\delta/2)\cdot\delta_{0}$ in each
iteration $C_{i,j}$ changes from $\bot$ to $(A\cdot B)_{i,j}$, and once it
changes, it never changes its value again.

By repeating the procedure $R=O(\frac{\log(n)}{\delta\cdot\delta_{0}})$ times,
the probability that in the end of the algorithm $C_{i,j}=\bot$ is upper
bounded by $\Pr[C_{i,j}=\bot]\leq(1-\delta_{0})^{R}<1/n^{3}$. This completes
the proof of the claim. ∎

We conclude the proof with the analysis of the running time of the algorithm.

#### Running Time:

The total running time of Algorithm 3 is essentially dominated by the running
time of Algorithm 2 multiplied by
$O\left(\frac{\log(n)}{\delta\cdot\delta_{0}}\right)$. Each iteration of
Algorithm 2 involves $O(t^{2})$ calls to ${\mathrm{A}LG}$ plus additional
$O(tn^{2})$ time. Therefore, the running time is
$O(t^{2}\log(n)T(n)/\delta\cdot\delta_{0})$. Since $t=O(\log 1/{\delta})$, and
$\delta_{0}=0.5^{O(\log^{3}(1/\delta))}$ the total running time of the
algorithm is is $2^{O(\log^{3}(1/\delta))}\cdot T(n)\log(n)$.

In particular, even for a slightly sub-constant
$\delta\geq\exp(-\log^{0.33}(n))$, our algorithm runs in time $T(n)\cdot
n^{o(1)}$.

This completes the proof of Theorem 1.3.

## References

  * [AGG+] Vahid R. Asadi, Alexander Golovnev, Tom Gur, Igor Shinkar, and Sathyawageeswar Subramanian.  Quantum Worst-Case to Average-Case Reductions for All Linear Problems, pages 2535–2567.
  * [AGGS22] Vahid R. Asadi, Alexander Golovnev, Tom Gur, and Igor Shinkar.  Worst-case to average-case reductions via additive combinatorics.  In STOC 2022, pages 00–00. ACM, 2022.
  * [AV21] Josh Alman and Virginia Vassilevska Williams.  A refined laser method and faster matrix multiplication.  In SODA 2021, pages 522–539. SIAM, 2021.
  * [BCRL79] Dario Bini, Milvio Capovani, Francesco Romani, and Grazia Lotti.  O(n2.7799) complexity for n × n approximate matrix multiplication.  Information Processing Letters, 8(5):234–235, 1979.
  * [BDCGL92] Shai Ben-David, Benny Chor, Oded Goldreich, and Michel Luby.  On the theory of average case complexity.  Journal of Computer and System Sciences, 44(2):193–219, 1992.
  * [BLR90] Manuel Blum, Michael Luby, and Ronitt Rubinfeld.  Self-testing/correcting with applications to numerical problems.  In STOC 1990, pages 73–83. ACM, 1990.
  * [BT06] Andrej Bogdanov and Luca Trevisan.  Average-case complexity.  Found. Trends Theor. Comput. Sci., 2(1):1–106, oct 2006.
  * [Cha02] Mei-Chu Chang.  A polynomial bound in Freiman’s theorem.  Duke Mathematical Journal, 113(3):399 – 419, 2002.
  * [CW90] Don Coppersmith and Shmuel Winograd.  Matrix multiplication via arithmetic progressions.  Journal of Symbolic Computation, 9(3):251–280, 1990.  Computational algebraic complexity editorial.
  * [DWZ23] Ran Duan, Hongxun Wu, and Renfei Zhou.  Faster matrix multiplication via asymmetric hashing.  In 64th IEEE Annual Symposium on Foundations of Computer Science, FOCS 2023, Santa Cruz, CA, USA, November 6-9, 2023, pages 2129–2138. IEEE, 2023.
  * [Gol11] Oded Goldreich, editor.  Studies in complexity and cryptography: miscellanea on the interplay between randomness and computation.  Springer-Verlag, Berlin, Heidelberg, 2011.
  * [HS23] Shuichi Hirahara and Nobutaka Shimizu.  Hardness self-amplification: Simplified, optimized, and unified.  In Proceedings of the 55th Annual ACM Symposium on Theory of Computing, STOC 2023, page 70–83. Association for Computing Machinery, 2023\.
  * [IL90] Russell Impagliazzo and Leonid A. Levin.  No better ways to generate hard NP instances than picking uniformly at random.  In 31st Annual Symposium on Foundations of Computer Science, St. Louis, Missouri, USA, October 22-24, 1990, Volume II, pages 812–821. IEEE Computer Society, 1990.
  * [Imp95] Russell Impagliazzo.  A personal view of average-case complexity.  Proceedings of Structure in Complexity Theory. Tenth Annual IEEE Conference, pages 134–147, 1995.
  * [Imp11] Russell Impagliazzo.  Relativized separations of worst-case and average-case complexities for NP.  In Proceedings of the 26th Annual IEEE Conference on Computational Complexity, CCC 2011, San Jose, California, USA, June 8-10, 2011, pages 104–114. IEEE Computer Society, 2011.
  * [Lev86] Leonid A. Levin.  Average case complete problems.  SIAM Journal on Computing, 15(1):285–286, 1986.
  * [LG14] François Le Gall.  Powers of tensors and fast matrix multiplication.  In Proceedings of the 39th International Symposium on Symbolic and Algebraic Computation, ISSAC ’14, page 296–303, New York, NY, USA, 2014\. Association for Computing Machinery.
  * [Pan78] V. Ya. Pan.  Strassen’s algorithm is not optimal trilinear technique of aggregating, uniting and canceling for constructing fast algorithms for matrix operations.  In 19th Annual Symposium on Foundations of Computer Science (SFCS 1978), pages 166–176, 1978.
  * [Rom82] Francesco Romani.  Some properties of disjoint sums of tensors related to matrix multiplication.  SIAM Journal on Computing, 11(2):263–267, 1982.
  * [Sch81] A. Schönhage.  Partial and total matrix multiplication.  SIAM Journal on Computing, 10(3):434–455, 1981.
  * [Sto10] Andrew James Stothers.  On the complexity of matrix multiplication.  2010\.
  * [Str69] Volker Strassen.  Gaussian elimination is not optimal.  Numerische Mathematik, 13(4):354–356, 1969.
  * [Str86] V. Strassen.  The asymptotic spectrum of tensors and the exponent of matrix multiplication.  In 27th Annual Symposium on Foundations of Computer Science (sfcs 1986), pages 49–54, 1986.
  * [Wil12] Virginia Vassilevska Williams.  Multiplying matrices faster than coppersmith-winograd.  In Proceedings of the Forty-Fourth Annual ACM Symposium on Theory of Computing, STOC ’12, page 887–898, New York, NY, USA, 2012. Association for Computing Machinery.


# An adaptive kernel estimator for the intensity function of spatio-temporal
point processes

Jonatan A<EMAIL_ADDRESS>[ Paula
<EMAIL_ADDRESS>[ Computer, Electrical and
Mathematical Sciences and Engineering Division, King Abdullah University of
Science and Technology (KAUST), Thuwal 23955-6900, Saudi Arabia presep=,
]e1,e2

###### Abstract

In spatio-temporal point pattern analysis, one of the main statistical
objectives is to estimate the first-order intensity function, i.e., the
expected number of points per unit area and unit time. This estimation is
usually carried out non-parametrically through kernel functions, where one of
the most frequent handicaps is the selection of kernel bandwidths prior to
estimation. This work presents an intensity estimation mechanism in which the
spatial and temporal bandwidths change at each data point in a spatio-temporal
point pattern. This class of estimators is called adaptive estimators, and
although there have been studied in spatial settings, little has been said
about them in the spatio-temporal context. We define the adaptive intensity
estimator in the spatio-temporal context and extend a partitioning technique
based on the bandwidths quantiles to perform a fast estimation. We demonstrate
through simulation that this technique works well in practice with the
partition estimator approximating the direct estimator and much faster
computation time. Finally, we apply our method to estimate the spatio-temporal
intensity of fires in the Amazonia basin.

Amazonia fires,

Bandwidth selection,

Intensity function,

Partitioning algorithm,

Spatio-temporal point process,

###### keywords:

, and

##  1 Introduction

In spatio-temporal point processes, one of the essential characteristics of a
given observation, that is, of a point pattern, is the first-order intensity,
which corresponds to the localised expected number of points in the
observation window and temporal interval (González et al., 2016). Knowing the
intensity of a point process is often extremely important for interpretation
since it represents the first moment, an analogy to the average of a
population (Baddeley, Rubak and Turner, 2015), and therefore it is a
fundamental descriptive characteristic. Additionally, this first-order
statistic is a key ingredient in the subsequent estimation of some higher-
order descriptors, such as those that measure the interaction between pairs of
points (González et al., 2016). Furthermore, the intensity of one process can
be related to the intensity of another process, obtaining a quantity commonly
known as relative risk, which is used to investigate the fluctuations of one
process compared to another (Davies, Marshall and Hazelton, 2018).

Kernel smoothing is a non-parametric technique classically used to estimate
some types of functions, such as intensities and probability density
functions. This technique is increasingly common for intensity estimation in
spatial and spatio-temporal statistics, especially in those cases where
additional information is unavailable to understand the distribution of points
in space or space-time (Fernando and Hazelton, 2014). One of the main
disadvantages of Kernel estimation is the need of prior knowledge of the
bandwidth or smoothing parameter of the kernel used for estimation. Larger
bandwidths provide smoother estimates and vice versa. Regardless of the data
dimensionality, this smoothing parameter is fundamental for the adequate
estimation of the intensity, and a wrong choice may have unfortunate
consequences (see e.g. Baddeley, Rubak and Turner, 2015, and references
therein). In the statistical literature, it is widespread to use fixed
bandwidth kernels, mainly due to the simplicity of the implementation and the
steps involved. However, this classical approach’s lack of spatial and
temporal adaptability often results in poor estimation, especially in highly
heterogeneous point processes with very complicated underlying features that
affect the intensity structure within the study region (Davies and Baddeley,
2018; Davies, Marshall and Hazelton, 2018). The fixed kernel density estimator
struggles to capture important finer details in crowded areas when much
smoothing is applied to control noise where data is sparse. Conversely, if
less smoothing is applied (to retain more of this detail), we can expect to
see spurious bumps generated by isolated points in low-density regions.

Therefore, a more intuitive approach emerges consisting of variable smoothing;
in this technique, the amount of smoothing is inversely related to the density
of the points. Known as adaptive, this smoothing has shown better bias levels
(Abramson, 1982; Hall and Marron, 1988) and lower integrated square error
(Diggle, Rowlingson and Su, 2005). There are two different ways to approach
variable bandwidth. The fitting can be performed at each location on the
function evaluation grid, i.e., the spatio-temporal grid in which the
intensity function is numerically evaluated. This technique is related to the
so-called locally adaptive or balloon estimators. Alternatively, the rescaled
bandwidth of the kernels can be associated with the data points; this
technique is called the point fitting, sample smoothing, or sample point
method (see e.g. Scott, 2015, and references therein). Following Davies and
Baddeley (2018), we focus on sample point estimators with the square root
methodology; this methodology assigns bandwidths proportional to the square
root of the pilot estimates (Abramson, 1982). To illustrate the idea of point
adaptation, Figure 1 provides an example of a spatio-temporal point pattern.
The ellipsoids correspond to a variable bandwidth in space and time for each
point.

Figure 1: (Right) Spatio-temporal point pattern with fixed kernel bandwidth in
space and time for each data point (drawn at one particular level as
translucent isosurfaces). (Left) Spatio-temporal point pattern with adaptive
kernel bandwidth spreading in space and time. In the latter image, the spatial
projection of the points is shown where an increase in time is interpreted as
moving up along the vertical axis.

We focus on the practical aspects of the methodology, combining theory with
implementation. We provide the code to apply the methods discussed here in the
newly developed software package kernstadapt (kernel adaptive estimation in
the spatio-temporal context). The package is freely available for the R
language on GitHub at https://github.com/jagm03/kernstadapt.

This article is organised as follows. Section 2 presents some basic concepts
related to spatio-temporal point processes, the first-order intensity
function, and its classical estimation through kernels. Section 3 explicitly
defines the expression for the space-time intensity estimator in its adaptive
version with two types of bandwidth, spatial and temporal; and the partition
algorithm presented for fast estimation. Section 4 shows how our method works
in terms of squared integrated errors and speed through simulation. In Section
5, we apply our method to a database of fire reports on Amazonia. Finally, we
conclude with specific comments and possible future lines of research in
Section 6.

##  2 Spatio-temporal point processes in a nutshell

We say that $X$ is a spatio-temporal point process whenever it is a random
countable subset of $\mathbb{R}^{2}\times\mathbb{R}$, for which
$|X\cap(A\times B)|<\infty$ for bounded $A\times
B\subseteq\mathbb{R}^{2}\times\mathbb{R}$; we use $|\cdot|$ to denote both set
cardinality and absolute value of a real number. On the other hand, a spatio-
temporal point pattern observed in $W\times T$ is a realisation of a spatio-
temporal point process in $\mathbb{R}^{2}\times\mathbb{R}$ to $W\times T$. We
denote $N(A\times B)$ as the number of points of a set $(A\times B)\cap X$,
where $A\subseteq W$ and $B\subseteq T$.

We say that a spatio-temporal point process $X$ is (spatio-temporally)
stationary when shifted processes $(\mathbf{u},v)+X$ have the same
distribution as the original process $X$ for any $(\mathbf{u},v)\in W\times
T$. $X$ is defined as (spatially) isotropic if, for any rotation ${\bf r}$
around the origin, every rotated point process ${\bf r}X=\\{({\bf
r}\mathbf{u},v):(\mathbf{u},v)\in X\\}$ has the same distribution as $X$ (see
e.g. Møller and Ghorbani, 2012; González et al., 2016).

###  2.1 Intensity

The expected value of the number of points governs the univariate
distributions of the points of $X$ in $W\times T$. A statistical analysis of
spatio-temporal point patterns usually begins estimating and modelling the
intensity function, i.e., the expected number of points per unit area per unit
time. When the intensity function $\lambda(\cdot)$ exists, the following
identity holds,

$\mathbb{E}[N(A\times B)]=\int_{A\times
B}\lambda(\mathbf{u},v)\mathrm{d}\mathbf{u}\mathrm{d}v.$

When $X$ is stationary, also referred to as $X$ being homogeneous (or
completely stationary, see Illian et al., 2008), then
$\lambda(\mathbf{u},v)\equiv\lambda>0$. This constant is referred to as the
intensity of $X$.

####  2.1.1 Separability

The separability of the intensity function of a spatio-temporal point process
is a pragmatic and common assumption in the spatio-temporal point process
literature (see, e.g., Gabriel and Diggle, 2009; Møller and Ghorbani, 2012;
Gabriel, Rowlingson and Diggle, 2013; González et al., 2016). When the
intensity function of a spatio-temporal temporal point process can be
factorised (almost everywhere) as

$\lambda(\mathbf{u},v)=\lambda_{1}(\mathbf{u})\lambda_{2}(v),$

where $\lambda_{1}(\cdot)$ and $\lambda_{2}(\cdot)$ are non-negative
functions, the process is referred to as first-order spatio-temporal
separable. Note that these functions are not unique. The separability
assumption simplifies the analysis and the estimation; however, it rarely
holds in practice, so it becomes essential to perform the estimation through
non-separable estimators. Several works have focused on statistical tests to
verify the separability assumption; see, e.g., Schoenberg (2004); Díaz-Avalos,
Juan and Mateu (2013); Fuentes-Santos, González-Manteiga and Mateu (2018), and
more recently, Ghorbani et al. (2021).

####  2.1.2 Kernel estimator with fixed bandwidths

For non-parametric estimation of the spatio-temporal intensity function, not
necessarily separable, it is common to use a kernel estimator (see, e.g.,
Fernando and Hazelton, 2014; González et al., 2016, and references therein),

$\hat{\lambda}\left(\mathbf{u},v\right)=\frac{1}{e\left(\mathbf{u},v\right)}\sum_{i=1}^{n}K^{\text{s}}_{\epsilon}\left(\mathbf{u}-\mathbf{u}_{i}\right)K^{\text{t}}_{\delta}\left(v-v_{i}\right),\qquad(\mathbf{u},v)\in
W\times T,$ (1)

where $K^{\text{s}}_{\epsilon}(\cdot)$ and $K^{\text{t}}_{\delta}(\cdot)$ are
bivariate and univariate kernels for space and time, $\epsilon,\delta>0$ are
the bandwidths, i.e., the smoothing parameters, and
$e\left(\mathbf{u},v\right)$ is an edge-correction factor included in the
estimation. The edge correction term intends to correct the bias that comes
from the estimation at the edges (see e.g. Diggle, 1985; Jones, 1993) or
preserve the estimator’s mass, i.e., to guarantee that $\int\hat{\lambda}=n$
(van Lieshout, 2011; Ghorbani, 2013). The uniform edge-correction factor can
be defined as

$e\left(\mathbf{u},v\right)=\int_{W}\int_{T}K^{\text{s}}_{\epsilon}\left(\mathbf{z}-\mathbf{u}\right)K^{\text{t}}_{\delta}\left(t-v\right)\mathrm{d}\mathbf{z}\mathrm{d}t.$
(2)

The bandwidth selection for the estimation of spatio-temporal (or any
dimension) kernels remains an open issue in statistics (Illian et al., 2008),
as there is no universal rule to set it; however, there are many methods to
choose a bandwidth in the spatial case. Some of those methods rely on
minimising error measures (Diggle, 1985; Berman and Diggle, 1989), in rules of
thumb (Scott, 2015; Baddeley, Rubak and Turner, 2015), on likelihood cross
validation criterion (Loader, 1999), or Campbell’s formula (Cronie and Van
Lieshout, 2018). In the specific spatio-temporal context, unfortunately,
bandwidths have not been so profoundly covered.

####  2.1.3 Computation by convolution and FFT

Let $\mathcal{N}$ a counting measure assigning 1 to every point of the point
pattern $X=\\{(\mathbf{u}_{i},v_{i})\\}_{i=1}^{n}$, and let
$K(\mathbf{u},v):=K^{\text{s}}_{\epsilon}\left(\mathbf{u}\right)K^{\text{t}}_{\delta}\left(v\right)$
for some fixed bandwidths $\epsilon$ and $\delta$ and $(\mathbf{u},v)\in
W\times T$, then the convolution of $K$ with the counting measure
$\mathcal{N}$ is given by

$\displaystyle(K\ast\mathcal{N})(\mathbf{u},v)$ $\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{2}\times\mathbb{R}}K(\mathbf{u}-\mathbf{u^{\prime}},v-v^{\prime})\mathrm{d}\mathcal{N}(\mathbf{u^{\prime}},v^{\prime})$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{n}K(\mathbf{u}-\mathbf{u}_{i},v-v_{i})$
$\displaystyle=$ $\displaystyle\tilde{\lambda}(\mathbf{u},v),$

where $\tilde{\lambda}(\mathbf{u},v)$ is the numerator in Eq.(1). For the edge
correction, consider the indicator function $\mathbf{1}_{W\times
T}(\mathbf{u},v)$ that assigns $1$ to every point of the region $W\times T$
and $0$ elsewhere, then we consider the convolution of this function with $K$

$\displaystyle(\mathbf{1}_{W\times T}\ast K)(\mathbf{u},v)$ $\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{2}\times\mathbb{R}}\mathbf{1}_{W\times
T}(\mathbf{u^{\prime}},v^{\prime})K(\mathbf{u}-\mathbf{u^{\prime}},v-v^{\prime})\mathrm{d}(\mathbf{u^{\prime}},v^{\prime})$
$\displaystyle=$ $\displaystyle\int_{W\times
T}K(\mathbf{u}-\mathbf{u^{\prime}},v-v^{\prime})\mathrm{d}(\mathbf{u^{\prime}},v^{\prime})$
$\displaystyle=$ $\displaystyle e(\mathbf{u},v),$

this idea goes back to Koch, Ohser and Schladitz (2003) and it was considered
in the planar case by Davies and Baddeley (2018).

A well-known numerical method used to compute the intensity fast and
efficiently is first to discretise the data to a very fine grid, and then use
the Fast Fourier transform to apply the convolution over the data with the
kernel to obtain the estimate. This method is known as the binned estimator,
and Silverman (1982) first proposed it in the univariate context. Let
$\mathcal{R}\supseteq W\times T$ be a three-dimensional rectangle enclosing
$W\times T$; if we partition $\mathcal{R}$ into $M_{1}\times M_{2}\times
M_{3}$ rectangular voxels with equal volume $V=m_{1}m_{2}m_{3}$; let
$\mathbf{c}_{ijk}$ be the centroid of the voxel $ijk$ ($i$th horizontal, $j$th
vertical and $k$th temporal axes position). Let $I=\\{1,\ldots,M_{1}\\}$,
$J=\\{1,\ldots,M_{2}\\}$ and $L=\\{1,\ldots,M_{3}\\}$ be the index sets; we
define the grid $C=\\{\mathbf{c}_{ijk}:i\in I,j\in J,k\in L\\}$; there are
several alternatives to discretise the counting measure $\mathcal{N}$ measure
as a collection of weights $\omega_{ijk}$. The simple binning involves
assigning a unit mass to the grid point closest to every data point, i.e.,
$\omega_{ijk}$ is the number of observations of $X$ falling in the voxel
$ijk$. On the other hand, the linear binning (see, e.g., Jones and Lotwick,
1983; Wand, 1994) assigns a fractional mass to the grid point depending on how
far is the data point. In what follows, we use simple binning for simplicity.
Once we have the weights calculated, an approximate form of the numerator of
Eq. (1) is

$\tilde{\lambda}^{*}\left(\mathbf{c}_{ijk}\right)=\sum_{\ell_{3}=1}^{M_{3}}\sum_{\ell_{2}=1}^{M_{2}}\sum_{\ell_{1}=1}^{M_{1}}K\left(\mathbf{c}_{ijk}-\mathbf{c}_{\ell_{1}\ell_{2}\ell_{3}}\right)\omega_{\ell_{1}\ell_{2}\ell_{3}},$

Note that $\omega_{\ell_{1}\ell_{2}\ell_{3}}=0$ outside of the set $I\times
J\times L$, so after a straightforward change of variable we have that

$\tilde{\lambda}^{*}\left(\mathbf{c}_{ijk}\right)=\sum_{\ell_{3}=1-M3}^{M_{3}-1}\sum_{\ell_{2}=1-M_{2}}^{M_{2}-1}\sum_{\ell_{1}=1-M1}^{M_{1}-1}\omega_{i-\ell_{1},j-\ell_{2},k-\ell_{3}}K\left(m_{1}\ell_{1},m_{2}\ell_{2},m_{3}\ell_{3}\right),$
(3)

and for the edge correction factor we have that

$e^{*}\left(\mathbf{c}_{ijk}\right)=V\sum_{\ell_{3}=1-M3}^{M_{3}-1}\sum_{\ell_{2}=1-M_{2}}^{M_{2}-1}\sum_{\ell_{1}=1-M1}^{M_{1}-1}\mathbf{1}_{W\times
T}\left(\mathbf{c}_{i-\ell_{1},j-\ell_{2},k-\ell_{3}}\right)K\left(m_{1}\ell_{1},m_{2}\ell_{2},m_{3}\ell_{3}\right).$

A noteworthy point is that the array
$\tilde{\lambda}^{*}\left(\mathbf{c}_{ijk}\right)$ and
$e^{*}\left(\mathbf{c}_{ijk}\right)$ correspond to the discrete convolution of
the $\omega$’s and the indicators with the kernel so that they can be computed
rapidly using the fast Fourier transform (FFT). This algorithm only requires
$O(M_{1}M_{2}M_{3}\log(M_{1})\log(M_{2})\log(M_{3}))$ operations outperforming
the direct estimation that requires $O(M_{1}^{2}M_{2}^{2}M_{3}^{2})$
operations (Wand, 1994).

Let $\mathscr{F}$ denote the discrete Fourier transform, $\mathscr{F}^{-1}$
its inverse and $f$ and $g$ two integrable functions, then, as a consequence
of the convolution theorem $f\ast
g=\mathscr{F}^{-1}(\mathscr{F}(f)\mathscr{F}(g))$. However, the theorem
requires that the involved functions be periodic; we employ a zero-padding
technique to ensure the periodicity. We consider larger arrays of sizes
$P_{M_{1}}\times P_{M_{3}}\times P_{M_{3}}$ where $P_{M_{i}}=2M_{i},i=1,2,3.$
The weights array $\omega_{\mathbf{0}}$ is an array of dimensions
$P_{M_{1}}\times P_{M_{3}}\times P_{M_{3}}$ filled with zeros, except for the
sub-array $\omega=(\omega_{ijk})_{i\in I,j\in J,k\in L}$ located in the lowest
vertex. In the same way, we define a zero-padded array $\mathbf{1}_{W\times
T,\mathbf{0}}$ for the edge correction; the schemes of these arrays are
displayed in Fig. 2.

$\omega_{\mathbf{0}}:=$ | ${\mathbf{0}}$${\mathbf{0}}$${\mathbf{0}}$${\mathbf{0}}$${\omega}$${\mathbf{0}}$${\mathbf{0}}$${\mathbf{0}}$ | $\mathbf{1}_{W\times T,\mathbf{0}}:=$ | ${\mathbf{0}}$${\mathbf{0}}$${\mathbf{0}}$${\mathbf{0}}$${\mathbf{1}_{W\times T}}$${\mathbf{0}}$${\mathbf{0}}$${\mathbf{0}}$
---|---|---|---
Figure 2: Visual representation of the zero-padded weight matrix and indicator
function, $\mathbf{0}$ represents a generic matrix of zeros.

We denote the reflected kernel evaluation array by $K_{\mathbf{0}}$; in this
array, the kernel values are padded in a wrap-around order, and its dimensions
are $P_{M_{1}}\times P_{M_{2}}\times P_{M_{3}}.$ Note that the dimensional
index for the horizontal axis ranges from
$0,\ldots,m_{1}(M_{1}-1),-m_{1}M_{1},\ldots-m_{1}$, the vertical axis from
$0,\ldots,m_{2}(M_{2}-1),-m_{2}M_{2},\ldots-m_{2}$ and the depth axis from
$0,\ldots,m_{3}(M_{3}-1),-m_{3}M_{3},\ldots-m_{3}$. Fig. 3 shows the padding
scheme at two of the 2D layers of $K_{\mathbf{0}}$ at values $0$ and $-m_{3}$
(the first and the last).

$K_{\mathbf{0}}:=\quad$ |
---|---
Figure 3:  Two layers of the reflected kernel evaluation array at depth
coordinates $0$ and $-m_{3}$.

Finally, we gather all the ingredients to estimate the intensity by using the
extended arrays Eq. (3) to obtain

$\tilde{\lambda}^{*}\left(\mathbf{c}_{ijk}\right)=\Re\left[\mathscr{F}^{-1}(\mathscr{F}(\omega_{\mathbf{0}})\mathscr{F}(K_{\mathbf{0}}))\right]_{ijk};\quad
i\in I,j\in J,k\in L,$

and for the edge correction

$e^{*}\left(\mathbf{c}_{ijk}\right)=\Re\left[\mathscr{F}^{-1}(\mathscr{F}(\mathbf{1}_{W\times
T,\mathbf{0}})\mathscr{F}(K_{\mathbf{0}}))\right]_{ijk};\quad i\in I,j\in
J,k\in L,$

where $\Re[\cdot]$ represents the real component of a complex number required
due to the numerical error (Davies and Baddeley, 2018). By definition
$\tilde{\lambda}^{*}\left(\mathbf{c}_{ijk}\right):=0,\forall\mathbf{c}_{ijk}\notin
W\times T$ and the final spatio-temporal fixed estimator is given by

$\hat{\lambda}\left(\mathbf{c}_{ijk}\right)=\frac{\tilde{\lambda}^{*}\left(\mathbf{c}_{ijk}\right)}{e^{*}\left(\mathbf{c}_{ijk}\right)}.$
(4)

##  3 Spatio-temporal adaptive estimator

Although the adaptive estimator using kernels is well known in the
multivariate context; i.e., when the coordinates to be smoothed have several
dimensions, it has not been explicitly defined or studied in the spatio-
temporal point processes framework. The Gaussian kernels make both the three-
dimensional and the spatio-temporal estimators algebraically equivalent,
although with different bandwidths (Wand and Jones, 1995). This case is
usually ignored in the multivariate context or considered only in anisotropic
settings. We define the spatio-temporal adaptive kernel estimator as follows,

$\hat{\lambda}_{\epsilon,\delta}\left(\mathbf{u},v\right)=\frac{1}{e_{\epsilon,\delta}\left(\mathbf{u},v\right)}\sum_{i=1}^{n}{K^{\text{s}}_{\epsilon(\mathbf{u}_{i})}\left(\mathbf{u}-\mathbf{u}_{i}\right)K^{\text{t}}_{\delta(v_{i})}\left(v-v_{i}\right)},\qquad(\mathbf{u},v)\in
W\times T,$ (5)

where

$e_{\epsilon,\delta}\left(\mathbf{u},v\right)=\int_{W}\int_{T}K^{\text{s}}_{\epsilon(\mathbf{u^{\prime}})}\left(\mathbf{u}-\mathbf{u^{\prime}}\right)K^{\text{t}}_{\delta(v^{\prime})}\left(v-v^{\prime}\right)\mathrm{d}\mathbf{u^{\prime}}\mathrm{d}v^{\prime},$
(6)

is the spatio-temporal edge correction inspired in the initially developed by
Marshall and Hazelton (2010). This uniform edge correction attempts to fix the
bias of the estimate at the edges of the study region (spatial and temporal),
which can produce potentially unreliable estimations. Other edge correction
alternatives, such as Diggle’s edge correction for space-time (in which mass
is preserved) (see e.g. González, Hahn and Mateu, 2019; Ghorbani et al.,
2021), are compatible with our estimator and the partitioning technique
presented in Section 3.1 to speed up the computation. The bandwidth functions
are defined as

$\epsilon(\mathbf{u})=\frac{\epsilon^{\star}}{\gamma^{\text{s}}}\sqrt{\frac{n}{\lambda^{\text{s}}(\mathbf{u})}},\quad\text{and}\quad\delta(v)=\frac{\delta^{\star}}{\gamma^{\text{t}}}\sqrt{\frac{n}{\lambda^{\text{t}}(v)}},$
(7)

where $\epsilon^{\star}$ and $\delta^{\star}$ stand for spatial and temporal
smoothing multipliers known as global bandwidths,
$\lambda^{\text{s}}(\mathbf{u})$ and $\lambda^{\text{t}}(v)$ are marginal
intensity functions in space and time, and $\gamma^{\text{s}}$ and
$\gamma^{\text{t}}$ are the geometric mean terms for the marginal intensities
evaluated in the points of the point pattern. This approach was proposed
originally by Abramson (1982) and requires pilot (fixed-bandwidth) estimates
of the marginal intensities by using the global bandwidths; the inclusion of
the geometric mean terms aims to free the bandwidth factor from the data
scale, as pointed out by Silverman (1986) and Davies and Hazelton (2010).

###  3.1 Fast estimation through a partitioning algorithm

To quickly calculate the adaptive estimator, our idea is to reduce the number
of bandwidths, which in the direct estimator is $n$, to a much lower number,
and therefore reduce the number of operations in that same proportion. For
this purpose, we go to group the bandwidths in bins. Approximating an adaptive
estimate through binning dates back to Sain and Scott (1996) in the univariate
case and to Sain (2002) in the multivariate case. They used piecewise
functions for the bandwidth, and they showed that their estimators are much
better from a practical point of view as they offer good approximations and
significantly reduce the number of operations while maintaining the mean
integrated square error (MISE).

To calculate the adaptive estimator given in Eq.(5), we follow the methodology
proposed by Davies and Baddeley (2018). This is based on the discretisation of
the bandwidths chosen for each point through the empirical quartiles of its
sampling distribution. Starting from Eq. (7), we obtain a set of spatial
bandwidths $\\{\epsilon_{1},\ldots,\epsilon_{n}\\}$ and a set of temporal ones
$\\{\delta_{1},\ldots,\delta_{n}\\}$ by using pilot intensities calculated in
the classical context. We consider the empirical $p$th quantiles,
$\hat{\epsilon}^{(p)}$ and $\hat{\delta}^{(p)}$ of the bandwidths together
with two quantile steps, $\xi_{1},\xi_{2}\in(0,1],$ such that
$C_{1}=\xi_{1}^{-1}$ and $C_{2}=\xi_{2}^{-1}$ are integers. We define the
bandwidth bins through the values
$\\{\hat{\epsilon}^{(0)},\hat{\epsilon}^{(\xi_{1})},\hat{\epsilon}^{(2\xi_{1})},\ldots,\hat{\epsilon}^{(1)}\\}$
and
$\\{\hat{\delta}^{(0)},\hat{\delta}^{(\xi_{2})},\hat{\delta}^{(2\xi_{2})},\ldots,\hat{\delta}^{(1)}\\}$,
and we place each observation $(\mathbf{u}_{i},v_{i})$ in one of the bins
shown in Table 1.

$\left(\begin{tabular}[]{ c c c c
}$\left[\hat{\epsilon}^{(0)},\hat{\epsilon}^{(\xi_{1})}\right]\times\left[\hat{\delta}^{(0)},\hat{\delta}^{(\xi_{2})}\right]$&$\left(\hat{\epsilon}^{(\xi_{1})},\hat{\epsilon}^{(2\xi_{1})}\right]\times\left[\hat{\delta}^{(0)},\hat{\delta}^{(\xi_{2})}\right]$&$\cdots$&$\left(\hat{\epsilon}^{(\\{C_{1}-1\\}\xi_{1})},\hat{\epsilon}^{(1)}\right]\times\left[\hat{\delta}^{(0)},\hat{\delta}^{(\xi_{2})}\right]$\\\\[8.53581pt]
$\left[\hat{\epsilon}^{(0)},\hat{\epsilon}^{(\xi_{1})}\right]\times\left(\hat{\delta}^{(\xi_{2})},\hat{\delta}^{(2\xi_{2})}\right]$&$\left(\hat{\epsilon}^{(\xi_{1})},\hat{\epsilon}^{(2\xi_{1})}\right]\times\left(\hat{\delta}^{(\xi_{2})},\hat{\delta}^{(2\xi_{2})}\right]$&$\cdots$&$\left(\hat{\epsilon}^{(\\{C_{1}-1\\}\xi_{1})},\hat{\epsilon}^{(1)}\right]\times\left(\hat{\delta}^{(\xi_{2})},\hat{\delta}^{(2\xi_{2})}\right]$\\\
$\vdots$&$\vdots$&$\ddots$&$\vdots$\\\
$\left[\hat{\epsilon}^{(0)},\hat{\epsilon}^{(\xi_{1})}\right]\times\left(\hat{\delta}^{(\\{C_{2}-1\\})\xi_{2}},\hat{\delta}^{(1)}\right]$&$\left(\hat{\epsilon}^{(\xi_{1})},\hat{\epsilon}^{(2\xi_{1})}\right]\times\left(\hat{\delta}^{(\\{C_{2}-1\\})\xi_{2}},\hat{\delta}^{(1)}\right]$&$\cdots$&$\left(\hat{\epsilon}^{(\\{C_{1}-1\\}\xi_{1})},\hat{\epsilon}^{(1)}\right]\times\left(\hat{\delta}^{(\\{C_{2}-1\\})\xi_{2}},\hat{\delta}^{(1)}\right]$\end{tabular}\right)$

Table 1:  Array of spatio-temporal discretised bandwidth bins.

The sets of bins generate a disjoint partition of the original point pattern
$X$ into $C_{1}\times C_{2}$ sets $Y_{ij}$ and

$X=\bigcup_{ij}Y_{ij}.$

Suppose that in each subset $Y_{ij},$ the intensity is estimated using a
bandwidth defined as the midpoint of the respective bin where it belongs. In
that case, the intensity can be approximated as

$\hat{\lambda}_{\epsilon,\delta}\left(\mathbf{u},v\right)\approx\sum_{i=1}^{C_{1}}\sum_{j=1}^{C_{2}}\hat{\lambda}^{*}_{\bar{\epsilon}_{i},\bar{\delta}_{j}}\left(\mathbf{u},v|Y_{ij}\right),$
(8)

where $\bar{\epsilon}_{i}$ and $\bar{\delta}_{j}$ represent the midpoints of
the $i$th spatial and $j$th temporal bins and
$\hat{\lambda}^{*}_{\bar{\epsilon}_{i},\bar{\delta}_{j}}\left(\mathbf{u},v|Y_{ij}\right)$
is a fixed-bandwidth estimate based on the sub-pattern $Y_{ij}$. Each of the
addends on the right-hand side of Eq. (8) is calculated by the procedure
described in Section 2.1.3, that is, by Fourier transforms and Eq. (4). Notice
that in each addend we include the corresponding edge correction factor.

##  4 Simulation study

In this section, we perform a simulation study to assess the performance of
the partition method described in Section 3 compared with the direct method,
which consists of evaluating the kernel at each data point and each point of
the grid. For a single simulation, we generate a spatio-temporal point pattern
as a realisation of a spatio-temporal log-Gaussian Cox process. Then we
estimate the spatial and temporal global bandwidths and compute the pilot
marginal intensities using the projection of point pattern coordinates to
space and time, respectively. Finally, we compute the adaptive intensity
through the direct and partitioning algorithms, the latter for different
binning schemes. Along with each estimation, we calculate the wall timing
clock.

We start generating the point patterns, which are observations of point
processes with stochastic intensities. Let $\\{Z(\mathbf{u},v)\\}$ be a real-
valued Gaussian random field with mean $\mu(\mathbf{u},v)$ and covariance
$\operatorname{Cov}(Z(\mathbf{u}_{i},v_{i}),Z(\mathbf{u}_{j},v_{j}))$, and let
$\Lambda(\mathbf{u},v)=\exp{(Z(\mathbf{u},v))}$ be the intensity. Then, a log-
Gaussian Cox process is a Poisson point process in which the intensity
function is a realisation of $\Lambda(\mathbf{u},v)$, say
$\lambda(\mathbf{u},v)$ (see e.g., Møller, Syversveen and Waagepetersen, 1998;
Gabriel, Rowlingson and Diggle, 2013; Diggle, 2013; Diggle et al., 2013;
González et al., 2016). In this case, we employ an non-separable stationary
covariance function of the form (Gneiting, Genton and Guttorp, 2006),

$\operatorname{Cov}(\mathbf{u},v)=\frac{1}{\psi^{2}(v^{2})}\phi\left(\frac{||\mathbf{u}||^{2}}{\psi(v^{2})}\right),\quad\mathbf{u}\in\mathbb{R}^{2},v\in\mathbb{R},$
(9)

where $\phi(r),r>0$, is a strictly monotone function, and $\psi(v),v>0$ is a
positive function with a completely monotone derivative. In this case we use
the so-called stable model, i.e., $\phi(r)=\sigma^{2}\exp(-cr)$ where $c$ is a
non-negative scale parameter for space and $\sigma^{2}$ is the spatio-temporal
variance, and

$\psi(v)=\left(1+av^{\alpha}\right)^{\beta},\quad\alpha,\beta\in(0,1],$

were $a$ is non-negative scale parameter for time, $\alpha$ is an smoothness
parameter and $\beta$ is the space-time interaction parameter. We fix three
sets of parameters for the simulation that encompass some first-order non-
separable point processes with different clustering degrees;

  1. Set 1.

$\sigma^{2}=2,c=0.5,a=0.8$ and $\alpha=\beta=0.1$

  2. Set 2.

$\sigma^{2}=0.5,c=0.5,a=0.8$ and $\alpha=\beta=0.1$

  3. Set 3.

$\sigma^{2}=0.5,c=0.5,a=0.8$ and $\alpha=\beta=0.8$

Figure 4 displays three realisations of spatio-temporal log-Gaussian Cox
processes coming from Set 1, 2 and 3, with their underlined intensity
functions, where the numbers of points come from a Poisson distribution with a
mean of 1000.

Figure 4: Spatial snapshots of three realisation of a spatio-temporal log-
Gaussian point pattern with their underlying intensity function having non-
separable Gneiting covariance functions. The irregular spatial window is
enclosed in the unit square.

We generate 1000 realisations of log-Gaussian Cox processes with roughly 1000
points in a region $W\times T$, where $W$ is an irregular polygon contained in
the unit square, and $T$ is the unit interval; the covariance of the Gaussian
random field is Gneiting type (see, Eq. (9)).

The next step in the simulation process is estimating the global spatial
bandwidths; we use the maximal smoothing principle (Terrell, 1990; Davies,
Marshall and Hazelton, 2018) that consists of upper bounding the optimal
bandwidths with respect to the asymptotic mean integrate squared error (MISE)
of the density estimates to obtain oversmoothing bandwidths. These bandwidths
are given by $\epsilon^{\star}\approx 1.085\hat{\sigma}n^{-1/6},$ where
$\hat{\sigma}$ is the minimum half sum between the axis-specific standard
deviations and the interquartile ranges normalised by 1.34 (Silverman, 1986).
For the global temporal bandwidths, we use the method of Sheather and Jones
(1991), which relies on the estimation of the density derivatives.

Now, for computing the adaptive smoothing bandwidth functions (Eq. (7)), we
employ two separated pilot estimates, one for space and one for time computed
using convolutions and FFT in two and one dimensions, respectively.

Finally, for every spatio-temporal point pattern, we estimate the intensity by
using the partitioning algorithm by combining the number of bins in space and
time, i.e., using all combinations of
$\delta,\beta\in\\{0.1,0.05,0.025,0.01\\}$. We use a spatio-temporal
resolution of $128\times 128\times 64$; the spatial resolution is a default in
the spatstat and sparr packages (Baddeley, Rubak and Turner, 2015; Davies,
Marshall and Hazelton, 2018), and the temporal resolution was chosen as a
compromise between the spatial axes.

As we wish to compare the performance of our estimator (by assigning different
amounts of spatio-temporal bins) we opt for the $L_{2}$-norm, which in the
densities context is known as integrated squared error (ISE), and it is
defined as (Scott, 2015)

$\text{ISE}[\hat{f}]:=\int_{W}\int_{T}\left(\hat{f}(\mathbf{u},v)-f(\mathbf{u},v)\right)^{2}\mathrm{d}\mathbf{u}\mathrm{d}v,$

where $\hat{f}$ is the estimate using the partitioning algorithm and $f$ is
the direct estimate, both of them expressed in on density scale, i.e.,
normalising so that they integrate one.

Figure 5 shows the ISE boxplots for each collection of estimates across
bandwidth partitions. As expected and in line with Davies and Baddeley (2018),
coarser space-time partitions lead to larger errors. However, it is worth
mentioning that even the coarsest partition, i.e., ten spatial bins and ten
temporal bins $(\delta=0.1,\beta=0.1)$, generates an excellent estimate in
terms of error. As expected, the estimates get closer to the direct estimate
as the partitions get finer.

Figure 5: Integrated square errors with respect to direct estimation of
intensity function using the partitioning algortithm. We consider combinations
of the spatial $(\delta)$ and temporal $(\beta)$ bins ranging from 10 to 100
bins. The mean number of points in each simulated point pattern is $\mu=1000$
and the resolution of the intensity arrays is $128\times 128\times 64$. We
truncate the most extremes outliers for better understanding.

We also observe how elapsed execution times and bins increase together (Figure
6). However, in general terms, these times are always minimal compared to the
time taken by direct estimation, implying that the binning procedure is always
a considerably better solution in terms of execution times. It should be noted
that when the number of spatial and temporal bins is maximum
$(\delta=0.01,\beta=0.01)$, the time is considerably longer, and the gains in
terms of error do not seem to be justified.

Figure 6: Wall clock timings relative to direct estimation using the
partitioning algortithm.

##  5 Application: Amazonia fires

All around the globe, fires can affect all biomes significantly producing
severe alterations to the primary ecosystems (Andela et al., 2019). Fire
activity in Amazonia varies considerably from month to month, likely
influenced by human activity and climate change. The deployment of preventive
and corrective measures can be aided by adequately describing the spatio-
temporal evolution of the number of fires in every region of the Amazonia
biome. Here we estimate the intensity function of the point pattern of
wildfire ignition points observed in the Amazonian biome. The area covers
almost the complete Amazon basin and contains the Amazon rainforest, a part of
the tropical rainforest, and some adjacent areas to the north and east.

The Amazon Dashboard (http://www.globalfiredata.org) tracks individual fires
in the Amazon region using an cluster approach and classify active fire
detections by fire type. Known as VIIRS, The Visible Infrared Imaging
Radiometer Suite 375 m thermal anomalies provides data from a sensor aboard
the joint NASA/NOAA Suomi National Polar-orbiting Partnership (Suomi NPP) and
NOAA-20 satellites. The data are updated on a daily basis and there are
several characteristics recorded with each single fire. For this analysis we
select new and active deforestation fires (starting within past 24 hours) and
have a sample of $59,910$ fires from 01/01/2021 to 10/10/2021 (284 days).

Figure 7, top panel, displays the locations of fires in the Amazon. Although
it seems that there are few points in the map, this is a visual delusion due
to two reasons. The first is that the Amazon area is so large that it is
difficult to appreciate what happens on such a small scale; the second is the
high concentrations of points in minimal areas. Many fires seem to occur near
each other, resulting in many places showing no fire activity at all while all
of this activity is concentrated in a few locations. The pattern of temporal
counts (Figure 7, top panel) allows us to appreciate that temporally, the
pattern also lacks uniformity, being the month of August particularly prone to
fires with much more pronounced variations than in any other month. This high
inhomogeneity in the fire pattern makes estimating the intensity very
difficult if a fixed bandwidth was chosen, which is why we resorted to the
variable bandwidth kernel estimator discussed in this work to carry out the
estimation.

Figure 7: Deforestation fires in the Amazonia biome from 01/01/2021 to
10/10/2021, wider lighter points corresponds to older fires, and the
background corresponds to the elevation of the Amazonia biome. The bottom
panel shows the time series of fire counts throughout the sampling days.

On the other hand, given that the initial coordinates of both the fires and
the map are given in longitude-latitude, we have decided to use the Albers
equivalent conic projection, which preserves the areas, to estimate the
intensity.

We first used the partitioning algorithm described in Section 3.1. We chose 38
groups for the spatial bandwidths $(n^{1/3})$ and 6 for the temporal ones
$(n^{1/6})$; we select the numbers of groups this way to ensure that the total
number of groups is not bigger than $\sqrt{n}\approx 245$ meaning that so that
the computation time is only about 228 times slower than fixed-bandwidth
smoothing but roughly 59682 $(n-\sqrt{n})$ times faster that the direct
estimation. The estimation took less than three minutes on a conventional
computer.

A sample with some snapshots of the estimated spatio-temporal intensity can be
seen in Figure 8. The two panels of the figure have different scales since
they correspond to periods with hugely different densities in magnitude; if we
chose a single scale, the drawings would probably flatten out and be visually
indistinguishable.

Figure 8: Snapshots of the intensity fires (in a log scale) in Amazon region
upon eight specific times during the study period during . This estimation was
carried out through the partitioning algorithm described in Section 3.1. The
upper row corresponds to the first semester of the year and the lower one to
the second.

As expected, the local changes in the number of fires per area from one moment
to another show that the intensity is not separable. We can also observe that,
although the majority of the fires are concentrated in the south-central part,
that is, part of Brazil, Bolivia, Paraguay and Peru; the south-eastern part
also has large concentrations of fires; the northwestern part is the least
prone, and this low number extends throughout the north of the biome. The
fires are exacerbated in number during the warmer months, reaching spots with
entire kilometres burning.

##  6 Discussion

In this paper, we have presented an adaptive estimator for the first-order
intensity of a spatio-temporal point process; this estimator is non-parametric
and is based on only one realisation of the point process to be computed.
Given the nature of the estimator, in which every point is equipped with a
particular bandwidth, a possible direct estimation is computationally
unfeasible, and we have proposed an algorithm to perform this estimation in a
fast way. This algorithm is an extension of the one presented by Davies and
Baddeley (2018) for the spatial case, and it bases on the discretisation of
the candidate bandwidths for each point in a point pattern. We have proven via
simulation that our algorithm performs well in integrated squared error and
computation time. Finally, we have applied our method to estimate the
intensity of deforestation fires in the Amazon.

An alternative computational approach could be reformulating the adaptive
smoothing for a point pattern in $\mathbb{R}^{2}\times\mathbb{R}$ into one
convolution in $\mathbb{R}^{5}$, and calculating it using the FFT approach.
Davies and Baddeley (2018) conceived the adaptive estimation problem in the
planar case $(\mathbb{R}^{2})$ to a convolution problem in three-dimensional
scale-space $(\mathbb{R}^{3})$ (Chaudhuri and Marron, 2000). This method
dramatically reduced computational costs when performing the adaptive
estimation. Unfortunately, when we talk about the spatio-temporal case, things
work well in theory but not so well in practice. The key idea would be the
following. We use the scale-space concept and define a five-dimensional kernel
in scale-space as

$\mathscr{K}(\mathbf{u},\epsilon,t,\delta)=K^{2}_{\exp(-\epsilon)}(\mathbf{u})K^{1}_{\exp(-\delta)}(t),$

where $K^{2}$ and $K^{1}$ represent the two- and one-dimensional Gaussian
kernel functions (Silverman, 1986). Consider once again a counting measure
$\mathcal{N}$ that assigns 1 to each of the points
$\\{\mathbf{u}_{i},\log{(\epsilon_{i})},t_{i},\log{(\delta_{i})}\\}_{i=1}^{n}$,
and then, consider the convolution of the kernel with the measure,

$\displaystyle(\mathscr{K}\ast\mathcal{N})(\mathbf{u},\epsilon,t,\delta)$
$\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{5}}\mathscr{K}(\mathbf{u}-\mathbf{u}_{0},\epsilon-\epsilon_{0},t-t_{0},\delta-\delta_{0})\mathrm{d}\mathcal{N}(\mathbf{u}_{0},\epsilon_{0},t_{0},\delta_{0})$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\mathscr{K}(\mathbf{u}-\mathbf{u}_{i},\epsilon-\log(\epsilon_{i}),t-t_{i},\delta-\log(\delta^{\prime}_{i}))$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{n}K^{2}_{\epsilon_{i}\exp(-\epsilon)}(\mathbf{u}-\mathbf{u}_{i})K^{1}_{\delta_{i}\exp(-\delta)}(t-t^{\prime}_{i}),$

and evaluating the convolution in the hyperplane $\epsilon=\delta=0$, we
obtain

$(\mathscr{K}\ast\mathcal{N})(\mathbf{u},0,t,0)=\hat{\lambda}^{*}_{\epsilon,\delta}(\mathbf{u},t).$

Analogously, the edge-correction term may be expressed as a convolution sliced
at the corresponding hyperplane. The previous expressions constitute five-
dimensional functions, complicating things in terms of computational
resources. For example, the padding step to make the function periodic must
add $M_{1}M_{2}M_{3}M_{4}M_{5}$ zeros to the kernel array and the counting
measure array. Consequently, the memory expense and the number of operations
to calculate the transform, its inverse, and the kernel’s product with the
counting measure increase exponentially. Thus, it causes a conventional
computer to exceed its memory limits quickly, and therefore, the algorithm is
useless in practice. Scott (2015) indicated this problem in several dimensions
and recommended a practical limit of three dimensions; the advice is not to go
beyond the planar case developed by Davies and Baddeley (2018).

One of the key points for the adaptive estimation of the intensity is the
selection of the pilot intensity and the global bandwidth to use Abramson’s
rule. In this work, we have chosen spatial and temporal pilot intensities
separated and estimated using the traditional method with fixed bandwidth,
that is, the method through Fourier transforms in $\mathbb{R}$ for the
temporal case and in $\mathbb{R}^{2}$ for the spatial case. However, using a
non-separable estimator such as Eq. 4, could also be considered , and then
integrating to obtain the marginal intensities. The method is not the fastest
computationally speaking but could eventually improve the bandwidths for each
data point.

Since the edge correction is a much smoother surface than the intensity itself
as in most regions it usually takes the value of one; Davies and Baddeley
(2018) estimate this edge correction using a coarser mesh than the one used
for the numerator of the intensity. They then use simple interpolation
techniques to scale the surface to the intensity scale; this idea stems from
the need for faster estimations. Although we have used the same meshes for
both the intensity and the edge correction factor in the present work,
exploring this possibility of coarser evaluation grids for the edge correction
in space-time could be interesting.

To conclude, the adaptive spatio-temporal estimator presented here allows us
to face the first-order intensity estimation problem, especially in point
patterns with high concentrations of points in certain areas, that is, those
highly heterogeneous, providing each point with its bandwidth. Although the
complexity of the estimator is high, we can alleviate the computational burden
through the partitioning algorithm we have introduced. In addition, we also
provide the R code to easily and fast analyse point patterns that can arise in
many scientific research fields.

## References

  * Abramson (1982) [author] Abramson, Ian S.I. S. (1982). On Bandwidth Variation in Kernel Estimates-A Square Root Law. The Annals of Statistics 10 1217–1223.
  * Andela et al. (2019) [author] Andela, NielsN., Morton, Douglas CD. C., Giglio, LouisL., Paugam, RonanR., Chen, YangY., Hantson, StijnS., Van Der Werf, Guido RG. R. and Randerson, James TJ. T. (2019). The Global Fire Atlas of individual fire size, duration, speed and direction. Earth System Science Data 11 529–552.
  * Baddeley, Rubak and Turner (2015) [author] Baddeley, AdrianA., Rubak, EgeE. and Turner, RolfR. (2015). Spatial Point Patterns: Methodology and Applications with R. Chapman & Hall Interdisciplinary Statistics Series. CRC Press, Boca Raton, Florida.
  * Berman and Diggle (1989) [author] Berman, MarkM. and Diggle, Peter J.P. J. (1989). Estimating weighted integrals of the second-order intensity of a spatial point process. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 51 81-92.
  * Chaudhuri and Marron (2000) [author] Chaudhuri, ProbalP. and Marron, J. S.J. S. (2000). Scale Space View of Curve Estimation. The Annals of Statistics 28 408–428.
  * Cronie and Van Lieshout (2018) [author] Cronie, OO. and Van Lieshout, M. N. MM. N. M. (2018). A non-model-based approach to bandwidth selection for kernel estimators of spatial intensity functions. Biometrika 105 455-462.
  * Davies and Baddeley (2018) [author] Davies, Tilman M.T. M. and Baddeley, AdrianA. (2018). Fast computation of spatially adaptive kernel estimates. Statistics and Computing 28 937–956.
  * Davies and Hazelton (2010) [author] Davies, Tilman M.T. M. and Hazelton, Martin L.M. L. (2010). Adaptive kernel estimation of spatial relative risk. Statistics in Medicine 29 2423-2437.
  * Davies, Marshall and Hazelton (2018) [author] Davies, T. M.T. M., Marshall, J. C.J. C. and Hazelton, M. L.M. L. (2018). Tutorial on kernel estimation of continuous spatial and spatiotemporal relative risk. Statistics in Medicine 37 1191–1221.
  * Díaz-Avalos, Juan and Mateu (2013) [author] Díaz-Avalos, CarlosC., Juan, P.P. and Mateu, J.J. (2013). Similarity measures of conditional intensity functions to test separability in multidimensional point processes. Stochastic Environmental Research and Risk Assessment 27 1193-1205.
  * Diggle (1985) [author] Diggle, Peter JP. J. (1985). A Kernel Method for Smoothing Point Process Data. Journal of the Royal Statistical Society: Series C (Applied Statistics) 34 138-147.
  * Diggle (2013) [author] Diggle, Peter J.P. J. (2013). Statistical Analysis of Spatial and Spatio-Temporal Point Patterns, third ed. Chapman & Hall Monographs on Statistics & Applied Probability. CRC Press, Boca Raton, Florida.
  * Diggle, Rowlingson and Su (2005) [author] Diggle, Peter J.P. J., Rowlingson, BarryB. and Su, Ting-liT.-l. (2005). Point process methodology for on-line spatio-temporal disease surveillance. Environmetrics 16 423-434.
  * Diggle et al. (2013) [author] Diggle, Peter J.P. J., Moraga, PaulaP., Rowlingson, BarryB. and Taylor, BenjaminB. (2013). Spatial and Spatio-Temporal Log-Gaussian Cox Processes: Extending the Geostatistical Paradigm. Statistical Science 28 542-563.
  * Fernando and Hazelton (2014) [author] Fernando, W. T. P. SarojinieW. S. and Hazelton, Martin L.M. L. (2014). Generalizing the spatial relative risk function. Spatial and Spatio-temporal Epidemiology 8 1-10.
  * Fuentes-Santos, González-Manteiga and Mateu (2018) [author] Fuentes-Santos, IsabelI., González-Manteiga, WenceslaoW. and Mateu, JorgeJ. (2018). A first-order, ratio-based nonparametric separability test for spatiotemporal point processes. Environmetrics 29 1-18.
  * Gabriel and Diggle (2009) [author] Gabriel, EdithE. and Diggle, Peter J.P. J. (2009). Second-order analysis of inhomogeneous spatio-temporal point process data. Statistica Neerlandica 63 43–51.
  * Gabriel, Rowlingson and Diggle (2013) [author] Gabriel, EdithE., Rowlingson, BarryB. and Diggle, Peter J.P. J. (2013). stpp: An R Package for Plotting, Simulating and Analyzing Spatio-Temporal Point Patterns. Journal of Statistical Software 53(2) 1-29.
  * Ghorbani (2013) [author] Ghorbani, MohammadM. (2013). Testing the weak stationarity of a spatio-temporal point process. Stochastic Environmental Research and Risk Assessment 27 517-524.
  * Ghorbani et al. (2021) [author] Ghorbani, MohammadM., Vafaei, NafisehN., Dvořák, JiříJ. and Myllymäki, MariM. (2021). Testing the first-order separability hypothesis for spatio-temporal point patterns. Computational Statistics & Data Analysis 161 107245\.
  * Gneiting, Genton and Guttorp (2006) [author] Gneiting, TilmannT., Genton, Marc G.M. G. and Guttorp, PeterP. (2006). Geostatistical Space-Time Models, Stationarity, Separability and Full Symmetry. In Statistical Methods for Spatio-Temporal Systems, first ed. (BärbelB. Finkenstädt, LeonhardL. Held and ValerieV. Isham, eds.). Monographs on Statistics and Applied Probability 151–175. Chapman & Hall/CRC, New York.
  * González, Hahn and Mateu (2019) [author] González, Jonatan A.J. A., Hahn, UteU. and Mateu, JorgeJ. (2019). Analysis of tornado reports through replicated spatiotemporal point patterns. Journal of the Royal Statistical Society: Series C (Applied Statistics) 69 3–23.
  * González et al. (2016) [author] González, Jonatan A.J. A., Rodríguez-Cortés, Francisco J.F. J., Cronie, OttmarO. and Mateu, JorgeJ. (2016). Spatio-temporal point process statistics: A review. Spatial Statistics 18 505-544.
  * Hall and Marron (1988) [author] Hall, PeterP. and Marron, J. S.J. S. (1988). Variable window width kernel estimates of probability densities. Probability Theory and Related Fields 80 37–49.
  * Illian et al. (2008) [author] Illian, J.J., Penttinen, P. A.P. A., Stoyan, H.H. and Stoyan, D.D. (2008). Statistical Analysis and Modelling of Spatial Point Patterns. Statistics in Practice. Wiley.
  * Jones (1993) [author] Jones, M. C.M. C. (1993). Simple boundary correction for kernel density estimation. Statistics and Computing 3 135-146.
  * Jones and Lotwick (1983) [author] Jones, M. C.M. C. and Lotwick, H. W.H. W. (1983). On the errors involved in computing the empirical characteristic function. Journal of Statistical Computation and Simulation 17 133-149.
  * Koch, Ohser and Schladitz (2003) [author] Koch, KarstenK., Ohser, JoachimJ. and Schladitz, KatjaK. (2003). Spectral theory for random closed sets and estimating the covariance via frequency space. Advances in Applied Probability 35 603–613.
  * Loader (1999) [author] Loader, CliveC. (1999). Local Regression and Likelihood. Springer, New York.
  * Marshall and Hazelton (2010) [author] Marshall, Jonathan C.J. C. and Hazelton, Martin L.M. L. (2010). Boundary kernels for adaptive density estimators on regions with irregular boundaries. Journal of Multivariate Analysis 101 949-963.
  * Møller and Ghorbani (2012) [author] Møller, J.J. and Ghorbani, MM. (2012). Aspects of second-order analysis of structured inhomogeneous spatio-temporal point processes. Statistica Neerlandica 66 472-491.
  * Møller, Syversveen and Waagepetersen (1998) [author] Møller, JesperJ., Syversveen, Anne RandiA. R. and Waagepetersen, Rasmus PlengeR. P. (1998). Log Gaussian Cox Processes. Scandinavian Journal of Statistics 25 451-482.
  * Sain (2002) [author] Sain, Stephan R.S. R. (2002). Multivariate locally adaptive density estimation. Computational Statistics & Data Analysis 39 165-186.
  * Sain and Scott (1996) [author] Sain, Stephan R.S. R. and Scott, David W.D. W. (1996). On Locally Adaptive Density Estimation. Journal of the American Statistical Association 91 1525-1534.
  * Schoenberg (2004) [author] Schoenberg, Frederic PaikF. P. (2004). Testing Separability in Spatial-Temporal Marked Point Processes. Biometrics 60 471–481.
  * Scott (2015) [author] Scott, David WD. W. (2015). Multivariate density estimation: theory, practice, and visualization, second ed. Wiley series in probability and statistics. John Wiley & Sons.
  * Sheather and Jones (1991) [author] Sheather, S. J.S. J. and Jones, M. C.M. C. (1991). A Reliable Data-Based Bandwidth Selection Method for Kernel Density Estimation. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 53 683–690.
  * Silverman (1982) [author] Silverman, B. W.B. W. (1982). Algorithm AS 176: Kernel Density Estimation Using the Fast Fourier Transform. Journal of the Royal Statistical Society: Series C (Applied Statistics) 31 93–99.
  * Silverman (1986) [author] Silverman, Bernard WB. W. (1986). Density estimation for statistics and data analysis. Chapman & Hall Monographs on Statistics & Applied Probability. CRC press, Boca Raton, Florida.
  * R Core Team (2020) [author] R Core Team (2020). R: A Language and Environment for Statistical Computing R Foundation for Statistical Computing, Vienna, Austria.
  * Terrell (1990) [author] Terrell, George R.G. R. (1990). The Maximal Smoothing Principle in Density Estimation. Journal of the American Statistical Association 85 470-477.
  * van Lieshout (2011) [author] van Lieshout, M. N. M.M. N. M. (2011). On Estimation of the Intensity Function of a Point Process. Methodology and Computing in Applied Probability 14 567-578.
  * Wand (1994) [author] Wand, M. P.M. P. (1994). Fast Computation of Multivariate Kernel Estimators. Journal of Computational and Graphical Statistics 3 433–445.
  * Wand and Jones (1995) [author] Wand, Matt PM. P. and Jones, M ChrisM. C. (1995). Kernel smoothing. Monographs on statistics and applied probability. Chapman & Hall, New York.