[width=0.82height=4.13cm]Figures/Final_flow_diagram.png

Figure : Overview of attack process.

##  1 Design, Implementation and Evaluation of the Mercenary Component

###  Attack design

The attack is centered on wrist motions that occur during the writing of
methodically selected English characters. Each character maps to an answer. A
sequence of characters are carefully interleaved with pauses during the
attack. For illustration, we assume the case of an exam for which the answer
options to each question are “A”, “B”, “C” and “D”. Take the case of a
question on this exam whose answer is “A”. To implement our attack, the
mercenary first reads the question and determines that the answer to this
question is “A”. The mercenary then has the option to either immediately
trigger the attack and send this answer to the beneficiary, or to continue on
to find the answer to the next question so as to accumulate a sequence of
answers. Assuming the mercenary decides to send the answer to the first
question, he/she executes the attack using the following three steps:

  1. [label=.]

  2. 1.

initiate the attack by holding the hand still on top of the table for
$(t_{1}\pm\varepsilon)$ seconds (we also call this an opening pause).

  3. 2.

write the letter “A” on the paper (or mimic the writing of letter “A”)111The
mercenary is free to write on the exam sheet or on scratch paper, or to simply
pretend to write by tracing out a character., and,

  4. 3.

hold the hand still on the table to execute another pause of
$(t_{2}\pm\varepsilon)$ seconds (we also call this a closing pause).

For this research we have set $t_{1}=12$ seconds, $t_{2}=5$ seconds and
$\varepsilon=2$ seconds following trial experiments in which users tried out
different options for these parameters. After executing the closing pause, the
mercenary is no longer required to stick to the protocol, and may move the
hands freely, flip the pages of the exam, undertake calculations, etc. The
system ignores all these actions until the mercenary initiates another opening
pause that signals the start of the second answer. On completing this pause,
the mercenary writes the character representing the second answer, and then
executes the closing pause similarly to what was done with letter “A”. This
process continues until the end of the exam, with the questions answered in
sequential order. Because data could be delayed or reordered during transit
across the network, answers sent out to the beneficiary ought to have meta
information (e.g., time-stamps or question numbers) to enable realigning the
order of answers if needed.

###  Data collection experiments for evaluation of mercenary component

Data collection for the mercenary module was done using a smart watch app on
the Android platform. The app recorded gyroscope and
accelerometer222Experiments showed accelerometer data to not perform well for
our attack, so we do not reference it any further and instead focus on
gyroscope data going forward. data that was then sent via bluetooth to a smart
phone paired with the watch. All watches used for our experiments were either
LG Urbane (W200) or Fossil Gen 5 (FTW 4025). The watches were set to sample
sensor data at a rate of 60 Hz. To execute the attacks, we recruited 15
mercenaries each of whom wore a smart watch paired with a phone. Each
mercenary was given a detailed description of the attack protocol and allowed
to practice before beginning attack execution. Practice times were on average
about 20-30 minutes before mercenaries were ready to participate in the
experiment. The data collection experiment varied depending on whether
training or testing data was being collected. To collect training data,
mercenaries simply carefully wrote each of the characters supported by our
attack. For example, to collect training data for a character such as “A”, the
mercenary wrote the letter “A” immediately after the press of a start button
on the interface. On completion of writing the letter, a stop button was
immediately pressed so gyroscope data that precisely delimits this letter was
saved. This process was repeated 29 more times for a total of thirty samples
of letter “A”. Other characters were treated similarly (details of the full
range of characters that we evaluated will be discussed in Section 1).
Collecting training data in this way enabled us to cleanly capture the data
for each character without the pause portions of the protocol complicating
character segmentation and potentially adding noise that could adulterate the
trained model. To collect testing data, each mercenary followed the full
attack process (earlier described in Section 1), writing answers, and pausing
per specifications. Mercenaries wore the watch on the hand being used to
write, meaning that this watch both recorded the wrist gestures and provided
time measurements needed for the protocol. In total, each of the 15
mercenaries wrote 50 answers during the attack (an equivalent of solving 50
questions on an exam).

###  Classification framework

The framework to classify characters (or exam answers) begins with the
detection of opening and closing pauses, followed by the classification of
sensor data located between these two pauses into one of the answer classes.
Pause detection uses a threshold-based approach. In particular an
instantaneous gyroscope measurement $g_{x}$ is flagged as potentially being
part of pause if $-T_{h}\leq g_{x}\leq T_{h}$, where $T_{h}$ is a mercenary-
specific threshold. $T_{h}$ obtained from data collected from a training
dataset of pauses provided by the mercenary. When a series of consecutive
points meet this condition, they are compared to the criteria for the closing
and opening pauses (recall Section 1) and marked as such if they meet this
criteria. Otherwise they are ignored.

Table : List of temporal and spectral features used to train our classifiers

!  Type Feature Names Temporal features Mean, Standard Deviation,
Interquartile Range, Absolute Energy, Mean Absolute Deviation, Standard Error
of the Mean, Mean change, Autocovariance, Longest Strike Above Mean, Variance,
Absolute Sum of Changes, Kurtosis, Sample Entropy, Autocorrelation, Mean
absolute change, Sum, Skewness, Quantile, Median, Longest Strike Below Mean
and Complexity invariant distance Spectral features Spectral Centroid,
Spectral Flatness, Spectral Kurtosis, Spectral Skewness, Spectral
Decrease,Spectral Spread, Spectral Rolloff, Spectral Slope

Following pause detection, gyroscope data located between an opening and
closing pause is input to the machine learning module which infers which
character the data maps to. For each of the $X$, $Y$ and $Z$ components of
gyroscope data, the module extracts 32 features. These features are a mix of
frequency domain and time domain features (see Table 1 for the list of
features). Thus a total of 96 features are extracted from each character
(i.e., from the full window of data between the opening and closing pauses).
We classified this data with two classifiers: a random forest (with 60 trees)
and logistic regression. For classification, we used the default parameters in
python’s Scikit learn library for all classifiers as we found them to perform
well.

###  Performance of classification framework

The pause detection module detected and classified all opening and closing
pauses correctly for all mercenaries (i.e., 100% accuracy). In the bulk of
this section, we describe the character classification performance which had a
much more complex dynamics.

####  Exploratory evaluation of attack symbols

To prepare the ground for presentation of the classification performance of
the character classification portion of the framework, we first make the
distinction between what we refer to as answer options and attack symbols that
represent the answer options. Typically the four answer options in a multiple
choice exam are $A$, $B$, $C$ and $D$ for the first, second, third and forth
answer choices respectively. The attack symbols, which might or might not be
the same as the answer options, are the characters which the mercenary writes
(or pretends to write with simulated pen strokes) to represent the answer
options. The data consumed by the machine learning algorithms is in essence
the sensor data generated by the mercenary’s wrist motions when writing these
attack symbols.

Table : Results from clustering of characters using the K-means algorithm. The
results show that B, C, E, and A are in 4 different clusters most of the time,
supporting their usage for our attack.

!  User Data Clusters Cluster 1 Cluster 2 Cluster 3 Cluster 4 1
B,H,J,O,q,S,y,Z,8 C,K,m E,I,W A,X 2 B,I,J,m,O,q,S,8 C,X A,D,E,H,W,y,Z K 3
B,K,X,Z C,J,m,O,q,S,W,y,8 E,H,I A 4 B,I,J,K,W,X C,q E,O,S,Z,8 A,H,m,y

In selecting which attack symbols to be used by our attack, we first conducted
an exploratory investigation of the classification performance of the answer
options $A$, $B$, $C$ and $D$ as attack symbols given their usability
benefits. For this investigation, we collected a separate exploratory dataset
early on in the study from a subset of mercenaries. We found $B$ and $D$ to be
very frequently confused with each other (results not shown due to space
limitations), likely due to their similar shape. We decided to drop “D” as an
attack symbol and keep the others.

####  Definitive set of attack symbols

To find a replacement for “D” which wont get confused with any of the 3
confirmed attack symbols, we went back to the drawing board and analyzed the
alphabets A-Z and the numbers 0 to 9. This analysis involved two steps: a
visual selection in which we discarded a character if it had another character
having a visually similar shape to it, and, a clustering process in which we
split the remaining characters into 4 clusters using the K-means algorithm. In
the visual step, we dropped a wide range of characters and kept only 19. To
give some insight into our selection criteria, the following are examples of
characters which we dropped in favor of others at this visual stage. We chose
“E” and dropped “F”, chose “Z” and dropped “2”, chose “S” and dropped “5”,
chose “I” and dropped “1”, etc. Having dropped a subset of the characters, we
had each of our mercenaries write the remaining symbols while we captured the
associated gyroscope data. We then extracted the earlier described 96 features
and clustered the data using K-means with K=4. Table 1 shows the results from
this clustering stage for 4 representative mercenaries in our exploratory
dataset. Observe that “B” and “C” are always in clusters #1 and #2
respectively while “A” is in cluster #4 75% of the time. We observed similar
patterns for other sets of mercenaries (results not shown here). The fact that
these 3 characters are always in separate clusters, supports their choice as
attack symbols and validates the separability. For the fourth symbol, we
naturally had to focus on cluster #3 since it does not contain any of our 3
chosen characters. Observe that “E” is consistently located in this cluster,
making it our final choice of attack symbol. The final set of attack symbols
used in our fully fledged data collection experiment (described in Section 1)
were A, B, C and E to respectively represent the answer options A, B, C and D.
This is a highly usable set that should only provide a minor headache to the
mercenary (only letter E would be unusual).

[width=0.45]Figures/15_user_results_new.png

Figure : Classification performance of all 15 mercenaries.

####  Attack results for definitive set of attack symbols

Figures 1 and 1 summarize the performance of the attack when we used this
final character-set. Figure 1 shows that the Logistic Regression classifier
always performed best, achieving 80% or higher accuracy for all 15
mercenaries. The figure reveals that the Random Forest performed worse, but
still had an accuracy of 70% or more for 11 of the 15 mercenaries. l0.32

[width=0.32]Figures/Final_confusion_matrix2.png

