# Roman Numeral Analysis with Graph Neural Networks: Onset-wise Predictions
from Note-wise Features

###### Abstract

Roman Numeral analysis is the important task of identifying chords and their
functional context in pieces of tonal music. This paper presents a new
approach to automatic Roman Numeral analysis in symbolic music. While existing
techniques rely on an intermediate lossy representation of the score, we
propose a new method based on Graph Neural Networks (GNNs) that enable the
direct description and processing of each individual note in the score. The
proposed architecture can leverage notewise features and interdependencies
between notes but yield onset-wise representation by virtue of our novel edge
contraction algorithm. Our results demonstrate that ChordGNN outperforms
existing state-of-the-art models, achieving higher accuracy in Roman Numeral
analysis on the reference datasets. In addition, we investigate variants of
our model using proposed techniques such as NADE, and post-processing of the
chord predictions. The full source code for this work is available at
https://github.com/manoskary/chordgnn

##  1 Introduction

Automatic Chord Recognition is one of the core problems in Music Information
Retrieval. The task consists of identifying the harmonies or chords present in
a musical piece. Various methods have been proposed to address this task using
either an audio or symbolic representation of the music [1]. In the symbolic
domain, most approaches focus on the related and arguably more complex problem
of Automatic Roman Numeral Analysis, which is a functional harmony analysis
problem that has its roots in musicological research of Western classical
music.

Roman Numeral Analysis is a notational system used in music theory to analyze
chord progressions and identify the relationship between chords in a given
key. In this system, each chord in a piece of music is assigned a Roman
numeral based on its position within the key’s scale. For example, in the key
of C major, the I chord is C major, the IV chord is F major, and the V chord
is G major. Roman Numerals are an important tool for understanding and
analyzing the harmonic structure of music, and they are a valuable resource
for musicians, composers, and arrangers alike.

In Music Information Retrieval, a lot of work has been done to automate Roman
Numeral analysis. However, current approaches still face significant
challenges. Some of these are related to the large chord symbol vocabulary. A
common way to address this problem is to divide a Roman Numeral into several
components (e.g., key, degree, inversion) and transform the analysis into a
multitask learning scenario. However, multitask approaches themselves face
challenges with interdependencies among tasks. Lastly, Roman Numeral analysis
faces a score representation problem related to existing models such as CNNs
whose inputs must be in fixed-sized chunks. Recent state-of-the-art approaches
follow an audio-inspired strategy, dividing a musical score into fixed-length
time frames ("windows") which are then processed by a Convolutional Recurrent
Neural Network (CRNN). However, such a representation is unnatural for scores
and has the added practical disadvantage of being time-limited (for example
regarding notes extending beyond the current window) and, due to the fixed-
length (in terms of score time) constraint, capturing varying amounts of
musically relevant context.

In this paper, we propose a new method for automatic Roman Numeral analysis
based on Graph Neural Networks that can leverage note-wise information to
address the score representation issue. Our model, ChordGNN, builds on top of
existing multitask approaches but introduces several novel aspects, including
a graph convolutional architecture with an edge contraction pooling layer that
combines convolution at the note level but yields the learned representation
at the onset level.

Our proposed method, ChordGNN, is evaluated on a large dataset of Western
classical music, and the experimental results demonstrate that it outperforms
existing state-of-the-art methods, in terms of the commonly used Chord Symbol
Recall measure. To address the interdependencies among tasks we investigate
the effect of post-processing and other proposed techniques such as NADE and
gradient normalization. Finally, we look at a qualitative musical example and
compare our model’s predictions with other state-of-the-art models.

##  2 Related Work

There is a big body of literature covering the topic of Automatic Chord
Recognition applied in the audio domain; however, in our work, we focus on the
problem of automatic Roman Numeral Analysis in the symbolic domain. It
consists of labeling the chords and harmonic progressions in a piece of music
using Roman Numerals, where each numeral represents a chord built on a
particular scale degree. Numerous approaches have tried to automate Roman
Numeral analysis or infer harmonic relations between chords. Notable work
includes statistical models such as Melisma [2], HMM-based models [3], and
grammar-based approaches [4].

Figure 1: A Roman Numeral analysis for two bars for four-part harmony in $C$
major. Capital letters stand for major quality and lowercase for minor
quality. The third chord has a dominant seven as its primary degree and the
dominant of $C$ major as its secondary degree. The $V^{6}_{5}$ indicates a
major with a seven quality in second inversion. The bass (lowest chord note)
of that chord is $F$ sharp, the root is $D$, and the local key is $C$ major.

In recent years, research has shifted towards a deep learning and data-driven
approach. Due to the large vocabulary of possible Roman Numerals, the problem
has been divided into several component subtasks, thus resulting in a
multitask learning setting [5]. As a multitask problem, a Roman Numeral is
characterized by the following components: the primary and secondary degree
(as illustrated in Figure 1), the local key at the time point of prediction,
the root of the chord, the inversion of the chord, and the quality (such as
major, minor, 7, etc.). Although the root can be derived from the other
components, it was pointed out by [6] that redundancy is assisting Roman
Numeral analysis systems to learn. An example of Roman Numerals and their
components can be viewed in Figure 1. Recent state-of-the-art approaches
decompose the numeral prediction task to the simultaneous prediction of those
6 components [5, 6, 7, 8, 9].

Most deep learning approaches to Roman Numeral analysis are inspired by work
in audio classification, cutting a score into fixed-size chunks (in terms of
some constant score time unit; e.g., a 32nd note) and using these as input to
deep models. Using this quantized time frame representation, [6] introduced a
CRNN architecture to predict Roman Numerals. Other work has continued to build
on the latter by introducing more tasks to improve performance such as the
AugmentedNet model [7], or introducing intra-dependent layers to inform in an
orderly fashion the prediction of one task with the previously predicted task,
such as the model introduced by [8]. Other architectures, such as the CSM-T
model, have demonstrated good results by introducing modular networks which
treat a score as a sequence of notes ordered first by onset and then by
pitch[9].

Figure 2: Different representations of the score excerpt shown in the middle.
Top: quantized time frame representation, bottom: graph representation.

Should a musicologist perform music analysis on a piece of music, they would
consider the individual notes existing in the score. Thus, a time frame
representation would come across as unnatural for symbolic music and in
particular for such an analysis task. In this paper, we present a method that
no longer treats the score as a series of quantized frames but rather as a
partially ordered set of notes connected by the relations between them, i.e.,
a graph. A visual comparison of the two representations is shown in Figure 2.
Recently, modeling scores as graphs has also been demonstrated to be
beneficial for problems such as expressive performance generation [10],
cadence detection [11], voice separation [12], or boundary detection [13].

Automatic Roman Numeral analysis, as a multitask problem, is mostly tackled
with hard parameter-sharing models. These models share part of the model
across all tasks as an encoder, and then the common embeddings are branched to
a classification model per task [6, 7, 8]. However, some approaches separate
tasks from this paradigm to a more modular or soft parameter sharing approach
[9].

In the field of multitask learning, a lot of research has been done on the
problem of conflicting gradients during backpropagation in hard parameter-
sharing models. Issues with multi-objective optimization have been early
addressed by Zhang et al. [14] and recent solutions have been proposed for the
multitask setting in the form of dynamic task prioritization [15], gradient
normalization [16], rotation matrices [17], or even game-theoretic approaches
[18]. In our work, we experimentally evaluate some of these techniques in the
multitask setting to investigate whether Roman Numeral analysis subtasks
conflict with each other (see Section 5.2).

##  3 Methodology

Figure 3: The proposed Architecture Chord-GNN

###  3.1 Roman Numeral Analysis

We already discussed, in Section 2, how Roman Numeral analysis can be viewed
as a multi-task problem. In this section, we describe in detail the additional
tasks introduced by [7] that we also use for training and prediction. First,
let us assume that the prediction can be broken down into specific time
points, and each time point is attributed to a unique onset in the score.

The Roman Numeral prediction can be viewed as a simultaneous prediction of the
local key, degree (primary and secondary), quality, inversion, and root. Each
one of these tasks is a categorical, multiclass classification problem.
However, [7] indicated that only three tasks would be sufficient for $~{}98\%$
of the Roman Numeral annotations in our dataset (detailed in Section 4.1).
These three tasks comprise the prediction of a restricted vocabulary of common
Roman Numeral symbols in combination with the local key and the inversion. We
refer to Roman Numeral prediction involving the 5 tasks as conventional RN,
and the combined prediction of key, inversion, and restricted RN vocabulary
alternative RN, as $RN_{alt}$, in accordance with [7].

Several other tasks have been introduced that have been shown to improve the
performance of related models [7]. These include the Harmonic Rhythm, which is
used to infer the duration of a Roman Numeral at a given time point; the
Tonicization task, a multiclass classification task that refers to a tonicized
key implied by the Roman Numeral label and is complementary to the local key;
the Pitch Class Sets task, which includes a vocabulary of different pitch
class sets, and the Bass task, which aims to predict the lowest note in the
Roman Numeral label.

###  3.2 Graph Representation of Scores

Our approach to automatic Roman Numeral analysis no longer treats the score as
a sequence of quantized time frames but rather as a graph, which permits us to
specify note-wise information such as pitch spelling, duration, and metrical
position. We use graph convolution to model interdependencies between notes.
We model our score generally following Karystinaios and Widmer [11], but we
opt for a heterogeneous graph convolution approach, i.e., including different
edge relations/types. Furthermore, we develop an edge contraction pooling
layer that learns onset-wise representations from the note-wise embeddings and
therefore yields a sequence.

After the edge contraction, we follow [6, 7, 8] by adding to the graph
convolution a sequence model for the hard-sharing part of our model, and
simple shallow multi-layer perceptron heads for each task. In essence, we
replace the CNN encoder that works on quantized frames of the score in
previous approaches, with a graph convolutional encoder followed by an edge
contraction layer. Our proposed architecture is shown in Figure 3.

The input to the GNN encoder is an attributed graph $G=(V,E,X)$ where $V$ and
$E$ denote its node and edge sets and $X$ represents the node feature matrix,
which contains the features of the notes in the score. For our model, we used
pitch spelling, note duration, and metrical position features.

Given a musical piece, the graph-building process creates a set of edges $E$,
with different relation types $\mathcal{R}$. A labeled edge $(u,r,v)$ of type
$r$ between two notes $u,v$ belongs to $E$ if the following conditions are
met:

  * •

notes starting at the same time:
$on(u)=on(v)\xrightarrow{}r=\textrm{onset}$

  * •

note starting while the other is sounding: $on(u)>on(v)\land on(u)\leq
on(v)+dur(v)\xrightarrow{}r=\textrm{during}$

  * •

note starting when the other ends:
$on(u)+dur(u)=on(v)\xrightarrow{}r=\textrm{follow}$

  * •

note starting after a time frame when no note is sounding:
$on(u)+dur(u)<on(v)\land\nexists v^{\prime}\in V,\;on(v^{\prime})<on(v)\land
on(v^{\prime})>on(u)+dur(u)\xrightarrow{}r=\textrm{silence}$

###  3.3 Model

In this section, we introduce and describe ChordGNN, a Graph Convolutional and
Recurrent Neural Network. The structure of the network is visually outlined in
Figure 3. ChordGNN uses heterogeneous graphSAGE [19] convolutional blocks
defined as:

$\displaystyle\begin{aligned}
\mathbf{h}_{\mathcal{N}_{r}(v)}^{(l+1)}&=\mathrm{mean}\left(\\{\mathbf{h}_{u}^{l},\forall
u\in\mathcal{N}_{r}(v)\\}\right)\\\
\mathbf{h}_{v_{r}}^{(l+1)}&=\sigma\left(W\cdot\mathrm{concat}(\mathbf{h}_{v}^{l},\mathbf{h}_{\mathcal{N}_{r}(v)}^{l+1})\right)\\\
\mathbf{h}_{v}^{(l+1)}&=\frac{1}{|\mathcal{R}|}\sum_{r\in\mathcal{R}}\mathbf{h}_{v_{r}}^{(l+1)}\end{aligned}$
(1)

where $\mathbf{h}_{v}^{(0)}=\mathbf{x}_{v}$ and $\mathbf{x}_{u}$ is the input
features for node $u$, $\mathcal{N}(u)$ are the neighbors of node $u$, and
$\sigma$ is a ReLU activation function. We name the output representations of
all nodes after graphSAGE convolution $H=\\{h^{(L)}_{u}\mid u\in V\\}$ where
$L$ is the total number of convolutional layers.

