# Constraining pseudo-Dirac neutrinos from a galactic core-collapse supernova

Manibrata Sen  Max-Planck-Institut für Kernphysik, Saupfercheckweg 1, 69117
Heidelberg, Germany

###### Abstract

Neutrinos can be pseudo-Dirac in nature – Majorana fermions behaving as Dirac
fermions for all practical purposes. In such a scenario, active and sterile
neutrinos are quasi-degenerate in mass, and hence oscillations between the
two, due to their tiny mass-squared difference $(\delta m^{2})$, can develop
only over very long baselines. Under this hypothesis, we analyze the neutrino
data from SN1987A, and find a mild preference for a non-zero mass-squared
difference. The same data can also be used to exclude values of $\delta
m^{2}\sim 10^{-20}{\rm eV}^{2}$ – the smallest constrained so far. We also
discuss how next-generation experiments like the DUNE and Hyper-Kamiokande can
probe this scenario for a future galactic supernova.

##  1 Introduction

The origin of neutrino masses has been eluding the particle physics community
for decades. Broadly, neutrino masses can either be categorized into Dirac, or
Majorana, according whether lepton number is a conserved symmetry of the
Standard Model (SM) or not. In fact, probing lepton number violation is one of
the very few ways to distinguish between these two hypothesis, since in the
relativistic limit, Dirac and Majorana neutrinos behave identically in
oscillation experiments. However, it is also possible that while lepton number
is violated in the SM, it happens in a controlled, soft manner, such that for
all practical purposes, neutrinos behave as Dirac fermions, although they are
originally Majorana fermions. These pseudo-Dirac (PD) fermions are
characterized by a maximal mixture of almost degenerate active and sterile
states. Such a tiny mass-squared difference $\delta m^{2}$ allows for the
possibility of oscillations between active-sterile neutrinos over very long
baselines.

Since these active-sterile oscillations take place over a long distance
$\propto\left(E_{\nu}/\delta m^{2}\right)$, the strongest bounds come from
astrophyical sources like solar neutrinos [1]($\delta m^{2}\lesssim
10^{-12}{\rm eV}^{2}$), atmospheric neutrinos [2]($\delta m^{2}\lesssim
10^{-4}{\rm eV}^{2}$), and even smaller values ($\delta m^{2}\sim 10^{-24}{\rm
eV}^{2}$) can be probed by a measurement of the diffuse supernova background
neutrinos [3]. Of particular interest to this scenario is a galactic core-
collapse supernova (SN), which emits almost all of its energy in the form of
${\cal O}({\rm MeV})$ neutrinos. The relatively small energy of the neutrinos,
coupled with the baseline of around ${\cal O}({\rm kpc})$, makes it a natural
laboratory to test oscillations driven by tiny $\delta m^{2}$.

In this proceeding based on [4], we will discuss how one can use neutrino
observation data from SN1987A by the Kamiokande-II (KII) collaboration [5],
the IMB collaboration [6] and the Baksan collaboration [7] to probe active-
sterile oscillations arising in the PD hypothesis. Using an unbinned
likelihood analysis, we find that the combined data from all the three
experiments have a slight preference for active-sterile oscillations, over the
no-oscillation hypothesis. On the other hand, one can use the same
observations to constrain a mass-squared difference, $\delta m^{2}\sim
10^{-20}{\rm eV}^{2}$, the smallest value probed so far. Finally, we also
discuss how future experiments like the Hyper-Kamiokande (HK) and the Deep
Underground Neutrino Experiment (DUNE) can utilize a galactic SN to constrain
even tinier values of $\delta m^{2}$, thereby providing stringent bounds on
the PD hypothesis.

##  2 Pseudo-Dirac Neutrinos

The SM can be minimally extended by adding three right handed neutrinos $(N)$.
After electroweak symmetry breaking, the neutrino mass matrix $(M_{\nu})$
consists of a Dirac mass term, given by $Yv/\sqrt{2}$ (where $v/\sqrt{2}$
denotes the Higgs vacuum expectation value, and $Y$ refers to the Yukawa
matrix), and a Majorana mass $(M_{N})$ term. The limit of soft lepton-number
violation is achieved by setting $M_{N}\ll Yv$, such that $M_{N}$ breaks the
degeneracy between the masses of the left and right handed states. In such a
limit, $M_{\nu}$ can be diagonalized by a $6\times 6$ unitary matrix,
consisting of two diagonal unitary matrices: the PMNS matrix $U$, and another
unitary matrix $U_{N}$ that diagonalizes the sterile sector [8]. A neutrino
flavor field $\nu_{\beta L}$ ($\beta=e,\mu,\tau$) can be written as a
superposition of an active and a sterile field, such that $\nu_{\beta
L}=U_{\beta k}(\nu_{ks}+i\,\nu_{ka})/\sqrt{2}$ for $k=1,2,3$. Note that due to
the smallness of $M_{N}$, the states are maximally mixed, and have quasi-
degenerate masses $m_{ks,ka}^{2}=m_{k}^{2}\pm\delta m_{k}^{2}/2$. Henceforth,
for simplicity, we will drop the subscript $k$, and assume that the mass
difference is the same for all states, given by $\delta m^{2}$.

