# Ethically aligned Deep Learning: Unbiased Facial Aesthetic Prediction

Michael Danner111 222Reutlingen , Thomas Weber22footnotemark: 2 , Leping
Peng333Hunan , Tobias Gerlach22footnotemark: 2 , Xueping Su444Xian , Matthias
Rätsch22footnotemark: 2 44footnotemark: 4
11footnotemark: 1 University of Surrey 22footnotemark: 2 Reutlingen University
33footnotemark: 3 Hunan University of Science and Technology
44footnotemark: 4 Xi’an Polytechnic University

###### Abstract

Facial beauty prediction (FBP) aims to develop a machine that automatically
makes facial attractiveness assessment. In the past those results were highly
correlated with human ratings, therefore also with their bias in annotating.
As artificial intelligence can have racist and discriminatory tendencies, the
cause of skews in the data must be identified. Development of training data
and AI algorithms that are robust against biased information is a new
challenge for scientists. As aesthetic judgement usually is biased, we want to
take it one step further and propose an Unbiased Convolutional Neural Network
for FBP. While it is possible to create network models that can rate
attractiveness of faces on a high level, from an ethical point of view, it is
equally important to make sure the model is unbiased. In this work, we
introduce AestheticNet, a state-of-the-art attractiveness prediction network,
which significantly outperforms competitors with a Pearson Correlation of
0.9601. Additionally, we propose a new approach for generating a bias-free CNN
to improve fairness in machine learning.

_K_ eywords Fairness in Machine Learning $\cdot$ Responsible Artificial
Intelligence $\cdot$ Discrimination Prevention $\cdot$ Facial Aesthetics
$\cdot$ Unconscious Bias

##  1 Introduction

In 2016 "Beauty.ai", a Hong-Kong based technology company hosted the first
international beauty contest judged by artificial intelligence [1] but the
results are heavily biased for example against dark skin [2]. “Machine
learning models are prone to biased decisions, due to biases in data-sets”[3].
Biased training data potentially leads to discriminatory models, as the data-
sets are created by humans or derived from human activities in the past, for
example hiring algorithms [4]. The reason for racist and discriminatory
tendencies must be identified. As the learning algorithms become more complex,
understanding why the decisions are made, or even how, prove to be nearly
impossible [5]. Therefore the development of non-biased training data and AI
algorithms (defined by the European Commission High-Level Expert Group on
Artificial Intelligence [6]) is a new and increasingly complex challenge for
scientists around the world. The specific field of aesthetic judgement is
especially vulnerable to being biased, as aesthetic judgement itself is
already a subjective rating [7]. Our data analysis has already proven that
people consider their own ethnicity to be more attractive than others [8],
this is the biggest bias in our research and within our data. With this
tendency, it becomes difficult to generate data to train a machine-learning
algorithm, which assesses a person’s attractiveness without bias.

##  2 State of the art

Machine learning has evolved in the past decades and stands out due to the
fact that the knowledge in the system is not provided by experts. Facial
Beauty Prediction (FBP), that is consistent with human perception, is a
significant visual recognition problem and a much-studied subject in recent
decades. Eisenthal et al. [9] and Kagian et al. [10] were among the first to
publish their research about automatic facial attractiveness predictors and
supervised learning techniques, based on the extraction of feature landmarks
on faces.

With the introduction of convolutional neural networks (CNN) and large-scale
image repositories, facial image and video tasks get more powerful [11, 12,
13]. Xie et al. [14] present the SCUT-FBP500 dataset, containing 500 asian
female subjects with attractiveness ratings. Since “FBP is a multi-paradigm
computation problem” the successor SCUT-FBP5500 [15] is introduced in 2018,
including an increased database of 5500 frontal faces with multiple
attributes: male/female, Asian/Caucasian, age, beauty score. Liang et al.
(2018) evaluated their database “using different combinations of feature and
predictor, and various deep learning methods” on AlexNet [11], ResNet-18 [16]
and ResNeXt-50 and achieved the Pearson Correlation PC: $0.8777$; mean average
error MAE: $0.2518$; root-mean-square error RMSE: $0.3325$ as a benchmark. It
was observed, that all deep CNN models are superior to the shallow predictor
with hand-crafted geometric feature or appearance feature [15].

##  3 Ethics

Since a large number of studies have shown that machines trained by humans
will also repeat human bias [17, 18, 19, 20]. Unfortunately, everyone is
biased, because a large part of human bias is unconscious and it is not easy
to be noticed, which was demonstrated in a Seattle press conference at the
University of Washington by psychologists who developed a new tool that
measures the unconscious bias [21]. Unconscious bias comes from educational
background, culture, attitudes, and stereotypes we pick up from the world we
live in.

An example of gender bias: many people have the stereotype that females are
worse at mathematics tasks and better at verbal tasks than males [22]. A Yale
study shows that male candidates are judged to be more talented and
experienced and hired more often [23]. Affinity bias refers to unconsciously
preferring familiar people, this is verified in our first part of our
experiment. A machine trained on a data-set labelled by Europeans judges
European faces higher than Asian faces and vice versa. And when you
unconsciously notice people’s appearances and associate it with their
personality, you might have beauty bias.

How to build a fair artificial intelligence if human beings are biased? Many
technology companies begin to realise this problem and take efforts, such as
Google [24] and Facebook [25]. Facebook makes unconscious bias training videos
widely available and Google educates employees by an unconscious bias training
program to reduce human bias during human-computer interaction.

Bias also exists in the field of image processing, for example, artificial
intelligence applications tag pictures of White American brides as “brides”,
“dresses”, and “weddings” while pictures of North Indian brides are tagged as
“performing arts” and “costumes” [26]. Angwin’s and Larson’s [27] analysis of
ethical bias has prompted research showing that the disparity can be addressed
if the algorithms focus on the fairness of outcomes.

What is the significance of our research? Although both Hume [28] and Kant
[29] believe that aesthetic judgement is only a subjective feeling, which
regarding the pleasure that we take from a beautiful object, and aesthetics
itself is neither right nor wrong nor moral, the behaviour caused by aesthetic
judgements are related to morality [30, 31]. Many studies show that unfair
recruitment cases encountered due to aesthetic judgements, [32, 33, 34], more
attractive people have higher incomes than less attractive individuals [35,
36, 37]. Our research provides a practical example of how to build a fair and
trustworthy AI.

##  4 Experiment

We conducted online surveys on image datasets beginning in 2013 on multiple
datasets where thousands of students and their relatives participated.
Starting 2017 we used the Asian-European-dataset SCUT-FBP to evaluate biased
annotations from Chinese and German universities. The results proved the
assumption that German students favour images of European women and vice versa
Chinese students rate Asian portraits higher. Within our research, we use
Convolutional Neural Networks (CNN) for facial aesthetic score prediction and
introduce AestheticNet as our approach.

Since the SCUT-FBP 5500 dataset is a small dataset for deep learning tasks, we
use data augmentation methods to enlarge the sample size of the training set
by generating GAN images with either Asian or European or mixed images as
input and new synthesised images as output. This augmentation method proves
superior to geometric transformations like cropping and rotating. All images
are preprocessed by normalisation methods to harmonise face pose, facial
landmark positions and image size. For our training data, we use the SCUT-FBP
dataset and labels enhanced with Chinese annotators to train an Asian-biased
network first. The second network is trained with the same dataset but
labelled by German students and their relatives.

We compare our method with other state-of-the-art approaches on the SCUT-
FBP500 datasets. As shown in Table 1 our AestheticNet therefore significantly
surpasses previous approaches, which is mainly due to the augmentation with
synthetic images and the optimisation of the previous approaches. In our best
experiment, we achieve a Pearson correlation of 0.9601, a normalised mean
average error of 3.896% and a normalised root mean squared error of 5.580%.
The results are normalised because there are different datasets with different
score ranges.

Table 1: Comparison of prediction accuracy on SCUT-FBP5500 | PC | nMAE (%) | nRMSE (%)
---|---|---|---
AlexNet [15] | 0.8298 | 7.345 | 9.548
AlexNet [38] | 0.8634 |  |
ResNet-18 [15] | 0.8513 | 7.045 | 9.258
ResNeXt-50 [15] | 0.8777 | 6.295 | 8.313
HMTNet[39] | 0.8783 | 6.2525 |
8.158 AaNet [40] | 0.9055 | 5.590 | 7.385
P-AaNet [40] | 0.8965 | 5.713 | 7.588
2M BeautyNet [41] | 0.8996 |  |
EfficientNetB3 based AestheticNet (ours) | 0.9011 | 5.841 | 7.663
VGG-Face based AestheticNet (ours) | 0.9363 | 4.400 | 6.261
AestheticNet (ours) | 0.9601 | 3.896 | 5.580
Figure 1: Unconscious bias towards ethnic aesthetic of either German or
Chinese annotators. Left: average aesthetic score on SCUT-FBP by German
annotators, middle: average aesthetic score labelled by Chinese students,
right: aesthetic scores on the Eurasian dataset annotated by German students.

Figure 1 shows the biased average score of our networks on the SCUT-FBP
dataset and the Eurasian dataset. Figure 2 illustrates the analysis on the
distribution of aesthetic score and age for Asians, Europeans and three mixed-
racial subgroups. In a fair machine, the distance of all lines would overlap
themselves.

Having a state-of-the-art aesthetic prediction network, we then train a third
CNN on the features from the Asian and German labelled networks to generate a
non-biased network. Therefore, the synthesised Eurasian dataset is used with a
categorical-cross-entropy-loss-function to converge the subgroup’s
intersection points of Figure 2.

Figure 2: Biased correlation between attractiveness, age and ethnicity by
German annotators. In an ethnical fair network the attractiveness would be the
same height for equal age groups.

##  5 Conclusion

Our two main contributions are AestheticNet and a new approach to bias-free
machine learning tools. In this work, we have proposed to augment the SCUT-FBP
dataset by synthesised GAN images and show that AestheticNet predicts facial
attractiveness with higher correlation then competitive approaches. Then we
utilise categorical-cross-entropy-loss learning strategies to minimise bias in
networks. Unbiased networks are an important step towards a future, where more
decisions are made by AI and therefore more lives are influenced by artificial
intelligence - unbiased decision making is the foundation of ethical and moral
values.

Bias-free decision making is a challenging problem in machine learning tasks,
yet it yields the great potential to be one of the most significant strengths
of an AI. We have shown a method to eliminate bias in facial attractiveness
prediction and this method can be transferred to multiple similar networks.

