# Optimal Feedback Control for Modeling Human-Computer Interaction

Florian Fischer  University of BayreuthBayreuthGermany florian.j.fischer@uni-
bayreuth.de ,  Arthur Fleig  University of BayreuthBayreuthGermany
<EMAIL_ADDRESS>,  Markus Klar  University of
BayreuthBayreuthGermany<EMAIL_ADDRESS>and  Jörg Müller
University of BayreuthBayreuthGermany<EMAIL_ADDRESS>
###### Abstract.

Optimal feedback control (OFC) is a theory from the motor control literature
that explains how humans move their body to achieve a certain goal, e.g.,
pointing with the finger. OFC is based on the assumption that humans aim to
control their body optimally, within the constraints imposed by body,
environment, and task. In this paper, we explain how this theory can be
applied to understanding Human-Computer Interaction (HCI) in the case of
pointing. We propose that the human body and computer dynamics can be
interpreted as a single dynamical system. The system state is controlled by
the user via muscle control signals, and estimated from observations. Between-
trial variability arises from signal-dependent control noise and observation
noise. We compare four different models from optimal control theory and
evaluate to what degree these models can replicate movements in the case of
mouse pointing. We introduce a procedure to identify parameters that best
explain observed user behavior. To support HCI researchers in simulating,
analyzing, and optimizing interaction movements, we provide the Python toolbox
OFC4HCI. We conclude that OFC presents a powerful framework for HCI to
understand and simulate motion of the human body and of the interface on a
moment by moment basis.

Optimal Control, OFC, Human-Computer Interaction, Modeling, Parameter Fitting,
Aimed Movements, Mouse Pointing, LQR, LQG, Second-order Lag, Minimum Jerk,
Intermittent Control

††copyright: rightsretained††journal: TOCHI††journalyear: 2022††journalvolume:
1††journalnumber: 1††article: 1††publicationmonth: 1††doi:
10.1145/3524122††journal: TOCHI††ccs: Human-centered computing HCI theory,
concepts and models

##  1\. Introduction

We address the problem of understanding, and modeling, how users control a
virtual end-effector when interacting with computers. Traditionally, the field
of Human-Computer Interaction (HCI) has concentrated on models such as Fitts’
Law (Fitts, 1954; Fitts and Peterson, 1964), predicting summary statistics of
the movement such as movement time. Recently, more attention has been paid to
modeling the underlying process by which the end-effector is controlled,
predicting not only movement time, but end-effector position, velocity, and
acceleration sequences, as well as applied forces (e.g., (Müller et al., 2017;
Fischer et al., 2021; Álvarez Martín et al., 2021; Do et al., 2021)).

We argue that, in order to understand how users control user representations
(e.g., mouse pointer) (Seinfeld et al., 2020), or virtual objects, the field
of Human-Computer Interaction needs to learn more from human motor control.
While human motor control mainly addresses the question of how humans control
the movement of their body, the theories developed there also apply to and can
be adapted to the question of how humans control the state of a computer,
e.g., movement of the mouse pointer.

In the field of human motor control, modern understanding of human movement is
based on the theory of optimal feedback control (OFC) (Todorov and Jordan,
2002; Diedrichsen et al., 2010). This theory understands the human body, and
possibly the environment the body is interacting with, as a dynamical system
that can be controlled, e.g., via muscle control signals. Body and environment
put constraints on this control, e.g., via the system dynamics and constant
and signal-dependent motor noise. The theory assumes that humans continuously
observe the state of their own body and the environment they are interacting
with, e.g., by processing visual and proprioceptive signals. Humans are
assumed to control their body optimally with respect to an internalized cost
function, while respecting the constraints given by the system dynamics and
motor noise.

We believe that the OFC framework enables a better connection between the
field of HCI and recent advances in neighboring scientific disciplines, such
as the study of human movement in motor control (Schmidt et al., 2018; Flash
and Hogan, 1985) and neuroscience (Shadmehr and Wise, 2005). However, OFC
models are not very well known in the field of HCI, yet. In particular, it has
not yet been shown whether these models, developed to model how humans control
their body, can be used to model how users behave during interaction.

The objective of this work is to examine the applicability of optimal feedback
control to HCI, using the example of mouse pointing. The contribution of this
paper is fourfold:

First, we propose a unifying optimal control framework for understanding
movement in interaction with computers. This framework allows to predict the
kinematics and dynamics of the entire movement trajectory, including, e.g.,
end-effector position, or muscle excitation.

Second, we present the first qualitative and quantitative evaluation to what
degree different optimal control models (either open- or closed-loop,
deterministic or stochastic) can replicate movements of the mouse pointer. To
the best of our knowledge, these models have not yet been evaluated
quantitatively regarding their ability to predict movement trajectories during
interaction. We also discuss the possibilities and limitations of the
presented models regarding their suitability for other HCI tasks such as
target tracking, path-following, or handwriting.

Third, we propose a generic parameter fitting process, which can be used to
identify the components of both the system dynamics and the cost function that
best explain observed user behavior, using any desired optimal control model.
For each of the presented models, we systematically analyze the individual
effects of the parameters and show how the proposed parameter fitting can be
used to explain typical differences between users and/or task conditions,
which would remain hidden when using summary statistics only.

Fourth, we provide OFC4HCI, an open-source toolbox accessible from our GitHub
repository111https://github.com/fl0fischer/OFC4HCI that contains the
underlying Python code of this paper. This toolbox includes easy-to-use
scripts for three main use cases: running simulations of human pointing
movements using any of the presented control methods, comparing the resulting
trajectories to data from the Pointing Dynamics Dataset, and optimizing the
parameters of a given control model. While the focus of this toolkit currently
is on (one- or multidimensional) pointing tasks, using the toolkit, extensions
to other HCI tasks such as target tracking, keyboard typing, or gesture-based
input methods are possible.

Our results suggest that stochastic OFC models are able to explain average
user behavior significantly better than models that only account for
simplified movement dynamics (second-order lag) or pure kinematic models (jerk
minimization). In addition, stochastic models such as LQG are able to fit the
distribution of entire trajectories, given a specific user and task condition.
Moreover, the fitting is significantly better than using the recently proposed
Intermittent Control (IC) model (Álvarez Martín et al., 2021) with respect to
both KL divergence (Kullback and Leibler, 1951) and the $2-$Wasserstein
distance (Olkin and Pukelsheim, 1982) serving as evaluation metrics. The
considered deterministic OFC model, which does not take into account any noise
terms, is able to predict average user behavior, given a slightly modified
cost function.

We strongly believe that a proper modeling of the underlying control process
can provide intuition to interface designers as to why users move the way they
do during interaction, and enables a deeper understanding of the impact of
parameters of the interface and input device on the process of interaction. In
the long term, such models could be used for automated optimization of the
parameters of interaction techniques and input devices. Models that work in
real-time could further be used in predictive interfaces, which anticipate
what the user wants to do and respond accordingly, such as pointing target
prediction (Asano et al., 2005). While we will focus on the example of mouse
pointing throughout this paper, the framework we present is generic and
suitable for a wide range of pointing devices, using, e.g., joysticks,
keyboards, pens, touch-based input, mid-air gestures, etc.
The paper is structured as follows:

In Section 2, we start with a short overview of existing models and methods
from the fields of Human-Computer Interaction, Human Motor Control, and
Optimal Control Theory. The proposed optimal control framework for Human-
Computer Interaction is then introduced in Section 3. The models presented in
this paper are evaluated against an existing dataset of one-dimensional
pointing movements, which is described in Section 4.1. The generic parameter
fitting process we use to identify the model parameters that best explain
observed user behavior is described in Section 4.2.

In Sections 5-10, different optimal control models are presented, analyzed,
and adapted to the case of mouse pointing. Moreover, the predicted movements
are compared against user data. Since this paper is also supposed to serve as
a tutorial to Optimal Feedback Control for HCI researchers and interaction
designers, we start with an analysis of the individual components of the OFC
framework before combining them into a final model. In Section 5, we start
with a basic model of movement dynamics, the second-order lag, which has been
used to describe the overall human-computer system dynamics (Müller et al.,
2017) and serves as a baseline for the presented optimal control models. The
idea of (open-loop) optimal control is introduced in Section 6, using the
minimum jerk model (Flash and Hogan, 1985). In Section 7, both movement
dynamics and the assumption of optimality are integrated into one closed-loop
OFC model, the Linear-Quadratic Regulator, which is based on the assumptions
of linear dynamics and quadratic costs (Todorov, 2005). From a didactic point
of view, it is important to develop a thorough understanding of deterministic
optimal feedback control before progressing to stochastic OFC (SOFC) models.
For this reason, we first start with the (substantially simpler) deterministic
case, which can be used to predict average human movement. In Sections 8 and
9, we extend this framework to the general stochastic case by adding different
sensory-input models along with Gaussian motor and sensory noise (thus denoted
as Linear-Quadratic Gaussian Regulator). We compare the stochastic OFC models
to a recently proposed Intermittent Control model (Álvarez Martín et al.,
2021), which is briefly described in Section 10.

Finally, both qualitative and quantitative comparisons between all considered
models are given in Section 11. Difficulties and limitations of the proposed
framework with regard to its applicability to other HCI tasks are discussed in
Section 12, together with some practical advice for HCI researchers, and
conclusions are drawn in Section 13.

##  2\. Related Work

In the field of Human-Computer Interaction, interaction is most commonly
understood as a sequence of discrete actions, which is reflected in the
classification of tasks, such as command selection or target acquisition
(Buxton and Baeker, 1987). In particular, movement, e.g., of the mouse
pointer, is often reduced to summary statistics. The most prominent example is
the dependency of movement time $\operatorname{MT}$ from distance $D$ and
width $W$ of targets, which is described by Fitts’ Law (Fitts, 1954; Fitts and
Peterson, 1964) as $\operatorname{MT}=a+b\operatorname{ID}$, with Index of
Difficulty (ID) usually defined as $\operatorname{ID}=\log_{2}(D/W+1)$
(MacKenzie, 1992). This affine relationship has shown to apply for a variety
of tasks, including reciprocal tapping (Fitts, 1954), mouse pointing and
dragging (Card et al., 1978; Gillan et al., 1990), eye-gazing (Ware and
Mikaelian, 1986; Jagacinski and Monk, 1985), reaching with a joystick
(Jagacinski and Monk, 1985; Card et al., 1978), and ellipse drawing (Mottet et
al., 1994). A very good explanation of the information theoretic
interpretation of Fitts’ Law has been provided by Gori et al. (Gori et al.,
2018).

While aggregated metrics of movement trajectories, e.g., movement variability
or movement offset (MacKenzie et al., 2001), have been used to evaluate task
accuracy since the early days of HCI (Jacob et al., 1994), predictive models
of movement kinematics and dynamics are less common. Exceptions include the
works of Williamson (Williamson, 2006; Williamson et al., 2009), which
introduce an information-theoretic model of interaction with a focus on the
amount of uncertainty that is apparent in different sensor and control
channels, and Müller et al. (Müller et al., 2017), in which three feedback
control models (without optimization) are compared regarding their ability to
model mouse pointer movements. However, the former model is originally
designed for the specific needs of brain-computer interfaces, particularly
inference of the user’s intention based on noisy signal channels, whereas the
latter models only describe the biomechanical apparatus, while high-level
factors affecting the movement trajectory such as concrete task requirements
or intrinsic motivations are neglected. Ziebart et al. (Ziebart et al., 2012)
explore the use of inverse optimal control models for pointing target
prediction. They do not make particular a priori assumptions about the
structure of the cost function. Instead, they use an inverse optimal control
approach to fit a generic function with a large number of parameters (36) to a
dataset of mouse pointer movements. While Ziebart et al. (Ziebart et al.,
2012) focus on the application of inverse optimal control to pointing target
prediction, in this paper we investigate the ability of optimal (feedback)
control models to model movement of the mouse pointer more quantitatively.
From an engineering perspective, several interaction techniques that take into
account the underlying end-effector kinematics have been proposed, including
cursor jumping (Asano et al., 2005; Murata, 1998), target expansion (McGuffin
and Balakrishnan, 2005), and increased cursor activation areas (Mott and
Wobbrock, 2014; Chapuis et al., 2009). These approaches are either based on
target likelihood estimates (Ziebart et al., 2012; Murata, 1998) or
extrapolate sensor data measured during runtime (Asano et al., 2005; McGuffin
and Balakrishnan, 2005; Mott and Wobbrock, 2014). Other methods compare
observed trajectories to a set of pre-defined templates in order to predict
the desired end-point (“kinematic template matching”) (Pasqual and Wobbrock,
2014). In general, these methods are restricted to the kinematic end-effector
level, i.e., they ignore the movement dynamics of the human body (which play a
crucial role for non-standard interaction techniques such as gesture-based
input) and cannot be used to model interaction with dynamic objects.

In addition to their functional use in HCI, movement dynamics have been a
research focus within the field of Motor Control for a long time. Various
models have been developed, all of which predict complete trajectories, e.g.,
end-effector position, velocity, and acceleration profiles over the entire
movement (e.g., (Bootsma et al., 2004; Bullock and Grossberg, 1988; Flash and
Hogan, 1985; Flash et al., 2013; Guiard, 1993; Jagacinski and Flach, 2003;
Lacquaniti et al., 1983; Mottet and Bootsma, 1999; Plamondon, 1998)).
Biomechanical and neural models, in addition, explain how these trajectories
are dynamically generated. This can be either done on the joint-, muscle-, or
neuronal level, incorporating quantities internal to the human body such as
joint angles, joint moments, muscle forces and activations, or neural
excitation signals (Rosenbaum et al., 1995; Uno et al., 1989; Nakano et al.,
1999; Kawato, 1993, 1996; Berret et al., 2011).

Many models of motor control are also capable of modeling the characteristic
between-trial variability that is typically observed in human movements. This
variability is mainly attributed to multiple sources of noise within the human
biomechanical and neural system, most of which can be modeled as additive or
multiplicative Gaussian random variables (Sutton and Sykes, 1967; Schmidt et
al., 1979; Slifkin and Newell, 1999; Jones et al., 2002; Todorov, 1998, 2005).
Signal-dependent noise terms, e.g., Gaussians with zero mean and with a
standard deviation that linearly depends on the magnitude of the muscle
control signal, are also considered responsible for the well-known speed-
accuracy trade-off in goal-directed human movements (Sutton and Sykes, 1967;
Schmidt et al., 1979; Harris and Wolpert, 1998). These noise terms have the
effect that larger control signals, which might increase the speed of the end-
effector, also result in larger deviations from the desired end-effector
position. Users thus face a trade-off between accurate achievement of the
desired goal and fast, but noisy movements.

