# Biarchetype analysis: simultaneous learning of observations and features
based on extremes

Aleix Alcacer
Department of Mathematics
Jaume I University
<EMAIL_ADDRESS>
&Irene Epifanio
Department of Mathematics
Jaume I University
<EMAIL_ADDRESS>
&Ximo Gual-Arnau
Department of Mathematics
Jaume I University

###### Abstract

A new exploratory technique called biarchetype analysis is defined. We extend
archetype analysis to find the archetypes of both observations and features
simultaneously. The idea of this new unsupervised machine learning tool is to
represent observations and features by instances of pure types (biarchetypes)
that can be easily interpreted as they are mixtures of observations and
features. Furthermore, the observations and features are expressed as mixtures
of the biarchetypes, which also helps understand the structure of the data. We
propose an algorithm to solve biarchetype analysis. We show that biarchetype
analysis offers advantages over biclustering, especially in terms of
interpretability. This is because byarchetypes are extreme instances as
opposed to the centroids returned by biclustering, which favors human
understanding. Biarchetype analysis is applied to several machine learning
problems to illustrate its usefulness. The source code and examples both in R
and Python are available at https://github.com/aleixalcacer/JA-BIAA.

_Keywords_ Archetype analysis, biclustering, prototype, unsupervised learning.

##  1 Introduction

Cluster analysis (CLA) is one of the most widely used tools in exploratory
data analysis. The idea of clustering is to make groups of observations in
such a way that each group contains similar observations that are different to
those of the rest of the groups. If the data consist of well-separated
clusters, appropriate clustering techniques can obtain, on the one hand, the
representative of each cluster (the mean or centroid of the cluster for the
popular $k$-means technique), and, on the other hand, the assignations of each
observation to one cluster, or a degree of belonging to each cluster for fuzzy
clustering techniques.

However, CLA is also used as a segmentation technique in the absence of well-
separated (clearly differentiated) clusters in data. Many times, data follow a
fan-spread pattern, i.e. features vary continuously across observations. The
centroids are located in the middle of the data cloud since data points have
to be covered in such a way that the distance between them and the assigned
centroid is minimized (see [Wu et al., 2016] about the relationship between
CLA and set partitioning). In those cases, where data can be viewed as a
superposition of various populations, it is of particular interest to use
Archetype Analysis (AA) for segmenting [Keller et al., 2019].

Instead of segmenting on the centroids, AA segments on the extremes. AA was
defined by [Cutler and Breiman, 1994]. The objective of AA is to represent the
observations by means of a convex combination of archetypes, which in turn are
convex combinations of observations. Archetypes or ‘pure types’ lie on the
boundary of the convex hull of the data and are therefore extreme profiles.
Being extreme instances rather than central instances makes human
understanding and interpretation of data easier [Thurau et al., 2012] since
human cognition prefers extreme opposites [Davis and Love, 2010]. An
illustrative example of this was analyzed in [Cabero and Epifanio, 2020],
where CLA and AA were compared and archetypes were much more informative and
understandable than centroids, because archetypes are further apart from each
other than centroids.

Biclustering is a data mining technique introduced by [Hartigan, 1972],
although it was popularized by [Cheng and Church, 2000], who applied it to
gene expression data analysis. In biclustering (also known as block
clustering, co-clustering, or two-mode clustering), rows (observations) and
columns (features) of a data matrix are simultaneously clustered. An excellent
overview of biclustering and fuzzy biclustering is found in [Ferraro et al.,
2021]. Biclustering is widely used in biological and medical applications
[Zhao et al., 2012], especially in gene expression data [Kerr et al., 2008,
Xie et al., 2019]. However, it is also applied in many other fields, such as
marketing [Dolnicar et al., 2012], psychology [Van Mechelen et al., 2004],
recommender systems [Forsati et al., 2013], sports [Kaiser, 2011, Shkedy et
al., 2016], website traffic [Koutsonikola and Vakali, 2009], and many other
pattern recognition applications, such as collaborative filtering, text
mining, multimedia data processing and retrieval, etc. [Zhao et al., 2012,
Henriques et al., 2015].

In recent years, there has been growing interest in AA. On the one hand, there
has been an increasing number of papers proposing efficient computational
methods to calculate AA [Mørup and Hansen, 2012, Chen et al., 2014, Bauckhage
et al., 2015, Mair et al., 2017], with applications in computer vision. On the
other hand, AA has been applied in other very diverse fields, such as,
climatology [Steinschneider and Lall, 2015, Su et al., 2017], ergonomics
[Epifanio et al., 2013, Alcacer et al., 2020], genetics [Shoval et al., 2012,
Thøgersen et al., 2013, Wang and Zhao, 2022], image processing [Zois et al.,
2017, Sun et al., 2017a, Sun et al., 2017b, Cabero and Epifanio, 2019],
machine learning problems [Mørup and Hansen, 2012, Ragozini and D’Esposito,
2015, Alcacer et al., 2021, Keller et al., 2021], market research [Porzio et
al., 2008], multi-document summarization [Canhasi and Kononenko, 2014],
nanotechnology [Fernandez and Barnard, 2015], neuroscience [Tsanousa et al.,
2015, Hinrich et al., 2016], sports [Eugster, 2012, Vinué and Epifanio, 2017,
Vinué and Epifanio, 2019] and sustainability [Thurau et al., 2012]. Finally,
other papers have proposed extensions and new methodologies derived from AA
with applications in a broad spectrum of fields: kernel AA [Mørup and Hansen,
2012], AA with missing data [Mørup and Hansen, 2012, Epifanio et al., 2020],
robust AA [Eugster and Leisch, 2011, Moliner and Epifanio, 2019], interval
archetypes [D’Esposito et al., 2012], archetypoid analysis (ADA) [Vinué et
al., 2015], functional AA [Epifanio, 2016b], data-driven prototype
identification [Ragozini et al., 2017], archetypal networks [Ragozini and
D’Esposito, 2015], probabilistic AA [Seth and Eugster, 2016b], AA for nominal
[Seth and Eugster, 2016a, Cabero and Epifanio, 2020] and ordinal observations
[Fernández et al., 2021], directional AA [Olsen et al., 2022], AA for shapes
[Epifanio et al., 2018], deep AA [Keller et al., 2021], and outlier detection
[Millán-Roures et al., 2018, Vinué and Epifanio, 2020, Cabero et al., 2021].
Nevertheless, no previous work has developed archetypal analysis for both rows
and columns simultaneously, which we refer to as biarchetype analysis (biAA),
co-archetype analysis or two-mode archetype analysis.

[Xie et al., 2019] reviews biclustering in biological and biomedical fields.
They point out the need to improve the interpretability of biclustering
results, and they also highlight that possible overlapping homogeneous
submatrices have to be identified. This clashes with the idea of CLA, whose
origin was to find separate (not overlapping) groups, but it is in the line
with the basis of AA. Moreover, biclustering of human gene expression data has
been used to identify phenotype–genotype associations in studies of common or
rare diseases. Note that archetypes themselves are phenotypes [Keller et al.,
2021]; in fact, archetypes have been used also to explain the evolutionary
development of biological systems [Tendler et al., 2015]. Therefore, putting
all this together, it seems that biAA could be a reasonable alternative to
biclustering in biology, as biAA could improve the interpretability of
results. Nevertheless, the fields of application of biAA are not just
restricted to biology; they would be the same as for biclustering, i.e. biAA
can be applied to many pattern recognition problems.

Our contributions consist of defining biAA for the first time, proposing a
computational method to calculate it, whose implementation is available in the
R package biaa https://github.com/aleixalcacer/biaa and the Python package
https://github.com/aleixalcacer/archetypes, showing how it works and the
advantages of using archetypes (extremes) rather than the centroids of
biclustering in an illustrative example, and finally, applying it to several
real data sets in different fields to demonstrate the usefulness of biAA in
various problems.

The outline of the paper is as follows: previous methodologies (CLA,
biclustering, fuzzy biclustering, AA) are reviewed in a common framework in
Sec. 2. In Sec. 3, biAA is defined and a computational procedure is proposed.
An illustrative example is used to exemplify biAA and compare it to
biclustering. In Sec. 4, our proposal is applied to three real data sets. Some
conclusions and ideas for future work are provided in Sec. 5.

##  2 Background

Matrix factorization is our common framework for describing the established
methods (as used in [Maurizio, 2001] for clustering) and our proposal. Let
${\bf X}_{n\times m}$ be a data matrix with $n$ observations and $m$
continuous features (they should be standardized in order to avoid problems if
they measure different dimensions). Let ${\bf\alpha}_{n\times k}$ and
${\bf\gamma}_{c\times m}$ be matrices with values in $[0,1]$. ${\bf\alpha}$ is
the membership matrix of the observations, while ${\bf\gamma}$ is the
membership matrix of the features. ${\bf Z}$ is the matrix of representative
instances that approximates ${\bf X}$. The objective is to minimize:
$\left\Arrowvert{\bf X}-{\bf\alpha Z\gamma}\right\Arrowvert^{2}$, with
different constraints, where $\left\Arrowvert.\right\Arrowvert$ stands for the
Frobenius norm.

###  2.1 Clustering

For clustering, ${\bf Z}$ is the matrix of centroids, which is computed by
${\bf Z}$ = $({\bf\alpha}^{\prime}{\bf\alpha})^{-1}{\bf\alpha}^{\prime}{\bf
X}{\bf\gamma}^{\prime}({\bf\gamma}{\bf\gamma}^{\prime})^{-1}$, where ′ denotes
transpose.

  * •

$k$-means clustering: The constraints are: $\sum_{g=1}^{k}\alpha_{ig}=1$ with
$\alpha_{ig}\in\\{0,1\\}$ for $i=1,\ldots,n$ and $\bf\gamma$ = ${\bf
I}_{m\times m}$ is the identity matrix of order $m$. The matrix ${\bf
Z}_{k\times m}$ =
$({\bf\alpha}^{\prime}{\bf\alpha})^{-1}{\bf\alpha}^{\prime}{\bf X}$ has the
centroids of each one of $k$ groups that partition the data set.

  * •

Fuzzy clustering: In soft clustering, each observation is assigned membership
to each group. The restrictions are: $\sum_{g=1}^{k}\alpha_{ig}=1$ with
$\alpha_{ig}\geq 0$ for $i=1,\ldots,n$ and $\bf\gamma$ = ${\bf I}_{m\times
m}$. Again, the matrix ${\bf Z}_{k\times m}$ has the centroids of each one of
$k$ groups.

  * •

Biclustering: This is also called double $k$-means with hard partitions by
[Maurizio, 2001], where algorithms to solve it are proposed. The constraints
are: $\sum_{g=1}^{k}\alpha_{ig}=1$ with $\alpha_{ig}\in\\{0,1\\}$ for
$i=1,\ldots,n$ and $\sum_{h=1}^{c}\gamma_{hj}=1$ with
$\gamma_{hj}\in\\{0,1\\}$ for $j=1,\ldots,m$. Now, the dimension of ${\bf Z}$
is ${k\times c}$, since there are $k$ groups for observations and $c$ groups
of variables.

  * •

Fuzzy biclustering: This is also called fuzzy double $k$-means by [Maurizio,
2001]. The constraints are now continuous: $\sum_{g=1}^{k}\alpha_{ig}=1$ with
$\alpha_{ig}\geq 0$ for $i=1,\ldots,n$ and $\sum_{h=1}^{c}\gamma_{hj}=1$ with
$\gamma_{hj}\geq 0$ for $j=1,\ldots,m$. [Ferraro et al., 2021] proposed
several algorithms for solving fuzzy double $k$-means with continuous data,
called FDkM and FDkMpf (Fuzzy Double $k$-Means with polynomial fuzzifiers),
whose Matlab implementations are available in [Ferraro et al., 2021].

