# Deriving the Dark Matter-Dark Energy Interaction Term in the Continuity
Equation from the Boltzmann Equation

Kevin J. Ludwick<EMAIL_ADDRESS>Department of Chemistry and Physics,
LaGrange College, LaGrange, GA 30240, USA  Holston Sebaugh
<EMAIL_ADDRESS>Department of Chemistry and Physics, LaGrange College,
LaGrange, GA 30240, USA

###### Abstract

Dark energy and dark matter are two of the biggest mysteries of modern
cosmology, and our understanding of their fundamental nature is incomplete.
Many parameterizations of couplings between the two in the continuity equation
have been studied in the literature, and observational data from the growth of
perturbations can constrain these parameterizations. Assuming standard general
relativity with a simple Yukawa-type coupling between dark energy and dark
matter fields in the Lagrangian, we use the Boltzmann equation to analytically
express and calculate the interaction kernel $Q$ in the continuity equation
and compare it to that of a typical parametrization. We arrive at a comparably
very small result, as expected. Since the interaction is a function of the
dark matter mass, other observational data sets can be used to constrain the
mass. This calculation can be modified to account for other couplings of the
dark energy and dark matter fields. This calculation required obtaining a
distribution function for dark energy that leads to an equation of state
parameter that is negative, which neither Bose-Einstein nor Fermi-Dirac
statistics can supply, and this is the main result of this paper. Treating
dark energy as a quantum scalar field, we use adiabatic subtraction to obtain
a finite analytic approximation for its distribution function that assumes the
FLRW metric and nothing more.

Introduction

It is now well known that some form of dark energy, comprising about 70% of
our universe, is responsible for cosmic acceleration, and that dark matter is
the next most prevalent non-luminous substance in our universe, comprising
about 25% of it. The fundamental nature of dark matter and dark energy and how
they interact with each and the Standard Model is uncertain. However,
interaction between dark matter and dark energy can be constrained by the
matter power spectrum 1809.00550 ; 1809.02411 ; 1812.05493 ; 1812.06854 ;
1812.05594  . Usually, an ad hoc parametrization for an interaction between
dark energy (DE) and dark matter (DM) as perfect fluids is assumed since we
don’t know of any fundamental coupling between them. Conservation of energy-
momentum implies the continuity equation:

$\nabla_{\alpha}T^{\mu\nu}=0~{}~{}~{}\rightarrow~{}~{}~{}~{}\Sigma_{i}(\dot{\rho_{i}}+3H(\rho_{i}+p_{i}))=0,$
(1)

where $\rho$ is the energy density, $p$ is the pressure, and $H$ is the Hubble
parameter, and the sum is over the DE and DM components of the universe for
late cosmological times since they dominate. If there is an interaction
between DE and DM, we would have

$\dot{\rho_{DM}}+3H(\rho_{DM}+p_{DM})=Q,~{}~{}~{}\dot{\rho_{DE}}+3H(\rho_{DE}+p_{DE})=-Q.$
(2)

Interaction between DE and DM can ameliorate the coincidence problem and the
$H_{0}$ tension, and different parameterizations for the interaction kernel
$Q$ have been widely studied in the literature 1908.03663 ; 9908023 ;
1812.06854 ; 1310.0085 ; 1603.08299 ; 1706.04953 ; 1809.06883 ; 1805.08252 ;
1905.08286 ; 1908.04281 ; 1501.06540 ; 1502.04030 ; 1501.03073 ; 1807.05541 ;
1905.10382 ; 0606520 ; 1507.00187 ; 0609597 ; 0610806 ; 0702015 ; 0801.4233 ;
9408025 . Influences for the parametrizations in the literature include
computational convenience, ease of data comparison, and analogies from field
theory. The uncertain nature of the field theory of dark energy helps motivate
the phenomenological approach to interaction kernels 1110.3045 ; 0706.3064 .
These parametrizations are usually convenient choices that allow for analytic
solutions of Friedmann’s equations. In this work, we are interested in what
$Q$ would look like if a coupling on a more fundamental level were assumed,
namely, a coupling as fields in the Lagrangian. Field couplings between DE and
DM have been studied in the literature 9908023 ; 1603.08299 ; 2006.04618 ;
1605.00996 ; 1910.02699 ; 1808.05015 ; 1411.3660 , but what we are interested
in is a way of expressing an interaction in terms of a distribution function
that only takes effect at 1st order and above in cosmological perturbation
theory, and this scheme is typically carried out by expressing the energy
density and pressure of fluid components in terms of a weighted integral of
the distribution function 9506072 . Field coupling interactions are then taken
into account via the collision term in the Boltzmann equation.

Typically, dynamical DE is modeled as a scalar field, and the usual way of
ensuring a slowly varying field that provides sufficient late-time cosmic
acceleration in accordance with observations is to let the mass be on the
order of $H_{0}$ or less. In this case, if DE were in thermodynamic
equilibrium with the rest of the contents of the early universe, it would be a
Bose-Einstein condensate 1411.0753 , and expressing the energy density and
pressure as weighted integrals of a distribution function would not be valid
since the majority of the energy would be in the ground state; a discrete sum
instead of a continuous integral would be needed. However, we assume here that
DE is not a condensate, and our methodology that follows serves more as a
proof of concept for DE and other general scalar field theories rather than a
thorough investigation of the nature of DE.

We let the scalar field $\phi$ represent DE and the scalar field $\psi$
represent DM. Assuming a Yukawa-type coupling (last term) in the action for
late times in which DE and DM dominate, the action is

$S=\int d^{4}x\sqrt{-g}\left[\frac{R}{16\pi
G}-\frac{1}{2}g^{\mu\nu}\nabla_{\mu}\phi\nabla_{\nu}\phi-V(\phi)-\frac{1}{2}\xi
R\phi^{2}-\frac{1}{2}g^{\mu\nu}\nabla_{\mu}\psi\nabla_{\nu}\psi-V(\psi)-\frac{1}{2}\chi
R\psi^{2}-\frac{1}{2}g\psi^{2}\phi\right],$ (3)

in which $\xi$ and $\chi$ are the non-minimal coupling parameters for DE and
DM respectively, and they are in general present formally due to
renormalization of the scalar fields. The Yukawa-type coupling parameter $g$
is assumed to be small and perturbative.

In what follows, instead of assuming a certain parametrization for $Q$, we
calculate what $Q$ should be for this Yukawa-type coupling by utilizing the
tree-level cross section with the Boltzmann equation. In theory, a weakly
interacting scalar field would have a pressure given by a weighted sum or
integral involving the Bose-Einstein statistics factor, which would result in
a positive pressure. However, DE pressure must be negative when modeled as a
perfect fluid according to the cosmic acceleration requirements. It is also
known that the correspondence between a scalar field theory and a perfect
fluid representation is not perfect 1201.1448 . In order to deal with this
incompatibility, we must find a suitable effective distribution function for
DE that results in a negative pressure. We use adiabatic expansion in the next
section to arrive at an approximate expression for the distribution function
for DE, and we calculate $Q$ in the section after that. Then we analyze our
result and compare it to observational constraints on a typical
parameterization in the literature, and then we conclude.

In this work, we assume the FLRW metric
$ds^{2}=-dt^{2}+a(t)^{2}(dx^{2}+dy^{2}+dz^{2})$, and we assume
$\hbar=c=k_{B}=1$.

Distribution Function for Dark Energy

In this section, we calculate an effective distribution function for dark
energy that is compatible with negative pressure. We do this by quantizing the
scalar field theory for dark energy and finding the renormalized expectation
value for the energy density obtained from the stress-energy tensor. We
regularize and renormalize using adiabatic subtraction. We then compare our
result with the expression for energy density from statistical mechanics to
identify the effective distribution function that we will later use in the
Boltzmann equation.

The dark energy fluid must have a negative pressure, but DE as a scalar field
would have the Bose-Einstein distribution function, which would result in
positive pressure according to the statistical mechanics expression for
pressure,

