Propensity Score Alignment of Unpaired Multimodal Data

Authors: johnny  xi
jana  osea
david  xu
jason  hartford
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 00:16:20

Abstract: Multimodal representation learning techniques typically require paired samples to learn shared representations, but collecting paired samples can be challenging in fields like biology, where measurement devices often destroy the samples. This paper presents an approach to address the challenge of aligning unpaired samples across disparate modalities in multimodal representation learning. We draw an analogy between potential outcomes in causal inference and potential views in multimodal observations, allowing us to leverage Rubin's framework to estimate a common space for matching samples. Our approach assumes experimentally perturbed samples by treatments, and uses this to estimate a propensity score from each modality. We show that the propensity score encapsulates all shared information between a latent state and treatment, and can be used to define a distance between samples. We experiment with two alignment techniques that leverage this distance---shared nearest neighbours (SNN) and optimal transport (OT) matching---and find that OT matching results in significant improvements over state-of-the-art alignment approaches in on synthetic multi-modal tasks, in real-world data from NeurIPS Multimodal Single-Cell Integration Challenge, and on a single cell microscopy to expression prediction task.

Main paper: 
 introduction: 
 Large-scale multimodal representation learning techniques such as CLIP [Radford et al., 2021] have lead to remarkable improvements in zero-shot classification performance and have enabled the recent success in conditional generative models. However, the effectiveness of multimodal methods hinges on the availability of paired samples-such as images and their associated captions-across data modalities. This reliance on paired samples is most obvious in the InfoNCE loss [Gutmann and Hyvärinen, 2010, van den Oord et al., 2018] used in CLIP [Radford et al., 2021] which explicitly learns representations to maximize the true matching between images and their captions.While paired image captioning data is abundant on the internet, paired multimodal data is often challenging to collect in scientific experiments. For instance, unpaired data are the norm in biology for technical reasons: RNA sequencing, protein expression assays, and the collection of microscopy Figure 1 : Visualization of propensity score matching for two modalities (e.g., Microscopy images and RNA expression data). We first train classifiers to estimate the propensity score for samples from each modalities; the propensity score reveals the shared information p(t|z i ), which allows us to re-pair the observed disconnected modalities. The matching procedure is then performed within each perturbation class based on the similarity bewteen the propensity scores.images for cell painting assays are all destructive processes. As such, we cannot collect multiple different measurements from the same cell, and can only explicitly group cells by their experimental condition. If we could accurately match unpaired samples across modalities, we could use the aligned samples as proxies for paired samples and apply existing multimodal learning techniques.In this paper, we formalize this setting by viewing each modality as a potential measurement, X (1) (Z) ∈ X (1) , X (2) (Z) ∈ X (2) , of the same underlying latent state Z ∈ Z, where we are only able to make a single measurement for each sample unit (e.g. an individual cell). The task is to reconcile (match) unpaired observations x (1) and x (2) with the same (or maximally similar) z. Estimating the latent, Z, is hopelessly underspecified without making unverifiable assumptions on the system, and furthermore, Z may still be sparse and high-dimensional, leading to inefficient matching. This motivates the need for approaches that use only the observable data.We identify two major challenges for this problem. First, measurements are often made in very different spaces X (1) and X (2) (e.g., pixel space and gene expression counts), which make defining a notion of similarity across modalities challenging. Second, the measurement process inevitably introduces modality-specific variation that can be impossible to disentangle from the relevant information (Z). For example in cell imaging, we would not want the matching to depend on irrelevant appearance features such as the orientation of the cell or the lighting of the plate.In this paper, we address these challenges by appealing to classical ideas from causal inference [Rubin, 1974] , in the case where we additionally observe some label t for each unit, e.g., indexing an experiment. By making the assumption that t perturbs the observations via their shared latent state, we identify an observable link between modalities with the same underlying z. Under conditions which we discuss in Section 2, the propensity score, defined as p(t|Z), is a transformation of the latent Z that satisfies three remarkable properties (Proposition 3.1): (1) it provides a common space for matching, (2) it is fully identifiable via classification on individual modalities, and (3) it maximally reduces the dimension of Z, retaining only the information revealed by the perturbations.The practical implementation of the methodology (as illustrated in Fig. 1 ) is then straightforward: we train two separate classifiers, one for each modality, to predict the treatment t applied to X (i) . We then match across modalities based on the similarity between the predicted probabilities (the propensity score) within each treatment group. This matching procedure is highly versatile and can be applied to match labeled observations between any modalities for which a classifier can be efficiently trained. However, since the same sample unit does not appear in both modalities, we cannot use naive bipartite matching. To address this, we use soft matching techniques to estimate the missing modality for each sample unit by allowing matching to multiple observations. We experiment with two recent matching approaches: shared nearest neighbours (SNN) matching [Lance et al., 2022, Cao and Gao, 2022] and optimal transport (OT) matching Villani [2009] .In our experiments, we find that OT matching with distances defined on the proposenity score leads to significant improvement on matching and a downstream cross-modality prediction task on both synthetic and real-world biological data. Notably, our prediction method, which leverages the soft matching to optimize an OT projected loss, outperforms supervised learning on the true pairs on CITE-seq data from the NeurIPS Multimodal Single-Cell Integration Challenge [Lance et al., 2022] . Finally, we applied our method to match single-cell expression data (from a PeturbSeq assay [Dixit et al., 2016] ) with single cell crops of image data [Fay et al., 2023] . We find improved generalization in predicting the distribution of gene expression from the cell imaging data in with unseen perturbations.Unpaired and Multimodal Data Learning from unpaired data has long been considered for image translation [Liu et al., 2017 , Zhu et al., 2017 , Almahairi et al., 2018] , and more recently for biological modality translation [Amodio and Krishnaswamy, 2018, Yang et al., 2021] . In particular, Yang et al. [2021] also takes the perspective of a shared latent variable for biological modalities. This setting has been studied more generally for multi-view representation learning [Gresele et al., 2020 , Sturma et al., 2023] for its identifiability benefits.Perturbations and Heterogeneity Many methods in biology treat observation-level heterogeneity as a nuisance dimension to globally integrate, even when cluster labels are observed [Butler et al., 2018 , Korsunsky et al., 2019 , Foster et al., 2022] . This is sensible when clusters correspond to noise rather than the signal of interest. However, it is well known in causal representation learning that heterogeneity-particularly heterogeneity arising from perturbations-has theoretical benefits in constraining the solution set [Khemakhem et al., 2020 , Squires et al., 2023 , Ahuja et al., 2023 , Buchholz et al., 2023 , von Kügelgen et al., 2023] . There, the benefits (weakly) increase with the number of perturbations, which is also true of our setting (Proposition 3.2). In the context of unpaired data, only Yang et al. [2021] explicitly leverage this heterogeneity in their method, while Ryu et al. [2024] treat it as a constraint in solving OT. Specifically, Yang et al. [2021] require their VAE representations to classify experimental labels in addition to reconstructing modalities, while our method is simpler, only requiring the classification objective. Notably, Yang et al. [2021] treat our objective as a regularizer, but our theory suggests that it is actually primarily responsible for the matching performance. Our experiment results coincide with the theoretical insights; requiring reconstruction, as in a VAE, led to worse matching performance with identical model architectures.Optimal Transport Matching OT is a common tool in single-cell biology. In cell trajectory inference, the unpaired samples are gene expression values measured at different time points in a shared (metric) space. OT matching minimizes this shared metric between time points [Schiebinger et al., 2019 , Tong et al., 2020] . Recent work [Demetci et al., 2022] extends this to our setting where each modality is observed in separate metric spaces by using the Gromov-Wasserstein distance, which computes the difference between the metric evaluated within pairs of points from each modality [Demetci et al., 2022] . In concurrent work, this approach was recently extended to ensure matching within experimental labels [Ryu et al., 2024] . In addition to these "pure" OT approaches, Gossi et al. [2023] use OT on contrastive learning representations, though this approach requires matched pairs for training, while Cao et al. [2022] use OT in the latent space of a multi-modal VAE.Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?Answer: [Yes] Justification: The abstract simply describes the methodology and our experimental results. 
