Additional authors: Rafael S. Pereira(LNCC,danielnascramos@gmail.com), Eduardo
Ogasawara (CEFET-RJ, eogasawara@ieee.org), Flavia Delicato (UFF,
fdelicato@gmail.com),Paulo Pires (UFF, paulo.f.pires@gmail.com), and Artur
Ziviani (LNCC<EMAIL_ADDRESS>
# Hyperspherical embedding for novel class classification

Rafael S. Pereira

Alexis Joly

Patrick Valduriez

Fabio Porto

DEXL-LNCC  Petropolis, Brazil<EMAIL_ADDRESS>Zenith-INRIA  Montpellier,
France<EMAIL_ADDRESS>Zenith-INRIA  Montpellier, France
<EMAIL_ADDRESS>DEXL-LNCC  Petropolis, Brazil<EMAIL_ADDRESS>
(12 December 2019)

##  1 abstract

Deep learning models have become increasingly useful in many different
industries. On the domain of image classification, convolutional neural
networks proved the ability to learn robust features for the closed set
problem, as shown in many different datasets, such as: MNIST, [9], FASHION
MNIST [20], CIFAR10 and CIFAR 100 [8], and IMAGENET [2]. These approaches use
deep neural networks with dense layers with softmax activation functions in
order to learn features that can separate classes in a latent space. However,
this traditional approach is not useful for identifying classes unseen on the
training set, known as the open set problem[21]. A similar problem occurs in
scenarios involving learning on small data. To tackle both problems, few shot
learning[19] has been proposed. In particular, metric learning learns features
that obey constraints of a metric distance in the latent space in order to
perform classification. However, while this approach proves to be useful for
the open set problem, the state-of-the-art implementation requires pairwise
training, where both positive and negative examples of similar images are
presented during the training phase, which limits the applicability of these
approaches in large data or large class scenarios given the combinatorial
nature of the possible inputs [17]. In this paper, we present a constraint
based approach applied to the representations in the latent space under the
normalized softmax loss, proposed by [18]. We experimentally validate the
proposed approach for the classification of unseen classes on different
datasets using both metric learning and the normalized softmax loss, on
disjoint and joint scenarios. Our results show that not only our proposed
strategy can be efficiently trained on larger set of classes, as it does not
require pairwise learning, but also present better classification results than
the metric learning strategies surpassing its accuracy by a significant
margin.

##  2 Introduction

Humans have the ability to identify many different types of objects [3]. Even
when we are not able to name a certain object, we can tell its differences
from a second object, which contributes to the identification of objects we
have never seen before and group them into classes based on prior knowledge.
Metric learning [7] is a well adopted approach that identifies novel classes
without fine tuning a model on these classes. The approach applies an
optimization strategy, which guarantees that the classes a model has seen
during optimization form disjoint clusters on the latent space according to a
certain metric distance. Some common approaches that use this strategy are:
the triplet loss [11]; constrative loss [5]; prototypical networks [13];
constellation loss [10]; and matching networks [16]. Another approach in
metric learning is called similarity learning, where the model receives pairs
of inputs and learns that they are similar if they belong to the same class or
dissimilar otherwise, as discussed in [14]. During inference on novel classes,
the former uses the distance between labeled points of the novel class the
model was not optimized upon to obtain a representation in the latent space
for the novel class and then calculate the distance between new points and
each class representation. In the latter approach, a similarity score is
calculated between every (class,query) pair in order to find the most similar
pair.

However while enforcing metric properties on the latent space leverages the
model knowledge to novel classes, it requires pairwise learning, which limits
the scalability of such approaches given the amount of possible pairs.

In this paper we take into account the normalized softmax loss(NSL), proposed
by [18], and present how it enforces a latent space that obeys the cosine
similarity. Based on this, we then present a methodology to apply the NSL to
the novel classes classification problem. Considering a trained artificial
neural network, we add a new neuron to its last layer and infer the weights
that connect the penultimate layer of the network to this neuron. The
connection and the new neuron are used to classify a novel class by using few
labeled samples of it. Our approach for the open set problem allows us to
classify new classes without fine-tuning the model , instead we use the same
network parameters the model was optimized upon to classify its seen classes
and only adding a new neuron along with its inferred connection. We evaluate
state-of-the-art approaches to solve the open set problem against our proposed
approach, both in the disjoint and joint scenarios, as defined in section 4.4,
for different datasets. The experimental results show that our approach
outperforms other metric learning strategies and it, additionally, induces a
more scalable training process, as it does not require pairwise learning,
leveraging the open set problem technique to deal with large datasets.

The remainder of this paper is structured as follows. First, it presents some
theoretical background at section 4. Our methodology and how to classify new
classes is described in section 5. Next, we present the results on the joint
and disjoint open set problem in section 6. Moreover, we present the use of
the NSL approach in a more complex dataset in the field of botany, in section
7. Finally, we conclude in section 8.

##  3 Problem Formalization

Let there be a dataset $D$ composed of a set of examples $x\in X$, such that
$X\sqsubseteq\mathbb{R}^{N}$, an $N$ is the number of features describing each
instance $x$. A classification task aims to build a function that takes the
features belonging to each example and project the object into a space
$\mathbb{R}^{k}$, where $k$ is the number of classes in $Y$.
We want to build a function $Y=\phi(X)$,$Y\sqsubseteq\mathbb{R}^{M}$, modelled
by a neural network, where features $F$ characterizing examples in $X$ sharing
the same class $y$ are clustered together according to a metric distance in
$Y$. We aim to use this metric property for the identification of novel
classes. Furthermore, we consider the mean point, in a given metric space, of
a set of points in a class to be referred as the class prototype [13], which
represents a class virtual centroid. To model the function as a neural network
and enforce the metric properties on $Y$, we define the network with two
modules. The first is a embedding module that takes $F$ as an input and
returns a $F^{\prime}$ representation in $Y$. This representation is fed to a
classification module composed of $k$ neurons, where $k$ is the total number
of classes $Y$. The modules are connected by a matrix of learnable weights $W$
of size (Mxk). To enforce the metric space, we constrain the optimization
procedure to guarantee that a vector of features extracted from $x$ share high
cosine similarity with a vector $w$ of weights for the class $y$ to which $x$
is classified in.

##  4 Preliminaries

###  4.1 The Open set problem

The classification problem can be formulated as a closed set or open set
problem. In the closed set problem context, the optimization process trains a
model to learn features that can identify the classes with samples present in
the training set. The approach does not require the identification of classes
not present in the training set. This is commonly tackled using the Softmax-
Crossentropy loss [6][12],[15]. In contrast, in the open set problem we are
interested in not only identifying the classes present in the training set,
but also to be able to use the model to classify new classes by exploiting
properties in the latent space yielded during optimization.

###  4.2 Metric Learning and the latent space

Deep neural networks learn a set of transformations and relations among the
inputs in order to obtain the desired output during optimization. The network
can be broken into two parts: the projection module, which takes the data and
transforms it into a representation in the latent space; and the data
representation processing, applied by the layer implementing the desired task.
The latter enforces the latent space to have some properties that are defined
by the task. For example, the regression task is defined as
$y=\sum_{i}w_{i}x_{i}+b$ and so given a fixed $w$ if $x\approx x^{\prime}$
then $y\approx y^{\prime}$. When discussing the classification problem, the
softmax crossentropy approach enforces the inequality
$w_{i}x_{i}+b_{i}>>w_{j}x_{i}+b_{j}$, where $i,j$ represent different classes
[18]. A proposed alternative that optimizes the latent space directly and can
enforce metric properties that allows the model to be used for novel classes
is know as metric learning.

The metric learning approach learns a set of features that obey a metric
distance on the latent space. The model can be optimized to learn a similarity
metric between pairs, as proposed in [14], or can enforce the latent space to
obey a predefined metric distance like euclidean distance or cosine
similarity. Some strategies, such as the constrative loss [5],learn on pairs
of data, while others learn using triplets like triplet loss [11].

Optimization on these approaches aim to obtain disjoint clusters for each
class of interest in the latent space, according to a predefined metric
distance. As a consequence of the approach, classification can be performed
for novel classes by using the representation of an anchor example and
calculating the metric distance between a query point and the anchor.

###  4.3 Normalized Softmax Loss properties

Proposed in [18], the NSL (Normalized Softmax Loss) is a modification of the
softmax loss that enforces a cosine similarity metric between classes on the
latent space. It enforces the features that are projected into the latent
space to be contained in a hypersphere ($N>3)$ where each region of the sphere
contains features belonging to a certain class. To discuss it, we first
present the usual loss function adopted in classification tasks called softmax
loss, which corresponds to the use of the softmax activation function along
with crossentropy, in Equation 1. The latter presents the value of the
activation to a class $i$ in a classification problem of $N$ classes.

$G_{i}(x)=\frac{e^{w_{i}x+b_{i}}}{\sum_{j=1}^{j^{N}}e^{w_{j}x+b_{j}}}$ (1)

Since $0=<G_{i}(x)<=1$, the optimization on this function aims to obtain $w$
and $b$ so that $w_{i}x_{i}+b_{i}>>w_{j}x_{i}+b_{j}$ for all $i\neq j$ in
order to maximize $G_{i}(x_{i})$ while minimizing $G_{i}(x_{j})$. This
approach has been extensively used in the context of the closed set problem
[15][12][6], but lacks generalization capability for the open set problem.

Given this, the NSL contribution lies on constraining the softmax function.
Those constraints come in the form of defining $|b|=0$, $|w|=1$ and $|x|=S$,
$S\in\mathbb{R}$. Those constraints enforce that the weights $|w|$ must be
contained in the surface of a hypersphere of radius one, while $x$ is
contained in the surface of a second hypersphere of radius $S$, both
hyperspheres are of dimension $N$ equal to the amount of neurons in the
penultimate layer.

Equation 2 presents what those constraints imply in the optimization process
and the nature of its extracted features. A visual example of the features and
weights of the NSL can be seen in figures 2 and 3.

$\displaystyle w_{i}x_{i}+b_{i}>>w_{j}x_{i}+b_{j}$ $\displaystyle
w_{i}x_{i}+0>>w_{j}x_{i}+0$
$\displaystyle|w_{i}||x_{i}|cos(\theta(w_{i},x_{i}))>>|w_{j}||x_{i}|cos(\theta(w_{j},x_{i}))$
(2) $\displaystyle Scos(\theta(w_{i},x_{i}))>>Scos(\theta(w_{j},x_{i}))$

Since $S$ is a hyperparameter, optimization happens on $\theta$. On its turn,
$cos(\theta)$ is a function ranging in the interval [-1,1]. In this context,
to maximize the left side of the equation (2) the angle between $x_{i}$ and
$w_{i}$ must be close to $0$, while to minimize the right side of the equation
(2), $\theta(w_{j},x_{i})$ must become closer to $\pi$. Those values could be
achieved if the loss value reached zero.

###  4.4 Disjoint and Joint scenarios

In this section, we briefly define the joint and disjoint scenarios, as well
as the seen and unseen classes for both scenarios.
First, consider the dataset $D$ with a set of classes $C$. From $C$ we
consider a subset of classes $C^{\prime}$ that a model will see during its
training phase, while there is a complement subset $C*$ the model did not see,
so that $C=C^{\prime}\cup C*$.

In the disjoint scenario, a neural network is optimized on $C^{\prime}$. After
optimization, we sample a subset of $C*$ and use the knowledge obtained from
$C^{\prime}$ to classify only examples from this subset. The seen classes are
the classes contained in $C^{\prime}$, while the unseen classes are the subset
of $C*$. Metrics are calculated on the test set and the method to perform this
transfer of knowledge are specific to each method. The ones we explore in this
work have their methods discussed in sections 5.1 and 5.2.

The joint scenario follows mostly the same proposed methodology as the
disjoint scenario. The major difference is that during inference we do not
discard the classes used during optimization, instead the model performs
inference between the union of $C^{\prime}$ and the subset of $C*$ as possible
answers. This scenario considers we wish to extend a model that is being used
to novel classes without retraining. The disjoint scenario considers that your
problem of interest may not have much data, and its preferable to optimize on
a related dataset with larger amounts of data.

##  5 Methodology

In this paper we aim to compare pairwise strategies, commonly used in metric
learning, against the normalized softmax loss approach for the open set
problem. In this manner we consider both the problem where during inference
seen and unseen classes are disjoint, as well as the scenario where the model
must identify both the seen and unseen classes together. To guarantee that the
extracted features for unseen classes are robust, the model does not consider
them to be the same as the seen classes.

In the first scenario, we optimize a model using $N$ classes on the training
set. After optimization occurs for a number of epochs (eg. 5) we evaluate the
model accuracy on the $N$ seen classes and the $M$ unseen classes on the test
set by using its training set to infer the weights for the novel classes.

The approach that implements the metric learning strategies is discussed in
section 5.1, while the one applying Normalized softmax loss is presented in
section 5.2

###  5.1 Classifying new classes

When tackling the open set problem, we are interested in optimizing models in
which the full knowledge the network obtains during optimization can be
exploited for classes outside of the training set. The usual softmax
crossentropy approach lacks the ability to extract features that obey this
property, as the weights between the penultimate layer and the classification
layer $w$ are as important as the representation in the latent space of the
penultimate layer $x$ as seen in equation 1, and the former is undefined for
novel classes. Usual approaches for classifying novel classes are explored in
metric learning [5],[11],[14],[16] where model optimization occurs by learning
a metric space that creates clusters based on a metric distance [5],[11],[16],
or by learning whether two examples are similar or dissimilar [14]. These
strategies allow classification to be performed for classes outside of the
training set as the representation in the latent space can be exploited via a
K nearest neighbours approach. However, all these strategies are based on
pairs or triplets, which limits the approaches scalability whenever optimizing
with larger amounts of data per class, or a large amount of classes in a small
data scenario.

In our work, to classify new classes, we will be adopting the proposal of [13]
that builds a single representation on the latent space for each class and
then builds a K nearest neighbour (KNN) model with the same metric distance
the model was optimized on. To perform predictions on the test set we extract
the sample representation in the latent space and feed it to the KNN model.

###  5.2 Classifying new classes via NSL

As already discussed in section 4.3, the NSL learns a metric space that
enforces cosine similarity. To minimize its loss function it must satisfy
Equation 2. In order to maximize the left side of equation 2, equation 3 must
be satisfied, whereas to minimize the right side of equation 2, equation 4
must be satisfied for all $i,j$ class pairs.

Given this, in the limit, when the loss value goes to zero, both equations: 3
and 4 must hold true. In the specific case where the angle between $w_{i}$ and
$x_{i}$ is close to zero then equation 5 is satisfied since
$\frac{x_{im}}{|x_{im}|}\approx w_{i}$ when equation 3 is satisfied where
$x_{im}$ is the feature vector of a example $m$ belonging to class $i$ and
$w_{i}$ is the weight vector that connects the penultimate layer to the neuron
associated with class $i$.

Given this relation, we propose the Equation 6 as a form of inferring the
weights that connect the penultimate layer of the neural network to a neuron
representing a novel class, by using labeled samples from this novel class in
a manner similar to [13]. A geometrical representation that shows the
relationship between these weights and the feature vectors of the latent space
is shown in Figures 2 and 3. One can see that even when the loss does not
reach zero, the center of the feature vector is aligned with its corresponding
class weights. We also show that Equation 6 has a small weight inference error
on non zero losses. In order to do so, we optimized different models on
subsets of the Cifar100 dataset. Then we calculated the cosine distance matrix
between the original class weights and the one from the inferred weights and
present the relative error of the reconstruction in Figure 4. In the discussed
equations, we use the notation $x_{im}$ as the example $m$ of class $i$.

$\frac{w_{i}\cdot x_{im}}{|w_{i}||x_{im}|}\geq 1$ (3) $\frac{w_{i}\cdot
x_{jn}}{|w_{i}||x_{jn}|}\geq-1$ (4)

$\frac{x_{im}\cdot x_{in}}{|x_{in}||x_{im}|}\geq 1$ (5)

$w_{i}=\frac{1}{N}\sum_{i=1}^{N}\frac{x_{i}|w_{i}|}{|x_{i}|}$ (6)

A dataflow depicting our approach to infer the weights for novel classes is
presented in Figure 1.

Figure 1: Diagram presenting the approach to infer weights for the decision
layer for new classes Figure 2: Embedding obtained on the cifar10 dataset when
using a latent space with two dimensions, each color represents a different
class. Inner points are the class weights while outer points come from the
training set, notice how classes from the outer circle are aligned with the
inner circle Figure 3: Embedding obtained on the Plantnet dataset discussed in
section 7 when using a latent space with two dimensions Figure 4: Relative
error for rebuilding the optimized weights by using equation 6 on the cifar
100 dataset when optimizing on different amount of classes on 30 different
runs

##  6 Results

###  6.1 Experimental setup

In this section we present the experimental setup. All experiments presented
for the FASHION MNIST and CIFAR datasets were performed using google
collaboratory and its GPU instances. Experiments using the Pl@ntNet dataset
were performed using a Dell PowerEdge R730 server, with 2 CPUs Intel (R) Xeon
(R) CPU E5-2690 v3 @ 2.60GHz; 768 GB of RAM; and running on a Linux CentOS
7.7.1908 kernel version 3.10.0-1062.4.3.e17.x86_64. The machine is equipped
with a single NVIDIA Pascal P100 GPU, with 16GB RAM. Implementations were
performed using Python 3.7 along with the Keras deep learning library.

###  6.2 How does the amount of samples affect the class prototype ?

A common scenario in which the identification of unseen classes appears is the
one where the amount of available data samples for classes of interest is
small or the cost of optimizing another model to include the new classes
becomes too costly. Therefore, all strategies discussed in this paper classify
new classes based on labeled examples without retraining. The influence of the
number of samples needed to perform classification tasks using the NSL
approach is shown in three different datasets in Figure 5 for a model
optimized for 30 epochs. We define class prototype in the same way as [13].

Figure 5: F1 score on the test set for three different datasets where we show
results when weights are: (a) Trained: the ones found during optimization; (b)
Single anchor: weights are inferred using a single random anchor example; and
(c) weights are inferred using the training set to build a class prototype.

As it can be seen on Figure 5, weights inferred using Equation 6 on the
training set (c) maintain the same model accuracy as by using the weights
obtained during optimization. We also can observe that model quality, when
inferring via a single example, decays in relation to the task complexity. On
the x axis we present how the class weights were obtained where class
prototype uses the whole training set to infer the weights according to our
methodology, while single anchor uses a single random example from the
training set for weight inference. Y axis presents the F1-Score on the test
set.

###  6.3 NSL architectural search on the closed set

In this section, we present in Table 1 an ablation study comparing NSL against
softmax loss applied to the same network architecture. We define the following
deep neural networks: an Inception block, as a single Inception module as
proposed in [15]; a VGG block is a composition of two Convolutional layers
with 3x3 kernel and a 2x2 MaxPooling layer; while a Resnet block is built with
a 3x3 convolution with a residual connection. For the number of filters we
increase them in the order $<32,64,128>$ values.

Architecture | Number of blocks | SL | NSL
---|---|---|---
Inception | 1 | 0.473 | 0.633
Inception | 2 | 0.47 | 0.746
Inception | 3 | 0.42 | 0.803
ResNet | 1 | 0.1 | 0.610
ResNet | 2 | 0.1 | 0.75
ResNet | 3 | 0.1 | 0.81
VGG | 1 | 0.175 | 0.638
VGG | 2 | 0.187 | 0.738
VGG | 3 | 0.643 | 0.823
Table 1: Analysis on the closed set for different architectural setups
optimized on softmax vs normalized softmax for the cifar10 dataset. SL and NSL
present model accuracy on the test set.

In table 1 we compared the normalized softmax loss vs softmax loss on
different architectures. We note that for all architecture variations as well
as variations in the network depth, the normalized softmax loss outperformed
the softmax loss. We also denote that on softmax loss the choice of base
architecture (Inception,Resnet,VGG) significantly influenced the model
accuracy as much as depth. In contrast, the normalized softmax loss was more
influenced by model depth in this dataset.

###  6.4 Evaluating the disjoint scenario

In this section we show the results of evaluating the model accuracy on the
test set from seen and unseen classes by employing a VGG based model with two
blocks.

To perform the experiment in a dataset with $K$ classes, we optimize the model
with $K-N$ classes, and use the network structure to classify the $N$ unseen
classes without considering the seen ones as possible answers. The approach to
do so was discussed in section 5. Results are presented in Tables 2 and 3.

N | NSL | Triplet | Constrative
---|---|---|---
2 | 0.852 | 0.663 | 0.700
3 | 0.725 | 0.570 | 0.551
4 | 0.629 | 0.406 | 0.422
5 | 0.545 | 0.296 | 0.328
Table 2: Model results for the cifar10 dataset. $N$ refers to the amount of unseen classes. NSL shows model accuracy on unseen classes when a model is optimized on $10-N$ classes and only identifies $N$. Triplet and constrative refers to the same problem using triplet and constrative loss, respectively. N | NSL | Triplet | Constrative
---|---|---|---
2 | 0.897 | 0.62 | 0.703
3 | 0.876 | 0.39 | 0.469
4 | 0.841 | 0.27 | 0.312
5 | 0.807 | 0.22 | 0.2
Table 3: Model results for the fashion MNIST dataset. $N$ refers to the amount
of unseen classes. The NSL column shows model accuracy on unseen classes when
a model is optimized on $10-N$ classes and only identifies $N$. Triplet and
constrative refer to the same problem using triplet and constrative loss,
respectively.

In t]Tables 2 and 3, we compared NSL against two metric learning strategies in
a disjoint settings. In the first line, we present a scenario we train with 8
random classes and evaluate on the other two. The second line trained with 7
and so forth. Our results show that in both datasets the NSL outperformed
those metric learning strategies for evaluating novel classes in a disjoint
scenario.

