an effective heating of the mechanical modes and a thermal contribution to
reduced cavity state which no longer vanishes at the respective decoupling
times of the optical and mechanical states. To see this explicitly we note
from Appendix D that under unitary evolution, the components of the (reduced)
cavity state when the mechanics is initially in a thermal state are given by,

$\displaystyle\rho_{p,q,m,n}^{cav}(t)=$ $\displaystyle
a_{p}a_{m}^{*}b_{q}b_{n}^{*}e^{-i(\omega-A)(p+q-m-n)t}e^{iB(m^{2}+n^{2}-p^{2}-q^{2})}e^{iD(pq-
mn)}e^{-\frac{1}{2}|K_{\hat{N}_{c}}|^{2}(\bar{n}_{s}+\frac{1}{2})(q-p+m-n)^{2}}$
(132) $\displaystyle\quad\times
e^{-\frac{1}{2}\left(|K_{\hat{N}_{s}}|^{2}(\bar{n}_{s}+\frac{1}{2})(p+q-m-n)^{2}-(C_{s}K_{\hat{N}_{s}}^{*}-C_{s}^{*}K_{\hat{N}_{s}})(p+q-m-n)\right)}\,.$

for arbitrary initial cavity states $\ket{\psi_{1}(0)}=\sum a_{m}\ket{m}$,
$\ket{\psi_{2}(0)}=\sum b_{m}\ket{m}$. This implies we can formally replace
the thermal occupation $\bar{n}_{\mu}$ by a time dependent term,

$\bar{n}_{\mu}(t)=\bar{n}_{\mu}(0)+\frac{\kappa_{\mu}^{\mathrm{dec}}(t)}{|K_{\hat{N}_{\mu}}|^{2}}\,,$
(133)

where

$\kappa_{\mu}^{\mathrm{dec}}(t)=\int_{0}^{t}\frac{8\Upsilon}{\hbar^{2}}x_{0,\mu}^{2}\Re[e^{-i\omega_{\mu}(s-t)}K_{\hat{N}_{\mu}}(s-t)]^{2}ds.$
(134)

For the modulation frequency $\Omega_{n}=(1-1/n)(\omega_{c}+\omega_{s})$ and
$\phi_{1}=\pi/2$, $\kappa_{\mu}^{\mathrm{dec}}(t_{q})$ can be evaluated to
zeroth order in $\gamma_{g}$ as

$\displaystyle\kappa_{\mu}^{\mathrm{dec}}(t_{q})=$
$\displaystyle\frac{4\Upsilon}{\hbar^{2}}x_{0,m}^{2}\frac{qn\pi}{\omega_{m}}\left(3k_{0}^{2}+\epsilon^{2}\frac{n^{2}(1+(2n-1)^{2})}{2(2n-1)^{2}}\right)$
(135) $\displaystyle\approx$ $\displaystyle\frac{\Upsilon
x_{0,m}^{2}}{\hbar^{2}\omega_{m}\gamma_{g}}D(t_{q}),$

for $\mu=c,s$.

##  Appendix C Evaluation of Determinant-based Entanglement Witnesses for
Continuous Variable States

In section V.1 it was noted that the idealised limiting state of the final
(entangled) optical fields will be of the form
$(1+i)\ket{\alpha_{1},\alpha_{2}}+(1-i)\ket{-\alpha_{1},-\alpha_{2}}$. While
many of the simplest continuous variable criteria fail to identify
entanglement in states of this form, Shchukin and Vogel [40] introduced an
infinite series of inequalities for which the violation of any is sufficient
to witness entanglement. These can be conveniently written as the determinant
of a matrix of moments of the state, where negativity implies the violation of
one set of inequalities. The simplest which is suitable for our purposes is
given by,

$\mathcal{W}_{1}(t)=\begin{vmatrix}1&\braket{a_{2}^{\dagger}}&\braket{a_{1}a_{2}^{\dagger}}\\\
\braket{a_{2}}&\braket{a_{2}^{\dagger}a_{2}}&\braket{a_{1}a_{2}^{\dagger}a_{2}}\\\
\braket{a_{1}^{\dagger}a_{2}}&\braket{a_{1}^{\dagger}a_{2}^{\dagger}a_{2}}&\braket{a_{1}^{\dagger}a_{1}a_{2}^{\dagger}a_{2}}\end{vmatrix}$
(136)

For convenience we move to the Heisenberg picture,
$\hat{O}(t)=\hat{U}^{\dagger}\hat{O}\hat{U}$, where in the symmetric cavity
limit the evolution operator is given explicitly by (13) and (27),

$\displaystyle U(t)=$ $\displaystyle
e^{i(-\omega_{0}t+A})(\hat{N}_{1}+\hat{N}_{2})e^{-i\omega_{c}\hat{b}_{c}^{\dagger}\hat{b}_{c}t}e^{-i\omega_{s}\hat{b}_{s}^{\dagger}\hat{b}_{s}t}e^{i\theta}e^{iB(\hat{N}_{1}^{2}+\hat{N}_{2}^{2})}e^{iD\hat{N}_{1}\hat{N}_{2}}D_{c}[-K_{\hat{N}_{c}}(\hat{N}_{2}-\hat{N}_{1})]D_{s}[C_{s}-K_{\hat{N}_{s}}(\hat{N}_{1}+\hat{N}_{2})],$
(137)

with $\omega_{s}=\omega_{m}\sqrt{1-4\gamma_{g}}$ and $\omega_{c}=\omega_{m}$.
Using the relation,

$D(C-K(\hat{N}_{2}+(-1)^{k}\hat{N}_{1}))=D(C)D(-K\hat{N}_{2})D((-1)^{k+1}K\hat{N}_{1})e^{+\frac{1}{2}(CK^{*}-C^{*}K)(\hat{N}_{2}+(-1)^{k}\hat{N}_{1})},$
(138)

Together with
$D_{\mu}^{\dagger}(X\hat{N}_{j})\hat{a}_{j}D_{\mu}(X\hat{N}_{j})=D_{\mu}(X)\hat{a}_{j}$,
and $e^{-X\hat{N}^{2}}\hat{a}e^{X\hat{N}^{2}}=e^{X(2\hat{N}+1)}\hat{a}$ we
find,

$\displaystyle\hat{a}_{1}(t)$
$\displaystyle=e^{i\phi_{0}}\hat{\mathcal{B}}_{1}e^{i(D\hat{N}_{2}+2B\hat{N}_{1})}\hat{a}_{1}$
(139) $\displaystyle\hat{a}_{2}(t)$
$\displaystyle=e^{i\phi_{0}}\hat{\mathcal{B}}_{2}e^{i(D\hat{N}_{1}+2B\hat{N}_{2})}\hat{a}_{2}$
$\displaystyle(\hat{a}_{1}\hat{a}^{\dagger}_{2})(t)$
$\displaystyle=\hat{\mathcal{B}}_{12}\hat{a}^{\dagger}_{2}e^{i(2B-D)(\hat{N}_{1}-\hat{N}_{2})}\hat{a}_{1}$
$\displaystyle(\hat{a}^{\dagger}_{j}\hat{a}_{j})(t)$
$\displaystyle=\hat{a}^{\dagger}_{j}\hat{a}_{j}$
$\displaystyle(\hat{a}_{1}\hat{a}^{\dagger}_{2}\hat{a}_{2})(t)$
$\displaystyle=\hat{a}_{1}(t)\hat{a}^{\dagger}_{2}\hat{a}_{2}=e^{iD}\hat{a}^{\dagger}_{2}\hat{a}_{1}(t)\hat{a}_{2}$

Where $\phi_{0}=-\omega_{0}t+A+B$ and,