If the $\delta m^{2}$ is much smaller than the mass-squared differences
associated with the solar and atmospheric differences, oscillations due to the
latter are averaged out over astrophysical distances, whereas those due to
$\delta m^{2}$ survive. The oscillation probability between $\nu_{\beta}$ and
$\nu_{\gamma}$, driven by $\delta m^{2}$, can be computed as,

$\displaystyle P_{\beta\gamma}=\frac{1}{2}\left(1+e^{-\left(\frac{L}{L_{\rm
coh}}\right)^{2}}\cos\left(\frac{2\pi L}{L_{\rm
osc}}\right)\right)\times\sum_{k}\left|U_{\beta k}\right|^{2}\left|U_{\gamma
k}\right|^{2}\,,$ (1)

where $L_{\rm osc}=(4\pi E_{\nu})/(\delta m^{2})$ indicates the oscillation
length for a neutrino of energy $E_{\nu}$. Neutrinos traversing astrophysical
baselines can be sensitive to decoherence due to separation of wavepackets
after a certain distance, given by $L_{\rm
coh}=(4\sqrt{2}E_{\nu}^{2}\sigma_{x})/(|\delta m^{2}|)$, where $\sigma_{x}$
measures the size of the neutrino wave-packet.

##  3 Application to neutrino data from SN1987A

The neutrino fluence from a SN at a distance $d$ can be approximated by the
alpha-fit spectra [11],

$\displaystyle\phi_{\beta}(E_{\nu})=\frac{{\cal E}_{\rm
tot}}{E_{0\beta}}\frac{1}{4\pi
d^{2}}\frac{(1+\alpha)^{1+\alpha}}{\Gamma(1+\alpha)}\left(\frac{E_{\nu}}{E_{0\beta}}\right)^{\alpha}e^{-(1+\alpha)\frac{E_{\nu}}{E_{0\beta}}},$
(2)

where $E_{0\beta}$ is the average energy for $\nu_{\beta}$, and ${\cal E}_{\rm
tot}$ gives the total energy emitted in ergs. Here, $\alpha=2.3$ determines
the width of the distributions, and we fix $d=50\,$kpc.

To test our hypothesis, we consider that the fluence in Eq. 2 gets processed
by the PD probability in Eq. 1, and fit it to the $\mathcal{O}(30)$ events
observed in KII, IMB and Baksan in total. We perform an unbinned extended
likelihood $(\mathscr{L})$ analysis for each experiment, by carefully adopting
the background treatment [9]. We take effectively two flavors $e$, and $x$,
which can be a linear combination of $\mu$ and $\tau$. We fit the parameters
$\\{{\cal E}_{\rm tot},E_{0e},E_{0x},\delta m^{2}\\}$, while
$\sigma_{x}=10^{-13}$ m is kept fixed [10].

Figure 1: Left: Plot showing marginalized $\Delta\chi^{2}$ as a function of
$\delta m^{2}$ for KII, IMB, Baksan, and the combined set of three
experiments. Right: Sensitivity plot on the $\delta m^{2}-\sigma_{x}$ plane
for DUNE and HK, along with the exclusion contours from SN1987A for
comparison.

In the left panel of Fig. 1, we demonstrate the results of our fit through
$\Delta\chi^{2}\equiv-2(\ln{\mathscr{L}}-\ln{\mathscr{L}_{\rm max}})$, which
is calculated as a function of $\delta m^{2}$, by marginalizing over the other
parameters. In this analysis, values of $\delta m^{2}<10^{-21}{\rm eV}^{2}$
correspond to no active-sterile oscillations, as the oscillation lengths are
longer than what is relevant for a galactic SN. We find that KII (light-blue
dashed) prefers a non-zero $\delta m^{2}$ at a $\Delta\chi^{2}\approx 1.1$
over the non-oscillation hypothesis. This is primarily because KII observed
zero events in the energy range between $(21{\rm MeV},31{\rm MeV})$, which is
better explained in the PD hypothesis due to the presence of an oscillation
minima. Similarly IMB (green dotted) and Baksan (orange dot-dashed) also
prefer the PD scenario by a $\Delta\chi^{2}\approx 1.7$, although there is
some tension between their best fit values.

Finally, the combined data from all three experiments prefer a non-zero value
of $\delta m^{2}=6.31\times 10^{-20}{\rm eV}^{2}$, disfavoring the no-
oscillation scenario with $\Delta\chi^{2}\approx 3$. This is primarily
dominated by the tension in data between KII and IMB. Active-sterile
oscillations due to $\delta m^{2}$ allow to reconcile this difference, by
predicting a broader spectra for IMB, while leading to a suppression of events
observed in KII. Furthermore, the data can be used to exclude values
$2.55\times 10^{-20}{\rm eV}^{2}\lesssim\delta m^{2}\lesssim 3.0\times
10^{-20}{\rm eV}^{2}$ with a $\Delta\chi^{2}>9$. This makes it the lowest
values of $\delta m^{2}$ constrained by experiments till date, and provides
the strongest constraints on the PD hypothesis. However, one needs to do a
more dedicated Monte-Carlo analysis, before assigning a statistical
significance to the exclusion region.

###  3.1 Sensitivity studies for a future galactic supernova