Besides the previous framework, there are some proposals of model-based
biclustering. In that case, it is supposed that data are generated by a
mixture distribution, as in [Govaert and Nadif, 2003], referred to as BMM
(Block Mixture Model). Instead of memberships, it returns the final posterior
probabilities for rows and columns, in addition to the mean and variance of
each co-cluster. This is implemented in the R package blockcluster [Bhatia et
al., 2017].

###  2.2 Archetype analysis

In AA, ${\bf Z}_{k\times m}={\bf\beta}_{k\times n}{\bf X}_{n\times m}$, where
$\sum_{l=1}^{n}\beta_{gl}=1$ with $\beta_{gl}\geq 0$ for $g=1,\ldots,k$, i.e.
the archetypes are mixture of the data. The other restrictions are:
$\sum_{g=1}^{k}\alpha_{ig}=1$ with $\alpha_{ig}\geq 0$ for $i=1,\ldots,n$ and
$\bf\gamma$ = ${\bf I}_{m\times m}$. Therefore, the objective function to
minimize subject to the previous constraints is:

$\begin{split}RSS=\left\Arrowvert{\bf X}-{\bf\alpha
Z}\right\Arrowvert^{2}=\left\Arrowvert{\bf X}-{\bf\alpha\beta
X}\right\Arrowvert^{2}=\\\
\sum_{i=1}^{n}\sum_{j=1}^{m}\left(x_{ij}-\sum_{g=1}^{k}\alpha_{ig}z_{gj}\right)^{2}=\\\
\sum_{i=1}^{n}\sum_{j=1}^{m}\left(x_{ij}-\sum_{g=1}^{k}\alpha_{ig}\left(\sum_{l=1}^{n}\beta_{gl}x_{lj}\right)\right)^{2}.\end{split}$
(1)

The ${\bf\alpha}$ coefficients determine how much each archetype contributes
to the approximation of each observation, i.e. $\alpha_{ig}$ is the weight of
the archetype $g$ for the $i$-th observation. Archetypes are built as mixtures
of observations weighted by ${\bf\beta}$ coefficients.

If $k$ = 1, the archetype coincides with the mean, but with $k$ $>$ 1, the
archetypes are located on the boundary of the convex hull of the data [Cutler
and Breiman, 1994]. Archetypes are not necessarily nested, so different $k$s
may reveal distinct structures of the data. Therefore, as happens in other
unsupervised statistical learning procedures, the selection of the number $k$
of prototypes has to be determined. If we have prior knowledge of the
arrangement of the data, $k$ can be selected based on this. Otherwise, we can
use a simple but effective heuristic method, the elbow criterion, which has
been used elsewhere [Cutler and Breiman, 1994, Eugster and Leisch, 2009]. The
elbow criterion consists of displaying the RSS for different $k$ values and
choosing the value $k$ as the position where the elbow is located. This method
is also used in clustering.

[Cutler and Breiman, 1994] proposed an alternating minimizing algorithm to
find the matrices ${\bf\alpha}$ and ${\bf\beta}$ that minimizes RSS. This
consists of alternating between estimating the best ${\bf\alpha}$ for given
archetypes ${\bf Z}$, and the optimum archetypes ${\bf Z}$ for given
${\bf\alpha}$. In each phase, convex least squares problems have to be solved.
They used a penalized version of the non-negative least squares algorithm
[Lawson and Hanson, 1974].

##  3 Biarchetype analysis

###  3.1 Definition

In biAA, biarchetypes are ${\bf Z}_{k\times c}={\bf\beta}_{k\times n}{\bf
X}_{n\times m}{\bf\theta}_{m\times c}$, where $\sum_{l=1}^{n}\beta_{gl}=1$
with $\beta_{gl}\geq 0$ for $g=1,\ldots,k$ and $\sum_{r=1}^{m}\theta_{rh}=1$
with $\theta_{rh}\geq 0$ for $h=1,\ldots,c$, i.e. the archetypes are mixture
of the data points and variables. There are $k$ archetypes for rows and $c$
for columns. The other restrictions are: $\sum_{g=1}^{k}\alpha_{ig}=1$ with
$\alpha_{ig}\geq 0$ for $i=1,\ldots,n$ and $\sum_{h=1}^{c}\gamma_{hj}=1$ with
$\gamma_{hj}\geq 0$ for $j=1,\ldots,m$. Therefore, the objective function to
minimize subject to the previous constraints is:

$\begin{split}RSS=\left\Arrowvert{\bf X}-{\bf\alpha
Z\gamma}\right\Arrowvert^{2}=\left\Arrowvert{\bf X}-{\bf\alpha\beta
X\theta\gamma}\right\Arrowvert^{2}=\\\
\sum_{i=1}^{n}\sum_{j=1}^{m}\left(x_{ij}-\sum_{g=1}^{k}\sum_{h=1}^{c}\alpha_{ig}z_{gh}\gamma_{hj}\right)^{2}=\\\
\sum_{i=1}^{n}\sum_{j=1}^{m}\left(x_{ij}-\sum_{g=1}^{k}\sum_{h=1}^{c}\alpha_{ig}\left(\sum_{l=1}^{n}\sum_{r=1}^{m}\beta_{gl}x_{lr}\theta_{rh}\right)\gamma_{hj}\right)^{2}.\end{split}$
(2)

As before, the ${\bf\alpha}$ coefficients determine how much each archetype
contributes to the approximation of each observation, i.e. $\alpha_{ig}$ is
the weight of the archetype $g$ for the $i$-th observation. Analogously, the
${\bf\gamma}$ coefficients determine how much each archetype contributes to
the approximation of each variable, i.e. $\gamma_{hj}$ is the weight of the
archetype $h$ for the $j$-th variable. Biarchetypes are built as mixtures of
observations and variables weighted by ${\bf\beta}$ and ${\bf\theta}$
coefficients, respectively.

####  3.1.1 Location of biarchetypes

In this section, we state some results that help in understanding the behavior
of the row and column vectors of the biarchetype matrix ${\bf Z}_{k\times c}$.

Let ${\bf X}_{n\times m}$ be a data matrix with $n$ observations and $m$
continuous features. We denote by ${\bf x}_{i}^{d}$, $i=1,\dots,n$ the row
vectors of the matrix ${\bf X}_{n\times m}$ (in this case observations) and by
${\bf x}_{i}^{f}$, $i=1,\dots,m$ the column vectors of the matrix ${\bf
X}_{n\times m}$ (in this case features). This notation will be the same for
all the matrices used.

The problem is to find a matrix ${\bf Z}_{k\times c}$ with $1\leq k\leq n$,
$1\leq c\leq m$, which is expressed as ${\bf Z}_{k\times
c}={\bf\beta}_{k\times n}{\bf X}_{n\times m}{\bf\theta}_{m\times c}$ and the
matrices $\alpha_{n\times k}$, ${\bf Z}_{k\times c}$ and $\gamma_{c\times m}$
minimize

$RSS=\left\Arrowvert{\bf X}_{n\times m}-\alpha_{n\times k}{\bf Z}_{k\times
c}\gamma_{c\times m}\right\Arrowvert^{2}.$

Now we distinguish several cases depending on the values of $k$ and $c$.

  * •

Case I: $k=1$ and $c=1$. In this case, $\alpha_{n\times
1}=(1,\dots,1)^{\prime}$ and $\gamma_{1\times m}=(1,\dots,1)$; then, the real
value ${\bf Z}_{1\times 1}$ that minimizes RSS is the mean of all entries in
matrix ${\bf X}_{n\times m}$, that is,

${\bf Z}_{1\times 1}=\frac{1}{nm}\sum_{i=1}^{n}\sum_{j=1}^{m}x_{ij}.$

  * •

Case II: $k=n$ and $c=m$. In this case we consider $\alpha_{n\times n}={\bf
I}_{n\times n}$ and $\gamma_{m\times m}={\bf I}_{m\times m}$. Then, by
choosing ${\bf Z}_{n\times m}={\bf X}_{n\times m}$ we obtain RSS = 0.

  * •

Case III: $1\leq k<n$ and $c=m$. We consider $\gamma_{m\times m}={\bf
I}_{m\times m}$ and $\theta_{m\times m}={\bf I}_{m\times m}$, then, the
problem now consists of minimizing

$RSS=\left\Arrowvert{\bf X}_{n\times m}-\alpha_{n\times k}{\bf Z}_{k\times
m}\right\Arrowvert^{2},$

which is a typical problem in AA and the location of the archetypes is
explained in [Cutler and Breiman, 1994] and reviewed in Sec. 2.2.

  * •

Case IV: $k=n$ and $1\leq c<m$. This is the case of finding only the
archetypes of features, i.e. we consider ${\bf\alpha}={\bf I}_{n\times n}$ and
${\bf\beta}={\bf I}_{n\times n}$. As the Frobenius norm of a matrix is the
same as the Frobenius norm of its transpose, $\left\Arrowvert{\bf X}-{\bf
Z\gamma}\right\Arrowvert^{2}$ = $\left\Arrowvert{\bf
X}^{\prime}-{\bf\gamma^{\prime}Z^{\prime}}\right\Arrowvert^{2}$, features can
adopt the role of observations when ${\bf X}$ is transposed. Then, the same
reasoning as in the preceding case entails a problem of AA.

  * •

Case V: $k=1$ and $1<c<m$. In this case ${\bf\alpha}={\bf I}_{n\times n}$ and
the problem of minimizing

$RSS=\left\Arrowvert{\bf X}_{n\times m}-(1,\dots,1)^{\prime}\left({\bf
Z}_{1\times c}\gamma_{c\times m}\right)\right\Arrowvert^{2},$

is satisfied if the mean value $\frac{1}{n}\sum_{i=1}^{n}{\bf
x}_{i}^{d}=(\bar{x}_{1}^{d},\dots,\bar{x}_{m}^{d})=\left({\bf Z}_{1\times
c}\gamma_{c\times m}\right)$ where
$\bar{x}_{j}^{d}=\frac{1}{n}\sum_{i=1}^{n}x_{ij}$. Note that each real number
$\bar{x}_{i}^{d}$ belongs to the convex hull of the components of the vector
${\bf Z}_{1\times c}$.

  * •

Case VI: $c=1$ and $1<k<n$. In this case RSS is minimized if
$\frac{1}{m}\sum_{i=1}^{m}{\bf
x}_{i}^{f}=(\bar{x}_{1}^{f},\dots,\bar{x}_{n}^{f})=\left({\bf
Z}^{\prime}_{k\times 1}\alpha^{\prime}_{n\times k}\right)$ where
$\bar{x}_{i}^{f}=\frac{1}{m}\sum_{j=1}^{m}x_{ij}$. Note that each real number
$\bar{x}_{i}^{f}$ belongs to the convex hull of the components of the vector
${\bf Z}_{k\times 1}$.

  * •

Case VII: $1<k<n$ and $1<c<m$.

Let us call ${\bf V}_{n\times c}={\bf X}_{n\times m}{\bf\theta}_{m\times c}$;
then, each ${\bf v}_{j}^{f}$ $(j=1,\dots,c)$ belongs to the convex hull
$C_{X}^{f}$ of the data ${\bf x}_{i}^{f}$ $(i=1,\dots,m)$. Moreover, since
${\bf Z}_{k\times c}={\bf\beta}_{k\times n}{\bf V}_{n\times c}$, each vector
${\bf z}_{j}^{d}$ $(j=1,\dots,k)$ belongs to the convex hull $C_{V}^{d}$ of
the vectors ${\bf v}_{i}^{d}$ $(i=1,\dots,n)$.

###### Proposition 1

Having fixed the matrix ${\bf\theta}_{m\times c}$, there is a matrix of
biarchetypes ${\bf Z}_{k\times c}$ such that each row vector ${\bf z}_{j}^{d}$
$(j=1,\dots,k)$ belongs to the boundary of the convex hull $C_{V}^{d}$.