Confusion matrix for classification performance of 15 mercenaries. To gain
insights into the distribution of mis-classifications we show the confusion
matrix obtained with the Logistic regression classifier in Figure 1. We only
show the confusion matrix for the best performing classifier due to space
limitations. The figure is an average of the confusion matrices obtained with
each of the 15 mercenaries. The figure confirms that the character “E”
provided a boost, as no single character is causing a disproportionate amount
of mis-classifications. The average accuracy rose to 86.3% with this new
character-set. A similar trend was seen with the Random Forest classifier
(results not shown here). Overall, an average accuracy of close to 90% for the
best classifier suggests a highly effective attack. That said, this answer-
level (as opposed to exam-level) behavior of the attack only paints a limited
picture of its behavior. In Section LABEL:model-1 (Page LABEL:model-1), we
more rigorously characterize the attack as we provide more context as to how a
classification performance in the above cited range relates to the
beneficiary’s odds of passing the exam.

##  2 Design, Implementation and Evaluation of the Beneficiary Component

We designed and studied two approaches for answer retrieval by the
beneficiary: a haptic approach and a visual approach. Below, we describe how
we designed these approaches as well as how they performed in our experiments.

###  Haptic answer retrieval mechanism

l0.32 [width=0.31]Figures/Vibration_flow_arrow_new.png Vibration-based system
for answer feedback to the beneficiary. Each arrow represents a single
vibration. The example represented on the figure is for the answer sequence C,
A, B.

####  Design and implementation

The haptic approach encodes the answers as vibrations. The answer options A,
B, C and D, are respectively represented by one, two, three and four
vibrations within a cluster. Figure 2 illustrates this approach for an answer
sequence of 3 answers, C, A and B. The first three arrows represent three
vibrations and hence the answer C. The definitions of the other two (A and B)
follow the same idea. Across an exam, the cluster numbers represent the
question being answered (i.e., the $i^{th}$ cluster represents the $i^{th}$
question on the exam). The approach relies on 4 parameters: $T_{1}$ (duration
of an individual vibration), $T_{2}$ (time interval between vibrations in a
cluster), $T_{3}$ (time interval between clusters) and $h$ (amplitude/strength
of the vibrations). We discuss our settings for these parameters in the
following subsection.

####  Performance evaluation

A key aspect of our performance evaluation was how to configure the 4
parameters without compromising the attack. The following considerations
guided our decisions: (1) $T_{3}$ ought to be markedly distinct from $T_{2}$
if the end-user is to discern a new cluster from a vibration within a cluster,
(2) $T_{1}$ has to be long enough for the user not to miss a vibration when it
happens, (3) $T_{2}$ should not be too short that the user fails to
distinguish between two consecutive vibrations within a cluster, (4) $h$ has
to be of such magnitude that the user is guaranteed to feel the vibrations,
but it should not be too strong that the proctor hears the audio from the
vibrations, (5) $T_{3}$ is not within the full control of the software
developer. For example, if the mercenary takes a very long time to answer a
question, $T_{3}$ could get arbitrarily long. One can only fix its minimum
value to avoid confusion with $T_{2}$. l0.30
[width=0.30height=0.24]Figures/vibration_watch_v8.png Distance at which smart
watch vibrations of various amplitudes get detected by a third party in the
room. Due to space limitations we only provide details of our sensitivity
analysis on the dynamics of $h$ since audio generated based on this parameter
has a very strong connection with students potentially being caught by the
proctor. For the other parameters, we only point out here that we ultimately
set $T_{1}=200ms$ and $T_{2}=1second$ as these values enabled all our 15
participants to reliably delineate between individual vibrations. We set the
minimum value of $T_{3}$ as 45 seconds. With these settings all our
participants successfully retrieved all the answers. Figure 2 shows results of
our sensitivity analysis on $h$. The figure reveals that the distance at which
the vibrations are heard, grows about linearly as the amplitude increases
towards 250 units. Its also clear from the plot that wearing the watch on the
wrist will always be the safer strategy. In our experiments we found an
amplitude of 70 units to work well as participants detected all vibrations and
got answers correct under this setting with the other variables set to the
earlier described values. Assuming the live attack uses this value of 70
units, the plot reveals that a beneficiary wearing the watch on the wrist
would need to sit further than 0.5 meters from the proctor’s location. If the
classroom arrangement is so packed and such an amount of room cannot be found,
the beneficiary might then have to opt for the visual feedback system
(discussed next).

###  Visual answer retrieval mechanism

####  Design and implementation

Figure 2 shows the visual interface that we designed for answer retrieval by
the beneficiary. The interface is obfuscated within an analog clock. The big
dots represent the question numbers as well as the answers to each question.
For example the red arrow points to the dot that represents Question 1.
Question numbers increase clockwise with the maximum being Question 12 for a
single page. The color of a dot represents the answer to the associated
question. The answer options “A”, “B”, “C” and “D” are respectively
represented by the colors red, green, blue and yellow, while a question that
is yet to be answered is represented as a purple dot. For example, the answer
to Question 1 is “B”, while Questions 7 to 12 are yet to be answered. l0.23
[width=0.24height=0.23]Figures/newscreen_v7.jpg Stealthy visual interface for
answer retrieval at beneficiary end of the attack. The figure shows the answer
status for a single page of 12 questions. For questions 13 and beyond, the
small dots are used as a counter of page numbers. For example, in the figure
the small dark grey dot (labelled Page 1 by the white arrow) implies that the
screen shows the first page of 12 questions. For the next page of 12
questions, the current dark grey dot reverts to the default color (light grey)
while the small dot next to it (clockwise) becomes dark grey. In general, the
screen always has one small dark grey dot which determines the page (or the
specific set of 12 questions). Since a typical exam would not have more than
just a couple of sets of 12 questions, the small dark grey dot would in
practice rarely go beyond the first 2-4 positions. Hence the beneficiary will
easily see its location given a single glance on the screen, and easily work
out the page number (and hence question numbers). To change the colors of the
dots as new answers arrive we use the Android color class to assign the
respective color constants to the various dot locations as needed.

####  Performance evaluation

We recruited a group of 15 participants to evaluate the visual feedback
system. The participants were first trained on how to retrieve the answers
given the colored dots. Each of them was then later given the task to read a
series of 15 answers. All participants found the scheme very intuitive and got
all the answers correct. Our conclusion from this was that given a minimal
amount of training, a beneficiary should in practice successfully use the
system to retrieve the correct answers.


# Symmetric Discrete Optimal Control and Deep Learning

Anthony M. Bloch
Department of Mathematics
University of Michigan
Ann Arbor MI 48109
<EMAIL_ADDRESS>Research partially supported by NSF grant
DMS-2103026, and AFOSR grants FA 9550-22-1-0215 and FA 9550-23-1-0400. Peter
E. Crouch
College of Engineering
University of Texas at Arlington
Arlington, TX
<EMAIL_ADDRESS>Tudor S. Ratiu
School of Mathematical Sciences
Ministry of Education Laboratory of Scientific Computing (MOE-LSC)
Shanghai Frontier Science Center of Modern Analysis
Shanghai Jiao Tong University
800 Dongchuan Road, 200240 China
Section de Mathématiques
Ecole Polytechnique Fédérale de Lausanne
1500 Lausanne, Switzerland
<EMAIL_ADDRESS><EMAIL_ADDRESS>Research partially supported by the
National Natural Science Foundation of China grant number 11871334 and by the
Swiss National Science Foundation grant NCCR SwissMAP.

(February 8, 2024
)

###### Abstract

We analyze discrete optimal control problems and their connection with back
propagation and deep learning. We consider in particular the symmetric
representation of the discrete rigid body equations developed via optimal
control analysis and optimal flows on adjoint orbits

In memory of Roger Brockett

Keywords: optimal control, discrete dynamics, deep learning, back propagation

##  1 Introduction

This paper presents some connections between discrete optimal control, deep
learning, and back propagation. This goes back to the work of Bryson and Ho
[1969] on discrete optimal control and, in particular, to the paper LeCun
[1988]. We show that the formalism developed in Bloch, Crouch, Marsden and
Ratiu [2002] has much in common with this latter paper. Other interesting
recent work on deep learning and optimal control includes Benning et. al.
[2019], E [2017], E [2019] and references therein.

We consider here firstly the discrete setting and discuss also the smooth
setting and double bracket equations. We consider in particular the discrete
symmetric rigid body equations developed in Bloch, Crouch, Marsden and Ratiu
[2002].

The key differences between this work and our earlier work Bloch, Crouch,
Marsden and Ratiu [2002] is, as appropriate for machine learning, a terminal
cost rather than a fixed endpoint condition and multiple sets of initial data.
The connections between deep learning and optimal control are essentially the
following: the controls are the weights, the layers are the discrete time
points, the training data or patterns are the initial data, back propagation
is solving the two point boundary problem, and the endpoint cost is the
primary cost.

This is related to our earlier work on an alternative formulation of the
$n$-dimensional rigid body equations and an associated set of discrete
equations called the symmetric representation of the discrete rigid body
equations; see Bloch, Crouch, Marsden and Ratiu [2002] and Guibout and Bloch
[2004]. Both the continuous equations and their discrete counterpart evolve on
a Cartesian product $G\times G$ of a Lie group $G$ rather than on its
cotangent bundle $T^{\ast}G$. One interesting way to derive the continuous
equations is by means of the (Pontryagin) maximum principle of optimal control
theory. Likewise, the discrete equations can be derived from discrete optimal
control theory. Extensions of the theory in the smooth setting may be found,
for example, in Bloch, Crouch, Marsden and Sanyal [2008], Gay-Balmaz, F. and
T. S. Ratiu [2011], and Bloch, Gay-Balmaz and Ratiu [2018].

It is also interesting to consider the subRiemannian setting of these problems
as we did with a view to the quantum setting in Sanyal, Mosely and Bloch
[2009].

The work described in this paper is in part expository and in part new. It
owes much to earlier work by and with Roger Brockett, see, e.g. Brockett
[1973], Brockett [1989], Brockett [1994], Bloch, Brockett and Ratiu [1990],
Bloch, Brockett and Ratiu [1992], Bloch, Brockett and Crouch [1997].

##  2 Discrete Optimal Control and Associated Deep Learning

We first consider a general class of discrete optimal control problems and we
follow with the special case of the discrete symmetric representation of the
rigid body equations

######  Definition 2.1.