This brings us to an important question – in the event of a future galactic
SN, where we are far better prepared to detect thousands of events, how will
the bounds on $\delta m^{2}$ improve? For the analysis, we consider DUNE and
HK, because of their large fiducial volume, and excellent energy
reconstruction. We assume the SN to take place at $10{\rm kpc}$. For the
analysis, we treat ${\cal E}_{\rm tot}$ and $E_{0}$ as free parameters, while
$\delta m^{2}$ is varied for a series of values of $\sigma_{x}$. This allows
us to get the sensitivity to both oscillations over large distances, as well
as to decoherence due to wave-packet separation.

Our results are depicted in the right panel of Fig. 1, which shows the $95\%$
sensitivity contours for DUNE (purple) and HK (green). We find that, for the
baseline chosen, the maximum sensitivity appears around $\delta m^{2}\sim
10^{-20}{\rm eV}^{2}$. For smaller values of $\delta m^{2}$, oscillation
baselines are longer than $10{\rm kpc}$, and hence sensitivity decreases.
However, for $\sigma_{x}\leq 10^{-14}\text{m}$, some sensitivity still exists
since a lower value of $\sigma_{x}$ cancels the effect of low $\delta m^{2}$.
On the other hand, for a larger $\delta m^{2}$, decoherence nullifies the
effect of rapid oscillations. Furthermore, for values of
$\sigma_{x}>10^{-13}\text{m}$, there is little sensitivity to values of
$\sigma_{x}$ since the coherence length remains larger than the relevant
baseline. We find that HK and DUNE results can be complementary, since HK can
probe smaller values of $\delta m^{2}$ due to its massive size, whereas DUNE
has more sensitivity to larger values of $\delta m^{2}$ due to its good energy
resolution. To compare with existing results, the region which can be explored
using data from SN1987A is shown in blue. Both DUNE and HK can easily probe a
large chunk of the excluded region from SN1987A, and therefore, will be able
to strengthen the results significantly.

##  4 Conclusion

The signatures of pseudo-Dirac neutrinos can be manifested through active-
sterile oscillations over long baselines. In this respect, we explored the
scenario that neutrinos are pseudo-Dirac, and studied the impact on the
neutrino fluence from SN1987A. Curiously, we found that the existing data from
SN1987A presents a mild preference for the pseudo-Dirac hypothesis.
Furthermore, using the same data, we excluded the tiniest values of mass-
squared differences considered so far, $2.55\times 10^{-20}{\rm
eV}^{2}\lesssim\delta m^{2}\lesssim 3.0\times 10^{-20}{\rm eV}^{2}$, with a
$\Delta\chi^{2}>9$. Finally, we used next generation experiments like DUNE and
HK to study the constraints that can be imposed on this scenario from a future
galactic core-collapse supernova. Needless to say, the occurrence of a
galactic supernova will open up a plethora of avenues for exploration in
neutrino physics.

## Acknowledgments

We thank Ivan Martinez-Soler and Yuber Perez-Gonzalez, in collaboration with
whom this work was performed. We would like to thank the Instituto de Fisica
Teorica (IFT UAM-CSIC) in Madrid for support via the Centro de Excelencia
Severo Ochoa Program under Grant CEX2020-001007-S, during the Extended
Workshop “Neutrino Theories”, where this manuscript was completed.

## References

## References

  * [1] A. de Gouvea, W. C. Huang and J. Jenkins, Phys. Rev. D 80 (2009), 073007 doi:10.1103/PhysRevD.80.073007 [arXiv:0906.1611 [hep-ph]].
  * [2] J. F. Beacom, N. F. Bell, D. Hooper, J. G. Learned, S. Pakvasa and T. J. Weiler, Phys. Rev. Lett. 92 (2004), 011101 doi:10.1103/PhysRevLett.92.011101 [arXiv:hep-ph/0307151 [hep-ph]].
  * [3] A. De Gouvêa, I. Martinez-Soler, Y. F. Perez-Gonzalez and M. Sen, Phys. Rev. D 102 (2020), 123012 doi:10.1103/PhysRevD.102.123012 [arXiv:2007.13748 [hep-ph]].
  * [4] I. Martinez-Soler, Y. F. Perez-Gonzalez and M. Sen, Phys. Rev. D 105 (2022) no.9, 095019 doi:10.1103/PhysRevD.105.095019 [arXiv:2105.12736 [hep-ph]].
  * [5] K. Hirata et al. [Kamiokande-II], Phys. Rev. Lett. 58 (1987), 1490-1493 doi:10.1103/PhysRevLett.58.1490
  * [6] R. M. Bionta, G. Blewitt, C. B. Bratton, D. Casper, A. Ciocio, R. Claus, B. Cortez, M. Crouch, S. T. Dye and S. Errede, et al. Phys. Rev. Lett. 58 (1987), 1494 doi:10.1103/PhysRevLett.58.1494
  * [7] E. N. Alekseev, L. N. Alekseeva, I. V. Krivosheina and V. I. Volchenko, Phys. Lett. B 205 (1988), 209-214 doi:10.1016/0370-2693(88)91651-6
  * [8] M. Kobayashi and C. S. Lim, Phys. Rev. D 64 (2001), 013003 doi:10.1103/PhysRevD.64.013003 [arXiv:hep-ph/0012266 [hep-ph]].
  * [9] F. Vissani, J. Phys. G 42 (2015), 013001 doi:10.1088/0954-3899/42/1/013001 [arXiv:1409.4710 [astro-ph.HE]].
  * [10] J. Kersten and A. Y. Smirnov, Eur. Phys. J. C 76 (2016) no.6, 339 doi:10.1140/epjc/s10052-016-4187-5 [arXiv:1512.09068 [hep-ph]].
  * [11] I. Tamborra, B. Muller, L. Hudepohl, H. T. Janka and G. Raffelt, Phys. Rev. D 86 (2012), 125031 doi:10.1103/PhysRevD.86.125031 [arXiv:1211.3920 [astro-ph.SR]].


