# Towards Universal Mesh Movement Networks

Mingrui Zhang1 Chunyang Wang1 Stephan Kramer1 Joseph G. Wallwork2
Siyi Li1 Jiancheng Liu4 Xiang Chen3 Matthew D. Piggott1
1Imperial College London 2University of Cambridge 3Noah’s Ark Lab, Huawei
4Michigan State University
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Solving complex Partial Differential Equations (PDEs) accurately and
efficiently is an essential and challenging problem in all scientific and
engineering disciplines. Mesh movement methods provide the capability to
improve the accuracy of the numerical solution without increasing the overall
mesh degree of freedom count. Conventional sophisticated mesh movement methods
are extremely expensive and struggle to handle scenarios with complex boundary
geometries. However, existing learning-based methods require re-training from
scratch given a different PDE type or boundary geometry, which limits their
applicability, and also often suffer from robustness issues in the form of
inverted elements. In this paper, we introduce the Universal Mesh Movement
Network (UM2N), which – once trained – can be applied in a non-intrusive,
zero-shot manner to move meshes with different size distributions and
structures, for solvers applicable to different PDE types and boundary
geometries. UM2N consists of a Graph Transformer (GT) encoder for extracting
features and a Graph Attention Network (GAT) based decoder for moving the
mesh. We evaluate our method on advection and Navier-Stokes based examples, as
well as a real-world tsunami simulation case. Our method outperforms existing
learning-based mesh movement methods in terms of the benchmarks described
above. In comparison to the conventional sophisticated Monge-Ampère PDE-solver
based method, our approach not only significantly accelerates mesh movement,
but also proves effective in scenarios where the conventional method fails.
Our project page is at https://erizmr.github.io/UM2N/.

##  1 Introduction

Various natural phenomena are modeled by Partial Differential Equations (PDE).
The accurate and efficient approximation of the solutions to these complex and
often nonlinear equations represents a fundamental challenge across all
scientific and engineering disciplines. To solve real-world PDEs, many
numerical methods require a computational mesh or grid to discretize the
spatial domain. The quality of this mesh significantly impacts the balance
between a numerical solution’s accuracy and the computational cost required to
obtain it. To maintain high-resolution everywhere in the domain is
computationally expensive. Many systems modeled by PDEs however are multi-
resolution in nature. For example, a small part of the system may be highly
dynamic, while other regions are quasi-stationary; alternatively, some
locations may be more important while others less important to the question
being considered by the calculation [1, 2, 3, 4]. While a uniform, high-
resolution mesh can be utilized to capture the dynamics or what is most
important accurately, this is often wasteful of finite, precious computational
resources and comes with sustainability implications [5].

An active body of research is focused on developing deep learning based
methods to accelerate the PDE-solving process by learning surrogate neural
solvers [6, 7, 8]. A learned solver can output PDE solutions directly given
the physical parameters, and boundary conditions of a PDE. However, these
methods encounter challenges such as not being able to guarantee physical
plausibility of the PDE solution, weak generalization ability, and low data
efficiency.

An alternative approach to improve the efficiency of a PDE solver is to
utilize mesh adaptation, which is a technique for distributing mesh according
to the requirements of numerical accuracy. Two main categories of mesh
adaptation techniques can be identified: $h$-adaptation and $r$-adaptation.
$h$-adaptation refines or coarsens the mesh resolution dynamically through
topological operations such as adding/deleting nodes and swapping element
edges/faces. In contrast, $r$-adaptation (or _mesh movement_) relocates or
moves mesh nodes, keeping the mesh topology and thus the number of elements
and vertices in the mesh unchanged [9]. These traditional mesh adaptation
techniques can help reduce PDE-solving costs, but the mesh adaptation process
itself may come at the cost of significant computational overhead.

Deep learning based methods have been proposed to accelerate mesh adaptation.
[10, 11, 12, 13, 14, 15]. Most previous methods focus on $h$-adaptive
refinement [11, 12, 13, 16]. Some existing works focus on learning error
indicators or Riemannian metrics (which account for element shape and
anisotropy, as well as size) [13, 14, 15] to guide the mesh adaptation
process. The learned indicator or metric is then fed into a traditional
remesher for mesh generation, or a mesh optimizer, which overall limits the
performance of these works to being no better than traditional methods,
especially in terms of efficiency. Reinforcement learning based methods show
potential to improve the mesh adaptation task, but are difficult to train with
low data efficiency [11, 16, 17]. Only a small number of works focus on
$r$-adaptation. [18] and [19] investigated supervised and unsupervised
learning based mesh movement methods. However, the proposed methods require
re-training from scratch given a different PDE or geometry, which limits their
applicability. In addition, the unsupervised learning method requires
significant GPU memory and suffers from long training times, which makes re-
training even more prohibitive for large scale problems.

Targeting the generalization ability and efficiency of mesh movement, here we
introduce the Universal Mesh Movement Network (UM2N), a two stage, deep-
learning based model that learns to move the computational mesh guided by an
optimal-transport approach driven by solutions of equations of Monge-Ampère
(MA) type. Given an underlying PDE to solve numerically, the process of MA
based mesh movement includes choosing a suitable monitor function (i.e., a
measure for the desired mesh density) and solving an auxiliary PDE whose
solution prescribes the movement of the mesh. The Monge-Ampère equation is an
example of such an auxiliary PDE, which results in mesh movement with
attractive theoretical properties [9], but whose solution comes at a
significant computational cost. Therefore, we decouple the underlying PDE
solve and the mesh movement process. The learning of the mesh movement process
is essentially learning to solve the auxiliary Monge-Ampère equation.

The proposed UM2N consists of a Graph Transformer based encoder and a Graph
Attention Network (GAT) [20] based decoder.

Element volume loss is selected within the training loss, instead of the
coordinates loss used in previous works [18, 19]. This offers the advantage
that element volume loss not only provides supervision signals but also
penalizes negative volumes (i.e., inverted elements), thereby reducing mesh
tangling and enhancing robustness. We construct a PDE-independent dataset by
generating random generic fields for model training. The trained model can be
applied in a zero-shot manner to move meshes with different sizes and
structure, for solvers applicable to different PDE types and boundary
geometries. We comprehensively demonstrate the effectiveness of our methods on
different cases including advection, Navier-Stokes examples, as well as a
real-world tsunami simulations on meshes with different sizes and structures.
Our UM2N outperforms existing learning-based mesh movement methods in terms of
the benchmarks described above. In comparison to the sophisticated Monge-
Ampère PDE-solver based method, our approach not only significantly
accelerates mesh movement, but also proves effective in scenarios where the
conventional method fails. The trained model can be directly integrated into
any mesh based numerical solvers for error reduction and acceleration in a
non-intrusive way, which can benefit various engineering applications
suffering trade-offs between accuracy and computational cost.

Our main contributions are listed as follows:

$\bullet$ We present the UM2N framework, which once trained can be applied in
a non-intrusive, zero-shot manner to move meshes with different sizes and
structures, for solvers applicable to different PDE types and boundary
geometries.

$\bullet$ We propose PDE-independent training strategies and mesh tangling
aware loss functions for our Graph Transformer and GAT-based architecture.

$\bullet$ We demonstrate the universal ability of UM2N on different PDEs and
boundary geometries: Advection, Navier-Stokes examples, and a real-world
tsunami simulation case.

##  2 Related Works

Neural PDE solver. Neural PDE surrogate models can significantly accelerate
the PDE solving process. There are three main groups of existing works.
Physics-informed Neural Networks (PINNs) [6, 21, 22] are the first group,
modeling the PDEs through implicit representations by neural networks. PINNs
require knowledge of the governing equations and train neural networks with
equation residuals from the PDEs with given boundary conditions and initial
conditions. Neural operators are another group of methods for learning to
solve PDEs which seek to learn a mapping from a function describing the
problem to the corresponding solution function [7, 8, 23, 24, 25]. The third
group is mesh-based PDE solvers with deep learning, including the utilization
of CNNs and GNNs for processing structured and unstructured meshes [26, 27,
15, 28, 29, 30, 31]. Our method can be considered to be related to the neural
operator approach, although it learns a mapping from monitor function values
to the potential field describing the mesh movement instead of the PDE
solution directly. In addition, the networks in the proposed model share
similar ideas to mesh-based PDE solvers for processing unstructured meshes.

Learning for mesh generation and adaptation. Deep learning methods have been
advanced for applications in mesh generation and adaptation. Two tasks need to
be differentiated here, the first involves modeling the ideal mesh resolution
or density, e.g. through learning a metric or error indicator. The second
involves the step that actually adapts the mesh in response, where we can
distinguish between $h$-adaptivity, in which the connectivity structure of the
mesh is changed to refine or coarsen the mesh, and $r$-adaptivity which
maintains the connectivity and only moves the vertices.

An example of the first task is MeshingNet [10] which employs a neural network
to establish the requisite local mesh density, subsequently utilized by a
standard Delaunay triangulation-based mesh generator. Similarly, the
determination of optimal local mesh density for the purpose of mesh refinement
is described in [12]. In [13], the authors propose a learning-based model to
determine optimal anisotropy to guide $h$-adaptation. [15] introduces a
neural-network-based approach for predicting the sizing field for use with an
external remesher. Instead of learning a metric, some works learn to
manipulate the mesh directly. [11, 17, 16] formulate the $h$-adaptation
process as a reinforcement learning paradigm, learning a policy that guides to
adjust the mesh elements.

$r$-adaptation (or mesh movement) is generally formulated in a fundamentally
different manner to $h$-adaptation. Leading $r$-adaptation methods are based
on a mapping between a fixed reference computational domain and the physical
domain encompassing the adapted mesh, which is established through the
solution of a nonlinear auxiliary PDE. [18] first proposed the use of a GNN
based mesh movement network. [19] investigated the unsupervised learning of
mesh movement for improving performance of neural PDEs solvers. [32]
introduced a learning based mesh movement method specially tailored for
computation fluid dynamics simulation of airfoil design.

Our work focuses specifically on $r$-adaptation, i.e., mesh movement methods.
The approaches proposed in previous work require re-training the model when
applied to different problems, or can only generalize within a limited range
of PDE parameters. We propose a mesh movement network that aims to be
universal across applications to different PDEs without retraining.

##  3 Preliminaries and Problem Statement

Monge-Ampère mesh movement. A mesh movement problem can be defined as the
search for the transformation between a fixed computational domain
$\Omega_{C}$, and a physical domain $\Omega_{P}$. Coordinates in these spaces
are denoted $\boldsymbol{\xi}$ and $\mathbf{x}$, respectively. The mesh
$\mathcal{H}_{C}$ defined in the computational domain is the original, often
uniform mesh. The mesh $\mathcal{H}_{P}$ defined in the physical domain is the
adapted (or moved) mesh.111In what follows, we abuse notation slightly by also
using $\boldsymbol{\xi}$ to refer the the (discrete) vertices of
$\mathcal{H}_{C}$, rather than the (continuous) coordinate field. Similarly,
we also use $\boldsymbol{x}$ to denote the vertices of $\mathcal{H}_{P}$. The
goal of mesh movement is to find a mapping
$\boldsymbol{x}=f(\boldsymbol{\xi})$ between the continuous coordinate fields
of $\Omega_{C}$ and $\Omega_{P}$, although in practice we deduce a mapping of
the discrete vertex sets of $\mathcal{H}_{C}$ and $\mathcal{H}_{P}$. While the
coordinate transformation can also be used to accommodate time-dependent
moving boundaries [33] or the optimisation of the shape of the domain, here we
focus on its application to achieve variable resolution such that the solution
$u_{\boldsymbol{x}}$ of a PDE solved on the adapted mesh has higher accuracy
than the solution $u_{\boldsymbol{\xi}}$ solved on the original mesh.

A monitor function $m$ over the spatial domain is used to specify the desired
mesh density, which can be based on various characteristics of the solution
field $u$ of the PDE. It prescribes where the mesh should be refined or
coarsened, _i.e._ large monitor values indicate where high mesh resolution is
required. Therefore, the goal of the mesh movement process can be rephrased as
finding a mapping so that $m$ is equidistributed over the adapted (_i.e._ ,
physical) mesh [34]:

$m({\boldsymbol{x}})\>\det(\boldsymbol{J})=\theta,$ (1)

where $\boldsymbol{J}$ is the Jacobian of the map
$\boldsymbol{x}=f(\boldsymbol{\xi})$ with respect to the computational
coordinates $\boldsymbol{\xi}$, $\det(\boldsymbol{J})$ denotes the determinant
of the Jacobian which corresponds to the relative change in volume under the
transformation, and $\theta$ is a normalization constant. By using concepts
from optimal transport theory [35], the problem can be constrained to have a
unique solution, with the deformation of the map expressed in terms of the
gradient of a scalar potential $\phi$ such that:

$\boldsymbol{x}(\boldsymbol{\xi})=\boldsymbol{\xi}+\nabla_{\boldsymbol{\xi}}\phi(\boldsymbol{\xi}).$
(2)

Substituting the additional constraint described by equation (2) into equation
(1) gives a nonlinear PDE of Monge-Ampère type:

$m({\boldsymbol{x}})\>\det(\boldsymbol{I}+\boldsymbol{H}(\phi(\boldsymbol{\xi})))=\theta,$
(3)

where $\boldsymbol{I}$ is the identity matrix and $\boldsymbol{H}(\phi)$ is
the Hessian of $\phi$, with derivatives taken with respect to
$\boldsymbol{\xi}$, i.e,
$\boldsymbol{H}(\phi)_{ij}=\frac{\partial^{2}\phi}{\partial\xi_{i}\partial\xi_{j}}$.
It should be noted that guarantees for existence and uniqueness, and
convergence of solution methods for the Monge-Ampère equation generally rely
on the domain being diffeomorphic with a convex domain (e.g. [36]). In
practice, the method is known to break down in some cases where the domain is
not simply-connected or has non-smooth boundaries.

The MA equation (3) is an auxiliary PDE which is purely associated with the
the mesh movement process, and is independent of the underlying PDE or
physical problem we wish to solve, i.e., the MA equation keeps the same form
for different PDEs. This decoupling benefits learning-based methods in terms
of generalization properties, i.e., a well-learned MA neural solver has the
potential to be applied to different problems with a small cost of fine-tuning
or even without re-training.

##  4 UM2N: Universal Mesh Movement Network

###  4.1 Framework overview.

Figure 1: Overview of Universal Mesh Movement Network.

The proposed UM2N framework is shown in Fig. 1. Given an input mesh, vertex
features and edge features are collected separately. The coordinates and
monitor function values are gathered from the vertices and input into a graph
transformer to extract embeddings. The embedding vector $\boldsymbol{z}$
obtained from the graph transformer encoder is then concatenated with the
extracted edge features $\boldsymbol{e}$ to serve as the input for the Graph
Attention Network (GAT) decoder. The decoder processes this combined input
along with a mesh query to ultimately produce an adapted mesh.

Graph Transformer encoder. Transformers [37] are a popular architecture which
have been successfully applied for neural PDE solvers [38, 24, 25]. They have
strong expressivity and can naturally handle irregular data structures such as
unstructured meshes. Their attention mechanism can capture both the local and
global information for a vertex of interest in a mesh. Here we use a graph
transformer encoder as a feature extractor. The input features include the
coordinates of the original mesh $\boldsymbol{\xi}$ and monitor values
$m_{\boldsymbol{\xi}}$. The coordinates serve as the positional encoding.
Features are concatenated and encoded into $\boldsymbol{Q}$ (query),
$\boldsymbol{K}$ (key) and $\boldsymbol{V}$ (value) as inputs using MLPs for a
self-attention module. The encoder output embedding $\boldsymbol{z}$ are fed
into a downstream GAT based deformer for mesh movement.

Graph Attention Network (GAT) based decoder. We choose GATs to construct our
decoder. Its attention mechanism can help constrain the mesh vertex movement
to within one hop of its neighbors which assists in alleviation of mesh
tangling issues. The decoder consists of $N$ GAT blocks. The first GAT block
takes a mesh query $\boldsymbol{\xi}$, and the embedding $\boldsymbol{z}$ from
the transformer encoder as well as the edge features $\boldsymbol{e}$ as
inputs. For each following $k$th block, the inputs consist of the coordinates
of the initial mesh $\boldsymbol{\xi}$, the coordinates of the intermediate
moved mesh $\boldsymbol{\xi}^{(k-1)}$ and extracted features from the previous
layer $\boldsymbol{f}^{(k-1)}$.

###  4.2 PDE-independent dataset and training.

In the existing work of neural mesh movement [18, 19], proposed models are
trained on solutions of PDEs, with the training data generated by solving one
specific type of PDE, limiting generalizability of the trained model. Aiming
here to train universal mesh movement networks, we construct a PDE-independent
training dataset
$\mathcal{D}=\\{\boldsymbol{d}=(\boldsymbol{\xi},m_{\boldsymbol{\xi}};\boldsymbol{V}_{r},\boldsymbol{T}_{r})\\}$,
where $\boldsymbol{\xi},m_{\boldsymbol{\xi}}$ denotes the original mesh and
monitor values as the model’s input, and
$\boldsymbol{V}_{r},\boldsymbol{T}_{r}$ as the pre-calculated (reference)
mesh’s vertex coordinates and elements as ground truth. To build the dataset,
we randomly generate generic solution fields, which are composed as the
summation of a random number of Gaussian distribution functions centred in
random locations, and with random widths in different directions to introduce
anisotropy, samples are shown in Figure A1; refer to Appx. B for a detailed
description. The generated fields can be interpreted as representing solutions
of any PDE. For mesh movement, a widely used method to translate the PDE
solution $u(\bf x)$ into a monitor function, is the Hessian based formula
[35]:

$m({\bf x})=1+\alpha\frac{\|H(u)\|}{\max\|H(u)\|},$ (4)

where $\|H(u)\|$ indicates the Frobenius norm of the Hessian of $u$ and
$\alpha$ is a user chosen constant. This choice results in small cells where
the curvature of the solution is high, with a ratio of $1+\alpha$ between the
largest and smallest cell volumes. Our network is trained on the values of the
monitor function only, which can be computed relatively cheaply from the PDE
solution. The expected movement of the mesh vertices is computed by solving
the Monge-Ampère equation (3) for each of the monitor functions based on the
randomly generated PDE solutions. After training, the model can then be
applied using only the monitor values as input, independent of the PDE that is
being solved. Other choices of formulae for the monitor function in terms of
$u$ may be appropriate in different cases, but again these can be applied
without re-training the model, e.g., the flow-past-a-cylinder and tsunami test
cases in the following section use a monitor function based on the gradient of
the solution.

###  4.3 Loss functions.

Given the dataset $\mathcal{D}$ defined above, the final objective is to find
model parameters $\boldsymbol{\theta}$, by minimizing the total loss
$L_{total}$. Thus, we can formalize the process as

$\operatorname*{arg\,min}_{\boldsymbol{\theta}}L_{total}(\boldsymbol{\theta}):=\lambda_{vol}L_{vol}(\boldsymbol{\theta})+\lambda_{cd}L_{cd}(\boldsymbol{\theta}).$
(5)

where $L_{vol},L_{cd}$ represents the element volume loss and Chamfer distance
loss respectively, which are defined later. $\lambda_{cd},\lambda_{vol}>0$
represent hyper-parameters balancing these two effects.

We denote the modified (adapted) mesh with model parameters
$\boldsymbol{\theta}$ as
$\mathcal{M}(\boldsymbol{\xi},m_{\boldsymbol{\xi}};\boldsymbol{\theta})=\\{\boldsymbol{V}_{a},\boldsymbol{T}_{a}\\}$,
where $\boldsymbol{V}_{a}=\\{\boldsymbol{x}_{i}|i=1,2,\ldots,n_{1}\\}$
represents the coordinates of vertices and
$\boldsymbol{T}_{a}=\\{\boldsymbol{t}_{i}|i=1,2,\ldots,n_{2}\\}$ represents
the elements with $\boldsymbol{x}_{i},\boldsymbol{t}_{i}$ as the $i$-th vertex
and element. For ease of presentation, we omit the $\boldsymbol{\theta}$ here.
Similarly, we define $\boldsymbol{y}_{i},\boldsymbol{q}_{i}$ as
$\boldsymbol{V}_{r},\boldsymbol{T}_{r}$’s $i$-th vertex and element.

Element volume loss. In contrast to the coordinate loss used in previous work
[18, 19], here we use the element volume loss for training. Element volume
loss is computed as the averaged volume difference between each element in the
adapted mesh ($\boldsymbol{T}_{a}$) and the reference mesh
($\boldsymbol{T}_{r}$) computed from the MA method. It is defined as

$L_{vol}(\boldsymbol{\theta})=\mathbb{E}_{(\boldsymbol{\xi},m_{\boldsymbol{\xi}};\boldsymbol{V}_{r},\boldsymbol{T}_{r})\in\mathcal{D}}\left[\frac{1}{|\boldsymbol{T}_{a}|}\sum_{i=1}^{n_{2}}|\text{Vol}(\boldsymbol{t}_{i})-\text{Vol}(\boldsymbol{q}_{i})|\right],$
(6)

where Vol is the function computing the volume of a given element
$\boldsymbol{t}$.

The element volume loss is inspired by the equidistribution relation in
equation (1). Intuitively, the equidistribution relation enforces deformed
mesh elements such that their volumes correspond to rescaling by the monitor
function value $m({\boldsymbol{x}})$. Therefore, the element volume loss is
utilised to encourage the model to learn to move the mesh so as to conform to
the equidistribution relation, under the guidance of MA method, given that the
MA method provides one good way to achieve this relation. In addition, the
element volume loss penalizes negative Jacobian determinants, which helps
prevent element inversion, i.e., mesh tangling.

Chamfer distance loss. Chamfer distance finds for all vertices in the first
mesh, the nearest vertex in a second mesh and sums the square of these
distances. It encourages the model to output a mesh which has a similar
spatial distribution of vertices to that of the reference mesh. To achieve
vertex-to-vertex alignment, the bidirectional Chamfer distance is utilized
without additional sampling. It is defined as

$L_{cd}(\boldsymbol{\theta})=\mathbb{E}_{(\boldsymbol{\xi},m_{\boldsymbol{\xi}};\boldsymbol{V}_{r},\boldsymbol{T}_{r})\in\mathcal{D}}\left[\frac{1}{|\boldsymbol{V}_{a}|}\sum_{\boldsymbol{x}_{i}\in\boldsymbol{V}_{a}}\min_{\boldsymbol{y}_{j}\in\boldsymbol{V}_{r}}\|\boldsymbol{x}_{i}-\boldsymbol{y}_{j}\|_{\ell_{2}}+\frac{1}{|\boldsymbol{V}_{r}|}\sum_{y_{j}\in\boldsymbol{V}_{r}}\min_{\boldsymbol{x}_{i}\in\boldsymbol{V}_{a}}\|\boldsymbol{x}_{i}-\boldsymbol{y}_{j}\|_{\ell_{2}}\right].$
(7)

##  5 Experiments

###  5.1 Experiment setups.

Unlike existing approaches to mesh movement, such as [18, 19], which re-train
their model on a case-by-case basis, our study aims to explore the universal
applicability of the proposed UM2N. Therefore, all examples shown in this
section are tested in a zero-shot generalization manner, i.e., without re-
training. In addition to PDE type, our training meshes are of small size
typically comprising $500$ vertices, while our test meshes have far more
vertices: Advection ($2,052$ vertices), Cylinder ($4,993$ vertices), Tsunami
($8,117$ vertices), which also tests the model’s scalability. There are
examples with more complicated settings, please refer to Appx. D.

Training. The training dataset consists of 600 randomly generated generic
solution fields and original meshes with 463 and 513 vertices. The model is
trained using the Adam optimizer. The training and experiments are performed
on an Nvidia RTX 3090 GPU.

Metrics. The main metric for evaluating mesh movement quality is underlying
PDE solution error reduction (ER) ratio. Here PDE error is approximated by the
difference between the solutions obtained on a coarse mesh and accurate
solutions which are obtained on a very high resolution mesh. The PDE error
reduction ratio is the difference between the PDE error from an original mesh
and that from an adapted mesh. Another important metric for mesh movement is
mesh tangling. Once mesh tangling happens, the mesh is invalid for a numerical
solver, which breaks the simulation.