Given the hidden representation $H$ of all nodes, and onset edges
$E_{\textrm{On}}=\\{(u,v)\mid on(u)=on(v)\\}$, the onset edge contraction
pooling is described by the following equations: first, we update the hidden
representation with a learned weight, $H^{\prime}=HW^{(\textrm{cpool})}$.
Subsequently we need to unify the representations for every node $u$, such
that $\forall
v\in\mathcal{N}_{\textrm{On}}(v),\;h^{(\textrm{cp})}_{u}=h^{(\textrm{cp})}_{v}$:

$h^{(\textrm{cp})}_{u}=h_{u}+\sum_{v\in\mathcal{N}_{\textrm{On}}(v)}h_{v}$ (2)

where, $h_{u}$ and $h_{v}$ belong to $H^{\prime}$. Subsequently, we filter the
vertices:

$V^{\prime}=\\{v\in V|\;\forall u\in V,\;(v,u)\in E_{\textrm{On}}\implies
u\notin V^{\prime}\\}$ (3)

Therefore, $H^{(cp)}=\\{h_{u}^{(cp)}\mid\forall u\in V^{\prime}\\}$ are the
representations obtained. Sorting the representations by the onset on which
they are attributed we obtain a sequence
$S=[h_{u_{1}}^{(cp)},h_{u_{2}}^{(cp)},\dots h_{u_{k}}^{(cp)}]$ such that
$on(u_{1})<on(u_{2})<\dots<on(u_{k})$.

The sequence $S$ is then passed through an MLP layer and 2 GRU layers. This
concludes the hard-sharing part of our model. Thereafter, an MLP head is
attached per task, as shown in Figure 3.

For training, we use the dynamically weighted loss introduced by [20]. The
total loss $\mathcal{L}_{tot}$ of our network is calculated as a weighted sum
of the individual losses for every task, where the weights are learned during
training:

$\mathcal{L}_{\textrm{tot}}=\sum_{t\in\mathcal{T}}\mathcal{L}_{t}*\frac{1}{2\gamma_{t}^{2}}+\log(1+\gamma_{t}^{2})$
(4)

where $\mathcal{T}$ is the set of tasks; $\mathcal{L}_{t}$ is the cross-
entropy loss relating to task $t$; the $\gamma_{t}$ are learned scalars that
give the weight for each task $t$; and the $\log$ expression is a
regularization term [20].

Figure 4: Post-processing of Roman Numeral predictions.

####  3.3.1 Post-processing

We enhance our model with a post-processing phase after the model has been
trained. The post-processing phase combines the logits of all tasks’
predictions by concatenating them and, then, feeds them to a single-layer
bidirectional LTSM block. Then, again the embeddings of the sequential block
are distributed to 11 one-layer MLPs, one for each task. The post-processing
block is sketched in Figure 4.

##  4 Experiments and Corpora

In the experiments, we compare our model, ChordGNN, with other recent models
for automatic Roman Numeral analysis. We run experiments with our model in the
exact same way as described in the paper [7], including the specific data
splits, so that our results are directly comparable to the figures reported
there. A detailed comparison of the results will be given in Table 1.
Furthermore, we develop variants of our model using proposed techniques such
as NADE [8], and post-processing of the chord predictions. We report a
configuration study of our model on the use of gradient normalization
techniques and NADE that should improve results on Multi-Task learning
scenarios and avoid common Multi-Task Learning problems such as conflicting
gradients. Lastly, we compare our model with the updated version v1.9.1 of the
state-of-the-art model Augmented-Net [21] and datasets.

###  4.1 Datasets

| Model | Key | Degree | Quality | Inversion | Root | RN | RN (Onset) | $\textrm{RN}_{\textrm{alt}}$
---|---|---|---|---|---|---|---|---|---
BPS | Micchi (2020) | 82.9 | 68.3 | 76.6 | 72.0 | - | 42.8 | - | -
CSM-T (2021) | 69.4 | - | - | - | 75.4 | 45.9 | - | -
AugNet (2021) | 85.0 | 73.4 | 79.0 | 73.4 | 84.4 | 45.4 | - | 49.3
ChordGNN (Ours) | 79.9 | 71.1 | 74.8 | 75.7 | 82.3 | 46.2 | 46.6 | 48.6
ChordGNN+Post (Ours) | 82.0 | 71.5 | 74.1 | 76.5 | 82.5 | 49.1 | 49.4 | 50.4
Full |  AugNet (2021) | 82.9 | 67.0 | 79.7 | 78.8 | 83.0 | 46.4 | - | 51.5
ChordGNN (Ours) | 80.9 | 70.1 | 78.4 | 78.8 | 84.8 | 48.9 | 48.4 | 50.4
ChordGNN+Post (Ours) | 81.3 | 71.4 | 78.4 | 80.3 | 84.9 | 51.8 | 51.2 | 52.9
Table 1: Model comparison on two different test sets, the Beethoven Piano
Sonatas (BPS), and the full test set. $RN$ stands for Roman Numeral,
$RN_{alt}$ for the alternative Roman Numeral computations discussed in Section
3.1. $RN(Onset)$ refers to onset-wise prediction accuracy, all other scores
use the CSR score (see Section 5). Note that model CSM-T reports Mode instead
of Quality.

For training and evaluation, we combined six data sources into a single "Full"
Dataset of Roman Numeral annotations in accordance with [7]: the Annotated
Beethoven Corpus (ABC) [22]; the annotated Beethoven Piano Sonatas (BPS)
dataset [5]; the Haydn String Quartets dataset (HaydnSun) [23]; the TAVERN
dataset [24]; a part of the When-in-Rome (WiR) dataset [25, 26]; and the Well-
Tempered-Clavier (WTC) dataset [25] which is also part of the WiR dataset.

Training and test splits for the full dataset were also provided by [7]. It is
worth noting that the BPS subset splits were already predefined in [5]. In
total, approximately 300 pieces were used for training, and 56 pieces were
used for testing, proportionally taken from all the different data sources. We
draw a distinction for the BPS test set, which includes 32 Sonata first
movements and for which we ran an additional experiment. The full test set
also includes the 7 Beethoven piano sonatas.

In addition to the above datasets, we include data augmentations identical to
the ones described in [7]: texturization and transposition. The texturization
is based on a dataset augmentation technique introduced by [27]. The
transposition augmentation boils down to transposing a score to all the keys
that lie within a range of key signatures that have up to 7 flats or sharps.
It should be noted that the augmentations are only applied in the training
split.

For our last experiment (to be reported on in Section 5.3 below), we add
additional data that were recently introduced by [21]. The additional data
include the annotated Mozart Piano Sonatas (MPS) dataset [28] for which we
also applied the aforementioned augmentations.

###  4.2 Configuration

For all our experiments, we train our network with the AdamW optimizer. We fix
our architecture with a hidden size of $256$, a learning rate of $0.0015$, a
weight decay of $0.005$, and a dropout of $0.5$ which is applied to each
learning block of our architecture.

##  5 Results

As an evaluation metric, we use Chord Symbol Recall (CSR) [29] where for each
piece, the proportion of time is collected during which the estimated label
matches the ground truth label. We apply the CSR at the 32nd note granularity
level, in accordance with [6, 7, 9].

###  5.1 Quantitative Results

In the first experiment, which compares our ChordGNN to existing state-of-the-
art approaches, we evaluate the full dataset, but also the annotated Beethoven
Piano Sonatas (BPS) [5] subset, which many previous approaches had also used.
The results are shown in Table 1. We present the CSR scores (where they are
applicable) for Local Key, Degree, Quality, Inversion, Root, conventional
Roman Numeral, and Alternative Roman Numeral (see Section 3). Furthermore, we
include the onset-wise accuracy score for our models’ conventional Roman
Numeral predictions.

On the BPS subset, we compare our model ChordGNN with the Micchi (2020) model
[6], the CSM-T (2021) model [9] and the AugmentedNet 2021 model [7]. Our
results on Roman Numeral prediction surpass all previous approaches. Note that
the AugmentedNet model exhibits higher prediction scores on the individual
Key, Degree, Quality, and Root tasks, which are used jointly for the
prediction of the Roman numeral. These results indicate that our model obtains
more meaningfully interrelated predictions, with respect to the Roman numeral
prediction, resulting in a higher accuracy score.

Moreover, we compare ChordGNN to AugmentedNet on the full test dataset. Our
model surpasses AugmentedNet with and without post-processing in all fields
apart from local key prediction and quality. Our model obtains up to $11.6\%$
improvement in conventional Roman Numeral prediction.

In both experiments, post-processing has been shown to improve both $RN$ and
$RN_{alt}$. However, ChordGNN without post-processing already surpasses the
other models.

Variant | RN | $\textrm{RN}_{\textrm{alt}}$
---|---|---
ChordGNN (Baseline) | $46.1\pm 0.003$ | $47.8\pm 0.007$
ChordGNN + WLoss | $\mathbf{48.9}\pm 0.001$ | $\mathbf{50.4}\pm 0.010$
ChordGNN + Rotograd | $45.5\pm 0.003$ | $47.1\pm 0.005$
ChordGNN + R-GradN | $45.2\pm 0.006$ | $46.7\pm 0.005$
ChordGNN + NADE | $48.2\pm 0.005$ | $49.9\pm 0.005$
Table 2: Configuration Study: Chord Symbol Recall on Roman Numeral analysis on
the full test set. $RN$ stands for Roman Numeral, $RN_{alt}$ refers to the
alternative Roman Numeral computations discussed in section 3.1. WLoss stands
for the dynamically weighted loss described in Section 3, and R-GradN stands
for Rotograd with Gradient Normalization. Every experiment is repeated 5 times
with the same ChordGNN model as Table 1 without post-processing. Figure 5: A
comparison between the human annotation, AugmentedNet, and ChordGNN on a
passage of Haydn’s string quartet op.20 No.3 movement 4. The red (wrong)
markings on Human Analysis and AugNet (2022) are from [21]

###  5.2 Configuration Study

For a systematic study of multitask training, we investigated the effects of
extension modules, gradient normalization techniques, and learnable weight
loss. In detail, we test 5 configurations using as baseline the ChordGNN model
(without post-processing) with standard CE loss and no weighing. Furthermore,
we test our proposed architecture using the dynamically weighted loss
described in Section 3.3 (same as the model in Table 1), Rotograd [17] and
GradNorm [16] for Gradient Normalization, and NADE [8]. The models are run on
the Full data set described above and averaged over five runs with random
initialization. The results, summarized in Table 2, suggest that using the
dynamically weighted loss yields better results compared to other methods such
as the Baseline or Gradient Normalization techniques. Furthermore, the
dynamically weighted loss is comparable to NADE but also more robust on
Conventional Roman Numeral prediction on our datasets.

###  5.3 Latest developments

Our last experiment focuses on specific developments that have very recently
been published in Nápoles López’s Ph.D. thesis [21]. In the thesis, three
additional tasks, related to predicting the components of a canonical
representation of the current chord, as implied by the Roman Numeral, were
proposed and the dataset was extended with the Annotated Mozart Piano Sonatas
(MPS) corpus [28], as mentioned in Section 4.1 above.

To test the relevance of these updates, we trained an adapted version of our
model, now with 11+3=14 individual tasks and including the Mozart data. It
turns out that the updated model improves significantly in performance,
achieving a $53.5$ CSR score on conventional Roman Numeral (compare this to
row "ChordGNN (Ours)" in Table 1). Furthermore, post-processing can improve
the results by up to two additional percentage points.111Unfortunately, we
cannot directly compare these numbers to [21], as their results are not
reported in comparable terms.

###  5.4 A Musical Example

In Figure 5, we look at a comparison between the human annotations,
AugmentedNet and Chord-GNN predictions (The musical excerpt is taken from
Nápoles López’s thesis [21], and the predictions relate to the new models
trained as described in the previous section.). Marked in red are false
predictions, and marked in yellow are correct predictions of the model with
wrong ground-truth annotations. Both models’ predictions are very similar to
the human analysis. However, our model correctly predicts the initial pickup
measure annotation. In measure 2, the ground truth annotation marks a tonic in
first inversion; however, the viola at that point is lower than the cello and
therefore the chord is actually in root position. Both models obtain a correct
prediction at that point. Subsequently, our model predicts a harmonic rhythm
of eighth notes, which disagrees with the annotator’s half-note marking.
Analyzing the underlying harmony in that passage, we can justify our model’s
choices.