# Non-Volatile Memory Accelerated Geometric Multi-Scale Resolution Analysis
Andrew Wood1, Moshik Hershcovitch2, Daniel Waddington2, Sarel Cohen3
Meredith Wolf4, Hongjun Suh5, Weiyu Zong1, Peter Chin1

1Boston University
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>2IBM Research
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>3Hasso Plattner Institute
<EMAIL_ADDRESS>4Williams College
<EMAIL_ADDRESS>5Seoul National University
<EMAIL_ADDRESS>
###### Abstract

Dimensionality reduction algorithms are standard tools in a researcher’s
toolbox. Dimensionality reduction algorithms are frequently used to augment
downstream tasks such as machine learning, data science, and also are
exploratory methods for understanding complex phenomena. For instance,
dimensionality reduction is commonly used in Biology as well as Neuroscience
to understand data collected from biological subjects. However, dimensionality
reduction techniques are limited by the von-Neumann architectures that they
execute on. Specifically, data intensive algorithms such as dimensionality
reduction techniques often require fast, high capacity, persistent memory
which historically hardware has been unable to provide at the same time. In
this paper, we present a re-implementation of an existing dimensionality
reduction technique called Geometric Multi-Scale Resolution Analysis (GMRA)
which has been accelerated via novel persistent memory technology called
Memory Centric Active Storage (MCAS). Our implementation uses a specialized
version of MCAS called PyMM that provides native support for Python datatypes
including NumPy arrays and PyTorch tensors. We compare our PyMM implementation
against a DRAM implementation, and show that when data fits in DRAM, PyMM
offers competitive runtimes. When data does not fit in DRAM, our PyMM
implementation is still able to process the data.

##  I Introduction

Dimensionality reduction is a set of algorithms designed to extract the
“useful” information that is often assumed to be embedded in a higher
dimensional signal. Dimensionality reduction techniques have been used with
great success in the field of machine learning [van2009dimensionality,
sakurada2014anomaly] and information retrieval [tsuge2001dimensionality,
venna2010information, kowalski2007information]. Some of the most popular
variants, such as the linear Singular Value Decomposition (SVD)
[golub1971singular] (and by extension Principal Component Analysis (PCA)
[ringner2008principal]), as well as nonlinear versions such as ISOMAP
[balasubramanian2002isomap] have even become standard items in a researcher’s
toolbox. Most recently, modern machine learning has pivoted to be more
“biologically” based: developing algorithms and models which ever increasingly
mimic the brain. Likewise, dimensionality reduction techniques have followed
the same trend. One technique based on the behavior of the cortex called
Geometric Multi-Scale Resolution Analysis (GMRA) [allard2012multi], learns
nonlinear manifolds in data by cleverly examining the local geometry of points
in the pointcloud.

However, dimensionality reduction algorithms such as GMRA have historically
been limited by the von-Neumann architectures they execute on: primarily
constrained by small system memory (DRAM) capacity. In contrast, modern
research has become dominated by “big data”, with datasets such as ImageNet
[deng2009imagenet], Cifar100 [krizhevsky2009learning], MusicNet
[thickstun2018invariances, trabelsi2017deep], WMT14 [axelrod2011domain] being
exponentially larger than any in history.

In the last few years, a notable hardware breakthrough has been the invention
of Intel Optane Persistent Memory Modules111Herein termed Optane-PM . Optane-
PM in particular communicates via the memory bus, circumventing bottlenecks
such as PCIe lanes, using the same interface to the CPU as DRAM. While there
are other Persistent Memory technologies, Optane-PM is the most mature product
on the market. Optane-PM is based on 3D XPoint (3DXP) technology and operates
at a cache-line granularity with a latency of around 300ns [spectra2020,
izraelevitz2019basic]. While this latency is slower than current DRAM
($\scriptstyle\mathtt{\sim}$100ns), it is 30x faster than the current state of
the art NVMe SSDs. Additionally, a single DIMM of Optane-PM can reach 512GB,
which is 8x larger than the available DRAM. Thus, the maximum Optane-PM
capacity of a commodity 2U server machine is 12TB - significantly more than
DRAM.