Let $N$ be a positive integer and $X_{0},X_{N}\in\mathbb{R}^{n}$ be given. Let
$f:\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$ and
$g:\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{+}:=[0,\infty)$ be
smooth functions. Denote points in $\mathbb{R}^{n}$ by $x$, points in
$\mathbb{R}^{m}$ by $u$. Let $E\subset\mathbb{R}^{m}$ be a control constraint
set and assume that $E=h^{-1}(0)$, where $h:\mathbb{R}^{m}\to\mathbb{R}^{l}$
is a submersion. Let $\langle\cdot,\cdot\rangle$ denote the pairing between
vectors in $\mathbb{R}^{p}$ given by the dot product; $p$ can be $n$, $m$, or
$l$.

Define the optimal control problem:

$\min_{u_{k}\in E}\sum^{N-1}_{k=0}g(x_{k},u_{k})$ (2.1)

subject to $x_{k+1}=f(x_{k},u_{k})$, with $x_{0}=X_{0}$ and $x_{N}=X_{N}$, for
$u_{k}\in E$.

######  Proposition 2.2.

A solution to problem (2.1) satisfies the following extremal solution
equations

$p_{k}=\frac{\partial H}{\partial x_{k}}(p_{k+1},x_{k},u^{*}_{k}),\qquad
x_{k+1}=\frac{\partial H}{\partial p_{k+1}}(p_{k+1},x_{k},u^{*}_{k}),$ (2.2)

where

$\displaystyle H(p_{k+1},x_{k},u_{k}):=\langle
p_{k+1},f(x_{k},u_{k})\rangle-g(x_{k},u_{k}),\quad k=0,\ldots,N,\quad
p_{N+1}=0,$ (2.3)

and $\partial/\partial x$, $\partial/\partial p$ denote the partial $x$\- and
$p$-gradients. In these equations, $u_{k}^{\ast}$ is determined as follows.
Define

$\hat{H}(p_{k+1},x_{k},u_{k},\sigma):=H(p_{k+1},x_{k},u_{k})+\left\langle\sigma,h(u_{k})\right\rangle$
(2.4)

for $\sigma\in\mathbb{R}^{l}$, $\sigma$ a Lagrange multiplier; then
$u^{*}_{k}$ and $\sigma$ are solutions of the equations

$\frac{\partial\hat{H}}{\partial
u_{k}}(p_{k+1},x_{k},u^{*}_{k},\sigma)=0,\quad h(u^{*}_{k})=0.$ (2.5)

###### Proof.

Use the discrete maximum principle. We wish to minimize
$\sum^{N-1}_{k=0}g(x_{k},u_{k})$ subject to the given discrete dynamics and
control constraint set $E$. To implement the constraints defining $E$, we
consider maximizing the augmented functional (which is independent of $p_{0}$)

$\displaystyle V(p_{0},\ldots,p_{N},x_{0},\ldots
x_{N-1},u_{0},\ldots,u_{N-1})$
$\displaystyle\qquad:=\sum^{N-1}_{k=0}\left(\langle
p_{k+1},f(x_{k},u_{k})-x_{k+1}\rangle+\left\langle\sigma,h(u_{k})\right\rangle-g(x_{k},u_{k})\right)$
$\displaystyle\qquad\,=\sum^{N-1}_{k=0}\left(-\langle
p_{k+1},x_{k+1}\rangle+\hat{H}(p_{k+1},x_{k},u_{k},\sigma)\right)$
$\displaystyle\qquad\,=\left\langle
p_{0},x_{0}\right\rangle+\sum^{N}_{k=0}\left(-\langle
p_{k},x_{k}\rangle+\hat{H}(p_{k+1},x_{k},u_{k},\sigma)\right),$

where we set $p_{N+1}=0$ for notational convenience since there is no
$x_{N+1}$. The critical points of $V$ are hence given by

$0=\delta V=\left\langle\delta
p_{0},x_{0}\right\rangle+\sum^{N}_{k=0}\left(-\langle p_{k},\delta
x_{k}\rangle-\langle\delta p_{k},x_{k}\rangle+\frac{\partial\hat{H}}{\partial
p_{k+1}}\delta p_{k+1}+\frac{\partial\hat{H}}{\partial x_{k}}\delta
x_{k}+\frac{\partial\hat{H}}{\partial u_{k}}\delta u_{k}\right)$

since $\delta x_{0}=0$ because $x_{0}=X_{0}$ is a given constant vector. This
gives the extremal equations (2.2) and (2.5) since

$\displaystyle\frac{\partial\hat{H}}{\partial
p_{k+1}}(p_{k+1},x_{k},u_{k},\sigma)$ $\displaystyle=\frac{\partial
H}{\partial p_{k+1}}(p_{k+1},x_{k},u_{k}),$
$\displaystyle\frac{\partial\hat{H}}{\partial
x_{k}}(p_{k+1},x_{k},u_{k},\sigma)$ $\displaystyle=\frac{\partial H}{\partial
x_{k}}(p_{k+1},x_{k},u_{k}),$

and $h(u_{k})=0$, for $k=0,\ldots,N$ and $p_{N+1}=0$. ∎

Note that for the algorithm described by equation (2.2) to make sense we need
to able to compute $p_{k+1}$ from the given implicit form. This follows if the
$n\times n$ matrix with entries $\frac{\partial^{2}H}{\partial y_{i}\partial
x_{j}}(y,x,u)$ for $i,j=1,\ldots,n,$ is invertible at every
$(y,x,u)\in\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}^{m}$. We need
such a condition to be satisfied for any given algorithm.

We assume that both $u^{*}_{k}$ and $\sigma$ are determined uniquely by
equations (2.5). Also note that $u^{*}_{k}=u^{*}_{k}(p_{k+1},x_{k})$,
$\sigma=\sigma(p_{k+1},x_{k})$. Using this hypothesis, we have the following
consequence.

######  Corollary 2.3.

Assume that the extremal equations (2.2) hold. Then

$\displaystyle dp_{k}$ $\displaystyle=\frac{\partial^{2}\hat{H}}{\partial
x^{2}_{k}}(p_{k+1},x_{k},u^{*}_{k})dx_{k}+\frac{\partial^{2}\hat{H}}{\partial
p_{k+1}\partial x_{k}}(p_{k+1},x_{k},u^{*}_{k})dp_{k+1}\,,$ $\displaystyle
dx_{k+1}$ $\displaystyle=\frac{\partial^{2}\hat{H}}{\partial p_{k+1}\partial
x_{k}}(p_{k+1},x_{k},u^{*}_{k})dx_{k}+\frac{\partial^{2}\hat{H}}{\partial
p^{2}_{k+1}}(p_{k+1},x_{k},u^{*}_{k})dp_{k+1}\,.$

We remark that the implicit advance map
$\Phi:(x_{k},p_{k})\mapsto(x_{k+1},p_{k+1})$ generated by the extremal
evolution (2.2) is symplectic, i.e.,

$\Phi^{*}(dx_{k+1}\wedge dp_{k+1})=dx_{k}\wedge dp_{k}.$ (2.6)

This is easily demonstrated by using Corollary 2.3. One can also derive
symplecticity directly from Hamilton’s phase space principle; see Marsden and
West [2001].

We note that solving the above two point boundary value in practice in highly
nontrivial. Various techniques have been employed including so-called shooting
methods (see .e.g. Bryson and Ho [1969]) and this is still an ongoing topic of
research.

Now we modify this problem to include an endpoint cost and a form suitable for
machine learning.

######  Definition 2.4.

Let $N,M$ be a positive integers and $X_{0},X_{N}^{a}\in\mathbb{R}^{n}$,
$a=1,\ldots,M$ be given. Let
$f:\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$,
$g:\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{+}:=[0,\infty)$,
and $\phi:\mathbb{R}^{n}\rightarrow\mathbb{R}^{+}$ be smooth functions. Let
$E\subset\mathbb{R}^{m}$ be a control constraint set and assume that
$E=h^{-1}(0)$, where $h:\mathbb{R}^{m}\to\mathbb{R}^{l}$ is a submersion.

Define the optimal control problem:

$\min_{u_{k}\in
E}\sum_{a=1}^{M}\sum^{N-1}_{k=0}g(x^{a}_{k},u_{k})+\sum_{a=1}^{M}\phi(x^{a}_{N})$
(2.7)

subject to $x^{a}_{k+1}=f^{a}(x_{k},u_{k}):=f(x_{k},u_{k},X_{0}^{a})$ with
$x_{0}=X^{a}_{0}$ and $x^{a}_{N}=X_{N}^{a}$, for $u_{k}\in E$,
$k=0,\ldots,N-1$, and $a=1,\ldots,M$.

The proof of the next proposition is analogous to that of Proposition 2.2.

######  Proposition 2.5.

A solution to problem (2.7) satisfies the following extremal solution
equations

$p^{a}_{k}=\frac{\partial H}{\partial
x_{k}}(p^{a}_{k+1},x_{k},u^{*}_{k}),\qquad x^{a}_{k+1}=\frac{\partial
H}{\partial p_{k+1}}(p^{a}_{k+1},x_{k},u^{*}_{k}),$ (2.8)

where $k=0,\ldots,N-1$, $p^{a}_{N+1}=0$ for all $a=1,\ldots,M$, and

$\displaystyle H(p_{k+1},x_{k},u_{k})=\sum_{a=1}^{M}\left(\langle
p^{a}_{k+1},f^{a}(x_{k},u_{k})\rangle-g(x^{a}_{k},u_{k})-\phi(x^{a}_{N})\right).$
(2.9)

In these equations, $u_{k}^{\ast}$ is determined as follows. Let

$\hat{H}(p_{k+1},x_{k},u_{k},\sigma):=\sum_{a=1}^{M}\left(\langle
p^{a}_{k+1},f^{a}(x_{k},u_{k})\rangle-g(x^{a}_{k},u_{k})\right)+\left\langle\sigma,h(u_{k})\right\rangle$
(2.10)

for $\sigma\in\mathbb{R}^{l}$, $\sigma$ a Lagrange multiplier. Then
$u^{*}_{k}$ and $\sigma$ are solutions of the equations

$\frac{\partial\hat{H}}{\partial
u_{k}}(p^{a}_{k+1},x^{a}_{k},u^{*}_{k},\sigma)=0,\quad h(u^{*}_{k})=0.$ (2.11)

In addition we have the endpoint condition

$p^{a}_{N}=\frac{\partial\phi(x^{a}_{N})}{\partial x^{a}_{N}},\quad
a=1,\ldots,M.$ (2.12)

