# From Average Embeddings To Nearest Neighbor Search111Research supported in
part by NSF grants (CCF-2008733, CCF-1617955, and CCF-1740833), and Simons
Foundation (#491119).

Alexandr Andoni
CS and DSI
Columbia University  David Cheikhi
Columbia University

###### Abstract

In this note, we show that one can use average embeddings, introduced recently
in [Nao20], to obtain efficient algorithms for approximate nearest neighbor
search. In particular, a metric $X$ embeds into $\ell_{2}$ on average, with
distortion $D$, if, for any distribution $\mu$ on $X$, the embedding is $D$
Lipschitz and the (square of) distance does not decrease on average (wrt
$\mu$). In particular existence of such an embedding (assuming it is
efficient) implies a $O(D^{3})$ approximate nearest neigbor search under $X$.
This can be seen as a strengthening of the classic (bi-Lipschitz) embedding
approach to nearest neighbor search, and is another application of data-
dependent hashing paradigm.

##  1 Introduction

In the $c$-Approximate Nearest Neighbor Search ($c$-ANN) problem, we
preprocess a dataset $P$ living in some metric space $(X,d_{X})$ for some
threshold $r>0$, such that, for a given query $q\in X$, we can efficiently
find a point in $p\in P$ at distance $d_{X}(p,q)\leq cr$, as long as there
exists a point $p^{*}\in P$ at distance $d_{X}(q,p^{*})\leq r$. In many
applications, the metric space is high-dimensional, with the $d$-dimensional
Euclidean space $\Re^{d}$ being the most well-studied. See, e.g., survey
[AIR18] for further background on high-dimensional ANN problem.

A common approach for solving the ANN problem is via Locality Sensitive
Hashing (LSH) [HIM12], which can be viewed as a random, oblivious space
partition. LSH mainly applies to the Euclidean ($\ell_{2}$) and Hamming
($\ell_{1}$) spaces. For example, for any fixed $\epsilon>0$, one can solve
ANN under $\ell_{2}$ with $c=\sqrt{1/\epsilon}+o(1)$ approximation,
$O(n^{\epsilon}d)$ query time and $O(n^{1+\epsilon}d)$ space (we call such
query/space parameters as ”efficient” henceforth).

Beyond $\ell_{1}$ and $\ell_{2}$ metrics, there has been much less progress
and the best bounds are usually far from understood. A long-standing approach
is via metric embeddings: a map from a metric $(X,d_{X})$ into an ”easier”
space, say, $(\Re^{d},\ell_{2})$, while approximately preserving all the
distances:

######  Definition 1.1 (bi-Lipschitz embedding into $\ell_{2}$).

A map $f:X\to\ell_{2}$ is an embedding into $\ell_{2}$, with distortion $D\geq
1$, iff for any $x,y\in X$ we have that

$\|f(x)-f(y)\|_{2}\leq d_{X}(x,y)\leq D\cdot\|f(x)-f(y)\|_{2}.$

Using such an embedding for a metric $(X,d_{X})$, one can immediately extend
the ANN under $\ell_{2}$ algorithms to ANN under $X$, with an extra factor $D$
approximation. Indeed, this approach yields best-known ANN algorithms for
spaces such as the edit distance [OR07]. However, there are also known lower
bounds for the best possible distortion, including in the case of edit
distance [KN06, KR06]. A possible extension of the embedding approach is to
embed a snowflake of the metric $(X,d_{X})$ which means embedding
$(X,d_{X}^{\omega})$ for some $\omega<1$, though oftentimes there are lower
bounds for this approach as well (albeit weaker), see, e.g., [AK10] for our
example of edit distance.

A recent qualitative development is that of data-dependent hashing/space
partitioning methods. For $\ell_{2}$ and $\ell_{1}$ spaces, [AINR14, AR15,
ALRW17] developed data-dependent LSH that gave a better approximation: e.g.,
$c=\sqrt{\tfrac{1+\epsilon}{2\epsilon}}+o(1)$ for $\ell_{2}$; also see [Ind98]
for an earlier example of data-dependent technique for $\ell_{\infty}$.

For metrics beyond $\ell_{1},\ell_{2}$, [ANN+18] developed a new data-
dependent approach, based on the following geometric property of the space
$(X,d_{X})$, termed cutting modulus:

######  Definition 1.2 (Cutting modulus).

For fixed $\epsilon>0$, the cutting modulus of the metric ${\cal
X}=(X,d_{X})$, denoted $\Xi({\cal X},\epsilon)$ is the infimum $D$ such that
the following holds for any $r>0$ and any integer $m>1$. Fix a pointset
$V=\\{x_{1},\ldots x_{m}\\}\subset X$, and consider a positively-weighted
graph $G=(V,E,w)$ where $w_{pq}>0\implies d_{X}(p,q)\leq r$ and $\sum_{p\in
V}\sum_{q\in V}w_{p,q}=1$ (i.e., it’s a distribution over the ordered pairs
$(p,q)$). Then, one of the following must hold:

  * •

Either there exists some point $z\in X$ such that $\sum_{p:d_{X}(z,p)\leq
Dr}\sum_{q\in V}w_{p,q}\geq 1/2$ (i.e., there is a ball of radius $Dr$
containing 1/2 mass of the pointset),

  * •

Or there exists a non-trivial cut $S\subset V$ with conductance
$\Phi(S)\leq\epsilon$.

The main result of [ANN+18] is that any metric $X$ with cutting modulus $D$
admits a $D$-ANN, albeit with a couple important caveats. First, without
further properties, the query algorithm is efficient only in the cell-probe
model, as it involves operating with an exponentially-sized graph $G$ (think a
net in $\cal X$). Nonetheless, as [ANN+18] show, one can use a slightly
stronger notion of cutting modulus to remove this caveat, for example
obtaining efficient ANN for $\ell_{p}$ and Schatten-$p$ norms, with
$O(p/\epsilon)$ approximation, for any $p\geq 2$. Second, the preprocessing
time is exponential and seems unavoidable for this approach.

There are currently few bounds on the cutting modulus, and most of them follow
from the theory of non-linear spectral gaps, whose systematic study was
started in [MN14]; see also [Nao14, MN+15]. Perhaps most striking is that the
cutting modulus of any $d$-dimensional norm $(\Re^{d},\|\cdot\|_{X})$ is
$O(\tfrac{\log d}{\epsilon})$, following the bound on the non-linear spectral
gap from [Nao17].

The work of [Nao20] introduced the notion of average embedding and proved a
connection to non-linear spectral gap bounds. An average embedding relaxes the
bi-Lipschitz embedding from above in that the contraction holds only on
average with respect to a fixed given measure (and hence the embedding is
”data dependent”):

######  Definition 1.3 (Average embedding).

For a metric ${\cal X}=(X,d_{X})$, a map $f:X\rightarrow\ell_{2}$ is an
$q$-average embedding, for a dataset $P\subset X$, with distortion $D$ if the
following is satisfied:

$\displaystyle\|f(p)-f(q)\|_{2}$ $\displaystyle\leq D\cdot
d_{X}(p,q)\quad\text{ for all $p,q\in X$}$ $\displaystyle\sum_{p,q\in
P}\|f(p)-f(q)\|_{2}^{q}$ $\displaystyle\geq\sum_{p,q\in P}d_{X}(p,q)^{q}.$

When $q=2$, we call it simply average embedding.

In particular, [Nao20] showed that the average embedding distortion for
$1/2$-snowflake of a $d$-dimensional normed space is $O(\sqrt{\log d})$.
Similarly, [Nao14, Nao20] show that the average embedding distortion for
$\ell_{p}$ is $O(p)$, for $p>2$. Both of these embeddings are non-
constructive, proven by duality.

#### Our result.

In this note, we show that one can use average embeddings directly to solve
ANN.

######  Theorem 1.4.

Suppose a metric $X$ has an average embedding $f$ with distortion $\sqrt{D}$,
where time to construct $f$ is $T_{P}$ and time to compute $f$ on a point $q$
is $T_{c}$. Then there is a $c$-ANN with query time $O(n^{\epsilon}\cdot
T_{c})$, space $n^{O(1)}$, and preprocessing time $n^{O(1)}\cdot T_{P}$, and
approximation $c=O(D^{3}/\epsilon)$.

Our approach avoids some of the aforementioned caveats of using cutting
modulus, as long as we can construct an efficient average embedding for a
given metric $X$. In particular, there is no more need for exponential time
preprocessing (again assuming efficient average embedding map).

At high-level, our algorithm is similar to the ANN approach via standard bi-
Lipschitz embedding: embed the space $X$ into $\ell_{2}$ and then use an
efficient ANN for $\ell_{2}$, though the actual algorithm is a little bit more
nuanced. The overall algorithm is ”data-dependent” because the average
embedding $f$ is data-dependent (note that $f$ from above may depend on the
dataset $P$). In Appendix A we also propose an explicit, efficient average
embedding for $\ell_{p}$, $p>2$, albeit we leave it as a conjecture to prove
its second property.

#### Independent work.

In independent work, [KNT21] showed a similar result, obtaining a tighter
bound of $O(D\log D)$-ANN. Their algorithm requires one more step: to
construct a weak average embedding out of a generic average embedding. They
also show how to compute the $1$-average embedding for $\ell_{p}$ with
distortion $O(p)$ (thus also resolving a related version of our conjecture),
and hence obtain an $O(p)$-ANN with polynomial-time preprocessing (removing
also the $O(\log p)$ from their generic reduction). Their algorithm thus
improves over ANN from [BG19, ANN+18].

##  2 ANN Algorithm from Average Embedding

We now describe our main ANN algorithm that uses an average embedding. Our
algorithm is solves a special case of the ANN problem, termed bounded ANN: the
$\beta$-bounded ANN is the problem where we are guaranteed that all distances
within the dataset $P$ verify $r\leq d_{X}(p,q)\leq\beta\cdot cr$. We note
that it is enough to solve $\beta$-bounded ANN problem, for $\beta=18$, due to
the result of [BG19].

Our algorithm uses an LSH scheme for $\ell_{2}$, formally defined as follows.

######  Definition 2.1.

(LSH) For a metric $X$, a family of hash functions $h:X\to U$ is
$(p_{1},p_{2},r,cr)$-LSH if, for all $x,y\in X$:

  * •

if $d_{X}(x,y)\leq r$, then $\Pr[h(x)=h(y)]\geq p_{1}$,

  * •

if $d_{X}(x,y)>cr$, then $\Pr[h(x)=h(y)]<p_{2}$.

The exponent of the LSH scheme is defined as $\rho=\tfrac{\log 1/p_{1}}{\log
1/p_{2}}$.

The metric $X=\ell_{2}$ admits $(p_{1},p_{2},r,cr)$-LSH, for $\rho=1/c$ with
$p_{1}$ arbitrarily close to 1 [DIIM04]. We also assume that $r=1$ henceforth
for simplicity.

Fix some $\lambda,w\geq 1$ tbd. We construct the data structure for a
$\beta$-bounded pointset $P$ recursively as a collection of randomized trees,
constructed iid. A tree recursively partitions the dataset until the size of
the dataset becomes constant (below think of $P$ is the dataset assigned to
the current node). There are two type of nodes:

  * •

If there exists $x_{0}\in P$ such that $|P\cap B_{X}(x_{0},\lambda
D))|>|P|/8$, we store $x_{0}$ and a point $p_{0}\in P\cap B_{X}(x_{0},\lambda
D)$. For the rest of the points, we create a single child inheriting
$P\setminus B_{X}(x_{0},\lambda D))$, and recurse the construction on it.

  * •

