§ INTRODUCTION

This technical note aims to introduce geometric controllers to roboticist for aligning 3D robots with non-constant 3D vector fields. This alignment entails the control of the robot's 3D attitude. We derive with excessive detail all the calculations needed for the analysis and implementation of the controllers.

The sources for the technical results and formal proofs employed throughout this note can be found in <cit.> for a concise introduction to $\text{SO}(3)$, [2] for Lie groups in robotics, and [3] for a more profound understanding of geometric controllers.

§ PRELIMINARIES, ROBOT KINEMATICS AND PROBLEM FORMULATION

We represent a robot in the 3D space with a Body-fixed Cartesian frame of coordinates $\mathcal{F}_B$. The position and attitude of the robot with respect to the Earth-fixed frame $\mathcal{F}_E$ are respectively described by the vector $p \in \mathbb{R}^3$ and the matrix

\begin{equation} \label{eq: R}
R = \begin{bmatrix}x^E & y^E & z^E\end{bmatrix},
\end{equation}

where $x^E, y^E, z^E \in S^2$ are the orthonormal column vectors represented in $\mathcal{F}_E$ which generate the vector basis of $\mathcal{F}_B$ (see <ref>). Indeed, the three more familiar Euler angles that describe the orientation of $\mathcal{F}_B$ with respect to $\mathcal{F}_E$ are encoded in the (rotation) matrix $R$ <cit.>. The main advantage of representing rotations with $R$, as opposed to Euler angles, lies in the fact that the attitude is uniquely defined globally; i.e., every physical orientation of our robot corresponds to a unique rotational matrix.

In classical mechanics, we define the motion of an object as the variation of its body-fixed reference frame over time concerning an inertial reference frame, i.e., we need to determine/measure the linear and angular velocities of the robot relative to $\mathcal{F}_E$. We always have the freedom to represent inertial variables in more practical (although no convenient sometimes) frames of coordinates. For example, an Inertial Measurement Unit (IMU) measures linear accelerations and angular velocities relative to $\mathcal{F}_E$, but the reading is given in $\mathcal{F}_B$ since the IMU is onboard the robot.

This figure illustrates a 3D unicycle robot moving at a linear velocity $v$ and rotating at an angular velocity $\omega$. These velocities are observed from the inertial frame $\mathcal{F}_E$ but represented in the non-inertial frame $\mathcal{F}_B$, which is fixed at the body of the robot. $\{x_e, y_e, z_e\}$ and $\{x, y, z\}$ are the vector bases of $\mathcal{F}_E$ and $\mathcal{F}_B$, respectively.

In our scenario, considering $v^E \in \mathbb{R}^3$ and $\omega^E \in \mathbb{R}^3$ as the linear and angular velocities of the robot measured from $\mathcal{F}_E$ respectively, we denote its representation in $\mathcal{F}_B$ as  $v = \begin{bmatrix}v_{x} & v_{y} & v_{z}\end{bmatrix}^\top  \in \mathbb{R}^3$ and $\omega = \begin{bmatrix}\omega_{x} & \omega_{y} & \omega_{z}\end{bmatrix}^\top  \in \mathbb{R}^3$. Thus, considering the rotation matrix $R$ as the linear transformation that, when applied to a vector represented in $\mathcal{F}_B$, returns its representation in $\mathcal{F}_E$, we have that

\begin{equation} \label{eq: R_transform}
v^E = R v, \quad \omega^E = R \omega,
\end{equation}

where we can easily extract the kinematics for the position $p$ given by

\begin{equation}
\dot p(t) = R(t)v.
\end{equation}

The evolution of the signal $R(t)$ depends on the actuation $\omega(t)$ and can be straightforwardly derived from the Rodrigues' formula for the rotation of an angle $\theta\in\mathbb{R}$ around a unit axis $l = \begin{bmatrix}l_x & l_y & l_z \end{bmatrix}^\top$, which can be written in terms of the matrix exponential as

\begin{equation} \label{eq: rod_exp}
R_l(\theta) = \exp(\theta S_l),
\end{equation}

where $S_l = \left[\begin{smallmatrix}0 & -l_z & l_y \\
l_z & 0 & -l_x \\
-l_y & l_x & 0
\end{smallmatrix}\right]$; therefore, we have that

\begin{equation}
\label{eq: Omega11}
\frac{\mathrm{d}}{\mathrm{dt}}R_l(\theta) = \dot\theta S_l \exp(\theta S_l) = \exp(\theta S_l) \dot\theta S_l = R_l(\theta)\Omega,
\end{equation}

\begin{equation}
\label{eq: Omega1}
\Omega = \begin{bmatrix} 0 & -\omega_{z} & \omega_{y} \\ \omega_{z} & 0 & -\omega_{x} \\ -\omega_{y} & \omega_{x} & 0 \end{bmatrix}.
\end{equation}

Now, we are ready to write the kinematics of our 3D robot as

\begin{equation}
\label{eq: kin}
\begin{cases}
\begin{array}{ll}
\dot p(t) = R(t) v \\
\dot R(t) = R(t) \Omega \\
\end{array},
\end{cases}
\end{equation}

where $\Omega$ through $\omega$ is our actuation over the robot since we decided that $v$ will be constant.

In order to focus on the attitude control of the robot, we will consider $v$ in (<ref>) constant; in particular, $v = \begin{bmatrix}1 & 0 & 0\end{bmatrix}^\top.$

Finally, we are ready to pose the following question.

Given $R(t)$ as the signal for the 3D attitude representation of our robot and $R_a(p,t)$ as the attitude target. What should be the actuation $\Omega(t)$ in (<ref>) through $\omega$ such that $R(t) \to R_a(p,t)$ as $t\to\infty$?

In particular, we will construct $R_a(p,t)$ from a 3D vector field, and first, we will cover the trivial case where $v = 0$. The motivation for considering $R_a(p,t)$ comes from the necessity of some robotic problems where a mobile robot needs to track a given non-constant direction given by a 3D vector field, e.g., in order to track a path [5] or to track the source of a scalar field [6].

§ BRIEF INTRODUCTION TO LIE GROUPS AND $\TEXT{SO}(3)$

A group is a set $G$ of elements equipped with a binary operation “$\cdot$” that combines two elements $a,b\in G$ and outputs another element $c$ in the same set $G$, i.e., $a\cdot b = c$. This operation must be associative and must count with an identity element $e$, i.e., $a\cdot e = a$, and for each element in $G$, there is an inverse so that the result of the binary operation is the identity, e.g., $a\cdot b = e$ where $b$, sometimes denoted $a^{-1}$, is the inverse of $a$.

A Lie group joins together the concepts of group and smooth manifold, a differentiable topological space [7]. Sometimes, the elements of a Lie group can be associated with a physical meaning, e.g., the orientation of a reference frame, and that is why it makes them interesting in robotics. The local properties of a smooth manifold allow a coordinate system around any point (an element of the group). This fact enables us to identify a vector space around every element of the group, also known as tangent space. For example, the Earth's surface is a manifold, and its flat maps are tangent spaces. The tangent space at the identity of the Lie group, endowed with a bilinear nonassociative product called the Lie bracket, is known as the Lie algebra. Indeed, a Lie group is almost completely determined by its Lie algebra. Thus, for many purposes, one may operate directly in the Lie algebra instead of with the elements of its Lie group.

In this note, we focus on the Lie group known as $\text{SO}(3)$. Each element of this group represents a rotation in 3D, and all possible rotations are present. In particular, each element $R$ of the group is a linear transformation $\mathbb{R}^3 \rightarrow \mathbb{R}^3$; indeed, the element $R$ is also represented and known as a rotation matrix. It is a Lie group since we can associate each element $R$ of $\text{SO}(3)$ with a point in a smooth manifold. Although it is beyond the scope of this note to show it, for the sake of completeness, we note that the manifold $\text{SO}(3) \simeq \mathrm{RP}^3$ and its fundamental group is $\pi\left(\text{SO}(3)\right) = \mathbb{Z}/2\mathbb{Z}$.

§.§ The special Orthogonal Group $\text{SO}(3)$

Let us define the group $\text{SO}(3)$ as

\begin{equation}
\label{eq: SO3}
\text{SO}(3) := \{R \in \mathbb{R}^{3 \times 3} \,|\, R^\top R = I,\, \det(R) = 1\}.
\end{equation}

Since $\text{SO}(3)$ is also a smooth manifold, i.e., it is a Lie group, there exists a tangent space at every point $R$, and the velocity of $R$, i.e., $\dot R$, belongs to such a space. We can calculate $\dot R$ by taking the time derivative of the first constraint in (<ref>) as

\begin{align}
\dot R^\top R + R^\top\dot R &= 0 \nonumber \\
\Leftrightarrow \dot R &= - R \dot R^\top R \nonumber \\
\Leftrightarrow \dot R &= R \Omega, \label{eq: RdotOmega}
\end{align}

where $\Omega$ has already been found in (<ref>), and we note that it is indeed skew-symmetric since $\dot R ^\top R = -(\dot R ^\top R)^\top$. Looking at the definition, the Lie algebra is the tangent space associated with $\text{SO}(3)$ at the identity element of the group, i.e., setting $R=I$ in (<ref>) we have that $\Omega$ is in the Lie algebra of $\text{SO}(3)$. We formally denote the tangent spaces of $\text{SO}(3)$ as

\begin{equation}
\mathrm{T}_R\text{SO}(3) := \{R\Omega \, | \, \Omega \in \mathfrak{so}(3)\},
\end{equation}

where $\mathfrak{so}(3)$ is the vector space of the Lie algebra. This vector space has dimension $3$ because each of its elements can be written as

\begin{equation}
\Omega = \omega_x \left[\begin{smallmatrix}0 & 0 & 0 \\
0 & 0 & -1 \\
0 & 1 & 0
\end{smallmatrix}\right] + \omega_y \left[\begin{smallmatrix}0 & 0 & 1 \\
0 & 0 & 0 \\
-1 & 0 & 0
\end{smallmatrix}\right] + \omega_z \left[\begin{smallmatrix}0 & -1 & 0 \\
1 & 0 & 0 \\
0 & 0 & 0
\end{smallmatrix}\right],
\end{equation}

where the coordinates can be identified as the three angular velocities in $\omega$. Note that there is a bijective map between $\mathfrak{so}(3)$ and $\mathbb{R}^3$, i.e., $\mathfrak{so}(3) \cong \mathbb{R}^3$. This fact allows us to (re)define this vector space as

\begin{equation}
\mathfrak{so}(3) := \{x^\wedge \in \mathbb{R}^{3 \times 3} \, | \, x \in \mathbb{R}^{3},\, x^\wedge = - (x^\wedge)^\top \},
\end{equation}

where we have introduced the linear map hat $\wedge : \mathbb{R}^3 \rightarrow \mathfrak{so}(3)$ as the transformation of a given vector $x \in \mathbb{R}^3$ to a skew-symmetric matrix $x^\wedge \in \mathbb{R}^{3\times 3}$ similarly as we have constructed $\Omega$ in (<ref>) from $\omega$. This linear map is interesting since we have the following identity

\begin{equation} \label{eq: cross_identity}
x^\wedge y = x \times y,
\end{equation}

for any $y \in \mathbb{R}^3$, where $\times$ denotes the standard cross product.

To finish the bijection between $\mathfrak{so}(3)$ and $\mathbb{R}^3$, for the recovery of $\omega$ from a given $\Omega$ we define the linear map vee $\vee : \mathfrak{so}(3) \rightarrow \mathbb{R}^3$ as the opposite of the hat map, i.e., $\omega = \Omega^\vee$.

§.§ Geodesics

Given two elements $q_a,q_b \in M$, where $M$ is a smooth manifold and $t_a,t_b \in \mathbb{R^+}$ with $t_a < t_b$, we let

\begin{align}
\mathcal{C}(q_a, q_b) := \{ \gamma(t): [t_a,t_b] \rightarrow M \; | \; & \gamma(t_a) = q_a, \gamma(t_b) = q_b, \nonumber \\ &\gamma \in C^1\},
\end{align}

be the set of $C^1$-curves $\gamma$ in $M$ which connects $q_a$ with $q_b$. We define the Riemannian length of a curve $\gamma$ as

\begin{equation} \label{eq: rieman_len}
\ell (\gamma) = \int_{t_a}^{t_b} \sqrt{\langle \; \dot\gamma(t), \dot\gamma(t)\; \rangle_{\gamma(t)}} \, \mathrm{d}t,
\end{equation}

where $\langle \cdot , \cdot \rangle_x : T_x M \times T_x M \rightarrow \mathbb{R}$ is an inner product on $T_x M$ of $M$. A Riemannian metric $\langle \cdot , \cdot \rangle$ of $M$ is a smoothly chosen inner product $\langle \cdot , \cdot \rangle_x$ on each of the tangent spaces $T_x M$ of $M$.

Among all the curves $\gamma \in \mathcal{C}(q_a, q_b)$, we are interested in those with the minimum length since they define a Rimannian distance between $q_a,q_b \in M$ as

\begin{equation} \label{eq: rieman_dis}
d_M(q_a,q_b) = \mathrm{inf}\{\ell(\gamma) \; | \; \gamma \in \mathcal{C}(q_a, q_b), \; q_a, q_b \in M\}.
\end{equation}

In the following example, we show that this general definition of distance matches with the particular case of Euclidean distance.

Consider $M = \{(x,y) \in \mathbb{R}^2\}$, and two points $q_a, q_b \in \mathbb{R}^2$ connected by a curve $\gamma(t) = (x(t), y(t)) \subset \mathbb{R}^2$ such that $\gamma(0) = q_a = (x_a, y_a)$ and $\gamma(\Delta t) = q_b = (x_b, y_b)$. Since we are in $\mathbb{R}^2$, let $\langle \cdot , \cdot \rangle$ be the well-known inner product between two vectors. In this scenario, it is straightforward to identify that the curve $\gamma(t)$ with minimum length is the one with constant velocity $\dot \gamma(t)$ for every $t \in [0,\Delta t]$, i.e., a straight line. Moreover, the length of such a curve as defined in (<ref>) yields