###### Proof.

Use the discrete maximum principle. We wish to minimize
$\sum_{a=1}^{M}\sum^{N-1}_{k=0}g(x^{a}_{k},u_{k})+\sum_{a=1}^{M}\phi(x^{a}_{N})$
subject to the given discrete dynamics and control constraint set $E$. To
implement the constraints defining $E$, we consider maximizing the augmented
functional (which is independent of $p_{0}$)

$\displaystyle V(p_{0},\ldots,p_{N},x_{0},\ldots
x_{N-1},u_{0},\ldots,u_{N-1})$
$\displaystyle\qquad:=\sum^{N-1}_{k=0}\sum^{M}_{a=1}\left(\langle
p^{a}_{k+1},f^{a}(x^{a}_{k},u_{k})-x^{a}_{k+1}\rangle+\left\langle\sigma,h(u_{k})\right\rangle-g(x^{a}_{k},u_{k})\right)-\sum_{a=1}^{M}(p^{a}(N)x^{a}(N)-\phi(x^{a}_{N}))$
$\displaystyle\qquad\,=\sum_{a-1}^{M}\sum^{N-1}_{k=0}\left(-\langle
p^{a}_{k+1},x^{a}_{k+1}\rangle+\hat{H}(p^{a}_{k+1},x^{a}_{k},u_{k},\sigma)\right)-\sum_{a=1}^{M}(p^{a}(N)x^{a}(N)-\phi(x^{a}_{N}))$
$\displaystyle\qquad\,=\sum_{a=1}^{M}\left(\left\langle
p^{a}_{0},x^{a}_{0}\right\rangle+\sum^{N}_{k=0}(-\langle
p^{a}_{k},x^{a}_{k}\rangle+\hat{H}(p^{a}_{k+1},x^{a}_{k},u_{k},\sigma))\right)-\sum_{a=1}^{M}(p^{a}(N)x^{a}(N)-\phi(x^{a}_{N}))$

where we set $p^{a}_{N+1}=0$ for notational convenience since there is no
$x_{N+1}$. The critical points of $V$ are hence given by