By default, Optane-PM can be used in three modes: block device mode in which
it acts as a faster storage device, Memory Mode in which it acts as high
capacity volatile memory cached by DRAM, and App Direct mode which allows
load/store access to the device which is persistent. For a program to use
Optane-PM in block device mode, the program must still pay an OS I/O penalty
to access data on the device. In Memory Mode, a program does not pay this
penalty, but must copy-off to additional storage devices if persistence is
required. App Direct mode treats the memory as separate from volatile memory
in the system. Optane-PM is exposed as a DAX (direct access) device. Once it
has been mapped into the process virtual address space, it can be accessed
directly through load/store. Under-the-hood, data is moved around in cache
lines and 256B blocks internally to the Optane-PM device. Although App Direct
mode is very powerful, it does require programs to be aware of this ”special”
memory and pay close attention to crash-consistency and data recovery
(normally requiring access to low level persistent memory and software
transaction support libraries).

While Optane-PM technology is still in its infancy, there are some early
adopters in databases (see [SAPHANA]). However, to the best of our knowledge,
we are the first to exploit this technology to accelerate dimensionality
reduction algorithms such as GMRA.

In this paper, we use components of MCAS [mcas2021], a key-value store
framework built from the ground up for persistent memory, to implement GMRA.
Specifically, we use a simplified MCAS library integrated into the Python
programming language called Python Micro-MCAS (PyMM [WaHeDi21PyMM]). PyMM uses
local persistent memory resources only, unlike the broader MCAS solution which
is a network-attached programmable key-value store.

One of the key features of PyMM is its support for native Python types such as
NumPy ndarray and PyTorch tensor. We use PyMM to re-implement GMRA. In total,
we were required to only 15 lines of code change (out of $1000+$ lines of
code) to integrate our existing implementation with PyMM. We then show that
despite the higher latency of Optane-PM , PyMM-GMRA is still comparable to its
DRAM only counterpart when data fits into DRAM. Additionally, short of writing
your own Python interface to persistent memory, PyMM provides a
straightforward solution when data is larger than available DRAM.
Specifically, we evaluate our performance on the MNIST [lecun1998mnist] and
Cifar10 [krizhevsky2009learning] datasets.

##  II Background and Related Work

###  II-A MCAS: Memory Centric Active Storage

Figure 1: The high level architecture of MCAS. MCAS abstracts Optane in DAX
mode to a key-value store api via the network which allows clients to
load/store memory via the system memory bus while maintaining data
persistence. Image taken from [mcas2021].

MCAS is a networked key-value store designed to interface with Optane devices
in DAX mode. While the MCAS architecture can be seen in Figure 1, it can
generally be thought of as an abstraction layer between the persistent memory
and a client application. A key feature of MCAS is the support of in-store
operations. This feature, called Active Data Objects (ADO)s, allows code to
run directly on the MCAS server without needing to transfer values to the
client or move data objects into DRAM. Thus, MCAS enabled systems allow client
applications to modify existing data objects, create new data objects, and
delete data objects without the need to pay runtime penalties for
copying/transferring data. Additionally, if persistence is enabled, all data
objects will live after program execution and after system reboots. Therefore,
programs can avoid paying additional I/O runtime penalties by using MCAS. We
note that programs will also need to add crash-consistency code to ensure data
recovery in the event of a system crash, this behavior is not enabled by
default.

MCAS currently supports pure Python ADOs which execute Python code directly on
the MCAS server. However, this api is not easy to use, and requires
considerable skill by the developers to write correct applications. Therefore,
a novel abstraction library called PyMM was developed to encapsulate MCAS
Python calls as well as support native Python types.

Figure 2: The PyMM architecture. PyMM exports micro-MCAS, a networkless
version of MCAS, to the Python programming language, and also provides support
for native Python types as well as NumPy arrays and PyTorch tensors. PyMM
abstracts away memory management from the client through the use of a shelf
datatype which can be seen above.

###  II-B PyMM: Python Micro-MCAS

As previously mentioned, PyMM is a new Python library that encapsulates the
MCAS api and provides a significantly easier interface. Specifically, PyMM is
built on top of a networkless MCAS implementation which is called Micro-MCAS.
Additionally, PyMM includes built-in support for native Python types such as
NumPy arrays and PyTorch tensors. While arbitrary datatypes can be stored on
an MCAS server, the client is responsible for writing the code that interfaces
their datatypes with the persistent memory such as writing metadata, etc. For
NumPy and PyTorch, most linear algebra operations would require custom code.
PyMM handles this for the client, allowing the client to use existing NumPy
algorithms as-is, and without major modification to their existing
NumPy/PyTorch code. The PyMM stack can be seen in Figure 2.

PyMM is built around a shelf datatype. A shelf represents a logical grouping
of data objects which are held in a specific set of memory resources. In our
case, the memory resource is Optane acting in DAX mode, enabling both high
memory speeds and persistence. Therefore, PyMM shelf objects will exist across
program executions and system resets. Data objects can also be exchanged
between DRAM or any other memory device which is directly attached to the CPU
via the memory bus, including GPUs. Shelf objects can be created in two modes:
_devdax_ (visible as a character device file) or _fsdax_ (visible as block
device files). In our case, we create shelf objects in fsdax mode.

In our experiments, we use PyMM’s built-in support for NumPy arrays which
correspond to the shelf type pymm.ndarray. This shelf type is polymorphic with
NumPy’s native ndarray type, meaning that all existing NumPy algorithms will
work (and operate in-place) with a shelf.

