# Towards a statistical mechanics of consciousness: maximization of number of
connections is associated with conscious awareness

R. Guevara Erra1, D. M. Mateos2, R. Wennberg3, J.L. Perez Velazquez${}^{2}*$
1 Laboratoire Psychologie de la Perception, CNRS and Université Paris
Descartes, Sorbonne Paris Cité, Paris, France.
2 Neuroscience and Mental Health Programme, Division of Neurology, Hospital
for Sick Children.
Institute of Medical Science and Department of Paediatrics, University of
Toronto, Toronto, Canada.
3 Krembil Neuroscience Centre, Toronto Western Hospital, University of
Toronto, Toronto, Canada.
*<EMAIL_ADDRESS>
###### Abstract

It has been said that complexity lies between order and disorder. In the case
of brain activity, and physiology in general, complexity issues are being
considered with increased emphasis. We sought to identify features of brain
organization that are optimal for sensory processing, and that may guide the
emergence of cognition and consciousness, by analysing neurophysiological
recordings in conscious and unconscious states. We find a surprisingly simple
result: normal wakeful states are characterised by the greatest number of
possible configurations of interactions between brain networks, representing
highest entropy values. Therefore, the information content is larger in the
network associated to conscious states, suggesting that consciousness could be
the result of an optimization of information processing. These findings
encapsulate three main current theories of cognition, as discussed in the
text, and more specifically the conceptualization of consciousness in terms of
brain complexity. We hope our study represents the preliminary attempt at
finding organising principles of brain function that will help to guide in a
more formal sense inquiry into how consciousness arises from the organization
of matter

##  1 Introduction

How consciousness arises from the organization of matter is a subject of
debate that spans several disciplines, from philosophy to physics. Multitude
of studies focus on the investigation of patterns of synchrony in brain
activity based on magnitudes of a variety of synchrony indices, and thus the
search for organising principles of brain function is today more crucial than
ever. We sought to identify global features of brain organization that are
optimal for sensory processing and that may guide the emergence of conscious
awareness. Our results provide a (very simple) answer to the question of what
the magnitudes of synchrony indices represent in terms of the structure of
brain activity. Neurophysiological recordings of brain activity demonstrate
fluctuating patterns of cellular interactions, variability that allows for a
wide range of states, or configurations of connections of distributed networks
exchanging information, and support the flexibility needed to process sensory
inputs. Recent years have seen a surge in the study of fluctuations in brain
coordinated activity, studies that have raised conceptual frameworks such as
that of metastable dynamics [1] and that have motivated interest in the
practical application of assessments of nervous system variability for
clinical purposes [2, 3]. A prominent question is how to describe the
organizing principles of this collective activity, which allow features
associated with consciousness to emerge. What is the optimal brain
organization that allows it to adequately process sensory stimuli and enable
the organism to adapt to its environment?

Previous studies have revealed values of different indicators of brain
coordinated activity, such as synchronization, associated with healthy and
pathologic states by comparison of baseline values and those in, for instance,
unconscious states like coma and epileptic seizures [4, 5, 6]. These
observations prompt the question of what physiological organization underlies
the specific values of the synchrony indices found in normal alert states and
other conditions; in other words, what is special about these values found in
conscious states? We believe that we have provided an answer to this question
in our work. We propose that there is a certain general organization of brain
cell ensembles that will be optimal for conscious awareness, that is, for the
processing of sensory inputs. As an extension of previous work [7] where it
was proposed that a general organising principle of natural phenomena is the
tendency toward maximal – more probable – distribution of energy/matter, we
venture that the brain organization optimal for conscious awareness will be a
manifestation of the tendency towards a widespread distribution of energy (or,
equivalently, maximal information exchange). Whereas we do not deal with
energy or information in our work, we instead focus on the number of
(micro)states, or combinations of connected signals derived from specific
types of neurophysiologic recordings. We use the term “information” in the
intuitive sense that normally permeates neuroscience: cell ensembles that are
functionally connected to process/exchange information; furthermore, the
equivalence between information exchange and energy transactions has been the
subject of several studies, more specifically in [7] (see also [8, 9]). The
question then becomes: how do we capture the nature of these organizations of
cell interactions?

We have followed the classic approach in physics when it comes to
understanding collective behaviours of systems composed of a myriad of units:
the assessment of the number of possible configurations, or microstates, that
the system can adopt. In our study we focus on the collective level of
description and assume that coordinated patterns of brain activity evolve due
to interactions of mesoscopic areas ( [10, 11]). Thus we use several types of
brain recordings in conscious and unconscious states, evaluating the number of
“connections” between these areas and the associated entropy and complexity.
We present evidence that conscious states result from higher entropy and
complexity in the number of configurations of pairwise connections. The number
of pairwise channel combinations is near the maximum of all possible
configurations when the individual is processing sensory inputs in a normal
manner (e.g. with open eyes). Our interpretation is that a greater number of
configurations of interactions allows the brain to optimally process sensory
information, fostering the necessary variability in brain activity needed to
integrate and segregate sensorimotor patterns associated with conscious
awareness.

##  2 Materials and Methods

### Electrophysiological recordings

Recordings were analysed from 9 subjects, using magnetoencephalography (MEG),
scalp electroencephalography (EEG) or intracranial EEG (iEEG). Three epilepsy
patients were studied with MEG; 1 epilepsy patient was studied with iEEG; 3
epilepsy patients were studied with simultaneous iEEG and scalp EEG; and 2
nonepileptic subjects were studied with scalp EEG. For the study of seizures
versus alert states, the 3 subjects with MEG recordings and the one with iEEG
were used. Details of the patients’ epilepsies and seizure types have been
presented in previous studies (MEG patients in [5]; iEEG patients in [12]).
For the study of sleep versus alert states, the 3 patients with combined iEEG
and scalp EEG have been described previously (patients 1, 3, 4 in [13]); the 2
subjects studied with scalp EEG alone had been investigated because of a
suspected history of epilepsy, but both were ultimately diagnosed with
syncope, with no evidence of epilepsy found during prolonged EEG monitoring.
In brief, the MEG seizure recordings were obtained in one patient with primary
generalized absence epilepsy, in one patient with symptomatic generalized
epilepsy, and in one patient with frontal lobe epilepsy. The iEEG seizure
recordings were obtained from a patient with medically refractory temporal
lobe epilepsy as part of the patient’s routine clinical pre-surgical
investigation. MEG recordings were obtained using a whole head CTF MEG system
(Port Coquitlam, BC, Canada) with sensors covering the entire cerebral cortex,
whereas iEEG electrodes were positioned in various locations including, in the
temporal lobe epilepsy patient, the amygdala and hippocampal structures of
both temporal lobes. EEG recordings were obtained using an XLTEK EEG system
(Oakville, ON, Canada). The details of the acquisitions varied from patient to
patient (e.g., acquisition rate varied from 200 to 625 Hz) and were taken into
consideration for the data analyses. The duration of the recordings varied as
well: for the seizure study, the MEG sample epochs were of 2 minutes duration
each, with total recording times of 30-40 minutes; the iEEG patient sample was
of 55 minutes duration. The sleep study data segments were each 2-4 minutes in
duration, selected from continuous 24-hour recordings.

### Data analysis

The only pre-processed data were those of the scalp EEG recordings. These were
processed using a Laplacian derivation [11], to avoid the potential effects of
the common reference electrode on synchronization [14] using the DSC algorithm
[15]. Initially a phase synchrony index was calculated from all possible
pairwise signal combinations, for which we use the standard procedure of
estimating phase differences between two signals from the instantaneous phases
extracted using the analytic signal concept via the Hilbert transform. To
compute the synchrony index, several central frequencies, as specified in the
text and figure legends, were chosen with a bandpass filter of 2 Hz on either
side, hence, for one value of the central frequency $f$, the bandpass is $f\pm
2$ Hz. The central frequencies were chosen according to the relevant
behavioural states and some analytical limitations– thus for the sleep studies
we choose 4 Hz (not lower because the extraction of the instantaneous phase
was not optimal for central frequencies lower than 4). To see whether similar
results were obtained with different frequencies, we chose others (see
figures) provided there was power at those values. The phase synchrony index
($R$) was calculated using a 1-second running window and was obtained from the
phase differences using the mean phase coherence statistic which is a measure
of phase locking and is defined as $R=\lvert\langle
e^{i\Delta\theta}\rangle\lvert$ where $\Delta\theta$ is the phase difference
between two signals. This analytical procedure has been described in great
detail elsewhere [4, 5, 14]. The calculation of the index R was done for all
possible signal pairs. The mean value of the $R$ index thus obtained was then
estimated for the desired time length. For instance, for the sleep recordings,
the time period was the whole episode, which was, as noted above, between 2
and 4 minutes. For the seizure recordings, the periods to obtain the mean
synchrony varied depending on the behavioural condition, for instance, in
figure 1C the whole ictal event (labelled ‘Sz’) and the initial 10-second
portion of it (’10 sec Sz’) were taken for the reason explained in the text
and figure legend. The calculation of the number of connected signals and the
entropy associated is described in the Results section. We note here that to
assess entropy we assume that the different pairwise configurations are
equiprobable, thus the entropy is reduced to the logarithm of the number of
states, $S=lnC$ (see notation in the main text). However, the estimation of
$C$ (the combinations of connections between diverse signals), is not feasible
due to the large number of sensors; for example, for 35 sensors, the total
possible number of pairwise connections is $[1442]=10296$, then if we find in
the experiment that, say, $2000$ pairs are connected, the computation of
$[102962000]$ has too large numbers for numerical manipulations, as they
cannot be represented as conventional floating point values in, for instance,
MATLAB. To overcome this difficulty, we used the well-known Stirling
approximation for large $n:ln(n!)=n~{}ln(n)\textendash n$. The Stirling
approximation is frequently used in statistical mechanics to simplify entropy-
related computations. Using this approximation, and after some basic algebra,
the equation for entropy reads, $S=N~{}ln(N/N-p)-p~{}ln(p/N-p)$, where $N$ is
the total number of possible pairs of channels and p the number of connected
pairs of signals in each experiment (see Results for details and notation).
Because this equation is derived from the Shannon entropy, it indicates the
information content of the system as well (37).

In addition to entropy, we used another measure of complexity, the Lempel-Ziv
(L-Z) complexity, based in the Kolmogorov deterministic complexity [16]. This
complexity measure the amount of non-redundant information in a string by
estimating the minimal size of the ”vocabulary” necessary to describe the
string [17]. Strings with high L-Z complexity require a large number of
different patterns (“words”) to be reproduced, while strings with low
complexity can be largely compressed with a few patterns employed to eliminate
redundancy with no loss of information. For this purpose, values of the matrix
$B$ (defined in Results) were placed in a one-dimensional vector and its L-Z
complexity determined.

##  3 Result

Guided by proposals that consciousness requires medium values of certain
features of cell assemblies, e.g. not too high or low synchrony or
correlations [18, 19, 3], or halfway between order and disorder [20], we chose
to quantify the number of possible configurations the brain can adopt in
different behavioural conditions. Our basic approach consists in the
estimation of the number of possible pairwise connections between recorded
brain signals. Signals included MEG, iEEG and scalp EEG recordings; details
can be found in Methods. We are limited to pairwise combinations of the
signals because of the manner in which phase synchrony is computed – as phase
differences between two signals – and we use phase synchronization as the
means to determine “connectivity” between the two signals. Once the number of
“connected” signals is known, we estimate the entropy of those pairwise
combinations. The results obtained with recordings acquired during conscious
states are compared with those acquired during unconscious states, which
included sleep (all stages) and epileptic seizures. To determine connectivity,
we use an accepted approach of computing a phase synchronization index
(details in Methods).