setting: 
 We consider the setting where there exist two potential views, X (e) ∈ X (e) from two different modalities indexed by e ∈ {1, 2}, and experiment t that perturbs a shared latent state of these observations. This defines a jointly distributed random variable (X (1) , X (2) , e, t), from which we observe only a single modality, its index, and label, {x(ei) i , e i , t i } n i=1 . 2We aim to match or estimate the samples from the missing modality, which corresponds to the realization of the missing random variable. Since t is observed, in practice we match observations within the same label class t.Formally, we assume each modality is generated by a common latent random variable Z as follows:t ∼ P T , Z (t) | t ∼ P (t) Z , U (e) ∼ P (e)U , U (e) ⊥ ⊥Z, U (e) ⊥ ⊥U (e ′ ) , X (e) | t = f (e) (Z (t) , U (e) ), (1) where t indexes the experimental perturbations, and we take t = 0 to represent a base environment. U (e) represents the modality-specific measurement noise that is unperturbed by t, and also independent across samples. The structural equations f e are deterministic after accounting for the randomness in Z and U : it represents the measurement process that captures the latent state. For example, in a microscopy image, this would be the microscope and camera that maps a cell to pixels.Comparison to Multimodal Generative Models Our setting is technically that of a multimodal generative model with latent perturbations. However, by focusing on matching rather than generation, we are able to make significantly weaker and more meaningful assumptions while still ensuring the theoretical validity of our method. Without the effects of the perturbation, our Eq. ( 1) is essentially the same as [Yang et al., 2021, Equation 1] in an abstract sense. However, in order to fit the generative model, it is required to formulate explicit models over f (e) and P (t) Z , which requires specifying the function class (e.g., continuous) and the space of Z (e.g., R d ) as assumptions, even in universal approximation settings. In contrast, since we will not directly fit the model Eq. ( 1), we do not make any technical assumptions about the generative model. Instead, we will make the following assumptions on the underlying data generating process itself.Key Assumptions Our theory makes the following assumptions about the data generating process.(A1) t ̸ ⊥ ⊥Z, and t⊥ ⊥U (e) . In words, t has a non-trivial effect on Z, but does not affect U (e) , implying that interventions target the common underlying process without affecting modalityspecific properties. For example, an intervention that affects the underlying cell state, but not the measurement noise of the individual modalities.(A2) Injectivity of f (e) : f (e) (z, u) = f (e) (z ′ , u ′ ) =⇒ (z, u) = (z ′ , u ′ ).In words, each modality captures enough information to distinguish changes in the underlying state.3 (A1) ensures that the conditional distribution t | X (e) is identical for e = 1, 2. (A2) then ensures that t | X (1) d = t | X (2) d = t | Z,which allows us to estimate the conditional distribution t | Z with observed data alone. Though sharp assumptions are required for the theory, versions replaced with approximate distributional equalities intuitively also allow for effective matchings when combined with our soft matching procedures in practice. A particular relaxation of (A1) when combined with OT matching is described in Appendix A.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: We included a limitations section that explains the limitations of the method.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.• If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
multimodal propensity scores: 
 Under (1), if Z were observable, an optimal matching can be constructed by simply matching the samples with the most similar z i . However, the prerequisite of inverting the model and disentangling Z is arguably more difficult than the matching problem itself. In particular, Z is unidentifiable without strong assumptions on Eq. ( 1) [Xi and Bloem-Reddy, 2023] , and even formulating the identifiability problem requires well-specification of the model as a prerequisite. We take an alternative approach that is robust to these problems, by using the perturbations t as an observable link to reveal information about Z. Specifically, we show that the propensity scoreEQUATIONis identifiable as a proxy for the latent Z under our assumptions of the data generating process. This is a consequence of the injectivity of f (e) , since it will be that π(Z) = π(X (e) ), e = 1, 2, indicating that we can compute it from either modality. Not only does the propensity score reveal shared information, classical causal inference theory [Rubin, 1974] states that it captures all information about Z that is contained in t, and does so minimally, in terms of having minimum dimension and entropy. Since t contains the only observable information that is useful for matching, the propensity score is hence an optimal compression of the observed information. We collect these observations into the following proposition.Proposition 3.1. In the model described by Eq. (1), further assume that f (e) are injective for e = 1, 2. Then, the propensity scores in either modality is equal to the propensity score given by Z,i.e., π(X (1) ) = π(X (2) ) = π(Z) as random variables. This impliesEQUATIONfor each e = 1, 2, where I is the mutual information. Furthermore, any other function b(Z) satisfyingI(t, Z | b(Z)) = 0 is such that π(Z) = f (b(Z)).The proof can be found in Appendix C. Practically, Proposition 3.1 shows that computing the propensity score on either modality is equivalent to computing it on the unobserved shared latent, which means that it is identifiable, and thus estimable, from the observations alone. Furthermore, the estimation does not require modified objectives or architectures for joint multimodal processing, instead they are simple and separate classification problems for each modality. Finally, t does not affect U (e) by assumption, and thus the propensity score, being a representation of the information in t, discards the modality-specific information that may be counterproductive to matching. Therefore, even if Z were observed, it may be sensible to match on its propensity score instead.Number of Perturbations Note that point-wise equality of the propensity score π(z 1 ) = π(z 2 ) does not necessarily imply equality of the latents z 1 = z 2 , due to potential non-injectivity of π. For example, consider t ∈ {0, 1}, then π(z) is a compression to a single dimension z → p(t = 1 | z).Intuitively, collecting data from more perturbations improves the amount of information contained in the label t. If the latent space is R d , the propensity score necessarily compresses information about Z (t) if the latent dimension exceeds the number of perturbations, echoing impossibility results from the causal representation learning literature [Squires et al., 2023] .Proposition 3.2. Let Z (t) ∈ R d . Suppose that P (t)Z has a smooth density p(z|t) for each t = 0, . . . , T . Then, if T < d, the propensity score π, restricted to its strictly positive part, is non-injective.The proof can be found in Appendix C. Note the above only states an impossibility result when T < d. More generally, it can be seen from the proof of Proposition 3.2 that the injectivity of the propensity score depends on the injectivity of the following expression in z:EQUATIONwhich then depends on the latent process itself. If the above mapping is non-injective, this represents a fundamental indeterminacy that cannot be resolved without making strong assumptions on point-wise latent variable recovery. As we have already established in Proposition 3.1, the propensity score contains the maximal shared information across modalities. Nonetheless, collecting data form a larger number of perturbations is clearly beneficial for matching, since g in Eq. ( 4) is injective if any of the subset of its entries are.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer: [Yes] Justification: Our assumptions are clearly stated, and the reasoning why they are required clearly justified, in (A1) and (A2).Guidelines:• The answer NA means that the paper does not include theoretical results.• All the theorems, formulas, and proofs in the paper should be numbered and crossreferenced. • All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 
