# Novel split quality measures for stratified multilabel Cross Validation with
application to large and sparse gene ontology datasets

Henri Tiittanen  Institute of Biotechnology, Helsinki Institute of Life
Sciences (HiLife), University of Helsinki, Helsinki, Finland  Liisa Holm
Institute of Biotechnology, Helsinki Institute of Life Sciences (HiLife),
University of Helsinki, Helsinki, Finland  Organismal and Evotionary Biology
Research Program, Faculty of Biosciences, University of Helsinki  Petri
Törönen  Institute of Biotechnology, Helsinki Institute of Life Sciences
(HiLife), University of Helsinki, Helsinki, Finland

##  1 Abstract

Multilabel learning is an important topic in machine learning research.
Evaluating models in multilabel settings requires specific cross validation
methods designed for multilabel data. In this article, we show that the most
widely used cross validation split quality measures do not behave adequately
with multilabel data that has strong class imbalance. We present improved
measures and an algorithm, optisplit, for optimising cross validations splits.
We present an extensive comparison of various types of cross validation
methods in which we show that optisplit produces more even cross validation
splits than the existing methods and that it is fast enough to be used on big
Gene Ontology (GO) datasets.

##  2 Introduction

Cross validation is a central procedure in machine learning, statistics and
other fields. It evaluates model performance by testing models on data points,
excluded from training dataset. In standard cross validation a dataset $D$ is
split randomly into $k$ non overlapping subsets (folds) $D_{i},i\in k$. A
model is trained for each fold $i$ on the data $\cup_{j\in k|j\neq i}D_{j}$
and evaluated on $D_{i}$. So one subset, $D_{i}$, is left out of the training
process and used as an evaluation data. The averaged result over all the folds
represents the final performance. Cross validation is typically used when the
amount of data is limited. When abundant amount / very large amounts of data
is available, one can also use standard training test splits in place of cross
validation. Cross validation is further used, besides traditional model
evaluation, also when various predictive models are combined in classifier
stacking [11] and when various parameters in predictive models are optimized.
The typically used random split approach assumes that the positive and
negative class distributions are balanced. If the class distributions are
imbalanced the resulting splits may not allow efficient learning. As an
example, suppose that in binary classification settings one randomly generated
fold contains all the data points in the data belonging to one of the classes.
Then the corresponding training set consisting of the rest of the folds does
not contain any data points of that class and the model cannot learn anything
about the class.

Stratified cross validation methods are variants of cross validation that
ensure that the class distributions of the folds are close to the class
distributions of the whole data. With class imbalanced data these stratified
methods ensure especially the distribution of smaller classes. In this work we
focus on stratified cross validation applied to multilabel classification.
Multilabel classification presents additional challenges for stratified cross
validation as each data point can belong to multiple classes simultaneously
(see figure 1). Here each class represents a separate task for training and
evaluation. The cross validation split should be formed so that the correct
class distributions are maintained for all classes and all folds at the same
time. This ensures that the classifier to be evaluated can be trained and
tested effectively with all of the created data splits. Random splitting has
been historically a popular choice on multilabel data but is has been shown to
lead to poor results [6].

Figure 1: Examples of stratified multilabel cross validation with splits of
different quality. 1) presents a target data with 24 data points and three
classes. This could be a subset of a sparse high dimensional dataset. 2) Now,
the data is splitted into four cross validation folds in two ways, one
representing a random split and the other representing a stratified split
optimized over all classes. The ideal distribution of the positive data points
in each fold would be 2/6 for the first class and 1/6 for the second and third
class so that they would follow the class distributions of the whole data.
Random split results in very unbalanced folds while the stratified split
follows closely the data distribution for all classes. 3) Here, the well
balanced stratified split is rearranged for increased clarity.

As dataset sizes are growing across application domains, multilabel and
extreme classification are of growing interest and importance [2]. Thus,
assessing model quality on these settings is getting more and more important.
Furthermore, the model assessment is critical step when different predictions
are combined in various classifier stacking methods [11]. However, we show in
Section 4.1 and Section 6 that the current stratified cross validation methods
use suboptimal split quality measures. Also, many existing methods are
impractically slow for extremely big classification datasets. Here we present
a new algorithm for generating multilabel cross validation datasets based on
optimising the global distribution of all classes.

The methods presented here are developed in the context of gene ontology (GO)
data [1] for automated protein function prediction task [12]. The GO datasets
used here are high dimensional and contain over half a million data points.
The high number of classes, often in the order of thousands to tens of
thousands, presents a challenge for analysis. Furthermore, the classes have a
hierarchical tree structure, which causes the class distributions to be
strongly imbalanced. So, as a result, most classes contain a very small number
of data points while some classes have a very large number of data points. The
high number of small classes, i.e. classes that have few positive data points,
makes the negative data abundant for most of the classes and results in
difficulties to have enough positive data points to train the models.

This article is organised as follows. In Section 3, we review the related
work. In Section 4, we show that the currently most widely used stratified
multilabel cross validation split quality measure is inadequate. We present
new better measures and demonstrate the behaviour of these first with
synthetic data splits. In Section 5, we present an algorithm for optimising
cross validation splits with respect to any selected split quality measures.
In Section 6, we present a comparison of the algorithms on a wide range of
real-world datasets and conclude that the new algorithm is the best practical
choice for getting balanced splits as measured by the new measures. Finally,
we conclude with a discussion and present ideas for future work in Section 7.

##  3 Related work

The most widely known method for generating cross validation splits with
balanced class distributions is Iterative Stratification (IS) [6]. IS works by
dividing the data points evenly into the folds one class at a time. It always
chooses the class with the fewest positive data points for the processing, and
breaks ties first by the largest number of desired data points and further
randomly. Smaller classes are more difficult to balance equally among the
folds so starting from them makes sure they get well distributed. Bigger
classes are easier to distribute so distributing them later is justifiable.
Iterative stratification has also been extended to consider second order
relationships between labels. This method is known as Second Order Iterative
Stratification (SOIS) [8].

The recently introduced Stratified Sampling (SS) algorithm [5] is designed to
produce balanced train/test splits for extreme classification data with a high
number of data points and dimensions. It should be faster to use than
iterative stratification variants and produce splits with better
distributions. This method calculates the proportion of each class in training
and test sets. Then for each data point a score is calculated over its
positive classes and the data points with the highest scores are redistributed
from one to the other partition. This method needs three parameters that have
to be adjusted according to the data and it does not produce cross validation
splits directly but training/test splits.

The partitioning method based on stratified random sampling (PMBSRS) [4] is
based on using the similarity of the label distributions between data points
to group them and then dividing them into the cross validation folds [4]. A
similarity score is defined as the product of the relative frequencies of the
positive labels present in each data point. Then the data points are ordered
by the score into a list which is cut into as many disjoint subsets S as is
the desired number of cross validation folds. Each cross validation fold is
then generated by randomly selecting items without replacement from each set
$s\in S$ so that each fold gets an equal proportion of the samples from each
$s$. The end result is that each fold contains elements with different scores.
As a non iterative method resembling the basic random sample method this is
computationally less expensive than the iterative methods. However, this does
not measure the split quality directly but is more aimed to ensure that all
folds contain an equal amount of differently sized classes.

##  4 Cross validation split quality measures

The current cross validation split quality measures are either evaluating the
quality of the folds directly or applying some model on the folds and
comparing the learning results. Here we focus on directly comparing the
quality of the folds in order to make the comparisons model independent. The
most commonly used measures in literature are the Labels Distribution (LD)
measure and the Examples Distribution (ED) measure [6].

Let $n$ be the number of data points, $k$ be the number of cross validation
folds and $q$ be the number of classes. $D\in\\{0,1\\}^{n\times q}$ denotes
the multilabel target set and $S_{j},j\in 1\dots k$ denotes the folds which
are disjoint subsets of $D$. The subsets of $D$ and $S_{j}$ containing
positive data points of label $i \in 1\dots q$ are denoted as $D^{i}$ and
$S_{j}^{i}$

We define the positive and negative frequencies for fold $j$ and label $i$ as
$p_{j}^{i}=|S_{j}^{i}|/|S_{j}|$ and $1-p_{j}^{i}$, respectively. Similarly for
the whole data, positive frequency as $d^{i}=|D^{i}|/|D|$ and negative
frequency as $1-d^{i}$.

Then, we define

$\textsc{LD}=\frac{1}{q}\sum_{i=1}^{q}\left(\frac{1}{k}\sum_{j=1}^{k}\left|\frac{p_{j}^{i}}{1-p_{j}^{i}}-\frac{d^{i}}{1-d^{i}}\right|\right)$
(1)

and

$\textsc{ED}=\frac{1}{k}\sum_{j=1}^{k}\left||S_{j}|-\frac{|D|}{k}\right|.$ (2)

Intuitively, LD measures how the distribution of the positive and negative
data points of each label in each subset compares to the distribution in the
whole data 111LD compares odds. Odds := positive to negative ratio. and ED
measures the deviation between the number of data points of each subset and
the desired number of data points. Since the exact equality of the fold sizes
is not generally important in practice, the ED score is merely useful in
checking that the fold size differences are sufficiently small compared to the
data size.

As noted in Section 4.1, it is especially important for training and
evaluation of models that the distributions of small classes in the cross
validation folds are well balanced. Small classes are hardest to split well,
since there are fewer ways to distribute the data points and even small
differences in fold distributions give big relative differences. A good split
quality measure should be able to correctly quantify the quality of the folds
even on small classes.