\begin{align*}
\int_{0}^{\Delta t} \sqrt{\left(\dot x(t)\right)^2 + \left(\dot y(t)\right)^2} \mathrm{d}t &= \sqrt{\left(\dot x \Delta t\right)^2 + \left(\dot y \Delta t\right)^2}\\
&= \sqrt{\left(x_b - x_a\right)^2 + \left(y_b - y_a\right)^2},
\end{align*}

which is the Pythagorean equation. Therefore, we have that the length of our chosen $\gamma$ is in fact the Euclidean distance between $\gamma(0)$ and $\gamma(\Delta t)$.

The straight line within a flat manifold of Example <ref> is the most illustrative case of geodesic. Formally speaking, the geodesics are curves that allow us to go from $q_a$ to $q_b$ while maintaining a constant velocity in $M$; i.e., they are zero-acceleration curves <cit.>; of course, the acceleration here refers to within the manifold strictly. Indeed, it is possible to show for every manifold $M$ that the curves with minimum length are always geodesics [8], so we define the length-minimizing geodesic as follows.

Let $M$ be a manifold with a Riemannian metric $\langle \cdot , \cdot \rangle$. The length-minimizing geodesic, denoted as $\gamma_g$, is a $C^1$-curve $\gamma : [t_a,t_b] \rightarrow M$ with the property that $\ell(\gamma_g) = d_M(\gamma_g(t_a),\gamma_g(t_b))$.

This figure illustrates, for two different manifolds, a pair of geodesics ($\gamma$ in blue color and $\gamma_g$ in red color)  that connects $q_a, q_b \in \mathbb{T}^2$ (left) and $x,y \in S^2$ (right). Note that $\gamma_g$ is the length-minimizing geodesic.

Let us note that the definition of length in (<ref>) and the geodesics has a potential physical meaning. For any system whose state $x\in M$, we can design $\langle  \dot\gamma(t), \dot\gamma(t) \rangle$ according to their dynamics to represent kinetic energy along $\gamma$. Hence, if $\gamma$ is a length-minimizing geodesic, then $\ell (\gamma)$ minimizes the kinetic energy integral. Moreover, as the geodesics are zero acceleration curves, they represent the trajectories of free particles within the manifold, where motion is entirely determined by the curvature of $M$. This concept aligns with the fundamental idea in general relativity, where particles move along geodesics within space-time, and the curvature of such a manifold is induced by the presence of a mass (gravity).

Even when $M \subset \mathbb{R}^3$, e.g., $\mathbb{T}^2$ in <ref>, visualizing the geodesics is not straightforward. For an illustrative example, consider the unit sphere $S^2 := \{x \in \mathbb{R}^3 \, | \, \|x\| = 1\}$ and two points $x,y \in S^2$, as in <ref>. The geodesic $\gamma_g \in \mathcal{C}(x, y)$ can be interpreted geometrically as the (not necessarily unique) shortest arc in the unit sphere that connects $x$ with $y$, whose length, the geodesic distance in $S^2$ as defined in (<ref>), can be derived (see Appendix <ref>) as

\begin{equation} \label{eq: dis_s2}
d_{S^2}(x, y) = \arccos(x^\top y),
\end{equation}

where $\arccos(\cdot)$ takes values between $0$ and $\pi$.

In $\text{SO}(3)$ the geodesic is not that straightforward to identify, but we can follow the same steps as in $S^2$ to find a geodesic distance. Let $\langle \Omega_1, \Omega_2 \rangle = \mathrm{tr}(\Omega_1^\top \Omega_2) =: \langle\langle \Omega_1, \Omega_2 \rangle\rangle$ be our Riemannian metric for $\text{SO}(3)$, being $\Omega_1, \Omega_2 \in \mathfrak{so}(3)$. Consider a constant angular velocity $\omega$ that connects $R_a,R_b \in \text{SO}(3)$ along $\gamma_g \in \mathcal{C}(R_a, R_b)$, i.e., $\dot \gamma_g(t) = \omega$, such that $\gamma(0) = R_a$ and $\gamma(\Delta t) = R_b$. Then, computing the length of such a curve as in (<ref>) yields

\begin{align} \label{eq: len_so3}
\int_{0}^{\Delta t} \sqrt{\langle \dot \gamma(t), \dot \gamma(t) \rangle_{\gamma(t)}} \, \mathrm{d}t \ &= \Delta t \sqrt{\mathrm{tr}({\omega^\wedge}^\top {\omega^\wedge})} \nonumber \\
&= \Delta t \|\Omega\|_F = \|\tau^\wedge\|_F,
\end{align}

where $\|.\|_F$ is the Frobenius norm and $\tau^\wedge = \Omega \Delta t$ is the element of $T_{R_a}\text{SO}(3)$ that connects $R_a$ with $R_b$, as we will show in the following subsection. Indeed, let us note the following identity: since $\Omega = \omega^\wedge$ is a skew-symmetric matrix, then

\begin{equation} \label{eq: frob_product}
\|\omega^\wedge\|_F = \sqrt{\mathrm{tr}({\omega^\wedge}^\top \omega^\wedge)} = \sqrt{2(\omega \cdot \omega)},
\end{equation}

where $\cdot$ denotes the standard scalar product between two vectors in the Euclidean space.

The length computed in (<ref>) is the geodesic distance in $\text{SO}(3)$, expressed in terms of the integrated velocity vector $\tau$. Ideally, we would prefer to express this distance in terms of $R_a$ and $R_b$ directly, as we did in Example <ref> for a flat manifold, but to do so we have to introduce first the maps that enable us to establish such a relationship between $\mathfrak{so}(3)$ and $\text{SO}(3)$.

§.§ Exponential and logarithmic maps

The calculation of the velocity of $R$, as formulated in (<ref>), leads to an ordinary differential equation (ODE) whose solution is the signal

\begin{equation} \label{eq: exp_0}
R(t) = R(0) \exp(\omega^\wedge t) = R(0) \exp(\tau^\wedge (t)),
\end{equation}

where we denote $\tau(t) = \omega t$ as the integrated rotation vector or tangent vector. Since $R(t), R(0) \in \text{SO}(3)$, the associative property of the group yields

$$\exp(\tau^\wedge (t)) = R(0)^T R(t) \in \text{SO}(3).$$

Therefore, $\exp(\tau^\wedge (t))$ maps elements $\tau^\wedge (t)$ of the Lie algebra to the group. This is known as the exponential map. Intuitively, as <ref> illustrates for $S^1$, this map wraps the tangent vectors over the manifold following the length-minimizing geodesic from $x$ to $y$.

The $S^1$ manifold is a unit circle (blue) in $\mathbb{C}$, which contains every complex number $x$ such that $x^*x=1$. Its Lie algebra $T_1S^1$ is the line of imaginary numbers $i\mathbb{R}$ (red), and any tangent space $T_xS^1$ is isomorphic to the line $\mathbb{R}$ (red). Tangent vectors (red segments) wrap the manifold creating $\gamma_g$. The $\exp(\cdot)$ and $\log(\cdot)$ maps (arrows) wrap and unwrap elements of $\mathbb{R}$ to/from elements of $S^1$. This figure has been adapted from the reference [2].

By taking the absolutely convergent Taylor series

\begin{equation} \label{eq: exp_taylor}
\exp(\tau^\wedge) = \sum_{k=0}^\infty \frac{(\tau^\wedge)^k}{k!}
\end{equation}

and computing it as in [2], we formulate the $\exp(\cdot)$ map as

\begin{equation} \label{eq: exp}
\exp(\tau^\wedge) = I + \frac{\sin\theta}{\theta} \tau^\wedge +\frac{1 - \cos\theta}{\theta^2} ({\tau^\wedge})^2,
\end{equation}

where $\theta = \|\tau\|$ is denoted as the integrated rotation angle. This closed form of the exponential map is the well-known Rodrigues' formula for rotation matrices, as previously introduced in (<ref>).

The inverse of $\exp(\cdot)$ is what we know as the logarithmic map $\log(\cdot)$ and its closed form can be obtained as follows. For the sake of clarity, let $R(0) = I$, i.e., $\tau^\wedge \in \mathfrak{so}(3)$, so we directly have that $\exp(\tau^\wedge(t)) = R(t)$. Therefore, we can consider the skew-symmetric/symmetric decomposition of $R$ to regroup the terms in (<ref>) as

\begin{equation}\label{eq: skew_symm}
R = \underbrace{\frac{\sin\theta}{\theta} \tau^\wedge}_{\frac{R - R^T}{2}} + \underbrace{ I +  \frac{1 - \cos\theta}{\theta^2} ({\tau^\wedge})^2}_{\frac{R + R^T}{2}}.
\end{equation}

Thus, we have that $\frac{\sin\theta}{\theta} \tau^\wedge = \frac{R - R^T}{2}$, which directly yields

\begin{equation} \label{eq: log}
\tau^\wedge = \log(R) =
\begin{cases}
\begin{array}{ll}
0_{3 \times 3} & R = I\\
\left( \frac{\theta}{2 \sin\theta}  \right) (R - R^\top ) & R \neq I\\
\end{array}.
\end{cases}
\end{equation}

Here, $\tau$ is unknown and we cannot compute $\theta$ from it directly, but note that taking the trace of $R$ in (<ref>) leads us to $\cos\theta = \frac{1}{2}(\mathrm{tr}(R) - 1)$, which restricting $\theta \in [0,\pi]$ gives the unique solution

\begin{equation} \label{eq: arccos}
\theta = \arccos\left(\frac{\mathrm{tr}(R) - 1}{2}\right).
\end{equation}

At the boundary $\theta = \pi$, $R$ is purely symmetric, i.e., $R = R^T$, so the skew matrix provides no information, and thus, the solution of $\log(\cdot)$ is not unique. Therefore, we have a unique $\log(\cdot)$ for all $\tau$ within the interior ball of radius $\pi$, i.e., when $\|\tau\| < \pi$. Taking this into consideration, we formally define the exponential and logarithmic linear maps as

\begin{align*}
\exp :& \; \Pi_{\mathfrak{so}(3)} \rightarrow \Pi_{\text{SO}(3)} \nonumber \\
\log :& \; \Pi_{\text{SO}(3)} \rightarrow \Pi_{\mathfrak{so}(3)},
\end{align*}

whose closed formed are given by (<ref>) and (<ref>), respectively. Here, $\Pi_{\mathfrak{so}(3)} := \{\tau^\wedge \in \mathfrak{so}(3) \, : \, \tau \in \mathbb{R}^{3},\, \|\tau\| < \pi\}$ and $\Pi_{\text{SO}(3)} := \{R \in \text{SO}(3) \, : \, \mathrm{tr}(R) \neq -1\}$.

When $\tau^\wedge \in T_{R_a} \text{SO}(3)$, we have that $R(0) = R_a$ in (<ref>), i.e., $\exp(\tau^\wedge(t)) = R_a^\top R(t)$. Computing the inverse of this expression, as we did for the trivial case of $R(0)=I$, leads us to $\tau^\wedge(t) = \log(R_a^\top R(t))$. This same scenario is illustrated for $S^1$ in <ref> when $i\theta \in T_xS^1$. Therefore, revisiting (<ref>), we can now formalize the geodesic distance in $\text{SO}(3)$ in terms of $R_a$ and $R_b$ as

\begin{equation} \label{eq: metric}
d_{\text{SO}(3)}(R_a, R_b) = \|\log(R_a^\top R_b)\|_F,
\end{equation}

where $R_a^\top R_b \in \Pi_{\text{SO}(3)}$, i.e., $\mathrm{tr}(R_a^\top R_b) \neq -1$. Note that ensuring $\mathrm{tr}(R_a^\top R_b) \neq -1$ is not straightforward, but there are particular situations where it becomes easier, e.g., see the following technical result.

Given $R_a(t), R_b(t) \subset \text{SO}(3)$ such that

\begin{equation} \label{eq: dso3_exp}
d_{\text{SO}(3)}(R_a(t), R_b(t)) \leq d_{\text{SO}(3)}(R_a(0), R_b(0)),
\end{equation}

where $t\in \mathbb{R}^+$. If $\mathrm{tr}(R_a(0)^\top R_b(0)) \neq -1$, then $\mathrm{tr}(R_a(t)^\top R_b(t)) \neq -1$ for all $t$.

From (<ref>) we have that $d_{\text{SO}(3)}(R_a(t), R_b(t))$ is bounded by the value at $t=0$.  On the other hand, from (<ref>) we have that $R$ is the furthest element in $\text{SO}(3)$ from $I$ when $\mathrm{tr}(R) = -1$. Therefore, the distance between $R_a(t)$ and $R_b(t)$ will not ever be the maximum one; hence, the claim.

§.§ Computing the exponential and logarithmic maps

In practice, numerous methods exist for computing these two maps, and its choice may depend on practical requirement details. Since in $\text{SO}(3)$ we have the Rodrigues' formula, the implementation of the exponential map reduces to computing $\theta = \|\tau\|$ and applying (<ref>). Nonetheless, when $\theta \rightarrow 0$, the denominators of the coefficients in (<ref>) also approach zero, so the computation of $\exp(\cdot)$ may be compromised numerically. In such cases, since we are dealing with small angles, it is convenient to truncate the series (<ref>) up to a suitable order.

When we are working with the logarithmic map, we also have to take special care about these singularities. Implementing $\log(\cdot)$ also reduces to compute $\theta$ as in (<ref>) and apply (<ref>), but when $\theta \rightarrow 0$ we have the same numerical problem. In that case, it is convenient to use the Taylor's series expansion