If there is no $x_{0}\in P$ such that $|P\cap B_{X}(x_{0},\lambda D)|>|P|/8$,
then we compute an average embedding $f:X\to\ell_{2}$ for $P$. Then we sample
a $(p_{1},p_{2},D,wD)$-LSH $h:\ell_{2}\to U$ and use it to partition the
pointset $P$. In particular, each part $k\in U$ generates a child node $k$
that inherits the points $p\in P$ with $h(f(p))=k$ (storing only non-empty
parts as usual).

For a given query point $q$, we traverse the tree as follows:

  * •

if the node is of the first type: if $d(x_{0},q)<\sqrt{D}$, return the
associated $p_{0}$, and otherwise recurse into the single child;

  * •

if the node is of the second type: branch into the child node $k=h(f(q))$.

  * •

in a leaf node: check the distance to all stored nodes and return a valid
answer if there exists one.

######  Theorem 2.2.

For any $\epsilon>0$, there is $\lambda,w$, and $c=O(D^{3}/\epsilon)$ such
that the tree constructed as above has a probability at least $n^{-\epsilon}$
to find a $cr$-near neighbor, assuming there’s an $r$-near neighbor.

Hence, we can build a $\beta$-bounded ANN which constructs $n^{\epsilon}$
trees as above, and uses space $n^{1+O(\epsilon)}$ for approximation
$c=O(D^{3}/\epsilon)$. The reduction from [BG19] then yields an algorithm for
the vanilla ANN.

###  2.1 Analysis: Proof of Theorem 2.2

First we note that nodes of the first type do not affect the correctness. So
in the rest we fix a node of the tree above, and suppose we are in the second
case: there is no $x_{0}\in X$ such that $|P\cap B_{X}(x_{0},\lambda
D)|>|P|/8$. We further assume that $P$ contains a near neighbor point $p^{*}$,
with $d(q,p^{*})\leq 1$. We first prove the following lemma bounding the
number of points that are ”close” to $q$ after the embedding $f$.

######  Lemma 2.3.

Consider a set $P$ satisfying:

  * •

$|B_{X}(x,\lambda D)\cap P|\leq p|P|$ for all $x\in P$;

  * •

For all points $x,y\in P$, $d_{X}(x,y)\leq\beta c$.

Consider a query $q$ at distance at most 1 from some point in $P$. Then, we
have that

$\left|B_{\ell_{2}}(f(q),wD)\cap f(P)\right|\leq
1-\dfrac{(1-p)\lambda^{2}}{4\beta^{2}c^{2}},$

as long as $(1-p)\lambda^{2}=\Omega(w^{2})$, and $\beta c\geq\lambda$.

###### Proof.

Let $\alpha=\tfrac{\left|B_{\ell_{2}}(f(q),wD)\cap f(P)\right|}{|P|}$. By the
second property of the average embedding for $P$:

$\displaystyle\sum_{x,y\in\mathcal{P}}d(x,y)^{2}$
$\displaystyle\leq\sum_{x,y\in\mathcal{P}}\|f(x)-f(y)\|^{2}$
$\displaystyle\implies\sum_{x,y\in\mathcal{P}}d(x,y)^{2}$
$\displaystyle\leq(1-\alpha^{2})\beta^{2}c^{2}D^{2}n^{2}+4\alpha^{2}w^{2}D^{2}n^{2}$
using $D$-Lipschitzness of $f$, and that for an $\alpha$ ratio of $x\in P$,
$\|f(x)-f(q)\|\leq wD$. $\displaystyle\implies n^{2}(1-p)\lambda^{2}D^{2}$
$\displaystyle\leq\beta^{2}c^{2}D^{2}n^{2}-\alpha^{2}(\beta^{2}c^{2}D^{2}-4w^{2}D^{2})n^{2}$
since $d_{X}(x,y)\geq\lambda D$ for a $1-p$ ratio of $x,y\in P$.
$\displaystyle\implies\alpha^{2}$
$\displaystyle\leq\frac{\beta^{2}c^{2}D^{2}-(1-p)\lambda^{2}D^{2}}{\beta^{2}c^{2}D^{2}-4w^{2}D^{2}}$
$\displaystyle\implies\alpha^{2}$
$\displaystyle\leq\frac{\beta^{2}c^{2}-(1-p)\lambda^{2}}{\beta^{2}c^{2}-4w^{2}}$
$\displaystyle\implies\alpha^{2}$
$\displaystyle\leq(1-\tfrac{(1-p)\lambda^{2}}{\beta^{2}c^{2}})(1+O(w^{2}/\beta^{2}c^{2}))$
$\displaystyle\implies\alpha^{2}$ $\displaystyle\leq
1-\tfrac{(1-p)\lambda^{2}-O(w^{2})}{\beta^{2}c^{2}}.$
$\displaystyle\implies\alpha$ $\displaystyle\leq
1-\tfrac{(1-p)\lambda^{2}}{4\beta^{2}c^{2}},$

for $(1-p)\lambda^{2}\geq\tfrac{1}{2}O(w^{2})$. ∎

We use the lemma from above with $p=1/8$, $\beta=18$, and $c=\lambda D$. After
the embedding $f$, it is guaranteed that $|f(P)\cap
B(f(q),wD)|>|P|\cdot(1-\Omega\left(\dfrac{1}{D^{2}}\right))$. Using an
$(p_{1},p_{2},D,wD)$-LSH in the embedded space, the fraction of points that
collides with the query is:

$p_{2}^{\prime}\leq
1-\Omega\left(\dfrac{1}{D^{2}}\right)+p_{2}\left(1-(1-\Omega\left(\dfrac{1}{D^{2}}\right))\right)\leq
1-(1-p_{2})\dfrac{\Omega(1)}{D^{2}}.$

Therefore, embedding the query and the dataset and applying an
$(p_{1},p_{2},D,wD)$-LSH is like applying an $(p_{1},p_{2}^{\prime},1,\lambda
D)$-LSH in the original space. Therefore, we pick