$\displaystyle 0=\delta V=$ $\displaystyle\left\langle\delta
p_{0},x_{0}\right\rangle+\sum^{N}_{k=0}\left(\sum_{a=1}^{M}(-\langle
p^{a}_{k},\delta x^{a}_{k}\rangle-\langle\delta
p^{a}_{k},x^{a}_{k}\rangle+\frac{\partial\hat{H}}{\partial p^{a}_{k+1}}\delta
p^{a}_{k+1}+\frac{\partial\hat{H}}{\partial x^{a}_{k}}\delta
x^{a}_{k})+\frac{\partial\hat{H}}{\partial u_{k}}\delta u_{k}\right)$
$\displaystyle+\sum_{a=1}^{M}\left(\delta p^{a}(N)x^{a}(N)+p^{a}(N)\delta
x^{a}(N)-\frac{\partial\phi(x^{a}(N)}{\partial x^{a}(N)}\delta
x^{a}(N)\right),$

since $\delta x_{0}=0$ because $x_{0}=X_{0}$ is a given constant vector. This
gives the extremal equations (2.8), (2.11) and (2.12) since

$\displaystyle\frac{\partial\hat{H}}{\partial
p_{k+1}}(p_{k+1},x_{k},u_{k},\sigma)$ $\displaystyle=\frac{\partial
H}{\partial p_{k+1}}(p_{k+1},x_{k},u_{k}),$
$\displaystyle\frac{\partial\hat{H}}{\partial
x_{k}}(p_{k+1},x_{k},u_{k},\sigma)$ $\displaystyle=\frac{\partial H}{\partial
x_{k}}(p_{k+1},x_{k},u_{k}),$

and $h(u_{k})=0$, for $k=0,\ldots,N$ and $p_{N+1}=0$ and $p_{N}$ is fixed ∎

#### Remark

1\. As described in Benning et. al. [2019], a common choice for $f$ is
$f(x,u)=\sigma(Kx+\beta)$, where $u=(K,\beta)$ and $sigma$ is the sigmoid
function. This is the so-called ResNet framework. We can, of course, consider
other problems of this type but here we are interested in a particular group
theoretic form.

2\. The form of the solution in Proposition 2.5 is very close to that of LeCun
[1988] and, at least on a superficial level, even more so in the rigid body
case.

###  2.1 The discrete symmetric rigid body

We now derive the discrete symmetric rigid body equations by considering
discrete optimal control on the special orthogonal group.

######  Definition 2.6.

Let $\Lambda$ be a positive definite diagonal matrix. Let
$\overline{Q}_{0},\overline{Q}_{N}\in\operatorname{SO}(n)$ be given and fixed.
Let

$\hat{V}(U_{0},\ldots,U_{N-1}):=\sum_{k=0}^{N-1}\operatorname{trace}(\Lambda
U_{k}),\quad U_{k}\in\operatorname{SO}(n).$ (2.13)

Define the optimal control problem

$\mathop{\rm
min}_{U_{k}\in\operatorname{SO}(n)}\hat{V}(U_{0},\ldots,U_{N-1})=\mathop{\rm
min}_{U_{k}\in\operatorname{SO}(n)}\sum_{k=0}^{N-1}\operatorname{trace}(\Lambda
U_{k})$ (2.14)

subject to dynamics and initial and final data

$Q_{r+1}=Q_{k}U_{r},\qquad Q_{0}=\overline{Q}_{0},\qquad
Q_{N}=\overline{Q}_{N}$ (2.15)

for $Q_{r},U_{r}\in\operatorname{SO}(n)$, $r=0,1,\ldots,N-1$.

######  Theorem 2.7.

A solution of the optimal control problem in Definition 2.6 satisfies the
optimal evolution equations

$Q_{k+1}=Q_{k}U_{k}\,,\qquad P_{k+1}=P_{k}U_{k}\,,\qquad k=0,\ldots,N-1,$
(2.16)

where $P_{k}\in\operatorname{SO}(n)$ is the discrete covector in the discrete
maximum principle and $U_{k}\in\operatorname{SO}(n)$ is defined by

$U_{k}\Lambda-\Lambda U^{T}_{k}=Q_{k}^{T}P_{k}-P_{k}^{T}Q_{k}\,.$ (2.17)

Equation (2.17) can be solved for $U_{k}$ under certain circumstances, as
discussed in Moser and Veselov [1991] and Cardoso and Leite [2001]; we discuss
this issue further below.

###### Proof.

Applying Proposition 2.2, we get

$H(P_{k+1},Q_{k},U_{k})=\operatorname{trace}(P_{k+1}^{T}Q_{k}U_{k})-\operatorname{trace}(\Lambda
U_{k})=\operatorname{trace}\left((P^{T}_{k+1}Q_{k}-\Lambda)U_{k}\right)\,.$
(2.18)

Note that

$\hat{V}(U_{0},\ldots,U_{N-1})=\sum_{k=0}^{N-1}\operatorname{trace}(\Lambda
U_{k})=\sum_{k=0}^{N-1}\operatorname{trace}(U^{T}_{k}\Lambda)=\sum_{k=0}^{N-1}\operatorname{trace}(Q_{k}\Lambda
Q_{k+1}^{T})$

is the Moser-Veselov functional Moser and Veselov [1991] and that this
functional is linear in the controls.

We need to find the critical points of $H(P_{k+1},Q_{k},U_{k})$ where
$U_{k}^{T}U_{k}=I$ since $U_{k}\in\operatorname{SO}(n)$. Thus, we need to
minimize a functional of the form $\operatorname{trace}(AU)$, $A$ fixed,
subject to $U^{T}U=I$. Set

$\tilde{V}(U):=\operatorname{trace}(AU)+\tfrac{1}{2}\operatorname{trace}\left(\Sigma(U^{T}U-I)\right)\,,$

where $U\in\operatorname{SO}(n)$ and $\Sigma=\Sigma^{T}$ is a $N\times N$
matrix of Lagrange multipliers. Then $\delta\tilde{V}(U)\cdot\delta
U=\operatorname{trace}(A\delta U+\Sigma U^{T}\delta U)=0$ implies $A+\Sigma
U^{T}=0$ where $U^{T}U=I$. Hence $\Sigma=-AU$. But since $\Sigma=\Sigma^{T}$
the extrema of our optimization problem are obtained when $AU=U^{T}A^{T}$.
Applying this observation to our case (see (2.18)), we have
$\nabla_{U_{k}}H=0$ when

$\left(P_{k+1}^{T}Q_{k}-\Lambda\right)U_{k}=U^{T}_{k}\left(Q_{k}^{T}P_{k+1}-\Lambda\right)\,,$

that is,

$U^{T}_{k}\Lambda-\Lambda
U_{k}=U^{T}_{k}Q^{T}_{k}P_{k+1}-P^{T}_{k+1}Q_{k}U_{k}$

or, equivalently,

$U_{k}\Lambda-\Lambda
U^{T}_{k}=-Q^{T}_{k}P_{k+1}U^{T}_{k}+U_{k}P^{T}_{k+1}Q_{k}\,.$ (2.19)

Also, by (2.2),

$\displaystyle P_{k}$
$\displaystyle=\nabla_{Q_{k}}H=\left(U_{k}P^{T}_{k+1}\right)^{T}=P_{k+1}U^{T}_{k}\,,\qquad
Q_{k+1}=\nabla_{P_{k+1}}H=Q_{k}U_{k}\,.$

Hence we obtain equations (2.16). Combining (2.16) with (2.19) we get

$U_{k}\Lambda-\Lambda U^{T}_{k}=P_{k}^{T}Q_{k}-Q_{k}^{T}P_{k}$ (2.20)

Now replace $P_{k}$ by $-P_{k}$ and $P_{k+1}$ by $-P_{k+1}$; thus (2.16)
remains unchanged but (2.20) is transformed to (2.17) which yields the stated
result. ∎

We now define the symmetric representation of the discrete rigid body
equations as follows:

$Q_{k+1}=Q_{k}U_{k},\qquad P_{k+1}=P_{k}U_{k},\qquad k=0,\ldots,N-1,$ (2.21)

where $U_{k}\in\operatorname{SO}(n)$ is defined by

$U_{k}\Lambda-\Lambda U^{T}_{k}=Q_{k}^{T}P_{k}-P_{k}^{T}Q_{k}.$ (2.22)

We will write this as

$J_{D}U_{k}=Q_{k}^{T}P_{k}-P_{k}^{T}Q_{k}$ (2.23)

where $J_{D}:\operatorname{SO}(n)\rightarrow\mathfrak{so}(n)$ (the discrete
version of the moment of inertia operator $J$) is defined by
$J_{D}U:=U\Lambda-\Lambda U^{T}$. Notice that the derivative of $J_{D}$ at the
identity in the direction $\Omega\in\mathfrak{so}(n)$ is
$J:\mathfrak{so}(n)\ni\Omega\mapsto\Omega\Lambda+\Lambda\Omega\in\mathfrak{so}(n)$,
the classical moment of inertia operator on $\mathfrak{so}(n)$. Since $J$ is
invertible, $J_{D}$ is a diffeomorphism from a neighborhood of the identity in
$\operatorname{SO}(n)$ to a neighborhood of $0$ in $\mathfrak{so}(n)$. Using
these equations, we have the algorithm $(Q_{k},P_{k})\mapsto(Q_{k+1},P_{k+1})$
defined by: compute $U_{k}$ from (2.22), compute $Q_{k+1}$ and $P_{k+1}$ using
(2.21). Note that the update map for $Q$ and $P$ is done in parallel.

Equation (2.22) can be solved for $U_{k}$ under certain circumstances, as
discussed above, in Moser and Veselov [1991], and in Cardoso and Leite [2001];
we come back later to this issue.

As discussed in Bloch, Crouch, Marsden and Ratiu [2002] these equations are
equivalent on a certain set to the discrete Moser-Veselov equations for the
classical rigid body if we identify $U$ with the body angular momentum. We
shall say more about Moser-Veselov as well as equivalence in the smooth
setting below.

We can now obtain the machine learning generalization of the discrete rigid
body equations.

######  Definition 2.8.

Let $\Lambda$ be a positive definite diagonal matrix. Let
$\overline{Q}_{0},\overline{Q}_{N}^{a}\in\operatorname{SO}(n)$,
$a=1,\ldots,M$, be given and fixed. Let

$\hat{V}(U_{0},\ldots,U_{N-1}):=\sum_{k=1}^{N-1}\operatorname{trace}(\Lambda
U_{k})+\sum_{a=1}^{M}\phi(\overline{Q}^{a}_{N})$ (2.24)

and let $\phi:\mathbb{R}^{n}\times\mathbb{R}^{n}\rightarrow\mathbb{R}^{+}$ be
a given smooth function. Define the optimal control problem

$\mathop{\rm
min}_{U_{k}\in\operatorname{SO}(n)}\hat{V}(U_{0},\ldots,U_{N-1})=\mathop{\rm
min}_{U_{k}\in\operatorname{SO}(n)}\sum_{k=1}^{N-1}\operatorname{trace}(\Lambda
U_{k})+\sum_{a=1}^{M}\phi(\overline{Q}^{a}_{N})$ (2.25)

subject to dynamics and initial data

$Q^{a}_{r+1}=Q^{a}_{r}U_{r},\qquad Q^{a}_{0}=\overline{Q}^{a}_{0},\qquad
r=0,\ldots,N-1,$ (2.26)

for $Q_{k},U_{k}\in\operatorname{SO}(n)$.

######  Theorem 2.9.

A solution of the optimal control problem in Defintion 2.8 satisfies the
optimal evolution equations for each $a=1,\ldots,M$,

$Q^{a}_{k+1}=Q^{a}_{k}U_{k}\qquad P^{a}_{k+1}=P^{a}_{k}U_{k}\,,\qquad
k=0,\ldots,N-1,\qquad Q^{a}_{0}=\overline{Q}^{a}_{0},$ (2.27)

where $P^{a}_{k}$ is the discrete covector in the discrete maximum principle
and $U_{k}$ is defined by

$U_{k}\Lambda-\Lambda
U^{T}_{k}=\sum_{{}^{M}a=1}((Q^{a}_{k})^{T}P^{a}_{k}-(P^{a}_{k})^{T}Q^{a}_{k})\,.$
(2.28)

with

$P^{a}_{N}=\frac{\partial\phi(Q^{a}_{N})}{\partial Q^{a}_{N}}\,.$ (2.29)

###### Proof.

We apply Proposition 2.5 with

$\displaystyle H(P_{k+1},Q_{k},U_{k})$
$\displaystyle=\sum_{a=1}^{M}\operatorname{trace}((P^{a}_{k+1})^{T}Q^{a}_{k}U_{k})+\sum_{a=1}^{M}\phi(Q^{a}_{N})-\operatorname{trace}(\Lambda
U_{k})$ (2.30)

Then the computation is as above for the optimal control setting.

∎

###  2.2 Classical Moser-Veselov equations

The dynamics above are related to the Moser-Veselov equations as discussed in
Bloch, Crouch, Marsden and Ratiu [2002].

The Moser-Veselov equations for the discrete rigid body go back to Veselov
[1988] and Moser and Veselov [1991]. Their work is closely related to the
development of variational integrators; see, e.g., Marsden, Pekarsky and
Shkoller [1999] and Kane, Marsden, Ortiz and West [2000]. Another approach to
integrating differential equations on manifolds is discussed in Crouch and
Grossman [1993]. See also Iserles, McLachlan and Zanna [1999], Budd and
Iserles [1999] and Bobenko and Suris [1999].

#### Review of the Moser-Veselov Discrete Rigid Body.

We briefly review the Moser and Veselov [1991] discrete rigid body equations.
Let $Q_{k}\in\operatorname{SO}(n)$ denote the rigid body configuration at time
$k$, let $\Omega_{k}\in\operatorname{SO}(n)$ denote the body angular velocity
at time $k$, and let $M_{k}$ denote the body angular momentum at time $k$.
These quantities are related by the Moser-Veselov equations

$\displaystyle\Omega_{k}$ $\displaystyle=Q_{k}^{T}Q_{k-1}$ (2.31)
$\displaystyle M_{k}$ $\displaystyle=\Omega_{k}^{T}\Lambda-\Lambda\Omega_{k}$
(2.32) $\displaystyle M_{k+1}$ $\displaystyle=\Omega_{k}M_{k}\Omega_{k}^{T}.$
(2.33)

These equations may be viewed as defining two different algorithms.

#### MV-Algorithm 1.

Define the step ahead map

$\left(Q_{k},Q_{k+1}\right)\mapsto\left(Q_{k+1},Q_{k+2}\right)$ (2.34)

as follows: compute $\Omega_{k+1}$ from (2.31), compute $M_{k+1}$ from (2.32),
compute $M_{k+2}$ from (2.33), compute $\Omega_{k+2}$ from (2.32) and then
compute $Q_{k+2}$ from (2.31).

#### Remark.

Given $M_{k}$, conditions under which equation (2.32) is solvable for
$\Omega_{k}$ are discussed in Moser and Veselov [1991] and Cardoso and Leite
[2001].

#### MV-Algorithm 2.

Define the map:

$\left(Q_{k},M_{k}\right)\mapsto\left(Q_{k+1},M_{k+1}\right)$ (2.35)

as follows: compute $\Omega_{k}$ from (2.32), compute $M_{k+1}$ from (2.33),
compute $\Omega_{k+1}$ from (2.32) and compute $Q_{k+1}$ from (2.31).

#### Discrete Variational Principle.

The Moser-Veselov equations (2.31)-(2.33) can be obtained by a discrete
variational principle, as was done in Moser and Veselov [1991]. This
variational principle has the general form described in discrete mechanics;
see, e.g., Marsden and Wendlandt [1997], Bobenko and Suris [1999], and Marsden
and West [2001]. Namely, stationary points of the functional

$\hat{S}=\sum_{k}\operatorname{trace}(Q_{k}\Lambda Q_{k+1}^{T})$ (2.36)

on sequences of orthogonal $n\times n$ matrices yield the Moser-Veselov
equations. This variational approach can be justified as in Marsden, Pekarsky
and Shkoller [1999].

As mentioned above we can prove that symmetric representation of the rigid
body equations and the Moser-Veselov equations are equivalent when restricted
to a suitable set. It is easy to see the following: suppose that we have a
solution $(Q_{k},P_{k})$ to the symmetric discrete rigid body equations We can
then produce a solution $(Q_{k+1},M_{k+1})$ of the Moser-Veselov equations:

$M_{k+1}=Q_{k}^{T}P_{k}-P_{k}^{T}Q_{k}$ (2.37)

will give us the required $M_{k+1}$ that does the job. We refer to Bloch,
Crouch, Marsden and Ratiu [2002] for the full proof of equivalence in the
discrete setting and we shall discuss below equivalence of the symmetric and
standard rigid body in the smooth setting.

##  3 Smooth Setting and Double Bracket Equations

These general ideas can also be recast in the smooth setting. We consider here
the rigid body analysis followed by an analysis of certain flows on adjoint
orbits.

###  3.1 Classical $n$-dimensional rigid body equations

Firstly we review the classical rigid body equations in $n$ dimensions for
completeness. Further details may be found in Bloch, Crouch, Marsden and Ratiu
[2002].

We use the following pairing (multiple of the Killing form) on
$\mathfrak{so}(n)$, the Lie algebra of $n\times n$ real skew matrices regarded
as the Lie algebra of the $n$-dimensional proper rotation group
$\operatorname{SO}(n)$:

$\left\langle\xi,\eta\right\rangle=-\frac{1}{2}\operatorname{trace}(\xi\eta).$
(3.1)

The factor of $1/2$ in (3.1) is to make this inner product agree with the
usual inner product on $\mathbb{R}^{3}$ when it is identified with
$\mathfrak{so}(3)$ in the following standard way: associate the $3\times 3$
skew matrix $\hat{u}$ to the vector $u$ by $\hat{u}\cdot v=u\times v$, where
$u\times v$ is the usual cross product in ${\mathbb{R}}^{3}$.

We use this inner product to identify the dual of the Lie algebra, namely
$\mathfrak{so}(n)^{\ast}$, with the Lie algebra $\mathfrak{so}(n)$.

We recall from Manakov [1976] and Ratiu [1980] that the left invariant
generalized rigid body equations on $\operatorname{SO}(n)$ may be written as

$\dot{Q}=Q\Omega;\qquad\dot{M}=[M,\Omega]\,,$ (3.2)

where $Q\in\operatorname{SO}(n)$ denotes the configuration space variable (the
attitude of the body), $\Omega=Q^{-1}\dot{Q}\in\mathfrak{so}(n)$ is the body
angular velocity, and

$M:=J(\Omega)=\Lambda\Omega+\Omega\Lambda\in\mathfrak{so}(n)$

is the body angular momentum. Here
$J:\mathfrak{so}(n)\rightarrow\mathfrak{so}(n)$ is the symmetric (with respect
to the inner product (3.1)), positive definite, and hence invertible, operator
defined by

$J(\Omega)=\Lambda\Omega+\Omega\Lambda,$

where $\Lambda$ is a diagonal matrix satisfying $\Lambda_{i}+\Lambda_{j}>0$
for all $i\neq j$. For $n=3$ the elements of $\Lambda_{i}$ are related to the
standard diagonal moment of inertia tensor $I$ by
$I_{1}=\Lambda_{2}+\Lambda_{3}$, $I_{2}=\Lambda_{3}+\Lambda_{1}$,
$I_{3}=\Lambda_{1}+\Lambda_{2}$.

The equations $\dot{M}=[M,\Omega]$ are readily checked to be the Euler-
Poincaré equations on $\mathfrak{so}(n)$ for the Lagrangian
$l(\Omega)=\frac{1}{2}\left\langle\Omega,J(\Omega)\right\rangle.$ This
corresponds to the Lagrangian on $T\operatorname{SO}(n)$ given by

$L(g,\dot{g})=\frac{1}{2}\left\langle
g^{-1}\dot{g},J(g^{-1}\dot{g})\right\rangle\,.$ (3.3)

We note that the dynamic rigid body equations on $\operatorname{SO}(n)$ and
indeed on any semisimple Lie group are integrable (Mishchenko and Fomenko
[1976]). A key observation in this regard, due to Manakov, was that one could
write the generalized rigid body equations as Lax equations with parameter:

$\frac{d}{dt}(M+\lambda\Lambda^{2})=[M+\lambda\Lambda^{2},\Omega+\lambda\Lambda],$
(3.4)

where $M=J(\Omega)=\Lambda\Omega+\Omega\Lambda$, as in §2. The nontrivial
coefficients of $\lambda$ in the traces of the powers of
$M+\lambda\Lambda^{2}$ then yield the right number of independent integrals in
involution to prove integrability of the flow on a generic adjoint orbit of
$\operatorname{SO}(n)$ (identified with the corresponding coadjoint orbit).
Useful references are Bogayavlenski [1994] and Fedorov and Kozlov [1995].)
Moser and Veselov [1991] show that there is a corresponding formulation of the
discrete rigid body equations with parameter.