estimation and matching: 
 For the remainder of the paper, we drop the notation e and use (x i , t i ) to denote observations from modality 1, and (x j , t j ) to denote observations from modality 2. Given a multimodal dataset with observations {(x i , t i )} n1 i=1 and {(x j , t j )} n2 j=1 , we wish to compute a matching matrix (or coupling) between the two modalities. We define a n 1 × n 2 matching matrix M where M ij represents the likelihood of x i being matched to x j . Since t is observed, we always perform matching only within observations with the same value of t, so that in practice we obtain a matrix M t for each t.Our method approximates the propensity scores by training separate classifiers that predicts t given x for each modality. We denote the estimated propensity score by π i and π j respectively, whereπ i ≈ π(x i ) = P (T = t | X (e) i = x i ).(5)This yields the transformed datasets {π i } n1 i=1 and {π j } n2 j=1 , where π i , π j are in the T dimensional simplex. We use this correspondence to compute a cross-modality distance function:EQUATIONIn practice, we typically compute the Euclidean distance in R T of the logit-transformed classification scores, but any metric over a bijective transformation of the propensity scores are also theoretically valid. Given this distance function, we use existing matching techniques to constructing a matching matrix. In our experiments, we found that OT matching gave the best performance, but we also evaluated Shared Nearest Neighbour matching; details of the latter can be found in Appendix B.The propensity score distance allows us to easily compute a cost function associated with transporting mass between modalities, c(x i , x j ) = d ′ (π i , π j ). Let p 1 , p 2 denote the uniform distribution over {π i } n1 i=1 and {π j } n2 j=1 respectively. Discrete OT aims to solve the problem of optimally redistributing mass from p 1 to p 2 in terms of incurring the lowest cost. Let C ij = c(x i , x j ) denote the n 1 × n 2 cost matrix. The Kantorovich formulation of optimal transport aims to solve the following constrained optimization problem:EQUATION)This is a linear program, and for n 1 = n 2 , it can be shown that the optimal solution is a bipartite matching between {π i } n1 i=1 and {π j } n2 j=1 . We refer to this as exact OT; in practice we add an entropic regularization term, resulting in a soft matching, that ensures smoothness and uniqueness, and can be solved efficiently using Sinkhorn's algorithm. Entropic OT takes the following form:EQUATION)whereH(M ) = -i,j M ij log(M ij ), the entropy of the joint distribution implied by M . This approach regularizes towards a higher entropy solution, which has been shown to have statistical benefits [Genevay et al., 2018] , but nonetheless for small enough λ serves as a computationally appealing approximation to exact OT.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?Answer: [Yes] Justification: Full code, complete with scripts to obtain synthetic or publicly available data for experimental datasets (1) and (2), is provided as supplementary material. Dataset (3) is not released publicly, but enough details are given in the paper, and in descriptions of the method, to reproduce on similar data. 
downstream tasks: 
 The matching matrix M can be seen as defining an empirical joint distribution over the samples in each modality. The OT approach in particular makes this explicit. Each row is proportional to the probability that each sample i from modality (1) is matched to sample j in modality (2), i.e., M i,j = P (x j |x i ). We can thus use M to obtain pseudosamples for any learning task that uses paired samples by (x i , xj ), where xj is obtained by sampling from the conditional distribution defined by M , or by a suitable conditional expectation, e.g., the barycentric projection (conditional mean) asE M [X j | X i = x i ] = j M i,j x j .In what follows, we describe a cross-modality prediction method based on both barycentric projection and stochastic gradients according to M i,j .We can use the matching matrix to design a method for cross-modality prediction/translation. The following MSE loss corresponds to constructing a prediction function f θ such that the barycentric projectionE M [f θ (X j ) | X i = x i ],under M minimizes the squared error for predicting x i :EQUATION)However, this requires evaluating f θ for all n 2 examples from modality (2) for each of the n 1 examples in modality (1). In practice, we can avoid this cost with stochastic gradient descent by sampling from modality (2) via M i• for each training example (1). To obtain an unbiased estimate of ∇ θ L, we need two independent samples from modality (2) for each sample from modality (1),EQUATIONBy taking two samples as in Eq. ( 10), we get an unbiased estimator of ∇L(θ), whereas a single sample would have resulted in optimizing an upper-bound on equation ( 9); for details, see Hartford et al. [2017] where a similar issue arises in the gradient of their causal effect estimator. We thus refer to prediction models trained via Eq. ( 10) as unbiased. 
