# A note on geometric algebras and control problems with SO(3)–symmetries

Jaroslav Hrdina, Aleš Návrat, Petr Vašík and Lenka Zalabová  JH, AN, PV:
Institute of Mathematics, Faculty of Mechanical Engineering, Brno University
of Technology, Technická 2896/2, 616 69 Brno, Czech Republic; LZ: Institute of
Mathematics, Faculty of Science, University of South Bohemia, Branišovská
1760, 370 05 České Budějovice, and Department of Mathematics and Statistics,
Faculty of Science, Masaryk University, Kotlářská 2, 611 37 Brno, Czech
Republic<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

We study the role of symmetries in control systems through the geometric
algebra approach. We discuss two specific control problems on Carnot groups of
step $2$ invariant with respect to the action of $SO(3)$. We understand the
geodesics as the curves in suitable geometric algebras which allows us to
assess a new algorithm for the local control.

###### Key words and phrases:

local control and optimality, Carnot groups, symmetries, sub–Riemannian
geodesics, geometric algebras

###### 2020 Mathematics Subject Classification:

15A67, 53C17

The first three authors were supported by the grant no. FSI-S-20-6187. Fourth
author is supported by the grant no. 20-11473S Symmetry and invariance in
analysis, geometric modeling and control theory from the Czech Science
Foundation. Finally, we thank the referee for valuable comments.

##  1\. Introduction

Geometric control theory uses geometric methods to control various mechanical
systems [13, 5]. We use the methods of sub–Riemannian geometry and Hamiltonian
concept [2, 1]. As a reasonable starting point, we consider mechanisms moving
in the plane, typically wheeled mechanisms like cars (with or without
trailers) or robotic snakes [8, 11]. The movement of a planar mechanisms is
always invariant with respect to the action of the Euclidean group $SE(2)$. As
the prototypes of planar mechanisms, we choose those consisting of the body in
the shape of a triangle and three legs connected to the vertices of the body
by joints of various types and combinations. Although such mechanisms have
almost the same shape, the configuration spaces may differ. In particular,
possible motions of the mechanism induce a specific filtration in the
configuration space. We present two examples that carry the filtration $(3,6)$
and $(4,7),$ respectively [11, 10].

Figure 1. Generalized trident snakes

To control the mechanisms locally, we consider the nilpotent approximations of
the original control systems [4]. Although the configuration spaces and their
approximations have the same filtration, the approximations form Carnot groups
that are generally endowed with more symmetries [15]. One gets the symmetries
generated by the right–invariant vector fields and there may be additional
symmetries acting non–trivially on the distribution. Our Carnot groups of
filtrations $(3,6)$ and $(4,7)$ carry subgroups of the symmetries isomorphic
to $SO(3)$ [12, 10]. This observation leads to the idea of the local control
in geometric algebra approach.

We reformulate the control problems in the concept of geometric algebras
$\mathbb{G}_{3}$ and $\mathbb{G}_{4}$ [17, 18, 20]. We use the natural
$SO(3)$–invariant operations in geometric algebras to reduce the set of
geodesics to a simpler set of curves in the geometric algebra [7]. Namely,
each geodesic is a linear combination of orthogonal vectors, and $SO(3)$ acts
on the geodesics by means of the action on the appropriate orthonormal system
of vectors. So it is sufficient to study geodesics for one fixed orthonormal
basis, i.e. we can study just geodesics in the moduli space over the action of
$SO(3)$.

We present the local control algorithm for finding geodesics passing through
the origin and an arbitrary point in its neighbourhood. The algorithm is based
on the use of rotors in order to relate two orthogonal bases. We provide an
efficient method to such comparison using geometric algebras. We illustrate
our algorithm on two specific examples.

##  2\. Nilpotent control problems

We focus on two control problems such that their symmetry groups contain
$SO(3)$ as their subgroups. The first system has the growth vector $(3,6)$ and
the other one has the growth vector $(4,7)$ [15, 10].

###  2.1. Control problems on Carnot groups of step $2$

By nilpotent control problems we mean the invariant control problems on Carnot
groups and we consider the Carnot groups $G$ of step $2$ with the filtration
$(m,n)$ [1, 16, 14]. If we denote the local coordinates by
$(x,z)\in\mathbb{R}^{m}\oplus\mathbb{R}^{n-m}$, we can model the corresponding
Lie algebra $\mathfrak{g}$ of vector fields

(1)
$\displaystyle\begin{split}X_{i}&={\partial_{x_{i}}}-\frac{1}{2}\sum_{l=1}^{n-m}\sum_{j=1}^{m}c_{ij}^{l}x_{j}{\partial_{z_{l}}}\
\ \ j=1,\dots,m\\\ \ X_{m+j}&={\partial_{z_{j}}}\ \ \hskip
82.51282ptj=1,\dots,m-n,\end{split}$

where $c_{jl}^{k}$ are the structure constants of the Lie algebra
$\mathfrak{g}$ and the symbol $\partial$ stands for partial derivative. We
discuss the related optimal control problem

(2) $\displaystyle\dot{q}(t)=u_{1}X_{1}+\cdots+u_{m}X_{m}$

for $t>0$ and $q$ in $G$ and the control
$u=(u_{1}(t),\dots,u_{m}(t))\in\mathbb{R}^{m}$ with the boundary condition
$q(0)=q_{1},$ $q(T)=q_{2}$ for fixed points $q_{1},q_{2}\in G$, where we
minimize the cost functional ${1\over
2}\int_{0}^{T}(u_{1}^{2}+\cdots+u_{m}^{2})\mathrm{d}t$. The solutions $q(t)$
then correspond to the sub–Riemannian geodesics, i.e. admissible curves
parametrized by a constant speed whose sufficiently small arcs are the length
minimizers.

We use the Hamiltonian approach to this control problem [1]. There are no
strict abnormal extremals for the step $2$ Carnot groups, so we focus on the
normal geodesics and address them just as geodesics. The left–invariant vector
fields $X_{i}$, $i=1,\dots,m$ form a basis of $TG$ and determine the
left–invariant coordinates on $G$. We define the corresponding left–invariant
coordinates $h_{i}$, $i=1,\dots,m$ and $w_{i}$, $i=1,\dots,n-m$ on the fibres
of $T^{*}G$ by $h_{i}(\lambda)=\lambda(X_{i})$ and
$w_{i}(\lambda)=\lambda(X_{m+i})$, for arbitrary $1$–forms $\lambda$ on $G$.
Thus we use $(x_{i},w_{i})$ as the global coordinates on $T^{*}G$.

The geodesics are exactly the projections of normal Pontryagin extremals, i.e.
the integral curves of the left–invariant normal Hamiltonian

(3) $\displaystyle H={1\over 2}(h_{1}^{2}+h_{2}^{2}+\cdots+h_{m}^{2}),$

on $G$. Assume that $\lambda(t)=(x_{i}(t),z_{i}(t),h_{i}(t),w_{i}(t))$ in
$T^{*}G$ is a normal extremal. Then the controls $u_{j}$ to the system (2)
satisfy $u_{j}(t)=h_{j}(\lambda(t))$ and the base system takes the form of

(4) $\displaystyle\begin{split}\dot{x}_{i}&=h_{i},\ \ \ i=1,\dots,m\\\
\dot{z}_{j}&=-\frac{1}{2}\sum_{i=1}^{m}c_{ik}^{j}h_{i}x_{k},\ \ \
j=1,\dots,n-m\end{split}$

for $q=(x_{i},z_{i})$. Using $u_{j}(t)=h_{j}(\lambda(t))$ and the equation
$\dot{\lambda}(t)=\vec{H}(\lambda(t))$ for the normal extremals, we write the
fiber system as

(5)
$\displaystyle\begin{split}\dot{h}_{i}&=-\sum_{l=1}^{m-n}\sum_{j=1}^{m}c_{ij}^{l}h_{j}w_{l},\
\ \ i=1,\dots,m,\\\ \dot{w}_{j}&=0,\ \ \hskip
71.13188ptj=1,\dots,n-m,\end{split}$

where $c_{ij}^{l}$ are the structure constants of the Lie algebra
$\mathfrak{g}$ for the basis $X_{i}$. The solutions $w_{i}$, $i=1,\dots,n-m$
are constants that we denote by

(6) $\displaystyle w_{1}=K_{1},\ \dots,w_{n-m}=K_{n-m}.$

If $K_{1}=\cdots=K_{n-m}=0$ then $h(t)=h(0)$ is a constant and the geodesic
$(x_{i}(t),z_{i}(t))$ is a line in $G$ such that $z_{i}(t)=0$. If at least one
of $K_{i}$ is non–zero, the first part of the fibre system (5) forms a
homogeneous system of ODEs $\dot{h}=-\Omega h$ with constant coefficients for
$h=(h_{1},\dots,h_{m})^{T}$ and the system matrix $\Omega$. Its solution is
given by $h(t)=e^{-t\Omega}h(0)$, where $h(0)$ is the initial value of the
vector $h$ at the origin.

###  2.2. Left–invariant control problem with the growth vector $(3,6)$

Let us consider three vector fields on $\mathbb{R}^{6}$ with the local
coordinates $(x_{1},x_{2},x_{3},z_{1},z_{2},z_{3})$ in the form

(7) $\displaystyle\begin{split}X_{1}&=\partial_{x_{1}}+{x_{3}\over
2}\partial_{z_{2}}-{x_{2}\over 2}\partial_{z_{3}},\\\
X_{2}&=\partial_{x_{2}}+{x_{1}\over 2}\partial_{z_{3}}-{x_{3}\over
2}\partial_{z_{1}},\\\ X_{3}&=\partial_{x_{3}}+{x_{2}\over
2}\partial_{z_{1}}-{x_{1}\over 2}\partial_{z_{2}}.\end{split}$

The only non–trivial Lie brackets are

(8) $\displaystyle\begin{split}X_{4}=[X_{1},X_{2}]=\partial_{z_{3}},\ \
X_{5}=[X_{1},X_{3}]=-\partial_{z_{2}},\ \
X_{6}=[X_{2},X_{3}]=\partial_{z_{1}}.\end{split}$

These six vector fields determine a step $2$ nilpotent Lie algebra
$\mathfrak{m}$ with the multiplication table given by Table 1.

$\mathfrak{m}$ | $X_{1}$ | $X_{2}$ | $X_{3}$ | $X_{4}$ | $X_{5}$ | $X_{6}$
---|---|---|---|---|---|---
$X_{1}$ | $0$ | $X_{4}$ | $X_{5}$ | $0$ | $0$ | $0$
$X_{2}$ | $-X_{4}$ | $0$ | $X_{6}$ | $0$ | $0$ | $0$
$X_{3}$ | $-X_{5}$ | $-X_{6}$ | $0$ | $0$ | $0$ | $0$
$X_{4}$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$
$X_{5}$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$
$X_{6}$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$
Table 1. Lie algebra $\mathfrak{m}$

There is a Carnot group $M$ such that the fields $X_{i}$, $i=1\dots,6$ are
left–invariant for the corresponding group structure. When identified with
$\mathbb{R}^{6}=\mathbb{R}^{3}\oplus\mathbb{R}^{3}$, the group structure on
$M$ reads as

(9)
$\displaystyle(x,z)\cdot(x^{\prime},z^{\prime})=(x+x^{\prime},z+z^{\prime}+{1\over
2}x\times x^{\prime})$

for $x=(x_{1},x_{2},x_{3})$ and $z=(z_{1},z_{2},z_{3})$, where $\times$ stands
for the vector product on $\mathbb{R}^{3}$. In particular,
$\mathcal{M}=\langle X_{1},X_{2},X_{3}\rangle$ forms a $3$–dimensional
left–invariant distribution on $M$. We define the left–invariant
sub–Riemannian metric $g_{M}$ on $\mathcal{M}$ by declaring
$X_{1},X_{2},X_{3}$ orthonormal.