\begin{equation*}
\left( \frac{\theta}{\sin\theta}  \right) = 1 + \frac{1}{6}\theta^2 + \frac{7}{360}\theta^4 + \mathcal{O}(\theta^6).
\end{equation*}

On the other hand, we have also seen that $\log(\cdot)$ is not well defined when $\theta = \pi$, but we can computationally deal with that as follows. When $\theta = \pi$, the symmetric part of (<ref>) yields $R = I + 2(\tau^\wedge)^2$, which can be rewritten as

\begin{equation}
\tau^\wedge {\tau^\wedge}^\top = \frac{1}{2}(R + I)
\end{equation}

by considering that $(\tau^\wedge)^2 = \tau^\wedge {\tau^\wedge}^\top - I$.
Thus, we can still determine $\tau$ up to a factor of $\pm 1$.

The reader can explore an implementation example in Python for both the exponential and logarithmic maps in the GitHub repository [9]. The method we have chosen in this project to compute the logarithmic map is similar to the approach presented in [10].

§.§ The adjoint map

The relationship between $\Omega \in \mathfrak{so}(3)$ and its representation in $\mathcal{F}_E$ as in (<ref>), denoted as $\Omega^E$, can be seen as operating with a different basis. Indeed, knowing that $R^{-1} = R^T$, we have that

\begin{equation} \label{eq: r_omega_rt}
\Omega^E = R \Omega R^{T},
\end{equation}

that is, we first change the basis from $\mathcal{F}_E$ to $\mathcal{F}_B$ with $R^T$, then we apply $\Omega$, and finally we come back to $\mathcal{F}_E$ with $R$. Of course, we note that $\Omega^E \in \mathfrak{so}(3)$ as it comes from a linear transformation. In the context of Lie groups, this linear transformation is the well-known adjoint map, defined as the linear map $\mathrm{Ad}_R : \mathfrak{so}(3) \rightarrow \mathfrak{so}(3)$ such that $\mathrm{Ad}_R(\Omega) = R \Omega R^T$. Indeed, taking the vee map of (<ref>) leads to (<ref>) by definition, i.e.,

\begin{equation} \label{eq: adjoint_vee}
{\mathrm{Ad}_R(\Omega)}^\vee = R \Omega^\vee = R\omega,
\end{equation}

since $(\Omega^E)^\vee = \omega^E$. This is a powerful identity of the adjoint map that we will extensively use later in our analysis.

Taking the derivative of the adjoint map at $R=I$ gives the adjoint operator $\mathrm{ad}_\Omega : \mathfrak{so}(3) \rightarrow \mathfrak{so}(3)$. We show next in Example <ref> that the definition of this operator can be formulated as

\begin{equation} \label{eq: ad_op_def}
\mathrm{ad}_\Omega = \frac{\mathrm{d}}{\mathrm{d}t} \mathrm{Ad}_{e^{t\Omega}}\bigg|_{t=0},
\end{equation}

where $e^{t\Omega} := \exp(t\Omega)$ is a length-minimizing geodesic curve, and that

\begin{equation} \label{eq: ad_op}
\mathrm{ad}_{\Omega_a} (\Omega_b) = [\Omega_a, \Omega_b] = \Omega_a \Omega_b - \Omega_b \Omega_a,
\end{equation}

where $[.\,,\,.] : \mathfrak{so}(3) \times \mathfrak{so}(3) \rightarrow \mathfrak{so}(3)$ is the Lie bracket operator of the Lie algebra. Note that this expression directly leads us to $\mathrm{ad}_{\Omega_a}(\Omega_b) = 0$ when $\Omega_a = \Omega_b$. Moreover, if the Lie algebra is $\mathfrak{so}(3)$, then (<ref>) yields $[\omega^\wedge_a, \omega^\wedge_b] = (\omega_a \times \omega_b)^\wedge$. Finally, considering the Taylor expansion of the exponential, note that the adjoint map can be written as

\begin{align} \label{eq: Ad_series}
\mathrm{Ad}_{\exp(\Omega_a)}(\Omega_b) &= \sum^\infty_{m=0}\frac{(\mathrm{ad}_{\Omega_a})^m}{m!}(\Omega_b) \nonumber\\
&= \Omega_b + [\Omega_a,\Omega_b] + \frac{1}{2}[\Omega_a,[\Omega_a,\Omega_b]] + \text{H.O.T.}
\end{align}

Consider two signals $R_a(t), R_b(t) \in \text{SO}(3)$ such that $R_a(0) = R_b(0) = I$, $\dot R_a(0) = \Omega_a$ and $\dot R_b(0) = \Omega_b \in \mathfrak{so}(3)$; e.g., the length-minimizing geodesics $R_a(t) = \exp(t\Omega_a)$ and $R_b(t) = \exp(t\Omega_b)$. Let us choose an arbitrary $s\in \mathbb{R}^+$ and define the curve
$C_s(t) = \mathrm{Ad}_{R_a(s)}(R_b(t)) = R_a(s) R_b(t) R_a(s)^\top.$
When $t=0$, note that $C_s(t)$ is contained in $\mathfrak{so}(3)$, and thus its derivative with respect to $t$ or $s$. Hence, the derivative at $t = s = 0$ yields

\begin{align} \label{eq: ex_lie_bracket}
\scalemath{0.9}{
\frac{\mathrm{d}}{\mathrm{d}s}}&
\scalemath{0.9}{
\left( R_a(s) \Omega_b R_a^\top(s) \right)\bigg|_{s=0} = \frac{\mathrm{d}}{\mathrm{d}s}\left( \exp({s\Omega_a}) \Omega_b {\exp({s\Omega_a})}^\top \right)\bigg|_{s=0}} \nonumber\\
&= \dot R_a(0) \Omega_b R_a^\top(0) + R_a(0) \Omega_b \dot R_a^\top(0) \nonumber \\
&= \Omega_a \Omega_b - \Omega_b \Omega_a \nonumber \\
&= [\Omega_a, \Omega_b] \in \mathfrak{so}(3),
\end{align}

thus $[\Omega_a, \Omega_b]$ is an interesting operation of $\mathfrak{so}(3)$. Indeed, from (<ref>) one can directly derive (<ref>) and (<ref>).

§.§ Adjoint invariant metric

We say that a Riemannian metric $\langle . , .\rangle$ is Ad-invariant when

\begin{equation} \label{eq: ad-inv}
\langle \mathrm{Ad}_R(\Omega_1) , \mathrm{Ad}_R(\Omega_2) \rangle = \langle \Omega_1 , \Omega_2 \rangle.
\end{equation}

For example, considering the metric $\langle\langle \Omega_1, \Omega_2 \rangle\rangle = \mathrm{tr}(\Omega_1^\top \Omega_2)$ we have that

\begin{align*}
\langle\langle \mathrm{Ad}_R(\Omega_1) , \mathrm{Ad}_R(\Omega_2) \rangle\rangle &= \mathrm{tr}(R\Omega_1^TR^T R \Omega_2 R^T) \\
&= \mathrm{tr}(\Omega_1^T \Omega_2) = \langle\langle \Omega_1 , \Omega_2 \rangle\rangle.
\end{align*}

The metrics with this property will be particularly interesting in our analysis because they lead to the following identity.

<cit.> Considering a Riemannian metric $\langle . , .\rangle$ of a smooth manifold $M$, if it is Ad-invariant then

\begin{equation}
\langle \mathrm{ad}_X(Y) , Z \rangle = - \langle Y , \mathrm{ad}_X(Z) \rangle,
\end{equation}

where $X,Y,Z \in T_x M$, $x\in M$.

Given the definition of adjoint operator as in (<ref>), we have that

\begin{align*}
\langle \mathrm{ad}_{X}(Y) , Z \rangle &= \langle \frac{\mathrm{d}}{\mathrm{d}t} \mathrm{Ad}_{e^{tX}}(Y)\bigg|_{t=0} , Z \rangle  \\
&= \frac{\mathrm{d}}{\mathrm{d}t}\langle  \mathrm{Ad}_{e^{tX}}(Y) , Z \rangle \bigg|_{t=0}.
\end{align*}

Here, since $\langle . , . \rangle$ is Ad-invariant, we can apply (<ref>), which yields

\begin{align*}
\langle  \mathrm{Ad}_{e^{tX}}(Y) , Z \rangle &= \langle \mathrm{Ad}_{e^{-tX}}(\mathrm{Ad}_{e^{tX}}(Y)) , \mathrm{Ad}_{e^{-tX}}(Z) \rangle\\
&= \langle  Y , \mathrm{Ad}_{e^{-tX}}(Z) \rangle.
\end{align*}

Therefore, we can directly conclude that

\begin{align*}
\langle \mathrm{ad}_{X}(Y) , Z \rangle &= \frac{\mathrm{d}}{\mathrm{d}t}\langle  \mathrm{Ad}_{e^{tX}}(Y) , Z \rangle \bigg|_{t=0} \\
&= \frac{\mathrm{d}}{\mathrm{d}t}\langle  Y , \mathrm{Ad}_{e^{-tX}}(Z) \rangle \bigg|_{t=0} \\
&= \langle Y , -\mathrm{ad}_{X}(Z) \rangle = - \langle Y , \mathrm{ad}_{X}(Z) \rangle.
\end{align*}

§.§ The Jacobian of the exponential map in $\text{SO}(3)$

Given $R(t)$ as a smooth curve of $\text{SO}(3)$, we now want to discover the relationship between the time derivative of $\tau^\wedge(t) = \log(R(t))$ and the velocity tensor of $R(t)$ from $\mathcal{F}_B$, denoted as $\Omega = \dot R(t) R(t)^\top$. One may show as in <cit.> that

\begin{equation} \label{eq: omega_tau_int}
\Omega = \int_0^1 \mathrm{Ad}_{e^{-\lambda \tau^\wedge(t)}}(\dot\tau^\wedge) \, \mathrm{d}\lambda,
\end{equation}

which, looking at the identity (<ref>), can be interpreted as the integration of $\dot\tau$ along the length-minimizing geodesic that connects $I$ with $\exp(-\tau^\wedge(t))$. Following the proof of <cit.> we can compute the solution of the integral in (<ref>), and it yields

\begin{equation} \label{eq: series_derv}
\dot\tau^\wedge = \sum^\infty_{n=0} \frac{(-1)^n \mathrm{B}_n}{n!} \mathrm{ad}_{\tau^\wedge(t)}^n (\Omega),
\end{equation}

where $\{\mathrm{B}_n\}$ are the Bernoulli numbers. In fact, <cit.> shows that in $\text{SO}(3)$ the matrix series (<ref>) converges to

\begin{align} \label{eq: tau_dot}
\dot\tau^\wedge &= \mathcal{B}_{\tau^\wedge(t)} (\Omega) \nonumber \\
&= \Omega + \left(\frac{1}{2}\mathrm{ad}_{\tau^\wedge(t)} + \frac{1 - \alpha(\theta)}{\theta^2}\mathrm{ad}^2_{\tau^\wedge(t)} \right) (\Omega),
\end{align}

where $\alpha(\theta) := (\theta/2)\cot(\theta/2)$. This expression will be useful for us thanks to the following identity

\begin{equation}
\Omega = \tau^\wedge(t) \implies \dot\tau^\wedge = \Omega,
\end{equation}

since $\mathrm{ad}_{\Omega}(\Omega) = 0$.

§ EXPONENTIAL STABILITY AND CONTROL IN $\TEXT{SO}(3)$

A function $f : \mathbb{R}^+ \rightarrow \mathbb{R}^n$ is said to converge to $r_0 \in \mathbb{R}^n$ exponentially fast if there exist constants $\alpha, k >0$ and a time $t_0 \in \mathbb{R}^+$, such that for all $t > t_0$

$$\|f(t) - r_0\| \leq \alpha \|f(t_0) - r_0\| e^{-k (t - t_0)}.$$

Hence, considering $f(t)$ as a curve $\gamma(t) \subset M$, we extend this concept of exponential convergence to curves within smooth manifolds as follows.

Given the signal $\gamma(t) \in M$, we say that the tracking of the desired signal $\gamma_d(t) \in M$ is exponentially stable for the set of initial conditions $\mathcal{U}_\epsilon(\gamma_d(t_0)) := \{ x \in M \,|\, d_M(x, \gamma_d(t_0)) \leq \epsilon\}$, where $\epsilon > 0$ and $t_0 \in \mathbb{R}^+$, if there exist constants $\alpha, k >0$ such that

\begin{equation} \label{eq: exp_stability}
d_M(\gamma(t), \gamma_d(t)) \leq \alpha d_M(\gamma(t_0), \gamma_d(t_0)) e^{-k(t-t_0)}
\end{equation}

for all $\gamma(t_0) \in \mathcal{U}_\epsilon(\gamma_d(t_0))$ and $t \geq t_0$.

Therefore, from (<ref>) we have that $\lim_{t\rightarrow + \infty}\gamma(t) = \gamma_d(t)$; in other words, we can solve Problem <ref> by designing a control law such that, given a local set of initial conditions $\mathcal{U}_\epsilon(R_a(t_0)) \subset \text{SO}(3)$, the tracking of $R_a(t)$ is exponentially stable for all $t>t_0$. For the sake of clarity, we have denoted $R(p(t),t) = R(t)$.

Let us begin by defining an error signal. Given the metric (<ref>) and a target attitude matrix

\begin{equation} \label{eq: Ra}
R_a := \left[x_a \; y_a \; z_a\right] \in \text{SO}(3),
\end{equation}

where $\{x_a, y_a, z_a\}$ is an orthonormal basis of $\mathbb{R}^3$, we define the attitude error as

\begin{equation} \label{eq: mu_r}
\mu_{R_e} := d_{\text{SO}(3)}(R_a, R) = \|\log(R_e)\|_F,
\end{equation}

where $R_e = R_a^\top  R \in \text{SO}(3)$ is denoted as the attitude error matrix. Using $R_e$ in the control input of (<ref>) yields the  following feedback system