Now, let us call ${\bf Y}_{k\times m}=\beta_{k\times n}{\bf X}_{n\times m}$;
then, each ${\bf y}_{j}^{d}$ $(j=1,\dots,k)$ belongs to the convex hull
$C_{X}^{d}$ of the data ${\bf x}_{i}^{d}$ $(i=1,\dots,n)$. Moreover, since
${\bf Z}_{k\times c}={\bf Y}_{k\times m}{\bf\theta}_{m\times c}$, each vector
${\bf z}_{j}^{f}$ $(j=1,\dots,c)$ belongs to the convex hull $C_{Y}^{f}$ of
the vectors ${\bf y}_{i}^{f}$ $(i=1,\dots,m)$.

###### Proposition 2

Having fixed the matrix ${\bf\beta}_{k\times n}$, there is a matrix of
biarchetypes ${\bf Z}_{k\times c}$ such that each column vector ${\bf
z}_{j}^{f}$ $(j=1,\dots,c)$ belongs to the boundary of the convex hull
$C_{Y}^{f}$.

Proof of Propositions 1 and 2 is detailed in Appendix A.

###### Example 1

This toy example illustrates the location of the biarchetypes for different
values of $k$ and $c$, for the following matrix $\begin{pmatrix}1&2&3&4&5\\\
6&7&8&9&10\\\ 11&12&13&14&15\\\ 16&17&18&19&20\\\ 21&22&23&24&25\\\
\end{pmatrix}$.

For $k$ = 1 and $c$ = 1, $\bf z$ = 13 (mean of all entries of the matrix
according to Case I). For $k$ = 1 and $c$ =2, $\bf z$ =
$\begin{pmatrix}11&15\end{pmatrix}$; here
$(\bar{x}_{1}^{d},\dots,\bar{x}_{5}^{d})=(11,12,13,14,15)$ and each real
number $\bar{x}_{i}^{d}$ belongs to the convex hull of the components of $\bf
z$ according to Case V. For $k$ = 2 and $c$ = 1, $\bf z$ =
$\begin{pmatrix}3\\\ 23\end{pmatrix}$ and, according to Case VI, each real
number $\bar{x}_{i}^{f}$ belongs to the convex hull of the components of the
vector $\bf z$. For $k$ = 2 and $c$ = 2, $\bf z$ = $\begin{pmatrix}1&5\\\
21&25\end{pmatrix}$. For this last case, RSS = 0, and, according to Case VII,
the vectors ${\bf z}_{j}^{d}$ and ${\bf z}_{j}^{f}$, $j=1,2$ are located at
the boundary of convex sets.

###### Example 2

In this example, we will generate the data from a multivariate random
distribution.

The covariance matrix for both rows and columns is:

$\displaystyle{\Sigma=\begin{pmatrix}1&0.8&\dots&0.8&0.8&0&0&\dots&0&0\\\
0.8&1&\dots&0.8&0.8&0&0&\dots&0&0\\\
\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0.8&0.8&\dots&1&0.8&0&0&\cdots&0&0\\\ 0.8&0.8&\dots&0.8&1&0&0&\cdots&0&0\\\
0&0&\dots&0&0&1&0.8&\cdots&0.8&0.8\\\ 0&0&\dots&0&0&0.8&1&\cdots&0.8&0.8\\\
\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&\dots&0&0&0.8&0.8&\cdots&1&0.8\\\ 0&0&\dots&0&0&0.8&0.8&\cdots&0.8&1\\\
\end{pmatrix}}$

And the mean also for both rows and columns is:

$\displaystyle\mu=(0\cdots 0)$

The data generated can be seen in Figure 1(a). and Figure 1(b), the data are
represented using the coefficients $\alpha$ and $\gamma$ of biAA with $k$ = 2
and $c$ = 2. These coefficients are mapped to the $x$ and $y$ axes of the
figure and the $z$ axis represents the value of each observation. The
biarchetypes are represented in black and, as seen, they are at the extremes
of the data.

(a) Simulated data as a heatmap.

(b) Simulated data in the biarchetype space.

Figure 1: Representations of data and the biarchetypal space for Example 2,
i.e. representation of the coefficients $\alpha$ and $\gamma$.

####  3.1.2 Selecting the number of biarchetypes

As in the case of AA, if there is no information available a priori, we can
use the elbow criterion, but in this case, we look for the elbow of a surface
instead of a curve. In biAA, we run biAA for different values of $k$ and $c$
and display their RSS values in a 3D plot. We select the point ($k$,$c$) where
the surface ”flattens”, i.e. $(k,c)$ is the point at which the RSS of the
following points $(k+1,c)$, $(k,c+1)$, and $(k+1,c+1)$ stops decreasing
drastically with respect to the RSS of the point $(k,c)$.

###  3.2 Algorithm

The following iterative method is proposed to solve biAA. It is based on
alternating minimization as in the AA algorithm by [Cutler and Breiman, 1994].

  1. 1.

Data preparation: To randomly initialize the matrices $\bf\alpha$,
$\bf\gamma$, $\bf\beta$ and $\bf\theta$, fulfilling the constraints in Sect.
3.1.

  2. 2.

Repeat until RSS is sufficiently small or the number of maximum iterations is
reached:

    1. (a)

Find the best $\bf\alpha$ (fixed $\bf\gamma$): solve $n$ convex least squares
problems using ${\bf X}^{\prime}=({\bf
Z\gamma})^{\prime}{\bf\alpha}^{\prime}$.

    2. (b)

Find the best $\bf\gamma$ (fixed $\bf\alpha$): solve $m$ convex least squares
problems using $\bf X$ = $({\bf\alpha Z}){\bf\gamma}$.

    3. (c)

Recalculate the biarchetypes, where + stands for the Moore-Penrose
pseudoinverse: ${\bf Z}={\bf\alpha}^{+}{\bf X\gamma}^{+}$.

    4. (d)

Find the best $\bf\beta$ (fixed $\bf\theta$): solve $k$ convex least squares
problems using ${\bf Z}^{\prime}=({\bf X\theta})^{\prime}{\bf\beta}^{\prime}$.

    5. (e)

Find the best $\bf\theta$ (fixed $\bf\beta$): solve $c$ convex least squares
problems using ${\bf Z}=({\bf\beta X}){\bf\theta}$.

    6. (f)

Recalculate the biarchetypes: ${\bf Z}={\bf\beta X\theta}$.

    7. (g)

Calculate the new RSS.

Note that the convex least squares problems can be solved as proposed by
[Cutler and Breiman, 1994], i.e. using a penalized least squares problem
[Lawson and Hanson, 1974]. The idea is, given a least squares problem ${\bf
A}_{n\times k}{\bf X}_{k\times m}={\bf B}_{n\times m}$, to add a row with
constant elements $C$ to $\bf A$ and $\bf B$, in order to obtain a new problem
${\bf A}_{(n+1)\times k}{\bf X}_{k\times m}={\bf B}_{(n+1)\times m}$, in such
a way that RSS would be:

$RSS=\sum_{i=1}^{n+1}\sum_{j=1}^{m}\left(b_{ij}-\sum_{h=1}^{k}a_{ih}x_{hj}\right)^{2}=$
(3)
$\displaystyle\text{$\sum_{j=1}^{m}\left(\sum_{i=1}^{n}\left(b_{ij}-\sum_{h=1}^{k}a_{ih}x_{hj}\right)^{2}+\left(b_{n+1,j}-\sum_{h=1}^{k}a_{n+1,h}x_{hj}\right)^{2}\right)=$}$
$\displaystyle\text{$\sum_{j=1}^{m}\left(\sum_{i=1}^{n}\left(b_{ij}-\sum_{h=1}^{k}a_{ih}x_{hj}\right)^{2}+\sum_{j=1}^{m}\left(C-\sum_{h=1}^{k}Cx_{hj}\right)^{2}\right)=$}$
$\displaystyle\text{$\sum_{j=1}^{m}\left(\sum_{i=1}^{n}\left(b_{ij}-\sum_{h=1}^{k}a_{ih}x_{hj}\right)^{2}+\sum_{j=1}^{m}C^{2}\left(1-\sum_{h=1}^{k}x_{hj}\right)^{2}\right)$}.$

Therefore, if value $C$ is high, the term
$C^{2}\left(1-\sum_{h=1}^{k}x_{hj}\right)^{2}$ forces the convexity of the
elements of $\bf X$ in eq. 3.

###  3.3 Illustrative example and comparison with biclustering

The following example illustrates the use of biAA and its advantages in
comparison with biclustering, especially when working with non-clustered data.
We consider the data of 45 students from Universitat Jaume I, who reported the
number of hours per week spent working on a subject at home over 17 weeks. The
complete description of the data can be found in [Epifanio, 2016a]. Missing
data are imputed by $mice$ [van Buuren and Groothuis-Oudshoorn, 2011]. The
data range from 0 to 10, with mean 4.94 and standard deviation 2.46. We apply
biAA with $k$ = 4 and $c$ = 3, since the elbow is found at those values (see
Fig. 2).

Figure 2: RSS for $k$ from 2 to 6 and $c$ from 2 to 5.

The $\gamma$ coefficients for biAA, BMM, and FDkMpf are shown in Table 1 (we
sort them in order to make the comments easier). FDkM was also applied, but
the solution is not valid since the same prototypes are obtained for all the
groups.

Table 1: The $\gamma$ coefficients for biAA, BMM, and FDkMpf of the illustrative example. 1 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 0
---|---|---|---|---|---|---|---|---
0.71 | 0.12 | 0.17 | 1 | 0 | 0 | 1 | 0 | 0
0.57 | 0.42 | 0.01 | 1 | 0 | 0 | 0.85 | 0.08 | 0.08
0 | 0.80 | 0.20 | 0 | 0.01 | 0.99 | 0 | 0.5 | 0.5
0.13 | 0.76 | 0.11 | 0 | 1 | 0 | 0 | 0.5 | 0.5
0.14 | 0.86 | 0 | 0 | 1 | 0 | 0.29 | 0.35 | 0.35
0 | 0.80 | 0.20 | 0 | 1 | 0 | 0 | 0.5 | 0.5
0.01 | 0.99 | 0 | 0 | 0 | 1 | 0.2 | 0.4 | 0.4
0 | 0.90 | 0.11 | 0 | 1 | 0 | 0 | 0.5 | 0.5
0 | 0.89 | 0.11 | 0 | 1 | 0 | 0 | 0.5 | 0.5
0 | 0.42 | 0.58 | 0 | 1 | 0 | 0 | 0.5 | 0.5
0.05 | 0.45 | 0.5 | 0 | 0.96 | 0.04 | 0 | 0.5 | 0.5
0 | 0.62 | 0.38 | 0 | 1 | 0 | 0 | 0.5 | 0.5
0.05 | 0 | 0.95 | 0 | 0 | 1 | 0 | 0.5 | 0.5
0 | 0.04 | 0.96 | 0 | 0 | 1 | 0 | 0.5 | 0.5
0 | 0.10 | 0.91 | 0 | 0 | 1 | 0 | 0.5 | 0.5
0.13 | 0 | 0.87 | 0 | 0 | 1 | 0 | 0.5 | 0.5

(From now on, we will use archetype or archetypal instead of biarchetype or
biarchetypal to simplify the language).

The feature similar to the first archetypal variable corresponds to the first
week. The second and third weeks are also similar, but with a temporal
gradation (0.71 and 0.57). Week 8, an intermediate week of the semester, is
similar to the second archetypal variable. Other intermediate weeks (4, 5, 6,
7, 9, and 10) are also similar. Week 15, a week at the end of the semester, is
similar to the third archetypal variable, as well as weeks 14, 16, and 17.
Weeks 11, 12, and 13 are explained as mixtures (nearly 50% -50%) of the second
and third archetypal variables. Note that the third week was also explained as
a mixture close to 50% -50% of the first and second archetypal variables. In
summary, the archetypal variables correspond to the profile of the beginning,
middle and end of the semester, respectively. The weeks in the transitions
between these temporal points are reflected as mixtures.