For data-intensive algorithms such as dimensionality reduction, fast system
memory with high capacity and persistence is required. For dimensionality
reduction algorithms, the entire dataset is usually required to exist in
memory at once, as global knowledge is needed. These algorithms commonly
construct complicated internal representations of the data which they then use
to perform dimensionality reduction. For instance, SVD and PCA require
computing the gram matrix ($A^{T}A$ or $AA^{T}$), ISOMAP and LDA build
diffusion maps, and GMRA uses approximate nearest neighbor strategies. These
complicated data structures are expensive to build, and consume memory
proportional to the number of data points. Therefore, dimensionality reduction
techniques require persistence: they only build these structures once for a
specific dataset, and then save that structure for future use/reuse. By using
MCAS, the I/O penalty for serializing these complicated structure to disk can
be avoided.

While there has been work on building these structures in a streaming fashion
[rachakonda2016memory, huang2020memory, mitliagkas2013memory,
ghashami2015frequent], these algorithms do not provide general case solutions
for all dimensionality reduction techniques, and they often are approximate
algorithms [huang2020memory, ghashami2015frequent]. By using PyMM (and
therefore MCAS), dimensionality reduction techniques are able to process large
datasets that were not otherwise possible without using streaming
approximations. By specifically using PyMM, these algorithms require the
fewest changes to their existing codebases in order to benefit from MCAS.

We note that a separate project, called PMDK [PMDK] also provides persistent
memory access to Optane via the memory bus. A C/C++ wrapper called PyNVM
[PyNVM] provides PMDK functionality to Python and allows for custom datatypes
to be compatible with PMDK. However, PyNVM does not use a subclassing strategy
and also does not support flat datatypes such as NumPy ndarrays or PyTorch
tensors. Thus, we chose to use PyMM and MCAS for our implementation.

Figure 3: A visualization of the linear approximation of a point using a basis
function at scale $j$, the approximation at a finer scale $j$+1 and difference
operator (geometric wavelet) from scale $j$ to scale $j+1$. This image was
taken from [allard2012multi].

###  II-C Geometric Multi-Scale Resolution Analysis (GMRA)

GMRA is a dimensionality reduction technique that is inspired by the cortex.
At the microscale, neurons in the cortex fire, which induces synchrony between
neurons. Neural synchrony at the macroscale produces patterns, which combined
with other firings, have long been believed to be the intermediary
representation of data [dicarlo2012does]. Macroscale firings can also be
viewed as dimensionality reduction: at the macroscale, due to neural
synchrony, at any given point only a subset of the neurons are firing.
Therefore, by viewing ever larger populations of neurons, it is believed that
the cortex computes lower-dimensional representation of data at different
scales in increasing layers of abstraction. GMRA mimics this behavior by
processing a point cloud at different scales to produce increasingly fine-
grained manifolds.

The GMRA algorithm contains three steps to compute these manifolds at
different scales:

  1. 1.

It computes a leveled tree decomposition of the manifold into dyadic cells.
Dyadic cells have the following properties:

    1. (a)

Each dyadic cell contains a subset of the pointcloud, and the entire subset
exists within a sphere with a fixed radius.

    2. (b)

The children of a dyadic cell divide the points contained in the parent into
disjoint subsets.

    3. (c)

The children of a dyadic cell cover the points inside the parent.

It is worth noting that due to these rules, each level of the tree covers the
entire point cloud, and all cells at a level are disjoint. The root node of
the tree contains the entire point cloud, and the final level of the tree
contains the finest-scale (smallest radii) grouping of the points.

  2. 2.

It computes a $d$-dimensional affine approximation for each dyadic cell. This
approximation represents the basis of each dyadic cell and is a linear
piecewise approximation (i.e. the SVD decomposition of the cell’s covariance).

  3. 3.

It computes a sequence of low-dimensonal affine difference operators that
encode the difference between subsequent levels of the tree (i.e. scales).
These difference operators allow efficient querying of the points by scale.

GMRA essentially searches for manifolds by grouping points according to their
local topology and produces piecewise linear approximations of the groups.
First, GMRA considers the entire point cloud as a group and computes a linear
approximation of the entire dataset using SVD [golub1971singular]. Then, it
iteratively partitions the point cloud into smaller groups by finding open
balls which contain the maximum number of points, where it computes linear
approximations for those groups, producing the next level of cells in the
tree. This cycle repeats until some halting criteria is met. To query the
tree, GMRA computes difference operators between levels of the tree which
allow a query to start at the root, lookup the embeddings for the query, and
then walk the tree, updating the embeddings using the linear difference
operators (a process which can be seen in Figure 3).

By computing the linear approximation of a cell using SVD of the covariance,
the approximation geometrically becomes a parallel plane that passes through
the mean of the cell and is approximately tangent to the local area of the
manifold in the cell. As seen in Figure 3, the linear approximation, called
the scaling function fit to each group can be queried by starting at scale 0
(the roughest scale), getting the approximation for the query at that scale,
and then applying the difference operator (wavelet correction) to get the
approximation at the next scale (finer scale). Since each level of the tree
represents the decomposition of the point cloud at a scale, by walking from
the root to the child at the appropriate level, GMRA produces low-dimensional
embeddings for each point at arbitrary scales. GMRA has been used to improve
performance of downstream tasks [tran2014geometric, chen2012fast] such as
classification and anomaly detection.