## References

  * [1] beauty.ai, “The First International Beauty Contest Judged by Artificial Intelligence.”
  * [2] S. Levin, “A beauty contest was judged by AI and the robots didn’t like dark skin,” Sept. 2016.  Section: Technology.
  * [3] S. Sharma, Y. Zhang, J. M. Ríos Aliaga, D. Bouneffouf, V. Muthusamy, and K. R. Varshney, “Data Augmentation for Discrimination Prevention and Bias Disambiguation,” in Proceedings of the AAAI/ACM Conference on AI, Ethics, and Society, AIES ’20, (New York, NY, USA), pp. 358–364, Association for Computing Machinery, Feb. 2020.
  * [4] M. Bogen, “All the Ways Hiring Algorithms Can Introduce Bias,” Harvard Business Review, May 2019.  Section: Hiring.
  * [5] N. Bostrom and E. Yudkowsky, “The Ethics of Artificial Intelligence,” Chapman and Hall, p. 21, 2018.
  * [6] European Commission High-Level Expert Group on Artificial Intelligence [AI HLEG], “A definition of AI,” tech. rep., 2019.
  * [7] S. H. Richmond, “The Beholder’s I: The Perception of Beauty and the Development of the Self,” Perception of Beauty, Oct. 2017.  Publisher: IntechOpen.
  * [8] T. Gerlach, M. Danner, L. Peng, A. Kaminickas, W. Fei, and M. Rätsch, “Who Loves Virtue as much as He Loves Beauty?: Deep Learning based Estimator for Aesthetics of Portraits:,” in Proceedings of the 15th International Joint Conference on Computer Vision, Imaging and Computer Graphics Theory and Applications, (Valletta, Malta), pp. 521–528, SCITEPRESS - Science and Technology Publications, 2020.
  * [9] Y. Eisenthal, G. Dror, and E. Ruppin, “Facial Attractiveness: Beauty and the Machine,” Neural Computation, vol. 18, pp. 119–142, Jan. 2006.
  * [10] A. Kagian, G. Dror, T. Leyvand, I. Meilijson, D. Cohen-Or, and E. Ruppin, “A machine learning predictor of facial attractiveness revealing human-like psychophysical biases,” Vision Research, vol. 48, pp. 235–243, Jan. 2008\.
  * [11] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “ImageNet classification with deep convolutional neural networks,” Communications of the ACM, vol. 60, pp. 84–90, May 2017.
  * [12] M. D. Zeiler and R. Fergus, “Visualizing and Understanding Convolutional Networks,” arXiv:1311.2901 [cs], Nov. 2013.  arXiv: 1311.2901.
  * [13] J. Deng, W. Dong, R. Socher, L. Li, Kai Li, and Li Fei-Fei, “ImageNet: A large-scale hierarchical image database,” in 2009 IEEE Conference on Computer Vision and Pattern Recognition, pp. 248–255, June 2009.  ISSN: 1063-6919.
  * [14] D. Xie, L. Liang, L. Jin, J. Xu, and M. Li, “SCUT-FBP: A Benchmark Dataset for Facial Beauty Perception,” arXiv:1511.02459 [cs], Nov. 2015.  arXiv: 1511.02459.
  * [15] L. Liang, L. Lin, L. Jin, D. Xie, and M. Li, “SCUT-FBP5500: A Diverse Benchmark Dataset for Multi-Paradigm Facial Beauty Prediction,” arXiv:1801.06345 [cs], Jan. 2018.  arXiv: 1801.06345.
  * [16] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” CoRR, vol. abs/1512.03385, 2015.
  * [17] F. Zuiderveen Borgesius, “Discrimination, artificial intelligence, and algorithmic decision-making,” Directorate General of Democracy © Council of Europe, p. 51, 2018.
  * [18] S. Stephens-Davidowitz, “The cost of racial animus on a black candidate: Evidence using Google search data,” Journal of Public Economics, vol. 118, pp. 26–40, Oct. 2014.
  * [19] K. Munger, “Tweetment Effects on the Tweeted: Experimentally Reducing Racist Harassment,” Political Behavior, vol. 39, pp. 629–649, Sept. 2017.
  * [20] E. Horvitz, “AI, people, and society,” Science, vol. 357, pp. 7–7, July 2017.
  * [21] Joel Schwarz, “Roots of unconscious prejudice affect 90 to 95 percent of people, psychologists demonstrate at press conference.”
  * [22] M. Johns, T. Schmader, and A. Martens, “Knowing Is Half the Battle: Teaching Stereotype Threat as a Means of Improving Women’s Math Performance,” Psychological Science, vol. 16, pp. 175–179, Mar. 2005\.
  * [23] C. A. Moss-Racusin, J. F. Dovidio, V. L. Brescoll, M. J. Graham, and J. Handelsman, “Science faculty’s subtle gender biases favor male students,” Proceedings of the National Academy of Sciences, vol. 109, pp. 16474–16479, Oct. 2012.
  * [24] Google LLC, “re:Work - Guide: Raise awareness about unconscious bias.”
  * [25] Facebook Inc., “Managing Bias.”
  * [26] S. Shankar, Y. Halpern, E. Breck, J. Atwood, J. Wilson, and D. Sculley, “No Classification without Representation: Assessing Geodiversity Issues in Open Data Sets for the Developing World,” arXiv:1711.08536 [stat], Nov. 2017.  arXiv: 1711.08536.
  * [27] J. A. Larson, Jeff, “Bias in Criminal Risk Scores Is Mathematically Inevitable, Researchers Say,” 2016.
  * [28] H. David, Of the Standard of Taste.  Essays Moral, Political and Literary, London: Longmans, Green and Co., 1898.  edited by T.H. Green and T. H. Grose.
  * [29] I. Kant, Critique of the Power of Judgment.  Cambridge University Press, 1 ed., Sept. 2000.
  * [30] X. Cui, Q. Cheng, W. Lin, J. Lin, and L. Mo, “Different influences of facial attractiveness on judgments of moral beauty and moral goodness,” Scientific Reports, vol. 9, Dec. 2019.
  * [31] J. Haidt, “The emotional dog and its rational tail: A social intuitionist approach to moral judgment.,” Psychological Review, vol. 108, no. 4, pp. 814–834, 2001.
  * [32] K. B. Maddox and J. M. Perry, “Racial Appearance Bias: Improving Evidence-Based Policies to Address Racial Disparities,” Policy Insights from the Behavioral and Brain Sciences, vol. 5, pp. 57–65, Mar. 2018.
  * [33] A. Mason, “Appearance, Discrimination, and Reaction Qualifications,” Journal of Political Philosophy, vol. 25, no. 1, pp. 48–71, 2017.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/jopp.12099.
  * [34] G. Beattie and P. Johnson, “Possible unconscious bias in recruitment and promotion and the need to promote equality,” Perspectives: Policy and Practice in Higher Education, vol. 16, pp. 7–13, Jan. 2012.
  * [35] P. Anýžová and P. Matějů, “Beauty still matters: The role of attractiveness in labour market outcomes,” International Sociology, vol. 33, pp. 269–291, May 2018.
  * [36] M. T. French, “Physical appearance and earnings: further evidence,” Applied Economics, vol. 34, pp. 569–572, Mar. 2002.
  * [37] J. Cawley, “The Impact of Obesity on Wages,” Journal of Human Resources, vol. XXXIX, no. 2, pp. 451–474, 2004.
  * [38] Y. Zhai, H. Cao, W. Deng, J. Gan, V. Piuri, and J. Zeng, “Beautynet: Joint multiscale CNN and transfer learning method for unconstrained facial beauty prediction,” Comput. Intell. Neurosci., vol. 2019, pp. 1910624:1–1910624:14, 2019.
  * [39] L. Xu, H. Fan, and J. Xiang, “Hierarchical Multi-Task Network For Race, Gender and Facial Attractiveness Recognition,” in 2019 IEEE International Conference on Image Processing (ICIP), pp. 3861–3865, Sept. 2019.  ISSN: 2381-8549.
  * [40] L. Lin, L. Liang, L. Jin, and W. Chen, “Attribute-aware convolutional neural networks for facial beauty prediction,” in Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence, IJCAI 2019, Macao, China, August 10-16, 2019 (S. Kraus, ed.), pp. 847–853, ijcai.org, 2019.
  * [41] J. Gan, F. Scotti, L. Xiang, Y. Zhai, C. Mai, G. He, J. Zeng, Z. Bai, R. D. Labati, and V. Piuri, “2m beautynet: Facial beauty prediction based on multi-task transfer learning,” IEEE Access, vol. 8, pp. 20245–20256, 2020.


# Decomposing Representations for
Deterministic Uncertainty Estimation