$\displaystyle\hat{\mathcal{B}}_{j}$
$\displaystyle=e^{+\frac{1}{2}(C_{s}K_{\hat{N}_{s}}^{*}-C_{s}^{*}K_{\hat{N}_{s}})}D_{c}(-(-1)^{j}K_{\hat{N}_{c}})D_{s}(-K_{\hat{N}_{s}})$
(140) $\displaystyle\hat{\mathcal{B}}_{j}^{2}$
$\displaystyle=e^{+(C_{s}K_{\hat{N}_{s}}^{*}-C_{s}^{*}K_{\hat{N}_{s}})}D_{c}(-(-1)^{j}2K_{\hat{N}_{c}})D_{s}(-2K_{\hat{N}_{s}})$
$\displaystyle\hat{\mathcal{B}}_{12}$
$\displaystyle=\hat{\mathcal{B}}_{1}\hat{\mathcal{B}}_{2}^{\dagger}=D_{c}(+2K_{\hat{N}_{c}})$
$\displaystyle\hat{\mathcal{B}}_{1}\hat{\mathcal{B}}_{2}$
$\displaystyle=e^{+(C_{s}K_{\hat{N}_{s}}^{*}-C_{s}^{*}K_{\hat{N}_{s}})}D_{s}(-2K_{\hat{N}_{s}})$

are operators acting on the mechanical modes. Here we consider two scenarios:
1) The mechanical modes initially in a coherent state $\beta_{\mu}$, and 2)
The mechanical modes initially in a thermal state $\rho_{\rm{th},\mu}$. For
each case the expectation value of the displacement operator is given by,

$\displaystyle\langle D_{\mu}(X)\rangle_{\beta}$
$\displaystyle=e^{-\frac{1}{2}|X|^{2}}e^{X\beta^{*}-X^{*}\beta}$ (141)
$\displaystyle\langle D_{\mu}(X)\rangle_{\rho_{th}}$
$\displaystyle=e^{-\frac{1}{2}|X|^{2}}e^{-(\Re[X]^{2}+\Im[X]^{2})\bar{n}_{\mu}}=e^{-|X|^{2}(\bar{n}_{\mu}+\frac{1}{2})}$

So we have

$\displaystyle\langle\mathcal{B}_{j}\rangle_{\beta}$
$\displaystyle=e^{-\frac{1}{2}\left(|K_{\hat{N}_{c}}|^{2}+|K_{\hat{N}_{s}}|^{2}\right)}e^{+\frac{1}{2}(C_{s}K_{\hat{N}_{s}}^{*}-C_{s}^{*}K_{\hat{N}_{s}})}e^{(-1)^{j}\tilde{\beta}_{c}+\tilde{\beta}_{s}},$
(142) $\displaystyle\langle\mathcal{B}_{j}^{2}\rangle_{\beta}$
$\displaystyle=e^{-2\left(|K_{\hat{N}_{c}}|^{2}+|K_{\hat{N}_{s}}|^{2}\right)}e^{+(C_{s}K_{\hat{N}_{s}}^{*}-C_{s}^{*}K_{\hat{N}_{s}})}e^{2[(-1)^{j}\tilde{\beta}_{c}+\tilde{\beta}_{s}]},$
$\displaystyle\langle\mathcal{B}_{12}\rangle_{\beta}$
$\displaystyle=e^{-2|K_{\hat{N}_{c}}|^{2}}e^{-2\tilde{\beta}_{c}},$
$\displaystyle\langle\mathcal{B}_{1}\mathcal{B}_{2}\rangle_{\beta}$
$\displaystyle=e^{-2|K_{\hat{N}_{s}}|^{2}}e^{+(C_{s}K_{\hat{N}_{s}}^{*}-C_{s}^{*}K_{\hat{N}_{s}})}e^{2\tilde{\beta}_{s}}$

where
$\tilde{\beta}_{\mu}=-(K_{\hat{N}_{\mu}}\beta_{\mu}^{*}-K_{\hat{N}_{\mu}}^{*}\beta_{\mu})$
and

$\displaystyle\langle\mathcal{B}_{j}\rangle_{\rho_{th}}$
$\displaystyle=e^{-\frac{1}{2}\left(|K_{\hat{N}_{c}}|^{2}(\bar{n}_{c}+\frac{1}{2})+|K_{\hat{N}_{s}}|^{2}(\bar{n}_{s}+\frac{1}{2})\right)}e^{+\frac{1}{2}(C_{s}K_{\hat{N}_{s}}^{*}-C_{s}^{*}K_{\hat{N}_{s}})},$
(143) $\displaystyle\langle\mathcal{B}_{j}^{2}\rangle_{\rho_{th}}$
$\displaystyle=e^{-2\left(|K_{\hat{N}_{c}}|^{2}(\bar{n}_{c}+\frac{1}{2})+|K_{\hat{N}_{s}}|^{2}(\bar{n}_{s}+\frac{1}{2})\right)}e^{+(C_{s}K_{\hat{N}_{s}}^{*}-C_{s}^{*}K_{\hat{N}_{s}})},$
$\displaystyle\langle\mathcal{B}_{12}\rangle_{\rho_{th}}$
$\displaystyle=e^{-2|K_{\hat{N}_{c}}|^{2}(\bar{n}_{c}+\frac{1}{2})},$
$\displaystyle\langle\mathcal{B}_{1}\mathcal{B}_{2}\rangle_{\rho_{th}}$
$\displaystyle=e^{-2|K_{\hat{N}_{s}}|^{2}(\bar{n}_{s}+\frac{1}{2})}e^{+(C_{s}K_{\hat{N}_{s}}^{*}-C_{s}^{*}K_{\hat{N}_{s}})}$

The final ingredient are the expectation values with respect to the cavity
states, which we take to be coherent with parameters $\alpha_{1}$ and
$\alpha_{2}$ for the two cavities, respectively. Using

$\langle
e^{iX\hat{N}_{j}}\hat{a}_{j}^{n}\rangle_{\alpha_{j}}=\alpha_{j}^{n}e^{-|\alpha_{j}|^{2}\left(1-e^{iX}\right)},$
(144)

we have,

$\displaystyle\langle\hat{a}_{1}(t)\rangle$
$\displaystyle=\alpha_{1}e^{i\phi_{0}}\langle\mathcal{B}_{1}\rangle
e^{-|\alpha_{1}|^{2}\left(1-e^{2iB}\right)}e^{-|\alpha_{2}|^{2}\left(1-e^{iD}\right)},$
(145) $\displaystyle\langle\hat{a}_{2}(t)\rangle$
$\displaystyle=\alpha_{2}e^{i\phi_{0}}\langle\mathcal{B}_{2}\rangle
e^{-|\alpha_{2}|^{2}\left(1-e^{2iB}\right)}e^{-|\alpha_{1}|^{2}\left(1-e^{iD}\right)},$
$\displaystyle\langle(\hat{a}_{1}\hat{a}^{\dagger}_{2})(t)\rangle$
$\displaystyle=\alpha_{1}\alpha_{2}^{*}\langle\mathcal{B}_{12}\rangle
e^{-|\alpha_{1}|^{2}\left(1-e^{-i(D-2B)}\right)}e^{-|\alpha_{2}|^{2}\left(1-e^{i(D-2B)}\right)}$
$\displaystyle\langle(\hat{a}_{j}^{\dagger}\hat{a}_{j})(t)\rangle$
$\displaystyle=|\alpha_{j}|^{2}$
$\displaystyle\langle(\hat{a}_{1}\hat{a}^{\dagger}_{2}\hat{a}_{2})(t)\rangle$
$\displaystyle=\alpha_{1}|\alpha_{2}|^{2}e^{iD}e^{i\phi_{0}}\langle\mathcal{B}_{1}\rangle
e^{-|\alpha_{1}|^{2}\left(1-e^{2iB}\right)}e^{-|\alpha_{2}|^{2}\left(1-e^{iD}\right)}$
$\displaystyle\langle(\hat{a}_{1}^{\dagger}\hat{a}_{1}\hat{a}^{\dagger}_{2}\hat{a}_{2})(t)\rangle$
$\displaystyle=|\alpha_{1}|^{2}|\alpha_{2}|^{2}.$

These expressions are generally too complicated to gain much insight. However,
by considering the form of the determinant (136) we observe that the phase
$\phi_{0}$ ultimately vanishes together with the imaginary exponential terms
in (142) and (143). Similarly, only the modulus of $\alpha_{1}$ and
$\alpha_{2}$ is important and so from symmetry of the set-up this implies that
the optimal choice of initial cavity states is
$\alpha_{1}=\alpha_{2}\equiv\alpha$ and $|\alpha|^{2}=N_{p}$. With these
restrictions, the determinant witness (136) can be written as,