\begin{equation} \label{eq: kin_fb}
\begin{cases}
\begin{array}{l}
\dot p = R v \\
\dot R = R \Omega(R_e) \\
\dot R_e = R_e \Omega_e
\end{array},
\end{cases}
\end{equation}

where $\Omega_e \in \mathfrak{so}(3)$ denotes the velocity tensor of the attitude error matrix, which may depend on the velocity tensor of the desired attitude matrix, defined as $\Omega_a := R_a^\top \dot R_a \in \mathfrak{so}(3)$. Note that by taking the time derivative of $R_e = R_a^\top  R$

\begin{align} \label{eq: R_e_dot}
\dot R_e &= R_a^{\top}\dot R + \frac{\mathrm{d}}{\mathrm{dt}}(R_a^{\top})R \nonumber\\
&= R_a^{\top} \dot R - R_a^{\top} \dot R_a R_a^{\top} R \nonumber\\
&= R_e(\Omega - R^{\top}\dot R_a R_a^{\top} R),
\end{align}

where we have considered that $\frac{\mathrm{d}}{\mathrm{dt}}(R^{\top}) = \frac{\mathrm{d}}{\mathrm{dt}}(R^{-1}) = R^{-1}\dot R R^{-1} = R^{\top}\dot R R^{\top}$. Thus, since $\dot R_e = R_e \Omega_e$, from (<ref>) we can formulate the angular velocity tensor of $R_e$ as

\begin{align} \label{eq: omega_e}
\Omega_e &= \Omega - R^{T}\dot R_a R_a^{T} R \nonumber\\
&= \Omega - R_e^{T} \Omega_a R_e \nonumber\\
&= \Omega - \mathrm{Ad}_{R_e^{T}}(\Omega_a),
\end{align}

which by taking the vee map while considering (<ref>) yields

\begin{align} \label{eq: omega_e_vee}
\Omega_e^\vee = \Omega^\vee - R_e^{T}\Omega_a^\vee,
\end{align}

therefore, we can express the velocity tensor of $R_e$ as a vector in $\mathbb{R}^3$ by taking the difference between $\Omega^\vee$ and $R_e^\top \Omega_a^\vee$.

Once we have introduced our error signal $\mu_{R_e(t)}$ and the feedback system (<ref>), we are in place to formalize a control law that will allow a robot modeled by (<ref>) to converge to a desired time-varying target attitude $R_a(t)$.

Consider the attitude target $R_a(t)$ and the feedback system (<ref>), for any initial condition $R(t_0)$ such that $\mathrm{tr}(R_e(t_0)) \neq -1 $, the tracking of $R_a(t)$ is exponentially stable under the control law

\begin{align} \label{eq: omega_control}
\Omega(R_e) &= - k_w \log(R_a^{T}R) + R^{T} \dot R_a R_a^{T} R \nonumber \\
&= - k_w \log(R_e) + \mathrm{Ad}_{R_e^\top } (\Omega_a),
\end{align}

where $k_w \in \mathbb{R}^+$ is a positive constant that modulates the aggressiveness of the proportional controller.

Given the Lyapunov function $V(R_e) = \mu_{R_e}^2/2$, we have that

\begin{equation} \label{eq: V}
V(R_e) = \frac{1}{2} \|\log(R_e)\|^2 = \frac{1}{2} \mathrm{tr}(\log(R_e)^\top \log(R_e)),
\end{equation}

which is well-defined for all $R_e$ if $\mathrm{tr}(R_e) \neq -1$. Taking the time derivative of this expression leads us to

\begin{align} \label{eq: V_dot_B}
\dot V(R_e) &= \mathrm{tr}(\log(R_e)^\top \frac{\mathrm{d}}{\mathrm{dt}}\log(R_e)) \nonumber\\
&= \mathrm{tr}(\log(R_e)^\top \mathcal{B}_{\log(R)}(\Omega_e)) \nonumber \\
&= \mathrm{tr}(\log(R_e)^\top  \Omega_e),
\end{align}

where we have used Lemma <ref> so that

\begin{align*}
\langle\langle \mathrm{ad}_{\log(R)}(\Omega) , \log(R) \rangle\rangle &= - \langle\langle \Omega , \mathrm{ad}_{\log(R)}(\log(R)) \rangle\rangle\\
\langle\langle \mathrm{ad}^2_{\log(R)}(\Omega) , \log(R) \rangle\rangle &= \langle\langle \Omega , \mathrm{ad}^2_{\log(R)}(\log(R)) \rangle\rangle,
\end{align*}

and the fact that $\mathrm{ad}_{\log(R)}(\log(R)) = 0$. Therefore, if we let

\begin{align} \label{eq: omega_e_kw}
\Omega_e = - k_w \log(R_e),
\end{align}

then we have that $\dot V(R_e) = - 2 k_w V(R_e)$. The solution of this ODE is $V(R_e(t)) = V(R_e(t_0)) e^{-2 k_w t}$, which is equivalent to say that $\mu_{R_e}(t) = \frac{1}{\sqrt{2}} \mu_{R_e}(t_0) e^{-k_w t}$. Hence, given the Definition <ref>, we have that the tracking of $R_a(t)$ is exponentially stable if $\mathrm{tr}(R_e(t_0)) \neq -1$ because of Lemma <ref>.

In fact, this proposition enables us to present the following formal solution for Problem  <ref>.

Consider the attitude target $R_a(p(t),t) = R_a(t)$ such that $\dot R_a(t)$ is known. If there exist a time $t_0 \in \mathbb{R}^+$ such that $\mathrm{tr}(R_e(t_0)) \neq -1$, then taking $\Omega(t)$ as in (<ref>) solves Problem <ref> for every $t \geq t_0$.

§ PRACTICAL PROBLEMS

§.§ Aligning to a vector field with an unknown time variation

Given a, possibly time-varying, vector field $m_d(t) \in S^2$, whose time variation is unknown but bounded, i.e., $\|\dot m_d(t)\| \leq \omega_d$, where $\omega_d \in \mathbb{R}^+$, our objective now is to align the heading direction of our robot, i.e., $\hat {\dot p} := \dot p/\|\dot p\|$, with $m_d(t)$. Considering Remark <ref> and $R(t) = \begin{bmatrix}x & y & z \end{bmatrix}$ with $x,y,z$ as defined in (<ref>), we have that $\hat {\dot p}^\top x = 1$, and $\hat {\dot p}^\top y = \hat {\dot p}^\top z = 0$. Thus, this problem is equivalent for us to minimize the geodesic distance (<ref>) between $x(t)$ and $m_d(t)$. Indeed, we define the target attitude set as follows.

Consider $x(t)\in S^2, R_a(t) \in SO(3)$ with $R_a(t) = \begin{bmatrix}m_d(t) & * & *\end{bmatrix}$ and a positive constant $\delta \in [0,\pi]$, the target attitude set is defined as

\begin{equation*}
\mathcal{M}_\delta := \{R \,|\, \arccos(x(t)^\top m_d(t)) \leq \delta\}.
\end{equation*}

Note that we leave the desired attitude in the YZ axis up to the controller or designer. We now formally introduce the problem of this subsection as follows.

Consider a robot modeled by (<ref>) and a target attitude set $\mathcal{M}_{\delta^*}$ such that $\|\dot m_d(t) \| \leq \omega_d$. Then, design a control signal $\Omega(t)$ such that, for a given $\delta^* \in [0,\pi]$ then $R(t) \to \mathcal{M}_{\delta^*}$ as $t\to\infty$.

Firstly, let us note that, considering $R_a(t)$ as in (<ref>), with $x_a(t) = m_d(t)$, we can split the attitude target velocity into

\begin{equation} \label{eq: omega_a_unk}
\Omega_a = \Omega_a^u + \Omega_a^k,
\end{equation}

where $\Omega_a^u$ is an unknown variation of $R_a(t)$ associated to the time variation of $m_d (t)$, i.e., $\|{\Omega_a^u}^\vee\| \leq \omega_d$, and $\Omega_a^k$ can be configured through the design of the, possibly time-varying but known, vectors $y_a(t)$ and $z_a(t)$. Since we do not know $\Omega_a^u$, we reformulate the control law (<ref>) as

\begin{equation} \label{eq: omega_control_know}
\Omega(R_e) = - k_\omega \log(R_e) + \mathrm{Ad}_{R_e^\top} (\Omega_a^k),
\end{equation}

which leads us to the following technical result.

Given the feedback system (<ref>), an attitude target $R_a(t)$, whose velocity $\Omega_a$ is given by (<ref>), and a constant angular velocity $\omega_d \in \mathbb{R}^+$ such that $\|{\Omega_a^u}^\vee\| < \omega_d$. If $\mathrm{tr}(R_e(t_0)) \neq -1$, then the control law in (<ref>) makes $R(t) \to M_{\delta^*}$ exponentially fast as $t\to\infty$ with $k_\omega = \sqrt{2}\frac{\omega_d}{\mu^*_{R_e}}$, where $\mu^*_{R_e} \leq \delta^*$.

Substituting the vee form of (<ref>) and (<ref>) into (<ref>), we obtain

\begin{align} \label{eq: omega_e_unk}
\Omega_e^\vee &= \Omega^\vee - {R_e^\top}^\vee \Omega_a^\vee \nonumber \\
&= - k_w \log(R_e)^\vee - {R_e^\top}^\vee {\Omega_a^u}^\vee.
\end{align}

Consider the Lyapunov function $V = \frac{1}{4}\mu_{R_e}^2$ and its time derivative, similarly as in (<ref>)

\begin{equation*}
\frac{1}{4}\frac{\mathrm{d}}{\mathrm{dt}} \mu_{R_e}^2 =  \frac{1}{2}\mathrm{tr}\left(\log(R_e)^\top \Omega_e\right) =  \log(R_e)^\vee \cdot \Omega_e^\vee.
\end{equation*}

Therefore, given $\Omega_e^\wedge$ as in (<ref>) and considering (<ref>), we have that

\begin{equation} \label{eq: to_max}
\log(R_e)^\vee \cdot \Omega_e^\vee = - \frac{1}{2}k_w \mu_{R_e}^2 - \log(R_e)^\vee \cdot R_e^\top (\Omega_a^u)^\vee,
\end{equation}

where $\frac{1}{2}\mu_{R_e}^2 = \log(R_e)^\vee \cdot \log(R_e)^\vee$. We upper bound the second term of (<ref>) thanks to

$$\log(R_e)^\vee \cdot R_e^\top (\Omega_a^u)^\vee \leq \frac{\mu_{R_e}}{\sqrt{2}}\omega_d,$$

where $\|\log(R_e)^\vee\| = \frac{\mu_{R_e}}{\sqrt{2}}$ because of (<ref>); thus, the time derivative of $\mu_{R_e}^2$ can be bounded as

\begin{align*}
\frac{1}{2}\frac{\mathrm{d}}{\mathrm{dt}} \mu_{R_e}^2 \leq \left(-k_\omega  + \sqrt{2}\frac{\omega_d}{\mu_{R_e}} \right) \mu_{R_e}^2.
\end{align*}

Therefore, if $\mathrm{tr}(R_e(t_0)) \neq -1$, as seen in Lemma <ref>, and knowing that $\mu_{R_e}(t_0)>0$ because it is a distance, then $\mu_{R_e}(t)$ converges exponentially fast as $t\to\infty$ to the interval $[0, \mu_{R_e}^*]$ for $k_\omega = \sqrt{2}\frac{\omega_d}{\mu^*_{R_e}}$, where $\mu^*_{R_e} \leq \delta^*$.  Let us see why this selection of $\mu^*_{R_e}$ is sufficient. First, note that $\mu_{R_e}$ is the length-minimizing geodesic distance, so it is equal to the integrated rotation angle $\theta$, i.e., $\log(R_e)^\vee = \mu_{R_e} l$, where $l\in\mathbb{R}^3$ is the unit vector along the rotation axis. Next, let us decompose the axis $x$ into components parallel and perpendicular to the axis $l$,

\begin{equation*}
x = x_\parallel + x_\perp = (x^\top l) l - l \times (l \times x).
\end{equation*}

Applying the rotation $\log(R_e)^\vee$ to this vector yields [13]

\begin{align*}
x_\text{rot} &= x_\parallel + x_{\perp \text{rot}} \\
&= x_\parallel + \cos(\mu_{R_e}) x_\perp + \sin(\mu_{R_e})l \times v_\perp.
\end{align*}

Therefore, considering $x_\text{rot} = m_d$ we have that

\begin{align*}
\cos(\delta) = x^\top m_d &= \|x_\parallel\|^2 + \cos(\mu_{R_e})\|x_\perp\|^2\\
&= \cos(\mu_{R_e}) + \|x_\parallel\|^2(1 - \cos(\mu_{R_e})),
\end{align*}

so $\cos(\delta) \geq \cos(\mu_{R_e})$ and consequently $\delta \leq \mu_{R_e}$. Hence, $\mu^*_{R_e} \leq \delta^*$ ensures that $\delta \leq \delta^*$.

Now we are ready to provide a solution to Problem <ref>.

Given the feedback system (<ref>) and attitudes $R_a(t), R(t_0)$ meeting the conditions of Lemma <ref>, then the control law (<ref>) solves Problem <ref>.

§.§ Distances between robots while they align with a common and fixed attitude

Even though the following result looks cumbersome because we know the convergence rate to the desired fixed attitude, the result aims to show how to exploit certain norm properties.

Given two robots $i, j$ with dynamics (<ref>), $||v|| =1$ and a desired attitude matrix $R_a \in SO(3)$, if $\Omega_{\{i,j\}}$ is designed as in (<ref>) so that the rate of exponential convergence to $R_a$ is $r \in \mathbb{R}^+$, then

