marginparsep has been altered.
topmargin has been altered.
marginparwidth has been altered.
marginparpush has been altered.
The page layout violates the ICML style. Please do not change the page layout,
or include packages like geometry, savetrees, or fullpage, which change it for
you. We’re not able to reliably undo arbitrary changes to the style. Please
remove the offending package(s), or layout-changing commands and try again.

Rethinking Positive Aggregation and Propagation of Gradients in Gradient-based
Saliency Methods

Anonymous Authors1

††footnotetext:  1Anonymous Institution, Anonymous City, Anonymous Region,
Anonymous Country. Correspondence to: Anonymous Author
<EMAIL_ADDRESS>
Preliminary work. Under review by the International Conference on Machine
Learning (ICML). Do not distribute.

###### Abstract

Saliency methods interpret the prediction of a neural network by showing the
importance of input elements for that prediction. A popular family of saliency
methods utilize gradient information. In this work, we empirically show that
two approaches for handling the gradient information, namely positive
aggregation, and positive propagation, break these methods. Though these
methods reflect visually salient information in the input, they do not explain
the model prediction anymore as the generated saliency maps are insensitive to
the predicted output and are insensitive to model parameter randomization.
Specifically for methods that aggregate the gradients of a chosen layer such
as GradCAM++ and FullGrad, exclusively aggregating positive gradients is
detrimental. We further support this by proposing several variants of
aggregation methods with positive handling of gradient information. For
methods that backpropagate gradient information such as LRP, RectGrad, and
Guided Backpropagation, we show the destructive effect of exclusively
propagating positive gradient information.

##  1 Introduction

Figure 1: Saliency maps derived from different aggregation methods (columns)
for different output class predictions (rows). All methods except GradCAM use
positive aggregation and are not class-sensitive.

The quest for understanding the basis of neural networks’ predictions is
gaining momentum as the need for interpretability intensifies especially in
sensitive application domains such as medicine, finance, policy, and law. One
approach for understanding the predictions is using saliency methods (aka
attribution methods) which assign each input element (e.g. pixels) its
importance for the corresponding output prediction. There exist a plethora of
proposed saliency methods, and thus raising the question “which explanation
method to trust?”. Visual evaluation of saliency maps has been shown to be
unreliable as the method can generate human interpretable saliency maps but
not explaining model behavior Adebayo et al. (2018); Kindermans et al. (2019);
Nie et al. (2018); Hooker et al. (2019). One way to evaluate saliency methods
is by checking whether they satisfy certain desirable properties. Some of
these properties are formulated as axioms Sundararajan et al. (2017); Lundberg
& Lee (2017), such as efficiency (completeness), dummy (null-player),
symmetry, sensitivity, and implementation invariance. It is therefore
plausible to theoretically show whether methods abide by the axioms, though
there is the caveat that while a method’s formulation satisfies axioms, the
required practical assumptions and approximations can break them Sundararajan
& Najmi (2019).

It is also possible to formulate desirable properties as experiments (sanity
checks) and empirically test whether a method violates those properties. Two
important properties that a saliency method is required to satisfy are the
sensitivity of the saliency map to the predicted output class Nie et al.
(2018), aka class-sensitivity, and the sensitivity of the saliency map to
model parameter randomization Adebayo et al. (2018), i.e. the saliency map for
the model before randomizing its weights should be different from the saliency
map after the model parameters are randomized. Not satisfying one of these two
properties signifies that the method cannot explain model behavior.

In this work, we empirically show that two common approaches for handling
gradient information in gradient-based saliency methods, methods that utilize
gradient information to explain the model, can result in violating class-
sensitivity and randomization-sensitivity. We show that in methods which
aggregate gradients in a layer, if this aggregation is done only on positive
gradients (e.g. via using a ReLU function) the desired properties are not
satisfied. We show this phenomenon for FullGrad Srinivas & Fleuret (2019) and
GradCAM++ Chattopadhay et al. (2018), and other positive aggregation methods
that we propose: positive aggregated GradCAM Selvaraju et al. (2020) for
various layers, positive aggregated Gradients for various layers, and a
version of both where results from all layers are aggregated. Another approach
that we investigate is inherent in methods that propagate gradient information
such as LRP Montavon et al. (2017); Bach et al. (2015) (and Excitation
Backpropagation Zhang et al. (2018) as it is equivalent to LRP-$\alpha 1\beta
0$), RectGrad Kim et al. (2020), and Guided Backpropagation Springenberg et
al. (2015). We show that, when only positive information is backpropagated,
the properties are not satisfied. We investigate these hypotheses empirically
via ablative experiments. For sensitivity to model parameter randomization, we
use the sanity checks of Adebayo et al. (2018). For evaluating class
sensitivity we use the class sensitivity metric proposed by Rebuffi et al.
(2020) and the pointing game ofZhang et al. (2018); Fong et al. (2019). We
also propose the restricted pointing game for evaluating class sensitivity.

Rebuffi et al. (2020) propose a general formulation for aggregation methods,
however, do not discuss how the choice of filtering functions (such as
$|\cdot|$ or $\operatorname{ReLU}$) affect them. It is generally observed that
these aggregations are not class sensitive. In this work, we reveal the
culprit and also show positive filtering affects the propagation methods.

##  2 Experimental Setup

The goal in experiments is to evaluate class-sensitivity and sensitivity to
model parameter randomization. For this purpose we use the following
experiments (we use the VGG-16 network Simonyan & Zisserman (2014)):
Pointing game. Zhang et al. (2018); Fong et al. (2019). Saliency maps are
generated for each class present in the image. For each class, if the maximum
value in the generated saliency map is inside the ground truth mask the method
has correctly pointed to the class. The accuracy is the ratio between
correctly pointed classes and all pointing trials. As we aim to study class
sensitivity we use a subset of PASCAL VOC07 Everingham et al. (2015) dataset
where there are _multiple object classes present in each image_ (as defined in
Zhang et al. (2018)).
Restricted pointing game. As we are using a dataset with multiple object
classes in each image, class-sensitive methods can generally achieve better
accuracy in pointing game. However, the major problem with the pointing game
experiment is that it does not disentangle class-sensitivity from
localization. Therefore, comparing different methods with each other is
difficult. We propose the restricted pointing game to tackle this problem. For
each image and all its object classes, we only generate the saliency map for
the max output class. The ground truth masks of all objects in the image are
all compared against the max saliency map. If the resulting accuracy for the
restricted pointing game is similar to the accuracy of the original pointing
game, it can show that the saliency maps for different object classes are
similar and are pointing to the same object class.
Class-sensitivity. Nie et al. (2018) showed that for Guided Backpropagation
and deconvolution methods, the generated saliency maps for different output
classes in the same image are visually indiscernible. Motivated by this
experiment, Rebuffi et al. (2020) proposed a quantitative evaluation for
class-sensitivity. For a given attribution method, saliency maps are generated
for the max and min classes, and the correlation between them is computed. We
report the average on all images in PASCAL VOC07.
Parameter randomization. Adebayo et al. (2018) proposed a sanity check
experiment for evaluating the methods’ sensitivity to model parameter
randomization. The similarity between saliency maps is measured by Spearman on
HOG. The model parameters are increasingly randomized by reinitializing
($\sim\mathcal{N}(0,0.01)$) the parameters of layers in a cascading fashion
from the last layer to the first. We use a set of 1000 images of ImageNet Deng
et al. (2009).

##  3 Positive Aggregation

Figure 2: Saliency maps for different output predictions (rows) using GradCAM
with positive aggregation (GradCAM_Pos) on different layers of VGG-16. Left
column being the final layer.

The locally connected structure of convolutions results in internal feature
maps that retain the spatial information of the input, and therefore an
attribution of contributions to these feature maps can be scaled to input
space. Several methods utilize the gradient of the output with respect to
activation units (GradCAM, GradCAM++) or biases (FullGrad) in a chosen layer.
The gradient information in the chosen layer is aggregated across the channel
dimension. In this section, we show how positive aggregation — using an
absolute function (FullGrad) or using $\operatorname{ReLU}$ (GradCAM++) before
the aggregation — can result in feature maps that recover salient objects in
the image without being related to the output prediction being explained.

We also show that simply adding an absolute function to the formulation of
GradCAM Selvaraju et al. (2020) before aggregation can result in the same
behavior. The positive aggregation allows for GradCAM to be applicable to
earlier layers as well, but we show that it is only recovering salient feature
information without any regard for the prediction. The original formulation of
GradCAM is as follows:

$S_{\text{GradCAM}}=\operatorname{ReLU}(\sum_{k}w_{k}A^{k})$ (1)

where $A^{k}$ is the activation map at channel $k$, $f$ is the output being
explaind, and

$w_{k}=\frac{1}{Z}\sum_{i}\sum_{j}\frac{\partial f}{\partial A_{ij}^{k}}$ (2)

where $Z$ is the number of activation units in $A^{k}$.

Figure 3: Saliency maps generated for different output predictions (rows)
using GradCAM++ on different layers (columns) of VGG-16. Earlier layers
presented on the right.

This method is well justified for the final convolutional layer, as such
linear approximation of output with respect to activations, which is the basis
for the aggregation of feature maps in GradCAM holds stronger for the last
layer than previous layers (in a network where CAM Zhou et al. (2016) is
applicable, linear aggregation exactly derives the output). However, if
GradCAM is computed for earlier layers, the resulting saliency maps are
arbitrary as reported in Selvaraju et al. (2020); Rebuffi et al. (2020).
During the aggregation phase, GradCAM++ uses $\operatorname{ReLU}$ on
gradients before the summation directly in the main formulation of the method
as follows:

$S_{\text{GradCAM++}}=\sum_{k}w_{k}A^{k}$ (3)

where,

$w_{k}=\sum_{i}\sum_{j}\alpha_{ij}^{k}\operatorname{ReLU}(\frac{\partial
f}{\partial A_{ij}^{k}})$ (4)

and $\alpha_{ij}^{k}$ is a coefficient computed for each activation unit
(please refer to Chattopadhay et al. (2018)).

Fullgrad’s Srinivas & Fleuret (2019) main formulation distributes the
contribution to the output prediction between all inputs to the network (input
and mid-level biases) and there is no $\operatorname{ReLU}$ or absolute
summation. However such formulation cannot be visualized as a saliency map.
Thus an aggregation phase is added for CNN visualizations, and an absolute
function is used within the post-processing $\psi$:

$S_{\text{FullGrad}}=\psi(\frac{\partial f}{\partial x}\odot x)+\sum_{l\in
L}\sum_{k}\psi(\frac{\partial f}{\partial A^{k}}\odot b^{l})$ (5)

where $x$ is the input image, and $b^{l}$ stands for the bias parameters of
layer $l$. The absolute operation is added to visualize only the magnitude of
importance while ignoring the sign. However, ignoring the sign results in
aggregating the gradients related to all features in the image (compare this
to CAM/GradCAM where the summation is done without $|\cdot|$ because the
linear combination of activations in the last layer recovers the output).

Table 1: Pointing game on aggregation methods
| Original | Restricted
---|---|---
GradCAM | 74.1 | 48.9
FullGrad | 58.4 | 50.5
GradCAM++ | 66.2 | 54.1
GradCAM++ B4 | 36.3 | 37.1
GradCAM++ B3 | 35.2 | 36.4
GradCAM++ B2 | 38.0 | 39.2
GradCAM++ B1 | 40.4 | 40.3
GradCAM_Pos | 60.4 | 58.1
GradCAM_Pos B4 | 37.6 | 37.8
GradCAM_Pos B3 | 31.8 | 31.6
GradCAM_Pos B2 | 40.6 | 40.7
GradCAM_Pos B1 | 40.7 | 40.6
GradMid | 63.9 | 56.3
GradMid B4 | 65.4 | 56.3
GradMid B3 | 60.0 | 53.8
GradMid B2 | 53.8 | 48.8
GradMid B1 | 51.8 | 44.5
Cumulative_GradCAM | 58.1 | 56.2
Cumulative_GradMid | 64.1 | 56.3