The human annotation suggests that the entire second half of the 2nd measure
represents a $vii^{o}$ chord. However, it should not be in the first
inversion, as the cello plays an F# as the lowest note (which is the root of
$vii^{o}$). The AugNet analysis faces the same issue, in contrast with the
predictions of ChordGNN. However, there are two conflicting interpretations of
the segment. First, the $vii^{o}$ on the third beat is seen as a passing chord
between the surrounding tonic chords, leading to a dominant chord in the next
measure. Alternatively, the $vii^{o}$ could already be part of a prolonged
dominant harmony (with passing chords on the offbeats) leading to the $V^{7}$.
The ChordGNN solution accommodates both interpretations as it doesn’t attempt
to group chords at a higher level, treating each eighth note as an individual
chord rather than a passing event. The other two solutions prefer the second
option.

##  6 Conclusion

In this paper, we presented ChordGNN, a model for automatic Roman Numeral
analysis in symbolic music, based on a note-level, graph-based score
representation. We showed that ChordGNN improves on other state-of-the-art
models, and that post-processing can further improve the accuracy of the
predictions. A configuration study suggests that gradient normalization
techniques or techniques for carrying prediction information across tasks are
not particularly beneficial or necessary for such a model.

Follow-up work will focus on strengthening the robustness of our models by
pre-training with self-supervised methods on large corpora. We believe that
such pre-training can be beneficial for learning helpful intrinsic musical
information. Such a step is crucial since more data improves predictions but
Roman Numeral annotations are hard to find or produce. Moreover, we aim to
enrich the number of tasks for joint prediction by including higher-level
analytical targets such as cadence detection and phrase boundary detection.
Finally, we aim to extend our method to the audio domain.

##  7 Acknowledgements

We gratefully acknowledge the musical analysis of the $vii^{o}$ passage in
Fig. 5 (Section 5.4) that was offered by an anonymous reviewer, and which we
took the liberty of adopting for our text. This work is supported by the
European Research Council (ERC) under the EU’s Horizon 2020 research &
innovation programme, grant agreement No. 101019375 (“Whither Music?”), and
the Federal State of Upper Austria (LIT AI Lab).

## References

  * [1] J. Pauwels, K. O’Hanlon, E. Gómez, M. Sandler _et al._ , “20 years of Automatic Chord Recognition from Audio,” in _Proceedings of the International Society for Music Information Retrieval Conference (ISMIR)_ , 2019\.
  * [2] D. Temperley, _The cognition of basic musical structures_. MIT press, 2004.
  * [3] C. Raphael and J. Stoddard, “Functional Harmonic Analysis Using Probabilistic Models,” _Computer Music Journal_ , vol. 28, no. 3, pp. 45–52, 2004.
  * [4] J. P. Magalhaes and W. B. de Haas, “Functional Modelling of Musical Harmony: an experience report,” _ACM SIGPLAN Notices_ , vol. 46, no. 9, pp. 156–162, 2011.
  * [5] T.-P. Chen, L. Su _et al._ , “Functional Harmony Recognition of Symbolic Music Data with Multi-task Recurrent Neural Networks.” in _Proceedings of the International Society for Music Information Retrieval Conference (ISMIR)_ , 2018.
  * [6] G. Micchi, M. Gotham, and M. Giraud, “Not all roads lead to Rome: Pitch representation and model architecture for automatic harmonic analysis,” _Transactions of the International Society for Music Information Retrieval (TISMIR)_ , vol. 3, no. 1, pp. 42–54, 2020.
  * [7] N. Nápoles López, M. Gotham, and I. Fujinaga, “AugmentedNet: A Roman Numeral Analysis Network with Synthetic Training Examples and Additional Tonal Tasks.” in _Proceedings of the International Society for Music Information Retrieval Conference (ISMIR)_ , 2021.
  * [8] G. Micchi, K. Kosta, G. Medeot, and P. Chanquion, “A deep learning method for enforcing coherence in Automatic Chord Recognition.” in _Proceedings of the International Society for Music Information Retrieval Conference (ISMIR)_ , 2021.
  * [9] A. P. McLeod and M. A. Rohrmeier, “A modular system for the harmonic analysis of musical scores using a large vocabulary,” in _Proceedings of the International Society for Music Information Retrieval Conference (ISMIR)_ , 2021\.
  * [10] D. Jeong, T. Kwon, Y. Kim, and J. Nam, “Graph Neural Network for Music Score Data and Modeling Expressive Piano Performance,” in _Proceedings of International Conference on Machine Learning (ICML)_ , 2019.
  * [11] E. Karystinaios and G. Widmer, “Cadence Detection in Symbolic Classical Music using Graph Neural Networks,” _Proceedings of the International Society for Music Information Retrieval Conference (ISMIR)_ , 2022.
  * [12] E. Karystinaios, F. Foscarin, and G. Widmer, “Musical Voice Separation as Link Prediction: Modeling a Musical Perception Task as a Multi-Trajectory Tracking Problem,” in _Proceedings of the International Joint Conference on Artificial Intelligence (IJCAI)_ , 2023.
  * [13] C. Hernandez-Olivan, S. R. Llamas, and J. R. Beltran, “Symbolic Music Structure Analysis with Graph Representations and Changepoint Detection Methods,” 2023.
  * [14] Z. Zhang, P. Luo, C. C. Loy, and X. Tang, “Facial Landmark Detection by Deep Multi-task Learning,” in _Proceedings of the European Conference on Computer Vision (ECCV)_ , 2014.
  * [15] M. Guo, A. Haque, D.-A. Huang, S. Yeung, and L. Fei-Fei, “Dynamic Task Prioritization for Multitask Learning,” in _Proceedings of the European Conference on Computer Vision (ECCV)_ , 2018.
  * [16] Z. Chen, V. Badrinarayanan, C.-Y. Lee, and A. Rabinovich, “GradNorm: Gradient Normalization for Adaptive Loss Balancing in Deep Multitask Networks,” in _Proceedings of the International Conference on Machine Learning (ICML)_ , 2018.
  * [17] A. Javaloy and I. Valera, “RotoGrad: Gradient Homogenization in Multitask Learning,” in _Proceedings of the International Conference on Learning Representations (ICLR)_ , 2022.
  * [18] A. Navon, A. Shamsian, I. Achituve, H. Maron, K. Kawaguchi, G. Chechik, and E. Fetaya, “Multi-task Learning as a Bargaining Game,” _Proceedings of the International Conference on Machine Learning (ICML)_ , 2022.
  * [19] W. Hamilton, Z. Ying, and J. Leskovec, “Inductive representation learning on large graphs,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [20] L. Liebel and M. Körner, “Auxiliary tasks in multi-task learning,” _arXiv preprint arXiv:1805.06334_ , 2018.
  * [21] N. Nápoles López, “Automatic roman numeral analysis in symbolic music representations,” Ph.D. dissertation, Schulich School of Music McGill University, December 2022.
  * [22] M. Neuwirth, D. Harasim, F. C. Moss, and M. Rohrmeier, “The Annotated Beethoven Corpus (ABC): A dataset of harmonic analyses of all Beethoven string quartets,” _Frontiers in Digital Humanities_ , vol. 5, p. 16, 2018\.
  * [23] N. Nápoles López, “Automatic Harmonic Analysis of Classical String Quartets from Symbolic Score,” Ph.D. dissertation, Master’s thesis, Universitat Pompeu Fabra, 2017.
  * [24] J. Devaney, C. Arthur, N. Condit-Schultz, and K. Nisula, “Theme and Variation Encodings with Roman Numerals (TAVERN): A new data set for symbolic music analysis,” in _Proceedings of the International Society for Music Information Retrieval Conference (ISMIR)_ , 2015.
  * [25] M. Gotham, D. Tymoczko, and M. S. Cuthbert, “The RomanText Format: A Flexible and Standard Method for Representing Roman Numeral Analyses.” in _Proceedings of the International Society for Music Information Retrieval Conference (ISMIR)_ , 2019.
  * [26] M. R. H. Gotham and P. Jonas, “The Openscore Lieder Corpus,” in _Proceedings of the Music Encoding Conference (MEC)_ , 2021.
  * [27] N. Nápoles López and I. Fujinaga, “Harmonic Reductions as a Strategy for Creative Data Augmentation,” in _Late-Breaking Demo at International Society for Music Information Retrieval Conference (ISMIR)_ , 2020\.
  * [28] J. Hentschel, M. Neuwirth, and M. Rohrmeier, “The Annotated Mozart Sonatas: Score, Harmony, and Cadence,” _Transactions of the International Society for Music Information Retrieval (TISMIR)_ , vol. 4, no. ARTICLE, pp. 67–80, 2021.
  * [29] C. Harte, “Towards automatic extraction of harmony information from music signals,” Ph.D. dissertation, Queen Mary University of London, 2010.


# InfoOT: Information Maximizing Optimal Transport

Ching-Yao Chuang
MIT
<EMAIL_ADDRESS>
Stefanie Jegelka
MIT
<EMAIL_ADDRESS>
David Alvarez-Melis
Microsoft Research
<EMAIL_ADDRESS>

###### Abstract

Optimal transport aligns samples across distributions by minimizing the
transportation cost between them, e.g., the geometric distances. Yet, it
ignores coherence structure in the data such as clusters, does not handle
outliers well, and cannot integrate new data points. To address these
drawbacks, we propose InfoOT, an information-theoretic extension of optimal
transport that maximizes the mutual information between domains while
minimizing geometric distances. The resulting objective can still be
formulated as a (generalized) optimal transport problem, and can be
efficiently solved by projected gradient descent. This formulation yields a
new projection method that is robust to outliers and generalizes to unseen
samples. Empirically, InfoOT improves the quality of alignments across
benchmarks in domain adaptation, cross-domain retrieval, and single-cell
alignment. The code is available at https://github.com/chingyaoc/InfoOT.

##  1 Introduction

Optimal Transport (OT) provides a general framework with a strong theoretical
foundation to compare probability distributions based on the geometry of their
underlying spaces [58]. Besides its fundamental role in mathematics, OT has
increasingly received attention in machine learning due to its wide range of
applications in domain adaptation [22, 47, 59], generative modeling [6, 11],
representation learning [41, 20], and generalization bounds [19]. The
development of efficient algorithms [23, 45] has significantly accelerated the
adoption of optimal transport in these applications.

Computationally, the discrete formulation of OT seeks a matrix, also called
transportation plan, that minimizes the total geometric transportation cost
between two sets of samples drawn from the source and target distributions.
The transportation plan implicitly defines (soft) correspondences across these
samples, but provides no mechanism to relate newly-drawn data points. Aligning
these requires solving a new OT problem from scratch. This limits the
applicability of OT, e.g., to streaming settings where the samples arrive in
sequence, or very large datasets where we can only solve OT on a subset. In
this case, the current solution cannot be used on future data. To overcome
this fundamental constraint, a line of work proposes to directly estimate a
mapping, the pushforward from source to target, that minimizes the
transportation cost [44, 50]. Nevertheless, the resulting mapping is highly
dependent on the complexity of the mapping function [29].

OT could also yield alignments that ignore the intrinsic coherence structure
of the data. In particular, by relying exclusively on pairwise geometric
distances, two nearby source samples could be mapped to disparate target
samples, as in Figure 1, which is undesirable in some settings. For instance,
when applying OT for domain adaptation, source samples with the same class
should ideally be mapped to similar target samples. To mitigate this, prior
work has sought to impose structural priors on the OT objective, e.g., via
submodular cost functions [3] or a Gromov-Wasserstein regularizer [56, 55].
However, these methods still suffer from sensitivity to outliers [40] and
imbalanced data [31, 53].