###  3.2 Smooth optimal control and the symmetric rigid body equations

Now we briefly review, see Bloch and Crouch [1996] and Bloch, Brockett and
Crouch [1997], two results which link the theory of optimal control with the
rigid body equations.

######  Definition 3.1.

Let $T>0$, $Q_{0},Q_{T}\in\operatorname{SO}(n)$ be given and fixed. Let the
rigid body optimal control problem be given by

$\mathop{\rm min}_{U\in\mathfrak{so}(n)}\frac{1}{4}\int_{0}^{T}\langle
U,J(U)\rangle dt$ (3.5)

subject to the constraint on $U$ that there be a curve
$Q(t)\in\operatorname{SO}(n)$ such that

$\dot{Q}=QU\qquad Q(0)=Q_{0},\qquad Q(T)=Q_{T}.$ (3.6)

######  Proposition 3.2.

The rigid body optimal control problem has optimal evolution equations

$\dot{Q}=QU\qquad\dot{P}=PU$ (3.7)

where $P$ is the costate vector given by the maximum principle.

The optimal controls in this case are given by

$U=J^{-1}(Q^{T}P-P^{T}Q).$ (3.8)

#### Remark.

The proof (see Bloch and Crouch [1996]) simply involves, as in the discrete
analysis above, writing the Hamiltonian of the maximum principle as

$H=\left\langle P,QU\right\rangle+\frac{1}{4}\left\langle
U,J(U)\right\rangle,$ (3.9)

where the costate vector $P$ is a multiplier enforcing the dynamics, and then
maximizing with respect to $U$ in the standard fashion (see, for example,
Brockett [1973]).

We refer to the equations (3.7) as the symmetric representation of the rigid
body equations. We can now recover the classical rigid body equations:

######  Proposition 3.3.

If $(Q,P)$ is a solution of (3.7), then $(Q,M)$ where $M=J(\Omega)$,
$\Omega=U$, and $\Omega=Q^{-1}\dot{Q}$ satisfies the rigid body equations
(3.2).

###### Proof.

Differentiating $M=Q^{T}P-P^{T}Q$ and using the equations (3.7) gives the
second of the equations (3.2). ∎

While in general there are no constraints on the costate vector
$P\in\mathfrak{gl}(n)$ one can consider the restriction of the extremal flows
to invariant submanifolds. This limits possible extremal trajectories that can
be recovered. For example this system restricts to a system on
$\operatorname{SO}(n)\times\operatorname{SO}(n)$. One can make other
assumptions on the costate vector. For example, suppose we assume a costate
vector $B$ such that $Q^{T}B$ is skew. Then it is easy to check that that the
extremal evolution equations become

$\dot{Q}=QJ^{-1}(Q^{T}B);\qquad\dot{B}=BJ^{-1}(Q^{T}B)\,,$ (3.10)

and that these equations restrict to an invariant submanifold defined by the
condition that $Q^{T}B$ is skew symmetric. These are the McLachlan-Scovel
equations (McLachlan and Scovel [1995]).

We can now generalize to the machine learning setting:

######  Definition 3.4.

Let $T>0$, $Q_{0},Q_{T}\in\operatorname{SO}(n)$ be given and fixed. Let the
rigid body optimal control problem be given by

$\mathop{\rm min}_{U\in\mathfrak{so}(n)}\frac{1}{4}\int_{0}^{T}\langle
U,J(U)\rangle dt+\sum_{a=1}^{M}\phi(Q^{a}_{T})$ (3.11)

subject to the constraint on $U$ that there be a curve
$Q(t)\in\operatorname{SO}(n)$ such that

$\dot{Q}^{a}=Q^{a}U\qquad Q^{a}(0)=Q^{a}_{0},\,a=1\dots M.$ (3.12)

######  Proposition 3.5.

The smooth machine learning symmetric rigid body flow has optimal evolution
equations

$\dot{Q}^{a}=Q^{a}U,\,\dot{P}^{a}=P^{a}U$ (3.13)

where $P$ is the costate vector given by the maximum principle

The optimal controls in this case are given by

$U=\sum_{a}J^{-1}((Q^{a})^{T}P^{a}-(P^{a})^{T}Q^{a}).$ (3.14)

and we have the endpoint conditions

$P^{a}_{T}=\frac{\partial\phi(Q^{a}_{T})}{\partial Q^{a}_{T}}$ (3.15)

###  3.3 Local equivalence of classical rigid body and symmetric rigid body
equations.

Above we saw that solutions of the symmetric rigid body equations can be
mapped to solutions of the rigid body system. As in Bloch, Crouch, Marsden and
Ratiu [2002] we can consider the converse question. Thus, suppose we have a
solution $(Q,M)$ of the standard left invariant rigid body equations. We seek
to solve for $P$ in the expression

$M=Q^{T}P-P^{T}Q.$ (3.16)

For the following discussion, it will be convenient to make use of the
operator norm on matrices. Recall that this norm is given by $\|A\|_{\rm
op}=\sup\left\\{\|Ax\|\mid\|x\|=1\right\\},$ where the norms on the right hand
side are the usual Euclidean space norms.

Since elements of $\operatorname{SO}(n)$ have operator norms bounded by $1$
and since the operator norm satisfies $\|AB\|_{\rm op}\leq\|A\|_{\rm
op}\|B\|_{\rm op}$, we see that if $M$ satisfies $M=Q^{T}P-P^{T}Q$, then
$\|M\|_{\rm op}\leq 2$. Therefore, $\|M\|_{\rm op}\leq 2$ is a necessary
condition for solvability of (3.16) for $P$.

######  Definition 3.6.

Let $C$ denote the set of $(Q,P)$ that map to $M$’s with operator norm equal
to 2 and let $S$ denote the set of $(Q,P)$ that map to $M$’s with operator
norm strictly less than 2. Also denote by $S_{M}$ the set of points $(Q,M)\in
T^{*}\operatorname{SO}(n)$ with $\|M\|_{\rm op}<2$. For the left invariant
system we trivialize
$T^{*}\operatorname{SO}(n)\cong\operatorname{SO}(n)\times\mathfrak{so}(n)^{\ast}$
by means of left translation to the identity and we identify
$\mathfrak{so}(n)^{\ast}$ with $\mathfrak{so}(n)$ using the Killing metric
(3.1).

Note that $C$ contains pairs $(Q,P)$ with the property that $Q^{T}P$ is both
skew and orthogonal.

Recall that $\sinh:\mathfrak{so}(n)\rightarrow\mathfrak{so}(n)$ is defined by
$\sinh\xi=\left(e^{\xi}-e^{-\xi}\right)/2$. One sees that indeed $\sinh$ takes
values in $\mathfrak{so}(n)$ by using, for example, its series expansion:

$\sinh\xi=\xi+\frac{1}{3!}\xi^{3}+\frac{1}{5!}\xi^{5}+\ldots.$

Recall from calculus that the inverse function $\sinh^{-1}(u)$ has a
convergent power series expansion for $|u|<1$ that is given by integrating the
power series expansion of the function $1/\sqrt{1+u^{2}}$ term by term. This
power series expansion shows that the map
$\sinh:\mathfrak{so}(n)\rightarrow\mathfrak{so}(n)$ has an inverse on the set
$U=\left\\{u\in\mathfrak{so}(n)\mid\|u\|_{\rm op}<1\right\\}$. We shall denote
this inverse by $\sinh^{-1}$, so $\sinh^{-1}:U\rightarrow\mathfrak{so}(n).$

######  Proposition 3.7.

For $\|M\|_{\rm op}<2$, the equation(3.16) has the solution

$P=Q\left(e^{\sinh^{-1}M/2}\right)$ (3.17)