###  6.5 Evaluating the joint scenario

In this section we present results when the novel classes must be integrated
into the classification process along with the classes used for optimization.
The model is optimized with $10-N$ classes and we evaluate its accuracy on
these and on the $N$ unseen classes, in a 10 possible classes scenario.
Results are presented in Tables 4 and 5

N | Type | NSL | Triplet | Constrative
---|---|---|---|---
2 | Seen | 0.559 | 0.195 | 0.205
3 | Seen | 0.578 | 0.128 | 0.182
4 | Seen | 0.620 | 0.08 | 0.166
5 | Seen | 0.641 | 0.05 | 0.215
2 | Unseen | 0.501 | 0.226 | 0.189
3 | Unseen | 0.433 | 0.156 | 0.176
4 | Unseen | 0.391 | 0.127 | 0.160
5 | Unseen | 0.357 | 0.146 | 0.176
Table 4: Model results for the cifar10 dataset, N refers to the amount of unseen classes, Seen shows model accuracy for the seen classes after optimization for 5 epochs and Unseen classes shows model performance to identify the new classes using NSL N | Type | NSL | Triplet | Constrative
---|---|---|---|---
2 | Seen | 0.854 | 0.708 | 0.657
3 | Seen | 0.842 | 0.705 | 0.440
4 | Seen | 0.860 | 0.688 | 0.284
5 | Seen | 0.878 | 0.690 | 0.19
2 | Unseen | 0.773 | 0.724 | 0.613
3 | Unseen | 0.771 | 0.716 | 0.390
4 | Unseen | 0.739 | 0.680 | 0.286
5 | Unseen | 0.719 | 0.691 | 0.181
Table 5: Model results for the fashion MNIST dataset. $N$ refers to the amount
of unseen classes. Seen shows model accuracy for the seen classes after
optimization for 5 epochs and Unseen classes shows model performance to
identify the new classes using NSL.

Tables 4 and 5 depicts the results of comparing our approach using $NSL$ with
metric learning strategies: Triplet loss and constrative loss on both cifar
and MNIST datasets, considering a joint scenario. The NSL outperformed both
approaches on the two evaluated datasets, including seen and unseen classes
predictions.

##  7 Case Study: The Pl@ntnet dataset

In order to assess on real world data, we evaluate the approach for the closed
and open set problems on the Plantnet dataset. Plantnet [1] is a mobile
application which allows users to take pictures of different plant species,
and by the use of deep learning techniques these pictures are processed
returning the most probable species each plant belongs to. The dataset becomes
challenging as available pictures have different levels of quality, as well as
having multiple species from many different parts of the world as shown in
figure 6. Given this we wish to evaluate how our approach fares in the closed
set compared to traditional crossentropy approaches, as well as evaluate it as
an open set problem. The problem becomes relevant to the evaluation of the
proposed approach given that the scenario is usually represented by a long
tail distribution, in which some classes are very common, while others are
rare and lack significant available training data.

