[a]David Kim for the Prometheus authors

# Prometheus: An Open-Source Neutrino Telescope Simulation

###### Abstract

The soon-to-be-realized, global network of neutrino telescopes will allow new
opportunities for collaboration between detectors. While each detector is
distinct, they share the same underlying physical processes and detection
principles. The full simulation chain for these telescopes is typically
proprietary which limits the opportunity for joint studies. This means there
is no consistent framework for simulating multiple detectors. To overcome
these challenges, we introduce Prometheus, an open-source simulation tool for
neutrino telescopes. Prometheus simulates neutrino injection and final state
and photon propagation in both ice and water. It also supports user-supplied
injection and detector specifications. In this contribution, we will introduce
the software; show its runtime performance; and highlight successes in
reproducing simulation results from multiple ice- and water-based
observatories.

##  1 Introduction

Figure 1:  Schematic showing the physical processes Prometheus models. (1),
Prometheus selects an interaction vertex within simulation volume, depicted
here by the lighter-colored region. (2), the final states of this interaction
are then propagated, accounting for energy losses and any daughter particles
which may be produced. (3), these losses are then converted to a number of
photons. (4), finally, these photons are then propagated until they either are
absorbed or reach an optical module.

The global network of neutrino telescopes, defined as gigaton-scale neutrino
detectors, has allowed us observe the Universe in new ways. The subset of this
network of telescope deployed on ice or water includes the IceCube Neutrino
Observatory [1] near the South Pole, proposed detectors ORCA and ARCA [2] in
the Mediterranean Sea (KM3NeT collaboration), and Baikal-GVD in Lake Baikal,
Russia [3] (BDUNT collaboration). Additionally, new experiments like P-ONE [4]
off the coast of Vancouver and TRIDENT [5] in the South China Sea are
underway, as well as an expansion for the IceCube Observatory [6, 7].

These telescopes all share many technological features. Each of these
detectors operate by detecting Cherenkov photons emitted by neutrino
interaction byproducts, and as such follow the same general simulation chain
illustrated in Fig. 1. The only proprietary step is the final detector
response, which occurs after an optical module (OM) has detected a photon. Yet
for many years we have lacked a simulation framework that takes advantage of
this similarity. Existing packages individually cover one or two of these
common steps, but until now there has been no easy way to simulate a particle
from injection to photon propagation.

Prometheus [8] looks to correct this by providing an integrated framework to
simulate these common steps for arbitrary detectors in ice and water, using a
combination of publicly available packages and those newly developed for this
work. Neutrino injection is handled by LeptonInjector[9], an event generation
recently developed by the IceCube Collaboration. Taus and muons are then
propagated by PROPOSAL [10]. Light yield simulation and photon propagation in
ice relies on PPC [11], while in water these steps are covered by Fennel [12]
and Hyperion, respectively.

Prometheus’s flexibility allows one to optimize detector configurations for
specific physics goals, while the common format allows one to develop
reconstruction techniques that may be applied across different experiments.
With the recent explosion in machine-learning research, it is now more
important than ever that we are able to rapidly implement and test new ideas
without relying on tools and data that may be proprietary to their
experiments.

The rest of this article is organized as follows. In Sec. 2 we outline the
format of Prometheus’s output and validate against published results. In Sec.
3 we present community work that employs Prometheus. In Sec. 4 we provide a
short example running Prometheus code. Finally, in Sec. 5 we conclude and
offer our future outlook.

##  2 Prometheus: Output and Validation

Figure 2: Event views for various detector geometries. This shows the events
created by either $\nu_{\mu}$ charged-current or $\nu_{e}$ charged-current
interactions in a variety of geometries of current and proposed neutrino
telescopes. Each black dot is an OM, while each colored dot indicates the
average time at which photons arrived at the OM; black indicates an earlier
arrival, orange indicates a later arrival, and purple an arrival in between.
Furthermore, the size of the colored spheres is proportional to the number of
photons that arrived at the OM. Detectors which appear against lighter blue
backgrounds—the top row—are ice-based, while those against the darker blue
backgrounds are water-based.

Prometheus outputs to Parquet [13] files that include two primary
fields—photons and mc_truth. photons contains information on photons that
produce hits in user-defined detection regions. This includes photon arrival
time, OM identification numbers, OM position, and the final-state particle
that produced the photon. Fig. 1 shows event displays for various detectors
generated using the information in photons. mc_truth includes information on
the injection like the interaction vertex; the initial neutrino type, energy,
and direction; and the final state types, energies, directions, and parent
particles. Users may also save the configuration file (see Sec. 4) as a json
file. This allows the user to resimulate events using the same parameters,
which is useful for comparing the same event across multiple detector
geometries.

The information stored in mc_truth allows us to compute effective areas when
combined with the weights from LeptonWeighter. Since it mainly depends on the
physics implemented in Prometheus—such as neutrino-nucleon interactions,
lepton range, and photon propagation—effective area serves as a reliable
indicator of our code’s performance when all simulation steps are properly
integrated. We can then validate our simulation by comparing these effective
areas to published data. Fig. 3 compares different experiments published
effective areas to our estimation using Prometheus simulations. It is worth
noting that the calculations for effective area rely on detector-specific cuts
and OM response, to which we have limited or no access. We can therefore
expect differences of $\mathcal{O}(10\%)$ from these missing detector details.

Figure 3:  Effective area computed using Prometheus with comparisons to
published results. We compare the $\nu_{\mu}$ effective areas computed with
Prometheus for IceCube, P-ONE3, ORCA, and ARCA for three different hit
requirements, denoted by different line styles, to published effective areas.
The IceCube effective area, taken from Ref. [14], is for
$\nu_{\mu}+\nu_{\tau}$ events which pass the SMT-8 trigger and agrees with our
calculation to within uncertainties. The ARCA [2] and ORCA [15] cases
effective areas are constructed with more complicated hit requirements. Still,
the scale and shape of the ORCA and ARCA effective areas and the Prometheus
effective areas agree within uncertainties despite the simplified selection
criterion. As of the publication of this proceeding, there is no published
effective area for P-ONE3.

##  3 Community Contributions

Figure 4:  [16, 17, 18] Event rates of triggers in different neutrino
telescopes compared to the run-times of various reconstruction methods.
Notably, sparse submanifold CNNs can process events well above standard
trigger rates in both ice- and water-based experiments. The CNN and maximum
likelihood method run-times are taken from [19]. Reproduced with permission of
the authors from Ref. [20].  Figure 5:  Angular reconstruction performance as
a function of the true neutrino energy. The angular resolution results are
binned by the true neutrino energy, with the median taken from each bin to
form the lines shown [20].

Developing a published reconstruction is often a long-term effort internal to
an experiment, meaning methods are slow to implement and difficult to compare.
Prometheus remedies this by allowing users to readily generate the large data
sets necessary to test the feasibility and relative performance of different
reconstruction methods. In this section we highlight two submissions at ICRC
2023 that utilize Prometheus in such a way. The first of these works is a
software-focused effort that looks to improve the speed of first
reconstructions at detectors and the second is a hardware-focused effort
concerning GPU alternatives for low power computing.

Ref. [20], proposes sparse-submanifold convolutions (SSCNNs) as an alternative
to the convolutional neural network (CNN) and traditional trigger-level event
reconstructions currently used in neutrino telescopes. Their directional and
energy algorithm is trained on a data set of 412892 events and tested on a
further 50000 events, cut from a set of 3 million generated by Prometheus. As
seen in Fig. 4, the SSCNN is capable of running at speeds comparable to the
neutrino telescope trigger rates. In Fig. 5, it achieves a median angular
resolution below $4^{\circ}$ for the highest-energy trigger-level events,
which matches or outperforms currently employed reconstructions. Their work
could be applied to improve on-site reconstructions and filtering for notable
events at any detector.

Ref. [21], is another Prometheus user submission that explores hardware
accelerators to improving event reconstruction. Specifically, they look at how
a Tensor Processing Unit (TPU) compatible algorithm could lower energy
consumption while running comparable operations to the current GPU-based ones.
Their model reaches a median angular resolution around $5^{\circ}$ above
$10^{3}$ GeV. Meanwhile, approximate peak total power consumption drops from
100W on the best performing GPU (Apple M1 Pro chip) down to 3W on the Google
Edge TPU. The portion of the power consumption directed to the ML accelerator
component drops from 15W to 2W. This submission illustrates Prometheus’s
utility in testing more speculative ideas and providing proof-of-concepts to
encourage further work, particularly in the direction of machine learning.

The findings of both works are generic to any ice or water neutrino telescope.
Their models are trained and tested on example detector geometries rather than
the geometry of any existing or proposed detector. Since Prometheus is fully
open-source, we hope it will facilitate not only easier iteration on these
techniques but also greater collaboration in sharing and adapting them.

##  4 Examples

In this section we will walk through producing a simulation of $\nu_{\mu}$
charged-current interactions in an ice-based detector. This example will use
mostly default values for injection parameters to show the essential steps in
running a simulation, after which we will show how the user can set these
parameters.

In this example, we set the number of events to simulate, the detector
geometry file, and—using the final state particles—the event type. With just
this, we can simulate events with Prometheus!

⬇

1import prometheus

2from prometheus import config, Prometheus

3

4config["run"]["nevents"] = 100

5geofile = f"{resource_dir}/geofiles/demo_ice.geo"

6config["detector"]["geo file"] = geofile

7injection_config["simulation"]["final state 1"] = "MuMinus"

8injection_config["simulation"]["final state 2"] = "Hadrons"

9

10p = Prometheus(config)

11p.sim()

As briefly shown here, the config dictionary is our primary interface for
configuring Prometheus. Key parameters not shown here include the output
directory; random state seed; and injection parameters like injection angle
and energy. Information on a detector’s medium, either ice or water, is stored
in its geometry file. The geofiles directory in the GitHub repository has
geometry files for all of the detectors in Fig. 2. Again, this only scratches
the surface of Prometheus’s capabilities. For a more thorough description of
all the options and features available see Ref. [8], which has more in-depth
examples for $\nu_{\mu}$ charged-current events in ice and $\bar{\nu}_{e}$
neutral-current events in water as well as directions on constructing a
detector, weighting events, and getting event rates.

##  5 Conclusion

In this submission we have introduced Prometheus as an open-source software
package for simulating neutrino telescopes. We have provided a brief example
for simulating events in an ice-based detector and highlighted two current
works that employ Prometheus. Prometheus’s flexibility of input for detector
geometry and injection parameters allows it to handle simulation for the full
range of existing and proposed telescopes in both water and ice.

As we have demonstrated in the highlighted community contributions, Prometheus
facilitates the implementation of new ideas without the need for proprietary
data, or for data on a scale not yet available. Via its particular application
to machine-learning models, Prometheus can be a key piece in accelerating the
development of faster, more efficient reconstructions for all detectors.

Finally, it is our hope that Prometheus opens the door for greater
collaboration within the community. By encouraging the sharing of methods and
simulated data sets, we hope work done by any one effort more quickly and
easily becomes progress for every group in the global neutrino telescope
network.