$p=\frac{1}{a^{4}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{k^{2}}{3\sqrt{k^{2}+a^{2}m^{2}}}f(k),$
(4)

where $a$ is the scale factor from the FLRW metric and $f(k)$ is the
distribution function. The local momentum $k$ is defined as the scale factor
$a(t)$ times the proper momentum. We use the definitions for stress-energy
components according to 9506072 , which account for the curvature of FLRW
space since the negative pressure of dark energy is a large-scale effect, and
we quantize the field below with this in mind as well.

To derive a suitable distribution function for scalar-field dark energy, we
find the vacuum expectation value of the energy density from the stress-energy
tensor in Einstein’s equation and equate it with the expression for the energy
density,

$\rho=\frac{1}{a^{4}}\int\frac{d^{3}k}{(2\pi)^{3}}\sqrt{k^{2}+a^{2}m^{2}}f(k),$
(5)

to obtain a distribution function $f_{\phi}$ for the dark energy scalar field
$\phi$.

Assuming a perfect fluid model for dark energy that may be non-minimally
coupled to the metric via the coupling $\xi$, the stress-energy tensor
component for the energy density is

$T_{00}=\rho_{\phi}=\frac{1}{2}\dot{\phi}^{2}+\frac{1}{2a^{2}}(\partial_{i}\phi)^{2}+V(\phi)+3\xi
H^{2}\phi^{2}+3\xi H\frac{\partial}{\partial t}\left(\dot{\phi}^{2}\right),$
(6)

where the ”dot” denotes a derivative with respect to the $t$ coordinate of the
metric. We use $V(\phi)=\frac{1}{2}m^{2}\phi^{2}$, which is compatible with
$w\leq-1/3$ for late cosmological times 1804.02987 . In order to find the
expectation value of $T_{00}$, we must quantize the field ParkerToms ;
1807.10361 :

$\phi=\frac{1}{\sqrt{2(2\pi a)^{3}}}\int
d^{3}k\left[A_{\vec{k}}e^{i\vec{k}\cdot\vec{x}}h_{\vec{k}}(t)+A^{\dagger}_{\vec{k}}e^{-i\vec{k}\cdot\vec{x}}h^{*}_{\vec{k}}(t)\right],~{}~{}~{}u\equiv\frac{1}{\sqrt{2(2\pi
a)^{3}}},$ (7)

where

$[A_{\vec{k}},A^{\dagger}_{\vec{k^{\prime}}}]=(2\pi)^{3}\delta(\vec{k}-\vec{k^{\prime}}),~{}~{}~{}A_{\vec{k}}\lvert
0\rangle=0,~{}~{}~{}\langle 0\rvert A^{\dagger}_{\vec{k}}=0.$ (8)

Using Eqs. (7) and (8) with Eq. (6), we obtain

$\displaystyle\langle T_{00}\rangle=\rho_{\phi}=(2\pi)^{3}\int d^{3}k\bigg{[}$
$\displaystyle\frac{1}{2}(\dot{u}^{2}|h_{\vec{k}}|^{2}+\dot{u}u(h_{\vec{k}}\dot{h}_{\vec{k}}^{*}+\dot{h}_{\vec{k}}h_{\vec{k}}^{*})+u^{2}\dot{h}_{\vec{k}}\dot{h}_{\vec{k}}^{*})+\frac{u^{2}}{2a^{2}}k^{2}|h_{\vec{k}}|^{2}+\frac{1}{2}u^{2}m^{2}|h_{\vec{k}}|^{2}$
$\displaystyle+3u^{2}\xi H^{2}|h_{\vec{k}}|^{2}+6\xi
H(u\dot{u}|h_{\vec{k}}|^{2}+\frac{u^{2}}{2}(h_{\vec{k}}\dot{h}_{\vec{k}}^{*}+\dot{h}_{\vec{k}}h_{\vec{k}}^{*})\bigg{]}.$
(9)

To render this expectation value finite, we regularize and renormalize it via
adiabatic subtraction. Minimizing the action in Eq. (3) with respect to $\phi$
results in the equation of motion

$\ddot{\phi}+3H\dot{\phi}+\frac{dV}{d\phi}+\xi R\phi=0,$ (10)

which assumes that the DM field $\psi$ and the DE field $\phi$ are not
dependent on each other since the equation is equivalent to the second part of
Eq. (2) with $Q=0$, and this is approximately correct since we the cross
section between DM and DE is small when DE has very small mass, which is what
we will come to shortly. And for any direct coupling between DE and DM in the
Lagrangian, one can treat the interaction term in the Lagrangian as
perturbative so that Eq. (10) is at least still valid at 0th order.
Substituting Eq. (7) into Eq. (10), one gets a perturbative solution in
adiabatic orders (i.e., orders of derivative of the metric, where a new
derivative introduces an extra factor of time coordinate $T$ in the
denominator) ParkerToms :

$h_{\vec{k}}=W^{-1/2}e^{-i\int^{t}Wdt^{\prime}}+O(T^{-2(n+1)}),$ (11)

where $W$ has non-zero terms for even adiabatic order only, and $W$ in the
equation is given to $2n$th adiabatic order. $W$ can be found iteratively
order by order from the relation

$W^{2}=k^{2}/a^{2}+m^{2}+(6\xi-3/4)(\dot{a}/a)^{2}+(6\xi-3/2)\ddot{a}/a+W^{1/2}\frac{d^{2}}{dt^{2}}W^{-1/2}.$
(12)

Using Eq. (11) in Eq. (Deriving the Dark Matter-Dark Energy Interaction Term
in the Continuity Equation from the Boltzmann Equation) and equating Eq.
(Deriving the Dark Matter-Dark Energy Interaction Term in the Continuity
Equation from the Boltzmann Equation) with Eq. (5) (a similar process to what
is done in 1804.07471 ), we find that the distribution function for the DE
field $\phi$ is

$\displaystyle f_{\phi}=$
$\displaystyle\frac{(2\pi)^{6}a^{4}}{\sqrt{k^{2}+a^{2}m^{2}}}\bigg{[}\frac{1}{2}\bigg{(}\dot{u}^{2}W^{-1}-\dot{u}uW^{-2}\dot{W}+u^{2}(\frac{1}{4}W^{-3}\dot{W}^{2}+W)\bigg{)}+\frac{u^{2}k^{2}}{2a^{2}}W^{-1}+\frac{1}{2}u^{2}m^{2}W^{-1}$
$\displaystyle+3u^{2}\xi H^{2}W^{-1}+6\xi
H(u\dot{u}W^{-1}-\frac{u^{2}}{2}W^{-2}\dot{W})\bigg{]}.$ (13)

The method of adiabatic subtraction ParkerToms  says that for a given quantity
$A$, the physical, finite expression for $\langle A\rangle$ is $\langle
A\rangle=\sum\limits_{i}^{\infty}\langle A\rangle_{i}-\langle
A\rangle_{\rm{divergent~{}orders}}$. For a slowly varying (i.e., adiabatic)
FLRW spacetime, each term in adiabatic order contributes less than the
previous term in adiabatic order ParkerToms . So for a given quantity $A$, we
can approximate by truncating the infinite sum to some sufficiently high
adiabatic order. It turns out that our calculation involving the Boltzmann
equation for the interaction term $Q$ will result in divergent terms for 6th
adiabatic order and lower; we truncate our expression to 8th order. We can
express each of the terms in Eq. (Deriving the Dark Matter-Dark Energy
Interaction Term in the Continuity Equation from the Boltzmann Equation) as a
sum of adiabatic orders, accordingly111Expressions for these quantities up to
8th order and other important expressions pertinent to this work are contained
within a Mathematica notebook posted on Google Drive at
https://drive.google.com/open?id=1Z-V6HgfN_c-yEP-MtHpV8_qtqzIT6ZCQ .:

$\displaystyle W=\omega^{(0)}+\omega^{(1)}+\omega^{(2)}+\dots,$ $\displaystyle
W^{-1}=(W^{-1})^{(0)}+(W^{-1})^{(1)}+(W^{-1})^{(2)}+\dots,$ $\displaystyle
W^{-2}\dot{W}=(W^{-2}\dot{W})^{(0)}+(W^{-2}\dot{W})^{(1)}+(W^{-2}\dot{W})^{(2)}+\dots,$
$\displaystyle
W^{-3}\dot{W}^{2}=(W^{-3}\dot{W}^{2})^{(0)}+(W^{-3}\dot{W}^{2})^{(1)}+(W^{-3}\dot{W}^{2})^{(2)}+\dots~{}.$
(14)

Using the distribution function in Eqs. (4) and (5), the observational
constraints for the present time $t_{0}$, $\rho_{\phi}(t_{0})=\rho_{DE0}$ and
$\frac{p_{\phi}}{\rho_{\phi}}\approx-1$, can specify $m$ and $\xi$. Using
$\rho_{DE0}=1.6077\times 10^{-10}$ eV4 from best-fit values from Planck Planck
and choosing a constant $w_{\phi}=-0.9$ and the corresponding $a(t)$ for the
epoch of DE domination, we obtain $m=1.57\times 10^{-54}$ eV and $\xi=0.176$.
A more precise approach, especially in light of the interaction between DE and
DM, would involve a global fit for $m$ and $\xi$ along with all other
parameters which may or may not assume a constant $w_{\phi}$ , but since we
expect the exchange between DE and DM to be small (and since doing a rigorous
numerical fit to data is not the purpose of this work), we use these values in
our final analysis.

Calculating the Interaction Kernel

In this section, we write the Boltzmann equation in terms of the continuity
equation Eq. (2) in order to identify the interaction kernel $Q$ in terms of
the collision term from Boltzmann’s equation. We then calculate the collision
term for a tree-level two-to-two conversion via the Yukawa-type coupling we
assumed in Eq. (3) in terms of the effective distribution function for dark
energy we calculated in the previous section.

The Liouville operator in FLRW space applied to the DM distribution function
$f_{\psi}$ is Bernstein

$\frac{\partial f_{\psi}}{\partial t}-\frac{\dot{a}}{a}k_{i}\frac{\partial
f_{\psi}}{\partial k_{i}}=C(f_{\psi}),$ (15)

where $C(f_{\psi})$ is the collision term, and $f_{\psi}$ is given by the
Bose-Einstein distribution function since DM is assumed to be a scalar field:

$f_{\psi}(k)=\frac{1}{e^{\sqrt{k^{2}+M^{2}}/T}-1}.$ (16)

Technically, the use of the Bose-Einstein distribution is only valid assuming
there is no interaction, and we could be more precise by using cosmological
perturbation theory to get a first-order correction to the distribution
function for dark matter. However, the collision term we consider will end up
being very small, so using the Bose-Einstein distribution here is valid. Also,
the collision term is for local interactions, so we use the local definitions
of the stress-energy components and the Bose-Einstein distribution function,
according to Bernstein ; KolbTurner . We will still utilize our distribution
function for dark energy derived using the large-scale definitions of the
stress-energy components out of necessity. In theory, one could consider the
calculation of a general collision term from an $S$-matrix in curved space
1908.06717  on large scales, but this is not necessary for our purposes here.

Applying $\int\frac{d^{3}k}{(2\pi)^{3}}$ to both sides of Eq. (15), we obtain

$\int\frac{d^{3}k}{(2\pi)^{3}}\dot{f_{\psi}}-\int\frac{d^{3}k}{(2\pi)^{3}}\frac{\dot{a}}{a}k\frac{\partial
f_{\psi}}{\partial k}=\int\frac{d^{3}k}{(2\pi)^{3}}C(f_{\psi}),$ (17)

where we have used $\frac{d}{dk}=\frac{k^{i}}{k}\frac{d}{dk_{i}}$ where
$k_{i}k^{i}\equiv k^{2}$. Using the number density expression,

$n_{\psi}=\int\frac{d^{3}k}{(2\pi)^{3}}f_{\psi},$ (18)

and integrating the second term on the lefthand side by parts (the boundary
term goes to $0$ since $f_{\psi}\rightarrow 0$ at infinity), we get the
typical form of the Boltzmann equation:

$\dot{n}_{\psi}+3\frac{\dot{a}}{a}n_{\psi}=\int\frac{d^{3}k}{(2\pi)^{3}}C(f_{\psi})\equiv
B,$ (19)

where we have defined $B$ for brevity’s sake. We will consider the interaction
$\psi_{1}\psi_{2}\leftrightarrow\phi_{3}\phi_{4}$ (assuming in Eq. (3) that
$V(\psi)=\frac{1}{2}M^{2}\psi^{2}+$ any higher-order terms which we will
ignore). Including the interaction’s symmetry factor of $1/2$, $B$ is given by
Bernstein ; KolbTurner

$\displaystyle B=$
$\displaystyle\int\frac{d^{3}k_{1}}{(2\pi)^{3}2\sqrt{k_{1}^{2}+M^{2}}}\frac{d^{3}k_{2}}{(2\pi)^{3}2\sqrt{k_{2}^{2}+M^{2}}}\frac{d^{3}k_{3}}{(2\pi)^{3}2\sqrt{k_{3}^{2}+m^{2}}}\frac{d^{3}k_{4}}{(2\pi)^{3}2\sqrt{k_{4}^{2}+m^{2}}}\times$
$\displaystyle(2\pi)^{4}\delta^{4}(k_{1}+k_{2}-k_{3}-k_{4})\left[f_{\psi}(k_{1})f_{\psi}(k_{2})-f_{\phi}(k_{3})f_{\phi}(k_{4})\right]|\mathcal{M}|^{2},$
(20)

which assumes the symmetry of $\mathcal{M}_{12\rightarrow
34}=\mathcal{M}_{34\rightarrow 12}$ and neglects the terms that have three
factors of $f$ since these terms are comparatively small.

Since $\rho=n\langle E\rangle$, we can rewrite Eq. (19) as

$\dot{\rho}_{\psi}+3H\rho_{\psi}=B\langle
E_{\psi}\rangle+n_{\psi}\frac{\partial\langle E_{\psi}\rangle}{\partial t},$
(21)

where

$\langle
E_{\psi}\rangle=\frac{\int\frac{d^{3}k}{(2\pi)^{3}}\sqrt{k^{2}+M^{2}}f_{\psi}(k)}{\int\frac{d^{3}k}{(2\pi)^{3}}f_{\psi}(k)}$
(22)

and

$\rho_{\psi}=\int\frac{d^{3}k}{(2\pi)^{3}}\sqrt{k^{2}+M^{2}}f_{\psi}(k).$ (23)

Since $\dot{\rho}_{\psi}+3H(\rho_{\psi}+p_{\psi})=Q$ from Eq. (2) and
$p_{\psi}\approx 0$ since we are dealing with cold dark matter, we see from
Eq. (21) that

$Q=B\langle E_{\psi}\rangle+n_{\psi}\frac{\partial\langle
E_{\psi}\rangle}{\partial t}.$ (24)

Also, for cold (i.e., non-relativistic) dark matter, $\langle
E_{\psi}\rangle=M$ KolbTurner , so we can express $Q$ as

$Q=MB.$ (25)

Alternatively, we can multiply both sides of Eq. (15) by $\sqrt{k^{2}+M^{2}}$,
apply $\int\frac{d^{3}k}{(2\pi)^{3}}$ to both sides, and integrate by parts
the second term on the lefthand side, and we obtain

$\dot{\rho}_{\psi}+3H(\rho_{\psi}+p_{\psi})=\int\frac{d^{3}k}{(2\pi)^{3}}C(f_{\psi})\sqrt{k^{2}+M^{2}}.$
(26)

From this, we have a general expression for the 2-to-2 interaction $Q$ between
any two perfect fluids. However, Eq. (25) is more convenient to use in our
case, and we use it in what follows.

The Lorentz invariant scattering amplitude for tree-level Yukawa-type
interaction between two scalar fields that appears in $B$ is

$|\mathcal{M}|^{2}=g^{4}\left(\frac{1}{M^{2}-t}+\frac{1}{M^{2}-u}\right)^{2}.$
(27)

where $s=-\eta_{\mu\nu}(k_{1}+k_{2})^{\mu}(k_{1}+k_{2})^{\nu}$,
$t=-\eta_{\mu\nu}(k_{1}-k_{3})^{\mu}(k_{1}-k_{3})^{\nu}$, and
$u=-\eta_{\mu\nu}(k_{1}-k_{4})^{\mu}(k_{1}-k_{4})^{\nu}$.

Assuming the DM mass is bigger than the DE mass ($M>m$), and working in the
center-of-mass frame, we can reduce the $f_{\phi}(k_{3})f_{\phi}(k_{4})$ term
of $B$ in Eq. (Deriving the Dark Matter-Dark Energy Interaction Term in the
Continuity Equation from the Boltzmann Equation) to

$f_{\phi}(k_{3})f_{\phi}(k_{4})~{}{\rm
term}=\frac{1}{8(2\pi)^{5}}\int_{\sqrt{M^{2}-m^{2}}}^{\infty}dk_{3}dk_{4}k_{4}^{3}\frac{k_{4}}{k_{3}^{2}+m^{2}}\sqrt{1-\frac{M^{2}}{k_{3}^{2}+m^{2}}}f_{\phi}(k_{3})f_{\phi}(k_{4})\mathcal{A}(k_{3}),$
(28)

where the bounds of the integral apply for both $k_{3}$ and $k_{4}$,
$\mathcal{A}(k)\equiv\int_{-1}^{1}d(\cos\theta)|\mathcal{M}|^{2}$,
$k_{3}\equiv|\vec{k_{3}}|$, and $k_{4}\equiv|\vec{k_{4}}|$. Similarly, the
$f_{\phi}(k_{1})f_{\phi}(k_{2})$ term of $B$ reduces to

$f_{\psi}(k_{1})f_{\psi}(k_{2})~{}{\rm
term}=\frac{1}{8(2\pi)^{5}}\int_{0}^{\infty}dk_{1}dk_{2}k_{2}^{3}\frac{k_{2}}{k_{1}^{2}+M^{2}}\sqrt{1-\frac{m^{2}}{k_{1}^{2}+M^{2}}}f_{\psi}(k_{1})f_{\psi}(k_{2})\mathcal{A}(k_{1}).$
(29)

$\mathcal{A}(k_{3})$ evaluates to

$\displaystyle\mathcal{A}(k_{3})=\frac{2g^{4}\left(2k_{3}\left(2k_{3}^{2}+m^{2}\right)\sqrt{-k_{3}^{2}-m^{2}+M^{2}}+\left(m^{4}+4k_{3}^{2}M^{2}\right)\tan^{-1}\left[\frac{2k_{3}\sqrt{-k_{3}^{2}-m^{2}+M^{2}}}{2k_{3}^{2}+m^{2}}\right]\right)}{k_{3}\left(2k_{3}^{2}+m^{2}\right)\sqrt{-k_{3}^{2}-m^{2}+M^{2}}\left(m^{4}+4k_{3}^{2}M^{2}\right)},$
(30)

and $\mathcal{A}(k_{1})$ is

$\displaystyle\mathcal{A}(k_{1})=2g^{4}\left(\frac{2}{4k_{1}^{2}M^{2}+\left(m^{2}-2M^{2}\right)^{2}}+\frac{\tan^{-1}\left[\frac{2k_{1}\sqrt{-k_{1}^{2}+m^{2}-M^{2}}}{2k_{1}^{2}-m^{2}+2M^{2}}\right]}{\sqrt{-k_{1}^{2}+m^{2}-M^{2}}\left(2k_{1}^{3}-k_{1}m^{2}+2k_{1}M^{2}\right)}\right).$
(31)

Putting all of this together, we can numerically calculate a finite output for
$Q$ as a function of $a(t)$, $m$, $M$, $\xi$, and $g$. We plot in Fig. (1) $Q$
as a function of $a$ for the epoch of DE domination, using a typical value for
DM mass, $M=1$ GeV. The magnitude of $Q$ is very small, indicating a very weak
interaction, as expected. The $f_{\phi}(k_{3})f_{\phi}(k_{4})$ term
contributes essentially $0$ to $Q$, and this is not unexpected given the very
low mass value for DE that we use based on the observational constraints
discussed above. For virtually any value of $m$, $-1\lesssim\xi\lesssim 1$,
$1$ GeV $\lesssim M\lesssim 1$ TeV, and $g\lesssim 1$eV, the value of $Q$ is
$10^{-100}$ eV${}^{5}\lesssim Q\lesssim 10^{-80}$ eV5 or smaller222A phantom
value of $w_{\phi}<-1$ is still observationally plausible even for a normal-
sign kinetic term for the scalar field theory 1804.02987 , which is what we
use in this paper. We checked for values of $w<-1$, such as $w=-2$, and the
range for $Q$ was the same..

Figure 1: We plot $Q$ (in eV5) vs $a$ for $m=1.57\times 10^{-54}$eV, $M=1$
GeV, $\xi=0.176$, $g=0.1$eV. Values of $m$ and $\xi$ chosen by satisfying
$\rho_{\phi}(t_{0})=\rho_{DE0}$ and $\frac{p_{\phi}}{\rho_{\phi}}=-0.9$ (as
constrained by Planck data and discussed above). The plot is virtually
independent of $\xi$ and $m$ (as discussed in the text).

For comparison, for a typical $Q$ parametrization, $Q=\beta H\rho_{DM}$ from
1812.06854  has a good fit to data for $|Q|\approx-10^{-45}$ eV5 evaluated at
the present-day time. So one can see that $Q$ for the tree-level 2-to-2
conversion via the Yukawa-type coupling is very small in comparison. More
accurately, we could apply weak equivalence principle constraints to severely
limit the magnitude of our coupling $g$ 1605.00996  and thus our magnitude of
$Q$, but since the weak equivalence principle is not used to constrain the
comparison value of $Q$ from 1812.06854 , we will ignore this for the sake of
comparison.

Conclusion

Many seemingly ad hoc parametrizations are used in the literature to model the
interaction between dark energy and dark matter since we do not know the
interaction between them on a fundamental level. In this work, we have used
adiabatic subtraction to obtain a finite analytic expression for the
distribution function for dark energy that is compatible with a negative fluid
pressure. Using this distribution function in the Boltzmann equation, we
analytically calculate the interaction $Q$ due to the tree-level $2$-to-$2$
conversion via our Yukawa-type coupling between dark energy and dark matter.
The interaction is very small compared to typical interaction found in the
literature.

More accurate treatment of our methodology is certainly possible (i.e., using
cosmological perturbation theory, doing a global fit to the matter power
spectrum and other data), and one can apply this method to other field
theories and couplings for DE and DM. This work serves as a proof of concept
that one can treat DE using the formalism of the Boltmann equation with a
collision term.

## References

  * (1) Ivan de Martino, Symmetry 10, 372 (2018). arXiv: 1809.00550 [astro-ph.CO].
  * (2) Azam Hussain, J. Astrophys. Astron. 40, 5, 43 (2019). arXiv:1809.02411 [astro-ph.CO].
  * (3) Supriya Pan, John D. Barrow, and Andronikos Paliathanasis, Eur. Phys. J. C 79, 2, 115 (2019). arXiv:1812.05493 [gr-qc].
  * (4) Weiqiang Yang, Narayan Banerjee, Andronikos Paliathanasis, and Supriya Pan, Phys. Dark Univ. 26, 100383 (2019). arXiv:1812.06854 [astro-ph.CO].
  * (5) Matteo Cataneo et al, Mon. Not. Roy. Astron. Soc. 488, 2, 2121 (2019). arXiv:1812.05594 [astro-ph.CO].
  * (6) Lloyd Knox and Marius Millea, Phys. Rev. D 101, 4, 043533 (2020). arXiv:1908.03663 [astro-ph.CO].
  * (7) Yang-Jie Yan, Deng Wang, and Xin-He Meng. arXiv:1801.00689 [astro-ph.CO].
  * (8) Luca Amendola, Phys. Rev. D 62, 043511 (2000). arXiv:astro-ph/9908023.
  * (9) Yu. L. Bolotin, A. Kostenko, O.A. Lemets, and D.A. Yerokhin, Int. J. Mod. Phys. D 24 no.03, 1530007 (2014). arXiv:1310.0085 [astro-ph.CO].
  * (10) B. Wang, E. Abdalla, F. Atrio-Barandela, and D. Pavon, Rept. Prog. Phys. 79 no.9, 096901 (2016). arXiv:1603.08299 [astro-ph.CO].
  * (11) Weiqiang Yang, Supriya Pan, and John D. Barrow, Phys. Rev. D 97, 043529 (2018). arXiv:1706.04953 [astro-ph.CO].
  * (12) Weiqiang Yang, Ankan Mukherjee, Eleonora Di Valentino, and Supriya Pan, Phys. Rev. D 98, 123527 (2018). arXiv:1809.06883 [astro-ph.CO].
  * (13) Weiqiang Yang et al., JCAP 09, 019 (2018). arXiv:1805.08252 [astro-ph.CO].
  * (14) Weiqiang Yang et al., JCAP 07, 037 (2019). arXiv:1905.08252 [astro-ph.CO].
  * (15) Eleonora Di Valentino, Alessandro Melchiorri, Olga Mena, and Sunny Vagnozzi. arXiv:1908.04281 [astro-ph.CO].
  * (16) Christian G. Boehmer, Nicola Tamanini, and Matthew Wright, Phys. Rev. D 91, 12, 123002 (2015). arXiv:1501.06540 [gr-qc].
  * (17) Christian G. Boehmer, Nicola Tamanini, and Matthew Wright, Phys. Rev. D 91, 12, 123003 (2015). arXiv:1502.04030 [gr-qc].
  * (18) C. van de Bruck and J. Morrice, JCAP 04 036 (2015). arXiv:1501.03073 [gr-qc].
  * (19) Linfeng Xiao, Rui An, Le Zhang, Bin Yue, Yidong Xu et al., Phys. Rev. D 99, 2, 023528 (2019). arXiv:1807.05541 [astro-ph.CO].
  * (20) Stephon Alexander, Marina Cortês, Andrew R. Liddle, João Magueijo, Robert Sims et al., Phys. Rev. D 100, 8, 083507 (2019). arXiv:1905.10382 [gr-qc].
  * (21) Sergio del Campo, Ramon Herrera, German Olivares, and Diego Pavón, Phys. Rev. D 74, 023501 (2006). arXiv:astro-ph/0606520.
  * (22) Sergio del Campo, Ramón Herrera, and Diego Pavón, Phys. Rev. D 91, 123539 (2015). arXiv:1507.00187 [gr-qc].
  * (23) Hao Wei and Shuang Nan Zhang, Phys. Lett. B 644, 7 (2007). arXiv:astro-ph/0609597.
  * (24) Luca Amendola, Gabriela Camargo Campos, and Rogerio Rosenfeld, Phys. Rev. D 75, 083506 (2007). arXiv:astro-ph/0610806.
  * (25) Zong-Kuan Guo, Nobuyoshi Ohta, and Shinji Tsujikawa, Phys. Rev. D 76, 023508 (2007). arXiv:astro-ph/0702015.
  * (26) Jian-Hua He and Bin Wang, JCAP 06 010 (2008). arXiv:0801.4233 [astro-ph].
  * (27) Christof Wetterich, Astron. Astrophys. 301, 321 (1995). arXiv:hep-th/9408025.
  * (28) Federico Marulli, Marco Baldi, and Lauro Moscardini, MNRAS 420, 2377 (2012). arXiv:1110.3045 [astro-ph.CO].
  * (29) Luca Amendola, Marco Baldi, and Christof Wetterich, Phys. Rev. D 78, 023015 (2008). arXiv:0706.3064 [astro-ph].
  * (30) Joseph P. Johnson and S. Shankaranarayanan. arXiv:2006.04618 [gr-qc].
  * (31) Guido D’Amico, Teresa Hamill, and Nemanja Kaloper, Phys. Rev. D 94, 10, 103526 (2016). arXiv:1605.00996 [hep-th].
  * (32) Ryotaro Kase and Shinji Tsujikawa, Phys. Rev. D 101, 6, 063511 (2020). arXiv:1910.02699 [gr-qc].
  * (33) Sandro M.R. Micheletti, Int. J. Mod. Phys. D 29, 2050057 (2020). arXiv:1808.05015 [gr-qc].
  * (34) André A. Costa, Lucas C. Olivari, and E. Abdalla, Phys. Rev. D 92, 103501 (2015). arXiv:1411.3660 [astro-ph.CO].
  * (35) Chung-Pei Ma and Edmund Bertschinger, Astrophys. J., 455 7 (1995). arXiv:astro-ph/9506072.
  * (36) Saurya Das and Rajat K. Bhaduri, Class. Quant. Grav. 32, 10, 105003 (2015). arXiv:1411.0753 [gr-qc].
  * (37) Valerio Faraoni, Phys. Rev. D 85, 024040 (2012). arXiv:1201.1448 [gr-qc].
  * (38) Kevin J. Ludwick, Phys. Rev. D 98 043519 (2018). arXiv:1804.02987 [physics.gen-ph].
  * (39) Leonard E. Parker and David T. Toms, Quantum Field Theory in Curved Spacetime: Quantized Field and Gravity, Cambridge University Press (2009).
  * (40) Antonio Ferreiro, Jose Navarro-Salas, and Silvia Pla, Phys. Rev. D 98 045015 (2018). arXiv:1807.10361 [gr-qc].
  * (41) Yohei Ema, Kazunori Nakayama, and Yong Tang, JHEP 09 135 (2018). arXiv:1804.07471 [hep-ph].
  * (42) P.A.R. Ade et al. (Planck Collaboration), Astron. Astrophys. 594, A13 (2016). arXiv:1502.01589 [astro-ph.CO].
  * (43) Jeremy Bernstein, Kinetic theory in the expanding universe, Cambridge University Press (1988).
  * (44) Edward W. Kolb and Michael S. Turner, The Early Universe, Westview Press (1990).
  * (45) Susobhan Mandal and Subhashish Banerjee. arXiv:1908.06717 [hep-ph].


# S++: A Fast and Deployable Secure-Computation Framework for Privacy-
Preserving Neural Network Training

Prashanthi Ramachandran 1 Shivam Agarwal 1 Arup Mondal 1
Aastha Shah 1 Debayan Gupta 1

###### Abstract

We introduce S++, a simple, robust, and deployable framework for training a
neural network (NN) using private data from multiple sources, using secret-
shared secure function evaluation. In short, consider a virtual third party to
whom every data-holder sends their inputs, and which computes the neural
network: in our case, this virtual third party is actually a set of servers
which individually learn nothing, even with a malicious (but non-colluding)
adversary.

Previous work in this area has been limited to just one specific activation
function: ReLU, rendering the approach impractical for many use-cases. For the
first time, we provide fast and verifiable protocols for all common activation
functions and optimize them for running in a secret-shared manner. The ability
to quickly, verifiably, and robustly compute exponentiation, softmax, sigmoid,
etc., allows us to use previously written NNs without modification, vastly
reducing developer effort and complexity of code. In recent times, ReLU has
been found to converge much faster and be more computationally efficient as
compared to non-linear functions like sigmoid or tanh. However, we argue that
it would be remiss not to extend the mechanism to non-linear functions such as
the logistic sigmoid, tanh, and softmax that are fundamental due to their
ability to express outputs as probabilities and their universal approximation
property. Their contribution in RNNs and a few recent advancements also makes
them more relevant.

## Introduction

Neural networks (NN) are used in areas ranging from image classification to
machine translation, and there are often multiple parties that contribute
possibly private data to the training process. However, training data
interaction and the resulting model may still leak a significant amount of
private data. Thus, there arises a need for secure neural networks, which can
help parties collaboratively train a neural network without revealing their
private input data.

Several approaches have been proposed to overcome possible data leakage, based
on secure multi-party computation (MPC). MPC is a powerful way to protect
sensitive training data which uses a range of cryptographic primitives to
allow multiple parties to compute a function without revealing the inputs of
any individual party (beyond what is implied by the output of the function
itself). However, previously proposed MPC-based schemes (Ohrimenko et al.
2016; Hunt et al. 2018; Mohassel and Zhang 2017; Wagh, Gupta, and Chandran
2018; Patra and Suresh 2020) do not support models that use exponentiation-
based activation functions such as logistic sigmoid, softmax, and tanh. These
protocols are largely restricted to ReLU and its variants, which might
restrict applicability in some specific cases (as described in the motivation
section).

In this paper, we propose S++, a three-party secure computation framework for
secure exponentiation, and exponentiation-based activation functions such as
logistic sigmoid, tanh, and softmax. This enables us to construct three-party
secure protocols for training and inference of several NN architectures such
that no single party learns any information about the data. S++ is an
efficient MPC-based privacy-preserving neural network training framework based
on (Wagh, Gupta, and Chandran 2018) for 3PC with the activation functions
logistic sigmoid, tanh, their derivatives, and softmax. In our setting, there
are $D$ (where $D$ can be arbitrary) data owners who wish to jointly train a
model over their data with the help of $3$-servers. In the setup phase, these
data owners create ”additive secret-shares” of their input data and send one
share each to $2$-servers. In the computation phase, the $3$-servers (the $2$
primary servers and the helper server) collectively run an interactive
protocol to train a neural network on the data owners’ data without learning
any information beyond the trained model.

### Contributions

We summarize our key contributions as follows:

  * •

We first propose a secure protocol for exponentiation in the three-party
setting. This protocol is based on SCALE MAMBA’s (Aly et al. 2020) protocol
for base 2 exponentiation. It also uses primitives from (Catrina and Saxena
2010). We modify these ideas for our three-party setting and additively
secret-shared data.

  * •

We describe novel secure protocols for the logistic sigmoid, softmax, and
tanh—popular activation functions that are significantly more complex than
ReLU (described by (Wagh, Gupta, and Chandran 2018))—along with their
derivatives with the help of the above-mentioned exponentiation protocol. The
inclusion of these protocols in a secure and private setting vastly increases
the practicality of the framework and enables people to convert their
protocols into secure ones without having to redesign the actual internal
structure of their NNs.

### Organization of the paper

The remainder of this paper is organized as follows: we first go over our
motivation for extending secure protocols to functions such as the logistic
sigmoid, tanh and softmax. Then, we look at recent work in this area, largely
in the realm of MPC. After that, we describe the notations used in our
protocols and explain some cryptographic primitives. In the Protocols section,
we delve into our architecture and the supporting protocols, specifically the
supporting 3-server protocols that serve as building blocks for our main
protocols. We describe our main protocols: exponentiation, logistic sigmoid,
tanh, their derivatives, and softmax. After that, we describe the theoretical
efficiency of our protocols and provide an evaluation of our experiments.
Lastly, we describe the future plans and directions for this work.

### Motivation

#### Logistic sigmoid and variants

The logistic sigmoid is a common activation function to introduce non-
linearity, where sigmoid$(x)=\displaystyle\frac{1}{1+e^{-x}}$. It essentially
squashes a real value in [0, 1] and has been used widely over the years for
its simplicity, especially in binary classification tasks with a maximum
likelihood loss function. It also has important variants like the tanh and
softmax:

  * •

Tanh as a rescaled logistic sigmoid: The tanh is a rescaling and stretching of
the logistic sigmoid that maps to [-1, 1] i.e., tanh$(x)$ =
2$\times$sigmoid$(2x)-1$. Extending the protocol to the tanh facilitates
applications in recurrent neural networks where it is widely used.

  * •

Logistic sigmoid as a special case of the softmax: The softmax is often used
for classification tasks, where
softmax$(z_{i})=\displaystyle\frac{e^{z_{i}}}{\sum_{i=1}^{k}e^{z_{i}}}$. The
logistic sigmoid happens to be a special case of this function, often used for
binary classification. Enabling softmax thus also extends the framework to
multi-class classification problems.

#### Why is ReLU not enough?

Despite ReLU being more computationally efficient than logistic sigmoid or
tanh (Krizhevsky, Sutskever, and Hinton 2017), we argue that it is worth
extending the secure computation protocol to popular functions particularly
involving exponentiation such as logistic sigmoid, tanh, and softmax for the
following reasons (with more details in the appendix):

  * •

For classification tasks, output layers usually use softmax or logistic
sigmoid because they make for more interpretable values in [0, 1] as opposed
to ReLU and variations that are unbounded. In fact, though previously missing
in the literature, (Asadi and Jiang 2020) provide theoretical proofs on the
universal approximation of softmax used in the output layer for multiclass
classification tasks as well.

  * •

ReLU’s unbounded nature limits its applicability in RNNs. Capturing long-term
dependencies becomes harder with an unbounded function due to exploding
gradients, as explored in (Pascanu, Mikolov, and Bengio 2013). Because of the
temporal correlations captured by RNNs (by forming connections between hidden
layers), long-term components add up and the norm of the gradient of the
hidden layers explode and ReLU’s unbounded nature exacerbates this problem.
There have been a few notable improvements to assuage this, but even in LSTMs,
the tanh often gives more consistent results (Kent and Salem 2019) and ReLU
tends to require greater fine-tuning.

  * •

’Leaky’ ReLU has been used to overcome ReLU’s vanishing gradient problem. (Xu,
Huang, and Li 2016) have also revised logistic sigmoid to a scaled version of
it and tanh to a ‘leaky’ tanh. The ‘leaky’ version penalizes the negative part
of the domain. They found it to be much faster than the tanh, it gave better
results than ReLU, and is almost identical to leaky ReLU. Such potential
improvements highlight the need for more secure protocols for such functions.

Though ReLU has become widely adapted in the past few years, there are still
some tasks wherein sigmoids are preferred over ReLU, and it is worthwhile to
extend such secure protocols for these as well.

### Recent Work

There have been many privacy-preserving machine learning (PPML) frameworks for
various situations, such as Decision Trees (Agrawal and Srikant 2000; Lindell
and Pinkas 2000), Linear Regression (Du and Atallah 2001; Sanil et al. 2004),
k-means clustering (Jagannathan and Wright 2005; Bunn and Ostrovsky 2007), SVM
classification (Yu, Vaidya, and Jiang 2006; Vaidya, Yu, and Jiang 2008), and
logistic regression (Slavkovic, Nardi, and Tibbits 2007). However, these
cannot generalize to a number of (complex and high accuracy) standard machine
learning applications.

To overcome this, SecureML (Mohassel and Zhang 2017) provided secure protocols
for linear regression, logistic regression, and neural networks with linear
activations, using a combination of arithmetic and garbled circuits for a
single semi-honest adversary in both the 2 and 3-server models. It also
provided secure protocols for logistic sigmoid and softmax in the 2-server
setting (although very briefly) and a method for truncating decimal numbers.
MiniONN (Liu et al. 2017) optimized the protocols of SecureML by reducing the
offline cost of matrix multiplications and increasing the online cost for the
2-server model in the semi-honest adversary setting. Chameleon (Riazi et al.
2018) and Gazelle (Juvekar, Vaikuntanathan, and Chandrakasan 2018) provided
secure inference protocols which are computationally secure against one semi-
honest adversary in the 3-server and 2-server models. Chameleon removes
expensive oblivious transfer protocols by using the third party as a dealer,
while Gazelle focuses on making the linear layers more communication efficient
by providing specialized packing schemes for additively homomorphic encryption
schemes.

SecureNN (Wagh, Gupta, and Chandran 2018) provides secure neural network
training protocols for non-linear activation functions in the 3-party setting
with up to one malicious corruption and shows that the honest-majority can
improve the performance by several orders of magnitude. For linear regression,
logistic regression and neural networks, the problem is even more challenging
as the training procedure computes many instances of non-linear activation
functions such as logistic sigmoid, tanh, and softmax that are expensive to
compute inside a 2 and 3-server.

SecureNN provides efficient protocols for computing the ReLU activation,
maxpool and their derivatives. It suggests that the general idea can extended
to other variants and piecewise linear approximations. However, (Mohassel and
Zhang 2017) show experimentally that low-degree polynomial approximations,
like piecewise linear approximations, are ineffective in terms of accuracy. It
can be proven that hard sigmoid, a piecewise linear approximation of the
logistic sigmoid, performs worse than logistic sigmoid, especially in the case
of linear regression (Maksutov 2018). SecureML also provides higher-order
approximations of logistic sigmoid and softmax using a combination of ReLU
functions in place of exponentiation. However, they do not provide a detailed
algorithm and find the running time in the offline phase to be high.

To overcome this research gap, we propose S++, based on (Wagh, Gupta, and
Chandran 2018), for the 3-server setting with efficient protocols for logistic
sigmoid and tanh and their derivatives, as well as the softmax. This can be
extended to the derivative of the softmax as well.

## Cryptographic Primitives

### Notation

In our work, we use 2-out-of-2 additive sharing over the even ring $Z_{L}$
where $L=2^{l}$. For our purposes, we use the same ring as (Wagh, Gupta, and
Chandran 2018), i.e., $Z_{2^{64}}$. The 2-out-of-2 additive shares are
represented by $\langle x\rangle_{0}^{L}$ and $\langle x\rangle_{1}^{L}$.
where L represents the ring $Z_{L}$ over which the value $x$ has been shared.

We use the notation $\langle a\rangle_{j}$ to denote Party $P_{j}$’s additive
secret share of the value a, and $\lfloor x\rfloor$ to denote the floor of the
value x. The notation $(p_{0},p_{1},\cdots,p_{n-1})$ is used to denote the
bits of a n-bit value, $p$.

Further, we use the following primitives from (Wagh, Gupta, and Chandran 2018)
in our work:

  1. 1.

Matrix Multiplication of two secret shared matrices represented by
$\mathcal{F}_{MatMul}(\\{P_{0},P_{1}\\}P_{2})$.

  2. 2.

Division of a two secret shared values represented by
$\mathcal{F}_{Division}(\\{P_{0},P_{1}\\}P_{2})$.

These protocols are described in the following section.

## Protocols of S++

### Architecture

S++ works with the $3$-server setting similar to SecureNN (Wagh, Gupta, and
Chandran 2018). In this setting, there are two primary servers and one helper
server. The two primary servers hold additive secret shares of the data. Let
$P_{0},P_{1},P_{2}$ be the three servers and $C$ be a set of $n$ participants
$\\{C_{1},C_{2},\dots C_{n}\\}$, where the $i^{th}$ party $C_{i}$ holds its
own private dataset, $\mathbb{D}_{i}$. $\mathcal{M}_{nn}$ is the neural
network model to be trained on the parties’ private data. The participants
$\\{C_{1},C_{2},\dots C_{n}\\}$ split their data into 2-out-of-2 additive
secret shares and give one share to each of the two servers $P_{0}$ and
$P_{1}$. The third server, $P_{2}$, similar to (Wagh, Gupta, and Chandran
2018), is crucial to the protocols, but does not hold any data. We use the
same integer representation for fixed point numbers as (Wagh, Gupta, and
Chandran 2018). However, to tackle any overflow that can render the
exponentiation protocol impractical for use, we also consider extending the
integer ring to $Z_{128}$. S++ provides security against one semi-honest
adversary and up to one malicious corruption.

Figure 1: Function Dependencies in S++.

### Supporting Protocols

#### Matrix Multiplication ($\mathcal{F}_{MatMul}(\\{P_{0},P_{1}\\}P_{2})$):

In S++, this function from (Wagh, Gupta, and Chandran 2018) has been used for
multiplying fixed point values (represented as unsigned integers in the
$Z_{L}$ ring) by considering the values as $1\times 1$ matrices.

Two parties, $P_{0}$ and $P_{1}$ are required to hold shares of $X\in
Z_{L}^{m\times n}$ and $Y\in Z_{L}^{n\times v}$. After running this
interactive secure protocol, $P_{0}$ gets $\langle X.Y\rangle^{L}_{0}$ and
$P_{1}$ gets $\langle X.Y\rangle^{L}_{1}$.

#### Division ($\mathcal{F}_{Division}(\\{P_{0},P_{1}\\}P_{2})$):

This function, also from (Wagh, Gupta, and Chandran 2018), is used to perform
secure division on secret shared values. In this protocol, two parties,
$P_{0}$ and $P_{1}$ are required to hold shares of values $x_{j}$ and $y_{j}$
respectively, for $j\text{in}\\{0,1\\}$. After running the protocol, they
obtain $\langle x/y\rangle^{L}_{0}$ and $\langle x/y\rangle^{L}_{1}$
respectively.

#### Taylor Series Expansion
($\mathcal{F}_{taylorExp}(\\{P_{0},P_{1}\\}P_{2})$; see algorithm 1):

In this protocol, we compute $e^{x}$, where $x$ is a secret fractional value.
We do so by computing the value of the Taylor expansion up to four terms. In
the end of the protocol, parties $P_{0}$ and $P_{1}$ obtain shares of $e^{x}$.

Algorithm 1 Taylor Expansion,
$\mathcal{F}_{taylorExp}(\\{P_{0},P_{1}\\},P_{2})$

Input: $P_{0},P_{1}$ hold $\\{\langle x\rangle_{0}^{L}\\}$ and $\\{\langle
x\rangle_{1}^{L}\\}$ respectively where $x<1$.

Output: $P_{0}$ and $P_{1}$ obtain $\langle y\rangle_{j}^{L}$ = $\langle
e^{x}\rangle_{j}^{L}$.

Common Randomness: $P_{0}$ and $P_{1}$ hold random shares of zero - $u_{0}$
and $u_{1}$ respectively.

1: Parties $P_{0}$ and $P_{1}$ compute $\langle c\rangle_{j}^{L}=j$.

2: Parties $P_{0}$ and $P_{1}$ compute $\langle
numerator\rangle_{j}^{L}=\langle x\rangle_{j}^{L}$ and $denominator=1$.

3: Now, Parties $P_{0}$ and $P_{1}$ compute $\langle c\rangle_{j}^{L}=\langle
c\rangle_{j}^{L}+\langle numerator\rangle_{j}^{L}$.

4: for $i=2,\ldots,4$ do

5: Parties $P_{0}$, $P_{1}$ and $P_{2}$ invoke
$\mathcal{F}_{MatMul}(\\{P_{0},P_{1}\\}P_{2})$ with inputs $\langle
numerator\rangle_{j}^{L}$ and $\langle x\rangle_{j}^{L}$ to obtain $\langle
numerator\rangle_{j}^{L}$ for $j\in\\{0,1\\}$.

6: $denominator=denomiator\times i$

7: Parties $P_{0}$ and $P_{1}$ compute $\langle c\rangle_{j}^{L}=\langle
c\rangle_{j}^{L}+\frac{\langle numerator\rangle_{j}^{L}}{denominator}$

8: end for

9: $P_{j}$ for $j\in\\{0,1\\}$ output $\langle c\rangle_{j}^{L}$ \+ $u_{j}$.

#### Exponentiation ($\mathcal{F}_{Exp}(\\{P_{0},P_{1}\\},P_{2})$); see
algorithm 2:

In this protocol, we compute $e^{x}$, where $x$ is secret. This protocol is
based on SCALE MAMBA’s (Aly et al. 2020) base-2 exponentiation protocol. It
uses primitives like $\mathcal{F}_{Trunc}$, $\mathcal{F}_{BitDecomp}$,
$\mathcal{F}_{PreMult}$ from (Catrina and Saxena 2010).

While (Aly et al. 2020) uses the polynomial $P_{1045}(X)$ from (Hart 1978) to
compute the exponentiation for the fractional part of a number, we use a
secure protocol to compute the output using the Taylor series expansion of
$e^{x}$. We introduce function $\mathcal{F}_{taylorExp}$ defined above for
this purpose.

This protocol first uses $\mathcal{F}_{Trunc}$ to split the input $a$ into its
integer ($a_{\text{int}}$) and fractional ($a_{\text{frac}}$) parts. After
that, it evaluates $e^{a_{\text{int}}}$ using primitives
$\mathcal{F}_{BitDecomp}$ and $\mathcal{F}_{PreMult}$ described above. It
further evaluates $e^{a_{\text{frac}}}$ using $\mathcal{F}_{taylorExp}$. It
then multiplies these two values to obtain
$e^{a_{\text{int}}+a_{\text{frac}}}=e^{a}$.

Algorithm 2 Exponentiation $\mathcal{F}_{Exp}(\\{P_{0},P_{1}\\},P_{2})$

Input: $P_{0}$ and $P_{1}$ hold $\langle x\rangle_{0}^{L}$ and $\langle
x\rangle_{1}^{L}$ (shares of a value x).

Output: $P_{0}$ and $P_{1}$ obtain $\langle y\rangle_{j}^{L}$ = $\langle
e^{x}\rangle_{j}^{L}$.

1: For $j\in\\{0,1\\}$, party $P_{j}$ calls
$\mathcal{F}_{Trunc}(\\{P_{0},P_{1}\\})$ to obtain $\langle a\rangle_{j}^{L}$,
which is the $j^{\text{th}}$ share of $\lfloor x\rfloor$.

2: For $j\in\\{0,1\\}$, party $P_{j}$ gets the fractional part of x by locally
computing:  $\langle b\rangle_{j}^{L}=\langle x\rangle_{j}^{L}-\langle
a\rangle_{j}^{L}$

3: For $j\in\\{0,1\\}$ party $P_{j}$ invokes
$\mathcal{F}_{BitDecomp}(\\{P_{0},P_{1}\\})$ to obtain $(\langle
c_{0}\rangle_{j}^{L},\langle c_{1}\rangle_{j}^{L},\cdots,\langle
c_{m-1}\rangle_{j}^{L})$: shares of $(x_{0},x_{1},\cdots,x_{m-1})$, where $x$
is a m-bit number.

4: for $i=0,1,\ldots,m-1$ do

5: $P_{j}$ for $j\in\\{0,1\\}$ computes $\langle
v_{i}\rangle_{j}^{L}=e^{2^{i}}.(\langle c_{i}\rangle_{j}^{L})+j-(\langle
c_{i}\rangle_{j}^{L})$.

6: end for

7: $P_{j}$ for $j\in\\{0,1\\}$ invokes
$\mathcal{F}_{PreMult}(\\{P_{0},P_{1}\\})$ to get $\langle m\rangle_{j}^{L}$.

8: $P_{j}$ for $j\in\\{0,1\\}$ invokes
$\mathcal{F}_{Taylor}(\\{P_{0},P_{1}\\})$ to get $\langle n\rangle_{j}^{L}$

9: Finally, $P_{j}$ for $j\in\\{0,1\\}$ invokes
$\mathcal{F}_{MatMul}(\\{P_{0},P_{1}\\}P_{2})$ with inputs $\langle
m\rangle_{j}^{L}$ and $\langle n\rangle_{j}^{L}$ to obtain share $\langle
m\rangle_{j}^{L}$ of:  $y=m\times n$.

### Main Protocols

Our main protocols include logistic sigmoid ($\mathcal{F}_{Sigmoid}$;
described in algorithm 3), tanh ($\mathcal{F}_{Tanh}$; described in algorithm
5), and their derivatives; see algorithms 4 and 6, as well as the softmax
function ($\mathcal{F}_{Softmax}$; described in algorithm 7).

#### Sigmoid, $\mathcal{F}_{Sigmoid}(\\{P_{0},P_{1}\\},P_{2})$; see algorithm
3):

This protocol can be used to securely compute logistic sigmoid. Parties
$P_{0}$, $P_{1}$ and $P_{2}$ invoke $\mathcal{F}_{Exp}$ with shares of the
input $x$ to obtain output shares of $e^{x}$. They then jointly compute shares
of $1+e^{x}$. They use (Wagh, Gupta, and Chandran 2018)’s secure division
protocol to obtain $\sigma(x)=\frac{e^{x}}{1+e^{x}}$.

Algorithm 3 Sigmoid, $\mathcal{F}_{Sigmoid}(\\{P_{0},P_{1}\\},P_{2})$

Input: $P_{0},P_{1}$ hold $\\{\langle x\rangle_{0}^{L}\\}$ and $\\{\langle
x\rangle_{1}^{L}\\}$ respectively.

Output: $P_{0},P_{1}$ get $\\{\langle\sigma(x)\rangle_{0}^{L}\\}$ and
$\\{\langle\sigma(x)\rangle_{1}^{L}\\}$ respectively where
$\sigma(x)=\frac{e^{x}}{1+e^{x}}$.

Common Randomness: $P_{0}$ and $P_{1}$ hold random shares of zero - $u_{0}$
and $u_{1}$ respectively.

1: Parties $P_{0}$, $P_{1}$ and $P_{2}$ invoke
$F_{Exp}(\\{P_{0},P_{1}\\},P_{2})$ with inputs $\langle x\rangle_{0}^{L}$ and
$\langle x\rangle_{1}^{L}$ and obtain $\langle a\rangle_{j}^{L}=\langle
e^{x}\rangle_{j}^{L}$.

2: Now, parties $P_{0}$ and $P_{1}$ compute $\langle b\rangle_{j}^{L}=\langle
a\rangle_{j}^{L}+j$.

3: $P_{0}$, $P_{1}$ and $P_{2}$ invoke
$\mathcal{F}_{Division}(\\{P_{0},P_{1}\\},P_{2})$ with inputs $\langle
a\rangle_{j}^{L}$ and $\langle b\rangle_{j}^{L}$ to obtain $\langle
c\rangle_{j}^{L}$ for $j\in\\{0,1\\}$.

4: $P_{j}$ for $j\in\\{0,1\\}$ output $\langle c\rangle_{j}^{L}$ \+ $u_{j}$.

#### Derivative of Sigmoid,
$\mathcal{F^{D}}_{Sigmoid}(\\{P_{0},P_{1}\\},P_{2})$; see algorithm 4):