The subset of the Plantnet data we used was obtained from [4] and it has a
total of 182 classes.

Figure 6: Distribution of the training dataset. Note how the data is highly
imbalanced as some species are very common while others are rare.

###  7.1 Experimental design

As it is clear from Figure 6, there is a high imbalance among classes in the
Plantnet dataset. Also, there are many classes in the training set which have
very small amounts of data. Since many plant species have few samples, we are
interested in exploring the performance of NSL, where a model is optimized
only on more common species and weights to classify uncommon species are
inferred, as discussed in section 5.2. To this end, we perform experiments by
optimizing the model only on classes where the number of samples is larger or
equal to $N=\\{200,100,50,25\\}$, which results in $C=\\{10,16,28,43\\}$,
where $C$ is the number of classes, and present the results for joint and
disjoint settings. Unseen classes will be selected randomly among those with
number of samples $M$, $M<N$, and results will be presented with the average
of 30 runs. All models are optimized on 100 epochs and weights that minimize
validation loss are used for inference. Model architecture is presented in
Figure 7. Pre-processing steps on the data only include normalizing it to
$[0,1]$ range and reshaping it to a $<96,96,3>$ shape. Models are compiled
with weighted crossentropy to take the class imbalance into account on the
loss functions, ensuring that solutions that output only the majority class
are penalized. For the open set tasks, we report balanced accuracy to better
take into account class imbalance.

Figure 7: Network architecture used to perform the Plantnet experiments

###  7.2 Results

In this section we present results on the model accuracy for both the disjoint
and joint open set problems when the model was optimized on different amount
of classes. As already discussed, we selected seen classes based on a filter
on the number of samples $\geq N$. Given this, the relation between the number
of samples and the number of classes is as follows:
$[25,43],[50,28],[100,16],[200,10]$.

Number of classes | NSL accuracy | SL accuracy
---|---|---
10 | 0.7349 | 0.1312
16 | 0.6617 | 0.0565
28 | 0.5382 | 0.0228
43 | 0.3974 | 0.0038
Table 6: Model accuracy on the test set optimized for 100 epochs on weighted
crossentropy

In table 6 we presented the model accuracy for the four different models that
will be used in the plantnet analysis. All models were optimized by receiving
the same amount of samples per epoch as well as number of epochs. We shown how
in the base configuration while our network mostly remained in the NULL
hypothesis using softmax loss it could actually decrease its loss using NSL.

####  7.2.1 Disjoint scenario

In this subsection we present the disjoint analysis for the Plantnet dataset.
We instantiate our base model without the last layer and then perform
classification on novel classes randomly sampled from the total unseen classes
by inferring the weights between the penultimate layer and the decision layer.
We present balanced accuracy on the model capacity on novel classes by
performing 30 runs for each value of $N$. Results are presented in figure 8.

Figure 8: Comparing models with the same architecture and optimized upon the
same amount of data and number of epochs, but with different amount of seen
classes on their ability to classify novel classes

In figure 8 we show a comparison between models trained on different amount of
classes with the same amount of data by presenting their ability to identify
novel classes in a disjointed scenario. Our results show that diversity of
classes seen during training allowed the model to become more robust for novel
classes as the model trained with 10 classes performed the worst of all for
novel classes. However its also important to note that optimizing on higher
amount of classes is a more complex problem, requiring more data to be seen by
the model, more updates or a more complex model to learn robust features for
all seen classes. This is shown on the curve obtained from the model optimized
upon 43 classes, which ranks the second worst. Our best result was seen on the
model with 28 classes, which had a high class diversity while also learning
robust features during training. .

####  7.2.2 Joint scenario

In this subsection, we present the joint analysis for the Plantnet dataset. We
present results where we instantiate the base model trained on $N$ classes.
Then, we add other $M$ unseen classes so that the model must classify between
$N+M$ classes. Weights for the $M$ classes are inferred as described in
section 5.2 and we report model accuracy for the overall model, as well as for
the $M$ classes. Results are presented in figures 9 and 10.

Figure 9: Analysis for the joint scenario showing the results for different
models on overall class architecture  Figure 10: Analysis for the joint
scenario showing the results for different models on the unseen classes

In figures 9 and 10 we presented our results for a jointed scenario evaluating
the overall model quality as new classes are added as well as showing the
balanced accuracy calculated only on the unseen classes. Our conclusions on
the disjointed scenario seen in figure 8 also hold in these results, as we can
see for both of these cases the model trained on 28 classes has the higher
balanced accuracy given the same number of total classes and the model
optimized on 43 classes shows the worse balanced accuracy. The lack of
diversity of the model optimized on 10 classes can be seen influencing its
quality as the number of novel classes increase in both scenarios.

##  8 Conclusion

In this paper, we presented how the normalized softmax loss can be employed on
the open set problem. We presented results on different datasets for both the
disjoint and joint open set problems and compared them to metric learning
strategies. We show that the NSL based approach demonstrates superior results
producing more robust features and implementing a less costly optimization
procedure, as it does not require pairwise training. Results on real world use
case evaluating on a subset of the Plantnet data shows how our approach can be
employed to identify classes unseen during optimization, with weights
associated to the classification of new data inferred by the approach.

##  9 Acknowledgments

The authors would like to thank Petrobras for supporting this work through the
project ”Development of an Intelligent software platform”. We would also like
to thank the INRIA-Brazil Associated Team cooperation project HPDaSc.

## References

  * [1] A. Affouard, H. Goëau, P. Bonnet, J.-C. Lombardo, and A. Joly.  Pl@ntnet app in the era of deep learning.  In ICLR: International Conference on Learning Representations, 2017\.
  * [2] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei.  ImageNet: A Large-Scale Hierarchical Image Database.  In Conference on Computer Vision and Pattern Recognition, CVPR09, 2009.
  * [3] C. Fields.  Editorial: How humans recognize objects: Segmentation, categorization and individual identification.  Frontiers in Psychology, 7:400, 2016.
  * [4] C. Garcin.  Projects Â· garcin camille / plantnet_dataset.
  * [5] R. Hadsell, S. Chopra, and Y. LeCun.  Dimensionality reduction by learning an invariant mapping.  In Proceedings of the 2006 IEEE Computer Society Conference on Computer Vision and Pattern Recognition - Volume 2, CVPR 2006, page 1735-1742, USA, 2006. IEEE Computer Society.
  * [6] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 770–778, 2016.
  * [7] M. KAYA and H. Å. BÄ°LGE.  Deep metric learning: A survey.  Symmetry, 11(9), 2019.
  * [8] A. Krizhevsky.  Learning multiple layers of features from tiny images.  Technical report, 2009.
  * [9] Y. LeCun, C. Cortes, and C. Burges.  Mnist handwritten digit database.  ATT Labs [Online]. Available: http://yann.lecun.com/exdb/mnist, 2, 2010.
  * [10] A. Medela and A. Picon.  Constellation loss: Improving the efficiency of deep metric learning loss functions for the optimal embedding of histopathological images.  Journal of Pathology Informatics, 11(1):38, 2020.
  * [11] F. Schroff, D. Kalenichenko, and J. Philbin.  Facenet: A unified embedding for face recognition and clustering.  In 2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 815–823, 2015.
  * [12] K. Simonyan and A. Zisserman.  Very deep convolutional networks for large-scale image recognition.  In Y. Bengio and Y. LeCun, editors, 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, 2015.
  * [13] J. Snell, K. Swersky, and R. Zemel.  Prototypical networks for few-shot learning.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30, pages 4077–4087. Curran Associates, Inc., 2017\.
  * [14] F. Sung, Y. Yang, L. Zhang, T. Xiang, P. H. S. Torr, and T. M. Hospedales.  Learning to compare: Relation network for few-shot learning.  In 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 1199–1208, 2018.
  * [15] C. Szegedy, Wei Liu, Yangqing Jia, P. Sermanet, S. Reed, D. Anguelov, D. Erhan, V. Vanhoucke, and A. Rabinovich.  Going deeper with convolutions.  In 2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 1–9, 2015.
  * [16] O. Vinyals, C. Blundell, T. Lillicrap, k. kavukcuoglu, and D. Wierstra.  Matching networks for one shot learning.  In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 29, pages 3630–3638. Curran Associates, Inc., 2016.
  * [17] F. Wang, X. Xiang, J. Cheng, and A. L. Yuille.  Normface: L2 hypersphere embedding for face verification.  In Proceedings of the 25th ACM International Conference on Multimedia, MM ’17, page 1041-1049, New York, NY, USA, 2017. Association for Computing Machinery.
  * [18] H. Wang, Y. Wang, Z. Zhou, X. Ji, D. Gong, J. Zhou, Z. Li, and W. Liu.  Cosface: Large margin cosine loss for deep face recognition.  In 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 5265–5274, 2018.
  * [19] Y. Wang and Q. Yao.  Few-shot learning: A survey.  CoRR, abs/1904.05046, 2019.
  * [20] H. Xiao, K. Rasul, and R. Vollgraf.  Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms.  CoRR, abs/1708.07747, 2017.
  * [21] R. Yoshihashi, W. Shao, R. Kawakami, S. You, M. Iida, and T. Naemura.  Classification-reconstruction learning for open-set recognition.  In 2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 4011–4020, 2019.


# Stochastic entropy production for dynamical systems with restricted
diffusion

Jonathan Dexter and Ian J. Ford  Department of Physics and Astronomy,
University College London, Gower Street, London WC1E 6BT, U.K.

###### Abstract

Modelling the evolution of a system using stochastic dynamics typically
implies a greater subjective uncertainty in the adopted system coordinates as
time progresses, and stochastic entropy production has been developed as a
measure of this change. In some situations the evolution of stochastic entropy
production can be described using an Itô process, but mathematical
difficulties can emerge if diffusion in the system phase space is restricted
to a subspace of lower dimension. This can arise if there are constants of the
motion, for example, or more generally when there are functions of the
coordinates that evolve without noise. We discuss such a case for an open
three-level quantum system modelled within a framework of Markovian quantum
state diffusion and show how the problem of computing the stochastic entropy
production in such a situation can be overcome. We go on to illustrate how a
nonequilibrium stationary state of the three-level system, with a constant
mean production rate of stochastic entropy, can be established under suitable
environmental couplings.

##  I Introduction

Entropy quantifies subjective uncertainty in the configuration of a system
when limited information is available. If the evolution of such a system is
modelled using stochastic dynamics, representing the effects of coupling to an
underspecified environment, then subjective uncertainty in the configuration
of the world (the system plus its environment) increases with time,
corresponding to a growth in the total entropy. The second law of
thermodynamics arises when the world is characterised or perceived at a coarse
grained level, but governed by underlying equations of motion with a
sufficient degree of deterministic chaos [1].

For system variables that evolve continuously according to a set of Markovian
stochastic differential equations (SDEs), or Itô processes, it may be shown
that the associated stochastic entropy production can also be described using
an SDE, but only if the noise terms satisfy certain requirements [2, 3, 4].
Problems arise from constraints imposed on the dynamics, for example the
existence of constants of the motion. In such cases, the matrix describing the
diffusion of system coordinates in their phase space becomes singular. There
are directions in the space in which diffusion does not take place and
mathematical difficulties in the evaluation of the stochastic entropy
production arise as a consequence.

The central aim of this paper is to show how to take such constraints on
diffusion into account when computing the stochastic entropy production. In
Section II we briefly discuss how an SDE for stochastic entropy production may
be derived from Markovian SDEs for a set of system coordinates. The
mathematical treatment of cases where diffusion is restricted is discussed in
Section III. In Section IV we consider the stochastic dynamics of the reduced
density matrix of an open three-level quantum system, subjected to
environmental disturbance manifested as raising and lowering operators
associated with transitions between the levels [5, 6]. The diffusion matrix is
singular but we demonstrate how the stochastic entropy production can still be
evaluated. We go on to compute the environmental component of mean stochastic
entropy production numerically in order to characterise equilibrium and
nonequilibrium stationary states of the system. Application to classical
systems is also possible. Our conclusions are given in Section V.