The geodesics of the control problem are the solutions to the control system
(4),(5), with $(m,n)=(3,6)$ and one can read the structure constants in Table
1. Hence, the fibre system is given by $w_{1}=K_{1},w_{2}=K_{2},w_{3}=K_{3}$
for the constants $K_{1},K_{2},K_{3}$ and $\dot{h}=-\Omega h$ for
$h=(h_{1},h_{2},h_{3})^{T}$ and

(10) $\displaystyle\Omega=\left(\begin{smallmatrix}0&K_{1}&K_{2}\\\
-K_{1}&0&K_{3}\\\ -K_{2}&-K_{3}&0\end{smallmatrix}\right).$

Its solution is given by the exponential $h(t)=e^{-t\Omega}h(0)$, where $h(0)$
is the initial value of the vector $h$ at the origin. We write an explicit
formula for the general solution in terms of eigenvectors of (10). If at least
one of the constants $K_{i}$ is non–zero, the kernel of $\Omega$, i.e.
zero–eigenspace, is one–dimensional, generated by the vector
$(K_{3},K_{2},K_{1})^{T}$. Its orthogonal complement corresponds to the sum of
eigenspaces to the eigenvalues $\pm iK$, where
$K:=\sqrt{K_{1}^{2}+K_{2}^{2}+K_{3}^{2}}$, and is generated by the vectors
$(-K_{1}K_{3},-K_{1}K_{2},K_{2}^{2}+K_{3}^{2})\pm i(K_{2},-K_{3},0)$. Thus
solution to the fibre system can be written as

(11) $\displaystyle
h(t)=(C_{1}\cos(Kt)-C_{2}\sin(Kt))v_{1}+(C_{1}\sin(Kt)+C_{2}\cos(Kt))v_{2}+C_{3}v_{3},$

where $v_{1},v_{2},v_{3}$ is the eigenspace–adapted real orthonormal basis

$v_{1}=\frac{1}{K\sqrt{K_{2}^{2}+K_{3}^{2}}}\begin{pmatrix}-K_{1}K_{3}\\\
K_{1}K_{2}\\\
K_{2}^{2}+K_{3}^{2}\end{pmatrix},v_{2}=\frac{1}{\sqrt{K_{2}^{2}+K_{3}^{2}}}\begin{pmatrix}-K_{2}\\\
-K_{3}\\\ 0\end{pmatrix},v_{3}=\frac{1}{K}\begin{pmatrix}K_{3}\\\ -K_{2}\\\
K_{1}\end{pmatrix}$

and $C_{1},C_{2},C_{3}$ are the constants that satisfy the level set condition
$H=1/2$, i.e. $\|h(t)\|=1$, that reads $C_{1}^{2}+C_{2}^{2}+C_{3}^{2}=1.$ Let
us note that the choice $C_{1}=C_{2}=0$ leads to the constant solutions that
are irrelevant as the control functions. Thus we assume that at least one of
the constants $C_{1},C_{2}$ is non–zero.

Let us emphasize that the base system (4) can be written in terms of a vector
product as follows

(12) $\displaystyle\begin{split}&\dot{x}=h,\\\ &\dot{z}={1\over 2}x\times
h\end{split}$

for vectors $x=(x_{1},x_{2},x_{3})^{T}$ and $z=(z_{1},z_{2},z_{3})^{T}$. One
obtains the general solution by substituting (11) for $h$ and by consequent
direct integration. We are interested in the solutions passing through the
origin, i.e. we impose the initial condition

(13) $\displaystyle x_{i}(0)$ $\displaystyle=0,z_{i}(0)=0,\quad i=1,2,3.$

However, it may be difficult to find the integration constants giving the
geodesics through a fixed target point.

###  2.3. Left–invariant control problem with the growth vector $(4,7)$

Let us consider four vector fields on $\mathbb{R}^{7}$ with the local
coordinates $(x,\ell_{1},\ell_{2},\ell_{3},y_{1},y_{2},y_{3})$ in the form

(14) $\displaystyle\begin{split}&Y_{0}=\partial_{x}-{\ell_{1}\over
2}{\partial_{y_{1}}}-{\ell_{2}\over 2}{\partial_{y_{2}}}-{\ell_{3}\over
2}{\partial_{y_{3}}},\\\ &Y_{1}=\partial_{\ell_{1}}+{x\over
2}\partial_{y_{1}},\ \ Y_{2}=\partial_{\ell_{2}}+{x\over 2}\partial_{y_{2}},\
\ Y_{3}=\partial_{\ell_{3}}+{x\over 2}\partial_{y_{3}}.\end{split}$

The only non–trivial Lie brackets are

(15) $\displaystyle\begin{split}Y_{4}=[Y_{0},Y_{1}]=\partial_{y_{1}},\ \
Y_{5}=[Y_{0},Y_{2}]=\partial_{y_{2}},\ \
Y_{6}=[Y_{0},Y_{3}]=\partial_{y_{3}}\end{split}.$

These seven fields determine a step $2$ nilpotent Lie algebra $\mathfrak{n}$
with the multiplication table given by Table 2.

$\mathfrak{n}$ | $Y_{0}$ | $Y_{1}$ | $Y_{2}$ | $Y_{3}$ | $Y_{4}$ | $Y_{5}$ | $Y_{6}$
---|---|---|---|---|---|---|---
$Y_{0}$ | $0$ | $Y_{4}$ | $Y_{5}$ | $Y_{6}$ | $0$ | $0$ | $0$
$Y_{1}$ | $-Y_{4}$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$
$Y_{2}$ | $-Y_{5}$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$
$Y_{3}$ | $-Y_{6}$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$
$Y_{4}$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$
$Y_{5}$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$
$Y_{6}$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$
Table 2. Lie algebra $\mathfrak{n}$

There is a Carnot group $N$ such that the fields $Y_{i}$, $i=1,\dots,7$ are
left–invariant for the corresponding group structure. The group structure on
$N$, when identified with
$\mathbb{R}^{7}=\mathbb{R}\oplus\mathbb{R}^{3}\oplus\mathbb{R}^{3}$, yields

(16)
$\displaystyle(x,\ell,y)\cdot(x^{\prime},\ell^{\prime},y^{\prime})=(x+x^{\prime},\ell+\ell^{\prime},y+y^{\prime}+{1\over
2}\ell\times\ell^{\prime}).$

for $\ell=(\ell_{1},\ell_{2},\ell_{3})$ and $y=(y_{1},y_{2},y_{3})$. In
particular, $\mathcal{N}=\langle Y_{0},Y_{1},Y_{2},Y_{3}\rangle$ forms a
$4$–dimensional left–invariant distribution on $N$. Moreover, there is a
natural decomposition

(17) $\displaystyle\mathcal{N}=\langle Y_{0}\rangle\oplus\langle
Y_{1},Y_{2},Y_{3}\rangle$

into a $1$–dimensional distribution and a $3$–dimensional involutive
distribution, both left–invariant. We define the left–invariant sub–Riemannian
metric $g_{N}$ on $\mathcal{N}$ by declaring $Y_{0}$, $Y_{1}$, $Y_{2}$,
$Y_{3}$ orthonormal.

The geodesics of the control problem are solutions to the control system
(4),(5), with $(m,n)=(4,7)$ and we read the structure constants in Table 2.
Hence, the first part of the fibre system (5) is given by
$w_{1}=K_{1},w_{2}=K_{2},w_{3}=K_{3}$, where $K_{1},K_{2},K_{3}$ are
constants. Second part of the fibre system takes form $\dot{h}=-\Omega h,$
where $h:=(h_{0},h_{1},h_{2},h_{3})^{T}$ and

(18) $\displaystyle\Omega=\left(\begin{smallmatrix}0&K_{1}&K_{2}&K_{3}\\\
-K_{1}&0&0&0\\\ -K_{2}&0&0&0\\\ -K_{3}&0&0&0\end{smallmatrix}\right).$

Its solution is given by $h(t)=e^{-t\Omega}h(0)$, where $h(0)$ is the initial
value of the vector $h$ at the origin and we write its explicit form in terms
of the eigenvectors of (18). If $K_{1}=K_{2}=K_{3}=0$ then $h(t)=h(0)$ is a
constant and the geodesic $(x(t),\ell_{i}(t),y_{i}(t))$ is a line in $N$ such
that $y_{i}=0$. If at least one of the constants $K_{i}$ is non–zero, the
kernel of $\Omega$, i.e. zero–eigenspace, is two–dimensional and is generated
by the vectors $(0,-K_{3},0,K_{1})^{T}$ and $(0,-K_{2},K_{1},0)^{T}$. Its
orthogonal complement corresponds to the sum of the eigenspaces to the
eigenvalues $\pm iK$, where $K:=\sqrt{K_{1}^{2}+K_{2}^{2}+K_{3}^{2}}$, and is
generated by the eigenvectors $(0,K_{1},K_{2},K_{3})^{T}\pm i(K,0,0,0)^{T}$.
Thus the solution to the vertical system for non–zero $K$ takes form

(19) $\displaystyle\begin{split}h_{0}&=K(C_{2}\cos(Kt)-C_{1}\sin(Kt))\\\
\bar{h}&=K(C_{2}\sin(Kt)+C_{1}\cos(Kt))r_{1}+Cr_{2}\end{split}$

where $\bar{h}=(h_{1},h_{2},h_{3})^{T}$ and $r_{1},r_{2}$ are
eigenspace–adapted real orthonormal vectors

$r_{1}=\frac{1}{K}\left(\begin{smallmatrix}K_{1}\\\ K_{2}\\\
K_{3}\end{smallmatrix}\right),\ \ \
r_{2}=\frac{1}{C}\left(C_{3}\left(\begin{smallmatrix}-K_{3}\\\ 0\\\
K_{1}\end{smallmatrix}\right)+C_{4}\left(\begin{smallmatrix}-K_{2}\\\ K_{1}\\\
0\end{smallmatrix}\right)\right)$

with the constants $C_{1},C_{2},C_{3},C_{4}$ and the normalization factor

$C=\sqrt{(C_{3}K_{3}+C_{4}K_{2})^{2}+K_{1}^{2}(C_{3}^{2}+C_{4}^{2})}.$

The level set condition $\|h(t)\|=1$ reads $C_{1}^{2}+C_{2}^{2}+C_{3}^{2}=1.$
Let us note that the choice $C_{1}=C_{2}=0$ leads to the constant solutions
that are irrelevant as the control functions. Thus we assume that at least one
of the constants $C_{1},C_{2}$ is non–zero.

The base system (4) takes the explicit form of

(20) $\displaystyle\begin{split}&\dot{x}=h_{0},\\\ &\dot{\ell}=\bar{h},\\\
&\dot{y}={1\over 2}(x\bar{h}-h_{0}\ell).\end{split}$

We are interested in the solutions passing through the origin, i.e. we impose
the initial condition

(21) $\displaystyle x(0)=0,\ell_{i}(0)=0,y_{i}(0)=0,\quad i=1,2,3.$

By substitution of (19), the system (LABEL:xl) can be directly integrated.
Again, it may be difficult to find the geodesics through a fixed target point.
In the section 4 we show how the symmetries of the system and the geometric
algebra approach is used for finding a geodesic towards a given point.

###  2.4. Symmetries of the control systems

Symmetries of the control system in question coincide with the symmetries of
the corresponding left–invariant sub–Riemannian structure
$(M,\mathcal{M},g_{M})$ and $(N,\mathcal{N},g_{N})$, respectively. These are
precisely the automorphisms on groups preserving the distributions and
sub–Riemannian metrics. The group $SO(3)$ acts on $\mathbb{R}^{3}$ and
preserves the vector product which implies the following statement.

######  Proposition 1.

For each $R\in SO(3)$, the map

(22) $\displaystyle\begin{split}(x,z)\mapsto(Rx,Rz)\end{split}$

maps the geodesics of the system from Section 2.2 starting at the origin to
the geodesics starting at the origin. For each $R\in SO(3)$ the map

(23) $\displaystyle\begin{split}(x,\ell,y)\mapsto(x,R\ell,Ry)\end{split}$

maps the geodesics of the system from Section 2.3 starting at the origin to
the geodesics starting at the origin.

###### Proof.