###  5.2 Benchmarking mesh movement.

Table 1: Quantitative results across mesh movement methods. ER indicates PDE
error reduction ratio. “Fail” indicates that mesh tangling happens during the
simulation.

Methods | Swirl | Cylinder | Helmholtz
---|---|---|---
ER (%) $\uparrow$ | Time (ms) $\downarrow$ | ER $\uparrow$ (%) | Time (ms) $\downarrow$ | ER (%) $\uparrow$ | Time (ms) $\downarrow$
MA [35] | 37.43 | 5615.78 | Fail | - | 17.95 | 3265.23
M2N [18] | 14.15 | 25.01 | Fail | - | 16.32 | 4.62
UM2N (Ours) | 35.93 | 17.66 | 61.29 | 19.25 | 17.08 | 4.86

11footnotetext: All time measures are in milliseconds.

In the following section, we benchmark the mesh movement using both non-
learned Monge-Ampère (MA) mesh movement implemented in _Movement_ [39] and
learning-based baseline M2N [18] over three different scenarios to compare
with UM2N. Table 1 shows the quantitative results in three different
scenarios, Swirl, Cylinder, and Helmholz. An additional qualitative result of
Tōhoku Tsunami simulation, is provided to show our model can handle highly
complex boundaries in real-world scenarios. Moreover, to give a fair
comparison, we re-trained the model in M2N on our training dataset before
evaluation.

Helmholtz. The Helmholtz case is a time-independent problem solving a 2D
Helmholtz equation as shown in (A2). The source terms of the equation are
generated use formula as shown in (A1). Our model achieves better error
reduction than M2N and comparable to MA method with significant less time
cost.

Swirl. The swirl case considers a time-dependent pure advection equation. In
this scenario, the initial tracer field takes a narrow ring shape and is
advected within a swirling velocity field, see Eq. (A3) and (A4) in Appx. C.
The velocity swirling direction is anti-clockwise in the first half of the
simulation and clockwise for the second half. Given this setting, as the
advection process is reversible, the final state should be consistent with the
initial state. This constitutes a problem of significant multi-resolution
complexity, as the extremely narrow ring requires a high resolution to
accurately resolve the dynamics at the sharply defined interfaces of both
inner and outer sides.

The results are shown in Fig. 2. As shown in Table 1, our UM2N shows a average
$35.93\%$ error reduction which clearly outperforms M2N’s $14.15\%$. It is
close to that of the full Monge-Ampère (MA) mesh movement approach which
achieves an error reduction of $37.43\%$, but with significantly lower
computational time.

Figure 2: Results for swirl case. The plot shows the end state of the swirl
simulation (bottom row) and corresponding meshes (top row). As the process is
reversible, the end state should be consistent to the initial state. It can be
observed that the solution obtained on high resolution almost recover the
initial state. The proposed UM2N outperform M2N and achieve comparable
performance to MA method. Quantitatively in the rightmost plot, the orange
line (UM2N) and blue line (MA) best suppress error accumulation along
timesteps.

Flow past cylinder. Simulation of flow past a cylinder in a channel is a
classic and challenging multi-resolution case of computational fluid dynamics,
solving the Navier-Stokes equations. The cylinder and the channel’s top-bottom
surfaces require a high resolution mesh to resolve the boundary layers [40].
Given appropriate Reynolds number, the well-known von Kármán vortex street
phenomena can be observed within the wake flow which also benefits from higher
resolution to resolve the time-varying evolution of vorticity structures. The
drag $C_{D}$ and lift $C_{L}$ coefficients for the cylinder are other
quantities of interest in this problem. Here we use a benchmark configuration
from [41], please refer to Appx.C.3 for more details.

The simulation is run for $8,000$ time steps of size $0.001$ s, i.e., $8$ s
physical time for the whole simulation. The original mesh is a uniform mesh
with triangular elements and $4,993$ vertices. The reference solutions are
obtained on a high resolution mesh with $260,610$ vertices.

Figure 3: Results for the flow past cylinder case in time. In the upper part,
we visualize $5$ snapshots of adapted meshes and vorticity intensity from
$4.00s$ to $4.20s$. Please refer to the full video in the supplementary
materials. It shows that our UM2N can output meshes which adapt to the
dynamics around the cylinder. It can be observed from the zoom-in mini-plots
at the top-left that UM2N reduces errors/noise in vorticity intensity compared
to the original mesh. In the lower part, the blue, orange and green lines show
the drag $C_{D}$ coefficients obtained on a high resolution fixed mesh, UM2N
adapted mesh and the original coarse mesh for the first 6 seconds (i.e., 6000
steps). It can be observed that the UM2N output mesh improves the accuracy of
$C_{D}$ in magnitude and periodicity compared to the original mesh.  Figure 4:
Analysis on wake flow vorticity for the cylinder case. The adapted mesh output
by UM2N (see top-left part of the plot). Qualitatively, the error maps at the
top-right part indicate that our UM2N reduce the vorticity difference to the
high resolution mesh result in general. The quantitative comparison results
along the probes are shown in the lower part of the figure. There are four
cross-section probes place on $x=[0.5,1.0,1.5,2.0]$ marked by the black dash
line in the middle left plot. The green line (our UM2N) is much more
consistent to the blue dash line (high resolution) compared to the orange line
(original mesh), which further verifies that our method can reduce the PDE
errors.

We investigate the dynamics in the area around the cylinder and in the wake
flow respectively. For the former, we use the measure of $C_{D}$ as shown in
Fig. 3. The results for $C_{L}$ is shown in Fig. A5 in Appx. D. For the wake
flow, we use cross-section probes at $x=[0.5,1.0,1.5,2.0]$ to measure the
intensity of vorticity spatially along the $y$-axis as shown in Fig. 4.

Qualitatively, it can be observed that the proposed method moves the mesh in a
manner that captures the fluid dynamics for the whole domain in general (see
upper part of Fig. 3). At the cylinder and top-bottom surface, the vertices
are moved to increase the resolution in order to resolve the boundary layers
as expected. Quantitatively, the drag $C_{D}$ and lift $C_{L}$ coefficients
computed on our moved mesh show a $61.29\%$ average error reduction compared
to that of original mesh as shown in lower part of Fig. 3. The proposed UM2N
approach not only improves the accuracy of these coefficients in magnitude but
also prevents their periodic variation from phase shifting. For the wake flow,
as shown in Fig. 4, we selected a single snapshot at $t=4.0s$ for qualitative
investigation. Our UM2N also shows an accuracy improvement for vorticity
intensity on data sampled from all four cross-section probes comparing to the
that on original mesh.

Note that we intended to compare our UM2N to M2N and the conventional full MA
method. However, both MA and M2N fail i.e., the solver diverges, within 30
time steps due to mesh tangling and thus provide insufficient data points for
comparison.

Figure 5: Qualitative results for the Tōhoku tsunami simulation. The boundary
of the simulation is generated based on the real coastline data (see the mesh
overlaid on the satellite map). We show 9 snapshots of the wave elevation from
80 steps of the tsunami simulation enhanced by our UM2N. In these
visualizations, red and blue hues indicate wave elevations that are,
respectively, above or below mean sea level. Regions exhibiting significant
elevation magnitudes necessitate increased resolution to accurately resolve
the underlying dynamics. Our method dynamically and robustly adjusts the mesh
to increase resolution at the wave front, thereby effectively tracking its
propagation. Figure 6: Handling complex boundaries. The Monge-Ampère method
generates inverted elements, resulting in mesh tangling near the coastline as
shown in the areas delineated by red rectangles. Our UM2N can smoothly handle
the elements near the coastline without mesh tangling.

Tōhoku tsunami simulation. Here we present the case of a simulation of the
2011 Tōhoku Tsunami, to show that our methods can be applied to real world
scenarios with boundaries describing highly complex geometries. The tsunami is
simulated by solving the shallow water equations using the Thetis framework
[42]. For details on the datasets and software used to set up the Tōhoku case
study, see Appx. B.2. Note that the mesh ($8,117$ vertices) has an arbitrarily
irregular boundary. It can be observed that the mesh moves in a manner that
helps track the wave propagation with enhanced resolution, as shown in Fig. 5.
Our model can also robustly handle the highly complex boundary without mesh
tangling near the coastline. As a comparison, the conventional Monge-Ampère
(MA) mesh movement approach struggles to handle the complex geometries, i.e.,
fails to converge when solving the non-linear MA equation, resulting in highly
tangled elements near the coastline as shown in Fig. 6. Therefore, our method
not only provides a more efficient approach but also show advantages dealing
with complex boundaries compared to the full MA mesh movement approach. Please
see the full video of this case in the supplementary materials.

###  5.3 Ablation study.

Volume loss vs coordinate loss. We denote the UM2N variant trained with
coordinate loss [18] as UM2N-coord. As shown in Table 2, UM2N-coord achieves
comparable error reduction ($31.21\%$) to the proposed UM2N ($35.93\%$) on the
Swirl case. However, it fails (i.e., encounters mesh tangling) on the Cylinder
case which has a more complex boundary geometry. This indicates that the
element volume loss is mesh-tangling aware, i.e., helps to prevent producing
inverted elements.

Table 2: Ablation Study for design choices of learning monitors and element
volume loss.

| Helmholtz | Swirl | Cylinder
---|---|---|---
Method | ER (%) $\uparrow$ | ER $\uparrow$ (%) | ER (%) $\uparrow$
UM2N-coord | 16.64 | 31.21 | Fail
UM2N-sol | 13.52 | -0.35 | 39.94
UM2N (Ours) | 17.08 | 35.93 | 61.29

Monitors vs PDE solutions as inputs. We denote the UM2N variant using the PDE
solution as input, instead of monitor function values, as UM2N-sol. The
UM2N-sol is trained on solutions of the Helmholtz PDE. The UM2N-sol shows a
much inferior performance ($-0.35\%$) compared to UM2N ($35.93\%$) on the
Swirl case as well as on the Cylinder case. This indicates that learning mesh
movement from monitors directly improves the generalisability of the model.

##  6 Conclusions

We introduce the Universal Mesh Movement Network (UM2N), which once trained on
our generated PDE-independent dataset, can be applied in a non-intrusive,
zero-shot manner to move meshes with different sizes and structure, for
solvers applicable to different PDE types and boundary geometries. Our method
demonstrates superior mesh movement results on Advection and Navier-Stokes PDE
examples as well as a real-world tsunami simulation case. Our UM2N outperforms
existing learning based method and achieves comparable PDE error reduction to
the costly Monge-Ampére (MA) PDE based approach. Our UM2N demonstrates
effectiveness on cases with complex boundary geometries where existing
learning and conventional MA based methods fail. Our method also shows
significant acceleration compared to MA based mesh adaptation. See Appx. F for
the discussions on limitations and broader impacts.