##  II Stochastic entropy production for Itô processes

We consider a set of coordinates
$\boldsymbol{x}\equiv(x_{1},x_{2},\cdots,x_{N})$ that specify the
configuration of a system, and model their evolution using Markovian
stochastic differential equations, or Itô processes:

$dx_{i}=A_{i}(\boldsymbol{x},t)dt+\sum_{j}B_{ij}(\boldsymbol{x},t)dW_{j},$ (1)

where the $dW_{j}$ are independent Wiener increments. We define [3]

$\displaystyle A_{i}^{{\rm irr}}(\boldsymbol{x},t)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[A_{i}(\boldsymbol{x},t)+\varepsilon_{i}A_{i}(\boldsymbol{\varepsilon}\boldsymbol{x},t)\right]=\varepsilon_{i}A_{i}^{{\rm
irr}}(\boldsymbol{\varepsilon}\boldsymbol{x},t),\qquad$ (2)

and

$A_{i}^{{\rm
rev}}(\boldsymbol{x},t)=\frac{1}{2}\left[A_{i}(\boldsymbol{x},t)-\varepsilon_{i}A_{i}(\boldsymbol{\varepsilon}\boldsymbol{x},t)\right]=-\varepsilon_{i}A_{i}^{{\rm
rev}}(\boldsymbol{\varepsilon}\boldsymbol{x},t),$ (3)

where $\varepsilon_{i}=1$ for variables $x_{i}$ with even parity under time
reversal symmetry (for example position) and $\varepsilon_{i}=-1$ for
variables with odd parity (for example velocity). The notation
$\boldsymbol{\varepsilon}\boldsymbol{x}$ represents
$(\varepsilon_{1}x_{1},\varepsilon_{2}x_{2},\cdots)$. Defining an $N\times N$
diffusion matrix
$\boldsymbol{D}(\boldsymbol{x})=\frac{1}{2}\boldsymbol{B}(\boldsymbol{x})\boldsymbol{B}(\boldsymbol{x})^{\mathsf{T}}$,
the Fokker-Planck equation for the probability density function (pdf)
$p(\boldsymbol{x},t)$ is

$\frac{\partial p}{\partial t}=-\sum_{i}\frac{\partial}{\partial
x_{i}}\left(A_{i}p\right)+\frac{\partial}{\partial x_{i}\partial
x_{j}}\left(D_{ij}p\right).$ (4)

The stochastic entropy production of the system and its environment,
associated with the stochastic motion described by Eq. (1), is a measure of
the difference in probability between pairs of time-reversed sequences of
events, and is defined by [7]

$\Delta s_{{\rm tot}}=\ln\left(\frac{{\rm Prob(forward\>trajectory)}}{{\rm
Prob(backward\,trajectory)}}\right).$ (5)

This is usually separated into system and environmental contributions:

$d\Delta s_{{\rm tot}}=d\Delta s_{\text{sys}}+d\Delta s_{\text{env}},$ (6)

with $d\Delta s_{\text{sys}}=-d\ln p(\boldsymbol{x},t)$. The evolution of the
environmental stochastic entropy production is governed by [4]:

$\begin{split}&d\Delta s_{\text{env}}=-\sum_{i}\frac{\partial
A_{i}^{\text{rev}}(\boldsymbol{x})}{\partial
x_{i}}dt+\sum_{i,j}\Biggl{\\{}\frac{D_{ij}^{-1}(\boldsymbol{x})}{2}\left(A_{i}^{\text{irr}}(\boldsymbol{x})dx_{j}+A_{j}^{\text{irr}}(\boldsymbol{x})dx_{i}\right)\\\
&-\frac{D_{ij}^{-1}(\boldsymbol{x})}{2}\left(\left(\sum_{n}\frac{\partial
D_{jn}(\boldsymbol{x})}{\partial
x_{n}}\right)dx_{i}+\left(\sum_{m}\frac{\partial
D_{im}(\boldsymbol{x})}{\partial
x_{m}}\right)dx_{j}\right)-\frac{D_{ij}^{-1}(\boldsymbol{x})}{2}\left(A_{i}^{\text{rev}}(\boldsymbol{x})A_{j}^{\text{irr}}(\boldsymbol{x})+A_{j}^{\text{rev}}(\boldsymbol{x})A_{i}^{\text{irr}}(\boldsymbol{x})\right)dt\\\
&+\frac{D_{ij}^{-1}(\boldsymbol{x})}{2}\left(A_{j}^{\text{rev}}(\boldsymbol{x})\left(\sum_{m}\frac{\partial
D_{im}(\boldsymbol{x})}{\partial
x_{m}}\right)+A_{i}^{\text{rev}}(\boldsymbol{x})\left(\sum_{n}\frac{\partial
D_{jn}(\boldsymbol{x})}{\partial x_{n}}\right)\right)dt\\\
&+\frac{1}{2}\sum_{k}\Biggl{[}D_{ik}(\boldsymbol{x})\frac{\partial}{\partial
x_{k}}\left(D_{ij}^{-1}(\boldsymbol{x})A_{j}^{\text{irr}}(\boldsymbol{x})\right)+D_{jk}(\boldsymbol{x})\frac{\partial}{\partial
x_{k}}\left(D_{ij}^{-1}(\boldsymbol{x})A_{i}^{\text{irr}}(\boldsymbol{x})\right)\\\
&-D_{ik}(\boldsymbol{x})\frac{\partial}{\partial
x_{k}}\left(D_{ij}^{-1}(\boldsymbol{x})\left(\sum_{n}\frac{\partial
D_{jn}(\boldsymbol{x})}{\partial
x_{n}}\right)\right)-D_{jk}(\boldsymbol{x})\frac{\partial}{\partial
x_{k}}\left(D_{ij}^{-1}(\boldsymbol{x})\left(\sum_{m}\frac{\partial
D_{im}(\boldsymbol{x})}{\partial
x_{m}}\right)\right)\Biggr{]}dt\Biggr{\\}}.\end{split}$ (7)

It may be shown that the average of $d\Delta s_{\text{sys}}$ over all possible
trajectories is related to the incremental change in Gibbs entropy of the
system: $d\langle\Delta s_{{\rm sys}}\rangle=dS_{G}$, to which boundary terms
should be added in certain circumstances [5]. Computing the environmental
stochastic entropy production, on the other hand, presents particular
difficulties if the diffusion matrix is singular, since the inverse matrix
$\boldsymbol{D}^{-1}$ is required in the above expression. This is the problem
we wish to address here..

##  III Defining dynamical and spectator variables

A singular diffusion matrix may be regarded as a consequence of having fewer
independent noise terms than the number of coupled Itô processes. For example,
diffusion might occur on a two dimensional surface within a three dimensional
phase space of system coordinates when the motion is described by three Itô
processes with only two independent Wiener increments. There is a direction at
each point in the phase space in which there is no diffusive current, which
makes the $3\times 3$ diffusion matrix singular. These directions lie parallel
to spatially dependent eigenvectors of the diffusion matrix with zero
eigenvalues, to be referred to as null eigenvectors. The obvious solution is
to establish a reduced set of stochastic differential equations that describe
the random evolution of, in this example, two coordinates on the surface with
the third related deterministically to the other two. We shall denote the
stochastically evolving coordinates as ‘dynamical’ and the remaining
coordinates as ‘spectators’.

It might be possible in simple cases to identify such a reduced set of
coordinates, perhaps by identifying a constant of the motion. However, as we
increase the dimensionality of the phase space and hence the size of the
diffusion matrix, the difficulties in doing so may become insurmountable. We
therefore require a more general treatment of situations with a singular
diffusion matrix.

Let us consider a system described by $N$ variables $x_{i}$, each of which
evolves stochastically according to

$dx_{i}=A_{i}dt+\sum_{j=1}^{M}B_{ij}dW_{j},$ (8)

where the $dW_{j}$ are $M$ independent Wiener increments. According to Itô’s
lemma [8], the differential of a function $f$ of these variables can be
written

$df=\sum_{i=1}^{N}\frac{\partial f}{\partial
x_{i}}dx_{i}+\sum_{i,j=1}^{N}\frac{\partial^{2}f}{\partial x_{i}\partial
x_{j}}D_{ij}dt,$ (9)

where the elements of the $N\times N$ diffusion matrix are
$D_{ij}=\frac{1}{2}\sum_{k=1}^{M}B_{ik}B_{jk}$. Consider such a function where
the stochastic terms in Eq. (9) vanish, i.e.

$\sum_{i=1}^{N}\sum_{j=1}^{M}\frac{\partial f}{\partial x_{i}}B_{ij}dW_{j}=0.$
(10)

Taking the square and using $dW_{i}dW_{j}=\delta_{ij}dt$ we find that

$\displaystyle\sum_{ijkl}\frac{\partial f}{\partial x_{i}}\frac{\partial
f}{\partial x_{k}}B_{ij}B_{kl}\delta_{jl}dt$
$\displaystyle=\sum_{ijk}\frac{\partial f}{\partial x_{i}}\frac{\partial
f}{\partial x_{k}}B_{ij}B_{kj}dt$ $\displaystyle=2\sum_{ik}\frac{\partial
f}{\partial x_{i}}\frac{\partial f}{\partial x_{k}}D_{ik}dt=0,$ (11)

and we therefore have deterministic (noise-free) evolution of $f$ if

$\left(\boldsymbol{\nabla}f\right)^{\mathsf{T}}\boldsymbol{D}\boldsymbol{\nabla}f=0.$
(12)

Such an outcome can arise if $\boldsymbol{\nabla}f$ is an eigenvector of
$\boldsymbol{D}$ with an eigenvalue equal to zero. A matrix is singular if one
or more of its eigenvalues are zero, so we have established that
$\boldsymbol{D}$ is singular if there exists a function $f$ of the stochastic
variables $x_{i}$ that evolves deterministically. For dynamics evolving in an
$N$ dimensional phase space under the influence of $M$ noises where $N>M$, we
conjecture that there will be $L=N-M$ such functions, each associated with one
of the null eigenvectors of $\boldsymbol{D}$. We shall see that this condition
allows us to recast the calculation of the stochastic entropy production to
overcome the problematic singularity of $\boldsymbol{D}$.

###  III.1 Identifying constants of motion

If the deterministic as well as stochastic terms in Eq. (9) vanish then $df=0$
and the function $f$ is a constant of the motion: the evolution of variables,
or coordinates, is constrained to a contour of constant $f$. We can therefore
write

$\boldsymbol{\nabla}f\cdot d\boldsymbol{x}=0,$ (13)

meaning that the infinitesimal vector $d\boldsymbol{x}$ specified by Eq. (8)
is tangential to a contour of $f$. Since $\boldsymbol{\nabla}f$ is also a null
eigenvector of the diffusion matrix, this constraint on $d\boldsymbol{x}$ is
conveniently identified by evaluating the eigenvectors of $\boldsymbol{D}$.

If there are $L=N-M$ constant functions of the coordinates under the dynamics
we can remove $L$ coordinates from the original $N$ dimensional phase space
leaving a reduced $M$ dimensional phase space. Without loss of generality, the
first $M$ coordinates in the set $\left\\{x_{m}\right\\}$ with $m=1,\cdots,M$
will be denoted _dynamical variables_ and the remaining $L$ coordinates
$\left\\{x_{l}\right\\}$ with _ $l=M+1,\cdots,N$ _are designated _spectator
variables_. The $L$ constants of the motion mean we should be able to write
the spectator variables as functions of the dynamical variables, namely
$\left\\{x_{l}\left(\left\\{x_{m}\right\\}\right)\right\\}$. The division into
dynamical and spectator variables is arbitrary, but, as we shall see, some
choices are more convenient for computing stochastic entropy production than
others.

In this labelling scheme, the top left $M\times M$ block of the diffusion
matrix $\boldsymbol{D}$ remains relevant to the stochastic entropy calculation
and will be non-singular, allowing us to use Eq. (7) to compute $\Delta
s_{{\rm env}}$ associated with the evolution, but with $i$ and $j$ ranging
only between 1 and $M$ rather than 1 and $N$. However, if the elements of this
matrix block depend on spectator variables, we need to take this into account
when performing derivatives with respect to the dynamical variables. To
emphasise this point, we write the components of
$\boldsymbol{A},\boldsymbol{D}$ and $\boldsymbol{D}^{-1}$ that appear in Eq.
(7) to show their dependence on the dynamical and spectator variables
explicitly:

$\begin{split}A_{i}&\left(\left\\{x_{m}\right\\},\left\\{x_{l}\left(\left\\{x_{m}\right\\}\right)\right\\}\right)\\\
D_{ij}&\left(\left\\{x_{m}\right\\},\left\\{x_{l}\left(\left\\{x_{m}\right\\}\right)\right\\}\right)\\\
D_{ij}^{-1}&\left(\left\\{x_{m}\right\\},\left\\{x_{l}\left(\left\\{x_{m}\right\\}\right)\right\\}\right),\end{split}$
(14)

and Eq. (7) requires us to evaluate derivatives of $A_{i}$, $D_{ij}$ and
$D_{ij}^{-1}A_{j}$ with respect to the dynamical variables
$\left\\{x_{m}\right\\}$.

We consider derivatives of $D_{ij}$ in the following, but the argument is
easily extend to other expressions. We begin by noting that Eq. (13) can be
separated according to dynamical and spectator variables such that

$\sum_{m=1}^{M}\alpha_{km}dx_{m}+\sum_{l=M+1}^{N}\alpha_{kl}dx_{l}=0,$ (15)

where $\alpha_{km}$ and $\alpha_{kl}$ are the dynamical and spectator
components, respectively, of the $k$th null eigenvector of $\boldsymbol{D}$,
with $k=1,\cdots,L$. Therefore

$\sum_{l=M+1}^{N}\alpha_{kl}dx_{l}=-\sum_{m=1}^{M}\alpha_{km}dx_{m}.$ (16)

Arranging the $\alpha_{km}$ as elements of a rectangular $L\times M$ matrix
$Q_{km}$ and the $\alpha_{kl}$ as elements of a square $L\times L$ matrix
$P_{kl}$ we have

$dx_{l}=-\sum_{k,m}P_{lk}^{-1}Q_{km}dx_{m}=\sum_{m=1}^{M}R_{lm}dx_{m},$ (17)

where $R_{lm}$ is an element of the $L\times M$ matrix
$\boldsymbol{R}=\boldsymbol{P}^{-1}\boldsymbol{Q}$. Next we write

$dD_{ij}=\sum_{m=1}^{M}\frac{\partial D_{ij}}{\partial
x_{m}}dx_{m}+\sum_{l=M+1}^{N}\frac{\partial D_{ij}}{\partial x_{l}}dx_{l},$
(18)

and by substituting $dx_{l}$ from Eq. (17) into Eq. (18), we arrive at the
following expression for the derivative of $D_{ij}$ with respect to the
dynamical coordinate $x_{m}$:

$\frac{dD_{ij}}{dx_{m}}=\frac{\partial D_{ij}}{\partial
x_{m}}+\sum_{l=M+1}^{N}\frac{\partial D_{ij}}{\partial x_{l}}R_{lm}.$ (19)

We were seeking and have identified an additional term on the right hand side.

To summarise, the described framework allows the computation of entropy
production in cases where the diffusion matrix is singular as a result of
constraints on the dynamics through $L$ constants of motion. The method
employs these constants of motion to reduce the dimensionality of the phase
space across which the system evolves such that the appropriately reduced
diffusion matrix is non-singular. Having carried out this transformation, Eq.
(7) can be used to compute entropy production with derivatives determined
according to Eq. (19).

###  III.2 Identifying deterministically evolving functions

We have considered a function $f$ of the stochastic variables
$\left\\{x_{i}\right\\}$ evolving according to

$df=\sum_{i}\frac{\partial f}{\partial
x_{i}}dx_{i}+\sum_{i,j}\frac{\partial^{2}f}{\partial x_{i}\partial
x_{j}}D_{ij}dt,$ (20)

and looked at situations where both the deterministic and stochastic terms in
Eq. (20) vanish, making the function $f$ a constant of motion of the dynamics.
However, it is only necessary for the stochastic terms to vanish for there to
be a restriction on the diffusive motion. We now consider the more general
case where

$df=\left(\sum_{i}\frac{\partial f}{\partial
x_{i}}A_{i}+\sum_{i,j}\frac{\partial^{2}f}{\partial x_{i}\partial
x_{j}}D_{ij}\right)dt$ (21)

is nonzero. Following earlier arguments, $\boldsymbol{\nabla}f$ is still a
null eigenvector of $\boldsymbol{D}$ and we take each null eigenvector to
correspond to a deterministically evolving function. To illustrate this,
consider a system evolving stochastically in two dimensions ($x_{1}$,$x_{2}$)
with a diffusion matrix that possesses a single null eigenvector. A function
$f(x_{1},x_{2},t)$ evolving deterministically according to Eq. (21) allows us
to reduce the number of variables needed to describe the motion from two to
one. The other becomes a spectator variable. To understand this better,
imagine that at time $t_{0}$ the function $f$ is given by

$f(x_{1},x_{2},t_{0})=c_{0},$ (22)

where $c_{0}$ is a constant defining a contour of $f$ in the space. Equation
(22) allows us to express spectator variable $x_{2}$ as a function of
dynamical variable $x_{1}$ at time $t_{0}$.

Between times $t_{0}$ and $t_{1}=t_{0}+dt$, $f$ changes deterministically by
an amount $df_{t_{0}}$ to define a new contour

$f(x_{1},x_{2},t_{1})=c_{1},$ (23)

according to which we can again express $x_{2}$ as a function of $x_{1}$ at
the later time. The evolution of the system is confined to a sequence of
contours of the deterministically evolving function $f$ such that we can
always express $x_{2}$ in terms of $x_{1}.$ We can parametrise the evolution
with a single coordinate and thereby employ a reduced diffusion matrix to
compute an entropy production. The argument easily generalises to an arbitrary
number of dimensions.

As in the previous section we need to consider how derivatives are modified
when we reduce the dimensionality of the phase space. We write
$\boldsymbol{\nabla}f=g\boldsymbol{\alpha}$ where $g$ is equal to
$|\boldsymbol{\nabla}f|$ and $\boldsymbol{\alpha}$ is a normalised null
eigenvector of $\boldsymbol{D}$, , with $|\boldsymbol{\alpha}|=1$, lying
perpendicular to the contour of $f$. Infinitesimal changes in coordinates
$\boldsymbol{x}$ within a timestep $dt$ are given by $d\boldsymbol{x}$, and
$\boldsymbol{x}$ is constrained to pass between points on specified contours
of $f$ at specified times, hence with only certain $d\boldsymbol{x}$ allowed.
The situation is illustrated in Fig. 1.

Figure 1: Illustration of the evolution of system coordinates between contours
defined by a deterministically evolving function $f(\boldsymbol{x},t)$,
constructed to be normal to a spatially dependent null eigenvector
$\boldsymbol{\alpha}$ of the diffusion matrix. There is a limited choice of
stochastic increments $d\boldsymbol{x}$, defined by angle $\theta$ and the
contours visited at the beginning and end of the timestep, which restricts the
diffusive evolution and complicates the computation of the stochastic entropy
production.

The component of $d\boldsymbol{x}$ in the direction normal to $f$ is given by

$d\boldsymbol{x}_{\perp}=|d\boldsymbol{x}|\cos\theta=|\boldsymbol{\alpha}||d\boldsymbol{x}|\cos\theta=\boldsymbol{\alpha}\cdot
d\boldsymbol{x},$ (24)

where the angle $\theta$ is shown in Fig. 1. We can also write

$df=|\boldsymbol{\nabla}f|d\boldsymbol{x}_{\perp}=|\boldsymbol{\nabla}f|\boldsymbol{\alpha}\cdot
d\boldsymbol{x}=g\boldsymbol{\alpha}\cdot d\boldsymbol{x},$ (25)

such that

$g\boldsymbol{\alpha}\cdot d\boldsymbol{x}=\left(\sum_{i}\frac{\partial
f}{\partial x_{i}}A_{i}+\sum_{i,j}\frac{\partial^{2}f}{\partial x_{i}\partial
x_{j}}D_{ij}\right)dt.$ (26)

For notational convenience we rewrite the term in brackets on the right hand
side as $gG$ so that

$\boldsymbol{\alpha}\cdot d\boldsymbol{x}=Gdt.$ (27)

For a set of null eigenvectors of $\boldsymbol{D}$ labelled by $k=1,...,L$
there are several deterministically evolving functions of the coordinates, and
we can specify relationships between increments in dynamical and spectator
variables such that

$\sum_{m=1}^{M}\alpha_{km}dx_{m}+\sum_{l=M+1}^{N}\alpha_{kl}dx_{l}=G_{k}dt.$
(28)

Following the reasoning in Eq. (17) we write

$\begin{split}dx_{l}&=-P_{lk}^{-1}Q_{km}dx_{m}+P_{lk}^{-1}G_{k}dt=R_{lm}dx_{m}+S_{l}dt,\end{split}$
(29)

with implied summation over repeated indices, and substituting Eq. (29) into
Eq. (18) we obtain

$dD_{ij}=\sum_{m=1}^{M}\frac{\partial D_{ij}}{\partial
x_{m}}dx_{m}+\sum_{l=M+1}^{N}\frac{\partial D_{ij}}{\partial
x_{l}}\left(R_{lm}dx_{m}+S_{l}dt\right),$ (30)

such that, as before, the required derivatives of the relevant elements of the
diffusion matrix are

$\frac{dD_{ij}}{dx_{m}}=\frac{\partial D_{ij}}{\partial
x_{m}}+\sum_{l}\frac{\partial D_{ij}}{\partial x_{l}}R_{lm}.$ (31)

We therefore find that contributions to derivatives with respect to dynamical
variables, where the expression in question also depends on spectator
variables, are the same whether $f$ is a constant of the motion or a
deterministic function of the dynamics. This is useful since, in general, we
are unlikely to be able to determine if the singularity of a diffusion matrix
is due to the existence of a constant or a deterministically evolving
function.

###  III.3 Simple example of restricted diffusive evolution

To illustrate the above reasoning, consider the two SDEs

$\displaystyle dx_{1}$ $\displaystyle=\left(1-x_{1}^{2}\right)dW$
$\displaystyle dx_{2}$ $\displaystyle=-\frac{1}{2}x_{2}dt-x_{1}x_{2}dW,$ (32)

which present a case of restricted diffusive evolution since there are two Itô
processes but only one noise. The $2\times 2$ diffusion matrix for the
$(x_{1},x_{2})$ phase space can be shown to be singular. The reversible
deterministic terms $A_{1,2}^{{\rm rev}}$ in both SDEs are zero. We proceed
first by regarding $x_{1}$ as the dynamical variable and $x_{2}$ as a
spectator and use Eq. (7) to write

$\begin{split}&d\Delta s_{\text{tot}}=-d\ln
p_{1}(x_{1},t)+\frac{1}{D_{11}}A_{1}^{\text{irr}}dx_{1}-\frac{1}{D_{11}}\frac{dD_{11}}{dx_{1}}dx_{1}\\\
&+\Biggl{[}D_{11}\frac{d}{dx_{1}}\left(\frac{1}{D_{11}}A_{1}^{\text{irr}}\right)-D_{11}\frac{d}{dx_{1}}\left(\frac{1}{D_{11}}\frac{dD_{11}}{dx_{1}}\right)\Biggr{]}dt,\end{split}$
(33)

where $p_{1}(x_{1},t)=\int p(x_{1},x_{2},t)dx_{2}$ and $p(x_{1},x_{2},t)$
satisfies the Fokker-Planck equation associated with Eqs. (32). Since
$A_{1}^{{\rm irr}}=0$ and $D_{11}=\frac{1}{2}(1-x_{1}^{2})^{2}$ we can take
derivatives unencumbered by implicit dependence on $x_{1}$ arising from
dependence on $x_{2}$.

However, we could just as well decide to compute the stochastic entropy
production by regarding $x_{2}$ as the dynamical variable and $x_{1}$ as the
spectator and write

$\begin{split}&d\Delta s_{\text{tot}}=-d\ln
p_{2}(x_{2},t)+\frac{1}{D_{22}}A_{2}^{\text{irr}}dx_{2}-\frac{1}{D_{22}}\frac{dD_{22}}{dx_{2}}dx_{2}\\\
&+\Biggl{[}D_{22}\frac{d}{dx_{2}}\left(\frac{1}{D_{22}}A_{2}^{\text{irr}}\right)-D_{22}\frac{d}{dx_{2}}\left(\frac{1}{D_{22}}\frac{dD_{22}}{dx_{2}}\right)\Biggr{]}dt,\end{split}$
(34)