##  6 Acknowledgements

We would like to thank all users who tested early versions of this software,
including—in no particular order—Miaochen Jin, Eliot Genton, Tong Zhu, Rasmus
Ørsøe, Savanna Coffel, and Felix Yu. The authors that developed Prometheus
were supported by Faculty of Arts and Sciences of Harvard University, the
Alfred P. Sloan Foundation, NSF under grants PLR-1600823, PHY-1607644,
Wisconsin Research Council with funds granted by the Wisconsin Alumni Research
Foundation, Australian Government through the Australian Research Council’s
Discovery Projects funding scheme (project DP220101727), and Lynne Sacks and
Paul Kim.

## References

  * [1] IceCube Collaboration, M. G. Aartsen et al. JINST 12 no. 03, (2017) P03012.
  * [2] KM3Net Collaboration, S. Adrian-Martinez et al. J. Phys. G 43 no. 8, (2016) 084001.
  * [3] A. D. Avrorin et al. Phys. Part. Nucl. 46 no. 2, (2015) 211–221.
  * [4] P-ONE Collaboration, M. Agostini et al. Nature Astron. 4 no. 10, (2020) 913–915.
  * [5] Z. P. Ye et al. arXiv:2207.04519.
  * [6] IceCube Collaboration, A. Ishihara PoS ICRC2019 (2021) 1031.
  * [7] IceCube-Gen2 Collaboration, M. G. Aartsen et al. J. Phys. G 48 no. 6, (2021) 060501.
  * [8] J. Lazar, S. Meighen-Berger, C. Haack, D. Kim, S. Giner, and C. A. Argüelles.
  * [9] IceCube Collaboration, R. Abbasi et al. Comput. Phys. Commun. 266 (2021) 108018.
  * [10] J. H. Koehne, K. Frantzen, M. Schmitz, T. Fuchs, W. Rhode, D. Chirkin, and J. Becker Tjus Comput. Phys. Commun. 184 (2013) 2070–2090.
  * [11] D. Chirkin, “ppc.” https://github.com/icecube/ppc, 2022.
  * [12] S. Meighen-Berger, “Fennel: Light from tracks and cascades,” 2022.  https://github.com/MeighenBergerS/fennel.
  * [13] “Documentation | Apache Parquet.” https://parquet.apache.org/docs/.  Accessed: 20123-02-17.
  * [14] IceCube Collaboration, A. Karle PoS ICRC2009 (2010) .
  * [15] KM3NeT Collaboration, G. de Wasseige, A. Kouchner, M. Colomer Molla, D. Dornic, and S. Hallmann PoS ICRC2019 (2020) 934.
  * [16] IceCube, MAGIC, VERITAS Collaboration, M. G. Aartsen et al. JINST 11 no. 11, (2016) P11009.
  * [17] B. Bakker, Trigger studies for the Antares and KM3NeT neutrino telescopes.  PhD thesis, University of Amsterdam, 2011.
  * [18] A. Avrorin et al. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 639 no. 1, (2011) 30–32. Proceedings of the Seventh International Workshop on Ring Imaging Cherenkov Detectors.
  * [19] M. Hünnefeld, Online Reconstruction of Muon-Neutrino Events in IceCube using Deep Learning Techniques.  PhD thesis, Technische Universität Dortmund, 2017.
  * [20] F. J. Yu, J. Lazar, and C. A. Argüelles PoS ICRC2023 1004\.
  * [21] M. Jin, Y. Hu, and C. A. Argüelles.


# Dynamic Gesture Recognition in Ultra-Range Distance for Effective Human-
Robot Interaction
††thanks:

1st Eran Bamani Beeri  School of Mechanical Engineering
Tel-Aviv University
Tel-Aviv, Israel 6997801
<EMAIL_ADDRESS>2nd Eden Nissinman  School of Mechanical
Engineering
Tel-Aviv University
Tel-Aviv, Israel 6997801
<EMAIL_ADDRESS>3rd Avishai Sintov  School of Mechanical Engineering
Tel-Aviv University
Tel-Aviv, Israel 6997801
<EMAIL_ADDRESS>
###### Abstract

This paper presents a novel approach for ultra-range gesture recognition,
addressing Human-Robot Interaction (HRI) challenges over extended distances.
By leveraging human gestures in video data, we propose the Temporal-
Spatiotemporal Fusion Network (TSFN) model that surpasses the limitations of
current methods, enabling robots to understand gestures from long distances.
With applications in service robots, search and rescue operations, and drone-
based interactions, our approach enhances HRI in expansive environments.
Experimental validation demonstrates significant advancements in gesture
recognition accuracy, particularly in prolonged gesture sequences.

###### Index Terms:

Human-Robot Collaboration, Human-Robot Interaction, Video Recognition, Ultra-
Range Gesture Recognition

##  I Introduction

In recent years, Gesture Recognition (GR) in Human-Robot Interaction (HRI) has
been advanced significantly. While research primarily focuses on hand gesture
recognition [1] using conventional cameras and Kinect devices [2], a critical
gap remains in recognizing gestures over distances beyond seven meters [3, 4].
Bridging this gap is crucial for seamless collaboration between humans and
robots, especially in scenarios like surveillance, drone operations, and
search and rescue missions.

Figure 1: Six human gestures are depicted with two images each, illustrating
the start and end of each gesture. The gestures, arranged from the top left,
include beckoning, stop, null, thumbs-up, pointing, and thumbs-down.

Ultra-range gesture recognition (URGR) aims to extend recognition capabilities
to 25 meters [5], yet considers only static images as the majority of gesture
recognition approaches. However, some gestures are inherently dynamic such as
waving and beckoning. Additionally, some dynamic gestures may appear as
different if only a single frame is used for inference. For instance, a
snapshot of a waving gesture may look like a stop gesture. In addition, URGR
often faces challenges like image degradation and noise. Video data offers
promise, capturing motion dynamics better [6, 7]. Deep learning shows
potential for recognizing gestures in videos [8, 9], aiding robots in
understanding gestures across distances.

In this paper, we introduce the Temporal-Spatiotemporal Fusion Network (TSFN)
for URGR in HRI which is effective within a range of 28 meters and operates in
real-time. The proposed TSFN model uses Temporal Convolutional Networks (TCN)
[11] for temporal dependencies and R(2+1)D convolutional networks [10] for
spatiotemporal feature extraction, advancing ultra-range human-robot
interaction.

##  II Methods

###  II-A Problem Formulation and Data Collection

GR across diverse environments and distances is essential for effective human-
robot communication. Let $\boldsymbol{V}_{i}\in\mathbb{R}^{T\times H\times
W\times 3}$ denote video sequences up to 28 meters, with $T$, $H$, and $W$
representing frame quantity, height and width, respectively. GR in a video
faces challenges due to temporal dynamics, occlusions, and lighting
variations. Given a temporal frame sequence
$\boldsymbol{V}=\\{v_{1},v_{2},...,v_{T}\\}$ and a set of gestures $G$, the
problem entails finding the function $F:V\rightarrow G$ mapping video
sequences to human gestures.

A diverse dataset $\mathcal{D}=\\{(X_{i},y_{i},d_{i})\\}_{i=1}^{N}$ is crucial
for accurate gesture recognition. It comprises input video $X_{i}$, gesture
label $y_{i}$, and distance $d_{i}\in[4,28]$ meters from the camera. Webcam-
collected data offers varied scenarios, improving model robustness. Challenges
include pose/environmental diversity and image quality degradation with
distance, leading to reduced resolution and motion blur. Annotations include
ground truth gesture labels and distance, aiding model training/evaluation.

###  II-B Models

Our TSFN model combines TCN and R(2+1)D models to create a robust GR system
from video sequences. TCN captures temporal dependencies and the output of a
TCN layer is given by

$Y_{\text{TCN}}=\sigma(W_{\text{TCN}}\ast V+b_{\text{TCN}}),\vspace{-0.1cm}$
(1)

where $W_{\text{TCN}}$ and $b_{\text{TCN}}$ are the weights and biases of the
TCN layer, the operator $\ast$ denotes convolution, and $\sigma$ is the
Sigmoid activation function. Video data necessitates 3D convolutions. A
standard 3D convolution applies a filter $W\in\mathbb{R}^{K\times d\times
k\times k}$, where $K$ is the number of output channels, $d$ is the temporal
depth, and $k\times k$ is the spatial dimension, given by:

$Y_{t,h,w}^{k}=\sum_{c=1}^{C}\sum_{i=1}^{d}\sum_{j=1}^{k}\sum_{m=1}^{k}W_{c,i,j,m}^{k}\cdot
V_{c,t+i,h+j,w+m}.$

In the TCN model, the 3D convolution is factorized into 1D temporal
convolution. R(2+1)D networks extend two-dimensional convolutional operations
into the temporal domain, separating spatial and temporal convolutions. The
R(2+1)D convolution is expressed as

$Y_{\text{R(2+1)D}}=\sigma(W_{\text{spatial}}\ast(\sigma(W_{\text{temporal}}\ast
V+b_{\text{temporal}}))+b_{\text{spatial}}),$ (2)

where $W_{\text{spatial}}$ and $W_{\text{temporal}}$ are the spatial and
temporal convolutional weights, respectively. Finally. the outputs from the
TCN and R(2+1)D layers are concatenated to form combined features

$F_{\text{combined}}=[Y_{\text{TCN}},Y_{\text{R(2+1)D}}].\vspace{-0.1cm}$ (3)

These features are then passed through three fully connected layers to produce
the final gesture classification:

$\text{O}=\sigma(W_{\text{fc}}\cdot
F_{\text{combined}}+b_{\text{fc}}),\vspace{-0.1cm}$ (4)

where $W_{\text{fc}}$ and $b_{\text{fc}}$ are the weights and biases of the
fully connected layer. This fusion leverages the strengths of both models to
enhance GR accuracy.

To train our TSFN model to address the challenges of recognizing human
gestures in a video, particularly when the subject is far from the camera, we
propose a novel loss function that takes into account the degradation and
other phenomena affecting the results. This novel loss function will be
designed to balance classification accuracy with robustness to distance-
related degradation and is given by

$\mathcal{L}=\mathcal{L}_{CE}+\alpha\cdot\mathcal{L}_{global}+\beta\cdot\mathcal{L}_{dist}+\gamma\cdot\mathcal{L}_{robust}$

, where $\mathcal{L}_{CE}$ is the cross-entropy loss and
$\mathcal{L}_{global}$ is the global context loss that ensures the global
context of the video is captured accurately. Loss $\mathcal{L}_{dist}$ is the
distance aware loss that takes into account the degradation effect due to
varying distances given by

$\mathcal{L}_{dist}=\frac{1}{N}\sum_{i=1}^{N}d_{i}\cdot\mathcal{L}_{CE}(v_{i},y_{i}).\vspace{-0.2cm}$

$\mathcal{L}_{dist}$ penalizes more easily for errors when the subject is
farther from the camera. Robustness loss $\mathcal{L}_{robust}$ enhances the
model’s robustness to various degradations by adding a regularization
component that minimizes the variance in predictions due to changes in
distance:

$\vspace{-0.3cm}\mathcal{L}_{robust}=\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{M}(\mathcal{L}_{CE}(v_{i}^{(j)},y_{i})-\bar{\mathcal{L}}_{CE}(y_{i}))^{2}.\vspace{-0.1cm}$

Component $\bar{\mathcal{L}}_{CE}$ is the mean cross-entropy loss for the true
label, $\alpha,\beta,\gamma$ are hyper-parameters that control the relative
importance of each loss component.

##  III Experimental Results

To validate the effectiveness of the proposed TSFN model, we conducted
experiments comparing its performance with state-of-the-art video recognition
models. Figure 1 illustrates the six gesture categories we focused on:
pointing, thumbs-up, thumbs-down, beckoning, stop, and null gestures. The
models were tested on a dataset of human gestures observed from 4 to 28
meters. We collected 15 samples of 4-second videos per meter, totaling 2,250
videos for each gesture, with 150 videos used as the test set. Key performance
metrics including accuracy, loss, and Mean Average Precision (mAP) are shown
in Table I. Figure 2 shows the accuracy of the TSFN in GR with respect to
different distances. The results clearly show the superiority of TSFN over
existing methods.

TABLE I: Comparison of Video Recognition Models Model |  | Accuracy (%) | Loss | mAP
---|---|---|---|---
ViViT | [13] | 91.2 | 0.29 | 0.87
TCN | [11] | 89.3 | 0.37 | 0.79
R(2+1)D | [10] | 88.1 | 0.39 | 0.77
Vanilla CNN |  | 78.4 | 0.51 | 0.70
CNN+LSTM | [12] | 83.7 | 0.42 | 0.74
TSFN |  | 96.1 | 0.12 | 0.92
Figure 2: Model performance vs. distance. The plot shows the accuracy of the
TSFN model in recognizing gestures at various distances.

##  IV Conclusion

TSFN redefines ultra-range gesture recognition, which is crucial for seamless
Human-Robot Collaboration. By fusing TCNs and R(2+1)D networks, TSFN excels in
understanding gestures over vast distances. Future work includes refining TSFN
for real-time applications and exploring its potential in dynamic environments
including underwater.

## References

  * [1] S. Ahmed, K. D. Kallu, S. Ahmed, and S. H. Cho, ”Hand gestures recognition using radar sensors for human-computer-interaction: A review,” Remote Sensing, vol. 13, no. 3, p. 527, 2021.
  * [2] K. K. Biswas and S. K. Basu, ”Gesture recognition using Microsoft Kinect®,” in Proceedings of the 5th International Conference on Automation, Robotics, and Applications, pp. 100–103, IEEE, 2011.
  * [3] K. Nickel and R. Stiefelhagen, ”Visual recognition of pointing gestures for human-robot interaction,” Image and Vision Computing, vol. 25, no. 12, pp. 1875–1884, 2007.
  * [4] J. P. Wachs, M. Kölsch, H. Stern, and Y. Edan, ”Vision-based hand-gesture applications,” Communications of the ACM, vol. 54, no. 2, pp. 60–71, 2011.
  * [5] E. Bamani, E. Nissinman, I. Meir, L. Koenigsberg, and A. Sintov, ”Ultra-Range Gesture Recognition using a web-camera in Human–Robot Interaction,” Engineering Applications of Artificial Intelligence, vol. 132, p. 108443, 2024.
  * [6] J. Materzynska, G. Berger, I. Bax, and R. Memisevic, ”The jester dataset: A large-scale video dataset of human gestures,” in Proceedings of the IEEE/CVF International Conference on Computer Vision Workshops, 2019, pp. 0–0.
  * [7] Y. M. Lui, ”Human gesture recognition on product manifolds,” The Journal of Machine Learning Research, vol. 13, no. 1, pp. 3297–3321, 2012.
  * [8] K. Yang, R. Li, P. Qiao, Q. Wang, D. Li, and Y. Dou, ”Temporal pyramid relation network for video-based gesture recognition,” in 2018 25th IEEE International Conference on Image Processing (ICIP), 2018, pp. 3104–3108.
  * [9] X. Bu, ”Human motion gesture recognition algorithm in video based on convolutional neural features of training images,” IEEE Access, vol. 8, pp. 160025–160039, 2020.
  * [10] D. Tran et al., ”A closer look at spatiotemporal convolutions for action recognition,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2018.
  * [11] C. Lea et al., ”Temporal convolutional networks: A unified approach to action segmentation,” in Computer Vision–ECCV 2016 Workshops: Amsterdam, The Netherlands, October 8-10 and 15-16, 2016, Proceedings, Part III, vol. 14, pp. 47-54, Springer International Publishing, 2016.
  * [12] J. Chen, J. Wang, Q. Yuan, and Z. Yang, ”CNN-LSTM model for recognizing video-recorded actions performed in a traditional Chinese exercise,” IEEE Journal of Translational Engineering in Health and Medicine, 2023.
  * [13] A. Arnab, M. Dehghani, G. Heigold, C. Sun, M. Lučić, and C. Schmid, ”ViViT: A video vision transformer,” in Proceedings of the IEEE/CVF International Conference on Computer Vision, pp. 6836–6846, 2021.


# Simultaneous Dimensionality Reduction: A Data Efficient Approach for
Multimodal Representations Learning