This protocol can be used to securely compute the derivative of the logistic
sigmoid. We use the idea that
$\frac{dS(x)}{d(x)}=\sigma^{\prime}(x)=\sigma(x)\cdot(1-\sigma(x)).$ In the
end of the protocol, parties $P_{0}$ and $P_{1}$ obtain the shares of
$\sigma^{\prime}(x)$.

Algorithm 4 Derivative of Sigmoid
$\langle{\mathcal{F^{D}}_{sigmoid}(x)}\rangle^{L}$

Input: $P_{0}$, $P_{1}$ have inputs $\langle{x}\rangle_{j}^{L}$ for $P_{j}$
such that $j\in\\{0,1\\}$.

Output: $P_{0}$, $P_{1}$ get $\langle{\sigma^{\prime}(c)}\rangle_{0}^{L}$ and
$\langle{\sigma^{\prime}(c)}\rangle_{1}^{L}$ where
$\sigma^{\prime}(c)=\sigma(x)(1-\sigma(x))$

Common Randomness: $P_{0}$ and $P_{1}$ hold random shares of zero - $u_{0}$
and $u_{1}$ respectively.

1: $P_{0}$, $P_{1}$ and $P_{2}$ invoke
$\mathcal{F}_{Sigmoid}{(\\{P_{0},P_{1}\\},P_{2})}$ with $P_{j}$ for
$j\in\\{0,1\\}$ having input $\langle{x}\rangle_{j}^{L}$ to learn
$\langle{a}\rangle_{j}^{L}=\langle{\sigma(x)}\rangle_{j}^{L}$.