with $p_{2}(x_{2},t)=\int p(x_{1},x_{2},t)dx_{1}$, $A_{2}^{{\rm
irr}}=-\frac{1}{2}x_{2}$ and $D_{22}=\frac{1}{2}x_{1}^{2}x_{2}^{2}$. Since
$D_{22}$ depends on the (current) spectator variable $x_{1}$ we have to employ
derivatives like

$\frac{dD_{22}}{dx_{2}}=\frac{\partial D_{22}}{\partial x_{2}}+R\frac{\partial
D_{22}}{\partial x_{1}},$ (35)

and identify the coefficient $R$ using the (single) null eigenvector of
$\boldsymbol{D}$, which may be shown to be proportional to
$(x_{1}x_{2},(1-x_{1}^{2}))^{\mathsf{T}}$. In this example the dynamics
preserve the value of the function $f(t)=x_{1}^{2}(t)+Kx_{2}^{2}(t)-1$ with
arbitrary constant $K$, as long as $f(0)=0$ is imposed as an initial condition
(namely the motion is confined to an ellipse). We can therefore use Eq. (16)
in the form $\alpha_{1}dx_{1}+\alpha_{2}dx_{2}=0$ where $\alpha_{i}$ is the
$i$th component of the null eigenvector. We then obtain

$x_{1}x_{2}dx_{1}=-\left(1-x_{1}^{2}\right)dx_{2},$ (36)

such that $R=-(1-x_{1}^{2})/(x_{1}x_{2})$. The computation of $\Delta
s_{\text{env}}$ using Eqs. (34) and 35) can then proceed.

The point we are making is that in cases of restricted diffusive evolution, we
can divide the stochastic variables arbitrarily into dynamical and spectator
sets. The implication is that some choices of the division might be more
convenient than others; in the case just considered it is more sensible to
regard $x_{2}$ as a spectator variable rather than $x_{1}$.

##  IV An open three-level quantum system with restricted diffusion

###  IV.1 SDEs and selection of spectator variables

We have developed the present framework for computing stochastic entropy
production because there are physical systems of interest where some of the
stochastically evolving variables are spectators. Specifically, we consider
the dynamics of an open quantum system characterised by the stochastic
evolution of its (reduced) density matrix $\rho$. The stochasticity is brought
about by coupling to the environment, as described elsewhere [9, 5, 6]. In
Appendix A it is shown how a Markovian stochastic Lindblad equation for the
evolution of the reduced density matrix of an open system can be derived
starting from the so-called Lindblad operators that specify the dynamical
effect of the environment on the system. Using this formalism we consider a
three-level bosonic system with environmental coupling characterised by the
three raising ($c_{1-3}$) and three lowering ($c_{4-6}$) Lindblad operators
given by

$\begin{split}&c_{1}=\begin{pmatrix}0&0&0\\\ 1&0&0\\\
0&0&0\end{pmatrix}\quad\\!c_{2}=\begin{pmatrix}0&0&0\\\ 0&0&0\\\
1&0&0\end{pmatrix}\quad\\!c_{3}=\begin{pmatrix}0&0&0\\\ 0&0&0\\\
0&1&0\end{pmatrix}\\\ \\\ &c_{4}=\begin{pmatrix}0&1&0\\\ 0&0&0\\\
0&0&0\end{pmatrix}\quad\\!c_{5}=\begin{pmatrix}0&0&1\\\ 0&0&0\\\
0&0&0\end{pmatrix}\quad\\!c_{6}=\begin{pmatrix}0&0&0\\\ 0&0&1\\\
0&0&0\end{pmatrix}\end{split}$ (37)

in a basis of kets $|1\rangle,|2\rangle$ and $|3\rangle$ corresponding to the
three levels. The SDE describing the dynamics of the system is given by

$\begin{split}d\rho&=\sum_{i=1}^{6}\left[\left(c_{i}\rho
c_{i}^{\dagger}-\frac{1}{2}\rho
c_{i}^{\dagger}c_{i}-\frac{1}{2}c_{i}^{\dagger}c_{i}\rho\right)dt\right.\\\
&\left.+\left(\rho c_{i}^{\dagger}+c_{i}\rho-{\rm
Tr}\left[\rho\left(c_{i}+c_{i}^{\dagger}\right)\right]\rho\right)dW_{i}\right].\end{split}$
(38)

A sketch of the inter-level transitions brought about by the $c_{i}$ is given
in Fig. 2.

Figure 2: A sketch of the dynamics of an open three-level quantum system with
system-environment coupling characterised by the set of six raising
($c_{1-3}$) and lowering ($c_{4-6}$) operators. By reducing the probabilities
of the transitions shown in yellow, through multiplying the Lindblad operators
$c_{1,3,5}$ in the stochastic Lindblad equation (38) by a weighting factor
$w<1$, a probability current through states $1\to 3\to 2\to 1$ characterised
by positive mean stochastic entropy production can be generated. For $w=1$ the
system would otherwise adopt an equilibrium state with zero current and zero
stochastic entropy production.

The reduced density matrix $\rho$ for the open quantum system is a complex,
$3\times 3$ Hermitian matrix with a unit trace, corresponding to eight degrees
of freedom. We therefore parametrise $\rho$ in terms of an eight dimensional
vector $\boldsymbol{x}$ evolving as

$d\boldsymbol{x}=\boldsymbol{A}dt+\boldsymbol{B}d\boldsymbol{W},$ (39)

where $\boldsymbol{A}$ is an eight dimensional vector, $\boldsymbol{B}$ is an
$8\times 6$ matrix and $d\boldsymbol{W}$ is a six dimensional vector of
independent Wiener increments. It is clear that with fewer noise terms than
SDEs, the diffusive motion will be restricted in some way.

Figure 3: Environmental component of stochastic entropy production as a
function of time in the three-level quantum system, computed for 25
trajectories with equally weighted Lindblads. The black line represents the
ensemble mean.

In order to proceed we employ the eight Gell-Mann matrices given by

$\begin{split}&\lambda_{1}=\begin{pmatrix}0&1&0\\\ 1&0&0\\\
0&0&0\end{pmatrix}\quad\lambda_{2}=\begin{pmatrix}0&-i&0\\\ i&0&0\\\
0&0&0\end{pmatrix}\quad\lambda_{3}=\begin{pmatrix}1&0&0\\\ 0&-1&0\\\
0&0&0\end{pmatrix}\\\ \\\ &\lambda_{4}=\begin{pmatrix}0&0&1\\\ 0&0&0\\\
1&0&0\end{pmatrix}\quad\lambda_{5}=\begin{pmatrix}0&0&-i\\\ 0&0&0\\\
i&0&0\end{pmatrix}\quad\lambda_{6}=\begin{pmatrix}0&0&0\\\ 0&0&1\\\
0&1&0\end{pmatrix}\\\ \\\ &\lambda_{7}=\begin{pmatrix}0&0&0\\\ 0&0&-i\\\
0&i&0\end{pmatrix}\quad\lambda_{8}=\frac{1}{\sqrt{3}}\begin{pmatrix}1&0&0\\\
0&1&0\\\ 0&0&-2\end{pmatrix}.\end{split}$ (40)

It is convenient to work with eight real variables $s,t,u,v,w,x,y$ and $z$,
defined by ${\rm Tr}\left(\lambda_{1}\rho\right)=s$, ${\rm
Tr}\left(\lambda_{2}\rho\right)=t$, ${\rm Tr}\left(\lambda_{3}\rho\right)=u$
etc, corresponding to the eight components of $\boldsymbol{x}$. We can then
write $\rho$ as

$\rho=\begin{pmatrix}\frac{u}{2}+\frac{\sqrt{3}z}{6}+\frac{1}{3}&\frac{s}{2}-\frac{it}{2}&\frac{v}{2}-\frac{iw}{2}\\\
\frac{s}{2}+\frac{it}{2}&-\frac{u}{2}+\frac{\sqrt{3}z}{6}+\frac{1}{3}&\frac{x}{2}-\frac{iy}{2}\\\
\frac{v}{2}+\frac{iw}{2}&\frac{x}{2}+\frac{iy}{2}&-\frac{\sqrt{3}z}{3}+\frac{1}{3}\end{pmatrix}.$

The SDEs now take the form

$\begin{pmatrix}ds\\\ dt\\\ du\\\ dv\\\ dw\\\ dx\\\ dy\\\
dz\end{pmatrix}=\begin{pmatrix}-2s\\\ -2t\\\ -3u\\\ -2v\\\ -2w\\\ -2x\\\
-2y\\\
-3z\end{pmatrix}dt+\left(\begin{smallmatrix}-sv+x&-s^{2}-u+\frac{\sqrt{3}z}{3}+\frac{2}{3}&-sx+v&-sv&-s^{2}+u+\frac{\sqrt{3}z}{3}+\frac{2}{3}&-sx\\\
-tv-y&-st&-tx+w&-tv&-st&-tx\\\
v\left(1-u\right)&s\left(1-u\right)&-x\left(u+1\right)&-uv&-s\left(u+1\right)&-ux\\\
-v^{2}-\frac{2\sqrt{3}z}{3}+\frac{2}{3}&-sv+x&-vx&u-v^{2}+\frac{\sqrt{3}z}{3}+\frac{2}{3}&-sv&s-vx\\\
-vw&-sw+y&-wx&-vw&-sw&t-wx\\\
-vx&-sx&-x^{2}-\frac{2\sqrt{3}z}{3}+\frac{2}{3}&s-vx&-sx+v&-u-x^{2}+\frac{\sqrt{3}z}{3}+\frac{2}{3}\\\
-vy&-sy&-xy&-t-vy&-sy+w&-xy\\\
\frac{v\left(-3z+\sqrt{3}\right)}{3}&\frac{s\left(-3z+\sqrt{3}\right)}{3}&\frac{x\left(-3z+\sqrt{3}\right)}{3}&-\frac{v\left(3z+2\sqrt{3}\right)}{3}&\frac{s\left(-3z+\sqrt{3}\right)}{3}&-\frac{x\left(3z+2\sqrt{3}\right)}{3}\end{smallmatrix}\right)\begin{pmatrix}dW_{1}\\\
dW_{2}\\\ dW_{3}\\\ dW_{4}\\\ dW_{5}\\\ dW_{6}\end{pmatrix}.$ (41)

The dynamics of $\rho$ are expressed in terms of $N=8$ Itô processes with
$M=6$ noise terms. We therefore expect the diffusion matrix to be singular
(this has been checked using Mathematica) and for there to exist $L=N-M=2$
deterministically evolving functions of the stochastic variables $s,\cdots,z$.
We choose to assign $y$ and $z$ as spectator variables allowing us to focus
instead on the six SDEs:

$\begin{pmatrix}ds\\\ dt\\\ du\\\ dv\\\ dw\\\
dx\end{pmatrix}=\begin{pmatrix}-2s\\\ -2t\\\ -3u\\\ -2v\\\ -2w\\\
-2x\end{pmatrix}dt+\left(\begin{smallmatrix}-sv+x&-s^{2}-u+\frac{\sqrt{3}z}{3}+\frac{2}{3}&-sx+v&-sv&-s^{2}+u+\frac{\sqrt{3}z}{3}+\frac{2}{3}&-sx\\\
-tv-y&-st&-tx+w&-tv&-st&-tx\\\
v\left(1-u\right)&s\left(1-u\right)&-x\left(u+1\right)&-uv&-s\left(u+1\right)&-ux\\\
-v^{2}-\frac{2\sqrt{3}z}{3}+\frac{2}{3}&-sv+x&-vx&u-v^{2}+\frac{\sqrt{3}z}{3}+\frac{2}{3}&-sv&s-vx\\\
-vw&-sw+y&-wx&-vw&-sw&t-wx\\\
-vx&-sx&-x^{2}-\frac{2\sqrt{3}z}{3}+\frac{2}{3}&s-vx&-sx+v&-u-x^{2}+\frac{\sqrt{3}z}{3}+\frac{2}{3}\end{smallmatrix}\right)\begin{pmatrix}dW_{1}\\\
dW_{2}\\\ dW_{3}\\\ dW_{4}\\\ dW_{5}\\\ dW_{6}\end{pmatrix}.$ (42)