$\displaystyle\mathcal{W}_{1}(t_{q})=$ $\displaystyle
N_{p}^{3}\left(1-e^{-4\kappa_{c}^{\mathrm{th}}-4N_{p}(1-\cos{(D-2B)})}-2e^{-\kappa_{c}^{\mathrm{th}}-\kappa_{s}^{\mathrm{th}}-2N_{p}(2-\cos
D-\cos 2B)}\right.$ (146)
$\displaystyle\left.+2e^{-3\kappa_{c}^{\mathrm{th}}-\kappa_{s}^{\mathrm{th}}-2N_{p}(3-\cos
D-\cos{2B}-\cos{(D-2B)})}\cos D\right),$

where $\kappa_{\mu}^{\mathrm{th}}=|K_{\hat{N}_{\mu}}|^{2}(\bar{n}_{\mu}+1/2)$.

We note, however, that the $B$ terms should not in principle contribute to the
entanglement, even though they have a non-trivial effect on the witness. An
analogous situation arises in the estimation of the gravitational field by a
single optomechanical cavity, where it is known that the quantum Fisher
information (QFI) is independent of the corresponding $B$ coefficient at the
decoupling time, however certain POVMs only saturate the QFI when this is an
integer multiple of $2\pi$ [21]. This is occurs when using homodyne
measurements, for example, and typically require $k_{0}\in\mathbb{Z}$.
Interestingly, we find from a numerical analysis, that the optimal values of
$B(t_{q})$ for the witness are integer multiples of $\pi$. Considering a
measurement time $t_{q}=qn\pi/(1/\omega_{c}+1/\omega_{s})$, with the
optomechanical coupling modulated at frequency
$\Omega_{n}=\frac{1}{2}(1-1/n)(\omega_{c}+\omega_{s})$, we find from
$B=B_{c}+B_{s}$ that,

$B(t_{q})=2qn\pi(1+3\gamma_{g})\left(k_{0}^{2}+\frac{n^{2}}{2(2n-1)}\epsilon^{2}\right).$
(147)

Thus in the large $n$ limit, we also see that approximately integer values of
$k_{0}$ and $\epsilon$ are favourable111111On the other hand, heterodyne
measurements achieve the same scaling as the QFI (up to a factor of $2$)
independent of this constraint, and so we may hypothesize that a heterodyne-
based witness could be more convenient from an experimental perspective.
Adopting these $B(t_{q})$ values leads to immediately to the witness in the
main text (28). The prefactor (generically $N_{p}^{3}$) should not be
attributed to the entanglement scaling estimated in section IV.2. Instead this
term acts like a signal gain, which should be offset against a proper noise
analysis of measurements of the moments used in (136).

When the initial mode temperatures are equal, we find that
$\kappa_{c}^{\mathrm{th}}(t_{q})=\kappa_{s}^{\mathrm{th}}(t_{q})$ to first
order in $\gamma_{g}$,

$\displaystyle\kappa_{\mu}^{\mathrm{th}}(t_{q})$
$\displaystyle\approx\frac{(2\pi
n\gamma_{g})^{2}\left(k_{0}^{2}(2n-1)^{2}+\left(n-1\right)^{2}n^{2}\epsilon^{2}\right)}{2(2n-1)^{2}}(\bar{n}_{\mu}+1/2)$
(148) $\displaystyle\approx\frac{1}{2}(\pi
n\gamma_{g})^{2}\left(4k_{0}^{2}+n^{2}\epsilon^{2}\right)(\bar{n}_{\mu}+1/2)$

where in the second line we have assumed $n$ large. Thus

$\mathcal{W}_{1}(t_{q})=4N_{p}^{3}e^{-4y}\left(e^{2y}\sinh^{2}(y)-z^{2}\right),$
(149)

where $y=2N_{p}z^{2}+\kappa^{\mathrm{th}}(t_{q})$ and $z=\sin(D(t_{q})/2)$.
Taking the derivative with respect to $z$, the minimum will satisfy (34),

$4z_{\min}^{2}+e^{2y_{\min}}=1+\frac{2z_{\min}}{4N_{p}z_{\min}+\kappa^{\mathrm{th}\prime}|_{z_{\min}}}.$
(150)

For $\kappa^{\mathrm{th}\prime}|_{z_{\min}}\approx 0$, a formal solution for
$z_{\min}$ can be written as,

$z_{\min}=\frac{\sqrt{1+2N_{p}-2W(e^{\frac{1}{2}+N_{p}+2\kappa^{\mathrm{th}}}N_{p})}}{2\sqrt{2N_{p}}}$
(151)

where $W(x)$ is the Lambert-$W$ function defined as the solution to
$W(x)e^{W(x)}=x$. For large $x$, we can approximate $W(x)$ asymptotically as

$W(x)=\ln\left(\frac{x}{W(x)}\right)\approx\ln\left(\frac{x}{\ln{x}}\right)$
(152)

then to first order in $1/N_{p}$, the term
$W(e^{\frac{1}{2}+N_{p}+2\kappa^{\mathrm{th}}}N_{p})=\frac{1}{2}+N_{p}-\frac{1}{2N_{p}}[1-4(N_{p}-1)\kappa^{\rm{th}}]$,
and so

$z_{\min}\approx\frac{\sqrt{1-4(N_{p}-1)\kappa^{\mathrm{th}}}}{2\sqrt{2}N_{p}}\approx\frac{\sqrt{1-4N_{p}\kappa^{\mathrm{th}}}}{2\sqrt{2}N_{p}}.$
(153)

which agrees with the simple analysis provided in the main text.

When $\kappa^{\prime}|_{z_{\min}}$ is significant, it is convenient to first
write $\kappa$ explicitly as a function of $z$. For $\epsilon=0$,
$\kappa^{\mathrm{th}}(z)\approx\Gamma^{\rm{th}}_{0}z^{2}$, with
$\Gamma_{0}^{\rm{th}}=\frac{2\bar{n}+1}{36k_{0}^{2}}$ (see (33)), which
immediately leads to (38) under the substitutions $N_{p}\rightarrow
N_{p}+\Gamma_{0}^{\rm{th}}/2$ and $\kappa^{\mathrm{th}}\rightarrow 0$.
However, for modulated coupling, where $\kappa^{\rm{th}}=\Gamma z^{4/3}$,
finding a formal solution to the minimum condition is difficult, and so we
instead resort to an expansion of (150) under the assumption
$z_{\min},N_{p}z_{\min}^{2}\ll 1$, and $N_{p}\gg 1$ (which are all consistent
with the expected experimental regimes). In order to preserve the correct
behaviour as $\kappa^{\rm{th}}\rightarrow 0$, it is convenient to first
multiply (150) by
$(3/2)(4N_{p}z_{\min}+\kappa^{\prime}|_{z_{\min}})/z_{\min}$, before expanding
the resulting expression to second order in $z_{\min}$. A slight rearrangement
leads to,

$4(6N_{p}^{2}+\Gamma^{3})w^{3}+20N_{p}\Gamma w^{2}+4\Gamma^{2}w-3=0$ (154)

where $\Gamma=\frac{1}{2}(\bar{n}+\frac{1}{2})\left(\frac{\pi
q\gamma_{g}}{\epsilon}\right)^{\frac{2}{3}}$, and we have made the
substitution $z_{\min}\rightarrow w^{3/2}$. The solution for $w$ up to first
order in $1/N_{p}$ is given by,

$w=\frac{1-v}{2N_{p}^{2/3}},\qquad\text{where
}v=\frac{5\Gamma}{9N_{p}^{1/3}}=\frac{5}{18}\frac{\bar{n}+\frac{1}{2}}{N_{p}^{1/3}}\left(\frac{\pi
q\gamma_{g}}{\epsilon}\right)^{\frac{2}{3}}.$ (155)

Transforming back to $z_{\min}$ leads to (42) in the main text. This
corresponds to an optimal entanglement time of,

$t_{1}=\left(\frac{2\sqrt{2}\pi^{2}}{N_{p}\omega_{m}^{3}\epsilon^{2}\gamma_{g}}\right)^{\frac{1}{3}}\sqrt{1-v}.$
(156)