##  III Our Approach

###  III-A GMRA Implementation

The GMRA tree decomposition step can be implemented in a few ways
[allard2012multi]:

  1. 1.

Use approximate nearest neighbors to construct a weighted graph where the
nodes represent points in the point cloud, and draw weighted edges from a node
to its $k$ nearest neighbors. The weight of this edge takes the form
$e^{\frac{||x_{i}-x_{j}||_{2}^{2}}{\sigma}}$ where $k$ and $\sigma$ are
hyperparameters. In practice, $k$ is taken between 10 and 50, and $\sigma$ is
typically the distance between $x_{i}$ and its
$\left\lfloor{\frac{k}{2}}\right\rfloor$ nearest neighbor.

By partitioning this weighted graph (e.g. METIS routines [karypis1998fast]),
the dyadic cells are constructed as the vertex groups left after each
separation, while the tree structure is given by the action of separation. For
instance, the partitioning algorithm is run recursively on all previous
partitions, yielding the dyadic cells and their parent-child relationships.

  2. 2.

Use the CoverTree algorithm [beygelzimer2006cover] to partition the data into
a level tree based on radii at different scales. The CoverTree algorithm
builds a graph where each node corresponds to a single point in the
pointcloud, and draws a directed edge $x_{i}\rightarrow x_{j}$ if $x_{j}$ is
within radius $2^{s}$ from $x_{i}$ (as measured by euclidean distance).

Since a CoverTree has a root node, the dyadic cells are constructed by
computing all vertices reachable within the radius $2^{j}$ at scale $j$ given
the root of the subtree. The procedure then recurses for all children of the
root reachable by radius $2^{j-1}$ for the next scale $j-1$. Note that in
order to use the CoverTree construction, the max scale must be known apriori
(to the construction) and can be an expensive operation. However, CoverTree
construction is the only approach which provides the theoretical bounds for
the GMRA algorithm.

  3. 3.

Iterated PCA. At the largest scale, compute the top $d$ principle components,
and then assign the points into two children based on the sign of their
$d+1$-st singular vector. The tree relationship is constructed by recursing on
the two children.

  4. 4.

Iterated $k$-means clustering. At the largest scale, partition the data into
$d$ clusters, and recurse on each cluster.

Our implementation currently supports the CoverTree construction. The original
implementation of GMRA is in MATLAB and uses custom C/C++ code which are
distributed in binary form. This has the unfortunate side effect of locking
the GMRA-MATLAB implementation to a specific, outdated version of MATLAB using
outdated (and in some cases, unsupported) third party libraries. Additionally,
GMRA-MATLAB will not interface with MCAS and therefore is constrained by DRAM.

Our implementation also makes use of custom C/C++ code, however our C/C++ is
distributed in source code form and is exported to Python using libtorch (the
C++ library of PyTorch [paszke2019pytorch]). Most importantly, our
implementation interfaces with PyMM and therefore can leverage its PyTorch and
NumPy support. Specifically, we implemented the tree construction using the
CoverTree algorithm in C++ for speed and parallelism, and then exported our
code to Python using Pybind11 [jakob2017pybind11]. Final wavelet construction
occurs in pure Python with full PyMM support. Our implementation demonstrated
PyMM’s seamless support for existing NumPy routines. For instance, we used
NumPy’s svd algorithm and scipy’s qr algorithm with pymm.ndarray datatypes on
persistent memory.

Name | # examples | dimensionality | # labels |
---|---|---|---|---
MNIST | 70k | 784 | 10
Cifar10 | 60k | 3072 | 10
TABLE I: A breakdown of the datasets. Note that while Cifar10 has fewer
examples, the dimensionality of each example is significantly higher.

Dataset | Data Loading DRAM (s) | Data Loading PyMM(s) | Wavelet Construction DRAM(s) | Wavelet Construction PyMM(s) | num trials
---|---|---|---|---|---
MNIST | $0.33\pm 2.02\times 10^{-3}$ | $0.53\pm 1.73\times 10^{-3}$ | $62.4\pm 1.88$ | $71.2\pm 1.55$ | 20
Cifar10 | $0.83\pm 2.58\times 10^{-3}$ | $1.2\pm 8.26\times 10^{-3}$ | $134.7\pm 5.07$ | $140.6\pm 5.24$ | 20
TABLE II: Top: Timing results of our experiments. The top two figures show the
timings from loading data and wavelet construction respectively. The bottom
two figures aggregate these results into the total time of using DRAM/PyMM and
the total running time of the experiments respectively. Note that while PyMM
performs much worse during the data loading phase, this phase comprises of 1%
of the total time spent using PyMM/DRAM. Bottom: Expected timing values
including standard deviations. MNIST took 209.35MB on the shelf and Cifar10
took 703.12MB on the shelf.

##  IV Experiments

We took two datasets, MNIST [lecun1998mnist] and Cifar10
[krizhevsky2009learning] for our experiments. Both datasets contain images,
but with different number of examples, dimensionality, and content. Most
notably the Cifar10 dataset contains larger, color images while the MNIST
dataset contains smaller, grayscale images. A breakdown of the datasets
(preprocessed and in original form) is shown Table I.