$p_{1}=1-\epsilon(1-p_{2}^{\prime})=1-\epsilon(1-p_{2})\dfrac{\Theta(1)}{D^{2}}.$

Such a choice is possible as long as $w\geq\Theta(D^{2}/\epsilon)$. As, from
above, we require $\lambda\geq\Omega(w)$, we have that
$\lambda\geq\Omega(D^{2}/\epsilon)$. Finally, we get that the approximation
ratio is $c=\lambda D=O(D^{3}/\epsilon)$. We obtain $\rho=\epsilon$ and hence
$n^{-\epsilon}$ success probability.

## References

  * [AINR14] Alexandr Andoni, Piotr Indyk, Huy L. Nguyen, and Ilya Razenshteyn.  Beyond locality-sensitive hashing.  In Proceedings of the ACM-SIAM Symposium on Discrete Algorithms (SODA), 2014.  Full version at http://arxiv.org/abs/1306.1547.
  * [AIR18] Alexandr Andoni, Piotr Indyk, and Ilya Razenshteyn.  Approximate nearest neighbor search in high dimensions.  2018\.  This survey accompanied the ICM 2018 talk of Piotr Indyk. arXiv preprint arXiv:1806.09823.
  * [AK10] Alexandr Andoni and Robert Krauthgamer.  The computational hardness of estimating edit distance.  SIAM J. Comput. (SICOMP), 39(6):2398–2429, 2010.  Previously in FOCS’07.
  * [ALRW17] Alexandr Andoni, Thijs Laarhoven, Ilya P. Razenshteyn, and Erik Waingarten.  Optimal hashing-based time-space trade-offs for approximate near neighbors.  In Proceedings of the ACM-SIAM Symposium on Discrete Algorithms (SODA), 2017.
  * [ANN+18] Alexandr Andoni, Assaf Naor, Aleksandar Nikolov, Ilya P. Razenshteyn, and Erik Waingarten.  Data-dependent hashing via nonlinear spectral gaps.  In Proceedings of the Symposium on Theory of Computing (STOC), 2018\.
  * [AR15] Alexandr Andoni and Ilya Razenshteyn.  Optimal data-dependent hashing for approximate near neighbors.  In Proceedings of the Symposium on Theory of Computing (STOC), 2015\.  Full version at http://arxiv.org/abs/1501.01062.
  * [BG19] Yair Bartal and Lee-Ad Gottlieb.  Approximate nearest neighbor search for $\ell_{p}$-spaces (2<p<$\infty$).  Theor. Comput. Sci., 757:27–35, 2019.
  * [DIIM04] Mayur Datar, Nicole Immorlica, Piotr Indyk, and Vahab Mirrokni.  Locality-sensitive hashing scheme based on p-stable distributions.  In Proceedings of the ACM Symposium on Computational Geometry (SoCG), 2004.
  * [HIM12] Sariel Har-Peled, Piotr Indyk, and Rajeev Motwani.  Approximate nearest neighbor: Towards removing the curse of dimensionality.  Theory of Computing, 1(8):321–350, 2012.
  * [Ind98] Piotr Indyk.  Faster algorithms for string matching problems: matching the convolution bound.  Proceedings of the Symposium on Foundations of Computer Science (FOCS), 1998.
  * [KN06] Subhash Khot and Assaf Naor.  Nonembeddability theorems via Fourier analysis.  Math. Ann., 334(4):821–852, 2006.  Preliminary version appeared in FOCS’05.
  * [KNT21] Deepanshu Kush, Aleksandar Nikolov, and Haohua Tang.  Near neighbor search via efficient average distortion embeddings.  In Proceedings of the ACM Symposium on Computational Geometry (SoCG), 2021.  To appear.
  * [KR06] Robert Krauthgamer and Yuval Rabani.  Improved lower bounds for embeddings into $L_{1}$.  In Proceedings of the ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 1010–1017, 2006.
  * [MN14] Manor Mendel and Assaf Naor.  Nonlinear spectral calculus and super-expanders.  Publications mathématiques de l’IHÉS, 119(1):1–95, 2014\.
  * [MN+15] Manor Mendel, Assaf Naor, et al.  Expanders with respect to hadamard spaces and random graphs.  Duke Mathematical Journal, 164(8):1471–1548, 2015.
  * [Nao14] Assaf Naor.  Comparison of metric spectral gaps.  Analysis and Geometry in Metric Spaces, 2(1), 2014.
  * [Nao17] Assaf Naor.  A spectral gap precludes low-dimensional embeddings.  In Proceedings of the ACM Symposium on Computational Geometry (SoCG), 2017.  Also at https://arxiv.org/abs/1611.08861.
  * [Nao20] Assaf Naor.  An average John theorem.  arXiv preprint 1905.01280, 2020.  To appear in Geometry & Topology.
  * [OR07] Rafail Ostrovsky and Yuval Rabani.  Low distortion embedding for edit distance.  J. ACM, 54(5), 2007.  Preliminary version appeared in STOC’05.

##  Appendix A Average Embedding for $\ell_{p}$

We consider the following average embedding from $\ell_{p}$ to $\ell_{2}$. Let

$h(x)_{i}=\operatorname{sign}(x_{i})|x_{i}|^{p/2}$

and

$f(x)=\dfrac{h(x)}{\|h(x)\|_{2}}\|x\|_{p}=h(x)\dfrac{\|x\|_{p}}{\|x\|_{p}^{p/2}}=h\left(\dfrac{x}{\|x\|_{p}}\right)\|x\|_{2}$

######  Theorem A.1.

The embedding $f:\ell_{p}\rightarrow\ell_{2}$ is $p+1$ Lipschitz.

###### Proof.

We notice that $h$ is the Mazur map from $\ell_{p}$ to $\ell_{2}$. The
embedding $f$ first normalizes its input, applies the Mazur map and the
rescale it. We will therefore make use of the fact that the Mazur map is
$\frac{p}{2}$-Lipscitz on the unit sphere.

$\displaystyle\|f(x)-f(y)\|_{2}$
$\displaystyle=\left\|\dfrac{h(x)}{\|h(x)\|_{2}}\|x\|_{p}-\dfrac{h(y)}{\|h(y)\|_{2}}\|y\|_{p}\right\|_{2}$
$\displaystyle=\left\|\|x\|_{p}\left(\dfrac{h(x)}{\|h(x)\|_{2}}-\dfrac{h(y)}{\|h(y)\|_{2}}\right)+\dfrac{h(y)}{\|h(y)\|_{2}}\left(\|x\|_{p}-\|y\|_{p}\right)\right\|_{2}$
$\displaystyle\leq\|x\|_{p}\left\|\dfrac{h(x)}{\|h(x)\|_{2}}-\dfrac{h(y)}{\|h(y)\|_{2}}\right\|_{2}+\left\|\dfrac{h(y)}{\|h(y)\|_{2}}\right\|_{2}\left|\|x\|_{p}-\|y\|_{p}\right|$
$\displaystyle\leq\|x\|_{p}\left\|\dfrac{h(x)}{\|h(x)\|_{2}}-\dfrac{h(y)}{\|h(y)\|_{2}}\right\|_{2}+\left\|x-y\right\|_{p}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\|x\|_{p}\left\|h\left(\dfrac{x}{\|x\|_{p}}\right)-h\left(\dfrac{y}{\|y\|_{p}}\right)\right\|_{2}+\left\|x-y\right\|_{p}$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}\dfrac{p}{2}\|x\|_{p}\left\|\dfrac{x}{\|x\|_{p}}-\dfrac{y}{\|y\|_{p}}\right\|_{p}+\left\|x-y\right\|_{p}$
$\displaystyle=\dfrac{p}{2}\left\|x-\dfrac{\|x\|_{p}}{\|y\|_{p}}y\right\|_{p}+\left\|x-y\right\|_{p}$
$\displaystyle=\dfrac{p}{2}\|x-y\|+\dfrac{p}{2}\left|\dfrac{\|x\|_{p}}{\|y\|_{p}}-1\right|\|y\|_{p}+\left\|x-y\right\|_{p}$
$\displaystyle\leq(p+1)\left\|x-y\right\|_{p}$

Where (a) uses the fact that $h(x)/\|h(x)\|_{2}=h(x/\|x\|_{p})$ and (b) uses
the fact that $h$ is $p/2$-Lipschitz on the unit sphere. ∎

To complete the proof that the map $f$ is average embedding we need the
following inequality, left as a conjecture.

######  Conjecture A.2.

For any $x_{1},\dots,x_{n}$, we can find $z$ such that $x\rightarrow f(x-z)$
verifies, for some universal constant $C>0$:

$\sum_{i,j}\|f(x_{i}-z)-f(x_{j}-z)\|_{2}^{2}\geq
C\sum_{i,j}\|x_{i}-x_{j}\|_{p}^{2}.$