Follows from the invariance of (LABEL:rr1) and (LABEL:xl) with respect to the
action of $R\in SO(3)$. ∎

##  3\. Geometric algebra

The construction of the universal real geometric algebra is well-known [17,
18, 20, 9]. We provide only a brief description in a special case
$\mathbb{G}_{m}$ that we use later. In general, geometric algebras are based
on symmetric bilinear forms of arbitrary signature. Here, we deal with the
real vector space $\mathbb{R}^{m}$ endowed with a positive definite symmetric
bilinear form $B$ only.

###  3.1. Geometric product

Let us consider a positive definite symmetric bilinear form $B$ on
$\mathbb{R}^{m}$ and the associated orthonormal basis $(e_{1},\dots,e_{m})$ ,
i.e.

$\displaystyle B(e_{i},e_{j})=\begin{cases}1&\text{ if }i=j\\\ 0&\text{ if
}i\neq j\end{cases}\;\text{ where }1\leq i,j\leq m.$

The Grassmann algebra $\Lambda(\mathbb{R}^{m})$ is an associative algebra with
the anti-symmetric outer product $\wedge$ defined by the rule

$\displaystyle e_{i}\wedge e_{j}+e_{j}\wedge e_{i}=0\;\text{ for }1\leq
i,j\leq m.$

The Grassmann blade of grade $r$ is $e_{A}=e_{i_{1}}\wedge\cdots\wedge
e_{i_{r}}$, where the multi-index $A$ is a set of indices ordered in the
natural way $1\leq i_{1}\leq\cdots\leq i_{r}\leq m,$ and we put
$e_{\emptyset}=1.$ Blades of grades $0\leq r\leq m$ form the basis of the
Grassmann algebra $\Lambda(\mathbb{R}^{m})$ and for the outer product we have

$\displaystyle e_{j}\wedge e_{A}=\begin{cases}e_{j}\wedge
e_{i_{1}}\wedge\cdots\wedge e_{i_{r}}&\text{ if }j\notin A,\\\ 0&\text{ if
}j\in A\end{cases}$

and $1\wedge e_{A}=e_{A}$. For the vectors from $\mathbb{R}^{m}$, the inner
product $e_{i}\cdot e_{j}=B(e_{i},e_{j})$ and the outer product $e_{i}\wedge
e_{j}$ lead to the so-called geometric product

$\displaystyle e_{i}e_{j}=e_{i}\cdot e_{j}+e_{i}\wedge e_{j},\;\;1\leq i,j\leq
m.$

The definitions of inner and geometric products then extend to blades of the
grade $r$ as follows. For the inner product we put $1\cdot e_{A}=0$ and

$\displaystyle e_{j}\cdot e_{A}=e_{j}\cdot(e_{i_{1}}\wedge\cdots\wedge
e_{i_{r}})=\sum_{k=1}^{r}(-1)^{k}B(e_{j},e_{i_{k}})e_{A\setminus\\{i_{k}\\}},$

where $e_{A\setminus\\{i_{k}\\}}$ is the blade of grade $r-1$ created by
deleting $e_{i_{k}}$ from $e_{A}$. This product is also called the left
contraction in literature. For the geometric product we define

$\displaystyle e_{j}e_{A}=e_{j}\cdot e_{A}+e_{j}\wedge e_{A}.$

These definitions extend linearly to the whole vector space
$\Lambda(\mathbb{R}^{m})$. Thus we get an associative algebra over this vector
space, the so-called real geometric algebra, denoted by $\mathbb{G}_{m}$. Note
that this algebra is naturally graded; the grade zero and grade one elements
are identified with $\mathbb{R}$ and $\mathbb{R}^{m}$, respectively. Finally,
we can define the norm of a blade as the magnitude of the blade
$|e_{A}|=\sqrt{e_{A}\cdot\tilde{e}_{A}}$. Note that $e_{A}\cdot\tilde{e}_{A}$,
where $e_{A}\neq 0$ is always positive in $\mathbb{G}_{m}$.

###  3.2. Objects

The vectors in $\mathbb{R}^{m}$ with the coordinates $(x_{1},\dots,x_{m})$ are
given by $x=x_{1}e_{1}+\cdots+x_{m}e_{m}$ and the square with respect to the
geometric product $x^{2}=x_{1}^{2}+\cdots+x_{m}^{2}\in\mathbb{R}$ coincides
with the square of the Euclidean norm of $x$. A vector $x$ represents a one-
dimensional subspace (line) $p$ in $\mathbb{R}^{m}$ given by the scalar
multiples of $x$ which in $\mathbb{G}_{m}$ is expressed by the formula $u\in
p\Leftrightarrow u\wedge x=0$. In the same way, a plane $\pi$ generated by two
vectors $x$ and $y$ is represented by $x\wedge y$ in the sense
$u\in\pi\Leftrightarrow u\wedge x\wedge y=0$. In general, any $r$-dimensional
subspace $V_{r}\subseteq\mathbb{R}^{m}$ is represented by a blade $A_{r}$ of
grade $r$ such that

(24) $\displaystyle V_{r}=NO(A_{r})=\\{x\in\mathbb{R}^{m}:x\wedge A_{r}=0\\}.$

Such a representation is called the outer product null space (OPNS)
representation in the literature. In particular, the whole space
$\mathbb{R}^{m}$ is represented by a blade of maximal grade, so called
pseudoscalar. Similarly one defines the inner product null space (IPNS)
representation $A^{*}_{m-r}$ of $V_{r}$ as a blade of grade $m-r$ such that
$x\in V_{r}\Leftrightarrow x\cdot A^{*}_{m-r}=0$. The OPNS and IPNS
representations are mutually dual with respect to the duality on
$\mathbb{G}_{m}$ defined by the multiplication by pseudoscalar, namely

$A^{*}=AI,$

where $A$ is a blade and $I$ is the pseudoscalar. Indeed, one can show that
$(x\wedge A)I=x\cdot(AI)$ for each vector $x\in\mathbb{R}^{m}$, in particular

$x\wedge A=0\Leftrightarrow x\cdot A^{*}=0.$

######  Remark 1.

OPNS representations of objects in $\mathbb{G}_{3}$ are summarized in Table 3.
For example, the plane generated by the vectors $u$, $v$ has the OPNS
representation $u\wedge v$. Its IPNS representation $(u\wedge v)^{*}$ is a
vector perpendicular to the plane. More specifically, for the pseudoscalar
$I=e_{1}\wedge e_{2}\wedge e_{3}=e_{1}e_{2}e_{3}$ we receive the usual vector
product in geometric algebra form as

(25) $\displaystyle u\times v=-(u\wedge v)I.$

grade | name | blades | dimension | objects
---|---|---|---|---
0 | scalars | 1 | 1 | numbers
1 | vectors | $e_{1},\>e_{2},\>e_{3}$ | 3 | lines
2 | bivectors | $e_{1}\wedge e_{2},\>e_{1}\wedge e_{3},\>e_{2}\wedge e_{3}$ | 3 | planes
3 | pseudoscalars | $e_{1}\wedge e_{2}\wedge e_{3}$ | 1 | volume forms
Table 3. Blades of geometric algebra $\mathbb{G}_{3}$

###  3.3. Transformations

Let us fix a vector $n\in\wedge^{1}\mathbb{R}^{m}\subset\mathbb{G}_{m}$ such
that $n\cdot n=n^{2}=1$ and
$x\in\wedge^{1}\mathbb{R}^{m}\subset\mathbb{G}_{m}$ arbitrary. The negative
conjugation $-nxn$ defines the reflection with respect to the hyperplane
orthogonal to $n$, because

$\displaystyle-nx^{\perp}n$ $\displaystyle=-(n\wedge
x^{\perp})n=(x^{\perp}\wedge n)n=x^{\perp}nn=x^{\perp},$ $\displaystyle-
nx^{||}n$ $\displaystyle=-x^{||},$

where $x=x^{||}+x^{\perp}$ is the orthogonal decomposition of $x$ with respect
to $n$. Conjugation preserves the grades of blades and is an outermorphism
$n(u_{1}\wedge\cdots\wedge u_{l})n=(nu_{1}n)\wedge\cdots\wedge(nu_{l}n)$ for
any vectors $u_{1},\dots,u_{l}$, thus minus the conjugation is the
(anti)outermorphism depending on the dimension $m$. Since each rotation is a
composition of two reflections, a rotation in $\mathbb{G}_{m}$ is represented
by the conjugation with respect to the geometric product of two vectors. To
find a rotor between vectors $x$ and $y$ we have a nice formula at hand.

######  Lemma 1.

Let $x$ and $y$ be the unit vectors in $\mathbb{G}_{m}$, i.e.
$x,y\in\wedge^{1}\mathbb{G}_{m}$, then the formula

(26) $\displaystyle R_{xy}=\widehat{1+yx},$

where the hat symbol stands for the normalization
$\widehat{u}={u}/{\sqrt{u\cdot u}},$ defines the rotation in the plane
$x\wedge y$ which maps vector $x$ to $y$ and acts trivially on $(x\wedge
y)^{*}$.

###### Proof.

Multiplication of two vectors
$x,y\in\wedge^{1}\mathbb{R}^{m}\subset\mathbb{G}_{m}$, such that
$x^{2}=y^{2}=1,$ defines a multivector

$yx=\cos(\theta)+\sin(\theta)\widehat{y\wedge x},$

where $\theta$ is the angle between $x$ and $y$[18, Lemma 4.2]. The
conjugation by such multivector $yx$ represents the rotation in the plane
$x\wedge y$ with respect to angle $2\theta$ in the positive way. Using
standard trigonometric formulas we can see by straightforward calculation that

$\displaystyle\begin{split}R_{xy}&=\widehat{1+yx}=\frac{1+\cos(\theta)+\widehat{y\wedge
x}\sin(\theta)}{\sqrt{(1+\cos(\theta))^{2}+\sin^{2}(\theta)}}=\frac{1+\cos(\theta)+\widehat{y\wedge
x}\sin(\theta)}{\sqrt{2+2\cos(\theta)}}\\\
&=\sqrt{\frac{1+\cos(\theta)}{2}}+\widehat{y\wedge
x}\sqrt{\frac{1-\cos^{2}(\theta)}{2(1+\cos(\theta))}}=\sqrt{\frac{1+\cos(\theta)}{2}}+\widehat{y\wedge
x}\sqrt{\frac{1-\cos(\theta)}{2}}\\\ &=\cos(\frac{\theta}{2})+(y\wedge
x)\sin(\frac{\theta}{2}).\end{split}$

So $R_{xy}$ is the rotation in the plane $\widehat{x\wedge y}$ in the positive
way about the angle between $x$ and $y$ so the vector $x$ goes to the vector
$y$.

Finally, $(x\wedge y)^{*}$ is orthogonal to $x$ and $y$ and the
straightforward computation $xy(x\wedge y)^{*}yx=xyyx(x\wedge y)^{*}=(x\wedge
y)^{*}$ proves the rest of the statement. ∎

######  Remark 3.4.

One can see that $yxy$ is an axial symmetry with respect to $y$ and thus
$x+yxy=2\cos(\theta)y$. We can compute the square of the norm

$\displaystyle(1+yx)(1+xy)=1+xy+yx+1=2+2\cos(\theta)$

and with the help of the geometric product, we compute

$\displaystyle(1+yx)x(1+xy)$
$\displaystyle=(x+y)(1+xy)=(x+2y+yxy)=(2+2\cos(\theta))y$

so the conjugation by (26) maps $x$ to $y$.

###  3.5. Rotor construction

Let $(x_{1},\dots,x_{m})$ and $(y_{1},\dots,y_{m})$ be a pair of bases of
$\mathbb{R}^{m}$ such that

  1. (1)

$x_{i}\cdot x_{j}=y_{i}\cdot y_{j}$ for all $i,j=1,\dots,m$, i.e. all the
scalar products are equal, and

  2. (2)

$x_{1}\wedge\cdots\wedge x_{m}=y_{1}\wedge\cdots\wedge y_{m}$, i.e. the
pseudoscalars are equal.