## References

  * Wallwork et al. [2020] Joseph G. Wallwork, Nicolas Barral, Stephan C. Kramer, David A. Ham, and Matthew D. Piggott.  Goal-oriented error estimation and mesh adaptation for shallow water modelling.  _SN Applied Sciences_ , 2(6):1053, 2020.  ISSN 2523-3971.  doi: 10.1007/s42452-020-2745-9.
  * Vlachas et al. [2022] Pantelis R Vlachas, Georgios Arampatzis, Caroline Uhler, and Petros Koumoutsakos.  Multiscale simulations of complex systems by learning their effective dynamics.  _Nature Machine Intelligence_ , 4(4):359–366, 2022.
  * Zhang et al. [2021] Chi Zhang, Massoud Rezavand, and Xiangyu Hu.  A multi-resolution sph method for fluid-structure interactions.  _Journal of Computational Physics_ , 429:110028, 2021.  ISSN 0021-9991.  doi: https://doi.org/10.1016/j.jcp.2020.110028.
  * Athanasopoulos et al. [2009] Michael Athanasopoulos, Hassan Ugail, and Gabriela González Castro.  Parametric design of aircraft geometry using partial differential equations.  _Advances in Engineering Software_ , 40(7):479–486, 2009.
  * CF2 [2023] The carbon footprint of computational research.  _Nature Computational Science_ , 3(8):659–659, 2023.
  * Raissi et al. [2019] Maziar Raissi, Paris Perdikaris, and George E Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  _Journal of Computational Physics_ , 378:686–707, 2019.
  * Lu et al. [2019] Lu Lu, Pengzhan Jin, and George Em Karniadakis.  Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators.  _arXiv preprint arXiv:1910.03193_ , 2019.
  * Li et al. [2021] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Fourier neural operator for parametric partial differential equations.  _International Conference on Learning Representations (ICLR)_ , 2021.
  * Huang and Russell [2011] Weizhang Huang and Robert D. Russell.  _Adaptive Moving Mesh Methods_.  Applied Mathematical Sciences. Springer Science+Business Media, LLC, 2 edition, 2011.
  * Zhang et al. [2020] Zheyan Zhang, Yongxing Wang, Peter K. Jimack, and He Wang.  MeshingNet: A new mesh generation method based on deep learning, 2020.
  * Yang et al. [2021] Jiachen Yang, Tarik Dzanic, Brenden K. Petersen, Jun Kudo, Ketan Mittal, Vladimir Z. Tomov, Jean-Sylvain Camier, Tuo Zhao, Hongyuan Zha, Tzanio V. Kolev, Robert W. Anderson, and Daniel Faissol.  Reinforcement learning for adaptive mesh refinement.  _CoRR_ , abs/2103.01342, 2021.  URL https://arxiv.org/abs/2103.01342.
  * Huang et al. [2021] Keefe Huang, Moritz Krügener, Alistair Brown, Friedrich Menhorn, Hans-Joachim Bungartz, and Dirk Hartmann.  Machine learning-based optimal mesh generation in computational fluid dynamics.  _arXiv preprint arXiv:2102.12923_ , 2021.
  * Fidkowski and Chen [2021] Krzysztof J. Fidkowski and Guodong Chen.  Metric-based, goal-oriented mesh adaptation using machine learning.  _Journal of Computational Physics_ , 426:109957, 2021.  ISSN 0021-9991.  doi: https://doi.org/10.1016/j.jcp.2020.109957.  URL https://www.sciencedirect.com/science/article/pii/S0021999120307312.
  * Tingfan et al. [2021] Wu Tingfan, Liu Xuejun, An Wei, Zenghui Huang, and Lyu Hongqiang.  A mesh optimization method using machine learning technique and variational mesh adaptation.  _Chinese Journal of Aeronautics_ , 2021.
  * Pfaff et al. [2021] Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gonzalez, and Peter W Battaglia.  Learning mesh-based simulation with graph networks.  _International Conference on Machine Learning (ICLR)_ , 2021.
  * Wu et al. [2023] Tailin Wu, Takashi Maruyama, Qingqing Zhao, Gordon Wetzstein, and Jure Leskovec.  Learning controllable adaptive simulation for multi-resolution physics.  In _The Eleventh International Conference on Learning Representations_ , 2023.  URL https://openreview.net/forum?id=PbfgkZ2HdbE.
  * Foucart et al. [2023] Corbin Foucart, Aaron Charous, and Pierre F.J. Lermusiaux.  Deep reinforcement learning for adaptive mesh refinement.  _Journal of Computational Physics_ , 491:112381, 2023.  ISSN 0021-9991.  doi: https://doi.org/10.1016/j.jcp.2023.112381.  URL https://www.sciencedirect.com/science/article/pii/S002199912300476X.
  * Song et al. [2022] Wenbin Song, Mingrui Zhang, Joseph G Wallwork, Junpeng Gao, Zheng Tian, Fanglei Sun, Matthew D Piggott, Junqing Chen, Zuoqiang Shi, Xiang Chen, and Jun Wang.  M2n: Mesh movement networks for pde solvers.  In _Advances in neural information processing systems (NeurIPS)_ , 2022.
  * Hu et al. [2024] Peiyan Hu, Yue Wang, and Zhi-Ming Ma.  Better neural PDE solvers through data-free mesh movers.  In _The Twelfth International Conference on Learning Representations_ , 2024.
  * Veličković et al. [2017] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio.  Graph attention networks.  _arXiv preprint arXiv:1710.10903_ , 2017.
  * Cai et al. [2021] Shengze Cai, Zhiping Mao, Zhicheng Wang, Minglang Yin, and George Em Karniadakis.  Physics-informed neural networks (pinns) for fluid mechanics: A review.  _Acta Mechanica_ , 37(12):1727–1738, 2021.
  * Han et al. [2018] Jiequn Han, Arnulf Jentzen, and Weinan E.  Solving high-dimensional partial differential equations using deep learning.  _Proceedings of the National Academy of Sciences_ , 115(34):8505–8510, 2018.
  * Li et al. [2023a] Zijie Li, Kazem Meidani, and Amir Barati Farimani.  Transformer for partial differential equations’ operator learning.  _Transactions on Machine Learning Research_ , 2023a.  ISSN 2835-8856.  URL https://openreview.net/forum?id=EPPqt3uERT.
  * Li et al. [2023b] Zijie Li, Dule Shu, and Amir Barati Farimani.  Scalable transformer for pde surrogate modeling.  _Advances in neural information processing systems (NeurIPS)_ , 36, 2023b.
  * Hao et al. [2023] Zhongkai Hao, Chengyang Ying, Zhengyi Wang, Hang Su, Yinpeng Dong, Songming Liu, Ze Cheng, Jun Zhu, and Jian Song.  Gnot: A general neural operator transformer for operator learning.  _International Conference on Machine Learning (ICML)_ , 2023.
  * Thuerey et al. [2018] Nils Thuerey, Konstantin Weissenow, Lukas Prantl, and Xiangyu Y. Hu.  Deep learning methods for reynolds-averaged navier–stokes simulations of airfoil flows.  _AIAA Journal_ , 2018.  URL https://api.semanticscholar.org/CorpusID:67855877.
  * Sanchez-Gonzalez et al. [2020] Alvaro Sanchez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter Battaglia.  Learning to simulate complex physics with graph networks.  _International Conference on Machine Learning (ICML)_ , 2020.
  * Brandstetter et al. [2022] Johannes Brandstetter, Daniel Worrall, and Max Welling.  Message passing neural pde solvers.  _International Conference on Learning Representations (ICLR)_ , 2022.
  * JANNY et al. [2023] Steeven JANNY, Aurèlien Bènèteau, Madiha Nadri, Julie Digne, and Nicolas THOME.  Eagle: Large-scale learning of turbulent fluid dynamics with mesh transformers.  _International Conference on Learning Representations_ , 2023.
  * Prantl et al. [2022] Lukas Prantl, Benjamin Ummenhofer, Vladlen Koltun, and Nils Thuerey.  Guaranteed conservation of momentum for learning particle-based fluid dynamics.  In _Conference on Neural Information Processing Systems_ , 2022.
  * Han et al. [2022] Xu Han, Han Gao, Tobias Pfaff, Jian-Xun Wang, and Li-Ping Liu.  Predicting physics in mesh-reduced space with temporal attention.  _International Conference on Learning Representations_ , 2022.
  * Yu et al. [2024] Jian Yu, Hongqiang Lyu, Ran Xu, Wenxuan Ouyang, and Xuejun Liu.  Flow2Mesh: A flow-guided data-driven mesh adaptation framework.  _Physics of Fluids_ , 36(3):037124, 03 2024.
  * Nunez Rattia et al. [2018] J.M. Nunez Rattia, J.R. Percival, S.J. Neethling, and M.D. Piggott.  Modelling local scour near structures with combined mesh movement and mesh optimisation.  _Journal of Computational Physics_ , 375:1220–1237, 2018.
  * Budd et al. [2009] Chris J Budd, Weizhang Huang, and Robert D Russell.  Adaptivity with moving grids.  _Acta Numerica_ , 18:111–241, 2009.
  * McRae et al. [2018] Andrew T. T. McRae, Colin J. Cotter, and Chris J. Budd.  Optimal-transport–based mesh adaptivity on the plane and sphere using finite elements.  _SIAM Journal on Scientific Computing_ , 40(2):A1121–A1148, 2018.  doi: 10.1137/16M1109515.
  * Sulman et al. [2011] Mohamed M. Sulman, J.F. Williams, and Robert D. Russell.  An efficient approach for the numerical solution of the monge–ampère equation.  _Applied Numerical Mathematics_ , 61(3):298–307, 2011.
  * Vaswani et al. [2017] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  _Advances in neural information processing systems (NeurIPS)_ , 30, 2017.
  * Cao [2021] Shuhao Cao.  Choose a transformer: Fourier or Galerkin.  _Advances in neural information processing systems (NeurIPS)_ , 34, 2021.
  * Wallwork et al. [2024] Joseph G. Wallwork, Stephan C. Kramer, Mingrui Zhang, and Davor Dundovic.  Movement, May 2024.  URL https://doi.org/10.5281/zenodo.11239647.
  * Schlichting and Gersten [2017] Hermann Schlichting and Klaus Gersten.  _Boundary-Layer Theory_.  Springer, Berlin Heidelberg, 9 edition, 2017.  ISBN 978-3-662-52917-1.  doi: 10.1007/978-3-662-52919-5_2.
  * Schäfer et al. [1996] M. Schäfer, S. Turek, F. Durst, E. Krause, and R. Rannacher.  _Benchmark Computations of Laminar Flow Around a Cylinder_ , pages 547–566.  Vieweg+Teubner Verlag, Wiesbaden, 1996.  ISBN 978-3-322-89849-4.  doi: 10.1007/978-3-322-89849-4_39.  URL https://doi.org/10.1007/978-3-322-89849-4_39.
  * Kärnä et al. [2018] T. Kärnä, S. C. Kramer, L. Mitchell, D. A. Ham, M. D. Piggott, and A. M. Baptista.  Thetis coastal ocean model: discontinuous galerkin discretization for the three-dimensional hydrostatic equations.  _Geoscientific Model Development_ , 11(11):4359–4382, 2018.  doi: 10.5194/gmd-11-4359-2018.  URL https://gmd.copernicus.org/articles/11/4359/2018/.
  * Rathgeber et al. [2016] Florian Rathgeber, David A Ham, Lawrence Mitchell, Michael Lange, Fabio Luporini, Andrew TT McRae, Gheorghe-Teodor Bercea, Graham R Markall, and Paul HJ Kelly.  Firedrake: automating the finite element method by composing abstractions.  _ACM Transactions on Mathematical Software (TOMS)_ , 43(3):1–27, 2016.
  * Alnæs et al. [2014] Martin S Alnæs, Anders Logg, Kristian B Ølgaard, Marie E Rognes, and Garth N Wells.  Unified form language: A domain-specific language for weak formulations of partial differential equations.  _ACM Transactions on Mathematical Software (TOMS)_ , 40(2):1–37, 2014.
  * Balay et al. [2024] Satish Balay, Shrirang Abhyankar, Mark F. Adams, Steven Benson, Jed Brown, Peter Brune, Kris Buschelman, Emil Constantinescu, Lisandro Dalcin, Alp Dener, Victor Eijkhout, Jacob Faibussowitsch, William D. Gropp, Václav Hapla, Tobin Isaac, Pierre Jolivet, Dmitry Karpeev, Dinesh Kaushik, Matthew G. Knepley, Fande Kong, Scott Kruger, Dave A. May, Lois Curfman McInnes, Richard Tran Mills, Lawrence Mitchell, Todd Munson, Jose E. Roman, Karl Rupp, Patrick Sanan, Jason Sarich, Barry F. Smith, Stefano Zampini, Hong Zhang, Hong Zhang, and Junchao Zhang.  PETSc/TAO users manual.  Technical Report ANL-21/39 - Revision 3.21, Argonne National Laboratory, 2024.
  * Geuzaine and Remacle [2009] Christophe Geuzaine and Jean-François Remacle.  Gmsh: A 3-d finite element mesh generator with built-in pre- and post-processing facilities.  _International Journal for Numerical Methods in Engineering_ , 79(11):1309–1331, 2009.  doi: https://doi.org/10.1002/nme.2579.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.2579.
  * Avdis et al. [2018] Alexandros Avdis, Adam S. Candy, Jon Hill, Stephan C. Kramer, and Matthew D. Piggott.  Efficient unstructured mesh generation for marine renewable energy applications.  _Renewable Energy_ , 116:842–856, 2018.  ISSN 0960-1481.  doi: https://doi.org/10.1016/j.renene.2017.09.058.  URL https://www.sciencedirect.com/science/article/pii/S0960148117309205.
  * Wessel and Smith [1996] Pål Wessel and Walter H. F. Smith.  A global, self-consistent, hierarchical, high-resolution shoreline database.  _Journal of Geophysical Research: Solid Earth_ , 101(B4):8741–8743, 1996.  doi: https://doi.org/10.1029/96JB00104.  URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/96JB00104.
  * Amante and Eakins [2022] Christopher Amante and Barry W Eakins.  ETOPO1 arc-minute global relief model: procedures, data sources and analysis, 2022.  Accessed: 2022-03-22.
  * Okada [1985] Yoshimitsu Okada.  Surface deformation due to shear and tensile faults in a half-space.  _Bulletin of the seismological society of America_ , 75(4):1135–1154, 1985.