The LD (Equation 1), is defined as the arithmetic mean of the differences
between the positive frequencies of a cross validation fold and of the whole
data over all classes and folds, using a transformation of the form
$f(x)=x/(1-x)$ for both frequencies. Since the data size is fixed, the
frequencies $p^{i}_{j}$ and $d^{i}$ are linear functions of the class sizes
$|S_{j}^{i}|$ and $|D^{i}|$, respectively. Since $x<x/(1-x),\forall x\in[0,1)$
and $f^{\prime}(x)=1/(1-x^{2})>1$, the transformed quantities grow faster than
the class size, resulting that the difference is greater for bigger classes
than it should be. That is, the same absolute difference of $p_{j}^{i}$ from
$d^{i}$ results in a larger contribution by bigger classes.222The contribution
is equal for the special case $p_{j}^{i}=d^{i}$ when the difference is zero.

Hence, as an improved alternative, we propose the relative Labels Distribution
(rLD) measure

$\text{rLD}=\frac{1}{q}\sum_{i=1}^{q}\left(\frac{1}{k}\sum_{j=1}^{k}\left|\frac{d^{i}-p_{j}^{i}}{d^{i}}\right|\right).$
(3)

Intuitively, rLD measures the linear difference of the ratios of the positive
frequencies between each fold and the whole data weighted by the class size.
Since the difference is calculated without the offending nonlinear
transformation, the measure is class size independent, while still having the
same main operating principle as LD of comparing the distributions
(frequencies) of folds to those of the whole data.

In addition, we present another measure that is insensitive to the class size,
the Delta-Class Proportion (DCP)

$\textsc{DCP}=\frac{1}{q}\sum_{i=1}^{q}\left|\frac{\max_{j\in 1\dots
k}(|S_{j}^{i}|)}{|S^{i}|}-\frac{1}{k}\right|,$ (4)

where the first part of the function represents the observed result and 1/k is
the positive frequency of a flat distribution. Compared to rLD this measure
does not measure the relative distributions of other folds than the largest.
Hence, DCP is a coarser measure that is used here for comparison purposes.

###  4.1 Comparison of split quality measures

Here we compare the behavior of the measures presented in Section 4 with a
synthetic dataset. We generate different errors to class distribution with
various splits of this data. Next, we show how different measures respond to
these errors.

A good synthetic data should have the following charasteristics: 1) it should
be large and sparse, 2) there should be high class imbalance. In order to
achieve these conditions, the data was constructed as follows: The data
consists of a binary target matrix $D\in\\{0,1\\}^{n\times q}$ with
$n=100000,q=100$, that is, 100000 data points and 100 classes. In order to
make the data class imbalanced, the positive class sizes of the data were set
to vary from $2k$ to $n/2$, so that $|S_{j}|=2k+\frac{j(n/2-2k)}{q},j\in
0\dots q-1$.

In order to ensure that the results are consistent, the synthetic data should
be divided into cross validation folds that are of varying quality: 1) perfect
(all class distributions are equal among folds), 2) with varying levels of
error in the distributions between folds.

We divided the data $D$ into 10 cross validation folds in three ways: 1) all
folds were equally sized, 2) there was a 20% increase in one fold, 20%
decrease in one fold, and the rest were equal and 3) the class was missing
from one fold while the other folds were equally sized. These are denoted as
Equal, Difference and One missing, respectively. These three alternatives
represent clearly different levels of error in data splitting, and measures
should separate them from each other. Furthermore, all the classes of $D$ had
the same fold distributions while the absolute class sizes were different.
Therefore, a good split quality measure should give a similar score for all
classes, irrespective of their size. In practical settings it is especially
important to correctly divide the smaller classes since the bigger classes are
naturally better distributed.

We evaluated these three splits splits on the synthetic data using the
measures LD, rLD and DCP. The results for each measure are presented in the
subplots of the Figure 2. In each plot, class size is plotted against the
class specific score given by the measure (smaller is better) for the
different folds. The results confirm that the widely used LD depends on the
class size. LD gives smaller values to smaller classes and higher values to
bigger classes even though their fold distributions are identical (see Figure
2 (a)). The behaviour of rLD and DCP on synthetic data, Figure 2 (b) and (c),
show that both rLD and DCP are not affected by the class size. Therefore, it
is recommendable to use them instead of LD for measuring cross validation
split quality.

(a)  (b)  (c)

Figure 2: Behaviour of cross validation split quality measures on synthetic
data with cross validation splits of varying quality. The cross validation
folds have the same distribution for each class, so the scores should be
equal. However, LD score depends on the class size, while rLD and DCP give the
same weight for all classes. Note that as DCP concentrates on the size of the
largest fold, it behaves differently on data where the class is entirely
missing from one fold giving it comparatively higher score than the other
measures.

##  5 Algorithm

In this section, we present a new general algorithm, optisplit, for optimising
cross validation splits with respect to any split quality measure that can
produce class specific scores. optisplit can also be used to generate standard
train-test splits by generating a cross validation split of size $n$ with fold
size of the desired test set and then forming the training set from the folds
$1\dots n-1$ and the test set from the fold $n$. Unlike some of the competing
methods, optisplit does not need any data specific parameters that have to be
adjusted for different datasets. The details of optisplit are presented
formally in Algorithm 1.

Intuitively, we start by randomly generating the initial $k$ subsets i.e.
cross validation folds. Let $M$ be a measure that computes the quality of the
folds with respect to a single class $c$. We calculate the score vector $M(c)$
for all classes $c\in C$ and calculate the global loss function $L=\sum_{c\in
C}M(c)$.

Let $L_{0}$ be the initial value of $L$. Let $A$ be the index of the class
with the highest $M$. We balance the class $A$ by moving data points from
folds with excess data points to folds without enough data points so that the
class distributions are balanced for the class $A$ (Function balance in
Algorithm 1). After processing the class $A$ the global loss $L$ is
recalculated. Let $L_{1}$ be the recalculated loss. If $L_{1}>L_{0}$, undo the
changes and move to the next worse class. Otherwise keep the modification and
continue to the next worse class. Continue this process for all classes for
many epochs until either the global loss does not improve any more or a
desired max iterations limit is reached. We only allow balancing operations
that lead to improvement in the global loss $L$. Balancing a class can change
the distribution of other classes. Thus it may be possible to balance later a
class that is skipped in the first epoch.