[hyperref] captionUnknown document class

[1]Adam J. Stewart [2]Lijun Liu

1]Department of Computer Science, University of Illinois at Urbana-Champaign,
201 North Goodwin Avenue, Urbana, 61801, IL, USA 2]Department of Earth Science
& Environmental Change, University of Illinois at Urbana-Champaign, 1301 West
Green Street, Urbana, 61801, IL, USA

# Hemispheric Dichotomy of Mantle Dynamics Revealed by Machine Learning

<EMAIL_ADDRESS>Yanchong Li<EMAIL_ADDRESS>Zebin Cao
<EMAIL_ADDRESS><EMAIL_ADDRESS>[  [

###### Abstract

Past efforts to interrogate the mantle’s contribution to Earth’s topography
suffer from inadequate in-situ measurements of true bathymetry and the
involvement of an empirical plate model whose presumed lithospheric density
profile interferes with the interpretation of other mantle forces. Here, we
introduce a machine learning algorithm that estimates the oceanic residual
topography based solely on surface and crustal attributes, providing a more
objective proxy for the topographic contribution from the mantle. Ablation
studies show that seafloor age is the most important factor, while other
properties help further improve the fit to bathymetry. The resulting residual
topography has notably smaller amplitudes than previous estimates, indicating
that more surface features can be explained by crustal properties than
previously thought. This exercise, designed to allow detection of long-
wavelength signals, uncovers a more prominent degree-one topographic component
than the degree-two pattern, both with much smaller amplitudes than previously
thought. We suggest that this result reflects a hemispheric contrast in mantle
dynamics, especially the two large low shear-wave velocity provinces (LLSVPs).
This inference is further strengthened by quantifying the respective
topographic contributions from lithospheric isostasy and dynamic topography.
We conclude that mantle convection beneath the Atlantic is more vigorous than
that beneath the Pacific, likely due to their different thermal-chemical
states.

Although the operation of plate tectonics on Earth’s surface is well accepted,
how the underlying mantle functions remains less clear. One important approach
to deciphering the style of mantle convection is through examining ocean
bathymetry. The systematic deepening of Earth’s oceans away from mid-ocean
ridges is generally considered isostatic subsidence of the lithosphere as it
ages. This process can be described through the plate cooling model (PSM
parsons1977analysis, GDH1 stein1992model, H13 hasterok2013heat), while
deviations from this subsidence trend are usually attributed to non-conductive
processes beneath the crust, a signal termed residual (observation) or dynamic
(modeling) topography.

Among the various mantle dynamic processes affecting topography, the existence
and amplitude of the degree-2 superswell due to mantle upwelling originated
from the two large low shear-wave velocity provinces (LLSVPs) are contentious.
Geoid modeling with buoyant LLSVPs favors large-amplitude (1–2 km) dynamic
uplift at this wavelength hager1985lower; yang2016dynamic, a prediction
challenged by recent in-situ measurement of oceanic residual topography
hoggard2020global. However, the latter study was also called into question by
its spatially limited data coverage and the adoption of the empirical plate
cooling model as a reference level, both of which may prohibit representation
of the true convective state of the underlying mantle yang2016dynamic;
yang2017oceanic. More recent works steinberger2016topography; davies2019earth,
both based on the extrapolated residual topography hoggard2016global, suggest
that incorporation of lithospheric density structures in the mantle convection
calculation could partially reconcile the above debate, although the
respective contributions from the lithosphere and the convective mantle to
surface topography still remain inconclusive. Another complication in such
exercises is the underestimated density of the continental mantle lithosphere,
which may also affect the degree-2 pattern of surface topography
wang2022topography; wang2022geoid.

To overcome the sparsity of in-situ measurements and to avoid a priori
assumptions about the density structure of oceanic lithosphere, we revisit the
problem of oceanic residual topography using a machine learning (ML)
algorithm, which effectively covers all ocean basins and utilizes data
constraints largely independent of the mantle. The result reveals a power
spectrum significantly different from previous studies, with a prominent
hemispheric rather than degree-2 pattern of residual topography and another
peak around 5–14 degrees. These features outline a different picture of mantle
convection and lithosphere dynamics from the traditional view.

##  1 Estimating Earth’s residual topography

Our goal is to obtain an objective representation of the global oceanic
residual topography with as little mantle contamination as possible. ML using
only surface and crustal data could potentially achieve this objective. By
training an ML model to associate the measurable quantities with the observed
bathymetry, we aim to uncover topographic features that are not determined by
shallow crustal properties.

Here, we consider all major structural features (seafloor age, layer
thickness, seismic velocity, and density) and vertical layers (water, ice,
sediment, crust, and mantle) of available data that could influence
bathymetry. In practice, we utilize constraints for observed topography from
ETOPO1 amante2009etopo1, seafloor age from multiple recent plate
reconstructions seton2020global; muller2019global; muller2016ocean;
muller2013seawater; muller2008age, and other structural features and layer
properties from CRUST1.0 laske2013update. All data are projected onto the same
1° resolution as CRUST1.0. Using the CRUST1.0 crustal type add-on, we filter
out data related to continental crust and oceanic plateaus. To prevent the ML
model from memorizing the long-wavelength residual topography whose dimension
is greater than any individual plate as well as to avoid spatial dependence
between data points, we divide the Earth into 7 major tectonic plates: Africa,
Antarctica, Australia, Eurasia, North America, South America, and Pacific
(Extended Data Fig. 1). The remaining plates are merged into the smallest
neighboring plate in order to achieve similarly sized units.

Our training and evaluation procedure uses grouped cross-validation. For each
tectonic plate, a model is trained from scratch on the other 6 plates, and
then used to make predictions on the 7th plate. This process is repeated for
all 7 plates until a prediction is made for every bathymetric data point on
Earth. Once predictions have been made for all data points, the predictions
are compared to the ground truth from the CRUST1.0 dataset. We compute and
report root mean square error (RMSE) and coefficient of determination (R2).
Three ML models—ridge regression (RR) tikhonov1977solutions, support vector
regression (SVR) drucker1996support, and multi-layer perceptrons (MLP)
haykin1994neural—are tested. More descriptions of the data processing and
training procedures are provided in Methods and Extended Data Table 2.

Table 1: Performance metrics for plate and ML models. Model | RMSE (km) | R2
---|---|---
HS | 1.051 | 0.214
PSM | 0.689 | 0.662
GDH1 | 0.634 | 0.714
H13 | 0.652 | 0.697
RR | 0.592 | 0.751
SVR | 0.564 | 0.773
MLP | 0.476 | 0.839

Regarding the performance of the three ML models in capturing bathymetric
information based on the same inputs (Table 1), MLP achieves the lowest RMSE
(0.476 km) and highest R2 (0.839). To further evaluate the effectiveness of ML
models in capturing topographic information, we compare their performance with
that of four different physical models, further described in Methods. As
expected, the half-space (HS) cooling model turcotte2002geodynamics has the
worst performance, with the highest RMSE and lowest R2 among all models. Among
the three plate models, GDH1 has the best performance. However, all machine
learning models we tested have significantly better performance due to their
ability to capture information from more crustal layers and features than the
plate models. Overall, ML model performance increases with increasingly
complex models, with MLP providing the best performance.

Figure 1: Residual topography derived from a machine learning model and a
plate model. a) MLP model of residual topography, overlain with point
measurements of oceanic residual topography hoggard2016global and hotspots.
Their linear correlation is shown in Fig. 2a. Purple contours show depth-
integrated shear-wave anomalies within the two LLSVPs. b) Low-pass filtered
MLP and c) GDH1 residual topography to degree 50. Blue and red contours show
larger and smaller LAB depth relative to its global average values.

The corresponding residual topography from the best-performing ML model (MLP)
and plate model (GDH1) are shown in Fig. 1. The raw output of the MLP residual
topography (Fig. 1a) contains multi-scale features, with peak topographic
values reaching  $\ignorespaces\sim$ 1.5 km. Comparison with the point
measurements of residual topography from Hoggard et al. hoggard2016global
displays a near 1:1 correlation (Fig. 2a). This confirms the fidelity of our
ML approach in accurately identifying residual topography. The fact that our
result covers all ocean basins makes it a more objective representation of
global residual topography, as interpolation of point data may introduce
arbitrary information to long-wavelength components yang2017oceanic. The
observation that nearly all hotspots reside in regions of positive residual
topography (Fig. 1), implying thinner lithosphere or mantle upwelling, further
verifies our result.