## Appendix

##  Appendix A Software used in the paper.

In this paper, PDEs are solved using finite element methods with _Firedrake_
[43]. Firedrake is written in Python, uses the _Unified Form Language (UFL)_
[44] domain-specific language to represent finite element forms, and
automatically generates C kernels for computational efficiency. Firedrake uses
_PETSc_ [45] for solving linear and nonlinear equations that result from
discretisation, as well as for its underlying unstructured mesh
representation.

Original meshes (i.e., meshes before adaptation) are generated by either
Firedrake or Gmsh [46]. Gmsh is a meshing software especially designed for
finite element methods, providing access to several meshing algorithms. In the
case of the Tōhoku tsunami simulation, the mesh of the Japan Sea is generated
using _qmesh_ [47] – a meshing tool specifically designed for coastal ocean
modelling applications, which manages GIS (Geoscientific Information
System)-based datasets to construct a valid input geometry used by gmsh to
generate the mesh.

The conventional mesh movement strategy used to generate training data is
implemented in _Movement_ [39], which is itself written in UFL and Firedrake.
Movement implements several mesh movement approaches, including two based on
solutions of Monge-Ampère type equations. In this paper, we use the
‘relaxation’ approach, which solves the nonlinear MA equation in an iterative
fashion by introducing a pseudo-timestep (see [35] for details).

##  Appendix B Dataset description.

###  B.1 Dataset generation pipeline.

We developed a pipeline for generating PDE-independent training data and
further research. Through our pipeline, generic fields are generated and can
be utilized as solutions of any PDE. The monitor values computed from the
generic field are used to guide the MA-based method to perform mesh adaption
on the original mesh.

We stack $N$ randomly generated Gaussian distributions to form the generic
field $u$. Concretely, they can be generated by the formula:

$u=\sum_{k=1}^{N}\exp\left(\frac{(x-\mu_{x})^{2}}{{\sigma_{x}}^{2}}+\frac{(y-\mu_{y})^{2}}{{\sigma_{y}}^{2}}\right),$
(A1)

where $\mu_{x}$ and $\mu_{y}$ denote the means of the distribution along two
orthogonal directions in a 2D domain, while $\sigma_{x}$ and $\sigma_{y}$
represent the respective standard deviations.

Figure A1: Samples of generated generic fields.

All training data are generated within a 2-D unit square domain. Equivalent
sizes $n\times n$ meshes ($n\in\\{18,20\\}$) are generated for model training.
Here we show samples of the generated fields in Fig. A1.

###  B.2 Evaluation data for the Tōhoku tsunami simulation.

As mentioned in Sec. A, the mesh of the Japan Sea is generated using qmesh
[47]. For this, we provide qmesh with coastline data from [48]. Bathymetry
(sea bed topography) data is taken from the ETOPO1 dataset [49]. A highly
idealised tsunami source condition is used, comprised of a sum of Okada
functions [50].

##  Appendix C Details of experimental settings

###  C.1 Helmholtz

The Helmholtz equation used in the paper is defined as:

$-\nabla^{2}u+u=f,\quad\nabla
u\cdot\overset{\rightarrow}{n}=0\;\;\text{on}\;\Gamma,$ (A2)

where $f$ is the source term and $\Gamma$ denotes the boundary.

###  C.2 Swirl

We use a similar setting to the unsteady 2D linear advection case in [17].

#### Initial Condition

$u_{0}=\exp\left(-\frac{1}{2\sigma^{2}}(\sqrt{(x-x_{0})^{2}+(y-y_{0})^{2}}-r_{0})^{2}\right).$
(A3)

The initial condition defines a ring centered at $(x_{0},y_{0})$ with an inner
radius $r_{0}$. The thickness of the ring is $3\sigma$. In this paper,
$r_{0}=0.2$, $(x_{0},y_{0})=(0.30,0.30)$, and $\sigma=0.05/3$.

#### Velocity Field

$\boldsymbol{c}(x,y,t)=\left(\frac{9}{10}a(t)\sin^{2}(\pi x)\sin(2\pi
y),-\frac{9}{10}a(t)\sin^{2}(\pi y)\sin(2\pi x)\right).$ (A4)

In this paper, $0<t<1$, if $t<0.5$ then $a(t)=1$, if $0.5<=t<1$, $a(t)=-1$,
i.e., the velocity direction is reversed at $t=0.5$

#### Monitor Function

$m({\bf
x})=1+\max\left(\alpha\frac{\|H(u)\|}{\max\|H(u)\|},\beta\frac{\|G(u)\|}{\max\|G(u)\|}\right)$
(A5)

The $\|H(u)\|$ and $\|G(u)\|$ are hessian norm and gradient norm respectively.
In practical, we set $\alpha=5$ and $\beta=10$.

###  C.3 Flow past a cylinder

#### Setting

The Reynolds number is 100, a quadratic velocity profile is imposed at the
left boundary as inflow condition. At the right boundary a zero pressure
outflow condition is applied. The cylinder and top-bottom surface are imposed
with non-slip boundaries.

##  Appendix D Additional experimental results.

#### Flow past multiple cylinders.

The Fig. A2 shows a challenge scenario with 5 cylinders in the domain. The
proposed UM2N can move mesh capturing the complex dynamics without mesh
tangling.

Figure A2: Multiple Cylinders

#### Flow past V-formation cylinders.

The Fig. A3 shows a bio-inspired V-formation setting. The proposed UM2N can
move mesh capturing the wake flows.

Figure A3: V-formation

#### Flow past cylinders and squares.

The Fig. A4 shows the proposed UM2N can handle another challenging scenario
with a mixture of cylinders and squares of different sizes.

Figure A4: Cylinder and Square

#### Lift coefficient ($C_{L}$) of flow past cylinder.

In Fig. A5, we visualized the lift coefficient over timestep across different
methods.

Figure A5: Results of $C_{L}$ coefficients. The blue, orange and green lines
show the drag $C_{L}$ coefficients obtained on a high resolution fixed mesh,
UM2N adapted mesh and the original coarse mesh for the first 6 seconds (i.e.,
6000 steps). It can be observed that the UM2N output mesh improves the
accuracy of $C_{L}$ in magnitude and periodicity compared to the original
mesh.

##  Appendix E Mesh movement enhanced simulation.

The Mesh Movement Enhanced Simulation Algorithm, as outlined in Alg. 1,
integrates our mesh adaption strategies into the simulation process to enhance
the quality and accuracy of computational solutions. The proposed MMES
procedure can be adapted into almost any time-dependent numerical simulation
processes non-intrusively.

Algorithm 1 Mesh Movement Enhanced Simulation

1:procedure Simulation($\mathcal{P}$, $\boldsymbol{\xi}_{\text{init}}$, $T$)

2: $\boldsymbol{\xi}_{t}\leftarrow\boldsymbol{\xi}_{\text{init}}$

3: while $t<T$ do

4: $u$ $\leftarrow$ PDESolver($\mathcal{P},\boldsymbol{\xi}_{t}$)

5: $m_{\boldsymbol{\xi}}\leftarrow$ MonitorFunc($u$, $\boldsymbol{\xi}_{t}$)

6: $\boldsymbol{\xi}_{t+1}\leftarrow\text{UM2N($\boldsymbol{\xi}_{t}$,
$m_{\boldsymbol{\xi}}$)}$

7: MeshIsValid $\leftarrow$ CheckMeshNotTangle($\boldsymbol{\xi}_{t+1}$)

8: if MeshIsValid then

9: $\boldsymbol{\xi}_{t}\leftarrow\boldsymbol{\xi}_{t+1}$

10: else

11: Break

12: end if

13: TimeStepUpdate()

14: end while

15:end procedure

For a given time-dependent PDE $\mathcal{P}$ and an initial mesh
$\boldsymbol{\xi}_{\text{init}}$, the algorithm evaluates the monitor function
$m$ based on the solution field of the PDE at each time step $t$. The monitor
value $m$ is then used by our proposed model to guide mesh movement. For
robustness, a mesh integrity check is performed before time stepping (a
tangled mesh will interrupt the procedure).

##  Appendix F Limitations and broader impacts.

Limitations. 1) The GAT to an extent constrains the mesh movement to lie in a
restricted range, which may be less effective in scenarios requiring
substantial large deformation within a single iteration of the mesh mover. 2)
There is no theoretical guarantee that mesh tangling will be completely
prevented. 3) Given a low-quality or extremely coarse original mesh, the
efficacy of the proposed method may be constrained by the inherent limitation
on the degree of freedom count, as it reallocates rather than increases the
number of mesh vertices. A promising future direction could be exploring the
integration of mesh movement with $h$-adaptation, potentially developing a
hybrid $hr$-adaptation approach that combines the advantages of both
techniques.

Broader Impacts. This research contributes to the field of AI for physical
simulation by introducing an advanced learning-based mesh movement method that
enhances the accuracy and efficiency of numerical simulations. The primary
application of this method lies in areas requiring high-fidelity simulations
of physical phenomena, such as geophysics, renewable energy, climate modeling,
and aerodynamics. This method can aid in better weather or hazards
forecasting, more efficient energy harvesting, fast prototyping of aircraft
design etc.


* [164] Iain Barclay, Alun Preece, Ian Taylor, and Dinesh Verma.  Towards traceability in data ecosystems using a bill of materials model, 2019.
  * [165] Protocol Labs.  Libp2p.  IPFS Docs, 2021.  https://docs.ipfs.tech/concepts/libp2p/.
  * [166] Md. Mainul Islam, Mpyana Mwamba Merlec, and Hoh Peter In.  A comparative analysis of proof-of-authority consensus algorithms: Aura vs clique.  In IEEE International Conference on Services Computing (SCC), pages 327–332, 2022.
  * [167] Cyril Naves Samuel, Severine Glock, François Verdier, and Patricia Guitton-Ouhamou.  Choice of ethereum clients for private blockchain: Assessment from proof of authority perspective.  In IEEE International Conference on Blockchain and Cryptocurrency (ICBC), pages 1–5, 2021.
  * [168] Jie Xu, Cong Wang, and Xiaohua Jia.  A survey of blockchain consensus protocols.  ACM Comput. Surv., 2023.  Just Accepted.
  * [169] Juhar Abdella, Zahir Tari, Adnan Anwar, Abdun Mahmood, and Fengling Han.  An architecture and performance evaluation of blockchain-based peer-to-peer energy trading.  IEEE Transactions on Smart Grid, 12(4):3364–3378, 2021.
  * [170] Golam Bashar, Graham Hill, Subroto Singha, Praneeth Marella, Gaby G. Dagher, and Jidong Xiao.  Contextualizing consensus protocols in blockchain: A short survey.  In 2019 First IEEE International Conference on Trust, Privacy and Security in Intelligent Systems and Applications (TPS-ISA), pages 190–195, 2019.


# Parameter-free version of Adaptive Gradient Methods for Strongly-Convex
Functions