Initially, we show that if we change the formulation of GradCAM by adding an
absolute function during aggregation, which we call GradCAM_Pos, the resulting
saliency maps for all layers show salient image features, however, it seems
that these are merely silhouettes of salient image features and we observe
that these features are not relevant to the predicted class. A visual example
is presented in Fig. 2, 1. The maps are computed for different layers of VGG.
We use the convolutional maps at different resolution blocks, and represent
them as GradCAM_Pos plus block number, where B1 for instance represents the
first block. It is observed that as we move toward earlier layers the
generated maps recover low-level features of objects such as edges. _The maps
on earlier layers should not be interpreted as ’high resolution’, as they are
not attribution or saliency maps anymore._ It is also reported by Rebuffi et
al. (2020) that earlier layers result in class-insensitive maps in aggregation
methods, however, we are considering the role positive aggregation plays.
_With positive aggregation, even in the final layer, maps highlight salient
image features but for both classes present in the image._ This observation is
further supported in the pointing game experiments in Table 1 as the original
and restricted games accuracies are similar. Class-sensitivity metric in Table
2 shows that for all layers, except the final layer, the correlation is close
to 1. For the final layer the correlation is lower but still _significantly
higher than GradCAM with absolute function (0.69 vs. 0.03), showing how
detrimental the effect of positive aggregation is._ The positive aggregation
in this case also has a significant effect on the method’s sensitivity to
randomization. Fig. 5 and 4 shows that _these method’s saliency maps do not
change after randomization, further supporting the claim that methods are
recovering image features._

Figure 4: Sanity checks for sensitivity to parameter randomization for
aggregation (Up) and propagation (Down) methods.

As stated, the ability to generate saliency maps that look interpretable on
earlier layers for GradCAM can be achieved via positive aggregation. GradCAM++
also uses positive aggregation and its visual examples for the last layer and
preceding layers are presented in Fig 3 and Fig. 1. The observations for
GradCAM_Pos are visible here as well. The results on the restricted pointing
game and original are equivalent when experiments are done on all layers
except the final layer. However, correlations (Fig. 2) are better than the
GradCAM_Pos case, but still significantly high compared to GradCAM. These
lower correlations are expected when looking at the visual results as there
are arbitrary highlighted areas in the maps as well, but the method still
highlights all the salient features in the image, which is confirmed by
pointing game experiments (Table 1). In parameter randomization in Fig. 4 and
5 _the method gets more insensitive when applied to earlier layers._

Table 2: Class-sensitivity metric on aggregation methods
| Spearman HOG
---|---
GradCAM | 0.03
FullGrad | 0.48
GradCAM++ | 0.37
GradCAM++ B4 | 0.23
GradCAM++ B3 | 0.44
GradCAM++ B2 | 0.58
GradCAM++ B1 | 0.65
GradCAM_Pos | 0.69
GradCAM_Pos B4 | 0.94
GradCAM_Pos B3 | 0.98
GradCAM_Pos B2 | 0.99
GradCAM_Pos B1 | 0.99
GradMid | 0.29
GradMid B4 | 0.35
GradMid B3 | 0.32
GradMid B2 | 0.34
GradMid B1 | 0.43
Cumulative_GradCAM | 0.92
CumulativeMid | 0.36

Fullgrad has one major difference in formulation compared to the other
methods. The method has an extra aggregation step where maps from all _layers_
are aggregated. It is expected that is aggregations compounds the destructive
effect of positive aggregation within each layer. In this section, we study
the effect of aggregation on all layers. First, we propose an all-layer
aggregation on GradCAM_Pos and call the method CumulativeGradCAM. As can be
seen in Table 2 and 1 the method is not class sensitive. The result seems to
be bounded by the worst-case early layer and best case final layer for
GradCAM_Pos. This also explains why FullGrad is less class sensitive than
GradCAM++. CumulativeGradCAM uses the model’s activations, therefore we also
propose a variant where we only aggregate gradients. When done on one layer we
call this as GradMid and when done on all layers we call it
CumuluativeGradMid. This is to show that _positive aggregation on gradients
alone without any weighting can also generate visually interpretable feature
maps, though similar to other cases the visualization only shows salient image
features irrespective of the model and the output._ The results for these all-
layer summation methods are presented alongside all the aforementioned methods
in Fig. 1. _It is also observed in Fig. 4 that Fullgrad is insensitive to
randomization considerably, and the CumulativeGradCAM is the least sensitive._

Figure 5: From left: Saliency maps for GradCAM++ B5, B3 and B1 layers.
followed by GradCAM_Pos B5, B3, B1. From top: Original image, saliency map on
original model, half randomized model and fully randomized model.

##  4 Positive Propagation

Table 3: Pointing game on propagation methods
| Original | Restricted
---|---|---
Guided BackProp | 48.7 | 47.5
Gradients | 52.7 | 45.1
RectGrad | 51.2 | 50.0
RectGrad_mod | 59.5 | 53.8
LRP-$\gamma=\infty$ | 50.6 | 50.0
LRP-0 | 51.1 | 44.2
Figure 6: Saliency maps generated for different output predictions (rows)
using different propagation methods.

Gradient information backpropagation towards the input from the the output
being explained is the underlying principle of a class of gradient-based
saliency methods. In this section we study the effect of rules that propagate
only positive gradient information. Considering a ReLU based neural network,
let $a_{i}^{l}$ denote an activation unit $i$ in layer $l$ and the flowing
gradient into $a_{i}^{l}$ be $R_{i}^{l+1}$ and let the gradient backpropagated
from $a_{i}^{l}$ be $R_{i}^{l}$. In normal gradient backpropagation
$R_{i}^{l}=\mathbb{I}(a_{i}^{l})R_{i}^{l+1}$, where $\mathbb{I}(.)$ is the
indicator function. Guided Backpropagation rule is defined as:

$R_{i}^{l}=\mathbb{I}(a_{i}^{l}R_{i}^{l+1}>0)R_{i}^{l+1}$ (6)

and RectGrad uses the following propagation rule:

$R_{i}^{l}=\mathbb{I}(a_{i}^{l}R_{i}^{l+1}>\tau)R_{i}^{l+1}$ (7)

where $\tau$ is a threshold value.

Layerwise Relevance Propagation (LRP) has several propagation rules and a
widely adopted rule is the LRP-$\alpha 1\beta 0$, which is equivalent Samek et
al. (2019) to Deep Taylor $\alpha 1\beta 0$ Montavon et al. (2017), LRP z+-
rule and Excitation-Backprop Zhang et al. (2018) methods. This rule is also
equivalent to LRP-$\gamma$ with $\gamma=\infty$. LRP-$\gamma$ is introduced to
favor the effect of positive contributions over negative ones and is defined
as follows:

$R_{j}=\sum_{k}\frac{a_{j}(w_{jk}+\gamma
w_{jk}^{+})}{\sum_{j}a_{j}(w_{jk}+\gamma w_{jk}^{+})}R_{k}$ (8)

where $R_{j}$ is the relevance of neuron $j$ for prediction and $R_{k}$ is
relevance of neuron $k$ in the next layer. The controlling parameter for
considering positive contributions is $\gamma$. As $\gamma$ is increased the
effect is more pronounced. In this work we consider the limits $\gamma=\infty$
and $\gamma=0$ (LRP-0).

Table 4: Class-sensitivity metric on propagation methods
| Spearman HOG
---|---
Guided BackProp | 0.99
Gradients | 0.44
RectGrad | 0.99
RectGrad_Mod | 0.72
LRP-$\gamma=\infty$ | 0.99
LRP-0 | 0.59

For all these methods, we show the effect of positive gradient information
backpropagation by running the experiments in section 2 on these methods and
their counterpart versions where there is no bias towards positive gradients.
The counterpart version for Guided Backpropagation is changing the rule to
$R_{i}^{l}=\mathbb{I}(a_{i}^{l}>0)R_{i}^{l+1}$ so that both positive and
negative information is propagated and evidently this is equivalent to
unmodified gradients. The counterpart version for RectGrad is achievable by
using an absolute function in the formulation so that large positive and
negative gradients can both flow backwards, therefore the rule is modified to
$R_{i}^{l}=\mathbb{I}(\left|a_{i}^{l}R_{i}^{l+1}\right|>\tau)R_{i}^{l+1}$,
which we call RectGrad_Mod. Setting $\gamma=0$ in LRP-$\gamma$ removes the
bias towards positive information backpropagation.

We first investigate the effect of positive backpropagation on the sensitivity
of the methods to model parameter randomization. Adebayo et al. (2018); Nie et
al. (2018) show that the Guided Backpropagation method is insensitive to
parameter randomization. Khakzar et al. (2019) report insensitivity to
randomization of RectGrad using sanity checks. Our sanity check experiment
reported in Fig. 4 further shows that _RectGrad and LRP- $\gamma=\infty$ are
also insensitive to the same extent as Guided Backpropagation_. As the
resulting saliency maps resemble the images before and after randomization, it
points to the fact that such _positive gradient information backpropagation
recovers salient features in the image regardless of the model_. The results
for the counterpart versions of all methods signify that these counterpart
methods are sensitive to randomization.

The pointing game experiment in Table 3 shows that the resulting accuracy on
the original is similar to accuracy on restricted pointing games are for the
positive propagation methods. _This signifies that these methods are pointing
to the same salient object class in the images for different output
predictions_. On the other hand, for the counterpart versions, the accuracy
drops in the restricted version, implying that indeed the saliency maps for
different outputs differ. However, it still seems that these counterpart
methods (e.g. Gradients) are also class-insensitive for many images (Fig. 6),
and effect which needs to be studied in future, especially when other works
Rebuffi et al. (2020); Nie et al. (2018) have contradicting statements in this
regard.

The class-sensitivity metric in Table 4 shows that there is a significant
similarity between saliency maps of different outputs for the positive
backpropagation methods. This further confirms the observation that _these
methods are recovering salient image features rather than explaining the
output prediction._ It is also observed that for counterpart methods, the
saliency maps change for different output predictions (Table 3, 4), though it
does not mean that they are class sensitive (Fig. 6). _Nevertheless, positive
propagation makes these methods class-insensitive._

##  5 Conclusion

In this work, we empirically showed that positive aggregation or propagation
of gradients in gradient-based saliency methods results in saliency maps that
recover salient image features regardless of the model and output prediction,
and changes the methods towards being class-insensitive and insensitive to
randomization.