This work presents _Information Maximization Optimal Transport_ (InfoOT), an
information-theoretic extension of the optimal transport problem that
generalizes the usual formulation by infusing it with _global_ structure in
form of mutual information. In particular, InfoOT seeks alignments that
maximize mutual information, an information-theoretic measure of dependence,
between domains. To do so, we treat the pairs selected by the transportation
plan as samples drawn from the joint distribution and estimate the mutual
information with kernel density estimation based on the paired samples [39].
Interestingly, this results in an OT problem where the cost is the log ratio
between the estimated joint and marginal distributions
$f_{XY}(x,y)/(f_{X}(x)f_{Y}(y))$. Empirically, we show that using a cost
combining mutual information with geometric distances yields better alignments
across different applications. Moreover, akin to Gromov-Wasserstein [37], the
mutual information estimator only relies on intra-domain distances, which
—unlike the standard OT formulation— makes it suitable for aligning
distributions whose supports lie in different metric spaces, e.g., supports
with different modalities or dimensionality [1, 24].

By estimating a joint density, InfoOT naturally yields a novel method for out-
of-sample transportation by taking an expectation over the estimated densities
conditioned on the source samples, which we refer to as _conditional
projection_. Typically, samples are mapped via a barycentric projection [26,
27], which corresponds to the weighted average of target samples, where the
weights are determined by the transportation plan. The barycentric projection
inherits the disadvantages of standard OT: sensitivity to outliers and failing
to generalize to new samples. In contrast, our proposed conditional projection
is robust to outliers and cross-domain class-imbalanced data (Figure 1 and 4)
by averaging over samples with importance sampling, where the weight is,
again, the ratio between the estimated joint and marginal densities.
Furthermore, this projection is well-defined even for unseen samples, which
widens the applicability of OT in streaming or large-scale settings where
solving OT for the complete dataset is prohibitive.

In short, this work makes the following contributions:

  * •

We propose InfoOT, an information-theoretic extension to the optimal transport
that regularizes alignments by maximizing mutual information;

  * •

We develop conditional projection, a new projection method for OT that is
robust to outliers and class imbalance in data, and generalizes to new
samples;

  * •

We evaluate our approach via experiments in domain adaptation, cross-domain
retrieval, and single-cell alignment.

Figure 1: Illustration of InfoOT on 2D point cloud. Compared to classic OT,
InfoOT preserves the cluster structure, where the source points from the same
cluster are mapped to the same target cluster. For projection estimation
(dashed lines), the new conditional projection improves over barycentric
projection with better outlier robustness and out-of-sample generalization.

##  2 Related Works

#### Optimal Transport

Optimal transport provides an elegant framework to compare and align
distributions. The discrete formulation, also called Earth Mover’s Distance
(EMD), finds an optimal coupling between empirical samples by solving a linear
programming problem [10]. To speed up the computation, Cuturi [23] propose the
Sinkhorn distance, an entropic regularized version of EMD that can be solved
more efficiently via the Sinkhorn-Knopp algorithm [33]. Compared to EMD, this
regularized formulation typically yields denser transportation plans, where
samples can be associated with multiple target points. Various extensions of
OT have been proposed to impose stronger priors, e.g., Alvarez-Melis et al.
[3] incorporate additional structure by leveraging a submodular transportation
cost, while Flamary et al. [27] induce class coherence through a group-
sparsity regularizer. The Gromov-Wasserstein (GW) distance [37] is a variant
of OT in which the transportation cost is defined upon intra-domain pairwise
distances. Therefore, GW has been adopted to align ‘incomparable spaces’ [2,
24] as the source and target domains do not need to lie in the same space.
Since the GW objective is no longer a linear program, it is typically
optimized using projected gradient descent [45, 51]. The Fused-GW, which
combines the OT and GW objectives, was proposed by Vayer et al. [55] to
measure graph distances.

#### Mutual Information and OT

The proposed InfoOT extends the standard OT formulation by maximizing a kernel
density estimated mutual information. Recent works [7, 32] also explore the
connection between OT and information theory. Liu et al. [36] consider a semi-
supervised setting for estimating a variant of mutual information, where the
unpaired samples are leveraged to minimize the estimation error. Ozair et al.
[41] replace the KL divergence in mutual information with Wasserstein distance
and develop a loss function for representation learning. In comparison, the
objective of InfoOT is to seek alignments that maximize the mutual information
while being fully unsupervised by parameterizing the joint densities with the
transportation plan. Another line of work also combines OT with kernel density
estimation [14, 38], but focuses on different applications.

##  3 Background on OT and KDE

#### Optimal Transport

Let $\\{x_{i}\\}_{i=1}^{n}\in{\mathcal{X}}^{n}$ and
$\\{y_{i}\\}_{i=1}^{m}\in{\mathcal{Y}}^{m}$ be the empirical samples and
$C\in{\mathbb{R}}^{n\times m}$ be the transportation cost for each pair, e.g,.
Euclidean cost $C_{ij}=\|x_{i}-y_{j}\|$. Given two sets of weights over
samples ${\mathbf{p}}\in{\mathbb{R}}_{+}^{n}$ and
${\mathbf{q}}\in{\mathbb{R}}_{+}^{m}$ where
$\sum_{i=1}^{n}{\mathbf{p}}_{i}=\sum_{i=1}^{m}{\mathbf{q}}_{i}=1$, and a cost
matrix $C$, Kantorovich’s formulation of optimal transport solves

$\displaystyle\min_{\Gamma\in\Pi({\mathbf{p}},{\mathbf{q}})}\langle\Gamma,C\rangle,\quad\Pi({\mathbf{p}},{\mathbf{q}})=\\{\Gamma\in{\mathbb{R}}_{+}^{n\times
m}|\gamma\mathbf{1}_{m}={\mathbf{p}},\gamma^{T}\mathbf{1}_{n}={\mathbf{q}}\\},$

where $\Pi({\mathbf{p}},{\mathbf{q}})$ is a set of transportation plans that
satisfies the flow constraint. In practice, the Sinkhorn distance [23], an
entropic regularized version of OT, can be solved more efficiently via the
Sinkhorn-Knopp algorithm. In particular, the Sinkhorn distance solves
$\min_{\Gamma\in\Pi({\mathbf{p}},{\mathbf{q}})}\langle\Gamma,C\rangle-\epsilon
H(\Gamma)$, where $H(\Gamma)=-\sum_{i,j}\Gamma_{ij}\log\Gamma_{ij}$ is the
entropic regularizer that smooths the transportation plan.

#### Kernel Density Estimation

Kernel Density Estimation (KDE) is a non-parametric density estimation method
based on kernel smoothing [42, 48]. Here, we consider a generalized KDE for
metric spaces $({\mathcal{X}},d_{\mathcal{X}})$ and
$({\mathcal{Y}},d_{\mathcal{Y}})$ [34, 43]. In particular, given a paired
dataset
$\\{x_{i},y_{i}\\}_{i=1}^{n}\in\\{{\mathcal{X}}^{n},{\mathcal{Y}}^{n}\\}$
sampled i.i.d. from an unknown joint density $f_{XY}$ and a kernel function
$K:{\mathbb{R}}\rightarrow{\mathbb{R}}$, KDE estimates the marginals and the
joint density as

$\displaystyle\hat{f}_{X}(x)=\frac{1}{n}\sum_{i}K_{h_{1}}\left(d_{\mathcal{X}}(x,x_{i})\right);\quad\hat{f}_{XY}(x,y)=\frac{1}{n}\sum_{i}K_{h_{1}}\left(d_{\mathcal{X}}(x,x_{i})\right)K_{h_{2}}\left(d_{\mathcal{X}}(y,y_{i})\right),$
(1)

Figure 2: Example of KDE. The bandwidth affect the smoothness.

where $K_{h}(t)=K(\frac{t}{h})/Z_{h}$ and the normalizing constant $Z_{h}$
makes equation 1 integrate to one. The bandwidth parameter $h$ controls the
smoothness of the estimated densities. Figure 2 illustrates an example of KDE
on 1D data. In this work, we do not need to estimate the normalizing constant
as only the ratio between joint and marginal densities
$\hat{f}_{XY}(x,y)/(\hat{f}_{X}(x)\hat{f}_{Y}(y))$ is considered while
estimating the mutual information. For all the presented experiments, we adopt
the Gaussian kernel:

$\displaystyle
K_{h}\left(d_{\mathcal{X}}(x,x^{\prime})\right)=\frac{1}{Z_{h}}\exp\left(-\frac{d_{\mathcal{X}}(x,x^{\prime})^{2}}{2h^{2}\sigma^{2}}\right),$

where $\sigma^{2}$ controls the variance. The Gaussian kernel has been
successfully adopted for KDE beyond the Euclidean space [34, 49], and we found
it to work well in our experiments. For simplicity, we also set
$h_{1}=h_{2}=h$ for all the experiments.

##  4 Information Maximizing OT

Optimal transport captures the geometry of the underlying space through the
ground metric in its objective. Additional information is not directly
captured in this metric —such as coherence structure— will therefore be
ignored when solving the problem. This is undesirable in applications where
this additional structure matters, for instance in domain adaptation, where
class coherence is crucial. As a concrete example, the cluster structure in
the dataset in Figure 1 is ignored by classic OT. Intuitively, the reason for
this issue is that the classic OT is too _local_ : the transportation cost
considers each sample separately, without respecting coherence across close-by
samples. Next, we show that mutual information estimated with KDE can
introduce global structure into OT maps.

###  4.1 Measuring Global Structure with Mutual Information

Formally, mutual information measures the statistical dependence of two random
variables $X,Y$:

$\displaystyle
I(X,Y)=\int_{\mathcal{Y}}\int_{\mathcal{X}}f_{X,Y}(x,y)\log\left(\frac{f_{XY}(x,y)}{f_{X}(x)f_{Y}(y)}\right)dxdy$

where $f_{XY}$ is joint density and $f_{X},f_{Y}$ are marginal probability
density functions. For paired datasets, various mutual information estimators
have been defined [8, 39, 46]. In contrast, we are interested in the inverse:
given _unpaired_ samples $\\{x_{i}\\}_{i=1}^{n},\\{y_{j}\\}_{i=1}^{m}$, can we
find alignments that maximize the mutual information?

#### Discrete v.s. Continuous.

An immediate idea is to treat the discrete transportation plan $\Gamma$ as the
joint distribution between $X$ and $Y$, and write the mutual information as
$\sum_{i,j}\Gamma_{ij}\log(nm\Gamma_{ij})=\log(nm)-H(\Gamma)$. In this case,
maximizing mutual information would be equivalent to minimizing the entropic
regularizer $H(\Gamma)$ introduced by [23]. For a finite set of samples, this
mutual information estimator is trivially maximized for _any_ one-to-one
mapping as then $H(\Gamma)=0$. Figure 3 (a) illustrates two one-to-one
mappings $\Gamma_{a}$ and $\Gamma_{b}$ between points sampled from multi-mode
Gaussian distributions, where $\Gamma_{a}$ preserves the cluster structure and
$\Gamma_{b}$ is simply a random permutation. They both maximize the mutual
information estimate above, yet $\Gamma_{a}$ is a better alignment with high
coherence. In short, directly using the transportation plan estimated from
finite samples as the joint distribution to estimate mutual information
between continuous random variables is problematic. In contrast, joint
distributions estimated with KDE tend to be smoother, such as $\Gamma_{a}$ in
Figure 3 (b). This suggests that KDE may lead to a better objective for the
alignment problem.

Figure 3: Measuring Structure with Mutual information. (a) The $\Gamma_{a}$
and $\Gamma_{b}$ are two one-to-one mappings where $\Gamma_{a}$ preserves the
cluster structure and $\Gamma_{b}$ is a random permutation; (b) The estimated
joint density of $\Gamma_{a}$ is more concentrated than the one of
$\Gamma_{b}$, which also leads to higher mutual information under KDE.

###  4.2 InfoOT: Maximizing Mutual Information with KDE

Instead of directly interpreting the OT plan as the joint distribution for the
mutual information, we use it to inform the definition of a different one. In
particular, we treat $\Gamma_{ij}$ as _the weight of pair $(x_{i},y_{j})$
within the empirical samples drawn from the unknown joint distribution with
density $f_{XY}$_. Intuitively, $\Gamma_{ij}$ defines what empirical samples
we obtain by sampling from the joint distribution. Given a transportation plan
$\Gamma$, the kernelized joint density in equation 1 can be rewritten as

$\displaystyle\hat{f}_{\Gamma}(x,y)$
$\displaystyle=\sum_{i}\sum_{j}\Gamma_{ij}K_{h}\left(d_{\mathcal{X}}(x,x_{i})\right)K_{h}\left(d_{\mathcal{Y}}(y,y_{j})\right).$
(2)