Let us remind that a complete flag $\\{V\\}$ in an increasing sequence of
subspaces of the vector space $\mathbb{R}^{m}$

$\\{0\\}\subset V_{1}\subset V_{2}\subset\cdots\subset V_{m}=\mathbb{R}^{m},$

such that $\text{dim}(V_{i})=i$. We use the complete flags to find the
explicit rotation $R$ such that $Rx_{i}\tilde{R}=y_{i}$ for all $i=1,\dots m$.
Our method can be summarized as follows:

  * •

We consider the complete flags $\\{V\\}$ and $\\{W\\}$ by setting
$V_{i}=\langle x_{1},\dots,x_{i}\rangle=NO(x_{1}\wedge\cdots\wedge x_{i})$ and
$W_{i}=\langle y_{1},\dots,y_{i}\rangle=NO(y_{1}\wedge\cdots\wedge y_{i})$,
respectively.

  * •

We map the complete flag $\\{V\\}$ to the complete flag $\\{W\\}$ inductively
in $m$ steps. In the $j^{\text{th}}$ step, we assume $V_{i}=W_{i}$ for $i>j$
and we find the rotation $R_{i}$ such that $R_{i}V_{i}\tilde{R}_{i}=W_{i}$ for
$i>j-1$.

Before we formulate the construction in detail, we need several technical
lemmas.

######  Lemma 2.

Let $(x_{1},\dots,x_{m})$ and $(y_{1},\dots,y_{m})$ be a pair of bases such
that

  1. (1)

$x_{i}\cdot x_{j}=y_{i}\cdot y_{j}$ for all $i,j=1,\dots,m$, and

  2. (2)

$x_{1}\wedge\cdots\wedge x_{i}=y_{1}\wedge\cdots\wedge y_{i}$ for all
$i=1,\dots,m$.

If $\\{V\\}$ and $\\{W\\}$ are the corresponding complete flags, respectively,
and if $V_{i}=W_{i}$ for all $i=1,\dots,m,$ then $x_{i}=y_{i}$ for all
$i=1,\dots,m$.

###### Proof.

The equality $x_{1}=y_{1}$ holds trivially from the assumptions. Then
$x_{2}\cdot x_{2}=y_{2}\cdot y_{2}$ reads that $x_{2},y_{2}$ are of the same
length, $x_{2}\cdot x_{1}=y_{2}\cdot y_{1}$ reads that the angles between
$x_{1},x_{2}$ and $y_{1},y_{2}$ are identical and $x_{1}\wedge
x_{2}=y_{1}\wedge y_{2}$ reads that they have the same orientation. Then
$V_{2}=W_{2}$ and $x_{1}=y_{1}$ imply $x_{2}=y_{2}$ and so on for all the
basis vectors. ∎

######  Lemma 3.

Let $(x_{1},\dots,x_{i},z)$ and $(y_{1},\dots,y_{i},z)$ be a pair of sets of
independent vectors such that $x_{1}\wedge\cdots\wedge x_{i}\wedge
z=y_{1}\wedge\cdots\wedge y_{i}\wedge z$. If $NO(x_{1}\wedge\cdots\wedge
x_{i})=NO(x_{1}\wedge\cdots\wedge x_{i})$ then $x_{1}\wedge\cdots\wedge
x_{i}=y_{1}\wedge\cdots\wedge y_{i}$.

###### Proof.

The independence of the sets of vectors imply $x_{1}\wedge\cdots\wedge
x_{i}\wedge z\neq 0$, $y_{1}\wedge\cdots\wedge y_{i}\wedge z\neq 0$. If
$NO(x_{1}\wedge\cdots\wedge x_{i})=NO(y_{1}\wedge\cdots\wedge y_{i})$ then
$x_{1}\wedge\cdots\wedge x_{i}=\beta x_{1}\wedge\cdots\wedge x_{i}$ for
$\beta\in\mathbb{R}$.

Then

$\displaystyle x_{1}\wedge\cdots\wedge x_{i}\wedge z$
$\displaystyle=y_{1}\wedge\cdots\wedge y_{i}\wedge z$
$\displaystyle(x_{1}\wedge\cdots\wedge x_{i}-y_{1}\wedge\cdots\wedge
y_{i})\wedge z$ $\displaystyle=0$
$\displaystyle(1-\beta)x_{1}\wedge\cdots\wedge x_{i}\wedge z$
$\displaystyle=0$

Finally $\beta=1$ and $x_{1}\wedge\cdots\wedge x_{i}=y_{1}\wedge\cdots\wedge
y_{i}$. ∎

######  Lemma 4.

Consider two complete flags $\\{V\\}$ and $\\{W\\}$ in $\mathbb{R}^{m}$ and
$i\leq m$ such that $V_{j}=W_{j}$ for $j>i$. The rotor $R_{i}$ between the
hyperplanes $V_{i}\oplus V_{i+1}^{\perp}$ and $W_{i}\oplus W_{i+1}^{\perp}$
constructed by the formula (26) maps $V_{i}$ to $W_{i}$.

###### Proof.

The property $V_{i}\subset V_{i+1}$ implies that $V_{i+1}^{\perp}\subset
V_{i}^{\perp}$ and thus $V_{i}\oplus V_{i+1}^{\perp}$ is a hyperplane equipped
with an orthogonal decomposition. Recall that $V_{i+1}=W_{i+1}$, and thus
$V^{\perp}_{i+1}=W^{\perp}_{i+1}$. Any rotation preserves an orthogonal
decomposition and thus $R_{i}$ acts as the identity on
$V^{\perp}_{i+1}=W^{\perp}_{i+1}$ , so it maps $V_{i}$ to $W_{i}$. ∎

We use all these Lemmas to provide a constructive proof of the following
theorem.

######  Theorem 3.6.

Let $(x_{1},\dots,x_{m})$ and $(y_{1},\dots,y_{m})$ be a pair of bases of
$\mathbb{R}^{m}$ such that

  1. (1)

$x_{i}\cdot x_{j}=y_{i}\cdot y_{j}$ for all $i,j=1,\dots,m$, and

  2. (2)

$x_{1}\wedge\cdots\wedge x_{m}=y_{1}\wedge\cdots\wedge y_{m}$.

Then we can construct a rotor $R$ such that $Rx_{i}\tilde{R}=y_{i}$ for all
$j,i=1,\dots m$.

###### Proof.

Let $\\{V\\}$ and $\\{W\\}$ be a pair of the corresponding complete flags
$V_{i}=NO(x_{1}\wedge\cdots\wedge x_{i})$ and
$W_{i}=NO(y_{1}\wedge\cdots\wedge y_{i})$. We construct a rotor $R=R_{1}\cdots
R_{m}$ mapping the complete flag $\\{V\\}$ to the complete flag $\\{W\\}$ so
that $V_{i}=W_{i}$ for all $i=1,\dots m$. The result on bases $x_{i},y_{i}$
then follows by Lemmas 2 and 3.

We define $R_{m}$ as the identity and proceed inductively. It follows from
Lemma 1 that there is a rotation $R_{m-1}$ between the hyperplanes
$V_{m-1}\oplus V_{m}^{\perp}\cong V_{m-1}$ and $W_{m-1}\oplus
W_{m}^{\perp}\cong W_{m-1}$ which maps the complete flag $\\{V\\}$ to the
complete flag $\\{RV\tilde{R}\\}$ in such a way that
$W_{m-1}=RV_{m-1}\tilde{R}$, where $R=R_{m-1}R_{m}$.

As the induction step, we consider the rotor $R=R_{j}\cdots R_{m}$ such that
$RV_{i}\tilde{R}=W_{i}$ for all indices $i\geq j$. According to Lemma 4, the
rotation $R_{j-1}$ between the hyperplanes
$(RV_{j-1}\tilde{R})\oplus(RV_{j}\tilde{R})^{\perp}$ and $W_{j-1}\oplus
W_{j}^{\perp}$ maps the complete flag $\\{RV\tilde{R}\\}$ to the complete flag
$\\{R_{j}RV\tilde{R}\tilde{R}_{j}\\}$ in such a way that
$W_{i}=R_{j}RV_{i}\tilde{R}\tilde{R}_{j}$ for all $i\geq j-1$.

After $m$ steps, the rotor $R=R_{1}\cdots R_{m}$ maps the complete flag
$\\{V\\}$ to the complete flag $\\{W\\}$ in such a way that $V_{i}=W_{i}$ for
all $i=1,\dots,m$ and so $Rx_{i}\tilde{R}=y_{i}$ for all $j,i=1,\dots,m$
because of Lemma 2. ∎

The explicit construction in the proof of Theorem 3.6 gives us the following
algorithm.

Calculate the rotor $R=R_{1}\dots R_{m}$

$x_{i}\cdot x_{j}=y_{i}\cdot y_{j}$ and $x_{1}\wedge\cdots\wedge
x_{m}=y_{1}\wedge\cdots\wedge y_{m}$

$y_{i}=Rx_{i}\tilde{R}$

for $m>i>0$ do

$V_{i}\leftarrow x_{1}\wedge\cdots\wedge x_{i}$

$W_{i}\leftarrow y_{1}\wedge\cdots\wedge y_{i}$

end for

$R\leftarrow Id$

for $m>i>0$ do

$V_{i}\leftarrow RV_{i}\tilde{R}$

$H_{V}\leftarrow V_{i}\wedge W_{i+1}^{*}$

$H_{W}\leftarrow W_{i}\wedge W_{i+1}^{*}$

$R_{i}\leftarrow\widehat{1+\hat{H}_{V}^{*}\hat{H}_{W}^{*}}$

$R\leftarrow R_{i}R$

end for

##  4\. Nilpotent control problems in GA approach

We use the symmetries of $SO(3)$ to define an equivalence relation on the set
of geodesics passing through the origin, see Proposition 1. We find a
convenient representative of any equivalence class and describe the moduli
space in the language of GA.

###  4.1. Geodesics of $(3,6)$

Since the vector product $x\times h$ coincides with minus the dual of wedge
product $x\wedge h$ according to (25), the horizontal system (LABEL:rr1) can
be written in the form

(27) $\displaystyle\begin{split}\dot{x}&=h,\\\ \dot{z}&=-{1\over 2}x\wedge
h\end{split}$

where $x\in\wedge^{1}\mathbb{R}^{3}$ represents a line and
$z\in\wedge^{2}\mathbb{R}^{3}$ represents a plane in $\mathbb{R}^{3}$. In this
way we see the geodesics as curves in the geometric algebra $\mathbb{G}_{3}$.

######  Proposition 2.

Each arc–length parameterized sub–Riemannian geodesic satisfying the initial
condition $x_{i}(0)=0,z_{i}(0)=0$, $i=1,2,3$ is equivalent to a curve in
$M\cong\wedge^{1}\mathbb{R}^{3}\oplus\wedge^{2}\mathbb{R}^{3}\subset\mathbb{G}_{3}$
and, up to the action of a suitable $R\in SO(3),$ it takes form

(28)
$\displaystyle\begin{split}q(t)&=x(t)+z(t)=\frac{D}{K}(1-\cos(Kt))e_{1}+\frac{D}{K}\sin(Kt)e_{2}+C_{3}te_{3}\\\
&-\frac{D^{2}}{2K^{2}}(Kt-\sin(Kt))e_{1}\wedge e_{2}\\\
&-\frac{C_{3}D}{2K^{2}}(Kt-2\sin(Kt)+Kt\cos(Kt))e_{3}\wedge e_{1}\\\
&+\frac{C_{3}D}{2K^{2}}(2-Kt\sin(Kt)-2\cos(Kt))e_{2}\wedge e_{3},\end{split}$

where $K>0$ and $D,C_{3}$ satisfy the level set equation $D^{2}+C_{3}^{2}=1.$

###### Proof.

The solution to the vertical system (11) can be rewritten as

(29)
$\displaystyle\begin{split}h(t)=D\sin(Kt)\bar{v}_{1}+D\cos(Kt)\bar{v}_{2}+C_{3}v_{3}\end{split},$