As regards the prototypical variables for biclustering methodologies, the
first three weeks have probabilities of 1 (or nearly 1 for FDkMpf) for the
first prototypical variable of BMM. Unlike biAA where gradation was found, the
probabilities (memberships) are nearly crisp classifications for BMM, not only
for the first prototypical variable, but for the rest as well. The
intermediate weeks 5, 6, 7, 9, 10, 11, 12, and 13 have probabilities of 1 for
the second prototypical variable, while the final weeks (14, 15, 16, and 17),
and the intermediate weeks 4 and 8 have probabilities of nearly 1 or 1 for the
third prototypical variable.

Note the difference with biAA. On the one hand, in BMM there is a lack of
gradation over time in the memberships (no mixture is found, but the
memberships are extremely high, nearly all ones), as if changes between
adjoining weeks were radical (as breaking jumps) rather than smooth.
Therefore, the information provided by biAA is richer. On the other hand,
there are two intermediate weeks (4 and 8) belonging to the third prototypical
variable corresponding to the end of semester weeks, which is not very
coherent. Therefore, the information provided by biAA is more reasonable.
Finally, the second and third prototypical variables are identical for FDkMpf,
with 50%-50% or close degrees of membership for weeks 4 to 17. Therefore, the
information returned by FDkMpf is poorer than that of BMM and biAA.

Table 2 displays the representative points ${\bf Z}$, archetypes or centroids
for biAA and biclustering, respectively (we sort them in order to make the
comments easier). The first archetype describes a student who works very few
hours per week throughout the semester. The second archetype represents a
student who studies very few hours throughout the semester, except at the end
of the semester, when they work for 9 hours per week. The third archetype
describes a student who works very few hours at the beginning of the semester
(1h per week), many hours during the semester (10h per week), and intermediate
hours (4h per week) at the end of the semester. The fourth archetype
represents a student who studies many hours throughout the whole semester.

Table 2: The archetypes and centroids for biAA, BMM, and FDkMpf of the toy example. 1.54 | 2.36 | 0.36 | 1.67 | 2.07 | 3.49 | 2.03 | 3.28 | 3.28
---|---|---|---|---|---|---|---|---
1 | 2 | 9 | 2.32 | 4.15 | 5.02 | 2.74 | 4.89 | 4.89
1 | 10 | 4 | 3.79 | 5.70 | 6.05 | 3.42 | 5.97 | 5.97
8 | 6.32 | 9.36 | 5.01 | 7.69 | 6.60 | 5.55 | 6.97 | 6.97

For BMM, the prototypes are not as pure as the archetypes. For example, there
is no great difference between centroids 2 and 3: centroid 3 studies only one
or one and a half hours more than centroid 2 per week. The centroids are not
as intuitively interpretable as archetypes. Centroid 1 corresponds to a
student who studies 2 or 3 hours throughout the semester; centroid 2 studies
2h per week at the beginning and 4 or 5h per week for the rest of the
semester; centroid 3 works 4 hours at the beginning and 6 hours per week for
the rest of the semester; while centroid 4 studies 5h per week at the
beginning of the semester and 7 or 8 hours throughout the rest of the
semester. It seems that centroids are limited to following a gradation
according to the total number of hours studied throughout the semester rather
than by differences in behavior throughout the semester. For FDkMpf, the
comments are similar, but in addition there is no difference between the
intermediate and final weeks. For example, the profiles of students 32 and 33,
who are similar to archetype 2 (with $\alpha$s of 0.84 and 0.88,
respectively), would not be reflected by the centroids of BMM or FDkMpf. They
belong to cluster 2 of BMM, with probabilities of 0.97 and 1, respectively.
But this does not say anything about how far (or in which direction) from
centroid 2 those students are. This happens because the goal of clustering is
to assign the data to groups, not to explain the structure of the data more
qualitatively.

##  4 Results and discussion

Like biclustering, biAA can be applied to a wide range of fields. In this
section, we will apply it to biology, document analysis and community
detection.

The code and data sets for reproducing the results including those in Sec. 3.3
are available at https://github.com/aleixalcacer/JA-BIAA.

###  4.1 Gene expression data

To show how biAA can be applied, we examine data from gene expression of
cutaneous melanoma used in [Bittner et al., 2000, Rocci and Vichi, 2008].
Instead of meticulously re-analyzing this data set, we use it to highlight the
salient features of biAA.

The aim of this study was to test the idea that molecular profiles generated
by cDNA microarrays could be used to differentiate between several subtypes of
cutaneous melanoma, a kind of skin cancer. mRNA was collected from the 31
cutaneous melanoma samples, and Cy5-labeled cDNA was created. All samples were
examined with the same reference probe, identified as Cy3. For each sample,
Cy5 and Cy3-labeled cDNA were combined and hybridized to a different melanoma
microarray. Red and green lasers were used to scan the hybridization array,
and the resulting image was then analyzed.

The same pre-processing was carried out as in [Bittner et al., 2000, Rocci and
Vichi, 2008]. Only 3613 cDNAs of the 8150 observations were classified as well
measured. Cy5/Cy3 expression ratios were computed for the accurately measured
genes. Ratios that were more than or equal to 50 and less than or equal to
0.02 were reduced to 50 and 0.02, respectively. A logarithmic scale was
applied to the derived ratios (base 2). The log ratios were adjusted so that
the median log-ratio for each experiment was equal to zero by subtracting the
median log-ratio within an experiment from all log-ratios for that experiment.
Since a single reference probe was utilized in all experiments, there was no
standardization between trials.

Using one minus the Pearson correlation coefficient of log-ratios as a measure
of dissimilarity between two experiments, [Bittner et al., 2000] applied the
average linkage hierarchical clustering on the 31 cutaneous melanoma samples
and obtained two clusters of 12 and 19 samples.

Regarding [Rocci and Vichi, 2008], they used double $k$-means to analyze the
same data set. In this case, the columns were centered and scaled to unit
variance, finding that the separation between two columns is proportional to
one minus the Pearson correlation coefficient. In particular, their analysis
indicates that samples 4 and 7 are members of the ‘19-samples’ group obtained
by [Bittner et al., 2000], i.e. the main cluster group. The membership of
these two samples in [Rocci and Vichi, 2008] differs from that obtained by
[Bittner et al., 2000].

In our case, we applied biarchetype analysis to the same data set as [Rocci
and Vichi, 2008], extracting three archetypes for the genes (rows) and two
archetypes for the melanoma samples (columns).

As can be seen in Fig. 3, regarding the melanoma samples, if we cluster the
data using the location of the maximum archetypal coefficient $\gamma$ (i.e.
the archetype that is most similar to the sample), we obtain two clusters of 8
and 23 samples. Regarding the archetypes of the genes, the first two
discriminate the two groups quite well, while the third archetype (or group of
genes) is not expressed for any melanoma group.

Figure 3: Representation of the most similar observations to each archetype.
The color represents the expression ratio of each gene for each sample.

If we compare our results to those obtained by [Bittner et al., 2000], four
samples are classified in different clusters. Samples 1, 4, 7 and 8 belong to
the main cluster group with biAA (see Fig. 4), unlike results provided by
[Bittner et al., 2000]. If we compare our results to those obtained by [Rocci
and Vichi, 2008], two samples are classified in different clusters, samples 1
and 8. Classification group of samples 4 and 7 is shared with biAA and results
by [Rocci and Vichi, 2008]. According to the $\gamma$ coefficients of biAA,
samples 4 and 7 are a nearly equal mixture between both archetypes, with the
values of the coefficient being 0.4 and 0.6 corresponding to the first and
second archetype for sample 4, and 0.45 and 0.55 corresponding to the first
and second archetype for sample 7. Therefore, samples 4 and 7 could be in the
border between both groups, which could explain the difference in
classification by different methods. However, the $\gamma$ coefficients for
sample 8 and 1 are 0.75 (0.25) and 1 (0) for the second (first) archetype,
respectively. In other words, sample 1 (and to lesser extent sample 8) should
definitely be in the main group according with biAA.

(a) Representation of the observations ($\alpha$ matrix) in their archetypal
space. Each corner represents an archetype.

(b) Representation of the variables ($\gamma$ matrix) in their archetypal
space. Each corner represents an archetype.

Figure 4: Representations of the archetypal spaces for Gene expression data.

###  4.2 Text documents

Another common use for biclustering is for clustering documents and words. In
this case, we have applied biAA to a subset of the 20 Newsgroups collection,
set up by [Lang, 1995]. Specifically, we have analyzed three topics:
_rec.autos_ , _rec.sport.hockey_ and _talk.politics.guns_.

For each document, the number of times each word is repeated in the document
has been stored in a count matrix, where each row represents a document, each
column a word, and the values indicate how many times each word is repeated in
each document.

In addition, a Tfidf transformer [Lang, 1995] was used to convert a count
matrix into a normalized tf or tf-idf representation. Tf stands for term
frequency, and tf-idf stands for term frequency multiplied by inverse document
frequency. This is a standard term weighting method used in information
retrieval, and it is also effective for classifying documents.

We have applied biarchetype analysis to this normalized matrix, obtaining
three archetypal profiles for documents and three archetypal profiles for
words.

In Fig. 5(a), it is clear that the three archetypes discriminate the three
groups of documents perfectly.

(a) Representation of the documents in their archetypal space ($\alpha$
values). The color represents the category of the document.

(b) Representation of the most similar words to each archetype (filtered using
a threshold over the coefficient matrix). The weight of each edge is the
cosine similarity between the two words in their archetypal space. This plot
was created using the _networkx_ Python package.

Figure 5: Results for the text documents example

Regarding the words, in Fig. 5(b) the words are represented as a graph. The
weight of each edge represents the similarity of the words in terms of the
archetypes (the gamma coefficients). The graph weights (i.e. the gamma
coefficients) split the words into three groups, where the words within each
group are related to one of the selected topics.

###  4.3 Community detection

Finally, we have also applied biAA to detect communities within the company
Enron. For that, we have studied the data set described in [Klimt and Yang,
2004], which contains a collection of emails between the company’s employees.

We have created an adjacency matrix between employees, containing 1 if one
employee has emailed another or 0 if the first one has never sent an email to
the second one.

After applying biAA with $k=6$ and $c=6$ to this adjacency matrix, we obtained
the results in Fig. 6(a). In the ‘senders’ part, the group $Z_{3}$ could be
omitted (represents employees who haven’t sent emails to a specific group).
The same occurs with the group $Z_{6}$ in the ‘recipients’ part of the
adjacency matrix.

If we analyze the employees from the point of view of who they send emails to,
we obtain the results shown in Fig. 6(b), in which we have removed
$Z_{3}$-like employees as they do not exchange emails with a specific group.
As can be seen, the dark blue cluster represents the Legal department, the
green one, the ETS department; the olive green one, the Gas/Energy department;
the pink one, the West Power departments; and the blue one, Enron’s top
management.

(a) The adjacency matrix ordered according to the archetypes obtained with
BiAA.

(b) Representation of employees from the point of view of who they send emails
to. The weight of each node is computed as in Fig. 5(b). The size of each
employee is proportional to how similar it is to its closest archetype and the
color of each one is determined by the closest archetype.

Figure 6: Results for the community detection example

##  5 Conclusion

In this work, we have proposed a new unsupervised machine learning technique:
biarchetype analysis. We have seen its usefulness in several problems in
different fields. We have compared the results of biAA and biclustering in an
illustrative example, showing not only the greater interpretability provided
by biAA, but also the greater coherence of the results.

In this work, biAA has been defined for continuous data. In future work, it
could be extended to other kinds of data, such as functional data, to which AA
was also extended [Epifanio, 2016b]. Note that biclustering analysis of time
series is used in many fields such as neuroscience [Castanho et al., 2022] and
engineering [Silva et al., 2022]; therefore, biAA could also be used for the
same problems. Biarchetypoid analysis could also be introduced in the same way
that archetypoid analysis was defined [Vinué et al., 2015], where biarchetypes
are not determined by mixtures of observations and features, but by concrete
elements of the data set. Just as archetype analysis is sensitive to outliers,
biAA is too. Robust biAA could be defined in the same way as robust AA was
[Moliner and Epifanio, 2019]. Likewise, biAA for missing data could be defined
as it was for AA [Epifanio et al., 2020], and it could be used in recommender
systems to find profiles of users and products, for instance. Another line of
future work would be to apply biAA to different fields where biclustering
analysis is applied, and to study more computational methods to calculate
biAA, especially for big data. Finally, biAA could also be easily extended to
high dimensions in a similar way to the decomposition proposed in [Tucker,
1966].