This residual topography filtered down to spherical harmonic degree 50 (Fig.
1b) retains all the information of interest for this study while removing
local noise. At this scale, the largest amplitude of all topographic signals
is below 1 km. This is in contrast to the plate model result which, when
filtered to the same wavelengths, displays a notably higher amplitude (
$\ignorespaces\sim$ 1.5 km) (Fig. 1c). A more striking difference comes from
their respective spectrum and correlation with the depth variation of the
lithosphere-asthenosphere boundary (LAB), a topic to be further elaborated
later (Fig. 2b,c). The LAB depth variation is calculated independent of any a
priori plate model, where we first construct a global reference LAB depth
model by averaging its observed depths for every age, and then map out the
local LAB deviation from this reference profile.

##  2 Key controlling factors on the results

In order to determine which features and layers contribute most to model
performance, we perform an ablation study. For each feature or layer, we train
an MLP model with that feature/layer removed from the dataset to see how it
impacts performance. Table 2 lists the performance metrics for each ablation.

Table 2: Performance metrics for ablation study.

(a) Feature ablation. Features | RMSE (km) | R2
---|---|---
No thickness | 0.581 | 0.760
No $v_{p}$ | 0.507 | 0.817
No $v_{s}$ | 0.489 | 0.830
No $\rho$ | 0.496 | 0.825
No age | 0.600 | 0.744
All features | 0.476 | 0.839

(b) Layer ablation. Layers | RMSE (km) | R2
---|---|---
No water | 0.483 | 0.834
No ice | 0.495 | 0.825
No sediments | 0.625 | 0.722
No crust | 0.542 | 0.791
No mantle | 0.574 | 0.766
All layers | 0.476 | 0.839

Among all structural features, age turns out to be the most important, with
the greatest impact on performance when removed. This may not be surprising
given the obvious control of age on bathymetry. One could also expect crustal
thickness and density to be important in affecting topography. Indeed,
thickness is the second most influential feature. However, density is the
second least important. In addition, P- and S-wave velocities contribute a
similarly small amount to model performance. The three least important
features may correspond to the general lack of lateral crustal variations
below oceanic basins.

The layer ablation test reveals that sediments are most impactful to model
performance when removed. In fitting seafloor bathymetry, traditional plate
models usually ignore this contribution by removing sediments through
isostatic unloading. We find that sediments serve as an important factor for
ML models to predict bathymetry. In theory, sediments and age are related,
with thicker sediment above older seafloors, as could explain their similarly
dominant role in model performance. However, they also possess independent
information from each other, as can be seen from their spatial mismatch and
non-linear correlation (Extended Data Fig. 3). This may explain why removing
them separately both notably decreases model performance.

Among other layers, crust that contains multiple features (crustal thickness,
seismic velocities, and density) is obviously important, but the mantle layer
with information on seismic velocity and density also plays a significant
role, and is worth further exploration. The water and ice layers containing
density and velocity have little impact on model performance. It should be
noted that ablation studies do not control for codependent variables. For
example, if crustal features are largely a function of age, we can explain why
removing the crust from the model does not have as large of an impact on model
performance as it would if we had also removed age.

We further test the effect of several published seafloor age datasets. In
order to determine how robust our models are against small perturbations in
the data, we report the performance of all models trained and evaluated on all
these seafloor age datasets in Extended Data Table 1. We find that different
seafloor age compilations have minor impacts on model performance.

##  3 Dichotomic residual topography

Topographic effects from the convective mantle usually correspond to
relatively long wavelengths, which we search for using low-pass filtered
residual topography down to degree 30. We observe that the spectra of the ML
model results are similar to that of the traditional GDH1 residual topography
at degrees higher than 5 (Fig. 2d), but they contrast sharply at longer
wavelengths where they become nearly anti-correlated. This reflects their
different capabilities in recognizing long-wavelength topographic features. We
suggest that the plate model, assuming age as the only control variable, fails
to capture other long-wavelength topographic signals with a shallow origin
such as sediment distribution, seismic velocity, and layer thickness.

When evaluated against the tomography-derived dynamic topography (Methods),
both the ML and GDH1 spectra show considerably less power at the longest
wavelengths (< degree 4) (Fig. 2d). However, except at degree 2, the pattern
of the power spectrum from the ML results correlates much better with that of
dynamic topography than with that of the GDH1 model, an observation lasting up
to degree 10 (Fig. 2d). While this echoes an outstanding issue that the plate-
model-based residual topography correlates poorly with dynamic topography at
the longest wavelengths davies2019earth, it also demonstrates that ML
algorithms may better isolate the topographic contribution from the convective
mantle.

Figure 2: Properties of MLP residual topography. a) Comparison of MLP residual
topography with point measurements of oceanic residual topography
hoggard2016global shows a linear relationship. b) Correlation of MLP and c)
GDH1 residual topography and LAB thickness variation, and their respective
spectra expressed in the square root of power with units of km. d) Spectrum of
residual (dynamic) topography and LAB thickness variation, all represented in
the square root of power with units of km. e) Degree-1 and f) degree-2
components of MLP residual topography.

The most striking difference between the MLP prediction and previous estimates
of residual topography (e.g., GDH1) is in their spectra, where the topographic
power at very long wavelengths (less than degree 5) is notably less than
previously thought (Fig. 2b–d). Among this wavelength window, the degree-1
component of the ML results stands out most prominently. In map view, this
corresponds to a pattern consisting of hemispheric uplift centered on the
Atlantic and a hemispheric subsidence centered on the Pacific (Fig. 2e). This
topographic component has an amplitude of  $\ignorespaces\sim$ 150 m.

In stark contrast, the commonly proposed degree-2 residual topography, mostly
based on modeled dynamic topography (Fig. 2d), has negligible power in the ML
spectra, especially the MLP prediction. This corresponds to a small amplitude
of < 100 m in map view (Fig. 2f). More surprisingly, the spatial pattern of
this degree-2 residual topography is nearly anti-correlated with that of the
LLSVPs, with subsidence occurring in the southern Atlantic and west-central
Pacific.

Further comparison of these MLP results with the spectra of LAB thickness
variation reveals positive correlation at degree-1 (Fig. 2b), due to the fact
that hemispheric uplift occurs within thick oceanic lithospheres around the
Atlantic-African region and subsidence within the Pacific (Fig. 1b, 2e). This
observation is further confirmed with the longitudinal distribution of
residual topography (Fig. 3a) and LAB depth variation (Fig. 3b), both of which
demonstrate dominant degree-1 power. This result directly contradicts
lithospheric isostasy and suggests that the hemispheric residual topography
originates from the convective mantle. We note that the same inference cannot
be concluded based on the GDH1 result (Fig. 2c). In contrast, the degree-2
patterns for both MLP and GDH1 results demonstrate a negative correlation with
LAB thickness variation (Fig. 2b,c), permitting but not requiring a possible
origin of lithospheric isostasy for this topographic component. Both these
findings are different from traditional views.

Figure 3: Longitudinal distribution of residual topography and LAB thickness.
a) MLP residual topography and b) LAB depth variation averaged over latitude
and within longitudinal windows ranging from 30° to 180°. The x-axis shows the
center longitude of each window. These stacked profiles progressively approach
a degree–1 sinusoidal curve (black dashed line). The purple dashed line is the
degree–1 curve for the estimated mantle dynamic topography (Fig. 4c). Major
mantle slabs and LLSVPs are marked.

Another peak within the spectrum of MLP residual topography occurs at degrees
5–14 (Fig. 2b). Geographically, these signals primarily occur around the
LLSVPs and above subducted slabs (Fig. 1a; 3), with both prominent uplift and
subsidence. From the longitudinal distribution, the uplift signals (relative
to the degree-1 background profile) within this wavelength range
(corresponding to the 30°–90° averaging windows) also display a hemispheric
contrast, with greater amounts of intermediate-wavelength uplift above the
African LLSVP than above the Pacific LLSVP (Fig. 3a). This further supports
the different upwelling vigor of the two LLSVPs. The subsidence pattern,
however, lacks such a contrast, with clear negative topography above major
subduction zones. The slightly less pronounced subsidence above the Farallon
slab than above the Pacific slab is consistent with their respective
subduction histories, with the former largely finished by now and its slabs
sinking mostly into the lower mantle.

##  4 Implications on mantle dynamics

Further understanding the origin of the MLP residual topography and associated
mantle dynamics requires a quantitative estimate on the respective
contributions from the lithosphere and the convective mantle. Here we derive
the former from the observed LAB depth variation of oceanic lithosphere via
thermal isostasy. Removing the lithospheric topographic contribution from the
MLP residual topography generates a plausible dynamic topography due to mantle
convection.

The deviation of LAB depth from its global reference profile varies by up to
$\ignorespaces\sim$ 150 km, with overall thicker lithosphere in the Atlantic
hemisphere than in the Pacific hemisphere (Fig. 4a). This observation is
confirmed from the longitudinal distribution of LAB depth variation, where a
prominent degree-1 pattern exists (Fig. 3b). We interpret this hemispheric
difference in lithospheric thickness as reflecting the corresponding upper
mantle temperatures, with the Atlantic side cooler than the Pacific side
dalton2014geophysical; bao2022relative. We also observe that the oceanic
lithosphere is thicker than the degree-1 baseline above former slabs and
thinner above LLSVPs (Fig. 3b), as is intuitive given the implied temperature
effect on lithospheric thickness.