Eslam Abdelaleem1 Ahmed Roman1 K. Michael Martini1,3 Ilya Nemenman1,2,3
1Department of Physics, Emory University, Atlanta, GA 30322
2Department of Biology, Emory University, Atlanta, GA 30322
3Initiative in Theory and Modeling of Living Systems, Emory University,
Atlanta, GA 30322
{eslam.abdelaleem, ahmed.hemdan.roman, karl.michael.martini,
<EMAIL_ADDRESS>Corresponding authorCurrently
<EMAIL_ADDRESS>at Department of Medical Oncology - Dana-Farber
Cancer Institute - Boston, MA, USA & Broad Institute of MIT and Harvard -
Cambridge, MA, USA & Harvard Medical School - Boston, MA, USA

###### Abstract

Current experiments frequently produce high-dimensional, multimodal
datasets—such as those combining neural activity and animal behavior or gene
expression and phenotypic profiling—with the goal of extracting useful
correlations between the modalities. Often, the first step in analyzing such
datasets is dimensionality reduction. We explore two primary classes of
approaches to dimensionality reduction (DR): Independent Dimensionality
Reduction (IDR) and Simultaneous Dimensionality Reduction (SDR). In IDR
methods, of which Principal Components Analysis is a paradigmatic example,
each modality is compressed independently, striving to retain as much
variation within each modality as possible. In contrast, in SDR, one
simultaneously compresses the modalities to maximize the covariation between
the reduced descriptions while paying less attention to how much individual
variation is preserved. Paradigmatic examples include Partial Least Squares
and Canonical Correlations Analysis. Even though these DR methods are a staple
of statistics, their relative accuracy and data set size requirements are
poorly understood. We introduce a generative linear model to synthesize
multimodal data with known variance and covariance structures to examine these
questions. We assess the accuracy of the reconstruction of the covariance
structures as a function of the number of samples, signal-to-noise ratio, and
the number of varying and covarying signals in the data. Using numerical
experiments, we demonstrate that linear SDR methods consistently outperform
linear IDR methods and yield higher-quality, more succinct reduced-dimensional
representations with smaller datasets. Remarkably, regularized CCA can
identify low-dimensional weak covarying structures even when the number of
samples is much smaller than the dimensionality of the data, which is a regime
challenging for all dimensionality reduction methods. Our work corroborates
and explains previous observations in the literature that SDR can be more
effective in detecting covariation patterns in data. These findings suggest
that SDR should be preferred to IDR in real-world data analysis when detecting
covariation is more important than preserving variation.

##  1 Introduction

Many modern experiments across various fields generate massive multimodal data
sets. For instance, in neuroscience, it is common to record the activity of a
large number of neurons while simultaneously recording the resulting animal
behavior (Stringer et al., 2019; Steinmetz et al., 2021; Urai et al., 2022;
Krakauer et al., 2017). Other examples include measuring gene expressions of
thousands of cells and their corresponding phenotypic profiles, or integrating
gene expression data from different experimental platforms, such as RNA-Seq
and microarray data (Clark et al., 2013; Zheng et al., 2017; Svensson et al.,
2018; Huntley et al., 2015; Lorenzi et al., 2018). In economics, important
variables such as inflation are often measured using combinations of
macroeconomic indicators as well as indicators belonging to different economic
sectors (Gosselin & Tkacz, 2001; Baillie et al., 2002; Freyaldenhoven, 2022;
Rudd, 2020). In all of these examples, an important goal is to estimate
statistical correlations among the different modalities.

Analyses usually begin with dimensionality reduction (DR) into a smaller and
more interpretable representation of the data. We distinguish two types of DR:
independent (IDR) and simultaneous (SDR) (Martini & Nemenman, 2023). In the
former, each modality is reduced independently, while aiming to preserve its
variation, which we call self signal. In the latter, the modalities are
compressed simultaneously, while maximizing the covariation (or the shared
signal) between the reduced descriptions and paying less attention to
preserving the individual variation. It is not clear if IDR techniques, such
as the Principal Components Analysis (PCA) (Hotelling, 1933), are well-suited
for extracting shared signals since they may overlook features of the data
that happen to be of low variance, but of high covariance (Colwell et al.,
2014). In particular, poorly sampled weak shared signals, common in high-
dimensional datasets, can exacerbate this issue. SDR techniques, such as
Partial Least Squares (PLS) (Wold et al., 2001) and Canonical Correlations
Analysis (CCA) (Hotelling, 1936), are sometimes mentioned as more accurate in
detecting weak shared signal (Chin & Newsted, 1999; Hair et al., 2011;
Pacharawongsakda & Theeramunkong, 2016). However, the relative accuracy and
data set size requirements for detecting the shared signals remain poorly
understood for both classes of methods.

In this study, we aim to assess the strengths and limitations of linear IDR,
represented by PCA, and linear SDR, exemplified by PLS and CCA, in detecting
weak shared signals. For this, we introduce a generative linear model that
captures key features of relevant examples, including noise, the self signal,
and the shared signal components. Using this model, we analyze the performance
of the methods in different conditions. Our goal is to assess how well these
techniques can (i) extract the relevant shared signal and (ii) identify the
dimensionality of the shared and the self signals from noisy, undersampled
data. We investigate how the signal-to-noise ratios, the dimensionality of the
reduced variables, and the method of computing correlations combine with the
sample size to determine the quality of the DR. We propose best practices for
achieving high-quality reduced representations with small sample sizes using
these linear methods.

##  2 Model

###  2.1 Relations to Previous Work

The extraction of signals from large-dimensional data sets is a challenging
task when the number of observations is comparable to or smaller than the
dimensionality of the data. The undersampling problem introduces spurious
correlations that may appear as signals, but are, in fact, just statistical
fluctuations. This poses a challenge for DR techniques, as they may retain
unnecessary dimensions or identify noise dimensions as true signals.

Here, we focus exclusively on linear DR methods. For these, the Marchenko-
Pastur (MP) distribution of eigenvalues of the covariance matrix of pure noise
derived using the Random Matrix Theory methods (Marchenko & Pastur, 1967) has
been used to introduce a cutoff between noise and true signal in real
datasets. However, recent work (Fleig & Nemenman, 2022) has shown that, when
observations are a linear combination of uncorrelated noise and latent low-
dimensional self signals, then the self signals alter the distribution of
eigenvalues of the sampling noise, questioning the validity of this naive
approach.

Moving beyond a single modality, Bouchaud et al. (2007) calculated the
singular value spectrum of cross-correlations between two nominally
uncorrelated random signals. However, it remains unknown whether the linear
mixing of self signals and shared signals affects the singular value spectra
of noise, and how all of these components combine to limit the ability to
detect shared signals between two modalities from data sets of realistic
sizes. Filling in this gap using numerical simulations is the main goal of
this paper, and analytical treatment of this problem will be left for the
future.

The linear model and linear DR approaches studied here do not capture the full
complexity of real-world data sets and state-of-the-art algorithms. However,
if sampling issues and self signals limit the ability of linear DR methods to
extract shared signals, it would be surprising for nonlinear methods to
succeed in similar scaling regimes on real data. Thus extending the previous
work to explicitly study the effects of linear mixtures of self signals,
shared signals, and noise on limitations of DR methods is likely to endow us
with intuition that is useful in more complex scenarios routinely encountered
in different domains of science.

Examples of scenarios with shared and self signals include inference of
dynamics of a system through a latent space (Creutzig et al., 2009; Chen et
al., 2022), where shared signals correspond to latent factors that are
relevant for predicting the future of the system from its past, while self
signals correspond to nonpredictive variation (Bialek et al., 2001). In
economics, shared and self signals correspond to diverse macroeconomic
indicators that are grouped into correlated distinct categories in structural
factor models (Forni & Gambetti, 2010; Gosselin & Tkacz, 2001; Rudd, 2020;
Baillie et al., 2002). In neuroscience, shared signals can correspond to the
latent space, by which neural activity affects behavior, while self signals
encode neural activity that does not manifest in behavior and behavior that is
not controlled by the part of the brain being recorded from (Sponberg et al.,
2015; Stringer et al., 2019; Natraj et al., 2022; Sani et al., 2021; Pang et
al., 2016; Urai et al., 2022; Krakauer et al., 2017).

Interestingly, in the context of the neural control of behavior, it was
noticed that SDR reconstructs the shared neuro-behavioral latent space more
efficiently and using a smaller number of samples than IDR (Sani et al.,
2021). Similar observations have been made in more general statistical
contexts (Chin & Newsted, 1999; Hair et al., 2011; Pacharawongsakda &
Theeramunkong, 2016), though the agreement is not uniform (Goodhue et al.,
2006; 2012; 2013). Because of this, most practical recommendations for
detecting shared signals are heuristic (Hair Jr et al., 2021), with widely
acknowledged, but poorly understood limitations and possible resolutions (Kock
& Hadaya, 2018). Our goal is to ground such rules in numerical simulations and
scaling arguments.

###  2.2 Linear Model with Self and Shared Signals

We consider a linear model with noise, $m_{\text{self,X}},m_{\text{self,Y}}$
self signals that are relevant to each modality independently, as well as
$m_{\text{shared}}$ shared signals that capture the interrelationships between
modalities. It results in $T$ observations of two high-dimensional
standardized observables, $X$ and $Y$ :

$\displaystyle\left[\tilde{X}\in\mathbb{R}^{N_{X}}\right]$
$\displaystyle=\underbrace{R_{X}}_{\text{Independent white
noise}}+\underbrace{U_{X}V_{X}}_{\text{Self-Signal for
X}}+\underbrace{PQ_{X}}_{\text{Shared-Signal}},$
$\displaystyle\left[\tilde{Y}\in\mathbb{R}^{N_{Y}}\right]$
$\displaystyle=\underbrace{R_{Y}}_{\text{Independent white
noise}}+\underbrace{U_{Y}V_{Y}}_{\text{Self-Signal for
Y}}+\underbrace{PQ_{Y}}_{\text{Shared-Signal}},$ (1) $\displaystyle
X=\tilde{X}/\sigma_{\tilde{X}},Y=\tilde{Y}/\sigma_{\tilde{Y}}.$ (2)

The observations of $X$ and $Y$ are linear combinations of the following: (a)
Independent white noise components $R_{X}$ and $R_{Y}$ with variances
$\sigma^{2}_{R_{X}}$ and $\sigma^{2}_{R_{Y}}$. (b) Self-signal components
$U_{X}$ and $U_{Y}$ residing in lower-dimensional subspaces
$\mathbb{R}^{m_{\mathrm{self},X}}$ and $\mathbb{R}^{m_{\mathrm{self},Y}}$ with
variances $\sigma^{2}_{U_{X}}$ and $\sigma^{2}_{U_{Y}}$. (c) Shared-signal
components $P$ in a shared lower-dimensional subspace
$\mathbb{R}^{m_{\mathrm{shared}}}$ with variance $\sigma^{2}_{P}$. These
components are projected into their respective high-dimensional spaces
$\mathbb{R}^{N_{X}}$ and $\mathbb{R}^{N_{Y}}$ using fixed quenched projection
matrices $V_{X}$, $V_{Y}$, $Q_{X}$, and $Q_{Y}$ with specified variances
$\sigma^{2}_{V_{X}}$, $\sigma^{2}_{V_{Y}}$, $\sigma^{2}_{Q_{X}}$, and
$\sigma^{2}_{Q_{Y}}$, all respectively. Entries in these matrices are drawn
from a Gaussian distribution with a zero mean and the corresponding variances.
Further, division by $\sigma_{\tilde{X}}$ and $\sigma_{\tilde{Y}}$
standardizes each column of the data matrices by their empirical standard
deviations. The total variance in the matrix $\tilde{X}$ can be calculated as
the sum of the variances of its individual components:
$\sigma^{2}_{\tilde{X}}=\sigma^{2}_{R_{X}}+m_{\mathrm{self},X}\times\sigma^{2}_{U_{X}}\sigma^{2}_{V_{X}}+m_{\mathrm{shared}}\times\sigma^{2}_{P}\sigma^{2}_{Q_{X}}$.
A similar calculation can be done for the total variance in $\tilde{Y}$.

We define self and shared signal-to-noise ratios
$\gamma_{\mathrm{self},X/Y},\gamma_{\mathrm{shared},X/Y}$ as the relative
strength of signals compared to background noise per component in each
modality. These definitions allow us to examine how easily self or shared
signals in each dimension can be distinguished from the noise.

$\gamma_{\mathrm{self},X/Y}=\frac{\sigma^{2}_{U_{X/Y}}\sigma^{2}_{V_{X/Y}}}{\sigma^{2}_{R_{X/Y}}},\quad\gamma_{\mathrm{shared},X/Y}=\frac{\sigma^{2}_{P}\sigma^{2}_{Q_{X/Y}}}{\sigma^{2}_{R_{X/Y}}}$
(3)

Our main goal is to evaluate the ability of linear SDR and IDR methods to
reconstruct the shared signal $P$ between $X$ and $Y$ while overlooking the
effects of the self signals $U_{X},U_{Y}$ on the statistics of the shared
ones.

##  3 Methods

We apply DR techniques to $X$ and $Y$ to obtain their reduced dimensional
forms $Z_{X}$ and $Z_{Y}$, respectively. $Z_{X},Z_{Y}$ are of sizes that can
range from $T\times 1$ to $T\times N_{X}$ and $T\times N_{Y}$, respectively.
As an IDR method, we use PCA (Hotelling, 1933). As SDR methods, we apply PLS
(Wold et al., 2001) and CCA (Hotelling, 1936; Vinod, 1976; Årup Nielsen et
al., 1998), including both normal and regularized versions of the latter. Each
of these methods focuses on specific parts of the overall covariance matrix

$C_{X,Y}=\begin{bmatrix}C_{XX}&C_{XY}\\\
C_{YX}&C_{YY}\end{bmatrix}=\begin{bmatrix}\frac{1}{T}X^{\top}X&\frac{1}{T}X^{\top}Y\\\
\frac{1}{T}Y^{\top}X&\frac{1}{T}Y^{\top}Y\end{bmatrix}.$ (4)

PCA aims to identify the most significant features that explain the majority
of the variance in $C_{XX}$ and $C_{YY}$, independently. PLS, on the other
hand, focuses on singular values and vectors that explain the covariance
component $C_{XY}$. Along the same lines, CCA aims to find linear combinations
of $X$ and $Y$ that are responsible for the correlation
($C_{XY}/\sqrt{C_{XX}C_{YY}}$) between $X$ and $Y$. See Appendix A.1 for a
detailed description of these methods.

For every numerical experiment, we generate training and test data sets
$(X_{\text{train}},Y_{\text{train}})$ and $(X_{\text{test}},Y_{\text{test}})$
according to Eqs. (2.2-2). We apply PCA, PLS, CCA, and regularized CCA (rCCA)
to the training to obtain the singular directions $W_{X_{\text{train}}}$ and
$W_{Y_{\text{train}}}$ for each method (see Appendix A.1). We then evaluate
the projections of the test data on these singular directions

$\displaystyle Z_{X_{\rm test}}$
$\displaystyle=X_{\text{test}}W_{X_{\text{train}}},$ $\displaystyle Z_{Y_{\rm
test}}$ $\displaystyle=Y_{\text{test}}W_{Y_{\text{train}}}.$ (5)

Finally, we evaluate the reconstructed correlations metric
$\mathcal{RC}^{\prime}$, which measures how well these singular directions
recover the shared signals in the data, corrected by the expected positive
bias due to the sampling noise, see Appendix A.2 for details.
$\mathcal{RC}^{\prime}=0$ corresponds to no overlap between the true and the
recovered shared directions, and $\mathcal{RC}^{\prime}=1$ corresponds to
perfect recovery.

##  4 Results

We perform numerical experiments to explore the undersampled regime,
$T\lesssim N_{X},N_{Y}$. We use $T=\\{100,300,1000,3000\\}$ samples,
$N_{X}=N_{Y}=1000$. We explore the case of one shared signal only,
$m_{\text{shared}}=1$ and we mask this shared signal by a varying number of
self signals and noise. We vary the number of retained dimensions,
($|Z_{X}|,|Z_{Y}|$), and explore how many of them are needed to recover the
shared signal in the noise and the self signal background with different SNR.

For brevity, we explore two cases: (1) One self-signal in $X$ and $Y$ in
addition to the shared signal ($m_{\text{self}}=1$); (2) many self-signals in
$X$ and $Y$. For both cases, we calculate the quality of reconstruction as the
function of the shared and the self SNR, $\gamma_{\text{shared}}$ and
$\gamma_{\text{self}}$. In all figures, we show $\mathcal{RC}^{\prime}$ for
severely undersampled (first row, $T=300$) and relatively well sampled (second
row, $T=3000$) regimes. We also show the value of $\mathcal{RC}_{0}$, the bias
that we removed from our reconstruction quality metric, for completeness, see
Appendix A.2 for details. Experiments at different parameter values can be
found in Appendix A.4.

Figure 1 shows that, in Case 1, when one dimension is retained in DR of $X$
and $Y$, PCA populates the compressed variable with the largest variance
signals and hence struggles to retain the shared signal when $\gamma_{\rm
self}>\gamma_{\rm shared}$, regardless of the number of samples. However, both
PLS and rCCA excel in achieving nearly perfect reconstructions. When $T\ll
N_{X}$, straightforward CCA cannot be applied (see A.1.3-A.1.4), but it too
achieves a perfect reconstruction when $T>N_{X}$.

Figure 1: Performance of PCA, PLS, CCA, rCCA, and noise in recovery of the
shared signal for $|Z_{X}|=|Z_{Y}|=1=m_{\text{self}}$. PCA struggles to detect
shared signals when they are weaker than the self signals. PLS and rCCA
demonstrate nearly perfect reconstruction. CCA displays no reconstruction in
the undersampled regime $T\ll N_{X}$, and it is nearly perfect for large $T$.
Figure 2: Same as Fig. 1, but for
$|Z_{X}|=|Z_{Y}|=2=m_{\text{self}}+m_{\text{shared}}$. Now there are enough
compressed variables for PCA to detect the shared signal. Other methods
perform similarly to Fig. 1, albeit the noise is larger.

In Fig. 2, we allow two dimensions in the reduced variables. For PCA, we
expect this to be sufficient to preserve both the self and the shared signals.
Indeed, PCA now works for all $\gamma$s and $T$, although with a slightly
reduced accuracy for large shared signals compared to Fig. 1. PLS and rCCA
continue to deliver highly accurate reconstructions. So does the CCA for
$T>N_{X}$. Spurious correlations, as measured by $\mathcal{RC}_{0}$ grow
slightly with the increasing dimensionality of $Z_{X}$, $Z_{Y}$ compared to
Fig. 1. This is expected since more projections must now be inferred from the
same amount of data.

Figure 3: Reconstruction results for $m_{\rm self}=30$, $m_{\rm shared}=1$,
and $|Z_{X}|=|Z_{Y}|=1$. PCA struggles to detect any shared signals when they
are even comparable to the self ones. PLS performance also degrades. CCA
displays its usual impotence at small $T$. Finally, rCCA demonstrates nearly
perfect reconstruction for all parameter values.

We now turn to $m_{\rm self}\gg m_{\rm shared}$. We use $m_{\text{shared}}=1$,
$m_{\text{self}}=30$ for concreteness. We expect that the performance of SDR
methods will degrade weakly, as they are designed to be less sensitive to the
masking effects of the self signals. In contrast, we expect IDR to be more
easily confused by the many strong self-signals, degrading the performance.
Indeed, Fig. 3 shows that PCA now faces challenges in detecting shared
signals, even when the self signals are weaker than in Fig. 1. Increasing $T$
improves its performance only slightly. Somewhat surprisingly, PLS performance
also degrades, with improvements at $T\gg N_{X}$. CCA again displays no
reconstruction when $T\ll N_{X}$, switching to near perfect reconstruction at
large $T$. Crucially, rCCA again shines, maintaining its strong performance,
consistently demonstrating nearly perfect reconstruction.

Since one retained dimension is not sufficient for PCA to represent the shared
signal when $\gamma_{\rm shared}\lesssim\gamma_{\rm self}$, we increase the
dimensionality of reduced variables $|Z_{X}|=|Z_{Y}|=m_{\text{self}}\gg
m_{\text{shared}}$), cf. Fig. 4. PCA now detects shared signals even when they
are weaker than the self-signals, $\gamma_{\rm shared}<\gamma_{\rm self}$, but
at a cost of the reconstruction accuracy plateauing significantly below 1. In
other words, when self and shared signals are comparable, they mix, allowing
for partial reconstruction. However, even at $T\gg N_{X}$, PCA cannot break
into the phase diagram’s lower right corner. Other methods perform similarly,
reconstructing shared signals over the same or wider ranges of sampling and
the SNR ratios than in Fig. 3. For all of them, the improvement comes at the
cost of the decreased asymptotic performance. The most distinct feature of
this regime is the dramatic effect of noise, where 30-dimensional compressed
variables can accumulate enough sampling fluctuations to recover correlations
that are supposedly nearly twice as high as the data actually has.