Deepak Gouda∗, Hassan Naveed∗, Salil Kamath∗
Georgia Institute of Technology
Atlanta, USA
{deepakgouda, hnaveed3<EMAIL_ADDRESS>

###### Abstract

The optimal learning rate for adaptive gradient methods applied to
$\lambda$-strongly convex functions relies on the parameters $\lambda$ and
learning rate $\eta$. In this paper, we adapt a universal algorithm along the
lines of Metagrad, to get rid of this dependence on $\lambda$ and $\eta$. The
main idea is to concurrently run multiple experts and combine their
predictions to a “master" algorithm. This master enjoys $O(d\log T)$ regret
bounds.

**footnotetext: These authors contributed equally to this work

_K_ eywords Online Convex Optimization $\cdot$ Adagrad $\cdot$ Metagrad
$\cdot$ Strongly-Convex functions

##  1 Introduction

Online convex optimization (OCO) is a widely recognized approach for
representing sequential decision-making processes. The OCO algorithm operates
in the following manner:

  * •

the learner selects an action $\mathbf{x}_{t}$

  * •

an adversary chooses a loss function $f_{t}(\cdot)$

  * •

the learner incurs a loss $f_{t}(\mathbf{x}_{t})$

  * •

the learner updates its action based on this loss

The aim is to minimize regret, which is defined as the disparity between the
cumulative loss suffered by the learner and that of the optimal action in
retrospect [1]. This is defined as:

$Regret=\sum_{t=1}^{T}f_{t}(\mathbf{x}_{t})-\min_{\mathbf{x}\in
D}\sum_{t=1}^{T}f_{t}(\mathbf{x})$

Several work has been done on this under various loss function settings. Among
the most frequent are convex functions, exponentially concave functions and
strongly convex functions. These leverage the slope of the function to alter
the gradient update step in order to obtain sublinear regret bounds. OGD on
convex functions achieves $O(\sqrt{T})$ regret with step size decreasing by
$1/\sqrt{T}$ [2], and a $O(\log(T))$ regret bound on strongly convex functions
with a step-size of order $1/T$ [3]. Similarly, for exponentially concave
functions, the ONS algorithm achieves a regret bound of $O(\log(T))$ as well
[3].

Our proposed algorithm combines a variant of the meta algorithm described in
[4] using strongly convex experts from a version of Adagrad described in [5].

##  2 Related Work

Recent work on adaptive gradients set a parameter specific learning rate to
help with sparse gradients, i.e, settings where we only have a few non-zero
features. A popular method was Adagrad, which decays the learning rate faster
for parameters with larger gradients [5]. This way, we have a larger step size
for weights that have smaller magnitude. It still has a regret bound of the
order ${O}(\sqrt{T})$ similar to work by [1], but this data dependent bound
demonstrates better performance in practice.

The Adagrad algorithm without modification did not guarantee a ${O}(\log(T))$
bound in strongly convex functions. This extension to the strongly convex case
was made by [6] by one major change. First, they decayed the step size on the
parameters using the sum of squared gradients, rather than taking an L2 norm
of them. They also had a decay rate for the constant added by Adagrad
numerical stability, but this did not result in any better bounds or
performance.

At another front, there has been a push to create universal algorithms that
can perform well for most function types. As previously mentioned, the type of
alogrithm and learning rate strictly divided the type of functions into three
main categories: convex, exp-concave, and strongly convex. It would be
preferable to have algorithms that work for not only these three cases at
once, but also everything in between. This influenced work on an algorithm
called MetaGrad by [7], which works for functions that satisfy the so called
Bernstien condition. It consists of one master algorithm which spawns several
experts. Each of these slaves has its own learning rate and employs a
surrogate loss function which closely mimics the Online Newton Step (ONS). The
master alogrithm combines the predictions from each of these slaves and passes
down its gradients to them after incurring a loss. This approach only requires
gradient information and achieves $O(\sqrt{T\log T\log T})$ regret on convex
and $O(d\log T)$ on exponentially concave functions.

However, the $O(d\log T)$ regret Metagrad implies on strongly convex functions
has an $O(d)$ gap from the optimal $O(\log T)$ regret possible. This inspired
work on Maler by [4]. The main idea was to run multiple learning algorithms
for each of the three function classes in parallel. The master algorithm
continued to track the best one on the fly. This spawned three times the
number of experts as Metagrad, but enjoyed optimal regret bound on all three
function classes. We borrow inspirations the learning algorithm employed by
Maler on the strongly convex case and extend it to handle sparse gradients.

##  3 Setup

###  3.1 Notation

We denote $C\subseteq\mathbb{R}^{d}$ to be our convex decision set. We let
$\mathbf{x}_{t}\in C$ be the learner’s decision at time step $t$. We denote
our loss function at time step $t$ to be
$f_{t}(\cdot):C\rightarrow\mathbb{R}$. We let $\mathbf{g}_{t}$ be the gradient
of $f_{t}(\mathbf{x}_{t})$. We use $\lVert\cdot\rVert$ to denote the 2-norm,
and $\langle\cdot,\cdot\rangle$ to denote the Euclidean inner product. We use
boldface letters to denote vectors.

###  3.2 Assumptions and Definitions

We include two main assumptions akin to those from previous studies [4].

Assumption 1:  The gradients of all loss functions $f_{1}(\cdot)\dots
f_{T}(\cdot)$ are bounded i.e. $\forall t\in[T],\forall\mathbf{x}\in C$
$\lVert\nabla f_{t}(\mathbf{x})\rVert\leq G$ for some $G\in\mathbb{R}$.
Assumption 2:  The diameter of the decision set $C$ is bounded i.e.
$\forall\mathbf{x}_{1},\mathbf{x}_{2}\in C$,
$\lVert\mathbf{x}_{1}-\mathbf{x}_{2}\rVert\leq D$ for some $D\in\mathbb{R}$.
Definition 1: Strong Convexity A function $f:C\rightarrow\mathbb{R}$ is
$\lambda$-strongly convex if $\forall\mathbf{x}_{1},\mathbf{x}_{2}\in C$,

$f(\mathbf{x}_{1})\geq f(\mathbf{x}_{2})+\langle\nabla
f(\mathbf{x}_{2}),\mathbf{x}_{1}-\mathbf{x}_{2}\rangle+\frac{\lambda}{2}\lVert\mathbf{x}_{1}-\mathbf{x}_{2}\rVert^{2}$

Assumption 3:  All loss functions $f_{1}(\cdot)\dots f_{T}(\cdot)$ are
$\lambda$-strongly convex.
Definition 2: A-weighted norm Let $A\in\mathbb{R}^{d\times d}$ and let
$\mathbf{x}\in\mathbb{R}^{d}$. We define the $A$-weighted norm to be the
matrix norm
$\lVert\mathbf{x}\rVert_{A}^{2}=\mathbf{x}^{\top}A\mathbf{x}=\sum_{i,j=1}^{d}A_{ij}\mathbf{x}_{i}\mathbf{x}_{j}$
[6].
Definition 3: Weighted Projection We define a weighted projection function
onto the convex set $C$ with respect to the $A$-weighted norm to be

$\Pi_{C}^{A}(\mathbf{x}_{1})=\arg\min_{\mathbf{x}\in
C}\lVert\mathbf{x}_{1}-\mathbf{x}\rVert^{2}_{A}$

##  4 Algorithm

The algorithm draws inspiration from the strongly convex experts in Maler. It
similarly follows a two hierarchical structure. At the lower level, we have a
set of experts, each with their own learning rates $\eta$. Because we are only
dealing with the strongly convex case, all of them use the same learning
algorithm (Algorithm 2). At the higher level, we have a meta algorithm which
tracks the best one using a tilted exponential weighted average [7].

Input: Grid of learning rates $\eta_{1},\eta_{2},\dots,\eta_{k}$ and weights
$\pi_{1}^{\eta_{1}},\pi_{1}^{\eta_{2}},\dots,\pi_{1}^{\eta_{k}}$

Get prediction $\mathbf{x}_{t}^{\eta}$ from Algorithm 2
Play
$\mathbf{x}_{t}=\frac{\sum_{\eta}\pi_{t}^{\eta}\eta\mathbf{x}_{t}^{\eta}}{\sum_{\eta}\pi_{t}^{\eta}\eta}$
Observe gradient $\mathbf{g}_{t}$ and send $\mathbf{x}_{t},\mathbf{g}_{t}$ to
all experts
For all $\eta$, update weights according to
$\pi_{t+1}^{\eta}=\frac{\pi_{t}^{\eta}\exp{\left(-s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})\right)}}{\sum_{\eta}\pi_{t}^{\eta}\exp{\left(-s_{t}^{\eta}(\mathbf{x}_{t}^{\eta}\right))}}$

Algorithm 1 Meta Algorithm

We define the surrogate loss function to be

$s_{t}^{\eta}(\mathbf{x})=-\eta(\mathbf{x}_{t}-\mathbf{x})^{\top}\mathbf{g}_{t}+\eta^{2}G^{2}\lVert\mathbf{x}_{t}-\mathbf{x}\rVert^{2}$
(1)

Note, by taking the second derivative of the surrogate loss function, we see
that it is $\eta^{2}G^{2}$ strongly convex.

Input: $\mathbf{v}_{0}=\mathbf{0}\in\mathbb{R}^{d}$,$\eta\in\mathbb{R}$,
$\mathbf{x}_{1}^{\eta}=0$, $\delta>0$

for  _$t=1\dots T$_ do

Send $\mathbf{x}_{t}^{\eta}$ to Algorithm 1
Receive $\mathbf{x}_{t},\mathbf{g}_{t}$ from Algorithm 1
$\nabla
s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})=\eta\mathbf{g}_{t}+2\eta^{2}G^{2}(\mathbf{x}_{t}^{\eta}-\mathbf{x}_{t})$
$\mathbf{v}_{t}=\mathbf{v}_{t-1}+\nabla
s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})\odot\nabla
s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})$
$A_{t}=\text{diag}(\mathbf{v}_{t})+\delta I$
$\mathbf{x}_{t+1}^{\eta}=\Pi_{C}^{A_{t}}\left(\mathbf{x}_{t}^{\eta}-\alpha
A^{-1}_{t}\nabla s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})\right)$, where
$\alpha=\frac{1}{4\eta^{2}}$

end for

Algorithm 2 SC-Adagrad expert algorithm

Experts: Each expert uses an adaptive algorithm with an update rule derived
from the Adagrad for strongly convex function [6]. However, in contrast to the
original Adagrad, this adaptation is on the gradients of the surrogate losses,
rather than the $\mathbf{g}_{t}$ incurred by a single learner. We maintain
$\lceil\frac{1}{2}\log_{2}T\rceil$ experts. For
$i=0,1,2,...\frac{1}{2}\log_{2}T$, each expert is configured with:

$\eta_{i}=\frac{2^{-i}}{5GD}$

and

$\pi_{1}^{\eta_{i}}=\frac{c}{3(i+1)(i+2)}$

These values for $\eta$, $\pi$ and $\alpha$ were carefully chosen to satisfy
certain criteria which will help in proving regret bounds.

##  5 Regret Analysis

###  5.1 Preliminaries

#### Lemma 1

[Lemma 3.2 from [6]] Let assumptions A1 and A2 hold, then for all $T\geq 1$
and $A_{t},\delta_{t}$ defined from algorithm 2 and $\nabla
s_{t}=s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})$, we have:

$\sum_{t=1}^{T}\langle\nabla s_{t},A^{-1}_{t}\nabla
s_{t}\rangle\leq\sum_{i=1}^{d}\log\Bigg{(}\frac{||\nabla
s_{1:T,i}||^{2}+\delta}{\delta}\Bigg{)}$

#### Lemma 2

$\forall t$, let $A_{t}$ be positive semi-definite diagonal matrices as which
follows the update rules in Algorithm 2. Similarly, let
$\alpha\geq\frac{1}{4\eta^{2}}$ for some constant $\eta$ and let
$\mu=2\eta^{2}G^{2}$. Then the following holds:

$\displaystyle A_{t}-A_{t-1}-2\alpha\text{diag}(\mu)\preceq 0$

and

$\displaystyle A_{1}-2\alpha\text{diag}(\mu)\preceq\delta$

Proof: In our update rule we have:

$\displaystyle A_{t}$ $\displaystyle=\text{diag}(v_{t})+\text{diag}(\delta)$
$\displaystyle=v_{t-1}+\nabla s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})\bullet\nabla
s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})+\delta$ $\displaystyle=A_{t-1}+\nabla
s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})\bullet\nabla
s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})$ $\displaystyle\leq A_{t-1}+{(G^{s})}^{2}$

where, $\nabla s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})\bullet\nabla
s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})$ has an upper bound of $(G^{s})^{2}$.
This can be rearranged as:

$\displaystyle A_{t}-A_{t-1}$ $\displaystyle\leq{(G^{s})}^{2}$
$\displaystyle\Rightarrow A_{t}-A_{t-1}-2\alpha\text{diag}(\mu)$
$\displaystyle\leq{(G^{s})}^{2}-2\alpha\text{diag}(\mu)$
$\displaystyle\leq{(G^{s})}^{2}-2\left(\frac{1}{4\eta^{2}}\right)(2\eta^{2}{(G^{s})}^{2})$
$\displaystyle\leq{(G^{s})}^{2}-{(G^{s})}^{2}=0$

This proves our the first inequality. For the second one, since we have
$\mathbf{v}_{0}=0$ in our algorithm, we observe that
$A_{1}=s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})+s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})+\text{diag}(\delta)$.
With the same setting of $\alpha$, we have:

$\displaystyle A_{1}$ $\displaystyle\leq{(G^{s})}^{2}+\delta$ $\displaystyle
A_{1}-2\alpha\text{diag}(\mu)$
$\displaystyle\leq{(G^{s})}^{2}-2\alpha\text{diag}(\mu)+\delta$
$\displaystyle\leq\delta$

This proves the second inequality.

###  5.2 Meta Regret

The regret of the meta algorithm is defined as the difference of the surrogate
losses of the meta learner and the surrogate losses of a given expert. It is
given by:

$\displaystyle
R_{T}^{Meta}=\sum_{t=1}^{T}s_{t}(\mathbf{x}_{t})-\sum_{t=1}^{T}s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})$
(5.2.1)

#### Lemma 3

[Lemma 1 in [4]] For every grid point $\eta$ we have:

$\sum_{t=1}^{T}s_{t}(\mathbf{x}_{t})-\sum_{t=1}^{T}s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})\leq
2\ln\bigg{(}\sqrt{3}\left(\frac{1}{2}\log_{2}T+3\right)\bigg{)}$

The proof is attached in the appendix.

###  5.3 Expert Regret

The regret of the expert algorithm is defined as the difference of the
surrogate losses of the expert and the surrogate losses of the best decision
in hindsight. It is given by:

$\displaystyle
R_{T}^{Expert}=\sum_{t=1}^{T}s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})-\sum_{t=1}^{T}s_{t}^{\eta}(\mathbf{x}^{*})$
(5.3.1)

Let us denote $\nabla s_{t}=\nabla s_{t}^{\eta}(\mathbf{x}^{\eta})$

#### Lemma 4

$\displaystyle(\mathbf{x}_{t}^{\eta}-\mathbf{x}^{*})^{\top}\nabla s_{t}$
$\displaystyle\leq\frac{||\mathbf{x}_{t}^{\eta}-\mathbf{x}^{*}||^{2}_{A_{t}}-||\mathbf{x}_{t+1}^{\eta}-\mathbf{x}^{*}||^{2}_{A_{t}}}{2\alpha}+\frac{\alpha}{2}\langle\nabla
s_{t},A_{t}^{-1}\nabla s_{t}\rangle$

Proof :
From the update rule, we have :

$\displaystyle||\mathbf{x}_{t+1}^{\eta}-\mathbf{x}^{*}||_{A_{t}}$
$\displaystyle=\left\|\Pi_{C}^{A_{t}}(\mathbf{x}_{t}^{\eta}-\alpha
A_{t}^{-1}\nabla s_{t}-\mathbf{x}^{*})\right\|^{2}_{A_{t}}$
$\displaystyle\leq\left\|\mathbf{x}_{t}^{\eta}-\alpha A_{t}^{-1}\nabla
s_{t}-\mathbf{x}^{*}\right\|^{2}_{A_{t}}$
$\displaystyle=||\mathbf{x}_{t}^{\eta}-\mathbf{x}^{*}||^{2}_{A_{t}}+\alpha^{2}||A_{t}^{-1}\nabla
s_{t}||^{2}_{A_{t}}-2\alpha(\mathbf{x}_{t}^{\eta}-\mathbf{x}^{*})^{\top}A_{t}(A_{t}^{-1}\nabla
s_{t})$ $\displaystyle\Rightarrow
2\alpha(\mathbf{x}_{t}^{\eta}-\mathbf{x}^{*})^{\top}\nabla s_{t}$
$\displaystyle\leq||\mathbf{x}_{t}^{\eta}-\mathbf{x}^{*}||^{2}_{A_{t}}-||\mathbf{x}_{t+1}^{\eta}-\mathbf{x}^{*}||_{A_{t}}+\alpha^{2}||A_{t}^{-1}\nabla
s_{t}||^{2}_{A_{t}}$
$\displaystyle\Rightarrow(\mathbf{x}_{t}^{\eta}-\mathbf{x}^{*})^{\top}\nabla
s_{t}$
$\displaystyle\leq\frac{||\mathbf{x}_{t}^{\eta}-\mathbf{x}^{*}||^{2}_{A_{t}}-||\mathbf{x}_{t+1}^{\eta}-\mathbf{x}^{*}||_{A_{t}}}{2\alpha}+\frac{\alpha}{2}\langle\nabla
s_{t},A_{t}^{-1}\nabla s_{t}\rangle$

#### Lemma 5

The expert regret is given by:

$\displaystyle
R_{T}^{Expert}=\sum_{t=1}^{T}s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})-s_{t}^{\eta}(\mathbf{x}^{*})\leq\frac{D^{2}d\delta}{2\alpha}+\frac{\alpha}{2}\sum_{i=1}^{d}\log\Bigg{(}\frac{||\nabla
s_{1:T,i}||^{2}+\delta}{\delta}\Bigg{)}$

Proof: We can start from the definition of strong convexity

$\displaystyle\sum_{t=1}^{T}s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})-s_{t}^{\eta}(\mathbf{x}^{*})$
$\displaystyle\leq\sum_{t=1}^{T}\langle\nabla
s_{t}^{\eta}(\mathbf{x}^{\eta}),\mathbf{x}_{t}-\mathbf{x}^{*}\rangle-||\mathbf{x}_{t}^{\eta}-\mathbf{x}^{*}||^{2}_{\text{diag}(\mu)}$
$\displaystyle\leq\sum_{t=1}^{T}\frac{||\mathbf{x}_{t}^{\eta}-\mathbf{x}^{*}||^{2}_{A_{t}}-||\mathbf{x}_{t+1}^{\eta}-\mathbf{x}^{*}||^{2}_{A_{t}}}{2\alpha}+\frac{\alpha}{2}\langle\nabla
s_{t},A_{t}^{-1}\nabla
s_{t}\rangle-||\mathbf{x}_{t}^{\eta}-\mathbf{x}^{*}||^{2}_{\text{diag}(\mu)}$
(Using Lemma 5.3)
$\displaystyle=\frac{||\mathbf{x}_{1}^{\eta}-\mathbf{x}^{*}||^{2}_{A_{1}}}{2\alpha}+\sum_{t=2}^{T}\frac{||\mathbf{x}_{t}^{\eta}-\mathbf{x}^{*}||^{2}_{A_{t}}-||\mathbf{x}_{t}^{\eta}-\mathbf{x}^{*}||^{2}_{A_{t-1}}}{2\alpha}-\frac{||\mathbf{x}_{T+1}^{\eta}-\mathbf{x}^{*}||^{2}_{A_{T}}}{2\alpha}$
$\displaystyle+\sum_{t=1}^{T}\frac{\alpha}{2}\langle\nabla
s_{t},A_{t}^{-1}\nabla
s_{t}\rangle-\sum_{t=1}^{T}||\mathbf{x}_{t}^{\eta}-\mathbf{x}^{*}||^{2}_{\text{diag}(\mu)}$
(Separating the first and last terms)
$\displaystyle\leq\frac{||x_{1}^{\eta}-\mathbf{x}^{*}||^{2}_{A_{1}-2\alpha\text{diag}(\mu)}}{2\alpha}+\sum_{t=2}^{T}\frac{||\mathbf{x}_{t}^{\eta}-\mathbf{x}^{*}||^{2}_{A_{t}-A_{t-1}-2\alpha\text{diag}(\mu)}}{2\alpha}$
$\displaystyle+\frac{\alpha}{2}\langle\nabla s_{t},A_{t}^{-1}\nabla
s_{t}\rangle$
$\displaystyle\leq\frac{||\mathbf{x}_{1}^{\eta}-\mathbf{x}^{*}||^{2}_{\delta}}{2\alpha}+\sum_{t=2}^{T}\frac{||\mathbf{x}_{t}^{\eta}-\mathbf{x}^{*}||^{2}_{0}}{2\alpha}+\frac{\alpha}{2}\langle\nabla
s_{t},A_{t}^{-1}\nabla s_{t}\rangle$ (Using Lemma 5.1)
$\displaystyle\leq\frac{||\mathbf{x}_{1}^{\eta}-\mathbf{x}^{*}||^{2}_{\delta}}{2\alpha}+\frac{\alpha}{2}\langle\nabla
s_{t},A_{t}^{-1}\nabla s_{t}\rangle$ (because second term $\leq$ 0)
$\displaystyle\leq\frac{D^{2}d\delta}{2\alpha}+\frac{\alpha}{2}\langle\nabla
s_{t},A_{t}^{-1}\nabla s_{t}\rangle$ (Using A1)
$\displaystyle\leq\frac{D^{2}d\delta}{2\alpha}+\frac{\alpha}{2}\sum_{i=1}^{d}\log\Bigg{(}\frac{||\nabla
s_{1:T,i}||^{2}+\delta}{\delta}\Bigg{)}$ (Using Lemma 5.1 )

###  5.4 Total Regret

The total regret is given by

$\displaystyle R_{T}$
$\displaystyle=\sum_{t=1}^{T}f_{t}(\mathbf{x}_{t})-f_{t}(\mathbf{x}^{*})$
$\displaystyle\leq\sum_{t=1}^{T}\mathbf{g}_{t}^{\top}(\mathbf{x}_{t}-\mathbf{x}^{*})-\frac{\mu}{2}||\mathbf{x}_{t}-\mathbf{x}^{*}||^{2}$
(from definition of strongly convex functions)
$\displaystyle=\sum_{t=1}^{T}\frac{-s_{t}^{\eta}(\mathbf{x}^{*})+\eta^{2}G^{2}||\mathbf{x}_{t}-\mathbf{x}^{*}||}{\eta}-\frac{\mu}{2}||\mathbf{x}_{t}-\mathbf{x}^{*}||^{2}$
(from the definition of surrogate loss function)
$\displaystyle\leq\sum_{t=1}^{T}\frac{s_{t}^{\eta}(\mathbf{x}_{t})-s_{t}^{\eta}(\mathbf{x}^{*})}{\eta}+\eta\sum_{t=1}^{T}G^{2}||\mathbf{x}_{t}-\mathbf{x}^{*}||-\frac{\mu}{2}\sum_{t=1}^{T}||\mathbf{x}_{t}-\mathbf{x}^{*}||^{2}$
$\displaystyle=\sum_{t=1}^{T}\frac{s_{t}^{\eta}(\mathbf{x}_{t})-s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})}{\eta}+\sum_{t=1}^{T}\frac{s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})-s_{t}^{\eta}(\mathbf{x}^{*})}{\eta}+\eta\sum_{t=1}^{T}G^{2}||\mathbf{x}_{t}-\mathbf{x}^{*}||-\frac{\mu}{2}\sum_{t=1}^{T}||\mathbf{x}_{t}-\mathbf{x}^{*}||^{2}$
$\displaystyle\leq\frac{2}{\eta}\log\left(\sqrt{3}\left(\frac{1}{2}\log(T)+3\right)\right)+\frac{1}{\eta}\frac{D^{2}d\delta}{2\alpha}+\frac{\alpha}{2}\sum_{i=1}^{d}\log\Bigg{(}\frac{||\nabla
s_{1:T,i}||^{2}+\delta}{\delta}\Bigg{)}$ (Using lemma 5.2, 5.3)
$\displaystyle+\eta
V_{T}^{S}-\frac{\mu}{2}\sum_{t=1}^{T}||\mathbf{x}_{t}-\mathbf{x}^{*}||^{2}$
(5.4.1)