The $1/n$ factor is dropped as the plan $\Gamma$ is already normalized
($\sum_{ij}\Gamma_{ij}=1$). Specifically, we replace the prespecified paired
samples in equation 1 with the ones selected by the plan $\Gamma$.

######  Definition 1 (Kernelized Mutual Information).

The KDE estimated mutual information reads

$\displaystyle\hat{I}_{\Gamma}(X,Y)$
$\displaystyle=\sum_{i,j}\Gamma_{ij}\log\frac{\hat{f}_{\Gamma}(x_{i},y_{j})}{\hat{f}(x_{i})\hat{f}(y_{j})}=\sum_{i,j}\Gamma_{ij}\log\frac{nm\cdot\sum_{k,l}\Gamma_{kl}K_{h}\left(d_{\mathcal{X}}(x_{i},x_{k})\right)K_{h}\left(d_{\mathcal{Y}}(y_{j},y_{l})\right)}{\sum_{k}K_{h}\left(d_{\mathcal{X}}(x_{i},x_{k})\right)\cdot\sum_{l}K_{h}\left(d_{\mathcal{Y}}(y_{j},y_{l})\right)}.$

The normalizing constant $Z_{h}$ in equation (1) cancels out while calculating
the ratio between joint and marginal probability densities. To maximize the
empirical mutual information $\hat{I}_{\Gamma}(X,Y)$, the plan has to map
close-by points $i,k$ to close-by points $j,l$. Maximizing this information
can be interpreted as an optimal transport problem:

$\displaystyle\textnormal{(InfoOT)
}\quad\max_{\Gamma\in\Pi({\mathbf{p}},{\mathbf{q}})}\hat{I}_{\Gamma}(X,Y)=\min_{\Gamma\in\Pi({\mathbf{p}},{\mathbf{q}})}\sum_{i,j}\Gamma_{ij}\cdot\log\left(\frac{\hat{f}(x_{i})\hat{f}(y_{j})}{\hat{f}_{\Gamma}(x_{i},y_{j})}\right).$
(3)

Instead of pairwise (Euclidean) distances, the transportation cost is now the
log ratio between the estimated marginal and joint densities. The following
lemma illustrates the asymptotic relation between the kernel estimated mutual
information and the entropic regularizer.

######  Lemma 2.

When $h\rightarrow 0$ and $K(\cdot)$ is the Gaussian kernel, we have
$\hat{I}_{\Gamma}(X,Y)\rightarrow-H(\Gamma)+\log(nm)$.

When the bandwidth $h$ goes to zero, the estimated density is the sum of delta
functions centered at the samples, and the estimated mutual information
degenerates back to the standard entropic regularizer [23].

Note that the formulation of InfoOT does not require the support of $X$ and
$Y$ to be comparable. Similar to Gromov-Wasserstein [37], InfoOT only relies
on intra-domain distances, which makes it an appropriate objective for
aligning distributions when the supports do not lie in the same metric space,
e.g., supports with different modalities or dimensionalities, as section 6.4
shows.

#### Fused InfoOT: Incorporating the Geometry.

When the geometry between domains is informative, the mutual information can
act as a regularizer that refines the alignment. Along with a weighting
parameter $\lambda$, we define the Fused InfoOT as

$\displaystyle\textnormal{(F-InfoOT)
}\;\min_{\Gamma\in\Pi({\mathbf{p}},{\mathbf{q}})}\left\langle\Gamma,C\right\rangle-\lambda\hat{I}_{\Gamma}(X,Y)=\min_{\Gamma\in\Pi({\mathbf{p}},{\mathbf{q}})}\sum_{i,j}\Gamma_{ij}\cdot\left(C_{ij}+\lambda\cdot\log\left(\frac{\hat{f}(x_{i})\hat{f}(y_{j})}{\hat{f}_{\Gamma}(x_{i},y_{j})}\right)\right).$

The transportation cost becomes the weighted sum between the pairwise
distances $C$ and the log ratio of joint and marginals densities. As Figure 1
illustrates, the mutual information regularizer excludes alignments that
destroy the cluster structure while minimizing the pairwise distances.
Practically, we found F-InfoOT suitable for general OT applications such as
unsupervised domain adaptation [27] and color transfer [26] where the geometry
between source and target is informative.

###  4.3 Numerical Optimization

As the transportation cost is dependent on $\Gamma$, the objective is no
longer linear in $\Gamma$ and cannot be solved with linear programming.
Instead, we adopt the projected gradient descent introduced in [45]. In
particular, Benamou et al. [9] show that the projection can be done by simply
solving the Sinkhorn distance [23] if the non-linear objective is augmented
with the entropic regularizer $H(\Gamma)$. For instance, we can augment
F-InfoOT as follows:

$\displaystyle\min_{\Gamma\in\Pi({\mathbf{p}},{\mathbf{q}})}\left\langle\Gamma,C\right\rangle-\lambda\hat{I}_{\Gamma}(X,Y)-\epsilon
H(\Gamma).$

In this case, the update of projected gradient descent reads

$\displaystyle\Gamma_{t+1}\leftarrow\operatorname*{arg\,min}_{\Gamma\in\Pi({\mathbf{p}},{\mathbf{q}})}\left\langle\Gamma,C-\lambda\nabla_{\Gamma}\hat{I}_{\Gamma_{t}}(X,Y)\right\rangle-\epsilon
H(\Gamma).$ (4)

The update is done by solving the sinkhorn distance [23], where the cost
function is the gradient to the objective of F-InfoOT. We provide a detailed
derivation of (4) in Appendix A.2.

#### Matrix Computation

Practically, the optimization can be efficiently computed with matrix
multiplications. The gradient with respect to the transportation plan $\Gamma$
is

$\displaystyle\frac{\partial\hat{I}_{\Gamma}(X,Y)}{\partial\Gamma_{ij}}$
$\displaystyle=\log\left(\frac{\hat{f}_{\Gamma}(x_{i},y_{j})}{\hat{f}(x_{i})\hat{f}(y_{j})}\right)+\sum_{k,l}\Gamma_{kl}\frac{K_{h}\left(d_{\mathcal{X}}(x_{i},x_{k})\right)K_{h}\left(d_{\mathcal{Y}}(y_{j},y_{l})\right)}{\hat{f}_{\Gamma}(x_{k},y_{l})}.$