The diffusion matrix formed from Eq. (42) via
$\boldsymbol{D}=\frac{1}{2}\boldsymbol{B}\boldsymbol{B}^{\mathsf{T}}$ is too
elaborate to obtain an analytical expression for its inverse. It is possible,
however, to compute an inverse numerically, remembering to append terms to any
derivatives in the matrix elements according to the procedure described in the
previous section. In actual fact, with this choice of spectator variables,
_no_ such additional terms are required, making this route very convenient.

###  IV.2 Equilibrium and nonequilibrium stationary states

Codes have been written to solve the SDEs for the dynamics of the reduced
density matrix and the evolution of the environmental __ component of the
stochastic entropy production. Calculating the system component $\Delta
s_{{\rm sys}}$ requires solution of a Fokker-Planck equation, which is
computationally demanding, but which does not add to the understanding we
develop regarding the _stationary_ states of the system. We comment further on
this point later.

The computational demands of solving Eq. (7) are considerable for the system
under investigation, so only limited ensembles of trajectories were generated.
Runs for environmental stochastic entropy production were typically executed
with a timestep of $dt=10^{-6}$, and the reduced density matrix was initiated
in the condition $s=t=u=v=w=x=y=z=0.1$ throughout.

Figure 3 illustrates the range of environmental stochastic entropy production
for 25 runs. The average environmental entropy production is nearly always
within one standard deviation of zero, which, given the limited statistics,
provides a strong indication that a zero mean rate of environmental stochastic
entropy production has been established in the stationary state. This is
precisely what is to be expected of an equilibrium state, and the mean rate of
system stochastic entropy production ought to be zero as well since it
corresponds in typical situations to the rate of change of Gibbs entropy [5].

We now investigate a nonequilibrium stationary state with a probability
current passing through the system phase space. Our interest in the three-
level system arises precisely because it is the simplest quantum system in
which such a state might be possible. We create a non-equilibrium stationary
state by breaking detailed balance and favouring a pattern of transitions
$|3\rangle\to|2\rangle\to|1\rangle\to|3\rangle$ around the system. We do this
by reducing the coupling strength associated with the Lindblad operators
linked to the opposite pattern. Specifically, we multiply Lindblads $c_{1}$,
$c_{2}$ and $c_{5}$ in (37) by a weighting factor $w<1$ and derive a modified
set of dynamical equations. Such a non-equilibrium stationary state is
expected to be associated with a positive mean rate of stochastic entropy
production. Moreover, we expect the strength with which the system is weighted
towards such a non-equilibrium stationary state to be related to the
irreversibility of its behaviour and to the degree of mean stochastic entropy
production.

Figure 4: Comparison of the environmental stochastic entropy production
computed with timesteps $dt=1/1000000$ and $1/500000$ for a weighting of
$w=0.2$. Averages were taken over 50 runs for each value of timestep. Bands
represent the standard error and solid lines the average. Figure 5:
Environmental stochastic entropy production computed for 10 trajectories with
$w=0.1$ and $dt=10^{-6}$. The black line represents the ensemble mean. Figure
6: Mean rate of environmental stochastic entropy production against $1-w$. The
latter characterises the degree to which detailed balance is broken: a value
$w=1$ corresponds to no breakage and the equilibrium state. Each average was
computed by using linear regression to produce best fit lines for 10
trajectories and then an average was performed over the gradient of these
lines, for a time interval $\Delta t=10^{-4}$. The error bars show the
standard error of these averages.

We indeed observe a mean positive rate of environmental stochastic entropy
production, within statistical errors. In Fig. 4 we check the accuracy of the
calculations for $w=0.2$ by comparing the production for two values of the
timestep and find them to be consistent. Figure 5 shows the environmental
stochastic entropy production for an ensemble of 10 runs with a weighting of
$w=0.1$.

In line with expectations, we see a constant mean rate of environmental
stochastic entropy production which we associate with the system being in a
nonequilibrium stationary state. The mean system stochastic entropy production
is not expected to make a contribution in a stationary state since the Gibbs
entropy of the system is then constant in time. Figure 6 shows the mean rate
of environmental stochastic entropy production for a selection of weightings.
There is clear indication of a relationship between the breakage of detailed
balance and the mean rate of environmental stochastic entropy production. For
the points $w=0.4$ and $w=0$ we have $\langle\Delta
s_{\text{env}}\rangle/\Delta t<0$, a result at odds with second law, but we
attribute this to statistical error, which could be achieved by, for example,
increasing the number of runs used to generate each point. Figure 6 provides
strong support that our approach is a means by which to quantify the
irreversibility of open quantum systems.

##  V Conclusions

We have employed Itô processes to model the dynamics and thermodynamics of a
system interacting with an environment in the absence of detailed information
about the exact configuration of either. Such an approach has frequently been
used in situations described by classical dynamics [10], and recently this has
been extended to quantum systems [5]. In both cases, difficulties arise when
there are fewer independent sources of noise than dimensions of the system
phase space. Diffusion is restricted and the diffusion matrix becomes
singular, which complicates the calculation of stochastic entropy production.

The solution to the problem is simply to eliminate degrees of freedom
(spectator coordinates) from the entropy calculation to account for the
existence of functions of the coordinates that evolve without noise. We have
described a general method for doing so and illustrated it for a particle
occupying a three-level quantum system, thermalised by an environment that
brings about transitions between the levels. The dynamics take place in an
eight dimensional phase space with only six noise terms, and the $8\times 8$
diffusion matrix is singular.

We have shown that a stationary equilibrium state of the system may be
established, corresponding to a zero mean rate of environmental stochastic
entropy production (and implicitly by a zero mean rate of system stochastic
entropy production as well). More interestingly, we have adjusted the
probabilities of the transitions induced by the environment to create a
nonequilibrium stationary state as well, where the particle is made to cycle
through the levels in a particular order. This state is characterised by
nonzero mean environmental stochastic entropy production.

Stochastic entropy production realises Boltzmann’s programme of linking
thermodynamics, and specifically entropy production, to the dynamical
evolution of system coordinates [11]. Calculating entropy production
quantifies the irreversibility of open system behaviour, addressing how
unlikely it is that reversals of sequences of events might be observed.
Refining this framework to accommodate special cases such as restricted
diffusion adds confidence that such an approach is the most appropriate tool
for quantifying irreversibility in open classical and quantum systems.

## References

  * I. J. Ford [2013] I. J. Ford,  _Statistical Physics: an Entropic Approach_ (Wiley, 2013).
  * R. E. Spinney and I. J. Ford [2012] R. E. Spinney and I. J. Ford, Nonequilibrium thermodynamics of stochastic systems with odd and even variables, Physical Review Letters 108, 170603 (2012).
  * Spinney and Ford [2012] R. E. Spinney and I. J. Ford, Entropy production in full phase space for continuous stochastic dynamics, Physical Review E 85, 051113 (2012).
  * Spinney [2012] R. E. Spinney,  _The use of stochastic methods to explore the thermal equilibrium distribution and define entropy production out of equilibrium_ , Ph.D. thesis, UCL (University College London) (2012).
  * Matos  _et al._ [2022] D. Matos, L. Kantorovich, and I. J. Ford, Stochastic entropy production for continuous measurements of an open quantum system, Journal of Physics Communications 6, 125003 (2022).
  * Clarke and Ford [2023] C. L. Clarke and I. J. Ford, Stochastic entropy production associated with quantum measurement in a framework of Markovian quantum state diffusion, arXiv:2301.08197 [quant-ph] (2023).
  * Seifert [2005] U. Seifert, Entropy production along a stochastic trajectory and an integral fluctuation theorem, Physical Review Letters 95, 040602 (2005).
  * Gardiner [2009] C. Gardiner,  _Handbook of Stochastic Methods_ , Vol. 4 (Springer Berlin, 2009).
  * Clarke [2021] C. L. Clarke,  _Irreversibility Measures in a Quantum Setting_ , Ph.D. thesis, UCL (University College London) (2021).
  * Seifert [2008] U. Seifert, Stochastic thermodynamics: principles and perspectives, European Physical Journal B 64, 423 (2008).
  * Cercignani [1998] C. Cercignani,  _Ludwig Boltzmann: the man who trusted atoms_ (Oxford, 1998).
  * Jacobs [2014] K. Jacobs,  _Quantum Measurement Theory and its Applications_ (Cambridge University Press, 2014).
  * Walls  _et al._ [2022] S. M. Walls, J. M. Schachter, H. Qian, and I. J. Ford, Stochastic quantum trajectories demonstrate the Quantum Zeno Effect in an open spin system, arXiv:2209.10626 [quant-ph] (2022).
  * Breuer and Petruccione [2002] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2002).

##  Appendix A SDEs for quantum state diffusion

The reduced density matrix $\rho$ describing an open quantum system is
considered to evolve stochastically in a timestep $dt$ according to quantum
maps specified by Kraus operators [12], the details of which are provided
here.

We consider $i$ pairs of Kraus operators:

$M_{i\pm}=\sqrt{\frac{P_{i}}{2}}\left(\mathbb{I}-\frac{1}{2P_{i}}c_{i}^{\dagger}c_{i}dt\pm\frac{1}{\sqrt{P_{i}}}c_{i}\sqrt{dt}\right),$
(43)

where $c_{i}$ are the Lindblad operators that specify the mode of coupling
between the system and environment. The meaning of the $P_{i}$ will become
clear shortly. A similar scheme has been outlined elsewhere [9, 5, 13, 6], in
which the $P_{i}$ do not appear.

The Kraus operators satisfy

$\displaystyle M_{i\pm}^{\dagger}M_{i\pm}$
$\displaystyle=\frac{P_{i}}{2}\left(\mathbb{I}\pm\frac{1}{\sqrt{P_{i}}}(c_{i}+c_{i}^{\dagger})\sqrt{dt}\right),$
(44)

so that $M_{i+}^{\dagger}M_{i+}+M_{i-}^{\dagger}M_{i-}=P_{i}\mathbb{I}$ and
the completeness relation [14]

$\sum_{i}\left(M_{i+}^{\dagger}M_{i+}+M_{i-}^{\dagger}M_{i-}\right)=\mathbb{I},$
(45)

holds if $\sum_{i}P_{i}=1$.

A set of $2i$ reduced density matrices, all of which are positive definite and
of unit trace, are reachable in the timestep starting from $\rho$, namely

$\rho^{\prime i\pm}=\frac{M_{i\pm}\rho M_{i\pm}^{\dagger}}{{\rm
Tr}\left(M_{i\pm}\rho M_{i\pm}^{\dagger}\right)}.$ (46)

The probabilities of making the transitions to these targets are given by

$\displaystyle p_{i\pm}$ $\displaystyle={\rm Tr}\left(M_{i\pm}\rho
M_{i\pm}^{\dagger}\right)$
$\displaystyle=\frac{P_{i}}{2}\left(1\pm\sqrt{\frac{dt}{P_{i}}}{\rm
Tr}\left[(c_{i}+c_{i}^{\dagger})\rho\right]\right),$ (47)

which satisfy $p_{i+}+p_{i-}=P_{i}$ and
$p_{i+}-p_{i-}=P_{i}\sqrt{\frac{dt}{P_{i}}}{\rm
Tr}\left[(c_{i}+c_{i}^{\dagger})\rho\right]$. $P_{i}$ is therefore the
probability that one of the $i$th pair of Kraus operators is selected for the
map: it is the probability that the environment disturbs the system in such a
way as to transform $\rho$ into either $\rho^{\prime i+}$ or $\rho^{\prime
i-}$. As a consequence of these dynamical rules, the reduced density matrix is
driven along a Brownian trajectory under the influence of the environment.
Since the Kraus operators reduce to a multiple of the identity as $dt\to 0$,
the Brownian path is continuous. No quantum jumps are allowed (and neither are
they necessary [5]).

The statistically averaged form of the new reduced density matrix is