Another well-established finding from human motor control refers to the amount
of information that is used when selecting a specific control signal. Several
experiments suggest that information that becomes available to the controller
during the movement, e.g., proprioceptive and/or visual signals regarding the
end-effector, are utilized to adjust control signals online and to account for
unexpected perturbations (Woodworth, 1899; Meyer et al., 1988; Todorov and
Jordan, 2002; Wang et al., 2001; Thoroughman and Shadmehr, 2000). This is
reflected by feedback control models of movement, which are able to explain
how users correct errors and handle disturbances during the movement. An early
closed-loop model has been provided by Crossman and Goodeve (Crossman and
Goodeve, 1983). They assume that users observe hand and target and adjust
their velocity as a linear function of the distance, as a first-order lag. A
physically more plausible extension of the first-order lag is the second-order
lag (Crossman and Goodeve, 1983; Langolf et al., 1976). These dynamics can be
interpreted as a spring-mass-damper system, where a constant force is applied
to the mass, such that the system moves to and remains at the target
equilibrium. Because of its simplicity and widespread use, we use this model
as a baseline, called  _2OL-Eq_. Other models of human movement include VITE
(Bullock and Grossberg, 1988) and the models of Plamondon (Plamondon and
Alimi, 1997).

The desired trajectory hypothesis (Kawato, 1999) assumes that whenever
disturbances occur (e.g., due to internal control noise or external
perturbations), feedback is used to push the end-effector towards a
predetermined, deterministic trajectory that results from a separated planning
phase. In contrast, Todorov and Jordan (Todorov and Jordan, 2002) have
demonstrated that deviations are corrected only if they interfere with the
task performance, i.e., deviations that are irrelevant for achieving the
desired goal remain ignored. This minimum intervention principle particularly
implies that all task-specific requirements (end-point position, movement
time, accuracy, etc.) need to be reflected by an internal formulation that the
controller has access to.

Optimal control models provide exactly this internal representation by
assuming that humans try to behave optimally with respect to a certain
internalized cost function. Flash and Hogan (Flash and Hogan, 1985) proposed
that humans aim to generate smooth movements by minimizing the jerk of the
end-effector. We call this model _MinJerk_ in the following. Although the
hypothesis that people aim to minimize jerk has been questioned, see, e.g.,
Harris and Wolpert (Harris and Wolpert, 1998), the minimum jerk model is one
of the most established models. For example, it has been successfully used by
Quinn and Zhai (Quinn and Zhai, 2018) to model the shape of gestures on a
word-gesture keyboard.

Most modern theories of motor control are based on optimal feedback control
(OFC), i.e., they combine the assumptions of optimality and continuously
perceived feedback for closed-loop control. Excellent overviews of recent
progress in OFC theory are provided by Crevecoeur et al. (Crevecoeur et al.,
2014) and Diedrichsen (Diedrichsen et al., 2010). An early approach that
models perturbed reach and grasp movements by using the minimum-jerk
trajectory on a moment by moment basis was presented by Hoff and Arbib (Hoff
and Arbib, 1993). A more general, more recent, and better known OFC model is
the Linear-Quadratic Gaussian Regulator (LQG) (Hoff, 1992; Loeb et al., 1990),
which was mainly used by Todorov to model human movement from a sensorimotor
perspective (Todorov, 1998; Todorov and Jordan, 2002; Todorov, 2005). In this
work, we will present and discuss the assumptions and limitations of this
model, and analyze its applicability to standard HCI tasks such as mouse
pointing.

An important limitation of the LQG model (and many other optimal control
models, e.g., (Flash and Hogan, 1985; Uno et al., 1989; Harris and Wolpert,
1998)) is that the exact movement time needs to be known in advance. One way
to circumvent this issue is to use infinite-horizon OFC (Jiang et al., 2011;
Qian et al., 2013; Li et al., 2018), i.e., to formulate the optimal control
problem on an infinite-time horizon. With such models, (quadratic) distance
and effort costs are usually applied continuously, resulting in an optimal
trajectory that consists of both a transient phase (where the end-effector is
moved towards the target) and a steady-state equilibrium (where the end-
effector is kept at the target). The movement time thus emerges implicitly
from the optimal control problem.

Another strand of literature that specifically deals with the duration of
movement has produced the _Cost of Time_ theory (Hoff, 1994; Shadmehr et al.,
2010; Berret and Jean, 2016). To account for the fact that humans value
earlier achievement more than later achievement, this theory assumes that time
is explicitly penalized with a certain cost function (usually hyperbolic or
sigmoidal).

Recently, methods from the field of Reinforcement Learning (RL) have gained
increased attention. These methods are also based on the principles of optimal
control theory, however, they do not require the system dynamics to be known
in terms of equations and formulas, but solely rely on sampling from an
environment that is usually assumed a black-box to the controller. For this
reason, they are generally applicable to arbitrarily complex systems including
highly non-linear dynamics and discontinuous cost functions (Sutton and Barto,
2018).

Cheema et al. (Cheema et al., 2020) have applied recent RL methods to predict
fatigue during mid-air movements, using a torque-actuated linked-segment model
of the upper limb. Building on this work, it has recently been shown that RL
applied to a more realistic upper-limb model allows to synthesize human arm
movements that follow both Fitts’ Law and the $\nicefrac{{2}}{{3}}$ Power Law
and can predict human behavior in mid-air pointing and path following tasks
(Fischer et al., 2021). Moreover, an extension to mid-air keyboard typing has
been proposed (Hetzel et al., 2021).

In theory, policy-gradient RL methods can also be applied to model interaction
on a muscular level, using state-of-the-art biomechanical models of the human
body (Tieck et al., 2018; Kidziński et al., 2018; Lee et al., 2019; Nakada et
al., 2018). However, the high complexity of the neuromuscular system has so
far imposed considerable restrictions to each of these approaches, including
the reduction of degrees of freedom (Tieck et al., 2018; Nakada et al., 2018;
Kidziński et al., 2018) and the omission of muscle activation dynamics (Nakada
et al., 2018; Lee et al., 2019). Most importantly, for most RL algorithms no
theoretical convergence guarantees exist, which complicates a profound
interpretation or replication of the resulting simulation results (Sutton and
Barto, 2018). For this reason, in this paper we focus on the well-understood
theory of optimal control, as this allows us to use convergence guarantees
more often, which makes us less reliant on intuition and experience. For
example, the LQR introduced in Section 7 is guaranteed to converge to the
optimal movement trajectory, given a fixed set of parameters. This is a
decisive advantage compared to pure RL-based methods, as it allows to compare
optimal trajectories for different task conditions, cost functions, and user
models.

In summary, the fundamental question of human movement coordination has
produced a vast literature and a deep understanding of the nature of human
movement. Given that almost every interaction of humans with computers
involves movement of the body, it is surprising that this field is little
known, and applied, in HCI. It is important to bear in mind, however, that
most of these theories intend to model movement of the human body per se. In
HCI, it is also relevant how users control the movement of user
representations (e.g., mouse pointers) (Seinfeld et al., 2020) and virtual
objects in the computer. Since the control of user representations and virtual
objects is mediated by input devices, operating systems, and programs,
requires high precision, and is often learnt very well, it is unclear how the
theory of human motor control can be applied to the HCI context. To our
knowledge, it has not yet been investigated whether the above optimal motor
control models can be applied to HCI tasks such as mouse movements. Adapting
and validating such models regarding their ability to model HCI tasks such as
pointing thus remains an open research question for HCI.

In order to leverage the strengths of recent motor control theory in the field
of HCI, we believe that a general optimal control framework for Human-Computer
Interaction is necessary, which can explain both how and why humans behave in
interaction with arbitrary interfaces on a continuous level. Such a framework
constitutes a natural extension of the principle of “designing interaction,
not interfaces” (Beaudouin-Lafon, 2004) by conceptualizing interaction based
on neuroscientific, psychological, and biomechanical insights within one
coherent and mathematically profound framework.

In the following section, we will introduce the optimal control framework for
Human-Computer Interaction and explain its main constituents using the example
of mouse pointing.

##  3\. Introducing the Framework

Forward
---
Model
Computation:
---
$\min_{u}J_{N}(x,u)$
Human
---
Controller
Task$J_{N}$ Human
---
Body Dynamics
Interface
---
Dynamics
Human-Computer System Dynamics Muscle
---
Control
Signals $u^{*}$
Optimal State Trajectory $x^{*}$
---
HumanComputer

(a) Open-Loop Optimal Control

Forward
---
Model
Computation:
---
$\min_{u}J_{N}(x,u)$
Control Strategy
---
$u_{n}=\pi(\hat{x}_{n})$
Human
---
Observer
Forward
---
Model
$-$ $+$ Human
---
Controller
Task$J_{N}$ Human
---
Body Dynamics
Interface
---
Dynamics
Human
---
Perception
Display
---
Human-Computer System DynamicsFeedback Muscle Control
---
Signal $u_{n}$
Efference
---
Copy
State $x_{n+1}$ Optimal State Trajectory $x^{*}$
---
Observation
---
$y_{n+1}$
State
---
Estimate $\hat{x}_{n+1}$
Before MovementDuring MovementAfter Movement

(b) Closed-Loop Optimal Control (Discrete-Time Formulation)

Figure 1. In our optimal control framework for Human-Computer Interaction, the
user is assumed to _control_ the state $x$ of the interactive system, which
incorporates both the body state (e.g., arm and finger position) and the
interface state (e.g., mouse pointer position and velocity), and which evolves
according to the Body and Interface Dynamics. We assume that the user computes
the controls $u$ through _optimization_ , i.e., by minimizing a cost function
$J_{N}$ (e.g., incorporating time or effort costs) that depends on the task.
(a) In an open-loop model, this calculation is only based on an internal
Forward Model of the Human-Computer System Dynamics. The optimal state
trajectory $x^{*}$ is obtained by applying the resulting muscle control
signals $u^{*}$ in one forward pass. The Forward Model does not have to
coincide with the System Dynamics.
(b) A closed-loop model takes into account effects that appear only after
execution. The key difference in the Computation block is that, instead of
optimal control signals $u^{*}$, it yields an optimal Control Strategy $\pi$
that is computed before movement onset. At each time step $n$, this Control
Strategy is used to map an arbitrary (estimated) state $\hat{x}_{n}$ to the
corresponding optimal control $u_{n}$. Based on the resulting state $x_{n+1}$,
an observation $y_{n+1}$ is obtained via Feedback, which incorporates
descriptions of both the Display and the Human Perception. The Human Observer
then compares this observed state to an expected state it computes using an
efference copy of the current control signal $u_{n}$ and the Forward Model.
Based on the resulting difference between expected and observed signals, an
internal state estimate $\hat{x}_{n+1}$ is computed and used to select the
next control $u_{n+1}$, and so on.

Fully described in caption and text.

Modern motor control theory assumes that humans aim at controlling their
movements optimally, given the constraints imposed by the body and
environment. Important constraints are imposed by physics, e.g., via Newton’s
second law, i.e., force equals mass times acceleration, and by the muscles,
which cannot create forces instantaneously, but need to build up muscle
activation (and thus force) over time. In the case of Human-Computer
Interaction, constraints are imposed not only by the human body, but also by
the input device, sensor, and processing within the computer. Further, the
human perceptual system does not have direct access to the state of the world,
but can only observe certain variables that depend on the state and needs to
build up an internal estimate of the true world state over time.

Since these properties are characteristic for almost any HCI task, we propose
a generic optimal control framework for Human-Computer Interaction, which
consists of four submodels that continuously interact with each other:

  * •

the Human-Computer System Dynamics, which describe the biomechanics of the
considered body parts as well as the dynamics of the resulting interaction
with the application interface via an input device,

  * •

the Human Controller, i.e., the decisive part of the brain, which selects the
new muscle control signals,

  * •

the Feedback, which models how environmental information is sensed by the
human,

  * •

and the Human Observer, a cognitive model for how perceived sensory signals
from the Feedback are processed and evaluated.

Since the computer operates in discrete time, we use discrete-time dynamics,
i.e., we consider time steps $n\in\\{0,\dots,N\\}$ up to a final step
$N\in\mathbb{N}$, with each time step corresponding to $h$ seconds. However,
the proposed framework is designed to be as general as possible and the
following explanation also applies to the continuous-time case.

Figure 1 illustrates the relationship between the four components of the
Human-Computer Interaction loop, specifically distinguishing between open-loop
and closed-loop optimal control models. In the following, we will give a brief
description of the proposed framework, with focus on the differences between
both variants. Subsequently, the four submodels are explained in detail in
separate subsections, introducing a more technical and mathematically rigorous
notation. Readers who are already familiar with the differences between open-
and closed-loop models can proceed directly to Section 3.2.

###  3.1. Open-loop vs. closed-loop models

Open-loop models, as depicted in Figure 1(a), cannot infer any information
from the System Dynamics after applying muscle control signals $u$. For this
reason, the Human Controller block in Figure 1(a) does not depend on the
output of the System Dynamics block, but depends only on an internal Forward
Model. In particular, the Feedback and Human Observer blocks are not part of
the generic open-loop framework at all. For this reason, open-loop models
allow for a strict separation between planning and execution phase of a
movement, i.e., trajectories $x^{*}$ that are optimal with respect to the
objective function $J_{N}$ can be obtained in a two-step process. First, an
(open-loop) optimal control problem (OCP) is solved (Computation block), i.e.,
a sequence of controls $u^{*}=(u_{n}^{*})_{n\in\\{0,\dots,N-1\\}}$ is found
such that $J_{N}(x,u)$ becomes minimal among all permissible control sequences
$u$, given an initial state (see Section 3.3 for more details). Second, the
resulting optimal control sequence is applied in a forward pass, i.e., the
system dynamics are evaluated at each time step $n\in\\{0,\dots,N-1\\}$ to
obtain the subsequent optimal state $x_{n+1}^{*}$.

Note that if the system dynamics are deterministic and the forward model
internally used to compute $u^{*}$ matches these dynamics, there is no
advantage in computing the optimal controls online, that is, only at the time
they are executed during the forward pass. However, there are a few scenarios
where the internal forward model used for optimization cannot fully predict
the actual behavior of the human body and interface dynamics, i.e., the
outcome of the Human-Computer System Dynamics block in Figure 1. This is
particularly the case

  * •

if the OCP is stochastic, e.g., in the case of motor noise,

  * •

if there is a mismatch between internal model and actual system dynamics, or

  * •

if unexpected disturbances occur that the internal model did not account for.