\begin{equation*}
\|p_{ij}(t) - p_{ij}(t_0)\| \leq \frac{2\sqrt{3}\pi}{r}.
\end{equation*}

Given $\dot p_{ij} = (R_i - R_j)v$, the definition of norm for linear applications leads to

\begin{equation} \label{eq: pij_rij}
\|\dot p_{ij}\| \leq \|R_i - R_j\|_F\|v\| = \|R_i - R_j\|_F.
\end{equation}

If we left multiply $(R_i - R_j)$ by $R_a R_a^\top$, the expression (<ref>) yields

\begin{align*}
\|R_i - R_j\|_F &\leq \|R_a\|_F\|R_a^\top R_i - R_a^\top R_j\|_F \nonumber\\
&= \sqrt{3}\|R_e^i - R_e^j\|_F
\end{align*}

Hence, since $\|R_i - R_j\|_F$ and (<ref>) are bounded equivalent metrics [14], we can consider

\begin{equation*}
\|R_e^i - R_e^j\|_F \leq \|\log({R_e^j}^\top R_e^i)\|_F,
\end{equation*}

which, by the triangle inequality, leads us to

\begin{align*}
\|\log({R_e^j}^\top R_e^i)\|_F &\leq \|\log(R_e^i)\|_F + \|\log(R_e^j)\|_F \nonumber\\
&= \mu_{R^i_e} + \mu_{R^j_e}.
\end{align*}

Therefore, since the tracking of $R_a$ is exponentially stable with convergence rate $r$, we have that

\begin{align*}
\|\dot p_{ij}(t)\| &\leq \sqrt{3}\left(\mu_{R^i_e}(t_0) + \mu_{R^j_e}(t_0)\right) e^{- r t} \nonumber\\
&\leq 2\sqrt{3}\pi \, e^{- r t},
\end{align*}

so that the relative position will be also stabilized exponentially fast. Indeed, we have that $p_{ij}(t) = p_{ij}(t_0) + \int_0^t  \dot p_{ij}(s) \mathrm{ds}$, so we can conclude that

\begin{align*}
\|p_{ij}(t) - p_{ij}(t_0)\| &= \left\|\int_0^t  \dot p_{ij}(s) \mathrm{ds}\right\| \leq \int_0^t  \|\dot p_{ij}(s)\| \mathrm{ds}\\ &\leq 2\sqrt{3}\pi \int_0^\infty  e^{- r s} \mathrm{ds} =  \frac{2\sqrt{3}\pi}{r}.
\end{align*}

§.§ Distance between two robots while they align with a non-totally known vector field

We finish this tutorial by showing how to compute the distance bound between two robots when they align with a time-varying vector field $m_d(t) \in \mathbb{R}^3$. However, the vector field is only known at time $t$ with $\left\|\dot m_d(t)\right\| \leq \omega_d \in\mathbb{R}^+$, so the robots cannot compensate their tracking with a feed-forward signal.

Given two robots $i, j$ with dynamics (<ref>), $||v|| =1$ and a desired attitude matrix $R_a \in SO(3)$. If $\Omega_{\{i,j\}}$ is given by (<ref>), $R_{i,j}(t_0) \in \mathcal{M}_{\delta^*}$ with $\delta^*$ sufficiently small and $\mu_{R_e^{\{i,j\}}}^*=\mu_{R_e}^*$ is designed as in Lemma <ref>, then

\begin{equation*}
\|p_{ij}(t) - p_{ij}(t_0)\| \leq\frac{2\sqrt{3}\mu_{R_e}^*}{k_w}.
\end{equation*}

Given $R_{ij} := R_j^\top R_i$, taking its time derivative yields

\begin{align*}
\dot R_{ij} = R_j^\top \dot R_i + \dot R^\top R_i
&= R_{ij}\Omega_i - \Omega_j R_{ij} \\
&= R_{ij} (\Omega_i - Ad_{R_{ij}^\top}(\Omega_j)) =  R_{ij}\Omega_{ij},
\end{align*}

and substituting $\Omega_{i}$ and $\Omega_{j}$ by the control law in (<ref>) leads us to

\begin{align} \label{eq: omega_ij}
\Omega_{ij} &= -k_w\left(\log(R_e^i) - Ad_{R_{ij}^\top}(\log(R_e^j))\right) \nonumber\\
&+ Ad_{R_i^\top R_a}(\Omega_a^k) - Ad_{R_{ij}^\top R_j^\top R_a}(\Omega_a^k)\nonumber\\
&= -k_w\left(\log(R_e^i) - Ad_{R_{ij}^\top}(\log(R_e^j))\right),
\end{align}

where $k_\omega^{\{i,j\}}=k_\omega$ because $k_\omega^{\{i,j\}}$ are designed as in Lemma <ref> and $\mu_{R_e^{\{i,j\}}}^*=\mu_{R_e}^*$. On the other hand, let us consider the Lyapunov function

\begin{align*}
V(R_e^{ij}) = \frac{1}{2}{\mu_{R_e^{ij}}}^2 = \frac{1}{2}\text{tr}\left(\log(R_e^{ij})^\top \log(R_e^{ij})\right),
\end{align*}

where $R_e^{ij} := (R_a^\top R_j)^\top R_a^\top R_i = R^{ij}$. Hence, following the same steps as in (<ref>), its time derivative leads us to

\begin{align} \label{eq: V_dot_small}
\dot V(R_e^{ij}) &= \text{tr}\left(\log(R_e^{ij})^\top \Omega_{ij}\right)\\
&= -k_w\text{tr}\left(\log(R_e^{ij})^\top \left(\log(R_e^i) - Ad_{R_{ij}^\top}\left(\log(R_e^j)\right)\right)\right)\nonumber.
\end{align}

Since $\mu_{R_e}^*$ is designed accordingly to Lemma <ref>, we have that if $\delta^*$ is sufficiently small then $\mu_{R_e^{\{i,j\}}}$ are too, i.e., $\mu_{R_e^{\{i,j\}}} \leq \mu_{R_e}^* \leq \delta^*$. Thus, substituting (<ref>) in (<ref>) yields

\begin{align*}
\Omega_{ij} &= -k_w\left(\log(R_e^i) - \sum^\infty_{m=0}\frac{(\mathrm{ad}_{\log(R_{ij})})^m}{m!}(\log(R_e^j))\right)\\
&= -k_w(\log(R_e^i) - \log(R_e^j)) + \text{H.O.T.},
\end{align*}

and considering the Baker-Campbell-Hausdorff formula [15] we have that $\log(R_e^{ij}) = \log(R_e^i) - \log(R_e^j) + \text{H.O.T.}$. Consequently, when $\delta^*$ is sufficiently small (<ref>) yields

\begin{align*}
\dot V(R_e^{ij}) &= -k_w \text{tr}\left(\log(R_e^{ij})^\top \log(R_e^{ij})\right) + \text{H.O.T.} \\ &= -k_w V(R_e^{ij}) + \text{H.O.T.},
\end{align*}

i.e., $\mu_{R_e^{ij}}(t) = \mu_{R_e^{ij}}(t_0) e^{-k_w t}$. Moreover, when $\delta^*$ is sufficiently small, we know from <cit.> that the higher order terms are not an issue concerning the exponential convergence. Therefore, following the same steps as in the proof of Lemma <ref> we have that

\begin{align*}
\|\dot p_{ij}(t)\| &\leq \sqrt{3} \mu_{R_e^{ij}} (t_0)e^{- k_w t} \\
&\leq \sqrt{3}\left(\mu_{R^i_e}(t_0) + \mu_{R^j_e}(t_0)\right)e^{- k_w t}\\
&\leq 2\sqrt{3}\mu_{R_e}^* e^{- k_w t}
\end{align*}

\begin{align*}
\|p_{ij}(t) - p_{ij}(t_0)\| &\leq 2\sqrt{3}\mu_{R_e}^* \int_0^\infty  e^{- k_w s} \mathrm{ds} \\
&= \frac{2\sqrt{3}\mu_{R_e}^*}{k_w}.
\end{align*}

Ideally, we would like to have (almost) global exponential stability but note that on the proof of Proposition <ref> we only show exponential stability when the distance between $R_i$ and $R_j$ is sufficiently small. Therefore, although this result is useful because it shows that $R_i = R_j$ is an expoentially stable equilibrium point, it does not enable us to bound the distance between robots when they start from distant points in $\text{SO}(3)$. Nonetheless, when we focus on 2D unicycle robots, Corollary <ref> shows that it is straightforward to demonstrate global exponential stability. Consequently, we will show in the simulations of the following section that $\mu_{R_e^{ij}}$ converges to zero exponentially fast.

Consider two 2D unicycle robots $i,j$ tracking a vector field $m_d(t)$ with the control law

$$\omega = -k_w\delta + \omega_d,$$

where $\omega$ is the angular heading velocity, $\delta$ is the heading angle error, and $\omega_d = \|\dot m_d(t)\|$. Define $\delta_{ij} := (\delta_j - \delta_i) \in (-2\pi,2\pi)$, so that $\delta_{ij} = 0$ means that robots $i,j$ share the same heading. Hence,

$$\dot \delta_{ij} = \omega_j - \omega_i = - k_w \delta_{ij},$$

i.e., $\delta_{ij}(t) = \delta_{ij}(0)e^{-k_w t}$, so $\delta_{ij}$ globally converges to zero exponentially fast.

For more details on Corollary <ref> we refer to [6]. In particular, it is shown that the distance

\begin{equation}
\|p_{ij}(t) - p_{ij}(0)\| \leq \frac{2\pi}{k_\omega}, \, \forall 1\leq i < j \leq N,
\end{equation}

for 2D unicycle robots on the same plane for all initial conditions of $\delta_{ij}$.

A unicycle robot moving at the constant speed $v = 0.5$ space unit/time unit along its $x$ body axis aligns with a time-varying attitude target $R_a(t) = [x_a \; y_a(t) \; z_a(t)]$, with $x_a,y_a,z_a \in S^2$ as in (<ref>). The time variation of $y(t)$ and $z(t)$ are given by the body-fixed angular velocity vector $\omega^k = [\pi,0,0]$ in radians/time unit, which is known. The control law for alignment is given by (<ref>), with $k_\omega = 1$. Note that, at $t=8$, $x_a$ instantly changes from $[-1,1,1]$ to $[1,0,0]$. The left plot shows the time evolution of the attitude error $\mu_{R_e}$. The right plot depicts the trajectory of the unicycle along with its attitude vectors $x$, $y$ and $z$ represented in red, green and blue, respectively.

A unicycle robot moving at the constant speed $v = 0.5$ space unit/time unit along its $x$ body axis aligns with a time-varying attitude target $R_a(t) = [x_a(t) \; y_a(t) \; z_a(t)]$, with $x_a,y_a,z_a \in S^2$ as in (<ref>). The time variation of $R_a(t)$ is given by the earth-fixed angular velocity vector $w(t) = R_a(t)^\top w^k + w^u$, where $\omega^k = [\pi,0,0]$ in radians/time unit and $\omega^u = [0,0,-w_d]$ with $w_d = \pi/14$ radians/time unit. The known variables are $w^k$ and $\|\omega^u\| = w_d$, so the control law for alignment is given by (<ref>), with $k_\omega = \sqrt{2}w_d/\mu_{R_e}^*$ and $\mu_{R_e}^* = \delta^* = 0.4$. The left plot shows the time evolution of the attitude error $\mu_{R_e}$ and $\delta$. The right plot depicts the trajectory of the unicycle along with its attitude vectors $x$, $y$ and $z$ represented in red, green and blue, respectively.

Two unicycle robots moving at the same constant speed $v = 0.5$ space unit/time unit along their $x$ body axis align with the time-varying attitude target $R_a(t) = [x_a(t) \; y_a(t) \; z_a(t)]$, with $x_a,y_a,z_a \in S^2$ as in (<ref>). The time variation of $R_a(t)$ is given by the earth-fixed angular velocity vector $w(t) = R_a(t)^\top w^k + w^u$, where $\omega^k = [\pi,0,0]$ in radians/time unit and $\omega^u = [0,0,-w_d]$ with $w_d = \pi/15$ radians/time unit. The known variables are $w^k$ and $\|\omega^u\| = w_d$, so the control law for alignment is given by (<ref>), with $k_\omega = \sqrt{2}w_d/\mu_{R_e}^*$ and $\mu_{R_e}^* = 0.5$ for both robots. The left plot shows the time evolution of the attitude error $\mu_{R_e}$ and $\|p_{ij}(t) - p_{ij}(t_0)\|$ for $t_0 = 0.5$. The right plot depicts the trajectory of the two unicycle robots along with their attitude vectors $x$, $y$ and $z$ represented in red, green and blue, respectively.

A robot swarm of $N=10$ unicycle robots moving at the same constant speed $v = 15$ space unit/time unit along their $x$ body axis seek the source of a scalar field $\sigma(p): \mathbb{R}^3 \rightarrow \mathbb{R}^3$ by aligning with the time-varying attitude target $R_a(t) = [m_d(t) \; y_a(t) \; z_a(t)]$, with $m_d,y_a,z_a \in S^2$ as in (<ref>), and where $m_d$ is an ascending direction calculated as in [6]. The time variation of $R_a(t)$ is given by the earth-fixed angular velocity vector $w(t) = R_a(t)^\top w^k + w^u$, where $\omega^k = [\pi,0,0]$ in radians/time unit and $w^u \in \mathbb{R}^3$ depends on the scalar field and the position of each agent, as it is deeply explained in [6]. Each robot of the swarm knows $w^k$ and $w_d = \pi/4$ radians/time unit, where $w_d \geq \|w^u\|$ when the centroid of the swarm $p_c(t)$ is in a certain set where $\|\nabla\sigma(p_c(t))\| > \epsilon$ for a given $\epsilon > 0$ <cit.>. The control law for alignment is given by (<ref>), with $k_\omega = \sqrt{2}w_d/\mu_{R_e}^*$ and $\mu_{R_e}^* = 0.4$ for all robots. The left plot shows the time evolution of the attitude error $\mu_{R_e}$ and $\|p_{i}(t) - p^*\|$, where $p^*$ is the source of the scalar field. The right plot depicts the $x$, $y$, and $z$ projections of the scalar field, along with the trajectory of all robots and their attitude vectors $x$, $y$, and $z$ represented in red, green, and blue, respectively.