experiments: 
 We present a comprehensive evaluation of our proposed methodology on three distinct datasets: (1) synthetic paired images, (2) single-cell CITE-seq dataset (simultaneous measurement of single-cell RNA-seq and surface protein measurements) [Stoeckius et al., 2017] , and (3) Perturb-seq and singlecell image data. In the first two cases, there is a ground-truth matching that we use for evaluation, but samples are randomly permuted during training. This allows us to exactly compute the quality of the matching in comparison to the ground truth. The final dataset is a more realistic setting where ground truth paired samples do not exist, and matching becomes necessary in practice. In this case, we compute distributional metrics to compare our proposed methodology against other baselines.Experimental Details All models for the experiments are implemented using Torch v2.2.2 [Paszke et al., 2017] and Pytorch Lightning v2.2.4 [Falcon and PyTorch Lightning Team, 2023]. The classifier used to estimate the propensity score is always a linear head on top of an encoder E i , which is specific to each modality and dataset. All models are saved at the optimal validation loss to perform subsequent matching. Shared nearest neighbours (SNN) is implemented using scikit-learn v1.4.0 [Pedregosa et al., 2011 ] using a single neighbour, and OT is implemented using the Sinkhorn algorithm as implemented in the pot v0.9.3 package [Flamary et al., 2021] .Both SNN and OT use the Euclidean distance as the metric. Whenever random variation can affect the results of the experiments, we report quantiles corresponding to variation from different random seeds. Additional experimental details are provided in Appendix D.Description of Baselines Our main baseline, which we evaluate against on all three datasets, is matching using representations learned by the multimodal VAE of Yang et al. [2021] , which is the only published method that is able to leverage perturbation labels for unpaired multimodal data (they refer to the labels as "prior information"). The standard multimodal VAE loss is a reconstruction loss based on encoder and decoders E i , D i for each modality, plus a latent invariance loss that aims to align the modalities in the latent space. In our setting, the multimodal VAE loss further includes an additional label classification loss from the latent space of each modality, i.e., encouraging the encoder to simultaneously learn P (t | E i (x i )). This additional objective, which acts as a regularizer for the multimodal VAE, is exactly the loss for our proposed method. To ensure a fair comparison, we always use the same architecture in the encoders E i of multimodal VAE and in our propensity score classifier. The performance differences between propensity score matching and multimodal VAE then represent the effects of the VAE reconstruction objective and latent invariance objectives. For additional baselines, we also compare against a random matching, where the samples are matched with equal weight within each perturbation as a sanity check. For datasets (1) and (2), we also compare against Gromov-Wasserstein OT (SCOT) [Demetci et al., 2022] computed separately within each perturbation. SCOT uses OT directly by computing a cost function derived based on pairwise distances within each modality, thus learning a local description of the geometry which can be compared between modalities. For the CITE-seq dataset, we also compare against matching using a graph-linked VAE, scGLUE [Cao and Gao, 2022] , where the graph is constructed from linking genes with the associated proteins.Evaluation Metrics We use the known ground truth matching to compute performance metrics on datasets (1) and ( 2). The trace and FOSCTTM [Liu et al., 2019] measure how much weight M places on the true pairing. However, this is not necessarily indicative of downstream performance as similar, but not exact matches are penalized equally to wildly incorrect matches. For this reason, we also measure the latent MSE for dataset (1) and the performance of a CITE-seq gene-to-protein predictive model based on the learned matching for dataset (2). For more details, see Appendix D.1.Data We followed the data generating process Eq. ( 1) with a latent variable Z encoding the coordinates of two objects. Perturbations represent different do-interventions on the different dimensons of Z. The difference between modalities corresponds to whether the objects are circular or square, and a fixed transformation of Z, while the modality-specific noise U controls background distortions.We used a convolutional neural network adapted from Yang et al. [2021] as the encoder. We report two evaluation metrics: (1) the trace metric, and (2) the MSE between the matched and the true latents. The latent MSE metric does not penalize close neighbours of the true match (i.e. examples for which ∥z i -z * i ∥ is small) as heavily as the trace metric. These "near matches" will typically still be useful on downstream multimodal tasks.In Table 1 , metrics are computed on a held out test set over 12 groups corresponding to interventions on the latent position, with approximately 1700 observations per group. A random matching, with weight 1/n, will hence have a trace metric of of 1/1700 ≈ 0.588 × 10 -3 . This implies, for example, that the median performance of PS+OT is approximately 31 times that of random matching. On both metrics, we found that propensity scores matched with OT (PS + OT) consistently outperformed other matching methods on both metrics.Data We used the CITE-seq dataset from the NeurIPS 2021 Multimodal single-cell data integration competition [Lance et al., 2022] , consisting of paired RNA-seq and surface level protein measurements over 45 cell types. In the absence of perturbations, we used the cell type as the observed label to classify and match within. Note the cell types are determined by consensus by pooling annotations from marker genes/proteins. In most cells, the annotations from each modality agreed, suggesting that the label is independent from the modality-specific noise. We used the first 200 principal components as the gene expression modality, and normalized (but otherwise raw) protein measurements as input.We used fully-connected MLPs as encoders. To assess matching, we report (1) the trace, and (2) the Fraction Of Samples Closer Than the True Match (FOSCTTM) ( [Demetci et al., 2022] , [Liu et al., 2019] ) (lower is better, 0.5 corresponds to random guessing). To evaluate against a downstream task, we also compared the performance of random and VAE matching procedures, as well as directly using the ground truth (M ii = 1), on predicting protein levels from gene expression. We trained a 2-layer MLP (the same architecture for all matchings) with both MSE loss and the unbiased procedure as described in Section 5 using pseudosamples sampled according to the matching matrix. We evaluated the predictive models against ground truth pairs by computing the prediction R 2 (higher is better) on a held-out, unpermuted, test set.In Table 1 , metrics are computed on a held-out test set averaged over 45 cell types with varying observation counts per group. While interpreting the average trace can be challenging due to group size variations, OT matching on PS consistently outperformed other methods both within and across groups. In these experiments, OT matching on PS was consistently the top performer, often followed by SNN matching on PS or OT matching on VAE embeddings.We present downstream task performance in Table 2 . Note that R 2 is computed using the sample average across possibly multiple cell types, which explains why random matching within each cell type results in non-zero R 2 (see Appendix D.1). We found that PS + OT matching outperforms other methods on this task. Surprisingly, the PS + OT prediction model performed even better on average than training with the standard MSE loss on ground truth pairings (though confidence intervals overlap). This highlights the potential benefit of soft (OT) matching as a regularizer, beyond that of simply reconciling most likely pairs: the soft matching effectively averages over modality specific variation from samples with similar latent states in a manner analogous to data augmentation (with an unknown group action).Data We collected PerturbSeq data (200 genes) and single-cell images of HUVEC cells with 24 gene perturbations and a control perturbation, resulting in 25 total labels across both modalities. As preprocessing, we embed the raw PerturbSeq counts into a 128-dimensional space using scVI [Lopez et al., 2018] and the cell images into a 1024-dimensional space using a pre-trained Masked Autoencoder [He et al., 2022 , Kraus et al., 2023] to train our gene expression and image classifiers.We used a fully connected 2-layer MLP as the encoder for both PerturbSeq and cell image classifiers. Similarly to the CITE-seq dataset, we evaluated the matchings based on downstream prediction of gene expression from (embeddings of) images. We used the unbiased procedure to minimize the projected loss Eq. ( 9) and evaluated on two held-out sets, one consisting of in-distribution samples from the 25 perturbations the classifier was trained on, and an out-ofdistribution set consisting of an extra perturbation not seen in training. In the absence of ground truth matching, we assessed three distributional metrics between the actual and predicted gene expression values within each perturbation: the L2 norm of the difference in means, the Kullback-Leibler (KL) divergence, and 1-Wasserstein distance (lower indicates better alignment). We report inverse cell-count weighted averages over each perturbation group. Each metric measures a slightly different aspect of fit-the L2 norm reports a first-order deviation, while the KL divergence is an empirical estimate of the deviation of the underlying predicted distribution, while the 1-Wasserstein distance measures deviations in terms of the empirical samples themselves.Note that matching is performed using classifiers trained on scVI embeddings, but the cross-modal predictions are generated in the original log transformed gene expression space (i.e. we predicted actual observations, not embeddings). We also evaluated distance measures on an out-of-distribution gene perturbation that was not used in either the matching or training of the translation model.We present KL divergence values for in-distribution and out-of-distribution in Table 2 . 4Additional metrics show similar patterns and can be found in Appendix D.4. OT + PS matching consistently outperforms its VAE counterpart both on in-distribution and out-of-distribution metrics, supporting our findings on the CITE-seq data to the case where ground truth pairs are not available.As in our Perturb-seq and cell imaging example, the ground truth matching is typically unknown in real problems. It is hence desirable to have an observable proxy of the matching performance Notice that validation cross-entropy inversely tracks the ground truth matching metrics, and thus can be used as a proxy in practical settings where the ground truth is unknown. The same pattern does not hold for the VAE [Yang et al., 2021] , which we suspect is because reconstruction is largely irrelevant for matching. as a validation during hyperparameter tuning. Figure 2 demonstrates that the propensity score validation loss (cross-entropy) empirically satisfies this role in our CITE-seq experiments, where lower validation loss corresponds to better matching performance, as if it were computed with the ground truth. By contrast, we found that the optimal VAE, in terms of matching, had higher validation loss. This empirically supports our intuition that the reconstruction loss minimization requires the VAE to capture modality specific information, i.e., the U (e) variables, which hinders its matching performance. 