## References

  * Adebayo et al. (2018) Adebayo, J., Gilmer, J., Muelly, M., Goodfellow, I., Hardt, M., and Kim, B.  Sanity checks for saliency maps.  In _Advances in Neural Information Processing Systems_ , 2018.
  * Bach et al. (2015) Bach, S., Binder, A., Montavon, G., Klauschen, F., Müller, K. R., and Samek, W.  On pixel-wise explanations for non-linear classifier decisions by layer-wise relevance propagation.  _PLoS ONE_ , 2015.  ISSN 19326203.  doi: 10.1371/journal.pone.0130140.
  * Chattopadhay et al. (2018) Chattopadhay, A., Sarkar, A., Howlader, P., and Balasubramanian, V. N.  Grad-CAM++: Generalized gradient-based visual explanations for deep convolutional networks.  In _Proceedings - 2018 IEEE Winter Conference on Applications of Computer Vision, WACV 2018_ , 2018.  ISBN 9781538648865.  doi: 10.1109/WACV.2018.00097.
  * Deng et al. (2009) Deng, J., Dong, W., Socher, R., Li, L.-J., Li, K., and Fei-Fei, L.  Imagenet: A large-scale hierarchical image database.  In _2009 IEEE conference on computer vision and pattern recognition_ , pp. 248–255. Ieee, 2009.
  * Everingham et al. (2015) Everingham, M., Eslami, S. M. A., Van Gool, L., Williams, C. K. I., Winn, J., and Zisserman, A.  The pascal visual object classes challenge: A retrospective.  _International Journal of Computer Vision_ , 111(1):98–136, January 2015.
  * Fong et al. (2019) Fong, R., Patrick, M., and Vedaldi, A.  Understanding deep networks via extremal perturbations and smooth masks.  In _Proceedings of the IEEE International Conference on Computer Vision_ , pp. 2950–2958, 2019.
  * Hooker et al. (2019) Hooker, S., Erhan, D., Kindermans, P.-J., and Kim, B.  A benchmark for interpretability methods in deep neural networks.  In _Advances in Neural Information Processing Systems_ , pp. 9737–9748, 2019.
  * Khakzar et al. (2019) Khakzar, A., Baselizadeh, S., Khanduja, S., Kim, S. T., and Navab, N.  Explaining neural networks via perturbing important learned features.  _arXiv preprint arXiv:1911.11081_ , 2019.
  * Kim et al. (2020) Kim, B., Seo, J., Jeon, S., Koo, J., Choe, J., and Jeon, T.  Why are Saliency Maps Noisy? Cause of and Solution to Noisy Saliency Maps.  2020\.  doi: 10.1109/iccvw.2019.00510.
  * Kindermans et al. (2019) Kindermans, P. J., Hooker, S., Adebayo, J., Alber, M., Schütt, K. T., Dähne, S., Erhan, D., and Kim, B.  The (Un)reliability of Saliency Methods.  In _Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)_. 2019.  doi: 10.1007/978-3-030-28954-6˙14.
  * Lundberg & Lee (2017) Lundberg, S. M. and Lee, S. I.  A unified approach to interpreting model predictions.  In _Advances in Neural Information Processing Systems_ , 2017.
  * Montavon et al. (2017) Montavon, G., Lapuschkin, S., Binder, A., Samek, W., and Müller, K. R.  Explaining nonlinear classification decisions with deep Taylor decomposition.  _Pattern Recognition_ , 2017.  ISSN 00313203.  doi: 10.1016/j.patcog.2016.11.008.
  * Nie et al. (2018) Nie, W., Zhang, Y., and Patel, A. B.  A theoretical explanation for perplexing behaviors of backpropagation-based visualizations.  In _35th International Conference on Machine Learning, ICML 2018_ , 2018.  ISBN 9781510867963.
  * Rebuffi et al. (2020) Rebuffi, S.-A., Fong, R., Ji, X., and Vedaldi, A.  There and Back Again: Revisiting Backpropagation Saliency Methods.  apr 2020.  URL http://arxiv.org/abs/2004.02866.
  * Samek et al. (2019) Samek, W., Montavon, G., Vedaldi, A., Hansen, L. K., and Muller, K.-R.  Explainable AI: Interpreting, Explaining and Visualizing Deep Learning.  _Lecture Notes in Computer Science_ , 2019.  doi: 10.1007/978-3-030-28954-6.
  * Selvaraju et al. (2020) Selvaraju, R. R., Cogswell, M., Das, A., Vedantam, R., Parikh, D., and Batra, D.  Grad-CAM: Visual Explanations from Deep Networks via Gradient-Based Localization.  _International Journal of Computer Vision_ , 2020.  ISSN 15731405.  doi: 10.1007/s11263-019-01228-7.
  * Simonyan & Zisserman (2014) Simonyan, K. and Zisserman, A.  Very deep convolutional networks for large-scale image recognition.  _arXiv preprint arXiv:1409.1556_ , 2014.
  * Springenberg et al. (2015) Springenberg, J. T., Dosovitskiy, A., Brox, T., and Riedmiller, M.  Striving for simplicity: The all convolutional net.  In _3rd International Conference on Learning Representations, ICLR 2015 - Workshop Track Proceedings_ , 2015.
  * Srinivas & Fleuret (2019) Srinivas, S. and Fleuret, F.  Full-Gradient Representation for Neural Network Visualization.  Technical report, 2019.
  * Sundararajan & Najmi (2019) Sundararajan, M. and Najmi, A.  The many shapley values for model explanation.  _arXiv preprint arXiv:1908.08474_ , 2019.
  * Sundararajan et al. (2017) Sundararajan, M., Taly, A., and Yan, Q.  Axiomatic attribution for deep networks.  In _34th International Conference on Machine Learning, ICML 2017_ , 2017.  ISBN 9781510855144.
  * Zhang et al. (2018) Zhang, J., Bargal, S. A., Lin, Z., Brandt, J., Shen, X., and Sclaroff, S.  Top-down neural attention by excitation backprop.  _International Journal of Computer Vision_ , 126(10):1084–1102, 2018.
  * Zhou et al. (2016) Zhou, B., Khosla, A., Lapedriza, A., Oliva, A., and Torralba, A.  Learning Deep Features for Discriminative Localization.  In _Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition_ , 2016.  ISBN 9781467388504.  doi: 10.1109/CVPR.2016.319.


# The geometry of financial institutions -
Wasserstein clustering of financial data

Lorenz Riess
University of Vienna and
Oesterreichische Nationalbank &Mathias Beiglboeck
University of Vienna &Johannes Temme
Oesterreichische Nationalbank &Andreas Wolf
Oesterreichische Nationalbank &Julio Backhoff-Veraguas
University of Vienna

###### Abstract

The increasing availability of granular and big data on various objects of
interest has made it necessary to develop methods for condensing this
information into a representative and intelligible map. Financial regulation
is a field that exemplifies this need, as regulators require diverse and often
highly granular data from financial institutions to monitor and assess their
activities. However, processing and analyzing such data can be a daunting
task, especially given the challenges of dealing with missing values and
identifying clusters based on specific features.

To address these challenges, we propose a variant of Lloyd’s algorithm that
applies to probability distributions and uses generalized Wasserstein
barycenters to construct a metric space which represents given data on various
objects in condensed form. By applying our method to the financial regulation
context, we demonstrate its usefulness in dealing with the specific challenges
faced by regulators in this domain. We believe that our approach can also be
applied more generally to other fields where large and complex data sets need
to be represented in concise form.

##  1 Introduction

The contribution of this article is to provide an algorithm which takes as
input a number of (discrete) probability distributions in a high dimensional
space which are then _clustered and represented_ as points in a metric space,
where the distance gauges differences between the individual distributions. A
crucial feature of the algorithm is that we allow for missing data, which may
be missing systematically. Specifically, individual distributions may be
missing all information in one or more coordinates.

Our original motivation to devise this type of algorithm stems from a
challenge faced by regulators of the financial industry. Data that are
delivered from financial institutions to the regulator consist of various
different formats, from highly aggregate data such as the total volume of the
balance sheet, down to very granular data about individual credits described
by characteristics such as volume, interest rate, etc. Two individual credits
might then be considered as similar if those characteristics take similar
numerical values, i.e. have small Euclidean distance when viewed as vectors in
$\mathbb{R}^{d}$. To build a distance between financial institutions, one can
view these institutions as probability distributions over a high dimensional
space, e.g. by considering the empirical measure of reported credits, and
determine the respective distance as a Wasserstein distance. One could then
interpret the landscape of all financial institutions as an ensemble of points
in the Wasserstein space, susceptible to familiar methods of clustering,
outlier detection, etc. A particular challenge which renders the problem more
complicated is the “missingness of data”. Data delivered by institutions may
have missing values. More crucially, data may be missing systematically, as
different institutions are required to deliver different data at varying
levels of detail and granularity.

The algorithm we propose as a ramification simultaneously clusters probability
measures with missing data and represents them as elements of a metric space.
The principal structure follows the idea of Lloyd’s algorithm for $k$-means
clustering and combines it with the concept of generalized Wasserstein
barycenters, recently introduced by Delon, Gozlan and Saint-Dizier [7]. A
classical approach to dealing with missing data would be to impute from
(weighted) nearest neighbors or a probability distribution in the same
cluster. However these types of imputation systematically skew results in the
sense that probability distributions with less reported data can appear as
being closer to other points which the imputation procedure is attempting to
mimic. This type of bias may be highly undesirable, e.g. when one is trying to
identify distributions that are “atypical”, where “atypical” could
specifically refer to the manner in which data are missing. We devise a
particular way of soft imputation which accounts for a random element in
filling up missing values, and in particular avoids the above mentioned bias.

For technical reasons we will present our algorithm in a more general form in
Section 3. That is, we cluster and perform soft imputation for points in an
abstract metric space rather than a Wasserstein space of probability measures.
This allows us to simplify the presentation and has the additional benefit
that we obtain a version of classical Euclidean $k$-means clustering with
missing data. This general perspective allows us in Section 7.1 to compare our
method with existing solutions for clustering in $\mathbb{R}^{d}$ with missing
data. As for the important case of probability distributions with missing data
our contribution is detailed in Section 5 and complemented experimentally in
Section 7.2 with results using actual (anonymized) data stemming from the
original motivation of regulation of financial institutions.

#### Related Literature

Clustering distributions using Wasserstein distances has been explored in
various works including [17, 25, 23, 12, 10]. A different financial context
was also considered in this area, including [11] by Horvath et al., for
clustering of market regimes. In [21], Staib et al. proposed the Wasserstein
$k$-means++ and its initialization strategy, generalizing the classical
$k$-means++ algorithm (introduced in [2]). However, to the best of our
knowledge, clustering distributions with missing coordinates has not been
investigated before. To address this problem, we rely heavily on the concept
of the generalized Wasserstein barycenter, introduced by Delon, Gozlan, and
Saint-Dizier [7], which extends the classical Wasserstein barycenter of Agueh
and Carlier [1]. For computational aspects related to optimal transport and
regularized Wasserstein distance, we refer to the work of Cuturi [5] and the
book [18] by Cuturi and Peyré. Additionally, we utilize the POT (0.8.2)
package (a Python package for optimal transport), as detailed in [8],
extensively for implementation purposes.

##  2 $k$-means clustering in metric spaces: a summary

Let $x_{1},\ldots,x_{N}$ be given points in a fixed metric space
$(\mathcal{X},d)$. The metric $k$-means problem is concerned with assigning
the data points into $k$ clusters. The clusters hereby are governed by $k$
barycenters, which will be some type of average or barycenter of the points in
the corresponding cluster. $k$-means was first introduced by MacQueen in [16].
The problem can be formalized as

$\displaystyle\textstyle\min\big{\\{}\sum_{i=1}^{N}d(x_{i},c_{a_{i}})^{2}:c_{1},\ldots,c_{k}\in\mathcal{X},~{}a\in\\{1,\ldots,k\\}^{N}\big{\\}}.$
(KM)

In this formulation $a$ denotes a vector of assignments, i.e. $a_{i}$
indicates the cluster membership of data point $x_{i}$. The points
$c_{1},\ldots,c_{k}$ serve as cluster barycenters. Furthermore a notion of
barycenter, being typically a function of some of the data points, is needed.
Then, problem (KM) can be tackled by the well-established Lloyd algorithm, see
also [14].

Following an initialization phase, two steps are iterated:

  1. i)

assignment step: given barycenters $c_{1},\ldots,c_{k}\in\mathcal{X}$, pick
for each $i\in\\{1,\ldots,N\\}$,

$\displaystyle\textstyle
a_{i}\in\operatorname*{arg\min}\\{d(x_{i},c_{j}):j\in\\{1,\ldots,k\\}\\},$ (1)

  2. ii)

barycenter step: given assignment $a$, update the barycenters, i.e. for each
$j\in\\{1,\ldots,k\\}$ set

$\displaystyle\textstyle c_{j}\in\text{barycenter}(\\{x_{i}:a_{i}=j\\}).$ (2)

Let us give two examples with particular choices of data space $\mathcal{X}$,
metric $d$, and barycenter operation, which will be of interest in the sequel:

######  Example 2.1.

Let $\mathcal{X}=\mathbb{R}^{d}$ and $d(x,y)=\lVert x-y\rVert_{2}$, the
Euclidean space and distance. Then (KM) simplifies to the classical $k$-means
problem in Euclidean space,

$\displaystyle\textstyle\min\big{\\{}\sum_{i=1}^{N}\lVert
x_{i}-c_{a_{i}}\rVert_{2}^{2}:c_{1},\ldots,c_{k}\in\mathbb{R}^{d},~{}a\in\\{1,\ldots,k\\}^{N}\big{\\}}.$
($k$-means)

The barycenter function is the Euclidean mean/average in this case.

######  Example 2.2.

Take $\mathcal{X}=\mathcal{P}_{2}(\mathbb{R}^{d})$, the space of all Borel
probability measures on $\mathbb{R}^{d}$ with finite second moments, and take
for $d$ the Wasserstein distance
$d(\mu,\nu):=W_{2}(\mu,\nu),~{}\mu,\nu\in\mathcal{P}_{2}(\mathbb{R}^{d})$
(recalled in (10) below). We will henceforth call
$\mu_{i}=x_{i},~{}i=1,\ldots,N$. Then (KM) turns into

$\displaystyle\textstyle\min\big{\\{}\sum_{i=1}^{N}W_{2}^{2}(\mu_{i},\nu_{a_{i}}):\nu_{1},\ldots,\nu_{k}\in\mathcal{P}_{2}(\mathbb{R}^{d}),~{}a\in\\{1,\ldots,k\\}^{N}\big{\\}}.$
(W-$k$-means)

As barycenter one takes the Wasserstein barycenter, i.e. the barycenter
updating step (2) reads as

$\displaystyle\textstyle\nu_{j}\in\operatorname*{arg\min}\big{\\{}\sum_{i:a_{i}=j}W_{2}^{2}(\mu_{i},\nu):\nu\in\mathcal{P}_{2}(\mathbb{R}^{d})\big{\\}}.$
(3)

In what follows we want to generalize the previous examples, especially
Example 2.2, to the case in which certain coordinates/marginals are missing
for certain observations, for an illustration see Figure 1 in Section 5.

For Lloyd’s algorithm, both a distance function and a notion of barycenter are
needed. In the case of a metric space one could always consider a Fréchet mean
for the barycenter operator, see [9]. In the case of missing coordinates,
however, the data does not come from one metric space but from several
different ones, therefore several distances will be considered. For the same
reason, in order to define a proper notion of barycenter, some generalization
will be needed. We will introduce a general approach to this in the next
section and discuss our specific examples in more detail in Example 3.1 and
Section 5.