§ SIMULATIONS

In this section, we validate numerically our results by conducting four different simulations. Firstly, in Figure <ref>, a unicycle robot aligns with an attitude target with a known time variation, which shows that Problem <ref> is solved by Theorem <ref>. In Figure <ref>, the same unicycle robot aligns with an attitude target with a non-totally known time variation, showing that Problem <ref> is solved by Theorem <ref>. Another robot is introduced in Figure <ref> to illustrate that the distance between two unicycle robots is exponentially stable and can be bounded as in Proposition <ref>. Finally, in Figure <ref>, we apply all the tools introduced in this work to the source seeking problem in 3D. In the last two simulations, we observe numerically that the exponential convergence of $\mu_{R_e^{ij}}$ to zero is global.

§ CONCLUSION

In this technical note, we introduced roboticists to the essential mathematical tools and techniques for analysing and implementing geometric controllers to align 3D robots with non-constant 3D vector fields. Additionally, we show how to apply this mathematical framework by exploring practical problems where the time variation of the vector field is unknown, and by presenting some technical results with detailed calculations. In future revisions of this work, we aim to extend the proof of Proposition <ref> for global exponential stability and to explore different practical problems that will enable us to introduce more advanced analysis tools.

§.§ The calculation of geodesic distance in $S^2$

Let $x,y \in S^2$ and an arc $\gamma(t) = (t, h(t)) \subset S^2$, where $h(t): \mathbb{R}^+ \rightarrow \mathbb{R}^2$. Consider that this arc connects $x$ with $y$, so $x = \gamma(t_x) = (t_x, h(t_x))$ and $y = \gamma(t_y) = (t_y, h(t_y))$ for a given $t_y \geq t_x$. Hence, given $\langle x , y \rangle = x \cdot y = x^\top y$ as a valid Riemannian metric for $S^2$, the definition of Riemannian length in (<ref>) yields

\ell(\gamma)  = \int \sqrt{\langle \dot\gamma(t) , \dot\gamma(t) \rangle_{\gamma(t)}} \mathrm{d}t = \int \|\dot\gamma(t)\| \mathrm{d}t,

where $\dot \gamma(t) = (1,\dot h(t))$, an so $\|\dot \gamma(t)\| = \sqrt{1 + \|\dot h(t)\|^2}$.

Note now that $t^2 + h(t)^\top h(t) = \langle \dot\gamma(t) , \dot\gamma(t) \rangle = 1$, so differentiating both sides with respect to $t$ yields

2t + 2h(t)^\top \dot h(t) = 0 \quad \Leftrightarrow \quad t + h(t)^\top \dot h(t) = 0.

Considering this identity in the Cauchy Schwartz inequality leads us to

\|\dot h(t)\|^2 \|h(t)\|^2 \geq \|h(t)^\top \dot h(t)\|^2 = t^2,

\|\dot h(t)\|^2 \geq \frac{t^2}{\|h(t)\|^2} \quad \Leftrightarrow \quad 1 + \|\dot h(t)\|^2 \geq \frac{t^2}{\|h(t)\|^2} + 1.

But $\|h(t)\|^2 = 1 - t^2$, thus

\sqrt{1 + \|\dot h(t)\|^2} \geq \sqrt{1 + \frac{t^2}{1 - t^2}} = \sqrt{\frac{1}{1 - t^2}},

and so

\begin{align*}
\ell(\gamma) = \int_{t_x}^{t_y} \sqrt{1 + \|\dot h(t)\|^2} \mathrm{d}t
&\geq \int_{t_x}^{t_y} \sqrt{\frac{1}{1 - t^2}} \mathrm{d}t
\end{align*}

This last integral represents the minimum length of an arc connecting $x$ to $y$. In fact, this is the definition of length-minimizing geodesic $\gamma_g$. Therefore, if we assume without loss of generality that $x = (a,s,0)$ and $y = (1,0,0)$, then we have

\begin{align*}
\ell(\gamma_g) = \int_{a}^{1} \sqrt{\frac{1}{1 - t^2}} \mathrm{d}t = -\arccos(t) \Big|_{a}^{1} &= \arccos(a) \\
&= \arccos(x^\top y).
\end{align*}

[1]
K. M. Lynch and F. C. Park, Modern Robotics: Mechanics, Planning, and Control, 1st ed.1em plus 0.5em minus 0.4emUSA: Cambridge University Press, 2017.

[2]
J. Solà, J. Deray, and D. Atchuthan, “A micro lie theory for state estimation in robotics,” 2021.

[3]
F. Bullo and A. D. Lewis, Geometric Control of Mechanical Systems, ser. Texts in Applied Mathematics.1em plus 0.5em minus 0.4emNew York-Heidelberg-Berlin: Springer Verlag, 2004, vol. 49.

[4]
H. A. Hashim, “Special orthogonal group so(3), euler angles, angle-axis, rodriguez vector and unit-quaternion: Overview, mapping and challenges,” 2021.

[5]
W. Yao, H. G. de Marina, B. Lin, and M. Cao, “Singularity-free guiding vector field for robot navigation,” IEEE Transactions on Robotics, vol. 37, no. 4, pp. 1206–1221, 2021.

[6]
A. Acuaviva, J. Bautista, W. Yao, J. Jimenez, and H. G. de Marina, “Resilient source seeking with robot swarms,” arXiv preprint, 2023.

[7]
Differentiable manifold. <https://en.wikipedia.org/wiki/Differentiable_manifold>.

[8]
M. Do Carmo, Riemannian Geometry, ser. Mathematics (Boston, Mass.).1em plus 0.5em minus 0.4emBirkhäuser, 1992.

[9]
Jesús B. V. (2023) Source Seeking 3D Simulator, GitHub repository. <https://github.com/jesusBV20/source_seeking_3D>.

[10]
R. Vemulapalli and R. Chellappa, “Rolling rotations for recognizing human actions from 3d skeletal data,” pp. 4471–4479, 2016.

[11]
J. H. Gallier, “Lecture notes in advanced geometric methods in computer science,” January 2023.

[12]
F. Bullo and R. M. Murray, Proportional Derivative (PD) Control on the Euclidean Group, ene 1995.

[13]
Derivation of the rodrigues' rotation formula. <https://en.wikipedia.org/wiki/Rodrigues'_rotation_formula#Derivation>.

[14]
D. Huynh, “EnglishMetrics for 3d rotations: Comparison and analysis,” EnglishJournal of Mathematical Imaging and Vision, vol. 35, no. 2, pp. 155–164, 2009.

[15]
N. Thanh and P. Qi, “Baker-campbell-hausdorff theorem,” 2023.

[16]
H. Khalil, Nonlinear Control, ser. Always Learning.1em plus 0.5em minus 0.4emPearson, 2014.

[17]
A. A. Huertos, “Control of multi-agent systems and their applications,” 2022. [Online]. Available: <https://www.swarmsystemslab.eu/team>


Further author information: (Send correspondence to Ali Vosoughi)
Ali Vosoughi: E-mail<EMAIL_ADDRESS>
# Large-Scale Extended Granger Causality for Classification of Marijuana Users
From Functional MRI

M. Ali Vosoughi  Department of Electrical and Computer Engineering, University
of Rochester, NY, USA  Axel Wismüller  Department of Electrical and Computer
Engineering, University of Rochester, NY, USA  Department of Imaging Sciences,
University of Rochester, NY, USA  Department of Biomedical Engineering,
University of Rochester, NY, USA  Faculty of Medicine and Institute of
Clinical Radiology, Ludwig Maximilian University, Munich, Germany

###### Abstract

It has been shown in the literature that marijuana use is associated with
changes in brain network connectivity. We propose large-scale Extended Granger
Causality (lsXGC) and investigate whether it can capture such changes using
resting-state fMRI. This method combines dimension reduction with source time-
series augmentation and uses predictive time-series modeling for estimating
directed causal relationships among fMRI time-series. It is a multivariate
approach, since it is capable of identifying the interdependence of time-
series in the presence of all other time-series of the underlying dynamic
system. Here, we investigate whether this model can serve as a biomarker for
classifying marijuana users from typical controls using 126 adult subjects
with a childhood diagnosis of ADHD form the Addiction Connectome Preprocessed
Initiative (ACPI) database. We use brain connections estimated by lsXGC as
features for classification. After feature extraction, we perform feature
selection by Kendall’s-tau rank correlation coefficient followed by
classification using a support vector machine. As a reference method, we
compare our results with cross-correlation, which is typically used in the
literature as a standard measure of functional connectivity. Within a cross-
validation scheme of 100 different training/test (90%/10%) data splits, we
obtain a mean accuracy range of [0.714, 0.985] and a mean Area Under the
receiver operating characteristic Curve (AUC) range of [0.779, 0.999] across
all tested numbers of features for lsXGC, which is significantly better than
results obtained with cross-correlation, namely mean accuracy of [0.728,
0.912] and mean AUC of [0.825, 0.969]. Our results suggest the applicability
of lsXGC as a potential biomarker for marijuana use.

###### keywords:

Resting-state fMRI, Large-Scale Extended Granger Causality, functional
connectivity, machine learning, support vector machine, Attention Deficit
Hyperactivity Disorder (ADHD), effect of marijuana use

##  1 INTRODUCTION

Marijuana has been used as an illegal drug for decades. Despite many concerns
about the adverse effects of repeated marijuana abuse on brain functionalities
and cognition, considerably few studies have been conducted to investigate the
effects of repeated marijuana use on the brain [1]. Currently, the diagnosis
of cannabis use disorder is based on the clinical evaluations of the symptoms
and behaviors. However, more objective biomarkers are of interest. To this
end, more recently, studies have investigated, if information can be extracted
non-invasively from brain activity data. Although these studies have shown
promising results, there is still scope for improvement, especially with
regards to using more meaningful connectivity analysis approaches [2].

Some evidence has demonstrated that frequent marijuana use affects the
connectivity of the brain [2, 1]. Biomarkers from resting-state functional MRI
(rs-fMRI) for frequent marijuana use can be derived using Multi-Voxel Pattern
Analysis (MVPA) techniques [3, 2]. MVPA is a machine-learning framework that
extracts differences in patterns of brain connectivity to discriminate between
connectivity profiles of individuals with neurological diseases/disorders and
healthy individuals. Most MVPA studies commonly use cross-correlation to
construct a functional connectivity profile. For example, one such study has
obtained an accuracy of 0.69 on the slow frequency bands (0.01-0.1 Hz) [2].
These results demonstrate that meaningful information can be learned from fMRI
data. However, correlation does not provide a measure of directed
connectivity. Hence, there may be more relevant information in the fMRI data
that is not being captured by cross-correlation.

Various methods have been proposed to obtain directional relationships in
multivariate time-series data, e.g., transfer entropy [4] and mutual
information [5]. However, as the multivariate problem’s dimensions increase,
computation of the density function becomes computationally expensive [6, 7].
Under the Gaussian assumption, transfer entropy is equivalent to Granger
causality [8]. However, the computation of multivariate Granger causality for
short time series in large-scale problems is challenging [9, 10]. We introduce
large-scale Extended Granger Causality (lsXGC) as a novel method that combines
dimension reduction with source time-series augmentation and uses predictive
time-series modeling for estimating directed causal relationships among fMRI
time-series [11]. In this work, we investigate, if changes in directed
connectivity manifest the frequent marijuana use among adults with a childhood
diagnosis of ADHD and if such directed measures are better able to
discriminate between frequent marijuana users and typical controls. To this
end, we apply lsXGC in the MVPA framework for estimating a measure of directed
causal interdependence between fMRI time-series.

This work is embedded in our group’s endeavor to expedite artificial
intelligence in biomedical imaging by means of advanced pattern recognition
and machine learning methods for computational radiology and radiomics, e.g. [
12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49,
50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69
].

##  2 DATA

This data respiratory of the Multimodal Treatment Study of ADHD (MTA) contains
a longitudinal study of participants for fourteen years pursued at multiple
sites in North America. The study is one of the largest conducted studies for
ADHD diagnosis and treatment, with primary results being published in 1999.
The Addiction Connectome Preprocessed Initiative (ACPI) data respiratory
contains datasets from a subset of MTA participants at six sites with and
without childhood ADHD, who were studied as part of a follow-up multimodal
brain imaging examination. The principal aim of the MTA study was to
investigate the effect of cannabis use among adults with a childhood diagnosis
of ADHD. The study was a 2 x 2 design of those with and without childhood ADHD
and those who did or did not regularly use cannabis [70].

The scan parameters of resting-state fMRI the ACPI dataset can be found on the
website of the project [70]. Preprocessed data from the MTA 1 dataset [70] of
the ACPI database was used. The dataset included 126 subjects, 101 males and
25 females, with ages between 21-27 years, 86 diagnosed with ADHD (68%), and
62 of them regularly used marijuana (49%). Preprocessing of the raw 4D rs-fMRI
data had been made using a Configurable Pipeline for the Analysis of
Connectomes (C-PAC) [71] and the Advanced Normalization Tools (ANTs) pipeline
and consisted in the removal of the first five fMRI volumes, anatomical
registration, tissue segmentation, functional registration in the Montreal
Neurological Institute (MNI) space, functional masking, temporal bandpass
filtering (0.01 - 0.1 Hz), motion correction, spatial smoothing, and various
nuisance corrections [72]. The MODL parcellation has been used as a pre-
computed atlas for regions definition and built using a form of online
dictionary learning [2].

##  3 METHODS