2: $P_{j}$ computes $\langle{a}\rangle_{j}^{L}$ =
$-1\times\langle{x}\rangle_{j}^{L}+j$ for $j\in\\{0,1\\}$

3: $P_{j}$ computes $\langle{b}\rangle_{j}^{L}$ = $\langle{a}\rangle_{j}^{L}$
for $j\in\\{0,1\\}$

4: $P_{0}$, $P_{1}$ and $P_{2}$ invoke
$\mathcal{F}_{MatMul}{(\\{P_{0},P_{1}\\},P_{2})}$ with $P_{j}$ for
$j\in\\{0,1\\}$ having inputs $\langle{a}\rangle_{j}^{L}$ and
$\langle{b}\rangle_{j}^{L}$ to learn
$\langle{c}\rangle_{j}^{L}=\langle{a}\rangle_{j}^{L}\times\langle{b}\rangle_{j}^{L}$.

5: $P_{j}$ for $j\in\\{0,1\\}$ output:
$\langle{c}\rangle_{j}^{L}+\langle{u}\rangle_{j}^{L}$

#### Tanh, $\mathcal{F}_{Tanh}(\\{P_{0},P_{1}\\},P_{2})$; see algorithm 5):

This protocol can be used to securely compute the tanh function. For this
protocol, we use the idea that $tanh(x)$ is related to $\sigma(x)$ as
$tanh(x)=2\sigma(2x)-1$.