Assuming a depth-averaged temperature perturbation of  $\ignorespaces\sim$ 150
°C near the base of the lithosphere turcotte2002geodynamics, a 150 km thick
layer within the lowermost lithosphere generates  $\ignorespaces\sim$ 0.7 km
of topographic relief. The pattern of the corresponding isostatic topography
variation resembles that of the LAB depth variation, with subsidence above the
thick Atlantic lithosphere and uplift above the thin Pacific lithosphere (Fig.
4b).

The resulting mantle dynamic topography, computed by removing the above
isostatic topography from the MLP residual topography, displays an even more
striking degree-1 pattern, as the hemispheric topography contrast is further
enhanced (Fig. 4c), given the positive correlation between the degree-1 mode
of MLP residual topography and LAB thickness variation (Fig. 2b). The peak-to-
trough difference in local dynamic topography reaches > 2 km, and that of its
degree-1 component is nearly twice that of the degree-1 MLP residual
topography. The overall pattern of dynamic topography is not affected by the
scaling from LAB depth variation to isostatic topography (Extended Data Fig.
2).

Figure 4: Shallow vs. deep origins of residual topography. a) LAB depth
variation relative to its global reference profile. b) Estimated contribution
to residual topography from the oceanic mantle lithosphere due to its
thickness variation. c) Estimated mantle dynamic topography by subtracting b
from the original MLP residual topography in Fig. 1b.

The robust pattern of degree-1 dynamic topography reveals that the Pacific
mantle generates less uplift than the Atlantic mantle does. We suggest that
this reflects that the African LLSVP is more dynamically unstable than the
Pacific LLSVP, as may be related to the lower intrinsic density of the former
yuan2022instability. The observation that dynamic uplift at both the
hemispheric scale and intermediate wavelengths (degrees 5–14) within the
Atlantic Ocean supports this inference. The more prominent uplift above the
African hemisphere is also consistent with the greater amount of depth-
integrated slow seismic anomalies within the African LLSVP relative to the
Pacific LLSVP (Fig. 1a).

This finding, however, is at odds with the reduced upper mantle temperature
dalton2014geophysical; bao2022relative and thicker oceanic lithosphere (Fig.
4a) below most of the Atlantic and Indian Oceans in comparison with the
Pacific Ocean. There may be two possible explanations for this. One
possibility links to the temporal evolution of the lower mantle structures,
where the African LLSVP is rapidly elevating in response to the converging
motion of the surrounding slab piles, many of which (Nazca, Farallon, Pacific,
etc.) reached the lowermost mantle during the Cenozoic peng2022quantifying.
This tends to generate more long-wavelength surface uplift than the Pacific
LLSVP that is farther away from lower mantle slabs. The latter, due to its
longer stability, could have accumulated more localized hot materials in the
upper mantle to erode the overlying lithosphere, resulting in thinner
lithosphere and milder dynamic uplift.

Another possibility concerns the hemispheric compositional difference of the
deep mantle. It has been recently proposed that the distinct geochemical
signature of the Pacific mantle from the African mantle suggests different
influences of the circum-Pacific subduction doucet2020distinct. The African
mantle is not only circumvented by multiple slab piles, but also filled with
delaminated continental lithosphere beneath much of the Atlantic Ocean
hu2018modification. Geodynamic modeling suggests that these continental keels,
mostly distributed at upper-mantle depth as seismically observed
hu2018modification, are buoyant, thus favoring mantle upwelling and surface
uplift peng2022fate. These compositional properties could explain the
prominent dynamic uplift and thick oceanic lithosphere in the Atlantic
hemisphere, as is largely absent in the Pacific side.

In summary, our new residual topography inference based on machine learning
reveals that the African mantle generates more dynamic uplift than the Pacific
mantle. We attribute this hemispheric dichotomy to the contrasting dynamic
states of two LLSVPs, where the African LLSVP is still shrinking in response
to surrounding subducted slabs, thus generating more uplift than the Pacific
side. The residence of cold but buoyant continental lithospheric materials
within the African upper mantle prevents lithospheric thinning while
maintaining its high topography. Consequently, this study highlights the
importance of lateral variation in the thermal-chemical configuration of the
deep mantle.

##  5 Methods

Below, we detail the steps taken to produce our ML results.

###  5.1 Data

The primary dataset used in this work is the CRUST1.0 model laske2013update.
This dataset provides measurements of 4 _features_ (boundary topography,
P-wave velocity, S-wave velocity, and density) across 9 _layers_ (water; ice;
upper, middle, and lower sediments; upper, middle, and lower crystalline
crust; and below the Moho). Observations are recorded on a 1° grid for the
entire Earth, with surface topography derived from ETOPO1 amante2009etopo1. We
also make use of the crustal type add-on during preprocessing, which divides
the Earth’s crust into 36 categories.

A number of seafloor age datasets produced by Müller et al. seton2020global;
muller2019global; muller2016ocean; muller2013seawater; muller2008age are used
in this study. The results reported above are produced using the most recent
dataset. In order to ensure that our results are robust across all datasets,
we also report results for older datasets in Extended Data Table 1.

A plate boundary model produced by Bird bird2003updated—specifically, a
version converted into an easy-to-use shapefile ahlenius2014world—is used
during cross-validation. We divide the Earth into 7 major tectonic plates:
Africa, Antarctica, Australia, Eurasia, North America, South America, and
Pacific. The remaining plates are merged into the smallest neighboring plate
in order to achieve similarly sized units.

###  5.2 Preprocessing

All datasets are downsampled to the same 1° resolution as CRUST1.0 and
spatially joined using GeoPandas kelsey_jordahl_2020_3946761. Using the
CRUST1.0 crustal type add-on, we filter out all data except for normal oceanic
crust (A1) and young oceanic crust (A0), including all continental crust and
oceanic plateaus. The depth to the boundary between lower sediments and upper
crystalline crust is chosen as the ground truth that our models attempt to
predict. All boundary topographies are converted to layer thicknesses, and the
thickness of the water layer is removed from the dataset since this is
approximately what we are trying to predict. Finally, all data is standardized
by subtracting the mean and dividing by the standard deviation. This
standardization is reverted before calculating any performance metrics.

###  5.3 Models

For baseline comparison, we implement the HS turcotte2002geodynamics, PSM
parsons1977analysis, GDH1 stein1992model, and H13 hasterok2013heat models
using the hyperparameters reported by the authors. We apply a sediment
correction factor, as outlined in the Supplementary Information. We experiment
with 3 different machine learning models. All machine learning models are
implemented in Scikit-learn scikit-learn.

The simplest model we try is linear regression (LR) kenney1962linear via
ordinary least squares. However, our robustness study reveals that outliers in
the seafloor age dataset result in wildly different predictions. We instead
focus on ridge regression (RR) tikhonov1977solutions, a variant of LR that
adds a regularization term to avoid overfitting.

Support vector regression (SVR) drucker1996support, a variant of support
vector machines (SVM) cortes1995support designed for regression, is chosen for
a more advanced model. Support vector models use a similar loss function as
linear regression, but ignore points within a small margin around the decision
boundary or prediction. They also allow non-linear solutions.

The most advanced model we experiment with is a multi-layer perceptron (MLP)
haykin1994neural, a type of neural network. MLPs are a series of matrix
multiplications interspersed with non-linear activation functions. They allow
the model to fit arbitrary non-linear functions, but require more training
data than other models.

###  5.4 Training and evaluation

Our training and evaluation procedure uses grouped cross-validation. All
preprocessed data is divided into 7 major tectonic plates as described above.
For each tectonic plate, a model is trained from scratch on the other 6
plates, and then used to make predictions on the 7th plate. This process is
repeated for all 7 plates until a prediction is made for every data point on
Earth. This plate tectonic-based grouped cross-validation procedure is
designed to avoid geospatial dependence between neighboring points seen in
prior work stewart2020understanding.

Once predictions have been made for all data points, the predictions are
compared to the ground truth from the CRUST1.0 dataset. We compute and report
root mean square error (RMSE) and coefficient of determination (R2).

###  5.5 Dynamic topography

We define dynamic topography as the topographic contribution from the
convective mantle beneath the lithosphere. We estimate the dynamic topography
by solving for the instantaneous global mantle flow using the spherical
finite-element code CitcomS zhong2008benchmark. The mantle density structure
is converted from a shear wave tomography S40RTS ritsema2011s40rts using an
empirical seismic-to-density conversion following a recent study
wang2022topography. The resulting density anomalies are further converted to
effective temperature perturbations, whose dynamic properties are solved for
in a thermomechanical model.