We now expand the exponential terms in $\mathcal{W}_{1}(t_{1})$ in powers of
$1/N_{p}$, ignoring overall terms proportional to any inverse powers of
$N_{p}$ (which is consistent with our final bound on $\Gamma$). This leads to,

$\mathcal{W}_{1}(z_{\min})\approx\frac{3}{8}-\frac{N_{p}}{4}+\frac{1}{2}N_{p}^{2/3}\Gamma-\frac{4}{9}N_{p}^{1/3}\Gamma^{2}+\frac{35\Gamma^{3}}{162}.$
(157)

In the limit $N_{p}\gg 1$ the first term can be neglected, and so we again
recover the asymptotic noiseless witness minimum of $-N_{p}/4$. The solution
for $\Gamma$ such that (157) is half the optimal value is then approximately,

$\Gamma\lesssim 0.33N_{p}^{1/3}.$ (158)

###  C.1 Decoherence

The comparatively slower scaling of the temperature bound is a direct result
of measuring close to (or at) the cavity-mechanics decoupling time. This
highlights the advantage of the fractional frequencies introduced in [21]. On
the other hand mechanical noise depends on integrating over the full motion
and therefore always contains a term independent of $\gamma_{g}$. In Appendix
B we saw that thermal decoherence effects can be formally included by
rewriting expectation values of the mechanical mode operators as (see equation
(133))

$\displaystyle\langle\mathcal{B}_{j}\rangle_{\rho_{th}}$
$\displaystyle=e^{-\frac{1}{2}\left(|K_{\hat{N}_{c}}|^{2}(\bar{n}_{c}+\frac{1}{2})+|K_{\hat{N}_{s}}|^{2}(\bar{n}_{s}+\frac{1}{2})+\kappa_{c}^{\mathrm{dec}}+\kappa_{s}^{\mathrm{dec}}\right)}e^{+\frac{1}{2}(C_{s}K_{\hat{N}_{s}}^{*}-C_{s}^{*}K_{\hat{N}_{s}})},$
(159) $\displaystyle\langle\mathcal{B}_{j}^{2}\rangle_{\rho_{th}}$
$\displaystyle=e^{-2\left(|K_{\hat{N}_{c}}|^{2}(\bar{n}_{c}+\frac{1}{2})+|K_{\hat{N}_{s}}|^{2}(\bar{n}_{s}+\frac{1}{2})+\kappa_{c}^{\mathrm{dec}}+\kappa_{s}^{\mathrm{dec}}\right)}e^{+(C_{s}K_{\hat{N}_{s}}^{*}-C_{s}^{*}K_{\hat{N}_{s}})},$
$\displaystyle\langle\mathcal{B}_{12}\rangle_{\rho_{th}}$
$\displaystyle=e^{-2|K_{\hat{N}_{c}}|^{2}(\bar{n}_{c}+\frac{1}{2})}e^{-2\kappa_{c}^{\mathrm{dec}}},$
$\displaystyle\langle\mathcal{B}_{1}\mathcal{B}_{2}\rangle_{\rho_{th}}$
$\displaystyle=e^{-2|K_{\hat{N}_{s}}|^{2}(\bar{n}_{s}+\frac{1}{2})}e^{-2\kappa_{s}^{\mathrm{dec}}}e^{+(C_{s}K_{\hat{N}_{s}}^{*}-C_{s}^{*}K_{\hat{N}_{s}})}$

where $\kappa_{\mu}^{\mathrm{dec}}(t)$ is defined in equation (134). We find,

$\mathcal{W}_{1}(t_{1})=N_{p}^{3}\left(1-e^{-4\kappa_{c}-4N_{p}(1-\cos
D)}-2e^{-\kappa_{c}-\kappa_{s}-2N_{p}(1-\cos
D)}+2e^{-3\kappa_{c}-\kappa_{s}-4N_{p}(1-\cos D)}\cos D\right)\,,$ (160)

where $\kappa_{\mu}=\kappa_{\mu}^{\mathrm{dec}}+\kappa_{\mu}^{\mathrm{th}}$.
To zeroth order in $\gamma_{g}$,
$\kappa_{c}^{\mathrm{dec}}(t_{q})\approx\kappa_{s}^{\mathrm{dec}}(t_{q})\gg\kappa_{s}^{\mathrm{th}}(t_{q})$
and so from (135) we can set
$\kappa_{\mu}(t_{q})=\kappa(z(t_{q}))\approx\Gamma_{\rm{dec}}z(t_{q})$, with
$\Gamma_{\rm{dec}}=\frac{2\Upsilon
x_{0,m}^{2}}{\hbar^{2}\omega_{m}\gamma_{g}}$. In this case the analysis
proceeds in a similar fashion to the above. Substituting for $\kappa(z)$, we
find the analogue of (154) as,

$2(12N_{p}^{2}+\Gamma_{\rm{dec}}^{4}+12N_{p}(1+\Gamma_{\rm{dec}}^{2}))z_{\min}^{2}+3\Gamma_{\rm{dec}}(6N_{p}+\Gamma_{\rm{dec}})z_{\min}+3(\Gamma_{\rm{dec}}-1)=0.$
(161)

Expanding to first orders in $1/N_{p}$, followed by $\Gamma_{\rm{dec}}$, the
solution is given by,

$z_{\min}=\frac{1}{2\sqrt{2}N_{p}}\left(1-\frac{3}{2\sqrt{2}}\Gamma_{\rm{dec}}\right).$
(162)

while to second order in $\Gamma_{\rm{dec}}$, the corresponding value of the
witness is,

$\mathcal{W}_{1}(z_{\min})=-\frac{N_{p}}{4}\left(1-2\sqrt{2}\Gamma_{\rm{dec}}+\frac{5}{2}\Gamma_{\rm{dec}}^{2}\right).$
(163)

Now, solving $\min\mathcal{W}_{1}=\min
a\mathcal{W}_{1|_{\Gamma_{\rm{dec}}=0}}$ gives,

$\Gamma_{\rm{dec}}\leq\frac{\sqrt{2}}{5}\left(2-\sqrt{5a-1}\right),$ (164)

so for $a=1/2$, $\Gamma_{\rm{dec}}\lesssim 0.219$, which is typically accurate
to within a little over $1\%$ of the (numerical) exact bound. In contrast, the
absolute bound for this decoherence model (i.e. when infinitesimal negative
values of $\mathcal{W}_{1}$ can be measured) is given by
$\Gamma_{\rm{dec}}\leq 1$. We then have,

$\gamma_{R}k_{B}T\lesssim\begin{cases}\frac{\hbar
Gm}{2d^{3}},\qquad\qquad&\text{Verify {any} entanglement}\\\ 0.11\frac{\hbar
Gm}{d^{3}},&\text{Achieve half the minimum possible value of the
witness.}\end{cases}$ (165)

##  Appendix D Logarithmic negativity for optical Fock state superpositions

In the case of finite dimensional bipartite systems, a convenient entanglement
measure for mixed states is given by the logarithmic logarithmic negativity,

$LN(\rho)=\log||\rho^{T_{2}}||_{1},$ (166)

where $\rho^{T_{2}}$ is the partial transposition of $\rho$ with respect to
subsystem one and $||\rho||_{1}=\mathrm{Tr}\sqrt{\rho^{\dagger}\rho}$ denotes
the trace norm [75]. Here we compare the entanglement behaviour when the
cavity fields are prepared in equally weighted superpositions of Fock states,
$\ket{\psi_{1,2}}=\frac{1}{\sqrt{2}}(\ket{0}+\ket{N})$. It is useful to first
give the general expression for the time dependent reduced joint cavity state
$\rho_{cav}=\mathrm{Tr}_{s,c}[\rho]$ when the mechanical modes are initially
in the thermal states $\rho_{s,th}\otimes\rho_{c,th}$. Writing the light
fields (which we assume to be pure) in the number basis, $\ket{\psi}_{1}=\sum
a_{m}\ket{m}$ and $\ket{\psi}_{2}=\sum b_{n}\ket{n}$, the reduced density
matrix after the evolution (13) can be found using (141) and the composition
of displacement operators
$\hat{D}(\alpha)\hat{D}(\beta)=e^{\frac{1}{2}(\alpha\beta^{*}-\alpha^{*}\beta)}\hat{D}(\alpha+\beta)$.
This is given by,