## Acknowledgments

This research was partially supported by the Spanish Ministry of Universities
(FPU grant FPU20/01825), Spanish Ministry of Science and Innovation
(PID2020-118763GA-I00 and PID2020-115930GA-I00) and UJI-B2020-22 from
Universitat Jaume I, Spain.

## References

  * [Alcacer et al., 2020] Alcacer, A., Epifanio, I., Ibáñez, M. V., Simó, A., and Ballester, A. (2020).  A data-driven classification of 3D foot types by archetypal shapes based on landmarks.  PLOS ONE, 15(1):e0228016.
  * [Alcacer et al., 2021] Alcacer, A., Epifanio, I., Valero, J., and Ballester, A. (2021).  Combining classification and user-based collaborative filtering for matching footwear size.  Mathematics, 9(7).
  * [Bauckhage et al., 2015] Bauckhage, C., Kersting, K., Hoppe, F., and Thurau, C. (2015).  Archetypal analysis as an autoencoder.  In Workshop New Challenges in Neural Computation, pages 8–15.
  * [Bhatia et al., 2017] Bhatia, P. S., Iovleff, S., and Govaert, G. (2017).  blockcluster: An R package for model-based co-clustering.  Journal of Statistical Software, 76(9):1–24.
  * [Bittner et al., 2000] Bittner, M., Meltzer, P., Chen, Y., Jiang, Y., Seftor, E., Hendrix, M., Radmacher, M., Simon, R., Yakhini, Z., Ben-Dor, A., et al. (2000).  Molecular classification of cutaneous malignant melanoma by gene expression profiling.  Nature, 406(6795):536–540.
  * [Cabero and Epifanio, 2019] Cabero, I. and Epifanio, I. (2019).  Archetypal analysis: an alternative to clustering for unsupervised texture segmentation.  Image Analysis & Stereology, 38:151–160.
  * [Cabero and Epifanio, 2020] Cabero, I. and Epifanio, I. (2020).  Finding archetypal patterns for binary questionnaires.  SORT, 44(1):39–66.
  * [Cabero et al., 2021] Cabero, I., Epifanio, I., Piérola, A., and Ballester, A. (2021).  Archetype analysis: A new subspace outlier detection approach.  Knowledge-Based Systems, 217:106830.
  * [Canhasi and Kononenko, 2014] Canhasi, E. and Kononenko, I. (2014).  Weighted archetypal analysis of the multi-element graph for query-focused multi-document summarization.  Expert Systems with Applications, 41(2):535 – 543.
  * [Castanho et al., 2022] Castanho, E. N., Aidos, H., and Madeira, S. C. (2022).  Biclustering fMRI time series: a comparative study.  BMC bioinformatics, 23(1):1–30.
  * [Chen et al., 2014] Chen, Y., Mairal, J., and Harchaoui, Z. (2014).  Fast and Robust Archetypal Analysis for Representation Learning.  In CVPR 2014 - IEEE Conference on Computer Vision & Pattern Recognition, pages 1478–1485.
  * [Cheng and Church, 2000] Cheng, Y. and Church, G. M. (2000).  Biclustering of expression data.  In Proceedings of the Eighth International Conference on Intelligent Systems for Molecular Biology, volume 8, pages 93–103.
  * [Cutler and Breiman, 1994] Cutler, A. and Breiman, L. (1994).  Archetypal Analysis.  Technometrics, 36(4):338–347.
  * [Davis and Love, 2010] Davis, T. and Love, B. (2010).  Memory for category information is idealized through contrast with competing options.  Psychological Science, 21(2):234–242.
  * [D’Esposito et al., 2012] D’Esposito, M. R., Palumbo, F., and Ragozini, G. (2012).  Interval Archetypes: A New Tool for Interval Data Analysis.  Statistical Analysis and Data Mining, 5(4):322–335.
  * [Dolnicar et al., 2012] Dolnicar, S., Kaiser, S., Lazarevski, K., and Leisch, F. (2012).  Biclustering: Overcoming data dimensionality problems in market segmentation.  Journal of Travel Research, 51(1):41–49.
  * [Epifanio, 2016a] Epifanio, I. (2016a).  Cargas de trabajo no presencial ECTS arquetípicas del estudiantado:¿cómo se reparten el trabajo semanalmente?  In Actas del Congreso Virtual: Avances en Tecnologías, Innovación y Desafíos de la Educación Superior ATIDES 2016, pages 367–376.
  * [Epifanio, 2016b] Epifanio, I. (2016b).  Functional archetype and archetypoid analysis.  Computational Statistics & Data Analysis, 104:24 – 34.
  * [Epifanio et al., 2020] Epifanio, I., Ibáñez, M. V., and Simó, A. (2020).  Archetypal analysis with missing data: see all samples by looking at a few based on extreme profiles.  The American Statistician, 74(2):169–183.
  * [Epifanio et al., 2018] Epifanio, I., Ibáñez, M. V., and Simó, A. (2018).  Archetypal shapes based on landmarks and extension to handle missing data.  Advances in Data Analysis and Classification, 12(3):705–735.
  * [Epifanio et al., 2013] Epifanio, I., Vinué, G., and Alemany, S. (2013).  Archetypal analysis: contributions for estimating boundary cases in multivariate accommodation problem.  Computers & Industrial Engineering, 64(3):757–765.
  * [Eugster and Leisch, 2009] Eugster, M. J. and Leisch, F. (2009).  From Spider-Man to Hero - Archetypal Analysis in R.  Journal of Statistical Software, 30(8):1–23.
  * [Eugster, 2012] Eugster, M. J. A. (2012).  Performance profiles based on archetypal athletes.  International Journal of Performance Analysis in Sport, 12(1):166–187.
  * [Eugster and Leisch, 2011] Eugster, M. J. A. and Leisch, F. (2011).  Weighted and robust archetypal analysis.  Computational Statistics & Data Analysis, 55(3):1215–1225.
  * [Fernandez and Barnard, 2015] Fernandez, M. and Barnard, A. S. (2015).  Identification of nanoparticle prototypes and archetypes.  ACS Nano, 9(12):11980–11992.
  * [Fernández et al., 2021] Fernández, D., Epifanio, I., and McMillan, L. F. (2021).  Archetypal analysis for ordinal data.  Information Sciences, 579:281–292.
  * [Ferraro et al., 2021] Ferraro, M. B., Giordani, P., and Vichi, M. (2021).  A class of two-mode clustering algorithms in a fuzzy setting.  Econometrics and Statistics, 18:63–78.
  * [Forsati et al., 2013] Forsati, R., Doustdar, H. M., Shamsfard, M., Keikha, A., and Meybodi, M. R. (2013).  A fuzzy co-clustering approach for hybrid recommender systems.  International Journal of Hybrid Intelligent Systems, 10(2):71–81.
  * [Govaert and Nadif, 2003] Govaert, G. and Nadif, M. (2003).  Clustering with block mixture models.  Pattern Recognition, 36(2):463–473.
  * [Hartigan, 1972] Hartigan, J. A. (1972).  Direct clustering of a data matrix.  Journal of the American Statistical Association, 67(337):123–129.
  * [Henriques et al., 2015] Henriques, R., Antunes, C., and Madeira, S. C. (2015).  A structured view on pattern mining-based biclustering.  Pattern Recognition, 48(12):3941–3958.
  * [Hinrich et al., 2016] Hinrich, J. L., Bardenfleth, S. E., Roge, R. E., Churchill, N. W., Madsen, K. H., and Mørup, M. (2016).  Archetypal analysis for modeling multisubject fMRI data.  IEEE Journal on Selected Topics in Signal Processing, 10(7):1160–1171.
  * [Kaiser, 2011] Kaiser, S. (2011).  Biclustering: methods, software and application.  PhD thesis, Ludwig-Maximilians-Universität München.
  * [Keller et al., 2021] Keller, S. M., Samarin, M., Arend Torres, F., Wieser, M., and Roth, V. (2021).  Learning extremal representations with deep archetypal analysis.  International journal of computer vision, 129(4):805–820.
  * [Keller et al., 2019] Keller, S. M., Samarin, M., Wieser, M., and Roth, V. (2019).  Deep archetypal analysis.  In Fink, G. A., Frintrop, S., and Jiang, X., editors, Pattern Recognition, pages 171–185, Cham. Springer International Publishing.
  * [Kerr et al., 2008] Kerr, G., Ruskin, H. J., Crane, M., and Doolan, P. (2008).  Techniques for clustering gene expression data.  Computers in Biology and Medicine, 38(3):283–293.
  * [Klimt and Yang, 2004] Klimt, B. and Yang, Y. (2004).  The Enron corpus: A new dataset for email classification research.  In Proceedings of the 15th European Conference on Machine Learning, page 217–226, Berlin, Heidelberg. Springer-Verlag.
  * [Koutsonikola and Vakali, 2009] Koutsonikola, V. A. and Vakali, A. (2009).  A fuzzy bi-clustering approach to correlate web users and pages.  IJ Knowledge and Web Intelligence, 1(1/2):3–23.
  * [Lang, 1995] Lang, K. (1995).  Newsweeder: Learning to filter netnews.  In Proceedings of the Twelfth International Conference on Machine Learning, pages 331–339.
  * [Lawson and Hanson, 1974] Lawson, C. L. and Hanson, R. J. (1974).  Solving Least Squares Problems.  Prentice Hall, Englewood Cliffs.
  * [Mair et al., 2017] Mair, S., Boubekki, A., and Brefeld, U. (2017).  Frame-based data factorizations.  In International Conference on Machine Learning, pages 2305–2313.
  * [Maurizio, 2001] Maurizio, V. (2001).  Double k-means clustering for simultaneous classification of objects and variables.  In Borra, S., Rocci, R., Vichi, M., and Schader, M., editors, Advances in Classification and Data Analysis, pages 43–52, Berlin, Heidelberg. Springer.
  * [Millán-Roures et al., 2018] Millán-Roures, L., Epifanio, I., and Martínez, V. (2018).  Detection of anomalies in water networks by functional data analysis.  Mathematical Problems in Engineering, 2018(Article ID 5129735):13.
  * [Moliner and Epifanio, 2019] Moliner, J. and Epifanio, I. (2019).  Robust multivariate and functional archetypal analysis with application to financial time series analysis.  Physica A: Statistical Mechanics and its Applications, 519:195 – 208.
  * [Mørup and Hansen, 2012] Mørup, M. and Hansen, L. K. (2012).  Archetypal analysis for machine learning and data mining.  Neurocomputing, 80:54–63.
  * [Olsen et al., 2022] Olsen, A. S., Høegh, R. M. T., Hinrich, J. L., Madsen, K. H., and Mørup, M. (2022).  Combining electro- and magnetoencephalography data using directional archetypal analysis.  Frontiers in Neuroscience, 16.
  * [Porzio et al., 2008] Porzio, G. C., Ragozini, G., and Vistocco, D. (2008).  On the use of archetypes as benchmarks.  Applied Stochastic Models in Business and Industry, 24:419–437.
  * [Ragozini and D’Esposito, 2015] Ragozini, G. and D’Esposito, M. R. (2015).  Archetypal networks.  In Proceedings of the 2015 IEEE/ACM International Conference on Advances in Social Networks Analysis and Mining 2015, pages 807–814, New York, NY, USA. ACM.
  * [Ragozini et al., 2017] Ragozini, G., Palumbo, F., and D’Esposito, M. R. (2017).  Archetypal analysis for data-driven prototype identification.  Statistical Analysis and Data Mining: The ASA Data Science Journal, 10(1):6–20.
  * [Rocci and Vichi, 2008] Rocci, R. and Vichi, M. (2008).  Two-mode multi-partitioning.  Computational Statistics and Data Analysis, 52:1984–2003.
  * [Seth and Eugster, 2016a] Seth, S. and Eugster, M. J. A. (2016a).  Archetypal analysis for nominal observations.  IEEE Trans. Pattern Anal. Mach. Intell., 38(5):849–861.
  * [Seth and Eugster, 2016b] Seth, S. and Eugster, M. J. A. (2016b).  Probabilistic archetypal analysis.  Machine Learning, 102(1):85–113.
  * [Shkedy et al., 2016] Shkedy, Z., Sengupta, R., and Perualila, N. J. (2016).  Identification of local patterns in the NBA performance indicators.  In Applied Biclustering Methods for Big and High-Dimensional Data Using R, pages 323–344. Chapman and Hall/CRC.
  * [Shoval et al., 2012] Shoval, O., Sheftel, H., Shinar, G., Hart, Y., Ramote, O., Mayo, A., Dekel, E., Kavanagh, K., and Alon, U. (2012).  Evolutionary trade-offs, Pareto optimality, and the geometry of phenotype space.  Science, 336(6085):1157–1160.
  * [Silva et al., 2022] Silva, M. G., Madeira, S. C., and Henriques, R. (2022).  Water consumption pattern analysis using biclustering: When, why and how.  Water, 14(12).
  * [Steinschneider and Lall, 2015] Steinschneider, S. and Lall, U. (2015).  Daily precipitation and tropical moisture exports across the Eastern United States: An application of archetypal analysis to identify spatiotemporal structure.  Journal of Climate, 28(21):8585–8602.
  * [Su et al., 2017] Su, Z., Hao, Z., Yuan, F., Chen, X., and Cao, Q. (2017).  Spatiotemporal variability of extreme summer precipitation over the Yangtze river basin and the associations with climate patterns.  Water, 9(11).
  * [Sun et al., 2017a] Sun, W., Yang, G., Wu, K., Li, W., and Zhang, D. (2017a).  Pure endmember extraction using robust kernel archetypoid analysis for hyperspectral imagery.  ISPRS Journal of Photogrammetry and Remote Sensing, 131:147 – 159.
  * [Sun et al., 2017b] Sun, W., Zhang, D., Xu, Y., Tian, L., Yang, G., and Li, W. (2017b).  A probabilistic weighted archetypal analysis method with Earth mover’s distance for endmember extraction from hyperspectral imagery.  Remote Sensing, 9(8):841.
  * [Tendler et al., 2015] Tendler, A., Mayo, A., and Alon, U. (2015).  Evolutionary tradeoffs, Pareto optimality and the morphology of ammonite shells.  BMC systems biology, 9(1):1–12.
  * [Thøgersen et al., 2013] Thøgersen, J. C., Mørup, M., Damkiær, S., Molin, S., and Jelsbak, L. (2013).  Archetypal analysis of diverse pseudomonas aeruginosa transcriptomes reveals adaptation in cystic fibrosis airways.  BMC Bioinformatics, 14:279.
  * [Thurau et al., 2012] Thurau, C., Kersting, K., Wahabzada, M., and Bauckhage, C. (2012).  Descriptive matrix factorization for sustainability: Adopting the principle of opposites.  Data Mining and Knowledge Discovery, 24(2):325–354.
  * [Tsanousa et al., 2015] Tsanousa, A., Laskaris, N., and Angelis, L. (2015).  A novel single-trial methodology for studying brain response variability based on archetypal analysis.  Expert Systems with Applications, 42(22):8454 – 8462.
  * [Tucker, 1966] Tucker, L. R. (1966).  Some mathematical notes on three-mode factor analysis.  Psychometrika, 31(3):279–311.
  * [van Buuren and Groothuis-Oudshoorn, 2011] van Buuren, S. and Groothuis-Oudshoorn, K. (2011).  mice: Multivariate imputation by chained equations in R.  Journal of Statistical Software, 45(3):1–67.
  * [Van Mechelen et al., 2004] Van Mechelen, I., Bock, H.-H., and De Boeck, P. (2004).  Two-mode clustering methods: a structured overview.  Statistical Methods in Medical Research, 13(5):363–394.
  * [Vinué and Epifanio, 2017] Vinué, G. and Epifanio, I. (2017).  Archetypoid analysis for sports analytics.  Data Mining and Knowledge Discovery, 31(6):1643–1677.
  * [Vinué and Epifanio, 2019] Vinué, G. and Epifanio, I. (2019).  Forecasting basketball players’ performance using sparse functional data.  Statistical Analysis and Data Mining: The ASA Data Science Journal, 12:534–547.
  * [Vinué and Epifanio, 2020] Vinué, G. and Epifanio, I. (2020).  Robust archetypoids for anomaly detection in big functional data.  Advances in Data Analysis and Classification, pages 1–26.
  * [Vinué et al., 2015] Vinué, G., Epifanio, I., and Alemany, S. (2015).  Archetypoids: A new approach to define representative archetypal data.  Computational Statistics & Data Analysis, 87:102 – 115.
  * [Wang and Zhao, 2022] Wang, Y. and Zhao, H. (2022).  Non-linear archetypal analysis of single-cell RNA-seq data by deep autoencoders.  PLoS computational biology, 18(4):e1010025.
  * [Wu et al., 2016] Wu, C., Kamar, E., and Horvitz, E. (2016).  Clustering for set partitioning with a case study in ridesharing.  In IEEE 19th International Conference on Intelligent Transportation Systems (ITSC), pages 1384–1388.
  * [Xie et al., 2019] Xie, J., Ma, A., Fennell, A., Ma, Q., and Zhao, J. (2019).  It is time to apply biclustering: a comprehensive review of biclustering applications in biological and biomedical data.  Briefings in Bioinformatics, 20(4):1450–1465.
  * [Zhao et al., 2012] Zhao, H., Wee-Chung Liew, A., Z Wang, D., and Yan, H. (2012).  Biclustering analysis for pattern discovery: current techniques, comparative studies and applications.  Current Bioinformatics, 7(1):43–55.
  * [Zois et al., 2017] Zois, E. N., Theodorakopoulos, I., and Economou, G. (2017).  Offline handwritten signature modeling and verification based on archetypal analysis.  In 2017 IEEE International Conference on Computer Vision (ICCV), pages 5515–5524.