If one of these assumptions holds (which usually is the case in practice), the
controller will benefit from any information it receives during execution, as
this allows to condition the choice of an individual control $u_{n}$ on the
true state in a closed-loop manner, instead of using a prior state estimate.
Note that the feedback loop immediately eliminates the strict separation
between planning and execution phase, which is prevalent for open-loop models.
Instead, the optimal control sequence $u^{*}$ needs to be computed online in
an iterative manner. In particular, all information available to the
controller at a time step $n\in\\{0,\dots,N-1\\}$ is used to compute the
muscle control signal $u_{n}^{*}$ at this time step, which is applied to the
Human Body Dynamics. In combination with the Interface Dynamics, this results
in a new system state $x_{n+1}$. The state (or partial information thereof,
see Section 3.4) is then fed back to the controller, which again selects the
next control $u_{n+1}^{*}$ based on this information, i.e., a sensorimotor
control loop between Human and Computer is established (see Figure 1(b)). In
particular, the optimal state trajectory $x^{*}$ does not only depend on the
control sequence $u$, but also vice versa. Since in such models, feedback is
given to the controller during execution, these models are often denoted as
optimal feedback control (OFC) models (Diedrichsen et al., 2010; Crevecoeur et
al., 2014).

It is important to note that with many OCP solution methods that take into
account feedback during execution (including the ones presented in this
paper), the actual optimization can be performed offline, i.e., before the
controls are applied to the actual system dynamics. Instead of a single
optimal control sequence $u^{*}$, such methods usually yield an optimal
Control Strategy, that is, a function $\pi:\mathcal{X}\rightarrow\mathcal{U}$
that maps an arbitrary state $x\in\mathcal{X}$ to a corresponding control
$u\in\mathcal{U}$ that is optimal starting from this state.

In the later sections, we will present both an open-loop model (Section 6) as
well as different variants of one of the most widely used OFC models (Sections
7-9).

###  3.2. Human-Computer System Dynamics

The Human-Computer System Dynamics form the basis of each optimal control
model of Human-Computer Interaction and consist of two parts: the Human Body
Dynamics and the Interface Dynamics.

_Human Body Dynamics._ Given a vector of neural muscle control signals $u$,
the Human Body Dynamics describe how these signals are transformed into joint
torques and accelerations. The corresponding joint postures are obtained via
integration. If required, kinematic models that map joint postures to world-
centered positions of arbitrary body parts can be included as well, e.g., to
get more realistic movement of the wrist or the index finger based on the
computed joint angles.

_Interface Dynamics._ In interaction with computers, the forces and
accelerations resulting from body dynamics are applied to the physical input
device, e.g., the mouse device. Inertial properties of this physical end-
effector determine the input device motion, which is sensed, filtered, and
mapped to a motion of the corresponding virtual end-effector. For example, in
the case of mouse pointing, the mouse pointer position might result from the
application of a pointing transfer function to the mouse device velocity
measured via optical sensors (Casiez et al., 2008; Casiez and Roussel, 2011).
The virtual end-effector is then used for interaction with graphical user
interfaces. These dynamics of the computer part, including physical properties
of the input device as well as visualizations and animations that appear from
interaction with buttons, sliders, etc., are summarized within the Interface
Dynamics block in Figure 1.

Combined, the Human-Computer System Dynamics yield the new state vector $x$,
which incorporates all relevant information about the current state of the
human body, of the physical and/or the virtual end-effector, and of the
applications. In the discrete-time formulation used in this paper, both the
state and control vectors $x_{n}$ and $u_{n}$ are given at time steps
$n\in\\{0,\dots,N\\}$ up to a final step $N\in\mathbb{N}$, with each time step
corresponding to $h$ seconds. The next state $x_{n+1}$ depends on the current
state $x_{n}$ and control $u_{n}$, which in most cases can be formalized as

(1) $x_{n+1}=f(x_{n},u_{n}),$

where $x_{0}$ is a given initial state.222This is closely related to a
continuous-time formulation based on differential equations, where the control
is assumed to be piecewise constant (i.e., it only changes once every $h$
seconds).

This system dynamics function $f$ can be either deterministic or stochastic.
In the deterministic case, starting from the current state
$x_{n}\in\mathbb{R}^{k}$ ($k\in\mathbb{N}$) and applying a control
$u_{n}\in\mathbb{R}^{m}$ ($m\in\mathbb{N}$), the subsequent state
$x_{n+1}\in\mathbb{R}^{k}$ is uniquely determined by the function
$f:\mathbb{R}^{k}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{k}$. In the
stochastic case, $x_{n+1}$ randomly emerges from a set of possible successor
states, according to a conditional probability distribution $p$, i.e.,
$x_{n+1}=f(x_{n},u_{n})\sim p(\cdot\mid x_{n},u_{n})$. Examples of stochastic
systems include, e.g., Body Dynamics with internal motor noise or Interface
Dynamics with noisy input signals.

It is important to note that the controller, which is described in Section
3.3, is agnostic to the partitioning of the system dynamics into effects
attributed to the Body Dynamics and effects attributed to the Interface
Dynamics. All the controller needs to know is the overall system dynamics $f$
(or an internal approximation of it), which maps an arbitrary state-control
pair $(x_{n},u_{n})$ to the subsequent state $x_{n+1}$ reached after $h$
seconds. Thus, an optimal control model of Human-Computer Interaction can be
instantiated in two ways. The first one is to include accurate submodels of
arbitrary granularity (e.g., a separate model for each muscle activation, arm
and hand dynamics, input device dynamics, and application dynamics), and
combine them along the interaction loop into one aggregated system dynamics
function $f$. However, the framework also allows to test whether some generic
dynamics, such as a spring-mass-damper system or simplified muscle activation
dynamics, are suitable to model the overall Human-Computer System Dynamics for
a given task setting. The focus of this paper will be on the latter approach,
since we aim to start with an easily understandable and well-established model
from the field of Human Motor Control, and test whether these dynamics are
applicable to the context of mouse pointing. We believe that this approach is
well suited for introducing optimal control methods to HCI without going too
much into (biomechanical) detail.

The system dynamics of all models considered in this paper are linear in both
the state and the control, i.e.,

(2) $x_{n+1}=f(x_{n},u_{n})=Ax_{n}+Bu_{n}.$

Here, the matrix $A$ describes how the dynamics evolve when no control is
exerted. The matrix $B$ describes how the control influences the system.

At first glance, this assumption seems to be very limiting, especially with
regards to the complexity of the human neuro- and biomechanical system, as
well as of most interaction methods and application GUIs, However, the tools
and methods proposed in this paper for the case of linear dynamics are also
beneficial for more complex models of interaction, which, for example, include
muscle-driven models of the human body or non-linear pointer acceleration
functions. Using a reference trajectory, it is always possible to linearize a
non-linear system around this particular trajectory in order to obtain a
linear system of the form (2), which locally approximates the non-linear one.
While linearization-based extensions of the considered optimal control models
to the non-linear case have been proposed (Li and Todorov, 2004; Todorov and
Li, 2005; Tassa et al., 2014), an application of these methods to typical HCI
tasks would be an interesting area for future work.

The main advantage of using linear dynamics is that, when combined with
quadratic costs and Gaussian noise, the resulting optimal control problem (see
Section 3.3) can be solved analytically and thus quickly and exactly. Finally,
the linear case is easier to understand and formalize and thus well suited for
the explanatory purposes of this paper.

In the case of mouse pointing, which usually requires only small movements of
the arm, the hand, and the input device, linearization around a single
trajectory, i.e., using constant system matrices $A$ and $B$ as in (2), is a
reasonable initial approach to model (moderate) mouse movements. Indeed, we
will show that linear system dynamics can account for many phenomena that are
characteristic in mouse pointing.

###  3.3. Human Controller

In general, various control sequences can produce the same movement
trajectory. For example, the arm can rest on the table or stay in the air, as
long as the mouse device is controlled appropriately. This is referred to as
the joint-redundancy problem (Bullock and Grossberg, 1988). For a large number
of degrees of freedom, e.g., motor signals that are applied to individual
muscles, the same goal can be achieved with different controls, raising the
question of which control is actually chosen by the central nervous system
(CNS) and why.333Moreover, in the case of muscle-driven simulations, the set
of feasible controls is relatively small compared to the total decision space,
which makes it even less clear how appropriate controls are internally found
(Valero-Cuevas, 2015). This fundamental question, however, cannot be answered
using movement dynamics alone. Instead, the optimal control framework has been
proposed to address this question (Friston, 2011; Kording, 2007; Todorov and
Jordan, 2002).

#### Optimal Control Problems (OCPs)

Optimal control methods make use of a specific cost function, which is to be
minimized. Previous approaches include minimization of either jerk (Flash and
Hogan, 1985; Hoff and Arbib, 1993), peak acceleration (Nelson, 1983), end-
point variance (Harris and Wolpert, 1998), duration (Artstein, 1980; Tanaka et
al., 2006), or torque-change (Uno et al., 1989), among others. Different
objectives can be combined in one cost function to model trade-offs, e.g.,
between accuracy and effort (Todorov, 1998; Li and Todorov, 2004), accuracy
and stability (Liu and Todorov, 2007), or jerk and movement time (Hoff, 1994).
Recently, it has been argued that several abilities associated with
intelligence such as knowledge, perception, or imitation naturally emerge from
behaving optimally with respect to an ultimate goal (Silver et al., 2021). For
goals that can be expressed by an adequate cost function, this particularly
suggests that the optimal control framework is able to explain intelligent
behavior.

In general, the (finite-horizon) discrete-time optimal control problem is
given by

(3c) $\displaystyle\begin{gathered}\textsl{Minimize}\quad
J_{N}(x,u)=g_{N}(x_{N})+\sum_{n=0}^{N-1}g(x_{n},u_{n})\\\ \textsl{with respect
to
}u=(u_{n})_{n\in\\{0,\dots,N-1\\}}\subset\mathcal{U}\subset\mathbb{R}^{m},\end{gathered}$
where $x=(x_{n})_{n\in\\{0,\dots,N\\}}\subset\mathcal{X}\subset\mathbb{R}^{k}$
satisfies (3f)
$\displaystyle\begin{gathered}x_{n+1}=f_{\text{int}}(x_{n},u_{n}),\enspace
n\in\\{0,\dots,N-1\\},\\\ x_{0}=\bar{x}_{0}\end{gathered}$ for some given
initial state $\bar{x}_{0}\in\mathcal{X}\subset\mathbb{R}^{k}$.