###  3.1 Large-scale Extended Granger Causality (lsXGC)

Large-scale Extended Granger Causality (lsXGC) has been developed based on 1)
the principle of original Granger causality, which quantifies the causal
influence of time-series $\mathbb{x_{s}}$ on time-series $\mathbb{x_{t}}$ by
quantifying the amount of improvement in the prediction of $\mathbb{x_{t}}$ in
presence of $\mathbb{x_{s}}$. 2) the idea of dimension reduction, which
resolves the problem of the tackling a under-determined system, which is
frequently faced in fMRI analysis, since the number of acquired temporal
samples usually is not sufficient for estimating the model parameters [9, 10].

Consider the ensemble of time-series $\mathcal{X}\in\mathbb{R}^{N\times T}$,
where $N$ is the number of time-series (Regions Of Interest – ROIs) and $T$
the number of temporal samples. Let
$\mathcal{X}=(\mathbb{x_{1}},\mathbb{x_{2}},\dots,\mathbb{x_{N}})^{\mathsf{T}}$
be the whole multidimensional system and $x_{i}\in\mathbb{R}^{1\times T}$ a
single time-series with $i=1,2,\dots,N$, where
$\mathbb{x_{i}}=(x_{i}(1),x_{i}(2),\dots,x_{i}(T))$. In order to overcome the
under-determined problem, first $\mathcal{X}$ will be decomposed into its
first $p$ high-variance principal components
$\mathcal{Z}\in\mathbb{R}^{p\times T}$ using Principal Component Analysis
(PCA), i.e.,

$\mathcal{Z}=W\mathcal{X},$ (1)

where $W\in\mathbb{R}^{p\times N}$ represents the PCA coefficient matrix.
Subsequently, the dimension-reduced time-series ensemble $\mathcal{Z}$ is
augmented by one original time-series $\mathbb{x_{s}}$ yielding a dimension-
reduced augmented time-series ensemble $\mathcal{Y}\in\mathbb{R}^{(p+1)\times
T}$ for estimating the influence of $\mathbb{x_{s}}$ on all other time-series.

Following this, we locally predict $\mathcal{X}$ at each time sample $t$, i.e.
$\mathcal{X}(t)\in\mathbb{R}^{N\times 1}$ by calculating an estimate
$\hat{\mathcal{X}}_{\mathbb{x_{s}}}(t)$. To this end, we fit an affine model
based on a vector of $m$ vector of m time samples of
$\mathcal{Y}(\tau)\in\mathbb{R}^{(p+1)\times 1}$($\tau=t-1,t-2,\dots,t-m$),
which is $\mathbb{y}(t)\in\mathbb{R}^{m.(p+1)\times 1}$, and a parameter
matrix $\mathcal{A}\in\mathbb{R}^{N\times m.(p+1)}$ and a constant bias vector
$\mathbb{b}\in\mathbb{R}^{N\times 1}$,

$\hat{\mathcal{X}}_{\mathbb{x_{s}}}(t)=\mathcal{A}\mathbb{y}(t)+\mathbb{b},~{}~{}t=m+1,m+2,\dots,T.$
(2)

Now $\hat{\mathcal{X}}_{\setminus{\mathbb{x_{s}}}}(t)$, which is the
prediction of $\mathcal{X}(t)$ without the information of $\mathbb{x_{s}}$,
will be estimated. The estimation processes is identical to the previous one,
with the only difference being that we have to remove the augmented time-
series $\mathbb{x_{s}}$ and its corresponding column in the PCA coefficient
matrix $W$.

The computation of a lsXGC index is based on comparing the variance of the
prediction errors obtained with and without consideration of $\mathbb{x_{s}}$.
The lsXGC index $f_{\mathbb{x_{s}}\xrightarrow{}\mathbb{x_{t}}}$ , which
indicates the influence of $\mathbb{x_{s}}$ on $\mathbb{x_{t}}$, can be
calculated by the following equation:

$f_{\mathbb{x_{s}}\xrightarrow{}\mathbb{x_{t}}}=\log{\frac{\mathrm{var}(e_{s})}{\mathrm{var}(e_{\setminus
s})}},$ (3)

where $e_{\setminus s}$ is the error in predicting $\mathbb{x_{t}}$ when
$\mathbb{x_{s}}$ was not considered, and $e_{s}$ is the error, when
$\mathbb{x_{s}}$ was used. In this study, we set $p=8$ and $m=4$.

###  3.2 Multi-voxel pattern analysis

In this study, brain connections served as features for classification and
were estimated by two methods, namely lsXGC and cross-correlation. Before
using high-dimensional connectivity feature vectors as input to a classifier,
feature selection was carried out to reduce the dimension of input features.

####  3.2.1 Feature selection

To reduce the number of features, feature selection was performed on each
training data set with k-fold cross-validation using Kendall’s Tau rank
correlation coefficient [73] and $10\%-90\%$ of test-to-train split ratio.
This approach quantifies the relevance of each feature to the task of
classification and assigns ranks by testing for independence between different
classes for each feature [73].

####  3.2.2 Classification

For the 100 iteration cross-validation scheme for classification, the data set
was divided into two groups: a training data set ($90\%$) and a test data set
($10\%$) in a way that the percentage of samples for each class was preserved.
Also, this was repeated with different numbers of features ranging from 5 to
175. A Support Vector Machine (SVM) [74] was used for classification between
control subjects and marijuana users. All procedures were implemented using
MATLAB 9.8 (MathWorks Inc., Natick, MA, 2020a) and Python 3.8.

##  4 RESULTS

Mean connectivity matrices that were extracted using lsXGC and cross-
correlation, are shown in Fig. 1 for both typical control and frequent
marijuana user cohorts. Different patterns are visible to the naked eye for
both methods. In the following, we quantitatively investigate the difference
between connectivity patterns of the two subject cohorts using an MVPA
approach.

Figure 1: Mean connectivity matrices: top left: mean connectivity matrix of
typical control subjects using lsXGC, top right: mean connectivity matrix of
frequent marijuana users using lsXGC, bottom left: mean connectivity matrix of
typical controls using cross-correlation, bottom right: mean connectivity
matrix of frequent marijuana users using cross-correlation. Note that the
different methods appear to extract different connectivity features and that
they appear to be slight differences in connectivity patterns between the
typical controls and the frequent marijuana users.

Classification results were evaluated using the Area Under the Receiver
Operator Characteristic Curve (AUC) and accuracy. An AUC = 1 indicates a
perfect classification, AUC = 0.5 indicates random class assignment.

(a) Mean AUC

(b) Mean accuracy

Figure 2: Plots comparing the performance of cross-correlation and the
proposed large-scale extended Granger causality (lsXGC). The shaded areas
represent the $95\%$ confidence interval. It is clearly demonstrated that
lsXGC outperforms cross-correlation for most numbers of selected features.

In this study, based on preliminary analyses, we chose 8 as the number of the
retained components of PCA in the lsXGC algorithm, and model order of 4 for
the multivariate vector autoregression function. From the plots of AUC and
accuracy results in Fig. 2, we can clearly see that lsXGC outperforms cross-
correlation for varying numbers of features. Across the wide range of
investigated numbers of features, the performance of lsXGC is consistently
good with its mean AUC within [0.779, 0.999] and its mean accuracy within
[0.714, 0.985]. On the other hand, cross-correlation performs quite poorly
compared to lsXGC with its mean AUC within [0.825, 0.969] and its mean
accuracy within [0.728, 0.912].

##  5 CONCLUSIONS

In this study, we use a recently developed method for brain connectivity
analysis, large-scale Extended Granger Causality (lsXGC), and apply it to
subjects of the ACPI data repository in order to classify individuals history
of the childhood ADHD diagnosis with frequent marijuana use from typical
controls by estimating a measure of directed causal interaction between
regional brain activities measured in resting-state fMRI. After constructing
connectivity matrices as characterizing features for brain network analysis,
we use Kendall’s tau rank correlation coefficient for feature selection and a
support vector machine for classification. We demonstrate that lsXGC compares
favorably to conventional cross-correlation analysis, as indicated by the
significantly higher AUC and accuracy values. The potential of lsXGC as an
effective biomarker for identifying frequent marijuana use is yet to be
validated in prospective clinical trials. However, our results suggest that
our approach outperforms the current clinical standard, namely cross-
correlation, at uncovering meaningful information from functional MRI data.

###### Acknowledgements.

This research was funded by Ernest J. Del Monte Institute for Neuroscience
Award from the Harry T. Mangurian Jr. Foundation. This work was conducted as a
Practice Quality Improvement (PQI) project related to American Board of
Radiology (ABR) Maintenance of Certificate (MOC) for Prof. Dr. Axel Wismüller.
This work is not being and has not been submitted for publication or
presentation elsewhere.