In practice, we assume an incompressible mantle that satisfies the Boussinesq
approximation. The 3D structure of the effective mantle viscosity is dependent
both on depth and on temperature wang2022topography. The convective mantle is
defined as that below the LAB, estimated as the base of the high shear-wave
heterogeneity. In practice, we take 300 km as the LAB depth below continents
and 200 km below oceans for the dynamic topography model, similar to previous
studies davies2019earth; wang2022topography. The surface dynamic topography is
obtained from the sub-lithospheric normal stress.

Data availability

The CRUST1.0 model used in our experiments is available at
https://igppweb.ucsd.edu/~gabi/crust1.html. All seafloor age datasets can be
found on EarthByte: https://www.earthbyte.org/category/resources/data-
models/seafloor-age/. The plate tectonic boundaries shapefile can be found on
GitHub: https://github.com/fraxen/tectonicplates. All predictions made by our
ML models are released under a CC-BY 4.0 license and can be found on Hugging
Face: https://huggingface.co/datasets/torchgeo/bathymetry.

Code availability

In order to ensure the reproducibility of the experiments presented in this
paper, all code used is released under an MIT license and can be found on
GitHub: https://github.com/adamjstewart/bathymetry.

Acknowledgments

The authors acknowledge the Texas Advanced Computing Center (TACC) at The
University of Texas at Austin for providing HPC resources that have
contributed to the research results reported within this paper. URL:
https://www.tacc.utexas.edu

This work made use of the Illinois Campus Cluster, a computing resource that
is operated by the Illinois Campus Cluster Program (ICCP) in conjunction with
the National Center for Supercomputing Applications (NCSA) and which is
supported by funds from the University of Illinois at Urbana-Champaign.

## References

  *   * (1) Parsons, B., Sclater, J.G.: An analysis of the variation of ocean floor bathymetry and heat flow with age. Journal of Geophysical Research 82(5), 803–827 (1977)
  * (2) Stein, C.A., Stein, S.: A model for the global variation in oceanic depth and heat flow with lithospheric age. Nature 359(6391), 123–129 (1992)
  * (3) Hasterok, D.: A heat flow based cooling model for tectonic plates. Earth and Planetary Science Letters 361, 34–43 (2013)
  * (4) Hager, B.H., Clayton, R.W., Richards, M.A., Comer, R.P., Dziewonski, A.M.: Lower mantle heterogeneity, dynamic topography and the geoid. Nature 313(6003), 541–545 (1985)
  * (5) Yang, T., Gurnis, M.: Dynamic topography, gravity and the role of lateral viscosity variations from inversion of global mantle flow. Geophysical Supplements to the Monthly Notices of the Royal Astronomical Society 207(2), 1186–1202 (2016)
  * (6) Hoggard, M.J., Czarnota, K., Richards, F.D., Huston, D.L., Jaques, A.L., Ghelichkhan, S.: Global distribution of sediment-hosted metals controlled by craton edge stability. Nature Geoscience 13(7), 504–510 (2020)
  * (7) Yang, T., Moresi, L., Müller, R.D., Gurnis, M.: Oceanic residual topography agrees with mantle flow predictions at long wavelengths. Geophysical Research Letters 44(21), 10–896 (2017)
  * (8) Steinberger, B.: Topography caused by mantle density variations: Observation-based estimates and models derived from tomography and lithosphere thickness. Geophysical Supplements to the Monthly Notices of the Royal Astronomical Society 205(1), 604–621 (2016)
  * (9) Davies, D.R., Valentine, A.P., Kramer, S.C., Rawlinson, N., Hoggard, M.J., Eakin, C.M., Wilson, C.R.: Earth’s multi-scale topographic response to global mantle flow. Nature Geoscience 12(10), 845–850 (2019)
  * (10) Hoggard, M.J., White, N., Al-Attar, D.: Global dynamic topography observations reveal limited influence of large-scale mantle flow. Nature Geoscience 9(6), 456–463 (2016)
  * (11) Wang, Y., Liu, L., Zhou, Q.: Topography and gravity reveal denser cratonic lithospheric mantle than previously thought. Geophysical Research Letters 49(1), 2021–096844 (2022)
  * (12) Wang, Y., Liu, L., Zhou, Q.: Geoid reveals the density structure of cratonic lithosphere. Journal of Geophysical Research: Solid Earth 127(8), 2022–024270 (2022)
  * (13) Amante, C., Eakins, B.W.: ETOPO1 arc-minute global relief model: Procedures, data sources and analysis (2009)
  * (14) Seton, M., Müller, R.D., Zahirovic, S., Williams, S., Wright, N.M., Cannon, J., Whittaker, J.M., Matthews, K.J., McGirr, R.: A global data set of present-day oceanic crustal age and seafloor spreading parameters. Geochemistry, Geophysics, Geosystems 21(10), 2020–009214 (2020)
  * (15) Müller, R.D., Zahirovic, S., Williams, S.E., Cannon, J., Seton, M., Bower, D.J., Tetley, M.G., Heine, C., Le Breton, E., Liu, S., et al.: A global plate model including lithospheric deformation along major rifts and orogens since the Triassic. Tectonics 38(6), 1884–1907 (2019)
  * (16) Müller, R.D., Seton, M., Zahirovic, S., Williams, S.E., Matthews, K.J., Wright, N.M., Shephard, G.E., Maloney, K.T., Barnett-Moore, N., Hosseinpour, M., et al.: Ocean basin evolution and global-scale plate reorganization events since Pangea breakup. Annual Review of Earth and Planetary Sciences 44(1), 107–138 (2016)
  * (17) Müller, R.D., Dutkiewicz, A., Seton, M., Gaina, C.: Seawater chemistry driven by supercontinent assembly, breakup, and dispersal. Geology 41(8), 907–910 (2013)
  * (18) Müller, R.D., Sdrolias, M., Gaina, C., Roest, W.R.: Age, spreading rates, and spreading asymmetry of the world’s ocean crust. Geochemistry, Geophysics, Geosystems 9(4) (2008)
  * (19) Laske, G., Masters, G., Ma, Z., Pasyanos, M.: Update on CRUST1.0—–a 1-degree global model of Earth’s crust. In: Geophysical Research Abstracts, vol. 15, p. 2658 (2013). EGU General Assembly Vienna, Austria
  * (20) Tikhonov, A.N., Arsenin, V.J., Arsenin, V.I., Arsenin, V.Y., et al.: Solutions of Ill-posed Problems. Vh Winston, Washington, D.C. (1977)
  * (21) Drucker, H., Burges, C.J., Kaufman, L., Smola, A., Vapnik, V.: Support vector regression machines. Advances in Neural Information Processing Systems 9 (1996)
  * (22) Haykin, S.: Neural Networks: a Comprehensive Foundation. Prentice Hall PTR, Hoboken, NJ (1994)
  * (23) Turcotte, D.L., Schubert, G.: Geodynamics. Cambridge University Press, Cambridge, UK (2002)
  * (24) Dalton, C.A., Langmuir, C.H., Gale, A.: Geophysical and geochemical evidence for deep temperature variations beneath mid-ocean ridges. Science 344(6179), 80–83 (2014)
  * (25) Bao, X., Lithgow-Bertelloni, C.R., Jackson, M.G., Romanowicz, B.: On the relative temperatures of earth’s volcanic hotspots and mid-ocean ridges. Science 375(6576), 57–61 (2022)
  * (26) Yuan, Q., Li, M.: Instability of the African large low-shear-wave-velocity province due to its low intrinsic density. Nature Geoscience 15(4), 334–339 (2022)
  * (27) Peng, D., Liu, L.: Quantifying slab sinking rates using global geodynamic models with data-assimilation. Earth-Science Reviews 230, 104039 (2022)
  * (28) Doucet, L.S., Li, Z.-X., Gamal El Dien, H., Pourteau, A., Murphy, J.B., Collins, W.J., Mattielli, N., Olierook, H.K.H., Spencer, C.J., Mitchell, R.N.: Distinct formation history for deep-mantle domains reflected in geochemical differences. Nature Geoscience 13(7), 511–515 (2020)
  * (29) Hu, J., Liu, L., Faccenda, M., Zhou, Q., Fischer, K.M., Marshak, S., Lundstrom, C.: Modification of the Western Gondwana craton by plume-lithosphere interaction. Nature Geoscience 11(3), 203–210 (2018)
  * (30) Peng, L., Liu, L., Liu, L.: The fate of delaminated cratonic lithosphere. Earth and Planetary Science Letters 594, 117740 (2022)
  * (31) Bird, P.: An updated digital model of plate boundaries. Geochemistry, Geophysics, Geosystems 4(3) (2003)
  * (32) Ahlenius, H.: World tectonic plates and boundaries (2014). https://github.com/fraxen/tectonicplates
  * (33) Jordahl, K., den Bossche, J.V., Fleischmann, M., Wasserman, J., McBride, J., Gerard, J., Tratner, J., Perry, M., Badaracco, A.G., Farmer, C., Hjelle, G.A., Snow, A.D., Cochran, M., Gillies, S., Culbertson, L., Bartos, M., Eubank, N., maxalbert, Bilogur, A., Rey, S., Ren, C., Arribas-Bel, D., Wasser, L., Wolf, L.J., Journois, M., Wilson, J., Greenhall, A., Holdgraf, C., Filipe, Leblanc, F.: GeoPandas. Zenodo. https://doi.org/10.5281/zenodo.3946761
  * (34) Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., Duchesnay, E.: Scikit-learn: Machine learning in Python. Journal of Machine Learning Research 12, 2825–2830 (2011)
  * (35) Kenney, J.F., Keeping, E.S.: Linear regression and correlation. Mathematics of Statistics 1, 252–285 (1962)
  * (36) Cortes, C., Vapnik, V.: Support-vector networks. Machine Learning 20(3), 273–297 (1995)
  * (37) Stewart, A., Liu, L., Peng, J.: Understanding ocean bathymetry using machine learning. In: AGU Fall Meeting Abstracts, vol. 2020, pp. 013–08 (2020)
  * (38) Zhong, S., McNamara, A., Tan, E., Moresi, L., Gurnis, M.: A benchmark study on mantle convection in a 3-d spherical shell using CitcomS. Geochemistry, Geophysics, Geosystems 9(10) (2008)
  * (39) Ritsema, J., Deuss, A., Van Heijst, H.J., Woodhouse, J.H.: S40RTS: A degree-40 shear-velocity model for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements. Geophysical Journal International 184(3), 1223–1236 (2011)