limitations: 
 Our methods are limited to settings where we have some signal to play the role of an experiment label, but we believe this is where these methods are most needed. Matching is impossible in general-e.g., if you tried to match modalities that have no shared information, it would clearly fail-but our theory formally articulates both where we expect this method to succeed and its limitations. Both (A1) and (A2) are strong assumptions, but the empirical results suggest the method is fairly robust to failures. 
conclusion: 
 This work presents a simple algorithm for aligning unpaired data from different modalities using propensity scores. The method is very general, requiring only a classifier to be trained on each modality, and demonstrates excellent matching performance, which we validate both theoretically and empirically. We also showcase the effectiveness of the matching algorithm in a downstream cross-modality prediction task, achieving better generalization compared to random matching, VAEbased matching, and even the ground truth matching on the evaluated dataset. This improved generalization over the ground truth may be attributed to implicitly enforcing invariance to modalityspecific information; a rigorous investigation of this phenomenon would be interesting for further investigation. 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained? Answer: [NA] Justification: Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: Relaxing Assumption 1 Consider the propensity score π(x (e,t) ) = P (t|X (e,t) = x (e,t) )(11)where we do not necessarily require U (e) ⊥ ⊥t | Z (t) , and thus we obtainEQUATIONsee the proof of Proposition 3.1 for details.Suppose that the two observed modalities are indeed generated by a shared {z i } n i=1 , but where the indices of modality 2 are potentially permuted, and with values differing by modality specific information:EQUATIONEQUATIONwhere j = π(i) denotes a permutation of the sample indices. Under (A1), we would be able to find some j such that π(x(1,t) i ) = π(x (2,t) j) for each i.Matching via OT can allow us to relax (A1) in a very particular way. Consider the simple case where t ∈ {0, 1}, so that π can be written in a single dimension, e.g., P (t = 1|X (e,t) = x (e,t) ) ∈ [0, 1]. In this case, exact OT is equivalent to sorting π(x(1,t) i) and π(x(2,t) j), and matching the sorted versions 1-to-1. Under (A1), the sorted versions will be exactly equal. A relaxed version of (A1) that would still result in the correct ground truth matching is to assume that t affects U (1) and U (2) differently, but that the difference is order preserving, or monotone. Denote (π(xEQUATIONi ) as the true pairing, noting that we use the same index i. We require the following:EQUATIONThis says that, even if π(1) i̸ = π (2)i , that their relative orderings will still coincide. Then, exact OT will still recover the ground truth matching. See Fig. 3 for a visual example of this type of monotonicity. For example, suppose that t is a chemical perturbation of a cell, and thus π(1) i , π(2) j can be seen as a measure of biological response to the perturbation, e.g., in a treated population, π i1 > π i2 indicates samples i 1 had a stronger response than sample i 2 , as perceived by the first modality indexed by i. Then, this monotonocity states that we should see the same π j1 > π j2 in the other modality as well, if the samples i 1 and i 2 truly corresponded to j 1 and j 2 .We can see the monotonicity requirement Eq. ( 14) as the monotonicity of the function with graph (πEQUATIONi ) ∈ [0, 1] 2 . In higher dimensions, we require that the "graph" satisfies the following cyclic monotonicity property [Villani, 2009] :Definition A.1. The collection {(π (1) i , π (2) i )} ni=1 is said to be c-cyclically monotone for some cost function c, if for any n = 1, . . . , N , and any subset of pairs (πEQUATIONImportantly, we define π n+1 = π 1 , so that the sequence represents a cycle.Note in our setting, the OT cost function is the Euclidean distance, c(x, y) = ∥x -y∥ 2 . It is known that the OT solution must satisfy cyclic monotonicity. Thus, if the true pairing is uniquely cyclically monotone, we can recover it with OT. However, we are unaware of common violations of (A1) that would satisfy cyclic monotonicity. Modality 2 Ground truth matchingP(t | z, u (2) ) P(t | z, u (1) ) P(t | z, u (2) )Figure 3 : OT matching allows for t to have different effects on the modality specific information, here u(1)i and u(2)i , as long as they can be written as transformations that preserve the relative order within modalities. Exact OT in 1-d always matches according to the relative ordering, and thus exhibits this type of "no crossing" behaviour shown in the figure on the left. The figure on the right shows a case where we would fail to correctly match across modalities because of the crossing shown in orange.Using the propensity score distance, we can compute nearest neighbours both within and between the two modalities. We follow Cao and Gao [2022] and compute the normalized shared nearest neighbours (SNN) between each pair of observations as the entry of the matching matrix. For each pair of observations (π(1) i , π(2) j ), we define four sets:• 11 ij : the k nearest neighbours of π (1) i amongst {π (1) i } n1 i=1 . π (1) i is considered a neighbour of itself. • 12 ij : the k nearest neighbours of π (2) j amongst {π (1) i } n1 i=1 . • 21 ij : the k nearest neighbours of π (1) i amongst {π (2) j } n2 j=1 . • 22 ij : the k nearest neighbours of π (2) j amongst {π (2) j } n2 j=1 . π (2)j is considered a neighbour of itself.(1) i and π(2) j correspond to the same underlying propensity score, their nearest neighbours amongst observations from each modality should be the same. This is measured as a set difference between 11 ij and 12 ij , and likewise for 21 ij and 22 ij . Then, a modified Jaccard index is computed as follows. DefineEQUATION)the sum of the number of shared neighbours measured in each modality. Then, we compute the following Jaccard distance to populate the unnormalized matching matrix:EQUATION)where notice that 4k =|11 ij | + |12 ij | + |21 ij | + |22 ij |, since each set contains k distinct neighbours,and thus 0 ≤ Mij ≤ 1, as with the standard Jaccard index. Then, we normalize each row to produce the final matching matrix:EQUATION)Note M ij is always well defined because π(1) i and π(2) j are always considered neighbours of themselves. Lemma B.1. Mij has at least one non-zero entry in each of its rows and columns for any number of neighbours k ≥ 1.Proof. We prove that J ij > 0 for at least one j in each i, which is equivalent to Mij > 0. Fix an arbitrary i. 21 ij by definition is the same set for every j. By the assumption of k ≥ 1 it is non-empty, so there exists π(2)j * ∈ 21 ij . Since π (2) j * is a neighbour of itself, we have π (2) j * ∈ 22 ij * , showing that J ij * > 0.The same reasoning applied to 11 and 12 also shows that J ij for at least one i in each j.C.1 Proof of Proposition 3.1Proof. Let x (e) denote the observed modality and z, u (e) be the unique corresponding latent values. By injectivity,EQUATIONfor e = 1, 2, since we assumed U (e) ⊥ ⊥t | Z. Since this holds pointwise, it shows that π(X (1 ) = π(X (2) ) = π(Z) as random variables. Now, a classical result of Rubin [1974] gives that Z⊥ ⊥t | π(Z), and that for any other function b (a balancing score) such that Z⊥ ⊥t | b(Z), we have π(Z) = g(b(Z)).The first property written in information theoretic terms yields,EQUATIONsince π(X (e) ) = π(Z (t) ) as random variables, as required.C.2 Proof of Proposition 3.2In what follows, we write π to be the restriction to its domain where it is strictly positive. The i-th dimension of the propensity score can be written asEQUATIONwhich, when restricted to be strictly positive, maps to the relative interior of the T -dimensional probability simplex. Consider the following transformation:EQUATIONwhere C = log(p(t = i)) -log(p(t = 0)) is constant in z, and that h(π(z)) 0 ≡ 0. Ignoring the constant first dimension, we can view h as an invertible map to R T . Under this convention, the map h • π : R d → R T is smooth (log is smooth, and the densities are smooth by assumption). Since it is smooth, it cannot be injective if T < d [Spivak, 2018] . Finally, since h is bijective, this implies that π cannot be injective.In the synthetic image and CITE-seq datasets, a ground truth matching is known, and we can evaluate the quality of the synthetic matching directly against the truth. In these cases, the dataset sizes are necessarily balanced, so that n = n 1 = n 2 . In each case, we evaluate the quality of our n × n matching matrix M , which we compute within samples with the same t. Our reported results are then averaged over each cluster. Note we randomize the order of the datasets before performing the matching to avoid pathologies.Trace Metric Assuming the sample indices correspond to the true matching, we can compute the average weight on correct matches, which is the normalized trace of M :EQUATION)As a baseline, notice that a uniformly random matching that assigns M ij = 1/n for each cell yields Tr(M ) = 1 and hence will obtain a metric of 1/n. This metric however does not capture potential failure modes of matching. For example, exactly matching one sample, while adversarially matching dissimiliar samples for the remainder also yields a trace of 1/n, which is equal to that of a random matching.Latent MSE On the image dataset, we have access to the ground truth latent values that generated the images, z = {z i } n i=1 . We compute matched latents as M z, the barycentric projection according to the matching matrix. Then, to evaluate the quality of the matching in terms of finding similar latents, we compute the MSE:EQUATION)FOSCTTM We do not have access to ground truth latents in the CITE-seq dataset, so use the Fraction Of Samples Closer Than the True Match (FOSCTTM) [Demetci et al., 2022 , Liu et al., 2019] as an alternative matching metric. First, we use M to project 2) . Then, we can compute a cross-modality distance as follows. For each point in x(1) , we compute the Euclidean distance to each point in x (1) , and compute the fraction of samples in x (1) that are closer than the true match. We also repeat this for each point in x (1) , computing the fraction of samples in x(1) in this case. That is, assuming again that the given indices correspond to the true matching, we compute:x (2) = {x j } n j=1 to x (1) = {x i } n i=1 as x(1) = M x (EQUATIONEQUATIONwhere notice that this evaluates M through the computation x(1) = M x (2) . As a baseline, we should expect a random matching, when distances between points are randomly distributed, to have an FOSCTTM of 0.5.We also trained a cross-modality prediction (translation) model f θ,M to predict CITE-seq protein levels from gene expression based on matched pseudosamples. LetEQUATION= {x j } denote protein and gene expression, respectively. We trained a simple 2-layer MLP minimizing either the standard MSE, using pairs (x i , xj ), xj ∼ M i• , or following the projected loss with unbiased estimates in Section 5. Each batch in general consists of samples from all t, but the xj sampling step occurs within the perturbation. Let x(1) test = {f θ,M (x j )}. We report the R 2 on a randomly held-out test set of ground truth pairs (again, consisting of samples from all t), which is defined as the following:EQUATIONwherex(1) test is the naive mean (over all perturbations) estimator which acts as a baseline.We train a cross-modality prediction model to predict gene expression from cell images based on matched pseudosamples in the same way as in CITE-seq, but only using the projected loss with unbiased estimates. Denote this model for a matching matrix M by f θ,M .Because we do not have access to ground truth pairs within each perturbation, we resort to distributional metrics. Let x (1) t = {x i,t } nt1 i=1 , x (2) t = {x j,t } nt2 j=1 denote gene expression and cell images in a held out test set respectively in perturbation t. Letx(1) t = {f θ,M (x j,t )} nt2 j=1 . We compute empirical versions of statistical divergencesEQUATIONwhere D is either the L2 norm of the difference in empirical mean, empirical Kullback-Leibler divergence or 1-Wasserstein distance. We report these weighted averages of D t over the perturbations t according to the number of samples in the modality of prediction interest.In this section we describe experimental details pertaining to the propensity score and VAE [Yang et al., 2021] . SCOT [Demetci et al., 2022] and scGLUE [Cao and Gao, 2022] are used according to tutorials and recommended default settings by the authors.The propensity score approach minimizes the standard cross-entropy loss for both modalities. The VAE includes, in addition to the standard ELBO loss (with parameter λ on the KL term), two cross-entropy losses based on classifiers from the latent space: one, weighted by a parameter α to classify t as in the propensity score, and another, weighted by a parameter β, that classifies which modality the latent point belongs to.We use the Adam optimizer with learning rate 0.0001 and one cycle learning rate scheduler. We follow Yang et al. [2021] and set α = 1, β = 0.1, but found that λ = 10 -9 (compared to λ = 10 -7 in Yang et al. [2021] ) resulted in better performance. We used batch size 256 in both instances and trained for either 100 epochs (image) or 250 epochs (CITE-seq).For the VAE and classifiers of experiment 3, we use an Adam optimizer with learning 0.001 and weight decay 0.001 and max epoch of 100 (PerturbSeq) and 250 (single cell images) using batch sizes of 256 and 2048 correspondingly. We follow similar settings as Yang et al. [2021] and implement α = 1 with λ = 10 -9 , and since we do not have matched data, β = 0. For the cross-modal prediction models in experiment 3, we use Stochastic Gradient Descent optimizer with learning rate 0.001 and weight decay 0.001 with max epochs 250 and batch size 256. We implement early stopping with delay of 50 epochs which we then checkpoint the last model to use for downstream tasksArchitecture For the synthetic image dataset, we use an 5-layer convolutional network (channels = 32, 54, 128, 256, 512) with batch normalization and leaky ReLU activations, with linear heads for classification (propensity score and VAE) and posterior mean and variance estimation (VAE). For the VAE, the decoder consists of convolutional transpose layers that reverse those of the encoder.For the CITE-seq dataset, we use a 5-layer MLP with constant hidden dimension 1024, with batch normalization and ReLU activations (adapted from the fully connected VAE in Yang et al. [2021] ) as both the encoder and VAE decoder. We use the same architecture for both modalities, RNA-seq (as we process the top 200 PCs) and protein.For the PerturbSeq classifier encoder, we use a 2-layer MLP architecture. Each layer consists of a linear layer with an output feature dimension of 64, followed by Rectified Linear Unit (ReLU) activation, Batch Normalization, and dropout (p=0.1). A final layer with Leaky ReLU activation that brings dimensionality to 128 before feeding into a linear classification head with an output feature dimension of 25.For the single-cell image encoder classifier, we use a proprietary Masked Autoencoder [Kraus et al., 2023] to generate 1024-dimensional embeddings. Subsequently, a 2-layer MLP is trained on these embeddings. Each MLP layer has a linear layer, Batch Normalization, and Leaky ReLU activation. The output feature dimensions of the linear layers are 512 and 256, respectively, and the latent dimension remains at 1024 before entering a linear classification head with an output feature dimension of 25.Optimal Transport We used POT [Flamary et al., 2021] to solve the entropic OT problem, using the log-sinkhorn solver, with regularization strength γ = 0.05. We also use the CITE-seq dataset from Lance et al. [2022] as a real-world benchmark (obtained from GEO accession GSE194122). These consist of paired RNA-seq and surface level protein measurements, and their cell type annotations over 45 different cell types. We used scanpy, a standard bioinformatics package, to perform PCA dimension reduction on RNA-seq by taking the first 200 principal components. The protein measurements (134-dimensional) was processed in raw form. For more details, see Lance et al. [2022] .PerturbSeq and Single Cell Image Data We collect single-cell PerturbSeq data (200 genes) and single-cell painting images in HUVEC cells with 24 gene perturbations and a control perturbation, resulting in 25 labels for matching across both modalities. The target gene perturbations are selected based on the 24 genes with the highest number of cells affected by the CRISPR guide RNAs targeting those genes. The PerturbSeq data is filtered to include the top 200 genes with the highest mean count, then normalized and log-transformed. The single-cell painting images are derived from multi-cell images, with each single-cell nucleus centered within a 32x32 pixel box. We use scVI Lopez et al. [2018] to embed the raw PerturbSeq counts into a 128-dimensional space before training the gene expression classifier. Similarly, we train our image classifier using 1024-dimensional embeddings obtained from a pre-trained Masked Autoencoder Kraus et al. [2023] , He et al. [2022] . Following matching, we perform cross-modality translation from the single-cell embeddings to the transformed gene expression counts.In • The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [NA] Justification: The paper describes methodology that may advance the sciences which can have many downstream societal impacts, but we do not otherwise see an immediate potential for societal impact. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).We will denote random variables by upper-case letters, and samples by their corresponding lower-case letter.Note that the injectivity is in the sense of f as a function of both u and z, which allows observations that have a shared z but differ by their value in u, and the function remains injective. For example, rotated images with the exact same content can have a shared z, but remain injective due to the rotation being captured in u.We computation of in-distribution metrics using random subsamples from the test set. The out-of-distribution metric was computed on a small dataset with a single perturbation and subsamples were not needed. 