Figure 5 now explores a regime when the dimensionality of the compressed
variables is enough to store both the self and the shared interactions at the
same time, $|Z_{X}|=|Z_{Y}|=m_{\text{self}}+m_{\text{shared}}=31$. With just
one more dimension than Fig. 4, PCA abruptly transitions to being able to
recover shared signals for all SNRs, albeit still saturating at a far from
perfect performance at large $T$. PLS, CCA, rCCA, and noise show behavior
remain similar to Fig. 4.

Figure 4: DR performance for
$|Z_{X}|=|Z_{Y}|=m_{\text{self}}>m_{\text{shared}}$). PCA now detects shared
signals even when they are weaker than the self signals. However, the quality
of reconstruction is significantly lower than in Fig. 2. PLS detects signals
in a larger part of the phase space, but also with a significant reduction in
quality, which improves with sampling. CCA has its usual problem for $T\ll
N_{X}$, and, like PLS, it has a significantly lower reconstruction quality
than in the regime in Fig. 3. rCCA is able to detect the signal in the whole
phase space, but again with worse quality. Finally, spurious correlations are
high, though they decrease with better sampling. Figure 5: PCA, PLS, CCA,
rCCA, and noise results when 31 dimensions are kept after reduction
($|Z_{X}|=|Z_{Y}|=m_{\text{self}}+m_{\text{shared}}$). PCA now can detect more
shared signals when they are weaker than the self signals (A1), however, with
a significantly lower quality compared to figure 2, but suddenly explores the
whole phase space, still with lower accuracy than Case 1. PLS, CCA, rCCA, and
noise show similar behavior to figure 4.

Our analysis suggests that there are three relevant factors that determine the
ability of DR to reconstruct shared signals. The first is the strength of the
shared and the self signals compared to each other and to noise. For brevity,
in the following analysis, we fix $\gamma_{\text{self}}$ and define the ratio
$\tilde{\gamma}=\gamma_{\text{shared}}/\gamma_{\text{self}}$ to represent this
effect. The second factor affecting the performance is the ratio between the
number of shared and self signals, denoted by
$\tilde{m}=m_{\text{shared}}/m_{\text{self}}$. The third factor is the number
of samples per dimension of the reduced variable, denoted by
$\tilde{q}=T/|Z|$.

In Fig. 6, we illustrate how these parameters influence the performance of DR,
$\mathcal{RC}^{\prime}$. Each subplot varies $\tilde{q}$, while holding $T$
constant and changing $|Z_{X}|$. We compare the results of PCA (representing
IDR) and rCCA (representing SDR). Each curve is averaged over 10 trials, with
error bars indicating 1 standard deviation around the mean, using algorithmic
parameters as described in Appendix A.3.

We see that the relative strength of signals, as represented by
$\tilde{\gamma}$, plays a significant role in determining which method
performs better. If the shared signals are larger (bottom) both approaches
work. However, for weak shared signals (top), SDR is generally more effective.
Further, the ratio between the number of shared and self signals, $\tilde{m}$,
also plays an important role. When $\tilde{m}$ is large (left), IDR is more
likely to detect the shared signal before the self signals, and it approaches
the performance of SDR. However, when $\tilde{m}$ is small, IDR is more likely
to capture the self signals before moving on to the shared signals, degrading
performance (right). Finally, not surprisingly, the number of samples per
dimension of the compressed variables, $\tilde{q}$, is also critical to the
success. If $\tilde{q}$ is small, the signal is drowned in the sampling noise,
and adding more retained dimensions hurts the DR process. This expresses
itself as a peak for SDR performance around $|Z_{X}|=m_{\text{shared}}$. For
IDR, the peak is around $|Z_{X}|=m_{\text{self}}+m_{\text{shared}}$, thus
requiring more data to achieve performance similar to SDR.

Figure 6: Performance of PCA (IDR) and rCCA (SDR) for different values of the
relevant parameters of the model: the number of samples per dimension of the
compressed variable ($\tilde{q}$), the strength of shared signals relative to
the self ones ($\tilde{\gamma}$), and the ratio of the number of shared to
self signal components ($\tilde{m}$), while fixing the number of samples
($T=1000$) and the number of shared dimensions ($m_{\rm shared}=10$). Note
that decreasing $\tilde{q}$ (left to right) corresponds to increasing the
dimension of the latent space $|Z_{X}|$ at a fixed number of samples $T$.

We observe that the performance of rCCA (SDR) is almost independent of
changing $\tilde{m}$ or $\tilde{\gamma}$, indicating that it focuses on shared
dimensions even if the latter is masked by self signals. The algorithm
crucially depends on $\tilde{q}$, where adding more dimensions (decreasing
$\tilde{q}$) than needed hurts the reduction. This is because, for a fixed
number of samples, the reconstruction of each dimension then gets worse. In
contrast, for PCA (IDR), the performance depends on all three relevant
parameters, $\tilde{q}$, $\tilde{m}$, and $\tilde{\gamma}$. At some parameter
combinations, the performance of IDR in reconstructing shared signals
approaches SDR. However, in all cases, SDR never performs worse than IDR on
this task.

##  5 Discussion

We proposed a generative linear model, which we believe captures multiple
desired features of multimodal data with shared and non-shared signals. The
model focused only on data with two measured modalities. However, while not a
part of this study, the model can be readily extended to accommodate more than
two modalities (e. g., $X_{i}=R_{i}+U_{i}V_{i}+PQ_{i}$ for $i=1,...,n$, where
$n$ represents the number of modalities). Then, methods such as Tensor CCA,
which can handle more than two modalities (Luo et al., 2015), can be used to
get insight into DR on such data.

We analyzed different DR methods on data generated by this model in different
parameter regimes. We showed that linear SDR methods are clearly superior to
their IDR counterparts for detecting shared signals. We can thus make a strong
practical suggestion that, whenever the goal is to reconstruct a low
dimensional representation of covariation between two components of the data,
IDR methods (PCA) should always be avoided in favor of SDR. Of the examined
SDR approaches, rCCA is a clear winner in all parameter regimes and should
always be preferred. This finding explains the results of, for example, Sani
et al. (2021) and others that SDR can recover joint neuro-behavioral latent
spaces with fewer latent dimensions and using fewer samples than IDR methods.
Further, our observation that SDR is always superior to IDR in the context of
our model corroborates the theoretical findings of Martini & Nemenman (2023),
who proved a similar result in the context of discrete data and a different
SDR algorithm, namely the Symmetric Information Bottleneck Friedman et al.
(2013). Collectively, these diverse investigations, linear and nonlinear,
theoretical, computational, and empirical, provide strong evidence that
generic (not just linear) SDR methods are likely to be more efficient in
extracting covariation than their IDR analogs.