where we denote $D=\sqrt{C_{1}^{2}+C_{2}^{2}},$ and the orthonormal vectors
$\bar{v}_{1},\bar{v}_{2}$ are obtained by the rotation of orthonormal vectors
$v_{1},v_{2}$ as

$\bar{v}_{1}=\frac{1}{\sqrt{C_{1}^{2}+C_{2}^{2}}}(-C_{1}v_{1}+C_{2}v_{2}),\;\bar{v}_{2}=\frac{1}{\sqrt{C_{1}^{2}+C_{2}^{2}}}(C_{2}v_{1}+C_{1}v_{2}).$

Thus the vectors $\bar{v}_{1},\bar{v}_{2},v_{3}$ are orthonormal with respect
to the Euclidean metric on $\mathbb{R}^{3}$. So, there is an orthogonal matrix
$R\in SO(3)$ that aligns vectors $\bar{v}_{1},\bar{v}_{2},v_{3}$ with the
standard basis of $\mathbb{R}^{3}$. Thus we get

$\bar{v}_{1}=Re_{1},\ \ \ \bar{v}_{2}=Re_{2},\ \ v_{3}=Re_{3},$

where $e_{1},e_{2}$ and $e_{3}$ are the elements of the standard Euclidean
basis of $\mathbb{R}^{3}$. According to (22), the rotor $R$ defines a
representative of a geodesic class $(R^{T}x(t),R^{T}z(t))$ which is a solution
to (27) for $h(t)=D\sin(Kt)e_{1}+D\cos(Kt)e_{2}+C_{3}e_{3}$. The solution (28)
then follows by a direct integration when the initial condition is applied.
Equation for the level set follows from the definition of $D$. ∎

The action of $SO(3)$ on $M\cong\mathbb{R}^{6}$ given by equation (22) defines
a moduli space $M/SO(3)$. We see $M$ as a subset of $\mathbb{G}_{3}$ and the
group $SO(3)$ is represented by rotors instead of matrices, which act on $M$
by conjugation. The action preserves the vector and bivector parts, inner
product, norm and dualization with respect to $*$. We can see the elements of
$M$ as the pairs consisting of lines and planes. The natural invariants are
the norms of lines’ directional vectors, norms of the planes’ normal vectors
and angles between these pairs of vectors. Square norm of the normal vector of
the plane $z^{*}\cdot z^{*}$ is $-z\cdot z$. Scalar product between the
directional vector of the line $x$ and the normal vector of the plane $z$ can
be rewritten as $(x\wedge z)^{*}$ because $(x\cdot z^{*})^{*}=x\wedge z$ and
$x\cdot z^{*}=(x\wedge z)^{*}$. Altogether, we consider three invariants

  * •

the square norm of the vector $x$, i.e. $x\cdot x$,

  * •

the square norm of the bivector $z$, i.e. $z\cdot z$,

  * •

the element $(x\wedge z)^{*}$ ,

where $\cdot$ coincides with the inner product on $\mathbb{G}_{3}$. In
particular, these invarians form a coordinate system on the moduli space
$M/SO(3)$.

######  Proposition 3.

Each geodesic starting at the origin defines a curve in the moduli space
$M/SO(3)$, which is determined by the invariants in the following way

(30) $\displaystyle\begin{split}x\cdot
x=&-\frac{2D^{2}}{K^{2}}(\cos(Kt)-1)+C_{3}^{2}t^{2}\\\ z\cdot
z=&-\frac{D^{2}}{4K^{4}}((4C_{3}^{2}K^{2}-4C_{3}^{2}-D^{2})\cos(Kt)^{2}\\\
&+2KC_{3}^{2}(2t(K-1)\sin(Kt)+t^{2}K-4)\cos(Kt)\\\
&-2Kt(4C_{3}^{2}+D^{2})\sin(Kt)+t^{2}(2C_{3}^{2}+D^{2})K^{2}+D^{2}+8C3^{2})\\\
(x\wedge
z)^{*}=&\>\frac{D^{2}{C_{3}}}{2K^{3}}((-2K+2)\cos(Kt)^{2}+(2K+2)\cos(Kt)+K^{2}t^{2}+Kt\sin(Kt)-4)\end{split}$

###### Proof.

Follows directly from (28). ∎

###  4.2. Geodesics of $(4,7)$

The base system (LABEL:xl) can be seen as a system in geometric algebra
$\mathbb{G}_{4}$

(31) $\displaystyle\begin{split}\dot{x}+\dot{\ell}&=h_{0}+\bar{h},\\\
\dot{y}&=-x\wedge\bar{h}-\ell\wedge h_{0},\end{split}$

where we assume that $x$ and $h_{0}$ are collinear with $e_{1}$ and
$\ell,\bar{h}$ in the subspace generated by $e_{2},e_{3},e_{4}.$ The form of
the second equation implies that $y$ is given by minus the wedge product of
$e_{1}$ and a vector from this subspace. Hence the solution $y(t)$ can be
viewed as a curve of planes in $\mathbb{G}_{4}.$

######  Proposition 4.

Each arc–length parameterized sub–Riemannian geodesic satisfying the initial
condition $x(0)=0,\ell_{i}(0)=0,y_{i}(0)=0,i=1,2,3$ is equivalent to a curve
in
$N\cong\wedge^{1}\mathbb{R}^{4}\oplus\wedge^{2}\mathbb{R}^{4}\subset\mathbb{G}_{4}$
and, up to the action of suitable $R\in SO(3),$ it takes form

(32)
$\displaystyle\begin{split}q(t)&=x(t)+\ell(t)+y(t)=(C_{1}\cos(Kt)+C_{2}\sin(Kt)-C_{1})e_{1}\\\
&+(C_{1}\sin(Kt)-C_{2}\cos(Kt)+C_{2})e_{2}+Cte_{3}\\\
&+\frac{1}{2}(C_{1}^{2}+C_{2}^{2})(tK-\sin(Kt))e_{1}\wedge e_{2}\\\
&+\frac{C}{2K}((2C_{1}-C_{2}Kt)\sin(Kt)-(C_{1}Kt+2C_{2})\cos(Kt)+2C_{2}-tC_{1}K)e_{1}\wedge
e_{3},\end{split}$

where $K>0$ and the constants $C_{1},C_{2},C$ satisfy the level condition
$K^{2}(C_{1}^{2}+C_{2}^{2})+C^{2}=1.$

###### Proof.

According to the vertical system (19), the vector $\bar{h}(t)$ lies in the
subspace generated by the vectors $r_{1},r_{2}$ for any $t$. Since the vectors
$r_{1}$ and $r_{2}$ are orthonormal, there is an orthogonal matrix $R\in
SO(3)$ that aligns these vectors with the second and third vector of the
standard basis of $\mathbb{R}^{3}$, i.e

$r_{1}=Re_{2},\ \ \ r_{2}=Re_{3}.$

Due to the symmetry of this system, see (23), this rotor defines a
representative of the geodesic class $(x(t),R^{T}\ell(t),R^{T}y(t))$ which is
the solution to the horizontal system (LABEL:xl) for

$\displaystyle h_{0}$ $\displaystyle=K(C_{2}\cos(Kt)-C_{1}\sin(Kt))$
$\displaystyle\bar{h}(t)$
$\displaystyle=K(C_{2}\sin(Kt)+C_{1}\cos(Kt))e_{1}+Ce_{2}$

or, equivalently, a curve in
$\mathbb{R}^{4}\oplus\Lambda^{2}\mathbb{R}^{4}\in\mathbb{G}_{4}$ given by the
solution of (31). By direct integration of this equation and by imposing the
initial conditions, we get the formula (32) for the solution. ∎

The action of $SO(3)$ on $N\cong\mathbb{R}^{7}$ given by (23), defines a
moduli space $N/SO(3)$. We see $N$ as a subset of $\mathbb{G}_{4}$ and the
group $SO(3)$ is represented by rotors instead of matrices, which act on $N$
by the conjugation. The action preserves the vector and bivector part, the
split $x+\ell$, inner product, norm and dualization with respect to $*$. The
orbits of this action are determined by natural invariants. For the same
reason as in the case of $(3,6)$ and due to the invariant split, we have three
invariants as follows

  * •

the value of the coordinate $x$,

  * •

the square of the norm of the vector $\ell$, i.e. $\ell\cdot\ell$,

  * •

the square of the norm of the bivector $y$, i.e. $y\cdot y$.

We need one more invariant for the dimensional reasons but the element
$(\ell\wedge y)^{*}$ is not scalar but vector. On the other hand, $(\ell\cdot
y)$ is a multiple of the vector $e_{1}$, so the value of $(\ell\cdot y)e_{1}$
is a scalar. As the last invariant we consider

  * •

the value of $(\ell\cdot y)e_{1}$.

These form the coordinate system on the moduli space $N/SO(3)$.

######  Proposition 5.

Each geodesic starting at the origin defines a curve in the moduli space
$N/SO(3),$ which is determined by the invariants in the following way

(33) $\displaystyle\begin{split}x&=C_{1}(\cos Kt-1)+C_{2}\sin Kt,\\\
\ell\cdot\ell&=(C_{1}\sin Kt+C_{2}(1-\cos Kt))^{2}+(Ct)^{2},\\\ (\ell\cdot
y)e_{1}&=\frac{1}{2}\Big{(}(C_{1}^{2}+C_{2}^{2})\big{(}C_{1}\sin
Kt+C_{2}(1-\cos Kt)\big{)}(Kt-\cos Kt)\\\ &+\frac{C^{2}}{K}t\big{(}C_{1}(2\sin
Kt-Kt\cos Kt-Kt)\\\ &+C_{2}(2-2\cos Kt-Kt\sin Kt)\big{)}\Big{)},\\\ y\cdot
y&=\frac{1}{4}\Big{(}(C_{1}^{2}+C_{2}^{2})^{2}(Kt-\cos
Kt)^{2}+\frac{C^{2}}{K^{2}}\big{(}C_{1}(2\sin Kt-Kt\cos Kt-Kt)\\\
&+C_{2}(2-2\cos Kt-Kt\sin Kt)\big{)}^{2}\Big{)}.\end{split}$

###### Proof.

Follows directly from (32). ∎

##  5\. Examples

In the sequel, we present two examples of controls based on the symmetries in
geometric algebra approach. We have the following scheme based on Algorithm
3.5.

  1. (1)

For the target point $q_{t}$ compute the invariants of the chosen particular
control system (2).

  2. (2)

Solve the system of non–linear equations (30) or (33) in the moduli space.

  3. (3)

Find the family of curves (28) or (32) going from the origin to the same point
$q_{o}$ that belongs to the same $SO(3)$ orbit of $q_{t}$.

  4. (4)

Find $R\in SO(3)$, such that $R(q_{o})=q_{t}$.

  5. (5)

Apply $R$ on the set of curves (28) or (32) to get a family of curves going
from the origin to the target point $q_{t}$.

The explicit calculations were acquired using a CAS system Maple similarly to
the paper [6].

###  5.1. Example in $(3,6)$

Our goal is to find the geodesic going from the origin to the target point

$q_{t}=(x_{t},z_{t})=2e_{1}-e_{2}+3e_{3}+e_{1}\wedge e_{2}-2e_{1}\wedge
e_{3}-2e_{2}\wedge e_{3}$

using the invariants (30) in the target point. We have

$x\cdot x=14,\>\>z\cdot z=-9,\>\>(x\wedge z)^{*}=3$

and together with the level set condition we get the system with the
invariants at $q_{t}$. We solve the system numerically in Maple and present
the solution with rounding up to four decimal digits

(34) $\displaystyle C_{3}$ $\displaystyle=0.7252,D=0.6885,K=0.9886$ (35)
$\displaystyle t$ $\displaystyle=5.0236$

Using the constants (34), we get the geodesic in the moduli space from the
origin to the point $q_{o}$ in the form