##  3 Clustering projected elements of metric spaces

Suppose we are working in a metric space $(\mathcal{X},d)$ and that there are
points $x_{1},\ldots,x_{N}\in\mathcal{X}$ which we want to assign to $k$
clusters. However, what we really observe are the points
$\tilde{x}_{i}:=\varphi_{i}(x_{i})$, where
$\varphi_{i}:\mathcal{X}\to\mathcal{X}_{i}$ is a known map into another metric
space $\mathcal{X}_{i}$. One can think for example of a projection onto some
coordinates in the case of $\mathcal{X}=\mathbb{R}^{d}$. Having observed
$\tilde{x}_{i}\in\mathcal{X}_{i}$ define for $y\in\mathcal{X}$

$\displaystyle\textstyle
d_{i}(\tilde{x}_{i},y):=d(\varphi_{i}^{-1}(\tilde{x}_{i}),y):=\inf\\{d(x,y):x\in\mathcal{X},~{}\varphi_{i}(x)=\tilde{x}_{i}\\},$
(4)

which serves as a type of dissimilarity measure of a point $y\in\mathcal{X}$
in the “full” metric space.

###  3.1 Problem and algorithm

Having this at hand, we introduce the problem

$\displaystyle\textstyle\min\big{\\{}\sum_{i=1}^{N}d_{i}(\tilde{x}_{i},c_{a_{i}})^{2}:c_{1},\ldots,c_{k}\in\mathcal{X},~{}a\in\\{1,\ldots,k\\}^{N}\big{\\}}.$
(5)

This says we want to find optimal cluster barycenters in $\mathcal{X}$, as
well as a vector optimally assigning each observed point to a cluster. It is
important to note that we are always looking for cluster barycenters in the
“full” space $\mathcal{X}$.

We suggest to tackle (5) by the following two steps which are iterated in a
Lloyd algorithm fashion:

  1. i)

given $c_{1},\ldots,c_{k}$, set
$a_{i}\in\operatorname*{arg\min}\\{d_{i}(\tilde{x}_{i},c_{j}):j\in\\{1,\ldots,k\\}\\},$

  2. ii)

given $a$, put
$c_{j}\in\operatorname*{arg\min}\\{\sum_{i:a_{i}=j}d_{i}(\tilde{x}_{i},y)^{2}:y\in\mathcal{X}\\}.$

Let us note that step ii does not necessarily admit minimizers. However, in
our two applications, i.e. NA $k$-means in Example 3.1 and NA Wasserstein
$k$-means, which is discussed in Section 5, minimizers always exist. For a
more precise discussion of this technical point see Appendix B.

###  3.2 Using clustering for imputation

After clustering we obtain a vector of assignments $a$ and barycenters
$c_{1},\ldots,c_{k}$. Define

$\displaystyle\textstyle
I_{f}:=\\{i:\varphi_{i}=Id_{\mathcal{X}}\\},~{}~{}~{}I_{m}:=\\{1,\ldots,N\\}\setminus
I_{f},$ (6)

i.e. respectively the set of indices being observed in $\mathcal{X}$ (i.e.
fully observed data) and the set of indices which are only observed after some
non-trivial map (i.e. with missing data). Let us also define the final
clusters, i.e. for each $j\in\\{1,\ldots,k\\}$ define

$\displaystyle\textstyle C_{j}:=\\{i:a_{i}=j\\}.$

We want to use the clusters to find, for a point $\tilde{x}_{i}$ with $i\in
I_{m}$, a probability measure on $\mathcal{X}$, i.e. an element of
$\mathcal{P}(\mathcal{X})$, which hopefully concentrates around the true
unobserved point $x_{i}$. We define for $i\in I_{m}$ the set of indices we use
for filling up the missing data of $\tilde{x}_{i}$, as

$\displaystyle\textstyle J_{i}:=I_{f}\cap C_{a_{i}}.$

This says that we want to use the “full” points in the same cluster as
$\tilde{x}_{i}$ in order to compensate the incomplete information that we have
about $\tilde{x}_{i}$. It can of course happen that $J_{i}=\emptyset$. In this
case we use the corresponding cluster barycenter $c_{a_{i}}$ to complement the
information about $\tilde{x}_{i}$. In the following we assume
$J_{i}\neq\emptyset$. Taking now some $x_{\ell}$ with $\ell\in J_{i}$ we set

$\displaystyle\textstyle
y^{i}_{\ell}\in\operatorname*{arg\min}\\{d(x,x_{\ell}):x\in\mathcal{X},~{}\varphi_{i}(x)=\tilde{x}_{i}\\},$
(7)

which will be one choice of “filling up”. We introduce the shorthand
$D^{i}_{\ell}:=d(y^{i}_{\ell},x_{\ell})=d(\varphi_{i}^{-1}(\tilde{x}_{i}),x_{\ell})$.
In order to determine the weights of a probability measure, we set
$p^{i}_{\ell}:=f(D^{i}_{\ell})$ with $f:[0,\infty)\to[0,\infty)$ being some
positive decreasing function fixed in advance, e.g. $f(x):=\exp(-x^{2})$. This
way $p^{i}_{\ell}\geq 0$ and also $\sum_{\ell\in J_{i}}p^{i}_{\ell}=1$ after
possibly renormalizing the weights by a positive constant. Having obtained the
weights, we define the probability measure which should represent a _randomly
reconstructed_ $x_{i}$ as

$\displaystyle\textstyle\theta_{i}:=\sum_{\ell\in
J_{i}}p^{i}_{\ell}\delta_{y^{i}_{\ell}}.$

To embed also the fully observed points, for $i\in I_{f}$ we set
$\theta_{i}:=\delta_{x_{i}}$, i.e. the Dirac delta concentrated on $x_{i}$.
For $i\in I_{m}$ with $J_{i}=\emptyset$, we set $\theta_{i}$ to be the Dirac
delta concentrated on some minimizer of
$d(\varphi_{i}^{-1}(\tilde{x}_{i}),c_{a_{i}})$. Thus we embed all observed
points $\tilde{x}_{1},\ldots,\tilde{x}_{N}$ in the same space, i.e.
$\mathcal{P}(\mathcal{X})$. In other words, we identify the possibly
unobserved point $x_{i}$ with a probability measure $\theta_{i}$.

In order to compare the hitherto constructed probability measures, as we will
need to do in Section 7, we define a (generalized) metric on
$\mathcal{P}(\mathcal{X})$, see Lemma A.1, through the cost induced by the
_product_ or _independent coupling_. That is, for
$\mu,\nu\in\mathcal{P}(\mathcal{X})$ we set $\rho(\mu,\nu):=0$ if $\mu=\nu$
and otherwise

$\displaystyle\textstyle\rho(\mu,\nu):=\int_{\mathcal{X}\times\mathcal{X}}d(x,x^{\prime})\mathrm{d}\mu\otimes\nu(x,x^{\prime}):=\int_{\mathcal{X}}\int_{\mathcal{X}}d(x,x^{\prime})\mathrm{d}\mu(x)d\nu(x^{\prime}).$
(8)

With a Wasserstein distance, or a similar notion of distance, imputed points
are biased to be closer than “fully” observed points. The metric $\rho$ is
designed to avoid this type of bias. Indeed, the choice in (8) formalises the
idea that missing values are not imputed in a deterministic sense but in a
_random_ or _soft_ fashion, as indicated in the introduction. The distance
between two randomly imputed values is then estimated as an _independent_
average of distances, following the rationale that the imputation for one
point does not inform the imputation for a different point.

######  Example 3.1 (NA $k$-means).

Corresponding to classical $k$-means, i.e. Example 2.1, consider
$\mathcal{X}=\mathbb{R}^{d}$. Let
$\varphi_{i}=P_{i}:\mathbb{R}^{d}\to\mathbb{R}^{d_{i}}$ be projections to some
of the coordinates. Then (5) becomes

$\displaystyle\textstyle\min\big{\\{}\sum_{i=1}^{N}\lVert
P_{i}(x_{i})-P_{i}(c_{a_{i}})\rVert_{2}^{2}:c_{1},\ldots,c_{k}\in\mathbb{R}^{d},~{}a\in\\{1,\ldots,k\\}^{N}\big{\\}}.$
(9)

The cluster assignment iteration steps assigns each point to the cluster to
whose barycenter it is closest to, by only considering the known coordinates.
The barycenter updating step is solved by
$\big{(}\sum_{i:a_{i}=j}P_{i}^{T}P_{i}\big{)}^{-1}\sum_{i:a_{i}=j}P_{i}^{T}P_{i}(x_{i}),$
which in each coordinate is the average of all the points in the cluster
having that coordinate. See Section 5 for details if the inverse does not
exit. By imputing according to the above discussion, we obtain measures in
$\mathcal{P}_{2}(\mathbb{R}^{d})$, on which we can then calculate pairwise
distances using the metric $\rho$ in (8). Experiments using this method may be
found in Section 7.1. It is worth mentioning that Chi et. al. in [4] consider
the same loss function, i.e. (9), but they propose a different algorithm for
clustering Euclidean points with missing values.

Before discussing our second example in Section 5, corresponding to the
generalization of Example 2.2, we recall the needed notions from optimal
transport theory.

##  4 Wasserstein distance and generalized Wasserstein barycenter: a summary

#### Optimal transport and Wasserstein distance

Let $\mu,\nu$ be probabilities on Polish spaces $X,Y$ resp. Denote by