where, $V_{T}^{S}=\sum_{t=1}^{T}G^{2}||x_{t}-x^{*}||^{2}$.

Let us denote

$\displaystyle E_{T}$
$\displaystyle=\frac{D^{2}d\delta}{2\alpha}+\frac{\alpha}{2}\sum_{i=1}^{d}\log\Bigg{(}\frac{||\nabla
s_{1:T,i}||^{2}+\delta}{\delta}\Bigg{)}=O\left(d\log\left(\frac{T}{\delta}\right)\right)$
(5.4.2) $\displaystyle A_{T}$
$\displaystyle=2\log\left(\sqrt{3}\left(\frac{1}{2}\log(T)+3\right)\right)+E_{T}$
(5.4.2)

Using (5.4.2) in (5.4.1), we have

$\displaystyle R_{T}$ $\displaystyle\leq\frac{A_{T}}{\eta}+\eta
V_{T}^{S}-\frac{\mu}{2}\sum_{t=1}^{T}||\mathbf{x}_{t}-\mathbf{x}^{*}||^{2}$
(5.4.3)

The optimal $\hat{\eta}$ to minimize $A_{T}+\eta V_{T}^{S}$ is ([4])

$\displaystyle\hat{\eta}=\sqrt{\frac{A_{T}}{V_{T}^{S}}}\geq\frac{1}{5GD}$
(5.4.4)

If $\hat{\eta}\leq\frac{1}{5GD}$, then by construction there exists a grid
point $\eta\in[\frac{\hat{\eta}}{2},\hat{\eta}]$. Thus.

$\displaystyle R_{T}$ $\displaystyle\leq\eta
V_{T}^{S}+\frac{A_{T}}{\eta}-\frac{\mu}{2}\sum_{t=1}^{T}||\mathbf{x}_{t}-\mathbf{x}^{*}||^{2}$
$\displaystyle\leq\hat{\eta}V_{T}^{S}+\frac{2A_{T}}{\hat{\eta}}-\frac{\mu}{2}\sum_{t=1}^{T}||\mathbf{x}_{t}-\mathbf{x}^{*}||^{2}$
$\displaystyle=3\sqrt{V_{T}^{S}A_{T}}-\frac{\mu}{2}\sum_{t=1}^{T}||\mathbf{x}_{t}-\mathbf{x}^{*}||^{2}$
(5.4.5)

If $\hat{\eta}\geq\frac{1}{5GD}$, by (5.4.4), we have

$\displaystyle V_{T}^{S}\leq 25G^{2}D^{2}A_{T}$ (5.4.6)

Using (5.4.4), (5.4.6) and $\hat{\eta}=\frac{1}{5GD}$ in (5.4.3), we have

$\displaystyle R_{T}$ $\displaystyle\leq
2\sqrt{A_{T}V_{T}^{S}}-\frac{\mu}{2}\sum_{t=1}^{T}||\mathbf{x}_{t}-\mathbf{x}^{*}||^{2}$
$\displaystyle\leq
10GDA_{T}-\frac{\mu}{2}\sum_{t=1}^{T}||\mathbf{x}_{t}-\mathbf{x}^{*}||^{2}$
(5.4.7)

Hence, combining (5.4.5) and (5.4.6), we have

$\displaystyle R_{T}$ $\displaystyle\leq
3\sqrt{V_{T}^{S}A_{T}}+10GDA_{T}-\frac{\mu}{2}\sum_{t=1}^{T}||\mathbf{x}_{t}-\mathbf{x}^{*}||^{2}$
(5.4.8)

Let $x=V_{T}^{S}$ and $y=A_{T}$. For $x,y,\gamma\geq 0$

$\displaystyle\sqrt{xy}$
$\displaystyle\leq\frac{\gamma}{2}x+\frac{1}{2\gamma}y$
$\displaystyle\Rightarrow 3\sqrt{V_{T}^{S}A_{T}}$
$\displaystyle\leq\frac{3\gamma}{2}V_{T}^{S}+\frac{3}{2\gamma}A_{T}$ (5.4.9)

Using (5.4.9) in (5.4.8), we have

$\displaystyle\noindent R_{T}$
$\displaystyle\leq\frac{3\gamma}{2}V_{T}^{S}+\frac{3}{2\gamma}A_{T}+10GDA_{T}-\frac{\mu}{2}\sum_{t=1}^{T}||x_{t}-x^{*}||^{2}$
$\displaystyle=\frac{3\gamma}{2}V_{T}^{S}+\frac{3}{2\gamma}A_{T}+10GDA_{T}-\frac{\mu}{2G^{2}}V_{T}^{S}$
$\displaystyle=\left(\frac{3\gamma}{2}-\frac{\mu}{2G^{2}}\right)V_{T}^{S}+\left(\frac{3}{2\gamma}+10GD\right)A_{T}$
$\displaystyle\text{Setting }\gamma$ $\displaystyle=\frac{\mu}{3G^{2}}\text{,
we have}$ $\displaystyle R_{T}$
$\displaystyle=\left(\frac{G^{2}}{\mu}+10GD\right)A_{T}$
$\displaystyle=\left(\frac{G^{2}}{\mu}+10GD\right)\left(2\log\left(\sqrt{3}\left(\frac{1}{2}\log(T)+3\right)\right)+\frac{D^{2}d\delta}{2\alpha}+\frac{\alpha}{2}\sum_{i=1}^{d}\log\Bigg{(}\frac{||\nabla
s_{1:T,i}||^{2}+\delta}{\delta}\Bigg{)}\right)$

Remark : In [6] the authors claim that

$\displaystyle\sum_{i=1}^{d}\log\Bigg{(}\frac{||\nabla
s_{1:T,i}||^{2}+\delta_{T,i}}{\delta_{1,i}}\Bigg{)}=O(\log T)$

And hence, the overall regret is $O(\log T)$ for strongly-convex functions.
Even with the monotonic-decay factor, the term

$\displaystyle\frac{||\nabla
s_{1:T,i}||^{2}}{\delta_{1,i}}\leq\frac{TG^{2}}{\delta_{1,i}}$

Hence, the regret term is bound by $O(d\log(T))$ rather than $O(\log(T))$.
This result hasn’t been explained in the paper which seems like a weakness on
the results. Hence, we take state the final regret bound as $O(d\log(T))$.

##  6 Conclusion and Future Work

By taking a Metagrad like approach to Adagrad, we can have created a parameter
free version of Strongly Convex Adagrad. This enjoys $O(d\log T)$ regret bound
in strongly convex functions. If every expert employs SAdam [8], the regret
bound can be improved to $O(logT)$ instead. This work can be extended to make
parameter free versions of other algorithms such as RMSProp or Adam.

## References

  * [1] Elad Hazan.  Introduction to online convex optimization, 2021.
  * [2] Martin Zinkevich.  Online convex programming and generalized infinitesimal gradient ascent.  In Proceedings of the Twentieth International Conference on International Conference on Machine Learning, ICML’03, page 928–935. AAAI Press, 2003.
  * [3] Elad Hazan, Adam Kalai, Satyen Kale, and Amit Agarwal.  Logarithmic regret algorithms for online convex optimization.  pages 499–513, 06 2006.
  * [4] Guanghui Wang, Shiyin Lu, and Lijun Zhang.  Adaptivity and optimality: A universal algorithm for online convex optimization, 2019.
  * [5] John Duchi, Elad Hazan, and Yoram Singer.  Adaptive subgradient methods for online learning and stochastic optimization.  Journal of Machine Learning Research, 12(61):2121–2159, 2011.
  * [6] Mahesh Chandra Mukkamala and Matthias Hein.  Variants of rmsprop and adagrad with logarithmic regret bounds, 2017.
  * [7] Tim van Erven and Wouter M. Koolen.  Metagrad: Multiple learning rates in online learning, 2016.
  * [8] Guanghui Wang, Shiyin Lu, Quan Cheng, Wei wei Tu, and Lijun Zhang.  Sadam: A variant of adam for strongly convex functions.  In International Conference on Learning Representations, 2020.

##  Appendix A Proof of Meta Regret

This proof follows closely from the proof of the meta regret bound in [4]. We
define a potential function

$\displaystyle\Phi(T)=\sum_{\eta}\pi_{1}^{\eta}\exp\left(-\sum_{t=1}^{T}s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})\right)$

We have

$\displaystyle\Phi(T+1)-\Phi(T)$
$\displaystyle=\sum_{\eta}\pi_{1}^{\eta}\exp\left(-\sum_{t=1}^{T}s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})\right)\exp\left(-s_{T+1}^{\eta}\left(\mathbf{x}_{T+1}^{\eta}\right)-1\right)$
$\displaystyle\leq\sum_{\eta}\pi_{1}^{\eta}\exp\left(-\sum_{t=1}^{T}s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})\right)\eta(\mathbf{x}_{T+1}-\mathbf{x}_{T+1}^{\eta})^{\top}\mathbf{g}_{t}$
$\displaystyle=(\mathbf{a}_{T}\mathbf{x}_{T+1}-\mathbf{b}_{T})^{\top}\mathbf{g}_{t}$

where

$\displaystyle\mathbf{a}_{T}$
$\displaystyle=\sum_{\eta}\pi_{1}^{\eta}\exp\left(-\sum_{t=1}^{T}s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})\right)\eta$
$\displaystyle\mathbf{b}_{T}$
$\displaystyle=\sum_{\eta}\pi_{1}^{\eta}\exp\left(-\sum_{t=1}^{T}s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})\right)\eta\mathbf{x}_{T+1}^{\eta}$

By definition,

$\displaystyle\mathbf{x}_{T+1}$
$\displaystyle=\frac{\sum_{\eta}\pi_{T+1}^{\eta}\eta\mathbf{x}_{T+1}^{\eta}}{\sum_{\eta}\pi_{T+1}^{\eta}\eta}$
$\displaystyle=\frac{\sum_{\eta}\pi_{1}^{\eta}\exp\left(-\sum_{t=1}^{T+1}s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})\right)\eta\mathbf{x}_{T}^{\eta}}{\sum_{\eta}\pi_{1}^{\eta}\exp(-\sum_{t=1}^{T}s_{t}^{\eta}(\mathbf{x}_{t}^{\eta}))\eta}$
$\displaystyle=\frac{\mathbf{b}_{T}}{\mathbf{a}_{T}}$

This implies that $\Phi(T+1)-\Phi(T)\leq 0$ and
$1=\Phi(0)\geq\Phi(1)\geq\dots\Phi(T)\geq 0$. Then for any $\eta$

$\displaystyle\ln\left(\pi_{1}^{\eta}\exp\left(-\sum_{t=1}^{T}s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})\right)\right)$
$\displaystyle\leq 0$
$\displaystyle\implies\ln\left(\frac{1}{\pi_{1}^{\eta}}\right)+\sum_{t=1}^{T}s_{t}^{\eta}(\mathbf{x}_{t}^{\eta})\geq
0$

By definition, $\pi_{1}^{\eta_{i}}=\frac{c}{3(i+1)(i+2)}$, where
$c=1+\frac{1}{1+k}$, $i\in\\{0,\dots,k\\}$ and
$k=\left\lceil\frac{1}{2}\log_{2}(T)\right\rceil$. Therefore,

$\displaystyle\ln\left(\frac{1}{\pi_{1}^{\eta}}\right)$
$\displaystyle=\ln\left(\frac{3(i+1)(i+2)}{c}\right)$
$\displaystyle\leq\ln\left(3(1+k)(2+k)\right)$
$\displaystyle\leq\ln\left(\left(\sqrt{3}\left(3+\frac{1}{2}\log_{2}(T)\right)\right)^{2}\right)$
$\displaystyle=2\ln\left(\sqrt{3}\left(3+\frac{1}{2}\log_{2}(T)\right)\right)$