For our experiments, we measured the runtime that the GMRA algorithm took to
process both the MNIST and Cifar10 datasets. We computed the CoverTree
structure outside of our expeiments, as it only needs to occur once, is shared
between the PyMM and DRAM versions, and only uses DRAM. We then serialized the
CoverTree to disk. To compare DRAM to PyMM, we ran the tree deserialization,
Dyadic cell computation, and wavelet computation using both a DRAM and PyMM
enabled version of GMRA. Since the PyMM api is similar to the raw NumPy api,
and supports NumPy api operations on PyMM-stored arrays, the differences
between the two implementations was minimal. Our experiments were run with
PyMM crash consistency disabled.

We ran the partial GMRA algorithm (CoverTree deserialization through wavelet
construction) 10 times each using DRAM and PyMM, and recorded timing splits
for each stage of the script. We note that changing between DRAM and PyMM has
minimal effect on the Dyadic cell construction (due to this operation always
occurring in DRAM and being minorly influenced by L2-L4 cache misses).
However, all major linear algebra operations occur during the wavelet
construction stage, and this operation will either occur in DRAM or PyMM as
well as storing the dataset itself. In our experiment, we included the cost of
moving the dataset from DRAM/disk to PyMM in the PyMM reported times.

In our experiments, we used a server equipped with two Intel Xeon Gold 6248
processors containing 80cpu cores running at 2.5GHz base clock. The server is
also equipped with 384GB of DDR4 DRAM and 786GB of Optane DC. This server also
has a NVIDIA Tesla m60 GPU, however we did not make use of it in our
experiments.

##  V Results and Discussion

In our experiments, we observed that our PyMM GMRA implementation performed
slower than the DRAM implementation. This behavior is expected, and results
can be seen in Table II. We were however surprised at how fast the PyMM
implementation was compared to how slow Optane memory is with respect to DRAM
speeds. In our experiments, we observed a slowdown of 1.6x on MNIST and a
1.45x slowdown on Cifar10 for loading data into PyMM versus DRAM. This finding
is surprising given the reported latency: Optane is 3x slower than DRAM. For
wavelet construction; the stage where PyMM has the greatest influence over
GMRA, we observed a slowdown of 1.14x on MNIST and a 1.04x slowdown on
Cifar10. Our timing results can also be seen in Table II.

While the slowdown is almost negligible for MNIST and Cifar10, it is worth
noting that for larger datasets, the slowdown can amount to hours of
additional runtime. However, there is a hidden cost that we did not measure in
our experiments: the penalty for serializing the GMRA wavelets to disk. This
I/O penalty only applies to the DRAM version since PyMM is persistent, and it
is unlikely that this serialization would close the runtime gap on large
datasets. However, for downstream tasks, the I/O penalty of using a DRAM
system would be nontrivial, since there is no guarantee that the lower
dimensional dataset would fit into DRAM, and that read speeds from persistent
storage such as a NVMe SSD or HDD would be significantly slower than load
speeds from PyMM.

Our experiments had crash consistency disabled. Crash consistency is
implemented via a software logging strategy in PyMM that occurs at every write
operation. While crash consistency is important for our algorithms, it induces
a copy operation every time a data object is written. Since this policy is in
software, this adds significant overhead to the runtime of any algorithm which
uses PyMM if crash consistency is enabled. We omitted crash consistency at
this point because our algorithm is not crash consistent. In order to make our
algorithm crash consistent, it would require additional code changes which
would create an unfair comparison between our PyMM implementation and a
(already crash inconsistent) DRAM implementation.

##  VI Future Work

One major improvement that we leave to future work is to integrate our custom
C/C++ tree structures into PyMM. Currently, our tree structures (CoverTree and
Dyadic cell tree) are implemented on DRAM and execute in series. In the
future, we would like to put both of these structures on persistent memory
using PyMM and operate on them in-place. A further direction is to add
parallelization to our CoverTree implementation. CoverTrees are highly
parallelizable, however we omitted this in our current implementation in order
to simplify the code.

Additionally, we envision a scenario where GMRA and other dimensionality
reduction techniques are supported by all the features available in MCAS.
Currently, PyMM can only access local storage devices, which caps the capacity
to how many Optane-PM drives can fit inside a single machine. Fully-fledged
MCAS on the other hand provides a network abstraction, allowing memory to be
shared across multiple machines. This would enable our algorithms to scale-out
across nodes and process datasets currently only possible by large Hadoop
clusters.

Our last avenue of future work is to utilize GPU resources. By transferring
data between persistent memory and GPU memory, we can further accelerate our
algorithms either in down stream tasks like using the lower-dimensional
datasets for learning, or by using GPU resources in our algorithms to speed up
calculations. CXL-attached GPU and FPGA accelerators combined with persistent
memory might also be an interesting direction for future research.

##  VII Conclusion

In conclusion, we have demonstrated that MCAS and PyMM provide a suitable
framework to implement data intensive algorithms such as GMRA. While
persistent memory has higher latency than DRAM when the dataset fits into
DRAM, PyMM-GMRA still provides competitive runtimes. When data does not fit
into DRAM, PyMM-GMRA is the only implementation which can process this data.
We demonstrate that with PyMM, data intensive algorithms have new
opportunities to escape the classical restrictions of von-Neumann
architectures and can be used more effectively in the world of big data.