Let $K_{X}$ and $K_{Y}$ be the kernel matrices where
$(K_{X})_{ij}=K_{h}\left(d_{\mathcal{X}}(x_{i}-x_{j})\right)$,
$(K_{Y})_{ij}=K_{h}\left((d_{\mathcal{Y}}(y_{i}-y_{j})\right)$. The gradient
has the following matrix form:

$\displaystyle\nabla_{\Gamma}\hat{I}_{\Gamma}(X,Y)=\log\left(K_{X}\Gamma
K_{Y}^{T}\oslash M_{X}M_{Y}^{T}\right)+K_{X}\left(\Gamma\oslash K_{X}\Gamma
K_{Y}^{T}\right)K_{Y}^{T}$

where $(M_{X})_{i}=\hat{f}(x_{i})$, $(M_{Y})_{i}=\hat{f}(y_{i})$ are the
marginal density vectors and $\oslash$ denotes element-wise division. The
gradient can be computed with matrix multiplications in
${\mathcal{O}}(n^{2}m+nm^{2})$.

##  5 Conditional Projection with InfoOT

Many applications of optimal transport involve mapping source points to a
target domain. For instance, when applying OT for domain adaptation, the
classifiers are trained on projected source samples that are mapped to the
target domain. When ${\mathcal{X}}={\mathcal{Y}}$, given a transportation plan
$\Gamma$, a _barycentric projection_ maps source samples to the target domain
by minimizing the weighted cost to target samples [27, 44]. The mapping is
equivalent to the weighted average of the target samples when the cost
function is the squared Euclidean distance $c(x,y)=\|x-y\|^{2}$:

$\displaystyle
x_{i}\mapsto\operatorname*{arg\,min}_{y\in{\mathcal{Y}}}\sum_{j=1}^{m}\Gamma_{ij}\|y-y_{j}\|^{2}=\frac{1}{\sum_{j=1}^{m}\Gamma_{ij}}\sum_{j=1}^{m}\Gamma_{ij}y_{j}.$
(5)

Despite its simplicity, the barycentric projection fails when (a) aligning
data with outliers, (b) imbalanced data, and (c) mapping new samples. For
instance, if sample $x_{i}$ is mostly mapped to an outlier $y_{j}$, then its
projection will be close to $y_{j}$. Similar problems occur when applying OT
for domain adaptation. If the size of a same class differs between domains,
false alignments would emerge due to the flow constraint of OT as Figure 4
illustrates, which worsen the subsequent projections.

Since the barycentric projection relies on the transportation plan to
calculate the weights, any new source sample requires re-computing OT to
obtain the transportation plan for it. This can be computationally prohibitive
in large-scale settings. In the next section, we show that the densities
estimated via InfoOT can be leveraged to compute the conditional expectation,
which leads to a new mapping approach that is both robust and generalizable.

###  5.1 Conditional Expectation via KDE

When treating the transportation plan as a probability mass function in the
right-hand side of (equation 5), the barycentric projection resembles the
conditional expectation $\mathbb{E}[Y|X=x]$. Indeed, the classic definition of
barycentric projection [5] is defined as the integral over the conditional
distribution. But, this again faces the issues discussed in section 4.
Instead, equipped with the densities estimated via KDE and InfoOT, the
conditional expectation can be better estimated with classical Monte-Carlo
importance sampling using samples from the marginal $P_{Y}$:

$\displaystyle x\mapsto\underset{{P_{Y|X=x}}}{\mathbb{E}}[y]=\underset{y\sim
P_{Y}}{\mathbb{E}}\left[\frac{f_{Y|X=x}(y)}{f_{Y}(y)}y\right]=\underset{y\sim
P_{Y}}{\mathbb{E}}\left[\frac{f_{XY}(x,y)}{f_{X}(x)f_{Y}(y)}y\right]\approx\frac{1}{Z}\frac{\hat{f}_{\Gamma}(x,y_{j})}{\hat{f}_{X}(x)\hat{f}_{Y}(y_{j})}y_{j}$
(6)

where
$Z=\sum_{j=1}^{m}\hat{f}_{\Gamma}(x,y_{j})/(\hat{f}_{X}(x)\hat{f}_{Y}(y_{j}))$
is the normalizing constant. Compared to the barycentric projection, the
_importance weight_ for each $y_{j}$ is the ratio between the joint and the
marginal densities. To distinguish the KDE conditional expectation with
barycentric projection, we will refer to the proposed mapping as _conditional
projection_.

#### Robustness against Noisy Data.

By definition in equation 2, the joint density $\hat{f}_{\Gamma}(x,y)$
measures the similarity of $(x,y)$ to all other pairs selected by the
transportation plan $\Gamma$. Even if $x$ is aligned with outliers or wrong
clusters, as long as the points near $x$ are mostly aligned with the correct
samples, the conditional projection will project $x$ to similar samples as
they are upweighted by the joint density in (6). This makes the mapping much
less sensitive to outliers and imbalanced datasets. See Figure 1 and Figure 4
for illustrations.

#### Out-of-sample Mapping.

The conditional projection is well-defined for any $x\in{\mathcal{X}}$, and
naturally generalizes to _new samples_ without recomputing the OT.
Importantly, the importance weight
$\hat{f}_{\Gamma}(x,y)/(\hat{f}_{X}(x)\hat{f}_{Y}(y))$ can be interpreted as a
similarity score between $(x,y)$, which is useful for retrieval tasks as
section 6.3 shows.

The conditional projection tends to cluster points together with larger
bandwidths that lead to more averaging. We found that using a smaller
bandwidth (e.g., $h=0.1$) for the conditional projection improves the
diversity of the projection when the dataset is less noisy, e.g., the data in
Figure 1. For noisy or imbalanced datasets, the same bandwidth used for
optimizing InfoOT works well. Note that analogous to Lemma 2, when the
bandwidth $h\rightarrow 0$, the conditional projection converges to the
barycentric projection, making the barycentric projection a special case of
the conditional projection (Figure 4).

##  6 Experiments

We now evaluate InfoOT with experiments in point cloud matching, domain
adaptation, cross-domain retrieval, and single-cell alignment. All the optimal
transport approaches are implemented or adopted from the POT library [28].
Detailed experimental settings and additional experiments can be found in the
appendix.

Figure 4: Projection under imbalanced samples. When the cluster sizes mismatch
between source and target, barycentric projection wrongly projects samples to
the incorrect cluster. In contrast, increasing the bandwidth of conditional
projection gradually improves the robustness and yields better projection.

###  6.1 Point Cloud Matching

We begin with a 2D toy example, where both source and target samples are drawn
from a Gaussian distribution with 2 modes, but the latter is rotated and has
two outliers added to it, as Figure 1 shows. We compare the behavior of
different variants of OT and mappings on this data. Perhaps not surprisingly,
standard OT maps the source points in the same cluster to two different target
clusters, overlooking the intrinsic structure of the data. In comparison, the
alignment of InfoOT retains the cluster structure. On the right hand side, the
barycentric projection maps two source points wrongly to the target outliers,
while the conditional projection is not affected by the outliers. Lastly, we
demonstrate an out-of-sample mapping with the conditional projection, where
newly sampled points are correctly mapped to clusters.

Figure 4 depicts an class-imbalanced setting, where the corresponding clusters
in source and target have different numbers of samples. Therefore, the
barycentric projection wrongly maps samples from the same source cluster to
different target clusters. When increasing the bandwidth in the conditional
projection, the smoothing effect of KDE gradually corrects the mapping and
yields more concentrated projections. In appendix B.1, we further demonstrate
that InfoOT improves the baselines in a color transfer task, where pixels are
treated as points in RGB space.

Figure 5: tSNE visualization of projections. We show the t-SNE visualization of projectrd source samples (circles) along with the target samples (triangles) on A$\rightarrow$C. Classes are indicated by colors. | OT | Sinkhorn | GL-OT | FGW | Linear | F-InfoOT |  F-InfoOT∗
---|---|---|---|---|---|---|---
MNIST-USPS
M$\rightarrow$U |  46.6$\pm$1.2 |  62.7$\pm$2.0 |  63.0$\pm$2.0 |  49.6$\pm$7.3 |  60.8$\pm$2.3 | 69.9$\pm$3.1 |  61.1$\pm$1.9
U$\rightarrow$M |  48.1$\pm$1.1 |  58.6$\pm$0.9 |  58.8$\pm$1.0 |  37.8$\pm$5.3 |  59.5$\pm$1.0 | 65.1$\pm$1.4 |  57.0$\pm$1.7
Office-Caltech
C$\rightarrow$D |  61.3$\pm$10.9 |  71.9$\pm$15.9 |  76.9$\pm$18.6 |  61.3$\pm$10.9 |  58.1$\pm$12.5 |  78.1$\pm$13.9 | 87.5$\pm$7.8
C$\rightarrow$W |  64.0$\pm$7.0 |  66.0$\pm$9.3 |  66.7$\pm$9.4 |  64.0$\pm$7.0 |  62.0$\pm$7.4 |  79.7$\pm$4.3 | 81.0$\pm$6.7
C$\rightarrow$A |  78.6$\pm$4.7 |  77.3$\pm$4.8 |  83.3$\pm$6.1 |  79.0$\pm$4.5 |  77.9$\pm$6.7 |  87.0$\pm$3.7 | 90.6$\pm$2.0
D$\rightarrow$W |  90.7$\pm$3.8 |  90.7$\pm$4.7 | 93.7$\pm$4.3 |  90.7$\pm$3.8 |  89.0$\pm$6.3 |  91.0$\pm$4.1 |  93.3$\pm$4.7
D$\rightarrow$A |  73.8$\pm$4.5 |  73.4$\pm$2.9 |  84.4$\pm$2.9 |  73.8$\pm$4.5 |  72.7$\pm$3.9 |  83.6$\pm$4.3 | 89.8$\pm$1.8
D$\rightarrow$C |  67.0$\pm$3.0 |  66.4$\pm$3.8 |  76.5$\pm$2.9 |  67.0$\pm$3.0 |  67.5$\pm$4.1 |  70.2$\pm$3.7 | 80.8$\pm$1.8
W$\rightarrow$D |  81.3$\pm$7.8 |  80.6$\pm$10.4 |  85.6$\pm$10.6 |  81.3$\pm$7.8 |  79.4$\pm$7.2 |  79.4$\pm$8.9 | 89.4$\pm$11.8
W$\rightarrow$C |  64.8$\pm$4.6 |  65.8$\pm$4.1 |  73.5$\pm$5.4 |  64.8$\pm$4.6 |  65.5$\pm$4.5 | 75.8$\pm$3.3 |  74.4$\pm$3.7
W$\rightarrow$A |  67.3$\pm$4.9 |  69.3$\pm$5.4 |  79.6$\pm$3.0 |  67.3$\pm$4.9 |  70.5$\pm$4.8 |  85.5$\pm$1.9 | 89.3$\pm$2.3
A$\rightarrow$D |  73.1$\pm$10.6 |  68.8$\pm$8.8 |  76.3$\pm$8.2 |  73.1$\pm$10.6 |  66.9$\pm$7.8 |  80.6$\pm$7.5 | 81.3$\pm$9.8
A$\rightarrow$W |  64.7$\pm$6.3 |  70.0$\pm$7.5 |  70.0$\pm$7.0 |  64.7$\pm$6.3 |  69.3$\pm$6.6 |  82.0$\pm$7.2 | 87.0$\pm$4.6
A$\rightarrow$C |  65.4$\pm$5.3 |  69.3$\pm$6.0 |  79.8$\pm$5.8 |  65.5$\pm$5.7 |  66.9$\pm$4.5 |  74.4$\pm$3.4 | 81.2$\pm$3.6
AVG |  71.0$\pm$2.5 |  72.4$\pm$3.7 |  78.9$\pm$4.5 |  71.0$\pm$2.5 |  70.5$\pm$2.4 |  80.6$\pm$3.3 | 85.6$\pm$3.3
Table 1: Optimal Transport for Domain Adaptation. The Fused-InfoOT with conditional projection (F-InfoOT∗) performs significantly better than the barycentric counterpart (F-InfoOT) and the other baselines when the dataset exhibit class imbalance, e.g., Office-Caltech.  | Office-Caltech | ImageCLEF
---|---|---
| P@1 | P@5 | P@15 | P@1 | P@5 | P@15
L2-NN |  70.0$\pm$13.9 |  62.9$\pm$14.0 |  53.6$\pm$12.1 |  80.4$\pm$10.5 |  77.6$\pm$9.5 |  71.7$\pm$9.3
Sinkhorn+NN |  69.0$\pm$11.6 |  65.1$\pm$6.4 |  58.0$\pm$6.9 | 81.9$\pm$10.2 |  81.1$\pm$10.6 |  79.2$\pm$10.4
FGW+NN |  69.6$\pm$11.7 |  65.6$\pm$6.6 |  58.5$\pm$7.0 |  81.9$\pm$10.4 |  81.2$\pm$10.6 |  79.3$\pm$10.3
F-InfoOT | 76.4$\pm$9.0 | 75.1$\pm$9.1 | 70.4$\pm$10.6 |  81.2$\pm$9.3 | 81.9$\pm$9.8 | 80.0$\pm$10.7
Table 2: Optimal Transport for Cross-Domain Retrieval. With conditional
projection, InfoOT is capable to perform alignment for unseen samples without
any modification.

###  6.2 Domain Adaptation

Next, we apply the fused version of InfoOT to two domain adaptation
benchmarks: MNIST-USPS and the Office-Caltech dataset [30]. The MNIST (M) and
USPS (U) are digit classification datasets, and the Office-Caltech dataset
contains 4 domains: Amazon (A), Dslr (D), Webcam (W) and Caltech10 (C), with
images labeled as one of 10 classes. For MNIST and USPS, the raw images are
directly used to compute the distances, while we adopt decaf6 features [25]
extracted from pretrained neural networks for Office-Caltech. Following
previous works on OT for domain adaptation [3, 27, 44], the source samples are
first mapped to the target, and 1-NN classifiers are then trained on the
projected samples with source labels. The barycentric projection is adopted
for all the baselines, while F-InfoOT is tested with both barycentric and
conditional projection.

Following [27], we present the results over 10 independent trials. In each
trial of Office-Caltech, the target data is divided into 90%/10% train-test
split, where OT and 1-NN classifiers are only computed on the training set.
For MNIST-USPS, only 2000 samples from the source and target training set are
used, while the original test sets are used. The strength of the entropy
regularizer $\epsilon$ is set to $1$ for every entropic regularized OT, and
the $\lambda$ of F-InfoOT is set to $100$ for all the experiments. The
bandwidth for each benchmark is selected from $\\{0.2,0.3,...,0.8\\}$ with the
_circular validation procedure_ [12, 44, 61] on M$\rightarrow$U and
A$\rightarrow$D, which is $0.4$ and $0.5$, respectively. We compare F-InfoOT
with Sinkhorn distance [23], group-lasso regularized OT [27], fused Gromov-
Wasserstein (FGW) [57], and linear-mapping OT [44]. For OTs involving intra-
domain distances such as F-InfoOT and FGW, we adopt the following class-
conditional distance for the source:
$\|x_{i}-x_{j}\|+5000\cdot\mathbb{1}_{f(x_{i})\neq f(x_{j})}$, where the
second term is a penalty on class mismatch [1, 60] and $f$ is the labeling
function. As Table 1 shows, F-InfoOT with barycentric projection outperforms
the baselines in both benchmarks, demonstrating that mutual information
captures the intrinsic structure of the datasets. In Office-Caltech, many
datasets exhibit the class-imbalance problem, which makes F-InfoOT with
conditional projection significantly outperform the barycentric projection and
the other baselines. Figure 5 visualizes the projected source and target
samples with tSNE [54]. The barycentric projection tends to produce one-to-one
alignments, which suffer from class-imbalanced data. In contrast, conditional
projection yields concentrated projections that preserves the class structure.

| scGEM | SNAREseq
---|---|---
| FOS, | Acc, | FOS. | Acc..
UnionCom | 0.210 | 58.2 | 0.265 | 42.3
MMD-MA | 0.201 | 58.8 | 0.150 | 94.2
SCOT (GW) | 0.190 | 57.6 | 0.150 | 98.2
InfoOT | 0.178 | 68.9 | 0.156 | 98.8

Table 3: Single-Cell Alignment Performance. Similar to GW, InfoOT also
performs well in cross-domain alignment.

Figure 6: scGEM alignment with InfoOT. The cell types are indicated by colors.

###  6.3 Cross-domain Retrieval

We now consider unsupervised cross-domain image retrieval, where given a
source sample, the algorithms have to determine the top-k similar target
samples. Given fixed source and target samples, this can be formulated as an
optimal transport problem, where the transportation plan $\Gamma_{ij}$ gives
the similarity score between the candidate source sample $x_{i}$ and target
samples $y_{j}$. Nevertheless, this formulation fails when new source samples
come. For standard OT, one has to solve the OT problem again to obtain the
alignment for new samples. In contrast, the importance weight
$\hat{f}_{\Gamma}(x,y)/(\hat{f}_{X}(x)\hat{f}_{Y}(y))$ defined in conditional
projection (6) naturally provides the similarity score between the candidate
$x$ and each target sample $y$. We test F-InfoOT on the Office-Caltech [30]
and ImageClef datasets [16], where we adopt the same hyperparameter for
Office-Caltech from the previous section. In the unsupervised setting, the in-
domain transportation cost for the source is the Euclidean distance instead of
the class-conditional distance. To compare with standard OTs, we adopt a
nearest neighbor approach for the baselines: (1) retrieve the nearest source
sample given an unseen sample, and (2) use the transportation plan of the
nearest source sample to retrieve target samples. Along with a simple nearest
neighbor retrieval baseline (L2-NN), the average top-k precision over 10
trials is shown in Table 2. The fuesed InfoOT significantly outperforms the
baselines on Office-Caltech across different choices of $k$.

###  6.4 Single Cell Alignment

Finally, we examine InfoOT in unsupervised alignment between incomparable
spaces with the single-cell multi-omics dataset from [24]. Recent techniques
allow to obtain different cellular features at the single-cell resolution [13,
17, 52]. Nevertheless, different features are typically collected from
different sets of cells, and aligning them is crucial for unified data
analysis. We examine InfoOT with the sc-GEM [18] and SNARE-seq [17] dataset
provided by [24] and follow the same data preprocessing steps, distance
calculation, and evaluation setup. Here, two evaluation metrics are
considered: “fraction of samples closer than the true match” (FOSCTTM) [35]
and the label transfer accuracy [15]. We compare InfoOT with UnionCom [15],
MMD-MA [35], and SCOT [24], where SCOT is an optimal transport baseline with
Gromov-Wasserstein distance. Similarly, the bandwidth for each dataset is
selected from $\\{0.2,0.3,...,0.8\\}$ with the _circular validation
procedure_. As Table 3 shows, InfoOT significantly improves the baselines on
the sc-GEM dataset, while being comparable on the SNARE-seq dataset,
demonstrating the applicability of InfoOT on cross-domain alignment. Figure 6
further visualizes the barycentric projection with InfoOT, where we can see
that cells with the same type are well aligned.

##  7 Conclusion

In this work, we propose InfoOT, an information-theoretic extension of optimal
transport. InfoOT produces smoother, coherent alignments by maximizing the
mutual information estimated with KDE. InfoOT leads to a new mapping method,
conditional projection, that is robust to class imbalance and generalizes to
unseen samples. We extensively demonstrate the applicability of InfoOT across
benchmarks in different modalities.

#### Acknowledgements

We thank Nicolo Fusi for useful conversations and feedback. This work was in
part supported by NSF BIGDATA IIS-1741341, NSF AI Institute TILOS, NSF CAREER
1553284. CC is supported by a IBM PhD Fellowship.

## References

  * Alvarez-Melis & Fusi [2020] David Alvarez-Melis and Nicolo Fusi.  Geometric dataset distances via optimal transport.  _Advances in Neural Information Processing Systems_ , 33:21428–21439, 2020.
  * Alvarez-Melis & Jaakkola [2018] David Alvarez-Melis and Tommi Jaakkola.  Gromov-wasserstein alignment of word embedding spaces.  In _Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing_ , pp. 1881–1890, 2018.
  * Alvarez-Melis et al. [2018] David Alvarez-Melis, Tommi Jaakkola, and Stefanie Jegelka.  Structured optimal transport.  In _International Conference on Artificial Intelligence and Statistics_ , pp. 1771–1780. PMLR, 2018.
  * Alvarez-Melis et al. [2019] David Alvarez-Melis, Stefanie Jegelka, and Tommi S Jaakkola.  Towards optimal transport with global invariances.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pp. 1870–1879. PMLR, 2019.
  * Ambrosio et al. [2005] Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré.  _Gradient flows: in metric spaces and in the space of probability measures_.  Springer Science & Business Media, 2005.
  * Arjovsky et al. [2017] Martin Arjovsky, Soumith Chintala, and Léon Bottou.  Wasserstein generative adversarial networks.  In _International conference on machine learning_ , pp. 214–223. PMLR, 2017.
  * Bai et al. [2020] Yikun Bai, Xiugang Wu, and Ayfer Özgür.  Information constrained optimal transport: From talagrand, to marton, to cover.  In _2020 IEEE International Symposium on Information Theory (ISIT)_ , pp. 2210–2215. IEEE, 2020.
  * Belghazi et al. [2018] Mohamed Ishmael Belghazi, Aristide Baratin, Sai Rajeshwar, Sherjil Ozair, Yoshua Bengio, Aaron Courville, and Devon Hjelm.  Mutual information neural estimation.  In _International conference on machine learning_ , pp. 531–540. PMLR, 2018.
  * Benamou et al. [2015] Jean-David Benamou, Guillaume Carlier, Marco Cuturi, Luca Nenna, and Gabriel Peyré.  Iterative bregman projections for regularized transportation problems.  _SIAM Journal on Scientific Computing_ , 37(2):A1111–A1138, 2015.
  * Bonneel et al. [2011] Nicolas Bonneel, Michiel Van De Panne, Sylvain Paris, and Wolfgang Heidrich.  Displacement interpolation using lagrangian mass transport.  In _Proceedings of the 2011 SIGGRAPH Asia conference_ , pp. 1–12, 2011.
  * Bousquet et al. [2017] Olivier Bousquet, Sylvain Gelly, Ilya Tolstikhin, Carl-Johann Simon-Gabriel, and Bernhard Schoelkopf.  From optimal transport to generative modeling: the vegan cookbook.  _arXiv preprint arXiv:1705.07642_ , 2017.
  * Bruzzone & Marconcini [2009] Lorenzo Bruzzone and Mattia Marconcini.  Domain adaptation problems: A dasvm classification technique and a circular validation strategy.  _IEEE transactions on pattern analysis and machine intelligence_ , 32(5):770–787, 2009.
  * Buenrostro et al. [2015] Jason D Buenrostro, Beijing Wu, Howard Y Chang, and William J Greenleaf.  Atac-seq: a method for assaying chromatin accessibility genome-wide.  _Current protocols in molecular biology_ , 109(1):21–29, 2015.
  * Canas & Rosasco [2012] Guillermo Canas and Lorenzo Rosasco.  Learning probability measures with respect to optimal transport metrics.  _Advances in Neural Information Processing Systems_ , 25, 2012.
  * Cao et al. [2020] Kai Cao, Xiangqi Bai, Yiguang Hong, and Lin Wan.  Unsupervised topological alignment for single-cell multi-omics integration.  _Bioinformatics_ , 36(Supplement_1):i48–i56, 2020.
  * Caputo et al. [2014] Barbara Caputo, Henning Müller, Jesus Martinez-Gomez, Mauricio Villegas, Burak Acar, Novi Patricia, Neda Marvasti, Suzan Üsküdarlı, Roberto Paredes, Miguel Cazorla, et al.  Imageclef 2014: Overview and analysis of the results.  In _International Conference of the Cross-Language Evaluation Forum for European Languages_ , pp. 192–211. Springer, 2014.
  * Chen et al. [2019] Song Chen, Blue B Lake, and Kun Zhang.  High-throughput sequencing of the transcriptome and chromatin accessibility in the same cell.  _Nature biotechnology_ , 37(12):1452–1457, 2019\.
  * Cheow et al. [2016] Lih Feng Cheow, Elise T Courtois, Yuliana Tan, Ramya Viswanathan, Qiaorui Xing, Rui Zhen Tan, Daniel SW Tan, Paul Robson, Yuin-Han Loh, Stephen R Quake, et al.  Single-cell multimodal profiling reveals cellular epigenetic heterogeneity.  _Nature methods_ , 13(10):833–836, 2016.
  * Chuang et al. [2021] Ching-Yao Chuang, Youssef Mroueh, Kristjan Greenewald, Antonio Torralba, and Stefanie Jegelka.  Measuring generalization with optimal transport.  _Advances in Neural Information Processing Systems_ , 34:8294–8306, 2021.
  * Chuang et al. [2022] Ching-Yao Chuang, R Devon Hjelm, Xin Wang, Vibhav Vineet, Neel Joshi, Antonio Torralba, Stefanie Jegelka, and Yale Song.  Robust contrastive learning against noisy views.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 16670–16681, 2022.
  * Conneau et al. [2017] Alexis Conneau, Guillaume Lample, Marc’Aurelio Ranzato, Ludovic Denoyer, and Hervé Jégou.  Word translation without parallel data.  _arXiv preprint arXiv:1710.04087_ , 2017.
  * Courty et al. [2017] Nicolas Courty, Rémi Flamary, Amaury Habrard, and Alain Rakotomamonjy.  Joint distribution optimal transportation for domain adaptation.  _Advances in Neural Information Processing Systems_ , 30, 2017.
  * Cuturi [2013] Marco Cuturi.  Sinkhorn distances: Lightspeed computation of optimal transport.  _Advances in neural information processing systems_ , 26, 2013.
  * Demetci et al. [2020] Pinar Demetci, Rebecca Santorella, Björn Sandstede, William Stafford Noble, and Ritambhara Singh.  Gromov-wasserstein optimal transport to align single-cell multi-omics data.  _BioRxiv_ , 2020.
  * Donahue et al. [2014] Jeff Donahue, Yangqing Jia, Oriol Vinyals, Judy Hoffman, Ning Zhang, Eric Tzeng, and Trevor Darrell.  Decaf: A deep convolutional activation feature for generic visual recognition.  In _International conference on machine learning_ , pp. 647–655. PMLR, 2014.
  * Ferradans et al. [2014] Sira Ferradans, Nicolas Papadakis, Gabriel Peyré, and Jean-François Aujol.  Regularized discrete optimal transport.  _SIAM Journal on Imaging Sciences_ , 7(3):1853–1882, 2014.
  * Flamary et al. [2016] R Flamary, N Courty, D Tuia, and A Rakotomamonjy.  Optimal transport for domain adaptation.  _IEEE Trans. Pattern Anal. Mach. Intell_ , 1, 2016.
  * Flamary et al. [2021] Rémi Flamary, Nicolas Courty, Alexandre Gramfort, Mokhtar Z. Alaya, Aurélie Boisbunon, Stanislas Chambon, Laetitia Chapel, Adrien Corenflos, Kilian Fatras, Nemo Fournier, Léo Gautheron, Nathalie T.H. Gayraud, Hicham Janati, Alain Rakotomamonjy, Ievgen Redko, Antoine Rolet, Antony Schutz, Vivien Seguy, Danica J. Sutherland, Romain Tavenard, Alexander Tong, and Titouan Vayer.  Pot: Python optimal transport.  _Journal of Machine Learning Research_ , 22(78):1–8, 2021.  URL http://jmlr.org/papers/v22/20-451.html.
  * Galanti et al. [2021] Tomer Galanti, Sagie Benaim, and Lior Wolf.  Risk bounds for unsupervised cross-domain mapping with ipms.  _J. Mach. Learn. Res._ , 22:90–1, 2021.
  * Gong et al. [2012] Boqing Gong, Yuan Shi, Fei Sha, and Kristen Grauman.  Geodesic flow kernel for unsupervised domain adaptation.  In _2012 IEEE conference on computer vision and pattern recognition_ , pp. 2066–2073. IEEE, 2012.
  * Hsu et al. [2015] Tzu Ming Harry Hsu, Wei Yu Chen, Cheng-An Hou, Yao-Hung Hubert Tsai, Yi-Ren Yeh, and Yu-Chiang Frank Wang.  Unsupervised domain adaptation with imbalanced cross-domain data.  In _Proceedings of the IEEE International Conference on Computer Vision_ , pp. 4121–4129, 2015.
  * Khan & Zhang [2022] Gabriel Khan and Jun Zhang.  When optimal transport meets information geometry.  _Information Geometry_ , pp. 1–32, 2022.
  * Knight [2008] Philip A Knight.  The sinkhorn–knopp algorithm: convergence and applications.  _SIAM Journal on Matrix Analysis and Applications_ , 30(1):261–275, 2008.
  * Li et al. [2020] Didong Li, Yulong Lu, Emmanuel Chevallier, and David B Dunson.  Density estimation and modeling on symmetric spaces.  _arXiv preprint arXiv:2009.01983_ , 2020.
  * Liu et al. [2019] Jie Liu, Yuanhao Huang, Ritambhara Singh, Jean-Philippe Vert, and William Stafford Noble.  Jointly embedding multiple single-cell omics measurements.  In _Algorithms in bioinformatics:… International Workshop, WABI…, proceedings. WABI (Workshop)_ , volume 143. NIH Public Access, 2019.
  * Liu et al. [2021] Yanbin Liu, Makoto Yamada, Yao-Hung Hubert Tsai, Tam Le, Ruslan Salakhutdinov, and Yi Yang.  Lsmi-sinkhorn: Semi-supervised mutual information estimation with optimal transport.  In _Joint European Conference on Machine Learning and Knowledge Discovery in Databases_ , pp. 655–670. Springer, 2021.
  * Mémoli [2011] Facundo Mémoli.  Gromov–wasserstein distances and the metric approach to object matching.  _Foundations of computational mathematics_ , 11(4):417–487, 2011.
  * Mokrov et al. [2021] Petr Mokrov, Alexander Korotin, Lingxiao Li, Aude Genevay, Justin M Solomon, and Evgeny Burnaev.  Large-scale wasserstein gradient flows.  _Advances in Neural Information Processing Systems_ , 34:15243–15256, 2021.
  * Moon et al. [1995] Young-Il Moon, Balaji Rajagopalan, and Upmanu Lall.  Estimation of mutual information using kernel density estimators.  _Physical Review E_ , 52(3):2318, 1995.
  * Mukherjee et al. [2021] Debarghya Mukherjee, Aritra Guha, Justin M Solomon, Yuekai Sun, and Mikhail Yurochkin.  Outlier-robust optimal transport.  In _International Conference on Machine Learning_ , pp. 7850–7860. PMLR, 2021.
  * Ozair et al. [2019] Sherjil Ozair, Corey Lynch, Yoshua Bengio, Aaron Van den Oord, Sergey Levine, and Pierre Sermanet.  Wasserstein dependency measure for representation learning.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Parzen [1962] Emanuel Parzen.  On estimation of a probability density function and mode.  _The annals of mathematical statistics_ , 33(3):1065–1076, 1962.
  * Pelletier [2005] Bruno Pelletier.  Kernel density estimation on riemannian manifolds.  _Statistics & probability letters_, 73(3):297–304, 2005.
  * Perrot et al. [2016] Michaël Perrot, Nicolas Courty, Rémi Flamary, and Amaury Habrard.  Mapping estimation for discrete optimal transport.  _Advances in Neural Information Processing Systems_ , 29, 2016.
  * Peyré et al. [2016] Gabriel Peyré, Marco Cuturi, and Justin Solomon.  Gromov-wasserstein averaging of kernel and distance matrices.  In _International Conference on Machine Learning_ , pp. 2664–2672. PMLR, 2016.
  * Poole et al. [2019] Ben Poole, Sherjil Ozair, Aaron Van Den Oord, Alex Alemi, and George Tucker.  On variational bounds of mutual information.  In _International Conference on Machine Learning_ , pp. 5171–5180. PMLR, 2019.
  * Redko et al. [2019] Ievgen Redko, Nicolas Courty, Rémi Flamary, and Devis Tuia.  Optimal transport for multi-source domain adaptation under target shift.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pp. 849–858. PMLR, 2019.
  * Rosenblatt [1956] Murray Rosenblatt.  Remarks on some nonparametric estimates of a density function.  _The annals of mathematical statistics_ , pp. 832–837, 1956.
  * Said et al. [2017] Salem Said, Lionel Bombrun, Yannick Berthoumieu, and Jonathan H Manton.  Riemannian gaussian distributions on the space of symmetric positive definite matrices.  _IEEE Transactions on Information Theory_ , 63(4):2153–2170, 2017.
  * Seguy et al. [2017] Vivien Seguy, Bharath Bhushan Damodaran, Rémi Flamary, Nicolas Courty, Antoine Rolet, and Mathieu Blondel.  Large-scale optimal transport and mapping estimation.  _arXiv preprint arXiv:1711.02283_ , 2017.
  * Solomon et al. [2016] Justin Solomon, Gabriel Peyré, Vladimir G Kim, and Suvrit Sra.  Entropic metric alignment for correspondence problems.  _ACM Transactions on Graphics (ToG)_ , 35(4):1–13, 2016.
  * Stoeckius et al. [2017] Marlon Stoeckius, Christoph Hafemeister, William Stephenson, Brian Houck-Loomis, Pratip K Chattopadhyay, Harold Swerdlow, Rahul Satija, and Peter Smibert.  Simultaneous epitope and transcriptome measurement in single cells.  _Nature methods_ , 14(9):865–868, 2017.
  * Tan et al. [2020] Shuhan Tan, Xingchao Peng, and Kate Saenko.  Class-imbalanced domain adaptation: An empirical odyssey.  In _European Conference on Computer Vision_ , pp. 585–602. Springer, 2020.
  * Van der Maaten & Hinton [2008] Laurens Van der Maaten and Geoffrey Hinton.  Visualizing data using t-sne.  _Journal of machine learning research_ , 9(11), 2008.
  * Vayer et al. [2018a] Titouan Vayer, Laetita Chapel, Rémi Flamary, Romain Tavenard, and Nicolas Courty.  Fused gromov-wasserstein distance for structured objects: theoretical foundations and mathematical properties.  _arXiv preprint arXiv:1811.02834_ , 2018a.
  * Vayer et al. [2018b] Titouan Vayer, Laetitia Chapel, Rémi Flamary, Romain Tavenard, and Nicolas Courty.  Optimal transport for structured data with application on graphs.  _arXiv preprint arXiv:1805.09114_ , 2018b.
  * Vayer et al. [2019] Titouan Vayer, Laetitia Chapel, Rémi Flamary, Romain Tavenard, and Nicolas Courty.  Optimal transport for structured data with application on graphs.  In _ICML 2019-36th International Conference on Machine Learning_ , pp. 1–16, 2019.
  * Villani [2009] Cédric Villani.  _Optimal transport: old and new_ , volume 338.  Springer, 2009.
  * Xu et al. [2020] Renjun Xu, Pelen Liu, Liyan Wang, Chao Chen, and Jindong Wang.  Reliable weighted optimal transport for unsupervised domain adaptation.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pp. 4394–4403, 2020.
  * Yurochkin et al. [2019] Mikhail Yurochkin, Amanda Bower, and Yuekai Sun.  Training individually fair ml models with sensitive subspace robustness.  _arXiv preprint arXiv:1907.00020_ , 2019.
  * Zhong et al. [2010] Erheng Zhong, Wei Fan, Qiang Yang, Olivier Verscheure, and Jiangtao Ren.  Cross validation framework to choose amongst models and datasets for transfer learning.  In _Joint European Conference on Machine Learning and Knowledge Discovery in Databases_ , pp. 547–562. Springer, 2010.

##  Appendix A Proofs

###  A.1 Proof of Lemma 2

  * Proof.

In the limit when $h\rightarrow 0$, the Gaussian kernel converges to

$\displaystyle K_{h}\left(t\right)$
$\displaystyle=\begin{cases}1/Z_{h}&\text{if $t=0$}\\\
0&\text{otherwise}.\end{cases}$

Therefore, the kernel $K_{h}(d_{\mathcal{X}}(x_{i},x_{k}))$ will only have
non-zero value when $x_{i}=x_{k}$, which implies that the kernelized mutual
information will converge as follows:

$\displaystyle\lim_{h\rightarrow 0}\hat{I}_{\Gamma}(X,Y)$
$\displaystyle=\lim_{h\rightarrow
0}\sum_{i,j}\Gamma_{ij}\log\frac{n^{2}\cdot\sum_{k,l}\Gamma_{kl}K_{h}\left(d_{\mathcal{X}}(x_{i},x_{k})\right)K_{h}\left(d_{\mathcal{Y}}(y_{j},y_{l})\right)}{\sum_{k}K_{h}\left(d_{\mathcal{X}}(x_{i},x_{k})\right)\cdot\sum_{l}K_{h}\left(d_{\mathcal{Y}}(y_{j},y_{l})\right)}.$
$\displaystyle=\sum_{i,j}\Gamma_{ij}\log\frac{n^{2}\cdot\Gamma_{ij}/Z_{h}^{2}}{1/Z_{h}^{2}}$
$\displaystyle=\sum_{i,j}\Gamma_{ij}\log\Gamma_{ij}+2\log(n)$
$\displaystyle=-H(\Gamma)+2\log(n).$

∎

###  A.2 Projected Gradient Descent

The classic mirror descent iteration is written as:

$\displaystyle
x_{t+1}\leftarrow\operatorname*{arg\,min}_{x}\left\\{\tau\langle\nabla
f(x_{t}),x\rangle+D(x\|x_{t})\right\\}.$

When $D(y||x)$ is the KL divergence:
$D_{\mathrm{KL}}(y||x)=\sum_{i}y_{i}\log\frac{y_{i}}{x_{i}}$, the update has
the following form:

$\displaystyle(x_{t+1})_{i}=e^{\log(x_{t})_{i}-\tau\nabla
f(x_{t})}=(x_{t})_{i}e^{-\tau\nabla f(x_{t})}.$

In our case, before the projection, the update reads

$\displaystyle\Gamma^{\prime}_{t+1}=\left(\Gamma_{t}\odot
e^{-\tau(C-\lambda\nabla_{\Gamma_{t}}\hat{I}_{\Gamma_{t}}(X,Y)-\epsilon\nabla
H(\Gamma_{t})})\right).$

Next, we solve the following projection w.r.t. KL metric:

$\displaystyle\Gamma_{t+1}=\operatorname*{arg\,min}_{\Gamma\in\Pi({\mathbf{p}},{\mathbf{q}})}D_{\mathrm{KL}}(\Gamma\|\Gamma^{\prime}_{t+1}).$

As [9] shows, the KL projection is equivalent to solving the entropic
regularized optimal transport problem, which is usually refer to the sinkhorn
distance [23]. Following [45], we set the stepsize $\tau=1/\epsilon$ to
simplify the iterations and reach the following update rule:

$\displaystyle\Gamma_{t+1}\leftarrow\operatorname*{arg\,min}_{\Gamma\in\Pi({\mathbf{p}},{\mathbf{q}})}\left\langle\Gamma,C-\lambda\nabla_{\Gamma}\hat{I}_{\Gamma_{t}}(X,Y)\right\rangle-\epsilon
H(\Gamma).$

##  Appendix B Additional Experiments

###  B.1 Color Transfer

Color transfer aims to transfer the colors of the target images into the
source image. Optimal transport achieves this by treating pixels as points in
the RGB space, and maps the source pixels to the target ones. Here, 500 pixels
are sampled from each image to compute the OT, then the barycentric projection
is applied to map all the source pixels to target. We compare fused InfoOT
with standard OT, Sinkhorn distance [23], and linear mapping estimation [44]
and show the results in Figure 7. We can see that InfoOT produces a sharper
results than the baselines while decently recovering the colors in the target
image.

Figure 7: Color Transfer via Optimal Transport. Fused InfoOT produces sharper
results while preserving the target color compared to the baselines.

###  B.2 Word Embedding Alignment

Here, we explore the possibility of applying InfoOT for unsupervised word
embedding alignment. We follow the setup in [2], where the goal is to recover
cross-lingual correspondences with word embedding in different languages. In
this case, the pairwise distance between domains might not be meaningful, as
the word embedding models are trained separately. Previous works suggest that
cross-lingual word vector spaces are approximately isometric, which makes
Gromov-Wasserstein an ideal choice due to its ability to align isometric
spaces. Here, we treat GW as the oracle, and show that InfoOT can perform
comparably to GW [2] and other baselines such as InvOT [4], Adv-NN [21], and
supervised PROCRUSTES. We report the results on the dataset of Conneau et al.
[21] in Table 4, where both GW and InfoOT are trained with 12000 words and
refined with Cross-Domain Similarity Scaling (CSLS) [21]. The entropy
regularizer is $0.0001$ and $0.02$ for GW and InfoOT, respectively. We can see
that InfoOT performs comparably with the baselines and GW, demonstrating its
applicability in recovering cross-lingual correspondence.

|  | EN-ES | EN-FR | EN-DE | EN-IT | EN-RU
---|---|---|---|---|---|---
| Supervision | $\rightarrow$ | $\leftarrow$ | $\rightarrow$ | $\leftarrow$ | $\rightarrow$ | $\leftarrow$ | $\rightarrow$ | $\leftarrow$ | $\rightarrow$ | $\leftarrow$
PROCRUSTES | 5K Words | 81.2 | 82.3 | 81.2 | 82.2 | 73.6 | 71.9 | 76.3 | 75.5 | 51.7 | 63.7
Adv-NN | None | 81.7 | 83.3 | 82.3 | 82.1 | 74.0 | 72.2 | 77.4 | 76.1 | 52.4 | 61.4
InvOT | None | 81.3 | 81.8 | 82.9 | 81.6 | 73.8 | 71.1 | 77.7 | 77.7 | 41.7 | 55.4
InfoOT (h=0.55) | None | 81.6 | 78.5 | 82.4 | 80.5 | 75.4 | 74.2 | 78.6 | 75.7 | 48.1 | 52.9
GW | None | 84.3 | 83.2 | 84.8 | 83.6 | 77.4 | 75.2 | 82.5 | 79.8 | 52.0 | 61.4
Table 4: Cross-lingual Word Alignment. The InfoOT achieves comparable
performance to GW, demonstrating its potential in recovering cross-lingual
correspondence.