## References

  * [1] Solowij, N., [Cannabis and cognitive functioning ], Cambridge University Press (2006).
  * [2] Dadi, K., Rahim, M., Abraham, A., Chyzhyk, D., Milham, M., Thirion, B., Varoquaux, G., Initiative, A. D. N., et al., “Benchmarking functional connectome-based predictive models for resting-state fMRI,” Neuroimage 192, 115–134 (2019).
  * [3] Norman, K. A., Polyn, S. M., Detre, G. J., and Haxby, J. V., “Beyond mind-reading: multi-voxel pattern analysis of fMRI data,” Trends in cognitive sciences 10(9), 424–430 (2006).
  * [4] Schreiber, T., “Measuring information transfer,” Physical review letters 85(2), 461 (2000).
  * [5] Kraskov, A., Stögbauer, H., and Grassberger, P., “Estimating mutual information,” Physical review E 69(6), 066138 (2004).
  * [6] Mozaffari, M. and Yilmaz, Y., “Online multivariate anomaly detection and localization for high-dimensional settings,” arXiv preprint arXiv:1905.07107 (2019).
  * [7] Mozaffari, M. and Yilmaz, Y., “Online anomaly detection in multivariate settings,” in [2019 IEEE 29th International Workshop on Machine Learning for Signal Processing (MLSP) ], 1–6, IEEE (2019).
  * [8] Barnett, L., Barrett, A. B., and Seth, A. K., “Granger causality and transfer entropy are equivalent for Gaussian variables,” Physical review letters 103(23), 238701 (2009).
  * [9] Vosoughi, M. A. and Wismüller, A., “Large-scale kernelized Granger causality to infer topology of directed graphs with applications to brain networks,” arXiv preprint arXiv:2011.08261 (2020).
  * [10] Wismüller, A., DSouza, A. M., Abidin, A. Z., and Vosoughi, M. A., “Large-scale nonlinear Granger causality: A data-driven, multivariate approach to recovering directed networks from short time-series data,” arXiv preprint arXiv:2009.04681 (2020).
  * [11] Wismüller, A., Foxe, J. J., Geha, P., and Saboksayr, S. S., “Large-scale extended Granger causality (lsxgc) for classification of autism spectrum disorder from resting-state functional MRI,” in [Medical Imaging 2020: Computer-Aided Diagnosis ], 11314, 113141Y, International Society for Optics and Photonics (2020).
  * [12] Nattkemper, T. W. and Wismüller, A., “Tumor feature visualization with unsupervised learning,” Medical Image Analysis 9(4), 344–351 (2005).
  * [13] Bunte, K., Hammer, B., Wismüller, A., and Biehl, M., “Adaptive local dissimilarity measures for discriminative dimension reduction of labeled data,” Neurocomputing 73(7-9), 1074–1092 (2010).
  * [14] Wismueller, A., Vietze, F., and Dersch, D. R., “Segmentation with neural networks,” in [Handbook of medical imaging ], 107–126, Academic Press, Inc. (2000).
  * [15] Leinsinger, G., Schlossbauer, T., Scherr, M., Lange, O., Reiser, M., and Wismüller, A., “Cluster analysis of signal-intensity time course in dynamic breast MRI: does unsupervised vector quantization help to evaluate small mammographic lesions?,” European radiology 16(5), 1138–1146 (2006).
  * [16] Wismüller, A., Vietze, F., Behrends, J., Meyer-Baese, A., Reiser, M., and Ritter, H., “Fully automated biomedical image segmentation by self-organized model adaptation,” Neural Networks 17(8-9), 1327–1344 (2004).
  * [17] Hoole, P., Wismüller, A., Leinsinger, G., Kroos, C., Geumann, A., and Inoue, M., “Analysis of tongue configuration in multi-speaker, multi-volume MRI data,” (2000).
  * [18] Wismüller, A., “Exploratory morphogenesis (xom): a novel computational framework for self-organization,” Ph. D. thesis, Technical University of Munich, Department of Electrical and Computer Engineering (2006).
  * [19] Wismüller, A., Dersch, D. R., Lipinski, B., Hahn, K., and Auer, D., “A neural network approach to functional MRI pattern analysis—clustering of time-series by hierarchical vector quantization,” in [International Conference on Artificial Neural Networks ], 857–862, Springer (1998).
  * [20] Wismüller, A., Vietze, F., Dersch, D. R., Behrends, J., Hahn, K., and Ritter, H., “The deformable feature map-a novel neurocomputing algorithm for adaptive plasticity in pattern analysis,” Neurocomputing 48(1-4), 107–139 (2002).
  * [21] Behrends, J., Hoole, P., Leinsinger, G. L., Tillmann, H. G., Hahn, K., Reiser, M., and Wismüller, A., “A segmentation and analysis method for MRI data of the human vocal tract,” in [Bildverarbeitung für die Medizin 2003 ], 186–190, Springer (2003).
  * [22] Wismuller, A., “Neural network computation in biomedical research: chances for conceptual cross-fertilization,” Theory in Biosciences (1997).
  * [23] Bunte, K., Hammer, B., Villmann, T., Biehl, M., and Wismüller, A., “Exploratory observation machine (xom) with kullback-leibler divergence for dimensionality reduction and visualization.,” in [ESANN ], 10, 87–92 (2010).
  * [24] Wismüller, A., Vietze, F., Dersch, D. R., Hahn, K., and Ritter, H., “The deformable feature map—adaptive plasticity for function approximation,” in [International Conference on Artificial Neural Networks ], 123–128, Springer (1998).
  * [25] Wismüller, A., “The exploration machine–a novel method for data visualization,” in [International Workshop on Self-Organizing Maps ], 344–352, Springer (2009).
  * [26] Wismüller, A., “Method, data processing device and computer program product for processing data,” (July 28 2009).  US Patent 7,567,889.
  * [27] Huber, M. B., Nagarajan, M., Leinsinger, G., Ray, L. A., and Wismüller, A., “Classification of interstitial lung disease patterns with topological texture features,” in [Medical Imaging 2010: Computer-Aided Diagnosis ], 7624, 762410, International Society for Optics and Photonics (2010).
  * [28] Wismüller, A., “The exploration machine: a novel method for analyzing high-dimensional data in computer-aided diagnosis,” in [Medical Imaging 2009: Computer-Aided Diagnosis ], 7260, 72600G, International Society for Optics and Photonics (2009).
  * [29] Bunte, K., Hammer, B., Villmann, T., Biehl, M., and Wismüller, A., “Neighbor embedding xom for dimension reduction and visualization,” Neurocomputing 74(9), 1340–1350 (2011).
  * [30] Meyer-Bäse, A., Lange, O., Wismüller, A., and Ritter, H., “Model-free functional MRI analysis using topographic independent component analysis,” International journal of neural systems 14(04), 217–228 (2004).
  * [31] Wismüller, A., “A computational framework for nonlinear dimensionality reduction and clustering,” in [International Workshop on Self-Organizing Maps ], 334–343, Springer (2009).
  * [32] Meyer-Base, A., Auer, D., and Wismueller, A., “Topographic independent component analysis for fMRI signal detection,” in [Proceedings of the International Joint Conference on Neural Networks, 2003. ], 1, 601–605, IEEE (2003).
  * [33] Meyer-Baese, A., Schlossbauer, T., Lange, O., and Wismüller, A., “Small lesions evaluation based on unsupervised cluster analysis of signal-intensity time courses in dynamic breast MRI,” International journal of biomedical imaging 2009 (2009).
  * [34] Wismueller, A., Lange, O., Auer, D., and Leinsinger, G., “Model-free functional MRI analysis for detecting low-frequency functional connectivity in the human brain,” in [Medical Imaging 2010: Computer-Aided Diagnosis ], 7624, 76241M, International Society for Optics and Photonics (2010).
  * [35] Meyer-Bäse, A., Saalbach, A., Lange, O., and Wismüller, A., “Unsupervised clustering of fMRI and MRI time series,” Biomedical Signal Processing and Control 2(4), 295–310 (2007).
  * [36] Huber, M. B., Nagarajan, M. B., Leinsinger, G., Eibel, R., Ray, L. A., and Wismüller, A., “Performance of topological texture features to classify fibrotic interstitial lung disease patterns,” Medical Physics 38(4), 2035–2044 (2011).
  * [37] Wismüller, A., Verleysen, M., Aupetit, M., and Lee, J. A., “Recent advances in nonlinear dimensionality reduction, manifold and topological learning.,” in [ESANN ], (2010).
  * [38] Meyer-Baese, A., Lange, O., Wismüller, A., and Hurdal, M. K., “Analysis of dynamic susceptibility contrast MRI time series based on unsupervised clustering methods,” IEEE Transactions on Information Technology in Biomedicine 11(5), 563–573 (2007).
  * [39] Wismueller, A., Behrends, J., Hoole, P., Leinsinger, G. L., Reiser, M. F., and Westesson, P.-L., “Human vocal tract analysis by in vivo 3d MRI during phonation: a complete system for imaging, quantitative modeling, and speech synthesis,” in [International Conference on Medical Image Computing and Computer-Assisted Intervention ], 306–312, Springer (2008).
  * [40] Wismüller, A., “Method and device for representing multichannel image data,” (Nov. 17 2015).  US Patent 9,189,846.
  * [41] Huber, M. B., Bunte, K., Nagarajan, M. B., Biehl, M., Ray, L. A., and Wismüller, A., “Texture feature ranking with relevance learning to classify interstitial lung disease patterns,” Artificial intelligence in medicine 56(2), 91–97 (2012).
  * [42] Wismuller, A., Meyer-Baese, A., Lange, O., Reiser, M. F., and Leinsinger, G., “Cluster analysis of dynamic cerebral contrast-enhanced perfusion MRI time-series,” IEEE transactions on medical imaging 25(1), 62–73 (2005).
  * [43] Twellmann, T., Saalbach, A., Muller, C., Nattkemper, T. W., and Wismuller, A., “Detection of suspicious lesions in dynamic contrast enhanced MRI data,” in [The 26th Annual International Conference of the IEEE Engineering in Medicine and Biology Society ], 1, 454–457, IEEE (2004).
  * [44] Otto, T. D., Meyer-Baese, A., Hurdal, M., Sumners, D., Auer, D., and Wismuller, A., “Model-free functional MRI analysis using cluster-based methods,” in [Intelligent Computing: Theory and Applications ], 5103, 17–24, International Society for Optics and Photonics (2003).
  * [45] Varini, C., Nattkemper, T. W., Degenhard, A., and Wismuller, A., “Breast MRI data analysis by lle,” in [2004 IEEE International Joint Conference on Neural Networks (IEEE Cat. No. 04CH37541) ], 3, 2449–2454, IEEE (2004).
  * [46] Huber, M. B., Lancianese, S. L., Nagarajan, M. B., Ikpot, I. Z., Lerner, A. L., and Wismuller, A., “Prediction of biomechanical properties of trabecular bone in mr images with geometric features and support vector regression,” IEEE Transactions on Biomedical Engineering 58(6), 1820–1826 (2011).
  * [47] Meyer-Base, A., Pilyugin, S. S., and Wismuller, A., “Stability analysis of a self-organizing neural network with feedforward and feedback dynamics,” in [2004 IEEE International Joint Conference on Neural Networks (IEEE Cat. No. 04CH37541) ], 2, 1505–1509, IEEE (2004).
  * [48] Meyer-Baese, A., Lange, O., Schlossbauer, T., and Wismuller, A., “Computer-aided diagnosis and visualization based on clustering and independent component analysis for breast MRI,” in [2008 15th IEEE International Conference on Image Processing ], 3000–3003, IEEE (2008).
  * [49] Wismüller, A., Meyer-Bäse, A., Lange, O., Schlossbauer, T., Kallergi, M., Reiser, M., and Leinsinger, G., “Segmentation and classification of dynamic breast magnetic resonance image data,” Journal of Electronic Imaging 15(1), 013020 (2006).
  * [50] Bhole, C., Pal, C., Rim, D., and Wismüller, A., “3d segmentation of abdominal ct imagery with graphical models, conditional random fields and learning,” Machine vision and applications 25(2), 301–325 (2014).
  * [51] Nagarajan, M. B., Coan, P., Huber, M. B., Diemoz, P. C., Glaser, C., and Wismüller, A., “Computer-aided diagnosis in phase contrast imaging x-ray computed tomography for quantitative characterization of ex vivo human patellar cartilage,” IEEE Transactions on Biomedical Engineering 60(10), 2896–2903 (2013).
  * [52] Wismüller, A., Meyer-Bäse, A., Lange, O., Auer, D., Reiser, M. F., and Sumners, D., “Model-free functional MRI analysis based on unsupervised clustering,” Journal of Biomedical Informatics 37(1), 10–18 (2004).
  * [53] Meyer-Baese, A., Wismueller, A., Lange, O., and Leinsinger, G., “Computer-aided diagnosis in breast MRI based on unsupervised clustering techniques,” in [Intelligent Computing: Theory and Applications II ], 5421, 29–37, International Society for Optics and Photonics (2004).
  * [54] Nagarajan, M. B., Coan, P., Huber, M. B., Diemoz, P. C., Glaser, C., and Wismüller, A., “Computer-aided diagnosis for phase-contrast x-ray computed tomography: quantitative characterization of human patellar cartilage with high-dimensional geometric features,” Journal of digital imaging 27(1), 98–107 (2014).
  * [55] Nagarajan, M. B., Huber, M. B., Schlossbauer, T., Leinsinger, G., Krol, A., and Wismüller, A., “Classification of small lesions on dynamic breast MRI: Integrating dimension reduction and out-of-sample extension into cadx methodology,” Artificial intelligence in medicine 60(1), 65–77 (2014).
  * [56] Yang, C.-C., Nagarajan, M. B., Huber, M. B., Carballido-Gamio, J., Bauer, J. S., Baum, T. H., Eckstein, F., Lochmüller, E.-M., Majumdar, S., Link, T. M., et al., “Improving bone strength prediction in human proximal femur specimens through geometrical characterization of trabecular bone microarchitecture and support vector regression,” Journal of electronic imaging 23(1), 013013 (2014).
  * [57] Wismüller, A., Nagarajan, M. B., Witte, H., Pester, B., and Leistritz, L., “Pair-wise clustering of large scale Granger causality index matrices for revealing communities,” in [Medical Imaging 2014: Biomedical Applications in Molecular, Structural, and Functional Imaging ], 9038, 90381R, International Society for Optics and Photonics (2014).
  * [58] Wismüller, A., Wang, X., DSouza, A. M., and Nagarajan, M. B., “A framework for exploring non-linear functional connectivity and causality in the human brain: mutual connectivity analysis (mca) of resting-state functional MRI with convergent cross-mapping and non-metric clustering,” arXiv preprint arXiv:1407.3809 (2014).
  * [59] Schmidt, C., Pester, B., Nagarajan, M., Witte, H., Leistritz, L., and Wismueller, A., “Impact of multivariate Granger causality analyses with embedded dimension reduction on network modules,” in [2014 36th Annual International Conference of the IEEE Engineering in Medicine and Biology Society ], 2797–2800, IEEE (2014).
  * [60] Wismüller, A., Abidin, A. Z., D’Souza, A. M., Wang, X., Hobbs, S. K., Leistritz, L., and Nagarajan, M. B., “Nonlinear functional connectivity network recovery in the human brain with mutual connectivity analysis (MCA): convergent cross-mapping and non-metric clustering,” in [Medical Imaging 2015: Biomedical Applications in Molecular, Structural, and Functional Imaging ], 9417, 94170M, International Society for Optics and Photonics (2015).
  * [61] Wismüller, A., Abidin, A. Z., DSouza, A. M., and Nagarajan, M. B., “Mutual connectivity analysis (MCA) for nonlinear functional connectivity network recovery in the human brain using convergent cross-mapping and non-metric clustering,” in [Advances in Self-Organizing Maps and Learning Vector Quantization ], 217–226, Springer (2016).
  * [62] Schmidt, C., Pester, B., Schmid-Hertel, N., Witte, H., Wismüller, A., and Leistritz, L., “A multivariate Granger causality concept towards full brain functional connectivity,” PloS one 11(4) (2016).
  * [63] Abidin, A. Z., Chockanathan, U., DSouza, A. M., Inglese, M., and Wismüller, A., “Using large-scale granger causality to study changes in brain network properties in the clinically isolated syndrome (CIS) stage of multiple sclerosis,” in [Medical Imaging 2017: Biomedical Applications in Molecular, Structural, and Functional Imaging ], 10137, 101371B, International Society for Optics and Photonics (2017).
  * [64] DSouza, A. M., Abidin, A. Z., Leistritz, L., and Wismüller, A., “Exploring connectivity with large-scale Granger causality on resting-state functional MRI,” Journal of neuroscience methods 287, 68–79 (2017).
  * [65] Chen, L., Wu, Y., DSouza, A. M., Abidin, A. Z., Wismüller, A., and Xu, C., “MRI tumor segmentation with densely connected 3d cnn,” in [Medical Imaging 2018: Image Processing ], 10574, 105741F, International Society for Optics and Photonics (2018).
  * [66] Abidin, A. Z., DSouza, A. M., Nagarajan, M. B., Wang, L., Qiu, X., Schifitto, G., and Wismüller, A., “Alteration of brain network topology in hiv-associated neurocognitive disorder: A novel functional connectivity perspective,” NeuroImage: Clinical 17, 768–777 (2018).
  * [67] Abidin, A. Z., Deng, B., DSouza, A. M., Nagarajan, M. B., Coan, P., and Wismüller, A., “Deep transfer learning for characterizing chondrocyte patterns in phase contrast x-ray computed tomography images of the human patellar cartilage,” Computers in biology and medicine 95, 24–33 (2018).
  * [68] DSouza, A. M., Abidin, A. Z., Chockanathan, U., Schifitto, G., and Wismüller, A., “Mutual connectivity analysis of resting-state functional MRI data with local models,” NeuroImage 178, 210–223 (2018).
  * [69] Chockanathan, U., DSouza, A. M., Abidin, A. Z., Schifitto, G., and Wismüller, A., “Automated diagnosis of HIV-associated neurocognitive disorders using large-scale Granger causality analysis of resting-state functional MRI,” Computers in Biology and Medicine 106, 24–30 (2019).
  * [70] MTA 1, “http://fcon_1000.projects.nitrc.org/indi/acpi/html/,” (2015).  Last accessed 20 August 2020.
  * [71] C-PAC, “https://fcp-indi.github.io/,” (2015).  Last accessed 20 August 2020.
  * [72] The INDI Team, “ACPI preprocessing pipelines,” (2015).  Last accessed 20 August 2020.
  * [73] Kendall, M. G., “The treatment of ties in ranking problems,” Biometrika 33(3), 239–251 (1945).
  * [74] Suykens, J. A. and Vandewalle, J., “Least squares support vector machine classifiers,” Neural processing letters 9(3), 293–300 (1999).