##  6 Extended data

Extended Data Fig. 1: Plate boundaries used in our study. Extended Data Fig.
2: Different scaling of isostatic topography contribution from oceanic mantle
lithosphere due to thickness variation (left) and the corresponding mantle
dynamic topography (right). a) Half and b) twice the peak-to-trough amplitude
of lithospheric residual topography as in Fig. 4b (0.9 km). The resulting
mantle dynamic topography on the right for both cases retains a similar
degree-1 dominant pattern as in Fig. 4c.

(a) Sediment thickness (km).

(b) Seafloor age (Ma).

Extended Data Fig. 3: Comparison between sediment thickness and seafloor age.
Although these features have a high degree of correlation in the Atlantic and
Indian Oceans, they have little correlation in the Western Pacific and Arctic
Oceans. Extended Data Fig. 4: Diagram of isostatic balance equation. Layers
include water, sediment, crust, and mantle. Dashed line represents depth of
compensation. Extended Data Table 1: Performance metrics for robustness study.

(a) HS model. Year | RMSE (km) | R2
---|---|---
2020 | 1.051 | 0.214
2019 | 1.030 | 0.231
2016 | 1.131 | 0.284
2013 | 1.154 | 0.259
2008 | 1.166 | 0.255

(b) PSM model. Year | RMSE (km) | R2
---|---|---
2020 | 0.689 | 0.662
2019 | 0.678 | 0.667
2016 | 0.729 | 0.703
2013 | 0.731 | 0.702
2008 | 0.739 | 0.701

(c) GDH1 model. Year | RMSE (km) | R2
---|---|---
2020 | 0.634 | 0.714
2019 | 0.627 | 0.715
2016 | 0.675 | 0.745
2013 | 0.676 | 0.746
2008 | 0.682 | 0.745

(d) H13 model. Year | RMSE (km) | R2
---|---|---
2020 | 0.652 | 0.697
2019 | 0.645 | 0.699
2016 | 0.692 | 0.732
2013 | 0.693 | 0.733
2008 | 0.699 | 0.733

(e) RR model. Year | RMSE (km) | R2
---|---|---
2020 | 0.592 | 0.751
2019 | 0.590 | 0.748
2016 | 0.609 | 0.792
2013 | 0.613 | 0.791
2008 | 0.613 | 0.794

(f) SVR model. Year | RMSE (km) | R2
---|---|---
2020 | 0.564 | 0.773
2019 | 0.560 | 0.773
2016 | 0.592 | 0.804
2013 | 0.601 | 0.799
2008 | 0.557 | 0.830

(g) MLP model. Year | RMSE (km) | R2
---|---|---
2020 | 0.476 | 0.839
2019 | 0.478 | 0.834
2016 | 0.513 | 0.853
2013 | 0.497 | 0.863
2008 | 0.462 | 0.883

Extended Data Table 2: Hyperparameter tuning experiments.

(a) Ridge Regression (RR). Parameter | Parameter Sweep | Optimal Parameter
---|---|---
$\alpha$ | $\left\\{0,10^{-3},10^{-2},10^{-1},10^{0},10^{1},10^{2},10^{3}\right\\}$ | $10^{2}$

(b) Support Vector Regression (SVR). Parameter | Parameter Sweep | Optimal Parameter
---|---|---
kernel | {linear, poly, rbf, sigmoid} | rbf
$\gamma$ | {scale, auto} | scale
$c_{0}$ | $\left\\{0,10^{-4},10^{-3},10^{-2},10^{-1}\right\\}$ | $10^{-4}$
$C$ | $\left\\{10^{-1},10^{0},10^{1},10^{2},10^{3}\right\\}$ | $10^{0}$
$\epsilon$ | $\left\\{0,10^{-2},10^{-1},10^{0},10^{1}\right\\}$ | $10^{-1}$

(c) Multi-Layer Perceptron (MLP). Parameter | Parameter Sweep | Optimal Parameter
---|---|---
activation | {logistic, tanh, relu} | tanh
solver | {lbfgs, sgd, adam} | adam
$\alpha$ | $\left\\{10^{-3},10^{-2},10^{-1},10^{0},10^{1}\right\\}$ | $10^{-1}$
hidden layers | {3, 4, 5, 6, 7} | 6
hidden size | $\left\\{2^{7},2^{8},2^{9},2^{10},2^{11}\right\\}$ | $2^{8}$
learning rate | $\left\\{10^{-6},10^{-5},10^{-4},10^{-3},10^{-2}\right\\}$ | $10^{-3}$

##  7 Supplementary Information

All plate models include a sediment correction factor that accounts for the
load of sediments on oceanic crust. This sediment correction factor can be
calculated as follows.

Extended Data Fig. 4 illustrates the isostatic balance equation we are trying
to solve. Column $t_{1}$ represents the current state, while column $t_{2}$
represents a fictional state where sediments have been removed. We set the
depth of compensation to the dashed line at the crust-mantle boundary in
column $t_{1}$. Following Airy isostasy, we assume that all densities remain
constant and that crustal thickness does not change during this correction.

Our goal is to solve for the uplift of the oceanic crust, $h_{m}$, after
removing the sediments. We do this using Archimedes’ Principle, which states
that the weight (and therefore mass) of these two columns must be equivalent
in order to achieve isostatic balance. We use the fact that oceanic crust
cancels out and that $h_{1}+h_{s}=h_{2}+h_{m}$.

$\displaystyle\begin{split}\rho_{w}h_{1}+\rho_{s}h_{s}+\cancel{\rho_{c}h_{c}}&=\rho_{w}h_{2}+\cancel{\rho_{c}h_{c}}+\rho_{m}h_{m}\\\
\rho_{w}h_{1}+\rho_{s}h_{s}&=\rho_{w}h_{2}+\rho_{m}h_{m}\\\
\rho_{w}h_{1}+\rho_{s}h_{s}&=\rho_{w}(h_{1}+h_{s}-h_{m})+\rho_{m}h_{m}\\\
\cancel{\rho_{w}h_{1}}+\rho_{s}h_{s}&=\cancel{\rho_{w}h_{1}}+\rho_{w}h_{s}-\rho_{w}h_{m}+\rho_{m}h_{m}\\\
\rho_{m}h_{m}-\rho_{w}h_{m}&=\rho_{s}h_{s}-\rho_{w}h_{s}\\\
h_{m}(\rho_{m}-\rho_{w})&=\rho_{s}h_{s}-\rho_{w}h_{s}\\\
h_{m}&=\frac{\rho_{s}h_{s}-\rho_{w}h_{s}}{\rho_{m}-\rho_{w}}\end{split}$ (1)

In our case, sediments are divided into upper, middle, and lower layers. Thus,
the equation becomes

$\displaystyle
h_{m}=\frac{\rho_{us}h_{us}+\rho_{ms}h_{ms}+\rho_{ls}h_{ls}-\rho_{w}h_{s}}{\rho_{m}-\rho_{w}},$
(2)

where $h_{s}=h_{us}+h_{ms}+h_{ls}$.