Algorithm 5 Tanh $\langle{\mathcal{F}_{Tanh}(x)}\rangle^{L}$

Input: $P_{0}$, $P_{1}$ have inputs $\langle{x}\rangle_{j}^{L}$ for $P_{j}$
such that $j\in\\{0,1\\}$.

Output: $P_{0}$, $P_{1}$ get $\langle{tanh(c)}\rangle_{0}^{L}$ and
$\langle{tanh(c)}\rangle_{1}^{L}$ where $tanh(c)=\frac{2}{1+e^{-2x}}-1$

Common Randomness: $P_{0}$ and $P_{1}$ hold random shares of zero - $u_{0}$
and $u_{1}$ respectively.

1: $P_{j}$ computes $\langle{a}\rangle_{j}^{L}$ =
$2\times\langle{x}\rangle_{j}^{L}$ for $j\in\\{0,1\\}$

2: $P_{0}$, $P_{1}andP_{2}$ invoke
$\mathcal{F}_{Sigmoid}{(\\{P_{0},P_{1}\\},P_{2})}$ with $P_{j}$ for
$j\in\\{0,1\\}$ having input $\langle{a}\rangle_{j}^{L}$ to learn
$\langle{p}\rangle_{j}^{L}=\langle{\sigma(a)}\rangle_{j}^{L}$.