Decision: This paper introduces a novel method for aligning unpaired multimodal data using label-based classifiers. The approach involves training classifiers to predict labels for each modality and using the resulting logits to calculate a similarity matrix between modalities. The method leverages a common latent space, influenced by a perturbation signal, to align samples across modalities. This alignment is further refined using optimal transport matching with entropic regularization. Empirical results on multiple datasets, including imaging and single-cell data, demonstrate that this method effectively improves sample alignment compared to existing techniques.

All reviewers agree that this work is interesting and important and the paper is well written and easy to follow. Experiments show that the proposed method outperforms the literature in terms of MSE and other specific metrics. I think this is an interesting paper that should be accepted for the conference.

Reviews: review: paper_summary: paper_summary: This paper presents a new method for pairing unpaired examples across different modalities using the labels of the examples. The method essentially trains a classifier for predicting labels for examples in each modality, then uses the classifier's logits across two modalities to calculate a similarity matrix. Empirical results show that the method is promising in aligning data from different modalities.
strength_weakness: strength_weakness: The method uses the information in the labels in order to match labels via classifiers. There are two drawbacks with this approach: (1) the quality of matching depends directly on the choice of the classifier and its capacity/complexity. (2) more importantly, in many multimodal representation learning settings the label signal is unavailable.
questions: questions: Is it possible that the use of labels as the sole signal to match examples from two modalities results in example clustering based on labels? i.e. examples from modality 1 that belong to class t are matched strongly against those from modality 2 that also belong to class t?
In that case, is the method still beneficial compared to not doing any matching (i.e. just using data from individual modalities)? Having a discussion on this in the paper will be useful.
limitations: limitations: There's a short section on limitations of the method, explaining (1) the reliance of the method on labels, and (2) the strong assumptions made in the method.
paper_summary: paper_summary: The paper proposes a novel matching method for unpaired multimodal (bimodal) data. It aligns unpaired samples across modalities using a propensity score. Based on additional treatment information, the propensity score is learned and used to align multimodal samples.
The method is evaluated on three datasets: one imaging dataset and one single-cell dataset (both with available matching ground truth); and one additional single-cell microscopy dataset without available ground truth.
strength_weakness: strength_weakness: - Final purpose of the method: To me, it was not fully clear what the final goal/purpose of the method is. Is it (1) the aligning of unpaired multimodal samples? Or (2), learning from unpaired, multimodal samples while aligning the samples? If (1) is the goal, what can we conclude from the improved results? If (2) is the goal, it would be interesting to see if the improved alignment/matching results in better downstream task/learning from multimodal data compared to other works (e.g., Yang et al., 2021) or just unimodal methods without any pairing of data.
- Evaluation: It is not straightforward to understand the used metrics (maybe also related to the point above) and relate the performance metrics to the quality of the method (besides a relative performance between different alignment methods). Maybe reporting some classification-based performances for alignment (is the assigned sample from the other modality the correct one) or some downstream task performance (what task can we solve using the trained method?) would help.
- Although the method is advertised as multimodal, no experiments or proofs use more than two modalities. Maybe the term bimodal would be more accurate? Otherwise, an experiment or derivation for more than two modalities would be appreciated.
questions: questions: - what are the labels/perturbations in the image-based dataset? I did not see the information in section 6 (but maybe I missed it. Apologies in that case)
- section 6.2: what is the reason for using the first 200 principal components for the gene expression modality (and not for the protein data as well)?
- does assumption A1 hold in the datasets used? Is there a way to check whether the assumption holds?
- there are some typos:
    - line 90: the our instead of the or our
    - line 168: form instead of from
    - line 227: against on instead of against or on