Our study also answers an open question in the literature surrounding the
effectiveness of SDR techniques. Specifically, there has been debate about
whether PLS, an SDR method, is effective at low sampling (Chin & Newsted,
1999; Hair et al., 2011; Goodhue et al., 2006; 2012). Our results show that
SDR is not necessarily effective in the undersampled regime. It works well
when the number of samples per retained dimension is high (even if the number
of samples per observed dimension is low), but only when the dimensionality of
the reduced description is matched to the actual dimensionality of the shared
signals.

Finally, our results can be used as a diagnostic test to determine the number
of shared versus self signals in data. As demonstrated in Fig. 6, total
correlations between $Z_{X}$ and $Z_{Y}$ obtained by applying PCA and rCCA
increase monotonically as the dimensionality of $Z$s increases, until this
dimensionality becomes larger than the signal dimensionality. For PCA, the
signal dimensionality is equal to the sum of the number of the shared and the
self signals, $m_{\text{shared}}+m_{\text{self}}$. For rCCA, it is only the
number of the shared signal. Thus increasing the dimensionality of the
compressed variables and tracking the performance of rCCA and PCA until they
diverge can be used to identify the number of self signals in the data,
provided that the data, indeed, has a low-dimensional latent structure. This
approach can be a valuable tool in various applications, where the
characterization of shared and self signals in complex systems can provide
insights into their structure and function.

In summary, we highlight a general principle that, when searching for a shared
signal between different modalities of data, SDR methods are preferable to IDR
methods. Additionally, the differences in performance between the two classes
of methods can tell us a lot about the underlying structure of the data.
Finally, for a limited number of samples, naive approaches, such as increasing
the number of compressed dimensions indefinitely to overcome the masking of
shared signals by self signals are infeasible. Thus, the use of SDR methods
becomes even more essential in such cases.

##  6 Limitations, and Future Work

While this work has provided useful insight, the assumptions made here may not
fully capture the complexity of real-world data. Specifically, our data is
generated by a linear model with random Gaussian features. It is unlikely that
real data have this exact structure. Therefore, there is a need for further
exploration of the advantages and limitations of linear DR methods on data
that have a low-dimensional, but nonlinear shared structure. This can be done
using more complex nonlinear generative models, such as nonlinearly
transforming the data generated by Eq. (2.2-2), or random feature two-layered
neural network models (Rocks & Mehta, 2022).

A different possible future research direction is to explore the performance
of nonlinear DR methods on data from generative models with a latent low-
dimensional nonlinear structure. Autoencoders and their variational extensions
are a natural extension of IDR to learn nonlinear reduced dimensional
representations (Hinton & Salakhutdinov, 2006; Kingma & Welling, 2014; Higgins
et al., 2016). Meanwhile, Deep CCA and its variational extensions (Andrew et
al., 2013; Wang et al., 2015; Chandar et al., 2016; Wang et al., ) should be
explored as a nonlinear version of SDR. Both of these types of methods can
potentially capture more complex relationships between the modalities and
improve the quality of the reduced representations, and it is not clear if the
SDR class of methods is always more efficient than the IDR one.

Further, our analysis depends on the choice of metric used to quantify the
performance of DR, and different choices should also be explored. For example,
to capture nonlinear correlations, mutual information can be utilized to
quantify the relationships between the reduced representations.

Despite the aforementioned limitations, we believe that our work provides a
compelling addition to the body of knowledge that SDR outperforms IDR in
detecting shared signals quite generally.

#### Acknowledgments

We thank Sean Ridout, Michael Pasek, and Arabind Swain for providing feedback
on the manuscript. IN thanks Philipp Fleig for useful discussions. This work
was funded, in part, by NSF Grant No. 2010524 and 2014173 and by the Simons
Investigator award to IN.

## References

  * Andrew et al. (2013) Galen Andrew, Raman Arora, Jeff Bilmes, and Karen Livescu.  Deep canonical correlation analysis.  In Sanjoy Dasgupta and David McAllester (eds.), _Proceedings of the 30th International Conference on Machine Learning_ , volume 28 of _Proceedings of Machine Learning Research_ , pp. 1247–1255, Atlanta, Georgia, USA, 17–19 Jun 2013. PMLR.
  * Baillie et al. (2002) Richard T Baillie, G Geoffrey Booth, Yiuman Tse, and Tatyana Zabotina.  Price discovery and common factor models.  _Journal of financial markets_ , 5(3):309–321, 2002.
  * Bialek et al. (2001) William Bialek, Ilya Nemenman, and Naftali Tishby.  Predictability, complexity, and learning.  _Neural computation_ , 13(11):2409–2463, 2001.
  * Bouchaud et al. (2007) Jean Philippe Bouchaud, Laurent Laloux, M Augusta Miceli, and Marc Potters.  Large dimension forecasting models and random singular value spectra.  _The European Physical Journal B_ , 55:201–207, 2007.
  * Bun et al. (2017) Joël Bun, Jean-Philippe Bouchaud, and Marc Potters.  Cleaning large correlation matrices: Tools from random matrix theory.  _Physics Reports_ , 666:1–109, 2017.  ISSN 0370-1573.  doi: https://doi.org/10.1016/j.physrep.2016.10.005.  Cleaning large correlation matrices: tools from random matrix theory.
  * Chandar et al. (2016) Sarath Chandar, Mitesh M Khapra, Hugo Larochelle, and Balaraman Ravindran.  Correlational neural networks.  _Neural Computation_ , 28(2):257–285, 2016.  doi: 10.1162/NECO˙a˙00801.
  * Chapman & Wang (2021) James Chapman and Hao-Ting Wang.  Cca-zoo: A collection of regularized, deep learning based, kernel, and probabilistic cca methods in a scikit-learn style framework.  _Journal of Open Source Software_ , 6(68):3823, 2021.
  * Chen et al. (2022) Boyuan Chen, Kuang Huang, Sunand Raghupathi, Ishaan Chandratreya, Qiang Du, and Hod Lipson.  Automated discovery of fundamental variables hidden in experimental data.  _Nature Computational Science_ , 2(7):433–442, 2022.
  * Chin & Newsted (1999) Wynne Chin and P. Newsted.  Structural equation modeling analysis with small samples using partial least square.  _Statistical Strategies for Small Sample Research_ , 01 1999.
  * Clark et al. (2013) Kenneth Clark, Bruce Vendt, Kirk Smith, John Freymann, Justin Kirby, Paul Koppel, Stephen Moore, Stanley Phillips, David Maffitt, Michael Pringle, et al.  The cancer imaging archive (tcia): maintaining and operating a public information repository.  _Journal of digital imaging_ , 26:1045–1057, 2013.
  * Colwell et al. (2014) Lucy J Colwell, Yu Qin, Miriam Huntley, Alexander Manta, and Michael P Brenner.  Feynman-hellmann theorem and signal identification from sample covariance matrices.  _Physical Review X_ , 4(3):031032, 2014.
  * Creutzig et al. (2009) Felix Creutzig, Amir Globerson, and Naftali Tishby.  Past-future information bottleneck in dynamical systems.  _Physical Review E_ , 79(4):041925, 2009.
  * Fleig & Nemenman (2022) Philipp Fleig and Ilya Nemenman.  Statistical properties of large data sets with linear latent features.  _Physical Review E_ , 106(1):014102, 2022.
  * Forni & Gambetti (2010) Mario Forni and Luca Gambetti.  The dynamic effects of monetary policy: A structural factor model approach.  _Journal of Monetary Economics_ , 57(2):203–216, 2010.
  * Freyaldenhoven (2022) Simon Freyaldenhoven.  Factor models with local factors—determining the number of relevant factors.  _Journal of Econometrics_ , 229(1):80–102, 2022.
  * Friedman et al. (2013) Nir Friedman, Ori Mosenzon, Noam Slonim, and Naftali Tishby.  Multivariate information bottleneck.  _arXiv preprint arXiv:1301.2270_ , 2013.
  * Goodhue et al. (2006) Dale L Goodhue, William Lewis, and Ron Thompson.  Pls, small sample size, and statistical power in mis research.  In _Proceedings of the 39th Annual Hawaii International Conference on System Sciences (HICSS’06)_ , volume 8, pp. 202b–202b, 2006.  doi: 10.1109/HICSS.2006.381.
  * Goodhue et al. (2012) Dale L Goodhue, William Lewis, and Ron Thompson.  Does pls have advantages for small sample size or non-normal data?  _MIS Quarterly_ , 36(3):981–1001, 2012.  ISSN 02767783.
  * Goodhue et al. (2013) Dale L Goodhue, Ron Thompson, and William Lewis.  Why you shouldn’t use pls: Four reasons to be uneasy about using pls in analyzing path models.  In _2013 46th Hawaii International Conference on System Sciences_ , pp. 4739–4748. IEEE, 2013.
  * Gosselin & Tkacz (2001) Marc-André Gosselin and Greg Tkacz.  Evaluating factor models: An application to forecasting inflation in canada.  Technical report, Bank of Canada, 2001.
  * Hair et al. (2011) Joe F Hair, Christian M Ringle, and Marko Sarstedt.  Pls-sem: Indeed a silver bullet.  _Journal of Marketing Theory and Practice_ , 19(2):139–152, 2011.  doi: 10.2753/MTP1069-6679190202.
  * Hair Jr et al. (2021) Joe Hair Jr, Joseph F Hair Jr, G Tomas M Hult, Christian M Ringle, and Marko Sarstedt.  _A primer on partial least squares structural equation modeling (PLS-SEM)_.  Sage publications, 2021.
  * Higgins et al. (2016) Irina Higgins, Loic Matthey, Arka Pal, Christopher Burgess, Xavier Glorot, Matthew Botvinick, Shakir Mohamed, and Alexander Lerchner.  beta-vae: Learning basic visual concepts with a constrained variational framework.  In _International conference on learning representations_ , 2016.
  * Hinton & Salakhutdinov (2006) Geoffrey E Hinton and Ruslan Salakhutdinov.  Reducing the dimensionality of data with neural networks.  _Science_ , 313(5786):504–507, 2006.  doi: 10.1126/science.1127647.
  * Hotelling (1933) Harold Hotelling.  Analysis of a complex of statistical variables into principal components.  _Journal of educational psychology_ , 24(6):417, 1933.
  * Hotelling (1936) Harold Hotelling.  Relations between two sets of variates.  _Biometrika_ , 1936.  doi: 10.1007/978-1-4612-4380-9˙14.
  * Huntley et al. (2015) Rachael P Huntley, Tony Sawford, Prudence Mutowo-Meullenet, Aleksandra Shypitsyna, Carlos Bonilla, Maria J Martin, and Claire O’Donovan.  The goa database: gene ontology annotation updates for 2015.  _Nucleic acids research_ , 43(D1):D1057–D1063, 2015.
  * Kingma & Welling (2014) Diederik P Kingma and Max Welling.  Auto-Encoding Variational Bayes.  In _2nd International Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014, Conference Track Proceedings_ , 2014.
  * Kock & Hadaya (2018) Ned Kock and Pierre Hadaya.  Minimum sample size estimation in pls-sem: The inverse square root and gamma-exponential methods.  _Information systems journal_ , 28(1):227–261, 2018.
  * Krakauer et al. (2017) John W Krakauer, Asif A Ghazanfar, Alex Gomez-Marin, Malcolm A MacIver, and David Poeppel.  Neuroscience needs behavior: correcting a reductionist bias.  _Neuron_ , 93(3):480–490, 2017.
  * Lorenzi et al. (2018) Marco Lorenzi, Andre Altmann, Boris Gutman, Selina Wray, Charles Arber, Derrek P Hibar, Neda Jahanshad, Jonathan M Schott, Daniel C Alexander, Paul M Thompson, et al.  Susceptibility of brain atrophy to trib3 in alzheimer’s disease, evidence from functional prioritization in imaging genetics.  _Proceedings of the National Academy of Sciences_ , 115(12):3162–3167, 2018.
  * Luo et al. (2015) Yong Luo, Dacheng Tao, Kotagiri Ramamohanarao, Chao Xu, and Yonggang Wen.  Tensor canonical correlation analysis for multi-view dimension reduction.  _IEEE transactions on Knowledge and Data Engineering_ , 27(11):3111–3124, 2015.
  * Marchenko & Pastur (1967) Vladimir Alexandrovich Marchenko and Leonid Andreevich Pastur.  Distribution of eigenvalues for some sets of random matrices.  _Matematicheskii Sbornik_ , 114(4):507–536, 1967.
  * Martini & Nemenman (2023) K Michael Martini and Ilya Nemenman.  Data efficiency, dimensionality reduction, and the generalized symmetric information bottleneck.  _arXiv preprint arXiv:2309.05649_ , 2023.
  * Natraj et al. (2022) Nikhilesh Natraj, Daniel B Silversmith, Edward F Chang, and Karunesh Ganguly.  Compartmentalized dynamics within a common multi-area mesoscale manifold represent a repertoire of human hand movements.  _Neuron_ , 110(1):154–174.e12, 2022.  ISSN 0896-6273.  doi: https://doi.org/10.1016/j.neuron.2021.10.002.
  * Pacharawongsakda & Theeramunkong (2016) Eakasit Pacharawongsakda and Thanaruk Theeramunkong.  A comparative study on single and dual space reduction in multi-label classification.  In Andrzej M.J. Skulimowski and Janusz Kacprzyk (eds.), _Knowledge, Information and Creativity Support Systems: Recent Trends, Advances and Solutions_ , pp. 389–400, Cham, 2016. Springer International Publishing.  ISBN 978-3-319-19090-7.
  * Pang et al. (2016) Rich Pang, Benjamin J Lansdell, and Adrienne L Fairhall.  Dimensionality reduction in neuroscience.  _Current Biology_ , 26(14):R656–R660, 2016.
  * Pedregosa et al. (2011) Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al.  Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830, 2011.
  * Rocks & Mehta (2022) Jason W Rocks and Pankaj Mehta.  Memorizing without overfitting: Bias, variance, and interpolation in overparameterized models.  _Physical review research_ , 4(1):013201, 2022.
  * Rudd (2020) Jeremy B Rudd.  Underlying inflation: Its measurement and significance.  2020\.
  * Sani et al. (2021) Omid G Sani, Hamidreza Abbaspourazad, Yan T Wong, Bijan Pesaran, and Maryam M Shanechi.  Modeling behaviorally relevant neural dynamics enabled by preferential subspace identification.  _Nature Neuroscience_ , 24(1):140–149, 2021.
  * Sponberg et al. (2015) Simon Sponberg, Thomas L Daniel, and Adrienne L Fairhall.  Dual dimensionality reduction reveals independent encoding of motor features in a muscle synergy for insect flight control.  _PLOS Computational Biology_ , 11(4):1–23, 04 2015.  doi: 10.1371/journal.pcbi.1004168.
  * Steinmetz et al. (2021) Nicholas A Steinmetz, Cagatay Aydin, Anna Lebedeva, Michael Okun, Marius Pachitariu, Marius Bauza, Maxime Beau, Jai Bhagat, Claudia Böhm, Martijn Broux, et al.  Neuropixels 2.0: A miniaturized high-density probe for stable, long-term brain recordings.  _Science_ , 372(6539):eabf4588, 2021.
  * Stringer et al. (2019) Carsen Stringer, Marius Pachitariu, Nicholas Steinmetz, Charu Bai Reddy, Matteo Carandini, and Kenneth D. Harris.  Spontaneous behaviors drive multidimensional, brainwide activity.  _Science_ , 364(6437):eaav7893, 2019.  doi: 10.1126/science.aav7893.
  * Svensson et al. (2018) Valentine Svensson, Roser Vento-Tormo, and Sarah A Teichmann.  Exponential scaling of single-cell rna-seq in the past decade.  _Nature protocols_ , 13(4):599–604, 2018.
  * Urai et al. (2022) Anne E Urai, Brent Doiron, Andrew M Leifer, and Anne K Churchland.  Large-scale neural recordings call for new insights to link brain and behavior.  _Nature neuroscience_ , 25(1):11–19, 2022.
  * Vinod (1976) Hrishikesh D Vinod.  Canonical ridge and econometrics of joint production.  _Journal of Econometrics_ , 1976.  doi: 10.1016/0304-4076(76)90010-5.
  * (48) Weiran Wang, Xinchen Yan2 Honglak Lee, and Karen Livescu.  Deep variational canonical correlation analysis.
  * Wang et al. (2015) Weiran Wang, Raman Arora, Karen Livescu, and Jeff Bilmes.  On deep multi-view representation learning.  In _International conference on machine learning_ , pp. 1083–1092. PMLR, 2015.
  * Wold et al. (2001) Svante Wold, Michael Sjöström, and Lennart Eriksson.  Pls-regression: a basic tool of chemometrics.  _Chemometrics and Intelligent Laboratory Systems_ , 58(2):109–130, 2001.  ISSN 0169-7439.  doi: https://doi.org/10.1016/S0169-7439(01)00155-1.  PLS Methods.
  * Zheng et al. (2017) Grace XY Zheng, Jessica M Terry, Phillip Belgrader, Paul Ryvkin, Zachary W Bent, Ryan Wilson, Solongo B Ziraldo, Tobias D Wheeler, Geoff P McDermott, Junjie Zhu, et al.  Massively parallel digital transcriptional profiling of single cells.  _Nature communications_ , 8(1):14049, 2017.
  * Årup Nielsen et al. (1998) Finn Årup Nielsen, Lars Kai Hansen, and Stephen C Strother.  Canonical ridge analysis with ridge parameter optimization.  _NeuroImage_ , 1998.  doi: 10.1016/s1053-8119(18)31591-x.

