Ich versichere hiermit, dass ich die Arbeit selbstständig verfasst habe, keine anderen als die
angegebenen Quellen und Hilfsmittel benutzt habe, die wörtlich oder inhaltlich übernommenen
Stellen als solche kenntlich gemacht habe und die Satzung des Karlsruher Instituts für Technologie
zur Sicherung guter wissenschaftlicher Praxis beachtet habe.

Karlsruhe, den 31. August 2015

Diese Bachelorarbeit befasst sich mit der Klassifikation menschlicher Ganzkörperbewegung mittels Hidden Markov Modellen (HMMs). Die menschlichen Bewegungen wurden mittels eines optisches Motion-Capture-Systems mit passiven Markern aufgenommen. Im Laufe der Arbeit werden zunächst Merkmale (features) erläutert, die die so aufgenommenen menschlichen Ganzkörperbewegungen beschreiben. Weiterhin wurden neue Merkmale aus den vorhandenen Rohdaten abgeleitet. Ein weiterer Schwerpunkt der Arbeit findet sich in der Diskussion von verschiedenen Ansätzen zur Lösung des Multi-Label-Klassifikationsproblems, also der Zuweisung von mehreren Klassen auf eine Bewegung. Hierbei wurden zwei grundsätzliche Ansätze erarbeitet: Beim ersten Ansatz wird das Multi-Label-Problem zu einem Multi-Class-Problem transformiert, welches sich dann einfach lösen lässt. Der zweite Ansatz beschäftigt sich mit Möglichkeiten, das Multi-Label-Problem ohne vorhergehende Transformation zu lösen, indem die Bewertungen einer unbekannten Bewegung durch die HMMs zu einer Gesamtvorhersage zusammengesetzt werden. Alle Aspekte des Klassifikators wurden anschließend auf einem aus $454$ Bewegungen bestehenden Datensatz evaluiert. Hierbei wurden verschiedene Parameter und Konfigurationen verglichen. Weiterhin fand ein Vergleich zwischen HMMs und Factorial Hidden Markov Modellen (FHMMs) statt. Die zwei zuvor genannten Ansätze wurden in zwei verschiedenen Systemen realisiert und quantitativ miteinander verglichen. Hierbei erkannte das erste System $98,02\%$ und das zweite System $93,39\%$ der menschlichen Bewegungen auf dem Testdatensatz korrekt.

CHAPTER: INTRODUCTION

Human whole-body motion plays an important role in many fields, including sports, medicine, entertainment,
computer graphics and robotics. Motion capture technologies are readily available to relatively
easily record vast amounts of data. Today, whole characters in movies
and video games are created using computer-generated imagery (CGI) and the recorded movements of an actor.
An impressive example is the completely computer-generated character Gollum from the Lord
of the Rings movies. In sports and medicine, motion data can be used for gait analysis of humans.
For example, stroke patients have been recorded using motion capture techniques to analyze their
disease-induced walking disorders. The insights from this analysis can be used to better understand
the symptoms of the patient and potentially allows the development of new rehabilitation therapies.
Gait analysis is also used in professional sports to analyze the performance of athletes and to optimize
their training. In humanoid robotics, human whole-body motion plays an important role
in the construction and control of biped humanoid robots. The use of motion in the field of robotics will be discussed
in-depth in the Related Works chapter.

Given the interest in human motion and the availability of the recording equipment, a growing number of motion
data is recorded. A number of human whole-body motion databases exist today to make this data available to artists, physicians and researchers.
Retrieval of information from these databases is often realized through so-called tags which are associated with a motion record.
A motion can potentially have many tags that describe it. For example, a motion of a human playing tennis
might be labeled with the tags tennis, forehand and right hand. These tags can be used
to query the database for motion data of interest. However, annotating a motion with these tags is usually
done by hand which is both an error-prone and slow process. The results are highly subjective since different
annotators will use different tags to label the same data. This may be because different annotators have a different
understanding of a tag or simply because they are not aware of the full set of available tags. Additionally,
labeling every motion by hand becomes infeasible as more data is recorded.

Autonomous motion recognition and classification can be used to automatically label new motions without
human involvement. This approach solves the two main problems of labeling by hand. Firstly, the classification algorithm
produces objective and reproducible results. Secondly, the system can be scaled to handle more and more
data by increasing the available computational resources. The objective of this thesis is to develop
a system that can be used to perform such autonomous classification of human whole-body motion.

This thesis is organized as follows: Chapter <ref> provides a brief overview of the
relevant literature and important authors in the field. Some fundamental concepts are introduced in chapter <ref> which forms
the basis for all following chapters. This includes an in-depth discussion of motion capture systems and
ways to represent this recorded data. The KIT Whole-Body Human Motion Database is introduced since
the system developed during this thesis will be used to classify motions stored in this database. Additionally,
the foundations of Hidden Markov Models (HMMs) and Factorial Hidden Markov Models (FHMMs) are discussed
since they play a key role for the devised classifier. Chapter <ref>
discusses possible features to discriminate motions. The discussion includes different representations
of motion data, the extraction of additional features from the data and important preprocessing steps
before the features can be used in a classifier. Chapter <ref> is concerned
with the problem of classification. In the chapter the previously discussed basics of Hidden Markov Models
are extended so that they can be used for motion recognition. The remainder of the chapter is devoted
to the use of multiple HMMs for the classification of human motion. The theoretical concepts
discussed in the previous chapters are then evaluated in chapter <ref>. Besides a description
of the used tools and dataset, the evaluation includes the selection of features, a comparison between different
HMM configurations as well as the evaluation of different classification approaches. Finally, the best
results from each area are used to perform an end-to-end evaluation of the entire system. Chapter <ref>
summarizes the work and describes possible improvements for future work.

CHAPTER: RELATED WORK

In robotics, a promising idea is to use human motion as an intuitive way to
instruct and program machines. This approach is commonly referred to as Programming by Demonstration
(PbD) [18, 7]. In PbD, a human instructor teaches a machine how to complete a
given task by performing the necessary steps themselves. The machine then observes the instructor
and attempts to also complete the task by imitating what it perceives.
Human motion can also be used to gain a better understanding of how different parts of the human body work together to complete a task or goal.
For example, observing how a human reacts to counter balance perturbations can potentially be used to transfer this knowledge to biped humanoid robots [39, 11].
Since motion plays a key role in humanoid robotics most of the hereafter mentioned literature has a background in this
area of research.

Different approaches exist for representing motions. Ogata et al. [43] used
Recurrent Neural Networks (RNNs) for interactive learning. In their work, RNNs were
used in a cooperative navigation task with a humanoid robot and a human partner. Taylor et al. [55, 53]
proposed Conditional Restricted Boltzmann Machines (CRBMs) to learn from and then generate human whole-body motion. The proposed
model is capable of generating continuous motion sequences (e.g. walking) and also allows to smoothly transition between them by adding higher-order layers to the model.
Nonlinear Oscillators were used by Nakanishi et al. [42] in a framework for learning
biped locomotion. Breazeal et al. [6] represented motion as a path through a directed weighted graph where
each node represents a pose. The edges of the graph define transitions between poses that are physically possible and safe.
Calinon et al. [16] used a Mixture Model of Gaussian and Bernoulli distributions (GMM/BMM)
to encode motion data. Yamane et al. [69] represented continuous motions in binary trees which can be used for
motion recognition and generation. Lastly, Hidden Markov Models (HMMs) have been a popular choice to represent human whole-body motion.
Since this work is concerned with Hidden Markov Models, the following paragraphs review some works in this domain in greater detail.

Takano at al. [60] developed a system for recognizing and generating human motion
for primitive nonverbal communication. Their approach uses a hierarchy of Hidden Markov Models for human
whole-body motion recognition and motion generation. The lower layer represents motion primitives, also referred to as proto symbols,
whereas the upper layer models the transitions between the motion primitives and therefore represents higher-level
interactions. The lower layer of the system was trained on joint angle data recorded with an optical motion capture system.
Multi-dimensional scaling was used to construct a multi-dimensional space of proto symbols (the proto symbol space)
on which the Hidden Markov Model in the upper layer was trained. The authors evaluated their approach by
recording a kickboxing match between two humans. One of the human subjects was then replaced with
a humanoid robot. The model trained on the recorded data was used by the robot to generate and perform motions in response to
the actions of its human counterpart.

Kulić et al. [36, 35, 37] proposed a system for learning, clustering
and hierarchy formation of human whole-body motion in humanoid robots. The authors used Hidden Markov Models and Factorial Hidden Markov Models to represent
motion as a sequence of motion primitives. Additionally, the system described by the authors is capable of on-line learning. This was achieved
by two essential properties of the devised system: sequential training of FHMMs and incremental hierarchical
formation of the motion primitives by clustering. The sequential training algorithm allowed Kulić et al. to initially encode an observed motion
into a simple Hidden Markov Model. As more and more data is observed, additional chains can be added and trained on-line, transforming the
HMM into an FHMM. Secondly, newly observed motions are dynamically organized into an hierarchical tree structure, the motion symbol tree.
This can be done efficiently by performing a tree search and placing the new motion into the node that is most similar.
Local clustering is performed to split groups into new subgroups as new knowledge is added. As a result, specialized motions
are placed at the leaves of the tree, whereas more generalized motions can be found near the root. Both properties allow
a humanoid robot to incrementally and efficiently build, organize and access knowledge during operation. The authors evaluated
their work with a database of recorded human whole-body motion data. The data set contained 28 motions for walking, 15 cheering motions,
7 dancing motions, 19 kicking motions, 14 punching motions, 13 sumo leg raise motions, 13 squatting motions, 13 throwing motions, and 15 bowing
motions. Each motion was represented in a humanoid model with 20 degrees of freedom. The results indicated that Factorial
Hidden Markov Models outperform single-chain Hidden Markov Models in their discriminative and generative

The work by Kulić et al. was extended in [56] and [34]. The authors used the motion symbol tree to perform efficient classification. This was achieved
by traversing the tree from the root and only recursively considering the subtree with the highest likelihood. As soon as a leave node is reached, the classification
is complete. The approach greatly decreased the computational cost since fewer comparisons are required in order to classify an unknown motion. The authors
further introduced the concept of a motion symbol graph. This directed graph allowed the authors to model likely transitions between motion symbols when observing continuous
motion. The motion symbol graph was used to predict motion patterns
and to generate motions in humanoid robots that consists of sequences of motion primitives.

In [54, 57, 58], Takano et al. proposed a system for mapping between motion symbols and word labels. Motions were encoded into HMMs
and the distances between all models was calculated. Like in the earlier work of Takano, the distance measures was used to construct a multi-dimensional space, the motion symbol space. Multiple word labels
associated with the motion primitive were encoded into a binary vector, which can be seen as a point in word label space. Finally, a linear mapping between the motion symbol space
and the word label space was learned using Canonical Correlation Analysis (CCA). CCA attempts to find a mapping in such a way that the correlation of the positions of motion symbols and
word labels is maximized. An advantage of this model is that it can be used to map from motion symbol to word label and
vice versa. This means that the system is both capable of motion classification given an unknown motion and motion retrieval given a query of word labels. In the latter case, since motion symbol
space and word label space are metric spaces, it is also possible to calculate the distance between a word label query and a motion symbol, making it easy to quantify the similarity.

An interesting extension to Hidden Markov Models was proposed in [66]. Wilson et al. used Parametric Hidden Markov Models (PHMMs)
to recognize parameterized gestures. An example of such a gesture is the movement of the hands that accompanies the speech “I saw a robot this big!”.
Here, this is a parameter of the gesture, namely the scalar size of the observed robot. The authors showed that traditional HMM-based recognition
cannot adequately model this spatial variance. Furthermore, HMMs do not allow to estimate the parameter (e.g. the size of the robot) from an unknown gesture.
The PHMM devised by the authors can solve both problems efficiently. It works by weighting a parameter vector and adding it to the mean
of the emission distribution of each hidden state. A modified version of the Baum-Welch algorithm was used to estimate the weights of the parameter vector. Recognition with PHMMs is complicated
by the fact that the parameter vector is unknown. This was solved by estimating the parameter vector (using an EM algorithm) for the observed sequence and each PHMM. The PHMM with the highest likelihood
was then selected. Furthermore, the authors extended PHMMs to a non-linear mapping from parameter vector to the means of the emission distributions. In this case, gradient ascent techniques were used to estimate the necessary parameters.

Herzog et al. [27] used a variation of PHMMs to recognize and imitate motions in humanoid robots. Their approach differed from the model proposed by Wilson et al. The basic idea
proposed by the authors is to use linear interpolation of HMMs that were trained on known parameters to generate a new HMM for new parameters. In their work the authors further discussed
how a humanoid robot can generate motions from such a model. The authors evaluated their approach on pointing and reaching motions and were able to show that, in those cases, PHMMs outperform traditional HMMs in classification.

Krüger et al. [32] used PHMMs for action recognition. Their work builds on the idea that
an action can be represented by a sequence of action primitives. The authors proposed a system that
used unsupervised segmentation to discover the action primitives. PHMMs were then used to encode and recognize them,
as well as synthesis motions with a desired effect (e.g. grabbing an object).

CHAPTER: BASICS

Classification of human whole-body motion first and foremost requires motion data. This chapter therefore
starts with a brief discussion of motion capture (section <ref>) and the Master Motor Map as
a framework for representing motion (section <ref>). Section <ref> gives an overview
of the KIT Whole-Body Human Motion Database, which plays an important role in this work since it
stores all motion data and also provides structures for labeling motions. Hidden Markov Models are
introduced in section <ref>, which are used to learn and recognize motions in this work. Lastly,
an extension of HMMs, Factorial Hidden Markov Models, are discussed (section <ref>).

§ MOTION CAPTURE

For recording motion data, the VICON MX motion capture system can be used. The system uses
passive optical markers that can be attached to both humans and objects. Cameras that are
positioned at multiple locations around the scene record the position of the markers within line of
sight. To do so, each camera features a ring of LEDs that surrounds its lens. The LEDs emit light in
the infrared spectrum, which is then reflected by the markers. Each camera records this reflected
light and (depending on the mode of operation) reports the 2D coordinates of the markers. The final 3D
coordinates for each marker are calculated by triangulation using the data from each camera [61].

Marker placement on the human body [41].

The motions are recorded by eight stationary and two portable
VICON T10 cameras. Each camera records with a sampling rate of 100. A total of 56 markers
are placed onto the human as depicted in figure <ref>.