limitations: limitations: The authors address some of the limitations.
paper_summary: paper_summary: This work addresses the problem of aligning unpaired samples from multimodal data, whereby the problem is to find sample $x_i$ from modality $i$ that is best “related” to sample $x_j$ from modality $j$, when those two samples come from the observation of the same phenomenon according to two different modalities (e.g. sensors).

This is an important problem that arise in many endeavors including biology, and that has been somehow overlooked in the literature of multimodal representation learning, whereby large quantities of paired data (images and their corresponding captions) are available.
In a nutshell, the idea proposed in this work is as follows. Provided the following assumptions are valid, that is there exists a common latent space $Z$, which is entangled with a “perturbation signal” or label $t$, but that it is independent of modality specific noise $U_{i,j}$, then it is possible to estimate a common space for matching samples through the “bridge” across modalities offered by the perturbation signal $t$.

In practice, this amounts to training a classifier to predict $t$ given one of the available modality, which allows the creation of a transformed dataset endowed with a Euclidean distance. Such distance can be used to compute a matching matrix $M$, whereby entry $(i,j)$ represents the likelihood of $x_i$ being matched to $x_j$. Matching matrix $M$ can be obtained efficiently through the lenses of optimal transport matching with entropic regularization.
It should be noted that $M$ can then be used to construct pseudo-samples for any learning task that requires paired samples $(x_i, x_j)$ by replacing them with $(x_i, \hat{x}_j)$, 