It must be noted that, while many studies use the words synchrony and
connectivity as synonymous, in reality phase synchrony analysis reveals only a
correlation between the phases of the oscillations between two signals, and
not a real connectivity which depends on several other factors (this matter
has been covered in detail in [21, 22]. Nevertheless, due to the unfeasibility
of an accurate, realistic estimation of connectivity which would necessitate
individual cell recordings from entire cell ensembles as well as structural
connectivity details, we use an accepted version under the assumption that the
phase relations may represent, at least, some aspect of a functional
connectivity. Hence, in order to evaluate interactions (“connections”), we
take each sensor/channel as one “unit”, and define a pair of signals as
“connected” if the phase synchrony index is larger than a threshold. The
threshold is determined for each individual, and is the average synchrony
index in the most normal alert state, the ‘awake eyes-open’ condition, when
the individual is fully alert and processing the sensorium in a regular
fashion. Because our data include the three recording methodologies
aforementioned, we have the opportunity to assess the reproducibility of the
results in various types of recordings. While we work at the signal level we
will make the reasonable assumption that the MEG and scalp EEG sensors record
cortical activity underlying those sensors and thus throughout the text we
will use the terms brain signals or brain areas/networks as synonymous. The
iEEG, obviously, records signals at the source level. Note that we are not
interested in the specific pattern of connectivity among brain sources/areas,
but rather in the global states. These points are further discussed in the
Discussion section.

Phase synchronization for specific frequencies (details in Methods) is
calculated for each pair of channels and a “connectivity” matrix $S$ is
obtained, whose entries are the average values of the synchrony index during a
certain time period for each pairwise configuration. From this matrix, a
Boolean connectivity matrix $B$ is calculated, with 0 entry if the
corresponding synchrony index is lower than a threshold, and 1 if higher. We
define two channels as “connected” if the corresponding entry in matrix $B$ is
1. Then we use the combinations of connected channels as a ‘complexity’
measure. The total number of possible pairs of channels given a specific
channel montage is given by the binomial coefficient $N=Nc!/2!~{}(Nc-2)!$
where $Nc$ is the total number of channels in the recording montage, normally
$144-146$ in the case of MEG sensors, and between 19 and 35 with iEEG and
scalp EEG. For instance, in our MEG recordings we have$Nc=144$, thus
$N=10,296$possible pairs of connected sensors are obtained. For each subject
we calculate p, the number of connected pairs of signals in the different
behavioural states, using the aforementioned threshold of the synchrony index
(which varies for each subject), and estimate $C$, the number of possible
combinations of those p pairs, using the binomial coefficient again:
$C=N!/p!(N-p)!$ where $N$ is the aforesaid value. In sum, all these
calculations represent the relatively simple combinatorial problem we are
trying to solve: given a maximum total of $N$ pairs of connected signals, in
how many ways can our experimental observation of p connected pairs (that is,
the number of 1’s in matrix $B$) be arranged. We then compute the entropy and
Lempel-Ziv complexity associated with those p values.

Figure 1 depicts the entropy ($S$), which, assuming equiprobable states, is
the logarithm of the number of states, $S=ln~{}C$ (see Methods for the
estimation of entropy using very large $C$ values) in three epileptic
patients. We note that this equation using the natural logarithm allows for a
calculation of both the Gibbs and the Shannon entropy, if needed, which differ
from $S$ by a constant multiplicative factor k (the Boltzmann constant) and
$1/(ln2)$ respectively. In reality the entropy estimation does not provide any
further information, as the main, crucial result is the number of
configurations, $C$. However, we have done it since it is a standard manner to
quantify the “complexity” of the number of microstates. The entropy data
points are graphed on the curve that represents the entropy of all points in
the binomial distribution, where the maximum number of configurations (that
is, maximum entropy) occurs in the middle of the graph. Note in that figure
that during conscious states, when patients are not having generalised
seizures with loss of awareness, the entropy is close to the maximum, whereas
entropy is lower (more distant from the top) for the seizure states. The
values during the seizures fall on the right-hand side of the graph because,
due to the higher synchrony during ictal (seizure) events such that the number
of coupled channels (note that the x-axis is the number of connected signals)
is larger than during interictal (between seizure) activity, there are fewer
pairwise configurations and thus lower entropy. This phenomenon seems
associated with the level of consciousness since when the seizures are not
generalised (Figure 1C and 1D), and the patients remain responsive and
conscious, the entropy values are similar to those of interictal (baseline)
activity. Where in the curve the data points are located depends of course on
the synchrony index. Because seizures had higher synchrony than interictal
periods (“baseline”), the number of coupled signals is greater and the number
of combinations is lower.

We observe similar trends in the case of sleep. Figure 2 depicts some
examples. Note how during wakefulness the entropy is closer to the maximum of
the curve, whereas the deeper the sleep stage, the more distant from the
maximum the values are. The deepest sleep stage, slow wave 3-4 (‘sws3-4’), has
consistently the lowest entropy. Interestingly, the entropy during REM sleep
is very close, in most cases, to the normal, alert state. This result may not
be as surprising as it sounds if we consider the mental activity during REM
episodes that are normally associated with dreams. It is worth noting too that
in recordings taken when the subjects had their eyes closed, the entropy is
much lower than during the eyes open condition, and sometimes it is as low as
during slow-wave 3-4 sleep. The results for two central frequencies (4 and 8
Hz) are shown in Figure 2A and 2C, to demonstrate that it is not always the
case of high synchrony having lower entropy; sometimes it is lower synchrony
(e.g., results at 8 Hz) that resulted in fewer channel combinations and thus
lower entropy as compared to fully alert states. In Figures 1 and 2 we have
shown results with iEEG and MEG (Figure 1) or iEEG and scalp EEG (Figure 2)
recordings to demonstrate that the same qualitative results are obtained with
these three recordings techniques. Thus these results do not depend on
recording methodology. Shown in Figure 3 is an example of the time course of
the entropy before and during an ictal event. Hence, we have demonstrated that
the specific values of the synchrony index R in fully alert states represent
the largest number of combinations (microstates, see Discussion) of pairwise
signal configurations. This method solves the potential problem of the
interpretation of the different R values obtained with various recordings
techniques. For example, in our experiments, average values in baseline
conditions are $0.248\pm 0.2$ for MEG, $0.428\pm 0.04$for iEEG, and $0.46\pm
0.05$ for scalp EEG, nevertheless, in our study the number of combinations are
“normalised” to the number of recording sensors and show a final result (the
entropy) that is independent of the structure and synchrony magnitudes of the
recording methodology. The main idea derived from these results is represented
in Fig. 4.

To further explore whether the decrease in entropy has a parallel with a
decrease in other forms of complexity, the Lempel-Ziv complexity of the number
of configurations was assessed. Tables 1 and 2 illustrate the Lempel-Ziv
complexity estimated for the $B$ matrices, where a complementary result can be
seen: unconscious states have lower values of complexity. A decrease of
complexity in the raw electrophysiological signals was obtained too (data not
shown), hence this may be a phenomenon observable at different levels of
description.

##  4 Discussion

Our attempts at seeking features of brain organization that allow for adequate
processing of sensory stimuli have provided evidence that a greater number of
possible configurations of interactions between brain networks is associated
with alert states, representing high entropy associated with the number of
those combinations, whereas lower entropy (and thus fewer combinations of
connections) is characteristic of either unconscious states or altered states
of alertness (eyes closed). This observation reflects a relatively simple
general organising principle at this collective level of description, which
results in the emergence of properties associated with consciousness.

With the advent of the ‘Big Data’ era and the related torrent of empirical
observations, the search for organising principles that result in the
emergence of biological phenomena seems more crucial than ever. We tried to
uncover relatively simple laws that capture the bounds in the global
organization of a biological system that enable it to become adaptable (i.e.,
responsive) to an environment, or, in neuroscientific terms, the features of
optimal brain organization (in terms of connections) that allow brains to
adequately process sensory stimuli. We focused on the global states and did
not investigate specific patterns of connectivity between brain areas as a
variety of other studies have assessed these inter-regional interactions in
conscious and unconscious states [23, 24]. The fact that our results are
similar, independent of recording methodology, demonstrates the robustness of
the analysis. On that note, we remark that while we have used the term
‘connectivity’, in reality the analysis reveals only correlation between
phases of oscillation, as already discussed in Results.

The present study can be considered an extension of previous work where it was
proposed that a general organising principle of natural phenomena is the
tendency toward a maximal, or more probable, distribution of energy [7], which
can be encapsulated by the notion of the maximization of information transfer
[25]. As well, the notions of information and energy exchange are conceptually
related: “the common currency paying for all biological information is energy
flow” [26, 9]. In the final analysis, information exchange implies energy
exchange, hence we interpret information exchange as energy redistribution as
proposed in [7], even though our study is focused not on energy considerations
but on the number of states. Other studies have assessed the importance of
brain synchronization to optimally transfer information [27]. We interpret our
observation that the number of pairwise channel combinations – that we take as
interactions/connections between brain networks – occurs near the maximum of
possible configurations in periods with normal alertness, as that greater
number of configurations of interactions represents the most probable
distribution of energy/information resulting in conscious awareness. The
configuration entropy we calculate measures the information content of the
functional network, and has been used in other works for the purpose of
quantifying information [28]. One somewhat surprising result is the low
entropy during the eyes-close condition. One could argue that having the eyes
closed does not change much the conscious state, and yet we observe reduced
entropy as compared with eyes-open condition. Hence our observations may
indicate not only states of consciousness but also optimality of sensory
processing – considering the great importance of visual processing in humans,
interrupting visual inputs should result in considerable changes in the
dynamics; for instance, an apparent alteration of brain dynamics upon eye
closure is the appearance of alpha waves most clearly in parieto-occipital
regions. Interesting too is the relatively high entropy associated with REM
episodes, perhaps a reflection of the partial awareness during dreams. Perhaps
the main difference between dream and awake states is psychological, but they
share similar brain dynamics.

It has been proposed that aspects of awareness emerge when certain levels of
complexity are reached [29]. It is then possible that the organization
(complexity) needed for consciousness to arise requires the maximum number of
configurations that allow for a greater variety of interactions between cell
assemblies because this structure leads to optimal segregation and integration
of information, two fundamental aspects of brain information processing [20] .
Our results help cast the study of consciousness and cognition into more of a
physics framework that may provide insight into simple principles guiding the
emergence of conscious awareness, and perhaps the proposed thermodynamics for
a network of connected neurons [18] can be extended to explain cognition. Some
classical studies [8] have already characterized biological order in terms of
functions of the number of states of a system. It is tempting to link our
observations with the typical chemical equilibrium that, despite being
composed of a myriad of microstates, when viewed at the macroscopic level
produces some useful laws, like the law of mass action describing chemical
balance, and thus the perspective we develop here may help guide research to
uncover organising principles in the neurosciences. In physics, microstates
that yield the same macrostate form an ensemble. A system tends to approach
the most probable state, maximising entropy under present constraints, and the
resulting macrostate will be represented by the maximum number of microstates.
Hence, the macrostate with higher entropy (see scheme in Fig. 4) we have
defined, composed of many microstates (the possible combinations of
connections between diverse networks, our $C$ variable defined in Results),
can be thought of as an ensemble characterised by the largest number of
configurations. Here we define an ensemble of microstates as all possible
configurations of connectivity leading to the same macrostate (having the same
number of connected pair of signals, $p$). The entropy of this macrostate,
given by the logarithm of the number of combinations (our $C$), is the number
of microstates that are compatible with the given macrostate (assuming all
microstates have same statistical weight). In neurophysiological terms, each
microstate represents a different connectivity pattern and thus is associated
with, in principle, different behaviours or cognitive processes. The
macrostate that we find associated with wakeful normal states (eyes open) is
the most probable because it has the largest entropy (largest number of
combinations of connections). While many elementary cellular microscopic
processes are far from equilibrium (e.g., ionic gradients), at the macroscopic
level the system tends towards equilibrium, as most natural phenomena remain
in near-equilibrium conditions [30]. At the same time, the ensemble of
microstates associated with normal sensory processing features the most varied
configurations and therefore offers the variability needed to optimally
process sensory inputs. In this sense, our results support current views on
the metastability of brain states [1] in that the states should not be too
stable for efficient information processing, hence the larger the number of
possible interactions, the more variability is possible.

Equally, the results are consistent with the global workspace theory [31] in
that the most widespread distribution of information leads to conscious
awareness. Furthermore, these observations relate as well to the information
integrated theory [32], in that consciousness increases in proportion to the
system’s repertoire of states, thus the more combinations possible, the more
states available, and here we can define states as configurations of
interactions. Additionally, the results support computational/theoretical
studies showing that patterns of organised activity arise from the
maximization of fluctuations in synchrony [33], and that transitions between
conscious states are achieved by just varying the probability of connections
in neural nets [34]. In general, our observations highlight the fundamental
importance of fluctuations in neuronal activity as the source of healthy brain
dynamics. More specifically, our results have a precise parallel with the work
of Hudetz et al. (2014), where they graph a dispersion index versus an
activation level (their figure 7B) and propose that consciousness resides at
the top of the curve, and anaesthetic states and seizures to lower and higher
activation levels respectively, as we show in Figures 1 and 2. Their
‘activation level’ could correspond to the number of signals that take part in
the combinations (our x axes in the graphs), and their ‘dispersion index’ to
our number of combinations (the y axes). Other studies have proposed as well
that consciousness requires medium values of certain features of cell
assemblies [32, 35].

In sum, along with others, we consider cognition/consciousness not a static
property but a dynamic process with constant flux of energy, or information
exchange [25]. Even though we have talked above about a macrostate, this
should not be taken as a fixed state, rather it contains dynamic processes
represented by the microstates, the re-arrangements of connections among brain
cell ensembles. The emergent features of cognitive phenomena that can be
termed “conscious” arise once an efficient web of connections endowed with
certain complexity appears. The fact that values of phase synchrony during
fully alert states gave us the largest entropy of the number of pairwise
signal combinations explains in part the neurophysiological organization
underlying these specific values of the quantification of brain synchrony.
Studies at this level of description may help to understand how consciousness
arises from organization of matter. In our view, consciousness can be
considered as an emergent property of the organization of the (embodied)
nervous system, especially a consequence of the most probable distribution
that maximizes information content of brain functional networks.

On a technical note, we remark that in our analysis we had to choose a
threshold to define “connectivity”, and therefore a suitable baseline had to
be chosen. This procedure is obviously biased in that a signal corresponding
to a baseline has to be ascribed as the one providing the threshold, and we
chose the signals corresponding to the (psychological) state that is most
suited for the purposes of adaptability: fully alert and receiving all sensory
inputs (awake with open eyes). By this choice, and according to our methods,
that signal is already ascribed as having maximum entropy. We tried another
less prejudiced method, using surrogates of the original signals, and then
computing the average synchrony index among the surrogate population (10
phase-randomised surrogates per original channel/signal) that is the threshold
to define connectivity. It turns out that the value of the magnitude of
synchrony of the surrogates is close to the one for the aforementioned
baseline chosen, so the results do not vary; nevertheless this new method
still assigns the largest entropy to the random signals (surrogates), so there
is still the assumption that the average synchrony of the stochastic signals
is a good approximation to define connections among brain networks. Our
current purpose is to completely avoid using thresholds of synchrony indices,
for which we are presently working on a scheme that assigns connectivity in
the time domain.

##  5 Conclusions and conjectures on the structure of brain-behaviour-
environment

It is tempting to speculate, based on these results and the conclusions of a
previous study [7] that there could be a universal logic ruling the evolution
of natural phenomena — biological and nonbiological— and the nervous system in
particular: patterns emerge from a central theme captured by maximising
information exchange. Because, in the final analysis, all exchange of
information implies exchange of energy, natural phenomena tend towards the
most probable distribution of energy, and thus the interactions among system
constituents tends to be maximized.

Because the brain functions to maintain a predictive model of the environment
(the reason the brain evolved is to model the environment, after all), then
perhaps the brain’s global configuration has to “copy” what is out there: and
out there energy distributes in all possible microstates (second principle of
thermodynamics). Then to process such variability in nature, the nervous
system should have same structure, and the result is the ‘inverted U’ that has
appeared in our analysis and has been theoretically proposed in other
publications [20], the top of the curve representing more possible
combinations to handle information/energy exchanges. On the other hand, in the
extremes of this curve we find fewer microstates, thus these are not optimal
situations to process the many microstates in the environment. The key then is
not to reach the maximum number of units interacting (which would be all-to-
all connections and thus only one possible microstate), but rather the largest
possible number of configurations allowed by the constraints. In a similar
fashion, it has been argued that the brain needs to show criticality because
natural phenomena possess critical dynamics [36]. Then, perchance,
consciousness can be considered as an emergent property of the organization of
the embodied nervous system submerged in an environment, consequence of the
most probable distribution of energy (information exchange) in the brain. In
this regard, consciousness (like biochemistry) may represent thus an optimal
channel for accessing sources of (free) energy.

## References

  * [1] JA. S. Kelso.  Dynamic patterns: The self-organization of brain and behavior.  MIT press, 1997.
  * [2] D. D Garrett, G. R. Samanez-Larkin, S. WS. MacDonald, U. Lindenberger, A. R McIntosh, and C. L Grady.  Moment-to-moment brain signal variability: a next frontier in human brain mapping?  Neuroscience & Biobehavioral Reviews, 37(4):610–624, 2013.
  * [3] V. Nenadovic, J. L. Perez Velazquez, and J. S. Hutchison.  Phase synchronization in electroencephalographic recordings prognosticates outcome in paediatric coma.  PloS one, 9(4):e94942, 2014.
  * [4] F. Mormann, K. Lehnertz, P. David, and C. E. Elger.  Mean phase coherence as a measure for phase synchronization and its application to the eeg of epilepsy patients.  Physica D: Nonlinear Phenomena, 144(3):358–369, 2000.
  * [5] L.. Garcia Dominguez, R. A. Wennberg, W. Gaetz, D. Cheyne, O. C. Snead, and J. L. Perez Velazquez.  Enhanced synchrony in epileptiform activity? local versus distant phase synchronization in generalized seizures.  The Journal of neuroscience, 25(35):8077–8084, 2005.
  * [6] D. C. Shields, J. W. Leiphart, D. L. McArthur, P. M. Vespa, M. L. Quyen, J. Martinerie, and J. R. Soss.  Cortical synchrony changes detected by scalp electrode electroencephalograph as traumatic brain injury patients emerge from coma.  Surgical neurology, 67(4):354–359, 2007.
  * [7] J. L. Perez Velazquez.  Finding simplicity in complexity: general principles of biological and nonbiological organization.  Journal of biological physics, 35(3):209–221, 2009.
  * [8] H. J. Morowitz.  Some order-disorder considerations in living systems.  The bulletin of mathematical biophysics, 17(2):81–86, 1955.
  * [9] E. Smith.  Thermodynamics of natural selection i: Energy flow and the limits on organization.  Journal of theoretical biology, 252(2):185–197, 2008.
  * [10] J.J. Wright and D.T.J. Liley.  Dynamics of the brain at global and microscopic scales: Neural networks and the eeg.  Behavioral and Brain Sciences, 19(02):285–295, 1996.
  * [11] P. L. Nunez.  Toward a quantitative description of large-scale neocortical dynamic function and eeg.  Behavioral and Brain Sciences, 23(03):371–398, 2000.
  * [12] J. L. Perez Velazquez, Garcia Dominguez L., V. Nenadovic, and R. A. Wennberg.  Experimental observation of increased fluctuations in an order parameter before epochs of extended brain synchronization.  Journal of biological physics, 37(1):141–152, 2011.
  * [13] R. Wennberg.  Intracranial cortical localization of the human k-complex.  Clinical Neurophysiology, 121(8):1176–1186, 2010.
  * [14] R. Guevara, J. L. Pérez Velazquez, V. Nenadovic, R. Wennberg, G Senjanović, and L. G. Dominguez.  Phase synchronization measurements using electroencephalographic recordings.  Neuroinformatics, 3(4):301–313, 2005.
  * [15] J. Kayser and C. E. Tenke.  Principal components analysis of laplacian waveforms as a generic method for identifying erp generator patterns: I. evaluation with auditory oddball tasks.  Clinical neurophysiology, 117(2):348–368, 2006.
  * [16] A. N. Kolmogorov.  Three approaches to the quantitative definition of information.  Problems of information transmission, 1(1):1–7, 1965.
  * [17] A. Lempel and J. Ziv.  On the complexity of finite sequences.  Information Theory, IEEE Transactions on, 22(1):75–81, 1976.
  * [18] G. Tkacik, T. Mora, O. Marre, D. Amodei, II Berry, J. Michael, and W. Bialek.  Thermodynamics for a network of neurons: Signatures of criticality.  arXiv preprint arXiv:1407.5946, 2014.
  * [19] V Nenadovic, J. S. Hutchison, L. G. Dominguez, H. Otsubo, M. P. Gray, R Sharma, J. Belkas, and J. L. Perez Velazquez.  Fluctuations in cortical synchronization in pediatric traumatic brain injury.  Journal of neurotrauma, 25(6):615–627, 2008.
  * [20] G. Tononi, G. M. Edelman, and O. Sporns.  Complexity and coherency: integrating information in the brain.  Trends in cognitive sciences, 2(12):474–484, 1998.
  * [21] J. L. Perez Velazquez, R. G. Erra, R. Wennberg, and L. G. Dominguez.  Correlations of cellular activities in the nervous system: physiological and methodological considerations.  In Coordinated Activity in the Brain, pages 1–24. Springer, 2009\.
  * [22] J. L. Perez Velazquez and M. Frantseva.  The Brain-Behavior Continuum: The Subtle Transition Between Sanity and Insanity.  World Scientific, 2011.
  * [23] G Dumermuth and D Lehmann.  Eeg power and coherence during non-rem and rem phases in humans in all-night sleep analyses.  European neurology, 20(6):429–434, 1981.
  * [24] K. Wang, M. L. Steyn-Ross, D.A. Steyn-Ross, M. T. Wilson, and J. W. Sleigh.  Eeg slow-wave coherence changes in propofol-induced general anesthesia: experiment and theory.  Frontiers in systems neuroscience, 8, 2014.
  * [25] Hermann Haken.  Information and self-organization: A macroscopic approach to complex systems.  Springer Science & Business Media, 2006.
  * [26] H. J Morowitz.  Energy flow in biology; biological organization as a problem in thermal physics.  1968\.
  * [27] A Buehlmann and G Deco.  Optimal information transfer in the cortex through synchronization.  PLoS Comput Biol, 6(9):e1000934, 2010.
  * [28] W. Bialek, F. Rieke, R de Ruyter van Steveninck, and D Warland.  Spikes: exploring the neural code.  MIT. Roddey, JC, Girish, B., & Miller, JP (2000). Assessing the performance of neural encoding models in the presence of noise. Journal of Computational Neuroscience, 8(95):112, 1997.
  * [29] M. Gell-Mann.  Consciousness, reduction, and emergence.  Annals of the New York Academy of Sciences, 929(1):41–49, 2001\.
  * [30] I. Prigogine.  Introduction to thermodynamics of irreversible processes.  New York: Interscience, 1967, 3rd ed., 1, 1967.
  * [31] B. J. Baars.  A cognitive theory of consciousness.  Cambridge University Press, 1993.
  * [32] G. Tononi.  An information integration theory of consciousness.  BMC neuroscience, 5(1):42, 2004.
  * [33] V. Vuksanović and P. Hövel.  Dynamic changes in network synchrony reveal resting-state functional networks.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 25(2):023116, 2015.
  * [34] D. W. Zhou, D. D. Mowrey, P. Tang, and Y. Xu.  Percolation model of sensory transmission and loss of consciousness under general anesthesia.  Physical review letters, 115(10):108103, 2015.
  * [35] A. A Fingelkurts, A. A. Fingelkurts, S. Bagnato, C. Boccagni, and G. Galardi.  Do we need a theory-based assessment of consciousness in the field of disorders of consciousness?  Frontiers in human neuroscience, 8:402, 2014.
  * [36] D. R. Chialvo.  Emergent complex neural dynamics.  Nature physics, 6(10):744–750, 2010.

| State | L-Z complexity
---|---|---
Subject #1 | Baseline | 0.7
| Seizure | 0.2
Subject #2 | Baseline | 0.61
| Seizure | 0.39
Subject #3 | Baseline | 0.86
| Seizure | 0.71
Subject #4 | Baseline | 0.61
| Seizure | 0.61
Table 1: Values of Lempel-Ziv (L-Z) complexity derived from the string of connections (details in Methods) in conscious (baseline) and unconscious (seizure) states, in four patients. Note lower complexity during seizures in patients $\\#~{}1-3$; patient $\\#~{}4$ (Figure 1D) did not have fully generalised seizures. | State | iEEG L-Z complexity | scalp L-Z complexity
---|---|---|---
Subject #1 | Alert eyes open | 0.81 | 1.01
| SWS 2 | 0.42 | 1.0
| SWS 3-4 | 0.018 | 0.0
| REM | 0.73 | 0.96
Subject #2 | Alert eyes open | 0.94 | 1.06
| Alert eyes close | 0.94 | 0.85
| SWS 3-4 | 0.33 | 0.22
Subject #3 | Alert eyes open | N/A | 1.07
| SWS 2 | N/A | 1.07
| SWS 3-4 | N/A | 1.07
| REM | N/A | 0.97
Subject #4 | Alert eyes open | 0.88 | 1.0
| Alert eyes closed | 0.55 | 1.0
| SWS 2 | 0.57 | 9.8
| SWS 3 | 0.05 | 0.0
| REM | 0.6 | 1.06
Subject #5 | Alert eyes open | N/A | 0.97
| SWS 1 | N/A | 1.08
| SWS 2 | N/A | 1.04
| SWS 3-4 | N/A | 0.8
| REM | N/A | 1.12
Table 2: Values of Lempel-Ziv (L-Z) complexity derived from the string of
connections in different sleep stages. When both recordings were obtained from
a subject, both iEEG and scalp EEG data were analysed (subjects $\\#~{}3$ and
$\\#~{}5$ did not have iEEG recordings). The L-Z complexity is consistently
lower, regardless of recording methodology, in the deepest sleep stage (SWS
3-4). These results parallel those of entropy estimations shown in figure 2

Figure 1: Graphs representing the entropy of the number of pairwise
configurations of signals in epileptic patients during conscious (baseline)
and unconscious (generalised seizure) states. A, derived from MEG recordings
in a patient with primary generalised epilepsy, shows entropy associated with
a normal alert period (baseline, ‘Base’) and a generalised absence seizure
(‘Sz’), estimated from synchrony values at two central frequencies (defined in
Methods) of 5 and 12 Hz. The curve in this and other graphs here and in Figure
2 represents the possible entropy values of all possible numbers of pairwise
combinations, yielding an inverted U or, for very large numbers, a Gaussian.
Note that here as well as in all generalised seizures analysed, the entropy
values associated with alert, baseline conditions were closer to the maximum
(top of the curve) than those associated with the seizures. B, entropy values
of two seizures and their corresponding baseline (‘Base’) activity (computed
using a time period of 30-40 minutes before the ictus) in a patient with
secondary (symptomatic) generalised epilepsy (MEG recordings). C, derived from
iEEG recordings in a patient with temporal lobe epilepsy, shows the entropy
during the alert state (‘baseline’), during the initial 10 seconds of the
seizure when the patient was still responsive and alert (‘10 sec Sz’), and
during the rest of the seizure when it became generalised and the patient was
unresponsive (‘Sz’). Note that when the ictus has not yet generalised, the
entropy is similar to that of normal alertness. D, another example of a non-
generalised seizure in a patient with frontal lobe epilepsy (MEG recordings).
Figure 2:  Same graph types as in Figure 1, using sleep recordings. In each
subject, data samples were of 2-4 minutes duration during wakefulness with
eyes open (‘awake’) or eyes closed (‘Eyecl’), and sleep stages slow-wave 2
(‘Sws2’), slow-wave 3-4 (‘Sws3-4’) and rapid eye movement (‘REM’). A, results
derived from iEEG recordings in a subject investigated with bilateral frontal
and temporal electrodes and simultaneous scalp EEG. Entropy estimated from
synchrony values at two central frequencies of 4 and 8 Hz. As occurred in the
patient recordings shown in Figure 1, the baseline, alert state (in this case
labelled ‘awake’) is closer to the top of the curve, having greater entropy.
The deepest sleep stage, slow wave 3-4 (‘sws3-4’), has the lowest entropy. B,
same subject but using the scalp EEG recordings for the calculations, showing
similar trend (evaluated at central frequency of 4 Hz). C, results derived
from a different subject investigated with right frontal electrodes and
simultaneous scalp EEG, with synchrony evaluated at two central frequencies of
4 and 8 Hz using the iEEG signals. Note that the eyes closed (‘Eyecl’)
condition has lower entropy than that of the normal alert state with open
eyes. Depending on the frequency of analysis, the entropy in ‘Eyecl’ falls
toward the left or right side of the curve; e.g., at 8 Hz, because the
synchrony is higher (more coupled channels) due to alpha waves at 8-10 Hz and
the entropy is reduced due to fewer combinations of connections, as occurs
similarly during seizures (Figure 1). D, same subject but using the scalp EEG
signals, showing similar results to those obtained with iEEG (evaluated at
central frequency of 4 Hz). Figure 3:  Time course of the entropy of the
number of configurations of connected MEG signals before, during and after a
generalised absence seizure. MEG signal from one channel is shown at top, the
ictus (‘Sz’) occurs towards the end of the 2 minute recording. Notice the drop
in entropy during the seizure. Figure 4:  Proposed general scheme of the
relation between global brain connectivity and behavioural states. Normal
alertness resides at the top of the curve representing the number of
configurations of connections the system can adopt, or the associated entropy.
The maximisation of the configurations (microstates) provides the variability
in brain activity needed for normal sensorimotor action. Abnormal, or
unconscious states, are located farther from the top, and are characterised by
either large or small number of “connected” networks therefore exhibiting
lower number of microstates (hence lower entropy) that are not optimal for
sensorimotor processing.


11institutetext: University of British Columbia 22institutetext: Princeton
University
22email: {chunyinh<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Efficient Medical Image Assessment via Self-supervised Learning

Chun-Yin Huang 11 Qi Lei 22 Xiaoxiao Li 11

###### Abstract

High-performance deep learning methods typically rely on large annotated
training datasets, which are difficult to obtain in many clinical applications
due to the high cost of medical image labeling. Existing data assessment
methods commonly require knowing the labels in advance, which are not feasible
to achieve our goal of ‘knowing which data to label.’ To this end, we
formulate and propose a novel and efficient data assessment strategy,
EXponentiAl Marginal sINgular valuE ($\mathsf{EXAMINE}$) score, to rank the
quality of unlabeled medical image data based on their useful latent
representations extracted via Self-supervised Learning (SSL) networks.
Motivated by theoretical implication of SSL embedding space, we leverage a
Masked Autoencoder [8] for feature extraction. Furthermore, we evaluate data
quality based on the marginal change of the largest singular value after
excluding the data point in the dataset. We conduct extensive experiments on a
pathology dataset. Our results indicate the effectiveness and efficiency of
our proposed methods for selecting the most valuable data to label.

##  1 Introduction

Artificial intelligence (AI) such as deep learning has became a powerful tool
for medical image analysis. Its success relies on the availability of abundant
high quality dataset. However, medical images collected from different sources
vary in their quality due to the various imaging devices, protocols and
techniques. When trained with low-quality data, AI models can be compromised.
Furthermore, labeling medical images for AI training requires domain experts
and is usually costly and time consuming. Therefore, it is demanding to have
an automated framework to effectively assess and screen data quality before
data labeling and model training.

There are numerous definitions of data quality. Data is generally considered
to be of high quality if “fit for [its] intended uses in operations, decision
making and planning.” [16, 4, 5]. In the context of training an AI predictive
model, good data are the fuel of AI. Namely, data with better quality can help
obtain higher prediction accuracy. However, how to quantitatively assess
data’s quality for AI tasks is under-explored. Previous works [10, 6] mainly
propose to estimate data values in the context of supervised machine learning,
which requires knowledge of labels and repeated training of a target utility.
Such setting lacks practical value as data labels are typically not available
at the data preparation stage for data privacy, labeling cost, and
computational efficiency concerns. Differently, we aim to develop a cost-
effective scheme for data assessment in the context of unsupervised learning
to tackle the limitations of the existing methods, in which no labeling is
required during assessment.

A trending and powerful unsupervised representation learning strategy is self-
supervised learning (SSL). SSL solves auxiliary pretext tasks without
requiring labeled data to learn useful semantic representations. These pretext
tasks are created solely using the input features, such as predicting a
missing image patch [8], recovering the color channels of an image from
context [19], predicting missing words in texts [12], forcing the similarity
of the different views of images [1, 7], etc. Motivated by the recent
discovery that SSL could embed data into linearly separable representations
under proper data assumptions [13, 17], we show that ‘good’ and ‘bad’ data can
be distinguished by examining the change of the data representation matrices’
singular value by removing a certain data point.

In this work, we tackle a practically demanding yet challenging problem —
medical image assessment (also referred as data assessment in this work). To
this end, we develop a novel and efficient pipeline for medical image
assessment without knowing data labels. As shown in Fig 1, we propose a new
metric, EXponentiAl Marginal sINular valuE ($\mathsf{EXAMINE}$) score to
evaluate the value (or referred as quality) of individual data by first using
SSL to extract the features, and then calculate the value of the data using
Singular Value Decomposition (SVD). $\mathsf{EXAMINE}$ scores are useful in
indicating the essential data to be annotated, which can not only abundantly
reduce the effort in manual labeling but also mitigate the negative effect of
mislabeled data, and further improve the target model. Our chief contributions
are summarized as follows:
\- We are the first to show the feasibility of using an unsupervised learning
framework to assess medical data by utilizing SSL and SVD, which is a more
cost-efficient and practical method to evaluate data compared to previous
work.
\- $\mathsf{EXAMINE}$ can assess data without knowing the label, which reduces
annotation efforts and the chance of mislabeling.
\- We conduct experiments on the simulated medical dataset to demonstrate the
feasibility of using $\mathsf{EXAMINE}$ scores to distinguish data with
different qualities and show comparable performance to previous supervised
learning based works.

##  2 Preliminaries

###  2.1 Supervised-learning-based Data Assessment

The goal of data assessment is using a valuation function to map an input data
to a single value that indicate its quality. Supervised-learning-based data
assessments assume knowing a labeled dataset
$\mathcal{N}^{l}=\\{(x_{i},y_{i})|i\in[N],x_{i}\in\mathcal{X},y_{i}\in\mathcal{Y}\\}$
where $N$ is the number of the labeled data, an utility model
$f:\mathcal{N}^{l}\mapsto\mathcal{Y}$, a held-out labeled testing set
$\mathcal{N}^{t}=\\{(x^{\prime}_{i},y^{\prime}_{i})|i\in[M],x^{\prime}_{i}\in\mathcal{X},y^{\prime}_{i}\in\mathcal{Y}\\}$
where $M$ is the number of the testing data, and a value function
$V:(f,\mathcal{N}^{l},\mathcal{N}^{t})\mapsto\mathbb{R}$ (e.g., the accuracy
of $\mathcal{N}^{t}$ evaluated by $f$ that is trained on $\mathcal{N}^{l}$ ).
The simplest assessment metric is by performing leave-one-out (LOO) on the
training set and calculating the performance differences on the testing set.
The $i$-th data samples value is defined as:

$\displaystyle\phi_{i}^{\rm
LOO}=V_{f}(\mathcal{N}^{l})-V_{f}(\mathcal{N}^{l}\backslash\\{i\\}).$ (1)

A more advanced but computational costly approach is Data Shapley [6]. Shapley
value for data valuation resembles a game where training data points are the
players and the payoff is defined by the goodness of fit achieved by a model
on the testing data. Given a subset $S$, let $f_{S}(\cdot)$ be a model trained
on $\mathcal{S}$. Then Shapley value of a data point
$(x_{i},y_{i})\in\mathcal{N}$ is defined as:

$\displaystyle\phi_{i}^{\rm
SHAP}=\mathop{\sum}\limits_{{S\subseteq\mathcal{N}\backslash\left\\{{x_{i}}\right\\}}}\frac{{V_{f}\left({S\cup\left\\{{x_{i}}\right\\}}\right)-V_{f}\left(S\right)}}{{\left({\begin{array}[]{*{20}c}{\left|\mathcal{N}^{l}\right|-1}\\\
{\left|S\right|}\\\ \end{array}}\right)}},$ (4)

where $V_{f}(S)$ is the performance of the utility model $f$ trained on subset
$S$ of the data. Suppose each training of $f$ takes time $T$, the
computational complexity of Eq (1) and Eq (4) is $\mathcal{O}(TN)$ and
$\mathcal{O}(T2^{N})$111In practice, there are approximation methods for
calculating Shapley value, but the it still requires around
$\mathcal{O}(T\rm{poly}(N))$ [11]., respectively. Also, training a deep
utility function (e.g., neural networks) leads to a large $T$.

###  2.2 Formulation of Unsupervised-learning-based Data Assessment

#### Motivation story

Labeling is costly and time consuming in many medical imaging tasks. AI
developers may want to pay for labeling some data points to train a particular
machine learning model. In such a scenario, supervised-learning-based methods
(Sec 2.1) cannot fulfill the aim. Therefore, algorithms that can automatically
identify low quality data before labeling data are highly desired.

To address computational issue and demand for task/label-agnostic data
quality, we propose to conduct quantitative data quality assessment via
unsupervised learning. Different from the formulations of LOO (Eq. (1)) and
Shapley value (Eq. (4)), here we propose a new problem formulation to infer
$i$-th data’s quality by assigning it a value
$\phi_{i}:(\mathcal{X},i)\mapsto\mathbb{R}$ using unlabeled data only.

##  3 Our Method

###  3.1 Theoretical Implication

Our proposed $\mathsf{EXAMINE}$ is well motivated by representation theory of
SSL. We start by restating Theorem 3.1 proved in [13] that under proper
assumptions, the embedded space obtained by the reconstruction-based SSL
strategy forms a linearly separable space of the embedded feature and a
related task. Then, Remark 1 presents how we use Theorem 3.1 to guide the
design of $\mathsf{EXAMINE}$.

######  Theorem 3.1 (informal [13])

For two views of a data $X_{1},X_{2}\in\mathcal{X}$ and their classification
label $Y\in\mathbb{R}^{k}$. Under the class conditional independence
assumption, i.e., $X_{1}\perp X_{2}|Y$, for some $w\in\mathbb{R}^{m\times k}$
the representation $\psi^{*}:\mathcal{X}\mapsto\mathbb{R}^{m}$ that minimizes
a reconstruction loss
${\mathcal{L}(\psi)}=\mathbb{E}_{(X_{1},X_{2})}\left[\|X_{1}-\psi(X_{2})\|^{2}\right]$
satisfies

$\displaystyle w^{\top}\psi^{*}(X_{1})=\mathbb{E}[Y|X_{1}].$

###### Remark 1

Theorem 3.1 indicates two desired properties with ‘good’ data that
(approximately) satisfy class-conditional independence. First, the data will
have a good geometric property in the learned representation space, namely
they become clusters that are (almost) linearly separable (by $w$). Second,
the learned representation has variance (top singular value of its covariance
matrix) controlled by that of $\mathbb{E}[Y|X_{1}]$, which is very small
(since the label is almost determined entirely by the image itself). Such
properties on the reconstruction-based SSL embedding space are not satisfied
for ‘bad’ data. Therefore when observing data $X_{1}$ is noisy or with low
quality, $\psi^{*}(X_{1})$ tends to have higher variance.

Figure 1: Proposed pipeline for $\mathsf{EXAMINE}$ data assessment. (a) Using
the state-of-the art reconstruction-based SSL strategy, MAE [8] architecture
for pre-training an representation extractor (encoder). (b) $\mathsf{EXAMINE}$
first utilizes the pre-trained encoder to extract semantic features $f_{S}$
and $f_{S\backslash\\{i\\}}$ from input data $D_{S}$ and
$D_{S\backslash\\{i\\}}$, where $D_{S\backslash\\{i\\}}$ denotes input data
$D_{S}$ without data point $i$. The features then pass the SVD module to find
the largest singular values $\lambda_{S}$ and $\lambda_{S\backslash\\{i\\}}$.
The $\mathsf{EXAMINE}$ score of data point $i$ is defined as Eq. (5).

###  3.2 Data Assessment on Singular Value

As shown in the Fig.1(b), to assess data from a dataset $D_{S}\in R^{N\times
C}$, where $N$ is the number of data and $C$ is the dimension of data, we
denote the dataset without $i$-th data point as $D_{S\backslash\\{i\\}}\in
R^{(N-1)\times C}$. To begin with, we employ SSL and the unlabeled data to
train an encoder that is able to extract the low-dimensional semantic
information. We denote the representation of the SSL embedding space of
$D_{S}$ and $D_{S\backslash\\{i\\}}$ as $f_{S}$ and $f_{S\backslash\\{i\\}}$,
respectively. Lastly, we perform SVD on both feature representations, and use
the largest singular values ($\lambda_{S}$ and $\lambda_{S\backslash\\{i\\}}$)
as the assessment indicator, that is, removing a ‘good’ data point $i$ results
in small change in the top singular value of embedded data representation $f$
(explained in Sec. 3.1). Thus, the $\mathsf{EXAMINE}$ score is defined as

$\displaystyle\phi_{i}^{\rm
EXAMINE}=\exp\left(-(\lambda_{S}-\lambda_{S\backslash\\{i\\}})\right),$ (5)

where $\phi_{i}^{\rm EXAMINE}\in(0,1)$ and a larger $\phi_{i}^{\rm EXAMINE}$
indicates better data quality222$\lambda_{S}>\lambda_{S\backslash\\{i\\}}$ is
for sure given the properties of singular value.. Note that $\mathsf{EXAMINE}$
is also a leave-one-out strategy but evalated on the change of the largest
singular value.

We claim two advantages of using SVD for data assessment. First, by performing
the SVD-based evaluation, we do not need any knowledge about the corresponding
labels. This is not only the primary difference from previous methods, but
also a perfect fit to our problem set up - finding good data to be labeled.
Second, unlike previous methods (i.e., LOO and Data Shapely, see Sec. 2.1)
that rely on extensively training a new model for different data combinations,
our proposed method is efficient by performing SVD once for each data point
without additional model training after the SSL encoder has been trained
offline.

###  3.3 Forming Embedding Space using Masked Auto-encoding

As the raw medical images are high-dimensional and have spurious features
(e.g., density, light, dose) that are irrelevant to their labels, directly
applying SVD to them cannot capture task-related variance. Based on our theory
developed on reconstruction-based SSL (Sec. 3.1), we utilize a state-of-the-
art reconstructed-based strategy, Masked Auto-Encoder (MAE) [8] to learn
lower-dimensional semantic feature embedding. As shown in Fig. 1(a), MAE
utilizes state-of-the-art image classification framework, Vision Transformer
(ViT) [3], as the encoder for semantic feature extraction, and uses a lighter
version of ViT as decoder. It first divides an input image into patches,
randomly blocks a certain percentage of image patches, and then feeds them
into the autoencoder architecture. By blocking out a large amount of image
patches, the model is forced to learn a more complete representation. With the
aim of positional embedding and transformer architecture, MAE is able to
generalize the relationship between each image patch and obtain the semantic
information among the whole image, which achieves the state-of-the-art
performance in self-supervised image representation training. This also
reduces the correlation between spurious features and labels, compared to the
traditional dimension reduction methods [2].

##  4 Experiment

###  4.1 Experiment setup and dataset

We evaluate $\mathsf{EXAMINE}$ on a binary classification task for PCam [18],
a microscopic dataset (image size $96\times 96$) for identifying metastatic
tissue in histopathologic scans of lymph node sections. Since noise is usually
the main corruption in medical images, we add non-zero mean Gaussian noise to
a portion of the data to simulate real world scenario.

We split the dataset into four disjoint sets following the scale of [6]:
\- SSL Pre-Training Set and Assessed Set: 160,000 and 500 unlabeled data
points randomly sampled from PCam, respectively. We add 4 different level
($\mathcal{N}(\delta,\delta\times m)$, where $m$ is the mean of the dataset
and $\delta=\\{0.1,0.3,0.5,1\\}$) of noise to 60000 data points of SSL Pre-
Training Set and 400 data points in Assessed Set.

\- Clean Train Set and Validation Set: 100 and 20,000 labeled data points
randomly sampled from PCam. They are used to validate the data selection in an
example downstream task after obtaining $\mathsf{EXAMINE}$ scores
($\phi_{i}^{\rm EXAMINE}$).

The experiments are run on NVIDIA GeForce RTX 3090 Graphics card with PyTorch.
For MAE training, we select Cosine Annealing LR scheduler [14] and AdamW LR
optimizer [15] with weight decay 0.05 and momentum $\\{0.9,0.95\\}$. We train
MAE for 200 epochs using batch size 256 ,and set image size 72, patch size 8,
masking ratio $40\%$. As indicated in [1, 7, 8] that SSL training requires a
large amount of data, we begin with training a MAE on SSL Pre-Training Set. To
ensure the training stability, we first train MAE without any noisy data, and
finetune it afterward. This step is to distinguish our setting from detection
out of distribution samples. After pretraining the MAE encoder, we use the
frozen encoder layers as our backbone to extract the low dimensional
representations. All of our experiments are repeated five times with different
random seeds and we report the mean value of the five trials.

(a)

(b)

Figure 2: Proof of concept and comparison with baseline embedding methods.
Observe that $\mathsf{EXAMINE}$ scores ($\phi_{i}^{\rm EXAMINE}$) get lower
when noise level increases (a), and the $\mathsf{EXAMINE}$ scores
($\phi_{i}^{\rm EXAMINE}$) in (b) shows that reconstrucion-based SSL methods
perform better in separating different noise levels.

###  4.2 Proof of Concept with ‘Ground-truth’

To validate the correctness of ranking the samples in Assessed Set, we plot
the distributions of $\mathsf{EXAMINE}$ scores ($\phi_{i}^{\rm EXAMINE}$) at
different data corruption levels of noisy setting in Fig 2(a) Note that
$\phi_{i}^{\rm EXAMINE}\in(0,1)$ and data point $i$ with larger $\phi_{i}^{\rm
EXAMINE}$ indicates that it is considered a good data. Specifically,
$\mathsf{EXAMINE}$ score ($\phi_{i}^{\rm EXAMINE}$) approaching 1 indicates
the difference between the top singular values are small, thus it will be
considered good data. The $\mathsf{EXAMINE}$ scores ($\phi_{i}^{\rm EXAMINE}$)
for the high-quality data ($\delta=0$) are close to 1 and significantly higher
than the corrupted data. The $\mathsf{EXAMINE}$ scores ($\phi_{i}^{\rm
EXAMINE}$) of data with low-level corruption ($\delta=0.5$) are also separable
from those with high-level corruption($\delta=1$).

###  4.3 Comparison with Alternative Embedding Methods

We investigate the alternative feature encoders and compare their performance
with MAE. Specifically, we replace MAE with SimCLR [1] and BYOL [7], two
alternative SSL algorithms, and Autoencoder (AE) [9], a naïve reconstruction-
based embedding strategy. SimCLR learns embedding by enforcing the closeness
of an image and its augmented views while enlarging the distance from other
images in the dataset (or batch). BYOL regularizes the multi-views of an image
without sampling negative samples by training two similar networks (the online
network and the target network) simultaneously. We use the same strategy as
training MAE for these alternative encoders. Fig. 2(b) shows that using MAE
embedding to calculate $\mathsf{EXAMINE}$ scores ($\phi_{i}^{\rm EXAMINE}$)
provides the best separability for the clean data from the corrupted data. The
reason is that MAE best satisfies the theoretical conditions that support our
proposal (Theorem 3.1). AE is second to MAE, but separation boundaries are
less clear.

(a)

(b)

(c)

(d)

Figure 3: Comparison with baseline data valuation methods. Adding good data(a)
and removing bad data(d) should increase accuracy. Adding bad data(b) and
removing good data(c) should result in accuracy drop. We conclude that
$\mathsf{EXAMINE}$ can help model training by identifying good and bad data.

###  4.4 Comparison with Baseline Data Valuation Methods

We compare $\mathsf{EXAMINE}$ with supervised data valuation methods, LOO (Eq.
1 and Truncated Monte Carlo (TMC) version of Data Shapley (Eq. 4) [6], as well
as a baseline method that randomly assigns data values. All these methods are
applied on Assessed Set’s features extracted by the pre-trained MAE. We design
four experiments to evaluate how selecting data using the different data
assessment methods can affect the classification accuracy. We report the
averaged test accuracy on Validation Set using logistic regression models
(LRM).

Fig. 3(a) and Fig. 3(b) show the results of adding data for training. We start
with a LRM trained on the small Clean Training Set, and then add good/bad data
from Assessed Set following the descending/ascending orders of their data
values. Our results show that adding data with high $\mathsf{EXAMINE}$ score
($\phi_{i}^{\rm EXAMINE}$) achieves comparable accuracy curve as Data Shapley,
while adding data with low $\mathsf{EXAMINE}$ scores ($\phi_{i}^{\rm
EXAMINE}$) results in similar curve as Data Shapley in the beginning and
overall lies in between Data Shapley and LOO. This indicates that
$\mathsf{EXAMINE}$ score ($\phi_{i}^{\rm EXAMINE}$) is able to identify what
data to be labeled and added to training set. Fig. 3(c) and Fig. 3(d) show the
results of removing data for training. We first train LRM on Assessed Set, and
then remove good/bad data following the descending/ascending orders of their
data values. Our result shows that removing high/low $\mathsf{EXAMINE}$ score
($\phi_{i}^{\rm EXAMINE}$) data results in accuracy curve that is slightly
worse than Data Shapley. We would like to emphasize Data Shapley uses utility
function which requires labels to determine the data value, while
$\mathsf{EXAMINE}$ score ($\phi_{i}^{\rm EXAMINE}$) is calculated only on data
itself, which is more efficient in real-world scenario. Overall,
$\mathsf{EXAMINE}$ shows comparable(or at best slightly worser) data
assessment performance to Data Shapley without knowing the labels of the data.

In addition to successfully providing correct inspection on data quality, our
method significantly reduces the computational cost without requiring training
utility functions. The running time to obtain the data values using
$\mathsf{EXAMINE}$, LOO, and TMC Data Shapley for the whole Assessed Set under
our experiment setting are 23 seconds, 10 minutes, and 460 minutes,
respectively333The running time for LOO and Data Shapley can significantly
increase if we use a deep neural network as the utility model..

##  5 Discussion and Conclusion

We present a new and efficient unsupervised data evaluation method,
$\mathsf{EXAMINE}$ scores $\phi_{i}^{\rm EXAMINE}$, to assess data quality.
With the help of MAE encoder, we can map data to the provable low-dimensional
embedding space. The marginal differences on the largest singular value of
data representation matrices can effectively separate data at different
quality levels and achieve comparable performance with supervised data
valuation methods when considering a specific task. This work takes a novel
approach to promote AI in healthcare by identifying low quality data. We plan
to test on larger scale medical datasets and collect domain experts’
evaluations in the future.

## Acknowledgement

This work is supported in part by the Natural Sciences and Engineering
Research Council of Canada (NSERC) and NVIDIA Hardware Award.

## References

  * [1] Chen, T., Kornblith, S., Norouzi, M., Hinton, G.: A simple framework for contrastive learning of visual representations. In: International conference on machine learning. pp. 1597–1607. PMLR (2020)
  * [2] Chen, Y., Wei, C., Kumar, A., Ma, T.: Self-training avoids using spurious features under domain shift. Advances in Neural Information Processing Systems 33, 21061–21071 (2020)
  * [3] Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., et al.: An image is worth 16x16 words: Transformers for image recognition at scale. arXiv preprint arXiv:2010.11929 (2020)
  * [4] Fadahunsi, K.P., Akinlua, J.T., O’Connor, S., Wark, P.A., Gallagher, J., Carroll, C., Majeed, A., O’Donoghue, J.: Protocol for a systematic review and qualitative synthesis of information quality frameworks in ehealth. BMJ open 9(3), e024722 (2019)
  * [5] Fadahunsi, K.P., O’Connor, S., Akinlua, J.T., Wark, P.A., Gallagher, J., Carroll, C., Car, J., Majeed, A., O’Donoghue, J.: Information quality frameworks for digital health technologies: systematic review. Journal of medical Internet research 23(5), e23479 (2021)
  * [6] Ghorbani, A., Zou, J.: Data shapley: Equitable valuation of data for machine learning. In: Chaudhuri, K., Salakhutdinov, R. (eds.) Proceedings of the 36th International Conference on Machine Learning. Proceedings of Machine Learning Research, vol. 97, pp. 2242–2251. PMLR (09–15 Jun 2019)
  * [7] Grill, J.B., Strub, F., Altché, F., Tallec, C., Richemond, P., Buchatskaya, E., Doersch, C., Avila Pires, B., Guo, Z., Gheshlaghi Azar, M., et al.: Bootstrap your own latent-a new approach to self-supervised learning. Advances in Neural Information Processing Systems 33, 21271–21284 (2020)
  * [8] He, K., Chen, X., Xie, S., Li, Y., Dollár, P., Girshick, R.: Masked autoencoders are scalable vision learners. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 16000–16009 (2022)
  * [9] Hinton, G.E., Salakhutdinov, R.R.: Reducing the dimensionality of data with neural networks. science 313(5786), 504–507 (2006)
  * [10] Jia, R., Dao, D., Wang, B., Hubis, F.A., Hynes, N., Gürel, N.M., Li, B., Zhang, C., Song, D., Spanos, C.J.: Towards efficient data valuation based on the shapley value. In: The 22nd International Conference on Artificial Intelligence and Statistics. pp. 1167–1176. PMLR (2019)
  * [11] Jia, R., Sun, X., Xu, J., Zhang, C., Li, B., Song, D.: An empirical and comparative analysis of data valuation with scalable algorithms (2019)
  * [12] Kenton, J.D.M.W.C., Toutanova, L.K.: Bert: Pre-training of deep bidirectional transformers for language understanding. In: Proceedings of naacL-HLT. pp. 4171–4186 (2019)
  * [13] Lee, J.D., Lei, Q., Saunshi, N., Zhuo, J.: Predicting what you already know helps: Provable self-supervised learning. Advances in Neural Information Processing Systems 34 (2021)
  * [14] Loshchilov, I., Hutter, F.: Sgdr: Stochastic gradient descent with warm restarts. arXiv preprint arXiv:1608.03983 (2016)
  * [15] Loshchilov, I., Hutter, F.: Decoupled weight decay regularization. arXiv preprint arXiv:1711.05101 (2017)
  * [16] Redman, T.C.: Data driven: profiting from your most important business asset. Harvard Business Press (2008)
  * [17] Tosh, C., Krishnamurthy, A., Hsu, D.: Contrastive learning, multi-view redundancy, and linear models. In: Algorithmic Learning Theory. pp. 1179–1206. PMLR (2021)
  * [18] Veeling, B.S., Linmans, J., Winkens, J., Cohen, T., Welling, M.: Rotation equivariant cnns for digital pathology. In: International Conference on Medical image computing and computer-assisted intervention. pp. 210–218. Springer (2018)
  * [19] Zhang, R., Isola, P., Efros, A.A.: Colorful image colorization. In: European conference on computer vision. pp. 649–666. Springer (2016)


# Adaptive Model Predictive Control of Wheeled Mobile Robots ††thanks: The
authors are graduate students in the department of Mechanical Engineering at
the University of California at Berkeley, Berkeley, CA 94720, USA.
Email:{1nikhilps,2tamara_perreault,3tvoth,
<EMAIL_ADDRESS>
Nikhil Potu Surya Prakash1 Tamara Perreault2 Trevor Voth3 and Zejun Zhong4

###### Abstract

In this paper, a control algorithm for guiding a two wheeled mobile robot with
unknown inertia to a desired point and orientation using an Adaptive Model
Predictive Control (AMPC) framework is presented. The two wheeled mobile robot
is modeled as a knife edge or a skate with nonholonomic kinematic constraints
and the dynamical equations are derived using the Lagrangian approach. The
inputs at every time instant are obtained from Model Predictive Control (MPC)
with a set of nominal parameters which are updated using a recursive least
squares algorithm. The efficacy of the algorithm is demonstrated through
numerical simulations at the end of the paper.

##  I Introduction

Model Predictive Control (MPC) as a method for designing optimal control
inputs for dynamical systems has gained a lot of attention in the past few
years. With the increase in computational power, using MPC to obtain optimal
control inputs has become more feasible. Since MPC is a closed loop control
strategy, it can robustly and optimally compensate for minor disturbances and
parametric uncertainties which cannot be accounted for while computing open
loop optimal control inputs.

But when the dynamic model of the system is not known or when the model is
known but the physical parameters in the model are not known, the control
inputs from MPC do not necessarily achieve the desired objective. It is
therefore essential to first learn the model from the behaviour of the system
to various inputs and then use MPC to obtain optimal control inputs. In this
paper, an Adaptive Model Predictive Control (AMPC) framework for systems whose
physical model is known but its parameters are unknown is presented so that
the learning of parameters and obtaining the optimal control inputs is done
simultaneously. More specifically, a system in which the dynamical model is
linear in its physical parameters. Such a situation is very common when the
robot has to picks up and carry a payload with unknown mass and mass moment of
inertia like a warehouse robot which needs to pick up and transfer payloads of
different masses and sizes. For such tasks, a non-adaptive control law might
not do a good job and hence a necessity to learn more information about the
system arises.

Paper structure: In Section II, we discuss the nonlinear dynamics used to
model the wheeled robot with nonholonomic constraints. In Section III, we
present the Adaptive MPC Algorithm used to control the robot along a
trajectory and the parameter adaptation algorithm used to learn the parameters
of the dynamical system. In Section IV, we discuss the performance of the
controller as modeled by simulations.

##  II Dynamics of a wheeled mobile robot

In this section, the dynamic evolution of the states of the wheeled mobile
robot will be presented. The mechanical structure of the wheeled mobile robot
is such that it is driven as a differential drive with two wheels actuated by
motor torques at the rear and a smooth caster in the front keeps the robot
stable. The dynamics are based on the dynamic model of a unicycle with a
nonholonomic constraint that the admissible velocity of the robot is only
along its orientation, thereby implying that the velocity in the direction
perpendicular to its orientation is zero. The mathematical structure of the
robot is isomorphic to that of a knife edge actuated by a thrust for motion
control along its heading and a torque along the vector perpendicular to the
motion plane of the robot to control its orientation. Using the isomorphic
structure between the dynamics of the wheeled mobile robot and the knife edge,
the motor torques on the wheels of the robot can be easily calculated using a
one-to-one transformation. Hence to simplify the analysis and notation, a
dynamic knife edge model will be used in the rest of the paper. The dynamics
are described using the Lagrangian formulation.
The system in consideration is similar to a knife in structure but with more
assumptions on the constraints that the frictional capacity available to keep
the skate from moving in the direction perpendicular to its
orientation.Further we assume that the frictional force along the direction of
motion allows only pure rolling of the wheels and no slippage is allowed.
We define two frames, an inertial frame (with subscript A) and a body fixed
frame (with subscript B) in our analysis.The transformation between these two
frames is achieved through the yaw angle.
The kinetic energy of the system with states
$[x,y,\psi,\dot{x},\dot{y},\dot{\psi}]$ is given by

$T=\frac{1}{2}m\dot{x}^{2}+\frac{1}{2}m\dot{y}^{2}+\frac{1}{2}J\dot{\psi}^{2}$
(1)

where $x$ and $y$ are the coordinates of the robot with respect to an inertial
frame, $\psi$ is its heading angle and the rest are their respective
derivatives.
We assume that the motion of the knife edge is restricted to a plane with the
same gravitational potential throughout and hence there will not be any change
in the potential energy of the system as the knife edge moves. But this
analysis can easily be extended to a case where gravitational potential need
not be constant along the path. Since the other potentials remain constant,
the kinetic energy will be equal to the Lagrangian ($L:=T-V$) of the system.
The constraint on the motion is given by the condition that the velocity in
the direction $j_{B}$ $i.e.,$ along the direction perpendicular to its
orientation is $0$

$-\dot{x}sin\psi+\dot{y}cos\psi=0$ (2)

We control the system by giving two inputs on the system. A thrust $R$ in the
direction of the velocity and a moment $M$ to change the orientation of the
knife edge.
The equations of motion are derived using the Lagrangian formulation by
introducing Lagrangian multipliers ($\lambda$) since the system is
nonholonomic. One can even use the Appelian approach using the acceleration
energy to arrive at the equations of motion. The generalized forces on the
system can be derived by writing the virtual work ($\delta W$) done by the
control inputs on the system.

$\displaystyle\delta W=$ $\displaystyle(\delta xcos\psi+\delta
ysin\psi)(R-bv)$ (3) $\displaystyle+\delta\psi(M-c\dot{\psi})$

Where $v=\dot{x}cos\psi+\dot{y}sin\psi$ is the longitudinal velocity, $b$ and
$c$ are the drag coefficients corresponding to the linear and angular motions
respectively. Therefore the generalized forces $Q_{x}$ , $Q_{y}$ and
$Q_{\psi}$ are equal to $(R-bv)cos\psi$, $(R-bv)sin\psi$ and $M-c\dot{\psi}$
respectively.

$\displaystyle\frac{d}{dt}\frac{\partial{L}}{\partial{\dot{x}}}-\frac{\partial{L}}{\partial{x}}-\lambda
a_{1}=Q_{x}$ (4)
$\displaystyle\frac{d}{dt}\frac{\partial{L}}{\partial{\dot{y}}}-\frac{\partial{L}}{\partial{y}}-\lambda
a_{2}=Q_{y}$
$\displaystyle\frac{d}{dt}\frac{\partial{L}}{\partial{\dot{\psi}}}-\frac{\partial{L}}{\partial{\psi}}-\lambda
a_{3}=Q_{\psi}$

The coefficients ($a_{i}$) to the Lagrange multipliers are the same as the
coefficients of the states in the constraint equations. Substituting all the
quantities and simplifying we get.

$m\ddot{x}+\lambda sin\psi=(R-bv)cos\psi$ (5) $m\ddot{y}-\lambda
cos\psi=(R-bv)sin\psi$ (6) $J\ddot{\psi}=M-c\dot{\psi}$ (7)

From Newton’s second law, it can easily be seen that the lateral frictional
force acting on the skate is the same as the Lagrangian multiplier $\lambda$.

$m\ddot{x}cos\psi+m\ddot{y}sin\psi=R-bv$ (8)
$-m\ddot{x}sin\psi+m\ddot{y}cos\psi=\lambda$ (9)

By differentiating the constraint equation (2) we get

$-\ddot{x}sin\psi-\dot{x}\dot{\psi}cos\psi+\ddot{y}cos\psi-\dot{y}\dot{\psi}sin\psi=0$
(10)

Therefore

$m\dot{\psi}(\dot{x}cos\psi+\dot{y}sin\psi)=\lambda$ (11)

Here $\dot{x}cos\psi+\dot{y}sin\psi$ is the longitudinal velocity of the knife
edge. Denoting the longitudinal velocity by $v$, we get

$m\dot{\psi}v=\lambda$ (12)

Using the expression for longitudinal velocity, (8) can be simplified as

$m\dot{v}=R-bv$ (13)

We can combine all these equations into a simple state space form as follows

$\begin{bmatrix}\dot{x}\\\ \dot{y}\\\ \dot{\psi}\\\ \dot{v}\\\
\dot{\omega}\end{bmatrix}=\begin{bmatrix}vcos(\psi)\\\ vsin(\psi)\\\ \omega\\\
\frac{(R-bv)}{m}\\\ \frac{(M-c\dot{\psi})}{J}\end{bmatrix}$ (14)

The nonholonomic constraint that the robot cannot move in the direction
perpendicular to its orientation given by (2) is satisfied at every instant
and can be easily verified by substituting (14) in (2). To simplify the
computation process by approximating the infinite dimensional optimization
problem to a finite dimensional optimization problem, the dynamics will be
discretized as follows and will be used to model the behavior of the robot in
the rest of the paper. The discretized equations are obtained using Euler
discretization with a sampling interval $\Delta t$. A more complicated yet
better approximation can be obtained using first order hold or any other
similar higher order interpolations for the inputs between the sampling
points. The discretized dynamics are given by

$\begin{bmatrix}x_{t+1}\\\ y_{t+1}\\\ \psi_{t+1}\\\ v_{t+1}\\\
\omega_{t+1}\end{bmatrix}=\begin{bmatrix}x_{t}+v_{t}cos(\psi_{t})\Delta t\\\
y_{t}+v_{t}sin(\psi_{t})\Delta t\\\ \psi_{t}+\omega_{t}\Delta t\\\
v_{t}+\frac{(R-bv_{t})\Delta t}{m}\\\ \omega_{t}+\frac{(M-c\omega_{t})\Delta
t}{J}\\\ \end{bmatrix}$ (15)

An equivalent representation of (15) as following will be used in the Adaptive
MPC algorithm.

$\begin{bmatrix}x_{t+1}\\\ y_{t+1}\\\ \psi_{t+1}\\\ v_{t+1}\\\
\omega_{t+1}\end{bmatrix}=\begin{bmatrix}x_{t}+v_{t}cos(\psi_{t})\Delta t\\\
y_{t}+v_{t}sin(\psi_{t})\Delta t\\\ \psi_{t}+\omega_{t}\Delta t\\\
\alpha_{v}v_{t}+\beta_{v}R\\\ \alpha_{\omega}\omega_{k}+\beta_{\omega}M\\\
\end{bmatrix}$ (16)

Where

$\displaystyle\alpha_{v}=1-\frac{b\Delta t}{m}$ (17)
$\displaystyle\beta_{v}=\frac{\Delta t}{m}$
$\displaystyle\alpha_{\omega}=1-\frac{c\Delta t}{J}$
$\displaystyle\beta_{\omega}=\frac{\Delta t}{J}$

It can be easily seen that by treating
$\theta=[\alpha_{v},\beta_{v},\alpha_{\omega},\beta_{\omega}]^{T}$ as proxy
system parameters, the original system parameters $[m,b,J,c]^{T}$can be
calculated using a simple one-to-one transformation. This transformation (17)
makes the system dynamics a linear function of the system parameters. This
property will become handy to use existing algorithms during the
identification of the system parameters. These dynamics can be compactly
written as

$X_{t+1}=f(X_{t},u_{t},\theta)$ (18)

where $X_{k}\in\mathcal{X}\subseteq\mathbb{R}^{5}$ is the state vector and
$u_{k}\in\mathcal{U}\subseteq\mathbb{R}^{2}$ is the vector of control inputs
$[R,M]^{T}$.The equations in (15) will be used to describe the dynamics of the
robots.

Figure 1: A sample image of a wheeled mobile robot.

##  III Adaptive MPC

In the dynamic model derived in section II, we start with no knowledge of the
system parameters like the mass, mass moment of inertia and the drag
coefficients of the robot and estimate them as we go along with the knowledge
of the states of the system with a parameter adaptation algorithm and feed it
as parameters to an optimization problem at every instance to obtain the
optimal control inputs. We assume that we have full access to the states of
the robot.
A typical MPC problem is formulated as follows

$\displaystyle U^{*}_{t}(x(t))=$
$\displaystyle\underset{U_{t}}{\text{argmin}}\sum_{k=0}^{N-1}q(X_{t+k},u_{t+k})$
(19) subject to $\displaystyle X_{t}=X(t)$ $\displaystyle(c1)$ $\displaystyle
X_{t+k+1}=f(X_{t+k},u_{t+k},\theta)$ $\displaystyle(c2)$ $\displaystyle
u_{t+k}\in\mathcal{U}$ $\displaystyle(c3)$ $\displaystyle
U_{t}=\\{u_{t},u_{t+1},\dots,u_{t+N-1}\\}$ $\displaystyle(c4)$

where the constraint (c1) is the measurement, (c2) is the system model as a
function of parameters, (c3) is the admissible input set and (c4) is the set
of optimization variables. The control input at $t$ would be the first
component of $U^{*}_{t}(x(t))$ $i.e.,$ $u^{*}_{t}$. It is worth noting that
the constraints on the state are not included here as the system parameters
are not known and guaranteeing the states to be confined to a desired set
might not be possible. Instead one could use soft constraints for state
constraints
MPC algorithm in (19) is not implementable just by itself as the system
parameters are unknown. So, we start with a guess of the system parameters and
calculate the optimal control corresponding to the guess and update the guess
at every instant using a parameter adaptation algorithm. MPC and parameter
adaptation run simultaneously as shown in 2.The adaptive MPC structure is
formulated as follows with just a minor modification where $\theta$ is
replaced by $\hat{\theta}_{t}$ which is an output of the parameter adaptation
algorithm.

$\displaystyle U^{*}_{t}(x(t))=$
$\displaystyle\underset{U_{t}}{\text{argmin}}\sum_{k=0}^{N-1}q(X_{t+k},u_{t+k})$
(20) subject to $\displaystyle X_{t}=X(t)$ $\displaystyle
X_{t+k+1}=f(X_{t+k},u_{t+k},\hat{\theta}_{t})$ $\displaystyle
u_{t+k}\in\mathcal{U}$ $\displaystyle
U_{t}=\\{u_{t},u_{t+1},\dots,u_{t+N-1}\\}$

Figure 2: Adaptive MPC Algorithm

### Parameter Adaptation Algorithm (PAA):

A simple series parallel parameter adaptation algorithm from existing
literature will be used in tandem with MPC to update the parameters, but any
other PAA could be used depending on the complexity of the problem. The series
parallel structure for PAA is formulated in this subsection. Due to the
transformation shown in (17), dynamical equations containing the system
parameters can be written as an inner product of the system parameters and a
regressor vector containing functions of the states.

$\displaystyle v_{t+1}=\theta_{v}^{T}\phi_{v,t}$ (21)
$\displaystyle\omega_{t+1}=\theta_{\omega}^{T}\phi_{\omega,t}$

where $\theta_{v}=[\alpha_{v},\beta_{v}]^{T}$ and
$\theta_{\omega}=[\alpha_{\omega},\beta_{\omega}]^{T}$ At every instant with
the estimate of the system parameters available, an estimate of the states at
the next instance can be found using

$\displaystyle\hat{v}_{t+1}=\hat{\theta}_{v,t}^{T}\phi_{v,t}$ (22)
$\displaystyle\hat{\omega}_{t+1}=\hat{\theta}_{\omega,t}^{T}\phi_{\omega,t}$

The error between the actual states measured and the projected states at the
next instance can now be used to update the parameters using a series parallel
PAA as

$\displaystyle\hat{\theta}_{v,t+1}=\hat{\theta}_{v,t}+F_{v,t+1}\phi_{v,t}\epsilon_{v,t+1}^{0}$
(23) $\displaystyle
F_{v,t+1}=F_{v,t}-\frac{F_{v,t}\phi_{v,t}\phi_{v,t}^{T}F_{v,t}^{T}}{1+\phi_{v,t}^{T}F_{v,t}\phi_{v,t}}$
$\displaystyle\hat{\theta}_{\omega,t+1}=\hat{\theta}_{\omega,t}+F_{\omega,t+1}\phi_{v,t}\epsilon_{\omega,t+1}^{0}$
$\displaystyle
F_{\omega,t+1}=F_{\omega,t}-\frac{F_{\omega,t}\phi_{\omega,t}\phi_{\omega,t}^{T}F_{\omega,t}^{T}}{1+\phi_{\omega,t}^{T}F_{\omega,t}\phi_{\omega,t}}$

where $\epsilon_{v}^{0}(t+1)=v_{t+1}-\hat{v}_{t+1}$ and
$\epsilon_{\omega}^{0}(t+1)=\omega_{t+1}-\hat{\omega}_{t+1}$ are the
estimation errors and
$\hat{\theta}(t)=[\hat{\theta}_{v}(t)^{T}\hat{\theta}_{\omega}(t)^{T}]^{T}$.
These updated system parameters at every instance can now be used in the
Adaptive MPC structure in (20).

##  IV Numerical Experiments

In this section, numerical results of the algorithm developed in the previous
section tested on a mobile robot with mass m = 5 kg, drag coefficient
corresponding to linear motion b = 0.1 kgs-1, Moment of inertia J = 0.2 kgm2
and drag coefficient corresponding to angular motion c = 0.1 kgm2s-1 are
presented. A quadratic cost with all the matrices P,Q and R as identity
matrices of proper dimensions is considered. A total horizon M = 500 and a
prediction horizon N = 30 is considered to evaluate the algorithm. The robot
starting from the location (1,1) with a heading angle of 0 radians was asked
to reach the origin with a heading of o radians.

$\displaystyle U^{*}_{t}(x(t))=$
$\displaystyle\underset{U_{t}}{\text{argmin}}\frac{1}{2}e_{N}^{T}P_{N}e_{N}+\sum_{k=1}^{N}\frac{1}{2}e_{k}^{T}Qe_{k}+\frac{1}{2}u_{k}^{T}Ru_{k}$
(24) subject to $\displaystyle X_{t}=X(t)$ $\displaystyle
X_{t+k+1}=f(X_{t+k},u_{t+k},\hat{\theta}_{t})$ $\displaystyle
u_{t+k}\in\mathcal{U}$ $\displaystyle
U_{t}=\\{u_{t},u_{t+1},\dots,u_{t+N-1}\\}$

### Performance

As shown in the plots in the Appendix, using the Adaptive MPC controller, the
robot was able to successfully track the desired location without too much
deviation. Testing different system parameters also showed that the robot was
able to adapt and continue to successfully reach the goal.

Figure 3: Final robot path (all positions are in meters)

Figure 4: Optimal states of the robot with Adaptive MPC (All the states are in
SI units and the time is in deci seconds to match the iteration number)

Figure 5: Optimal input values with no knowledge of the system parameters
(Both the inputs are in SI units)

Figure 6: Optimal path of the robot with MPC with full knowledge of the system

Figure 7: Optimal states of the robot with MPC with full knowledge of the
system

Figure 8: Optimal inputs of the robot with MPC with full knowledge of the
system

## Conclusion

In this work, adaptive MPC was shown to be an effective strategy for optimal
control of a simple wheeled robot with unknown parameters. Utilizing a
nonlinear system model and introducing an adaptive element to the MPC to
control for unknown system parameters makes it more difficult to show
persistent feasibility. However, our simulation results look reasonably
accurate and have very interesting applications. The goal tracking problem can
easily be extended for trajectory tracking by modifying the cost function to
minimize the distance from the trajectory and soft constraints can also be
introduced on the states of the system.

## References

  * [1] Francesco Borrelli, Alberto Bemporad, Manfred Morari, ”Predictive control for linear and hybrid systems”, Cambridge University Press.
  * [2] Jason Kong, Mark Pfeiffer, Georg Schildbach, Francesco Borrelli, ”Kinematic and dynamic vehicle models for autonomous driving control design”, in IEEE Intelligent Vehicles Symposium (IV), 2015, pp 1094–1099.
  * [3] Veronica Adetola, Darryl DeHaan, Martin Guay, ”Adaptive model predictive control for constrained nonlinear systems”, in Systems and Control Letters, Volume 58, Issue 5, 2009, Pages 320–326.
  * [4] Hiroaki Fukushima, Tae-Hyoung Kim, Toshiharu Sugie,”Adaptive model predictive control for a class of constrained linear systems based on the comparison model”, in Automatica, Volume 43, Issue 2, 2007, Pages 301–308.
  * [5] Tae-Hyoung Kim, H. Fukushima and T. Sugie, ”Robust adaptive model predictive control based on comparison model,” 2004 43rd IEEE Conference on Decision and Control (CDC) (IEEE Cat. No.04CH37601), 2004, pp. 2041-2046 Vol.2, doi: 10.1109/CDC.2004.1430348.
  * [6] K. Pereida and A. P. Schoellig, ”Adaptive Model Predictive Control for High-Accuracy Trajectory Tracking in Changing Conditions,” 2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2018, pp. 7831-7837, doi: 10.1109/IROS.2018.8594267.
  * [7] Anthony M Bloch, ”Nonholonomic mechanics and control”, Springer, New York, NY, 2003.
  * [8] AM Bloch, PS Krishnaprasad, JE Marsden, RM Murray, ”Nonholonomic mechanical systems with symmetry”, in Archive for Rational Mechanics and Analysis, Vol 136, Issue 1, 1996, pp 21–99.
  * [9] Lennart Ljung, Keith Glover, ”Frequency domain versus time domain methods in system identification”, in Automatica, Vol. 17, Issue 1, 1981, pp. 71-86.
  * [10] MATLAB. version 9.10.0 (R2016b). Natick, Massachusetts: The MathWorks Inc.,2016.
  * [11] Zhang, Qinghua. “Using Wavelet Network in Nonparametric Estimation”, IEEE Transactions on Neural Networks 8, no. 2 , March 1997. https://doi.org/10.1109/72.557660.
  * [12] Sjöberg, Jonas, Qinghua Zhang, Lennart Ljung, Albert Benveniste, Bernard Delyon, Pierre-Yves Glorennec, Håkan Hjalmarsson, and Anatoli Juditsky. “Nonlinear Black-Box Modeling in System Identification: A Unified Overview”, Automatica 31, no. 12 (December 1995): 1691–-1724. https://doi.org/10.1016/0005-1098(95)00120-8.
  * [13] Ljung, Lennart, and Torkel Glad. Modeling of Dynamic Systems. Prentice Hall Information and System Sciences Series. Englewood Cliffs, NJ: PTR Prentice Hall, 1994.
  * [14] Ljung, Lennart. System Identification: Theory for the User. Second edition. Prentice Hall Information and System Sciences Series. Upper Saddle River, NJ: PTR Prentice Hall, 1999.
  * [15] Söderström, Torsten, and Petre Stoica. System Identification. Prentice Hall International Series in Systems and Control Engineering. New York: Prentice Hall, 1989.
  * [16] Lennart Ljung, ”Analysis of recursive stochastic algorithms”, IEEE transactions on automatic control 22.4, 1977, pp. 551–575.