$\displaystyle q$
$\displaystyle=(x,z)=(0.6965(1-\cos(0.9886t)))e_{1}+0.6965\sin(0.9886t)e_{2}+0.7252te_{3}$
$\displaystyle-(0.2425(0.9886t-\sin(0.9886t)))e_{1}\wedge e_{2}$
$\displaystyle+(0.2555(0.9886t-2\sin(0.9886t)+0.9886t\cos(0.9886t)))e_{1}\wedge
e_{3}$
$\displaystyle+(0.2555(2-0.9886t\sin(0.9886t)-2\cos(0.9886t)))e_{2}\wedge
e_{3}$

and at the time $t=5.0236$ we reach the point

(36)
$\displaystyle\begin{split}q_{o}&=0.5216e_{1}-0.6741e_{2}+3.643e_{3}-1.439e_{1}\wedge
e_{2}+2.082e_{1}\wedge e_{3}+1.611e_{2}\wedge e_{3}\end{split}$

We are looking for the rotor which maps the multivector $q_{o}$ on the
multivector $q_{t}$. We consider the complete flags

$\displaystyle\\{0\\}\subset NO(x_{t})\subset NO(x_{t}\wedge z_{t}^{*})\subset
NO(z_{t}\wedge z_{t}^{*})\cong\mathbb{R}^{3},$ $\displaystyle\\{0\\}\subset
NO(x_{o})\subset NO(x_{o}\wedge z_{o}^{*})\subset NO(z_{o}\wedge
z_{o}^{*})\cong\mathbb{R}^{3}.$

We set $R_{m}=R_{3}=\text{id}$ and map the plane $x_{o}\wedge z_{o}^{*}$ to
the plane $x_{t}\wedge z_{t}^{*}$ by the rotor $R_{m-1}=R_{2}$ according to
the formula (26). Explicitly,

$R_{2}:=0.1334+0.7083e_{1}\wedge e_{2}-0.5483e_{1}\wedge
e_{3}-0.4242e_{2}\wedge e_{3}$

and we can map the multivector $q_{o}$ on the multivector
$q_{s}=(x_{s},z_{s})=R_{2}q_{o}\tilde{R_{2}}$ in such a way that $x_{s}$ and
$z_{s}$ lie in the plane $x_{o}\wedge z_{o}^{*}$. Explicitly,

$q_{s}=(x_{s},z_{s})=-2.8510e_{1}+2.3208e_{2}-0.6956e_{3}-2.641e_{1}\wedge
e_{2}+1.2523e_{1}\wedge e_{3}+0.6767e_{2}\wedge e_{3}.$

Finally, we map the plane $x_{s}\wedge(x_{s}\wedge z_{s})^{*}$ to the plane
$x_{t}\wedge(x_{t}\wedge z_{t})^{*}$ by rotor

$R_{m-2}=R_{1}=0.3727+0.1716e_{1}\wedge e_{2}+0.6863e_{1}\wedge
e_{3}-0.6005e_{2}\wedge e_{3}.$

Altogether, we found the rotor $R=R_{1}R_{2}R_{3}$ and, when applied on (28),
we got a geodesic going from the origin to the point $q_{t}$ in the form

$\displaystyle q$
$\displaystyle=(x,z)=(-0.5302+0.4673t+0.5302\cos(0.9886t)-0.05136\sin(0.9886t))e_{1}$
$\displaystyle-(0.008(22.124+36.347t-22.124\cos(0.9886t)+76.628\sin(0.9886t)))e_{2}$
$\displaystyle+(-0.4156\cos(0.9886t)+0.4723t+0.4156-0.3267\sin(0.9886t))e_{3}$
$\displaystyle-(0.2(-1.5244+0.7536t\sin(0.9886t)+0.4086\sin(0.9886t)+1.5244\cos(0.9886t)$
$\displaystyle-0.5923t\cos(0.9886t)+0.1884t)){e_{1}\wedge e_{2}}$
$\displaystyle+(-0.3184t+0.1299-0.2223t\cos(0.9886t)-0.1299\cos(0.9886t)$
$\displaystyle-0.06419t\sin(0.9886t)+0.547\sin(0.9886t)){e_{1}\wedge e_{3}}$
$\displaystyle+(-0.389+0.1923t\sin(0.9886t)-0.1358t+0.0186t\cos(0.9886t)$
$\displaystyle+0.389\cos(0.9886t)+0.1186\sin(0.9886t))e_{2}\wedge e_{3}.$

In Figure 2 we present the trajectories $(x_{1},x_{2},x_{3})$ and
$(z_{1},z_{2},z_{3})$, respectively.

Figure 2. Trajectories $(x_{1},x_{2},x_{3})$ and $(z_{1},z_{2},z_{3})$

###  5.2. Example in $(4,7)$

Our goal is to find the geodesic going from the origin to the target point

$q_{t}=(x_{t},\ell_{t},y_{t})=e_{1}+2e_{2}+e_{3}+3e_{4}-e_{1}\wedge
e_{2}+2e_{1}\wedge e_{3}+2e_{1}\wedge e_{4}$

using the invariants (33) at the target point. We compute

$x=1,\>\>\ell\cdot\ell=14,\>\>y\cdot y=-9,\>\>(\ell\cdot y)e_{1}=-6,$

and together with the level set condition, we get the system with the
invariants at $q_{t}$. We solve the system numerically in Maple and we present
the solution with the constants rounded up to four decimal digits as follows

(37) $\displaystyle C$
$\displaystyle=0.6126,C_{1}=-0.7816,C_{2}=-0.5324,K=0.8358,$ (38)
$\displaystyle t$ $\displaystyle=6.0748.$

Using the constants (37) we get a geodesic in the moduli space from the origin
to the point $q_{o}$ in the form

$\displaystyle\begin{split}q=(x,\ell,y)&=(-0.7816\cos(0.8358t)-0.5324\sin(0.8358t)+0.7816)e_{1}\\\
&+(-0.7816\sin(0.8358t)+0.5324\cos(0.8358t)-0.5324)e_{2}\\\
&+0.6126te_{3}+0.4471(0.8358t-\sin(0.8358t))e_{1}\wedge e_{2}\\\
&+0.3665\big{(}(0.4449t-1.563)\sin(0.8358t)+(0.6533t+1.065)\cos(0.8358t)\\\
&-1.065+0.6533t\big{)}e_{1}\wedge e_{3}\end{split}$

and at the time $t=6.0748$, we reach the point

(39)
$\displaystyle\begin{split}q_{o}&=(x_{o},\ell_{o},y_{o})=e_{1}+0.3878e_{2}+3.722e_{3}+2.688e_{1}\wedge
e_{2}+1.332e_{1}\wedge e_{3}.\end{split}$

We are looking for the rotor which maps the multivector $q_{o}$ on the
multivector $q_{t}$. We shall consider the complete flags starting with the
line $NO(\ell)$ and ending with the space $NO(\ell\wedge y)$. To find the
middle one, we can use the projection of the line $NO(\ell)$ onto the plane
$NO(\ell\wedge y)$. Thus we get

$\displaystyle\\{0\\}\subset NO(\ell_{t})\subset
NO(\ell_{t}\wedge(\ell_{t}\wedge
y_{t}^{*})^{*})=NO(\ell_{t}\wedge(\ell_{t}\cdot y_{t}))$ $\displaystyle\subset
NO(\ell_{t}\wedge y_{t})\subset NO(y_{t}\wedge y_{t}^{*})\cong\mathbb{R}^{4},$
$\displaystyle\\{0\\}\subset NO(\ell_{o})\subset
NO(\ell_{o}\wedge(\ell_{o}\wedge
y_{o}^{*})^{*})=NO(\ell_{o}\wedge(\ell_{o}\cdot y_{o}))$ $\displaystyle\subset
NO(\ell_{o}\wedge y_{o})\subset NO(y_{o}\wedge y_{o}^{*})\cong\mathbb{R}^{4}.$

First, we map the hyperplane $\ell_{o}\wedge y_{o}$ to the hyperplane
$\ell_{t}\wedge y_{t}$ by the rotor $R_{3}$ according to the formula (26). We
obtain

$\displaystyle R_{3}$ $\displaystyle:=0.4863-0.4335e_{2}\wedge
e_{4}-0.7587e_{3}\wedge e_{4}.$

The next step is to map the hyperplane $NO(\ell_{o}\wedge(\ell_{o}\cdot
y_{o})\wedge(\ell_{o}\wedge y_{o})^{*})$ to the hyperplane
$NO(\ell_{t}\wedge(\ell_{t}\cdot y_{t})\wedge(\ell_{t}\wedge y_{t})^{*})$ by
the rotor $R_{2}$ according to the formula (26). Explicitly,

(40) $\displaystyle R_{2}$ $\displaystyle:=0.0387+0.9993e_{2}\wedge e_{3}.$

Finally, we map the hyperplane
$NO(\ell_{o}\wedge(\ell_{o}\wedge(\ell_{o}\wedge y_{o}^{*})^{*})^{*})$ to the
hyperplane $NO(\ell_{t}\wedge(\ell_{t}\wedge(\ell_{t}\wedge
y_{t}^{*})^{*})^{*})$ by the rotor $R_{1}$ according to the formula (26). It
turns out that $R_{1}=1.$ Altogether, $R=R_{1}R_{2}R_{3}$ and, when applied on
(5.2), we get a geodesic going from the origin to the point $q_{t}$ as

$\displaystyle q=(x,\ell,y)$
$\displaystyle=(-0.7816\cos(0.8358t)-0.5324\sin(0.8358t)+0.7816)e_{1}$
$\displaystyle+(0.5261\sin(0.8358t)-0.3584\cos(0.8358t)+0.3946t+0.3584)e_{2}$
$\displaystyle+(-0.4749\sin(0.8358t)+0.3235\cos(0.8358t)+0.1235t-0.3235)e_{3}$
$\displaystyle+((0.2513+0.1542t)\cos(0.8358t)+(-0.06803+0.1050t)\sin(0.8358t)$
$\displaystyle-0.2514-0.09731t)e_{1}\wedge e_{2}$
$\displaystyle+((0.07868+0.04827t)\cos(0.8358t)+(-0.3871+0.03287t)\sin(0.8358t)$
$\displaystyle-0.07869+0.2753t)e_{1}\wedge e_{3}$
$\displaystyle+((0.2880+0.1767t)\cos(0.8358t)+(0.1203t-0.6112)\sin(0.8358t)$
$\displaystyle+0.3342t-0.2880)e_{1}\wedge e_{4}.$

In Figure 3 we present trajectories $x$, $(\ell_{2},\ell_{3},\ell_{4})$ and
$(y_{1},y_{2},y_{3})$, respectively.

Figure 3. Trajectories $x$, $(\ell_{2},\ell_{3},\ell_{4})$ and
$(y_{1},y_{2},y_{3})$

##  6\. Conclusion

We presented the use of geometric algebra for the control systems invariant
with respect to the orthogonal transformations. The main contribution of GA
lies in a construction of the rotor between two bases of a vector space based
only on algebraic computations in a chosen GA. This allows us to use the
geometric objects effectively and, analogously to quaternions, the
implementations are faster than the usual computations with matrices. We
assessed an algorithm and illustrated its use on two particular examples with
filtration $(3,6)$ corresponding to a trident snake robot control, and $(4,7)$
corresponding to the control of a trident snake with flexible leg. All
calculations were acquired using Maple packages Clifford [19] and
DifferentialGeometry [3].

###  6.1. Acknowledgements

This work does not have any conflicts of interest. The first three authors
were supported by the grant no. FSI-S-20-6187. Fourth author is supported by
the grant no. 20-11473S Symmetry and invariance in analysis, geometric
modeling and control theory from the Czech Science Foundation. Finally, we
thank the referee for valuable comments.