where the second term is obtained through the application of the matching matrix to a structural equation $f_\theta$, which is a parametric function with learnable parameters.

Several experiments on synthetic data, and a dataset whereby ground truth pairing is available, as well as experiments where no ground truth is available, indicate that the proposed method, in combination with optimal transport matching, provides an improvement over state of the art matching methods.
strength_weakness: strength_weakness: * The proposed method relies on strong assumptions, which are only superficially discussed. Relaxation of the conditional independence assumption A1 is discussed in appendix A, but only for the case of a small perturbation set $t \in \{ 0,1 \}$, which allows the authors conclude that under exact optimal transport, this is equivalent to operate on order preserving effects of $t$. The validity of the method also  relies on A2, which allows computing propensity score from data alone, and not the shared latent space $Z$.

* The overall perturbation method, as well as the existence of a shared latent representation $Z$ is only vaguely described. I have looked up the key reference [Yang et al, 2021], but the lack of detail is uncomfortable. It is possible to get an intuition of how the method works, but I would prefer such important technical details to be spelled out clearly.

* There is no mention of the scalability of the proposed approach. Building a matching matrix can be a daunting task for some domains whereby a large amount of unpaired training data is available. In this case, the matching matrix could be either a big square matrix, or a very skinny and tall matrix that could pose some computational challenges. It would also be important to discuss about what would be the best operating scenario for the proposed method in terms of the number of samples per modality. In very unbalanced scenarios, would the classifiers using one or the other modality really behave similarly as per Equation (5)?
questions: questions: * Would it be possible to explain (and eventually add to the appendix) in more detail the setup from [Yang et al, 2021] which is used in this work? More precisely, it would be great to have the definition of the multimodal autoencoder, the existence of a shared latent space, and the details of the perturbations $t$ and their effects on $Z$. This would really help formalizing better the use of $t$ as a bridge between modalities.

* Can you provide a thorough discussion about the scalability of the proposed approach, the impact of imbalance in the size of the multimodal data (e.g. $n_1 >> n_2$), and the possible application to other domains with more abundant data?

* Can you elaborate more on the baselines you compare your method to? Lines 235-239 are not detailed enough to fully appreciate the differences w.r.t. the proposed method.

======= POST REBUTTAL UPDATE =======
Dear authors, thank you for the rebuttal. I have raised my score.
limitations: limitations: Yes, the authors acknowledge that their work relies on strong assumptions, but suggest that empirical results indicate their method is robust to failures. Nevertheless, I could not find the discussion on the experimental results that clearly indicates such robustness.

Comments: title: Thank you for your rebuttal
comment: Thank you so much!
Dear authors,
thank you for the thorough rebuttal to my observations and questions.

I appreciate the explanation on the assumptions made in this work, and how do they compare to the ones done in the state of the art.

I also appreciate the discussion on the obtained results, which surpass currently available methods.

I will raise my score to reflect these observations.
Dear authors,

There is a long line of work discussing truly multimodal methods, e.g. [1,2,3], which is important as many new problems arise when trying to scale methods from 2 to $M$ modalities.

But I acknowledge your efforts in explaining and answering my questions. I will raise my score.


--- 

### References
[1] Wu and Goodman, "Multimodal Generative Models for Scalable Weakly-Supervised Learning", Neurips 2018

[2] Shi et al., "Variational mixture-of-experts autoencoders for multi-modal deep generative models", Neurips 2019

[3] Sutter et al., "Generalized Multimodal ELBO", ICLR 2021
Likewise, thank you so much for the time you put into your review! Please let us know if there's anything else you need to know, and if you have any other concerns that are preventing you from raising your score.
Dear authors,

Thanks for your rebuttals and replies to my questions.