3: $P_{j}$ for $j\in\\{0,1\\}$ output:
$2\times\langle{P}\rangle_{j}^{L}-j+\langle{u}\rangle_{j}^{L}$

#### Derivative of Tanh, $\mathcal{F^{D}}_{Tanh}(\\{P_{0},P_{1}\\},P_{2})$;
see algorithm 6):

This protocol can be used to securely compute the derivative of the tanh
function. For this protocol, we use the idea that $tanh^{\prime}(x)$ is
related to $\sigma^{\prime}(x)$ as $tanh^{\prime}(x)=4\sigma^{\prime}(2x)$.

Algorithm 6 Derivative of Tanh $\langle{\mathcal{F^{D}}_{Tanh}(x)}\rangle^{L}$

Input: $P_{0}$, $P_{1}$ have inputs $\langle{x}\rangle_{j}^{L}$ for $P_{j}$
such that $j\in\\{0,1\\}$.

Output: $P_{0}$ and $P_{1}$ get $\langle{\tanh^{\prime}(c)}\rangle_{0}^{L}$
and $\langle{\tanh^{\prime}(c)}\rangle_{1}^{L}$ respectively.

Common Randomness: $P_{0}$ and $P_{1}$ hold random shares of zero - $u_{0}$
and $u_{1}$ respectively.