##  Appendix A Appendix

###  A.1 Linear Dimensionality Reduction Methods

####  A.1.1 PCA

PCA is a widely used linear IDR method that aims to find the orthogonal
principal directions, such that a few of them explain the largest possible
fraction of the variance within the data. PCA decomposes the covariance matrix
of the data matrix $X$, $C_{XX}=\frac{1}{T}X^{\top}X$, into its eigenvectors
and eigenvalues through singular value decomposition (SVD). The SVD yields
orthogonal directions, represented by the vectors $w^{(i)}_{X}$, that capture
the most significant variability in the data. In most numerical
implementations (Pedregosa et al., 2011), these directions are obtained
consecutively, one by one, such that the dot product between any two
directions is zero $w_{X}^{(i)}\cdot w_{X}^{(j)}=\delta_{ij}$. The
eigenvectors $w_{X}^{(i)}$ are obtained as the best solution to the
optimization problem:

$w^{*(i)}_{X}=\underset{w^{(i)}_{X}}{\mathrm{arg\max}}\frac{{w^{(i)}_{X}}^{\top}{X^{(i)}}^{\top}X^{(i)}w^{(i)}_{X}}{{w^{(i)}_{X}}^{\top}w^{(i)}_{X}}.$
(6)

Here $X^{(i)}$ is the $i$th deflated matrix where $X^{(1)}$ is the original
matrix, and for every subsequent $i+1$, the matrix is deflated by subtracting
the projection of $X$ on the obtained weights:
$X^{(i+1)}=X-\Sigma_{s=1}^{i}Xw_{(s)}w_{(s)}^{\top}$. The eigenvectors are
sorted in decreasing order according to their corresponding eigenvalues, and
the first $k$ eigenvectors $w_{X}^{(i=1:k)}$ are selected to form the
projection matrix $W_{X}$. The obtained vectors determine the size of the
reduced form $Z_{X}$, where $|Z_{X}|=k$ is the number of vectors retained from
the decomposition of $X$. The vectors $w_{X}^{(i)}$ are then stacked together
to form the projection matrix $W_{X}$. The low-dimensional representation
$Z_{X}$ is then obtained by multiplying the original data matrix $X$ with this
projection matrix, resulting in the reduced data matrix $Z_{X}=XW_{X}$.
Similar treatment is done for $Y$ in order to obtain $Z_{Y}=YW_{Y}$

One of the main advantages of PCA is its simplicity and efficiency. However,
one of the drawbacks of this method is that it performs DR for $X$ and $Y$
independently, and one then searches for relations between $Z_{X}$ and $Z_{Y}$
by regressing one on the other. Thus obtained low-dimensional descriptions may
capture variance but not the covariance between the two datasets.

####  A.1.2 PLS

PLS, or Partial Least Squares, performs SDR by finding the shared signals that
explain the maximum covariance between two sets of data (Wold et al., 2001).
PLS performs the SVD of the covariance matrix $C_{XY}=\frac{1}{T}X^{\top}Y$
(or equivalently $C_{YX}=\frac{1}{T}Y^{\top}X$). The left and right singular
vectors ($w^{*(i)}_{X},w^{*(i)}_{Y}$) are obtained consecutively pair by pair
such that $w_{X}^{(i)}\cdot w_{Y}^{(j)}=\delta_{ij}$. They are solutions of
the optimization problem:

$(w^{*(i)}_{X},w^{*(i)}_{Y})=\underset{w^{(i)}_{X},w^{(i)}_{Y}}{\mathrm{arg\max}}\frac{{w^{(i)}_{X}}^{\top}{X^{(i)}}^{\top}Y^{(i)}w^{(i)}_{Y}}{\sqrt{({w^{(i)}_{X}}^{\top}w^{(i)}_{X})({w^{(i)}_{Y}}^{\top}w^{(i)}_{Y})}}$
(7)