Haiwen Huang Joost van Amersfoort Yarin Gal
OATML Group
Department of Computer Science
University of Oxford
Oxford, United Kingdom
{haiwen.huang2, joost.van.amersfoort<EMAIL_ADDRESS>

##  1 Introduction

Uncertainty estimation is a key component in any deployed machine learning
system: when the system encounters unfamiliar or noisy input, it needs to
refrain from making a decision and instead defer to an expert. There has been
great progress in uncertainty estimation in recent years, particularly in a
sub-field called deterministic uncertainty estimation, where only a single
forward pass is necessary to accurately estimate uncertainty [1, 2, 3].

One way to evaluate uncertainty estimation is using “out-of-distribution”
(OoD) detection, where it is expected that data points far away from the
training distribution have higher uncertainty than points similar to the
training distribution. A concrete way to evaluate this is by distinguishing
between two data distributions, one used for training and one never seen
during training, using uncertainty. Recently a challenging version of this
benchmark was introduced: distinguishing CIFAR-10 from CIFAR-100 and vice
versa. This was called “near” OoD detection to distinguish from “far” OoD
detection [4], for example in CIFAR-10 vs SVHN. Near OoD detection is
significantly more difficult as the near OoD data are more semantically
similar to the training data. For example, CIFAR-10 and CIFAR-100 both come
from the 80-million tiny images dataset, but contain no class overlap [5].

It was shown that in the near OoD setting, standard methods such as
Mahalanobis distance [1] and DDU [2] do not perform well [6]. We show that the
proposed fix in [6] does not work well in the far OoD setting, and propose a
simple solution that does work. We identify two types of features:
discriminative (class-relevant) and non-discriminative (class-irrelevant)
features. Intuitively, one can consider the discriminative features sufficient
for distinguishing near OoD data, while non-discriminative features are used
for far OoD detection. In Table 1, we show that when using a representation
that is a mix of these features, uncertainty estimators based on feature
density that perform well on one OoD detection benchmark have poor performance
on the other.

This is because when the OoD data mainly differ in one specific type of
representations, the other type can still influence the uncertainty. To solve
this, we propose to decompose the learned representations into two proposed
types of representations and integrate the uncertainties estimated on them
separately. Through experiments, we demonstrate that we can greatly improve
the performance and the interpretability of the uncertainty estimation.

##  2 Motivating Observations

First, we use Mahalanobis distance for OoD detection [1] and its variants [2,
4, 6] to demonstrate the influence of class-relevance on uncertainty
estimation. Since our analysis is focused on the feature maps, it can be
extended to any other representation-based uncertainty estimation methods.

Given a pretrained feature extractor
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ and dataset
$\\{x_{i},y_{i}\\}_{i=1}^{N}$ where input $x_{i}\in\mathbb{R}^{n}$, and label
$y_{i}\in\\{1,\cdots,C\\}$, the Mahalanobis distance for OoD detection first
models the distribution of the features as a class-conditional multivariate
Gaussian distribution with a single covariance matrix shared across classes.
That is, for each class $c$, it models the distribution of features as
$p(z|y=c)\sim\mathcal{N}(\boldsymbol{\hat{\mu}}_{c},\Sigma)$. It then
calculates class-wise Mahalonobis distance (i.e. the same quantity as in the
exponent of a multi-variate Gaussian distribution) of a feature representation
$z=f(x)$ and uses the maximum across classes as the distance:

$M(x)=\max_{c}-(f(x)-\boldsymbol{\hat{\mu}}_{c})^{T}\hat{\Sigma}^{-1}(f(x)-\boldsymbol{\hat{\mu}}_{c}),$
(1)

where
$\hat{\Sigma}=\frac{1}{N}\sum_{c=1}^{C}\sum_{i:y_{i}=c}\left(f(\boldsymbol{x}_{i})-\boldsymbol{\hat{\mu}}_{c}\right)\left(f(\boldsymbol{x}_{i})-\boldsymbol{\hat{\mu}}_{c}\right)^{T},\;\boldsymbol{\hat{\mu}}_{c}=\sum_{i:y_{i}=c}f(\boldsymbol{x}_{i}).$
Recently, many works [2, 4] model the feature distribution using a class-
dependent covariance matrix and found it beneficial. In this case, the
uncertainty metric is
$M(x)=\max_{c}[-\left(f(x)-\boldsymbol{\hat{\mu}}_{c}\right)^{T}\boldsymbol{\hat{\Sigma}}_{c}^{-1}\left(f(x)-\boldsymbol{\hat{\mu}}_{c}\right)-\log((2\pi)^{m}\operatorname{det}\mathbf{\hat{\Sigma}}_{c})]$
where
$\mathbf{\hat{\Sigma}}_{c}=\frac{1}{|\\{i:y_{i}=c\\}|}\sum_{i:y_{i}=c}(f(\boldsymbol{x}_{i})-\boldsymbol{\hat{\mu}}_{c})(f(\boldsymbol{x}_{i})-\boldsymbol{\hat{\mu}}_{c})^{T}.$

Variants of Mahalanobis distance were proposed to improve on some specific
tasks. Notably, marginal Mahalanobis distance (Marginal Maha) [7] calculates
the mean and covariance without using class-relevant information.
Specifically, in Equation (1),
$\Sigma=\frac{1}{N}\sum_{i=1}^{N}\left(f(\boldsymbol{x}_{i})-\boldsymbol{\mu}\right)\left(f(\boldsymbol{x}_{i})-\boldsymbol{\mu}\right)^{T},\;\boldsymbol{\mu}_{c}=\boldsymbol{\mu}=\frac{1}{N}\sum_{i=1}^{N}f(\boldsymbol{x}_{i}).$
Another variant is called Relative Mahalanobis distance (Relative Maha) [6]
$M_{rel}(x)=M(x)-M_{marg}(x)$ where $M_{rel}$ and $M_{marg}$ refers to
Relative Maha and Marginal Maha, respectively. It aims to mitigate the
influence of the class-irrelevant information. As we can see from the
definitions, the key difference among the three variants is how the class-
relevant information is (not) used. In terms of the degree of class-relevance:
Relative Maha $>$ Maha $>$ Marginal Maha.

To empirically compare these three methods, we use two OoD detection
benchmarks: CIFAR10 [5] vs. SVHN [8] and CIFAR10 vs. CIFAR100. The OoD
detection results can be found in Table 1. We can see that on CIFAR10 vs.
SVHN, all the methods can perform relatively well with Maha being the top.
While on CIFAR10 vs. CIFAR100, there is a sharp difference among the methods:
Relative Maha shows significantly better performance than others while
Marginal Maha almost completely mixes the two datasets (AUROC=50% means random
guess). Note that although the compared methods are well-known, none of them
has been tested simultaneously on both benchmarks. Our experiments demonstrate
that none of these methods can beat every other method on both benchmarks.

Table 1: AUROC (%) comparison of variants of Mahalanobis distance calculated
on features at penultimate layer (pre-logit).

Method | CIFAR10 vs SVHN | CIFAR10 vs CIFAR100
---|---|---
Maha [1] | 97.96 | 75.75
Marginal Maha [9] | 96.07 | 60.10
Relative Maha [6] | 95.43 | 90.97

##  3 Decomposing Representations

Considering the difference among the three compared methods in the last
section, we hypothesize that the differences in OoD detection performances can
be attributed to the different degrees of class relevance. Therefore, we
propose decomposing the discriminative and non-discriminative information from
the learned representations. This allows us to estimate uncertainty on each
type of representations separately and provides a new perspective to interpret
the uncertainty estimation.

Formally, an ideal decomposition of the features $z$ into the discriminative
part $z_{d}$ and non-discriminative part $z_{n}$, i.e., $z=(z_{d};z_{n})$,
should have the following properties:

  * •

$I(z_{d};z_{n}|y)=0$ (the two decomposed parts should be conditionally
independent given the labels $y$);

  * •

$I(z_{n};y)=0$ ($z_{n}$ is non-discriminative, i.e., independent from the
labels $y$);

  * •

$z_{d}=\operatorname*{arg\,min}_{S}I(z_{d};z)$ where $S=\\{z_{d}\mid
I(z_{d};y)=I(z;y)\\}$ ($z_{d}$ includes and only includes all the
discriminative information in $z$).

To perform the decomposition, a simple way is to use Principal Component
Analysis (PCA) to transform the features and then take the first $d$ principal
components (PCs) as $z_{d}$ and the remaining PCs as $z_{n}$. However, since
PCA is a linear method and does not use any label information, there is no
guarantee that PCA can achieve the ideal decomposition. Alternatively, we can
use the independent cross-entropy (iCE) loss proposed by [10]. Specifically,
given a fixed pretrained feature extractor
$f:\mathbb{R}^{M}\rightarrow\mathbb{R}^{D}$ which maps from inputs to
representations, we first perform a transformation of the features $z=f(x)$
using an invertible function $F:\mathbb{R}^{D}\rightarrow\mathbb{R}^{D}$. The
invertibility is to make sure we keep all the information of the features
after the transformation. Our aim then is to extract discriminative features
$z_{d}$ into the first $d$ dimensions of $F(z)$ and non-discriminative
features $z_{n}$ into the other dimensions. As a convenient choice, we set
$d=C$, i.e., the number of classes of the training dataset. $z_{d}$ then
serves as logits for the classification. We then uses another linear map
$D:\mathbb{R}^{D-d}\rightarrow\mathbb{R}^{C}$ to map $z_{n}$ to a
$C-$dimensional logits. In this way, we can now use the iCE loss to encourage
the decomposition of discriminative and non-discriminative features. The iCE
loss is defined as:

$\min_{\theta}\max_{\phi}\mathcal{L}_{iCE}\left(\theta,\phi\right)=\sum_{i=1}^{C}-y_{i}\log[\text{softmax}[F_{\theta}(z)]]_{i}+\sum_{i=1}^{D-C}y_{i}\log[\text{softmax}[D_{\phi}\left(F_{\theta}(z)_{[C+1:D]}\right)]]_{i}.$
(2)

We point out that for the second loss term, the minimization is concerning a
lower bound on $I(y;z_{n})$, while the maximization aims to tighten the bound
[10]. Therefore, the iCE loss aims to maximize $I(y;z_{s})$ (first term) and
minimize $I(y;z_{n})$ (second term) at the same time to achieve the second and
third properties in our definition of decomposition (note the minimization of
$I(z;z_{d})$ is also achieved during optimization [11]). The first
independence property is not explicitly enforced in our implementation, but
can still be achieved through the optimization. We discuss in the appendix B.

When the decomposition is finished, we can get a new uncertainty estimate
based on the uncertainty calculated on the decomposed features. Because
ideally, $\log p(z)=\log p(z_{d};z_{n})=\log p(z_{d})+\log p(z_{n})$ (since
$I(z_{d};z_{n})=0$), the new uncertainty estimate can simply be
$M^{\prime}(z):=1/2M(z_{d})+1/2M(z_{n})$, where $M$ is an uncertainty
estimator on the features, e.g., Maha. We can also design a KL-based dataset
distance metric based on the decomposition:
$d(D_{in},D_{out})=\mathcal{D}_{KL}[p_{in}(M(z))\|p_{out}(M(z))]$, where
$p_{in}$ and $p_{out}$ are the distributions of the OoD scores (e.g. Maha
distance) on $\mathcal{D}_{in}$ and $\mathcal{D}_{out}$ respectively.
Depending on the features we use (i.e., $z_{d}$ or $z_{n}$), we can define
$d_{dis}$ and $d_{nondis}$ correspondingly. Details can be seen in the
appendix B.3.

##  4 Experiments

Table 2: AUROC (%) comparison of different methods calculated on features at
penultimate layer (pre-logit). $D_{in}$ = CIFAR10.

Method | Datasets
---|---
$D_{out}$ = SVHN |  $D_{out}$ =CIFAR100 |  $D_{out}$ = SVHN $\cup$ CIFAR100
$d_{dis}=3.75$ | $d_{dis}=4.55$ | $d_{dis}=2.08$
$d_{nondis}=4.09$ | $d_{nondis}=1.15$ | $d_{nondis}=1.73$
Maha | 97.96 | 75.75 | 86.86
Marginal Maha | 96.07 | 60.10 | 78.08
Relative Maha | 95.43 | 90.97 | 93.20
PCA top score | 94.59 | 90.20 | 92.39
PCA bottom score | 98.23 | 83.77 | 91.00
PCA top score + bottom score | 97.95 | 89.90 | 93.92
Discriminative score | 96.14 | 90.46 | 93.30
Non-discriminative score | 98.08 | 76.14 | 87.11
Dis score + Non-dis score | 98.41 | 89.95 | 94.18

#### Main Result

In Table 2, we show the comparison of different methods on the three OoD
benchmarks. The scoring function we use is Maha. From Table 2, we can see
that:

  1. 1.

Our dataset distance metric is a good indicator of the OoD types and the AUROC
differences. SVHN is far from CIFAR10 in both features but $d_{nondis}$ is
slightly larger. So while either type of features can achieve high AUROC for
detecting SVHN, non-discriminative features are slightly better. For CIFAR100,
$d_{dis}$ is much larger than $d_{nondis}$. So when using discriminative
features to detect CIFAR100, we can see considerable improvement.

  2. 2.

Discriminative score performs best on discriminative features, and non-
discriminative score performs best on non-discriminative features. Moreover,
the sum of the two scores combines the best of both features, yielding close
to best performances on both benchmarks and best performance on the union
benchmark.

  3. 3.

PCA is a strong baseline for the decomposition. Note that this conclusion can
change when we are not using the feature extractor of a discriminatively
trained neural network.

(a)

(b)

Figure 1: AUROCs of OoD detection using discriminative and non-discriminative
features (y-axis) against the averaged Maha distances on the two features for
each class (x-axis). Every two dots with the same x-axis coordinate belong to
the same class in either SVHN or CIFAR100.

###  4.1 Interpreting Uncertainty

The dataset distance under different features also provides us with a tool to
interpret the uncertainty estimates. In Figure 1, we conduct a class-wise
analysis for a fine-grained understanding of the OoD classes. Specifically, we
calculate the AUROC using discriminative and non-discriminative scores to
detect the 110 classes of SVHN $\cup$ CIFAR100 from CIFAR10 (one class at a
time). This gives us two AUROC scores for each class, dis AUROC and nondis
AUROC. We then plot them against the dataset distance between CIFAR10 and each
class. In this way, we can interpret the OoD detection performances (y-axis)
by looking at dataset distances (x-axis).

We highlight three types of OoD: (1) Classes with much higher non-
discriminative distances than discriminative ones, e.g., SVHN classes. For
these classes, non-discriminative features are more suitable for OoD
detection. (2) Classes with small discriminative distances, e.g., streetcar.
The small discriminative distances are usually due to similar categories in
CIFAR10, e.g., automobile and truck for the streetcar. These classes typically
also have small non-discriminative distances. Usually, discriminative
distances are still larger, so more suitable for OoD detection. (3) Classes
with large discriminative distances, e.g., shark. These classes typically also
have large non-discriminative distances. The OoD performances using the two
kinds of distances are usually similar.

##  5 Discussion

In this work, we showed that current state-of-the-art uncertainty estimation
methods could not consistently outperform other methods across different
benchmarks. We solved this problem by decomposing the features. Specifically,
we summed the OoD scores calculated separately on the discriminative and non-
discriminative features and achieved consistently high performance across
different types of benchmarks. Our decomposition can also help us interpret
the uncertainty by looking at the uncertainty estimates on the decomposed
features.

## References

  * [1] Kimin Lee, Kibok Lee, Honglak Lee, and Jinwoo Shin.  A simple unified framework for detecting out-of-distribution samples and adversarial attacks.  In Proceedings of the 32nd International Conference on Neural Information Processing Systems, NIPS’18, page 7167–7177, Red Hook, NY, USA, 2018. Curran Associates Inc.
  * [2] Jishnu Mukhoti, Andreas Kirsch, Joost van Amersfoort, Philip H. S. Torr, and Yarin Gal.  Deterministic neural networks with appropriate inductive biases capture epistemic and aleatoric uncertainty.  CoRR, abs/2102.11582, 2021.
  * [3] Joost van Amersfoort, Lewis Smith, Andrew Jesson, Oscar Key, and Yarin Gal.  Improving deterministic uncertainty estimation in deep learning for classification and regression, 2021.
  * [4] Jim Winkens, Rudy Bunel, Abhijit Guha Roy, Robert Stanforth, Vivek Natarajan, Joseph R. Ledsam, Patricia MacWilliams, Pushmeet Kohli, Alan Karthikesalingam, Simon Kohl, A. Taylan Cemgil, S. M. Ali Eslami, and Olaf Ronneberger.  Contrastive training for improved out-of-distribution detection.  CoRR, abs/2007.05566, 2020.
  * [5] Alex Krizhevsky.  Learning multiple layers of features from tiny images.  Technical report, 2009.
  * [6] Jie Ren, Stanislav Fort, Jeremiah Liu, Abhijit Guha Roy, Shreyas Padhy, and Balaji Lakshminarayanan.  A Simple Fix to Mahalanobis Distance for Improving Near-OOD Detection.  arXiv:2106.09022 [cs], June 2021.  arXiv: 2106.09022.
  * [7] Ryo Kamoi and Kei Kobayashi.  Why is the Mahalanobis Distance Effective for Anomaly Detection?  arXiv:2003.00402 [cs, stat], April 2020.  arXiv: 2003.00402.
  * [8] Yuval Netzer, Tao Wang, Adam Coates, Alessandro Bissacco, Bo Wu, and Andrew Y. Ng.  Reading digits in natural images with unsupervised feature learning.  In NIPS Workshop on Deep Learning and Unsupervised Feature Learning 2011, 2011.
  * [9] Ryo Kamoi and Kei Kobayashi.  Why is the Mahalanobis Distance Effective for Anomaly Detection?  arXiv:2003.00402 [cs, stat], April 2020.  arXiv: 2003.00402.
  * [10] Joern-Henrik Jacobsen, Jens Behrmann, Richard Zemel, and Matthias Bethge.  Excessive invariance causes adversarial vulnerability.  In International Conference on Learning Representations, 2019.
  * [11] Naftali Tishby, Fernando C Pereira, and William Bialek.  The information bottleneck method.  arXiv preprint physics/0004057, 2000.
  * [12] Alex Krizhevsky, Geoffrey Hinton, et al.  Learning multiple layers of features from tiny images.  Tech Report, 2009.
  * [13] Sergey Zagoruyko and Nikos Komodakis.  Wide residual networks.  In BMVC, 2016.
  * [14] Jens Behrmann, Will Grathwohl, Ricky T. Q. Chen, David Duvenaud, and Joern-Henrik Jacobsen.  Invertible Residual Networks.  In International Conference on Machine Learning, pages 573–582. PMLR, May 2019.  ISSN: 2640-3498.
  * [15] Arthur Gretton, Kenji Fukumizu, Choon Teo, Le Song, Bernhard Schölkopf, and Alex Smola.  A kernel statistical test of independence.  In J. Platt, D. Koller, Y. Singer, and S. Roweis, editors, Advances in Neural Information Processing Systems, volume 20. Curran Associates, Inc., 2008.
  * [16] Qing Wang, Sanjeev R Kulkarni, and Sergio Verdú.  Divergence estimation for multidimensional densities via $k$-nearest-neighbor distances.  IEEE Transactions on Information Theory, 55(5):2392–2405, 2009\.

##  Appendix A Additional Experiments

###  A.1 Toy experiments

We consider both low-dimensional (2D) simulation and high-dimensional (128D)
simulation. For the low dimensional simulation, we generate two Gaussian
distributions with mean (3,0) and (-3,0) and tied covariance matrix
$\Sigma=0.5I_{2}$. The two OoD distributions are also Gaussian distributions
with mean (3,1.4) and (1.6, 0) and the same covariance matrix
$\Sigma=0.3I_{2}$. We then use Equation (1) to calculate the Mahalanobis
distance on the chosen features (full, discriminative, non-discriminative).
The results in Figure 2 show that Maha calculated on the full features can
only reach sub-optimal performance (94% AUROC) compared to on the decomposed
features for the corresponding type of OoD data (98% AUROC).

For the high-dimensional simulation, we generate ten 128D Gaussian
distributions
$\\{\boldsymbol{z}_{c}\sim\mathcal{N}(\boldsymbol{\mu}_{c},\Sigma)\\}_{c=1}^{10}$
with mean $\boldsymbol{\mu}_{c}=10\boldsymbol{e}_{c}$, diagonal covariance
matrix $\Sigma=I_{128}$. Here $\boldsymbol{e}_{c}$ is the standard basis. We
also consider two types of OoD: (1) OoD along discriminative directions:
$\boldsymbol{z}\sim\mathcal{N}(\boldsymbol{\mu}_{k},\Sigma)$ where
$\boldsymbol{\mu}_{k}=5\boldsymbol{e}_{k}$, $k\in\\{1,2,...,10\\}$.; (2) OoD
along non-discriminative directions:
$\boldsymbol{z}\sim\mathcal{N}(\boldsymbol{\mu}_{k^{\prime}},\Sigma)$ where
$\boldsymbol{\mu}_{k^{\prime}}=10\boldsymbol{e}_{c}+5\boldsymbol{e}_{l}$,
$c\in\\{1,2,...,10\\},11\leq l\leq 128$. The remaining process of calculating
Mahalanobis distance is the same as the low-dimensional simulation. When using
Mahalanobis distance on the full feature, the OoD detection performance is
$83.95\%$ AUROC. While we we use Mahalanobis distance on the decomposed
features, we can achieve $99.13\%$ AUROC on OoD data along discriminative
directions and $84.41\%$ along non-discriminative directions. This shows that
the uncertainty estimation is dominated by the non-discriminative
representations. This is similar to the CIFAR10 vs CIFAR100 benchmark (see
Table 2).

(a) Mahalanobis distance on full features.
AUROCred = 0.94,
AUROCgreen=0.94.

(b) Mahalanodis distance on discriminative features.
AUROCred= 0.98,
AUROC${}_{green}=0.36.$

(c) Mahalanodis distance on non-discriminative features.
AUROCgreen = 0.98,
AUROC${}_{red}=0.34.$

Figure 2: A 2D example showing the effect of decomposition. The orange and
blue points are two in-distribution classes, green points are OoD from non-
discriminative directions and red points are OoD from discriminative
directions. Contour lines are colored according to the Mahalanobis distance to
the in-distribution points (darker means higher). We can see that Maha on the
full feature weighs on both discriminative and non-discriminative directions
simultaneously (and in our case, equally). This makes Maha on full features
((a)) suboptimal in detecting OoD data compared to Maha calculated on either
direction only ((b), (c)).

###  A.2 Histogram comparison on image datasets

To further understand the dataset distances under two different features, in
Figure 3, we show the histograms of the two scores (log likelihoods)
calculated on the three datasets. From the histogram, we can see that SVHN is
far from CIFAR10 in both histograms, but the absolute values of non-
discriminative log likelihoods are larger, so they can be easier to detect.
For CIFAR100, there is a sharp difference between two features: the
discriminative log likelihoods are more flat and far from CIFAR100, while the
non-discriminative log likelihoods are very close to CIFAR10. So
discriminative scores are a better choice to detect CIFAR100. These
observations are aligned with the AUROC results of OoD detection and the
dataset distances in Table 2.

###  A.3 Setup of experiments on image datasets

For the experiments on image datasets, we use CIFAR10 [12] as the in-
distribution datasets and SVHN [8], CIFAR100 [12], and the union of the two
datasets as three different benchmarks. This is to demonstrate better how OoD
datasets’ differences influence the OoD detection performances using different
methods.

We use WRN-28-10 [13] trained on CIFAR10 as our pretrained model to extract
the penultimate features for different methods. For the decomposition using
iCE loss, the invertible transformation is implemented as a 4-layer invertible
residual networks [14] with linear residual function and 0.9 spectral
coefficient. The decomposition takes 3000 iterations using SGD optimizer, and
we get 96.12% test accuracy for discriminative logits and 11.2% test accuracy
for non-discriminative logits. For each method, we repeat the experiments
using five differently randomly initialized models and report the mean.

(a)

(b)

Figure 3: OoD detection using uncertainty estimated on different features.
CIFAR10 is the in-distribution dataset, and SVHN and CIFAR100 are two
different OoD datasets. Notably, when using non-discriminative features,
CIFAR100 is much closer to the in-distribution than SVHN.

##  Appendix B Discussion of the Decomposition

###  B.1 Definition of the decomposition

In Section 3, we listed three properties that an ideal decomposition should
possess. For demonstration purposes, we can write them in a slightly redundant
way:

  * •

$I(z_{d};z_{n}|y)=0$ (the two decomposed parts should be conditionally
independent given the labels $y$);

  * •

$I(z_{n};y)=0$ ($z_{n}$ are non-discriminative, i.e., independent from labels
$y$);

  * •

$I(z_{d};y)=I(z;y)$ ($z_{d}$ includes all the discriminative information in
$z$);

  * •

$z_{d}=\operatorname*{arg\,min}_{S}I(z_{d};z)$ where $S=\\{z_{d}\mid
I(z_{d};y)=I(z;y)\\}$ ($z_{d}$ only includes all the discriminative
information in $z$).

We first show that we can derive the third property from the first two. To
start with, from the first two properties we have
$I(z_{d};z_{n})=I(z_{n};y)+I(z_{n};z_{d}\mid y)-I(z_{n};y\mid
z_{d})=-I(z_{n};y\mid z_{d})\leq 0$, thus $I(z_{d};z_{n})=0$. Then since
$I(z_{d};z_{n};y)=I(z_{n};y)-I(z_{n};y\mid
z_{d})=I(z_{d};z_{n})-I(z_{d};z_{n}\mid y)$, we can derive that $I(z_{n};y\mid
z_{d})=I(z_{d};z_{n};y)=0$. Finally, combined with the chain rule of mutual
information, $I(z;y)=I(z_{d};y)+I(z_{n};y\mid z_{d})$, we can derive the third
property $I(z;y)=I(z_{d};y)$. Similarly, when the first and third properties
are satisfied, we can also derive the second property.

The fourth property is to ensure that $z_{d}$ does not include more
information than discriminative part. For example, simply setting $z_{d}=z$
and $z_{n}=0$ could satisfy the first three properties. Note that since
$I(z_{d};y)\leq I(z;y)$, to achieve the fourth property, we may use the
following optimization scheme:

$z_{d}=\operatorname*{arg\,min}I(z_{d};z)-I(z_{d};y).$

We can achieve this by adding a cross entropy loss on the logits ($z_{d}$)
based on the information bottleneck method [11].

We can also show that when the last three properties are satisfied, the first
property will be satisfied as well. As a sketch proof, assume the last three
properties are already satisfied for $z_{d},z_{n}$, if $I(z_{d};z_{n}|y)>0$,
we can then simply set $z_{d}^{\prime}=z_{d}|z_{n}$ 111Note
$I(z_{d};z_{n}|y)=I(z_{d};z_{n})$ when the second and third properties hold.
This is because $I(z_{n};y|z_{d})=I(z;y)-I(z_{d};y)=0$, thus
$I(z_{d};z_{n};y)=I(z_{n};y)-I(z_{n};y|z_{d})=0$. And since we also have
$I(z_{d};z_{n};y)=I(z_{d};z_{n})-I(z_{d};z_{n}|y)$, we now get
$I(z_{d};z_{n}|y)=I(z_{d};z_{n})$. , then
$I(z_{d}^{\prime};z)=I(z_{d}^{\prime};z_{n})+I(z_{d}^{\prime};y)+H(z_{d}^{\prime}|y,z_{n})=I(z_{d};y)+H(z_{d}|y,z_{n})<I(z_{d};z_{n})+I(z_{d};y)+H(z_{d}|y,z_{n})=I(z_{d};z)$,
while $I(z_{n};y)=0$ and $I(z_{d}^{\prime};y)=I(z;y)$ still hold. This leads
to a contradiction against the last property. Intuitively, when $z_{d}$ only
includes discriminative information and $z_{n}$ only includes the rest
information of $z$, $z_{d}$ and $z_{n}$ are conditionally independent given
labels $y$.

###  B.2 Implementation of the decomposition

The relationship among the properties of the decomposition gives us two
options of implementing the decomposition. One implementation is to maximize
$I(z_{d};y)$ (and minimize $I(z_{d};z)$) and minimize $I(z_{n};y)$
simultaneously using iCE loss, as detailed in Section 3 of the paper. From our
discussion above, this will also lead to independence of $z_{d}$ and $z_{n}$.
Alternatively, we can also put cross entropy loss only on $z_{d}$ and add an
independence regularizer (e.g., HSIC [15]) on $z_{d}$ and $z_{n}$. This can
enforce the first and third properties in Section 3 (or first, third, and
fourth in A.1), and as we demonstrated, the second property will be then
satisfied automatically. This can be a future direction for exploration.

It is worth pointing out that throughout our work, we are focusing on the
scenario where we are given a pre-trained feature extractor $f$, so $z=f(x)$
is fixed for each $x$. In practice, the properties of $z$ also play a big role
in the success of the decomposition. If the feature extractor $f$ cannot
extract semantically meaningful features, the decomposition will also fail.
When the feature extractor is bi-Lipschitz, we will have $I(z;y)=I(x;y)$. In
most cases, we can simply think $I(x;y)$ as an upper bound for $I(z_{d};y)$ in
the decomposition, i.e., we want $z_{d}$ to extract all the discriminative
information in the input $x$.

###  B.3 Dataset distance metric based on decomposition

We can design a new dataset distance metric by measuring the difference in the
distributions of uncertainty estimates. The decomposition gives such
divergence a clear meaning: dataset distances in either discriminative or non-
discriminative directions.

To make the scores of the two features comparable, we first normalize the
discriminative (dis) and non-discriminative (non-dis) scores of the test data
using the mean and standard deviation of the dis and non-dis scores calculated
on the training data. To reflect the difference between the two distributions
of scores $p_{in}(M(z))$ and $p_{out}(M(z))$, we can compute a
Kullback–Leibler (KL) divergence
$D_{KL}[p_{in}\|p_{out}]=\sum_{k=1}^{K}p_{in,k}\log\frac{p_{in,k}}{p_{out,k}}$.
In our case, since we only have the samples drawn from the two distributions
$p_{in},p_{out}$, we use an estimator of the KL-divergence based on k-nearst-
neighbor distances [16]. When we use the decomposed discriminative and non-
discriminative features, we can also define
$d_{dis}=D_{KL}[p_{in}(M(z_{d}))\|p_{out}(M(z_{d}))]$ and
$d_{nondis}=D_{KL}[p_{in}(M(z_{n}))\|p_{out}(M(z_{n}))]$. We show $d_{dis}$
and $d_{nondis}$ in Table 2.


# The von Neumann entropy for the Pearson correlation matrix:
A test of the entropic brain hypothesis

H. Felippe1, A. Viol2, D. B. de Araujo3, M. G. E. da Luz4, F. Palhano-Fontes3,
H. Onias3, E. P. Raposo5, G. M. Viswanathan1,6
1Department of Physics, Federal University of Rio Grande do Norte, Natal–RN,
59078-970, Brazil
2Cognitive Neuroscience, Scuola Internazionale Superiore di Studi Avanzati,
Trieste, 34136, Italy
3Brain Institute, Federal University of Rio Grande do Norte, Natal–RN,
59078-970, Brazil
4Departamento de Física, Universidade Federal do Paraná, Curitiba-PR,
81531-980, Brazil
5Laboratório de Física Teórica e Computacional, Departamento de Física,
Universidade Federal de Pernambuco, Recife–PE, 50670-901, Brazil
6National Institute of Science and Technology of Complex Systems, Federal
University of Rio Grande do Norte, Natal–RN, 59078-970, Brazil

(August 27, 2024)

###### Abstract

We address the general problem of how to estimate an entropy given a Pearson
correlation matrix. Most methods currently in use inject a degree of
arbitrariness due to the thresholding of correlations. Here we propose an
entirely objective method of entropy estimation that requires no thresholding.
Let $\bm{R}$ be an $N\times N$ Pearson correlation matrix. We define the
matrix $\bm{\rho}=\bm{R}/N$ and prove that $\bm{\rho}$ satisfies all the
properties required of the density operator. Hence, the von Neumann entropy
$S=-\operatorname{tr}\,(\bm{\rho}\log\bm{\rho})$ can be directly calculated
from the Pearson correlation matrix. To demonstrate the generality and power
of the method, we estimate the entropy of functional correlations of the human
brain. We find that the entropy increases for altered brain states under
serotonergic agonist action, as predicted by the entropic brain hypothesis.

An important open problem in physics concerns how to estimate the entropy from
the Pearson correlation matrix of an arbitrary data set alb02 ; ana09 ; man16
. This problem has wide and immediate applicability to diverse phenomena,
e.g., correlated asset prices alm19 ; cha21 , diffusion processes gom08 ;
hol12 , and brain networks bul09 ; san19 ; nic20 . The Shannon entropy solves
the simpler problem of how to calculate the entropy given a probability
distribution shannon . However, attempts to extend this formalism to estimate
the entropy of a correlation matrix have required admittedly ad hoc methods,
such as the thresholding of the correlation matrix to generate a new symmetric
and unweighted adjacency matrix consisting only of 0’s and 1’s lat01 ; yan18 .
The resulting adjacency matrix represents a simple graph, which can then be
studied using topological metrics. These methods have led to significant
advances in diverse fields, examples of which include applications to market
volatility in financial crisis kum12 , and to scale-free and small-world brain
networks egu05 ; bas06 . Nevertheless, they have drawbacks that motivate the
search for better approaches rub11 .

The main limitation of the existing methods is that the thresholding operation
inevitably leads to loss of information kuk20 . Consider, for example, the set
of all possible $2\times 2$ Pearson correlation matrices. Since a correlation
coefficient can assume any real value in the closed interval $[-1,1]$, there
are an infinite number of possible $2\times 2$ Pearson matrices. In contrast,
there are only two corresponding adjacency matrices of size $2\times 2$
because a simple graph with 2 nodes is either connected or disconnected.
Significant information can thus be lost due to thresholding. Moreover, the
results are not always adequately robust with respect to the choice of
thresholding lan13 ; gar15 . Insofar as unweighted matrices are used, these
methods are only indirectly measuring the entropy of correlations, generally
via topological metrics such as network motifs tag14 , degree distribution
vio17 , and geodesic diversity vio19 . There is thus a growing and vexing
demand for a better method of estimating entropy.

Here, we show that the problem of how to calculate the entropy of a
correlation matrix is akin to one solved long ago by John von Neumann von27
and Lev Landau lan65  for the density operator

$\displaystyle\bm{\rho}=\sum_{i=1}^{M}p_{i}\,|\psi_{i}\rangle\langle\psi_{i}|$
(1)

that describes an ensemble of $M$ state vectors $|\psi_{i}\rangle$ with
assigned classical probabilities $p_{i}$ ohy93 . Given a density operator
$\bm{\rho}$ with $N$ eigenvalues $\lambda_{j}$, the von Neumann entropy is
defined according to

$\displaystyle
S(\bm{\rho})=-\operatorname{tr}\,(\bm{\rho}\log\bm{\rho})=-\sum_{j=1}^{N}\lambda_{j}\log\lambda_{j}\leavevmode\nobreak\
.$ (2)

It is well known that a linear operator $\bm{\rho}$ must satisfy three
necessary and sufficient conditions for the von Neumann entropy to be well
defined: $\bm{\rho}$ must (i) be Hermitian, (ii) have unit trace, and (iii) be
positive semidefinite (i.e., the eigenvalues must be non-negative). Given an
arbitrary Pearson correlation matrix, we will show below how to construct a
corresponding density matrix. This procedure will thereby allow the
calculation of the von Neumann entropy from the Pearson matrix.

Figure 1: The construction of a density matrix from a Pearson correlation
matrix. (a) The system of interest (say, the brain) is partitioned into $N$
interacting components or nodes (regions). The $i$-th and $j$-th regions are
shown highlighted. (b) The correlation between regions $i$ and $j$ is then
calculated. For each $i$ region there are $n$ data points $x_{i}^{(k)}$ where
$k=1,\dots,n$. In the example of the brain that we study here, the
$x_{i}^{(k)}$ represent fluctuations of the blood-oxygen-level-dependent
(BOLD) signal for region $i$. The correlation $R_{ij}$ between regions $i$ and
$j$ is given by Eq. (3), and lies in the interval $-1\leq R_{ij}\leq 1$. (c)
For each correlation coefficient $R_{ij}$, we calculate $\rho_{ij}=R_{ij}/N$,
yielding the matrix $\bm{\rho}=\bm{R}/N$. This matrix $\bm{\rho}$ fully
satisfies the properties of a density operator, i.e., (i) hermiticity, (ii)
unit trace, and (iii) positive semidefiniteness. Finally, we calculate the von
Neumann entropy from this density matrix, via (2). Panel (a) is adapted from
lyn06 .

We first review the definition of the Pearson correlation matrix. Consider $N$
data sets, each with a population of $n$ measurements. Let
$\bm{x}_{i}=(x_{i}^{(1)},\dots,x_{i}^{(n)})$ denote the sequence of $n$
measurements of the $i$-th data set. Let $\langle\leavevmode\nobreak\
\cdot\leavevmode\nobreak\ \rangle$ denote the average taken over the
population of $n$ points, so that
$\bm{X}_{i}=\bm{x}_{i}-\langle\bm{x}_{i}\rangle$ represents deviations from
the mean $\langle\bm{x}_{i}\rangle$ and
$\sigma_{i}^{2}=\langle\bm{X}_{i}^{2}\rangle$ is the variance of $\bm{x}_{i}$.
The correlation between the $i$-th and $j$-th data sets is then given by the
Pearson correlation coefficient

$\displaystyle R_{ij}=\frac{1}{n}\leavevmode\nobreak\
\sum_{k=1}^{n}\left[\frac{{x}_{i}^{(k)}-\langle\bm{x}_{i}\rangle}{\sigma_{i}}\right]\left[\frac{{x}_{j}^{(k)}-\langle\bm{x}_{j}\rangle}{\sigma_{j}}\right]=\frac{\langle\bm{X}_{i}\bm{X}_{j}\rangle}{\sigma_{i}\sigma_{j}}\leavevmode\nobreak\
.$ (3)

The Pearson correlation matrix is the $N\times N$ matrix $\bm{R}$ with entries
$R_{ij}$ given by (3).

Given a Pearson correlation matrix $\bm{R}$, the entropy is often estimated by
using a different (but related) object as a proxy, viz., a binary adjacency
matrix $\bm{A}$. The latter is defined by imposing a threshold value $\xi$
over the entries of $\bm{R}$. Specifically, one sets $A_{ij}=1$ whenever
$|R_{ij}|\geq\xi$ for $i\neq j$, and $A_{ij}=0$ otherwise liu08 . The matrix
$\bm{A}$ defines a graph with $N$ nodes and undirected and unweighted edges.
The entropy estimation, then, becomes a matter of computing the Shannon
entropy of certain probability distributions embedded in the graph structure
of $\bm{A}$, as explained earlier. We note that, in network analysis, the von
Neumann entropy has previously been used with adjacency matrices
passerini-2009 ; estrada-2014 , however it has never been applied directly to
Pearson correlation matrices.

Numerous properties of correlation matrices can be extracted from their
thresholded counterparts, partly due to the information that remains stored in
the topology of simple graphs vio21 . Hence, these methods have become pivotal
and ubiquitous in the study of complex systems kum12 ; egu05 ; bas06 . But
there are yet other drawbacks to this approach, in addition to the limitations
previously mentioned rub11 ; kuk20 ; lan13 ; gar15 . Not only can thresholding
lead to loss of information, it can even inject spurious complexity into
genuinely random networks can20 . For these reasons, such methods require
caution when applied to real-world networks nic20 .

In order to overcome these difficulties, we develop a way to estimate the
entropy directly from the Pearson correlation matrix, without the intermediary
step requiring an adjacency matrix. Our method rests on the following
mathematical foundation: if $\bm{R}$ is an arbitrary $N\times N$ Pearson
correlation matrix, then

$\bm{\rho}={\bm{R}\over N}$ (4)

satisfies all the conditions for a density operator: (i) hermiticity, (ii)
unit trace, and (iii) positive semidefiniteness. To see this, observe that
$\bm{R}$ has real and symmetric entries, from which follows claim (i). Next,
the diagonal entries are given by

$\displaystyle
R_{ii}=\frac{\langle\bm{X}_{i}^{2}\rangle}{\sigma_{i}^{2}}=1\leavevmode\nobreak\
,$ (5)

so that $\operatorname{tr}\bm{R}=N$, leading to claim (ii). For (iii), it
suffices to show that $\bm{R}$ is positive semidefinite. In other words, we
must show that $\bm{v}\cdot\bm{R}\,\bm{v}\geq 0$ for all
$\bm{v}\in\mathbb{R}^{N}$. In fact, the positive semidefiniteness of Pearson
correlation matrices is well known. Observe that

$\displaystyle\sum_{ij}^{N}v_{i}\left(\sum_{k}^{n}X_{i}^{(k)}X_{j}^{(k)}\right)v_{j}=\sum_{k}^{n}\left|\sum_{i}^{N}v_{i}X_{i}^{(k)}\right|^{2}\geq
0\leavevmode\nobreak\ .$ (6)

We now state our main result. Given an arbitrary Pearson correlation matrix,
the entropy can be calculated from Eq. (2), with $\bm{\rho}$ given by (4). We
illustrate our general method in Figure 1. We remark in passing that many
other correlation matrices (e.g., Spearman, Kendall) will satisfy analogous
conditions, so long as they consist of real symmetric entries bou12 .

Next we give illustrative examples of maximal and minimal entropies for $N$
data sets. Consider first the scenario where all $\bm{x}_{i}$ and $\bm{x}_{j}$
are statistically independent, i.e., uncorrelated:
$\langle\bm{X}_{i}\bm{X}_{j}\rangle=\langle\bm{X}_{i}^{2}\rangle\,\delta_{ij}$.
Thus $R_{ij}=\delta_{ij}$, and, as expected, the density operator admits the
(completely mixed state) representation $\bm{\rho}=\bm{I}/N$, where $\bm{I}$
is the $N\times N$ identity matrix. The von Neumann entropy then assumes the
maximal value of $S=\log N$, as it should.

Compare the above with the other extreme, i.e., when all pairs are perfectly
autocorrelated, so that $\bm{X}_{i}/\sigma_{i}=\bm{X}_{j}/\sigma_{j}$ for all
$i,\leavevmode\nobreak\ j$. Then, clearly, $R_{ij}=1$, so that
$\bm{R}^{2}=N\bm{R}$. The density operator is thus idempotent:
$\bm{\rho}^{2}=\bm{\rho}$. It follows immediately that

$\operatorname{tr}\,(\bm{\rho}\log\bm{\rho})=\operatorname{tr}\,(\bm{\rho}\log\bm{\rho^{2}})=2\operatorname{tr}\,(\bm{\rho}\log\bm{\rho})\leavevmode\nobreak\
,$

hence $\operatorname{tr}\,(\bm{\rho}\log\bm{\rho})=0$. Therefore, the von
Neumann entropy assumes $S=0$, which is the absolute minimum possible value
(corresponding to a pure state configuration).

We thus have a robust and direct approach to obtain the entropy of correlation
matrices. Crucially, at no point have we introduced a priori estimates, hence
the resulting entropy is free from arbitrary or extraneous assumptions (e.g.,
choice of threshold $\xi$). Of course, one can additionally use thresholding
to obtain a suitable $\bm{A}$, but even then our protocol can still be useful.
For example, it can be used to obtain an exact expression for the entropy of
an idealized ranked set of correlated data (see Supplemental Material supp ).
Such results might then yield a null hypothesis for subsequent statistical
inference or to undertake more rigorous hypothesis testing.

As a practical demonstratation of the method’s power and wide applicability,
we calculate the entropy of Pearson correlation matrices obtained from
functional magnetic resonance imaging (fMRI) of the human brain smi12 .
Specifically, we test the hypothesis of increased brain entropy while under
the influence of a psychedelic car14 . For this purpose, we chose the
indigenous beverage ayahuasca lun11 , a decoction with rapid antidepressant
effects partially mediated by the serotonergic agonist N,N-dimethyltryptamine
(DMT) rib03 ; oso15 ; pal19 .

A total of 9 healthy human volunteers were assigned to a neuroimaging session
both before and after ayahuasca intake. We have a set of $N=104$ fMRI time
series of the BOLD signal of the human brain from each subject at both
conditions. See, again, the Supplemental Material supp  for details on the
method of parcellating the brain into regions of interest, of acquiring the
fMRI images, and of obtaining the BOLD time series.

Figure 2: Functional connectivity of the brain of a volunteer (Subject 4) in
the ayahuasca condition. (a) Matrix with entries $\rho_{ij}=R_{ij}/N$, where
$R_{ij}$ is the Pearson correlation of the $i$-th and $j$-th BOLD fMRI time
series. (b) The matrix $\bm{\rho}\log\bm{\rho}$. Up to a sign, the trace of
$\bm{\rho}\log\bm{\rho}$ gives the von Neumann entropy of $\bm{\rho}$.

We next calculated the density operators. The $18$ correlation matrices were
individually scaled according to Eq. (4) using $N=104$. Figure 2(a) shows an
example of the functional connectivity network (represented as a matrix) from
a single subject in the ayahuasca condition. Every matrix entry $\rho_{ij}$ is
shown $(i,j=1,\dots,104)$, expressing the (scaled) correlation between pairs
of time series $\bm{x}_{i}$ and $\bm{x}_{j}$. In 2(b), the matrix
$\bm{\rho}\log\bm{\rho}$ is also shown in its entirety. Up to a sign, the sum
of its diagonal entries is precisely the von Neumann entropy of $\bm{\rho}$.

Figure 3 shows the normalized entropy of all $18$ matrices $\bm{\rho}$,
grouped in terms of subjects (from 1 to 9) for both experimental conditions
(before and after ayahuasca ingestion). With the exception of 2 individuals
(Subjects 1 and 9), we find positive entropy variation $\Delta S\equiv
S_{\textrm{after}}-S_{\textrm{before}}$ among subjects (dotted lines).
However, the low number of individuals demands careful analysis in regards to
statistical significance.

Figure 3:  Entropy before and after ayahuasca. The plot shows the von Neumann
entropy $S$ of the 9 subjects (horizontal axis), both before
($S_{\textrm{before}}$) and after ($S_{\textrm{after}}$) ayahuasca ingestion,
normalized by $S_{\textrm{max}}=\log\left(104\right)$. Entropy variation
$\Delta S=S_{\textrm{after}}-S_{\textrm{before}}$ was found positive for 7
subjects (dotted lines). These findings are broadly consistent with the
entropic brain hypothesis, although one must be careful due to the small
number of subjects involved.

The relevance of our method’s findings is best seen through the lens of the
entropic brain hypothesis (EBH) car14 , a recent approach to cognitive
neuroscience grounded on notions of criticality kin06 ; chi10  and prompted by
neuroimaging research with psychedelic compounds (e.g., psilocybin, LSD) dra12
; car12 ; car16 . The EBH asserts that the quality of conscious states is
dependent upon the entropy of key parameters of brain function. Specifically,
it predicts an increase of entropy when the brain is under the influence of a
psychedelic substance.

Through various indirect measures of spontaneous brain activity, the EBH was
shown consistent with previously conducted works, for instance, the enhanced
repertoire of network motifs post-psilocybin infusion tag14 , the greater
sample entropy (negative predictability) of fMRI time courses following LSD
administration leb16 , as well as with the broadening of functional patterns
due to ayahuasca intake vio17 ; vio19 . In this work, similarly, we found
increases in brain entropy among individuals in the ayahuasca condition.
Specifically, the entropy of correlations embedded in the functional
connectivity of the human cortex (i.e., correlated BOLD fMRI time series)
increases in 7 out of 9 subjects undergoing the effects of ayahuasca. Our test
of the entropic brain hypothesis requires careful interpretation because of
the small sample size. Indeed, the main goal of this work is to show the
generality and usefulness of the proposed method based on the von Neumann
entropy, as a proof of concept. It should be clear to the reader that the
proposed method can be used to perform rigorous hypothesis testing on
arbitrarily large data sets, for example.

From a methodological standpoint, it is instructive to compare our results
with those of Viol et al. in vio17  and vio19 : following the ayahuasca
consumption, the latter found an enhanced neighborhood diversity (as measured
by the geodesic entropy) of the nodes within functional networks; whereas the
former measured a greater entropy of the network’s degree distribution.
Crucially, both works applied a range of threshold values to examine the
properties of numerous adjacency matrices. This was done in an effort to
compensate for the inherent limitations of the thresholding procedure. In
contrast, the method here proposed uses no thresholding in its entropy
calculation, thus evading the drawbacks associated with thresholding rub11 ;
kuk20 ; lan13 ; gar15 ; can20 . Calculation of the von Neumann entropy is also
fast, for it uses a single correlation matrix instead of many adjacency
matrices. Moreover, the von Neumann entropy has a solid mathematical
foundation ohy93 .

The concept of entropy has found growing importance within the general area of
biomedicine, particularly in neuroscience. Indeed, entropic features of the
brain activity have recently served as the basis for new models of brain
function and structure car19 ; her20 ; san21 . On such grounds, the present
protocol should be considered a relevant new addition to the existing tools to
analyze functional correlations in the human brain. For instance, it has
corroborated the view that psychedelics induce brain entropy increase.
Finally, in a broader context we believe our approach can be successfully
applied to many other complex systems, such as risk management in financial
markets nam11 , and the dynamics of disease spreading and information routing
boc06 ; iacopini18 .

We thank M. Copelli and S. Ribeiro for discussions. AV thanks the support of
PRIN grant 20174TPEFJ “TRIPS”. This work was supported by the Brazilian
agencies CNPq (Grants No. 304532/2019-3, No. 305062/2017-4 and No.
302051/2018-0), CAPES, and FACEPE.

## References

  * (1) R. Albert and A. Barabási, Rev. Mod. Phys. 74, 47 (2002).
  * (2) K. Anand and G. Bianconi, Phys. Rev. E 80, 045102(R) (2009).
  * (3) M. De Domenico and J. Biamonte, Phys. Rev. X 6, 041062 (2016).
  * (4) A. Almog and E. Shmueli, Sci. Rep. 9 10832 (2019).
  * (5) A. Chakraborti, Hrishidev, K. Sharma, and H. K. Pharasi, J. Phys. Complex. 2, 015002 (2021).
  * (6) J. Gómez-Gardeñes and V. Latora, Phys. Rev. E 78, 065102(R) (2008).
  * (7) P. Holme and J. Saramäki, Phys. Rep. 519, 97 (2012).
  * (8) E. Bullmore and O. Sporns, Nat. Rev. Neurosci. 10, 186 (2009).
  * (9) F. A. N. Santos, E. P. Raposo, M. D. Coutinho-Filho, M. Copelli, C. J. Stam, and L. Douw, Physical Review E 100, 032414 (2019)
  * (10) C. Nicolini, G. Forcellini, L. Minati, and A. Bifone, Neuroimage 211, 116603 (2020).
  * (11) C. Shannon, Bell System Tech. J. 27 623 (1948).
  * (12) V. Latora and M. Marchiori, Phys. Rev. Lett. 87, 198701 (2001).
  * (13) X. Yan, L. G. S. Jeub, A. Flammini, F. Radicchi, and S. Fortunato, Phys. Rev. E 98, 042304 (2018).
  * (14) S. Kumar and N. Deo, Phys. Rev. E 86, 026101 (2012).
  * (15) V. M. Eguíluz, D. R. Chialvo, G. A. Cecchi, M. Baliki, and A. V. Apkarian, Phys. Rev. Lett. 94, 018102 (2005).
  * (16) D. S. Bassett and E. Bullmore, The Neuroscientist 12, 512 (2006).
  * (17) M. Rubinov and O. Sporns, Neuroimage 56, 2068 (2011).
  * (18) V. Kukreti, H. K. Pharasi, P. Gupta, and S. Kumar, Front. Phys. 8, 323 (2020).
  * (19) N. Langer, A. Pedroni, and L. Jäncke, PLOS ONE 8, 1 (2013).
  * (20) K. A. Garrisona, D. Scheinostb, E. S. Finnc, X. Shenb, R. T. Constablebd, Neuroimage 118, 651 (2015).
  * (21) E. Tagliazucchi, R. L. Carhart-Harris, R. Leech, D. Nutt, and D. R. Chialvo, Hum. Brain. Mapp. 35, 5442 (2014).
  * (22) A. Viol, F. Palhano-Fontes, H. Onias, D. B. de Araujo, and G. M. Viswanathan, Sci. Rep. 7, 7388 (2017).
  * (23) A. Viol, F. Palhano-Fontes, H. Onias, D. B. de Araujo, P. Hövel, and G. M. Viswanathan, Entropy 21, 128 (2019).
  * (24) J. von Neumann, Thermodynamik quantenmechanischer Gesamtheiten, Nachr. Ges Wiss. Göttingen (1927).
  * (25) L. D. Landau, The Damping Problem in Wave Mechanics, in Collected Papers of L.D. Landau, edited by D. T. Haar (Pergamon, 1965), pp. 8-18.
  * (26) M. Ohya and D. Petz, Quantum Entropy and Its Use (Springer-Verlag, Berlin, 1993).
  * (27) Y. Liu, M. Liang, Y. Zhou, Y. He, Y. Hao, M. Song, C. Yu, H. Liu, Z. Liu, and T. Jiang, Brain 131, 945 (2008)
  * (28) F. Passerini and S. Severini,  Int. J. Agent Tech. Syst. 1 58 (2009).
  * (29) E. Estrada, J. A. de la Pẽna, and N. Hatano,  Linear Algebra Appl. 443, 235 (2014).
  * (30) A. Viol, V. Vuksanović, and P. Hövel, Physica A 561, 125233 (2021).
  * (31) G. T. Cantwell, Y. Liu, B. F. Maier, A. C. Schwarze, C. A. Serván, J. Snyder, and G. St-Onge, Phys. Rev. E 101, 062302 (2020).
  * (32) P. J. Lynch and C. C. Jaffe, CC BY 2.5 (2006).
  * (33) K. Boudt, J. Cornelissen, and C. Croux, Stat. Comput. 22, 471 (2012).
  * (34) See Supplemental Material at `PRL-link` for a detailed description of the experimental procedure, data acquisition, and mathematical analysis. It includes Refs. [35–43] and [50–54].
  * (35) L. Guttman, The quantification of a class of attributes: a theory and method of scale construction, in The prediction of personal adjustment, edited by P. Horst, P. Wallin, L. Guttman, F. B. Wallin, J. A. Clausen, R. Reed, and E. Rosenthat (Social Science Research Council, New York, 1941), pp. 319-348.
  * (36) L. Guttman, The principal components of scale analysis, in Measurement and prediction, edited by S. A. Stouffer, L. Guttman, E. A. Suchman, P. F. Lazarsfeld, S. A. Star, and J. A. Clausen (Princeton University Press, Princeton, 1950), pp. 312-361.
  * (37) R. Zwick, Psychometrika 52, 512 (1987).
  * (38) C. P. Davis-Stober, J.-P. Doignon, and R. Suck, Front. Psychol. 6, 1767 (2015).
  * (39) M. Kendall, Rank Correlation Methods, 5th ed. (Oxford University Press, Oxford, 1990).
  * (40) L. Piao and Z. Fu, Sci. Rep. 6, 36759 (2016).
  * (41) L. A. Goodman, J. Am. Stat. Assoc. 70, 755 (1975).
  * (42) D. J. Hand, J. R. Statist. Soc. A 159, 445 (1996).
  * (43) C. H. Coombs, L. C. Coombs, and J. C. Lingoes, Stochastic cumulative scales, Theory construction and data analysis in the behavioral sciences, edited by S. Shye (Jossey-Bass, San Francisco, 1978), pp. 280–298.
  * (44) S. M. Smith, Neuroimage 62, 1257 (2012).
  * (45) R. L. Carhart-Harris, R. Leech, P. J. Hellyer, M. Shanahan, A. Feilding, E. Tagliazucchi, D. R. Chialvo, and D. Nutt, Front. Hum. Neurosci. 8, 20 (2014).
  * (46) L. E. Luna, Indigenous and Mestizo use of Ayahuasca. An Overview, in The Ethnopharmacology of Ayahuasca, edited by R. G. dos Santos (Transworld Research Network, Tivandrum, 2011), pp 1-21.
  * (47) J. Riba, M. Valle, G. Urbano, M. Yritia, A. Morte, and M. J. Barbanoj, J. Pharmacol. Exp. Ther. 306, 73 (2003).
  * (48) F. L. Osório, R. F. Sanches, L. R. Macedo, R. G. dos Santos, J. P. Maia-de-Oliveira, L. Wichert-Ana, D. B. de Araujo, J. Riba, J. A. Crippa, and J. E. Hallak, Rev. Bras. Psiquiatr. 37, 13 (2015).
  * (49) F. Palhano-Fontes, D. Barreto, H. Onias, K. C. Andrade, M. M. Novaes, J. A. Pessoa, S. A. Mota-Rolim, F. L. Osório, R. Sanches, R. G. dos Santos et al., Psychol. Med. 49, 655 (2019).
  * (50) A. P. Association, Diagnostic and Statistical Manual of Mental Disorders: DSM-IV-TR (American Psychiatric Association, 2000).
  * (51) D. B. de Araujo, S. Ribeiro, G. A. Cecchi, F. M. Carvalho, T. A. Sanchez, J. P. Pinto, B. S. de Martinis, J. A. Crippa, J. E. C. Hallak, and A. C. Santos, Hum. Brain Mapp. 33, 255 (2012).
  * (52) M. Jenkinson, C. F. Beckmann, T. E. Behrens, M. W. Woolrich, and S. M. Smith, Neuroimage 62 782 (2012).
  * (53) M. Brett, I. S. Johnsrude, and A. M. Owen, Nat. Rev. Neurosci. 3, 243 (2002).
  * (54) H. Onias, A. Viol, F. Palhano-Fontes, K. C. Andrade, M. Sturzbecher, G. Viswanathan, and D. B. de Araujo, Epilepsy Behav. 38, 71 (2014).
  * (55) O. Kinouchi and M. Copelli, Nat. Phy. 2, 348 (2006).
  * (56) D. R. Chialvo, Nat. Phys. 6, 744 (2010).
  * (57) R. L. Carhart-Harris, D. Erritzoe, T. Williams, J. M. Stone, L. J. Reed, A. Colasanti, R. J. Tyacke, R. Leech, A. L. Malizia, K. Murphy et al., PNAS 109, 2138 (2012).
  * (58) R. L. Carhart-Harris, S. Muthukumaraswamy, L. Roseman, M. Kaelen, W. Droog, K. Murphy, E. Tagliazucchi, E. E. Schenberg, T. Nest, C. Orban et al., PNAS 113, 4853 (2016).
  * (59) A. V. Lebedev, M. Kaelen, M. Lövdén, J. Nilsson, A. Feilding, D. J. Nutt, and R. L. Carhart-Harris, Hum. Brain. Mapp. 37, 3203 (2016).
  * (60) R. L. Carhart-Harris and K. J. Friston, Neuropharmacology 71, 316 (2019).
  * (61) R. Herzog, P. A. M. Mediano, F. E. Rosas, R. Carhart-Harris, Y. S. Perl, E. Tagliazucchi, and R. Cofre, Sci. Rep. 10, 17725 (2020).
  * (62) C. Sanz, C. Pallavicini, F. Carrillo, F. Zamberlan, M. Sigman, N. Mota, M. Copelli, S. Ribeiro, D. Nutt, R. Carhart-Harris, and E. Tagliazucchi, Conscious. Cogn. 87, 103070 (2021).
  * (63) A. Namaki, A. Shirazi, R. Raei, and G. Jafari, Physica A 390, 3835 (2011).
  * (64) S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D. -U. Hwang, Phys. Rep. 424, 175 (2006).
  * (65) I. Iacopini, S. Milojević, and V. Latora, Phys. Rev. Lett. 120, 048301 (2018).

Supplementary Material for:
The von Neumann entropy for the Pearson correlation matrix:
A test of the entropic brain hypothesis
H. Felippe1, A. Viol2, D. B. de Araujo3, M. G. E. da Luz4, F. Palhano-Fontes3,
H. Onias3, E. P. Raposo5, G. M. Viswanathan1,6

1Department of Physics, Federal University of Rio Grande do Norte, Natal–RN,
59078-970, Brazil
2Cognitive Neuroscience, Scuola Internazionale Superiore di Studi Avanzati,
Trieste, 34136, Italy
3Brain Institute, Federal University of Rio Grande do Norte, Natal–RN,
59078-970, Brazil
4Departamento de Física, Universidade Federal do Paraná, Curitiba-PR,
81531-980, Brazil
5Laboratório de Física Teórica e Computacional, Departamento de Física,
Universidade Federal de Pernambuco, Recife–PE, 50670-901, Brazil
6National Institute of Science and Technology of Complex Systems, Federal
University of Rio Grande do Norte, Natal–RN, 59078-970, Brazil
(Dated: August 27, 2024)

##  I I. Guttman Scaling

###  I.1 A. The von Neumann entropy for the Homogeneous Guttman-scalable
Pearson correlation matrix

Here we discuss an idealized type of data, following the Guttman scaling
guttman-1941 ; guttman-1950 . In the homogeneous case zwick-1987  (see below),
many exact results are known for the corresponding Pearson correlation matrix
$\bm{R}_{hG}$ zwick-1987 ; stober-2015 . As we will show, this allows to write
an analytical expression for the von Neumann entropy $S_{hG}$ of
$\bm{R}_{hG}$. Although real data rarely fits in with a perfect Guttman scale,
the present $S_{hG}$ could play the role of a null hypothesis model zwick-1987
. Indeed, it may be useful for testing hierarchies as well as features of
correlations in distinct processes involving $N$ populations kendall-1990 ;
piao-2016 ; in our case all with the same size $n$ and described by
($i=1,\ldots,N$)

${\bm{x}}_{i}=(x_{i}^{(1)},\ldots,x_{i}^{(n)})\leavevmode\nobreak\ .$ (S1)

Moreover, Guttman scaling is a very instructive example of how the von Neumann
entropy portrays data with strong but not absolute correlations.

Next we address (homogeneous) Guttman-scalable data as usually done in the
literature zwick-1987 ; goodman-1975 . Possible ways to relate such standard
framework to data in the form of Eq. (S1) is briefly considered in Sec. IB.
For our purposes, we concentrate on dichotomous or binary events (or “items”,
a more common terminology in the area of statistical measurement theory
zwick-1987 ; hand-1996 ). Thus, to each properly defined event $i$ we can
associate two distinct numerical values, 0 and 1. For a given instance or
realization pattern for the collection of events, the numerical value 1 (0)
for the event $i$ means it has (has not) taken place in that pattern.

The idea underlying the Guttman scaling, also known as cumulative scaling or
scalogram analysis, is to be able to create a table or indicator matrix —
which should not be confused with the correlation matrix $\bm{R}$ itself —
assembling the data in a very particular manner. If such construction is
possible for the data set, the indicator matrix would display a triangular
structure, with all entries in the lower (upper) part equal to 1 (0), see
Table 1. Hence, within a given pattern of outcomes, there is a priority of
events. If the event $i$ has occurred, then event $i-1$ has also occurred. The
Table 1 explains in more detail the data organization.

Figure S1: The von Neumann entropy, $S_{hG}$, for a Pearson correlation matrix resulting from an ideal homogeneous Guttman-scalable data set versus its population (or items) number $N$. The entropy $S_{hG}$ tends to a constant value for $N$ going to infinity (see the inset). For comparison, it also shows the von Neumann entropy for a completely uncorrelated data set, for which $S=\ln[N]$. Table 1:  The data table (or indicator matrix) organized in a specific ordering of events $i$, with $i=1,\ldots,N$. Along each row — representing a particular configuration of events occurrence: a pattern — the entry 1 (0) means that the corresponding event has (has not) taken place. The depicted structure illustrates a perfect Guttman scaling: for the pattern $i$, only the events up to $i-1$ do take place. This ascribes a ranking or degree of prevalence among the distinct events within a specific pattern. ${\mathcal{N}}_{i}$ is the number of times the pattern $i$ is found in the full data set collection. Pattern | Event 1 | Event 2 | Event 3 | Event 4 | Event 5 | … | Event $N$ | Number of occurrences
---|---|---|---|---|---|---|---|---
1 | 0 | 0 | 0 | 0 | 0 | … | 0 | ${\mathcal{N}}_{1}$
2 | 1 | 0 | 0 | 0 | 0 | … | 0 | ${\mathcal{N}}_{2}$
3 | 1 | 1 | 0 | 0 | 0 | … | 0 | ${\mathcal{N}}_{3}$
4 | 1 | 1 | 1 | 0 | 0 | … | 0 | ${\mathcal{N}}_{4}$
5 | 1 | 1 | 1 | 1 | 0 | … | 0 | ${\mathcal{N}}_{5}$
6 | 1 | 1 | 1 | 1 | 1 | … | 0 | ${\mathcal{N}}_{6}$
⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | 0 | ⋮
$N+1$ | 1 | 1 | 1 | 1 | 1 | … | 1 | ${\mathcal{N}}_{N+1}$

Two extra assumptions for data already observing the Guttman arrangement are
known as homogeneity conditions zwick-1987 . They can be stated as the
following. (a) There is no full correlation between any two populations
${\bm{x}}_{i^{\prime\prime}}$ and ${\bm{x}}_{i^{\prime}}$. Actually, in many
contexts this is not a really restrictive hypothesis, either because such
coincidence is rare or, if it is not, for $N$ large enough we can just
eliminate the akin ${\bm{x}}_{i}$’s from the data set without loosing relevant
statistical significance. (b) The probability of obtaining the different
patterns of outcomes indicated in Table 1 (see also Sec. IB) are the same,
i.e., the ${\mathcal{N}}_{i}$’s in Table 1 are essentially all equal.
Arguably, this may seem a too restricting imposition stober-2015 . But perhaps
a bit surprising, it has been verified numerically that certain results (like
principal component analysis and certain properties of covariance matrices)
for homogeneous data do not depart very much from the case where the
${\mathcal{N}}_{i}$’s display a normal distribution (provided the populations
have ample variability zwick-1987 ; stober-2015 ).

Therefore, it is correct to affirm that data obeying the homogeneous Guttman
scaling is not ubiquitous. But the relevant fact for our purposes is that,
precisely because of that, the Pearson correlation matrix $\bm{R}_{hG}$ can be
derived analytically zwick-1987 . Furthermore, its exact eigenvalues
$\tilde{\lambda}_{i}$’s have been obtained in zwick-1987 . Since the
eigenvalues ${\lambda}_{i}$’s of $\bm{\rho}_{hG}$ are simply
$\tilde{\lambda}_{i}/N$, then from zwick-1987  we get

$\lambda_{i}=\frac{1+N^{-1}}{i\,(i+1)}\leavevmode\nobreak\ .$ (S2)

Notice that the $\lambda_{i}$’s depend only on a single parameter, the
population number $N$. So, finally we find

$\displaystyle S_{hG}$
$\displaystyle=-(1+N^{-1})\sum_{i=1}^{N}\,\frac{1}{i\,(i+1)}\ln\left[\frac{1+N^{-1}}{i\,(i+1)}\right]$
$\displaystyle=\ln[N]-(1-N^{-2})\ln[N+1]$
$\displaystyle\quad+2\,(1+N^{-1})\sum_{i=1}^{N-1}\frac{1}{i\,(i+2)}\ln[i+1]\leavevmode\nobreak\
.$ (S3)

Since it is well known that the series $\sum_{i=1}^{\infty}\ln[i]/i^{2}$
converges, it is not difficult to prove that $S_{hG}$ also converges as
$N\rightarrow\infty$. Also, due to the peculiar way items are correlated in a
homogeneous Guttman-scalable data collection, after a rapid increase with $N$,
the entropy $S_{hG}$ tends to a plateau as $N$ increases. Figure S1
illustrates the general behavior of $S_{hG}(N)$, also comparing it with the
von Neumann entropy for a Pearson correlation matrix of a totally uncorrelated
data set, for which $S=\ln[N]$.

Table 2:  Consider that for a data set of $N=3$ and $n=4$, the twelve $x_{i}^{(j)}$’s assume $M_{0}=6$, $M_{1}=3$, $M_{2}=2$, $M_{3}=1$ values, respectively, in the numerical intervals $\Delta_{0}$, $\Delta_{1}$, $\Delta_{2}$ and $\Delta_{3}$. Then, taking as the events the instances (irrespective of $i$ and $j$) in which the variables are within $\Delta_{1}$, $\Delta_{2}$ and $\Delta_{3}$, one can construct an indicator matrix exactly satisfying a homogeneous Guttman scaling. Indeed, the three, two and one “1”’s in the columns Event 1, Event 2 and Event 3, correspond to $M_{1}$, $M_{2}$ and $M_{3}$, whereas the total of six “0”’s spread in the whole table are due to $M_{0}$. Pattern | Event 1: Interval $\Delta_{1}$ | Event 2: Interval $\Delta_{2}$ | Event 3: Interval $\Delta_{3}$ | Occurrence of the pattern
---|---|---|---|---
1 | 0 | 0 | 0 | 1
2 | 1 | 0 | 0 | 1
3 | 1 | 1 | 0 | 1
4 | 1 | 1 | 1 | 1

###  I.2 B. Simple possible ways to frame a general data set in terms of
“events” for the Guttman scaling description

An important question is certainly how a given data set, given as in Eq. (S1),
can be described in terms of events and thus to allow an indicator matrix
representation. Below we summarize two possible procedures to do so.

The first is when we can associate phenomenological behavior to the data.
Then, we can define an event by classifying appropriate (i.e., easier to
characterize in terms of ranking of predominance) information which can be
inferred once we know the explicit values of the $x_{i}^{(j)}$’s. For example,
suppose whenever $x_{i}^{(j)}$ is above a threshold $\gamma$, we can conclude
it has been a high electrical pulse crossing the population $i$ at the
measurement step $j$. Then, an event could be the existence (1) or not (0) of
an intense pulse. Further, suppose that a determined biochemical reaction
takes place when $\alpha<x_{i}^{(j)}<\beta$. Such reaction could also
constitute another event, with the values 0 and 1 being ascribed conforming to
the numerical intervals for the $x_{i}^{(j)}$’s. Importantly, if
$\langle{\bm{x}}_{i}\rangle<\gamma<\alpha$, in principle we may expect an
order of priority between the high pulse and the biochemical reaction. In this
way, we might try to organize the events as in Table 1 according to the
numerical sequences in Eq. (S1). Of course, such an approach is heavily based
on the concrete understanding of the physical, chemical, biological, etc.,
aspects underlying the problem studied. In the absence of this type of
knowledge, the protocol is not feasible.

A second possibility is to consider an algorithmically oriented method, not
relying on any a priori qualitative information about the process. But then,
admittedly it has a higher chance of not satisfying the perfect homogeneous
Guttman scaling stober-2015 ; coombs-1978 . For the full data set, let us
define $a$ and $b$ as the minimum and maximum values assumed by the entire
collection of $x_{i}^{(j)}$’s. Thus, we divide $[a,b]$ into $N+1$ intervals
$\Delta_{0},\Delta_{1},\ldots,\Delta_{N}$. Here, $\Delta_{0}$ must represent
the most frequent numerical interval for the data set, with the other
$\Delta_{i}$’s presenting a decreasing order of occurrence with $i$. The
events to be considered in the indicator matrix are the occurrence of the
intervals $\Delta_{1},\Delta_{2},\ldots,\Delta_{N}$ ($\Delta_{0}$ is not
explicitly included in the events table, but it contributes with the “0”’s in
it). For $M_{i}$ the number of times the interval $\Delta_{i}$ appears in the
data set, we have $\sum_{i=0}^{N}M_{i}=n\times N$. All the incidences
$M_{i}$’s ($i=0,1,\ldots,N$) of the events $\Delta_{i}$’s should thus be
combined to generate a structure like that in Table 1. An illustration of how
to proceed in the particular case the data allows an exact homogeneous Guttman
scaling is depicted in Table 2. Finally, given that there is a certain freedom
in choosing the intervals $\Delta_{i}$’s, this must be done in such way to get
an indicator matrix as close as possible to the pattern in Table 1.

##  II II. Data Acquisition

A total of 9 healthy adults (5 women) with no history of neurological or
psychiatric disorders (assessed by DSM-IV structured interview dsm ),
volunteered to ingest 120–200 mL (2.2 mL/kg of body weight) of ayahuasca known
to contain 0.8 mg/mL of its main psychoactive compound N,N-dimethyltryptamine
(DMT) and 0.21 mg/mL of harmine, a $\beta$-carboline which allows for the
entrance of DMT into the bloodstream dra12 . The volunteers were assigned with
two task-free fMRI sessions: the first one prior to the ayahuasca ingestion,
and the second one at 40 minutes after the brew intake (when the acute effects
of ayahuasca start to emerge, lasting for approximately 4 hours). Both
sessions required the participants to maintain an awake resting state. The
experimental procedure was approved by the Ethics and Research Committee of
the University of São Paulo at Ribeirão Preto (No. 14672/2006).

The fMRI images were acquired in a 1.5 T scanner (Siemens, Magneton Vision),
using an EPI-BOLD sequence to obtain 150 volumes with parameters TR = 1700 ms,
TE = 66 ms, FOV = 220 mm, matrix $64\times 64$, and voxel dimensions of 1.72
mm $\times$ 1.72 mm $\times$ 1.72 mm. The images were preprocessed in the FSL
software jen12 , and consisted of slice-timing, head motion, and spatial
smoothing corrections (Gaussian kernel, FWHM = 5 mm). A total of 9 regressors
were used within a general linear model (GLM): 6 regressors to movement
correction; 1 to white matter signal; 1 to cerebrospinal fluid; and 1 to
global signal. The images were normalized to the standard anatomical space of
the Montreal Neurological Institute (MNI152 template) bret02 .

The preprocessed fMRI images were then parcellated into 110 anatomical regions
of interest (ROIs) in accordance with the Harvard-Oxford cortical and
subcortical atlas (see Table 3). Due to acquisition limitations, six regions
were excluded from further analysis. Thus, for each one of the remaining 104
ROIs, we averaged the BOLD signal of the regions’s associated voxels,
resulting in time courses listing a sequence of $n=150$ data points. To reduce
confounders in the signal, a maximum overlap discrete wavelet transform
(MODWT) was applied to the series in order to select the typical frequency
range (0.01–0.1 Hz) of the resting state signal oni14 . Finally, in possession
of the $N=104$ time series $\bm{x}_{i}$, we pairwise-correlated them using the
Pearson correlation coefficient of Eq. (3) from the main text, resulting in
the $104\times 104$ Pearson correlation matrix $\bm{R}$. We thus ended up with
a total of $18$ functional connectivity networks: two conditions (before and
after ayahuasca) for each of the 9 subjects involved.

Table 3: Brain regions in accordance with the Harvard-Oxford cortical atlas. Asterisks (∗) are assigned to the six excluded regions. Ordered pairs $(\ell,r)$ refers to the left ($\ell$) and right ($r$) hemispheres of the corresponding region. For example, in reference to Eq. (3) of the main text, the correlation between the left Thalamus ($\ell=97$) and right Hippocampus ($r=108$) is given by the Pearson coefficient $R_{97,108}$. Index | Region
---|---
$(1,2)$ | Frontal Pole
$(3,4)$ | Insular Cortex
$(5,6)$ | Superior Frontal Gyrus
$(7,8)$ | Middle Frontal Gyrus
$(9,10)$ | Inferior Frontal Gyrus, pars triangularis
$(11,12)$ | Inferior Frontal Gyrus, pars opercularis
$(13,14)$ | Precentral Gyrus
$(15,16)$ | Temporal Pole
$(17,18)$ | Superior Temporal Gyrus, anterior division
$(19,20)$ | Superior Temporal Gyrus, posterior division
$(21,22)$ | Middle Temporal Gyrus, anterior division
$(23,24)$ | Middle Temporal Gyrus, posterior division
$(25,26)$ | Middle Temporal Gyrus, temporooccipital part
$(27,28)$ | Inferior Temporal Gyrus, anterior division*
$(29,30)$ | Inferior Temporal Gyrus, posterior division
$(31,32)$ | Inferior Temporal Gyrus, temporooccipital part
$(33,34)$ | Postcentral Gyrus
$(35,36)$ | Superior Parietal Lobule*
$(37,38)$ | Supramarginal Gyrus, anterior division
$(39,40)$ | Supramarginal Gyrus, posterior division
$(41,42)$ | Angular Gyrus
$(43,44)$ | Lateral Occipital Cortex, superior division
$(45,46)$ | Lateral Occipital Cortex, inferior division
$(47,48)$ | Intracalcarine Cortex
$(49,50)$ | Frontal Medial Cortex
$(51,52)$ |  Juxtapositional lobule cortex (form. supp. motor cortex)
$(53,54)$ | Subcallosal Cortex
$(55,56)$ | Paracingulate Gyrus
$(57,58)$ | Cingulate Gyrus, anterior division
$(59,60)$ | Cingulate Gyrus, posterior division
$(61,62)$ | Precuneous Cortex
$(63,64)$ | Cuneal Cortex
$(65,66)$ | Frontal Orbital Cortex
$(67,68)$ | Parahippocampal Gyrus, anterior division
$(69,70)$ | Parahippocampal Gyrus, posterior division
$(71,72)$ | Lingual Gyrus
$(73,74)$ | Temporal Fusiform Cortex, anterior division*
$(75,76)$ | Temporal Fusiform Cortex, posterior division
$(77,78)$ | Temporal Occipital Fusiform Cortex
$(79,80)$ | Occipital Fusiform Gyrus
$(81,82)$ | Frontal Operculum Cortex
$(83,84)$ | Central Opercular Cortex
$(85,86)$ | Parietal Operculum Cortex
$(87,88)$ | Planum Polare
$(89,90)$ | Heschl’s Gyrus (includes H1 and H2)
$(91,92)$ | Planum Temporale
$(93,94)$ | Supracalcarine Cortex
$(95,96)$ | Occipital Pole
$(97,104)$ | Thalamus
$(98,105)$ | Caudate
$(99,106)$ | Putamen
$(100,107)$ | Pallidum
$(101,108)$ | Hippocampus
$(102,109)$ | Amygdala
$(103,110)$ | Accumbens