$\displaystyle\langle\rho^{\prime}\rangle$
$\displaystyle=\sum_{i}\left(p_{i+}\frac{M_{i+}\rho M_{i+}^{\dagger}}{{\rm
Tr}\left(M_{i+}\rho M_{i+}^{\dagger}\right)}+p_{i-}\frac{M_{i-}\rho
M_{i-}^{\dagger}}{{\rm Tr}\left(M_{i-}\rho M_{i-}^{\dagger}\right)}\right)$
$\displaystyle=\sum_{i}\left(M_{i+}\rho M_{i+}^{\dagger}+M_{i-}\rho
M_{i-}^{\dagger}\right),$ (48)

which corresponds in notation with the standard Kraus map [12].

The average evolution having been established, we now derive an SDE that
describes the stochastic evolution of the reduced density matrix. The possible
increments in the timestep are $d\rho^{i\pm}=\rho^{\prime i\pm}-\rho$ where

$\displaystyle d\rho^{i\pm}$ $\displaystyle=\frac{1}{P_{i}}\left(c_{i}\rho
c_{i}^{\dagger}-\frac{1}{2}\rho
c_{i}^{\dagger}c_{i}-\frac{1}{2}c_{i}^{\dagger}c_{i}\rho\right)dt$
$\displaystyle-$ $\displaystyle\frac{1}{P_{i}}\left(\rho
c_{i}^{\dagger}+c_{i}\rho-{\rm
Tr}\left[\rho\left(c_{i}+c_{i}^{\dagger}\right)\right]\rho\right){\rm
Tr}\left[\rho\left(c_{i}+c_{i}^{\dagger}\right)\right]dt$ $\displaystyle\pm$
$\displaystyle\frac{1}{\sqrt{P_{i}}}\left(\rho c_{i}^{\dagger}+c_{i}\rho-{\rm
Tr}\left[\rho\left(c_{i}+c_{i}^{\dagger}\right)\right]\rho\right)\sqrt{dt}.$
(49)

The average increment in the reduced density matrix is then

$\langle d\rho\rangle=\sum_{i}\left(p_{i+}d\rho^{i+}+p_{i-}d\rho^{i-}\right),$
(50)

which leads after some manipulation to

$\langle d\rho\rangle=\sum_{i}\left(c_{i}\rho c_{i}^{\dagger}-\frac{1}{2}\rho
c_{i}^{\dagger}c_{i}-\frac{1}{2}c_{i}^{\dagger}c_{i}\rho\right)dt.$ (51)

The terms in the sum are the dissipators associated with each Lindblad
operator, in the form in which they appear in the standard Lindblad equation
[12].

Next we compute the variance of $d\rho$. For clarity, we write

$d\rho^{i\pm}=\frac{1}{P_{i}}A_{i}dt\pm\frac{1}{\sqrt{P_{i}}}B_{i}\sqrt{dt},$
(52)

where

$\displaystyle A_{i}$ $\displaystyle=c_{i}\rho c_{i}^{\dagger}-\frac{1}{2}\rho
c_{i}^{\dagger}c_{i}-\frac{1}{2}c_{i}^{\dagger}c_{i}\rho$
$\displaystyle-\left(\rho c_{i}^{\dagger}+c_{i}\rho-{\rm
Tr}\left[\rho\left(c_{i}+c_{i}^{\dagger}\right)\right]\rho\right){\rm
Tr}\left[\rho\left(c_{i}+c_{i}^{\dagger}\right)\right],$ (53)

and $B_{i}=\rho c_{i}^{\dagger}+c_{i}\rho-{\rm
Tr}\left[\rho\left(c_{i}+c_{i}^{\dagger}\right)\right]\rho$. We then construct
the average of $\left(d\rho^{i\pm}-\langle d\rho\rangle\right)^{2}$ to lowest
order in $dt$. The $A_{i}$ and $\langle d\rho\rangle$ terms do not contribute
because they are already of order $dt$ and we get

$\displaystyle\sum_{i}\left(p_{i+}\left(d\rho^{i+}-\langle
d\rho\rangle\right)^{2}+p_{i-}\left(d\rho^{i-}-\langle
d\rho\rangle\right)^{2}\right)$
$\displaystyle=\sum_{i}\left(p_{i+}\frac{1}{P_{i}}B_{i}^{2}dt+p_{i-}\frac{1}{P_{i}}B_{i}^{2}dt\right)=\sum_{i}B_{i}^{2}dt,$
(54)

suggesting that the Itô process governing the evolution of $\rho$ is

$\displaystyle d\rho$ $\displaystyle=\langle d\rho\rangle+\sum_{i}B_{i}dW_{i}$
$\displaystyle=\sum_{i}\Big{[}\left(c_{i}\rho c_{i}^{\dagger}-\frac{1}{2}\rho
c_{i}^{\dagger}c_{i}-\frac{1}{2}c_{i}^{\dagger}c_{i}\rho\right)dt$
$\displaystyle+\left(\rho c_{i}^{\dagger}+c_{i}\rho-{\rm
Tr}\left[\rho\left(c_{i}+c_{i}^{\dagger}\right)\right]\rho\right)dW_{i}\Big{]},$
(55)

and this is the form employed in Eq. (38).

The Kraus operators in Eq. (43) that underpin our framework are specified in
terms of $P_{i}$ parameters that are to be interpreted as the probabilities of
selection of one of the $i$th pair of Kraus operators for the transformation
of the reduced density matrix. By design, these parameters do not appear in
the SDE, but it is worth discussing how they might be chosen for use in a
Monte Carlo simulation, for example. If there are $N$ Lindblads, a possible
choice could be $P_{i}=1/N$. This is not completely satisfactory, however,
since it would give equal selection weight to the Kraus operators irrespective
of the degree of coupling between the system and its environment through each
Lindblad operator. It would be better if we ascribe a zero probability of
selection to a Lindblad with zero coupling strength.

A solution would be to make $P_{i}$ dependent on the norm of the matrix
representing the Lindblad operator. We could use the Frobenius norm
$\bigparallel
c_{i}\bigparallel_{\mathrm{F}}=\sqrt{\mathrm{Tr}(c_{i}^{\dagger}c_{i})}$, for
example, and employ probabilities of selection

$P_{i}=\frac{\bigparallel c_{i}\bigparallel_{\mathrm{F}}}{\sum_{i}\bigparallel
c_{i}\bigparallel_{\mathrm{F}}}.$ (56)

Adding a null Lindblad $c_{i}=0$ to the set would not change the Kraus
operators representing the other Lindblads, and furthermore such a null
Lindblad would have zero probability of being selected. The matter is somewhat
academic since the $P_{i}$ do not affect the derived form of the stochastic
dynamics, by design, but it does have impact on the elegance of the framework
that emerges for modelling the quantum state diffusion.


* [399] Heath B Henninger, Shawn P Reese, Andrew E Anderson, and Jeffrey A Weiss.  Validation of computational models in biomechanics.  Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine, 224(7):801–812, 2010.
  * [400] Xue Ying.  An overview of overfitting and its solutions.  In Journal of physics: Conference series, volume 1168, page 022022. IOP Publishing, 2019.
  * [401] Benyamin Ghojogh and Mark Crowley.  The theory behind overfitting, cross validation, regularization, bagging, and boosting: tutorial.  arXiv preprint arXiv:1905.12787, 2019.
  * [402] Peter Bühlmann and Sara Van De Geer.  Statistics for high-dimensional data: methods, theory and applications.  Springer Science & Business Media, 2011.
  * [403] Junhua Ding, Xin-Hua Hu, and Venkat Gudivada.  A machine learning based framework for verification and validation of massive scale image data.  IEEE Transactions on Big Data, 7(2):451–467, 2017.
  * [404] Kexin Pei, Linjie Zhu, Yinzhi Cao, Junfeng Yang, Carl Vondrick, and Suman Jana.  Towards practical verification of machine learning: The case of computer vision systems.  arXiv preprint arXiv:1712.01785, 2017.
  * [405] Weiming Xiang, Patrick Musau, Ayana A Wild, Diego Manzanas Lopez, Nathaniel Hamilton, Xiaodong Yang, Joel Rosenfeld, and Taylor T Johnson.  Verification for machine learning, autonomy, and neural networks survey.  arXiv preprint arXiv:1810.01989, 2018.
  * [406] Xiaowei Huang, Marta Kwiatkowska, Sen Wang, and Min Wu.  Safety verification of deep neural networks.  In Computer Aided Verification: 29th International Conference, CAV 2017, Heidelberg, Germany, July 24-28, 2017, Proceedings, Part I 30, pages 3–29. Springer, 2017.
  * [407] Nicholas Carlini, Anish Athalye, Nicolas Papernot, Wieland Brendel, Jonas Rauber, Dimitris Tsipras, Ian Goodfellow, Aleksander Madry, and Alexey Kurakin.  On evaluating adversarial robustness.  arXiv preprint arXiv:1902.06705, 2019.
  * [408] Wenjie Ruan, Xiaowei Huang, and Marta Kwiatkowska.  Reachability analysis of deep neural networks with provable guarantees.  arXiv preprint arXiv:1805.02242, 2018.
  * [409] Xiaowei Huang, Daniel Kroening, Wenjie Ruan, James Sharp, Youcheng Sun, Emese Thamo, Min Wu, and Xinping Yi.  A survey of safety and trustworthiness of deep neural networks: Verification, testing, adversarial attack and defence, and interpretability.  Computer Science Review, 37:100270, 2020.
  * [410] Sergio Conti, Stefan Müller, and Michael Ortiz.  Data-driven problems in elasticity.  Archive for Rational Mechanics and Analysis, 229(1):79–123, 2018.
  * [411] Andrius Vabalas, Emma Gowen, Ellen Poliakoff, and Alexander J Casson.  Machine learning algorithm validation with a limited sample size.  PloS one, 14(11):e0224365, 2019.
  * [412] Christoph Bergmeir and José M Benítez.  On the use of cross-validation for time series predictor evaluation.  Information Sciences, 191:192–213, 2012.
  * [413] Arthur E Hoerl and Robert W Kennard.  Ridge regression: Biased estimation for nonorthogonal problems.  Technometrics, 12(1):55–67, 1970.
  * [414] Twan Van Laarhoven.  L2 regularization versus batch and weight normalization.  arXiv preprint arXiv:1706.05350, 2017.
  * [415] Ilya Loshchilov and Frank Hutter.  Decoupled weight decay regularization.  arXiv preprint arXiv:1711.05101, 2017.
  * [416] Hui Zou and Trevor Hastie.  Regularization and variable selection via the elastic net.  Journal of the royal statistical society: series B (statistical methodology), 67(2):301–320, 2005.
  * [417] Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov.  Dropout: a simple way to prevent neural networks from overfitting.  The journal of machine learning research, 15(1):1929–1958, 2014.
  * [418] Karl A Kalina, Lennart Linden, Jörg Brummund, Philipp Metsch, and Markus Kästner.  Automated constitutive modeling of isotropic hyperelasticity based on artificial neural networks.  Computational Mechanics, 69(1):213–232, 2022.
  * [419] David González, Francisco Chinesta, and Elías Cueto.  Thermodynamically consistent data-driven computational mechanics.  Continuum Mechanics and Thermodynamics, 31(1):239–253, 2019.
  * [420] David Stutz, Matthias Hein, and Bernt Schiele.  Disentangling adversarial robustness and generalization.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 6976–6987, 2019.
  * [421] Dongxian Wu, Shu-Tao Xia, and Yisen Wang.  Adversarial weight perturbation helps robust generalization.  Advances in Neural Information Processing Systems, 33:2958–2969, 2020.
  * [422] Aditi Raghunathan, Sang Michael Xie, Fanny Yang, John C Duchi, and Percy Liang.  Adversarial training can hurt generalization.  arXiv preprint arXiv:1906.06032, 2019.
  * [423] Davis Blalock, Jose Javier Gonzalez Ortiz, Jonathan Frankle, and John Guttag.  What is the state of neural network pruning?  Proceedings of machine learning and systems, 2:129–146, 2020.
  * [424] Knut A Meyer and Fredrik Ekre.  Thermodynamically consistent neural network plasticity modeling and discovery of evolution laws.  engrXiv preprint https://doi.org/10.31224/2961, 2023.
  * [425] Christoph Torens, Franz Juenger, Sebastian Schirmer, Simon Schopferer, Theresa D Maienschein, and Johann C Dauer.  Machine learning verification and safety for unmanned aircraft-a literature study.  In AIAA Scitech 2022 Forum, page 1133, 2022.