The matrices ${X^{(i)}},Y^{(i)}$ are deflated in a similar manner to PCA
A.1.1. The singular vectors are sorted in the decreasing order of their
corresponding singular values, and the first $k$ vectors are selected to form
the projection matrices $(W_{X},W_{Y})$. The obtained vectors determine the
size of the reduced form $(Z_{X},Z_{Y})$, where $|Z_{X}|=|Z_{Y}|=k$ is the
number of vectors retained. The vectors $(w_{X}^{(i)},w_{Y}^{(i)})$ are then
stacked together to form the projection matrices $(W_{X},W_{Y})$ respectively.
The low-dimensional representations $(Z_{X},Z_{Y})$ are obtained by projecting
the original data matrices $(X,Y)$ onto these projection matrices:
$Z_{X}=XW_{X}$, and $Z_{Y}=YW_{Y}$.

In summary, PLS performs simultaneous reduction on both datasets, maximizing
the covariance between the reduced representations $Z_{X}$ and $Z_{Y}$. This
property makes PLS a powerful tool for studying the relationships between two
datasets and identifying the underlying factors that explain their joint
variability.

####  A.1.3 CCA

CCA is another SDR method, which aims to find the directions that explain the
maximum correlation between two datasets Hotelling (1936). However, unlike
PLS, CCA obtains the shared signals by performing SVD on the correlation
matrix $\frac{C_{XY}}{\sqrt{C_{XX}}\sqrt{C_{YY}}}$. The singular vectors
($w^{*(i)}_{X},w^{*(i)}_{Y}$) are obtained consecutively pair by pair such
that $w_{X}^{(i)}\cdot w_{Y}^{(j)}=\delta_{ij}$. CCA enforces the
orthogonality of $w_{X}^{(i)},w_{Y}^{(i)}$ independently as well, such that
$w_{X}^{(i)}\cdot w_{X}^{(j)}=w_{Y}^{(i)}\cdot w_{Y}^{(j)}=\delta_{ij}$. The
singular vectors are obtained by solving the optimization problem:

$(w^{*(i)}_{X},w^{*(i)}_{Y})=\underset{w^{(i)}_{X},w^{(i)}_{Y}}{\mathrm{arg\max}}\frac{{w^{(i)}_{X}}^{\top}{X^{(i)}}^{\top}Y^{(i)}w^{(i)}_{Y}}{\sqrt{({w^{(i)}_{X}}^{\top}{X^{(i)}}^{\top}X^{(i)}w^{(i)}_{X})({w^{(i)}_{Y}}^{\top}{Y^{(i)}}^{\top}Y^{(i)}w^{(i)}_{Y})}}.$
(8)

Like in PLS A.1.2, the matrices ${X^{(i)}},Y^{(i)}$ are deflated in a similar
manner. In addition, the first $k$ singular vectors
$(w^{*(i)}_{X},w^{*(i)}_{Y})$ are stacked together to form the projection
matrices $(W_{X},W_{Y})$, which then are used to obtain the reduced data
matrices $Z_{X}=XW_{X}$, and $Z_{Y}=YW_{Y}$.

One of the key differences between PLS and CCA is that while both perform SDR,
CCA also simultaneously performs IDR implicitly. Indeed, it involves
multiplication of $C_{XY}$ by $C_{XX}^{-1/2}$ on the left and $C_{YY}^{-1/2}$
on the right, which, in turn, requires finding singular values of the $X$ and
the $Y$ data matrices independently.

####  A.1.4 Regularized CCA - rCCA

While CCA is a useful method for finding the maximum correlating features
between two sets of data, it does have some limitations. Specfically, in the
undersampled regime, where $T\leq\max(N_{X},N_{Y})$, the matrices $C_{XX}$ and
$C_{YY}$ are singular and their inverses do not exist. Using the pseudoinverse
to solve the problem can lead to numerical instability and sensitivity to
noise. Regularized CCA (rCCA) (Vinod, 1976; Årup Nielsen et al., 1998)
overcomes this problem by adding a small regularization term to the covariance
matrices, allowing them to be invertible. Specifically, one tales

$\displaystyle\tilde{C}_{XX}$ $\displaystyle=C_{XX}+c_{X}I_{X},$ (9)
$\displaystyle\tilde{C}_{YY}$ $\displaystyle=C_{YY}+c_{Y}I_{Y},$ (10)

where $\tilde{C}_{XX},\tilde{C}_{YY}$ are the new regularized matrices,
$c_{X},c_{Y}>0$ are small regularization parameters and $I_{X},I_{Y}$ are
identity matrices with sizes $N_{X}\times N_{X},N_{Y}\times N_{Y}$
respectively.

This original implementation of rCCA resulted in correlation matrices with
diagonals not equal to one. Thus, a better implementation uses a different
form of regularization (Årup Nielsen et al., 1998) by adding the
regularization parameters $c_{X}$ and $c_{Y}$ individually to the equations as
an affine combination (i. e., $\sum_{i}^{n}c_{i}=1$) as the following:

$\displaystyle\tilde{C}_{XX}$
$\displaystyle=\frac{1}{T}(c_{X_{1}}w_{X}^{\top}X^{\top}Xw_{X}+c_{X_{2}}w_{X}^{\top}w_{X})$
(11) $\displaystyle\tilde{C}_{YY}$
$\displaystyle=\frac{1}{T}(c_{Y_{1}}w_{Y}^{\top}Y^{\top}Yw_{Y}+c_{Y_{2}}w_{Y}^{\top}w_{Y}).$
(12)

This results in the regularized equations for $X$ and $Y$ to be:

$\displaystyle\tilde{C}_{XX}$
$\displaystyle=\frac{1}{T}\big{(}(1-c_{X})w_{X}^{\top}X^{\top}Xw_{X}+c_{X}w_{X}^{\top}w_{X}\big{)}$
(13) $\displaystyle\tilde{C}_{YY}$
$\displaystyle=\frac{1}{T}\big{(}(1-c_{Y})w_{Y}^{\top}Y^{\top}Yw_{Y}+c_{Y}w_{Y}^{\top}w_{Y}\big{)},$
(14)

where $c_{X}$ and $c_{Y}$ are the regularization parameters, with values
between 0 and 1. Writing the regularization conditions in this form is in fact
a convex interpolation problem between PLS and CCA, which is a more robust
solution and does not suffer from shortening the length of correlations due to
the added regularization. As a result, this implementation of rCCA achieves
the best accuracy among all other methods.

###  A.2 Assessing Success and Sampling Noise Treatment

To assess the success of DR, we calculated the ratio between the total
correlation between $Z_{X_{\rm test}}$ and $Z_{Y_{\rm test}}$, defined as in
Eq. (3), and the total correlation between $X$ and $Y$, which we input into
the model. Specifically, we take the total correlation as the Frobenius norm
of the correlation matrix, $||{A}||_{F}=\sqrt{\sum_{i}{\sigma_{i}^{2}(A)}}$,
where $\sigma(A)$ are the singular values of the matrix $A$. Therefore, the
metric of the quality of the DR is

$\mathcal{RC}=\frac{||\text{Corr}(Z_{X_{\rm test}},Z_{Y_{\rm
test}})||_{F}}{||\text{Corr}(P,P)||_{F}}=\frac{||\text{Corr}(Z_{X_{\rm
test}},Z_{Y_{\rm test}})||_{F}}{m_{\rm shared}},$ (15)

where Corr stands for the correlation matrix between its arguments, and we use
$||\text{Corr}(P,P)||_{F}=m_{\text{shared}}$ as the total shared correlation
that one needs to recover. Statistical fluctuations aside, $\mathcal{RC}$
should vary between zero (bad reconstruction of the shared variables) and one
(perfect reconstruction).

Figure 7: The resulting correlations are averages of all the points in the
phasespace, then averaged over 10 different realizations of the matrices. The
error bars are for two standard deviations around the mean

In many real-world applications, the number of available samples, $T$, is
often limited compared to the dimensionality of the data, $N_{X}$ and $N_{Y}$.
This undersampling can introduce spurious correlations. We are not aware of
analytical results to calculate the effects of the sampling noise on
estimating singular values in the model in Eq. (2.2) (Bun et al., 2017). Thus,
to estimate the effect of the sampling noise, we adopt an empirical approach.
Specifically, we generate two random matrices, $Z_{X_{\text{random}}}$ and
$Z_{Y_{\text{random}}}$, of sizes $T\times|Z_{X}|$ and $T\times|Z_{Y}|$,
respectively. We then calculate the correlation between these matrices,
denoted as $\mathcal{RC}_{0}$, for multiple such trials using the metric in
Eq. (15). For random $Z_{X_{\text{random}}}$ and $Z_{Y_{\text{random}}}$,
$\mathcal{RC}$ should be zero. However, Fig. 7 shows that, especially for
large dimensionalities of the compressed variables and small $T$, the sampling
noise results in a significant spurious $\mathcal{RC}_{0}>0$, which may even
be larger than 1! Crucially, $\mathcal{RC}_{0}$ does not fluctuate around its
mean across trials, so that the sampling bias is narrowly distributed.

To compensate for this sampling bias, we subtract it from the reconstruction
quality metric,

$\mathcal{RC}^{\prime}=\mathcal{RC}-\mathcal{RC}_{0}.$ (16)

It is this $\mathcal{RC}^{\prime}$ that we plot in all Figures in this paper
as the ultimate metric of the reconstruction quality. While subtracting the
bias is not the most rigorous mathematically, it provides a practical approach
for reducing the effects of the sampling noise.

###  A.3 Implementation

We used Python and the scikit-learn (Pedregosa et al., 2011) library for
performing PCA, PLS, and CCA, while the cca-zoo (Chapman & Wang, 2021) library
was used for rCCA. For PCA, SVD was performed with default parameters. For
PLS, the PLS Canonical method was used with the NIPALS algorithm. For both PLS
and CCA, the tolerance was set to $10^{-4}$ with a maximum convergence limit
of 5000 iterations. For rCCA, regularization parameters were set as
$c_{1}=c_{2}=0.1$. All other parameters not explicitly here were set to their
default values.

All figures shown in this paper were averaged over 10 independent realizations
of $R_{X},R_{Y},U_{X},U_{Y},P$, while fixing the projection matrices
$V_{X},V_{Y},Q_{X},Q_{Y}$. We then performed an additional round of averaging
everything over 10 realizations of the projection matrices themselves. The
simulations were parallelized and run on Amazon Web Services (AWS) servers of
instance types ml.c5.2xlarge.

###  A.4 Extended Figures

Here we provide results of simulations similar to the main text, but with a
wider range of $T$.

Figure 8: Performance of PCA, PLS, CCA, and rCCA in detecting shared signals
with one self signal and one dimension kept after DR. Figure 9: Performance of
PCA, PLS, CCA, and rCCA in detecting shared signals with one self signal and
two dimensions kept after DR. Figure 10: Performance of PCA, PLS, CCA, and
rCCA in detecting shared signals with 30 self signals and one dimension kept
after DR. Figure 11: Performance of PCA, PLS, CCA, and rCCA in detecting
shared signals among 30 self signals and with 30 dimensions kept after DR.
Figure 12: Performance of PCA, PLS, CCA, and rCCA in detecting shared signals
among 30 self signals and with 31 dimensions kept after DR.