##  Appendix A Proof of Propositions 1 and 2

Proof of Proposition 1.

Let us denote ${\bf H}_{n\times m}=\alpha_{n\times k}{\bf Z}_{k\times
c}\gamma_{c\times m}$ and ${\bf R}_{n\times c}=\alpha_{n\times k}{\bf
Z}_{k\times c}$; then, each ${\bf r}_{j}^{d}$ $(j=1,\dots,n)$ belongs to the
convex hull $C_{Z}^{d}$ of ${\bf z}_{i}^{d}$ $(i=1,\dots,k)$. Moreover, since
${\bf H}_{n\times m}={\bf R}_{n\times c}{\gamma}_{c\times m}$, each vector
${\bf h}_{j}^{f}$ $(j=1,\dots,m)$ belongs to the convex hull $C_{R}^{f}$ of
the vectors ${\bf r}_{i}^{f}$ $(i=1,\dots,c)$.

Having fixed ${\bf Z}_{k\times c}$; that is, having fixed ${\bf z}_{j}^{d}$
$(j=1,\dots,k)$, $RSS=\left\Arrowvert{\bf X}_{n\times m}-{\bf H}_{n\times
m}\right\Arrowvert^{2}$ is minimized with respect to $\alpha^{\prime}s$ and
$\gamma^{\prime}s$ by choosing $\forall i=1,\dots,n$, ${\bf h}_{i}^{d}$ to be
the point in the convex set defined by $C_{R}^{f}$ and $C_{Z}^{d}$ that is
closest to ${\bf x}_{i}^{d}$.

Suppose without loss of generality that ${\bf z}_{1}^{d}$ is strictly interior
to $C_{V}^{d}$; then choose $t$ such that ${\bf z}(t)={\bf z}_{j}^{d}+t({\bf
z}_{1}^{d}-{\bf z}_{j}^{d})$ is on the boundary of $C_{V}^{d}$. Since the
matrix ${\bf\theta}_{m\times c}$ is considered fixed, the convex hull $C$ of
${\bf z}(t),{\bf z}_{2}^{d},\dots,{\bf z}_{k}^{d}$ contains the convex hull
$C_{Z}^{d}$ of ${\bf z}_{1}^{d},{\bf z}_{2}^{d},\dots,{\bf z}_{k}^{d}$;
therefore, we obtain a larger set over which to minimize
$RSS=\left\Arrowvert{\bf X}_{n\times m}-{\bf H}_{n\times
m}\right\Arrowvert^{2}$.

Proof of Proposition 2.

Now, let us denote ${\bf S}_{k\times m}={\bf Z}_{k\times c}\gamma_{c\times
m}$; then, each ${\bf s}_{j}^{f}$ $(j=1,\dots,m)$ belongs to the convex hull
$C_{Z}^{f}$ of ${\bf z}_{i}^{f}$ $(i=1,\dots,c)$. Moreover, since ${\bf
H}_{n\times m}=\alpha_{n\times k}{\bf S}_{k\times m}$, each vector ${\bf
h}_{j}^{d}$ $(j=1,\dots,n)$ belongs to the convex hull $C_{S}^{d}$ of the
vectors ${\bf s}_{i}^{d}$ $(i=1,\dots,k)$.

As before, having fixed ${\bf Z}_{k\times c}$; that is, having fixed ${\bf
z}_{j}^{d}$ $(j=1,\dots,k)$, $RSS=\left\Arrowvert{\bf X}_{n\times m}-{\bf
H}_{n\times m}\right\Arrowvert^{2}$ is minimized with respect to
$\alpha^{\prime}s$ and $\gamma^{\prime}s$ by choosing $\forall i=1,\dots,n$,
${\bf h}_{i}^{f}$ to be the point in the convex set defined by $C_{S}^{d}$ and
$C_{Z}^{f}$ that is closest to ${\bf x}_{i}^{d}$.

Suppose without loss of generality that ${\bf z}_{1}^{f}$ is strictly interior
to $C_{Y}^{f}$; then choose $t$ such that ${\bf z}(t)={\bf z}_{j}^{f}+t({\bf
z}_{1}^{f}-{\bf z}_{j}^{f})$ is on the boundary of $C_{Y}^{f}$. Since the
matrix ${\bf\beta}_{k\times n}$ is considered fixed, the convex hull $C$ of
${\bf z}(t),{\bf z}_{2}^{f},\dots,{\bf z}_{c}^{f}$ contains the convex hull
$C_{Y}^{f}$ of ${\bf z}_{1}^{f},{\bf z}_{2}^{f},\dots,{\bf z}_{c}^{f}$;
therefore we obtain a larger set over which to minimize
$RSS=\left\Arrowvert{\bf X}_{n\times m}-{\bf H}_{n\times
m}\right\Arrowvert^{2}$.


# Fairness With Minimal Harm:
A Pareto-Optimal Approach For Healthcare

Natalia Martinez† Martin Bertran† Guillermo Sapiro
Department of Electrical and Computer Engineering
Duke University
<EMAIL_ADDRESS>

###### Abstract

Common††$\dagger$ Equal contributions fairness definitions in machine learning
focus on balancing notions of disparity and utility. In this work, we study
fairness in the context of risk disparity among sub-populations. We are
interested in learning models that minimize performance discrepancies across
sensitive groups without causing unnecessary harm. This is relevant to high-
stakes domains such as healthcare, where non-maleficence is a core principle.
We formalize this objective using Pareto frontiers, and provide analysis,
based on recent works in fairness, to exemplify scenarios were perfect
fairness might not be feasible without doing unnecessary harm. We present a
methodology for training neural networks that achieve our goal by dynamically
re-balancing subgroups risks. We argue that even in domains where fairness at
cost is required, finding a non-unnecessary-harm fairness model is the optimal
initial step. We demonstrate this methodology on real case-studies of
predicting ICU patient mortality, and classifying skin lesions from
dermatoscopic images.