$\displaystyle\rho_{cav}\left(t\right)=$
$\displaystyle\sum_{m,n,p,q}M_{mpnq}\ket{m}\\!\\!\bra{p}\ket{n}\\!\\!\bra{q}$
(167)

where,

$\displaystyle M_{mpnq}=$
$\displaystyle\,a_{m}a_{p}^{*}b_{n}b_{q}^{*}e^{-i(\omega_{0}t-A)(m+n-p-q)}e^{iB(m^{2}+n^{2}-p^{2}-q^{2})}e^{iD(mn-
pq)}e^{-\frac{1}{2}|K_{\hat{N}_{c}}|^{2}(\bar{n}_{c}+\frac{1}{2})(n-m+p-q)^{2}}$
(168) $\displaystyle\quad\times
e^{-\frac{1}{2}\left(|K_{\hat{N}_{s}}|^{2}(\bar{n}_{s}+\frac{1}{2})(m+n-p-q)^{2}-(C_{s}K_{\hat{N}_{s}}^{*}-C_{s}^{*}K_{\hat{N}_{s}})(m+n-p-q)\right)},$

and $\bar{n}_{\mu}=1/(e^{\hbar\omega_{\mu}/k_{B}T}-1)$ is the average thermal
occupation of the transformed mechanical modes. Note that as entanglement is
invariant under local operations, we can equivalently consider the state
$\tilde{\rho}_{cav}(t)$ defined by the coefficients,

$\displaystyle\tilde{M}_{mpnq}=$
$\displaystyle\,a_{m}a_{p}^{*}b_{n}b_{q}^{*}e^{iD(mn-
pq)}e^{-\frac{1}{2}|K_{\hat{N}_{c}}|^{2}(\bar{n}_{c}+\frac{1}{2})(n-m+p-q)^{2}}e^{-\frac{1}{2}|K_{\hat{N}_{s}}|^{2}(\bar{n}_{s}+\frac{1}{2})(m+n-p-q)^{2}},$
(169)

where the local actions on the optical states in (13) have been ignored.
Further simplifications can also be made if we restrict our attention to short
times around the optimal measurement time, $t_{q}-\Delta t\leq t\leq
t_{q}+\Delta t$. This is consistent with our expectation of a small
measurement window, and for $\bar{n}_{c}=\bar{n}_{s}=\bar{n}$, an optimal
strategy is to choose the modulation frequency
$\Omega_{n}=(1-1/n)(\omega_{c}+\omega_{s})/2$, with
$t_{q}=qn\pi(1/\omega_{c}+1/\omega_{s})$. The difference between
$|K_{\hat{N}_{s}}(t)|^{2}$ and $|K_{\hat{N}_{c}}(t)|^{2}$ is then proportional
to $\gamma_{g}\Delta t$, and so to lowest order in $\gamma_{g}$, we can assume
that

$\displaystyle\tilde{M}_{mpnq}\approx$
$\displaystyle\,a_{m}a_{p}^{*}b_{n}b_{q}^{*}e^{iD(mn-
pq)}e^{-\kappa((m-p)^{2}+(n-q)^{2})},$ (170)

where
$\kappa=\kappa_{c}^{\mathrm{th}}=|K_{\hat{N}_{c}}|^{2}(\bar{n}+\frac{1}{2})$.

For initial states described by
$a_{m}=b_{m}=(\delta^{0}_{m}+\delta^{N}_{m})/\sqrt{2}$, we can map the partial
transpose of the cavity state,

$\displaystyle\tilde{\rho}^{T_{2}}_{cav}\left(t\right)=$
$\displaystyle\sum_{m,n,p,q\in\\{0,N\\}}\tilde{M}_{mpqn}\ket{m}\\!\\!\bra{p}\ket{n}\\!\\!\bra{q},$
(171)

to its $4\times 4$ matrix representation
$(\tilde{\rho}^{T_{2}}_{cav}(t))_{(m,n),(p,q)}=\tilde{M}_{mpqn}$, where the
indices run over $(j,k)=(0,0),(N,0),(0,N),(N,N)$. The trace norm is then found
from the sum of the eigenvalues of
$(\tilde{\rho}^{T_{2}}_{cav}(t))^{\dagger}(\tilde{\rho}^{T_{2}}_{cav}(t))$,
leading to,

$LN(\rho_{cav}(t))\approx\left\\{\begin{array}[]{ll}0&:\,\sinh(\kappa
N^{2})\geq\sin\left(\frac{DN^{2}}{2}\right)\\\ \log\left[e^{-\kappa
N^{2}}\left(\cosh(\kappa
N^{2})+\sin\left(\frac{DN^{2}}{2}\right)\right)\right]&:\,\rm{else}\end{array}\right..$
(172)

This expression shows that the logarithmic negativity is non-vanishing only in
regions where $\sinh(\kappa N^{2})<\sin\left(\frac{DN^{2}}{2}\right)$. These
define the available measurement windows in which entanglement can be detected
by measurements on the optical fields. In particular, for small arguments of
the $\sinh$ and $\sin$ functions, the condition takes the form $\kappa\leq
D/2$ which coincides with equation (32). A repeat of the analysis leading to
(49) is straightforward, though in this case
$\kappa^{\rm{th}}_{\max}\approx\sinh^{-1}(1)/N^{2}$ and so from (45),

$\delta
t\approx\frac{\sqrt{2\sinh^{-1}(1)}}{\omega_{m}k_{0}N\sqrt{\bar{n}+1/2}}.$
(173)

Note that as one can expect the optimal measurement time now scales as
$t_{q}\sim 1/N^{2}$ (see Appendix A), for both coherent states and Fock state
superpositions the measurement window is inversely proportional to the square
root of $t$, i.e. its behaviour is relatively insensitive to the initial
state.

Figure 4:  The logarithmic negativity of initial cat states
$\ket{\psi_{1,2}}=\frac{1}{\sqrt{2}}(\ket{0}+\ket{10})$ for optomechanical
couplings modulated at the fractional frequency
$\Omega_{q}=(1-1/n)(\omega_{c}+\omega_{s})/2$, $n=1000$ (and phase
$\phi_{1}=\pi/2$) as a function of the time-difference to the measurement time
$t_{1}=n\pi(1/\omega_{c}+1/\omega_{s})$, in units of
$\tau_{m}=2\pi/\omega_{m}$. Here we choose the more optimistic parameters,
$\omega_{m}=2\pi\times 10^{2}$Hz and $m=2.4\times 10^{-8}$kg, leading to
$\gamma_{g}=5\times 10^{-13}$. Plot a) shows $LN$ for different values of the
thermal occupation of the mechanical system $\bar{n}$ and plot b) shows $LN$
for different decoherence rates quantified by the parameter
$\Gamma_{\mathrm{dec}}$ for $\bar{n}=10$.

The above analysis can be repeated to include decoherence via the substitution
$\kappa_{\mu}\rightarrow\kappa_{\mu}=\kappa^{\rm{th}}_{\mu}+\kappa^{\rm{dec}}_{\mu}$
in (169). Just as with the thermal contributions, for small times around
$t_{q}$ the difference between $\kappa_{s}^{\mathrm{dec}}$ and
$\kappa_{c}^{\mathrm{dec}}$ is proportional to $\gamma_{g}$, thus (170) is
again fulfilled in this regime, which immediately recovers equation (172) with
$\kappa=\kappa_{c}^{\mathrm{th}}+\kappa_{c}^{\mathrm{dec}}$.

Figure 4 shows plots of the logarithmic negativity close to $t_{1}$ for the
analytic result above. These agree with numerical evaluations to within a
relative error of less than $10^{-4}$ everywhere up to regions close to the
boundary of the measurement windows (where the numerical errors becomes
significant). For the case of decoherence, we find that the measurement window
closes as $\Gamma_{\mathrm{dec}}$ approaches $1$.

##  Appendix E Displacement and Variance of the Mechanical Modes

The time dependent position operators for the c.o.m. and stretched modes can
be calculated as in Appendix B. We find,