Here, $f_{\text{int}}$ denotes the dynamics of the internal model used for
optimization. In Figure 1, this corresponds to the Forward Model block within
the Human Controller. In most optimal control models, including those
considered in this work, the internal model is assumed to be exact, i.e., it
matches the actual system dynamics (corresponding to the Human-Computer System
Dynamics block in Figure 1), which are analogously described by some function
$f$. The objective function $J_{N}(x,u)$ that we want to minimize consists of
some _terminal cost_ function
$g_{N}:\mathbb{R}^{k}\rightarrow\left[0,\infty\right[$ and the sum of _running
costs_ $g:\mathbb{R}^{k}\times\mathbb{R}^{m}\rightarrow\left[0,\infty\right[$
accumulated over $N$ time steps. These might be chosen dependent on the task
under consideration. For example, in a tracking task, the distance between
end-effector and target could be penalized in each step, whereas in a steering
task, large costs might be applied whenever one of the bounds is reached. The
sets $\mathcal{X}$ and $\mathcal{U}$ can be used to restrict the states and
controls that are permissible at each time step. In the following, however, we
will set $\mathcal{X}=\mathbb{R}^{k}$ and $\mathcal{U}=\mathbb{R}^{m}$, i.e.,
we do not impose any restrictions. The initial state $x_{0}=\bar{x}_{0}$ is
assumed to be given, and the (unique) optimal solution $(x,u)$ to an OCP
(assuming that it exists) is denoted by $(x^{*},u^{*})$ in the following.

For deterministic OCPs, both the internal model and the actual system dynamics
are given by deterministic functions
$f_{\text{int}}:\mathbb{R}^{k}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{k}$
and $f:\mathbb{R}^{k}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{k}$,
respectively. For stochastic OCPs, these dynamics are replaced by conditional
probability distributions (see Section 3.2). It is important to note that
albeit in stochastic OCPs, the concrete successor state $x_{n+1}$ resulting
from the application of a hypothetical control $u_{n}$ in the state $x_{n}$ is
not available to the controller during optimization, the underlying transition
probability density function $p(\cdot\mid x_{n},u_{n})$ usually is.444In the
case of unknown transition dynamics $p$ (or $f$, in the deterministic case),
the controller would need to rely on sampling transitions from the environment
in order to be able to estimate the expected future costs of different
controls. This problem is addressed in Model-Free Reinforcement Learning.
Stochastic OCPs are thus capable of modeling the between-trial variability
that typically occurs in human movement.

The methods to find the optimal solution of an OCP depend on its problem
structure, i.e., the properties of the cost function (e.g., differentiability,
convexity), the system dynamics (e.g., linearity, stochasticity), and the
permissible state space $\mathcal{X}$ and control space $\mathcal{U}$. Often,
the optimal solution can only be determined approximately, using numerical
methods such as Multiple Shooting (Bock and Plitt, 1984), Direct Collocation
(Betts, 2010), or Reinforcement Learning (Sutton and Barto, 2018). However, in
some cases an explicit solution formula exists that yields the exact (and
unique) optimal control sequence $u^{*}$. In this paper, we will focus on the
most widely known class of OCPs that allows for such an analytical solution
scheme: those with linear system dynamics and convex, quadratic costs.

###  3.4. Feedback & Human Observer

In the closed-loop case of the optimal control framework for Human-Computer
Interaction (Figure 1(b)), the Feedback block accounts for the fact that
usually not all information included in the state $x$ are (immediately)
available to the user. First, the visual output, which is shown on the
Display, is created based on the respective state components. This information
is then sensed and processed by the Human Perception, which describes how
visual, proprioceptive, and/or auditive signals are perceived and integrated
into the stream of observations $y=(y_{n})_{n\in\\{0,\dots,N-1\\}}$ the
controller can condition on. The same holds for information on the own body
state, which is directly obtained from the Human Body Dynamics, e.g., via
proprioceptive input signals.

In general, these observations might be delayed, noisy, or incomplete. To
decide for an appropriate control $u_{n+1}$ at time step $n+1$, thus an
estimate $\hat{x}_{n+1}$ of the true current state $x_{n+1}$ is required. This
estimate is computed by the Human Observer, which compares the observed state
$y_{n+1}$ to an expected state it computes using an efference copy of the most
recent control signal $u_{n}$ and the Forward Model.555In this paper, we
assume perfect system knowledge, i.e., the forward model consists of the same
system dynamics and perception functions as used in the actual interaction
loop. Based on the resulting difference between expected and observed signals,
an internal state estimate $\hat{x}_{n+1}$ is computed, which is then used by
the Human Controller to select the next muscle control signal $u_{n+1}$,
resulting in the above discussed closed interaction loop.

###  3.5. Applications of the Proposed Framework

In this paper, we analyze and compare several optimal control models of
interaction. It is important to note that not all of the considered models
exploit the complexity of the generic models depicted in Figure 1. For
example, the closed-loop LQR model includes a trivial perceptual model in that
it assumes that the controller has complete access to the true system state
$x$ immediately. Another example is the open-loop MinJerk model, which does
not include task-specific system dynamics.

Each optimal control model, however, yields a continuous representation of all
relevant quantities of the interactive system. In contrast to summary
statistics such as Fitts’ Law, this allows to simulate and predict complete
movement trajectories on both the kinematic and the dynamic level. It also
allows to analyze the effects of the control $u$ and of different cost terms
incorporated in the objective function $J_{N}$, on the human body and the
interface (e.g., user representations, buttons, or sliders).

Most importantly, the modularity of the proposed framework enables high
flexibility and generalizability. For example, it is possible to analyze the
effects of different input devices and/or GUIs on movement trajectories and
control sequences, using the same description of the human biomechanical and
perceptual system. Additionally, a given interface can be evaluated for
different tasks such as pointing, dragging, steering, etc., by modulating the
internal objective function accordingly. The resulting continuous
representations can then be evaluated with respect to different metrics, e.g.,
remaining distance to target (Diedrichsen et al., 2010), effort (Shadmehr et
al., 2016), fatigue (Cheema et al., 2020), movement time (Tanaka et al.,
2006), etc.

Finally, the framework can be used to reverse-engineer the internal objective
function (inverse optimal control) as well as properties of the human
biomechanical system (system identification), such that the resulting
trajectories best fit some experimentally observed user trajectories. Before
we explain how to identify such model-specific parameters using a data-driven
parameter fitting procedure (see Section 4.2), we give a brief overview of the
experimental data we use to evaluate the presented models in this paper.

##  4\. User Trajectories and Parameter Fitting

###  4.1. The Pointing Dynamics Dataset

For the evaluations in this paper, we use the Pointing Dynamics Dataset. Task,
apparatus, and experiment are described in detail in (Müller et al., 2017).
The dataset contains the mouse trajectories for a reciprocal mouse pointing
task in 1D for ID 2, 4, 6, and 8 (12 participants, 8 task conditions, 7732
trajectories in total). We use the raw, unfiltered position data in our
parameter fitting process to avoid artifacts from the filtering process. In
this section, we explain how we pre-process this data for the purposes of this
paper.

Pointing experiments both in the reciprocal and discrete Fitts’ paradigm
introduce reaction times as experimental artifacts. In real mouse usage, users
first decide on a pointing target themselves, and then start moving the mouse.
In this sense, the pointing process can be considered initiated as soon as the
pointer begins to move. In contrast, in the experimental paradigm used in
(Müller et al., 2017), the next trial started as soon as the user clicked the
mouse in the previous trial. The target given to the user appeared at that
instant. This introduces a potential confound in the starting time of each
trial. The beginning of each trial can be partly attributed to belong to the
end of the previous trial, and partly to a reaction time adjusting to the new
target. During this time, velocity and acceleration of the pointer are close
to zero. This reaction time shows considerable variation both within and
between participants.

Because the focus of this paper is not on modeling reaction times, we remove
them from each mouse movement. To this end, we drop all frames before the
velocity reaches $1\%$ of its maximum/minimum value (depending on the movement
direction) for the first time in each trial.666 For improved temporal
alignment of the individual trajectories and to remove outliers that sometimes
occur at the beginning of a movement, we additionally assume that the
acceleration remains positive/negative (depending on the movement direction)
for at least 40ms after the initial time. If no time step is found that
satisfies both criteria, we discard the entire movement. This was the case
only for a single movement (ID 2, participant 5).

Since the deterministic optimal control models considered in this paper can
only predict average user behavior, we compute mean trajectories for each
user, task condition, and direction from the raw data, resulting in 192 mean
trajectories. This is done as follows.

First, we remove outlier trials, where at any time step the position was more
than three standard deviations from the respective mean. This was the case for
397 trajectories in total, i.e., $5.1\%$ of all trials. We found this to be
necessary as the averaging process is highly sensitive to outliers. In
particular, delayed movement onsets, which, e.g., might have occurred due to a
lack of attention of the participant, would inject a high bias into the
computed statistics.

Second, in order to make trials of different length comparable, we assume that
the pointer would not move after the mouse click. Given a set of trials to be
averaged (with reaction times removed as described above), movements shorter
than the longest one are extended by their last position, zero velocity, and
zero acceleration to achieve the same length. To avoid unnecessarily long
trajectories for conditions where a few of the recorded trials were of
exceptional length, we additionally remove trials with duration longer than
three standard deviations from the mean duration of the respective condition,
before extending the remaining trajectories to maximum length. This was the
case for 87 trajectories in total (i.e., $1.1\%$ of all trials), with a
maximum of two trials removed per user, task condition, and direction.
Finally, we average the resulting trajectories on a frame by frame basis. We
also compute the respective sample covariance matrices at each time step to
capture the between-trial-variability observed from user data.

###  4.2. Parameter Identification

In this section, we present a method to identify the parameters of a given
instantiation of the optimal control framework introduced in Section 3 that
best explain experimentally observed user behavior. More specifically, for a
given interaction model and a given dataset of user trajectories, we aim to
find the model-specific parameter values such that the resulting trajectories
approximate a subset of user trajectories (e.g., all trajectories of a
specific user for some task condition) as closely as possible. In the
following, we will apply this procedure to each of the presented models, using
the Pointing Dynamics Dataset as reference data. Since only stochastic models
can account for the between-trial variability typically observed in human
movements, we need to distinguish between the deterministic and the stochastic
case in the following.

In the deterministic case, we use the squared error between predicted and
observed mouse pointer position, summed over time (sum squared error (SSE)) as
a measure of distance between simulation and mean user trajectories. Given a
model with parameter vector $\Lambda$, where $p_{n}^{\Lambda}$ denotes the
position time series of the resulting simulation trajectory, and given a mean
user position time series $p_{n}^{\text{USER}}$, the goal is to find the
parameter vector $\Lambda^{*}$ such that the loss function of the parameter
fitting process

(4)
$\mathcal{L}(\Lambda)=\text{SSE}(\Lambda)=\sum_{n=0}^{N}\left(p_{n}^{\Lambda}-p_{n}^{\text{USER}}\right)^{2}$

takes its minimum in $\Lambda^{*}$. This is done for each mean trajectory,
resulting in 192 optimal parameter vectors.

Minimizing Equation (4) with respect to $\Lambda$ can be considered a least-
squares problem (Björck, 1996), where each function evaluation of
$\mathcal{L}(\Lambda)$ requires computation of the respective model simulation
trajectory to obtain $p_{n}^{\Lambda}$. In the case of optimal control models,
this particularly implies that an OCP must be solved to obtain
$p_{n}^{\Lambda}$ for a given $\Lambda$, that is, the complete parameter
fitting process consists of two nested optimizations.

Stochastic models, in contrast, yield a sequence of distributions of the
state, such as the distributions of pointer position and velocity, over
multiple trials. These distributions can be used to sample individual
trajectories. We measure the “similarity” between two distributions using the
$2-$Wasserstein distance (often denoted as Earth mover’s distance) (Olkin and
Pukelsheim, 1982). Given two normal distributions $\rho_{1}$ and $\rho_{2}$
with means $\mu_{1}$ and $\mu_{2}$ and covariance matrices $\Sigma_{1}$ and
$\Sigma_{2}$, respectively, the $2-$Wasserstein distance can be written as

(5)
$W_{2}\left(\rho_{1},\rho_{2}\right)=\left(\|\mu_{1}-\mu_{2}\|_{2}^{2}+\text{tr}(\Sigma_{1})+\text{tr}(\Sigma_{2})-2\text{tr}\left((\Sigma_{1}\Sigma_{2})^{\frac{1}{2}}\right)\right)^{\frac{1}{2}},$

and can be interpreted as the amount of work required to transform the
probability distribution $\rho_{1}$ into the probability distribution
$\rho_{2}$ (and vice versa). In the following, we will use this formula to
measure the distance between the simulation state distribution of a model with
parameter vector $\Lambda$, $\rho^{\Lambda}$, and the empirically observed
state distribution, $\rho^{\text{USER}}$, at some time step
$n\in\\{0,\dots,N\\}$, i.e., $\rho_{1}=\rho_{n}^{\Lambda}$ and
$\rho_{2}=\rho_{n}^{\text{USER}}$.

One advantage of this measure is that it is only based on the relative
distance between the two means and covariance matrices of the two
distributions, independent of the magnitude of these quantities. This is in
contrast to the KL divergence (Kullback and Leibler, 1951), which increases as
the variance of the reference distribution decreases. In the special case
where both distributions have the same, diagonal covariance matrix, the
$2-$Wasserstein distance corresponds to the Euclidean distance between the
means of both distributions.

As a measure for the similarity between complete sequences of distributions
(e.g., of mouse pointer positions and velocities), we use the mean Wasserstein
distance (MWD) over time:

(6)
$\mathcal{L}(\Lambda)=\text{MWD}(\Lambda)=\frac{1}{N+1}\sum_{n=0}^{N}W_{2}\left(\rho_{n}^{\Lambda},\rho_{n}^{\text{USER}}\right).$

In both the deterministic and the stochastic case, we solve the (outer)
optimization problem of minimizing $\mathcal{L}(\Lambda)$ with respect to
$\Lambda$ using differential evolution (Storn and Price, 1997), which is a
simple, gradient-free global optimization algorithm suitable for continuous
parameter spaces. This algorithm has proven to yield robust and reliable
results even for ill-conditioned problems (Auger et al., 2009).

Of course, more efficient optimization methods, e.g., gradient-based ones, are
always desirable, and algorithmic differentiation is a promising step forward
in that regard. The main question here is the applicability of algorithmic
differentiation in the case of iteratively alternating between control and
estimation problems, as is required for the considered case of LQG with
signal-dependent noise (see Section 8.1). Pursuing this endeavor, however,
might very well enable real-time predictions of parameter effects on the
entire interaction loop.

DE$\text{MODEL}_{\Lambda}$$J_{N}$$\mathcal{L}$$\Lambda^{0}$$\Lambda^{*}$$\Lambda$
$p^{\Lambda}$ / $\rho^{\Lambda}$
---
$\text{SSE}(\Lambda)$ /
---
$\text{MWD}(\Lambda)$
$x^{\Lambda^{*}}$ / $\chi^{\Lambda^{*}}$
---
Deterministic Parameter Fitting |  Stochastic Parameter Fitting

Figure 2. Starting with an initial parameter vector $\Lambda=\Lambda^{0}$, the
differential evolution (DE) algorithm obtains the loss $\mathcal{L}(\Lambda)$
(sum squared error (SSE) value in the deterministic and mean Wasserstein
distance (MWD) in the stochastic case) for the currently considered parameter
vector $\Lambda$. To do this, it calls $\text{MODEL}_{\Lambda}$, which
computes the resulting model trajectory sequence $x^{\Lambda}$ (or a sequence
of state distributions $\chi^{\Lambda}$). In case of optimal control models,
this requires an inner optimization with respect to a model-specific objective
function $J_{N}$. The resulting position time series $p^{\Lambda}$ included in
$x^{\Lambda}$ (or the sequence of position-velocity distributions
$\rho^{\Lambda}$ included in $\chi^{\Lambda}$) is used to compute the loss
$\mathcal{L}(\Lambda)$ of the current parameter vector $\Lambda$, which in
turn is used by the DE algorithm. After obtaining the loss for a few requested
parameter vectors $\Lambda$, the DE algorithm chooses the next parameter
vectors $\Lambda$ until $\Lambda^{*}$ with minimum loss is found. Finally,
$\Lambda^{*}$ is returned along with the respective optimal trajectory
$x^{\Lambda^{*}}$ (or the optimal sequence of state distributions
$\chi^{\Lambda^{*}}$).

Workflow diagram of the proposed parameter fitting process, both for the
deterministic and the stochastic case. Details are described in the caption
and in the text.

Figure 2 gives an overview of our parameter identification process for both
the deterministic and the stochastic case. The parameter boundaries for all
models introduced below are given in Table B.1 in the Appendix. Descriptions
of how discrete parameters are relaxed in order to optimize them using
standard continuous optimization methods are given in the respective model
sections.

##  5\. Pointing as a Dynamical System: The Second-Order Lag

One of the basic models of mouse pointer dynamics is the second-order lag,
which has been used as a baseline in several papers, including (Müller et al.,
2017; Álvarez Martín et al., 2021). We therefore also include it as a
baseline. The parameters $\Lambda$ of the model are the stiffness of the
spring $k>0$ and the damping factor $d>0$. In the setting described below, the
mass is a redundant parameter, and we thus set it to 1. In continuous time, we
denote the position of the mouse pointer as $y(t)$, and its first and second
derivatives with respect to time (i.e., velocity and acceleration) as
$\dot{y}(t)$ and $\ddot{y}(t)$, respectively. The behavior is then described
by the second-order lag equation:

(2OL) $\ddot{y}(t)=u(t)-ky(t)-d\dot{y}(t).$

Figure 3. Mouse pointer model with spring and damper.

The 2OL-Eq can be interpreted as a spring-mass-damper model, with the mouse
cursor attached to one edge of the screen via a spring and a damper.

An intuitive illustration of these dynamics is given in Figure 3: Assuming
that the mouse pointer is fixed at one edge of the screen via a spring, the
parameters $k$ and $d$ correspond to the stiffness and the damping of this
spring, respectively, and the control value $u$ can be interpreted as the
force acting on the mouse pointer. In particular, the pointer acceleration
$\ddot{y}$ is assumed to be directly proportional to the control $u$ (apart
from the damping and stiffness terms), i.e., (2OL) defines a (linear)
dynamical system of second order. A control flow diagram of the model is shown
in Figure B.1 in the Appendix.

Given a target position $T\in\mathbb{R}$, it can be shown that for the
particular choice $u\equiv kT$, with $k,d>0$, the position $y(t)$ approaches
$T$ for large enough $t$, independent of the initial position $y(0)=y_{0}$ and
initial velocity $\dot{y}(0)=\dot{y}_{0}$ (Hannsgen, 1987). More precisely,
the state $(y,\dot{y})=(T,0)$, i.e., the desired target $T$ is reached and the
velocity is zero, is an equilibrium, meaning that once reached, that state
will forever be maintained. The resulting trajectory, which is uniquely
determined given $y_{0}$, $k$, $d$, and $T$, is often referred to as second-
order lag trajectory, and can be used to model mouse pointing movements
towards a given target $T$. Since the control $u(t)$ is constant in time and
converges towards the equilibrium state, we denote this variant of (2OL) as
2OL-Eq in the following.

Following the general notation of (2), we derive a discrete-time version of
(2OL), with a step size of two milliseconds, i.e., $h=0.002$, which
corresponds to the mouse sensor sampling rate. Considering our example case of
1D pointing tasks, in which the mouse can only be moved horizontally, the
state $x_{n}$ encodes the horizontal position and velocity of the pointer,
denoted by $p_{n}\in\mathbb{R}$ and $v_{n}\in\mathbb{R}$, respectively, i.e.,

(7) $x_{n}=\left(p_{n},v_{n}\right)^{\top}\in\mathbb{R}^{2}.$

Using the forward Euler method777While we could also use the exact solution
here, the (fairly good) approximation via forward Euler yields matrices that
are more suitable for our explanatory purposes., we obtain the matrices $A$
and $B$ for (2) as

(8) $A=\begin{pmatrix}1&h\\\ -hk&1-hd\end{pmatrix},\quad B=\begin{pmatrix}0\\\
h\end{pmatrix}.$

This model can easily be extended to 2D or 3D pointing tasks by augmenting
$x_{n}$ and $u_{n}$ with the respective components for the additional
dimensions.

###  5.1. Analysis of Parameters

Since the 2OL model is an important baseline for mouse pointing dynamics, in
this section, we provide an analysis and intuitive understanding of influence
of the model parameters on model behavior.

While the convergence of 2OL-Eq towards a target $T$ of fixed width $W>0$ can
be easily shown under the assumptions described above, both the time until
this target is reached first (i.e., the time at which the remaining distance
to target is smaller than $W$) and the transient behavior (i.e., specific
characteristics of the trajectory until this time) largely depend on the
parameters $k$ and $d$. In the following, we analyze the effects of the
stiffness $k$ and the damping ratio $\zeta$, which is defined as

(9) $\zeta=\frac{d}{2\sqrt{k}},$

as this is easier to interpret than the actual damping parameter $d$. Note
that given two of the three parameters $k$, $d$, and $\zeta$, the remaining
one (in this case $d$) is uniquely determined by the others and can be easily
computed.

Figure 4.  Position, velocity, and acceleration time series of typical 2OL-Eq
trajectories with target shown as the orange box. Top: Effect of stiffness
parameter $k$ with fixed damping ratio $\zeta=1$ (red dashed: $k=15$, green
dash-dotted: $k=25$, blue dotted: $k=60$). Bottom: Effect of damping ratio
$\zeta$ with fixed stiffness parameter $k=40$ (red dashed: $\zeta=0.5$, green
dash-dotted: $\zeta=1$, blue dotted: $\zeta=2$).
The stiffness parameter $k$ mainly affects the instantaneous initial
acceleration and thus the speed at which the target is approached. The
parameter $\zeta$ determines the relative amount of damping, with $\zeta<1$
(underdamped) resulting in oscillations around the target and $\zeta>1$
(overdamped) leading to trajectories that reach the target later.

Fully described in caption and text.

The position, velocity, and acceleration time series of typical 2OL-Eq
trajectories are shown in Figure 4. The initial and target values stem from an
ID 4 task condition from the Pointing Dynamics Dataset. The most
characteristic feature of 2OL-Eq trajectories is the large positive
acceleration at the beginning of the movement. This is due to the model being
second-order, i.e., the control $u$ is proportional to the acceleration
$\ddot{y}$ (apart from the damping and stiffness terms). The velocity profile
is typically left-skewed, since the deceleration phase is considerably longer
than the acceleration phase. As a consequence, the peak velocity is reached
relatively early during the movement.

Given a constant damping ratio $\zeta$, the stiffness parameter $k$ mainly
affects the initial acceleration, and consequently the peak velocity and the
time at which the target is reached first. As can be seen in the top row of
Figure 4, a large stiffness (blue dotted line) leads to a high initial
acceleration and peak velocity, and thus the target is reached earlier than
with lower stiffness values (red dashed line). For damping ratios $\zeta<1$
(red dashed line in the bottom row of Figure 4), i.e., the damping $d$ is
small compared to the stiffness $k$, oscillations occur in the trajectories,
leading to multiple peaks in velocity and acceleration time series and to
overshooting in the position time series (the so-called underdamped case). For
$\zeta=0$, the pointer does not even converge towards the target, but
oscillates indefinitely (not shown). If $\zeta>1$ (the so-called overdamped
case, blue dotted line in the bottom row of Figure 4), the pointer converges
towards the target slowly, without oscillations. If $\zeta=1$, the trajectory
is critically damped, which means that the pointer reaches (and stays at) the
equilibrium (i.e., the target $T$) in minimum time.

###  5.2. Results of Parameter Fitting

For each combination of participant, task, and direction, we identify the
parameters $\Lambda=(k,\zeta)$ that best explain the corresponding mean
trajectory from the Pointing Dynamics Dataset, using the deterministic
parameter fitting process described in Section 3. The loss function is the SSE
on position, see (4). Note that the results obtained from our parameter
fitting do not exactly match those presented in (Müller et al., 2017), since
we apply a different pre-processing to the experimental user trajectories and
optimize the parameters with respect to the positional error only.

Figure 5. While the position time series visually matches the observed user
data quite well, 2OL-Eq yields a much less symmetric velocity and acceleration
profile during the surge (here: up to 0.36s). This is due to the assumption of
constant equilibrium control, which results in a physically implausible
instantaneous peak acceleration.

[Position, velocity, and acceleration time series of 2OL-Eq vs. user data (ID
4 task)]Position, velocity, and acceleration time series of both the Pointing
Dynamics Dataset (Participant 3, ID 4 (765px distance, 51px width), right
direction) and the corresponding 2OL-Eq simulation trajectory.

The fitted trajectory for a representative ID 4 task condition is shown in
Figure 5. As discussed in (Müller et al., 2017), the main differences between
model and human behavior are the less symmetric velocity profile and the large
initial accelerations produced by the 2OL-Eq model. In particular, the user
trajectories exhibit velocity profiles that are close-to-symmetric and bell-
shaped, at least for the initial ballistic movement towards the target (the
“surge” (Müller et al., 2017)), which is consistent with previous findings
(Morasso, 1981). The differences can be explained with the physical
interpretation of the 2OL-Eq as a spring-mass-damper system: Since $u$ is
constant in this model, as the system is released, the spring instantaneously
accelerates the system with a force that is proportional to the extension of
the spring. Because human muscles cannot build up force instantaneously
(Schmidt et al., 2018), this behavior is not physically plausible.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6. Parameters of 2OL-Eq, optimized for the mean trajectories of all
participants, tasks, and directions, grouped by participants (left) and by ID
(middle), and as overall distribution (right). While the optimal stiffness
parameter $k$ considerably differs between participants, the damping ratio
$\zeta$ is mainly affected by the task ID. Interestingly, all resulting
simulation trajectories are underdamped (mean $\zeta=0.71$), with $\zeta$
increasing as ID increases.

Fully described in caption and text.

In Figure 6, the optimal values of $k$ and $\zeta$ are given for all
participants, tasks, and movement directions, both grouped by participants
(left) and by ID (right). Interestingly, different behavior between individual
users is mainly captured by different stiffness values $k$. Participant 2, for
instance, seems to be characterized best by a large stiffness (between $46.8$
and $107.6$, with mean $78.3$), while the trajectories of participant 9 are
best explained by a considerably lower stiffness (between $10.3$ and $18.3$,
with mean $14.4$) in the 2OL-Eq model. In contrast, the damping ratio $\zeta$
seems to center around $0.7$-$0.76$, independent of the participant. Instead,
it is mainly influenced by the index of difficulty of the task. Lower IDs tend
to result in a lower damping ratio, with trajectories of ID 2 tasks being
considerably more underdamped than others. This is consistent with previous
findings (Guiard, 1993; Bootsma et al., 2004; Müller et al., 2017) and might
be explained by the reciprocal nature of the considered pointing task, where
participants alternately moved between two given targets (which we denote as
initial and target position for a given movement direction) without dwell
time.

In summary, the stiffness $k$ mostly seems to account for movement strategies
that are characteristic of specific users, whereas the damping ratio $\zeta$
mainly differs between indices of task difficulty.

###  5.3. Discussion

The main shortcomings of the 2OL-Eq as a model of mouse pointer movements are
the unrealistically high initial acceleration and the resulting skewed
velocity profile. This is mainly due to the assumption of equilibrium control,
while the literature suggests that the motor control signals are actively
changed during the movement (Bizzi et al., 1992; Georgopoulos et al., 1982;
Todorov, 1998). From a conceptual standpoint, the 2OL-Eq only describes the
passive dynamics of the mouse pointer as a differential equation. It does not
separately model the user’s “brain” or intention in controlling these
dynamics. In particular, it does not describe what the user is trying to
achieve. This can be explained by optimal control models.

##  6\. Pointing as Optimal Open-Loop Control: The Minimum Jerk Model

An elementary model of aimed movements that assumes that users behave
optimally according to an internal cost function is the minimum jerk model by
Flash and Hogan (Flash and Hogan, 1985). This model, which we will refer to as
_MinJerk_ in the following, assumes that the objective of users is to generate
smooth movements by minimizing the jerk of the end-effector, i.e., the time
derivative of the end-effector’s acceleration, while reaching the target
exactly at a prescribed movement time with zero final velocity and
acceleration. Within HCI, this model has been successfully used by Quinn and
Zhai (Quinn and Zhai, 2018) to model the shape of gestures on a word-gesture
keyboard.

The model assumes that the movement is controlled open-loop, and thus cannot
explain how users would correct for disturbances or inaccurate execution.
Similar to 2OL-Eq, the choice of parameters and boundary values already
determines the complete trajectory. However, there are a few important
differences to 2OL-Eq. First, MinJerk does not only require information about
the initial, but also about the final state (as a third-order model, initial
and final position, velocity, and acceleration need to be specified). Second,
the overall movement time, which is denoted by $N_{MJ}$ in the following,
needs to be known in advance. This is in contrast to the discrete-time
formulation of 2OL-Eq.

In discrete time, minimizing jerk corresponds to minimizing the differences
between subsequent accelerations. While in principle, the MinJerk optimal
control problem could be transformed into a closed-loop discrete-time system
similar to (2) (but with time-dependent matrices $A_{n}$ and $B_{n}$, see
(Hoff and Arbib, 1993)), here we make use of the analytical solution of the
original continuous-time problem, and afterwards discretize the resulting 5th-
degree polynomial with respect to time. For arbitrary initial state
$x_{0}=\left(\bar{p}_{0},\bar{v}_{0},\bar{a}_{0}\right)$ and final state
$x_{N_{MJ}}=\left(\bar{p}_{N_{MJ}},\bar{v}_{N_{MJ}},\bar{a}_{N_{MJ}}\right)$,
where the third component is the respective end-effector acceleration, the
(discrete-time) MinJerk trajectory is given by

(10a) $\displaystyle
x_{n}=\sum_{i=0}^{5}c_{i}\left(\frac{n}{N_{MJ}}\right)^{i},$ (10b)
$\displaystyle\begin{pmatrix}c_{0}\\\ c_{1}\\\ c_{2}\\\ c_{3}\\\ c_{4}\\\
c_{5}\end{pmatrix}=\begin{pmatrix}1&0&0&0&0&0\\\ 0&t_{f}&0&0&0&0\\\
0&0&0.5t_{f}^{2}&0&0&0\\\ -10&-6t_{f}&-1.5t_{f}^{2}&10&-4t_{f}&0.5t_{f}^{2}\\\
15&8t_{f}&1.5t_{f}^{2}&-15&7t_{f}&-t_{f}^{2}\\\
-6&-3t_{f}&-0.5t_{f}^{2}&6&-3t_{f}&0.5t_{f}^{2}\end{pmatrix}\begin{pmatrix}\bar{p}_{0}\\\
\bar{v}_{0}\\\ \bar{a}_{0}\\\ \bar{p}_{N_{MJ}}\\\ \bar{v}_{N_{MJ}}\\\
\bar{a}_{N_{MJ}}\end{pmatrix},$

where $t_{f}=N_{MJ}h$ denotes the final time, and $h$ is the fixed time
interval between two consecutive time steps.

###  6.1. Extension to Complete Trajectories

Figure 7. Visually, MinJerk trajectories show a relatively good fit to
observed user data, as long as the duration parameter $N_{MJ}$ is chosen to
cover only the surge (black solid line; here, $N_{MJ}=223$ corresponds to
0.446s). However, the MinJerk model predicts zero velocity at the end of the
surge, which results in a considerably worse overall fit of user trajectories
that exhibit clear submovements. If $N_{MJ}$ is set as total movement duration
$N$, the resulting simulation trajectory (dash-dotted blue line) does not fit
the data at all.

[Position, velocity, and acceleration time series of MinJerk vs. user data (ID
4 task)]Position, velocity, and acceleration time series of both the Pointing
Dynamics Dataset (Participant 3, ID 4 (765px distance, 51px width), right
direction) and the corresponding MinJerk simulation trajectories for
$N_{MJ}<N$ and $N_{MJ}=N$.

The MinJerk model has been derived from data of an experiment that did not
involve any corrective submovements (Flash and Hogan, 1985). However, in mouse
pointing tasks with large ID, submovements occur regularly. If MinJerk is used
for modeling of the entire movement, i.e., until the final time step $N$, the
fit is thus very poor (see blue dash-dotted lines in Figure 7). Instead of a
quick movement towards the target with corrective submovements, as in our
data, the model predicts a slow, smooth movement, reaching the target only at
the final time.

A much better fit is obtained by using MinJerk only for the first, rapid
movement towards the target (the “surge”), and assuming that the pointer does
not move afterwards. To this end, we define the (extended) MinJerk trajectory
as follows. For $n\leq N_{MJ}$, $x_{n}$ corresponds to the minimum jerk
polynomial (10) with $\bar{p}_{0}$, $\bar{v}_{0}$, and $\bar{a}_{0}$ taken
from user data, $\bar{p}_{N_{MJ}}=T$, and
$\bar{v}_{N_{MJ}}=\bar{a}_{N_{MJ}}=0$. For $n>N_{MJ}$, the trajectory is
constantly extended by the final state of the MinJerk polynomial, i.e.,
$x_{n}=x_{N_{MJ}}=(T,0,0)^{\top}$.

Figure 8.  Position, velocity, and acceleration time series of different
MinJerk trajectories with target shown as the orange box. The parameter
$N_{MJ}$ determines the end time of the symmetric and smooth jerk-minimizing
movement, after which the trajectory is constantly extended by its final
position value and zero velocity and acceleration (red dashed: $N_{MJ}=50$,
green dash-dotted: $N_{MJ}=75$, blue dotted: $N_{MJ}=N=100$).

Fully described in caption and text.

The effect of the parameter $N_{MJ}$ in our MinJerk model is shown in Figure
8. Varying this parameter allows to model variable peak velocities and
accelerations. However, with $N_{MJ}$ being a pure scaling parameter of the
trajectory, the velocity profile is always bell-shaped and the acceleration
profile N-shaped. Moreover, the target center is reached at time step
$n=N_{MJ}$ by definition, i.e., the model naturally cannot account for
corrections that typically occur after the surge. As illustrated in Figure 7
for the same ID 4 task as above (black solid lines), this can result in a
considerably worse overall fit of the trajectory, at least for movements that
consist of several submovements.

###  6.2. Results of Parameter Fitting

Similar to the parameter fitting process for 2OL-Eq, we identify the optimal
duration parameter $N_{MJ}$ with respect to positional SSE for each
participant, task, and movement direction, using the respective mean
trajectory as our reference. To improve convergence of the used optimization
algorithm, we relax this parameter by allowing for continuous values of
$N_{MJ}$ in (10) and compute the discrete-time MinJerk states $x_{n}$ for
$n\leq\left\lceil N_{MJ}\right\rceil$.

(a)

(b)

(c)

Figure 9. Duration parameter $N_{MJ}$ of MinJerk, optimized for the mean
trajectories of all participants, tasks, and directions. Left: Absolute
parameter values, grouped by participants. Middle: Parameter values relative
to the total movement duration $N$, grouped by target width. Right: Parameter
values relative to the total movement duration $N$, grouped by distance to
target.

Fully described in caption and text.

The optimal values of $N_{MJ}$ grouped by participants are shown in the left
plot of Figure 9. Different users can be characterized by different values for
$N_{MJ}$. Interestingly, the user-specific effects match those observed for
the stiffness parameter $k$ in the 2OL-Eq model fairly well. Indeed, there is
an inverse-linear relationship between $k$ and $N_{MJ}$, as illustrated in
Figure B.2 in the Appendix. Comparing the effects of $k$ (Figure 4, top row)
and $N_{MJ}$ (Figure 8) on the respective model trajectories, however, this is
not very surprising. Both a higher stiffness $k$ and a lower surge duration
$N_{MJ}$ lead to a faster movement that reaches the target earlier (note that
participant 9 moved considerably slower (average duration: 1.7s) than the rest
of the participants (average duration: 0.95s)). Differences between the
effects of these two parameters are mainly related to the initial
acceleration, which scales with $k$ in 2OL-Eq, but is always fixed in MinJerk,
and to the skewness of the velocity profile, which is only affected by $k$.

In contrast, the average surge duration does not differ much between different
task IDs (not shown). However, the relative time spent in the surge, i.e.,
$\nicefrac{{N_{MJ}}}{{N}}$, clearly increases as the target width increases,
while it is unaffected by the distance between initial and target center (see
middle and right plot of Figure 9). This is consistent with previous findings,
which suggest that the skewness of the velocity profile and thus the duration
of corrective submovements is mainly determined by the target size, while the
distance to target has a much smaller effect on the relative duration of the
initial ballistic movement (Thompson et al., 2007; Bohan et al., 2003; Müller
et al., 2017).

###  6.3. Discussion

The minimum jerk model can explain the shape of the initial ballistic movement
(the surge) towards the target very well. However, both initial and terminal
conditions need to be known in advance. The same holds for the overall
movement time, unless it is identified through a parameter fitting process,
using experimentally observed user data. It should also be noted that it is
difficult to explain conceptually why users should aim to minimize the jerk of
the movement (see, e.g., Harris and Wolpert (Harris and Wolpert, 1998)).

The main limitation of the model is that it is a pure kinematic model. That
is, the trajectory of the mouse pointer is uniquely defined given the initial
and terminal conditions. The underlying reasons for the movement, such as the
acting forces, are not explained. In particular, the model does not involve
any explanation of the underlying biomechanics of the user, not even as a
point-mass model such as 2OL-Eq. Due to its deterministic nature, it cannot
account for the between-trial variability typically observed in user movements
(see red dashed lines in Figure 7). Further, as an open-loop model, the
movement trajectory is completely specified at the beginning of the movement,
and in its standard form, the model cannot react to perturbations or
inaccuracies in the movement. In order to explain how users react to visual or
proprioceptive feedback, models based on optimal feedback control are
required.

##  7\. Pointing as Optimal Feedback Control: The LQR

In general, the optimal control problem given by (3) is very difficult to
solve, since no solution method is known that guarantees convergence towards
the global optimum without imposing (fairly strong) assumptions on the costs
and system dynamics (e.g., convexity, continuous differentiability, etc.). One
important subclass of problems where an analytic solution method exists is the
case of linear dynamics and quadratic costs. The solution in this case is
given by the Linear-Quadratic Regulator (LQR).

These optimal control problems usually are of the following form:

(11c) $\displaystyle\begin{gathered}\textsl{Minimize}\quad
J_{N}(x,u)=\sum_{n=0}^{N}x_{n}^{\top}Q_{n}x_{n}+\sum_{n=0}^{N-1}u_{n}^{\top}R_{n}u_{n}\\\
\textsl{with respect to
}u=(u_{n})_{n\in\\{0,\dots,N-1\\}}\subset\mathbb{R}^{m},\end{gathered}$ where
$x=(x_{n})_{n\in\\{0,\dots,N\\}}\subset\mathbb{R}^{k}$ satisfies (11f)
$\displaystyle\begin{gathered}x_{n+1}=Ax_{n}+Bu_{n},\quad
n\in\\{0,\dots,N-1\\},\\\ x_{0}=\bar{x}_{0}\end{gathered}$ for some given
initial state $\bar{x}_{0}\in\mathbb{R}^{k}$.

As before, $x_{n}$ is the state of the human-computer system, $u_{n}$ is the
control (e.g., muscle excitation), the matrix $A$ describes the dynamics of
the human-computer system if no control is applied, i.e., $u\equiv 0$, and $B$
describes how the control influences the system (e.g., the force generated by
muscles). The matrices $Q_{n}$ and $R_{n}$ can be interpreted as coefficients
or weights for the state and control costs, respectively, where, e.g., the
former formalizes that users aim to reach the target and the latter formalizes
that they aim to do so with minimal effort. Note that in our case, the
controls $u_{n}$ are one-dimensional ($m=1$), i.e., the matrix $R_{n}$ only
consists of a single entry.

Regarding the optimal control framework for Human-Computer Interaction, the
minimization in Equation (11c) corresponds to the Human Controller block in
Figure 1(b), and Equation (11f) corresponds to the Human-Computer System
Dynamics block. The Feedback and Human Observer blocks in Figure 1(b) are
assumed trivial, i.e., observation and internal state estimate both equal the
actual state of the system ($y_{n+1}=\hat{x}_{n+1}=x_{n+1}$). This is clearly
different from the open-loop case depicted in Figure 1(a), where the controls
are independent from the state estimates.

In the following, we use the muscle model and the cost function that have been
used by Todorov (Todorov, 2005) in the case of the (stochastic, and therefore
significantly more complex) Linear-Quadratic Gaussian regulator, which we
discuss in Section 8. To understand the stochastic case, it is useful to first
consider the deterministic case, which we introduce in this section. The
simplified second-order muscle model that we use has been proposed by van der
Helm (van der Helm and Rozendaal, 2000), and obtains control signals as input
$u(t)$ and yields the forces applied to the end-effector $F(t)$ as output:

(12) $F(t)=u(t)-\tau_{1}\tau_{2}\ddot{F}(t)-(\tau_{1}+\tau_{2})\dot{F}(t),$

with time constants $\tau_{1}$, $\tau_{2}>0$. Throughout this section, we use
$\tau_{1}=\tau_{2}=0.04$.

A discrete-time approximation of these dynamics is obtained by the Forward
Euler method with time interval $h>0$,

(15)
$\displaystyle\begin{gathered}f_{n+1}=f_{n}+\frac{h}{\tau_{2}}\left(g_{n}-f_{n}\right),\\\
g_{n+1}=g_{n}+\frac{h}{\tau_{1}}\left(u_{n}-g_{n}\right),\end{gathered}$

where $f_{n}$ and $g_{n}$ denote the muscle activation (corresponding to
force) and excitation at time step $n$, respectively. Recall that $h=0.002$
corresponds to the two millisecond sampling rate of the mouse sensor.
Following Todorov (Todorov, 2005), we assume a unit mass of 1 kg of the hand-
mouse system.

Combining these muscle dynamics with a second-order integrator, this leads to
the following system dynamics matrices:

(16) $A=\begin{pmatrix}1&h&0&0&0\\\ 0&1&h&0&0\\\
0&0&1-\frac{h}{\tau_{2}}&\frac{h}{\tau_{2}}&0\\\
0&0&0&1-\frac{h}{\tau_{1}}&0\\\ 0&0&0&0&1\end{pmatrix},\quad
B=\begin{pmatrix}0\\\ 0\\\ 0\\\ \frac{h}{\tau_{1}}\\\ 0\end{pmatrix}.$

Here, the state $x_{n}=(p_{n},v_{n},f_{n},g_{n},T)^{\top}\in\mathbb{R}^{5}$
(i.e., $k=5$ in (11f)) consists of the pointer position $p_{n}$ and velocity
$v_{n}$, as well as muscle force $f_{n}$ and muscle excitation $g_{n}$.
Moreover, the fixed target $T$ is included in the state for technical reasons.
Note that, in contrast to the 2OL-Eq model, the controls $u_{n}$ do not equal
a fixed, target-dependent value (that is, we do not prescribe pure equilibrium
control), but are chosen to minimize the cost $J_{N}$.

In principle, the cost matrices $Q_{n}$ and $R_{n}$ can be chosen freely.
Based on (Todorov, 2005), we derive them from the following assumptions:

  * •

Users aim to minimize the distance between pointer and target.

  * •

Users aim to stay inside the target after reaching it.

  * •

Users aim to minimize the effort required to fulfill the task.

Ideally, no distance costs should occur within the target, which is assumed to
be a box of width $W$. However, a fundamental limitation of the LQR setting is
that cost terms need to be quadratic in the states and controls. Therefore,
costs that are zero everywhere within the target are unfortunately infeasible
within the LQR setting. To approximate such costs, we construct the distance
costs such that we have lower costs inside the target and higher costs
outside. In particular, we penalize the squared remaining Euclidean distance
between the end-effector position $p_{n}$ and the desired target position $T$,
which is given as

(17) $D_{n}^{2}=|p_{n}-T|^{2}.$

To create an incentive to keep the end-effector inside the target once it is
reached, the squared velocity $v_{n}^{2}$ and the squared force $f_{n}^{2}$
(which can be interpreted as acceleration, since unit mass is assumed) are
penalized as well, weighted with $\omega_{v},\omega_{f}\geq 0$. All these cost
terms are quadratic with respect to the state $x_{n}$, i.e., a positive semi-
definite matrix $Q_{n}$ can be found such that

(18)
$x_{n}^{\top}Q_{n}x_{n}=D_{n}^{2}+\omega_{v}v_{n}^{2}+\omega_{f}f_{n}^{2}.$

At the same time as minimizing the distance to the target, we assume that
users aim at minimizing their effort. This assumption is well-established in
motor control theory, motivated by both neuroscientific findings and
mathematical requirements (Todorov and Jordan, 2002; Li and Todorov, 2004;
Guigon et al., 2007).888In optimal control problems, penalizing the controls
$u$ acts as a regularizer, i.e., it constrains the subspace of optimal
solutions, which often results in a unique optimal solution. Moreover, we
assume that the effort cost matrices are constant in time, i.e., $R_{n}=R$
holds for all $n\in\\{0,\dots,N-1\\}$, and normalized with respect to the
duration of the movement for better comparability between conditions. In
particular, we choose

(19) $R=\frac{\omega_{r}}{N-1},$

where the weight parameter $\omega_{r}>0$ determines how much effort the user
is willing to invest to reduce the distance to the target more quickly.

Figure 10. If state costs are only applied at the final time step $N$ (blue
dashed line), the LQR cannot capture the correction phase, which is highly
pronounced for tasks with sufficiently small target. With continuous state
costs (black solid line), LQR trajectories visually exhibit a good fit,
although the peak acceleration is slightly higher. In contrast to MinJerk, the
end of the surge phase does not need to be prescribed, but emerges implicitly
from the model.

[Position, velocity, and acceleration time series of the LQR vs. user data (ID
4 task)]Position, velocity, and acceleration time series of both the Pointing
Dynamics Dataset (Participant 3, ID 4 (765px distance, 51px width), right
direction) and the corresponding LQR simulation trajectories, each with
continuous state costs and with terminal state costs only. The control time
series of both LQR variants are shown as well.

While in (Todorov, 2005), the state costs (18) were only applied at the final
time $n=N$, which must be known in advance, in the deterministic LQR case
without neuromotor noise, this does not work. As can be seen in Figure 10, for
tasks with fairly small targets that require a considerable correction phase,
the resulting LQR trajectory (dashed blue lines) does not resemble the
observed user behavior at all. The main reason for this is that distance costs
that only occur at a single time step do not create an incentive to reach this
state earlier than necessary. The constant incurrence of control costs adds to
this, resulting in an optimal policy that chooses relatively small controls
during most of the movement, while shortly before the final time step $N$,
larger controls are applied to reach the target “just in time” with low
velocity and acceleration.

This problem can be addressed in two ways. Either neuromotor noise can be
included, as it is done in (Todorov, 2005) and in the next section. Or, in the
deterministic case, behavior more similar to humans can be achieved by
assuming that the state costs from Equation (18) are applied during the entire
movement. This clearly creates an incentive to move the pointer towards the
target from the start, in order to reduce the sum of the distance costs.

In summary, the objective function of our LQR model is given by

(20)
$J_{N}^{\text{(LQR)}}(x,u)=\sum_{n=0}^{N}\left(D_{n}^{2}+\omega_{v}v_{n}^{2}+\omega_{f}f_{n}^{2}\right)+\frac{\omega_{r}}{N-1}\sum_{n=0}^{N-1}u_{n}^{2}.$

We assume that the user computes the _optimal_ control, denoted by
$u^{*}_{n}$, in a feedback manner, based on the current state (i.e., the model
is closed-loop). It has been proven that for these kinds of problems the
optimal control $u^{*}_{n}$ depends linearly on the state (Dorato and Levis,
1971), i.e.,

(21) $u_{n}^{*}=\pi(x_{n})=-L_{n}x_{n}$

holds for some uniquely determined feedback gain matrices $L_{n}$. Recall that
in our case, the control $u_{n}$ is one-dimensional. Since the state $x_{n}$
is a vector in $\mathbb{R}^{5}$, $L_{n}$ is thus a $1\times 5$ matrix.

Given the matrices $A$ and $B$ as well as the cost function
$J_{N}^{\text{(LQR)}}$, these feedback gain matrices can be computed once
before movement onset (i.e., at the planning stage) by iteratively solving the
corresponding _Discrete Riccati Equation_ (see Appendix; details are given in
(Todorov, 1998, Theorem 7)). This results in a very fast on-line computation
of the optimal controls, yet taking into account the most recent state
observations.

###  7.1. Analysis of Parameters

Figure 11.  Position, velocity, acceleration, and control time series of
typical LQR trajectories with target shown as the orange box. Top: Effect of
velocity cost weight $\omega_{v}$ with fixed cost weights $\omega_{f}=0$,
$\omega_{r}=~{}$5e-3 (red dashed: $\omega_{v}=0.01$, green dash-dotted:
$\omega_{v}=0.05$, blue dotted: $\omega_{v}=0.1$). Middle: Effect of force
cost weight $\omega_{f}$ with fixed cost weights $\omega_{v}=0$,
$\omega_{r}=~{}$5e-3 (red dashed: $\omega_{f}=~{}$1e-4, green dash-dotted:
$\omega_{f}=~{}$1e-3, blue dotted: $\omega_{f}=~{}$1e-2). Bottom: Effect of
effort cost weight $\omega_{r}$ with fixed cost weights $\omega_{v}=0.01$,
$\omega_{f}=~{}$1e-4 (red dashed: $\omega_{r}=~{}$5e-4, green dash-dotted:
$\omega_{r}=~{}$5e-3, $\omega_{r}=0.05$).

Fully described in caption and text.

In Figure 11, the individual effects of the cost function weights
$\omega_{v}$, $\omega_{f}$, and $\omega_{r}$ are shown. A higher velocity cost
weight $\omega_{v}$ (top plots, blue dashed lines) results in a lower peak
velocity, as expected, since velocity is penalized quadratically. Keeping the
remaining parameters constant, this leads to a less symmetric velocity
profile, as higher velocities towards the end of the movement are necessary to
compensate for the lower peak. Moreover, the target is reached later. Note
that this only occurs as long as the velocity cost weight is not dominant.
Otherwise, for large enough $\omega_{v}$, there is no incentive to reach the
target at all.

Similar effects (peak values, symmetry) can be observed in the acceleration
profile for the force cost weight $\omega_{f}$ (since the forces are applied
to a unit mass and thus can be interpreted as acceleration) and in the control
profile for the effort cost weight $\omega_{r}$. Moreover, large force cost
weights lead to a constant, positive velocity at the end of the movement
(middle plots, velocity profile), i.e., the pointer moves across the target
instead of staying inside (middle plots, position profile). In contrast, the
magnitude of the effort cost weight $\omega_{r}$ mainly affects the duration
of the surge (bottom plots, position profile), while the target is still
reached with relatively low velocity and acceleration for moderate values of
$\omega_{r}$.

###  7.2. Results of Parameter Fitting

Analogously to the parameter fitting process for the 2OL-Eq and MinJerk
models, we identify optimal values for the cost weights $\omega_{v}$,
$\omega_{f}$, and $\omega_{r}$ collected in the paramter vector $\Lambda$ for
each mean trajectory. Since these parameters only define the objective
function $J_{N}^{\text{(LQR)}}$ of the OCP, and the system matrices $A$ and
$B$ are uniquely determined given the above fixed values of $m$, $\tau_{1}$,
and $\tau_{2}$, the parameter fitting for the LQR can be regarded as an
inverse optimal control problem (Jin et al., 2021).999Including the system
dynamics parameters $m$, $\tau_{1}$, and $\tau_{2}$ in the parameter fitting
process did not improve the fit to observed user behavior. As usual, we use
the bi-level approach described in Section 4.2, i.e., at each iteration of the
parameter fitting method, the OCP subject to the parameter vector $\Lambda$ is
solved as described above.

####  7.2.1. Qualitative Results

The resulting optimal LQR trajectory (for the same ID 4 task considered for
the previous models) is shown as the black solid line in Figure 10. The
position time series is approximated fairly well at first glance, even better
than by the MinJerk trajectory (see Figure 20; quantitative comparisons
between all models are given in Section 11.2). Peak velocity and acceleration
are also very close to the values of the respective user trajectory, albeit
the maximum acceleration is higher and the timing of the minimum and maximum
acceleration does not match exactly. However, it is important to note that the
duration of the surge phase was not explicitly built into the LQR
model101010In theory, it would be possible to include prior knowledge about
submovements and thus induce time-dependent behavior by choosing the time-
dependent cost matrices $Q_{n}$ and $R_{n}$ appropriately., but emerges
naturally from the interplay of the optimal parameters $\omega_{v}$,
$\omega_{f}$, and $\omega_{r}$. In contrast, in the MinJerk model, the
duration of the surge phase needs to be either known in advance or determined
using the parameter fitting process.

####  7.2.2. Quantitative Results

(a)

(b)

(c)

(d)

(e)

(f)

Figure 12. Parameters of the LQR model, optimized for the mean trajectories of
all participants, tasks, and directions, grouped by participants (left) and by
ID (right). For better visibility, the optimal values of $\omega_{f}$ and
$\omega_{r}$ are plotted on a logarithmic scale.

Fully described in caption and text.

All optimal parameter values, grouped by both user and task ID, are shown in
Figure 12 (note the logarithmic scale for $\omega_{f}$ and $\omega_{r}$).

Each of the three parameters exhibit a large between-user variability.
Interestingly, the between-user effects show similar trends and
characteristics for all cost weights. For example, trajectories of participant
9 are characterized by very large weights (i.e., relatively small distance
costs, since the weights affect every term in the cost function but the
distance costs), while the behavior of participant 2 is explained best by
considerably lower weights (i.e., relatively high distance costs). Moreover,
these characteristic differences in the behavior of individual users are
similar to those identified for the surge duration $N_{MJ}$ in MinJerk (see
Figure 9 (left)). In contrast, the task difficulty does not seem to have a
clear effect on the optimal values of the cost weights. Hence, our findings
suggest that the identified optimal parameters rather encode user-specific
characteristics than the task under consideration.

###  7.3. Discussion

The LQR combines beneficial features of both 2OL-Eq (movement duration emerges
from the model) and MinJerk (smooth movements with “close-to-bell-shaped”
velocity profiles, as observed in many user studies). With a continuous
penalization of remaining distance to target, velocity, and force during the
entire movement, the LQR model is able to explain average user behavior in
terms of position, velocity, and acceleration profiles. However, whether the
modified cost terms are plausible from a biomechanical and neuroscientific
perspective is debatable (Berret et al., 2011; Richardson and Flash, 2002).
Moreover, as the MinJerk model, due to its deterministic nature, the LQR model
cannot account for the considerable between-trial variability, which is
typically observed in user movements.

In the next section, we will thus introduce signal-dependent control noise,
which in combination with continuous effort costs and only terminal distance,
velocity, and force costs, allows to replicate user behavior similarly well as
the proposed LQR variant, with the additional benefit of obtaining a
distribution of optimal trajectories.

##  8\. Pointing as Optimal Feedback Control Subject to Signal-Dependent Motor
Noise: The LQG

While the LQR modelvisually captures typical average user behavior relatively
well (see Figure 10; for a quantitative comparison between different models,
see Section 11), it has one major drawback: the proposed optimal trajectory is
necessarily deterministic and thus cannot account for the large variability
observed among multiple trials of the same participant for the same task
(Dhawale et al., 2017; Todorov and Jordan, 2002).

An extension that allows to consider different sources of variability, which
might occur during the sensorimotor control loop, is the Linear-Quadratic
Gaussian Regulator (LQG) (Hoff, 1992; Loeb et al., 1990; Todorov, 2005). It
belongs to the class of Stochastic Optimal Feedback Control (SOFC) models, as
it takes into account the variance resulting from various noise terms.

In the following, we will present and analyze the LQG model introduced by
Todorov for reaching movements (Todorov, 2005). In comparison to the LQR model
from Section 7, this model includes signal-dependent Gaussian control noise
and an observation model with additive Gaussian noise. Moreover, all state
costs (i.e., distance, velocity, and force costs) are only applied at the
final time step $N$.

###  8.1. Linear-Quadratic Gaussian Regulator (LQG)

Starting from the linear-quadratic control problem (11), we extend the
deterministic LQR model presented in the previous section with stochastic
noise terms. In principle, there are two main types of noise terms that can be
included in the system dynamics: additive noise and multiplicative
noise.111111In this paper, we only consider white noise, i.e., the noise terms
that occur at different time steps are assumed to be independent. However, the
LQG model could be extended to incorporate temporarily correlated noise by
augmenting the state space accordingly (for details, see (Todorov, 1998,
Section 3.4.3)). While the former introduces the same level, or variability,
to the system at each time step, the variance of multiplicative noise depends
on the system variables themselves. In the considered case of signal-dependent
(multiplicative) control noise, a larger magnitude of the applied control
$u_{n}$ thus results in a higher uncertainty about the subsequent state
$x_{n+1}$. From a neuroscientific perspective, this dependency can be
justified by the empirically observed effect of neural motor commands on the
variance of the resulting motor-neuronal firing (Sutton and Sykes, 1967;
Schmidt et al., 1979; Clamann, 1969; Jones et al., 2002). In particular,
signal-dependent noise can account for well-established phenomena such as
cosine tuning, muscle synergies, smooth movements, and the trade-off between
speed (or duration) and the end-point accuracy of a movement (Todorov, 2002;
Harris and Wolpert, 1998; Todorov, 2004).

Introducing the control noise level $\sigma_{u}>0$ and the sequence
$(\eta_{n})_{n\in\\{0,\dots,N-1\\}}$ of (univariate) Gaussian random variables
$\eta_{n}\sim\mathcal{N}(0;1)$, the resulting discrete-time system dynamics
can be written as

(24)
$\displaystyle\begin{gathered}x_{n+1}=Ax_{n}+(1+\sigma_{u}\eta_{n})Bu_{n},\quad
n\in\\{0,\dots,N-1\\},\\\
x_{0}\sim\mathcal{N}(\bar{x}_{0},\Sigma_{0}),\end{gathered}$

where the initial state $x_{0}$ is drawn from a multivariate Gaussian
distribution with given mean $\bar{x}_{0}$ and covariance matrix $\Sigma_{0}$.
The specific values for $\bar{x}_{0}$ and $\Sigma_{0}$ can be extracted from
data, see Section 8.2. Following (Todorov, 2005), we do not include additive
control noise, albeit this would be easily possible in the proposed framework.

As in the previous section, the state $x_{n}\in\mathbb{R}^{5}$ contains
pointer position, velocity, force, and muscle excitation, as well as the
target position. We assume that the controller, which needs to decide for a
control $u_{n}\in\mathbb{R}$ at each time step $n\in\\{0,\dots,N-1\\}$, does
not have complete access to the current state of the system $x_{n}$. The main
reason is that usually, not all information stored in $x_{n}$ is observable.
This means that the control $u_{n}$ at time step $n\in\\{0,\dots,N-1\\}$ may
depend only on some observation $y_{n}\in\mathbb{R}^{l}$ ($l\in\mathbb{N}$),
but not on the true state $x_{n}$. Moreover, the observation can be noisy,
i.e., errors during observation may occur. More specifically, to maintain
linear model dynamics, we assume that the (noisy) observation $y_{n}$ is
linear in $x_{n}$, i.e.,

(25) $y_{n}=H_{n}x_{n}+G\xi_{n}.$

For simplicity, we assume that position, velocity, and force values can be
observed immediately in global coordinates, while muscle excitation and target
position cannot be directly observed, i.e.,

(26) $H_{n}=\begin{pmatrix}1&0&0&0&0\\\ 0&1&0&0&0\\\ 0&0&1&0&0\end{pmatrix},$

which results in

(27) $H_{n}x_{n}=\left(p_{n},v_{n},f_{n}\right)^{\top}.$

The error that occurs during observation is modeled by the additive noise term
$G\xi_{n}$. The $3\times 3$ matrix $G$ determines the individual noise levels.
As in (Todorov, 2005), it is given by

(28) $G=\sigma_{s}~{}diag(0.02,0.2,1),$

where $\sigma_{s}>0$ is a scaling parameter and the constants were chosen in
order to reflect the different magnitudes between position, velocity, and
force (i.e., acceleration). The vector $\xi_{n}$ is a three-dimensional
Gaussian random variable, i.e., $\xi_{n}\sim\mathcal{N}(0;I_{3})$, where
$I_{3}$ denotes the $3\times 3$ identity matrix. In particular, the
observation $y_{n}$ is a three-dimensional vector.

The objective function is related to the one used in the LQR model, (20),
i.e., a weighted combination of distance, velocity, and force costs, plus
effort costs. The two differences are: (i) distance, velocity, and force costs
incur only at the final time step $N$, while only the effort costs incur
continuously throughout the movement, and (ii) since state and control are not
deterministic due to noise, we use the expected value of these terms, denoted
by $\mathbb{E}[\cdot]$. The objective function is thus given by

(29)
$J_{N}^{\text{(LQG)}}(x,u)=\mathbb{E}\left[D_{N}^{2}+\omega_{v}v_{N}^{2}+\omega_{f}f_{n}^{2}+\frac{\omega_{r}}{N-1}\left(\sum_{n=0}^{N-1}u_{n}^{2}\right)\right].$

In total, this results in the following stochastic optimal control problem:

(30c) $\displaystyle\begin{gathered}\textsl{Minimize}\quad
J_{N}^{\text{(LQG)}}(x,u)\\\ \textsl{with respect to
}u=(u_{n})_{n\in\\{0,\dots,N-1\\}}\subset\mathbb{R}^{m},\end{gathered}$ where
$y=(y_{n})_{n\in\\{0,\dots,N-1\\}}\subset\mathbb{R}^{l}$ satisfies (30e)
$\displaystyle\begin{gathered}y_{n}=H_{n}x_{n}+G\xi_{n},\quad
n\in\\{0,\dots,N-1\\},\end{gathered}$ and
$x=(x_{n})_{n\in\\{0,\dots,N\\}}\subset\mathbb{R}^{k}$ satisfies (30h)
$\displaystyle\begin{gathered}x_{n+1}=Ax_{n}+(1+\sigma_{u}\eta_{n})Bu_{n},\quad
n\in\\{0,\dots,N-1\\},\\\
x_{0}\sim\mathcal{N}(\bar{x}_{0},\Sigma_{0}).\end{gathered}$

We specifically use the dimensions $m$, $l$, and $k$ throughout to underline
the generality and easy expandability of the LQG model. In our case, we recall
that $m=1$, $l=3$, and $k=5$.

Since the true state $x_{n}$ at time step $n$ is not available to the
controller, it needs to compute internal state estimates $\hat{x}_{n}$ based
on the information available. Under the LQG assumptions (linear dynamics,
quadratic costs, Gaussian noise), the state estimates $\hat{x}_{n}$ can be
computed using a linear estimator:

(31) $\hat{x}_{n+1}=A\hat{x}_{n}+Bu_{n}+K_{n}(y_{n}-H_{n}\hat{x}_{n}),\quad
n\in\\{0,\dots,N-1\\}.$

Here, $y_{n}$ and $H_{n}\hat{x}_{n}$ denote the obtained and the expected
sensory input at time step $n$, respectively (recall the Human Observer block
in Figure 1(b)), where initially, at $n=0$, we set $\hat{x}_{0}=\bar{x}_{0}$.
The matrix $K_{n}$, which is to be determined, specifies to which degree the
observed differences between these quantities should be taken into account for
the computation of the subsequent state estimate $\hat{x}_{n+1}$. As a trivial
example, consider $K_{n}=0$, which corresponds to the open-loop case, i.e., no
sensory information is used by the controller. When it comes to choosing a
“good” matrix $K_{n}$, the above LQG assumptions allow to analytically derive
the matrix that is optimal in the sense that it minimizes the objective
function $J_{N}^{\text{(LQG)}}$. The resulting optimal estimator is known as
the Kalman Filter.121212In the case of state-dependent observation noise,
i.e., if $G$ is replaced by $G_{n}(x_{n})$, the Kalman Filter is still optimal
among all linear estimators (Todorov, 1998).

Moreover, it can be shown that, similar to the LQR case, the optimal closed-
loop solution $u^{*}=(u_{n}^{*})_{n\in\\{0,\dots,N-1\\}}$ is linear in the
state estimate $\hat{x}_{n}$, i.e., there exist unique feedback gain matrices
$L_{n}$, $n\in\\{0,\dots,N-1\\}$ such that

(32) $u_{n}^{*}=\pi(\hat{x}_{n})=-L_{n}\hat{x}_{n}$

holds. As in the previous section, $L_{n}$ is a $1\times 5$ matrix. For more
details, we refer to (Todorov, 2005).

It is important to note that in the considered case of signal-dependent
control noise, the feedback gain matrices $L_{n}$ and the Kalman gain matrices
$K_{n}$ non-trivially depend on each other. More precisely, each feedback gain
matrix $L_{n}$ depends on the subsequent Kalman filter matrices
$(K_{i})_{i\in\\{n+1,\dots,N-1\\}}$, and each Kalman filter matrix $K_{n}$
depends on the previous feedback gain matrices
$(L_{i})_{i\in\\{0,\dots,n-1\\}}$. Thus, the feedback gain matrices $L_{n}$
can only be computed backward in time, starting at the final time step $n=N$,
whereas the Kalman filter matrices $K_{n}$ can only be computed forward in
time, starting at the first time step $n=0$. In particular, one must be given
in order to optimally determine the other. Fortunately, iterating alternately
between the two optimizations results in a coordinate descent algorithm, which
can be shown to converge towards a (local) minimum (Todorov,
2005).131313Numerically, it has been shown that this algorithm even converges
towards the global minimum (Todorov, 2005). This iterative computation of the
matrices $K_{n}$ and $L_{n}$ can be done offline, i.e., before movement onset.

As suggested by Todorov (Todorov, 2005), we therefore initialize the Kalman
gain matrices to zero, i.e., we first compute the optimal open-loop control
strategy matrices $L_{n}$, and iteratively compute $K_{n}$ and $L_{n}$ until
the resulting objective function value converges towards its optimum. In
particular, we terminate the iterative optimization procedure in the $i$-th
step if the relative improvement of $J_{N}^{\text{(LQG)}}$ falls below a
predefined threshold $\epsilon_{J}>0$, i.e.,

(33)
$\left|\frac{\left(J_{N}^{\text{(LQG)}}\right)_{i+1}-\left(J_{N}^{\text{(LQG)}}\right)_{i}}{\left(J_{N}^{\text{(LQG)}}\right)_{i}}\right|\leq\epsilon_{J},$

or after a maximum number of iterations (set to 20; we found that this value
sufficed to obtain parameters that did not change considerably afterwards). In
our implementation, we use $\epsilon_{J}=~{}$1e-3.

We do not make use of the adaptive Kalman filter, as it was proposed for an
LQG model designed for via-point movements (Todorov, 1998). The adaptive
Kalman filter computes Kalman gain matrices $K_{n}$ that explicitly depend on
the observations $(y_{i})_{i\in\\{0,\dots,n-1\\}}$ received so far, i.e., this
needs to be done online. Moreover, since these observations are stochastic,
the matrices $K_{n}$ (and thus the expected behavior) differ between several
runs. In addition, it was empirically observed that the adaptive Kalman filter
might be unstable (Todorov, 2005). In contrast, the non-adaptive Kalman filter
from (Todorov, 2005) computes matrices $K_{n}$ that only depend on information
available before movement onset. The resulting a priori expectations over the
state sequences can be directly used in the loss function of the parameter
fitting process to compare the expected outcome of different parameter
vectors.

In principle, using an adaptive filter only during runtime (with $L_{n}$
optimized with respect to the non-adaptive filter) could further reduce the
expected total costs. However, this effect has shown to be minor (Todorov,
2005), which is why we refrain from using an adaptive filter at all.

###  8.2. Analysis of Parameters

The above mentioned a priori expectations over state trajectories, which are
used to compute the non-adaptive Kalman filter, allow to compare the
stochastic results of the LQG controller between different sets of parameters.
In total, there are 5 parameters we aim to optimize:

  * •

the (terminal) velocity cost weight $\omega_{v}$,

  * •

the (terminal) force cost weight $\omega_{f}$,

  * •

the effort cost weight $\omega_{r}$,

  * •

the signal-dependent control noise level $\sigma_{u}$, and

  * •

the observation noise level $\sigma_{s}$.

The first three parameters correspond to those from Section 7.1. In
particular, setting all other parameters as well as the initial variance to
zero and using exact initial state estimates $\hat{x}_{0}=\bar{x}_{0}=x_{0}$
and $\hat{\Sigma}_{0}=\Sigma_{0}=0$, the stochastic optimal control problem
(30) equals the deterministic optimal control problem (11), i.e., the LQR is a
special case of the LQG. However, in contrast to the LQR model from Section 7,
distance, velocity, and force costs are only applied at the final time step
$N$. The effect of the individual cost weight parameters is thus much smaller.

In principle, the same fixed noise level parameters $\sigma_{u}$ and
$\sigma_{s}$ could be used for all participants and task conditions. However,
choosing them from literature is difficult since their effects strongly depend
on the assumptions and system dynamics of the considered model.

Figure 13.  Typical LQG trajectories with target shown as the orange box, as
well as selected entries of the corresponding feedback gain matrices, and
control time series, with noise levels $\sigma_{u}=0.2$, $\sigma_{s}=0.5$.
Top: Effect of velocity cost weight $\omega_{v}$ with fixed cost weights
$\omega_{f}=0$, $\omega_{r}=~{}$1e-7 (red dashed: $\omega_{v}=0$, green dash-
dotted: $\omega_{v}=1$, blue dotted: $\omega_{v}=20$). Middle: Effect of force
cost weight $\omega_{f}$ with fixed cost weights $\omega_{v}=0$,
$\omega_{r}=~{}$1e-7 (red dashed: $\omega_{f}=0$, green dash-dotted:
$\omega_{f}=1$, blue dotted: $\omega_{f}=100$). Bottom: Effect of effort cost
weight $\omega_{r}$ with fixed cost weights $\omega_{v}=\omega_{f}=2$ (red
dashed: $\omega_{v}=~{}$5e-8, green dash-dotted: $\omega_{v}=~{}$5e-7, blue
dotted: $\omega_{v}=~{}$5e-6).

Fully described in caption and text.

As shown in Figure 13 (top plots), an increase in the velocity cost weight
$\omega_{v}$ results in a higher peak velocity, which is attained earlier.
This can be explained by the corresponding decrease of the end-point velocity,
which effectively reduces the total costs. This is also visible from the
entries of the optimal feedback gain matrices $L_{n}$, which are depicted in
Figure 13 as well. The magnitude of the position component141414We recall that
to compute the resulting control $u_{n}$, each entry of the matrix $L_{n}$ is
multiplied by the corresponding entry of the state $x_{n}$, see (32). In our
case, $L_{n}$ is a $1\times 5$ matrix, which allows for easy interpretation of
the matrix entries. For example, the first entry of $L_{n}$ is multiplied by
the first entry of $x_{n}$, which is the position, and hence referred to as
the _position component_ of $L_{n}$. of $L_{n}$ determines the importance of
correcting the error between end-effector and target at time step $n$, whereas
the magnitude of the velocity component of $L_{n}$ determines the importance
of adjusting the end-effector velocity. For the LQG model, an increase in the
velocity cost weight $\omega_{v}$ results in a later peak of the velocity
component, i.e., the controller would correct for a large velocity that occurs
very shortly before the end of the movement, in order to reduce the terminal
velocity costs. In the position component of the feedback gain matrices, a
second, earlier peak occurs for moderate values ($\omega_{v}=1$, green dash-
dotted line). This suggests that it is particularly important to the
controller to eliminate deviations from the target both at the end of the
surge phase (peak at 0.212s for the considered case) and at the end of the
movement (peak at 0.68s). For large velocity cost weights ($\omega_{v}=20$,
blue dotted line), the earlier peak becomes more dominant and remains the only
peak, i.e., end-point errors that occur shortly before the end of the movement
are not corrected anymore, as this would result in a larger end-point velocity
and thus a higher total expected cost. In addition, the (relative) length of
the correction phase increases as the terminal velocity cost weight increases.
It is important to note that these effects only occur under both observation
noise and signal-dependent control noise, i.e., for $\sigma_{u}$,
$\sigma_{s}>0$.

Terminal force costs only, i.e., without terminal velocity costs, cannot
account for the typically observed corrective movements following the surge
(see Figure 13, middle plots). In particular, there is no incentive to reduce
the velocity towards the end of the movement. Instead, the acceleration is
reduced, which is visible from the constant velocity at the end of the
movement for large $\omega_{f}$ (green dash-dotted and blue dotted lines).
Similarly to the velocity cost weight $\omega_{v}$, the force cost weight
$\omega_{f}$ also affects the time period at which deviations from the desired
terminal position and velocity are corrected. As $\omega_{f}$ increases, this
period shifts forward, i.e., late-occurring deviations are not corrected
anymore.

For positive velocity and force cost weights $\omega_{v}$ and $\omega_{f}$, an
increase of the effort cost weight $\omega_{r}$ leads to a trajectory that
reaches the target later, with lower peak velocity (see Figure 13, bottom
plots). This is intuitive, since higher effort costs reduce the magnitude of
the optimal control signals, resulting in lower accelerations and velocities.
The terminal velocity is close to zero for any moderate $\omega_{r}$.

Note that the shape of the optimal control sequences differs considerably from
the deterministic LQR model (see Figure 11). In the LQG model, the control
usually attains its maximum at the beginning, then linearly decreases towards
its minimum, and increases again towards zero. This is mainly due to the
velocity and the acceleration being penalized only at the final time step $N$,
which allows to reach a higher peak velocity and acceleration (achieved
through larger control signals at the beginning of the movement; also note the
large (non-normalized) positional feedback gain values in bottom left plot).
Under the assumption of signal-dependent control and constant observation
noise, the control is very close to zero during the correction phase (which,
for ID $\geq 4$, makes up a considerable part of the movement). In Figure 13
(bottom right plot), this is shown for moderate control and observation noise
levels $\sigma_{u}=0.2$, $\sigma_{s}=0.5$.

If the control noise $\sigma_{u}$ is set large enough, similar effects can be
observed without any observation noise, i.e., $\sigma_{s}=0$ (not shown). In
contrast, large observation noise levels $\sigma_{s}$ cannot compensate a
missing control noise. This is intuitively plausible – since the controller is
aware that there is no control noise, the deterministic system states
resulting from the closed-loop system become perfectly predictable through the
internal (correct) forward model, i.e., there is no need to rely on noisy
observations at all.

Figure 14.  Top: Effect of the signal-dependent noise level $\sigma_{u}$ on
the variance time series of the LQG position and velocity profiles (red
dashed: $\sigma_{u}=3$, green dash-dotted: $\sigma_{u}=4$, blue dotted:
$\sigma_{u}=5$) and on the time until the target is reached; all simulations
with observation noise level $\sigma_{s}=0.5$. Bottom: Effect of the
observation noise level $\sigma_{s}$ on the variance time series of the LQG
position and velocity profiles (red dashed: $\sigma_{s}=0.02$, green dash-
dotted: $\sigma_{s}=0.2$, blue dotted: $\sigma_{s}=2$) and on the time until
the target is reached; all simulations with signal-dependent noise level
$\sigma_{u}=4$.

Fully described in caption and text.

The effects of the control noise level $\sigma_{u}$ (with $\sigma_{s}=0.5$)
and the observation noise level $\sigma_{s}$ (with $\sigma_{u}=4$) are shown
in Figure 14. In the left and in the middle plots, the variance of the
resulting LQG position and velocity profiles is plotted against time. The
black solid lines show a representative variance profile from the Pointing
Dynamics Dataset. Visually, the combination of signal-dependent control noise
and constant observation noise can explain the characteristic variance
profiles relatively well. Under signal-dependent noise, the application of
smaller controls in the second half of the movement (during the correction
phase, which follows the surge) results in a decreasing movement variability
towards the end of the movement. This is in accordance with the two-phase
positional variance profiles that are typically observed in aimed movements
(Gutman and Gottlieb, 1992; Lai et al., 2005; Gori and Rioul, 2020). From an
information-theoretic perspective, the decaying rate of these profiles can be
explained by a user-specific channel capacity (Gori and Rioul, 2020). In the
considered LQG model, the idea of user-specific control and observation noise
levels affecting the (expected) end-effector variability provides a different
interpretation. However, very large and physically implausible noise levels
($\sigma_{u}\approx 1-4$) are required to account for the substantial variance
observed in the position and velocity profiles of mouse pointing movements.
The velocity variance profiles, which are typically bimodal (Gutman and
Gottlieb, 1992; Darling et al., 1988), can also be replicated by the noise
model of the LQG, albeit the second peak is usually less pronounced in the
simulation.

The effect of the signal-dependent noise level $\sigma_{u}$ on the (average)
time until the target is reached is depicted in the top right plot of Figure