$\displaystyle\textstyle\Pi(\mu,\nu):=\\{\pi\in\mathcal{P}(X\times
Y):pr_{X\\#}\pi=\mu,~{}pr_{Y\\#}\pi=\nu\\}$

the set of all measures on the product space having $\mu$ and $\nu$ as
marginals. The Kantorovich problem for a cost function $c:X\times
Y\to\mathbb{R}_{+}$, introduced in [13], is

$\displaystyle\textstyle\inf_{\pi\in\Pi(\mu,\nu)}\int_{X\times Y}c(x,y)\
\mathrm{d}\pi(x,y).$ (KP)

We specialise to the case $X=Y=\mathbb{R}^{d}$ and $c(x,y)=\lVert
x-y\rVert_{p}^{p}$ for $p\in[1,\infty)$. The Wasserstein distance on the space
$\mathcal{P}_{p}(\mathbb{R}^{d})$ of probability measures with finite
$p$-moment is then defined as

$\displaystyle\textstyle
W_{p}(\mu,\nu):=\big{(}\inf_{\pi\in\Pi(\mu,\nu)}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\lVert
x-y\rVert_{p}^{p}\ \mathrm{d}\pi(x,y)\big{)}^{1/p}.$ (10)

#### Wasserstein barycenter

A notion of averaging probabilities which has recently received significant
attention is the notion of Wasserstein barycenters, introduced by Agueh and
Carlier [1]. A Wasserstein barycenter of
$\mu_{1},\ldots,\mu_{n}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ with weights
$\lambda_{1},\ldots\lambda_{n}\geq 0$, is a solution of

$\displaystyle\textstyle\inf\\{\sum_{i=1}^{n}\lambda_{i}W_{2}^{2}(\mu_{i},\nu):\nu\in\mathcal{P}_{2}(\mathbb{R}^{d})\\}.$

Since we want to generalize the algorithm introduced in Sections 2-3 to the
setting of distributions we need to have a variant of the Wasserstein
barycenter, which is applicable if only some coordinates of the measures are
known. An appropriate generalized Wasserstein barycenter was introduced and
studied by Julien, Gozlan and Saint-Dizier in [7]. To formulate it, let
probability measures
$\mu_{1},\ldots,\mu_{n}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ be given and
consider linear maps $P_{i}:\mathbb{R}^{d}\to\mathbb{R}^{d_{i}}$. In our
intended applications, these will be projections onto some of the coordinates.
A generalized Wasserstein barycenter of the push-forwarded measures
$P_{1}\\#\mu_{1},\ldots,P_{n}\\#\mu_{n}$ with weights
$\lambda_{1},\ldots,\lambda_{n}\geq 0$, is a solution of

$\displaystyle\textstyle\inf\\{\sum_{i=1}^{n}\lambda_{i}W_{2}^{2}(P_{i}\\#\mu_{i},P_{i}\\#\nu):\nu\in\mathcal{P}_{2}(\mathbb{R}^{d})\\}.$
(11)

To solve the problem it is useful to translate it back to a classical
Wasserstein barycenter problem. By associating to $P_{i}$ a matrix
$P_{i}\in\mathbb{R}^{d_{i}\times d}$, define
$A:=\sum_{i=1}^{n}\lambda_{i}P_{i}^{T}P_{i}$ and assume in the following that
$A$ is invertible. Then set
$\tilde{\mu_{i}}:=(A^{-1/2}P_{i}^{T})\\#\mu_{i}\in\mathcal{P}_{2}(\mathbb{R}^{d})$
for $i=1,\ldots,n$, and consider the conventional Wasserstein barycenter
problem

$\displaystyle\textstyle\inf\Big{\\{}\sum_{i=1}^{n}\lambda_{i}W_{2}^{2}(\tilde{\mu_{i}},\tilde{\nu}):\tilde{\nu}\in\mathcal{P}_{2}(\mathbb{R}^{d})\Big{\\}}.$
(12)

Then, $\nu$ is a solution to (11) if and only if $A^{1/2}\\#\nu$ is a solution
to (12), thanks to [7, Proposition 3.1]. Concerning computational aspects,
especially in the discrete case, efficient algorithms for calculating
Wasserstein barycenters have already been developed; see the work [6] by
Cuturi and Doucet, as well as [15] by Álvarez, del Barrio, Cuesta-Albertos and
Matrán, to name but two.

##  5 NA Wasserstein $k$-means

#### Algorithm description

Having the two operations needed in Section 3 at hand, we can discuss the case
of $\mathcal{X}=\mathcal{P}_{2}(\mathbb{R}^{d})$ in detail. Suppose we want to
cluster probability measures
$\mu_{1},\ldots,\mu_{N}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ into $k$ clusters
but we only observe the push forwards
$P_{i}\\#\mu_{i}\in\mathcal{P}_{2}(\mathbb{R}^{d_{i}})$ under projections
$P_{i},~{}i=1,\ldots,N$. Hence we have $\varphi_{i}(\cdot):=P_{i}\\#(\cdot)$.
Then problem (5) reads as

$\displaystyle\textstyle\inf\\{\sum_{i=1}^{N}W_{2}^{2}(P_{i}\\#\mu_{i},P_{i}\\#\nu_{a_{i}}):\nu_{1},\ldots,\nu_{k}\in\mathcal{P}_{2}(\mathbb{R}^{d}),~{}a\in\\{1,\ldots,k\\}^{N}\\}.$
(NA-W-$k$-means)

We can tackle problem (NA-W-$k$-means) by the suggested iterations in Section
3.1, which read as

  1. i)

update cluster assignments given barycenters $\nu_{1},\ldots,\nu_{k}$, i.e.
for each $i\in\\{1,\ldots,N\\}$ set

$\displaystyle\textstyle
a_{i}\in\operatorname*{arg\min}\\{W_{2}(P_{i\\#}\mu_{i},P_{i\\#}\nu_{j}):j\in\\{1,\ldots,k\\}\\},$
(13)

  2. ii)

update cluster barycenters given an assignment $a$, i.e. for each
$j\in\\{1,\ldots,k\\}$ set

$\displaystyle\textstyle\nu_{j}\in\operatorname*{arg\min}\big{\\{}\sum_{i:a_{i}=j}W_{2}^{2}(P_{i\\#}\mu_{i},P_{i\\#}\nu):\nu\in\mathcal{P}_{2}(\mathbb{R}^{d})\big{\\}}.$

Again, we want to find a vector $a$ of optimally assigning each measure to one
of the $k$ clusters, as well as optimal barycenters, or in this case
generalized Wasserstein barycenters, $\nu_{1},\ldots,\nu_{k}$. Note that once
more we are looking for barycenters on the full space, i.e. having all
coordinates.

It can happen that a cluster ends up consisting of measures which all miss the
same coordinate, i.e. for some $j\in\\{1,\ldots,k\\}$ we have
$\bigcap_{i:a_{i}=j}\ker P_{i}\neq\\{0\\}.$ In order not to run into numerical
problems in that case, one can adapt step ii: Assume at iteration $t$ the
barycenters are $\nu^{(t)}_{1},\ldots,\nu^{(t)}_{k}$ and a new assignment
vector $a^{(t+1)}$ has been calculated. Introduce a weight
$\lambda^{(t)}\in[0,1)$ and set

$\displaystyle\textstyle\nu_{j}^{(t+1)}\in\operatorname*{arg\min}\big{\\{}(1-\lambda^{(t)})\sum_{i:a_{i}^{(t+1)}=j}W_{2}^{2}(P_{i}\\#\mu_{i},P_{i}\\#\nu)+\lambda^{(t)}W_{2}^{2}(\nu_{j}^{(t)},\nu):\nu\in\mathcal{P}_{2}(\mathbb{R}^{d})\big{\\}}.$

Experimentally we have observed that $\lambda^{(t)}:=\frac{1}{(t+1)^{1/2}}$
works well. Having defined the two steps to be iterated, we can present the
main computational contribution of this work, Algorithm 1:

Data: $N$ observed measures $P_{1}\\#\mu_{1},\ldots,P_{N}\\#\mu_{N}$, number
of cluster $k$, maximum number of iterations $T$, weighting schedule
$(\lambda^{(t)})_{t=0}^{T}$

Result: barycenters $\nu_{1},\ldots,\nu_{k}$, assignment vector
$a\in\\{1,\ldots,K\\}^{N}$

$t:=0$, initialize $\nu_{1}^{(0)},\ldots,\nu_{k}^{(0)}$, and $a^{(0)}$

while  _$t <T$_ do

for  _$i=1,\ldots,N$_ do

$\displaystyle\textstyle
a_{i}^{(t+1)}:=\operatorname*{arg\min}\\{W_{2}^{2}(P_{i}\\#\mu_{i},P_{i}\\#\nu_{j}^{(t)}):j\in\\{1,\ldots,k\\}\\}$
(14)

for  _$j=1,\ldots,k$_ do

set $\nu_{j}^{(t+1)}$ as solution of
$\displaystyle\textstyle\min\big{\\{}(1-\lambda^{(t)})\sum_{i:a_{i}^{(t+1)}=j}W_{2}^{2}(P_{i}\\#\mu_{i},P_{i}\\#\nu)+\lambda^{(t)}W_{2}^{2}(\nu_{j}^{(t)},\nu):\nu\in\mathcal{P}_{2}(\mathbb{R}^{d})\big{\\}}$

if  _$a^{(t)}==a^{(t+1)}$_ then

break

$t:=t+1$

Algorithm 1 NA Wasserstein $k$-means algorithm

In (14) we use the entry of the old assignment vector $a_{i}^{(t)}$ if it is a
minimizer. We can formulate a simple convergence result whose proof is
reported in Appendix A. Let

$\displaystyle\textstyle
L(\nu_{1},\ldots,\nu_{K},a):=\sum_{i=1}^{N}W_{2}^{2}(P_{i}\\#\mu_{i},P_{i}\\#\nu_{a_{i}}).$

######  Theorem 5.1.

Given initialization
$\nu_{1}^{(0)},\ldots,\nu_{k}^{(0)}\in\mathcal{P}_{2}(\mathbb{R}^{d}),~{}a\in\\{1,\ldots,k\\}^{N}$,
Algorithm 1 strictly decreases $L$ until it terminates after finitely many
steps, for any choice of $(\lambda^{(t)})_{t\in\mathbb{N}}$ (even for
$T=\infty$). If in every iteration each cluster has a measure with all
coordinates, then Algorithm 1 with $\lambda^{(t)}=0$ yields a local minimum of
$L$.

#### A toy example

To illustrate the algorithm, consider six measures on the plane, of which one
misses the vertical coordinate (denoted by vertical lines). Each measure
consists of three support points with equal weight. The measures are
visualized by different colors and depicted in the left of Figure 1. If we
think of clustering these measures into three clusters, it is quite clear
which clusters we want. By applying the introduced algorithm this is indeed
the case and leads to the results depicted in Figure 1. To obtain it we
calculated free support barycenters in each step of the algorithm with a fixed
support size of three. For the blue and brown barycenters each support point
is precisely the average of support points of the two measures in that
clusters. For the pink barycenter the situation is slightly different. In the
$y$ coordinate it inherits the values of the red measure, while as in the $x$
coordinate the support points are averages of the $x$ coordinates of a red and
rose support point, respectively.

Figure 1: (left): 6 measures on the plane of which one misses the vertical
coordinate of its points. (right) Clustering of the 6 measures into 3
clusters.

#### Imputation via NA Wasserstein $k$-means

Let us look how the imputation procedure according to Section 3.2 works in the
case where $\mathcal{X}=\mathcal{P}_{2}(\mathbb{R}^{d})$. Having clustered
measures$P_{i}\\#\mu_{i}\in\mathcal{P}_{2}(\mathbb{R}^{d_{i}}),~{}i=1,\ldots,N$,
we obtain barycenters
$\nu_{1},\ldots,\nu_{k}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ and a vector of
assignments $a\in\\{1,\ldots,k\\}^{N}$. For $i\in I_{m}$, see (6) for the
definition, we then set

$\displaystyle\textstyle\boldsymbol{\mu_{i}}:=\sum_{\ell\in
J_{i}}p_{\ell}^{i}\boldsymbol{\delta_{\tau_{i}\\#(P_{i}\\#\mu_{i}\otimes
P_{i}^{C}\\#\mu_{\ell})}}$

to obtain a measure
$\boldsymbol{\mu_{i}}\in\mathcal{P}_{2}(\mathcal{P}_{2}(\mathbb{R}^{d}))$,
i.e. a random measure. Here, by $P_{i}^{C}$ we denote the projection on the
missing coordinates of $P_{i}\\#\mu_{i}$, and by $\tau_{i}$ we denote the
unique permutation such that
$\tau_{i}\circ(P_{i},P_{i}^{C})=Id_{\mathbb{R}^{d}}$. For $i\in I_{f}$ we set
$\boldsymbol{\mu_{i}}:={\delta_{\mu_{i}}}$. For $i\in I_{m}$ and $\lvert
J_{i}\rvert>1$ we suggest the weights

$\displaystyle\textstyle
p_{\ell}^{i}\propto\exp\big{(}-\frac{\lambda}{2\sigma_{i}^{2}}W_{2}^{2}(P_{i}\\#\mu_{i},P_{i}\\#\mu_{\ell})\big{)},~{}~{}~{}\ell\in
J_{i}.$

Here, $\lambda>0$ and the variance $\sigma_{i}^{2}$ is used to standardize the
distances, i.e.

$\displaystyle\textstyle\sigma_{i}^{2}:=\frac{1}{\lvert
J_{i}\rvert-1}\sum_{\ell\in
J_{i}}W_{2}^{2}(P_{i}\\#\mu_{i},P_{i}\\#\mu_{\ell}).$

If $\lvert J_{i}\rvert=1$ it is clear what to do, as we only have one measure
to fill up with.

After imputation we obtain random measures
$\boldsymbol{\mu_{1}},\ldots,\boldsymbol{\mu_{N}}\in\mathcal{P}_{2}(\mathcal{P}_{2}(\mathbb{R}^{d}))$.
To calculate their distances, we use metric $\rho$ from (8), i.e., for
$\boldsymbol{\mu},\boldsymbol{\nu}\in\mathcal{P}(\mathcal{P}_{2}(\mathbb{R}^{d}))$,
$\rho(\boldsymbol{\mu},\boldsymbol{\nu})=0$, if
$\boldsymbol{\mu}=\boldsymbol{\nu}$, and

$\displaystyle\textstyle\rho(\boldsymbol{\mu},\boldsymbol{\nu})=\int_{\mathcal{P}_{2}(\mathbb{R}^{d})}\int_{\mathcal{P}_{2}(\mathbb{R}^{d})}W_{2}(\mu,\nu)\mathrm{d}\boldsymbol{\mu}(\mu)\mathrm{d}\boldsymbol{\nu}(\nu),~{}if~{}\boldsymbol{\mu}\neq\boldsymbol{\nu}.$
(15)

#### Limitations

For the dataset of our application in Section 5 the method developed above is
sufficient. However, for other applications it might be desirable to
investigate the use of the _regularized Wasserstein distance_ or _entropic
optimal transport_ in our method for reduced computational complexity. For
this, a notion of a _generalized regularized Wasserstein barycenter_ , as the
NA-Wassertein $k$-means algorithm relies heavily on the generalized
Wasserstein barycenter has to be developed which we intend to do in future
work.

##  6 Evaluation method

Our clustering/imputation method will be evaluated on Euclidean and
distributional data. Hence we present here in abstract language the setting in
which the method will be tested. Suppose we have elements $x_{1},\ldots,x_{N}$
in a metric space $(\mathcal{X},d_{\mathcal{X}})$, which we collect in the set
$X:=\\{x_{1},\ldots,x_{N}\\}$. We only observe data via the maps
$\varphi_{i}:\mathcal{X}\to\mathcal{X}_{i}$ from Section 3, i.e.
$\tilde{x}_{i}=\varphi_{i}(x_{i})$. We collect the observed points in
$(\tilde{x}_{1},\ldots,\tilde{x}_{N})\in\mathcal{X}_{NA}$, where we define
$\mathcal{X}_{NA}:=\mathcal{X}_{1}\times\cdots\times\mathcal{X}_{N}.$ Let us
summarize this by saying that we only observe the values of
$(x_{1},\dots,x_{N})$ under the map $h:\mathcal{X}^{N}\to\mathcal{X}_{NA}$,
i.e. $h_{i}(x_{1},\ldots,x_{N})=\varphi_{i}(x_{i})=\tilde{x}_{i}$.

The goal is to reconstruct the metric structure of $X$ as accurately as
possible. This means that we want to find a _good_ metric space
$(\mathcal{Y},d_{\mathcal{Y}})$ and a reconstruction map
$R:\mathcal{X}_{NA}\to\mathcal{Y}^{N}$.

Before finding a way of comparing the reconstructed points
$R((\tilde{x}_{1},\ldots,\tilde{x}_{N}))$ to the original points $X$, we also
want to consider another piece of information, namely that some observations
of $X$ might be more important for us than others. This is motivated by our
main practical application of reconstructing a banking landscape, where we
think of the importance of a bank being linked to attributes such as the size
of that bank as measured, e.g. by its total assets. It is then more important
to reconstruct the characteristics of a very large bank than those of a much
smaller bank. We model this by assigning weights proportional to the size of
the bank. Formally, we define a probability measure $\mu_{X}$ on $X$ by
setting $\mu_{X}(x_{i}):=p_{i}\geq 0$, where we assume
$\sum_{i=1}^{N}p_{i}=1$. Naturally we define $\mu_{Y}$, a probability measure
on $\mathcal{Y}$, via
$\mu_{Y}(R_{i}(\tilde{x}_{1},\ldots,\tilde{x}_{N})):=p_{i}$.

To compare the metric measure spaces $(X,d_{\mathcal{X}},\mu_{X})$ and
$(R((\tilde{x}_{1},\ldots,\tilde{x}_{N})),d_{\mathcal{Y}},\mu_{Y})$, we use a
_Gromov-Wasserstein distance_ (see [22, 8] for resp. theory and
implementation):

$\displaystyle\textstyle
GW_{2}((\mathcal{X},d_{\mathcal{X}},\mu),(\mathcal{Y},d_{\mathcal{Y}},\nu))^{2}:=\inf_{\pi\in\Pi(\mu,\nu)}\int_{\mathcal{X}^{2}\times\mathcal{Y}^{2}}\lvert
d_{\mathcal{X}}(x,x^{\prime})-d_{\mathcal{Y}}(y,y^{\prime})\rvert^{2}\mathrm{d}\pi(x,y)\mathrm{d}\pi(x^{\prime},y^{\prime}).$

##  7 Experimental results

###  7.1 Reconstructing points from a Gaussian Mixture Model

We start by comparing our clustering/reconstruction method to existing
imputation methods in the Euclidean case, i.e. $\mathcal{X}=\mathbb{R}^{d}$,
and evaluate the results using the Gromov-Wasserstein distance.

We simulate data from a Gaussian Mixture Model
$\gamma:=\sum_{j=1}^{k}\alpha_{j}\mathcal{N}(\mu_{j},\Sigma_{j})$, with
$\sum_{j=1}^{k}\alpha_{j}=1$. To simulate the importance of points, i.e.
$\mu_{X}(x_{i})$ from Section 6, we sample from
$\text{Lognormal}(\mu,\sigma^{2})$ to weigh each observation proportional to
its sample value. In the simulation study we set $k=5$, $d=10$, and we always
sample $N=500$ points. To simulate missing values we use several missingness
structures, i.e. several choices for the map $h$ from Section 6. Before
creating missing values, we adjust each dimension by its mean and standard
deviation. In order to define $h$ from Section 6, we calculate for each
dimension the $25\%$ quantile, draw a uniform random variable $h_{j}\sim
U(0,2/k)$ for each cluster $j$ and choose a correlation parameter
$\beta\in[0,1]$. We then define $h$ depending on $\beta$ by the following
properties: for each cluster $j\in\\{1,\ldots,k\\}$ set $100\min\\{(\beta
h_{j}+(1-\beta)\alpha_{j}),0.95\\}\%$ to NA, and for each dimension set the
values below the previously calculated $25\%$ quantile to NA.

This means we are left with choosing the parameters
$\alpha_{1},\ldots,\alpha_{5},\mu_{1},\ldots,\mu_{5},\Sigma_{1},\ldots,\Sigma_{5}$
for the Gaussian Mixture model, $\beta$ for the map $h$, and $\mu,\sigma^{2}$
for the lognormal which will govern the importance of sampled points. We set
$\mu=20,\sigma=1.5$ as with these parameters the lognormal distribution models
total assets of a financial institution reasonably well (empirically). We let
$\beta$ vary between $0$ and $1$ with steps $0.1$. For the Gaussian mixture
model we want to simulate the weights of the normal distributions, the means
and the covariance, as non-informative as possible. Therefore, the
$\alpha_{j}$’s are uniformly sampled from the probability simplex. The
$\mu_{j}$’s are sampled uniformly on the cube $[-10,10]^{d}$ and the
$\Sigma_{j}$’s are sampled from a Wishart distribution with parameters
$(d,Id_{d})$.

For each choice of $\beta$ we sample the explained parameters $100$ times,
create data, standardize it and calculate pairwise distances between data.
Then, we create missing values according to $h$ and compare different
imputation procedures presented below to reconstruct the points and then
calculate pairwise distances between the imputed points in the case if missing
values are imputed by points. In the case of imputing the points by measures,
we use the distance $\rho$ from Section 3.2 to calculate pairwise distances.
Afterwards we calculate the corresponding Gromov Wasserstein distance between
original pairwise distances and imputed pairwise distances with the
corresponding weights which come from the $Lognormal(\mu,\sigma^{2})$. We look
at the following imputation techniques:

  * •

NA $k$-means: this is our method coming from Example 3.1,

  * •

NA $k$-means-m: this corresponds to using NA $k$-means for clustering and
imputation, but after obtaining the measures $\theta_{i}$ we calculate their
expectation and use this as imputed point,

  * •

mean imputation: impute each missing value by the average value of this
attribute,

  * •

median imputation: impute each missing value by the median value of this
attribute,

  * •

multiple imputation: a method [20, 3] trying out several values for each
missing value,

  * •

KNN: K-nearest-neighbour imputation which imputes each missing value by a
weighted average of points close in the coordinates which are not missing,
where we choose $K=4$,

  * •

LR: impute by regressing on points, with missing values and using predictions
as imputation.

In Table 1 we report the estimated Gromov Wasserstein Distances $\pm$ standard
errors based on sampling $100$ times for each scenario of missing values.

Table 1: Gromov Wasserstein Distance for Euclidean simulation.

settings | NA $k$-means. | NA $k$-means-m | mean imp. | median imp. | multiple imp. | KNN | LR
---|---|---|---|---|---|---|---
$~{}\pm$

We see that, in the case of Euclidean points with missing values, the
introduced method outperforms classical imputation methods when considering
Gromov-Wasserstein as evaluation measure. Indeed, whether we apply _NA
$k$-means_ directly, or perform another averaging step, i.e.  _NA $k$-means-
m_, these two columns seem to outperform the remaining five for all settings
of missing value creation. In Appendix C we use another way to evaluate the
method, namely the Rand index, when viewing NA $k$-means solely as clustering
algorithm and not as a method to reconstruct a metric structure. Also, further
settings of missing value generations can be found there to accompany our
results.

###  7.2 Reconstructing financial institutions

As mentioned in the introduction, our main motivation for developing the
proposed method comes from clustering financial institutions based on granular
loan data they are obliged to report to the central bank. Applying the
proposed algorithm to such real data in a practical empirical example, we
consider 500 financial institutions. For each of them we include information
on 100 loans in 7 dimensions, i.e. each loan is described by up to 7
attributes. Among the attributes considered, we include e.g. credit volume,
time to maturity, and interest rate. Out of the $500$ institutions, $3$ just
reported $2$ attributes, $9$ institutions reported $3$, $56$ reported $4$,
$110$ reported $5$, $134$ reported $6$ and $188$ reported all attributes. We
cluster the institutions into $15$ clusters using the NA Wasserstein $k$-means
algorithm proposed in Section 5. Afterwards we use the imputation method
described in the same section. Using the distance $\rho$ on
$\mathcal{P}_{2}(\mathcal{P}_{2}(\mathbb{R}^{7}))$, defined in (15), we
calculate pairwise distances between the financial institutions. In order to
get a visual representation of the landscape of financial institutions, and
the clustering results, we use _multidimensional scaling_ to represent the
matrix of pairwise distances $(\rho(\theta_{i},\theta_{j}))_{i,j=1}^{500}$ as
points in $\mathbb{R}^{3}$ with a similar matrix of pairwise distances.

Figure 2: Visualization of pairwise distances.

In Figure 2 we see the representation of a banking landscape and clusters.
Each color describes a cluster, and each point a bank. We observe that
clusters are not strictly separated, due to multidimensional scaling, since we
want to embed points via pairwise distances from a very high-dimensional into
a very low-dimensional space. The visualisation may help analysts to further
their understanding which financial institutions are similar to each other and
also identify financial institutions which behave significantly different to
others with respect to their loan structure. Detecting potential outliers is
of course of particular interest.

##  8 Discussion

We have proposed a novel method for clustering elements in a metric space that
are only observable under maps to potentially smaller metric spaces. The
method can also reconstruct the metric structure of the original elements. In
a simulation study in Euclidean space, our proposed method, boiling down to a
new NA $k$-means algorithm, outperforms classical imputation methods when
using the Gromov-Wasserstein distance as evaluation criterion. A particular
variant of our method, the NA Wasserstein $k$-means algorithm, clusters
probability measures that are only observable as push-forwards under
projections. This algorithm is novel, and to our knowledge, there are no
references in the literature dealing with partially observed distributions. In
upcoming work with a central bank we will study, much more extensively than in
Section 7 and from an enhanced economic point of view, the application to
financial data for regulatory purposes.

## References

  * [1] M. Agueh and G. Carlier.  Barycenters in the Wasserstein space.  SIAM Journal on Mathematical Analysis, 43(2):904–924, 2011.
  * [2] D. Arthur and S. Vassilvitskii.  K-means++: The advantages of careful seeding.  In Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA ’07, page 1027–1035, USA, 2007. Society for Industrial and Applied Mathematics.
  * [3] M. J. Azur, E. A. Stuart, C. Frangakis, and P. J. Leaf.  Multiple imputation by chained equations: what is it and how does it work?  International Journal of Methods in Psychiatric Research, 20(1):40–49, 2011.
  * [4] J. T. Chi, E. C. Chi, and R. G. Baraniuk.  $k$-pod: A method for $k$-means clustering of missing data.  The American Statistician, 70(1):91–99, jan 2016.
  * [5] M. Cuturi.  Sinkhorn distances: Lightspeed computation of optimal transport.  In C. J. C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Q. Weinberger, editors, Advances in Neural Information Processing Systems, volume 26. Curran Associates, Inc., 2013.
  * [6] M. Cuturi and A. Doucet.  Fast computation of wasserstein barycenters.  In E. P. Xing and T. Jebara, editors, Proceedings of the 31st International Conference on Machine Learning, volume 32 of Proceedings of Machine Learning Research, pages 685–693, Bejing, China, 22–24 Jun 2014\. PMLR.
  * [7] J. Delon, N. Gozlan, and A. Saint-Dizier.  Generalized Wasserstein barycenters between probability measures living on different subspaces.  To appear in Annals of Applied Probability, 2021.
  * [8] R. Flamary, N. Courty, A. Gramfort, M. Z. Alaya, A. Boisbunon, S. Chambon, L. Chapel, A. Corenflos, K. Fatras, N. Fournier, L. Gautheron, N. T. Gayraud, H. Janati, A. Rakotomamonjy, I. Redko, A. Rolet, A. Schutz, V. Seguy, D. J. Sutherland, R. Tavenard, A. Tong, and T. Vayer.  Pot: Python optimal transport.  Journal of Machine Learning Research, 22(78):1–8, 2021.
  * [9] M. Fréchet.  Les éléments aléatoires de nature quelconque dans un espace distancié.  In Annales de l’institut Henri Poincaré, volume 10, pages 215–310, 1948.
  * [10] N. Ho, X. L. Nguyen, M. Yurochkin, H. H. Bui, V. Huynh, and D. Phung.  Multilevel clustering via wasserstein means.  In Proceedings of the 34th International Conference on Machine Learning - Volume 70, ICML’17, page 1501–1509. JMLR.org, 2017.
  * [11] B. Horvath, Z. Issa, and A. Muguruza.  Clustering market regimes using the wasserstein distance.  arXiv:2110.11848, 2021.
  * [12] V. Huynh, N. Ho, N. Dam, X. Nguyen, M. Yurochkin, H. Bui, and D. Phung.  On efficient multilevel clustering via wasserstein distances.  Journal of Machine Learning Research, 2021.
  * [13] L. Kantorovich.  On the translocation of masses.  C. R. (Doklady) Acad. Sci. URSS (N.S.), 37:199–201, 1942.
  * [14] S. Lloyd.  Least squares quantization in pcm.  IEEE Transactions on Information Theory, 28(2):129–137, 1982.
  * [15] P. C. Álvarez Esteban, E. del Barrio, J. Cuesta-Albertos, and C. Matrán.  A fixed-point approach to barycenters in wasserstein space.  Journal of Mathematical Analysis and Applications, 441(2):744–762, 2016.
  * [16] J. MacQueen.  Some methods for classification and analysis of multivariate observations.  Proc. 5th Berkeley Symp. Math. Stat. Probab., Univ. Calif. 1965/66, 1, 281-297., 1967.
  * [17] G. I. Papayiannis, G. N. Domazakis, D. Drivaliaris, S. Koukoulas, A. E. Tsekrekos, and A. N. Yannacopoulos.  On clustering uncertain and structured data with wasserstein barycenters and a geodesic criterion for the number of clusters.  Journal of Statistical Computation and Simulation, 91(13):2569–2594, mar 2021.
  * [18] G. Peyré and M. Cuturi.  Computational optimal transport. With applications to data sciences.  Found. Trends Mach. Learn., 11(5-6):1–262, 2018.
  * [19] W. M. Rand.  Objective criteria for the evaluation of clustering methods.  Journal of the American Statistical Association, 66:846–850, 1971\.
  * [20] D. B. Rubin.  Multiple imputation for nonresponse in surveys. Hoboken, NJ: John Wiley & Sons, reprint of the 1987 original edition, 2004.
  * [21] M. Staib and S. Jegelka.  Wasserstein k-means++ for cloud regime histogram clustering.  In Proceedings of the Seventh International Workshop on Climate Informatics: CI 2017, 2017.
  * [22] K.-T. Sturm.  The space of spaces: curvature bounds and gradient flows on the space of metric measure spaces.  to appear in Memoirs AMS, 2020.
  * [23] I. Verdinelli and L. Wasserman.  Hybrid Wasserstein distance and fast distribution clustering.  Electronic Journal of Statistics, 13(2):5088 – 5119, 2019.
  * [24] C. Villani.  Topics in optimal transportation, volume 58 of Grad. Stud. Math. Providence, RI: American Mathematical Society (AMS), 2003.
  * [25] Y. Zhuang, X. Chen, and Y. Yang.  Wasserstein $k$-means for clustering probability distributions.  In A. H. Oh, A. Agarwal, D. Belgrave, and K. Cho, editors, Advances in Neural Information Processing Systems, 2022.

##  Appendix A Proofs

###### Proof of Theorem 5.1..

Let us first show that the loss function decreases monotonically. Calculate,

$\displaystyle\textstyle(1-\lambda^{(t+1)})L(\nu_{1}^{(t)},\ldots,\nu_{k}^{(t)},a^{(t)})$
$\displaystyle\textstyle=(1-\lambda^{(t+1)})\sum_{i=1}^{N}W_{2}^{2}(P_{i}\\#\mu_{i},P_{i}\\#\nu^{(t)}_{a_{i}^{(t)}})$
$\displaystyle\textstyle\overset{(\star)}{\geq}(1-\lambda^{(t+1)})\sum_{i=1}^{N}W_{2}^{2}(P_{i}\\#\mu_{i},P_{i}\\#\nu^{(t)}_{a_{i}^{(t+1)}})$
$\displaystyle\textstyle=\sum_{j=1}^{k}\Big{[}(1-\lambda^{(t+1)})\sum_{i:a_{i}^{(t+1)}=j}W_{2}^{2}(P_{i}\\#\mu_{i},P_{i}\\#\nu_{j}^{(t)})+\lambda^{(t+1)}\underbrace{W_{2}^{2}(\nu_{j}^{(t)},\nu_{j}^{(t)})}_{=0}\Big{]}$
$\displaystyle\textstyle\geq\sum_{j=1}^{k}\Big{[}(1-\lambda^{(t+1)})\sum_{i:a_{i}^{(t+1)}=j}W_{2}^{2}(P_{i}\\#\mu_{i},P_{i}\\#\nu_{j}^{(t+1)})+\lambda^{(t+1)}\underbrace{W_{2}^{2}(\nu_{j}^{(t)},\nu_{j}^{(t+1)})}_{\geq
0}\Big{]}$
$\displaystyle\textstyle\geq(1-\lambda^{(t+1)})\sum_{i=1}^{N}W_{2}^{2}(P_{i}\\#\mu_{i},P_{i}\\#\nu_{a_{i}^{(t+1)}}^{(t+1)})$
$\displaystyle\textstyle=(1-\lambda^{(t+1)})L(\nu_{1}^{(t+1)},\ldots,\nu_{k}^{(t+1)},a^{(t+1)}).$

Since $1-\lambda^{(t+1)}>0$ this shows the monotonicity. In inequality
$(\star)$ we have equality only if nothing changes in the assignment step,
i.e. there is equality if and only if

$\displaystyle\textstyle a_{i}^{(t)}=a_{i}^{(t+1)},~{}\forall
i\in\\{1,\ldots,N\\}.$

In that case the algorithm has converged and if there is some
$i\in\\{1,\ldots,N\\}$ such that $a_{i}^{(t)}\neq a_{i}^{(t+1)}$ we have a
strict inequality, since we only change the assignment in Algorithm 1 if $L$
decreases due to the new assignment.

Since only a finite number of possible assignments is possible because we are
aiming to cluster finitely many measures into finitely many clusters, the
algorithm strictly decreases $L$ and converges after a finite number of steps,
finishing the proof. ∎

Let us also prove that $\rho$, defined in (8), is a distance, or to be more
concrete just a generalized distance, i.e. we allow for $\rho$ to take the
value infinity.

######  Lemma A.1.

Let $(\mathcal{X},d)$ be a metric space and define on
$\mathcal{P}(\mathcal{X})$ the map
$\rho:\mathcal{P}(\mathcal{X})\times\mathcal{P}(\mathcal{X})\to[0,\infty]$ by

$\displaystyle\textstyle\rho(\mu,\nu):=\begin{cases}\textstyle\int_{\mathcal{X}\times\mathcal{X}}d(x,x^{\prime})\mathrm{d}(\mu\otimes\nu)(x,x^{\prime}),&\text{if
}\mu\neq\nu\\\ \textstyle 0,&\text{if }\mu=\nu.\end{cases}$

Then, $\rho$ is a generalized metric on $\mathcal{P}(X)$.

###### Proof.

Symmetry of the generalized metric is clear due to the symmetry of the
underlying distance $d$. Concerning the triangle inequality take
$\mu,\nu,\theta\in\mathcal{P}(\mathcal{X})$ which we assume to be different
from each other (otherwise there is nothing to prove). Furthermore, take three
independent random variables $X\sim\mu,Y\sim\nu,Z\sim\theta$. Then,

$\displaystyle\textstyle\rho(\mu,\theta)$ $\displaystyle=$
$\displaystyle\textstyle\int_{\mathcal{X}\times\mathcal{X}}d(x,x^{\prime})\mathrm{d}(\mu\otimes\theta)(x,x^{\prime})=\mathbb{E}[d(X,Z)]\leq$
$\displaystyle\leq$
$\displaystyle\textstyle\mathbb{E}[d(X,Y)+d(Y,Z)]=\rho(\mu,\nu)+\rho(\nu,\theta),$

which proves the triangle inequality for $\rho$. Concerning definiteness, if
$\mu=\nu$, we have $\rho(\mu,\nu)=0$ by definition.

Suppose now $\mu\neq\nu$, but $\rho(\mu,\nu)=0$. Furthermore, take two
independent random variables $X\sim\mu,~{}Y\sim\nu$. We then have,

$\displaystyle\textstyle
0=\int_{\mathcal{X}\times\mathcal{X}}d(x,x^{\prime})\mathrm{d}\mu\otimes\nu(x,x^{\prime})=\mathbb{E}[d(X,Y)].$

From this we conclude $d(X,Y)=0$ almost surely, which translates to $X=Y$
almost surely. This implies $\mu=\nu$ which is a contradiction. Therefore,
whenever $\mu\neq\nu$ we have $\rho(\mu,\nu)>0$. By contraposition,
$\rho(\mu,\nu)=0$ implies $\mu=\nu$, which concludes the prove of
definiteness. ∎

##  Appendix B Existence of “barycenters”

In this section we want to discuss the existence of minimizers in the general
barycenter updating step of Section 3.1. Thus a metric space $(\mathcal{X},d)$
and finitely many continuous maps into other metric spaces
$\varphi_{i}:\mathcal{X}\to\mathcal{X}_{i}$ are given. We observe the image
points $\tilde{x}_{i}=\varphi_{i}(x_{i})$ for $i=1,\ldots,n$. The question is:
under which conditions is

$\displaystyle\textstyle\inf\\{\sum_{i=1}^{n}\lambda_{i}d(\varphi_{i}^{-1}(\tilde{x}_{i}),y)^{2}:y\in\mathcal{X}\\}$
(16)

attained? Here, $\lambda_{i}>0,\sum_{i=1}^{n}\lambda_{i}=1$ denote some
weights. Since in this general formulation the maps $\varphi_{i}$ can be very
arbitrary, minimizers do not have to exist. However, if one of the maps
$\varphi_{i}$ is the identity on $\mathcal{X}$, we can guarantee existence
under mild assumptions.

######  Lemma B.1.

Assume there exists a metrizable topology $\tau$ on $\mathcal{X}$ which is
weaker than the topology induced by the metric $d$, and such that

  1. i)

$\varphi_{i}^{-1}(\tilde{x}_{i})$ is closed in $\tau$ for all
$i\in\\{1,\ldots,n\\}$,

  2. ii)

balls, i.e. $B_{d}(y,c):=\\{x\in\mathcal{X}:d(x,y)\leq c\\}$, are $\tau$
compact for all $y\in\mathcal{X},c>0$,

  3. iii)

$d(\cdot,\cdot)$ is lower semi continuous w.r.t. the product topology
$\tau\times\tau$,

  4. iv)

there exists $j\in\\{1,\ldots,n\\}$ such that $\varphi_{j}=Id_{\mathcal{X}}$.

Then (16) admits a minimizer.

###### Proof.

Let us first prove that for $d(\varphi_{i}^{-1}(\tilde{x}_{i}),y)$ there
exists $z\in\varphi_{i}^{-1}(\tilde{x}_{i})$ such that

$\displaystyle\textstyle
d(\varphi_{i}^{-1}(\tilde{x}_{i}),y)=\inf\\{d(x,y):x\in\varphi_{i}^{-1}(\tilde{x}_{i})\\}=d(z,y).$
(17)

To this end, take a sequence $\\{x_{n}\\}_{n\geq
1}\subset\varphi_{i}^{-1}(\tilde{x}_{i})$ such that

$\displaystyle\textstyle
d(x_{n},y)\searrow\inf\\{d(x,y):x\in\varphi_{i}^{-1}(\tilde{x}_{i})\\}.$

We then set $c:=d(x_{1},y)$, so $\\{x_{n}\\}_{n\geq 1}\subset B_{d}(y,c)$.
Since $B_{d}(y,c)$ is $\tau$-compact we obtain a subsequence
$(x_{n_{k}})_{k\geq 1}$ such that $x_{n_{k}}\overset{\tau}{\to}z\in
B_{d}(y,c)$. Since $\\{x_{n_{k}}\\}_{k\geq
1}\subset\varphi_{i}^{-1}(\tilde{x}_{i})$ and
$\varphi_{i}^{-1}(\tilde{x}_{i})$ is closed in $\tau$, we also have
$z\in\varphi_{i}^{-1}(\tilde{x}_{i})$. By lower semicontinuity of $d$ w.r.t.
$\tau\times\tau$ we obtain

$\displaystyle\textstyle
d(z,y)\leq\liminf_{k\to\infty}d(x_{n_{k}},y)=\inf\\{d(x,y):x\in\varphi_{i}^{-1}(\tilde{x}_{i})\\},$

which proves the existence of $z\in\varphi_{i}^{-1}(\tilde{x}_{i})$ such that
$d(z,y)=d(\varphi_{i}^{-1}(\tilde{x}_{i}),y)$.

Let us next prove that for a set $A=\varphi_{i}^{-1}(\tilde{x}_{i})$ the map
$x\mapsto d(A,x)$ is almost lower semi continuous w.r.t. $\tau$. Precisely,
let us show that for a sequence $\\{x_{n}\\}_{\geq 1}\subset\mathcal{X}$ such
that $x_{n}\overset{\tau}{\to}x\in\mathcal{X}$ and $d(x_{n},x)\leq c$ we have

$\displaystyle\textstyle\liminf_{n\to\infty}d(A,x_{n})\geq d(A,x).$ (18)

Without loss of generality, assume the left hand side of (18) is finite, as
otherwise there is nothing to prove. Therefore, we may assume the existence of
a constant $c^{\prime}>0$ such that $d(A,x_{n})\leq c^{\prime}$. By what we
already proved, for $n\geq 1$ there exists $z_{n}\in A$ such that
$d(A,x_{n})=d(z_{n},x_{n})$. We then have

$\displaystyle\textstyle d(z_{n},x)\leq
d(z_{n},x_{n})+d(x_{n},x)=d(A,x_{n})+d(x_{n},x)\leq c^{\prime}+c.$

By $\tau$ compactness of $B_{d}(x,c+c^{\prime})$ and since $A$ is closed in
$\tau$, there exists $z\in A\cap B_{d}(x,c+c^{\prime})$ such that
$z_{n}\overset{\tau}{\to}z$ up to switching to a subsequence. By lower semi
continuity of $d$ w.r.t. $\tau\times\tau$, i.e. assumption iii, we have

$\displaystyle\textstyle\liminf_{n\to\infty}d(A,x_{n})=\liminf_{n\to\infty}d(z_{n},x_{n})\geq
d(z,x)\geq d(A,x),$

which proves (18).

To prove the existence of minimizers for (16) we set

$\displaystyle\textstyle
V:=\inf\\{\sum_{i=1}^{n}\lambda_{i}d(\varphi_{i}^{-1}(\tilde{x}_{i}),y)^{2}:y\in\mathcal{X}\\}$

and take a sequence $\\{y_{n}\\}_{n\geq 1}\subset\mathcal{X}$ such that
$\sum_{i=1}^{n}\lambda_{i}d(\varphi_{i}^{-1}(\tilde{x}_{i}),y_{n})^{2}\searrow
V$. We assume without loss of generality that $d(\tilde{x}_{j},y_{n})\leq 2V$,
since $\varphi_{j}=Id_{\mathcal{X}}$, so we have $\\{y_{n}\\}_{n\geq 1}\subset
B_{d}(\tilde{x}_{j},2V)$. By compactness of $B_{d}(\tilde{x}_{j},2V)$ w.r.t.
$\tau$ we have the existence of a subsequence $\\{y_{n_{k}}\\}_{k\geq
1}\subset B_{d}(\tilde{x}_{j},2V)$ such that
$y_{n_{k}}\overset{\tau}{\to}y^{\star}$ for some $y^{\star}\in
B_{d}(\tilde{x}_{j},2V)$. We can now use the proved property (18) to obtain

$\displaystyle\textstyle
V=\liminf_{k\to\infty}\sum_{i=1}^{n}\lambda_{i}d(\varphi_{i}^{-1}(\tilde{x}_{i}),y_{n_{k}})^{2}\geq\sum_{i=1}^{n}\lambda_{i}d(\varphi_{i}^{-1}(\tilde{x}_{i}),y^{\star})^{2}.$

Therefore, we have found $y^{\star}\in\mathcal{X}$ such that

$\displaystyle\textstyle\sum_{i=1}^{n}\lambda_{i}d(\varphi_{i}^{-1}(\tilde{x}_{i}),y^{\star})^{2}=\inf\\{\sum_{i=1}^{n}\lambda_{i}d(\varphi_{i}^{-1}(\tilde{x}_{i}),y)^{2}:y\in\mathcal{X}),$

finishing the proof. ∎

In the proof we also showed the existence of minimizers in (7), i.e. points
used for “filling up”, under the same assumptions. This is precisely (17).

Let us discuss the assumptions of Lemma B.1. Note that assumptions i, ii and
iii are fulfilled in our two important examples of
$\mathcal{X}=\mathbb{R}^{d}$ and
$\mathcal{X}=\mathcal{P}_{2}(\mathbb{R}^{d})$. Indeed, for
$\mathcal{X}=\mathbb{R}^{d}$ this is clear, whereas for
$\mathcal{X}=\mathcal{P}_{2}(\mathbb{R}^{d})$ the role of the metric $\tau$ is
played by the metric induced by weak convergence of probability measures. For
details check Chapter $7$, _The Metric Side of Optimal Transportation_ , in
[24]. As for Assumption iv, this is not a big restriction in our two examples,
since this basically (in practical terms) means that for each cluster we need
at least one fully observed data point. Even if this is not the case, we can
still include the old barycenter with a small weight to enforce the existence
of minimizers, as done in Section 5 and Algorithm 1. This also helps
numerically, as otherwise infinite solutions for the generalized mean in the
Euclidean case, and the generalized Wasserstein barycenter in the case of
probability measures, may exist.

##  Appendix C Further simulation experiments

Let us discuss some additional simulation experiments to evaluate the proposed
method. In Section 7.1 we evaluated our method as an imputation method, but
without the imputation step it can also be viewed as a clustering algorithm.
In the case of Euclidean data it can be viewed as a clustering algorithm for
points with missing values. Thus, we can also compare it to classical
$k$-means after imputing the missing data with the more standard methods for
imputing missing values we used in 7.1. To recall, the standard methods for
imputing missing values we consider are mean imputation, median imputation,
multiple imputation, K-nearest neighbour imputation, and imputation through
linear regression. However, in the case of clustering Euclidean points with
missing data we also consider the $k$-pod method which was introduced in [4].
These authors consider the same loss function as we do in Example 3.1, i.e.
(9), however they propose a different algorithm. To perform the comparison, we
use the same simulation from Section 7.1, i.e. sampling from a Gaussian
Mixture Model with $k=5$ clusters in $10$ dimensions and different settings of
missing values. We use precisely the same simulated data, and for each of
those we apply classical $k$-means to cluster the points after imputing them
with a standard imputation method. For our proposed method, i.e. NA $k$-means,
and for the $k$-pod method, no imputation step is needed as these directly
cluster points with missing values. To evaluate the clustering results we use
the so-called Rand index, also known as the Rand score, introduced in [19]. It
compares a baseline clustering (in our case, the labels from the normal
distribution of the Gaussian Mixture Model to which an observation belongs) to
a clustering of the data using an imputation and/or clustering algorithm. The
Rand score takes values in $[0,1]$, with $0$ meaning nothing is clustered the
same, and $1$ meaning the clusterings coincide except for renaming of
clusters. In particular, a higher Rand score corresponds to a better
clustering method. In Table 2 we report the corresponding Rand scores $\pm$
standard errors for these simulations.

Table 2: Rand scores for Euclidean simulation.

settings | NA $k$-means | mean imp. | median imp. | multiple imp. | KNN | LR |  $k$-pod
---|---|---|---|---|---|---|---
$\pm$

We can see that the introduced method NA $k$-means either performs best or
lies within the standard error of the best performing method for each setting
of missing values when considering it solely as clustering algorithm for
points with missing data. The naive approaches of mean and median imputation
are outperformed.

#### Generalizing our results

In our simulation scenarios we have two parameters concerning the creation of
missing values in the data. One is the _quantile_ below which we remove values
for each dimension. The second one is $h_{j}$, which is a uniform random
variable per cluster which is then correlated with the cluster size to decide
what percentage of missing values per cluster is generated additionally. In
Section 7.1 and in the previous paragraph, the _quantile_ is chosen as $25\%$
quantile and we set $h_{j}\sim U(0,2/k)$ where $k$ is the number of clusters,
i.e. $5$, in the considered simulations. Essentially, this fixes the share of
missing values per simulated dataset. However, we want to vary these two
parameters in order to show the robustness of our method. In Table 3 and Table
4 we report the corresponding results analogous to Table 1 and Table 2,
respectively, when choosing $5$ different choices of the _quantile_ and
$h_{j}$.

Table 3: Gromov Wasserstein distance for Euclidean simulation with varying
share of missing values.

settings | NA $k$-means. | NA $k$-means-m | mean imp. | median imp. | multiple imp. | KNN | LR
---|---|---|---|---|---|---|---
$10\%$ quantile and $h_{j}\sim U(0,0.1)$
$~{}\pm$
$20\%$ quantile and $h_{j}\sim U(0,0.2)$
$~{}\pm$
$30\%$ quantile and $h_{j}\sim U(0,0.3)$
$~{}\pm$
$40\%$ quantile and $h_{j}\sim U(0,0.4)$
$~{}\pm$
$50\%$ quantile and $h_{j}\sim U(0,0.5)$
$~{}\pm$

From Table 3 we conclude that the proposed method NA $k$-means outperforms all
the other considered methods independent of the parameters needed for the
missing value generation, when considering the Gromov Wasserstein distance as
evaluation measure.

Table 4: Rand scores for Euclidean simulation with various shares of missing
values.

settings | NA $k$-means | mean imp. | median imp. | multiple imp. | KNN | LR |  $k$-pod
---|---|---|---|---|---|---|---
$10\%$ quantile and $h_{j}\sim U(0,0.1)$
$\pm$
$20\%$ quantile and $h_{j}\sim U(0,0.2)$
$\pm$
$30\%$ quantile and $h_{j}\sim U(0,0.3)$
$\pm$
$40\%$ quantile and $h_{j}\sim U(0,0.4)$
$\pm$
$50\%$ quantile and $h_{j}\sim U(0,0.5)$
$\pm$

In Table 4 we can see that also w.r.t. the Rand score, when NA $k$-means is
solely viewed as clustering procedure for Euclidean points with missing
values, it competes with the best performing imputation + clustering
techniques.

##  Appendix D An example highlighting NA $k$-means strength

We give a toy example highlighting a scenario in which the proposed method
works well, and which hopefully helps building some intuition. Therefore, we
work in dimension 2 and consider Euclidean points with missing values. In
particular, we consider three clusters, the first and third one with $200$
points each and the second one with $100$ points, coming from the three
Gaussians

$\displaystyle\textstyle\mathcal{N}\Big{(}\Big{(}\begin{matrix}3\\\
3\end{matrix}\Big{)},\Big{(}\begin{matrix}1&0\\\
0&1\end{matrix}\Big{)}\Big{)},~{}~{}~{}\mathcal{N}\Big{(}\Big{(}\begin{matrix}0\\\
-1.5\end{matrix}\Big{)},\Big{(}\begin{matrix}1&0.6\\\
0.6&1\end{matrix}\Big{)}\Big{)},~{}~{}~{}\mathcal{N}\Big{(}\Big{(}\begin{matrix}-3\\\
3\end{matrix}\Big{)},\Big{(}\begin{matrix}1&0\\\
0&1\end{matrix}\Big{)}\Big{)}.$

In order to model _structural missingness of data_ , we assume that we only
observe the $y$ coordinate of the points if it is positive, so missing values
are created for non-positive $y$ coordinates. This can be seen in Figure 3.

Figure 3: (left): The true points with true labels. (right) The observed
points with the true labels.

In order to cluster in this situation, most methods for imputation lead to a
bias by shifting up the points with the unobserved $y$ coordinates.
Especially, as in this case there is essentially only one cluster with points
having negative $y$ coordinates and they are the majority of that cluster. NA
$k$-means however, does not introduce this bias and therefore recovers the
true cluster membership much better. The corresponding rand scores in this
example are in Table 5.

Table 5: Rand scores.

NA $k$-means | mean imp. | median imp. | multiple imp. | KNN | LR |  $k$-pod
---|---|---|---|---|---|---

In Figure 4 below we present the corresponding observed clusterings, i.e. we
plot the true points and color them with the labels from the corresponding
clustering method.

Figure 4: (first row, left to right): NA $k$-means, mean imputation, median
imputation, multiple imputation. (second row, left to right): KNN, LR,
$k$-pod. The reader is encouraged to zoom in.

We can see that most methods solely base their clustering on the values of the
$x-$coordinate. The reason is that, by imputing the points coming from the
turquoise cluster upwards, these points end up having $y-$values similar to
the $y-$values of the points with no missing values. The clustering coming
from NA $k$-means recovers the original cluster membership of the points the
best.