The time complexity of the Algorithm 1 consists of processing $m$ classes, for
each calculating the loss (complexity: $\mathcal{O}\left(nm\right)$, finding
the class with the highest loss ($\mathcal{O}\left(m\right)$) and
redistributing the elements ($\mathcal{O}\left(n\right)$). Therefore, the
total time complexity is $\mathcal{O}\left(nm^{2}\right)$. In practice,
optisplit could be also easily used on top of another possibly faster method
to fine tune the results.

Note, that in the accompanying practical implementation the classes that have
more positive than negative data points are balanced with respect to the
negative distribution. This is not important for GO data but could be useful
in some other applications.

Input : measure $M$,

target data $D\in\\{0,1\\}^{n\times m}$,

number of cross validation folds $k$,

max number of epochs $\mathrm{max\\_epochs}$

Output : List of cross validation splits

Function  _ stratified_split(_D, k_)_

$\mathrm{folds}_{0}\leftarrow\textnormal{{make\\_random\\_split}}(D,k)$

$\mathrm{epoch}\leftarrow 0$

while  _true_ do

$\mathrm{max\\_offset}\leftarrow 0$

$L_{0}\leftarrow M(D,\mathrm{folds}_{0})$

for  _$i\in[m]$_ do

$A\leftarrow\operatorname*{arg\,max}{L_{0}}[\mathrm{max\\_offset}]$

$\mathrm{folds}_{1}\leftarrow\textnormal{{balance}}(A)$

$L_{1}\leftarrow M(D,\mathrm{folds}_{1})$

if  _$L_{0}\leq L_{1}$_ then

$\mathrm{max\\_offset}\leftarrow\mathrm{max\\_offset}+1$

continue

end if

$L_{0}\leftarrow L_{1}$

$\mathrm{folds}_{0}\leftarrow\mathrm{folds}_{1}$

end for

if  _$L_{0}=l_{1}\textnormal{{ or }}\mathrm{epoch}=\mathrm{max\\_epochs}$_
then

break

end if

end while

return _$\mathrm{folds}_{0}$_

Algorithm 1 Generating stratified cross validation splits

##  6 Cross validation experiments

In Section 4, we showed that LD is not a good measure for evaluating cross
validation split quality. Most existing multilabel stratified cross validation
methods found in the literature (see Section 3) are optimised for producing
splits with good LD or similar class size dependent scores. In this Section,
we will examine the performance of the existing methods, namely, SOIS, IS, SS
and PMBSRS as well as our own optisplit with respect to the proposed new
measures rLD and DCP. We will show that if one wants to get cross validation
splits that are good with respect to rLD and DCP, optisplit is the best option
available since it can be used to directly optimise them.

In the following experiments we set $k=5$. All the results presented are
averages over 10 runs with different random initialisations. The experiments
were run using Python 3.6.8 on a machine with AMD opteron-6736 1.4GHz.
Implementations of optisplit and the experiments presented in this article are
available at https://github.com/xtixtixt/optisplit.

We used iterative stratification implementations from the popular skmultilearn
library [7] for first order and second order iterative stratification
experiments. The Stratified Sampling [5] implementation used was the one
provided with the article. Note that the SS implementation produces train test
splits not cross validation splits. In order to compare it to the rest of the
methods we have split the data by recursively splitting it to approximately
1/k sized test sets. Thus the method is run 5 times.

The methods compared here can be divided into three categories. IS and SOIS
are iterative stratification based methods, SS and optisplit are optimisation
based methods and PMBSRS is a random split based method.

We optimised optisplit with respect to both rLD and DCP to compare the effect
of the cost function to the outcome.

###  6.1 Datasets

We used a wide range of diverse datasets: bibtex, delicious and mediamill from
the MULAN dataset collection [10] that have been used to evaluate earlier
similar methods. Datasets CC (cellular component), MF (molecular function) and
BP (biological process) are our own GO subset datasets used in protein
function prediction (see [9, 12] for more info). These and are considerably
bigger and sparser than the MULAN datasets used here. The dataset Wiki10-31K
[3] is a large and very sparse extreme classification dataset. Here, classes
without any positive or negative data points are excluded. Detailed properties
of the datasets are presented in Table 1

Data | Data size | Labels | Density | Min | 25% | 50% | 75% | Max
---|---|---|---|---|---|---|---|---
bibtex | 7395 | 159 | 0.0151 | 51 | 61 | 82 | 130 | 1042
delicious | 16015 | 983 | 0.0193 | 21 | 58 | 105 | 258 | 6495
mediamill | 43907 | 101 | 0.0433 | 31 | 93 | 312 | 1263 | 33869
CC | 577424 | 1688 | 0.0077 | 5 | 66 | 225 | 891 | 577410
MF | 637552 | 3452 | 0.0028 | 11 | 61 | 150 | 498 | 637533
BP | 666349 | 11288 | 0.0028 | 4 | 41 | 123 | 493 | 666338
Wiki10-31K | 20762 | 30938 | 0.0006 | 2 | 2 | 3 | 6 | 16756
Table 1: Properties of the datasets used in evaluations. CC, MF and BP are GO subsets. Min and Max represent the minimum and maximum class sizes of the datasets. Columns 25%, 50% and 75% are the corresponding percentiles of the class sizes. Dataset | Method | ED | LD | DCP | rLD | Runtime (s)
---|---|---|---|---|---|---
bibtex | $\textsc{optisplit}_{\text{rLD}}$ | 27 | 0.0004 | 0.0073 | 0.0234 | 5
$\textsc{optisplit}_{\textsc{DCP}}$ | 38 | 0.0005 | 0.0068 | 0.0315 | 5
SOIS | 16 | 0.0005 | 0.0143 | 0.0425 | 5
IS | 17 | 0.0007 | 0.0206 | 0.0604 | 1
| SS | 57 | 0.0007 | 0.0173 | 0.0465 | 4
| Random | 0 | 0.0022 | 0.0564 | 0.1693 | 1
| PMBSRS | 2 | 0.0022 | 0.0558 | 0.1660 | 1
mediamill | $\textsc{optisplit}_{\text{rLD}}$ | 53 | 0.0005 | 0.0053 | 0.0187 | 10
$\textsc{optisplit}_{\textsc{DCP}}$ | 7 | 0.0005 | 0.0047 | 0.0176 | 11
SOIS | 1 | 0.0003 | 0.0205 | 0.0610 | 77
IS | 1 | 0.0008 | 0.0280 | 0.0854 | 50
| SS | 36 | 0.0009 | 0.0068 | 0.0231 | 31
| Random | 0 | 0.0019 | 0.0379 | 0.1142 | 1
| PMBSRS | 2 | 0.0019 | 0.0386 | 0.1150 | 1
delicious | $\textsc{optisplit}_{\text{rLD}}$ | 35 | 0.0010 | 0.0223 | 0.0666 | 75
$\textsc{optisplit}_{\textsc{DCP}}$ | 28 | 0.0010 | 0.0215 | 0.0772 | 76
SOIS | 16 | 0.0012 | 0.0458 | 0.1357 | 381
IS | 13 | 0.0013 | 0.0489 | 0.1461 | 11
| SS | 75 | 0.0007 | 0.0221 | 0.0625 | 40
| Random | 0 | 0.0015 | 0.0507 | 0.1515 | 1
| PMBSRS | 2 | 0.0015 | 0.0512 | 0.1525 | 1
CC | $\textsc{optisplit}_{\text{rLD}}$ | 193 | 8.1053 | 0.0065 | 0.0230 | 2762
$\textsc{optisplit}_{\textsc{DCP}}$ | 78 | 8.1195 | 0.0062 | 0.0248 | 2872
SOIS | 5 | 5.7623 | 0.0305 | 0.0894 | 204823
IS | 1 | 5.6075 | 0.0448 | 0.1320 | 97120
| SS | 182 | 8.8831 | 0.0133 | 0.0416 | 1118
| Random | 0 | 10.2802 | 0.0455 | 0.1342 | 1
| PMBSRS | 1 | 6.3606 | 0.0448 | 0.1332 | 17
MF | $\textsc{optisplit}_{\text{rLD}}$ | 607 | 3.2147 | 0.0064 | 0.0229 | 4980
$\textsc{optisplit}_{\textsc{DCP}}$ | 89 | 3.1782 | 0.0063 | 0.0252 | 5209
SOIS | 84 | 3.0536 | 0.0490 | 0.1450 | 53646
IS | 1 | 2.7503 | 0.0490 | 0.1451 | 20442
| SS | 656 | 4.7935 | 0.0129 | 0.0400 | 1018
| Random | 0 | 6.2255 | 0.0493 | 0.1465 | 1
| PMBSRS | 1 | 4.8213 | 0.0498 | 0.1480 | 17
BP | $\textsc{optisplit}_{\text{rLD}}$ | 612 | 2.1053 | 0.0161 | 0.0516 | 59436
$\textsc{optisplit}_{\textsc{DCP}}$ | 173 | 2.0951 | 0.0156 | 0.0576 | 52412
SS | 279 | 2.9140 | 0.0282 | 0.0857 | 2734
Random | 0 | 2.896 | 0.0568 | 0.1664 | 1
| PMBSRS | 0 | 2.4250 | 0.0567 | 0.1662 | 21
Wiki10-31K | $\textsc{optisplit}_{\text{rLD}}$ | 1678 | 0.0002 | 0.2579 | 0.7563 | 3065
$\textsc{optisplit}_{\textsc{DCP}}$ | 367 | 0.0002 | 0.2068 | 0.8033 | 3053
SS | error | error | error | error | error
Random | 0 | 0.0002 | 0.3008 | 0.9316 | 1
| PMBSRS | 2 | 0.0002 | 0.3013 | 0.9323 | 1
Table 2: Performances of the algorithms evaluated on diverse datasets. The
scores are means over 10 runs. Bold font highlights the best DCP and rLD
results for each dataset. Error marks case where the method failed to run.
Notice that we ran slower methods (SOIS and IS) only with smaller datasets.

###  6.2 Results

Scores of the split quality measures for all datasets and methods are
presented in Table 2 with the following exceptions: IS and SOIS results are
not presented for the biggest datasets BP and Wiki10-31K because their runtime
was prohibitively high. Wiki10-31K results are not presented for SS because
the implementation used produced an error when run on that particular data.
For comparison purposes we have also presented scores for random split
(Random).

The results show that optisplit performs better than previous methods with
respect to rLD and DCP scores when optimising with either of those. The
runtimes of optisplit are also competitive when compared to other top
performing methods. Generally, iterative stratification based methods perform
quite well but are unusably slow on bigger datasets. Random split based
methods are fast but produce poor quality folds compared to more advanced
methods. Optimisation based methods (optisplit and SS) usually give best
results and their runtimes are in the middle of iterative stratification based
and random split based methods.

We can see that DCP and rLD are very correlated, the ordering of the methods
is similar with respect to both measures and optimising optisplit with respect
to DCP produces nearly as good rLD results as optimising rLD directly.
However, since rLD measures the folds more thoroughly i.e. it does not just
concentrate on the biggest fold it seems to be a better practical choice than
DCP in most cases.

Note that optisplit does not attempt to produce exactly equally sized splits.
This results in quite high ED (Example Distribution) scores compared to some
other methods. This should not be a problem in practical machine learning
settings.

For completeness, we have included LD evaluations in Table 2. As is to be
expected from the results presented in Section 4, we can see that the method
ordering is often considerably different with respect to LD scores. In smaller
and less imbalanced datasets LD gives results more in line with rLD and DCP.
For bigger and more imbalanced datasets, when LD weakness gets more
pronounced, the results differ more significantly. There, LD favours iterative
stratification based methods and gives random split based PMBSRS noticeably
better score than to Random, in contrast to rLD or DCP evaluations.

##  7 Discussion and future work

In this article, we have shown that the most widely used multilabel cross
validation split quality measure, LD, does not measure split quality correctly
when used on unequally sized classes. In response, we have presented new
measures with better properties and have presented a new general method,
optisplit, for generating and optimising multilabel stratified cross
validation splits. We have compared optisplit to existing methods and found
that it produces better quality cross validation folds with respect to the new
measures than the previous methods and scales well for GO sized datasets. We
note for future work that optisplit could be made faster by calculating the
loss only for the classes that have been modified in the previous balancing
operation. For sparse data that should allow it to be used even on
considerably larger datasets. Also, optisplit uses now a greedy hill-climbing
approach for optimising the target function. However, a Monte Carlo /
simulated annealing based version could achieve even better performance.

## Acknowledgements

This work was funded by NNF20OC0065157 of the Novo Nordisk Foundation.
Computations were partly done using resources belonging to Biocenter Finland’s
Bioinformatics platform.

## References

  * [1] M. Ashburner, C. A. Ball, J. A. Blake, D. Botstein, H. Butler, J M. C., A. P. Davis, Kara Dolinski, S. S. Dwight, J. T. Eppig, et al.  Gene ontology: tool for the unification of biology.  Nature genetics, 25(1):25–29, 2000.
  * [2] S. Bengio, K Dembczynski, T. Joachims, M. Kloft, and M. Varma.  Extreme Classification (Dagstuhl Seminar 18291).  Dagstuhl Reports, 8(7):62–80, 2019.
  * [3] K. Bhatia, K. Dahiya, H. Jain, P. Kar, A. Mittal, Y. Prabhu, and M. Varma.  The extreme classification repository: Multi-label datasets and code, 2016\.
  * [4] F. Charte, A. Rivera, M. J. del Jesus, F. Herrera, A. Troncoso, H. Quintián, and E. Corchado.  On the impact of dataset complexity and sampling strategy in multilabel classifiers performance.  In Hybrid Artificial Intelligent Systems, pages 500–511, Cham, 2016\. Springer International Publishing.
  * [5] M Merrillees and L Du.  Stratified Sampling for Extreme Multi-Label Data.  arXiv e-prints, page arXiv:2103.03494, March 2021.
  * [6] K. Sechidis, G. Tsoumakas, and I. Vlahavas.  On the stratification of multi-label data.  In Machine Learning and Knowledge Discovery in Databases, pages 145–158, Berlin, Heidelberg, 2011. Springer Berlin Heidelberg.
  * [7] P. Szymański and T. Kajdanowicz.  A scikit-based Python environment for performing multi-label classification.  ArXiv e-prints, February 2017.
  * [8] P. Szymański and T. Kajdanowicz.  A network perspective on stratification of multi-label data.  In Proceedings of the First International Workshop on Learning with Imbalanced Domains: Theory and Applications, volume 74 of Proceedings of Machine Learning Research, pages 22–35. PMLR, 2017.
  * [9] P. Törönen, A. Medlar, and L. Holm.  Pannzer2: a rapid functional annotation web server.  Nucleic acids research, 46(W1):W84–W88, 2018.
  * [10] G. Tsoumakas, E. Spyromitros-Xioufis, J. Vilcek, and I. Vlahavas.  Mulan: A java library for multi-label learning.  Journal of Machine Learning Research, 12:2411–2414, 2011.
  * [11] David H. Wolpert.  Stacked generalization.  Neural Networks, 5(2):241–259, 1992.
  * [12] N. Zhou, Y. Jiang, T. R. Bergquist, A. J. Lee, B. Z. Kacsoh, A. W. Crocker, K. A. Lewis, G. Georghiou, H. N. Nguyen, M. N. Hamid, et al.  The cafa challenge reports improved protein function prediction and new functional annotations for hundreds of genes through experimental screens.  Genome biology, 20(1):1–23, 2019.


# Statistical Consequences of Dueling Bandits

Nayan Saxena
Department of Statistical Sciences
University of Toronto
<EMAIL_ADDRESS>
&Pan Chen
Department of Computer Science
University of Toronto
<EMAIL_ADDRESS>
&Emmy Liu
Department of Computer Science,
Cognitive Science Program
University of Toronto
<EMAIL_ADDRESS>

###### Abstract

Multi-Armed-Bandit frameworks have often been used by researchers to assess
educational interventions, however, recent work has shown that it is more
beneficial for a student to provide qualitative feedback through preference
elicitation between different alternatives, making a dueling bandits framework
more appropriate. In this paper, we explore the statistical quality of data
under this framework by comparing traditional uniform sampling to a dueling
bandit algorithm and find that dueling bandit algorithms perform well at
cumulative regret minimisation, but lead to inflated Type-I error rates and
reduced power under certain circumstances. Through these results we provide
insight into the challenges and opportunities in using dueling bandit
algorithms to run adaptive experiments.

##  1 Introduction

Rapid growth in online learning has provided scientists and researchers a new
digital platform to conduct randomised experiments with students in real-world
settings. This data can be leveraged by education researchers to adaptively
assign students carefully tailored educational material, and further improve
the quality of their content through careful exploratory assessment of the
assigned conditions [19, 23]. To achieve these goals, one of the most popular
problem frameworks is the Multi-Armed Bandit (MAB) problem which focuses on
optimal assignment of conditions based on numerical reward signals from the
participant (student) whose goal is to trade-off between exploration and
exploitation of conditions [3, 20]. Historically, MAB algorithms have been
used in industrial settings to leverage user feedback and adaptively display
more popular advertisements, websites and produce content [8, 13], whereas in
research settings these algorithms have been used for clinical trials [14, 21,
4], robot control [15, 16], and by behavioral and social scientists for crowd-
sourcing experiments [1, 5, 2].

While MAB algorithms focus entirely on quantitative feedback (scalar reward),
it has been shown that human motor learning is maximised when subjected to
both qualitative and quantitative feedback [12] and recent work suggests that
quantitative metrics might not be the best indicators of true human
preferences [18]. Therefore, a more qualitative approach for preference
elicitation from human participants should be more ideal to leverage student
feedback, which can in turn be achieved through preference based MAB (dueling
bandit) algorithms like Double Thompson Sampling (DTS) [7, 24], which focuses
on participants simply choosing between two presented alternatives.

Through this paper, we take a closer look at leveraging these dueling bandit
algorithms for adaptive experimentation by empirically investigating the
statistical properties of the DTS algorithm against uniform random
experimental condition selection. Our study is primarily motivated by prior
work assessing the statistical properties of MAB algorithms in an educational
setting [19], and builds upon recent work exploring the challenges that come
with conducting hypothesis tests to analyze data from adaptive experiments
using bandit algorithms [22], such as proposing strategies for modifying MAB
algorithms to trade-off reward and power [25] and improving coverage of
confidence intervals [10].

In summary, our main contributions are:

  1. 1.

We emphasise the conceptual significance of quantifying the differences in
statistical power, regret and false positive rate between uniform sampling and
Double Thompson Sampling

  2. 2.

Through simulation experiments we show that for a reasonable decrease in
power, it may be advantageous to run experiments with dueling bandits when the
number of arms is small and the expected effect size is large

  3. 3.

By applying the same analysis on the real-world Microsoft Learning to Rank
(LTR) dataset we confirm that using DTS results in lowered power, regret and a
higher proportion of people being assigned to the better alternative.

##  2 Preliminaries

###  2.1 Multi-armed Bandits

As a precursor to dueling bandits, we consider a central problem in sequential
decision making and adaptive experiment design known as the multi-armed bandit
(MAB) problem [20]. The problem setup consists of $k$ arms (actions) such that
for every action $a_{i}\in\mathcal{A}$ the corresponding probability of
yielding a success is $p_{i}\in[0,1]$. With no a priori information about
these success probabilities, the central goal for an agent is to maximise the
overall number of successes by pulling these arms, or performing these
actions, and ultimately settling on one over a fixed period of time denoted by
$t$. Assume that under this framework the expected reward of arm $a_{i}$ is
given by $\mu(a_{i})$ and $\mu^{*}=\operatorname*{arg\,max}\mu(a_{i})$, for
$i\in\\{1\ldots k\\}$. The objective of an agent is to then minimise the
cumulative regret,

$\mathcal{R}^{\text{\tiny MAB}}_{t}=\sum_{i=1}^{t}[\mu^{*}-\mu_{a_{(i)}}]$

Here the agent is allowed to only choose one action, $a_{(i)}$ at a given time
step. However, in many cases during adaptive experimentation, it may be
difficult to frame the result of an action as a success or failure. In
particular, in many cases what experimenters are interested in is simply which
arm is preferred to the others, in which case pairwise preferences may be more
appropriate for the subject to make. To accommodate pairwise preferences we
now consider a generalised form of the MAB problem also known as the dueling
bandit problem.

###  2.2 Dueling Bandits

The dueling bandit problem can be characterised as a special case of the
popular Multi-armed Bandit (MAB) problem that focuses on pairwise comparisons
between actions at every iteration [26]. The problem setup is similar to the
MAB problem, except at every iteration the agent chooses two actions
$a_{m},a_{n}\in\mathcal{A}$ and performs a comparison before choosing one
action that they prefer 111The traditional dueling bandit framework allows
comparisons between the same actions (self-dueling) where $m=n$, but
throughout this paper we assume $m\neq n$ as it is unlikely for users to be
presented with comparisons between the same actions during adaptive
experimentation .. Throughout this paper we specifically utilise the Double
Thompson Sampling (DTS) algorithm, which is an adaptation of the popular
Thompson Sampling algorithm in the regular MAB context [24]. A comprehensive
outline of the DTS procedure can be found in Appendix A.

Formally, the probability of one arm winning over another is given by
${P}(a_{m}\succ a_{n})$ which we abbreviate as $P_{mn}\equiv{P}(a_{m}\succ
a_{n})$. For each trial $i$ the outcome of each comparison is binary
$x_{i}\sim\text{Bernoulli}(P_{mn})$ where the probability of one arm winning
is formally given by,

${P}(a_{m}\succ a_{n})=\Delta(a_{m},a_{n})+0.5$

Here, $\Delta(a_{m},a_{n})\in[-0.5,0.5]$ is the difference measure between the
two actions, abbreviated as $\Delta_{mn}\equiv\Delta(a_{m},a_{n})$, such that
$a_{m}\succ a_{n}\Longleftrightarrow\Delta_{mn}>0$. It is further assumed that
the probabilities representing user preferences are stationary over time and
modelled by the unknown preference matrix denoted by
$\mathbf{P}=[P_{mn}]_{k\times k}$ where the entries satisfy $P_{mn}+P_{nm}=1$.

The goal of an agent operating in this environment is to choose the optimal
winning arm when conducting pairwise comparisons and ultimately minimise
regret. Motivated by classical voting theory, there are different notions of a
winning arm like the Borda winner [11] and Copeland winner [27, 9]. In this
paper we focus on using the notion of a Condorcet winner [6] and strong regret
[24] as described below.

##### Condorcet Winner

The Condorcet winner is a single action $a^{*}$ that beats all other actions
in a pairwise comparison, such that $\Delta_{a^{*}a}>0$ for all $a\neq a^{*}$.
The Copeland winner, found by maximising the normalised Copeland score
$\frac{1}{k-1}\sum_{m\neq n}\mathbbm{1}[P_{mn}>0.5]$ is a Condorcet criterion,
which means that it always finds the Condorcet winner if one exists. It should
be noted that it is possible for no Condorcet winner to exist, in which case
the same definition can be extended to obtain a set of Copeland winners which
always exist. [27, 24].

##### Strong Regret

Suppose $a^{(i)}_{m}$ and $a^{(i)}_{n}$ are actions chosen at a given timestep
$i$ and $a_{o}$ is the optimal action. Then under this framework the overall
objective for an agent is to minimise the cumulative strong regret,

$\mathcal{R}^{\text{\tiny DUEL}}_{t}=\sum_{i=1}^{t}[\Delta_{om}+\Delta_{on}]$

##  3 Methods

###  3.1 Experimental Setup

In order to examine the impact of dueling bandit algorithms as compared to
uniform random assignment during adaptive experimentation, we looked at
several conditions relevant to experimenters by varying the number of arms and
the effect sizes between pairs of arms. For each set of conditions 5000
simulations were carried out, each corresponding to a real-world study, where
data at each timestep would correspond to data collected from a participant
exposed to a treatment arm. Furthermore, to assess the statistical quality of
each algorithm, we considered the statistical power in comparing pairs of
arms, the regret accumulated, and the false positive rate. During its
operation, each simulation assigned hypothetical participants using either the
DTS algorithm, or uniformly at random and was initialized using a preference
matrix $P\in\mathbbm{R}^{n\times n}$, where $n$ is the number of arms. For
each given pair of arms $(i,j)$, the effect size was calculated (Cohen’s w of
0.1, 0.3 or 0.5) of the comparison based on the difference (0.05, 0.15, or
0.25) between $P_{ij}$ and $P_{ji}$, noting that $P_{ij}+P_{ji}=1$. Finally,
in order to compute the false positive rate, preference matrices with zero
effect size for all comparisons between pairs of arms were also used. These
experimental conditions are also summarized in Table 1. Data and code are
publicly available at https://bit.ly/35J2xH0.

###  3.2 Experimental Dataset

In addition to simulated data, we also present experiments using data
collected from users in a real-world context. In particular, we perform
simulations using the Microsoft Learning to Rank (LTR) dataset which presents
pairs of search queries and documents, along with the relevance ratings of the
document [17]. As A/B testing is commonly used to test ranking and
recommendation algorithms, this is an appropriate context in which we
benchmark outcomes of such experiments. We used an implicit preference matrix
over all 136 rankers derived from this dataset similar to [27], and then
randomly sampled Condorcet and non-Condorcet submatrices from this larger
matrix in order to simulate experiments comparing different rankers. The
sample size is set as ${n\choose 2}\times m$, where arm $i$ and arm $j$ are
the pair of arms with the smallest effect size, and $m$ is the number of
participants needed to achieve expected statistical power of 0.8 for effect
size 0.1. This simulates the real-world scenario where the number of required
participants is unknown. To further evaluate the long-term performance of the
two assignment methods, we ran the simulations with sample sizes up to 10
times more than the initial sample size and tracked the same set of metrics as
in the main analysis.

Condition |  | Description
---|---|---
Sampling type |  | Double Thompson Sampling
|  | Uniform Random Sampling
Number of arms |  | $n\in\\{3,5\\}$
Sample sizes | Uniform effect sizes | Non-zero effect sizes: $s=10mn$ where $m$ is the number of participants
|  | needed to achieve expected statistical power of 0.8 and $n$ is
|  | the number of paired comparisons in that simulation
| Zero effect size | Same sample sizes as in the non-zero effect case calculated using the
|  | number of participants needed to achieve expected statistical power of
|  | 0.8 for effect size 0.3.
| Learning to Rank dataset | Different effect sizes: $s=10mn$, where $m$ is the number of participants
|  | needed to achieve expected statistical power of 0.8 for 0.1 effect size,
|  | and n is the number of paired comparisons in that simulation.
Effect sizes | Simulated dataset | None (0), Small (0.1), Medium (0.3), Large (0.5). Same effect size between
|  | each pair of arms. The winning arm in each pair is randomly assigned
|  | in non-zero cases.
| Learning to Rank dataset | For each pair of arms $i$ and $j$, the effect size is in $[0,1)$.

Table 1: Summary of conditions varied across simulations

###  3.3 Analysis

After each simulation, we shift our focus to analyzing several outcomes
relevant to experimenters: (a) statistical power to detect an effect between
each pair of arms, (b) false positive rate in detecting effects when none
exist, (c) regret over time, and (d) percentage of participants assigned to
the Condorcet winner when one exists. For simulations where there was some
difference between pairs of arms, a Chi-squared contingency test was performed
for each pair of arms, with significance level $\alpha=0.05$, which is the
standard across multiple domains including educational experiments. For each
simulation, regret over time, reward over time, and percentage of participants
assigned to the Condorcet winner was recorded and aggregated. Finally, to
calculate the false positive rate, we considered a pair of arms as a false
positive if the comparison between them reached significance, given there
should be no effect between these two arms.

##  4 Results

###  4.1 Synthetic Data

####  4.1.1 Conditions that differ in terms of effectiveness

##### Average Power Over Time

In educational settings, different conditions, such as sample solutions, may
have a different impact on students’ learning efficiency and engagement with
the educational resources. However, one condition that outperforms the others
is often not statistically significant until testing is conducted with many
students involved. Yet, different ways of conducting the experiments may
require a different number of participants to reach significant effects
amongst conditions.

Figure 1: Average power over time for DTS algorithm (blue) and Uniform
assignment (red) with synthetic data.

We modeled the means of the power of different pairs of conditions, as this
measures how well each assignment did in finding statistically significant
effects among all the conditions. It can be observed in Figure 1 that uniform
sampling consistently reached the 0.8 power threshold with fewer participants
as compared to DTS, in both 3 and 5-condition settings. This pattern was also
observed when final power was recorded for both uniform sampling and DTS, as
seen in Figure 2, thus indicating that DTS is more susceptible to inflated
Type-II error rates (reduced power). This shows that when there are
insufficient students to test among multiple variants, uniform sampling is
more often a better way of assignment given the effects among different
conditions are of the highest priority.

Figure 2: Comparison of final power between DTS algorithm (blue) and Uniform
assignment (red) for simulated dataset where the dotted line indicates effect
detected.

##### Proportion of Condorcet Winners

In educational experiments, given that there is one unknown best condition, we
examined how uniform sampling and double Thompson sampling fulfill the goal of
providing students with the best condition by calculating the proportion of
trials in which a student would be assigned to the Condorcet winner.

Figure 3: Proportion of trials with Condorcet winners for DTS algorithm (blue)
and Uniform assignment (red) with synthetic data. Here, the dotted line
indicates expected assignment.

From Figure 3, it can be observed that, on average, DTS assigned a higher
number of Condorcet winners to students who participated in the experiments.
Across both 3 and 5-condition simulations, we observe that DTS assigned the
Condorcet winner to almost every student, whereas with uniform sampling the
number of students receiving the Condorcet winner is in line with the
statistical expectation.

##### Cumulative Strong Regret

Apart from the proportion of Condorcet winners, we used cumulative strong
regret to measure the different experiences uniform sampling and DTS give to
students. As seen in Figure 4, a noticeable difference in cumulative strong
regret between uniform sampling and DTS can be seen, which grew more evident
with increase in effect size from 0.5 to 0.1, and number of arms from 3 to 5.
Across all simulations, DTS performed better by accumulating less regret as
compared to uniform sampling.

Figure 4: Comparison of cumulative strong regret over time between DTS
algorithm (blue) and Uniform assignment (red) for simulated dataset.

####  4.1.2 Conditions that are equally-effective

##### False Positive Rate

To simulate scientific settings like clinical trials, along with behavioral
and social sciences where little to no difference between arms is commonly
observed, we measured the false positive rate when conditions are equally
effective. False positives here mean that at least one comparison between arms
produces a significance value $<0.05$ while the effect size between each pair
of arms is 0. From Figure 5, we observe that the false positive rates for
uniform sampling remain consistent across 3-condition and 5-condition
scenarios and are relatively lower when compared to DTS algorithm.

Figure 5: False positive rate for DTS algorithm (blue) and Uniform assignment
(red) with simulated data.

This indicates that when there is little to no difference between conditions
the DTS algorithm will result in a higher Type-I error as compared to uniform
sampling. These issues are especially problematic for scientific research,
since erroneously believing an intervention is better than a control may
result in lack of reproducibility.

###  4.2 Learning to Rank Dataset

In real settings the effect sizes between any two conditions are not limited
to the set $\\{0.1,0.3,0.5\\}$, therefore we ran the same experiments on the
LTR dataset. Overall, our results on the real-world dataset were consistent
with our simulation studies. We noticed that DTS provided the best condition
to students more often and accrued relatively lower average power as seen in
Figure 6. It should further be noted that in a few trials DTS did not assign
the Condorcet winner to the majority of students possibly because the
differences between the best condition and other arms is miniscule, thus
requiring more participants for the algorithm to learn the best condition.

Figure 6: Proportion of trials with Condorcet winners alongside average power
over time for DTS algorithm (blue) and Uniform assignment (red) with LTR
dataset.

Consistent with our results on synthetic data, it can further be observed in
Figure 6 that uniform sampling reached the 0.8 power threshold with fewer
participants as compared to DTS across both 3 and 5-condition settings. This
pattern was also observed when final power was recorded for both uniform
sampling and DTS, as seen in Figure 7, thus reinforcing the finding that DTS
is more susceptible to reduced power even in real world settings.

Figure 7: Comparison of final power between DTS algorithm (blue) and Uniform
assignment (red) for LTR dataset where the dotted line indicates effect
detected.

Figure 8: Comparison of cumulative strong regret over time between DTS
algorithm (blue) and Uniform assignment (red) for LTR dataset

Finally, we also observed that DTS does a much better job at minimising regret
than uniform sampling on the LTR data set as seen in Figure 8. Overall, we
note that both uniform sampling and DTS are useful algorithms in educational
contexts depending on the available resources and potential intervention
goals. Uniform sampling is more appropriate in cases where finding significant
effects is of higher priority and the number of students is small, whereas DTS
is more advantageous when better student experience is more important through
optimal allocation.

##  5 Conclusion

Adaptive experimentation through bandits presents an important opportunity for
researchers, as it allows more participants to receive the probable best
intervention, compared to uniform sampling. Although the benefits and
drawbacks of experimentation with the Thompson sampling algorithm have been
explored, in some cases a dueling bandit paradigm may be more appropriate. To
explore the statistical quality of data during randomised experimentation
under this framework we conducted an analysis of statistical power, regret,
and false positive rates comparing traditional uniform sampling to the DTS
algorithm. Our findings indicate that the DTS algorithm shows excellent
performance with respect to cumulative regret minimisation but in the case
when there is little to no difference between arms, DTS leads to higher Type-I
error rates and reduced power. This makes it increasingly difficult to collect
data that enables statistical hypothesis tests between conditions, while also
balancing this with reward maximization and minimizing the false positive
rate. These findings pave the way for researchers to better leverage machine
learning for adaptive experimentation, while enabling sound statistical
inference. They also present a new opportunity for future work to explore ways
of incorporating the DTS algorithm in adaptive experiments where both reward
optimisation and statistical inference are important, possibly through
algorithm-driven dynamic interpolation between DTS and uniform exploration for
better performance.

## References

  * Abraham et al. [2013] I. Abraham, O. Alonso, V. Kandylas, and A. Slivkins.  Adaptive crowdsourcing algorithms for the bandit survey problem.  In _Conference on learning theory_ , pages 882–910. PMLR, 2013.
  * Athey et al. [2017] S. Athey, S. Wager, et al.  Efficient policy learning.  _arXiv preprint arXiv:1702.02896_ , 78, 2017.
  * Auer [2002] P. Auer.  Using confidence bounds for exploitation-exploration trade-offs.  _Journal of Machine Learning Research_ , 3(Nov):397–422, 2002.
  * Bather [1980] J. Bather.  Randomised allocation of treatments in sequential trials.  _Advances in Applied Probability_ , pages 174–182, 1980.
  * Berkowitz et al. [2018] J. P. Berkowitz, C. A. Depken II, and J. M. Gandar.  One-armed bandits and behavioral economics.  In _Handbook of Behavioral Industrial Organization_. Edward Elgar Publishing, 2018.
  * Black et al. [1958] D. Black et al.  The theory of committees and elections.  1958\.
  * Busa-Fekete and Hüllermeier [2014] R. Busa-Fekete and E. Hüllermeier.  A survey of preference-based online learning with bandit algorithms.  In _International Conference on Algorithmic Learning Theory_ , pages 18–39. Springer, 2014.
  * Chakrabarti et al. [2008] D. Chakrabarti, R. Kumar, F. Radlinski, and E. Upfal.  Mortal multi-armed bandits.  _Advances in neural information processing systems_ , 21:273–280, 2008.
  * Copeland [1951] A. Copeland.  A ‘reasonable’social welfare function, seminar on mathematics in social sciences, university of michigan.  _Cited indirectly from its mention by Luce and Raiffa (1957)_ , page 358, 1951.
  * Deshpande et al. [2018] Y. Deshpande, L. Mackey, V. Syrgkanis, and M. Taddy.  Accurate inference for adaptive linear models.  In _International Conference on Machine Learning_ , pages 1194–1203. PMLR, 2018.
  * Jamieson et al. [2015] K. Jamieson, S. Katariya, A. Deshpande, and R. Nowak.  Sparse dueling bandits.  In _Artificial Intelligence and Statistics_ , pages 416–424. PMLR, 2015.
  * Kilduski and Rice [2003] N. C. Kilduski and M. S. Rice.  Qualitative and quantitative knowledge of results: effects on motor learning.  _American Journal of Occupational Therapy_ , 57(3):329–336, 2003.
  * Kohli et al. [2013] P. Kohli, M. Salek, and G. Stoddard.  A fast bandit algorithm for recommendation to users with heterogenous tastes.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 27, 2013.
  * Kuleshov and Precup [2014] V. Kuleshov and D. Precup.  Algorithms for multi-armed bandit problems.  _arXiv preprint arXiv:1402.6028_ , 2014.
  * Matikainen et al. [2013] P. Matikainen, P. M. Furlong, R. Sukthankar, and M. Hebert.  Multi-armed recommendation bandits for selecting state machine policies for robotic systems.  In _2013 IEEE International Conference on Robotics and Automation_ , pages 4545–4551. IEEE, 2013.
  * McGuire et al. [2018] S. McGuire, P. M. Furlong, C. Heckman, S. Julier, D. Szafir, and N. Ahmed.  Failure is not an option: Policy learning for adaptive recovery in space operations.  _IEEE Robotics and Automation Letters_ , 3(3):1639–1646, 2018.
  * Qin and Liu [2013] T. Qin and T. Liu.  Introducing LETOR 4.0 datasets.  _CoRR_ , abs/1306.2597, 2013.  URL http://arxiv.org/abs/1306.2597.
  * Radlinski et al. [2008] F. Radlinski, M. Kurup, and T. Joachims.  How does clickthrough data reflect retrieval quality?  In _Proceedings of the 17th ACM conference on Information and knowledge management_ , pages 43–52, 2008.
  * Rafferty et al. [2019] A. Rafferty, H. Ying, and J. Williams.  Statistical consequences of using multi-armed bandits to conduct adaptive educational experiments.  _JEDM| Journal of Educational Data Mining_ , 11(1):47–79, 2019.
  * Robbins [1952] H. Robbins.  Some aspects of the sequential design of experiments.  _Bulletin of the American Mathematical Society_ , 58(5):527–535, 1952.
  * Sui and Burdick [2014] Y. Sui and J. Burdick.  Clinical online recommendation with subgroup rank feedback.  In _Proceedings of the 8th ACM conference on recommender systems_ , pages 289–292, 2014.
  * Villar et al. [2015] S. S. Villar, J. Bowden, and J. Wason.  Multi-armed bandit models for the optimal design of clinical trials: benefits and challenges.  _Statistical science: a review journal of the Institute of Mathematical Statistics_ , 30(2):199, 2015.
  * Williams et al. [2018] J. J. Williams, A. N. Rafferty, D. Tingley, A. Ang, W. S. Lasecki, and J. Kim.  Enhancing online problems through instructor-centered tools for randomized experiments.  In _Proceedings of the 2018 CHI Conference on Human Factors in Computing Systems_ , pages 1–12, 2018.
  * Wu and Liu [2016] H. Wu and X. Liu.  Double thompson sampling for dueling bandits.  _arXiv preprint arXiv:1604.07101_ , 2016.
  * Yao et al. [2020] J. Yao, E. Brunskill, W. Pan, S. Murphy, and F. Doshi-Velez.  Power-constrained bandits.  _arXiv preprint arXiv:2004.06230_ , 2020.
  * Yue et al. [2012] Y. Yue, J. Broder, R. Kleinberg, and T. Joachims.  The k-armed dueling bandits problem.  _Journal of Computer and System Sciences_ , 78(5):1538–1556, 2012.
  * Zoghi et al. [2015] M. Zoghi, Z. Karnin, S. Whiteson, and M. De Rijke.  Copeland dueling bandits.  _arXiv preprint arXiv:1506.00312_ , 2015.

## Appendix

## A. Double Thompson Sampling


# Extending the capabilities of vectorial ptychography to circular-polarizing
materials such as cholesteric liquid crystals

Patrick Ferrand  Aix-Marseille Univ, CNRS, Centrale Méditerranée, Institut
Fresnel, F-13013 Marseille, France<EMAIL_ADDRESS>Michel Mitov
Univ. Côte d’Azur, CNRS, Institut de Physique de Nice, INPHYNI, F-06200 Nice,
France

###### Abstract

The problem of imaging materials with circular polarization properties is
discussed within the framework of vectorial ptychography. We demonstrate, both
theoretically and numerically, that using linear polarizations to investigate
such materials compromises the unicity of the solution provided by this
computational method. To overcome this limitation, an improved measurement
approach is proposed, which involves specific combinations of elliptical
polarizations. The effectiveness of this strategy is demonstrated by numerical
simulations and experimental measurements on cholesteric liquid crystals
films, which possess unique polarization properties. With the help of Pauli
matrices algebra, our results highlight the technique’s ability to discern
between different types of circular polarizers, uniform vs. non-uniform, and
determine their handedness.

††journal: opticajournal

##  1 Introduction

Vectorial ptychography is a recent imaging technique that uses phase retrieval
algorithms to provide quantitative maps of Jones matrices [1]. It is used in
optical microscopy to study specimens that strongly affect the phase and
polarization of transmitted light, and it has a robust reference-free
experimental scheme [2]. This variant of optical ptychography [3] has been
successfully applied to various materials, ranging from natural substances
like biomineral calcareous shells [4, 5] or biological tissues [6] to advanced
optical components such as holographic polarization-controlled metasurfaces
[7]. In previous studies, vectorial ptychography has relied on linear
polarization states for both illumination and detection [2]. While this
approach has been effective in addressing various challenging situations, it
may have limitations when dealing with materials that exhibit strong circular-
polarization properties. Such materials are often encountered in chiral
molecular assemblies, which can make reconstruction difficult.

This letter presents an extension to the capabilities of vectorial
ptychography, demonstrating how it can be applied to circular-polarizing
materials. Initially, we provide theoretical evidence of the
underdetermination introduced by linear polarizations when investigating such
materials. Subsequently, we propose an enhanced measurement scheme that relies
on combinations of elliptical polarizations. We validate this approach through
numerical simulations, where the results are analyzed with the help of Pauli
matrices algebra. Finally, experimental vectorial ptychography measurements
are conducted on cholesteric liquid crystal (CLC) films, which exhibit
specific polarization properties.

##  2 Theory

###  2.1 Principle of vectorial ptychography

In vectorial ptychography, the recorded intensity at the $j$-th scanning
position is represented as the square modulus of the far field

$I_{jkl}(\mathbf{q})=\left|\mathcal{F}\left(\psi_{jkl}(\mathbf{r})\right)\right|^{2}$
(1)

where where $\psi_{jkl}$ denotes the scalar exit field after analysis for the
$k$-th polarized probe and $l$-th polarization analysis [1]. The operator
$\mathcal{F}$ represents a propagation operator, while $\mathbf{r}$ and
$\mathbf{q}$ are the direct and reciprocal space coordinates, respectively.
The exit field can be expressed as

$\psi_{jkl}(\mathbf{r})=\mathbf{h}_{l}^{\mathrm{t}}\mathbf{J}(\mathbf{r}-\mathbf{r}_{j})\mathbf{p}_{k}(\mathbf{r})$
(2)

where $\mathbf{h}_{l}$ is the polarization analysis operator, "t" denotes the
transpose operator, $\mathbf{J}(\mathbf{r}-\mathbf{r}_{j})$ represents the
Jones matrix map of the laterally shifted investigated object, and
$\mathbf{p}_{k}(\mathbf{r})$ corresponds to the vectorial field distribution
of the polarized illumination probe.

The success of the iterative ptychography reconstruction depends on the
changes observed in the recorded intensity patterns $I_{jkl}(\mathbf{q})$
during spatial ($j$) and polarization ($k,l$) scanning [2]. Previous works
have mainly focused on linear polarizations for both the illumination (at
angle $\alpha_{k}$) and analysis (at angle $\theta_{l}$), given by

$\mathbf{p}^{\mathrm{lin}}_{k}\propto\begin{bmatrix}\cos\alpha_{k}\\\
\sin\alpha_{k}\end{bmatrix}\quad\textrm{and}\quad\mathbf{h}^{\mathrm{lin}}_{l}\propto\begin{bmatrix}\cos\theta_{l}\\\
\sin\theta_{l}\end{bmatrix}.$ (3)

While this measurement scheme has been employed to investigate a wide range of
optical properties [4], the specific case of circular-polarizing materials has
not been studied.

###  2.2 Circular-polarizing materials

Circular-polarizing materials possess the property of transforming incident
light into circularly polarized light, in transmission and/or in reflection.
Additionally, if a circular polarization of a particular handedness (left or
right) can be transmitted with the same handedness, the circular polarizer is
referred to as homogeneous [8]. In this study, it is important to note that we
will adopt a convention for which the term "left" refers to a polarized field
that circulates in a counterclockwise direction when observed from the
detector’s viewpoint.

###  2.3 Linear polarization scheme

For the sake of simplicity, let us consider a material that behaves as a
homogeneous left-circular polarizer in transmission, along with transmittance
properties $T(\mathbf{r})$ that are independent of polarization. This material
can be described by the following Jones matrix,

$\mathbf{J}(\mathbf{r})=\frac{1}{2}T(\mathbf{r})\begin{bmatrix}1&-i\\\
+i&1\end{bmatrix}.$ (4)

Under a linear polarization scheme (Eq. 3), it can be demonstrated that the
exit field, as given by Eq. 2, becomes

$\psi^{\mathrm{lin}}_{jkl}(\mathbf{r})\propto
T(\mathbf{r}-\mathbf{r}_{j})e^{i(\theta_{l}-\alpha_{k})}.$ (5)

It is important to note that the polarization properties contribute to the
exit field solely as a spatially homogeneous phase factor in Eq. 5.
Consequently, the recorded intensity remains unchanged according to Eq. 1,
regardless of the combination of polarizations employed. In a more general
context, it can be shown that different types of circular-polarizing
materials, whether left or right, homogeneous or inhomogeneous, would yield
the same set of diffracted intensities $I_{jkl}(\mathbf{q})$ and, therefore,
would be indistinguishable using this measurement scheme, compromising the
unicity of the ptychographic reconstruction.

###  2.4 Improved polarization scheme

The effective resolution of this ambiguity can be achieved by utilizing a
wider range of polarizations. By substituting the operators $\mathbf{p}_{k}$
and $\mathbf{h}_{l}$ with the following expressions

$\mathbf{p}^{\mathrm{ell}}_{k}\propto\begin{bmatrix}\cos\alpha_{k}\\\
-i\sin\alpha_{k}\end{bmatrix}\quad\textrm{and}\quad\mathbf{h}^{\mathrm{ell}}_{l}\propto\begin{bmatrix}\cos\theta_{l}\\\
-i\sin\theta_{l}\end{bmatrix},$ (6)

which correspond to general elliptical polarizations with an azimuth of 0° or
90° and ellipticities defined by $\tan\alpha_{k}$ and $\tan\theta_{l}$,
respectively, the resulting exit field described by Eq. 2, in the example
considered earlier, becomes

$\psi^{\mathrm{ell}}_{jkl}(\mathbf{r})\propto
T(\mathbf{r}-\mathbf{r}_{j})(\cos\alpha_{k}-\sin\alpha_{k})(\cos\theta_{l}+\sin\theta_{l}).$
(7)

As it will be illustrated later, the concept can be easily grasped by
considering the Poincaré sphere representation for polarizations. It involves
investigating multiple points on a sphere’s meridian for both illumination and
analysis, involving spherical, elliptical, and linear polarizations. This
situation allows us to overcome the barrier imposed by linear polarizations
positioned on the sphere’s equator. This modification introduces a clear
signature of the circular polarization properties, resulting in a change in
the amplitude of the exit field as given by Eq. 7, depending on the specific
combination of polarizations. Furthermore, it can be demonstrated that the
amplitude factors in Eq. 7 take different algebraic forms, facilitating the
unequivocal identification of the type of circular polarizer. The relevance of
this latter approach has been tested on both numerical and experimental
datasets.

##  3 Materials and methods

###  3.1 Numerical simulations

Numerical datasets were generated to simulate an experiment of vectorial
ptychography at a wavelength $\lambda=635$ nm using 50-µm-diameter cropped
Gaussian-shaped illumination probes under which an object is raster-scanned
with a step size of 7 µm in both directions, resulting in an $11\times 11$
grid. Far-field intensity patterns were computed as if they were captured at
an infinite distance, within a numerical aperture of 0.3, and recorded on a
camera sensor with dimensions of $122\times 122$ pixels. To replicate the
effects of shot noise on the sensor, a Poisson random number generator was
employed. Typically, each frame had a total photon count of $10^{6}$. The
datasets were processed using the vectorial ptychographic iterative algorithm
detailed in a previous publication [1]. The algorithm employed random
distributions (modulus and phase) for the initial guesses and was run for 30
iterations with known probes.

###  3.2 CLC films

Cholesteric liquid crystals (CLC) films were experimentally investigated at
$\lambda=635$ nm. These materials are characterized by a helical structure of
pitch $p$ that produces a specific optical response for wavelengths within a
bandwidth centered at $\lambda_{0}=\bar{n}p$, at normal incidence, where
$\bar{n}$ is the mean refractive index. The reflected light is circularly
polarized. When unpolarized light is incident on a CLC, a maximum of 50% of
light is reflected (matching the helix handedness), and 50% of the light is
transmitted (circularly polarized with the inverse handedness) [9]. In the
following, this case will be referred to as "Bragg film", in analogy with
X-ray diffraction. On the contrary, if the wavelength $\lambda$ is outside the
bandwidth, the polarization rule is ineffective and the material will be
referred to as "off-Bragg". Wacker-Chemie GmbH provided us with CLCs
polysiloxane-based oligomers. Details regarding their chemical and physical
properties, as well as relevant references, are given in a prior paper [10].
Their glass transition temperature $T_{\textrm{g}}$ ranges from 40 to 55°C.
The helix is left-handed. Thin films were produced between two plain glass
substrates and annealed at 140°C in their viscous state. They can be vitrified
by quenching below $T_{\textrm{g}}$, and their cholesteric structure is
preserved at normal temperature in a solid state. The Silicon-Green chemical
was used to fabricate the off-Bragg film (half height bandwidth 430–500 nm).
It was annealed for 90 min without the top substrate, resulting in a blue-
shift in the reflection band [11]. The thickness measures $2.6\pm 0.1$ µm. The
Bragg film (half height bandwidth 550–680 nm) is a 13:87 wt. % mixture of
Silicon-Blue and Silicon-Red compounds. It was annealed for 10 min with both
substrates present. The thickness measures $15\pm 2$ µm.

###  3.3 Vectorial ptychography measurements

Measurements were carried out at $\lambda=635$ nm on an optical setup
described previously [2], adapted here by inserting two quarter waveplates
properly oriented, one before the object, one after, as shown in Fig. S1 in
Supplement 1. This modification provides a simple solution to upgrade the
linear polarization scheme to the improved one proposed in this work. Thus,
the mechanical control of the polarization, formerly of the angles of linear
polarizations (Eqs. 3) acts now as a control of the polarization ellipticities
(Eqs. 6). Measurements were carried out with an illumination probe of
effective diameter 100 µm. The far-field was collected through a numerical
aperture of 0.4 and recorded by a camera of effective dimensions $320\times
240$ pixels. Object reconstructions were performed by means of 500 iterations
of a conjugate-gradient algorithm described in a previous work [12], allowing
the joint estimation of the three probes together with the Jones maps of the
object, with a pixel size of about $0.73\times 0.97$ µm2.

##  4 Results and discussion

Figure 1: Description of the simulated numerical object. (a) Jones maps
$\mathbf{J}(\mathbf{r})$. Inset show the colour encoding for complex values in
the complex plane. (b) Pauli coefficients maps. (c) Values of the complex
Pauli coefficients represented in the complex plane. Values have been averaged
in the square dotted area shown on the maps of panel b. Scale bars are 20 µm.

We performed a simulation using an object described by the Jones matrix
presented in Eq. 4. The transmittance properties $T(\mathbf{r})$ of the object
corresponded to a uniform region intersected by vertical and horizontal
stripes with lower transmittance. The resulting Jones maps are depicted in
Fig. 1a. To better identify the presence of circular polarization properties
while preserving all information, we expanded them as a sum

$\mathbf{J}(\mathbf{r})=\sum_{n=0}^{3}C_{n}(\mathbf{r})\
\boldsymbol{\sigma}_{n}.$ (8)

Here, $\boldsymbol{\sigma}_{0}$ denotes the identity matrix, and
$\boldsymbol{\sigma}_{1}$, $\boldsymbol{\sigma}_{2}$, and
$\boldsymbol{\sigma}_{3}$ represent the Pauli matrices [8]

$\boldsymbol{\sigma}_{1}=\begin{bmatrix}1&0\\\ 0&-1\end{bmatrix};\
\boldsymbol{\sigma}_{2}=\begin{bmatrix}0&1\\\ 1&0\end{bmatrix};\
\boldsymbol{\sigma}_{3}=\begin{bmatrix}0&-i\\\ i&0\end{bmatrix}.$ (9)

The corresponding maps of the four complex Pauli coefficients
$C_{n}(\mathbf{r})$ are displayed in Fig. 1b. Their average values, obtained
near the center of the image, are illustrated in the complex plane (Fig. 1c).
These values highlight the homogeneous left-circular polarizer properties of
the simulated object, where $C_{0}$ and $C_{3}$ are equal, positive, and real,
while $C_{1}=C_{2}=0$.

Figure 2: Simulations of measurements with the linear polarization scheme. (a)
Combinations of polarizations for illuminations
($\mathbf{p}^{\mathrm{lin}}_{k}$, solid circle) and analyses
($\mathbf{h}^{\mathrm{lin}}_{l}$, empty circle) shown on the Poincaré sphere.
$S_{1}$, $S_{2}$ and $S_{3}$ are the standard Stokes coefficients [8]. RCP and
LCP stand for right and left circular polarization, respectively. (b)
Illustration of these polarizations. The colour indicates the ellipticity,
given by the latitude on the Poincaré sphere. (c-g) Pauli coefficient maps
calculated from retrieved Jones maps, for several reconstruction, run
independently, represented with the same convention as in Figs. 1b and c.
Scale bars are 20 µm.

First, we conducted a simulation to replicate a measurement using the linear-
polarization scheme. The specific angles employed were $\alpha_{1}=0$°,
$\alpha_{2}=60$°, $\alpha_{3}=120$°, along with $\theta_{1}=0$°,
$\theta_{2}=60$°, $\theta_{3}=120$° [7] in Eqs. 3, as visually depicted in
Fig. 2a,b. Multiple independent algorithm runs were executed, yielding diverse
solutions. Figs 2c-e provide an overview by displaying maps of the Pauli
coefficients for each solution. Alongside the correct solution (Fig. 2c),
which indicated a homogeneous left-circular polarizer, several alternative
solutions emerged. These alternative solutions comprised of either homogeneous
right-circular polarizers (Fig. 2d) or inhomogeneous ones (Fig.2e), thereby
confirming the theoretical ambiguity highlighted. Remarkably, the discovered
solutions extended beyond circular-polarizing characteristics, exhibiting
different properties. For example, Fig. 2f corresponded to a horizontal
quarter waveplate, while the last case, Figs. 2g, represented more intricate
combinations.

Figure 3: Simulations of measurements with the improved polarization scheme.
(a) Combinations of polarizations for illuminations
($\mathbf{p}^{\mathrm{ell}}_{k}$, solid circle) and analyses
($\mathbf{h}^{\mathrm{ell}}_{l}$, empty circle) shown on the Poincaré sphere.
(b) Corresponding polarizations. (c) Pauli coefficient maps calculated from
retrieved Jones maps, as obtained systematically for any reconstruction run
independently, represented with the same convention as in Figs. 1b and c.
Scale bars are 20 µm.

Subsequently, the simulation was executed using an improved polarization
scheme. The chosen angles were $\alpha_{1}=15$°, $\alpha_{2}=105$°,
$\alpha_{3}=45$°, in conjunction with $\theta_{1}=0$°, $\theta_{2}=-45$°,
$\theta_{3}=45$° within Eqs. 6, visually represented in Fig. 3a. These angles,
chosen empirically, represented a combination of linear, elliptical, and
circular polarization states, allowing to probe efficiently a large variety of
optical properties. Once again, multiple runs of the algorithm were performed,
resulting in consistent convergence towards the correct solution, now. The
obtained solution is illustrated in Fig. 3b, further validating the
effectiveness of this scheme.

Finally, Bragg and off-Bragg CLC films were investigated using vectorial
ptychography with the improved polarization scheme. Fig. 4a illustrates the
Jones and Pauli coefficients maps obtained from the reconstruction algorithm.
For the off-Bragg film (Fig. 4a), the Pauli coefficients maps show that only
the $C_{0}(\mathbf{r})$ map has non-zero values. This confirms that the film
has an optical response that is insensitive to polarization. On the contrary,
the Bragg film (Fig. 4b) displays different maps of Pauli coefficients, with
$C_{0}(\mathbf{r})=C_{3}(\mathbf{r})\neq 0$, while
$C_{1}(\mathbf{r})=C_{2}(\mathbf{r})=0$. This confirms that the film functions
as a uniform left-circular polarizer. Furthermore, variations in the Pauli
coefficient allow to report a texturing of the films (polygonal texture) that
has been discussed previously [11]. Addressing this aspect is beyond the scope
of the present study and will be the subject of a future article. It is worth
emphasizing the clear advantage of spatial resolution and intrinsic phase-
imaging capabilities offered by vectorial ptychography [4]. In contrast, other
interferometric methods typically provide spatially averaged measurements
[13].

Figure 4: Investigation of two CLC films. Pauli coefficient maps. (a) Off-
Bragg film. (b) Bragg film. Scale bars are 40 µm.

##  5 Conclusion

In conclusion, this letter presents a novel approach to vectorial
ptychography, expanding its capabilities to circular-polarizing materials such
as CLCs. By a measurement scheme based on elliptical polarizations, the
technique overcomes the limitations of linear polarizations and enables the
accurate characterization of circular polarizers at a micrometer-scale
resolution. This progress in quantitative imaging has implications for various
fields, including the study of chiral molecular assemblies and other materials
or advanced components with strong circular-polarization properties.

Funding We acknowledge funding from the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation program (Grant
agreement no. 724881).

Acknowledgments We thank Dr. E. Hanelt from Wacker-Chemie GmbH (Munich,
Germany) for providing us with cholesteric oligomers.

Disclosures The authors declare no conflicts of interest.

Data Availability Statement The datasets are available upon reasonable request
from the corresponding author.

Supplemental document See Supplement 1 for supporting content.

## References

  * [1] P. Ferrand, M. Allain, and V. Chamard, Optics Letters 40, 5144 (2015).
  * [2] P. Ferrand, A. Baroni, M. Allain, and V. Chamard, Optics Letters 43, 763 (2018).
  * [3] T. Wang, S. Jiang, P. Song, R. Wang, L. Yang, T. Zhang, and G. Zheng, Biomedical Optics Express 14, 489 (2023).
  * [4] A. Baroni, V. Chamard, and P. Ferrand, Physical Review Applied 13, 054028 (2020).
  * [5] J. Duboisset, P. Ferrand, A. Baroni, T. A. Grünewald, H. Dicko, O. Grauby, J. Vidal-Dupiol, D. Saulnier, L. M. Gilles, M. Rosenthal, M. Burghammer, J. Nouet, C. Chevallard, A. Baronnet, and V. Chamard, Acta Biomaterialia 142, 194 (2022).
  * [6] X. Dai, S. Xu, X. Yang, K. C. Zhou, C. Glass, P. C. Konda, and R. Horstmeyer, Biomedical Optics Express 13, 1457 (2022).
  * [7] Q. Song, A. Baroni, R. Sawant, P. Ni, V. Brandli, S. Chenot, S. Vézian, B. Damilano, P. de Mierry, S. Khadir, P. Ferrand, and P. Genevet, Nature Communications 11, 2651 (2020).
  * [8] R. A. Chipman, W. S. T. Lam, G. Young, W. S. T. Lam, and G. Young, _Polarized Light and Optical Systems_ (CRC Press, 2018).
  * [9] M. Mitov, Adv. Mater. 24, 6260 (2012).
  * [10] A. Scarangella, V. Soldan, and M. Mitov, Nature Commun. 11, 4108 (2020).
  * [11] G. Agez, R. Bitar, and M. Mitov, Soft Matter 7, 2841 (2011).
  * [12] A. Baroni, M. Allain, P. Li, V. Chamard, and P. Ferrand, Opt. Exp. 27, 8143 (2019).
  * [13] A. Sanchez-Castillo, S. Eslami, F. Giesselmann, and P. Fischer, Opt. Exp. 22, 31227 (2014).

myrefs