$\displaystyle\hat{x}_{c}(t)=$
$\displaystyle\hat{x}_{\omega_{c}}(t)+2x_{0,c}\Re[e^{-i\omega_{s}t}K_{\hat{N}_{c}}](\hat{N}_{1}-\hat{N}_{2})$
(174) $\displaystyle\hat{x}_{s}(t)=$
$\displaystyle\hat{x}_{\omega_{s}}(t)+2x_{0,s}\left(\Re[e^{-i\omega_{s}t}C_{s}]-\Re[e^{-i\omega_{s}t}K_{\hat{N}_{s}}](\hat{N}_{1}+\hat{N}_{2})\right)$

where,

$\displaystyle\hat{x}_{\omega_{\mu}}(t)=$ $\displaystyle
x_{0,\mu}\left(e^{i\omega_{\mu}t}\hat{b}^{{\dagger}}_{\mu}+e^{-i\omega_{\mu}t}\hat{b}_{\mu}\right)$
(175) $\displaystyle=$
$\displaystyle\hat{x}_{\mu}(0)\cos\omega_{\mu}t+\frac{1}{m\omega_{m}}\hat{p}_{\mu}(0)\sin\omega_{\mu}t$

The stretched mode gives information of the separation of the two mechanical
elements via $\hat{x}_{2}-\hat{x}_{1}=\sqrt{2}\hat{x}_{s}\equiv\delta\hat{x}$.
For thermal states, we have the separation expectation values is simply,

$d+\langle\delta\hat{x}\rangle=d+2\sqrt{2}x_{0,s}\left(\Re[e^{-i\omega_{s}t}C_{s}]-\Re[e^{-i\omega_{s}t}K_{\hat{N}_{s}}](\langle\hat{N}_{1}\rangle+\langle\hat{N}_{2}\rangle)\right)$
(176)

Similarly, the variance in the mode operators is given by,

$\displaystyle(\Delta\hat{x}_{c})^{2}=$
$\displaystyle(\Delta\hat{x}_{\omega_{c}})^{2}+4x_{0,c}^{2}\Re[e^{-i\omega_{s}t}K_{\hat{N}_{c}}]^{2}\left(\Delta(\hat{N}_{1}-\hat{N}_{2})\right)^{2}$
(177) $\displaystyle(\Delta\hat{x}_{s})^{2}=$
$\displaystyle(\Delta\hat{x}_{\omega_{s}})^{2}+4x_{0,s}^{2}\Re[e^{-i\omega_{s}t}K_{\hat{N}_{s}}]^{2}\left(\Delta(\hat{N}_{1}+\hat{N}_{2})\right)^{2}$

Now, for a thermal state,
$(\Delta\hat{x}_{\omega_{\mu}})^{2}=(2\bar{n}_{\mu}+1)x_{0,\mu}^{2}$, thus the
variance in separation is,

$(\Delta(d+\delta\hat{x}))^{2}=2x_{0,s}^{2}\left[1+2\bar{n}_{s}+4\Re[e^{-i\omega_{s}t}K_{\hat{N}_{s}}]^{2}\left(\Delta(\hat{N}_{1}+\hat{N}_{2})\right)^{2}\right].$
(178)


# Exploring Alternatives to Softmax Function