##  1 Introduction

Machine learning algorithms play an important role in decision making in
society. When these algorithms are used to make high-impact healthcare
decisions such as deciding whether or not to classify a skin lesion as benign
or not, or predicting mortality for intensive care unit patients, it is
paramount to guarantee that these decisions are accurate and unbiased with
respect to sensitive attributes such as gender or ethnicity. A model that is
trained naively may not have these properties by default (Barocas and Selbst
(2016)). It is desirable in these critical applications to impose some
fairness criteria. There are several lines of work in fairness in machine
learning, such as Demographic Parity (Louizos et al. (2015); Zemel et al.
(2013); Feldman et al. (2015)), Equality of Odds, Equality of Opportunity
(Hardt et al. (2016); Woodworth et al. (2017)), or Calibration (Pleiss et al.
(2017)). These notions of fairness are appropriate in many scenarios, but in
domains where quality of service is paramount, such as healthcare, we argue
that it is necessary to strive for models that are as close to fair as
possible without introducing unnecessary harm to any subgroup (Ustun et al.
(2019)).

In this work, we measure discrimination (unfairness) in terms of difference in
predictive risks across sub-populations defined by our sensitive attributes.
This concept has been explored in other works such as Calders and Verwer
(2010); Dwork et al. (2012); Feldman et al. (2015); Chen et al. (2018); Ustun
et al. (2019). We examine the subset of models from our hypothesis class that
have the best trade-offs between sub-population risks, and select from this
set the one with the smallest risk disparity gap. This is in contrast to
common post-hoc correction methods like the ones proposed in Hardt et al.
(2016); Woodworth et al. (2017) where randomness is potentially added to the
decisions of all sub-populations. While this type of approach diminishes the
accuracy-disparity gap, it does so by potentially introducing randomness into
the final decision (i.e., with some probability, disregard classifier output
and produce arbitrary decision based solely on sensitive label), which leads
to performance degradation in most common risk metrics. Since our proposed
methodology does not require test-time access to sensitive attributes, and can
be applied to any standard classification or regression task, it can also be
used to reduce risk disparity between outcomes, acting as an adaptive risk
equalization loss compatible with unbalanced classification scenarios.

##### Main Contributions

We formalize the notion of no-unnecessary-harm fairness using Pareto
optimality, a state of resource allocations from which it is impossible to
reallocate without making one subgroup worse. We show that finding a Pareto-
fair classifier is equivalent to finding a model in our hypothesis class that
is both Pareto optimal with respect to the sub-population risks (no
unnecessary harm) and minimizes risk disparity. This notion is already
amenable to non-binary sensitive attributes. We analyze Pareto fairness on an
illustrative example, and compare it to alternative approaches. We provide an
algorithm that promotes fair solutions belonging to the Pareto front; this
algorithm can be applied to any standard classifier or regression task.
Finally, we show how our methodology performs on real tasks such as predicting
ICU mortality rates in the MIMIC-III dataset from hospital notes Johnson et
al. (2016), and classifying skin lesions in the HAM10000 dataset Tschandl et
al. (2018).

##  2 Problem Statement

Consider we have access to a dataset
$\mathcal{D}=\\{(x_{i},y_{i},a_{i})\\}_{i=1}^{n}$ containing $n$ independent
triplet samples drawn from a joint distribution $(x_{i},y_{i},a_{i})\sim
P(X,Y,A)$ where $x_{i}\in\mathcal{X}$ are our input features (e.g., images,
tabular data, etc.), $y_{i}\in\mathcal{Y}$ is our target variable, and
$a_{i}\in\mathcal{A}$ indicates group membership or sensitive status (e.g.,
ethnicity, gender); our input features $X$ may or may not explicitly contain
$A$.

Let $h\in\mathcal{H}$ be a classifier from our hypothesis class $\mathcal{H}$
trained to infer $y$ from $x$, $h:\mathcal{X}\rightarrow\mathcal{Y}$; and a
loss function $\ell:\mathcal{Y}\times\mathcal{Y}\rightarrow\mathbb{R}$. We
define the class-specific risk of classifier $h$ on subgroup $a$ as
$R_{a}(h)=\mathbb{E}_{X,Y|A=a}[\ell(Y,h(X))]$. The risk discrimination gap
between two subgroups $a,a^{\prime}\in\mathcal{A}$ is measured as
$\Gamma_{a,a^{\prime}}(h)=|R_{a}(h)-R_{a^{\prime}}(h)|$, and we define the
pairwise discrimination gap as
$\vec{\Gamma}_{\mathcal{A}}(h)=\\{\Gamma_{a,a^{\prime}}(h)\\}_{a,a^{\prime}\in\mathcal{A}}$.
Our goal is to obtain a classifier $h\in\mathcal{H}$ that minimizes this gap
without causing unnecessary harm to any particular group. To formalize this
notion, we define:

######  Definition 2.1.

Pareto front: The set of Pareto front classifiers is defined as
$\mathcal{P}(\mathcal{H},\mathcal{A})=\\{h\in\mathcal{H}:\not\exists
h^{\prime}\in\mathcal{H},\not\exists
a^{\prime}\in\mathcal{A}:[R_{a^{\prime}}(h^{\prime})<R_{a^{\prime}}(h)]\land[R_{a}(h^{\prime})\leq
R_{a}(h)],\forall a\in\mathcal{A}\\}$

######  Definition 2.2.

Pareto-fair classifier and Pareto-fair vector: A classifier $h^{*}$ is an
optimal no-harm classifier if it minimizes the discrimination gap among all
Pareto front classifiers,
$h^{*}=\operatorname*{arg\,min}\limits_{h\in\mathcal{P}(\mathcal{H},\mathcal{A})}||\vec{\Gamma}_{\mathcal{A}}(h)||_{\infty}$.
The Pareto-fair vector is defined as
$\bm{r}^{*}\in\mathbb{R}^{|\mathcal{A}|}:$
$\bm{r}^{*}=\\{R_{a}(h^{*})\\}_{a\in\mathcal{A}}$.

The Pareto front defines the best achievable trade-offs between population
risks $R_{a}(h)$, while the Pareto-fair classifier gives the trade-off with
least disparity. Building on analysis done by Chen et al. (2018); Domingos
(2000), it is possible to decompose the risk in bias, variance and noise for
some given loss functions. The noise represents the smallest achievable risk
for infinitely large datasets. If it differs between sensitive groups, zero-
discrimination (perfect fairness) can be only achieved by introducing bias or
variance, hence doing harm. Figure 1 shows a scenario where the Pareto-front
does not intersect the equality of risk line for the case of two sensitive
groups $a\in\\{0,1\\}$ and a binary output variable $y\in\\{0,1\\}$. Here the
level of noise between subgroups differs, and the Pareto-fair vector
$\bm{r}^{*}$ would not be achieved by either a naive classifier (minimizes
expected global risk), or a classifier where low-occurrence subgroups are
over-sampled (re-balanced naive classifier).

Figure 1: Example of binary target ($y$) and sensitive ($a$) variable. Bottom
left figure shows conditional distributions of observation variable ($x$)
$p(x|a)$, while top left shows distribution of $y$ as a function of $x,a$,
$p(y|x,a)$, these distributions are simple piece-wise constants with two
levels $\rho^{low}$ and $\rho^{high}$, $\Delta$ represents the gap between
level transitions for each group. Middle image shows the Pareto front between
group risks, since noise levels are not the same across subgroups, perfect
fairness is not attainable. Equality of Risk requires pure degradation of
service for group $a=1$. Both Naive and Re-balanced Naive classifiers do not
attain the best possible no-harm classifier. Pareto-fair point is shown in
green. Rightmost figure shows the trade-offs attainable between discrimination
and mean risks, an equivalent problem to the risk trade-off figure.

##  3 Methods

Any loss function that is monotonically decreasing with $R_{a}(h),\forall a$
is minimized by a classifier in the Pareto front. Since we want to minimize
the discrimination gap, we will build an adaptive loss function that is shares
this property with the following form:

$\begin{array}[]{cc}\phi(h;\bm{\mu},c)&=\sum\limits_{a\in\mathcal{A}}R_{a}(h)+\mu_{a}(R_{a}(h)-c)^{+2},\end{array}$
(1)

with $c<\min_{a}R_{a}(h)$ and $\bm{\mu}=\\{\mu_{a}\\}_{a\in\mathcal{A}}$. It
can be shown that for convex Pareto sets (with respect to risk vectors
$\bm{r}=\\{R_{a}(h)\\}_{a\in\mathcal{A}}$), there exist a set of
$\bm{\mu}^{*}$ such that the Pareto-fair vector $\bm{r}^{*}$ is a unique
solution. In Algorithm 1 we jointly search for $\bm{\mu}^{*}$ ,$\bm{r}^{*}$,
and $h^{*}$.

Given: $h_{\theta},\mathcal{D},\gamma>0,\xi\in(0,1),\zeta\in(0,1),\text{lr},$

$\bm{\mu}\leftarrow\bm{0},\;\bm{\mu}^{*}\leftarrow\bm{0},\;,c\leftarrow
0,\;h^{*}\leftarrow h_{\theta}$

while  _Improving_ do

$h_{\theta}\leftarrow\text{SGDwithEarlyStopping}(h_{\theta},\phi,\bm{\mu},c,\mathcal{D})$

// Optimize adaptive loss with fixed $\bm{\mu}$ on dataset $\mathcal{D}$

$r^{\text{val}}\leftarrow\text{EvaluateRisk}(h_{\theta},\mathcal{D}^{\text{val}})$

if  _$||\vec{\Gamma}_{\mathcal{A}}(h)||_{\infty} <\Gamma^{*}$ and
$\bm{r}^{\text{Val}}$ is not dominated by previous validation risks _ then

$h^{*}\leftarrow
h_{\theta},\;\Gamma^{*}\leftarrow||\vec{\Gamma}_{\mathcal{A}}(h)||_{\infty},\;c_{\text{old}}\leftarrow
c,\;c\leftarrow\min_{a}\frac{r_{a}^{\text{Val}}}{k}$

$\bm{\mu}^{*}\leftarrow\bm{\mu}\cdot\frac{(\bm{r}^{\text{Val}}-c_{\text{old}})^{+}}{(\bm{r}^{\text{Val}}-c)^{+}},\;a^{\prime}\leftarrow\operatorname*{arg\,max}_{a}r_{a}^{\text{Val}}$

else

lr $\leftarrow$ $\zeta$ lr,
$\bm{\mu}\leftarrow\bm{\mu}^{*},\;\gamma\leftarrow\gamma\xi,\;h_{\theta}\leftarrow
h^{*}$

end if

$\mu_{a^{\prime}}\leftarrow(1+\gamma)\mu_{a^{\prime}}$

end while

Return: $h^{*}$

// Exit loop due to excessive iterations or no improvement in fairness

Algorithm 1 ParetoFairOptimization

##  4 Experiments and Results

We evaluate our method on mortality prediction and skin lesion classification
and show empirically how we reduce accuracy and risk disparity without
unnecessary harm. We compare it against a naive classifier, class-rebalancing
(equally sampled sub-groups), the post-processing framework presented in Hardt
et al. (2016), and the Disparate Mistreatment framework of Zafar et al.
(2017).

###  4.1 Predicting Mortality in Intensive Care Patients