1: $P_{0}$ and $P_{1}$ compute
$\langle{a}\rangle_{j}^{L}=2\times\langle{x}\rangle_{j}^{L}$.

2: $P_{0}$, $P_{1}$, $P_{2}$ invoke $\mathcal{F^{D}}_{sigmoid}$ with
$\langle{a}\rangle_{j}^{L}$ to obtain shares $\langle{b}\rangle_{j}^{L}$.

3: $P_{0}$ and $P_{1}$ output:
$4.\langle{b}\rangle_{j}^{L}+\langle{u}\rangle_{j}^{L}$

#### Softmax, $\mathcal{F}_{Softmax}(\\{P_{0},P_{1}\\},P_{2})$; see algorithm
6):

This protocol can be used to securely compute the softmax function. The two
parties compute the numerator $e^{z_{i}}$ and the denominator
$\sum_{i=1}^{k}e^{z_{i}}$ and invoke the $\mathcal{F}_{Division}$ function to
perform secure division.

Algorithm 7 Softmax, $\mathcal{F}_{Softmax}(\\{P_{0},P_{1}\\},P_{2})$

Input: $P_{0},P_{1}$ hold $\\{\langle z_{i}\rangle_{0}^{L}\\}_{i\in[k]}$ and
$\\{\langle z_{i}\rangle_{1}^{L}\\}_{i\in[k]}$ respectively.

Output: $P_{0},P_{1}$ get $\\{\langle
s_{max}(z_{i})\rangle_{0}^{L}\\}_{i\in[k]}$ and $\\{\langle
s_{max}(z_{i})\rangle_{1}^{L}\\}_{i\in[k]}$ respectively where
$s_{max}(z_{i})=\frac{e^{z_{i}}}{\sum_{i=1}^{k}e^{z_{i}}}$.

Common Randomness: $P_{0}$ and $P_{1}$ hold random shares of zero - $u_{0}$
and $u_{1}$ respectively.

1: for $i=1,2,\ldots,k$ do

2: Parties $P_{0}$, $P_{1}$ and $P_{2}$ invoke
$F_{Exp}(\\{P_{0},P_{1}\\},P_{2})$ with inputs $\\{\langle
z_{i}\rangle_{0}^{L}\\}_{i\in[k]}$ and $\\{\langle
z_{i}\rangle_{1}^{L}\\}_{i\in[k]}$ and $P_{0}$, $P_{1}$ obtain shares $\langle
c_{i}\rangle_{0}^{L}$, $\langle c_{i}\rangle_{1}^{L}$ resp. of $c_{i}^{L}$ =
$e^{z_{i}}$.

3: end for

4: for $j\in\\{0,1\\}$, $P_{j}$ calculates $S_{j}=\sum_{i=1}^{k}\langle
c_{i}\rangle_{j}^{L}$.

5: for $i=1,2,\dots,k$ do

6: Parties $P_{0}$, $P_{1}$ and $P_{2}$ invoke
$\mathcal{F}_{Division}(\\{P_{0},P_{1}\\},P_{2})$ with inputs $\langle
c_{i}\rangle_{j}^{L}$ and $\langle S\rangle_{j}^{L}$.

7: Parties $P_{j}$ for $j\in\\{0,1\\}$ output $\\{\langle
S_{max}(z_{i})\rangle_{j}^{L}\\}_{i\in[k]}+u_{j}$.

8: end for

## Empirical Evaluation

### Implementation and Experimental Setup

Our system is implemented in C++ using standard libraries. To be able to
execute machine learning protocols, we use fixed-point representation in the
integer ring $\mathbb{Z}_{64}$. This allows us to compute precise values using
our machine learning protocols. The fixed-arithmetic used in our
implementation is the same as (Wagh, Gupta, and Chandran 2018) which borrowed
the idea from (Mohassel and Zhang 2017).

Exponentiation is complicated and may overflow the integer ring. So, we
consider the ring $\mathbb{Z}_{128}$ as proposed by (Aly et al. 2020). So far,
we have implemented the Taylor series exponentiation (refer to step 7 of
Algorithm 1) which works for values less than 1. We provide benchmarks for our
logistic sigmoid, tanh and softmax functions, their derivatives and the Taylor
series exponentiation in the secure setting with input vectors of varying
sizes. These benchmarks include the execution time and the total communication
(in MB) of these protocols.

We run our protocols in a LAN setting on an Intel i5 7th Gen with 8GB RAM. The
average bandwidth measured was 40Mb/s for download and 9.2Mb/s for upload and
the average ping was 12ms.

### Experimental Results

We show the average time taken by the protocols in Table 1.

Protocol | Dimension | Time(s) | Comm.(mb)
---|---|---|---
| 64x16 | 0.08 | 0.025
$\mathcal{F}_{Exp}$ | 128x128 | 2.134 | 0.393
| 576x20 | 0.882 | 0.276
| 64x16 | 0.252 | 2.58
$\mathcal{F}_{Sigmoid}$ | 128x128 | 5.631 | 41.288
| 576x20 | 2.615 | 29.03
| 64x16 | 0.275 | 2.58
$\mathcal{F}_{tanh}$ | 128x128 | 5.32 | 41.288
| 576x20 | 2.613 | 29.03
| 64x16 | 0.324 | 2.58
$\mathcal{F}_{Softmax}$ | 128x128 | 5.438 | 41.288
| 576x20 | 2.617 | 29.03
| 64x16 | 0.464 | 2.597
$\mathcal{F^{D}}_{sigmoid}$ | 128x128 | 8.033 | 41.55
| 576x20 | 4.121 | 29.214
| 64x16 | 0.383 | 2.58
$\mathcal{F^{D}}_{tanh}$ | 128x128 | 4.465 | 41.288
| 576x20 | 2.84 | 29.03
| 64x16 | 0.032 | 0.005
$\mathcal{F}_{taylorExp}$ | 128x128 | 0.092 | 0.079
| 576x20 | 0.427 | 0.055
Table 1: Benchmarks in LAN setting on i5 7th Gen

## Conclusion and Future Work