Kunal Banerjee, Vishak Prasad C, Rishi Raj Gupta∗, Karthik Vyas∗, Anushree H∗,
Biswajit Mishra
{kunal.banerjee, vishak.prasad.c<EMAIL_ADDRESS>
Intel Corporation
Bangalore, Karnataka, India

###### Abstract

Softmax function is widely used in artificial neural networks for multiclass
classification, multilabel classification, attention mechanisms, etc. However,
its efficacy is often questioned in literature. The log-softmax loss has been
shown to belong to a more generic class of loss functions, called spherical
family, and its member log-Taylor softmax loss is arguably the best
alternative in this class. In another approach which tries to enhance the
discriminative nature of the softmax function, soft-margin softmax (SM-
softmax) has been proposed to be the most suitable alternative. In this work,
we investigate Taylor softmax, SM-softmax and our proposed SM-Taylor softmax,
an amalgamation of the earlier two functions, as alternatives to softmax
function. Furthermore, we explore the effect of expanding Taylor softmax up to
ten terms (original work proposed expanding only to two terms) along with the
ramifications of considering Taylor softmax to be a finite or infinite series
during backpropagation. Our experiments for the image classification task on
different datasets reveal that there is always a configuration of the SM-
Taylor softmax function that outperforms the normal softmax function and its
other alternatives.

††∗Work done during internship at Intel Corporation.

Keywords: softmax, spherical loss, function approximation, classification.

##  1 Introduction

Softmax function is a popular choice in deep learning classification tasks,
where it typically appears as the last layer. Recently, this function has
found application in other operations as well, such as the attention
mechanisms [8]. However, the softmax function has often been scrutinized in
search of finding a better alternative [9, 1, 6, 5, 4].

Specifically, Vincent et al. explore the spherical loss family in [9] that has
log-softmax loss as one of its members. Brebisson and Vincent further work on
this family of loss functions and propose log-Taylor softmax as a superior
alternative than others, including original log-softmax loss, in [1].

Liu et al. take a different approach to enhance the softmax function by
exploring alternatives which may improve the discriminative property of the
final layer as reported in [6]. The authors propose large-margin softmax (LM-
softmax) that tries to increase inter-class separation and decrease intra-
class separation. LM-softmax is shown to outperform softmax in image
classification task across various datasets. This approach is further
investigated by Liang et al. in [5], where they propose soft-margin softmax
(SM-softmax) that provides a finer control over the inter-class separation
compared to LM-softmax. Consequently, SM-softmax is shown to be a better
alternative than its predecessor LM-softmax [5].

In this work, we explore the various alternatives proposed for softmax
function in the existing literature. Specifically, we focus on two contrasting
approaches based on spherical loss and discriminative property and choose the
best alternative that each has to offer: log-Taylor softmax loss and SM-
softmax, respectively. Moreover, we enhance these functions to investigate
whether further improvements can be achieved. The contributions of this paper
are as follows:

  * •

We propose SM-Taylor softmax – an amalgamation of Taylor softmax and SM-
softmax.

  * •

We explore the effect of expanding Taylor softmax up to ten terms (original
work [1] proposed expanding only to two terms) and we prove higher order even
terms in Taylor’s series are positive definite, as needed in Taylor softmax.

  * •

We explore ramifications of considering Taylor softmax to be a finite or
infinite series during backpropagation.

  * •

We compare the above mentioned variants with Taylor softmax, SM-softmax and
softmax for image classification task.

It may be pertinent to note that we do not explore other alternatives such as,
dropmax [4], because it requires the true labels to be available; however,
such labels may not exist in other tasks where softmax function is used, for
example, attention mechanism [8]. Consequently, dropmax cannot be considered
as a drop-in replacement for softmax universally and hence we discard it.

The paper is organized as follows. Section 2 elaborates on the softmax
function and its several alternatives explored here. Experimental results are
provided in Section 3. Section 4 concludes the paper and shares our plan for
future work.

##  2 Alternatives to softmax

In this section, we provide a brief overview of the softmax function and its
alternatives explored in this work.

###  2.1 Softmax

The softmax function $sm:\mathbb{R}^{K}\to\mathbb{R}^{K}$ is defined by the
formula:

$sm(\mathbf{z})_{i}={\frac{e^{z_{i}}}{\sum_{j=1}^{K}e^{z_{j}}}}{\text{ for
}}i=1,\ldots,K{\text{ and }}\mathbf{z}=(z_{1},\dotsc,z_{K})\in\mathbb{R}^{K}$
(1)

To clarify, the exponential function is applied to each element $z_{i}$ of the
input vector $\mathbf{z}$ and the resulting values are normalized by dividing
by the sum of all the exponentials. The normalization guarantees that the
elements of the output vector $sm(\mathbf{z})$ sum up to 1.

###  2.2 Taylor softmax

The Taylor softmax function as proposed by Vincent et al. [9] uses second
order Taylor series approximation for $e^{z}$ as $1+z+0.5z^{2}$. They then
derive the Taylor softmax as follows:

$Tsm(\mathbf{z})_{i}={}{\frac{1+z_{i}+0.5z_{i}^{2}}{\sum_{j=1}^{K}{1+z_{j}+0.5z_{j}^{2}}}}\text{
for }i=1,\ldots,K{\text{ and
}}\mathbf{z}=(z_{1},\dotsc,z_{K})\in\mathbb{R}^{K}$ (2)

Moreover, the second order approximation of $e^{z}$ as $1+z+0.5z^{2}$ is
positive definite, and hence it is suitable to represent a probability
distribution of classes [1]. Again, it has a minimum value of 0.5, so the
numerator of the equation 2 never becomes zero, that enhances numerical
stability.

We explore higher order Taylor series approximation of $e^{z}$ (as $f^{n}(z)$)
to come up with an $n^{th}$ order Taylor softmax.

${f^{n}(z)}={}\sum_{i=0}^{n}\frac{z^{i}}{i!}$ (3)

Thus, the Taylor softmax for order $n$ is

$Tsm^{n}(\mathbf{z})_{i}={}{\frac{f^{n}(z_{i})}{\sum_{j=1}^{K}{f^{n}(z_{j})}}}\text{
for }i=1,\ldots,K{\text{ and
}}\mathbf{z}=(z_{1},\dotsc,z_{K})\in\mathbb{R}^{K}$ (4)

It is important to note that $f^{n}(z)$ is always positive definite if $n$ is
even. We will prove it by the method of induction.
Base case: We have already shown that $f^{n}(z)$ is positive definite for
$n=2$ in Section 2.
Induction hypothesis: $f^{n}(z)$ is positive definite for $n=2k$.
Induction step: We will prove that it holds for $n=2(k+1)=2k+2$, where $k$ is
an integer starting from 1.
We denote $f^{2k+2}(z)=S(k+1)$, so

$S(k+1)=\sum_{i=0}^{2k+2}\frac{z^{i}}{i!}$ $\displaystyle
S(k+1)=\sum_{i=0}^{2k}\frac{z^{i}}{i!}+\frac{z^{2k+1}}{(2k+1)!}+\frac{z^{2k+2}}{(2k+2)!}$
$\displaystyle\text{Let us consider this series with }p\in\mathbb{R}\text{ and
}p>1$ $\displaystyle
S(k+1,p)=\sum_{i=0}^{2k}\frac{z^{i}}{i!}+\frac{z^{2k+1}}{(2k+1)!}+\frac{z^{2k+2}}{(2k+2)!p}$
$\displaystyle\text{Clearly, }S(k+1)>S(k+1,p)\text{ and }$ $\displaystyle
S(k+1,p)=\sum_{i=0}^{2k-1}\frac{z^{i}}{i!}+\frac{z^{2k}}{(2k)!}(\frac{(4-p)k+2-p}{2(2k+1)}+\frac{(z+k+1)^{2}}{(2k+1)(2k+2)})$
$\displaystyle>\sum_{i=0}^{2k-1}\frac{z^{i}}{i!}+\frac{z^{2k}}{(2k)!}(\frac{(4-p)k+2-p}{2(2k+1)})$
$\displaystyle\text{ if we select }p<\frac{4k+2}{k+1}\text{ then
}\frac{(4-p)k+2-p}{2(2k+1)}>0$ $\displaystyle\text{if we set
}q=\frac{2(2k+1)}{(4-p)k+2-p}\text{ then the expression becomes}$
$\displaystyle
S(k+1,p)>\sum_{i=0}^{2k-1}\frac{z^{i}}{i!}+\frac{z^{2k}}{(2k)!q}=S(k,q)$
$\displaystyle\text{We go further to prove S(k,q) < S(k-1,r)},\text{ it
requires}$ $\displaystyle
q<\frac{4k-2}{k}\Leftrightarrow\frac{2(2k+1)}{(4-p)k+2-p}<\frac{4k-2}{k}$
$\displaystyle\text{which is true if }p<\frac{4k+2}{k+1}$
$\displaystyle\text{Hence, }S(k+1)>S(k,p)>S(k-1,q)...>S(1,t)$
$\displaystyle\text{and }S(1,t)>0\text{ for }t=\frac{5}{3}$
$\displaystyle\text{so, }S(k+1)>0$

$\blacksquare$

The actual back propagation equation for Taylor softmax cross entropy loss
function ($L$) is

$\frac{\partial{L}}{\partial{z}}_{i}={}{\frac{f^{n-1}(z_{i})}{\sum_{j=1}^{k}f^{n}(z_{j})}}-y_{i}\frac{f^{n-1}(z_{i})}{f^{n}(z_{i})}$
(5)

Instead of using equation 5, we used softmax like equation 6 for
backpropagation. For very large $n$ (i.e., as $n$ tends to infinity),
equations 5 and 6 are equivalent, we denote this variation as Taylor_inf. This
equation 5 is corresponding to negative log likelihood loss function of the
Taylor softmax probabilities with a regularizer R(z) defined by equation 7; it
is because of the regularization effect this method performs better.

${\frac{\partial{L}}{\partial{z}}}_{i}={}Tsm^{n}(\mathbf{z})_{i}-y_{i}$ (6)
$R(\mathbf{z})={\log{\frac{Tsm(\mathbf{z})}{sm(\mathbf{z})}}}$ (7)

###  2.3 Soft-margin softmax

Soft-margin (SM) softmax [5] reduces intra-class distances but enhances inter-
class discrimination, by introducing a distance margin into the logits. The
probability distribution for this, as described in [5], is as follows:

$SMsm(\mathbf{z})_{i}=\frac{e^{z_{i}-m}}{{\sum_{j\neq
i}^{K}e^{z_{j}}}+e^{z_{i}-m}}\text{ for }i=1,\ldots,K{\text{ and
}}\mathbf{z}=(z_{1},\dotsc,z_{K})\in\mathbb{R}^{K}$ (8)

###  2.4 SM-Taylor softmax

SM-Taylor softmax uses the same formula as given in equation 8 while using
equation 3, for some given order $n$, to approximate $e^{z}$.

##  3 Experimental Results

In this section, we share our results for image classification task on MNIST,
CIFAR10 and CIFAR100 datasets, where we experiment on the softmax function and
its various alternatives. Note that our goal was not to reach the state of the
art accuracy for each dataset but to compare the influence of each
alternative. Therefore, we restricted ourselves to reasonably sized standard
neural network architectures with no ensembling and no data augmentation.

Table 1: Topologies for different datasets MNIST | CIFAR10 | CIFAR100
---|---|---
{Conv[3x3,64]}x4 | {Conv[3x3,32],BN}x2 | Conv[3x3,384]
MaxPool[2x2,2] | MaxPool[2x2,1] | MaxPool[2x2,1],DO
{Conv[3x3,64]}x3 | Dropout[0.2] | Conv[1x1,384]
MaxPool[2x2,2] | {Conv[3x3,64],BN}x2 | Conv[2x2,384]
{Conv[3x3,64]}x3 | MaxPool[2x2,1] | {Conv[2x2,640]}x2
MaxPool[2x2,2] | Dropout[0.3] | MaxPool[2x2,1],DO
FC[256] | {Conv[3x3,128],BN}x2 | Conv[3x3,640]
FC[10] | MaxPool[2x2,1] | {Conv[2x2,768]}x3
| Dropout[0.4] | Conv[1x1,768]
| FC[128],BN | {Conv[2x2,896]}x2
| Dropout[0.5] | MaxPool[2x2,1],DO
| FC[10] | Conv[3x3,896]
|  | {Conv[2x2,1024]}x2
|  | MaxPool[2x2,1],DO
|  | Conv[1x1,1024]
|  | Conv[2x2,1152]
|  | MaxPool[2x2,1],DO
|  | Conv[1x1,1152]
|  | MaxPool[2x2,1],DO
|  | FC[100]

The topology that we have used for each dataset is given in Table 1. The
topology for MNIST is taken from [5]; we experimented with the topologies for
CIFAR10 and CIFAR100 given in [5] as well to make comparison with the earlier
work easier – however, we could not reproduce the accuracies mentioned in [5]
with the prescribed neural networks. Consequently, we adopted the topology for
CIFAR10 mentioned in [2] and for CIFAR100, we borrow the topology given in
[3]; in both cases, no data augmentation was applied. The abbreviations used
in Table 1 are explained below: (i) Conv[MxN,K] – convolution layer with
kernel size MxN and K output channels, we always use stride of 1 and padding
as “same” for convolutions; (ii) MaxPool[MxN,S] – maxpool layer with kernel
size MxN and stride of S; (iii) FC[K] – fullyconnected layer with K output
channels, an appropriate flatten operation is invoked before the
fullyconnected layer that we have omitted for brevity; (iv) BN – batchnorm
layer with default initialization values; (v) Dropout[R] – dropout layer with
dropout rate R; (vi) DO – dropout layer with rate 0.5, note that CIFAR100
topology uses a uniform rate for all its dropouts; (vii) {layer1[,layer2]}xN –
this combination of layer(s) is repeated N times. In all these topologies, we
replace the final softmax function by each of its alternatives in our
experiments.

Table 2: SM-softmax accuracies for different datasets Dataset | 0 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9
---|---|---|---|---|---|---|---|---|---|---
MNIST | 99.46 | 99.42 | 99.45 | 99.48 | 99.52 | 99.46 | 99.54 | 99.46 | 99.54 | 99.47
CIFAR10 | 87.09 | 87.10 | 87.29 | 87.15 | 87.33 | 87.30 | 87.33 | 87.22 | 87.12 | 87.25
CIFAR100 | 48.28 | 48.03 | 48.06 | 48.11 | 47.82 | 47.68 | 48.95 | 48.03 | 47.96 | 48.02

Table 2 shows the effect of varying soft margin $m$ on accuracy for the three
datasets. We vary $m$ from 0 to 0.9 with a step size of 0.1, as prescribed in
the original work [5]. We note that $m$ set to 0.6 provided the best accuracy
for all the datasets considered, although there are other values which provide
the same best accuracy for MNIST and CIFAR10. Hence, for simplicity, we fix
$m$ to 0.6 for all further experiments.

Table 3: Comparison among softmax and its alternatives Dataset | Variants | Accuracy | 2 | 4 | 6 | 8 | 10
---|---|---|---|---|---|---|---
MNIST | softmax | 99.41 |  |  |  |  |
| Taylor | 99.65 | 99.54 | 99.59 | 99.50 | 99.65 | 99.51
| Taylor_inf | 99.62 | 99.54 | 99.60 | 99.59 | 99.62 | 99.47
| SM-softmax | 99.54 |  |  |  |  |
| SM-Taylor | 99.67 | 99.67 | 99.59 | 99.63 | 99.47 | 99.45
CIFAR10 | softmax | 86.87 |  |  |  |  |
| Taylor | 87.29 | 86.86 | 87.06 | 87.17 | 87.29 | 87.29
| Taylor_inf | 87.46 | 87.46 | 87.37 | 87.34 | 87.00 | 87.38
| SM-softmax | 87.33 |  |  |  |  |
| SM-Taylor | 87.47 | 87.47 | 86.86 | 87.08 | 87.08 | 87.27
CIFAR100 | softmax | 48.57 |  |  |  |  |
| Taylor | 49.94 | 44.70 | 49.24 | 49.94 | 49.84 | 49.04
| Taylor_inf | 49.81 | 44.62 | 47.31 | 49.81 | 46.69 | 45.97
| SM-softmax | 48.95 |  |  |  |  |
| SM-Taylor | 49.95 | 44.77 | 49.95 | 49.56 | 49.69 | 48.11

Table 3 compares the softmax function and its various alternatives with
respect to image classification task for the different datasets. For Taylor
softmax (Taylor) and its variant where we consider equation 6 from Section 2
while doing gradient calculation (Taylor_inf), we look into Taylor series
expansion of orders 2 to 10 with a step size of 2. For SM-Taylor softmax, we
follow the same expansion orders while keeping soft margin fixed at 0.6. For
these three variants, we choose the order which gives the best accuracy and
mention it again in the column labeled “Accuracy”. As can be seen from Table
3, there is always a configuration for SM-Taylor softmax (namely,
$m=0.6,~{}order=2$ for MNIST and CIFAR10, and $m=0.6,~{}order=4$ for CIFAR100)
that outperforms other alternatives.

Figure 1: Plot of training loss vs epochs for MNIST dataset. Figure 2: Plot of
training loss vs epochs for CIFAR10 dataset. Figure 3: Plot of training loss
vs epochs for CIFAR100 dataset.

The plots for training loss vs epochs for MNIST, CIFAR10 and CIFAR100 are
given in Figure 1, Figure 2 and Figure 3, respectively, It may be pertinent to
note that in Figure 1, we see fluctuation in the training loss for the softmax
function, whereas the plot is comparatively smoother for all its alternatives.

##  4 Conclusion and Future Work

Softmax function can be found in almost all modern artificial neural network
models whose applications range from image classification, object detection,
language translation to many more, However, there has been a lot of research
dedicated to finding a better alternative to this popular softmax function.
One approach explores the loss functions belonging to the spherical family,
and proposes log-Taylor softmax loss as arguably the best loss function in
this family [9] Another approach that tries to amplify the discriminative
nature of the softmax function, proposes soft-margin (SM) softmax as the most
appropriate alternative. In this work, we investigate Taylor softmax, soft-
margin softmax and our proposed SM-Taylor softmax as alternatives to softmax
function. Moreover, we study the effect of expanding Taylor softmax up to ten
terms, in contrast to the original work that expanded only to two terms, along
with the ramifications of considering Taylor softmax to be a finite or
infinite series during gradient computation. Through our experiments for the
image classification task on different datasets, we establish that there is
always a configuration of the SM-Taylor softmax function that outperforms the
original softmax function and its other alternatives.

In future, we want to explore bigger models and datasets, especially, the
ILSVRC2012 dataset [7] and its various winning models over the years. Next we
want to explore other tasks where softmax is used, for example, image caption
generation [10] and language translation [8], and check how well do the
softmax alternatives covered in this work perform for the varied tasks.
Ideally, we would like to discover an alternative to softmax that can be
considered as its drop-in replacement irrespective of the task at hand.

## References

  * [1] Alexandre Brébisson and Pascal Vincent  “An Exploration of Softmax Alternatives Belonging to the Spherical Loss Family”  In _ICLR_ , 2016
  * [2] Jason Brownlee  “How to Develop a CNN From Scratch for CIFAR-10 Photo Classification” Accessed: 2020-06-21, https://machinelearningmastery.com/how-to-develop-a-cnn-from-scratch-for-cifar-10-photo-classification/
  * [3] Djork-Arné Clevert, Thomas Unterthiner and Sepp Hochreiter  “Fast and Accurate Deep Network Learning by Exponential Linear Units (ELUs)”  In _ICLR_ , 2016  URL: http://arxiv.org/abs/1511.07289
  * [4] Hae Beom Lee et al.  “DropMax: Adaptive Variationial Softmax”  In _NeurIPS_ , 2018
  * [5] Xuezhi Liang et al.  “Soft-Margin Softmax for Deep Classification”  In _ICONIP_ , 2017, pp. 413–421
  * [6] Weiyang Liu, Yandong Wen, Zhiding Yu and Meng Yang  “Large-Margin Softmax Loss for Convolutional Neural Networks”  In _ICML_ , 2016, pp. 507–516
  * [7] Olga Russakovsky et al.  “ImageNet Large Scale Visual Recognition Challenge”  In _International Journal of Computer Vision (IJCV)_ 115.3, 2015, pp. 211–252  DOI: 10.1007/s11263-015-0816-y
  * [8] Ashish Vaswani et al.  “Attention is All you Need”  In _NeurIPS_ , 2017, pp. 5998–6008
  * [9] Pascal Vincent, Alexandre Brébisson and Xavier Bouthillier  “Efficient Exact Gradient Update for training Deep Networks with Very Large Sparse Targets”  In _NeurIPS_ , 2015, pp. 1108–1116
  * [10] Kelvin Xu et al.  “Show, Attend and Tell: Neural Image Caption Generation with Visual Attention”  In _ICML_ 37, JMLR Workshop and Conference Proceedings, 2015, pp. 2048–2057