We analyze clinical notes from adult ICU patients at the Beth Israel Deaconess
Medical Center (MIMIC-III dataset) Johnson et al. (2016) to predict patient
mortality. We follow the pre-processing methodology outlined in Chen et al.
(2018) and use tf-idf statistics on the $10,000$ most frequent words in
clinical notes as input features. Fairness is measured with respect to age
(under/over 55 years old), ethnicity, and outcome. We used a fully connected
neural network with two 2048-unit hidden layers trained with Brier score (BS)
loss. Table 1 shows accuracy and BS of all tested methodologies. We observe
from columns (PF BS) and (PF Acc) that our model has the best accuracy and BS
discrepancies. We can reduce accuracy disparities further by applying Hardt
post processing (HPF Acc), at the cost of reducing overall performance. Note
that, as expected, it is better to apply this post processing on our method
than on the rebalanced naive classifier (HPF Acc vs HRen Acc). This
illustrates the goal of the proposed Pareto-fair paradigm: develop the fairest
algorithm with no unnecessary harm, and if (e.g., due to policy) the resulting
fairness level needs to be improved, use other post-processing techniques on
the Pareto-fair classifier such as Hardt et al. (2016).

Out/Age/Race | Ratio | Na Acc | ReN Acc | Zafar Acc | PF Acc | HReN Acc | HPF Acc | ReN BS | PF BS
---|---|---|---|---|---|---|---|---|---
A/A/NW | 5.7% |  99.1$\pm$0.4% |  86.3$\pm$1.5% |  93.0$\pm$1.4% |  83.4$\pm$2.6% |  76.3$\pm$1.9% |  71.6$\pm$3.2% |  0.2 $\pm$ 0.02 |  0.25 $\pm$ 0.03
A/A/W | 13.3% |  98.8$\pm$0.5% |  86.3$\pm$1.1% |  90.0$\pm$1.3% |  83.2$\pm$1.5% |  76.7$\pm$1.6% |  71.8$\pm$1.9% |  0.2 $\pm$ 0.01 |  0.25 $\pm$ 0.02
A/S/NW | 12.9% |  97.5$\pm$0.6% |  76.5$\pm$1.7% |  81.8$\pm$1.7% |  71.4$\pm$3.0% |  76.4$\pm$2.2% |  71.3$\pm$3.1% |  0.31 $\pm$ 0.02 |  0.36 $\pm$ 0.03
A/S/W | 56.7% |  97.9$\pm$0.3% |  79.0$\pm$0.6% |  77.4$\pm$0.7% |  74.6$\pm$1.6% |  76.2$\pm$2.2% |  72.1$\pm$1.2% |  0.28 $\pm$ 0.01 |  0.34 $\pm$ 0.02
D/A/NW | 0.4% |  23.4$\pm$9.4% |  76.1$\pm$8.5% |  47.7$\pm$9.6% |  78.6$\pm$6.1% |  66.6$\pm$9.9% |  74.1$\pm$9.2% |  0.36 $\pm$ 0.06 |  0.34 $\pm$ 0.04
D/A/W | 0.9% |  32.6$\pm$3.7% |  80.1$\pm$3.3% |  60.5$\pm$6.3% |  83.3$\pm$4.2% |  66.4$\pm$2.4% |  73.6$\pm$4.3% |  0.29 $\pm$ 0.04 |  0.28 $\pm$ 0.03
D/S/NW | 1.8% |  21.4$\pm$2.2% |  66.9$\pm$2.4% |  48.2$\pm$2.0% |  73.3$\pm$2.9% |  64.8$\pm$2.1% |  71.2$\pm$2.9% |  0.42 $\pm$ 0.02 |  0.37 $\pm$ 0.03
D/S/W | 8.3% |  23.4$\pm$2.2% |  67.4$\pm$1.9% |  57.1$\pm$2.2% |  72.5$\pm$3.6% |  66.2$\pm$2.9% |  72.4$\pm$3.5% |  0.42 $\pm$ 0.02 |  0.37 $\pm$ 0.04
Sample Mean | - | 89.5$\pm$0.2% |  78.9$\pm$0.7% |  78.0$\pm$0.7% |  75.7$\pm$1.1% | 75.1$\pm$1.8% |  71.9$\pm$1.2% | 0.29 $\pm$ 0.01 |  0.33 $\pm$ 0.01
Group Mean | 12.5% |  61.8$\pm$1.5% |  77.3$\pm$1.3% |  69.4$\pm$1.3% | 77.5$\pm$0.7% |  71.2$\pm$1.1% | 72.3$\pm$1.0% | 0.31 $\pm$ 0.01 |  0.32 $\pm$ 0.01
Discrepancy | 56.3% |  81.1$\pm$2.5% |  22.5$\pm$2.5% |  49.3$\pm$4.4% | 17.1$\pm$2.6% |  18.6$\pm$3.2% | 13.1$\pm$5.6% |  0.24 $\pm$ 0.03 | 0.16 $\pm$ 0.03
Table 1: Result summary table on MIMIC dataset. Sensitive attributes are the
combination of Outcome (Alive or Deceased), Age (Adult or Senior) and Race
(White or Non-White). Target attribute is Outcome. We report Brier scores (BS)
and Accuracies (Acc) for Naive (Na), Rebalanced Naive (ReN), Zafar (Zafar),
and Pareto-Fair (PF) classifiers; the prefix (H) indicates that Hardt post-
processing was applied to control for Age and Race. Zafar was applied on top
of the embeddings of the Naive classifier, as that produced the best results.

####  4.1.1 Skin Lesion Classification

The HAM10000 dataset Tschandl et al. (2018) collects over $10,000$
dermatoscopic images of skin lesions over a diverse population, lesions are
classified in 7 categories. We used a pretrained DenseNet121 Huang et al.
(2017) as our base classifier, and measured fairness with respect to diagnosis
class, casting balanced risk minimization as a particular use-case for Pareto
fairness. Table 2 shows our empirical results.

Groups | Ratio | PF Acc | ReN Acc | Na Acc | PF BS | ReN BS | Na BS
---|---|---|---|---|---|---|---
akiec | 2.5% | 51.9% | 55.6% | 3.7% | 0.741 | 0.671 | 1.289
bcc | 2.7% | 56.7% | 76.7% | 56.7% | 0.549 | 0.341 | 0.613
bkl | 7.2% | 59.5% | 58.2% | 36.7% | 0.62 | 0.606 | 0.931
df | 0.5% | 66.7% | 33.3% | 0.0% | 0.536 | 0.898 | 1.721
nv | 81.0% | 83.5% | 90.9% | 96.0% | 0.241 | 0.128 | 0.054
vasc | 1.3% | 71.4% | 85.7% | 0.0% | 0.36 | 0.246 | 1.604
mel | 4.8% | 53.8% | 48.1% | 32.7% | 0.586 | 0.657 | 0.941
Sample Mean | - | 78.6% |  84.8% | 83.6% | 0.308 | 0.213 | 0.234
Group Mean | 14.3% | 63.4% | 64.1% | 32.3% | 0.519 | 0.507 | 1.022
Discrepancy | 80.4% |  31.7% | 57.5% | 96.0% | 0.501 | 0.769 | 1.667
Table 2:  Group-specific risks and accuracies for Naive (Na), Rebalanced Naive
(ReN), and Pareto-Fair (PF) classifiers. Lesion Type was chosen as both target
and sensitive variable, they are classified as Actinic keratoses and
intraepithelial carcinoma (akiec), basal cell carcinoma (bcc), benign
keratosis-like lesions (bkl), dermatofibroma (df), melanoma (mel), melanocytic
nevi (nv) and vascular lesions (vasc). The Pareto-fair classifier exhibits the
lowest discrepancy on Acc and BS while still being Pareto-optimal.

##  5 Discussion

Here we explore the problem of reducing the risk disparity gaps in the most
ethical way possible (i.e., minimizing unnecessary harm). We provide an
algorithm that can be used with any standard classification or regression loss
to bridge risk disparity gaps without introducing unnecessary harm. We show
its performance on two real-world case studies; we take advantage of the fact
that our method does not require test-time access to sensitive attributes to
frame balanced classification as a fairness problem. In future work, we wish
to analyze if we can automatically identify high-risk sub-populations as part
of the learning process and attack risk disparities as they arise, rather than
relying on preexisting notions of disadvantaged groups or populations. We
believe that no-unnecessary-harm notions of fairness are of great interest for
several applications, especially so on domains such as healthcare.

## References

  * Barocas and Selbst (2016) Solon Barocas and Andrew D Selbst.  Big data’s disparate impact.  _Calif. L. Rev._ , 104:671, 2016.
  * Calders and Verwer (2010) Toon Calders and Sicco Verwer.  Three naive bayes approaches for discrimination-free classification.  _Data Mining and Knowledge Discovery_ , 21(2):277–292, 2010.
  * Chen et al. (2018) Irene Chen, Fredrik D Johansson, and David Sontag.  Why is my classifier discriminatory?  In _Advances in Neural Information Processing Systems_ , pages 3539–3550, 2018.
  * Domingos (2000) Pedro Domingos.  A unified bias-variance decomposition.  In _Proceedings of 17th International Conference on Machine Learning_ , pages 231–238, 2000.
  * Dwork et al. (2012) Cynthia Dwork, Moritz Hardt, Toniann Pitassi, Omer Reingold, and Richard Zemel.  Fairness through awareness.  In _Proceedings of the 3rd innovations in theoretical computer science conference_ , pages 214–226. ACM, 2012.
  * Feldman et al. (2015) Michael Feldman, Sorelle A Friedler, John Moeller, Carlos Scheidegger, and Suresh Venkatasubramanian.  Certifying and removing disparate impact.  In _Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pages 259–268. ACM, 2015.
  * Hardt et al. (2016) Moritz Hardt, Eric Price, Nati Srebro, et al.  Equality of opportunity in supervised learning.  In _Advances in neural information processing systems_ , pages 3315–3323, 2016.
  * Huang et al. (2017) Gao Huang, Zhuang Liu, Laurens Van Der Maaten, and Kilian Q Weinberger.  Densely connected convolutional networks.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 4700–4708, 2017.
  * Johnson et al. (2016) Alistair EW Johnson, Tom J Pollard, Lu Shen, H Lehman Li-wei, Mengling Feng, Mohammad Ghassemi, Benjamin Moody, Peter Szolovits, Leo Anthony Celi, and Roger G Mark.  Mimic-iii, a freely accessible critical care database.  _Scientific data_ , 3:160035, 2016.
  * Louizos et al. (2015) Christos Louizos, Kevin Swersky, Yujia Li, Max Welling, and Richard Zemel.  The variational fair autoencoder.  _arXiv preprint arXiv:1511.00830_ , 2015.
  * Pleiss et al. (2017) Geoff Pleiss, Manish Raghavan, Felix Wu, Jon Kleinberg, and Kilian Q Weinberger.  On fairness and calibration.  In _Advances in Neural Information Processing Systems_ , pages 5680–5689, 2017.
  * Tschandl et al. (2018) Philipp Tschandl, Cliff Rosendahl, and Harald Kittler.  The ham10000 dataset, a large collection of multi-source dermatoscopic images of common pigmented skin lesions.  _Scientific data_ , 5:180161, 2018.
  * Ustun et al. (2019) Berk Ustun, Yang Liu, and David Parkes.  Fairness without harm: Decoupled classifiers with preference guarantees.  In _International Conference on Machine Learning_ , pages 6373–6382, 2019.
  * Woodworth et al. (2017) Blake Woodworth, Suriya Gunasekar, Mesrob I Ohannessian, and Nathan Srebro.  Learning non-discriminatory predictors.  _arXiv preprint arXiv:1702.06081_ , 2017.
  * Zafar et al. (2017) Muhammad Bilal Zafar, Isabel Valera, Manuel Rodriguez, Krishna Gummadi, and Adrian Weller.  From parity to preference-based notions of fairness in classification.  In _Advances in Neural Information Processing Systems_ , pages 229–239, 2017.
  * Zemel et al. (2013) Rich Zemel, Yu Wu, Kevin Swersky, Toni Pitassi, and Cynthia Dwork.  Learning fair representations.  In _International Conference on Machine Learning_ , pages 325–333, 2013.