rational mappings and complex dynamics1997ISSN 0022-2518Indiana Univ. Math.
J.463897–932LinkReview MathReviews@article{Russakovskii–Shiffman, author =
{Russakovskii, Alexander}, author = {Shiffman, Bernard}, title = {Value
distribution for sequences of rational mappings and complex dynamics}, date =
{1997}, issn = {0022-2518}, journal = {Indiana Univ. Math. J.}, volume = {46},
number = {3}, pages = {897\ndash 932}, url =
{https://doi.org/10.1512/iumj.1997.46.1441}, review = {\MR{1488341}}}
RullgårdHansPolynomial amoebas and convexityPh.D. ThesisDepartment of
Mathematics, Stockholm University2001@thesis{Rull-thesis, author =
{Rullg{\aa}rd, Hans}, title = {Polynomial amoebas and convexity}, type =
{Ph.D. Thesis}, school = {Department of Mathematics, Stockholm University},
date = {2001}} SturmfelsBerndTevelevJeniaElimination theory for tropical
varieties2008ISSN 1073-2780Math. Res. Lett.153543–562LinkReview
MathReviews@article{Sturmfels-Tevelev, author = {Sturmfels, Bernd}, author =
{Tevelev, Jenia}, title = {Elimination theory for tropical varieties}, date =
{2008}, issn = {1073-2780}, journal = {Math. Res. Lett.}, volume = {15},
number = {3}, pages = {543\ndash 562}, url =
{https://doi.org/10.4310/MRL.2008.v15.n3.a14}, review = {\MR{2407231}}}
TevelevJeniaCompactifications of subvarieties of tori2007ISSN 0002-9327Amer.
J. Math.12941087–1104LinkReview MathReviews@article{Tevelev, author =
{Tevelev, Jenia}, title = {Compactifications of subvarieties of tori}, date =
{2007}, issn = {0002-9327}, journal = {Amer. J. Math.}, volume = {129}, number
= {4}, pages = {1087\ndash 1104}, url =
{https://doi.org/10.1353/ajm.2007.0029}, review = {\MR{2343384}}} ViroOlegFrom
the sixteenth Hilbert problem to tropical geometry2008ISSN 0289-2316Jpn. J.
Math.32185–214LinkReview MathReviews@article{Viro-Hilbert, author = {Viro,
Oleg}, title = {From the sixteenth {H}ilbert problem to tropical geometry},
date = {2008}, issn = {0289-2316}, journal = {Jpn. J. Math.}, volume = {3},
number = {2}, pages = {185\ndash 214}, url =
{https://doi.org/10.1007/s11537-008-0832-6}, review = {\MR{2465564}}}


# Pinpointing the Memory Behaviors of DNN Training
††thanks: This work is supported by the National Key R&D Program of China
(under Grant No. 2017YFB1003103) and the Science Fund for Creative Research
Groups of the National Natural Science Foundation of China (under Grant No.
61521092).

Jiansong Li†‡, Xiao Dong∗, Guangli Li†‡, Peng Zhao¶, Xueying Wang†‡, Xiaobing
Chen†‡,
Xianzhi Yu¶, Yongxin Yang†‡, Zihan Jiang†‡, Wei Cao†, Lei Liu†, Xiaobing
Feng†‡ †State Key Laboratory of Computer Architecture, Institute of Computing
Technology, CAS, Beijing, China
‡University of Chinese Academy of Sciences, Beijing, China
∗Youtu Lab, Tencent, Shanghai, China¶Huawei Technology Co., Ltd, Beijing,
China
{lijiansong, liguangli, wangxueying, chenxiaobing, yangyongxin, jiangzihan,
caowei, liulei<EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The training of deep neural networks (DNNs) is usually memory-hungry due to
the limited device memory capacity of DNN accelerators. Characterizing the
memory behaviors of DNN training is critical to optimize the device memory
pressures. In this work, we pinpoint the memory behaviors of each device
memory block of GPU during training by instrumenting the memory allocators of
the runtime system. Our results show that the memory access patterns of device
memory blocks are stable and follow an iterative fashion. These observations
are useful for the future optimization of memory-efficient training from the
perspective of raw memory access patterns.

###### Index Terms:

DNN Training, Memory Access Patterns, Workload Characterizing

##  I Introduction

Nowadays, DNNs trend to go deeper and wider for higher accuracy [1, 2]. The
statistics from OpenAI show that the computation and memory requirements of
DNNs are doubled every 3.4 months in recent years [3]. However, the memory
capacity constraint of DNN accelerators limits the upper bound on the size of
DNNs that can be offloaded to and trained. For instance, the DRAM size of
Nvidia GPU with the latest Ampere architecture is 40GB [4]. Even such a large
memory capacity fails to satisfy the demand of some DNNs. For example, the
typical Inception-V4 [5] requests up to 45GB of device memory to keep the
entire DNN on the GPU during training [6]. It is important to reduce the
memory pressures of DNN training. Characterizing the memory behaviors of DNN
training is critical to optimize the memory usage of DNN workloads.

##  II Benchmarking Method

To work around the memory pressures of DNN training, we carefully pinpoint the
memory behaviors (including malloc, free, read, write) of each device memory
block by manually instrumenting the memory allocators of PyTorch’s runtime
system. We firstly conduct a case study of the trivial MLP to show the memory
access patterns of DNN training. To see where the device memory is spent and
which kinds of memory content is the bottleneck, we also make a detailed
memory occupation breakdown for the typical linear and non-linear [6] DNNs
with different batch size and layer structures.

Figure 1:  Layer topology of MLP. Note that star means the mat_mul operator,
plus is add_bias, $f$ is ReLU activation. Figure 2:  Gantt chart of the first
five iterations in MLP training on the Nvidia Titan X Pascal GPU.

##  III Experimental Results

Iterative Memory Access Patterns. Due to space limit, we only use a trivial
MLP to show the memory access patterns of DNN training (our observations are
also appliable to other DNNs). Its layer topology is shown in Fig. 1. The
shape of $W_{0}$ is (2, 12288), $b_{0}$: (12288), $W_{1}$: (12288, 2),
$b_{1}$: (2). In Fig. 2, each rectangle represents the memory access to a
certain device memory block of GPU. Width of the rectangle is the elapsed time
from the allocation time to the free time of current device memory block,
which can be regarded as the lifetime of current memory block. Height of the
rectangle denotes the size of current device memory block. We observe that
there are obvious iterative memory access patterns in the first five rounds of
MLP training. Besides, the overlap between rectangles along the y-axis
represents the overlap of the live ranges of different memory blocks. The
blank space between two rectangles along the y-axis represents the memory
fragments of GPU. We can also observe that there are fewer memory fragments
during MLP training.

(a) CDF [7].

(b) Violin plot [8].

Figure 3:  CDF and violin plot of the memory block access intervals in MLP
training on the Nvidia Titan X Pascal GPU.

Furthermore, Fig. 3 shows the distributions of these memory behaviors. The
access time interval (abbr $ATI$) is the elapsed time between two adjacent
memory access to the same memory block. As presented in Fig. 3(b), the $ATI$s
of most memory behaviors range from 10us to 25us, and their distributions are
relatively concentrated. Besides, Fig. 3(a) shows that the $ATI$s of 90% of
the memory behaviors are less than 25us. Suppose we want to reduce the device
memory pressure by swapping data back and forth between the host CPUs and the
device GPUs [9, 10], the $ATI$ of current device memory block is $T$. Assume
the memory bandwidth from device GPU to host CPU is $B_{d2h}$ and that of host
to device is $B_{h2d}$. Then the maximum value of memory swapping size $S$
without sacrificing the runtime performance should meet the following
equation:

$\displaystyle\frac{S}{B_{d2h}}+\frac{S}{B_{h2d}}\leq T,$ (1) $\displaystyle
S\leq\frac{T}{1/B_{d2h}+1/B_{h2d}}$

We measured the memcpy bandwidth between the device GPU and host CPU with the
bandwidthTest tool from CUDA SDK samples [11]. The pinned memory transfer
bandwidth from host to the device is 6.3GB/s, from device to host is 6.4GB/s.
So the maximum value of memory swapping size without sacrificing the runtime
performance is about 79.37KB (by equation 1,
$S\leq\frac{25us}{1/6.4GB/s+1/6.3GB/s}=79.37KB$). This is just a drop in the
bucket for reducing the memory footprint of MLP training. It seems that memory
swapping is promising less.

Figure 4:  Pair-wise $ATI$ and the corresponding device memory block size of
each memory behavior during MLP training on the Nvidia Titan X Pascal GPU.
Note that the x-axis is the index of each memory behavior.

The detailed memory behaviors of Fig. 4, including $ATI$s along with the
corresponding memory block size in MLP training, show that the $ATI$s of most
memory access behaviors are negligible, which is consistent with the
observation from Fig. 3. However, there are some outliers whose access time
interval is larger than 0.8s and the corresponding memory block size is larger
than 600MB. For example, the $ATI$ and the corresponding memory block size of
the red marked outlier is 840211us and 1200MB, respectively. In this case, the
boundary data swapping size without sacrificing the runtime performance
between the host and device is about 2.54GB (by equation 1,
$S\leq\frac{0.8s}{1/6.4GB/s+1/6.3GB/s}=2.54GB$), which is much larger than
1200MB. This indicates that these outlier memory behaviors with high $ATI$s
and large memory block size ought to be the focus of attention. They are the
major contributors in terms of reducing the memory pressure of DNN training.
In the future, we plan to build an automatic cost model to sift out these
memory access behaviors to reduce the device memory pressure during training.

Figure 5:  Memory occupation breakdown of typical DNN training on the Nvidia
Titan X Pascal GPU.

(a) CIFAR-100 (32x32).

(b) CIFAR-100 (32x32).

Figure 6:  Memory ocupation breakdown of linear DNN (AlexNet) with different
batch size.

(a) ImageNet (224x224).

(b) ImageNet (224x224).

Figure 7:  Memory ocupation breakdown of non-linear DNN (ResNet) with
different layer structures (ResNet18, 34, 50, 101, 152).

Device Memory Occupation Breakdown. During training, the memory consumption
falls into three categories from the perspective of device memory storage
contents, i.e., input data, parameters and intermediate results [12]. We
observe from Fig. 5 that, for most DNNs, parameters only account for a small
fraction of the total memory footprint during training. This indicates that
weight pruning or quantization techniques are not efficient for reducing the
memory pressures of DNN training [13, 14]. In training, the intermediate
results are the primary contributor to the device memory footprint. To show
the detailed effects of datasets, batch size and model layer structures over
the memory footprint, we explored the enumeration space of these involving
factors. Due to space limit, we only report the data of one linear DNN
(AlexNet [15]) and one non-linear DNN (ResNet with different numbers of
residual layer blocks [16]) with CIFAR-100 and ImageNet dataset under
different batch size. As shown in Fig 6, for AlexNet, with the growth of batch
size, those intermediate results gradually dominate the device memory
consumption and the occupation of parameters is gradually weakened. Meanwhile,
the impact of input data increases slightly. This observation is also
applicable to the non-linear ResNet with different numbers of residual layer
blocks (Fig. 7).

##  IV Conclusion & Future Directions

From our characterizing results, we highlight the following observations. The
memory behaviors of DNN training illustrate obvious iterative patterns (Fig.
2). Besides, the distributions of the $ATI$s of most memory behaviors are
relatively concentrated. Especially, the $ATI$s of most memory access
behaviors are negligible, which can be filtered out by the bandwidth between
host and device memory (Fig. 3). Most importantly, those outlier memory
behaviors with high $ATI$ and large memory block size ought to be the focus of
attention (Fig. 4). They are the major contributors to reduce the memory
pressure of DNN training.

Based on these observations, and inspired by the swapping-based memory
pressure reduction schemes [9, 10], we plan to propose a more general approach
that takes the memory access patterns as input to automatically address the
device memory pressure issues of DNN training with small runtime overhead.

## References

  * [1] T. B. Brown, B. Mann, N. Ryder, M. Subbiah, J. Kaplan, P. Dhariwal, A. Neelakantan, P. Shyam, G. Sastry, A. Askell, S. Agarwal, A. Herbert-Voss, G. Krueger, T. Henighan, R. Child, A. Ramesh, D. M. Ziegler, J. Wu, C. Winter, C. Hesse, M. Chen, E. Sigler, M. Litwin, S. Gray, B. Chess, J. Clark, C. Berner, S. McCandlish, A. Radford, I. Sutskever, and D. Amodei, “Language Models are Few-Shot Learners,” 2020.
  * [2] W. Fedus, B. Zoph, and N. Shazeer, “Switch Transformers: Scaling to Trillion Parameter Models with Simple and Efficient Sparsity,” 2021\.
  * [3] OpenAI, “AI and Compute,” https://openai.com/blog/ai-and-compute, 2021\.
  * [4] NVIDIA, “Nvidia ampere architecture,” https://www.nvidia.com/en-us/data-center/nvidia-ampere-gpu-architecture/, 2020, [Online; accessed 20-Sept-2020].
  * [5] C. Szegedy, S. Ioffe, V. Vanhoucke, and A. A. Alemi, “Inception-v4, inception-resnet and the impact of residual connections on learning,” in _Proceedings of the Thirty-First AAAI Conference on Artificial Intelligence_ , ser. AAAI’17. AAAI Press, 2017, p. 4278–4284.
  * [6] D. Yang and D. Cheng, “Efficient gpu memory management for nonlinear dnns,” in _Proceedings of the 29th International Symposium on High-Performance Parallel and Distributed Computing_ , ser. HPDC ’20. New York, NY, USA: Association for Computing Machinery, 2020, p. 185–196. [Online]. Available: https://doi.org/10.1145/3369583.3392684
  * [7] S. Kokoska and C. Nevison, _The Binomial Cumulative Distribution Function_. New York, NY: Springer New York, 1989, pp. 50–51. [Online]. Available: https://doi.org/10.1007/978-1-4613-9629-1_5
  * [8] J. L. Hintze and R. D. Nelson, “Violin Plots: A Box Plot-Density Trace Synergism,” _The American Statistician_ , vol. 52, no. 2, pp. 181–184, 1998\. [Online]. Available: http://dx.doi.org/10.2307/2685478
  * [9] C.-C. Huang, G. Jin, and J. Li, “SwapAdvisor: Pushing Deep Learning Beyond the GPU Memory Limit via Smart Swapping,” in _Proceedings of the Twenty-Fifth International Conference on Architectural Support for Programming Languages and Operating Systems_ , ser. ASPLOS ’20. New York, NY, USA: Association for Computing Machinery, 2020, p. 1341–1355. [Online]. Available: https://doi.org/10.1145/3373376.3378530
  * [10] J. Ren, S. Rajbhandari, R. Y. Aminabadi, O. Ruwase, S. Yang, M. Zhang, D. Li, and Y. He, “ZeRO-Offload: Democratizing Billion-Scale Model Training,” 2021.
  * [11] NVIDIA., “CUDA Samples,” https://docs.nvidia.com/cuda/cuda-samples/index.html, 2020.
  * [12] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner, “Gradient-based learning applied to document recognition,” in _Proceedings of the IEEE_ , vol. 86, no. 11, 1998, pp. 2278–2324. [Online]. Available: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.42.7665
  * [13] M. Denil, B. Shakibi, L. Dinh, M. A. Ranzato, and N. de Freitas, “Predicting parameters in deep learning,” in _Advances in Neural Information Processing Systems 26_ , C. J. C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Q. Weinberger, Eds. Curran Associates, Inc., 2013, pp. 2148–2156. [Online]. Available: http://papers.nips.cc/paper/5025-predicting-parameters-in-deep-learning.pdf
  * [14] W. Niu, X. Ma, S. Lin, S. Wang, X. Qian, X. Lin, Y. Wang, and B. Ren, “Patdnn: Achieving real-time dnn execution on mobile devices with pattern-based weight pruning,” in _Proceedings of the Twenty-Fifth International Conference on Architectural Support for Programming Languages and Operating Systems_ , ser. ASPLOS ’20. New York, NY, USA: Association for Computing Machinery, 2020, p. 907–922. [Online]. Available: https://doi.org/10.1145/3373376.3378534
  * [15] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” in _Advances in Neural Information Processing Systems 25_ , F. Pereira, C. J. C. Burges, L. Bottou, and K. Q. Weinberger, Eds. Curran Associates, Inc., 2012, pp. 1097–1105. [Online]. Available: http://papers.nips.cc/paper/4824-imagenet-classification-with-deep-convolutional-neural-networks.pdf
  * [16] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2016, pp. 770–778.


# Variational Symplectic Accelerated Optimization on Lie Groups

Taeyoung Lee, Molei Tao, and Melvin Leok Taeyoung Lee, Mechanical and
Aerospace Engineering, George Washington University, Washington, DC 20052.
<EMAIL_ADDRESS>Tao, Mathematics, Georgia Institute of Technology, Atlanta,
GA 30332<EMAIL_ADDRESS>Leok, Mathematics, University of
California–San Diego, La Jolla, CA 92093<EMAIL_ADDRESS>research was
partially supported by NSF under the grants CNS-1837382, CMMI-1760928,
DMS-1813635, IIP-1747760, DMS-1847802, ECCS-1936776, and by AFOSR
FA9550-18-1-0288 and DoD 13106725.

###### Abstract

There has been significant interest in generalizations of the Nesterov
accelerated gradient descent algorithm due to its improved performance
guarantee compared to the standard gradient descent algorithm, and its
applicability to large scale optimization problems arising in deep learning. A
particularly fruitful approach is based on numerical discretizations of
differential equations that describe the continuous time limit of the Nesterov
algorithm, and a generalization involving time-dependent Bregman Lagrangian
and Hamiltonian dynamics that converges at an arbitrarily fast rate to the
minimum. We develop a Lie group variational discretization based on an
extended path space formulation of the Bregman Lagrangian on Lie groups, and
analyze its computational properties with two examples in attitude
determination and vision-based localization.

##  I Introduction

Nesterov’s accelerated gradient descent algorithm [1] was introduced in 1983,
and it exhibits the convergence rate of $\mathcal{O}(1/k^{2})$ when applied to
a convex objective function, which is faster than the $\mathcal{O}(1/k)$
convergence rate of standard gradient descent methods. It is shown in [2] that
this rate of convergence is optimal for the class of first-order gradient
methods. This improved rate of convergence over the standard gradient method
is referred to as acceleration, and there is a great interest in developing
systematic approaches to the construction of efficient accelerated
optimization algorithms, driven by potential applications in deep learning.

A continuous time limit of the Nesterov algorithm was studied in [3], whose
flow converges to the minimum at $\mathcal{O}(1/t^{2})$, and this was
generalized in [4] using a time-dependent Bregman Lagrangian and Hamiltonian
to obtain higher-order convergence of $\mathcal{O}(1/t^{p})$ for arbitrary
$p\geq 2$. However, it has been shown that discretizing Bregman dynamics is
not trivial as common discretizations fail to achieve the higher convergence
rate guaranteed in the continuous time limit. As such, there have been several
attempts to construct accelerated optimization algorithms using geometric
structure-preserving discretizations of the Bregman dynamics [5].

A natural class111Note that other classes of discretization methods exist,
such as those based on splitting (e.g., [6, 7]) and composition (e.g., [8]),
and such approaches also arise in variational discretization [9]. of geometric
numerical integrators [6] for discretizing such Lagrangian or Hamiltonian
systems is variational integrators [9, 10]. They are constructed by a discrete
analogue of Hamilton’s variational principle, and therefore, their numerical
flows are symplectic. They also satisfy a discrete Noether’s theorem that
relates symmetries with momentum conservation properties, and further exhibit
excellent exponentially long-time energy stability. One complication is that
such methods are typically developed for autonomous Lagrangian and Hamiltonian
systems on the Euclidean space. To address this, variational integrators have
been developed on a Lie group [11], and time-adaptive Hamiltonian variational
integrators have been proposed [12].

In this paper, we focus on the optimization problem to minimize an objective
function defined on an a Lie group. Optimization on a manifold or a Lie group
appears in various areas of machine learning, engineering, and applied
mathematics [13, 14], and respecting the geometric structure of manifolds
yields more accurate and efficient optimization schemes, when compared to
methods based on embeddings in a higher-dimensional Euclidean space with
algebraic constraints, or using local coordinates.

In particular, we formulate a Bregman Lagrangian system on a Lie group, and we
further discretize it using the extended Lie group variational integrator to
construct an intrinsic accelerated optimization scheme, which inherits the
desirable properties of variational integrators while also preserving the
group structure. Compared with [12] where the evolution of the stepsize is
prescribed, the proposed scheme adaptively adjusts the stepsize according to
the extended variational principle at the cost of increased computational
load. The resulting computational properties of the proposed approach are
analyzed with two examples in attitude determination and vision-based
localization, where it is observed that the scheme exhibits an interesting
convergence of the adaptive stepsize, and the variational discretization
provides robustness against the choice of stepsize, which is exploited in the
numerical experiments to improve computational efficiency. We also present
benchmark studies against other discretization schemes applied to the Bregman
dynamics, and other accelerated optimization schemes on a Lie group [7].

##  II Extended Lagrangian Mechanics

This section presents Lagrangian mechanics for non-autonomous systems on a Lie
group. It is referred to as extended Lagrangian mechanics as the variational
principle is extended to include reparamerization of time [9]. These are
developed in both of continuous-time and discrete-time formulations. The
latter yields a Lie group variational integrator [11], which will be applied
to accelerated optimization using the Bregman Lagrangian in the next section.

Consider an $n$-dimensional Lie group $\mathsf{G}$. Let $\mathfrak{g}$ be the
associated Lie algebra, or the tangent space at the identity, i.e.,
$\mathfrak{g}=\mathsf{T}_{e}\mathsf{G}$. Consider a left trivialization of the
tangent bundle of the group
$\mathsf{T}\mathsf{G}\simeq\mathsf{G}\times\mathfrak{g}$,
$(g,\dot{g})\mapsto(g,L_{g^{-1}}\dot{g})\equiv(g,\xi)$ More specifically, let
$\mathsf{L}:\mathsf{G}\times\mathsf{G}\rightarrow\mathsf{G}$ be the left
action defined such that $\mathsf{L}_{g}h=gh$ for $g,h\in\mathsf{G}$. Then the
left trivialization is a map
$(g,\dot{g})\mapsto(g,L_{g^{-1}}\dot{g})\equiv(g,\xi)$, where
$\xi\in\mathfrak{g}$, and the kinematics equation can be written as

$\displaystyle\dot{g}=g\xi.$ (1)

Further, suppose $\mathfrak{g}$ is equipped with an inner product
$\left\langle\cdot,\cdot\right\rangle$, which induces an inner product on
$\mathsf{T}_{g}\mathsf{G}$ via left trivialization. For any
$v,w\in\mathsf{T}_{g}\mathsf{G}$, $\left\langle
w,v\right\rangle_{\mathsf{T}_{g}\mathsf{G}}=\left\langle\mathsf{T}_{g}\mathsf{L}_{g^{-1}}v,\mathsf{T}_{g}\mathsf{L}_{g^{-1}}w\right\rangle_{\mathfrak{g}}$.
Given the inner product, we identify $\mathfrak{g}\simeq\mathfrak{g}^{*}$ and
$\mathsf{T}_{g}\mathsf{G}\simeq\mathsf{T}^{*}_{g}\mathsf{G}\simeq
G\times\mathfrak{g}^{*}$ via the Riesz representation. Throughout this paper,
the pairing is also denoted by the dot product $\cdot$. Let
$\mathbf{J}:\mathfrak{g}\rightarrow\mathfrak{g}^{*}$ be chosen such that
$\left\langle\mathbf{J}(\xi),\zeta\right\rangle$ is positive-definite and
symmetric as a bilinear form of $\xi,\zeta\in\mathfrak{g}$. Define the metric
$\langle\\!\langle\cdot,\cdot\rangle\\!\rangle:\mathfrak{g}\times\mathfrak{g}\rightarrow\mathbb{R}$
with
$\langle\\!\langle\xi,\zeta\rangle\\!\rangle=\left\langle\mathbf{J}(\xi),\zeta\right\rangle$.
This serves as a left-invariant Riemmanian metric on $\mathsf{G}$. Also
$\|\xi\|^{2}=\langle\\!\langle\xi,\xi\rangle\\!\rangle$ for any
$\xi\in\mathfrak{g}$. The adjoint operator is denoted by
$\mathrm{Ad}_{g}:\mathfrak{g}\rightarrow\mathfrak{g}$, and the ad operator is
denoted by $\mathrm{ad}_{\xi}:\mathfrak{g}\rightarrow\mathfrak{g}$. See, for
example [15] for detailed preliminaries.

###  II-A Continuous-Time Extended Lagrangian Mechanics

Consider a non-autonomous (left-trivialized) Lagrangian
$L(t,g,\xi):\mathbb{R}\times\mathsf{G}\times\mathfrak{g}\rightarrow\mathbb{R}$
on the extended state space. The corresponding extended path space is composed
of the curves $(c_{t}(a),c_{g}(a))$ on $\mathbb{R}\times\mathsf{G}$
parameterized by $a>0$. To ensure that the reparameterized time increases
monotonically, we require $c^{\prime}_{t}(a)>0$. For a given time interval
$[t_{0},t_{f}]$, the corresponding interval $[a_{0},a_{f}]$ for $a$ is chosen
such that $t_{0}=c_{t}(a_{0})$ and $t_{f}=c_{t}(a_{f})$. For any path
$(c_{t}(a),c_{g}(a))$ over $[a_{0},a_{f}]$ in the extended space, the
associated curve is

$\displaystyle g(t)=c_{g}(c_{t}^{-1}(t)),$ (2)

on $\mathsf{G}$ over the time interval $[t_{0},t_{f}]$. For a given extended
path, define the extended action integral as

$\displaystyle\mathfrak{G}(c_{t},c_{g})=\int_{t_{0}}^{t_{f}}L(t,g,\xi)\bigg{|}_{g(t)=c_{g}(c_{t}^{-1}(t))}dt,$
(3)

where the Lagrangian is evaluated on the associated curve (2), and $\xi$
satisfies the kinematics equation (1).

Taking the variation of $\mathfrak{G}$ with respect to the extended path, we
obtain the Euler–Lagrange equation according to the variational principle in
the extended phase space. As discussed in [9, Sec. 4.2.2], the resulting
Euler–Lagrange equations depends only on the associated curve (2), not on the
extended path $(c_{t},c_{g})$ itself, and the variational principle does not
dictate how the curve should be reparameterized.

Further, the resulting Euler–Lagrange equation share the exactly same form as
(unextended) Lagrangian mechanics for the associated curve. As such, the
Euler–Lagrange equation for non-autonomous Lagrangian
$L(t,g,\xi):\mathbb{R}\times\mathsf{G}\times\mathfrak{g}\rightarrow\mathbb{R}$
can be written as

$\displaystyle\frac{d}{dt}\\!\left(\frac{\partial
L}{\partial\xi}\right)-\mathrm{ad}^{*}_{\xi}\frac{\partial
L}{\partial\xi}-\mathsf{T}^{*}_{e}\mathsf{L}_{g}(\mathbf{D}_{g}L)=0,$ (4)

where $\mathbf{D}_{g}$ stands for the differential with respect to $g$ (see
[16, Sec. 8.6.3] for derivation of the above equation for autonomous
Lagrangians).

Introducing the Legendre transform $\mu=\frac{\partial
L}{\partial\xi}\in\mathfrak{g}^{*}$, and assuming that it is invertible, the
Euler–Lagrange equation can be rewritten as

$\displaystyle\dot{\mu}-\mathrm{ad}^{*}_{\xi}\mu-\mathsf{T}^{*}_{e}\mathsf{L}_{g}(\mathbf{D}_{g}L)=0.$
(5)

###  II-B Extended Lie Group Variational Integrator

Variational integrators are geometric numerical integration schemes that can
be viewed as discrete-time mechanics derived from a discretization of the
variational principle for Lagrangian mechanics [9]. The discrete-time flows of
variational integrators are symplectic and they exhibit a discrete analogue of
Noether’s theorem. This provides long-term structural stability in the
resulting numerical simulations. For Lagrangian mechanics evolving on a Lie
group, the corresponding Lie group variational integrators were developed in
[11].

Here, we develop extended Lie group variational integrators by discretizing
the extended variational principle presented above, following the general
framework of [9]. The extended discrete path space is composed of the sequence
$\\{(t_{k},g_{k})\\}_{k=0}^{N}$ on $\mathbb{R}\times\mathsf{G}$, satisfying
$t_{k+1}>t_{k}$. Next, the discrete kinematics equation is chosen to be

$\displaystyle g_{k+1}=g_{k}f_{k},$ (6)

for $f_{k}\in\mathsf{G}$ representing the relative update over a single
timestep. The discrete Lagrangian
$L_{d}(t_{k},t_{k+1},g_{k},f_{k}):\mathbb{R}\times\mathbb{R}\times\mathsf{G}\times\mathsf{G}\rightarrow\mathbb{R}$
is chosen such that the following extended discrete action sum

$\displaystyle\mathfrak{G}_{d}(\\{(t_{k},g_{k})\\}_{k=0}^{N})=\sum_{k=0}^{N-1}L_{d}(t_{k},t_{k+1},g_{k},f_{k}),$
(7)

approximates (3).

###### Proposition 1

The discrete path $\\{(g_{k},f_{k})\\}_{k=0}^{N-1}$ that extremizes the
discrete action sum (7) subject to fixed endpoints satisfies the following
discrete Euler–Lagrange equation,

$\displaystyle\mathsf{T}^{*}_{e}\mathsf{L}_{g_{k}}(\mathbf{D}_{g_{k}}L_{d_{k}})-\mathrm{Ad}^{*}_{f_{k}^{-1}}(\mathsf{T}^{*}_{e}\mathsf{L}_{f_{k}}(\mathbf{D}_{f_{k}}L_{d_{k}}))$
$\displaystyle+\mathsf{T}^{*}_{e}\mathsf{L}_{f_{k-1}}(\mathbf{D}_{f_{k-1}}L_{d_{k-1}})=0,$
(8)
$\displaystyle\mathbf{D}_{t_{k}}L_{d_{k-1}}+\mathbf{D}_{t_{k}}L_{d_{k}}=0,$
(9)

which together with the discrete kinematic equation (6) defines an extended
Lie group variational integrator.

###### Proof:

From (6), $\delta f_{k}=-g_{k}^{-1}(\delta
g_{k})g_{k}^{-1}g_{k+1}+g_{k}^{-1}\delta g_{k+1}$. Since $\delta g_{k}$ can be
written as $\delta g_{k}=g_{k}\eta_{k}$ for $\eta_{k}\in\mathfrak{g}$,

$\displaystyle f_{k}^{-1}\delta
f_{k}=-\mathrm{Ad}_{f_{k}^{-1}}\eta_{k}+\eta_{k+1}.$ (10)

Take the variation of (7) and substitute (10) to obtain

$\displaystyle\delta\mathfrak{G}_{d}=\sum_{k=0}^{N-1}$
$\displaystyle\mathsf{T}^{*}_{e}\mathsf{L}_{g_{k}}(\mathbf{D}_{g_{k}}L_{d_{k}})\cdot\eta_{k}$
$\displaystyle+\mathsf{T}^{*}_{e}\mathsf{L}_{f_{k}}(\mathbf{D}_{f_{k}}L_{d_{k}})\cdot(-\mathrm{Ad}_{f_{k}^{-1}}\eta_{k}+\eta_{k+1})$
$\displaystyle+\mathbf{D}_{t_{k}}L_{d_{k}}\cdot\delta
t_{k}+\mathbf{D}_{t_{k+1}}\mathbf{D}_{d_{k}}\cdot\delta t_{k+1}.$

Since the endpoints are fixed, we have $\eta_{0}=0$ and $\delta t_{0}=0$.
Therefore in the above expression, the range of summation for the terms paired
with $\eta_{k}$ and $\delta t_{k}$ can be reduced to $1\leq k\leq N-1$. Also,
using $\eta_{N}=0$ and $\delta t_{N}=0$, for the other terms paired with
$\eta_{k+1}$ and $\delta t_{k+1}$, the terms can be reindexed by reducing the
subscripts by one and summed over the same range. According to the variational
principle, $\delta\mathfrak{G}_{d}=0$ for any $\eta_{k}$ and $\delta t_{k}$,
which yields (8) and (9). ∎

The most notable difference compared to the continuous-time counterpart is
that in addition to the discrete Euler–Lagrange equation (8), we have the
additional equation (9) for the evolution of the discrete time. This is
because the discrete action sum $\mathfrak{G}_{d}$ depends on the complete
extended path $\\{(t_{k},g_{k})\\}_{k=1}^{N}$. Whereas the continuous-time
action $\mathfrak{G}$ is only a function of the associated curve (2).

The discrete Euler–Lagrange equation for the discrete time (9) is associated
with the energy. Define the discrete energy to be

$\displaystyle E^{+}_{k}$ $\displaystyle=-\mathbf{D}_{t_{k+1}}L_{d_{k}},$ (11)
$\displaystyle E^{-}_{k}$ $\displaystyle=\mathbf{D}_{t_{k}}L_{d_{k}}.$ (12)

Then, (9) can be rewritten as

$\displaystyle E^{+}_{k-1}=E^{-}_{k},$ (13)

which reflects the evolution of the discrete energy. When the discrete
Lagrangian is autonomous, (13) implies the conservation of discrete energy,
thereby yielding a symplectic-energy-momentum integrator [17].

To implement (8) and (9) as a numerical integrator, it is more convenient to
introduce the extended discrete Legendre transforms,
$\mathbb{F}^{\pm}L_{d_{k}}:\mathbb{R}\times\mathbb{R}\times\mathsf{G}\times\mathsf{G}\rightarrow\mathbb{R}\times\mathbb{R}\times\mathsf{G}\times\mathfrak{g}^{*}$
as

$\displaystyle\mathbb{F}^{+}L_{d_{k}}(t_{k},t_{k+1},g_{k},f_{k})$
$\displaystyle=(t_{k+1},E_{k+1},g_{k+1},\mu_{k+1}),$ (14)
$\displaystyle\mathbb{F}^{-}L_{d_{k}}(t_{k},t_{k+1},g_{k},f_{k})$
$\displaystyle=(t_{k},E_{k},g_{k},\mu_{k}).$ (15)

where

$\displaystyle\mu_{k}$
$\displaystyle=-\mathsf{T}^{*}_{e}\mathsf{L}_{g_{k}}(\mathbf{D}_{g_{k}}L_{d_{k}})+\mathrm{Ad}^{*}_{f_{k}^{-1}}(\mathsf{T}^{*}_{e}\mathsf{L}_{f_{k}}(\mathbf{D}_{f_{k}}L_{d_{k}})),$
(16) $\displaystyle\mu_{k+1}$
$\displaystyle=\mathsf{T}^{*}_{e}\mathsf{L}_{f_{k}}(\mathbf{D}_{f_{k}}L_{d_{k}}),$
(17)

and $E_{k+1}$ and $E_{k}$ are given by (11) and (12), respectively.

The resulting discrete flow map is defined by
$\mathbb{F}^{+}L_{d_{k}}\circ(\mathbb{F}L_{d_{k}})^{-1}$. More specifically,
for given $(t_{k},E_{k},g_{k},\mu_{k})$, (12) and (16) are solved together for
$t_{k+1},f_{k}$ with the constraint $t_{k+1}>t_{k}$. Then,
$(E_{k+1},g_{k+1},\mu_{k+1})$ are computed by (11), (6), and (17),
respectively. This yields the discrete flow map
$(t_{k},E_{k},g_{k},\mu_{k})\rightarrow(t_{k+1},E_{k+1},g_{k+1},\mu_{k+1})$
consistent with (8) and (9). While the flow map is expressed in terms of $E$
for convenience, the initial value of $E_{0}$ is often selected by choosing
the initial timestep $h_{0}$ and calculating the corresponding value of
$E_{0}$ through (12). This inherits the desirable properties of variational
integrators, and the group structure is also preserved through (6).

##  III Bregman Lagrangian Systems on $\mathsf{G}$

Let $\mathsf{f}:\mathsf{G}\rightarrow\mathbb{R}$ be a real-valued smooth
function on $\mathsf{G}$. We focus on the optimization problem:

$\displaystyle\min_{g\in\mathsf{G}}\mathsf{f}(g).$ (18)

A variational accelerated optimization scheme for the above problem was
developed in [7], where the Nesterov accelerated gradient (NAG) descent on a
finite-dimensional vector space was intrinsically generalized to a Lie group.
In this section, we introduce an intrinsic formulation of Bregman Lagrangian
dynamics [4], which encompasses a larger class of accelerated optimization
scheme, including NAG. More importantly, the continuous dynamics guarantees
polynomial convergence rates up to an arbitrary order.

###  III-A Continuous-Time Bregman Dynamics

The Bregman Lagrangian
$L(t,g,\xi):\mathbb{R}\times\mathsf{G}\times\mathfrak{g}\rightarrow$ is

$\displaystyle L(t,g,\xi)=\frac{t^{\lambda
p+1}}{2p}\|\xi\|^{2}-Cpt^{(\lambda+1)p-1}\mathsf{f}(g),$ (19)

where
$\|\xi\|^{2}=\langle\\!\langle\xi,\xi\rangle\\!\rangle=\left\langle\mathbf{J}(\xi),\xi\right\rangle$,
for $p,C>0$, and $\lambda\geq 1$. When $\mathsf{G}=\mathbb{R}^{n}$ and
$\lambda=1$, this recovers the Bregman Lagrangian for vector spaces [4], and
it yields the continuous-time limit of Nesterov’s accelerated gradient descent
for $p=2$ [18]. Also, in case $p=3$, it corresponds to the continuous-time
limit of Nesterov’s accelerated cubic-regularized Newton’s method [19]. When
$\mathsf{G}$ is considered as a Riemannian manifold, this corresponds to the
$p$-Bregman Lagrangian in [20]. The additional term $\lambda$ accounts for the
sectional curvature and diameter of the manifold [21].

The left-trivialized derivative of the objective function is

$\displaystyle\nabla_{\mathsf{L}}\mathsf{f}(g)=\mathsf{T}_{e}^{*}\mathsf{L}_{g}(\mathbf{D}_{g}\mathsf{f}(g)).$
(20)

Applying (4) to (19), the corresponding Euler–Lagrange equations are given
below.

###### Proposition 2

The Euler–Lagrange equations corresponding to the Bregman Lagrangian (19) are

$\displaystyle\frac{d\mathbf{J}(\xi)}{dt}+\frac{\lambda
p+1}{t}\mathbf{J}(\xi)-\mathrm{ad}^{*}_{\xi}\mathbf{J}(\xi)+Cp^{2}t^{p-2}\nabla_{\mathsf{L}}\mathsf{f}(g)=0,$
(21)

and (1). Further, the corresponding continuous flow locally converges to the
minimizer $g^{*}$ of $\mathsf{f}$ with the rate given by

$\displaystyle\mathsf{f}(g(t))-\mathsf{f}(g^{*})\in\mathcal{O}(t^{-p}),$ (22)

when $\mathsf{f}$ is geodesically convex.

###### Proof:

We have

$\displaystyle\frac{\partial L}{\partial\xi}$ $\displaystyle=\frac{t^{\lambda
p+1}}{p}\mathbf{J}(\xi)$

Substituting this into (4) and using (20),

$\displaystyle\frac{t^{\lambda
p+1}}{p}\frac{d\mathbf{J}(\xi)}{dt}+\frac{(\lambda p+1)t^{\lambda
p}}{p}\mathbf{J}(\xi)-\frac{t^{\lambda
p+1}}{p}\mathrm{ad}^{*}_{\xi}\mathbf{J}(\xi)$
$\displaystyle+Cpt^{(\lambda+1)p-1}\nabla_{\mathsf{L}}\mathsf{f}(g)=0.$

Dividing both sides by $\frac{t^{\lambda p+1}}{p}$ yields (21). The
convergence property is established by [20, Theorem 3.2]. ∎

Therefore, the optimization problem on $\mathsf{G}$ can be addressed by
numerically integrating (21) from an initial guess. However, it has been
observed that a naïve discretization is not able to match the polynomial
convergence rate established in [4]. Further, we need a guarantee that the
discrete trajectory evolves on the Lie group.

These two challenges can be addressed by applying a Lie group variational
integrator, as their structure-preserving properties provides long-term
numerical stability, and preservation of the group structure. In the
subsequent section, we derive Lie group variational integrators for the
Bregman Lagrangian system.

###  III-B Lie Group Variational Integrator for Bregman Dynamics

Let $h_{k}=t_{k+1}-t_{k}$ and $t_{k,k+1}=t_{k}+h_{k}/2$. We consider the
following form of the discrete Lagrangian

$\displaystyle L_{d}(t_{k},t_{k+1},g_{k},f_{k})$
$\displaystyle=\frac{\phi(t_{k,k+1})}{h_{k}}T_{d}(f_{k})-\frac{h_{k}}{2}\theta(t_{k})\mathsf{f}(g_{k})$
$\displaystyle\quad-\frac{h_{k}}{2}\theta(t_{k+1})\mathsf{f}(g_{k}f_{k}),$
(23)

where $T_{d}(f_{k}):\mathsf{G}\rightarrow\mathbb{R}$ is chosen such that it
approximates $T(f_{k})\approx h_{k}^{2}\|\xi_{k}\|^{2}/2$, and
$\phi,\theta:\mathbb{R}\rightarrow\mathbb{R}$ are

$\displaystyle\phi(t)$ $\displaystyle=\frac{t^{\lambda p+1}}{p},$ (24)
$\displaystyle\theta(t)$ $\displaystyle=Cpt^{(\lambda+1)p-1}.$ (25)

The corresponding variational integrators are presented as follows.

###### Proposition 3

The discrete-time Euler–Lagrange equations, or the Lie group variational
integrator for the discrete Lagrangian (23) corresponding to the Bregman
Lagrangian (19) are given by

$\displaystyle\mu_{k}$
$\displaystyle=\frac{\phi_{k,k+1}}{h_{k}}\mathrm{Ad}^{*}_{f_{k}^{-1}}(\mathsf{T}^{*}_{e}\mathsf{L}_{f_{k}}(\mathbf{D}_{f_{k}}T_{d_{k}}))+\frac{h_{k}\theta_{k}}{2}\nabla_{L}\mathsf{f}_{k},$
(26) $\displaystyle\mu_{k+1}$
$\displaystyle=\mathrm{Ad}^{*}_{f_{k}}(\mu_{k}-\frac{h_{k}\theta_{k}}{2}\nabla_{L}\mathsf{f}_{k})-\frac{h_{k}\theta_{k+1}}{2}\nabla_{\mathsf{L}}\mathsf{f}_{k+1},$
(27) $\displaystyle E_{k}$
$\displaystyle=\frac{\phi^{\prime}_{k,k+1}}{2h_{k}}T_{d_{k}}-\frac{h_{k}\theta^{\prime}_{k}}{2}\mathsf{f}_{k}$
$\displaystyle\quad+\frac{\phi_{k,k+1}}{h_{k}^{2}}T_{d_{k}}+\frac{\theta_{k}}{2}\mathsf{f}_{k}+\frac{\theta_{k+1}}{2}\mathsf{f}_{k+1},$
(28) $\displaystyle E_{k+1}$
$\displaystyle=-\frac{\phi^{\prime}_{k,k+1}}{2h_{k}}T_{d_{k}}+\frac{h_{k}\theta^{\prime}_{k+1}}{2}\mathsf{f}_{k+1}$
$\displaystyle\quad+\frac{\phi_{k,k+1}}{h_{k}^{2}}T_{d_{k}}+\frac{\theta_{k}}{2}\mathsf{f}_{k}+\frac{\theta_{k+1}}{2}\mathsf{f}_{k+1},$
(29)

together with (6).

###### Proof:

These can be derived by substituting (23) into (16), (17), (12), and (11),
respectively. ∎

As discussed at the end of Section III, these provide symplectic and momentum-
preserving discrete time flow maps. Since these corresponds to a
discretization of the Bregman Lagrangian system, they can be considered as a
geometric numerical integrator for (21), or utilized as an optimization
algorithm on $\mathsf{G}$. If $T_{d}(f_{k})=T_{d}(f_{k}^{-1})$, then the
discrete Lagrangian is self-adjoint, and the above integrator is symmetric and
therefore at least second-order accurate.

##  IV Optimization on $\mathsf{G}$

In this section, we present both of the continuous Bregman Lagrangian system
and the Lie group variational integrator for several Lie groups.

###  IV-A Euclidean Space $\mathbb{R}^{n}$

Suppose $\mathsf{G}=\mathbb{R}^{n}$, with the additive group action, and the
inner product is chosen to be $\left\langle x,y\right\rangle=x^{T}y$ for any
$x,y\in\mathbb{R}^{n}$. Let $\mathbf{J}(\dot{x})=I_{n\times n}\dot{x}$, and
$\lambda=1$.

From (21), the continuous Euler–Lagrange equation is given by

$\displaystyle\ddot{x}+\frac{p+1}{t}\dot{x}+Cp^{2}t^{p-2}\nabla\mathsf{f}(x)=0,$
(30)

which recovers the differential equation derived in [4].

Next, we develop variational integrators. The discrete kinematics equation (6)
is rewritten as $x_{k+1}=x_{k}+\Delta x_{k}$ for $\Delta
x_{k}\in\mathbb{R}^{n}$. The kinetic energy term in (23) is chosen as

$\displaystyle T_{d}=\frac{1}{2}\|\Delta x_{k}\|^{2}.$ (31)

According to 3, we obtain the discrete Euler–Lagrange equations as follows.

###### Proposition 4

When $G=\mathbb{R}^{n}$, the variational integrator for the discrete Bregman
Lagrangian (23) is given by

$\displaystyle v_{k}$ $\displaystyle=\frac{\phi_{k,k+1}}{h_{k}}\Delta
x_{k}+\frac{h_{k}\theta_{k}}{2}\nabla\mathsf{f}_{k},$ (32) $\displaystyle
v_{k+1}$
$\displaystyle=v_{k}-\frac{h_{k}\theta_{k}}{2}\nabla\mathsf{f}_{k}-\frac{h_{k}\theta_{k+1}}{2}\nabla\mathsf{f}_{k+1},$
(33)

and (28), (29) with (31).

These are implicit as (32) and (28) should be solved together for $\Delta
x_{k}$ and $h_{k}$. One straightforward approach is fixed-point iteration. For
a given $h_{k}$, (32) can be solved explicitly for $\Delta_{k}$, which yields
$x_{k+1}$. Then, (28) can be solved for $h_{k}$. These procedure are iterated
until $h_{k}$ converges.

###  IV-B Three-Dimensional Special Orthogonal Group $\mathsf{SO(3)}$

Next, consider $\mathsf{SO(3)}=\\{R\in\mathbb{R}^{3\times
3}\,|\,R^{T}R=I_{3\times 3},\,\mathrm{det}(R)]=1\\}$. Its Lie algebra is
$\mathfrak{so}(3)=\\{S\in\mathbb{R}^{3\times 3}\,|\,S^{T}=-S\\}$ with the
matrix commutator as the Lie bracket. This is identified with $\mathbb{R}^{3}$
through the hat map $\hat{\cdot}:\mathbb{R}^{3}\rightarrow\mathfrak{so}(3)$
defined such that $\hat{x}\in\mathfrak{so}(3)$ and $\hat{x}y=x\times y$ for
any $x,y\in\mathbb{R}^{3}$. The inverse of the hat map is denoted by the vee
map $\vee:\mathfrak{so}(3)\rightarrow\mathbb{R}^{3}$. The inner product is
given by

$\displaystyle\langle\hat{\eta},\hat{\xi}\rangle_{\mathfrak{so}(3)}=\frac{1}{2}\mbox{tr}\negthickspace\left[\hat{\eta}^{T}\hat{\xi}\right]=\eta^{T}\xi=\left\langle\eta,\xi\right\rangle_{\mathbb{R}^{3}}.$

The metric is chosen as

$\displaystyle\langle\mathbf{J}(\hat{\eta}),\hat{\xi}\rangle_{\mathfrak{so}(3)}=\mbox{tr}\negthickspace\left[\hat{\eta}^{T}J_{d}\hat{\xi}\right]=\eta^{T}J\xi=\left\langle
J\eta,\xi\right\rangle_{\mathbb{R}^{3}},$ (34)

where $J\in\mathbb{R}^{3\times 3}$ is a symmetric, positive-definite matrix,
and $J_{d}=\frac{1}{2}\mbox{tr}\negthickspace\left[J\right]I_{3\times
3}-J\in\mathbb{R}^{3\times 3}$. Further,

$\displaystyle\mathrm{ad}_{\eta}\xi$ $\displaystyle=\eta\times\xi,$
$\displaystyle\quad\mathrm{ad}^{*}_{\eta}\xi$ $\displaystyle=\xi\times\eta,$
$\displaystyle\mathrm{Ad}_{F}\eta$ $\displaystyle=F\eta,$
$\displaystyle\quad\mathrm{Ad}^{*}_{F}\eta$ $\displaystyle=F^{T}\eta.$

Consider

$\displaystyle L(t,R,\Omega)=\frac{t^{p+1}}{2p}\Omega\cdot J\Omega-
Cpt^{2p-1}\mathsf{f}(R).$

From (21), the Euler–Lagrange equations are given by

$\displaystyle
J\dot{\Omega}+\frac{p+1}{t}J\Omega+\hat{\Omega}J\Omega+Cp^{2}t^{p-2}\nabla_{\mathsf{L}}\mathsf{f}(R)=0,$
(35) $\displaystyle\dot{R}=R\hat{\Omega}.$ (36)

Next, we derive variational integrators. The kinematics equation is written as

$\displaystyle R_{k+1}=R_{k}F_{k},$ (37)

for $F_{k}\in\mathsf{SO(3)}$. Similar with [11], the angular velocity is
approximated with
$\hat{\Omega}_{k}\approx\frac{1}{h_{k}}R_{k}^{T}(R_{k+1}-R_{k})=\frac{1}{h_{k}}(F_{k}-I_{3\times
3})$. Substituting this into (34),

$\displaystyle T_{d}(F_{k})$
$\displaystyle=\mbox{tr}\negthickspace\left[(I_{3\times
3}-F_{k})J_{d}\right],$ (38)

which satisfies $T_{d}(F_{k})=T_{d}(F_{k}^{T})$.

###### Proposition 5

When $\mathsf{G}=\mathsf{SO(3)}$, the Lie group variational integrator for the
discrete Bregman Lagrangian (23) with (38) is given by

$\displaystyle\mu_{k}$
$\displaystyle=\frac{\phi_{k,k+1}}{h_{k}}(F_{k}J_{d}-J_{d}F_{k}^{T})^{\vee}+\frac{h_{k}\theta_{k}}{2}\nabla_{\mathsf{L}}\mathsf{f}_{k},$
(39) $\displaystyle\mu_{k+1}$
$\displaystyle=F_{k}^{T}\mu_{k}-\frac{h_{k}\theta_{k}}{2}\nabla_{\mathsf{L}}\mathsf{f}_{k}-\frac{h_{k}\theta_{k}}{2}\nabla_{\mathsf{L}}\mathsf{f}_{k+1},$
(40)

together with (29), (37), (28), and (38).

###### Proof:

Let $\delta F_{k}=F_{k}\hat{\chi}_{k}$. The derivative of (38) is

$\displaystyle\mathbf{D}_{F_{k}}T_{d_{k}}\cdot\delta
F_{k}=\mbox{tr}\negthickspace\left[-F_{k}\hat{\chi}_{k}J_{d}\right]=(J_{d}F_{k}-F_{k}^{T}J_{d})^{\vee}\cdot\chi,$

where the last equality is from the identity,
$\mathrm{tr}[-\hat{x}A]=x\cdot(A-A^{T})^{\vee}$ for any $x\in\mathbb{R}^{3}$
and $A\in\mathbb{R}^{3\times 3}$. Thus,
$\mathsf{T}^{*}_{I}\mathsf{L}_{F_{k}}(\mathbf{D}_{F_{k}}T_{d_{k}})=(J_{d}F_{k}-F_{k}^{T}J_{d})^{\vee}$.
Substituting this into (26) and (27) yields (39) and (40), respectively. ∎

To implement these, (40) and (28) should be solved together for $h_{k}$ and
$F_{k}$. For a given $h_{k}$, computational approaches to solve (39) for
$F_{k}$ are presented in [22, Sec 3.3.8]. When $J=I_{3\times 3}$, or
equivalently when $J_{d}=\frac{1}{2}I_{3\times 3}$, (39) can be solved
explicitly to obtain

$\displaystyle F_{k}=\exp\left(\frac{\sin^{-1}\|a\|}{\|a\|}\hat{a}\right),$
(41)

where
$a=\frac{h_{k}}{\phi_{k,k+1}}(\mu_{k}-\frac{h_{k}\theta_{k}}{2}\nabla_{\mathsf{L}}\mathsf{f}_{k})\in\mathbb{R}^{3}$.
This can replace (39).

###  IV-C Product of $\mathbb{R}^{n}$ and $\mathsf{SO(3)}$

Suppose $\mathsf{G}=\mathsf{SO(3)}\times\mathbb{R}^{n}$. As it is the direct
product of $\mathsf{SO(3)}$ and $\mathbb{R}^{n}$, the variation of the action
sum is decomposed into two parts of $\mathsf{SO(3)}$ and $\mathbb{R}^{n}$.
Therefore, the continuous Euler–Lagrange equations on
$\mathsf{SO(3)}\times\mathbb{R}^{n}$ are given by (30) and (35), after
replacing $\nabla\mathsf{f}(x)$ of (30) with $\nabla_{x}\mathsf{f}(R,x)$, and
replacing $\nabla_{\mathsf{L}}f(R)$ of (35) with
$\mathsf{T}^{*}_{I}\mathsf{L}_{R}(\mathbf{D}_{R}\mathsf{f}(R,x))$.

Similarly, the corresponding Lie group variational integrators are also given
by (32), (33), (39), and (40), in addition to the energy equations (28) and
(29) with

$\displaystyle T_{d_{k}}(F_{k},\Delta x_{k})=\frac{1}{2}\|\Delta
x_{k}\|^{2}+\mbox{tr}\negthickspace\left[(I_{3\times 3}-F_{k})J_{d}\right].$

##  V Numerical Examples

###  V-A Optimization on $\mathsf{SO(3)}$

Consider the objective function given by

$\displaystyle\mathsf{f}(R)$
$\displaystyle=\frac{1}{2}\|A-R\|^{2}_{\mathcal{F}}=\frac{1}{2}(\|A\|^{2}_{\mathcal{F}}+3)-\mbox{tr}\negthickspace\left[A^{T}R\right],$
(42)

where $\|\cdot\|_{\mathcal{F}}$ denotes the Frobenius norm, and
$A\in\mathbb{R}^{3\times 3}$. Optimization of the above function appears in
the least-squares estimation of attitude, referred to as Wahba’s problem [23].
Let the singular value decomposition of $A=USV^{T}$ for a diagonal
$S\in\mathbb{R}^{3\times 3}$ and $U,V\in\mathsf{O}(3)$. The optimal attitude
is explicitly given by $R^{*}=U\mathrm{diag}[1,1,\mathrm{det}(UV)]V^{T}$. The
left-trivialized gradient is
$\nabla_{\mathsf{L}}\mathsf{f}(R)=(A^{T}R-R^{T}A)^{\vee}$.

####  V-A1 Order of Convergence

(a) convergence with respect to $t$

(b) convergence with respect to $k$

Figure 1: Convergence rate of LGVI in 5 for varying $p$

First, we check if the theoretical order of convergence guaranteed by 2 is
achieved by the discrete Euler–Lagrange equations presented in 3. The elements
of the matrix $A$ in (42) are randomly chosen from the uniform distribution on
$[0,1]$. The initial guess of $R_{0}$ is chosen such that the initial error is
$0.9\pi$ in terms of the Euler-axis rotation. Lie group variational
integrators (LGVI) in 5 are simulated with fixed $J=I_{3\times 3}$, $C=1$, and
$h_{0}=0.1$ for varying $p\in\\{2,4,6,8\\}$. Since $J=I_{3\times 3}$, (39) is
replaced by (41). The remaining implicit equation (28) is solved for $h_{k}$
via the Matlab equation solver, lsqnonlin with the tolerance of $10^{-4}$. The
initial guess for $h_{k}$ is provided by $h_{k-1}$.

The resulting convergence rate represented by $\mathsf{f}-\mathsf{f}^{*}$ over
$t_{k}$ is illustrated in Figure 1.(a), where the empirical convergence rate
computed by manual fitting are also marked. It is shown that LGVI empirically
achieved the order of convergence greater than the theoretical guarantee of
$\mathcal{O}(t^{-p})$. It has been reported that naïve discretizations of
Bregman Lagrangian systems are not able to match the theoretical convergence
rate, or it might cause numerical instability [4, 5]. These results suggest
that LGVIs do not suffer from these discretization issues, and their
performance are consistent with the continuous-time analysis.

Next, given that the step size $h_{k}$ is adjusted adaptively according to
(28) and (29), it is likely that numerical simulation with higher $p$ requires
a smaller step size. In fact, the average step sizes are given by $6.15\times
10^{-2}$, $6.50\times 10^{-3}$, $4.89\times 10^{-4}$ and $1.21\times 10^{-5}$,
respectively for $p\in\\{2,4,6,8\\}$. To examine the effects of the step size
variations, the convergence with respect to the discrete time step is
illustrated in Figure 1.(b). It turns out that all of four cases of $p$
exhibit the similar order of long-term convergence, approximately
$\mathcal{O}(k^{-2.3})$. This is not surprising, as Nesterov [2] showed that
for every smooth first-order method, there exists a convex, $L$-smooth
objective function, such that the rate of convergence is bounded from below by
$\mathcal{O}(k^{-2})$, but it does not preclude the possibility of faster
rates of convergence for strongly convex functions.

However, the case of higher $p$ benefits from faster initial convergence, and
as a result, the terminal error for $p=4$ is more than 400 times smaller than
that of $p=2$.

####  V-A2 Effects of Initial Step Size

(a) convergence with respect to $t$

(b) evolution of step size $h_{k}$

Figure 2: Convergence rate of LGVI in 5 for varying $h_{0}$.

As discussed at the end of Section III, the extended LGVI requires choosing
the initial step size $h_{0}$. Here, we study the effects of $h_{0}$ in the
convergence. More specifically, the order is fixed to $p=4$, and the initial
step size is varied as $h_{0}\in\\{0.001,0.05,0.01,0.1,0.4\\}$. The
corresponding results are illustrated at Figure 2. Interestingly, in Figure
2.(a), the convergence with respect to $t$ is not much affected by the initial
step size $h_{0}$. Next, Figure 2.(b) presents the time-evolution of the step
size, and it is shown that the step size computed by (28) decreases at the
approximate order of $\mathcal{O}(t^{-1.6})$ for all cases. This might have
been caused by the fact that the forcing term in (35) increases over time.
Another notable feature is that after a certain period, the step sizes tend to
converge. More specifically, the step size initialized by $h_{0}=0.001$
converges to $1.8\times 10^{-4}$ when $t>10$, which is joined by the case of
$h_{0}=0.005$ later. It is expected that the next case for $h_{0}=0.01$ would
follow the similar trend if the simulation time is increased. This implies a
certain stability property of the extended LGVI in the step size. Furthermore,
observe that for the wide range of variations of step sizes presented in
Figure 2.(b), the convergence in Figure 2.(a) is fairly consistent, which
suggests that the LGVI is robust to the choice of the step size.

####  V-A3 Comparison with Other Discretizations of Bregman Euler–Lagrange
Equation

(a) convergence with respect to $t$

(b) orthogonality error of $R_{k}$

Figure 3: Comparison with other discretization schemes for Bregman
Euler–Lagrange equation

Next, we compare LGVI with other discretization schemes applied to (35) and
(36). Three methods are considered, namely the splitting approach introduced
in [7] applied to the proposed continuous dynamics (abbreviated as SPLT), a
4-th order fixed-step Runge–Kutta method (RK4), and a variable stepsize
Runge–Kutta method (RK45) implemented by the Matlab ode45 function with the
tolerance of $10^{-8}$. More precisely, the evolution of SPLT over step size
$h$ is written as $\phi_{h/2}\circ\psi_{h}\circ\phi_{h/2}$, where $\phi_{t}$
is the exact flow map of (36) with fixed $\Omega$, and $\psi_{t}$ is the exact
$t$-time flow map of (35) with fixed $R$ and $J=I_{3\times 3}$.

The goal of this comparison is not to claim that a certain method is superior
to the other methods. Rather, it is to identify the numerical properties of
LGVI compared with others. Having stated that, LGVI is implicit, and (28) is
solved by a general purpose nonlinear solver, instead of a numerical solver
tailored for (28). As a consequence, LGVI is substantially slower than the
three explicit methods, to the extent that the comparison is not meaningful.

Instead, for a more interesting comparison, we exploit the property of LGVI
providing consistent results for a wide range of step sizes, and we only
utilize (39) and (40) with a fixed prescribed step size. The resulting scheme,
denoted by ELGVI, is explicit as shown in (41). Overall ELGVI is quite
comparable with SPLT, but it benefits from a bit faster initial convergence,
especially when $p$ is larger and $h$ is smaller. One particular case for
$p=6$ and $h=0.001$ is illustrated in Figure 3.(a). With regard to RK4 and
RK45, their convergence is almost identical to ELGVI, but as presented in
Figure 3.(b), those methods do not preserve the orthogonality of the rotation
matrix, which is problematic. Whereas, both of LGVI and SPLT conserve the
structure of rotation matrices. Next, the computation time with Intel Core i7
3.2GHz, averaged for 10 executions, are 0.0727, 0.0258, 0.3847, and 1.1476
seconds for ELGVI, SPLT, RK4, and RK45, respectively. It is expected that RK4
requires more computation time as the gradient should be evaluated four times
per a step, and it seems that the time-adaptive RK45 algorithm requires more
frequent evaluations of the gradient.

####  V-A4 Comparison with Other Optimization Schemes on Lie Groups

Figure 4: Comparison with other accelerated optimization schemes on Lie groups

Finally, we compare ELGVI with other optimization schemes on Lie groups. In
particular, we consider variationally accelerated Lie-group methods based on
the NAG variational principle and operating splitting [7], referred to as Lie-
NAG-SC and Lie-NAG-C, which are conformally symplectic and group-structure
preserving. Note that Lie-NAG-C corresponds to SPLT with $p=2$.

Four cases are considered as marked in Figure 4 for varying $p$ and $h$.
Compared with Lie-NAG-C, ELGVI exhibits faster convergence at a higher order.
This does not contradict Nesterov’s oracle lower bound: the continuous Bregman
dynamics with $p>2$ should be discretized by smaller steps as $t$ increases,
and therefore, the asymptotic order of convergence is still
$\mathcal{O}(1/k^{2})$ as illustrated above. However, since ELGVI uses a fixed
stepsize, the initial error can decay faster than inverse quadratic, and
depending on the level of accuracy required, we can take the advantage of it
by employing early stopping. On the other hand, Lie-NAG-SC demonstrates
exponential convergence asymptotically when applied to strongly convex
functions. Overall, if moderate stopping criteria are employed, ELGVI may be
preferred, as they exhibit the fastest initial decay of the cost function.

###  V-B Optimization on $\mathsf{SO(3)}\times\mathbb{R}^{3}$

Next, we present an optimization problem on
$\mathsf{SO(3)}\times\mathbb{R}^{3}$ to estimate the position and the attitude
of a camera using the KITTI vision benchmark dataset [24]. This is to verify
the performance of ELGVI for a non-convex function in a higher-dimensional Lie
group, with more relevance to engineering practice. More specifically, we
consider $N=516$ distinct features on a single image frame, where their 2D
pixel coordinates in the image plane, and the actual 3D location in the world
coordinates are given by $p^{i}\in\mathbb{R}^{3}$ and
$P^{i}\in\mathbb{R}^{4}$, respectively as homogeneous coordinates. Assuming
that the camera calibration matrix $K\in\mathbb{R}^{3\times 3}$ is also known,
we wish to estimate the pose $(R,x)\in\mathsf{SO(3)}\times\mathbb{R}^{3}$ of
the camera.

This is formulated as an optimization problem to minimize the reprojection
error, which is the discrepancy between the actual pixel location of the
features and the features projected to the image plane by the current estimate
of $(R,x)$ [25]. For example, let $\tilde{p}^{i}\in\mathbb{R}^{3}$ be the
homogeneous coordinates for the feature corresponding to $P^{i}$ projected to
the image plane by $(R,x)$. From the perspective camera model,

$\displaystyle\lambda\tilde{p}^{i}=K[R,x]P^{i},$

for $\lambda>0$. The corresponding reprojected pixel is determined by the
dehomogenization of $\tilde{p}^{i}$, namely
$H^{-1}(\tilde{p}^{i})\in\mathbb{R}^{2}$ corresponding to the first two
elements of $\tilde{p}^{i}$ divided by the last element. The objective
function is the sum of the reprojection error given by

$\displaystyle\mathsf{f}(R,x)=\sum_{i=1}^{N}\|H^{-1}(p^{i})-H^{-1}(\tilde{p}^{i})\|^{2}.$
(43)

Figure 5 presents the optimization results by ELGVI, which are comparable to
the benchmark examples presented for $\mathsf{SO(3)}$. However, the terminal
phase is relatively noisy, partially because the gradients of (43) are
evaluated numerically with a finite-difference rule. Figure 6 illustrates the
reprojected features before and after the optimization.

Figure 5: Optimization on $\mathsf{SO(3)}\times\mathbb{R}^{3}$: convergence
with respect to $k$

(a) Initial guess $(R_{0},x_{0})$

(b) Optimized $(R^{*},x^{*})$

Figure 6: Reprojection error: the red $+$ markers denote the key points
detected, and the yellow $+$ markers represent the key points projected by the
estimated pose. The paired features are connected by solid lines.

##  VI Conclusions

In this paper, we proposed a Lie group variational integrator for the Bregman
Lagrangian dynamics on Lie groups, to construct an accelerated optimization
scheme. The variable stepsize prescribed by the extended variational principle
exhibits an interesting convergence property, and the variational
discretization is robust to the initial stepsize. It would be interesting to
explore the role of variable time-stepping in geometric discretizations of the
Bregman dynamics especially compared with Hamiltonian variational integrators.

## References

  * [1] Y. Nesterov, “A method of solving a convex programming problem with convergence rate $\mathcal{O}(1/k^{2})$,” _Soviet Mathematics Doklady_ , vol. 27, no. 2, pp. 372–376, 1983.
  * [2] ——, _Introductory lectures on convex optimization: A basic course_ , 2004\.
  * [3] W. Su, S. Boyd, and E. J. Candes, “A differential equation for modeling nesterov’s accelerated gradient method: Theory and insights,” _The Journal of Machine Learning Research_ , vol. 17, no. 1, pp. 5312–5354, 2016.
  * [4] A. Wibisono, A. C. Wilson, and M. I. Jordan, “A variational perspective on accelerated methods in optimization,” _proceedings of the National Academy of Sciences_ , vol. 113, no. 47, pp. E7351–E7358, 2016.
  * [5] M. Betancourt, M. I. Jordan, and A. C. Wilson, “On symplectic optimization,” _arXiv preprint arXiv:1802.03653_ , 2018.
  * [6] E. Hairer, C. Lubich, and G. Wanner, _Geometric Numerical Integration: Structure-preserving algorithms for ordinary differential equations_ , 2nd ed. Berlin: Springer-Verlag, 2006.
  * [7] M. Tao and T. Ohsawa, “Variational optimization on Lie groups, with examples of leading (generalized) eigenvalue problems,” in _International Conference on Artificial Intelligence and Statistics_ , 2020, pp. 4269–4280.
  * [8] M. Tao, H. Owhadi, and J. E. Marsden, “Nonintrusive and structure preserving multiscale integration of stiff odes, sdes, and hamiltonian systems with hidden slow dynamics via flow averaging,” _Multiscale Modeling & Simulation_, vol. 8, no. 4, pp. 1269–1324, 2010.
  * [9] J. Marsden and M. West, “Discrete mechanics and variational integrators,” in _Acta Numerica_. Cambridge University Press, 2001, vol. 10, pp. 317–514.
  * [10] M. Leok and J. Zhang, “Discrete Hamiltonian variational integrators,” _IMA J. Numer. Anal._ , vol. 31, no. 4, pp. 1497–1532, 2011.
  * [11] T. Lee, M. Leok, and N. McClamroch, “Lie group variational integrators for the full body problem,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 196, pp. 2907–2924, May 2007.
  * [12] V. Duruisseaux, J. Schmitt, and M. Leok, “Adaptive Hamiltonian variational integrators and symplectic accelerated optimization,” _arXiv preprint arXiv:1709.01975_ , 2021.
  * [13] J. Hu, X. Liu, Z.-W. Wen, and Y.-X. Yuan, “A brief introduction to manifold optimization,” _Journal of the Operations Research Society of China_ , vol. 8, no. 2, pp. 199–248, 2020.
  * [14] P.-A. Absil, R. Mahony, and R. Sepulchre, _Optimization Algorithms on Matrix Manifolds_. Princeton University Press, 2008.
  * [15] J. Marsden and T. Ratiu, _Introduction to Mechanics and Symmetry_ , 2nd ed., ser. Texts in Applied Mathematics. Springer-Verlag, 1999, vol. 17.
  * [16] T. Lee, M. Leok, and N. McClamroch, _Global Formulation of Lagrangian and Hamiltonian Dynamics on Manifolds_. Springer, 2018.
  * [17] C. Kane, J. Marsden, and M. Ortiz, “Symplectic-energy-momentum preserving variational integrators,” _Journal of Mathematical Physics_ , vol. 40, no. 7, pp. 3353–3371, 1999.
  * [18] Y. Nesterov, “Smooth minimization of non-smooth functions,” _Mathematical programming_ , vol. 103, no. 1, pp. 127–152, 2005.
  * [19] ——, “Accelerating the cubic regularization of Newton’s method on convex problems,” _Mathematical Programming_ , vol. 112, no. 1, pp. 159–181, 2008\.
  * [20] V. Duruisseaux and M. Leok, “A variational formulation of accelerated optimization on Riemannian manifolds,” _arXiv preprint arXiv:2101.06552_ , 2021.
  * [21] F. Alimisis, A. Orvieto, G. Bécigneul, and A. Lucchi, “A continuous-time perspective for modeling acceleration in Riemannian optimization,” in _International Conference on Artificial Intelligence and Statistics_. PMLR, 2020, pp. 1297–1307.
  * [22] T. Lee, “Computational geometric mechanics and control of rigid bodies,” Ph.D. dissertation, University of Michigan, 2008.
  * [23] G. Wahba, “A least squares estimate of satellite attitude, Problem 65-1,” _SIAM Review_ , vol. 7, no. 5, p. 409, 1965.
  * [24] A. Geiger, P. Lenz, C. Stiller, and R. Urtasun, “Vision meets robotics: The KITTI dataset,” _International Journal of Robotics Research (IJRR)_ , 2013\.
  * [25] Y. Ma, S. Soatto, J. Kosecka, and S. S. Sastry, _An invitation to 3-D vision: from images to geometric models_. Springer Science & Business Media, 2012, vol. 26.