All recorded motion data is stored in the C3D format. The C3D format is a binary file format under
public domain and is considered an industry standard. Besides storing marker coordinates in 3D space,
it also allows to store information about the human (e.g. body size and weight), the experiment setup
(e.g. marker positions) as well as additional data (e.g. data from additional sensors like force sensors) [15].

§ MASTER MOTOR MAP

Master Motor Map (MMM) [59, 1] is a framework for representation, mapping and reproduction of human motions on humanoid robots. The fundamental goal of MMM is to map and unify different motions performed by different
humans and recorded with different motion capture systems to the MMM reference model. Motions represented
under this reference model can then be converted to different outputs, e.g. to map a human motion onto
a humanoid robot like ARMAR-III [2]. The architecture of the MMM framework is depicted in
figure <ref>. The framework includes different command-line and graphical user interface (GUI) tools
and is open source[<http://h2t.anthropomatik.kit.edu/752.php>].

Architecture of the Master Motor Map framework [59].

At the core of the framework is the MMM reference model. It consists of a model of the human body with a
normalized height and weight, as well as kinematic and dynamic properties.
These properties are based on the research conducted by Winter et al. [67, 68]
and Buchholz et al. [4].
The kinematics of the MMM model consist of 104 degrees of freedom (DoF): 6 DoF cover the model pose,
23 DoF are assigned to each hand, and the remaining 52 DoF are distributed on arms, legs, head, eyes
and body. The reference coordinate system in every joint is chosen in such a way that the $x$-axis points to the right of the model,
the $y$-axis to the front and the $z$-axis upwards. If a joint has multiple DoF it is
split into multiple joints with a single DoF each. The MMM reference model also specifies upper and
lower limits for each joint. It is important to note that not all joints in the model must be used
to represent motion. For example, the movement of individual fingers might not be of interest
when recording a whole-body walking motion. In this case, the unspecified joints will simply remain
in their initial positions [59].

The command-line tool MMMConverter can be used to convert data recorded with a motion capture system to the
MMM reference model, i.e. reconstruct joint angles of the MMM model from motion data. This is accomplished by placing virtual markers onto the reference model and finding
a mapping from the position of the physical markers (as recorded by the motion capture system) to
virtual markers. The optimization problem can be solved efficiently by minimizing the distance between the position
of physical and virtual markers for each frame. Details on this mapping procedure are given in [59].
The converted motion is then stored in a XML-based file format. Such a XML file can contain multiple motions, which
is useful for scenes where a human interacts with objects or scenes with multiple humans in them. Each motion
is referenced by a unique name and consists of two parts: a preamble and the actual motion data.
The preamble specifies a model (e.g. the MMM reference model or a model of an object)
and can contain additional information about the human or object (e.g. body size and weight). The
motion data is encoded in a list of frames. Each frame has a relative time step in seconds (the first
frame starts at time step $t = 0$) and contains values for all properties that are specified by the model.
For example, each frame of human motion under the MMM reference model contains the root position ($x$, $y$, $z$ coordinates),
the root rotation (roll, pitch, yaw angles) and a list of joint angles. Additionally, the velocity and
acceleration information for each of the above properties as well as dynamic data (e.g. center of mass, angular momentum) can be stored for each frame [40].

The GUI tool MMMViewer can be used to visualize motions.
The whole motion can be played back or each frame can be inspected individually. The camera can be moved freely
to view the motion from different angles. Additionally, the joint angles of the currently visible frame
are displayed. Figure <ref> shows the tool during a visualization.

Visualization of a walking motion in MMMViewer.

§ KIT WHOLE-BODY HUMAN MOTION DATABASE

The KIT Whole-Body Human Motion Database[<https://motion-database.humanoids.kit.edu>]
contains motion data of both humans and objects that have been recorded using a marker-based approach
as described in section <ref>. Each entry in the database has a unique ID, belongs to a project,
and references the subjects and objects that participated in the recording. When recording motions,
multiple trials are usually performed. For each trial, the raw motion data is stored in the C3D format (see
section <ref>) and uploaded. Recorded data of additional sensors (e.g. force measurements for push recovery)
as well as video footage can be uploaded as well. The database system automatically converts the
C3D files to a subject-independent representation under the MMM model (see section <ref>)
and, optionally, estimates dynamic properties like the center of mass. Log files grant insight into
the conversion process. Furthermore, the database is capable of storing data related to subjects and
objects. Size, weight, gender and other anthropometric measurements can be stored in the record for
each subject. For objects, a 3D model of the object alongside a custom description can be
saved [41].

The Motion Description Tree (excerpt) [41].

Each motion is classified within the Motion Description Tree. The tree consists of a hierarchical declaration
of tags describing motion types (e.g. walk, kick, run) and additional tags for other properties like the
direction of a movement (e.g. left, right, forward, backward). The tree is organized in such a way
that the parent of a node has a broader semantic meaning than its children. For example, the tags clap
and bow are both child nodes of the tag gesticulation. An excerpt of the Motion
Description Tree is depicted in figure <ref>. An important property of
this classification approach is that a motion can be associated with an arbitrary number of nodes of
the Motion Description Tree. For example, a motion of a human that trips while walking to the left
with high speed can be categorized using the following tags:

* locomotion $\rightarrow$ bipedal $\rightarrow$ walk,

* speed $\rightarrow$ fast,

* direction $\rightarrow$ left,

* perturbation $\rightarrow$ result $\rightarrow$ failing, and

* perturbation $\rightarrow$ source $\rightarrow$ passive.

The whole tree is managed by the KIT Whole-Body Human Motion Database and can be extended if necessary [41].

The database can be accessed through a web interface or an application programming interface (API).
The web interface and the API are available publicly. For each motion the raw files as well as the processed files can
be downloaded. For convenience, bulk download options are available. The web interface is also used to
modify existing or upload new motions. These operations are restricted to registered accounts. The API
allows direct access to the database. This allows the integration of the database into existing
tools. The API is build on top of the Internet Communications Engine (Ice) [28].
Ice is a remote procedure call (RPC) framework and allows for easy integration with a wide variety
of platforms and programming languages [41].

At the time of this writing, the database contains $4\,457$ motions performed by $49$ different subjects. All motions
in total have a length of approximately $9$ hours and $20$ minutes, with the average length of a recording being
approximately $7.56$ seconds.

§ HIDDEN MARKOV MODELS

A Hidden Markov Model (HMM) [19, 47] is a statistical model popular for learning
sequential data. This is due to the fact that HMMs have the ability to have some degree of invariance
to local warping (compression and stretching) of the time axis [3]. The methods
discussed in this section are applicable to all forms of sequential data. However, since this work deals
with temporal sequences, this section and all following chapters use notation and phrases that imply temporal
sequences. Concretely, a temporal sequence of length $T$ is denoted by $\vec{o}_1, \vec{o}_2, \ldots, \vec{o}_T$,
where each $\vec{o}_t$ is a multi-dimensional observation. The following discussion of HMMs and
the underlaying concepts are based on Bishop et al. [3].

[->,scale=1.,auto=left,thick,every node/.style=circle,thick,minimum size=1.1cm]

(x1) at ( 0,0) ;
[observation] (x2) at ( 3,0) $\vec{o}_{t-1}$;
[observation] (x3) at ( 7,0) $\vec{o}_t$;
[observation] (x4) at (11,0) $\vec{o}_{t+1}$;
(x5) at (14,0) ;

/→in x1/x2, x2/x3, x3/x4, x4/x5
() – (→);

Graphical representation of a first-order Markov chain. Each observation $\vec{o}_t$ is only dependent on the previous
observation [3].

To understand Hidden Markov Models, it is helpful to first consider a simpler Markov model: the
Markov chain. In a Markov chain (of first order), given a sequence $\vec{o}_1, \ldots, \vec{o}_T$, the conditional
probability of an observation $\vec{o}_t$ (that is the observation at time $t$) is assumed to be independent
of all past (and, of course, future) observations except for observation $\vec{o}_{t-1}$. The conditional
distribution in a such a model is given by

\begin{equation}
p(\vec{o}_t \mid \vec{o}_1, \ldots, \vec{o}_{t-1}) = p(\vec{o}_t \mid \vec{o}_{t-1}).
\end{equation}

Consequently, the joint distribution is given by

\begin{equation}
p(\vec{o}_1, \ldots, \vec{o}_T) = p(\vec{o}_1) \prod\limits_{t=2}^T{p(\vec{o}_t \mid \vec{o}_{t-1})}.
\end{equation}

A graphical representation of a first-order Markov chain is depicted in figure <ref>.
The assumption that an observation is only dependent on its previous observation is rather strong.
This can be easily seen by considering an example: If one attempts to predict the weather for the next
hour by only considering the current weather situation instead of using the data of the last
24 hours, the prediction would be severely limited. The assumption can be relaxed by generalizing
the Markov chain to be of $M$-th order. Here, each observation in a sequence is dependent on the
past $M$ observations. However, in such a model the number of parameters grows exponentially with M,
so that this approach becomes impractical for large values of $M$.

[->,scale=1.,auto=left,thick,every node/.style=circle,thick,minimum size=1.1cm]

(z1) at ( 0,2) ;
[hidden] (z2) at ( 3,2) $\vec{z}_{t-1}$;
[hidden] (z3) at ( 7,2) $\vec{z}_t$;
[hidden] (z4) at (11,2) $\vec{z}_{t+1}$;
(z5) at (14,2) ;

(x1) at ( 0,0) ;
[observation] (x2) at ( 3,0) $\vec{o}_{t-1}$;
[observation] (x3) at ( 7,0) $\vec{o}_t$;
[observation] (x4) at (11,0) $\vec{o}_{t+1}$;
(x5) at (14,0) ;

/→in z1/z2, z2/z3, z3/z4, z4/z5
() – (→);
/→in z2/x2, z3/x3, z4/x4
() – (→);

Graphical representation of a Hidden Markov Model. Each observation $\vec{o}_t$ (in blue) is conditioned
on the state of its respective hidden variable $\vec{z}_t$ (in red). The hidden variables form a first-order
Markov chain [3].

To solve this problem, hidden (sometimes also referred to as latent) variables are
introduced. Concretely, each observation variable $\vec{o}_t$ is conditioned on the state of its hidden variable
$\vec{z}_t$. The hidden variables form a first-order Markov chain. Such a model is known as a state space model,
which is visualized in figure <ref>. The joint distribution for this model is given by

\begin{equation}
p(\vec{o}_1, \ldots, \vec{o}_T, \vec{z}_1, \ldots, \vec{z}_T) = p(\vec{z}_1) \prod\limits_{t=2}^T{p(\vec{z}_t \mid \vec{z}_{t-1})}
\prod\limits_{t=1}^T{p(\vec{o}_t \mid \vec{z}_t)}.
\end{equation}

An important property of this model is that any pair of observed variables $\vec{o}_i$ and $\vec{o}_j$
are connected via the hidden variables. It can be shown that the predictive distribution
$p(\vec{o}_{t+1} \mid \vec{o}_1,\ldots,\vec{o}_t)$ for observation $\vec{o}_{t+1}$ is dependent on all past
observations $\vec{o}_1,\ldots,\vec{o}_t$. This model is therefore not
constrained by the strong independence assumption of a Markov chain. If the hidden variables in a
space state model as described above are discrete, the Hidden Markov Model is obtained.

Formally, a Hidden Markov Model is defined by the set of parameters that govern the model [3, 47]:

* $K$ is the number of the individual hidden states. The appropriate number of states depends on the problem
at hand. The current state at time $t$ is denoted by $q_t$ and the set of possible states is $S = \{s_1, \ldots, s_K\}$.
Notice the relationship of $K$ and the hidden variables: Each of the $K$ states describes a possible value
for the hidden variables $\vec{z}_t$. One way to represent this is through the $1$-of-$K$ coding scheme, hence $\vec{z}_t \in \{0,1\}^K$ and
$\norm{\vec{z}_t}_1 = 1$ for each $t$. For example, the notation $q_t = s_2$ is equivalent to $\vec{z}_t = (0, 1, 0, \ldots, 0)$.

* $\vec{A} = (a_{ij}) \in [0,1]^{K \times K}$ is a matrix that consists of transition probabilities. Concretely,
each entry defines the probability to transition to state $s_j$ given the current state $s_i$:

\begin{equation}
a_{ij} = p(q_{t+1}=s_j \mid q_t=s_i).
\end{equation}

Since $\vec{A}$ is a probability distribution, it must hold that $\forall j: \sum\nolimits_k{a_{kj}}=1$. By setting
$a_{ij}=0$, it is possible to “disable” that specific transition from $s_i$ to $s_j$.

* $\vecsym{\pi}=(\pi_{1}, \ldots, \pi_{K}) \in [0,1]^K$ is the initial state distribution where

\begin{equation}
\pi_{i} = p(q_1=s_i).
\end{equation}

It must hold that $\norm{\vecsym{\pi}}_1 = 1$.

* $\phi$ describes the parameters of the conditional distributions of the observed variables:

\begin{equation}
p(\vec{o}_t \mid \vec{z}_t, \phi).
\end{equation}

These probabilities are known as emission probabilities and can be given by different
distributions. For example, if the observed values are discrete, a conditional probability table
can be used. For observations with continuous values, a Gaussian distribution is a often a good choice.
Other distributions are possible and picking an appropriate distribution depends on the observations.

Since $K$ is already encoded by the shape of $\vec{A}$, an HMM is fully described by the following set of parameters:
$\theta = \{\vec{A}, \vecsym{\pi}, \phi\}$. In this work, an HMM with parameters $\theta$ is denoted by $\lambda_\theta$.

Three fundamental problems can be identified when working with HMMs:

* The evaluation problem,

* the decoding problem, and

* the optimization problem.

The problems and their description are all based on the work of Rabiner et al. [47].
The first problem, the evaluation problem, is concerned with calculating the probability of a given
sequence under a given model. Formally, given a sequence $\vec{O} = (\vec{o}_1,\ldots,\vec{o}_T)$, how can
$p(\vec{O} \mid \lambda_\theta)$ be calculated efficiently. This can also be viewed as scoring how well
a model matches the given observations. The forward-backward algorithm [8, 14]
can be used to solve this problem (strictly speaking, only the forward pass is necessary to solve this first problem).
The second problem, the decoding problem, is concerned with finding the state sequence of the hidden
variables. Formally, given a sequence $\vec{O} = (\vec{o}_1,\ldots,\vec{o}_T)$ and a model $\lambda_\theta$, find a sequence
$\vec{Z} = (\vec{z}_1, \ldots, \vec{z}_T)$ of hidden states that is optimal. Different criteria of an optimal
state sequence exist, e.g. choosing the states that are individually most likely. The most popular criterion
is to find the single best state sequence. This is equivalent to maximizing $p(\vec{Z} \mid \vec{O}, \lambda_\theta)$, which is solved
efficiently by the Viterbi algorithm [62].
The third problem, the optimization problem, is concerned with adjusting the parameters of the model.
Formally, given a sequence $\vec{O} = (\vec{o}_1,\ldots,\vec{o}_T)$, find parameters $\theta$ such that $p(\vec{O} \mid \lambda_\theta)$
is maximized. Solving this problem corresponds with learning the parameters, that is “training” the model.
The Baum-Welch algorithm [12] solves this problem efficiently. However, since
Baum-Welch is a specific case of the expectation maximization algorithm (EM algorithm), it does
not necessarily find a global maximum.

On a final note, another important property of HMMs is that they are generative models. This means that a model
that has been trained on some data can be used to generate new samples. This is especially interesting
for human motions and humanoid robots. Here, a motion can be learned by observation and later be reproduced in a robot
by sampling from the model [60].

§ FACTORIAL HIDDEN MARKOV MODELS

A sever limitation of HMMs is that they cannot represent a lot of information about the history of a time
sequence. Factorial Hidden Markov Models (FHMMs) are a generalization of HMMs and offer a way to overcome this limitation.
For example, representing $30$ bit of information about the history requires $2^{30}$ hidden states in a
standard HMM whereas an FHMM can represent the same information with only $30$ binary state variables [23].
The discussion in this section is based on the work of Ghahramani et al. [23].

[->,scale=1.,auto=left,thick,every node/.style=circle,thick,minimum size=1.1cm]

(z11) at ( 0,6) ;
[hidden1] (z12) at ( 3,6) $\vec{z}^{(1)}_{t-1}$;
[hidden1] (z13) at ( 7,6) $\vec{z}^{(1)}_t$;
[hidden1] (z14) at (11,6) $\vec{z}^{(1)}_{t+1}$;
(z15) at (14,6) ;

(z21) at ( 0,4) ;
[hidden2] (z22) at ( 3,4) $\vec{z}^{(2)}_{t-1}$;
[hidden2] (z23) at ( 7,4) $\vec{z}^{(2)}_t$;
[hidden2] (z24) at (11,4) $\vec{z}^{(2)}_{t+1}$;
(z25) at (14,4) ;

(z31) at ( 0,2) ;
[hidden3] (z32) at ( 3,2) $\vec{z}^{(3)}_{t-1}$;
[hidden3] (z33) at ( 7,2) $\vec{z}^{(3)}_t$;
[hidden3] (z34) at (11,2) $\vec{z}^{(3)}_{t+1}$;
(z35) at (14,2) ;

(x1) at ( 0,0) ;
[observation] (x2) at ( 3,0) $\vec{o}_{t-1}$;
[observation] (x3) at ( 7,0) $\vec{o}_t$;
[observation] (x4) at (11,0) $\vec{o}_{t+1}$;
(x5) at (14,0) ;

/→in z11/z12, z12/z13, z13/z14, z14/z15
() – (→);
/→in z21/z22, z22/z23, z23/z24, z24/z25
() – (→);
/→in z31/z32, z32/z33, z33/z34, z34/z35
() – (→);
/→in z12/x2, z13/x3, z14/x4
() edge [bend left=60] (→);
/→in z22/x2, z23/x3, z24/x4
() edge [bend left=40] (→);
/→in z32/x2, z33/x3, z34/x4
() – (→);

Graphical representation of a Factorial Hidden Markov Model. Each observation $\vec{o}_t$ (in blue) is conditioned
on the state of all its respective hidden variables $\vec{z}^{(1)}_t$ (in red), $\vec{z}^{(2)}_t$ (in green) and
$\vec{z}^{(3)}_t$ (in purple). The hidden variables $\vec{z}^{(m)}_t$ form a first-order Markov chain each. [23].

In an FHMM, the current state is generalized by letting the state be represented by a collection of $M$
state variables:

\begin{equation}
\vec{z}_t = (\vec{z}^{(1)}_t,\ldots,\vec{z}^{(M)}),
\end{equation}

where each state variable $\vec{z}^{(m)}_t$ can take on $K$ different values. Each state
variable is constrained in such a way that it evolves according to its own dynamics and is therefore
uncoupled from the other state variables:

\begin{equation}
p(\vec{z}_t \mid \vec{z}_{t-1}) = \prod\limits_{m=1}^M{p(\vec{z}^{(m)}_t \mid \vec{z}^{(m)}_{t-1}}).
\end{equation}

This can be seen as $M$ independent first-order Markov chains that all contribute to the observation. The distribution of the observed variable
$\vec{o}_t$ is conditional on the states of all hidden variables $\vec{z}^{(1)}_t, \ldots, \vec{z}^{(M)}_t$
for each time step $t$. A simple way to represent this dependency for continuous observations is a
multivariate Gaussian. Concretely, given that each $\vec{z}^{(m)}_t$ uses the $1$-of-$K$ coding scheme as described
in section <ref> and that each observation $\vec{o}_t \in \mathbb{R}^D$, the conditional distribution
is given by

\begin{equation}
p(\vec{o}_t \mid \vec{z}_t) = \abs{\vecsym{\Sigma}}^{-\frac{1}{2}} (2 \pi)^{-\frac{D}{2}} \exp\left(-\frac{1}{2}
(\vec{o}_t-\vecsym{\mu}_t)^T \vecsym{\Sigma}^{-1} (\vec{o}_t-\vec{\mu}_t)\right),
\label{eq:fhmm-gaussian}
\end{equation}

\begin{equation}
\vecsym{\mu}_t = \sum\limits_{m=1}^M{\vec{W}^{(m)}\vec{z}^{(m)}_t}.
\label{eq:fhmm-mean}
\end{equation}

Here, each $\vec{W}^{(m)}$ is a $D \times K$ matrix that weights the contributions to the mean for each
$\vec{z}^{(m)}_t$, $\vecsym{\Sigma}$ is the $D \times D$ covariance matrix and $\abs{\cdot}$ denotes the matrix determinant.
Also note that the scalar $\pi$ in equation <ref> is not to be confused with the initial probability vector
$\vecsym{\pi}$ from
section <ref>. In words,
at each time step, the state of all chains are weighted, summed and output through an expectation function
(here equation <ref>) to produce the observation [37].

Like an HMM, a Factorial Hidden Markov Model is defined by the set of parameters that govern the model:
$\theta = \{\vec{A}^{(1)}, \ldots, \vec{A}^{(M)}, \vecsym{\pi}^{(1)}, \ldots, \vecsym{\pi}^{(M)}, \phi\}$. This is a simple extension of the parameters of a standard normal HMM: The transition probabilities $\vec{A}^{(1)}, \ldots, \vec{A}^{(M)} \in \mathbb{R}^{K \times K}$ and initial probabilities
$\vecsym{\pi}^{(1)}, \ldots, \vecsym{\pi}^{(M)} \in \mathbb{R}^K$ must be given for each of the $M$ Markov chains. $\phi$ still defines
the necessary parameters for the emission probability distribution, which is $\phi = \{\vec{W}^{(1)}, \ldots, \vec{W}^{(M)}, \vecsym{\Sigma}\}$ for
a Gaussian FHMM as described above.

A problem with FHMMs is learning their parameters. This is because although at each time step the hidden variables are marginally
independent, they become conditionally dependent given the observation sequence. This can be easily seen by considering
equation <ref> and <ref> that makes the mean and therefore the entire Gaussian
a function of all states. As a result, exact inference becomes infeasible. Concretely, the backward-forward algorithm
used in the E step of the Baum-Welch algorithm has time complexity $\mathcal{O}(T M K^{M+1})$, where $T$ is the length of the sequence,
$K$ is the number of states and $M$ is the number of Markov chains. Note however that the M step for FHMMs is
completely tractable and can therefore be calculated exactly. To work around the infeasibility of inference, several approximations
of the E step have been proposed: Ghahramani et al. devised inference using Gibbs sampling,
completely factorized variational inference and structured variational inference.
A fourth approach, the generalized backfitting algorithm, was described by Jacobs et al. [30].

CHAPTER: FEATURES

As already mentioned in section <ref>, motions can be recorded using an optical marker-based
motion capture system. The following section discusses different approaches to represent such motions
and describes possible features that are used to recognize and classify them in later chapters.

§ MARKER REPRESENTATION

A natural and obvious way to represent the recorded data is in 3-dimensional Cartesian space. For each
time sample $t$, the system records the location of each marker $n$:

\begin{equation}
\vec{r}_t^{(n)} = (x_t^{(n)},y_t^{(n)},z_t^{(n)}) \in \mathbb{R}^3.
\end{equation}

A complete motion or observation sequence $\vec{O}$ is then represented by the all marker locations for all sampled time steps. A way to write this
is to “unroll” all marker locations for a given time step $t$ into the $t$-th row of an observation matrix:

\begin{equation}
\vec{O}_{cartesian} = \begin{bmatrix}
x_1^{(1)} & y_1^{(1)} & z_1^{(1)} & \dots  & x_1^{(N)} & y_1^{(N)} & z_1^{(N)} \\
x_2^{(1)} & y_2^{(1)} & z_2^{(1)} & \dots  & x_2^{(N)} & y_2^{(N)} & z_2^{(N)} \\
\vdots    & \vdots    & \vdots    & \ddots & \vdots    & \vdots    & \vdots    \\
x_T^{(1)} & y_T^{(1)} & z_T^{(1)} & \dots  & x_T^{(N)} & y_T^{(N)} & z_T^{(N)}
\end{bmatrix} \in \mathbb{R}^{T \times 3N},
\end{equation}

where $T$ is the number of time samples and $N$ is the number of markers. A visualization of a motion and
the respective marker locations is depicted in figure <ref>.

Five key frames of a squatting motion with the location of the individual markers visualized. The subject looks to the left and is depicted from the side.

A problem with this representation is that an absolute coordinate system is used. Consider for example
two running motions. Assume that in the first case the subject moves towards a defined point and in the second
case turns 45 degrees and repeats the motion almost identically. However,
since an absolute coordinate system is used, the values in the observation matrix $\vec{O}$ will
be very different for the two almost identical motions. The same problem occurs if the start location of two motions is offset.
Again, similar motions will have different marker positions in the observation matrix. In short, the representation in an absolute
Cartesian coordinate system is neither invariant to translation nor rotation. A coordinate system that is relative to the
recorded subject is desirable to allow for robust motion recognition.

§ JOINT ANGLE REPRESENTATION

The Master Motor Map (MMM, [59]) framework (see also section <ref>) uses a relative coordinate system.
This is achieved by mapping the position of the physical markers onto virtual markers
on a reference model. To do so, the squared error between the physical and virtual markers is minimized by varying the pose
of the subject (as defined by its position and rotation in space as well as its joint angles) while maintaining the constraints of the reference
model. The optimization problem is solved by the reimplementation of the Subplex algorithm as provided by the NLOpt
library [39]. Figure <ref> depicts the kinematics and shows the location and labels of all joints.

The kinematics of the MMM reference model [59].

Some joints have multiple degrees of freedom (DoF). Take, for example, the body lower neck (BLN) joint that
has 3 DoF (to convince yourself that this is indeed the case, nod, shake your head and move your head from shoulder to shoulder).
In robotics, this is usually handled by splitting a joint with multiple degrees of freedom into multiple joints with a single
DoF each. In the case of the exemplary BLN joint, this means that the joints BLNx, BLNy and BLNz will replace it.
After this step, each joint has a single DoF and can therefore be represented as a scalar value that defines
the joint angle in radians. The complete joint configuration at time step $t$ can therefore be written as

\begin{equation}
\vecsym{\theta}_t = (\theta_t^{(1)},\theta_t^{(2)},\ldots,\theta_t^{(N)}) \in \mathbb{R}^N,
\end{equation}

where $N$ is the number of joints with a single DoF each. Similar to a representation in Cartesian space, a complete observation sequence
can then be written as:

\begin{equation}
\vec{O}_{mmm} = \begin{bmatrix}
\theta_1^{(1)} & \theta_1^{(2)} & \dots  & \theta_1^{(N)} \\
\theta_2^{(1)} & \theta_2^{(2)} & \dots  & \theta_2^{(N)} \\
\vdots    	   & \vdots    		& \ddots & \vdots   	  \\
\theta_T^{(1)} & \theta_T^{(2)} & \dots  & \theta_T^{(N)}
\end{bmatrix} \in \mathbb{R}^{T \times N},
\end{equation}

where $T$ denotes the number of time samples.

It is important to stress again that the the joint angle representation is relative to the subject.
However, some important information is lost: the position of the subject in space. To overcome this,
the absolute root position and the root rotation at each time step $t$ are included in the MMM framework as well:

\begin{align}
\vec{r}_t^{(root)} &= (x_t,y_t,z_t) \in \mathbb{R}^3, \\
\vecsym{\theta}_t^{(root)} &= (\theta_t^{(roll)}, \theta_t^{(pitch)}, \theta_t^{(yaw)}) \in \mathbb{R}^3.
\end{align}

However, these properties are once again given in an absolute coordinate system, suffering from the same
problem described above. Luckily, this is resolved easily by an affine transformation of the coordinate system
translating it such that the root position at $t = 0$ starts at the origin and rotating it such
that the $y$ axis points away from the front of the subject.
The following two equations describe the necessary calculations for each time step $t$:

\begin{align}
\Delta\vec{r}_t^{(root)} &= \vec{r}_t^{(root)} - \vec{r}_0^{(root)}, \\
\vec{\hat{r}}_t^{(root)} &= \vec{R}^{-1}~\Delta\vec{r}_t^{(root)}.
\label{eq:normalization-root}
\end{align}

The first equation describes the translation, whereas the second equation describes the rotation. A
rotation matrix for roll, pitch and yaw angles is given in [17]:

\begin{equation}
\vec{R} =
\begin{bmatrix}
\cos\alpha \cos\beta & \cos\alpha \sin\beta \sin\gamma - \sin\alpha \cos\gamma & \cos\alpha \sin\beta \cos\gamma + \sin\alpha \sin\gamma \\
\sin\alpha \cos\beta & \sin\alpha \sin\beta \sin\gamma + \cos\alpha \cos\gamma & \sin\alpha \sin\beta \cos\gamma - \cos\alpha \sin\gamma \\
-\sin\beta & \cos\beta \sin\gamma & \cos\beta \cos\gamma
\end{bmatrix},
\end{equation}

where $\alpha := \theta_0^{(yaw)}$, $\beta := \theta_0^{(pitch)}$ and $\gamma := \theta_0^{(roll)}$.

Figure <ref> plots the root position of two motions where the subject runs forward. The
same movements are depicted once before any normalization has been applied and once after normalization. Notice that without normalization,
the features are neither translation nor rotation invariant. This can be especially well seen in figure <ref>:
Although the subject is running only forward, the movement is split between the $x$ and $y$ components. This is because the subject moves
at an approximately 45 degree angle between the $x$ and $y$ axis of the absolute coordinate system. After normalization (figure <ref>),
the movement direction happens only in the direction of the $y$ axis of the transformed and now relative coordinate system.
Hence rotation and translation invariant features are obtained after normalization.

[The first running motion without normalization.]
[The second running motion without normalization.]
[The first running motion with normalization.]
[The second running motion with normalization.]

Comparison between the root positions of two different running motions before and after normalization.

The root rotation must be normalized as well to make it comparable. Since the rotation is given in
angles, the normalization is straightforward:

\begin{equation}
\vecsym{\hat{\theta}}_t^{(root)} = \vecsym{\theta}_t^{(root)} - \vecsym{\theta}_0^{(root)}.
\end{equation}

Notice that this work assumes that the roll, pitch and yaw angles are not limited to the interval $[-\pi,\pi]$.
If necessary, this is easily achieved by correcting overflows by adding $2\pi$ to all following angles
(and similarly subtracting $2\pi$ for underflows).

Another interesting usage of the MMM reference model is that it allows normalization of the marker positions.
Recall that the marker positions were previously given in an absolute coordinate system (see section <ref>).
However, since the initial pose of the subject is known in the MMM reference model, this information can
be used to normalize the marker positions as well. This is done analogously to the normalization of
the root position for the position of each marker (see equation <ref>).

§ DERIVED FEATURES

Multiple additional features are computed under the MMM reference model. Firstly, an obvious extension
is to calculate the velocities and accelerations of all features that describe positions in Cartesian coordinate
space. Given that the velocity is the first derivative of the position and the acceleration is the first
derivative of the velocity, both properties are easily calculated by approximating the respective derivatives:

\begin{equation}
\vec{v}_t = \frac{\vec{r}_{t+1} - \vec{r}_{t-1}}{2 \Delta t} \qquad \text{and} \qquad \vec{a}_t = \frac{\vec{v}_{t+1} - \vec{v}_{t-1}}{2 \Delta t},
\end{equation}

where $\vec{v}$ denotes the velocity, $\vec{a}$ the acceleration and $\Delta t$ is the time difference
between two subsequent samples (which is assumed to be equidistant over all samples). The velocity
and acceleration must be normalized as well. The normalization works similarly to
equation <ref> but normalizes each sample with the current pose of the respective
segment instead of normalizing each sample with the initial root pose.
An interesting modification of the velocities and accelerations is to reduce them to simple scalar values
by using their norm instead, e.g. the Euclidean one.

Comparison between the root position (left) and center of mass (right) of the same bowing motion.

Secondly, since the MMM reference model provides additional information about the subject, more advanced
features are computed as well. Two interesting dynamic properties are the center of mass (CoM) and
the angular momentum [44].
The CoM is conceptually similar to the root position discussed earlier in the sense
that it describes the position of a subject in 3-dimensional Cartesian space. However, while the root position is always
at a fixed point on the reference model, the center of mass is the barycenter of the subject. More concretely,
the center of mass is the average over the CoM positions of body segments weighted by their respective masses.
To clarify this, consider a motion where the subject performs a deep bow. A bowing motion is interesting in this case
since the lower body remains mostly fixed in space while the upper body moves down. Figure <ref> plots
the root position and the CoM of such a motion. Notice that the $z$ component of the root position
remains approximately constant although the upper body moves down during the bow. The $z$ component of the CoM on the other
hand first decreases while the subject bows down and then increases again as the subject comes back up.
Now consider the $y$ component of the root position: It decreases as the subject bows down since the hip of the subject moves backwards. Compare
this to the $y$ component of the CoM that instead increases as the subject bows down since that shifts the center of mass forward.
Like the root position, the CoM is given in an absolute coordinate system. The
normalization to a coordinate system relative to the subject is done analogous to the normalization
of the root position. The velocity and acceleration of the CoM are computed as well.

The angular momentum is a physical measure for the rotational configuration of an object or a system in 3D space.
In the MMM framework, the the angular momentum is calculated with respect to the center of mass at
each time step $t$ in all three spatial directions. The whole-body angular momentum is calculated as follows:

\begin{equation}
\vec{L}_t = \sum_{i=1}^M \left( m^{(i)} (\vec{r}_t^{(i)} \times \vec{v}_t^{(i)}) + \vec{I}_t^{(i)} \vecsym{\omega}_t^{(i)} \right) \in \mathbb{R}^3 ,
\end{equation}

\begin{equation}
\vec{r}_t^{(i)} = \vec{r}_t^{(CoM_i)} - \vec{r}_t^{(CoM)} \qquad\text{and}\qquad \vec{v}_t^{(i)} = \vec{\dot{r}}_t^{(CoM_i)} - \vec{\dot{r}}_t^{(CoM)}.
\end{equation}

The first part of the sum considers the angular momenta created by the orbital rotation of each segment around the whole-body center of mass. For each segment $i \in \{1,\ldots,M\}$,
$m^{(i)}$ describes its mass, $\vec{r}_t^{(i)}$ its position at time step $t$ w.r.t. the CoM and $\vec{v}_t^{(i)}$ its velocity at time step $t$ w.r.t. the CoM. The cross product is denoted as $\times$.
The second part of the sum takes the spin of each segment into account by computing the product of its inertia tensor $\vec{I}_t^{(i)}$ and its angular velocity $\vecsym{\omega}_t^{(i)}$.
The velocity $\vec{v}_t^{(i)}$ and the difference in CoM $\vec{r}_t^{(i)}$ must be normalized as previously described.

Thirdly, the position of body segments are used as features as well. Consider for example a waving
motion. In this case, the positions of the hands are an interesting feature. Similarly, the position
of the feet are interesting for other motions, e.g. a kick. The positions of the extremities must be
normalized. The velocities and accelerations are computed as previously described.

§ SMOOTHING

The described features can be noisy or contain errors. For example, noise is introduced during the
recording process. In addition, approximating the derivative for the computation of velocities and accelerations
can amplify errors and inaccuracies in the recorded data. Smoothing is used to reduce the impact
of these interferences. A signal can be smoothed using a wide variety of different filters. For a more complete
discussion, see [49].

In this work, only a simple filter is briefly discussed: the moving average or sliding window filter.
In such a filter of length $W$, the mean of the surrounding $W$ data points is used instead of the individual data point:

\begin{equation}
\vec{\hat{x}}_t = \frac{1}{W+1} \sum\limits_{j=-W/2}^{W/2}{\vec{x}_{t+j}},
\end{equation}

where $\vec{x}_t$ denotes a (potentially multi-dimensional) data point at time step $t$ and $\vec{\hat{x}}$ is the
smoothed version thereof. Including future samples into the average avoids introducing a time delay in the signal.
Notice that averaging over future samples is possible if the smoothing is applied off-line.

The components of a subject's root acceleration during a running motion before (left) and
after (right) a moving average filter with $W=3$ was applied.

Figure <ref> compares the root acceleration of
a subject before and after applying a moving average filter. Notice that the the original signal
is very noisy. The smoothed signal maintains the overall structure of the signal but reduces the
amount of jitter. Smoothing is a useful preprocessing step before feeding the features into
a model.

§ SCALING

Feature scaling is another preprocessing step. Take for example the joint angles and
the root position from the previous sections. The joint angles are physically constrained to a very narrow
value range, whereas the root position can potentially grow very large if the subject travels a large
distance from the start position. It should be obvious from this example that features are on very different scales.
This difference in scale becomes a problem if $k$-means clustering is used to initialize an HMM's emission
distribution parameters. If the data is on very different scales across dimensions, $k$-means will not
find clusters that properly fit the data because the same distance measure is minimized across dimensions.
This, in turn, results in bad estimates of the emission distribution parameters which results in vanishing probabilities
and numerical instabilities during inference. To counter this, feature scaling is performed. A very simple but effective strategy is to scale the features
such that each feature's values are in the same range, e.g. $[-1,1]$. This is achieved by applying
the following equation to each individual feature $x$ over all $T$ time steps:

\begin{equation}
\hat{x}_t = 2 \cdot \frac{x_t - \min\limits_{i\in\{1,\dots,T\}}{x_i}}{\max\limits_{i\in\{1,\dots,T\}}{x_i} - \min\limits_{i\in\{1,\dots,T\}}{x_i}} - 1.
\end{equation}

Notice that each feature needs to be scaled over all samples, not per-sample. Furthermore, if feature
scaling is used, the computation of the minimum and maximum are computed on the training data.
New samples (e.g. from the test dataset when evaluating or unknown motions when used productively)
are then simply scaled using the previously computed values. Otherwise the features that were
used to train the model and the features that are used to recognize unknown motions would end up on a different scales.

CHAPTER: CLASSIFICATION

This chapter discusses methods to recognize and classify motions represented by features as described
in the previous chapter. In the first section, the general concepts of Hidden Markov Models introduced
in chapter <ref> are concretely discussed for the case of motion recognition. In motion
recognition, the goal is to encode a motion into a Hidden Markov Model and compute a measure that describes
how likely an unknown motion is under the model. In the second section, motion recognition is extended
to perform multi-class classification. In contrast to motion recognition, the goal is now
to assign an unknown motion exactly one class out of a (potentially large) set of possible classes.
Lastly, methods for performing multi-label classification are introduced. In contrast to the
multi-class classification task, an unknown motion can have many labels that assign it to multiple classes.

To avoid confusion and to emphasize the distinction between multi-class (where one label assigns each motion one class)
and multi-label (where multiple labels assign each motion multiple classes), multi-class classification
is referred to as single-label classification hereinafter.

§ MOTION RECOGNITION

Hidden Markov Models are a popular choice for encoding human whole-body motions [60, 35, 37].
This section describes some properties of HMMs in depth and discusses properties and problems
that are especially relevant when dealing with motions.

§.§ Emission Distribution

Recall that the emission distribution models the observed data. Since in this case motions are observed, and all previously discussed features
are continuous, the emission distribution must also be continuous. Typically, a Gaussian distribution
or a mixture model thereof is used to model this case [37]. Since this work uses a
multivariate Gaussian distribution, the following discussion focuses on this distribution.

A multivariate Gaussian or normal distribution is defined by two parameters, its mean vector $\vecsym{\mu} \in \mathbb{R}^D$ and its
covariance matrix $\vecsym{\Sigma} \in \mathbb{R}^{D\times D}$, where $D$ is the dimension of the feature vector.
The probability density function (pdf) is then given by

\begin{equation}
f(\vec{x}) = \abs{\vecsym{\Sigma}}^{-\frac{1}{2}} (2 \pi)^{-\frac{D}{2}} \exp\left(-\frac{1}{2} (\vec{x}-\vecsym{\mu})^T \vecsym{\Sigma}^{-1} (\vec{x}-\vecsym{\mu})\right).
\end{equation}

In an HMM, the emission of each state $k$ is governed by its mean vector $\vecsym{\mu}_k$ and its covariance matrix
$\vecsym{\Sigma}_k$. When dealing with motions, the covariance matrices are often constrained to be diagonal to avoid numerical problems [35, 36].

§.§ Topologies

An important property of a Hidden Markov Models is that it uses hidden states. The transition probabilities
between the $K$ states are given by the transition matrix $\vec{A} \in [0,1]^{K \times K}$, while $\vecsym{\pi} \in [0,1]^K$
defines the start probabilities for each state. By constraining the transition matrix (and as a result the start probabilities),
different topologies can be realized. This is easily done by initializing the transition matrix and the
start probabilities with some entries set to zero. During training, all probabilities that were initially
set to zero will remain at zero [47].

Illustration of a $4$-state ergodic HMM (left) and a $4$-state left-to-right HMM with with $\Delta = 2$ (right) [47].

If the transition matrix is not constrained, a transition from any given state to every other state can occur.
Such an HMM is usually referred to as fully connected or ergodic. Another popular topology is
the left-to-right topology or Bakis topology [5]. In such a topology,
the states are thought to be aligned sequentially from left to right. At each state, only a transition to a
state that is right of the current state or
a self-transition is allowed. The model can thus only be traversed from left to right, hence the name.
In a left-to-right model, the start probabilities are set to $\vecsym{\pi} = (1, 0, \ldots, 0)$ while
the transition matrix takes the following shape:

\begin{equation}
\vec{A} =
\begin{bmatrix}
a_{1,1} & a_{1,2} & \dots  & a_{1,(K-1)}     & a_{1,K} \\
0 	    & a_{2,2} & \dots  & a_{2,(K-1)}     & a_{2,K} \\
\vdots  & \vdots  & \ddots & \vdots           & \vdots \\
0	    & 0		  & \dots  & a_{(K-1),(K-1)} & a_{(K-1),K} \\
0 	    & 0       & \dots  & 0               & a_{K,K} \\
\end{bmatrix}.
\end{equation}

The left-to-right constraint can thus be written as:

\begin{equation}
a_{i,j} = 0, \qquad j < i.
\end{equation}

To avoid skipping too many states while traversing from left to right, an additional constraint is introduced:

\begin{equation}
a_{i,j} = 0, \qquad j > i + \Delta.
\end{equation}

$\Delta$ limits the maximum number of allowed skips [47].
The left-to-right topology is frequently used for the recognition of motions [60, 35].
Both, the ergodic and the left-to-right topology are visualized in figure <ref>. Note
that, due to the variate of possible constraints for the transition matrix, other topologies are possible which are not discussed here.

§.§ Parameter Initialization

An interesting problem that arises is how to initialize the values of the transition matrix and the start probabilities as well as the
means and covariance matrices of the emission distributions. Since the Baum-Welch algorithm
does not necessarily converge to a global maximum but rather to a local one (see chapter <ref>), a proper initial estimate is important
to increase the chances of finding the global maximum during training.
According to [47], the start probabilities and transition matrix can either be initialized by a random
(while maintaining the constraint that the respective probabilities must sum to one)
or uniform estimation. Note that the constraints imposed by the topology choice must also be maintained.

The initialization of the mean vectors and covariance matrices of the emission distribution is more complicated. While a random
initialization is possible, it is beneficial to perform an initial estimate of the underlaying distribution.
A popular approach works as follows: The basic idea is to find $K$ clusters that correspond to the $K$ states of the Hidden Markov Model.
The mean  $\vecsym{\mu}_k$ and covariance matrix $\vecsym{\Sigma}_k$ for each state $k$ are then estimated over the respective
$k$-th cluster [47].
The necessary clustering can be performed by hand. Alternatively an unsupervised clustering algorithm like the k-means algorithm
can be used to automatically cluster the data. The k-means algorithm works by alternately assigning data samples
to a cluster given the current parameters and then updating the cluster center means such that the distance from the previously
associated samples is minimized. Lloyd's algorithm is often used to solve the clustering
problem efficiently [33].

§.§ Training and Recognition

After deciding on the hyperparameters of the model (namely the number of states, the topology
and the initialization method for the parameters of the model), the model can be trained. The training of an
HMM is performed efficiently using the Baum-Welch algorithm [12]: For each training sequence
$\vec{O}$ the parameters of the model are updated by first computing the expected likelihood given the current parameters (expectation step)
and then updating the parameters such that the expected quantity from the expectation step is maximized (maximization step). When performed iteratively, the likelihood
is maximized. The E and M steps are repeated until convergence or until a fixed number of iterations
have been performed.

A common problem when using the Baum-Welch algorithm are numerical instabilities. This is due
to the fact that the probabilities during the forward or backward pass can become extremely small since they
are multiplied at each time step. This becomes worse as the sequences grow longer. Since motions are usually
at least a couple of hundred samples long, this becomes a very real problem when training HMMs on motions.
The problem is mitigated by either using a scaling technique [47]
or by adapting the algorithm such that it uses the logarithm of the probabilities instead [38].

Another important property is that the training is unsupervised. This means that no target value
like a label is necessary to learn the parameters of the HMM. However, if multiple HMMs are used to classify
motions into classes, supervised learning becomes important. This will be discussed in the next section.
Finally, HMMs can be trained off-line and on-line. In off-line training, the HMM is trained
only once and the parameters are kept fixed even if previously unseen observation sequences become available. Most
discussion in the literature assume off-line learning, e.g. [47].
In on-line learning on the other hand, the HMMs are trained incrementally as new data becomes available.
Kulić et al. [34] describe such a system. This work only considers the off-line approach.

After the model has been trained, recognition is performed by calculating the likelihood under the
model $\lambda$ for an unknown observation sequence $\vec{O}$:

\begin{equation}
p(\vec{O} \mid \lambda).
\end{equation}

This is done efficiently by the forward algorithm [47]. Since the forward algorithm is used during training
as well, the same underflow issues as discussed earlier apply. Notice that the likelihood can become larger
than $1$ by definition, hence $0 \leq p(\vec{O} \mid \lambda) < \infty$. Since the likelihood can become both very small (that is, very close to zero but not negative)
and very large, the logarithmic likelihood (loglikelihood) is usually computed and presented.
A strongly positive loglikelihood thus indicates a motion that has been strongly recognized by the model whereas
a strongly negative loglikelihood indicates that the motion has not been recognized by the model at at. Finding
such a decision boundary will be discussed later in this chapter.

§.§ Extension to Factorial Hidden Markov Models

The previously discussed concepts apply equally to Factorial Hidden Markov Models. However, three additional
considerations are of interest: the number of Markov chains as an additional hyperparameter, efficient training of the FHMM
and computation of the likelihood under the model.

Firstly, the number of chains is an important hyperparameter since it directly controls the complexity
of time series that an FHMM can represent. However it also comes at the cost of increasing the computational
complexity of both the training and the evaluation given an unknown observation sequence (see chapter <ref>).
In the literature, few chains are usually used. For example, Kulić et al. [37] use
only two chains in their work. This does not seem like a lot at first. However, assume that an FHMM with
$15$ states and $2$ chains is used. This results in $15^2 = 225$ possible state combinations. Compared to
an HMM with $15$ states, even such a relatively small FHMM can already represent vastly more
history than its HMM counterpart.

Exact training of FHMMs is computationally very expensive and can be, depending
on the parameters, even infeasible. Four different approximations were already briefly discussed in chapter <ref>.
However, Kulić et al.[37] proposed a fifth approach that the authors used to train FHMMs
on motion data. The authors further showed that their approach is at least as good as and in most cases
even better than the exact inference algorithm when using it to train FHMMs on motion data. This makes
their approach especially interesting for this work. The algorithm devised by the authors works as follows:
The FHMM is trained sequentially. This means that each FHMM starts with a single chain. Inference
is then performed using the standard Baum-Welch algorithm on the training data. For the next chain $m$,
the residual error between the previously trained chains and the $n$-th training sample $\vec{O}^{(n)} \in \mathbb{R}^{T \times D}$ is computed for
each time step $t \in \{1,\ldots,T\}$:

\begin{equation}
\vec{e}_t^{(n)} = \frac{1}{W} \left(\vec{o}_t^{(n)} - \sum\limits_{i=1}^{m-1}{W \vec{c}_t^{(i)}}\right) \in \mathbb{R}^D,
\end{equation}

where $\vec{e}_t^{(n)}$ is the residual error between the frame at time step $t$ of the $n$-th training
sample and the summed contributions of the previous $m-1$ chains. Each chain's contribution
is weighted by $W = 1/M$, where $M$ is the number of all chains.
Finally, the contribution for each chain $m$ at time step $t$ is computed as follows:

\begin{equation}
\vec{c}_t^{(m)} = \sum\limits_{k=1}^{K}{\vecsym{\mu}_k^{(m)}\gamma_{t,k}^{(m)}} \in \mathbb{R}^D,
\end{equation}

where $K$ denotes the number of states, $\vecsym{\mu}_k^{(m)}$ is the $D$-dimensional mean vector of the
emission distribution of the already trained
chain $m$ in state $k$. Furthermore,  $\gamma_{t,k}^{(m)}$ denotes the probability that state $k$ in chain $m$ is active
at time step $t$. Algorithm <ref> summarizes the described training procedure.

initialize first chain

train first chain on time series $\vec{O}^{(n)} = (\vec{o}_t^{(n)})$ using the Baum-Welch algorithm

$m = 2,\ldots,M$
initialize next chain $m$

compute residual errors $\vec{e}_1^{(n)}, \ldots, \vec{e}_T^{(n)}$

train chain $m$ on time series $\vec{E}^{(n)} = (\vec{e}_t^{(n)})$ using the Baum-Welch algorithm

The sequential training algorithm for FHMMs in pseudo code [37].

A noteworthy and convenient property of the sequential training algorithm is that it uses procedures that are
already used when training HMMs. More concretely, the training of each chain is done by an unmodified
version of the well-known Baum-Welch algorithm. Computing $\gamma_{t,k}^{(m)}$ for the residual error
is achieved just as easily by using the standard forward-backward algorithm.

Lastly, the likelihoods under the FHMM must be calculated as well. This is done using the
exact algorithm given in [23]. In the more concrete case of the sequential training algorithm,
the necessary means and covariances are computed as follows:

\begin{equation}
\vecsym{\mu}_{k_1,\dots,k_M}    = W \sum\limits_{m=1}^{M}{\vecsym{\mu}_{k_m}^{(m)}} \qquad\text{and}\qquad
\vecsym{\Sigma}_{k_1,\dots,k_M} = W^2 \sum\limits_{m=1}^{M}{\vecsym{\Sigma}_{k_m}^{(m)}}.
\end{equation}

Notice that this must be done for each possible combination of states across all $M$ chains, as indicated
by the index $k_1,\ldots,k_M$. Take, for example, an FHMM with $3$ states and $2$ chains. This would result in $3^2 = 9$ combinations, with the following mean vectors:

\begin{align}
\vecsym{\mu}_{1,1} =& W \left(\vecsym{\mu}_1^{(1)} + \vecsym{\mu}_1^{(2)} \right), &
\vecsym{\mu}_{1,2} =& W \left(\vecsym{\mu}_1^{(1)} + \vecsym{\mu}_2^{(2)} \right), &
\vecsym{\mu}_{1,3} =& W \left(\vecsym{\mu}_1^{(1)} + \vecsym{\mu}_3^{(2)} \right), \\
\vecsym{\mu}_{2,1} =& W \left(\vecsym{\mu}_2^{(1)} + \vecsym{\mu}_1^{(2)} \right), &
\vecsym{\mu}_{2,2} =& W \left(\vecsym{\mu}_2^{(1)} + \vecsym{\mu}_2^{(2)} \right), &
\vecsym{\mu}_{2,3} =& W \left(\vecsym{\mu}_2^{(1)} + \vecsym{\mu}_3^{(2)} \right), \\
\vecsym{\mu}_{3,1} =& W \left(\vecsym{\mu}_3^{(1)} + \vecsym{\mu}_1^{(2)} \right), &
\vecsym{\mu}_{3,2} =& W \left(\vecsym{\mu}_3^{(1)} + \vecsym{\mu}_2^{(2)} \right), &
\vecsym{\mu}_{3,3} =& W \left(\vecsym{\mu}_3^{(1)} + \vecsym{\mu}_3^{(2)} \right).
\end{align}

The $9$ covariance matrices would be calculated analogously. The emission distribution for each
state combination is then given by a multivariate Gaussian distribution with mean $\vecsym{\mu}_{k_1,\dots,k_M}$ and
covariance $\vecsym{\Sigma}_{k_1,\dots,k_M}$ [30, 37].

§ SINGLE-LABEL CLASSIFICATION

In the previous chapter motion recognition was discussed. This discussion is now extended to a
classification problem. In a classification problem an unknown observation sequence, in this case a motion,
must be assigned to a finite set of classes $\mathcal{L}$. For example, assume that the set of classes
is run, jump and kick.
The goal is then to find the correct label $y \in \mathcal{L} = \{\text{run}, \text{jump}, \text{kick}\}$ for an
unknown motion $\vec{O}$ by making a prediction $p$. The classification is correct if $y = p$. Figure <ref>
provides an overview of the classification process. This process is used throughout this work, not just
for single-label classification.

[->, node distance=2cm, thick]
motion = [rectangle, rounded corners, minimum width=2cm, minimum height=1cm,text centered, draw=red!100, fill=red!25]
hmm = [rectangle, minimum width=2cm, minimum height=1cm, text centered, draw=green!100, fill=green!25]
ll = [rectangle, rounded corners, minimum width=2cm, minimum height=1cm,text centered, draw=purple!100, fill=purple!25]
dm = [rectangle, minimum width=1cm, minimum height=7.3cm, text centered, draw=orange!100, fill=orange!25]
pred = [rectangle, rounded corners, minimum width=2cm, minimum height=1cm,text centered, draw=blue!100, fill=blue!25]
(motion) [motion] $\vec{O}$;

(hmm_dot) [right=2cm of motion] $\ldots$;
(hmm_2) [hmm, above=0.5cm of hmm_dot] $\lambda_2$;
(hmm_1) [hmm, above of=hmm_2] $\lambda_1$;
(hmm_M1) [hmm, below=0.5cm of hmm_dot] $\lambda_{M-1}$;
(hmm_M) [hmm, below of=hmm_M1] $\lambda_M$;

(ll_dot) [right=3cm of hmm_dot] $\ldots$;
(ll_2) [ll, above=0.5cm of ll_dot] $p(\vec{O} \mid \lambda_2)$;
(ll_1) [ll, above of=ll_2] $p(\vec{O} \mid \lambda_1)$;
(ll_M1) [ll, below=0.5cm of ll_dot] $p(\vec{O} \mid \lambda_{M-1})$;
(ll_M) [ll, below of=ll_M1] $p(\vec{O} \mid \lambda_M)$;

(dm) [dm, right=2cm of ll_dot] 90Decision Maker;

(pred) [pred, right=1.3cm of dm] prediction;

(motion) |- (hmm_1);
(motion) |- (hmm_2);
(motion) |- (hmm_M1);
(motion) |- (hmm_M);

(hmm_1) – (ll_1);
(hmm_2) – (ll_2);
(hmm_M1) – (ll_M1);
(hmm_M) – (ll_M);

(ll_1.east) – (dm.west|-ll_1.center);
(ll_2.east) – (dm.west|-ll_2.center);
(ll_M1.east) – (dm.west|-ll_M1.center);
(ll_M.east) – (dm.west|-ll_M.center);

(dm) – (pred);

Overview of the classification process. For an unknown observation sequence $\vec{O}$ (in red), the likelihood
under each HMM (in green) is computed. Each HMM corresponds to a single class. The likelihoods (in purple)
are then fed into a decision maker (in orange) that computes the prediction (in blue).

If only a single label is assigned to a motion, the classification is straightforward. Instead
of using words, $y$ is encoded by natural numbers where each number corresponds to a class:

\begin{equation}
y \in \{1, \ldots, M\},
\end{equation}

where $M$ is the number of classes that need to be recognized. The mapping between word and class can
be chosen arbitrarily as long as it is bijective. In the above example, a possible mapping would be
$\text{run} \rightarrow 1$, $\text{jump} \rightarrow 2$ and $\text{kick} \rightarrow 3$.

Additionally, assume that $M$ HMMs have been trained per class: $\lambda_1, \ldots, \lambda_M$.
The first HMM was only trained with motions of class $1$, the second one only with motions of class $2$ and so on.
In contrast to the learning algorithm described in the previous chapter, the labels of the training data
are now used to select the appropriate HMM during training. While the training of the HMM itself is still unsupervised,
the selection happens in supervised fashion.

After training, the prediction of an unknown motion $\vec{O}$ must be computed. To do so, the likelihoods
of the motion under each individual HMM are evaluated. The last step is then performed by
what this work refers to as a decision maker. The job of a decision maker is to find a mapping from the likelihoods of the HMMs
to the prediction $p$. In case of the single-label classification problem, the prediction is done by a
very simple decision maker that chooses the label that corresponds to the HMM with the maximum likelihood:

\begin{equation}
p = \argmax_{m \in \{1,\ldots,M\}}{p(\vec{O} \mid \lambda_m)}.
\label{eq:single-classification}
\end{equation}

§ MULTI-LABEL CLASSIFICATION

Let us now extend the classification problem to a problem where multiple labels can be assigned to a
motion. This can be motivated by recalling the Motion Description Tree (see chapter <ref>):
A motion is specified by the leave nodes of the tree, e.g.

* locomotion $\rightarrow$ bipedal $\rightarrow$ walk,

* speed $\rightarrow$ fast,

* direction $\rightarrow$ left,

* perturbation $\rightarrow$ result $\rightarrow$ failing, and

* perturbation $\rightarrow$ source $\rightarrow$ passive.

It should be obvious from this example that most motions cannot be adequately described by a single class.
However, if a motion can have multiple labels, the simple classification described by equation <ref>
does not work anymore since it only computes the single label with the maximum likelihood.

To work around this problem, two different approaches can be identified: problem transformation using
the power set method or the use of more advanced decision makers. In the latter case, the decision maker
needs to truly handle multi-label classification. This can either be done by using the binary relevance method
or algorithm modification [48].

§.§ Power Set Method

A first attempt to solve this problem is to simply treat each possible combination of
labels as a single class. More formally, if $\mathcal{L}$ is the set of all labels,
compute the power set as the substitute label set:

\begin{equation}
\hat{\mathcal{L}} = \mathcal{P}(\mathcal{L}) = \{\mathcal{U} \mid \mathcal{U} \subseteq \mathcal{L}\} \setminus \emptyset.
\end{equation}

If each element in $\hat{\mathcal{L}}$ is then mapped to a natural number, the multi-label problem
has been transformed into a single-label problem. This means that equation <ref>
can be used to predict the substitute label with the maximum likelihood. Since the substitute label
represents multiple labels internally, reversing the transformation after the classification step solves the multi-label problem
using the already discussed approach for single-label classification. In the literature,
this idea is usually referred to as the power-set method [10, 48].

The power set method, however, has several downsides. Firstly, the number
of possible label combinations is $2^M-1$, where $M$ is the number of classes. Recall that for each class a Hidden
Markov Model must be trained. However, having sufficient training data now becomes a problem. If,
for example, $200$ walking motions exist but of those $200$ motions only a single a has the labels walk
and fast whereas all others have the labels walk and slow, the HMM for the fast walking
motion would only be trained on a single training sample.
This is inefficient since the HMM would benefit from the samples of the similar slow walking motion.
Secondly, to perform classification, the likelihood under each model must be calculated.
It can be easily seen that, for the worst case, this requires $\mathcal{O}(2^M)$ computations of the loglikelihood
for each unknown motion,
which makes this approach infeasible for even moderately large $M$.

§.§ Binary Relevance Method

Instead of treating all possible combinations of labels as a single label, an approach where a motion can truly
have multiple labels is interesting. One way to achieve this is to encode the label as a binary

\begin{equation}
\vec{y} \in \{0, 1\}^M,
\end{equation}

where $y_m$ is set to $1$ if the $m$-th label is active and set to $0$ otherwise.
To give an explanatory example, assume that the classes walk, fast and slow exist.
The label of a fast walking motion would then be encoded as $\vec{y}=(1,1,0)$ whereas motion
where the subject walks slowly would would be encoded as $\vec{y}=(1,0,1)$.

The training then works as follows:
$M$ HMMs are initialized, one for each class. Given a training sequence $\vec{O}$, use it to train
all HMMs which correspond to the set classes in the sequence's label $\vec{y}$. If, for example,
$\vec{y}=(1,0,1)$, train $\lambda_1$ and $\lambda_3$ on the same sequence. An advantage over the previously described method
is that only $M$ instead of $2^M-1$ HMMs (worst case) must be kept in memory. Additionally, the utilization of the available training
data is vastly improved. In the above example, the HMM that corresponds to the walk class now
benefits from both, the slow and fast walks. In the previously described approach,
both would only be considered in isolation. However, this can also become a problem
if classes are too generic. Take for example motion of a throw performed with the left hand and
a kick with the left foot. If both would be associated with the same class left, it is unlikely
that common patterns can be learned properly. Instead, hand-left and foot-left can be used
to distinguish the two, which is a more reasonable choice in this case.

However, it is unclear how the decision maker can be realized in this case. The previous approach
where the class with the maximum likelihood was selected cannot be used anymore. One way to perform classification in this case is to
use some fixed value as the decision boundary. Let $\bar{x}_m$ be such a decision boundary for the $m$-th label.
The multi-label classification is then achieved by computing all likelihoods of the unknown motion $\vec{O}$ under each
HMM and selecting all that are equal to or greater than the respective decision boundary:

\begin{align}
p_m = \begin{cases}
1, & \text{if } p(\vec{O} \mid \lambda_m) \geq \bar{x}_m \\
0, & \text{otherwise}\end{cases}
\qquad m = 1,\ldots,M,
\end{align}

where $y_m$ denotes the $m$-th element of the prediction $\vec{p}$. This approach is
promising, but unfortunately it all depends on a good choices of $\bar{x}_m$.

Instead of trying to find good values for all $\bar{x}_m$ by hand, it would be preferable if the decision boundaries could be determined automatically.
Since the approach described in this work is supervised, the association between the likelihoods and the
corresponding class is actually known for the training dataset. It seems like a good idea to use this knowledge to learn a decision
boundary from the training data instead of finding it by hand. This has the advantage that the model easily adapts to the current
dataset at hand without the need to fine-tune the decision boundary manually. The problem can indeed be seen
as binary classification: For each class, the features for the binary classifier are the likelihoods which
should be classified into two half-spaces that correspond to the labels $0$ and $1$ respectively. Since in this case, the
feature space is $1$-dimensional (just a scalar likelihood under one specific model), the hyperplane that separates the two
half-spaces is a simple point and corresponds to the decision boundary $\bar{x}_m$ described previously.
This means that any binary classifier like Logistic
Regression [31] or Support Vector Machines [26]
can be used to learn a mapping from the likelihoods to $y_m$ for each class.

Instead of only considering the $m$-th likelihood, in a slightly more advanced version of this basic idea
the binary classifier considers all likelihoods for its binary decision. This approach is commonly referred to as the binary relevance method [48].
More concretely, multi-label motion classification is realized with this method as follows:

* Given the training data $\vec{O}^{(1)}, \ldots, \vec{O}^{(N)}$ and respective labels $\vec{y}^{(1)}, \ldots, \vec{y}^{(M)}$, calculate the likelihoods for each motion under each model:

\begin{equation}
\label{eq:likelihood-mat}
\vec{x}^{(n)} = (p(\vec{O}^{(n)} \mid \lambda_1), p(\vec{O}^{(n)} \mid \lambda_2), \ldots, p(\vec{O}^{(n)} \mid \lambda_M)),
\qquad n \in \{1,\ldots,N\}.
\end{equation}

* Train $M$ binary classifiers. For each classifier, the features are $\vec{x}^{(1)}, \ldots, \vec{x}^{(N)}$. However, the labels are varied across classifiers, with the $m$-th classifier using labels $y^{(1)}_m, \ldots, y^{(N)}_m$.

* Given an unknown motion $\vec{O}$, calculate the likelihoods $\vec{x}$ under each of the $M$ models analogously to
equation <ref>. Afterwards, perform $M$ predictions using the previously trained classifiers to
obtain the complete multi-label prediction:

\begin{equation}
\vec{p} = (h^{(1)}(\vec{x}), \ldots, h^{(M)}(\vec{x})),
\end{equation}

where $h^{(m)}(\vec{x}) \in \{0,1\}$ denotes a function that computes the prediction of the $m$-th classifier.

In this case, the decision maker uses multiple binary classifiers internally to learn a mapping from
the likelihoods of the HMMs to the multi-label prediction using the binary relevance method. It is important
to stress again that the ability to use any binary classifier for multi-label predictions is a big
advantage of this method since a wide and well-understood variety of such classifiers exist.

However, this approach also has a downside: Each class is considered
in isolation. This means that information is lost since the correlation between classes in the label
vector $\vec{y}$ carries information as well. More concretely, it is likely that certain label combinations correspond
to certain patterns in the likelihoods. Such patterns cannot be detected by the approach
described above, since each classifier only “sees” a single class [48].

§.§ Modified Algorithms

Some learning algorithms have been modified to support multi-label classification “out of the box”.
In such a case, the multi-label classification is straightforward: Similarly to the previous approaches,
the likelihoods under each model are calculated for each of the $N$ training motions and the likelihood
vectors $\vec{x}^{(1)},\ldots,\vec{x}^{(N)}$ are obtained (compare equation <ref>).
In contrast to the binary relevance method, the classifier can now be trained on the entire label
vector $\vec{y}$ instead of training individual classifiers on the individual classes. This potentially
allows the learning algorithm to find patterns between the likelihoods and the classes since the correlation
between classes can now also be considered. In this case the decision maker is thus simply a classifier
that is capable of learning multi-label classification.

Multiple algorithms exist that have been adopted to the multi-label problem. Decision Trees
have been extended to work in multi-label classification problems [64]. With this
extension, Random Forests [13] can also be used to perform multi-label
classification since a Random Forest is an ensemble classifier that uses multiple Decision Trees
internally. While this work focuses on the mentioned algorithms, it should be noted that other algorithms have been modified
as well to support multi-label predictions, e.g. AdaBoost [51] or $k$-Nearest Neighbors [70].

CHAPTER: EVALUATION

The evaluation of the previously discussed concepts is challenging. This is due to the fact
that a huge number of system parameters can be identified: A set of features must be selected,
the hyperparameters of the Hidden Markov Models must be chosen, a comparison between FHMMs and
HMMs is necessary and a variety of different classifiers that perform the mapping from likelihoods to final
labels are available which also potentially have many hyperparameters that need fine-tuning. It should
be clear that evaluating all those DoFs simultaneously is impossible due to the large number of
combinations. Instead, the problem is split into smaller problems that are evaluated and optimized
individually. Only a small subset of the best parameter choices is then used in the next problem, and so on.

The structure of this chapter is based on this idea: To get started, a brief overview of the tools used and
developed for this work is given. The next section describes the used evaluation dataset. For the first evaluation step, feature selection is performed to find
a “good” subset of the wide variety of available features (see chapter <ref>).
Next, Hidden Markov Models and their set of hyperparameters are evaluated (see chapter <ref>). The discussion is continued
by a comparison between Hidden Markov Models and Factorial Hidden Markov Models (see chapter <ref>).
Different approaches to predict the labels are evaluated (see chapter <ref>).
Finally, an end-to-end evaluation of two classifier systems is conducted using the results from the previous sections.

§ TOOLS

A toolkit was developed to help to evaluate this work. The toolkit is mostly written in Python with some
modules using C and C++ either for speed or for interoperability. Additionally, the excellent scikit-learn
framework [46] was used throughout the toolkit. The toolkit
is used via multiple command-line scripts. This allows to easily automate long-running tasks and handles
the wide variety of different actions and parameters efficiently.
Areas of responsibility are split into 3 main modules which are each briefly covered in the remainder of this section.

§.§  Module

The  module provides classes and functions to load motion from different file formats
(MMM and C3D files). Furthermore, the derived features described in chapter <ref> can be computed as well.
The functionality of the Simox[<http://simox.sourceforge.net>] and
MMMTools[<https://gitlab.com/mastermotormap/mmmtools>] libraries are used to aid in the computation of more advanced
features like the angular momentum. SWIG[<http://swig.org>] is used to bridge
between Python and the C++ API of the two mentioned libraries[<https://gitlab.com/cmandery/pySimox>][<https://gitlab.com/cmandery/pyMMM>]. The module also provides functionality
to normalize features to obtain rotation and translation invariant representations. Features can also be smoothed
by applying a moving average filter and transformed to have similar scales across features. In short,
the  module implements everything that was discussed in chapter <ref>.

Another important responsibility of the  module is the handling of datasets. A dataset
is defined by a simple manifest. The manifest is realized as a JSON file that contains
an array of folders. Each folder can be label with multiple classes to perform supervised learning.
A dataset can be loaded from such a JSON manifest file. Multiple file formats (MMM and C3D) can be merged
into a combined dataset. A dataset can be checked for inconsistencies (e.g. missing files), a report
can be generated (e.g. the number of samples and the labels present in the dataset), selected
features can be plotted to visualize the data and the entire dataset can be exported into a single file
(different file formats are supported). Especially the export functionality is important since the
computation of the derived features and the normalization is computationally expensive. All evaluation
steps therefore load such an exported dataset instead of computing the dataset from scratch every time.

Lastly, the  module also provides functionality that helps with splitting the dataset
into two separate (and disjunct) datasets: one for training, one for testing. This is done using
stratified $k$-fold for multi-label data as proposed by Sechidis et al. [52]. The basic idea of a $k$-fold
evaluation is to split the dataset into $k$ disjunct subsets or folds. $k$ rounds can then be performed by
using $k-1$ folds for training and the remaining fold for testing. In a stratified $k$-fold, each
fold attempts to have the same distribution of data across classes as the entire dataset.

§.§  Module

The  module implements everything necessary to train and evaluate Hidden Markov Models.
The basic HMM functionality is provided by the hmmlearn library[<https://github.com/hmmlearn/hmmlearn>].
Other libraries were considered as well, namely pomegranate[<https://github.com/jmschrei/pomegranate>] and
GHMM[<http://ghmm.org>]. However, hmmlearn proved to be the most reliable and fastest

For this evaluation, a fork[<https://github.com/matthiasplappert/hmmlearn>] of the hmmlearn library
was created and used. It fixes some issues and, more importantly, implements functionality to use
and train Factorial Hidden Markov Models. The implementation allows to compute exact likelihoods
as described by Ghahramani et al. [23]. Additionally, it supports the sequential
training algorithm developed by Kulić et al. [37] (see chapter <ref>)

The module also provides functionality to combine the individual HMMs into an ensemble of HMMs.
It allows training of the individual HMMs on labeled training data as well as computing the loglikelihoods
of unknown motions under each model (see chapter <ref> and <ref>). The training
and evaluation of likelihoods is performed in parallel for each model to maximize performance.

Lastly, the  module provides utility functions to compute the initial transition matrix
and start probabilities using different topologies (see chapter <ref>). Functions to estimate the initial means and covariance matrices
using different strategies and constraints as described in chapter <ref> are included as well.

§.§  Module

The  module is a collection of smaller components. The most important one are the
decision makers. A decision maker is a classifier that performs the mapping from the likelihoods
as calculated by the HMM ensemble to the binary label vector. As described in chapter <ref>,
different approaches exist. The  module implements a decision maker that always picks
the maximum likelihood (chapter <ref>) and one that uses a fixed decision boundary
(chapter <ref>). Additionally, decision makers
that use Logistic Regression and Support Vector Machines are implemented using the binary relevance method (chapter <ref>).
Decision Trees and Random Forests are used with adapted learning algorithms to handle the multi-label problem directly
(chapter <ref>). The decision makers make especially heavy use of the scikit-learn framework.

§ DATASET

Each motion in the dataset was recorded using the motion capture system described in chapter <ref>. The raw positions of the
motion markers were stored in a C3D file. Each motion was recorded using the KIT reference marker set consistent of $56$
markers placed at well-defined anatomical locations and a sampling rate
of 100. The C3D files are converted to the MMM reference model using $40$ joints (each with a single DoF) and a non-linear optimization
algorithm (see chapter <ref>). The joint angles as well as the root position and root rotation were stored in an XML file.
The C3D and XML files are used to extract and compute all features described in chapter <ref>.
For all features, the normalized or unnormalized form can be used. During normalization, the root rotation
of the first frame was constrained such that the roll and pitch angles were set to zero. This is because
a small error in either of the two angles can cause an improperly rotated coordinate system. The features can also be smoothed using a moving average filter with $W=3$
(compare chapter <ref>). Table <ref> lists all $29$ available features without differentiating
between normalized/unnormalized and smoothed/not smoothed.

Dimension 		       Feature Name     Description
40     joint_pos	  		   				    angles of all $40$ joints

40     joint_vel	    	  				        velocities of each joint

1     joint_vel_norm 	   				        Euclidean norm of the joint velocities

40     joint_acc 			   			        acceleration of each joint

1     joint_acc_norm 	   				        combined Euclidean norm of the joint accelerations

3     root_pos 			    					    root position of the subject in Cartesian space

3     root_vel 			    					    directed root velocity

1     root_vel_norm 	 					    Euclidean norm of the root velocity

3     root_acc 			   					    directed root acceleration

1     root_acc_norm 	 					    Euclidean norm of the directed root acceleration

3     root_rot 			   					    pitch, roll and yawn angles of the subject's root

1     root_rot_norm 	 					    Euclidean norm of the root rotation

12     extremities_pos  				    position of the hands and feet in Cartesian space

12     extremities_vel	 				    directed velocities of the hands and feet

4     extremities_vel_norm  					    Euclidean norm of the directed extremity velocities per hand/foot

12     extremities_acc		    				    directed accelerations of the hands and feet

4     extremities_acc_norm     				    Euclidean norm of the directed extremity velocities per hand/foot

3     com_pos 			      				    position of the center of mass in Cartesian space

3     com_vel 			      				    directed velocity of the center of mass

1     com_vel_norm 	   				    Euclidean norm of the CoM's velocity

3     com_acc 			  	 				    directed acceleration of the center of mass

1     com_acc_norm 		      				    Euclidean norm of the CoM's acceleration

3     angular_momentum 	    				    whole-body angular momentum in $x$, $y$ and $z$ direction

1     angular_momentum_norm   				    Euclidean norm of the whole-body angular momentum

168     marker_pos 			    				    position of all $56$ markers in Cartesian space

168     marker_vel 			    				    directed velocities of all markers

1     marker_vel_norm 	      				    Euclidean norm of all markers' velocities

168     marker_acc 			    				    directed acceleration of all markers

1     marker_acc_norm 	     				    Euclidean norm of all markers' accelerations

List of all $29$ available features with a total of $702$ dimensions.

A total of $454$ motions were selected from the motion database. $10$ different human subjects performed
the motions. Of those $10$ subjects, $3$ are female and $7$ are male. Each motion was performed by a single
human subject. If a motion involves an object (e.g. playing the guitar), the object was not physically present
but instead imagined by the subject. The only exception to this is a set of upwards walking motions where
actual stairs were used since it is quite hard to climb imaginary ones. Each recording in the dataset
contains only a single motion. Additionally, each subject begins the motion in an upward-standing pose similar to the zero-pose of the MMM model
and finishes the motion in the same pose. For periodic motions such as walking, stirring or dancing,
a fixed number of repetitions were selected across the dataset. The motions in the dataset were manually labeled
with $49$ different labels. Table <ref> lists all labels and the number of samples in the dataset that are assigned to each label.
Since motions can have multiple labels, table <ref> lists all label combinations
that are present in the evaluation dataset.

2cMotions     2cAttributes
(r)1-2 (l)3-4
Samples     Label     Samples     Label
180     walk 		       31     turn-right

35     push-recovery     28     turn-left

41     run 		       9     speed-fast

28     kick             10     speed-normal

10     throw 	      10     speed-slow

10     bow            84     direction-forward

5     squat 		      14     direction-backward

10     punch 	      15     direction-left

10     stomp           12     direction-right

25     jump            15     direction-upward

11     golf            35     direction-circle

20     tennis          27     direction-slalom

15     wave            10     bend-right

11     play-guitar     10     bend-left

10     play-violin     18     counter-clockwise

11     stir            17     clockwise

11     wipe            23     foot-right

11     dance           15     foot-left

45     hand-right

48     hand-left

5     hand-both

5     deep

5     slight

10     high

10     low

5     putting

6     drive

10     smash

10     forehand

6     waltz

5     chachacha

List of all $49$ different labels used in the dataset.

Samples     Label Combination     Samples     Label Combination
(r)1-2 (l)3-4
10     walk, direction-forward, speed-normal     5     stomp, foot-right

10     walk, direction-forward, speed-slow     5     stomp, foot-left

9      walk, direction-forward, speed-fast     5     jump, direction-forward

31     walk, turn-right     5     jump, direction-backward

28     walk, turn-left     5     jump, direction-right

10     walk, bend-left     5     jump, direction-upward

10     walk, bend-right     5     jump, direction-left

18     walk, direction-circle, counter-clockwise     5     golf, putting

17     walk, direction-circle, clockwise     6     golf, drive

27     walk, direction-slalom      5     tennis, hand-right, smash

10     walk, direction-upward      5     tennis, hand-right, forehand

9     push-recovery, direction-forward     5     tennis, hand-left, smash

9     push-recovery, direction-backward     5     tennis, hand-left, forehand

7     push-recovery, direction-right     5     wave, hand-right

10     push-recovery, direction-left     5     wave, hand-left

41     run, direction-forward      5     wave, hand-both

8     kick, foot-right      5     play-guitar, hand-right

5     kick, foot-right, high     6     play-guitar, hand-left

5     kick, foot-right, low      5     play-violin, hand-right

5     kick, foot-left, high      5     play-violin, hand-left

5     kick, foot-left, low      5     stir, hand-right

5     throw, hand-right      6     stir, hand-left

5     throw, hand-left      5     wipe, hand-right

5     bow, deep      6     wipe, hand-left

5     bow, slight      6     dance, waltz

5     squat      5     dance, chachacha

5     punch, hand-right

5     punch, hand-left

All $54$ label combinations used in the dataset.

§ FEATURE SELECTION

In feature selection, the goal is to find a subset of features that are relevant to the problem at hand.
The selection process is a crucial first step. This is due to the fact that any classifier can only compute
predictions from the features it receives. If those features do not contain the relevant features, the classifier
obviously cannot produce optimal results. Additionally, a reduction of the feature dimensionality has
several advantages: Firstly, the computational complexity is usually proportional to the dimensionality
of the feature space. Secondly, if the feature space has very high dimensionality but the dataset is
relatively small, overfitting becomes more likely. Lastly, reducing the feature
set can also bring interesting insights to better understand a problem at hand [22].

Unfortunately, the optimal feature set cannot be computed directly by evaluating all possible subsets.
This is due to the fact that, for a feature set with $N$ different features, $2^N-1$ non-empty subsets exist.
As previously described, a set of $N=29$ different features can be identified in this work with
results in $2^{29}-1 = 536\,870\,911$ subsets. Evaluating each of these subsets is infeasible.
In the literature, three popular heuristics can be identified to solve this problem:
the filter method, the wrapper method and the embedded method.
The filter method selects features by ranking them with correlation coefficients. The wrapper method
evaluates the usefulness of a subset using a classifier. The embedded method works similarly to the
wrapper method in the sense that it uses a predictor. However, while the wrapper methods relies
on an external measure to assess the quality of the feature set, the embedded method relies on feature
selection that is inherent to the learning algorithms. To give a concrete example, Support Vector Machines
with L1 regularization can be used to perform feature selection using the embedded method. This is possible
because L1 regularization yields sparse weights. The selection can then be done
by simply picking all features that correspond to non-zero weights without the need of an externally defined measure [22].

In this work, feature selection is performed using the wrapper method with backwards elimination.
In backwards elimination, one starts with the full feature set. Each feature is temporarily excluded from the set once and a model
is trained on each of the subsets. For each subset, the model is evaluated on the test set and a measure
is computed. The feature with the least effect on the measure is then removed from the feature set. This process
is repeated until the feature set is empty (or until a stopping criterion is reached).
The feature selection was performed on the dataset described in section <ref>. Of the $29$ available features, the features
marker_pos, marker_vel and marker_acc were initially excluded. This is due
to the very high dimensionality of the three mentioned features ($168$ dimensions each) that resulted in
numerical problems. The three excluded features will be revisited later. Additionally, to avoid too much complexity, all features
under consideration were normalized, smoothed (using the previously described moving average with $W=3$) and scaled to be
in range $[-1, 1]$. The reasoning behind
this decision is that it seems very unlikely that unnormalized features perform well since they are,
as previously discussed, neither invariant to translation nor rotation. Another concern was that
the unnormalized features can cause overfitting. The smoothed features were
selected as the default for similar reasons since it seems unlikely that an error-prone signal with strong jitter
performs better than a slightly smoothed representation thereof. Scaling is important to ensure a good
initialization when using $k$-means clustering. Since the dataset contains $49$ different
classes, $49$ HMM models with $5$ states each and the left-to-right topology with $\Delta = 1$ were
used. Each model was trained using the Baum-Welch algorithm for $10$ iterations. The transition matrix
and start probabilities were initialized uniformly whereas $k$-means clustering was used to initialize the
means and covariance matrices. The covariance matrix was further constrained to be diagonal and the diagonal
elements were constrained to be larger than $0.0001$ to avoid numerical instabilities. The entire dataset was
initially shuffled to break correlations between nearby training samples. A stratified $3$-fold was used to
train and evaluate the HMMs three times per feature subset. In each round, the likelihoods under each model were
calculated for the test split. The likelihoods were additionally split into a set of positive and negative likelihoods.
The positive set contains all likelihoods of the motions that should have been recognized by the respective model whereas the
negative set contains the likelihoods of the motions that should have been rejected. Since $49$ different HMMs are used this results in
$98$ sets of likelihoods: $49$ sets of positive instances and $49$ sets of negative instances. The mean and standard deviation
were then calculated per set over all three rounds, resulting in $98$ means of the likelihoods and their respective standard deviations.

A possible way to measure the results in each round is to compute the distance between the distribution of positive and negative likelihoods
on a per-class basis. The Wasserstein metric [24, 65]
can be used to compute the distance between two Gaussian distributions:

\begin{equation}
m_i = \sqrt{|\bar{\mu}_{pos_i} - \bar{\mu}_{neg_i}| + \left(\bar{\sigma}_{pos_i}^2 + \bar{\sigma}_{neg_i}^2 - 2 \sqrt{\bar{\sigma}_{pos_i}^2 \bar{\sigma}_{neg_i}^2}\right)},
\end{equation}

where $i$ denotes the $i$-th class, $\bar{\mu}_{pos_i}$ is the mean over all positive
and $\bar{\mu}_{neg_i}$ the mean over all negative likelihoods for the $i$-th class. $\bar{\sigma}_{pos_i}$ and
$\bar{\sigma}_{neg_i}$ denote the respective standard deviations. It was further assumed that the positive and negative likelihood distributions are
Gaussian. This assumption is not formally proven here but seems likely to hold given the central limit theorem.

After computing $m_1, \ldots, m_{49}$,
the individual distance measures must be combined into a single measure. Since all
labels are equally important, the median is used here. This ensures that
the distance between positive and negative distribution is balanced across all $49$ classes. The total dimension of the features
should be considered as well. If, for example, roughly the same distance has been computed for a feature set with $100$ and
$10$ dimensions respectively, the feature set with only $10$ dimensions should be preferred. A simple way to achieve
this is to compute a measure per feature. Based on these considerations, the following combined measure is devised:

\begin{equation}
m = \frac{1}{D} \tilde{m},
\label{eq:wasserstein-metric}
\end{equation}

where $D$ denotes current number of feature dimensions and $\tilde{m}$ is the median over individual
Wasserstein metrics $m_1, \ldots, m_{49}$. Other measures
were considered (e.g. AICc and the Mahalanobis metric) but proofed to not work as well
as the Wasserstein-based metric.

Round     Score     Dimension     Deleted Feature
1       5 899      198     —

2       6 615      158     joint_acc

3       9 009      118     joint_vel

4      10 686       78     joint_pos

5      12 602       66     extremities_acc

6      15 508       54     extremities_vel

7      17 176       53     com_vel_norm

8      18 740       49     extremities_vel_norm

9      20 856       46     angular_momentum

10      23 644       42     extremities_acc_norm

11      26 253       41     root_vel_norm

12      32 415       38     com_vel

13      35 563       35     com_acc

14      38 544       32     root_acc

15      40 214       31     joint_acc_norm

16      42 362       30     joint_vel_norm

17      44 416       29     com_acc_norm

18      46 866       28     root_acc_norm

19      48 849       27     marker_acc_norm

20      50 754       26     angular_momentum_norm

21      52 949       25     marker_vel_norm

22      54 781       22     com_pos

23      51 928       19     root_pos

24      54 525       16     root_vel

25      48 524       13     root_rot

26     127 598        1     extremities_pos
3cRemaining Feature        root_rot_norm

The results of the feature selection using backward elimination. The table lists all elimination
rounds, the feature that was deleted in each round, the measure that was computed in each round without the deleted
feature and the dimension of the combined features without the deleted feature.

Table <ref> lists the results per elimination round. The best feature was found in the last round: root_rot_norm. However, the individual
likelihoods reveal that the single feature does not capture enough information to recognize all motion
types. For example, the stir motion had a mean loglikelihood of $-12\,046 \pm 43\,080$, which
clearly indicates that this motion cannot be recognized using only the normalized norm of the root rotation.
Some direction labels achieve similarly bad results. The label clockwise and counter-clockwise
for example achieved a loglikelihood score of $-311 \pm 67$ and $-275 \pm 103$ respectively. However, the
mean loglikelihoods for motions that should not be recognized by the models are $401 \pm 701$
and $405 \pm 725$ respectively. This feature set was therefore discarded as an outlier.

The second best result was achieved in round $22$ with the following set of features: root_pos,
root_vel, extremities_pos, root_rot and root_rot_norm.
This is a more reasonable result. The velocity of the subject's root certainly is important
to detect the direction the subject moves in and to decide if the motion is of stationary or dynamic nature.
The position of the extremities proved to be an especially important feature since it survived all elimination
rounds until the last one. This also makes sense: A lot of motions in the dataset involve the subject's hands
or feet. Examples are stirring, waving, throwing a ball, playing the guitar or violin and kicks. The
position of the extremities presumably also help with identifying which hand or foot was used to perform
the motion. Including the root rotation and its norm makes sense as well: Since the root position
and the root rotation are both included, it should be possible to adequately determine the subject's position
and rotation in space from the feature set. This is obviously important information. An interesting
observation is that the joint angles are not that important to recognize motions. The joint angles
were already discarded in the fourth round. This is presumably due to the relatively low
information contents per feature. While some joints might carry important information, this
information is also redundantly present in other features like the position of the extremities. Since
other feature sets achieved similar scores, they will be revisited in the last
section during the end-to-end evaluation.

All things considered, the selected features seem like a good choice. This claim is backed by the data in table <ref> (at
the very end of this section)
which lists the mean likelihoods and their standard deviation per label over the entire test dataset
for HMMs that were trained using the previously mentioned feature set from elimination round $22$.
The likelihoods of positive samples are consistently high with a reasonable standard deviation. Additionally, the models
only respond to the motions that they were trained to recognize. The only exception is the model
that represents walking motions, which seems likely to produce false positives. This is presumably due
to the fact that very different walking motions (e.g. walking and making a 90 degree turn vs. walking
forward) were used to train the same model. The high standard deviation of the negative
likelihoods can be explained by the fact that an extremely wide variety of very different motions were
all combined into a single score. In other words: The positive likelihoods were only computed for motions
of the same type whereas the negative likelihoods combines motion of all other types into a single score.

Score     Dimension     Additional Feature
54 781            22                     —

4 456            181                    marker_pos

4 459            181                    marker_vel

NaN                181                    marker_acc

Comparison between the baseline feature set and the previously excluded features. Each of the previously excluded
features is added once to the baseline set and the performance is evaluated.

Recall that three features were excluded from the previous feature selection process: marker_pos,
marker_vel and marker_acc. This was necessary because of numerical problems when including these features in the initial feature set during backwards elimination.
More concretely, the probabilities during inference would become so small that an underflow occurred even though
the HMM implementation used during the evaluation already uses logarithmic probabilities to counter the problem. However,
these features are now revisited. To do so, the best feature set as previously determined was used as a baseline.
Each of the three left out features was then added once. The achieved scores (using the same measure
and training procedure as before) are listed in table <ref>.

It is quite obvious from the data that the marker positions, velocities and accelerations do not
contribute anything to the feature set; instead the scores worsen significantly. Notice that the marker
accelerations even cause numerical instabilities due to vanishing probabilities, hence a score cannot be computed in this case.
The three temporarily excluded features are thus excluded permanently from further consideration.

Class      Positive Loglikelihood      Negative Loglikelihood
walk                    $12\,858 \pm 4\,924$           $10\,095 \pm 16\,747$

push-recovery           $84\,073 \pm 26\,271$          $-830\,206 \pm 1\,314\,193$

run                     $13\,655 \pm 3\,937$           $-161\,784 \pm 344\,406$

kick                    $24\,898 \pm 5\,235$          $-207\,653 \pm 422\,938$

throw                   $25\,122 \pm 6\,796$          $-4\,982\,845 \pm 8\,584\,547$

bow                     $31\,251 \pm 4\,014$          $-5\,946\,674 \pm 9\,198\,336$

squat                   $32\,149 \pm 6\,640$          $-2\,198\,971 \pm 3\,600\,540$

punch                   $35\,073 \pm 11\,606$         $-1\,828\,241 \pm 3\,125\,022$

stomp                   $34\,982 \pm 3\,633$          $-3\,409\,349 \pm 5\,592\,487$

jump                    $24\,124 \pm 6\,256$          $-293\,847 \pm 612\,549$

golf                    $37\,559 \pm 8\,758$          $-1\,880\,819 \pm 3\,414\,369$

tennis                  $34\,408 \pm 3\,933$          $-303\,574 \pm 551\,566$

wave                    $58\,801 \pm 7\,035$          $-14\,784\,535 \pm 21\,892\,345$

play-guitar             $49\,155 \pm 4\,275$          $-3\,531\,256 \pm 6\,955\,093$

play-violin             $55\,510 \pm 3\,780$          $-6\,878\,584 \pm 10\,443\,505$

stir                    $28\,269 \pm 52\,856$         $-5\,574\,138 \pm 8\,338\,055$

wipe                    $48\,898 \pm 4\,734$          $-8\,223\,779 \pm 13\,225\,863$

dance                   $33\,981 \pm 5\,709$          $-178\,554 \pm 259\,631$

turn-right              $18\,467 \pm 3\,158$          $-58\,574 \pm 111\,112$

turn-left               $18\,594 \pm 7\,066$          $-36\,449 \pm 69\,685$

speed-fast              $20\,673 \pm 3\,640$          $-333\,766 \pm 531\,406$

speed-normal            $27\,002 \pm 8\,508$          $-393\,349 \pm 618\,222$

speed-slow              $28\,689 \pm 2\,372$          $-410\,014 \pm 559\,746$

direction-forward       $20\,472 \pm 17\,043$         $-159\,039 \pm 328\,429$

direction-backward      $54\,408 \pm 24\,864$         $-720\,136 \pm 1\,467\,702$

direction-left          $86\,904 \pm 52\,238$         $-949\,005 \pm 1\,478\,147$

direction-right         $50\,982 \pm 24\,431$         $-1\,271\,632 \pm 1\,879\,763$

direction-upward        $27\,477 \pm 4\,773$           $-194\,048 \pm 366\,744$

direction-circle        $17\,183 \pm 6\,479$           $-24\,067 \pm 75\,746$

direction-slalom        $23\,776 \pm 4\,992$           $-38\,560 \pm 106\,968$

bend-right              $28\,950 \pm 3\,030$          $-191\,518 \pm 237\,383$

bend-left               $23\,457 \pm 4\,346$          $-206\,328 \pm 262\,216$

counter-clockwise       $21\,281 \pm 5\,297$          $-50\,763 \pm 115\,171$

clockwise               $20\,271 \pm 8\,250$            $-45\,160 \pm 83\,697$

foot-right              $25\,384 \pm 5\,914$           $-229\,056 \pm 472\,367$

foot-left               $31\,525 \pm 2\,765$          $-1\,583\,832 \pm 2\,753\,615$

hand-right              $39\,043 \pm 11\,234$          $-457\,165 \pm 819\,820$

hand-left               $40\,720 \pm 11\,209$          $-380\,057 \pm 723\,547$

hand-both               $62\,305 \pm 4\,304$          $-8\,376\,056 \pm 12\,155\,675$

deep                    $29\,047 \pm 2\,347$          $-5\,173\,123 \pm 7\,942\,443$

slight                  $31\,561 \pm 6\,827$          $-3\,272\,445 \pm 6\,239\,772$

high                    $29\,758 \pm 4\,651$          $-1\,173\,850 \pm 2\,088\,698$

low                     $27\,388 \pm 3\,102$          $-1\,653\,420 \pm 2\,928\,772$

putting                 $39\,531 \pm 4\,220$          $-4\,056\,969 \pm 6\,825\,540$

drive                   $36\,932 \pm 11\,987$         $-1\,846\,573 \pm 3\,369\,606$

smash                   $35\,313 \pm 4\,018$           $-328\,343 \pm 665\,057$

forehand                $31\,963 \pm 4\,084$           $-826\,177 \pm 1\,417\,883$

waltz                   $36\,618 \pm 10\,136$          $-396\,325 \pm 634\,596$

chachacha               $29\,362 \pm 8\,173$           $-797\,359 \pm 1\,213\,353$

The mean positive and negative loglikelihoods under each of the $49$ classes for the feature set selected during round $22$.

§ HIDDEN MARKOV MODELS

As discussed in chapter <ref>, Hidden Markov Models can be used to perform motion
recognition. However, HMMs have multiple hyperparameters that need to be selected, namely the number
of states and the topology. Additionally, different initialization strategies can be used. As mentioned
in chapter <ref>, proper initialization is a crucial step before training.
Moreover, a comparison between the performance of HMMs and FHMMs makes sense. FHMMs have another
important hyperparameter, the number of chains.

Throughout this section, the experimental setup from the previous feature selection section was used again.
More concretely, the same dataset was used. To make the results comparable between this and the previous section,
the same permutation was used to shuffle the dataset. $3$-fold cross-validation was used with $49$ HMMs, one for each
class. Each HMM was trained for $10$ iterations. The best feature set from the previous section was used throughout
the following experiments: root_pos, root_vel, extremities_pos, root_rot and root_rot_norm.
All features were normalized, smoothed and scaled. The same Wasserstein-based metric (equation <ref>)
from the previous chapter was used to score the results. The default parameters of the HMMs were as follows: $5$ states, left-to-right topology
with $\Delta = 1$, uniform initialization of the transition & start probabilities and $k$-means initialization
of the Gaussian emission distribution parameters with the covariance matrices constrained to be diagonal. Notice that these are only defaults. Each experiment
will vary either one or multiple of these parameters, which will be clearly stated in each subsection.

§.§ Hyperparameters

The optimal number of states and topology are found using grid search. In grid search, each possible
combination of hyperparameters is used to train a model on the training dataset which is then evaluated
on the test dataset. A measure is computed per combination and the best combination is selected. Such a search
is feasible in this case since only two different hyperparameters are evaluated.

In theory any natural number can be used for the number of states and countless different topologies are possible.
In practice however, it makes sense to limit the number of states $K$ to be between $5$ and $20$ states
when recognizing human motions [35, 37]. In this evaluation $K~\in~\{3,\ldots,20\}$,
resulting in $18$ different values for the number of states.
The topologies (chapter <ref>) under consideration are the following:

* fully-connected,

* left-to-right without a $\Delta$ constraint,

* left-to-right with $\Delta = 1$, and

* left-to-right with $\Delta = 2$.

The left-to-right topologies proved to be especially well-suited for motion recognition [37]. Different
variations of this topology are evaluated by varying the $\Delta$ parameter. For completeness, the fully-connected topology is considered
as well.

Hence $18$ different values for the number of states and $4$ different topologies were under evaluation.
This results in a total of $18~\cdot~4~=~72$ combinations, making the grid search easily feasible.
The previously described experimental setup was used, with all parameters fixed and set to their default values except
for the number of states and the topology. The results are depicted in figure <ref>.

The results of the hyperparameter grid search. The four different topologies are plotted individually. For each
topology the number of states are plotted against the respectively achieved score.

As can be seen from the results, it is indeed desirable to have a low number of states. Good choices
are between $5$ and $8$ states, depending on the topology. Another interesting observation is that
the left-to-right topology with $\Delta = 1$ is the only topology that performed consistently
without large jumps. However, the best result was achieved by the left-to-right topology without
a $\Delta$ constraint and with $8$ states. The result support the claim that a low number of states
is preferable for motion recognition and that the left-to-right topology is a good choice when
dealing with motions. However, except for a few cases, the performance across topologies was similar.

§.§ Parameter Initialization

Different strategies can be used to initialize the parameters of a Hidden Markov Model before starting
the training (see chapter <ref>). To quickly recap, the transition
matrix and the start probabilities can either be initialized uniformly or by randomization. The implementation
of a uniform initialization is trivial. The randomization can be easily achieved as well: The transition
and start probabilities are initialized uniformly. Next, each element from the transition matrix is multiplied by a pseudo-random
number (in this case the random number was drawn uniformly from the interval $[0,1]$). Notice that for each
multiplication a new random number must be used. After randomization, the probabilities need to be normalized
such that they sum to $1$. This can be easily done by summing over each row of the transition matrix
and dividing each element in that row by the the respective sum. The randomization of the start probabilities vector works
analogously. This randomization automatically ensures that the desired topology is maintained.

The means and covariance matrices can either be randomized or estimated using the $k$-means clustering algorithm.
The randomization of the mean vectors is straightforward: $D$ pseudo-random numbers (in this case drawn uniformly from $[-1,1]$ to match
the feature scaling) are combined to form the mean vector ($D$ denotes the number of features). This is repeated for each state to obtain $K$ random mean vectors.
The randomization of the covariance matrix is a bit more complicated since a covariance matrix
is by definition symmetric and positive semi-definite. However, this can be easily achieved as well:
First, initialize a randomized matrix $\vec{R}$ of the desired dimension (that is $D \times D$).
A symmetric and positive semi-definite matrix is then be obtained as follows [45]:

\begin{equation}
\vecsym{\Sigma} = \vec{R} \vec{R}^T.
\end{equation}

Again, this is repeated $K$ times to obtain one covariance matrix per state.

In contrast, the initialization using $k$-means is done as follows: First $K$ clusters are found in the training
dataset using the $k$-means algorithm as implemented by scikit-learn (the implementation uses Lloyd's algorithm [33]).
More concretely, all training observation sequences are “stacked” vertically. The clustering algorithm is then run
on this stacked matrix. Also recall that
$K$ is the number of states, in this case $K=5$. Next, the mean vectors of each state are simply