Previous work in secure computation of neural networks has been able to
achieve reasonable efficiency for many real world applications. However, the
lack of options in the choice of protocols that these implementations can run
can limit the capacity of activities the user is able to perform in the secure
setting.

Through the addition of the exponentiation function, it becomes possible to
add a large number of secure activation functions that were previously
impossible to compute. We have shown logistic sigmoid, tanh, softmax and their
derivatives as examples of such functions. We implement these functions with
similar overheads as the previously proposed activation functions.

In its current form, the exponentiation function overflows for small values
making it impractical for usage in regular neural network. This also makes us
unable to use the functions that call the exponentiation function, as part of
their execution, in practical settings.

Exponentiation protocols in a fixed-point MPC setting have been provided by
(Aly and Smart 2019) and implemented by (Aly et al. 2019). Ideas from these
protocols can be used to avoid using a workaround that would involve
increasing the size of shares to accommodate the overflow.

The future work on S++ will be focused on dealing with the overflow such that
the exponentiation function can handle the magnitude of the values that it
would encounter in a regular neural network. Once this is solved, the other
functions that we have proposed, that use exponentiation, will automatically
become suitable for practical usage.

> ## References
>
>   * Agrawal and Srikant (2000) Agrawal, R.; and Srikant, R. 2000.  Privacy-
> preserving data mining.  In _Proceedings of the 2000 ACM SIGMOD
> international conference on Management of data_ , 439–450.
>   * Aly et al. (2020) Aly, A.; Cong, K.; Cozzo, D.; Keller, M.; Orsini, E.;
> Rotaru, D.; Scherer, O.; Scholl, P.; Smart, N.; Tanguy, T.; et al. 2020.
> SCALE–MAMBA v1. 10: Documentation .
>   * Aly et al. (2019) Aly, A.; Orsini, E.; Rotaru, D.; Smart, N. P.; and
> Wood, T. 2019.  Zaphod: Efficiently combining lsss and garbled circuits in
> scale.  In _Proceedings of the 7th ACM Workshop on Encrypted Computing &
> Applied Homomorphic Cryptography_, 33–44.
>   * Aly and Smart (2019) Aly, A.; and Smart, N. P. 2019.  Benchmarking
> privacy preserving scientific operations.  In _International Conference on
> Applied Cryptography and Network Security_ , 509–529. Springer.
>   * Asadi and Jiang (2020) Asadi, B.; and Jiang, H. 2020.  On Approximation
> Capabilities of ReLU Activation and Softmax Output Layer in Neural Networks.
> _arXiv preprint arXiv:2002.04060_ .
>   * Bunn and Ostrovsky (2007) Bunn, P.; and Ostrovsky, R. 2007.  Secure two-
> party k-means clustering.  In _Proceedings of the 14th ACM conference on
> Computer and communications security_ , 486–497.
>   * Catrina and Saxena (2010) Catrina, O.; and Saxena, A. 2010.  Secure
> Computation with Fixed-Point Numbers.  In _International Conference on
> Financial Cryptography and Data Security_. Springer.
>   * Du and Atallah (2001) Du, W.; and Atallah, M. J. 2001.  Privacy-
> preserving cooperative scientific computations.  In _csfw_ , 0273. Citeseer.
>   * Hart (1978) Hart, J. F. 1978.  Computer Approximations.  In _Krieger
> Publishing Co., Inc._ Springer.
>   * Hunt et al. (2018) Hunt, T.; Song, C.; Shokri, R.; Shmatikov, V.; and
> Witchel, E. 2018.  Chiron: Privacy-preserving machine learning as a service.
> _arXiv preprint arXiv:1803.05961_ .
>   * Jagannathan and Wright (2005) Jagannathan, G.; and Wright, R. N. 2005.
> Privacy-preserving distributed k-means clustering over arbitrarily
> partitioned data.  In _Proceedings of the eleventh ACM SIGKDD international
> conference on Knowledge discovery in data mining_ , 593–599.
>   * Juvekar, Vaikuntanathan, and Chandrakasan (2018) Juvekar, C.;
> Vaikuntanathan, V.; and Chandrakasan, A. 2018.  $\\{$GAZELLE$\\}$: A low
> latency framework for secure neural network inference.  In _27th
> $\\{$USENIX$\\}$ Security Symposium ($\\{$USENIX$\\}$ Security 18)_,
> 1651–1669.
>   * Kent and Salem (2019) Kent, D.; and Salem, F. 2019.  Performance of
> three slim variants of the long short-term memory (LSTM) layer.  In _2019
> IEEE 62nd International Midwest Symposium on Circuits and Systems (MWSCAS)_
> , 307–310. IEEE.
>   * Krizhevsky, Sutskever, and Hinton (2017) Krizhevsky, A.; Sutskever, I.;
> and Hinton, G. E. 2017.  Imagenet classification with deep convolutional
> neural networks.  _Communications of the ACM_ 60(6): 84–90.
>   * Le, Jaitly, and Hinton (2015) Le, Q. V.; Jaitly, N.; and Hinton, G. E.
> 2015.  A simple way to initialize recurrent networks of rectified linear
> units.  _arXiv preprint arXiv:1504.00941_ .
>   * Lindell and Pinkas (2000) Lindell, Y.; and Pinkas, B. 2000.  Privacy
> preserving data mining.  In _Annual International Cryptology Conference_ ,
> 36–54. Springer.
>   * Liu et al. (2017) Liu, J.; Juuti, M.; Lu, Y.; and Asokan, N. 2017.
> Oblivious neural network predictions via minionn transformations.  In
> _Proceedings of the 2017 ACM SIGSAC Conference on Computer and
> Communications Security_ , 619–631.
>   * Maksutov (2018) Maksutov, R. 2018.  Deep study of a not very deep neural
> network. Part 2: Activation functions.  https://towardsdatascience.com/deep-
> study-of-a-not-very-deep-neural-network-part-2-activation-functions-
> fd9bd8d406fc.  Accessed: 2020-11-06.
>   * Mohassel and Zhang (2017) Mohassel, P.; and Zhang, Y. 2017.  SecureML: A
> system for scalable privacy-preserving machine learning.  In _2017 IEEE
> Symposium on Security and Privacy (SP)_ , 19–38. IEEE.
>   * Ohrimenko et al. (2016) Ohrimenko, O.; Schuster, F.; Fournet, C.; Mehta,
> A.; Nowozin, S.; Vaswani, K.; and Costa, M. 2016.  Oblivious multi-party
> machine learning on trusted processors.  In _25th $\\{$USENIX$\\}$ Security
> Symposium ($\\{$USENIX$\\}$ Security 16)_, 619–636.
>   * Pascanu, Mikolov, and Bengio (2013) Pascanu, R.; Mikolov, T.; and
> Bengio, Y. 2013.  On the difficulty of training recurrent neural networks.
> In _International conference on machine learning_ , 1310–1318.
>   * Patra and Suresh (2020) Patra, A.; and Suresh, A. 2020.  BLAZE: Blazing
> Fast Privacy-Preserving Machine Learning.  _arXiv preprint arXiv:2005.09042_
> .
>   * Riazi et al. (2018) Riazi, M. S.; Weinert, C.; Tkachenko, O.; Songhori,
> E. M.; Schneider, T.; and Koushanfar, F. 2018.  Chameleon: A hybrid secure
> computation framework for machine learning applications.  In _Proceedings of
> the 2018 on Asia Conference on Computer and Communications Security_ ,
> 707–721.
>   * Sanil et al. (2004) Sanil, A. P.; Karr, A. F.; Lin, X.; and Reiter, J.
> P. 2004.  Privacy preserving regression modelling via distributed
> computation.  In _Proceedings of the tenth ACM SIGKDD international
> conference on Knowledge discovery and data mining_ , 677–682.
>   * Slavkovic, Nardi, and Tibbits (2007) Slavkovic, A. B.; Nardi, Y.; and
> Tibbits, M. M. 2007.  ” Secure” Logistic Regression of Horizontally and
> Vertically Partitioned Distributed Databases.  In _Seventh IEEE
> International Conference on Data Mining Workshops (ICDMW 2007)_ , 723–728.
> IEEE.
>   * Vaidya, Yu, and Jiang (2008) Vaidya, J.; Yu, H.; and Jiang, X. 2008.
> Privacy-preserving SVM classification.  _Knowledge and Information Systems_
> 14(2): 161–178.
>   * Wagh, Gupta, and Chandran (2018) Wagh, S.; Gupta, D.; and Chandran, N.
> 2018.  SecureNN: Efficient and Private Neural Network Training.  _IACR
> Cryptology ePrint Archive_ 2018: 442.
>   * Xu, Huang, and Li (2016) Xu, B.; Huang, R.; and Li, M. 2016.  Revise
> saturated activation functions.  _arXiv preprint arXiv:1602.05980_ .
>   * Yu, Vaidya, and Jiang (2006) Yu, H.; Vaidya, J.; and Jiang, X. 2006.
> Privacy-preserving svm classification on vertically partitioned data.  In
> _Pacific-asia conference on knowledge discovery and data mining_ , 647–656.
> Springer.
>

## Appendix

Universal approximation of softmax: (Asadi and Jiang 2020) provide theoretical
proofs on the universal approximation of ReLU used in hidden layers of deep
networks and softmax used in the output layer for multiclass classification
tasks. They find that such a proof for the softmax is missing in the
literature despite it being widely used for these tasks, and provide their own
to prove the power of the softmax used in output layers to approximate
indicator functions.

ReLU and gradients: (Pascanu, Mikolov, and Bengio 2013) delve into the
exploding and vanishing gradients of activation functions that tend to cause
problems in recurrent neural networks because of the temporal correlations
they capture (by forming connections between hidden layers). Long-term
components add up and the norm of the gradient of the hidden layers explode.
ReLU’s unbounded nature exacerbates this problem. There have been a few
notable improvements to assuage this, namely gradient clipping, identity
initialization (Le, Jaitly, and Hinton 2015), and gated RNNs like LSTMs and
GRUs. However, even in LSTMs, the tanh often gives more consistent results
((Kent and Salem 2019)) and ReLU tends to require greater fine-tuning.

The Leaky Tanh: (Xu, Huang, and Li 2016) propose a ‘leaky’ tanh that penalizes
the negative part i.e., instead of the regular tanh function, the penalized
tanh behaves like this:

$penalized\\_tanh(x)=\begin{dcases}tanh(x)&x>0\\\ a.tanh(x)&otherwise\\\
\end{dcases}$ (1)

where $a\in(0,1)$. They also found this variation to be two times faster than
the tanh.