###### Proof.

Notice that $M=e^{\sinh^{-1}M/2}-e^{-\sinh^{-1}M/2}\,.$ ∎

######  Proposition 3.8.

The sets $C$ and $S$ are invariant under the double rigid body equations.

###### Proof.

Notice that the operator norm is invariant under conjugation; that is, for
$Q\in\operatorname{SO}(n)$ and $M\in\mathfrak{so}(n)$, we have
$\|QMQ^{-1}\|_{\rm op}=\|M\|_{\rm op}.$ This is readily checked from the
definition of the operator norm. Recall that under the identification of the
dual $\mathfrak{so}(n)^{\ast}$ with the space $\mathfrak{so}(n)$, the
coadjoint action agrees with conjugation. Thus, the map
$f:\mathfrak{so}(n)\rightarrow\mathbb{R}$; $M\mapsto\|M\|_{\rm op}$ is a
Casimir function and so is invariant under the dynamics. In particular, its
level sets are invariant and so the sets $S$ and $C$ are invariant. ∎

#### The Hamiltonian form of the symmetric rigid body equations.

Recall that the classical rigid body equations are Hamiltonian on
$T^{*}\operatorname{SO}(n)$ with respect to the canonical symplectic structure
on the cotangent bundle of $\operatorname{SO}(n)$. The following result gives
the corresponding theorem for the symmetric case. The proof is given in Bloch,
Crouch, Marsden and Ratiu [2002]

######  Proposition 3.9.

Consider the Hamiltonian system on the symplectic vector space
$\mathfrak{gl}(n)\times\mathfrak{gl}(n)$ with the symplectic structure

$\Omega_{\mathfrak{gl}(n)}(\xi_{1},\eta_{1},\xi_{2},\eta_{2})=\frac{1}{2}\operatorname{trace}(\eta_{2}^{T}\xi_{1}-\eta_{1}^{T}\xi_{2})$
(3.18)

where $(\xi_{i},\eta_{i})\,,i=1,2$ are elements of
$\mathfrak{gl}(n)\times\mathfrak{gl}(n)$ and Hamiltonian

$H(\xi,\eta)=-\frac{1}{8}\operatorname{trace}\left[\left(J^{-1}(\xi^{T}\eta-\eta^{T}\xi)\right)\left(\xi^{T}\eta-\eta^{T}\xi\right)\right].$
(3.19)

The corresponding Hamiltonian system leaves
$\operatorname{SO}(n)\times\operatorname{SO}(n)$ invariant and induces on it,
the flow of the symmetric representation of the rigid body system.

Note that the above Hamiltonian is equivalent to the standard rigid body
Hamiltonian $H=\frac{1}{4}\left\langle J^{-1}M,M\right\rangle,$ as in Ratiu
[1980].

###  3.4 Optimality on adjoint orbits and learning

These general ideas can also be extended to a variational problem on the
adjoint orbits of compact Lie groups as in Bloch and Crouch [1996].

Let $\mathfrak{g}$ be a complex semisimple Lie algebra, $\mathfrak{g}_{u}$ its
compact real form, and $G_{u}$ the corresponding compact group. In this case a
natural drift free control system on an orbit of $G_{u}$ takes the form

$\dot{x}=[x,u]$ (3.20)

We remark that we formulate the problem in this generality for convenience,
but the most useful case to bear in mind is the algebra $\mathfrak{su}(n)$ of
skew-Hermitian matrices or the algebra $\mathfrak{so(}n)$ of skew symmetric
matrices (the intersection of the compact and normal real forms of the the
algebra $\mathfrak{sl}(n,\mathbb{C})$). Orbits in this case are similarity
orbits under the group action.

We then consider the following generalization of the functional suggested by
Brockett [1994] (we shall return to Brockett’s precise problem shortly):

$\eta(x,u)=\int_{0}^{t_{f}}1/2||u||^{2}-V(x)dt$ (3.21)

where $||\cdot||=<\cdot,\cdot>^{1/2}$ is the norm induced on
$\mathfrak{g}_{u}$ by the negative of the Killing form $\kappa(\cdot,\cdot)$
on $\mathfrak{g}$ and V is a smooth function on $\mathfrak{g}_{u}$. The
pairing between vectors $x$ in $\mathfrak{g}$ and dual vectors $p$ in
$\mathfrak{g}^{*}$ may be written $<p,x>=-\kappa(x,p)$.

We have

######  Theorem 3.10.

The equations of the maximum principle for the variational problem with
functional 3.21 subject to the dynamics 3.20 are

$\displaystyle\dot{x}$ $\displaystyle=$ $\displaystyle[x,[p,x]]$
$\displaystyle\dot{p}$ $\displaystyle=$ $\displaystyle[p,[p,x]]-V_{x}\,.$
(3.22)

Proof. The Hamiltonian is given by

$H(x,p,u)=<p,[x,u]>-1/2||u||^{2}+V(x)\,.$ (3.23)

Hence

${\partial H\over\partial u}=-<[x,p],\cdot>-<u,\cdot>$

and thus the optimal control is given by

$u^{*}=[p,x]$ (3.24)

Substituting this into $H$ we find the Hamiltonian evaluated along the optimal
trajectory is given by

$H^{*}(p,x)=-1/2<x,[p,[p,x]]>+V(x)$ (3.25)

Computing

$\dot{x}=\left({\partial H^{*}\over\partial p}\right)^{T}$

and

$\dot{p}=-\left({\partial H^{*}\over\partial x}\right)^{T}$

gives the result. $\blacksquare$

A particularly interesting special case of this problem is that of Brockett
[1994] where we have

######  Corollary 3.11.

The equations of the maximum principle for the variational problem 3.21
subject to equations 3.20 with $V(x)=-\tfrac{1}{2}\|[x,n]\|^{2}$ are

$\displaystyle\dot{x}$ $\displaystyle=$ $\displaystyle[x,[p,x]$
$\displaystyle\dot{p}$ $\displaystyle=$ $\displaystyle[p,[p,x]]-[n,[n,x]]\,.$
(3.26)

The proof of the corollary follows immediately, setting
$V(x)=\tfrac{1}{2}\left\langle x,[n,[n,x]]\right\rangle$. Note that with this
functional the equations lie naturally on an adjoint orbit. In addition, these
equations are interesting in that the optimal flow may be related to the
integrable Toda lattice equations (see Bloch [1990], Bloch, Brockett and Ratiu
[1992] and Brockett [1994].)

The smooth machine learning version of this problem considers for smooth
functions $\phi^{a}:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}^{+}$,
$a=1,\dots,M$

$\eta(x,u)=\tfrac{1}{2}\int_{0}^{T}\|u\|^{2}-V(x)dt+\sum_{a=1}^{M}\phi(x^{a}_{T})$
(3.27)

with $\dot{x^{a}}=[x^{a},u],\,a=1\dots M$

As before we now have

######  Corollary 3.12.

The equations of the maximum principle for the variational problem subject to
equations 3.20 with $V(x)=-1/2||[x,n]||^{2}$ are

$\displaystyle\dot{x}^{a}$ $\displaystyle=$
$\displaystyle[x^{a},[p^{a},x^{a}]]$ $\displaystyle\dot{p}^{a}$
$\displaystyle=$ $\displaystyle[p^{a},[p^{a},x]]-[n,[n,x^{a}]]\,.$ (3.28)

with

$p^{a}_{T}=\frac{\partial\phi(x^{a}_{T})}{\partial x^{a}_{T}}$ (3.29)

In particular we would like to consider $\phi(x^{a}_{T})=<x^{a}_{T},n>$.

Then we can see the final solution tends to the double bracket equation
$\dot{x}=[x,[x,n]]$ and the second term in the costate equation goes to zero.

One can then write a discrete version of these equations using an appropriate
discrete algorithm and following our formalism above. This will be considered
in future work.

##  4 Conclusions

In this paper we have discussed discrete optimal control systems and related
them to equations for machine learning. In particular, we have considered the
symmetric formulation of the rigid body equations, both discrete and smooth,
and discussed double bracket equations.

We note also that the analysis here can be extended to other systems, such as
the full Toda dynamics. We intend to discuss such extensions and numerical
aspects, as well as the subRiemannian and quantum settings, in a future
publication.

#### Acknowledgement.

We acknowledge the inspiration and guidance of Roger Brockett over many years
and the wonderful collaborations we had together. We would also like to thank
Fred Leve for his support for the meeting in which Roger was honored and for
his support for research in nonlinear control in general, and to thank Maani
Ghaffari for his valuable comments on the manuscript.