## References

  * [1] A. Agrachev, D. Barilari, U. Boscain,  A comprehensive introduction to sub-Riemannian geometry. From the Hamiltonian viewpoint, Cambridge Studies in Advanced Mathematics, Vol. 181, Cambridge University Press (2020)
  * [2] A. A. Agrachev and Y. L. Sachkov, Control theory from the geometric viewpoint, Encyclopaedia of Mathematical Sciences, Control Theory and Optimization, Springer (2006)
  * [3] I. M. Anderson and C. G. Torre, New symbolic tools for differential geometry, gravitation, and field theory,  J. Math. Phys. 53 (2015) 1–12 http://dx.doi.org/10.1063/1.3676296
  * [4] A. Bellaiche,  The tangent space in sub-Riemannian geometry, Sub-Riemannian Geometry (1996) 1–78
  * [5] A.M. Bloch, Nonholonomic Mechanics and Control, Springer (2015)
  * [6] J. Hrdina, A. Návrat, L. Zalabová, Symmetries in geometric control theory using Maple. Math. Comput. Simul. 190 (2021) 474–493
  * [7] J. Hrdina, A. Návrat, L. Zalabová, On symmetries of Sub–Riemannian structure with growth vector (4,7),Annali di Matematica (2022). https://doi.org/10.1007/s10231-022-01242-6
  * [8] J. Hrdina, A. Návrat, P. Vašík, R. Matoušek, Geometric control of the trident snake robot based on CGA. Adv. Appl. Clifford Algebras (2017) 27: 633. https://doi.org/10.1007/s00006-016-0693-7.
  * [9] J. Hrdina, P. Vašík, Notes on differential kinematics in conformal geometric algebra approach. In Advances in Intelligent Systems and Computing; Springer: Cham, Switzerland, Volume 378 (2015) 363–374. https://doi.org/10.1007/978-3-319-19824-8_30
  * [10] J. Hrdina, L. Zalabová, Local geometric control of a certain mechanism with the growth vector (4,7), J. Dyn. Control Syst. 26 (2020) 199–216 https://doi.org/10.1007/s10883-019-09460-7
  * [11] M. Ishikawa,  Trident snake robot: locomotion analysis and control.  IFAC Symposium on Nonlinear Control Systems 6 (2004) 895–900 https://doi.org/10.1016/S1474-6670(17)31339-3
  * [12] F. Jean,  Control of Nonholonomic Systems: From Sub–Riemannian Geometry to Motion Planning.  Springer (2014)
  * [13] V. Jurdjevic, Geometric Control Theory (Cambridge Studies in Advanced Mathematics. Cambridge: Cambridge University Press (1996)
  * [14] A. Montanari, G. Morbidelli, On the sub–Riemannian cut locus in a model of free two-step Carnot group. Calc. Var. 56(36) (2017) https://doi.org/10.1007/s00526-017-1149-1
  * [15] O. Myasnichenko,  Nilpotent (3, 6) sub–Riemannian problem.  Journal of Dynamical and Control Systems, 8(4) (2002) 573–597https://doi.org/10.1023/A:1020719503741
  * [16] L. Rizzi, U. Serres, On the cut locus of free, step two Carnot groups, Proc. Amer. Math. Soc. 145 (2017) 5341–5357 https://doi.org/10.1090/proc/13658
  * [17] D. Hildenbrand, Foundations of Geometric Algebra Computing, Springer Verlag (2013)
  * [18] Ch. Perwass, Geometric Algebra with Applications in Engineering, Springer Verlag (2009)
  * [19] R. Ablamowicz, B. Fauser, Mathematics of Clifford - a Maple package for Clifford and Graßmann algebras, Adv. Appl. Clifford Algebr. 15, (2005) 157–181. https://doi.org/10.1007/s00006-005-0009-9
  * [20] P. Lounesto, Clifford Algebra and Spinors. 2nd edn. CUP, Cambridge (2006).


# Look-ups are not (yet) all you need
for deep learning inference

Calvin McCarter
<EMAIL_ADDRESS>
&Nicholas Dronen*
<EMAIL_ADDRESS>
Work done prior to employment at Amazon

###### Abstract

Fast approximations to matrix multiplication have the potential to
dramatically reduce the cost of neural network inference. Recent work on
approximate matrix multiplication proposed to replace costly multiplications
with table-lookups by fitting a fast hash function from training data. In this
work, we propose improvements to this previous work, targeted to the deep
learning inference setting, where one has access to both training data and
fixed (already learned) model weight matrices. We further propose a fine-
tuning procedure for accelerating entire neural networks while minimizing loss
in accuracy. Finally, we analyze the proposed method on a simple image
classification task. While we show improvements to prior work, overall
classification accuracy remains substantially diminished compared to exact
matrix multiplication. Our work, despite this negative result, points the way
towards future efforts to accelerate inner products with fast nonlinear
hashing methods.

##  1 Introduction

To reduce the computational cost of dense matrix multiplications (or, more
generally, inner products) in neural network (NN) inference, recent research
efforts have sought to develop cheap approximations with minimal sacrifice of
classification accuracy. Methods based on distillation (Bucilua et al.,, 2006;
Hinton et al.,, 2015) and pruning (Han et al.,, 2015; Janowsky,, 1989) have
proposed to reduce the number of stored parameters and the number of multiply-
accumulate (MAC) operations, while methods based on scalar quantization (Wu et
al.,, 2016) have proposed to reduce the size of stored parameters and the cost
of each MAC. A recent method, Maddness (Blalock & Guttag,, 2021), instead
proposed to approximate MACs by replacing them entirely with lookup-accumulate
(LAC) operations. This reduced computational cost by involving fewer LACs (and
less memory) than required by MACs, and by utilizing hardware multiplexers,
which require fewer transistors than hardware multipliers.

Maddness breaks down a matrix multiply $\bm{A}\bm{B}$ (where $\bm{A}$ contains
inputs / activations and $\bm{B}$ contains weights) into an assemblage of
$\bm{a}^{\top}\bm{b}$ inner products, and approximates each of them by
building on (Gray,, 1984) and product quantization (PQ) (Gray,, 1984; Jegou et
al.,, 2010). PQ enables fast inner products of the form $\bm{a}^{\top}\bm{b}$
when a training set of sampled $\bm{a}$s is given in advance. For paired
vectors $\bm{a}$ and $\bm{b}$ of dimension $D$, the PQ algorithm partitions
their dimensions to form $C\leq D$ pairs of subvectors (and subspaces). The
inner product is then computed as the sum of approximate partial inner
products over subspaces. To do this, the $K$-Means algorithm is run over each
subspace of the training set of $\bm{a}s$ to obtain $C$ subspace-specific sets
of prototypes. Since all $\bm{b}$s (corresponding to columns of NN weight
matrices) are also known in advance, all possible inner products between
subspace prototypes and $\bm{b}$ subvectors are precomputed offline and stored
in $C$ lookup tables. (Each subspace-specific lookup table is referred to as a
codebook, so $C$ is the number of codebooks.) To compute (approximate) inner
products for a newly-seen $\bm{a}$ at inference time, PQ identifies, for each
subspace, the nearest subspace prototype (represented by its index) for each
subvector of $\bm{a}$. This mapping from subvectors to categorical variables,
represented as $\log_{2}(K)$-bit integer indices where $K$ is the number of
prototypes per subspace, is called the encoding function. Then, the computed
inner product is simply the sum of partial inner products obtained from
indexing into the appropriate rows in the precomputed lookup tables.

The primary contribution of Maddness addressed the computational bottleneck of
PQ, with a fast new encoding function. Maddness proposed a new hashing-based
encoder, based on balanced binary regression trees, that runs efficiently on
modern processors when $K=16$; this encoder was learned from the training set
of $\bm{a}s$ to minimize reconstruction error. Also, rather than choosing
prototypes via $K$-Means, Maddness optimized a least-squares objective, which
also minimized reconstruction error; this objective conditioned on the
training set of $\bm{a}$s, as well as the (already learned) encoder’s
resulting indices from the training samples. Finally, lookup tables were
quantized to 8-bit block floating point format, with fast rounded 8-bit
summation over subspaces.

While proposed to accelerate entire NNs, Maddness was validated for a large
variety of models and settings but only on the final classification layers. In
each of these settings, the authors varied the number of codebooks, the key
tuning parameter which trades off accuracy versus efficiency. For a NN
feedforward layer with $M$ output dimensions and Float32 weights, the weight
matrix requires $4DM$ bytes of storage. With Maddness, the lookup table
requires $KCM=16CM$ bytes of storage.

In this work, we propose a new method which builds on Maddness, and analyze
its applicability to whole NN inference. We dub our proposed method, Inference
Targeted LookUp-based Matrix Multiplication (ITLUMM) 111Spelled backward, this
is MMULTI.. First, we intelligently partition the inner dimension rather than
naively partitioning it into contiguous subspaces. Second, we directly
optimize the lookup table rather than prototypes, taking into account the
subsequent nonlinear activation function (or classification loss function).
Furthermore, for accelerating entire neural networks, we propose incremental
fine-tuning of linear layers.

We apply our approach to feedforward-only image classification on MNIST (LeCun
et al.,, 2010) and CIFAR-10 (Krizhevsky et al.,, 2009). We found that, while
Itlumm improves upon Maddness, the accuracy-efficiency tradeoff is not yet
good enough for practical use. Our analysis of results suggests that an
improved encoding function should be a key focus of future work.

##  2 Method

Here we describe Itlumm, which comprises intelligent subspace partitioning
(Section 2.1), model-aware optimization of the lookup tables (Section 2.2),
and fine-tuning of full networks (Section 2.3).

###  2.1 Intelligent Subspace Partitioning

For each subvector, Maddness performs exactly 4 comparisons in order to map
the subvector into one of 16 lookup-table rows. We would therefore prefer a
subvector with maximal mutual information among its elements, such that this
small number of comparisons provides maximal information about the entire
subvector. Thus, we should partition the vector $\bm{a}$ such that mutually-
informative dimensions are placed together in the same subvector. To do this,
we find a permutation of the $D$ dimensions of $\bm{a}$, then slice the
permuted dimensions into $C$ contiguous subspaces.

This problem is closely related to Optimized Product Quantization (OPQ) (Ge et
al.,, 2013) which rotates vectors before partitioning them, but in our case
our rotation matrix must be a permutation matrix. Directly optimizing for the
optimal permutation is challenging, so we try two different approximate
solutions for this problem.

####  2.1.1 OPQ-based Partitioning

We first run OPQ to obtain a fully-dense rotation matrix $\bm{R}$. Then, we
find a permutation matrix $\bm{Q}$ that approximates the effect of $\bm{R}$,
by solving the following maximum weight matching problem:

$\displaystyle\max_{\bm{Q}}\sum_{i}\sum_{j}\bm{R}_{ij}\bm{Q}_{ij}.$ (1)

We solve this using the Hungarian algorithm (Kuhn,, 1955). The resulting
bipartite graph matches dimensions of the original vector to their “closest
matching” dimensions in the OPQ-rotated vector. Once we have a permuted
ordering of the dimensions, we split the reordered dimensions into $C$
contiguous equally-sized chunks.

####  2.1.2 Squared-correlation Hierarchical Clustering

Here, we instead try to directly identify a permutation of the dimensions that
places correlated dimensions closely together. We first compute the squared
correlation matrix $R^{2}$, then perform hierarchical/agglomerative clustering
on it to produce a dendrogram of the dimensions Müllner, (2013). Next, we
convert the dendrogram into a permutation by finding the ordering of the leaf
nodes such that the distance between successive leaves is minimal. Finally, we
split the reordered dimensions into $C$ contiguous equally-sized chunks.

###  2.2 Model-Aware Lookup-Table Optimization

The prototype-optimization objective function of Maddness sought to maximize
the ability to reconstruct $\bm{A}$, given the prototypes and the encoder’s
output (i.e. indices into the lookup-tables). In other words, given a row
vector $\bm{a}^{\top}$ and the corresponding encoder output $\bm{g}^{\top}$
(the concatenation of $C$ one-hot 16-dimensional vectors), the prototype
matrix $\bm{P}$ was optimized so that
$\bm{a}^{\top}\approx\bm{g}^{\top}\bm{P}$. To do this, given a training matrix
$\bm{A}$, its corresponding encoder output matrix $\bm{G}$, and the original
$K$-means prototype matrix $\bm{P}_{0}$, Maddness optimized the following:

$\displaystyle\arg\min_{\bm{P}}\|\bm{A}-\bm{GP}\|^{2}+\lambda\|\bm{P}-\bm{P}_{0}\|^{2},$
(2)

where $\lambda$ is a regularizer. Finally, Maddness computed the lookup-table
via $\bm{T}=\bm{P}\bm{B}$.

We note that this procedure ignores the fact that for deep learning inference,
we also know the already-learned weight matrix $\bm{B}$. We also observe that
it aims at reconstructing the layer inputs $\bm{A}$, rather than the matrix
product $\bm{AB}$. In fact, in deep learning we actually care about not the
matrix product $\bm{AB}$, but rather the layer output $\sigma(\bm{AB})$. The
function $\sigma(\cdot)$ is the composition of the bias term with a
nonlinearity (either an elementwise activation function such as ReLU or a row-
wise nonlinearity for classification such as softmax). Optimizing quantization
for downstream performance has shown great success in previous works Guo et
al., (2020); May et al., (2019).

Therefore, in Itlumm we take advantage of our knowledge of both the model
weights and the subsequent nonlinearity. Furthermore, for computationally
efficiency at training time, we directly optimize the lookup table, rather
than the prototypes. This leads to the following objective:

$\displaystyle\arg\min_{\bm{T}}\|\sigma(\bm{AB})-\sigma(\bm{GT})\|^{2}+\lambda\|\bm{T}-\bm{P}_{0}\bm{B}\|^{2}.$
(3)

###  2.3 Fine-tuning for Acceleration of Full Neural Networks

Replacing a neural network layer with approximate hashing-based lookups may
degrade accuracy in two ways. Layer replacement may alter the distribution of
the layer’s outputs without destroying information; inference accuracy would
still be degraded since subsequent layers would undergo domain shift. In this
case, it would in principle be possible to retrain subsequent layers to adjust
to this shift and regain lost accuracy, assuming subsequent layers are
sufficiently high-capacity. Alternatively, layer replacement may actually
destroy information contained in the layer’s inputs, and subsequent layers
would be unable to recover original accuracy.

To accelerate full networks, while ameliorating the first cause of degraded
accuracy, we propose incremental layer replacement. For an $L$-layer NN with
layers indexed as $l=1,\dots,L$, we start with the input layer $1$ and proceed
towards the classification layer $L$. At each step $l$, we first replace the
layer by collecting its inputs for the training data. We then freeze layers
$1,\dots,l$ (now using our non-differentiable lookups) and fine-tune the
weights of layers $l+1,\dots,L$ (still using exact matrix multiplication) on
the training data.

##  3 Experiments

###  3.1 Approximating Softmax Classifiers

We first repeat the softmax classifier experiments of Blalock & Guttag, (2021)
using the final dense classifier layers of VGG-like neural networks trained on
CIFAR-100. In both these networks, the final activations are 512-dimensional,
so the $\bm{A}$ matrices are $10{,}000\times 512$ in both cases, while the
$\bm{B}$ matrices are $512\times 10$ and $512\times 100$, respectively.
Because MADDNESS was Pareto optimal compared to its prior works, we compare
only it to our improved methods. Also, because the inference runtimes and
memory usage are identical between Maddness and Itlumm, depending only on the
number of codebooks, we show results in terms of number of codebooks rather
than the runtime speedup (which relies on assumptions regarding hardware
support and software implementation). For each experimental setting, we plot
the relative accuracy, the classification accuracy of a given setting divided
by the accuracy obtained from using exact matrix multiplication.

Our results are shown in Figure 1. Overall, we see that all settings of Itlumm
(depicted in orange, green, and red) are superior to Maddness (depicted in
blue). We also see the effect of the partitioning heuristic in Figure 1(A).
For Itlumm we fix the prototype-optimization objective to be the KL-
divergence; we also plot the accuracy of Maddness as a reference. We see that
the OPQ heuristic provides essentially no improvement over naive (PQ)
partitioning. Meanwhile, $R^{2}$-clustering partitioning provides a modest
improvement. Meanwhile, we see the effect of the prototype-optimization
objective in Figure 1(B & C). We see that, regardless of the choice of
partitioning heuristic, minimizing the KL-divergence (KLD) is better than
minimizing mean-squared error (MSE).

(A) | (B) | (C)
---|---|---
|  |
Figure 1: Comparison of Maddness and Itlumm accuracy vs efficiency tradeoff
for CIFAR-100 classifier layer. (A) Effect of Itlumm partitioning strategy.(B)
Effect of nonlinearity function when Itlumm uses $R^{2}$-based partitioning.
(C) Effect of nonlinearity function when Itlumm uses OPQ-based partitioning.

###  3.2 Accelerating Full Networks

Our goal is to replace all of a network’s matrix multiplications with lookup
tables. To that end, we first evaluated our proposed approach on a simple
multi-layer perceptron (MLP) trained on the MNIST dataset. The network has 4
layers, with 30 neurons per layer. We began with an ablation study to measure
the effect that substituting lookup for multiplication has on accuracy in each
layer, then evaluated the effect of said substitution on every layer.

Figure 2: Effect of Itlumm on MNIST accuracy when converting a single layer to
use Itlumm instead of matrix multiplication. The original network is a 4-layer
MLP. The accuracy of the original network, with matrix multiplication in every
layer, is the upper bound (black dotted line). Any configuration with fewer
codebooks than the performance breakeven point is faster than matrix
multiplication (blue dashed line).

Figure 2 shows the accuracy of a network as a function of the number of
codebooks when a given layer’s matrix multiplication has been replaced by a
lookup table. Higher layers in the network, those closer to the output, are
more robust to the error introduced by a lookup tables. However, the
performance breakeven point is between 2 and 3 codebooks, and the degradation
of accuracy is quite severe. Because of the compounding nature of error in
neural networks, it’s likely that replacing all of the layers in the network
with Itlumm will result in even more severe degradation.

Number of codebooks | Accuracy | Faster
---|---|---
1 | 36.1 | Yes
2 | 36.3 | Yes
4 | 52.0 | No
8 | 70.5 | No
16 | 84.9 | No
Table 1: Effect of Itlumm on MNIST accuracy when converting every layer to use
Itlumm instead of matrix multiplication. A configuration in the table is
faster than matrix multiplication when its number of codebooks is less than
the breakeven point in Figure 2.

Indeed, Table 1 shows that replacing all of the matrix multiplications in this
network with Itlumm results in accuracy worse than that due to replacing any
single layer. Since the performance breakeven point is 2 codebooks,
accelerating the entire network Itlumm requires a loss of more than 60%
accuracy.

##  4 Discussion

As expected, ITLUMM’s model-aware lookup-table optimization improves upon
Maddness. Meanwhile, optimized partitioning provides little improvement,
surprisingly so given the substantial benefits of OPQ in maximum inner-product
search (MIPS) (Norouzi & Fleet,, 2013; Ge et al.,, 2013). However, we are not
able to obtain the full benefits of OPQ because the applied rotation must be a
permutation. Furthermore, in the MIPS setting, input data are expected to have
redundancy/correlation and a wide range of variances, while NN activations are
regularized (e.g. with dropout) Srivastava et al., (2014) to avoid this.

We suggest that future research focus on improving the Maddness hash function,
which appears to be an accuracy bottleneck. One could exploit knowledge of
weight matrix $\bm{B}$ while learning hash function parameters, as done for NN
scalar quantization Stock et al., (2020). More substantially, a differentiable
hash function could improve inference accuracy and speed up conversion of full
networks.

##  5 Conclusion

We proposed Itlumm, an approach for accelerating matrix multiplication by
replacing multiplications with look-ups, which improves upon previous work. In
the context of deep learning inference, where model weights are known, our
approach advanced classification accuracy beyond prior work in Maddness
(Blalock & Guttag,, 2021). Our approach and analysis on full neural networks
informs the community on the current state of hashing-based acceleration, and
points the way to potential future advances.

## References

  * Blalock & Guttag, (2021) Blalock, Davis, & Guttag, John. 2021.  Multiplying Matrices Without Multiplying.  Pages 992–1004 of: Meila, Marina, & Zhang, Tong (eds), Proceedings of the 38th International Conference on Machine Learning.  Proceedings of Machine Learning Research, vol. 139.  PMLR.
  * Bucilua et al., (2006) Bucilua, Cristian, Caruana, Rich, & Niculescu-Mizil, Alexandru. 2006.  Model compression.  Pages 535–541 of: Proceedings of the 12th ACM SIGKDD international conference on Knowledge discovery and data mining.
  * Ge et al., (2013) Ge, Tiezheng, He, Kaiming, Ke, Qifa, & Sun, Jian. 2013.  Optimized product quantization.  IEEE transactions on pattern analysis and machine intelligence, 36(4), 744–755.
  * Gray, (1984) Gray, Robert. 1984.  Vector quantization.  IEEE Assp Magazine, 1(2), 4–29.
  * Guo et al., (2020) Guo, Ruiqi, Sun, Philip, Lindgren, Erik, Geng, Quan, Simcha, David, Chern, Felix, & Kumar, Sanjiv. 2020.  Accelerating large-scale inference with anisotropic vector quantization.  Pages 3887–3896 of: International Conference on Machine Learning.  PMLR.
  * Han et al., (2015) Han, Song, Mao, Huizi, & Dally, William J. 2015.  Deep compression: Compressing deep neural networks with pruning, trained quantization and huffman coding.  arXiv preprint arXiv:1510.00149.
  * Hinton et al., (2015) Hinton, Geoffrey, Vinyals, Oriol, Dean, Jeff, et al. 2015\.  Distilling the knowledge in a neural network.  arXiv preprint arXiv:1503.02531, 2(7).
  * Janowsky, (1989) Janowsky, Steven A. 1989.  Pruning versus clipping in neural networks.  Physical Review A, 39(12), 6600.
  * Jegou et al., (2010) Jegou, Herve, Douze, Matthijs, & Schmid, Cordelia. 2010.  Product quantization for nearest neighbor search.  IEEE transactions on pattern analysis and machine intelligence, 33(1), 117–128.
  * Krizhevsky et al., (2009) Krizhevsky, Alex, Hinton, Geoffrey, et al. 2009\.  Learning multiple layers of features from tiny images.
  * Kuhn, (1955) Kuhn, Harold W. 1955.  The Hungarian method for the assignment problem.  Naval research logistics quarterly, 2(1-2), 83–97.
  * LeCun et al., (2010) LeCun, Yann, Cortes, Corinna, & Burges, Christopher J. 2010.  MNIST handwritten digit database. 2010.  URL http://yann. lecun. com/exdb/mnist, 7(23), 6.
  * May et al., (2019) May, Avner, Zhang, Jian, Dao, Tri, & Ré, Christopher. 2019.  On the downstream performance of compressed word embeddings.  Advances in neural information processing systems, 32.
  * Müllner, (2013) Müllner, Daniel. 2013.  fastcluster: Fast hierarchical, agglomerative clustering routines for R and Python.  Journal of Statistical Software, 53, 1–18.
  * Norouzi & Fleet, (2013) Norouzi, Mohammad, & Fleet, David J. 2013.  Cartesian k-means.  Pages 3017–3024 of: Proceedings of the IEEE Conference on computer Vision and Pattern Recognition.
  * Srivastava et al., (2014) Srivastava, Nitish, Hinton, Geoffrey, Krizhevsky, Alex, Sutskever, Ilya, & Salakhutdinov, Ruslan. 2014.  Dropout: a simple way to prevent neural networks from overfitting.  The journal of machine learning research, 15(1), 1929–1958.
  * Stock et al., (2020) Stock, Pierre, Joulin, Armand, Gribonval, Rémi, Graham, Benjamin, & Jégou, Hervé. 2020.  And the Bit Goes Down: Revisiting the Quantization of Neural Networks.  Pages 1–11 of: ICLR 2020-Eighth International Conference on Learning Representations.
  * Wu et al., (2016) Wu, Jiaxiang, Leng, Cong, Wang, Yuhang, Hu, Qinghao, & Cheng, Jian. 2016.  Quantized convolutional neural networks for mobile devices.  Pages 4820–4828 of: Proceedings of the IEEE conference on computer vision and pattern recognition.