## References

  * Abraham and Marsden [1978] Abraham, R. and J. E. Marsden [1978], Foundations of Mechanics, Second edition. Revised and enlarged, with the assistance of Tudor Ratiu and Richard Cushman, Benjamin/Cummings Publishing Co., Inc., Advanced Book Program, Reading (1978) Addison-Wesley. Reprinted by Perseus Press, 1997.
  * Anderson [1983] Anderson, H. [1983], Rattle: A velocity version of the shake algorithm for molecular dynamics calculations, Journal of Computational Physics, 52, 24–34.
  * Benning et. al. [2019] Benning, M., E. Celledoni, M. Ehrhardt, B. Owren and C. B Schönlieb [2019], Deep learning as optimal control problems: models and numerical methods, J. Computational Dynamics, 6, 171–198.
  * Bloch [1990] Bloch, A. M. [1990], Steepest descent, linear programming and Hamiltonian flows, Contemp. Math. AMS, 114, 77–88.
  * Bloch et al. [1997] Bloch, A.M., R. Brockett and P. Crouch [1997], Double bracket equations and geodesic flows on symmetric spaces, Comm. Math. Phys., 187, 357–373.
  * Bloch, Brockett and Ratiu [1990] Bloch, A. M., R. W. Brockett and T. S. Ratiu [1990], A new formulation of the generalized Toda Lattice equations and their fixed point analysis via the momentum map, Bull. Amer. Math. Soc., 23, 477–485.
  * Bloch, Brockett and Ratiu [1992] Bloch, A. M., R. W. Brockett and T. S. Ratiu [1992], Completely integrable gradient flows, Comm. Math. Phys., 147, 57–74.
  * Bloch and Crouch [1996] Bloch, A.M. and P. Crouch [1996], Optimal control and geodesic flows, Systems Control Lett., 28, n 2, 65–72.
  * Bloch et al. [1998] Bloch, A. M., P. Crouch, J. E. Marsden and T. S. Ratiu [1998], Discrete rigid body dynamics and optimal control, Proc. CDC, 37, 2249–2254.
  * Bloch et al. [2000] Bloch, A. M., P. E. Crouch, J. E. Marsden and T. S. Ratiu [2000], An almost Poisson structure for the generalized rigid body equations, Proc. IFAC Workshop on Lagrangian and Hamiltonian Method in Nonlinear Control (2000), 92–97.
  * Bloch et al. [2002] Bloch, A. M., P. Crouch, J. E. Marsden and T. S. Ratiu [2002], The symmetric representation of the rigid body equations and their discretization, Nonlinearity, 15, 1309–1341.
  * Bloch et al. [2008] Bloch, A. M., P. E. Crouch, J. E. Marsden and A. K. Sanyal [2008], Optimal control and geodesics on matrix Lie groups, Foundations of Computational Mathematics, 8, 469–500.
  * Bloch et al. [1990] Bloch, A. M., H. Flaschka and T. S. Ratiu [1990], A convexity theorem for isospectral manifolds of Jacobi matrices in a compact Lie algebra, Duke Math. J., 61, 41–65.
  * Bloch, Gay-Balmaz and Ratiu [2018] Bloch, A. M., F. Gay-Balmaz and T. S. Ratiu [2016], The Clebsch representation in optimal control and low rank integrable systems, Computation and Combinatorics in Dynamics, Stochastics and Control, The Abel Symposium, 13, Rosendal, Norway, August 2016, Celledoni, E., Di Nunno, G., Ebrahimi-Fard, K., Munthe-Kaas, H. Z. (eds.), 129–158, Springer-Verlag, 2018.
  * Bobenko and Suris [1999] Bobenko, A. I. and Y. B. Suris [1999], Discrete time Lagrangian mechanics on Lie groups with an application to the Lagrange top, Comm. Math. Phys 204, 147–188.
  * Bogayavlenski [1994] Bogayavlenski, O. [1994], Breaking Solitons, Springer Verlag.
  * Brockett [1973] Brockett, R. [1973], Lie theory and control systems defined on spheres, SIAM J. Appl., 122/123/124, 761–777.
  * Brockett [1989] Brockett, R. [1989], Least squares matching problems, Linear Algebra App Math., 25, 213–225.
  * Brockett [1994] Brockett, R. [1994], The double bracket equation as the solution of a variational problem, Hamiltonian and Gradient Flows, Algorithms and Control, Amer. Math. Soc., Providence, RI., 69–76.
  * Bryson and Ho [1969] Bryson, A. E. and Ho, Y.-C, [1969] Applied Optimal Control, Blaisdell Pub.
  * Budd and Iserles [1999] Budd, C. and A. Iserles [1999], Geometric integration: numerical solution of differential equations on manifolds, R. Soc. Lond. Philos. Trans. Ser. A Math. Phys. Eng. Sci., 357, 945–956.
  * Cardoso and Leite [2001] Cardoso, J. and F. Leite [2001], The Moser-Veselov equation, Linear Algebra and its Applications, 360, 237–248.
  * Crouch and Grossman [1993] Crouch, P. and R. Grossman [1993], Numerical integration of ordinary differential equations on manifolds, J. Nonlinear Sci., 3, 1–33.
  * Deift et al. [1992] Deift, P., L. C. Li and C. Tomei [1992], Loop groups, discrete versions of some classical integrable systems, and rank 2 extensions, Mem. Amer. Math. Soc., 100.
  * E [2017] E., W. [2017], A proposal on machine learning via dynamical systems Comm. Math Sci., 5, 1–11.
  * E [2019] E., W. [2017], A mean-field optimal control formulation of deep learning. Research in Mathematical Sciences ., 10. .
  * Faybusovich [1988] Faybusovich, L. E. [1988], Explicitly solvable nonlinear optimal control problems, Internat. J. Control, 48, 2507–2526.
  * Fedorov and Kozlov [1995] Fedorov, Yu. N. and V. V. Kozlov [1995], Various aspects of $n$-dimensional rigid body dynamics, Amer. Math. Soc. Transl., 168, 141–171.
  * Gay-Balmaz, F. and T. S. Ratiu [2011] Gay-Balmaz, F. and T. S. Ratiu [2011], Clebsch optimal control formulation in mechanics, Journal of Geometric Mechanics, 3(1), 41–79.
  * Guibout and Bloch [2004] Guibout, V. and A.M. Bloch [2004], A discrete maximum principle for solving optimal control problems, Proc. 43rd CDC., 43, 1806–1811.
  * Holm et al. [1986] Holm, D. D., J. E. Marsden and T. S. Ratiu [1986], The Hamiltonian structure of continuum mechanics in material, inverse material, spatial, and convective representations, Les Presses de L’Université de Montrèal., Séminaire de Mathématiques supérieurs, 100, 11–122.
  * Iserles et al. [1999] Iserles, A., R. McLachlan and A. Zanna [1999], Approximately preserving symmetries in the numerical integration of ordinary differential equations, European J. Appl. Math., 10, 419–445.
  * Jurdjevic [1997] Jurdjevic, V. [1997], Geometric Control Theory, Cambridge Studies in Advanced Mathematics, 52, Cambridge University Press.
  * Kane et al. [1999] Kane, C, J. E. Marsden and M. Ortiz [1999], Symplectic energy-momentum integrators, J. Math. Phys., 40, 3353–3371.
  * Kane et al. [2000] Kane, C, J. E. Marsden, M. Ortiz and M. West [2000], Variational integrators and the Newmark algorithm for conservative and dissipative mechanical systems, Int. J. Num. Math. Eng., 49, 1295–1325.
  * Kucera [1972] Kucera, V. [1972], A contribution to matrix quadratic equations, IEEE Trans. on Automatic Control, 17, 344-347.
  * Lancaster and Rodman [1980] Lancaster, P. and L. Rodman [1980], Existence and uniqueness theorems for the algebraic Riccati equation, International Journal of Control, 32, 285–309.
  * Lancaster and Rodman [1995] Lancaster, P. and L. Rodman [1995], Algebraic Riccati Equations, Oxford University Press, New York.
  * LeCun [1988] LeCun, Y. [1988], A theoretical framework for back-propagation, Proceedings of the 1988 Connectionist Models Summer School, Lourestky, D., G. Hinton and T. Sejnowski (eds.), 21–28, CMU, Pittsburgh, Morgan Kaufmann.
  * Leimkuhler and Skeel [1994] Leimkuhler, B. and R. Skeel [1994], Symplectic numerical integrators in constrained Hamiltonian systems, Journal of Computational Physics, 112, 117–125.
  * Lewis and Simo [1995] Lewis, D. and J. C. Simo [1995], Conserving algorithms for the dynamics of Hamiltonian systems on Lie groups, J. Nonlinear Sci., 4, 253–299.
  * Lewis, Ratiu, Simo and Marsden [1992] Lewis, D., T. S. Ratiu, J. C. Simo and J. E. Marsden [1992], The heavy top: a geometric treatment, Nonlinearity, 5, 1–48.
  * Manakov [1976] Manakov, S. V. [1976], Note on the integration of Euler’s equations of the dynamics of an $n$-dimensional rigid body, Funct. Anal. and its Appl., 10, 328–329.
  * Marsden et al. [1998] Marsden, J. E., G. W. Patrick and S. Shkoller [1998], Multisymplectic geometry, variational integrators and nonlinear PDEs, Comm. Math. Phys., 199, 351–395.
  * Marsden et al. [1999] Marsden, J. E., S. Pekarsky and S. Shkoller [1999], Discrete Euler–Poincaré and Lie–Poisson equations, Nonlinearity, 12, 1647–1662.
  * Marsden and Ratiu [1999] Marsden, J. E. and T. S. Ratiu [1999], Introduction to Mechanics and Symmetry, Springer-Verlag, Texts in Applied Mathematics, 17, 1994, Second Edition, 1999.
  * McLachlan and Scovel [1995] McLachlan, R. I. and C. Scovel [1995], Equivariant constrained symplectic integration, J. Nonlinear Sci., 5, 233–256.
  * Marsden and Wendlandt [1997] Marsden, J. E. and J. M. Wendlandt [1997], Mechanical systems with symmetry, variational principles and integration algorithms, Current and Future Directions in Applied Mathematics, Alber, M., B. Hu and J. Rosenthal (eds.), Birkhäuser, 219–261.
  * Marsden and West [2001] Marsden, J. E. and M. West [2001], Discrete variational mechanics and variational integrators, Acta Numerica 10, 357–514.
  * Mishchenko and Fomenko [1976] Mishchenko, A. S. and A. T. Fomenko [1976], On the integration of the Euler equations on semisimple Lie algebras, Sov. Math. Dokl., 17, 1591–1593.
  * Mishchenko and Fomenko [1978] Mishchenko, A. S. and A. T. Fomenko [1978], Generalized Liouville method of integration of Hamiltonian systems, Funct. Anal. Appl., 12, 113–121.
  * Moser and Veselov [1991] Moser, J. and A. P. Veselov [1991], Discrete versions of some classical integrable systems and factorization of matrix polynomials, Comm. Math. Phys., 139, 217–243.
  * Ratiu [1980] Ratiu, T. [1980], The motion of the free $n$-dimensional rigid body, Indiana U. Math. J., 29, 609–627.
  * Reich [1994] Reich, S. [1994], Momentum preserving symplectic integrators, Physica D, 76, 375–383.
  * Sanyal, Mosely and Bloch [2009] Sanyal, A. K., C. Mosely and A. M. Bloch [2009], Analytical and numerical solution of a sub-Riemannian optimal control problem with applications of quantum spin systems, Communications in Information and Systems, 9, 59–76.
  * Suris [2001] Suris, Yu B. [2001] Integrable discretizations of some cases of the rigid body dynamics, J. Nonlinear Math. Phys. 8(4), 534–560.
  * Veselov [1988] Veselov, A. P. [1988], Integrable discrete-time systems and difference operators, Funct. Analysis and Appl., 22, 83–94.
  * Wendlandt and Marsden [1997] Wendlandt, J. M. and J. E. Marsden [1997], Mechanical integrators derived from a discrete variational principle, Physica D, 105, 130–162.