# R$\acute{\text{E}}$nyi Differential Privacy in the Shuffle Model: Enhanced
Amplification Bounds

###### Abstract

The shuffle model of Differential Privacy (DP) has gained significant
attention in privacy-preserving data analysis due to its remarkable tradeoff
between privacy and utility. It is characterized by adding a shuffling
procedure after each user’s locally differentially private perturbation, which
leads to a privacy amplification effect, meaning that the privacy guarantee of
a small level of noise, say $\epsilon_{0}$, can be enhanced to
$O(\epsilon_{0}/\sqrt{n})$ (the smaller, the more private) after shuffling all
$n$ users’ perturbed data. Most studies in the shuffle DP focus on proving a
tighter privacy guarantee of privacy amplification. However, the current
results assume that the local privacy budget $\epsilon_{0}$ is within a
limited range. In addition, there remains a gap between the tightest lower
bound and the known upper bound of the privacy amplification. In this work, we
push forward the state-of-the-art by making the following contributions.
Firstly, we present the first asymptotically optimal analysis of
R$\acute{\text{e}}$nyi Differential Privacy (RDP) in the shuffle model without
constraints on $\epsilon_{0}$. Secondly, we introduce hypothesis testing for
privacy amplification through shuffling, offering a distinct analysis
technique and a tighter upper bound. Furthermore, we propose a DP-SGD
algorithm based on RDP. Experiments demonstrate that our approach outperforms
existing methods significantly at the same privacy level. †† This work is
partially support by JST CREST JPMJCR21M2, JSPS KAKENHI Grant Number
JP22H00521, JP22H03595, JP21K19767, JST/NSF Joint Research SICORP JPMJSC2107.

Index Terms—  R$\acute{\text{e}}$nyi Differential Privacy, Shuffle Model,
Privacy Amplification, Hypothesis Testing

##  1 Introduction

R$\acute{\text{e}}$nyi Differential Privacy (RDP) [1, 2] is a variant of
Differential Privacy (DP) [3, 4], which has become the preferred and more
sophisticated framework for ensuring privacy protection. While DP is widely
regarded as the standard framework for privacy protection, RDP offers enhanced
privacy guarantees and finer control over privacy levels by introducing the
R$\acute{\text{e}}$nyi order parameter. Although advanced composition theorems
exist for DP [3] to quantify privacy leakage, it should be noted that these
bounds may not always be tight.

The interest in the shuffle model has been driven by its privacy amplification
effect [5, 6]. This effect becomes notably pronounced because there is often a
trade-off between utility and individual privacy protection through Local
Differential Privacy (LDP) [7]. The primary focus in this area is to achieve
more precise RDP bounds for privacy amplification [8, 9, 10, 11]. Figure 1
illustrates the shuffle model, where the choice between adaptive and non-
adaptive approaches is determined by the query’s dependency on previous
queries.

However, it is worth noting that the existing results not only impose
limitations on the local privacy budget $\epsilon_{0}$ but also leave a
considerable gap to reach the lower bound of RDP. This suggests the need for
further advancements and exploration in this area to bridge the existing gap.

In this work, we provide the first asymptotically optimal analysis for RDP in
the shuffle model without any additional constraints on $\epsilon_{0}$. By
transforming the original problem into a simpler task of measuring the
distance between two multinomial distributions [10], we improve the RDP bound
using hypothesis testing [12, 13]. To provide a clear and convenient
comparison, we outline both previous findings and our own results in Table 1.
Despite the advancements and continuous improvements in RDP bounds under the
shuffle model, previous studies have imposed constraints on the privacy budget
$\epsilon_{0}$ (such as $0<\epsilon_{0}<1/2$ [8]). In contrast, our method
eliminates such restrictions and is suitable for all $\epsilon_{0}>0$.
Furthermore, we introduce a DP-SGD algorithm based on RDP, leveraging the RDP
bounds obtained from our asymptotically optimal analysis. The numerical
results demonstrate that our proposed algorithm significantly surpasses
existing methods in terms of both preserving privacy and maintaining utility.

Fig. 1: The shuffle model with $\epsilon_{0}$-LDP users Methods |  Upper
Bound |  Lower
Bound [9]
---|---|---
Erlingsson et al. (SODA, 2019) [8] | $O(\frac{e^{6\epsilon_{0}}\lambda}{n})$ | $O(\frac{e^{\epsilon_{0}}\lambda}{n})$
Girgis et al. (CCS, 2021) [9] | $O(\frac{e^{2\epsilon_{0}}\lambda}{n})$
Feldman et al. (FOCS, 2022) [10] | $O(\frac{64e^{\epsilon_{0}}\lambda}{n})$
Feldman et al. (SODA, 2023) [11]
This work | $O(\frac{2e^{\epsilon_{0}}\lambda}{n})$ |
Table 1: Privacy amplification of RDP via shuffling with order $\lambda$. For
the sake of brevity, we have not listed the constraints on $\epsilon_{0}$
found in previous studies.

##  2 Preliminaries and Notations

This section gives essential terminology, definitions and properties related
to differential privacy. Due to space limitations, proofs will be included in
the appendix.

######  Definition 1 (Pure Differential Privacy).

A randomized algorithm $\mathcal{R}:\mathcal{D}\rightarrow\mathcal{S}$
satisfies $\epsilon$-DP if for all neighboring datasets
$D_{0},D_{1}\in\mathcal{D}^{n}$,

$\mathbb{P}(\mathcal{R}(D_{0})\in S)\leq
e^{\epsilon}\mathbb{P}(\mathcal{R}(D_{1})\in S).$ (1)

Here $D_{0}$ and $D_{1}$ are considered neighboring if they differ by exactly
one record. If the randomized mechanism is localized, it is called Local
Differential Privacy (LDP).

######  Definition 2.

($R\acute{e}nyi$ Divergence) For two random variables $U$ and $V$, the
$R\acute{e}nyi$ divergence of $U$ and $V$ of order $\lambda>1$ is defined as:

$D^{\lambda}(U\|V)=\frac{1}{\lambda-1}\log\mathop{\mathbb{E}}\limits_{x\sim
V}\left[\left(\frac{U(x)}{V(x)}\right)^{\lambda}\right].$ (2)

Introduced in [1], R$\acute{\text{e}}$nyi differential privacy (RDP) can be
defined based on R$\acute{\text{e}}$nyi divergence.

######  Definition 3 (R$\acute{\text{e}}$nyi Differential Privacy).

A randomized mechanism $\mathcal{R}$ is said to be
$(\lambda,\epsilon(\lambda))$-RDP if for all neighbouring pairs $D_{0}\sim
D_{1}$, it holds that
$D^{\lambda}(\mathcal{R}(D_{0})\|\mathcal{R}(D_{1}))\leq\epsilon.$

Finally, we establish the framework of privacy protection algorithm under
consideration in this paper. The notation $[n]$ represents the set of natural
numbers from $1$ to $n$. For a domain $\mathcal{D}$, let
$\mathcal{R}^{(i)}:\mathcal{S}^{(1)}\times\mathcal{S}^{(2)}\times\cdots\times\mathcal{S}^{(i-1)}\times\mathcal{D}\rightarrow\mathcal{S}^{(i)}$
($i$ $\in[n]$) be a sequence of algorithms such that
$\mathcal{R}^{(i)}(z_{1:i-1},\cdot)$ is an $\epsilon_{0}$-LDP randomizer for
all values of auxiliary inputs
$z_{1:i-1}\in\mathcal{S}^{(1)}\times\mathcal{S}^{(2)}\times\cdots\times\mathcal{S}^{(i-1)}$
where $\mathcal{S}^{(i)}$ is the range space of $\mathcal{R}^{(i)}$. Let
$\mathcal{A}_{R}:\mathcal{D}\rightarrow\mathcal{S}^{(1)}\times\mathcal{S}^{(2)}\times\cdots\times\mathcal{S}^{(n)}$
represent the algorithm applied to the given dataset
$x_{1:n}\in\mathcal{D}^{n}$. The algorithm sequentially computes
$z_{i}=\mathcal{R}^{(i)}(z_{1:i-1},x_{i})$ for $i\in[n]$ and outputs
$z_{1:n}$. We refer to $\mathcal{A}_{R}(\mathcal{D})$ as an $\epsilon_{0}$-LDP
adaptive process. Alternatively, if we first uniformly sample a permutation
$\pi:[n]\rightarrow[n]$, and then sequentially compute
$z_{i}=\mathcal{R}^{(i)}(z_{1:i-1},x_{\pi_{i}})$ for $i\in[n]$, we say it is a
shuffled process and denote it as $\mathcal{A}_{R,S}(\mathcal{D})$. Here,
$\pi_{i}=\pi(i)$ represents the position of $i$ after permutation.

For the sake of brevity and convenience of notation, we omit $D$ and use
$\mathcal{A}_{R}$, $\mathcal{A}_{R,S}$ to represent the adaptive process and
the shuffled adaptive process, respectively.

######  Proposition 1.

(Feldman et al. [10]) For a domain $\mathcal{D}$, let $\mathcal{A}_{R}$ be the
$\epsilon_{0}$-LDP adaptive process and $\mathcal{A}_{R,S}$ be the related
shuffled $\epsilon_{0}$-LDP adaptive process. Assume $X_{0}$ =
$(x^{0}_{1},x_{2},\ldots,x_{n})$ and $X_{1}=(x^{1}_{1},x_{2},\ldots,x_{n})$ be
two neighbouring datasets such that for all $j\neq
1,x_{j}\notin\\{x^{0}_{1},x^{1}_{1}\\}$. Suppose that there exists a positive
value $p\in(0,1]$ such that for all
$i\in[n],x\in\mathcal{D}\backslash\\{x^{0}_{1},x^{1}_{1}\\}$ and
$z_{1:i-1}\in\mathcal{S}^{(1)}\times\mathcal{S}^{(2)}\times\cdots\times\mathcal{S}^{(i-1)}$,
there exists a distribution $LO^{(i)}(z_{1:i-1},x)$ such that

$\begin{split}\mathcal{R}^{(i)}(z_{1:i-1},x)=&\frac{p}{2}\mathcal{R}^{(i)}(z_{1:i-1},x^{0}_{1})+\frac{p}{2}\mathcal{R}^{(i)}(z_{1:i-1},x^{1}_{1})\\\
&+(1-p)LO^{(i)}(z_{1:i-1},x).\end{split}$ (3)

Then there exists a randomized postprocessing algorithm $f$ such that
$\mathcal{A}_{s}(X_{0})$ is distributed identically to
$f(A+\Delta,C-A+1-\Delta)$ and $\mathcal{A}_{R,S}(X_{1})$ is distributed
identically to $f(A+1-\Delta,C-A+\Delta)$, where
$p=e^{-\epsilon_{0}},\Delta\sim
Bern(\frac{e^{\epsilon_{0}}}{e^{\epsilon_{0}}+1}),C\sim Bin(n-1,p),A\sim
Bin(C,1/2)$.

The starting point of this paper is Proposition 1, which transforms the
original problem into a simpler task of analyzing a non-adaptive protocol. It
is worth noting that Proposition 1 mentions the joint distribution of $A$ and
$C$, which corresponds to the multinomial distribution
$Multinom(n-1;p/2,p/2,1-p)$. Here, $A$ and $C-A$ represent the number of 0s
and 1s, respectively.

##  3 Privacy amplification by shuffling based on multinomial distribution

###  3.1 The Exact RDP Bound for the Shuffle Model

In order to provide our tighter exact closed-form bound, we employ a
combination of maximum likelihood and hypothesis testing methods, leveraging
their strengths to derive a more accurate and precise result.

######  Theorem 1.

Let $P=(A,C-A+\Delta)$ and $Q=(A+\Delta,C-A)$, where
$p=e^{-\epsilon_{0}},C\sim Bin(n-1,p),A\sim Bin(C,1/2),\Delta\sim
Bern(\frac{1}{1+p})$. Then

$D^{\lambda}(P\|Q)=\frac{1}{\lambda-1}\log\int_{0}^{1}|h^{\prime}(x)|^{1-\lambda}dx.$

Here, $h(\alpha)$ can be obtained as follows.

$\displaystyle h(\alpha)=1-\alpha-
c_{0}\sum\limits_{v=0}^{n-1}\mathbb{P}\left(B_{1}(v)<A\leq
B_{2}(v)|C=v\right)\mathbb{P}(C=v),$

where

$\alpha(t)=\mathbb{P}(\Delta=0)\mathbb{P}(\frac{A}{C-A+1}>t)+\mathbb{P}(\Delta=1)\mathbb{P}(\frac{A+1}{C-A}>t),$

$g(\alpha)=\inf\limits_{t}\\{t:\alpha(t)\leq\alpha\\},$
$c_{0}=\frac{e^{-\epsilon_{0}}-1}{e^{-\epsilon_{0}}+1}$,
$B_{1}(v)=\frac{g(\alpha)v-1}{g(\alpha)+1},B_{2}(v)=\frac{g(\alpha)v-1}{g(\alpha)+1}+1$,
$\mathbb{P}(C=v)=\binom{n-1}{v}e^{-v\epsilon_{0}}(1-e^{-\epsilon_{0}})^{n-1-v},v=0,1,\cdots,n-1$,
$P(A=k|C=v)=\binom{v}{k}(\frac{1}{2})^{k}$, $k=0,1,\cdots,v.$

Theorem 1 gives exact RDP bound for the shuffled output, however, since RDP is
symmetric, we also need to provide result of $D^{\lambda}(Q\|P)$ in a similar
way. Without causing any ambiguity, we use notations $P$ and $Q$ to represent
the same value in Theorem 1 in the following context. In fact, we can obtain
the R$\acute{\text{e}}$nyi divergence between two multinomial distributions by
numerical computation [10].

######  Corollary 1.

The $\epsilon_{0}$-LDP shuffled adaptive process satisfies
$(\lambda,\max\\{D^{\lambda}(P\|Q),D^{\lambda}(Q\|P)\\}$-RDP, where $P,Q$ and
$D^{\lambda}(\cdot\|\cdot)$ are defined in Theorem 1.

###  3.2 The Asymptotic RDP Bound for the Shuffle Model

Although Corollary 1 provides a RDP bound for the shuffle model, it is
excessively complicated and lacks intuitive understanding. In the following
context, we provide an asymptotic RDP bound. This approach not only simplifies
computations but also facilitates the extension to other privacy definitions
based on the relationship between divergences and Gaussian Differential
Privacy (GDP) [13].

######  Lemma 1.

Assume $\xi$ = $(n_{0},n_{1},n_{2})^{\prime}$ is a random variable which obeys
multinomial distribution with parameters $(n-1;\frac{p}{2},\frac{p}{2},1-p)$,
then $\xi$ approximately follows the multivariate normal distribution
$N(\tilde{\boldsymbol{\mu}},\tilde{\boldsymbol{\Sigma}})$ as
$n\rightarrow\infty$, where
$\tilde{\boldsymbol{\mu}}=\left(\frac{(n-1)p}{2},\frac{(n-1)p}{2},(n-1)(1-p)\right)^{\prime}$
and covariance matrix of $\xi$ is

$\tilde{\boldsymbol{\Sigma}}=(n-1)\left(\begin{array}[]{ccc}\frac{p}{2}(1-\frac{p}{2})&-\frac{p^{2}}{4}&-\frac{p(1-p)}{2}\\\
-\frac{p^{2}}{4}&\frac{p}{2}(1-\frac{p}{2})&-\frac{p(1-p)}{2}\\\
-\frac{p(1-p)}{2}&-\frac{p(1-p)}{2}&p(1-p)\\\ \end{array}\right).$

Lemma 1 provides the asymptotic normality of the multinomial distribution,
which is closely related to GDP. The Berry-Esseen type central limit theorem
[14] ensures that the theoretical convergence rate is at most
$O\left(\frac{1}{\sqrt{n}}\right)$ while the numerical calculation confirms a
convergence rate of approximately $O\left(\frac{1}{n}\right)$, which gives an
upper bound.

######  Theorem 2.

For a domain $\mathcal{D}$, the shuffled $\epsilon_{0}$-LDP adaptive process
approximately satisfies $\frac{2e^{\epsilon_{0}/2}}{\sqrt{n-1}}$-GDP.

According to the fact that $\mu$-GDP implies
$(\lambda,\frac{1}{2}\mu^{2}\lambda)$-RDP [13], we have the asymptotic RDP
bound of the shuffled output.

######  Corollary 2.

For a domain $\mathcal{D}$, the shuffled $\epsilon_{0}$-LDP adaptive process
approximately satisfies $(\lambda,\frac{2e^{\epsilon_{0}}\lambda}{n-1})$-RDP
for any $\lambda\geq 2$.

###  3.3 Comparison of RDP Bounds under the Shuffle Model

Bounds on RDP for privacy amplification via shuffling were initially
introduced by Erlingsson et al. [8]. Girgis et al. [9] gave both an upper
bound and a lower bound of RDP for $\epsilon_{0}\geq 0$ and any integer
$\lambda\geq 2$. That is, the RDP of the shuffled $\epsilon_{0}$-LDP adaptive
process is upper-bounded by

$\epsilon(\lambda)\leq\frac{1}{\lambda-1}\log\left(e^{\lambda^{2}\frac{(e^{\epsilon_{0}}-1)^{2}}{\bar{n}}}+e^{\epsilon_{0}\lambda-\frac{n-1}{8e^{\epsilon_{0}}}}\right),$
(4)

where $\bar{n}=\lfloor\frac{n-1}{2e^{\epsilon_{0}}}\rfloor+1$. And the RDP of
the shuffled $\epsilon_{0}$-LDP adaptive process is lower-bounded by

$\epsilon(\lambda)\geq\frac{1}{\lambda-1}\log\left(1+\frac{\lambda(\lambda-1)}{2}\frac{(e^{\epsilon_{0}}-1)^{2}}{ne^{\epsilon_{0}}}\right).$
(5)

The exponential term $e^{\epsilon_{0}\lambda-\frac{n-1}{8e^{\epsilon_{0}}}}$
in the upper bound comes from the Chernoff bound, it goes to $0$ rapidly as
$n$ increases. If we omit this term, the upper bound is
nearly$\frac{2}{n-1}\frac{\lambda}{\lambda-1}e^{\epsilon_{0}}(e^{\epsilon_{0}}-1)^{2},$
which is worse than our simplified bound in Corollary 2 by a multiplicative
factor of $\frac{\lambda}{\lambda-1}(e^{\epsilon_{0}}-1)^{2}$. Although the
RDP bound has been improved to $O(\frac{64e^{\epsilon_{0}}\lambda}{n})$ [11],
our proposed RDP bound shows a significantly superior performance.

Fig. 2: RDP as a function of $\epsilon_{0}$ for $\lambda=4$ and $n=10^{4}$

Fig. 3: RDP as a function of $\lambda$ for $\epsilon_{0}=2$ and $n=10^{4}$

Figure 3 illustrates that our RDP bound provides a significantly tighter bound
compared to the existing work for a fixed value of $\lambda=4$. Moreover, it
demonstrates that for both $\epsilon_{0}>1$ and $\epsilon_{0}<1$, our RDP
bound and the lower bound are in close proximity. Similarly, Figure 3 shows
that our RDP bound and the lower bound are close for each $\lambda\geq 2$,
with a fixed $\epsilon_{0}=2$.

##  4 Application and Experiments

Stochastic gradient descent (SGD) is a crucial algorithm for empirical risk
minimization (ERM), which aims to minimize a parameterized function given by
$\mathcal{L}(\boldsymbol{\theta})=\sum_{i=1}^{n}\ell(\boldsymbol{\theta},x_{i}),$
where $\boldsymbol{\theta}\in\mathbb{R}^{d}$. Several studies have focused on
a differentially private variant of Stochastic Gradient Descent (SGD) [15, 16,
17]. Additionally, researchers have investigated a deep learning version of
SGD tailored for the popular MNIST handwritten digit dataset [2].

In our experiments, we consider a scenario with $m$ clients, where each client
has $n/m$ samples. Our approach, which incorporates Principle Components
Analysis (PCA), achieves an impressive accuracy of 98.18% after approximately
50 epochs in the non-private scenario. This accuracy result is consistent with
the findings of a standard neural network trained on the same MNIST dataset
[18]. By employing this methodology, we can effectively train a simple
classifier that achieves high accuracy in recognizing handwritten digits from
the MNIST dataset. Details of the experiment setup can be found in Table 2.

To the best of our knowledge, the best RDP bound for the shuffled noisy SGD
with $\epsilon_{0}$-LDP adaptive process is listed in [9, 10], while Figure 3
and 3 show that the privacy bound in this work is tighter. Furthermore, our
technique is based on Laplace mechanism and can be applied to stochastic
gradient descent with batch size $m$. Figure 4 demonstrates that our method
achieves the closest accuracy (96.41%) to the true result (98.18%) when the
same total privacy budget (RDP) is applied. This result is natural because a
better RDP bound allows for more precise local computations, leading to more
accurate results.

Table 2: Experiment setting for the shuffled SGD on the MNIST dataset Parameters/Setting | Value | Explanation
---|---|---
Activation function | ReLU |
Output layer | Softmax |
Loss function | Cross-entropy |
Input layer | 60 variables | 60 PCA components
$C$ | $10$ | Clipping bound
$\epsilon_{0}$ | $[0.1,2]$ | Privacy budget
$\eta$ | $0.05$ | Step size
$m$ | $100$ | Batch size
$n$ | $60,000$ | Sample size
$T$ | $50$ | Epoch count

######  Proposition 2 (Composition theorem of GDP [13]).

The $k$-fold composition of $\mu_{i}$-GDP mechanisms is
$\sqrt{\mu_{1}^{2}+\cdots+\mu_{k}^{2}}$.

######  Proposition 3 (Laplace mechanism [3]).

For a function $g:D\rightarrow\mathbb{R}^{d}$, let $l_{1}$ sensitivity be
defined as $\Delta(g)=\max_{D_{0}\sim D_{1}}\|g(D_{0})-g(D_{1})\|_{1}$, then
for any $\epsilon_{0}>0$, the noisy output $h(D)=g(D)+Lap(\Delta(g)/\epsilon)$
satisfies $\epsilon_{0}$-LDP.

######  Theorem 3.

Algorithm 1 approximately satisfies
$(\lambda,\frac{2Te^{\epsilon_{0}}\lambda}{m-1})$-RDP.

###### Proof.

In each epoch, the algorithm is consisted of two main steps splitting and
shuffling, let

$\displaystyle\mathcal{R}^{(i)}(z_{1:i-1},D_{\pi(i)})$ $\displaystyle=$
$\displaystyle\tilde{\boldsymbol{\theta}}_{i}=\tilde{\boldsymbol{\theta}}_{i-1}(z_{1:i-1})-\eta_{i}(\nabla\ell(\tilde{\boldsymbol{\theta}}_{i-1}(z_{1:i-1}),D_{\pi(i)})+\boldsymbol{b}_{i}),$

then the output of Algorithm 1 can be seen as post processing of the shuffled
$m$ blocks. Since $l_{1}$ sensitivity of each
$\mathcal{R}^{(i)}(z_{1:i-1},\cdot)$ is $\frac{2C}{m}$, then it is
$\epsilon_{0}$-LDP according to Proposition 3. Combined with Theorem 2 and
Proposition 2, Algorithm 1 approximately satisfies
$\frac{2\sqrt{T}}{\sqrt{m-1}}e^{\frac{\epsilon_{0}}{2}}$-GDP. Since $\mu$-GDP
implies $(\lambda,\frac{1}{2}\mu^{2}\lambda)$-RDP, the proof is completed. ∎

Algorithm 1 The Shuffled noisy SGD for $\epsilon_{0}$-LDP users

0: $X=(x_{1},\ldots,x_{n})$, $\mathcal{L}(\boldsymbol{\theta},x)$,
$\boldsymbol{\theta}_{0}$, $\eta$, $T$, $\epsilon_{0}$, $d$, $m$, $C$

0: $\hat{\boldsymbol{\theta}}_{T,m}$

1: Split $[n]$ into $m$ disjoint subsets $S_{1},\cdots,S_{m}$ with equal size
$n/m$

2: Choose arbitrary initial point $\hat{\boldsymbol{\theta}}_{0,m}$

3: for each $t\in[T]$ do

4: $\tilde{\boldsymbol{\theta}}_{0}=\hat{\boldsymbol{\theta}}_{t-1,m}$

5: Choose a random permutation $\pi$ of $[m]$

6: for each $i\in[m]$ do

7: $\boldsymbol{b}_{i}\sim Lap(0,\frac{2C}{\epsilon_{0}m}\boldsymbol{I}_{d})$

8: for each $j\in S_{\pi(i)}$ do

9:
$\boldsymbol{g}_{i}^{j}=\nabla\ell(\tilde{\boldsymbol{\theta}}_{i-1},x_{j})$

10:
$\tilde{\boldsymbol{g}}_{i}^{j}=\boldsymbol{g}_{i}^{j}/\max(1,\|\boldsymbol{g}_{i}^{j}\|_{1}/C)$

11: end for

12:
$\tilde{\boldsymbol{\theta}}_{i}=\tilde{\boldsymbol{\theta}}_{i-1}-\eta(\frac{1}{m}\sum_{j}\tilde{\boldsymbol{g}}_{i}^{j}+\boldsymbol{b}_{i})$

13: end for

14: $\hat{\boldsymbol{\theta}}_{t,m}=\tilde{\boldsymbol{\theta}}_{m}$

15: end for

16: return $\hat{\boldsymbol{\theta}}_{T,m}$

Fig. 4: Comparison of train accuracy at the same privacy level of
$(\lambda,0.008\lambda$)-RDP. Non-Private refers to the standard SGD without
adding noise.

##  5 Conclusion

In conclusion, this study focused on R$\acute{\text{e}}$nyi Differential
Privacy in the shuffle model for privacy protection. We provide an analysis of
RDP without any constraints on the privacy budget for the first time.
Furthermore, a more efficient DP-SGD algorithm based on RDP was proposed,
outperforming existing methods within the same privacy budget.

## References

  * [1] Ilya Mironov,  “Rényi differential privacy,”  in 2017 IEEE 30th Computer Security Foundations Symposium. IEEE, 2017, pp. 263–275.
  * [2] Martin Abadi, Andy Chu, Ian Goodfellow, H Brendan McMahan, Ilya Mironov, Kunal Talwar, and Li Zhang,  “Deep learning with differential privacy,”  in Proceedings of the 2016 ACM SIGSAC Conference on Computer and Communications Security, 2016, pp. 308–318.
  * [3] Cynthia Dwork, Aaron Roth, et al.,  “The algorithmic foundations of differential privacy,”  Foundations and Trends® in Theoretical Computer Science, vol. 9, no. 3–4, pp. 211–407, 2014.
  * [4] Cynthia Dwork, Frank McSherry, Kobbi Nissim, and Adam Smith,  “Calibrating noise to sensitivity in private data analysis,”  Journal of Privacy and Confidentiality, vol. 7, no. 3, pp. 17–51, 2016.
  * [5] Borja Balle, James Bell, Adrià Gascón, and Kobbi Nissim,  “The privacy blanket of the shuffle model,”  in Advances in Cryptology–CRYPTO 2019: 39th Annual International Cryptology Conference, Santa Barbara, CA, USA, August 18–22, 2019, Proceedings, Part II 39. Springer, 2019, pp. 638–667.
  * [6] Albert Cheu, Adam Smith, Jonathan Ullman, David Zeber, and Maxim Zhilyaev,  “Distributed differential privacy via shuffling,”  in Advances in Cryptology–EUROCRYPT 2019: 38th Annual International Conference on the Theory and Applications of Cryptographic Techniques, Darmstadt, Germany, May 19–23, 2019, Proceedings, Part I 38. Springer, 2019, pp. 375–403.
  * [7] Xingxing Xiong, Shubo Liu, Dan Li, Zhaohui Cai, and Xiaoguang Niu,  “A comprehensive survey on local differential privacy,”  Security and Communication Networks, vol. 2020, pp. 1–29, 2020\.
  * [8] Úlfar Erlingsson, Vitaly Feldman, Ilya Mironov, Ananth Raghunathan, Kunal Talwar, and Abhradeep Thakurta,  “Amplification by shuffling: From local to central differential privacy via anonymity,”  in Proceedings of the Thirtieth Annual ACM-SIAM Symposium on Discrete Algorithms. SIAM, 2019, pp. 2468–2479.
  * [9] Antonious M Girgis, Deepesh Data, Suhas Diggavi, Ananda Theertha Suresh, and Peter Kairouz,  “On the renyi differential privacy of the shuffle model,”  in Proceedings of the 2021 ACM SIGSAC Conference on Computer and Communications Security, 2021, pp. 2321–2341.
  * [10] Vitaly Feldman, Audra McMillan, and Kunal Talwar,  “Hiding among the clones: A simple and nearly optimal analysis of privacy amplification by shuffling,”  in 2021 IEEE 62nd Annual Symposium on Foundations of Computer Science. IEEE, 2022, pp. 954–964.
  * [11] Vitaly Feldman, Audra McMillan, and Kunal Talwar,  “Stronger privacy amplification by shuffling for rényi and approximate differential privacy,”  in Proceedings of the 2023 Annual ACM-SIAM Symposium on Discrete Algorithms. SIAM, 2023, pp. 4966–4981.
  * [12] Peter Kairouz, Sewoong Oh, and Pramod Viswanath,  “The composition theorem for differential privacy,”  in International Conference on Machine Learning. PMLR, 2015, pp. 1376–1385.
  * [13] Jinshuo Dong, Aaron Roth, and Weijie J Su,  “Gaussian differential privacy,”  Journal of the Royal Statistical Society Series B: Statistical Methodology, vol. 84, no. 1, pp. 3–37, 2022.
  * [14] Iosif Pinelis and Raymond Molzon,  “Optimal-order bounds on the rate of convergence to normality in the multivariate delta method,”  Electronic Journal of Statistics, vol. 10, no. 1, pp. 1001, 2016\.
  * [15] Raef Bassily, Adam Smith, and Abhradeep Thakurta,  “Private empirical risk minimization: Efficient algorithms and tight error bounds,”  in 2014 IEEE 55th annual symposium on foundations of computer science. IEEE, 2014, pp. 464–473.
  * [16] Shuang Song, Kamalika Chaudhuri, and Anand D Sarwate,  “Stochastic gradient descent with differentially private updates,”  in 2013 IEEE Global Conference on Signal and Information Processing. IEEE, 2013, pp. 245–248.
  * [17] Antonious Girgis, Deepesh Data, Suhas Diggavi, Peter Kairouz, and Ananda Theertha Suresh,  “Shuffled model of differential privacy in federated learning,”  in International Conference on Artificial Intelligence and Statistics. PMLR, 2021, pp. 2521–2529.
  * [18] Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner,  “Gradient-based learning applied to document recognition,”  Proceedings of the IEEE, vol. 86, no. 11, pp. 2278–2324, 1998.

##  A appendix

First, we introduce a powerful tool called $f$-DP based on hypothesis testing
. Readers can refer to the reference [13] for more detailed information.

Let $U$ and $V$ denote the probability distributions of $M(D_{0})$ and
$M(D_{1})$, respectively. We consider a rejection rule $0\leq\phi\leq 1$, with
type I and type II error rates defined as

$\alpha_{\phi}=\mathbb{E}_{U}[\phi],\quad\beta_{\phi}=1-\mathbb{E}_{V}[\phi].$
(6)

######  Definition A.1.

(Trade-off function) For any two probability distributions $U$ and $V$ on the
same space $\Omega$, the trade-off function $T(U,V):[0,1]\rightarrow[0,1]$ is
defined as

$T(U,V)(\alpha)=\inf\\{\beta_{\phi}:\alpha_{\phi}\leq\alpha\\},$

where the infimum is taken over all measurable rejection rules $\phi$, and
$\alpha_{\phi}=\mathbb{E}_{U}(\phi)$ and
$\beta_{\phi}=1-\mathbb{E}_{V}(\phi)$.

######  Definition A.2.

($f$-differential privacy, $f$-DP) Let $f$ be a trade-off function, a
mechanism $M$ is said to be $f$-differentially private if

$T(M(D_{0}),M(D_{1}))(\alpha)\geq f(\alpha),$ (7)

for all neighboring data sets $D_{0}$ and $D_{1}$ and $0\leq\alpha\leq 1$.
Under the assumption of no ambiguity, we assume that the trade-off function is
a function of $\alpha$ and denote the abbreviated form of equation (7) as
$T(M(D_{0}),M(D_{1}))\geq f$.

######  Lemma A.1.

[13] If a mechanism $M$ is $f$-DP and $\inf\\{x\in[0,1]:f(x)=0\\}=1$, then it
is $(\lambda,\ell_{\lambda}^{R\acute{e}nyi}(f))$-RDP for any $\lambda>1$.
Here,
$\ell_{\lambda}^{R\acute{e}nyi}(f)=\frac{1}{\lambda-1}\log\int^{1}_{0}|f^{\prime}(x)|^{1-\lambda}dx$.

Here are some useful facts about $f$-DP and details can be found in [13] .

######  Fact 1.

If a mechanism $M$ is $f$-DP, then its post-processing $Proc\circ M$ is also
$f$-DP.

######  Fact 2.

A $f$-DP mechanism is called $\mu$-GDP if $f$ can be obtained by
$f=T(N(0,1),N(\mu,1))=\Phi(\Phi^{-1}(1-\alpha)-\mu)$. In other words, A
mechanism is $\mu$-GDP if and only if it is $(\epsilon,\delta(\epsilon))$-DP
for all $\epsilon\geq 0$, where

$\delta(\epsilon)=\Phi(-\frac{\epsilon}{\mu}+\frac{\mu}{2})-e^{\epsilon}\Phi(-\frac{\epsilon}{\mu}-\frac{\mu}{2}),$

where $\Phi(\cdot)$ is cumulative distribution function of standard normal
distribution $N(0,1)$.

######  Fact 3.

The $k$-fold composition of $\mu_{i}$-GDP mechanisms is
$\sqrt{\mu_{1}^{2}+\cdots+\mu_{k}^{2}}$.

######  Lemma A.2.

Given observation $(a,b)$, the likelihood function between $H_{0}:P=(A,C-A+1)$
and $H_{1}:Q=(A+1,C-A)$ is
$\Lambda=\frac{\mathcal{L}_{H_{1}}}{\mathcal{L}_{H_{0}}}=\frac{\mathbb{P}(Q=(a,b))}{\mathbb{P}(P=(a,b))}=\frac{a}{b}$,
where $C\sim Bin(n-1,p),A\sim Bin(C,1/2)$.

###### Proof.

Since $C\sim Bin(n-1,p)$, we can obtain that
$\mathbb{P}(C=k)=\binom{n-1}{k}p^{k}(1-p)^{n-1-k},k=0,1,\cdots,n-1$. Then
$\mathbb{P}(Q=(a,b))=\mathbb{P}(C=a+b-1)\mathbb{P}(A=a-1|C=a+b-1)=\mathbb{P}(C=a+b-1)\mathbb{P}(A=b|C=a+b-1)\\\
=\frac{a}{b}\mathbb{P}(P=(a,b)),$ which indicates that
$\Lambda=\frac{\mathcal{L}_{H_{1}}}{\mathcal{L}_{H_{0}}}=\frac{a}{b}$. ∎

######  Lemma A.3.

(Neyman-Pearson lemma) For a hypothesis test between null hypothesis
$H_{0}$:The sample distribution is $P$ and alternative hypothesis $H_{1}$:The
sample distribution is $Q$, let
$\Lambda=\frac{\mathcal{L}_{H_{1}}}{\mathcal{L}_{H_{0}}}$, then the optimal
reject region for $\omega\in\Omega$(sample space) is defined as follows:

$\phi(\omega)=\left\\{\begin{aligned} 1&,&\Lambda>t,\\\ 0&,&\Lambda\leq
t.\end{aligned}\right.$ (8)

######  Lemma A.4.

Assume $\boldsymbol{\mu}_{0}$,$\boldsymbol{\mu}_{1}\in\mathbb{R}^{d}$,
$\boldsymbol{I}_{d}\in\mathbb{R}^{d\times d}$ is an identity matrix and
$\boldsymbol{\Sigma}\in\mathbb{R}^{d\times d}$ is a positive definite matrix,
then

$T(N(\boldsymbol{0},\boldsymbol{I}_{d}),N(\boldsymbol{\mu},\boldsymbol{I}_{d}))=\Phi(\Phi^{-1}(1-\alpha)-\sqrt{\boldsymbol{\mu}^{\prime}\boldsymbol{\mu}}),$

and

$\displaystyle
T(N(\boldsymbol{\mu}_{0},\boldsymbol{\Sigma}),N(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}))$
$\displaystyle=\Phi(\Phi^{-1}(1-\alpha)-\sqrt{(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{0})^{\prime}\boldsymbol{\Sigma}^{-1}(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{0})}),$

where $\alpha\in[0,1]$ represents type I error.

###### Proof.

Notice that likelihood function between
$H_{0}:N(\boldsymbol{0},\boldsymbol{I}_{d})$ and
$H_{1}:N(\boldsymbol{\mu},\boldsymbol{I}_{d})$ is

$\frac{L_{H_{1}}}{L_{H_{0}}}=\frac{\prod\limits_{j=1}^{d}(2\pi)^{\frac{d}{2}}e^{-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^{\prime}(x-\boldsymbol{\mu})}}{\prod\limits_{j=1}^{d}(2\pi)^{\frac{d}{2}}e^{-\frac{1}{2}\boldsymbol{x}^{\prime}\boldsymbol{x}}}=e^{-\frac{1}{2}\boldsymbol{\mu}^{\prime}\boldsymbol{\mu}+\boldsymbol{\mu}^{\prime}\boldsymbol{x}},$
(9)

according to Lemma A.3, the reject rule is defined as follows:

$\phi(\omega)=\left\\{\begin{aligned}
1&,&\boldsymbol{\mu}^{\prime}\boldsymbol{x}>t,\\\
0&,&\boldsymbol{\mu}^{\prime}\boldsymbol{x}\leq t.\end{aligned}\right.$ (10)

Hence, type I error
$\alpha(t)=\mathbb{P}(\boldsymbol{\mu}^{\prime}\boldsymbol{x}>t)$, where
$x\sim N(0,\boldsymbol{I}_{d})$, which indicates that
$\boldsymbol{\mu^{\prime}}\boldsymbol{x}\sim
N(0,\boldsymbol{\mu}^{\prime}\boldsymbol{\mu})$. Therefore,

$\alpha(t)=\mathbb{P}(\frac{\boldsymbol{\mu}^{\prime}\boldsymbol{x}}{\sqrt{\boldsymbol{\mu}^{\prime}\boldsymbol{\mu}}}>\frac{t}{\sqrt{\boldsymbol{\mu}^{\prime}\boldsymbol{\mu}}})=1-\Phi(\frac{t}{\sqrt{\boldsymbol{\mu}^{\prime}\boldsymbol{\mu}}}),$
(11)

which indicates that
$t=\sqrt{\boldsymbol{\mu}^{\prime}\boldsymbol{\mu}}\Phi^{-1}(1-\alpha)$. In
addition, we can calculate type II error

$\displaystyle\beta(t)$
$\displaystyle=\mathbb{P}(\boldsymbol{\mu}^{\prime}(\boldsymbol{x}+\boldsymbol{\mu})\leq
t)$ (12)
$\displaystyle=\mathbb{P}(\frac{\boldsymbol{\mu}^{\prime}\boldsymbol{x}}{\sqrt{{\boldsymbol{\mu}^{\prime}\boldsymbol{\mu}}}}\leq\frac{t}{\sqrt{\boldsymbol{\mu}^{\prime}\boldsymbol{\mu}}}-\sqrt{\boldsymbol{\mu}^{\prime}\boldsymbol{\mu}}),$
$\displaystyle=\Phi(\frac{t}{\sqrt{\boldsymbol{\mu}^{\prime}\boldsymbol{\mu}}}-\sqrt{\boldsymbol{\mu}^{\prime}\boldsymbol{\mu}})=\Phi(\Phi^{-1}(1-\alpha)-\sqrt{\boldsymbol{\mu}^{\prime}\boldsymbol{\mu}}).$

Based on simple calculation, we can find that

$\displaystyle
T(N(\boldsymbol{\mu}_{0},\boldsymbol{\Sigma}),N(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}))$
$\displaystyle=T(N(\boldsymbol{0},\boldsymbol{I}_{d}),N(\boldsymbol{\Sigma}^{-1/2}(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{0}),\boldsymbol{I}_{d})$
$\displaystyle=\Phi(\Phi^{-1}(1-\alpha)-\sqrt{(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{0})^{\prime}\boldsymbol{\Sigma}^{-1}(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{0})}).$

∎

When $d=1$, we have $T(N(0,1),N(\mu,1))=\Phi(\Phi^{-1}(1-\alpha)-\mu)$, which
is consistent with the conclusion in Fact 2.

###  A.1 Proof of Theorem 1

######  Lemma A.5.

Let $P=(A,C-A+1)$ and $Q=(A,C-A)$, where $p=e^{-\epsilon_{0}},C\sim
Bin(n-1,p),A\sim Bin(C,1/2),$ then

$T(P,Q)\geq h.$

Here, $\alpha(t),\beta(t)$ and $h(\alpha)$ can be obtained as follows.

$\alpha(t)=\sum\limits_{v=0}^{n-1}\mathbb{P}(A>\frac{tv+t}{tv-1})\mathbb{P}(C=v),$
$g(\alpha)=\inf\limits_{t}\\{t:\alpha(t)\leq\alpha\\},$
$\beta(t)=h(\alpha)=1-\alpha-\sum\limits_{v=1}^{n-1}\mathbb{P}(A=\lceil\frac{g(\alpha)v-1}{g(\alpha)+1}\rceil)\mathbb{P}(C=v).$

###### Proof.

Consider the reject rule in Lemma A.3, then type I error

$\displaystyle\alpha(t)$ $\displaystyle=\mathbb{P}(\frac{A}{C-A+1}>t)$
$\displaystyle=\sum\limits_{v=0}^{n-1}\mathbb{P}(A>\frac{tv+t}{t+1})\mathbb{P}(C=v),$

and type II error

$\displaystyle\beta(t)$ $\displaystyle=\mathbb{P}(\frac{A+1}{C-A}\leq t)$
$\displaystyle=\sum\limits_{v=0}^{n-1}\mathbb{P}(A\leq\frac{tv-1}{t+1})\mathbb{P}(C=v)$
$\displaystyle=1-\alpha(t)-\sum\limits_{v=0}^{n-1}\mathbb{P}(\frac{tv-1}{t+1}<A\leq\frac{tv-1}{t+1}+1)\mathbb{P}(C=v).$

Since $\alpha(t)$ is non-increasing, then for each $\alpha$, there exists
unique $t$ such that $t=g(\alpha)$, where $g(\alpha)$ =
$\mathop{\inf}\limits_{t}$ $\\{t:\alpha(t)\leq\alpha\\}$. Hence,
$\beta(t)=h(\alpha)=1-\alpha-\sum\limits_{v=0}^{n-1}\mathbb{P}(A=\lceil\frac{g(\alpha)v-1}{g(\alpha)+1}\rceil)\mathbb{P}(C=v)$.
∎

#### Proof of Theorem 1

###### Proof.

Similar in proof of Lemma A.5, we can obtain that

$\displaystyle\alpha(t)$
$\displaystyle=\mathbb{P}(\frac{A+\Delta}{C-A+1-\Delta}>t)$
$\displaystyle=\mathbb{P}(\Delta=0)\mathbb{P}(\frac{A}{C-A+1}>t)+\mathbb{P}(\Delta=1)\mathbb{P}(\frac{A+1}{C-A}>t),$

$\displaystyle\beta(t)$
$\displaystyle=\mathbb{P}(\frac{A+1-\Delta}{C-A+\Delta}\leq t)$
$\displaystyle=\mathbb{P}(\Delta=0)\mathbb{P}(\frac{A+1}{C-A}\leq
t)+\mathbb{P}(\Delta=1)\mathbb{P}(\frac{A}{C-A+1}\leq t).$

Let $g(\alpha)$ = $\inf\limits_{t}$ $\\{t:\alpha(t)\leq\alpha\\}$ and
substitute the probability function of the $A$ and $C$ into the expression, we
can obtain that

$\displaystyle\beta(t)$ $\displaystyle=h(\alpha)$ (13)
$\displaystyle=1-\alpha-
c_{0}\sum\limits_{v=1}^{n-1}\mathbb{P}\left(B_{1}(v)<A\leq
B_{2}(v)|C=v\right)\mathbb{P}(C=v),$

where $c_{0}=\frac{e^{-\epsilon_{0}}-1}{e^{-\epsilon_{0}}+1}$,
$B_{1}(v)=\frac{g(\alpha)v-1}{g(\alpha)+1},B_{2}(v)=\frac{g(\alpha)v-1}{g(\alpha)+1}+1$,
$\mathbb{P}(C=v)=\binom{n-1}{v}e^{-v\epsilon_{0}}(1-e^{-\epsilon_{0}})^{n-1-v},v=0,1,\cdots,n-1$,
$P(A=k|C=v)=\binom{v}{k}(\frac{1}{2})^{k}$, $k=0,1,\cdots,v$. Then it follows
that the mechanism $M$ satisfies $h$-DP. By combining this with Lemma A.1, the
proof is completed. ∎

###  A.2 Proof of Lemma 1

Consider the following two multinomial distributions $P_{1}$ and $Q_{1}$,
where $P_{1}:=(n_{0}+1,n_{1},n_{2})$, $Q_{1}:=(n_{0},n_{1}+1,n_{2})$ and
$(n_{0},n_{1},n_{2})\sim MultiNom(n-1;\frac{p}{2},\frac{p}{2},1-p)$,
$p=e^{-\epsilon_{0}}$. It is easy to check that $T(P,Q)=T(P_{1},Q_{1})$, where
$P$ and $Q$ are distributions defined in Lemma A.2. In addition, if $P_{2}$
and $Q_{2}$ are distributions defined in Theorem 1, then we need to calculate
$T(P_{2},Q_{2})$, where $P_{2}=(n_{0}+1-\Delta,n_{1}+\Delta,n_{2})$ and
$Q_{2}=(n_{0}+\Delta,n_{1}+1-\Delta,n_{2})$, $\Delta\sim Bern(p)$. The
following lemma is a basic property of multinomial distribution.

######  Lemma A.6.

Assume $\xi$ = $(n_{0},n_{1},n_{2})^{\prime}$ is a random variable which obeys
$MultiNom(n-1;\frac{p}{2},\frac{p}{2},1-p)$, then expectation of $\xi$ is
$\tilde{\boldsymbol{\mu}}=(\frac{(n-1)p}{2},\frac{(n-1)p}{2},(n-1)(1-p))^{\prime}$
and covariance matrix of $\xi$ is

$\tilde{\boldsymbol{\Sigma}}=(n-1)\left(\begin{array}[]{ccc}\frac{p}{2}(1-\frac{p}{2})&-\frac{p^{2}}{4}&-\frac{p(1-p)}{2}\\\
-\frac{p^{2}}{4}&\frac{p}{2}(1-\frac{p}{2})&-\frac{p(1-p)}{2}\\\
-\frac{p(1-p)}{2}&-\frac{p(1-p)}{2}&p(1-p)\\\ \end{array}\right)$

This lemma can be obtained from simple calculation directly. Next, we prove
that the multinomial distribution approximately obeys the multivariate normal
distribution.

#### Proof of Lemma 1

###### Proof.

Since $\xi=(n_{0},n_{1},n_{2})^{\prime}$ obeys
$MultiNom(n-1;p_{0},p_{1},p_{2})$ with $p_{0}=p/2,p_{1}=p/2$ and $p_{2}=1-p$,
the characteristic function

$\phi_{\xi}(t_{0},t_{1},t_{2})=E(e^{it_{0}n_{0}+it_{1}n_{1}+it_{2}n_{2}}),$
(14)

then

$\displaystyle\phi_{\xi}(t_{0},t_{1},t_{2})$
$\displaystyle=\sum\limits^{n-1}_{n_{0},n_{1},n_{2}=0}\frac{n!}{n_{0}!n_{1}!n_{2}!}p_{0}^{n_{0}}p_{1}^{n_{1}}p_{2}^{n_{2}}e^{it_{0}n_{0}+it_{1}n_{1}+it_{2}n_{2}}$
(15) $\displaystyle=(p_{0}e^{it_{0}}+p_{1}e^{it_{1}}+p_{2}e^{it_{2}})^{n-1}.$

Let
$\tilde{\xi}=(\frac{n_{0}-(n-1)p_{0}}{\sqrt{n-1}},\frac{n_{1}-(n-1)p_{1}}{\sqrt{n-1}},\frac{n_{2}-(n-1)p_{2}}{\sqrt{n-1}})^{\prime}$,
then we have $\phi_{\tilde{\xi}}(t_{0},t_{1},t_{2})=$

$(p_{0}e^{\frac{it_{0}}{\sqrt{n-1}}}+p_{1}e^{\frac{it_{1}}{\sqrt{n-1}}}+p_{2}e^{\frac{it_{2}}{\sqrt{n-1}}})^{n-1}\\\
e^{-\sqrt{n-1}i(p_{0}t_{0}+p_{1}t_{1}+p_{2}t_{2})}.$

Based on Taylor expansion of $e^{y}=1+y+\frac{y^{2}}{2}+o(y^{2})$ near point
$0$, we can derive that

$\displaystyle
p_{0}e^{\frac{it_{0}}{\sqrt{n-1}}}+p_{1}e^{\frac{it_{1}}{\sqrt{n-1}}}+p_{2}e^{\frac{it_{2}}{\sqrt{n-1}}}$
$\displaystyle=$ $\displaystyle
1+i\frac{a}{\sqrt{n-1}}-\frac{b}{2(n-1)}+o(\frac{1}{n-1}),$

where
$a=\sum\limits_{j=0}^{2}p_{j}t_{j},b=\sum\limits_{j=0}^{2}p_{j}t_{j}^{2}.$ In
addition, since $\log(1+y)=1-\frac{y^{2}}{2}+o(y^{2})$ when $y\rightarrow 0$,
we can obtain that $ln(\psi_{\tilde{\xi}}(t))\rightarrow-\frac{b-a^{2}}{2}$
and

$\displaystyle-\frac{b-a^{2}}{2}$
$\displaystyle=-\frac{1}{2}[\sum\limits_{j=0}^{2}p_{j}(1-p_{j})t_{j}^{2}-\sum\limits_{j=0,j\neq
k}^{2}\sum\limits_{k=0}^{2}p_{j}p_{k}t_{i}t_{k}]$
$\displaystyle=-\frac{1}{2}t^{\prime}\tilde{\boldsymbol{\Sigma}}t,$

which is the logarithm of the characteristic function of multivariate normal
distribution $N(\boldsymbol{0},\tilde{\boldsymbol{\Sigma}})$, where
$\tilde{\boldsymbol{\Sigma}}$ is defined in Lemma A.6. According to the
uniqueness of the characteristic function, $\tilde{\xi}\sim
N(\boldsymbol{0},\tilde{\boldsymbol{\Sigma}})$ when $n$ is sufficiently large,
which indicates that $\xi$ approximately follows
$N(\tilde{\boldsymbol{\mu}},\tilde{\boldsymbol{\Sigma}})$. ∎

###  A.3 Proof of Theorem 2

###### Proof.

According to Proposition 1, the key component is to measure distance between
two distributions $P_{2}$ and $Q_{2}$, where
$P_{2}=(n_{0}+1-\Delta,n_{1}+\Delta,n_{2})$ and
$Q_{2}=(n_{0}+\Delta,n_{1}+1-\Delta)$, $(n_{0},n_{1},n_{2})\sim
MultiNom(n-1;\frac{p}{2},\frac{p}{2},1-p)$, $\Delta\sim Bern(p)$ and
$p=e^{-\epsilon_{0}}$. Let $P=(n_{0}+1,n_{1},n_{2})$ and
$Q=(n_{0},n_{1}+1,n_{2})$, we first prove that

$T(P,Q)=T(Q,P)\geq\Phi\left(\Phi^{-1}(1-\alpha)-\frac{2}{\sqrt{n-1}}e^{\frac{\epsilon_{0}}{2}}\right)$
(16)

approximately. If the inequality (16) holds, then

$T(P_{2},Q_{2})=\mathbb{P}(\Delta=0)T(P,Q)+\mathbb{P}(\Delta=1)T(Q,P),$

which indicates that

$T(P_{2},Q_{2})=(\mathbb{P}(\Delta=0)+\mathbb{P}(\Delta=1))T(P,Q)=T(P,Q)$ (17)

by using the fact that $T(P,Q)=T(Q,P)$.
We will now prove that the inequality (16) holds. Combined with Lemma 1, the
key is to prove the following formula,

$T(N(\tilde{\boldsymbol{\mu}}_{0},\tilde{\boldsymbol{\Sigma}}),N(\tilde{\boldsymbol{\mu}}_{1},\tilde{\boldsymbol{\Sigma}}))\geq\Phi\left(\Phi^{-1}(1-\alpha)-\frac{2}{\sqrt{n-1}}e^{\frac{\epsilon_{0}}{2}}\right),$
(18)

where
$\tilde{\boldsymbol{\mu}}_{0}=\mathbb{E}(P),\tilde{\boldsymbol{\mu}}_{1}=\mathbb{E}(Q)$,
and $\tilde{\boldsymbol{\Sigma}}=Cov(P)=Cov(Q)$. However,
$\tilde{\boldsymbol{\Sigma}}$ is not irreversible because there exists a
linear relationship between sub-vectors: $n_{0}+n_{1}+n_{2}=n-1$. Hence, we
need to combine postprocessing properties of $f$-DP and the fact that sub-
vector of a multivariate normal variable is still a multivariate normal
variable.
Since $P=(n_{0}+1,n_{1},n_{2})^{\prime}$ and
$Q=(n_{0},n_{1}+1,n_{2})^{\prime}$, we can use sub-vector
$P_{2}=(n_{0}+1,n_{1})^{\prime}$ and $Q_{2}=(n_{0},n_{1}+1)^{\prime}$ to
construct $P$ and $Q$ by formula $n_{2}=n-1-n_{0}-n_{1}$. Then,
$P=Proc(P_{2})$ and $Q=Proc(Q_{2})$, which indicates that $T(P,Q)\leq
T(P_{2},Q_{2})$ by Fact 1.
From another perspective, it is obvious that $P_{2}=Proc(P)$ and $Q=Proc(Q)$
if we take $Proc(\cdot)$ as a projection, then $T(P,Q)\geq T(P_{2},Q_{2})$.
Notice that $P_{2}$ obeys $N(\boldsymbol{\mu}_{0},\boldsymbol{\Sigma})$ and
$Q_{2}$ obeys $N(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma})$, where
$\boldsymbol{\mu}_{0}=(\frac{(n-1)p}{2}+1,\frac{(n-1)p}{2})^{\prime},\boldsymbol{\mu}_{1}=(\frac{(n-1)p}{2},\frac{(n-1)p}{2}+1)^{\prime}$
and

$\boldsymbol{\Sigma}=(n-1)\left(\begin{array}[]{cc}\frac{p}{2}(1-\frac{p}{2})&-\frac{p^{2}}{4}\\\
-\frac{p^{2}}{4}&\frac{p}{2}(1-\frac{p}{2})\\\ \end{array}\right).$

After simple algebraic operation, we can obtain that

$\boldsymbol{\Sigma}^{-1}=\frac{1}{n-1}\left(\begin{array}[]{cc}\frac{2-p}{p(1-p)}&\frac{1}{1-p}\\\
\frac{1}{1-p}&\frac{2-p}{p(1-p)}\\\ \end{array}\right),$

and

$\displaystyle(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{0})^{\prime}\boldsymbol{\Sigma}^{-1}(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{0})$
(19)
$\displaystyle=(-1,1)\frac{1}{n-1}\left(\begin{array}[]{cc}\frac{2-p}{p(1-p)}&\frac{1}{1-p}\\\
\frac{1}{1-p}&\frac{2-p}{p(1-p)}\\\
\end{array}\right)\left(\begin{array}[]{c}-1\\\ 1\\\ \end{array}\right)$
$\displaystyle=\frac{4}{(n-1)p}.$

Then $T(P,Q)=\Phi(\Phi^{-1}(1-\alpha)-\frac{2}{\sqrt{(n-1)p}})$, and the
inequality (16) can be obtained immediately by using Lemma A.4. ∎


Classification by Boosting Differences in Input Vectors††thanks: By Ninan
Sajeeth Philip, email: nspp@iucaa.ernet.in–LABEL:lastpage

# Classification by Boosting Differences in Input Vectors††thanks: By Ninan
Sajeeth Philip, email<EMAIL_ADDRESS>
N. S. Philip1, A. Mahabal2, S. Abraham.1, R. Williams2, S.G. Djorgovski2,3, A.
Drake2, C Donalek2 and M. Graham2
1St. Thomas College Kozhencheri  Kerala  India
2 Caltech  1200 E California Bl  Pasadena  CA  91125  USA
3 Distinguished Visiting Professor  King Abdulaziz University  Jeddah  Saudi
Arabia

(2011; Received 27th August 2024)

###### Abstract

There are many occasions when one does not have complete information in order
to classify objects into different classes, and yet it is important to do the
best one can since other decisions depend on that. In astronomy, especially
time-domain astronomy, this situation is common when a transient is detected
and one wishes to determine what it is in order to decide if one must follow
it. We propose to use the Difference Boosting Neural Network (DBNN) which can
boost differences between feature vectors of different objects in order to
differentiate between them. We apply it to the publicly available data of the
Catalina Real-Time Transient Survey (CRTS) and present preliminary results. We
also describe another use with a stellar spectral library to identify spectra
based on a few features. The technique itself is more general and can be
applied to a varied class of problems.

###### keywords:

methods: data analysis, techniques: photometric, techniques: spectroscopic,
stars: general

††volume: 00

##  1 Introduction

It is common to have fragmented and fragmentary evidence when the sources of
information are varied as well as unequal in strength. In time-domain
astronomy, just like in forensic science, one tries to piece together the
information in order to obtain a verdict. An object has been observed that was
either not there before, or was much fainter. The only extra information
available is from sporadic past observations of the same area in archival
surveys that answer questions like: ‘was this seen at such and such radio
frequency?’, ‘was it seen in the Sloan Digital Sky Survey?’, ‘what is the
distance to the nearest galaxy?’ etc. Depending on the location of the object,
and its nature, the resulting information can be sparse and very different
from source to source. The Catalina Real-Time Transient Survey
(CRTS111http://crts.caltech.edu) publishes the transients it finds in real-
time. The events, broadcast as VOEvent packets, are ingested by
Skyalert222http://www.skyalert.org (Williams et al., 2009) where a portfolio
is gathered for each object by annotating the initial information using
programs that query individual surveys and archives to answer, as best as
possible, a pre-determined set of questions as mentioned earlier. If one
thinks of all such bits associated with an individual object as its feature,
the entire input data is just a vector of features (with many, often 50-80
percent, individual features for each vector missing). Making sense of this
dataset is not trivial, but classification is important (Mahabal et al.
(2012), Djorgovski et al. (2011a)). Classifying the features into classes that
are as unambiguous as possible, and with minimal number of false detections
needs a probabilistic approach that can make full use of all available prior
information.

The present work summarizes one possible method to handle such situations
where the sparseness of data is diverse and patchy whereby making it difficult
for existing machine learning techniques to handle them. It uses Bayesian
belief update rule that can be used to progressively update the belief in the
outcome based on plausible evidence. We apply the method to the CRTS
transients to generate quick classification probabilities which can be revised
as more data become available. We also show how it can be generalized to be
used with a stellar spectral library (Valdes et al., 2004) making it
especially useful in future when IFUs make available large number of spectra
simultaneously.

##  2 Catalina Real-Time Survey

The Catalina Real-Time Transient Survey (Drake et al. (2009), Djorgovski et
al. (2011b), Mahabal et al. (2011)) use data from the Catalina Sky Survey
(CSS333http://www.lpl.arizona.edu/css/) for near-Earth objects and potential
planetary hazard asteroids (NEO/PHA), conducted by Edward Beshore, Steve
Larson, and their collaborators at the University of Arizona. CRTS looks for
astrophysical transient and variable objects using real-time processing and
carries out characterization, and distribution of these events, as well as
follow-up observations of selected events.

Optical transients (OTs) are detected as sources displaying significant
changes in brightness, as compared to the baseline comparison images and
catalogs with significantly fainter limiting magnitudes. Data cover time
baselines from 10 min to several years. The detected transients are published
electronically in real time, using a variety of mechanisms. One of the methods
is a Skyalert stream where additional data on the transients is gathered by
harvesting archival datasets as well as information on the proximity of the
source in what is termed as passive follow-up.

It is this comprehensive dataset that we make use of with DBNN in order to try
to predict the nature of the transient. The training set is based upon
classification by human experts. Besides the extremely sparse matrix visible
to DBNN (described in the next section) the humans can use such aids as the
historic light curve of the transient which a human neural network can use to
trivially discriminate between SNe (single hump) versus a CV (multiple
brightenings over a few years) in many cases. Such features are being
incorporated into separate tools. Some derived characteristics from such aids
will be incorporated into DBNN in the future. The different techniques will
also be incorporated in to a fusion network.

##  3 Difference Boosting Neural Network (DBNN)

According to the Bayesian Theorem (von Toussaint, 2011), it is possible to
start with an initial belief about the probability of occurrence of an event
even when there is no compelling evidence and later update this belief
periodically as new evidence is found. In the context of this paper, the
initial belief is called prior and the likelihood for an observation to cause
an event is called evidence. Then according to Bayesian theorem, popularly
called the Bayes rule, the updated belief, or the posterior, is the product of
the prior and the evidence normalized over all possibilities. Mathematically,
this may be written as:

$P(A|B)=\frac{P(B|A)\times P(A)}{\sum_{i}P(B|A_{i})\times P(A_{i})}$ (1)

where $P(A|B)$ is the updated belief that the event A is caused by the
occurrence of event B and $P(B|A)$ is the likelihood that event B may cause A.
P(A) is the prior, the initial belief, which is independent of whether B is
observed or not. One may call it the probability that such an event may be
found even when none of the different types of evidence is seen. The index $i$
is used to normalize all the possible events including $A$ that could have
produced $B$.

The beauty of the Bayesian rule is that it allows sequential updating of the
belief or confidence in an outcome as more and more evidence arrives. Each
time an update is made, the computed new belief becomes the new prior for the
next update, thus essentially making the situation conducive for systems where
one has to deal with a diverse set of inputs and come up with plausible causes
that created them. The second advantage is that, all these estimates are based
on the statistical distribution of the likelihoods and hence the final
evidence is the probability and can be directly used as the confidence one may
have in the prediction of the Bayesian classifier.

We have used a Bayesian Classifier named Difference Boosting Neural Network
(Philip & Joseph, 2000) for this study. The DBNN distributes the evidence in a
feature space so that every feature is associated with its likelihood as
learned through a process called learning (Philip, 2009, 2010). It also
estimates the prior for each evidence in individual cases that maximizes the
prediction accuracy on the data used for training. This data is referred to as
training data. Assuming that we have a large set of training data, it is
possible to have a reasonable estimate of the prior and the likelihood even
for very complex cases. After the training process, the estimated likelihoods
and priors are saved for future use.

Training is usually followed by a testing cycle in which the classifier is
tested with a fresh set of data that are similar but never used in the
training process. This is to test for adequate learning of the classifier, in
which case, the accuracy obtained on the training data and that on the test
data will be comparable. If that is not the case, one has to add the failed
examples also into the training data and update the likelihoods taking them
into consideration. This results in a continued learning process that Bayesian
classifiers can very efficiently handle. The second advantage of this
progressive learning process is that the system asymptotically converges to a
global optimum as the learning progresses. This similarity of Bayesian rule to
human learning made Laplace comment that it is the mathematical equivalent of
common sense.

Preparing the data in a format required for use by the classifier is here
referred to as preprocessing. This is somewhat similar to preprocessing and
data reduction that are familiar to astronomers. The preprocessing ensures
that the input features are ordered in some fashion and an appropriate label
is used to represent the associated class of each example in the data. Usually
the number of input features may be fixed and the names of the classes all
known. But the new situation narrated above has no such restriction. One can
have a dynamic situation where new features are to be incorporated. For
managing this new situation, we have adopted a sequence coding method in which
the presence or absence of a feature is marked by a 1 or 0 in the form of a
chain of binary numbers. New entries are appended to the sequence when new
evidence becomes available, allowing dynamical growth of the feature vector.

Num. Code | Class
---|---
1 | Cataclysmic Variable
2 | Supernova
3 | other
5 | Blazar Outburst
6 | Active Galactic Nucleus Variability
7 | UVCeti Variable
8 | Asteroid
9 | Variable
10 | Mira Variable
11 | High Proper Motion Star
12 | Comet
16 | Nova
Table 1: The 12 different CRTS transient classes identified with numeric
labels and common names.

The preprocessed data may have many missing entries depending on the list of
evidence available per feature vector. The present scheme makes it possible
for us to use any additional information available about an entry to make more
meaningful prediction about its nature. The number of new entries and the
learning based on them dynamically vary as learning progress. However, because
all updating is done statistically, the subtle uncertainties and errors in the
features do not affect the decision making process. The DBNN code is able to
directly train and test on the preprocessed data and we describe the
application of the algorithm to CRTS and stellar spectral studies in the
following section.

Real $\rightarrow$ | 1 | 2 | 3 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 16 | Total
---|---|---|---|---|---|---|---|---|---|---|---|---|---
Predicted |  |  |  |  |  |  |  |  |  |  |  |  |
1 | 273 | 14 | 4 | 2 | 1 | 4 | 3 | 3 | 6 | 1 | 0 | 0 | 311
2 | 15 | 447 | 7 | 3 | 6 | 2 | 4 | 5 | 1 | 4 | 1 | 2 | 497
3 | 1 | 1 | 47 | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 51
5 | 0 | 0 | 0 | 68 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 68
6 | 0 | 1 | 1 | 1 | 144 | 1 | 1 | 3 | 0 | 1 | 0 | 0 | 153
7 | 0 | 0 | 0 | 0 | 0 | 33 | 0 | 0 | 0 | 0 | 0 | 0 | 33
8 | 0 | 0 | 0 | 0 | 0 | 0 | 5 | 0 | 0 | 0 | 0 | 0 | 5
9 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 15 | 0 | 0 | 0 | 0 | 15
10 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 9 | 0 | 0 | 0 | 10
11 | 2 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 51 | 0 | 0 | 53
12 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 5 | 0 | 5
16 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
Total | 291 | 463 | 59 | 75 | 152 | 40 | 13 | 26 | 17 | 57 | 6 | 2 | 1201
Table 2: The confusion matrix of the class labels: humans (Horizontal axis)
against predicted (Vertical).

##  4 Results

The CRTS data that was used for this study has 12 distinct classification
labels enumerated in Table 1 with their numerical class labels followed by the
common names. The input sequence used for this study had 39 input features
(including detection magnitudes, colors from archival photometry, distance to
nearest star/galaxy etc.) of which 50 – 80 percent were missing in some cases.
None of the cases had all 39 inputs. This forms only a subset of the data that
human experts usually use for discrimination. We have not yet translated all
the available information that human experts use into machine recognizable
format (for example, archival light curves). This is something we want to do
in the near future.

Despite these limitations, the classifier agrees with human experts in more
than 90% of the cases as the graduated training progresses and it learns about
the variety in the evidence. This is shown in the confusion matrix Table 2.
The confusion matrix is a convenient representation showing all relevant
information such as how many objects in each class were correctly identified
and into which classes they were incorrectly labeled etc. It may also be used
to have a rough estimate of the number density of the different classes by
normalizing the totals in each class (last row) by the total number of objects
(last column of the last row).

Num. Code | Object Type | Total Objects | Completeness (%) | Contamination (%)
---|---|---|---|---
1 | Cataclysmic Variable | 291 | 93.8 | 12.2
2 | Supernova | 463 | 96.5 | 10.1
3 | Other | 59 | 79.7 | 7.8
5 | Blazar Outburst | 75 | 90.7 | 0.0
6 | Active Galactic Nucleus Variability | 152 | 94.7 | 5.9
7 | UVCeti Variable | 40 | 82.5 | 0.0
8 | Asteroid | 13 | 38.5 | 0.0
9 | Variable | 26 | 57.7 | 0.0
10 | Mira Variable | 17 | 52.9 | 10
11 | High Proper Motion Star | 57 | 89.5 | 3.7
12 | Comet | 6 | 83.3 | 0.0
16 | Nova | 2 | 0.0 | 0.0
Table 3: Completeness and contamination details of predictions by the DBNN classifier based on Table 2 for the different type of objects in CRTS survey are shown. Num. Code | Object Type | Total Objects | Completeness (%) | Contamination (%)
---|---|---|---|---
1 | Cataclysmic Variable | 61 | 63.9 | 35.0
2 | Supernova | 74 | 74.3 | 40.9
3 | Other | 14 | 0.0 | 100
5 | Blazar Outburst | 9 | 55.6 | 16.7
6 | Active Galactic Nucleus Variability | 26 | 76.9 | 48.7
7 | UVCeti Variable | 4 | 25.0 | 0.0
8 | Asteroid | 1 | 0.0 | 0.0
9 | Variable | 3 | 0.0 | 0.0
10 | Mira Variable | 2 | 0.0 | 0.0
11 | High Proper Motion Star | 5 | 0.0 | 0.0
12 | Comet | 1 | 0.0 | 0.0
16 | Nova | 0 | 0.0 | 0.0
Table 4: Classification of transients from one of the test samples. Possible
reasons for relatively lower completeness are detailed in the text.

The performance of the classifier was quantified in terms of contamination and
completeness. The fraction of the total number of objects in a class that were
correctly recognized by the classifier is referred to as completeness and the
fraction that comes as contaminants into a class due to incorrect labelling by
the classifier is referred to as contamination. Both measures are indicative
of how reliable the classifier is. Table 3 shows the completeness and
contamination of the classifier predictions for the CRTS data taking the
verdict of the human expert as reference.

As explained earlier, the classifier evaluates the Bayesian posterior
probability that can be considered synonymous to the true confidence the
classifier predicts. A graphical representation of the completeness and
confidence measures gives a more intuitive picture of the classifier
performance (Fig. 1). It may be noted that for most of the failed cases, the
confidence was low indicating that they are not as reliable as the rest. It is
advisable to use this as a guide to determine how safe it is to rely on low
confidence predictions and to draw a cut-off line which states that anything
below is unreliable. Quantification of this is in progress as more data come
in.

We have also carried out tests with training the method on a smaller subset
and testing on feature vectors never seen before. Due to the sheer variety in
features, and the fact that a large majority of them are missing, the
performance worsens (Table 4). One other important reason for this is that a
few of the classes, like asteroid, comet and high proper motion stars (HPM)
are not static and hence the archives do not have any information about them
(at the discovery location). To incorporate that, additional features will
have to be introduced. Even for many of the other classes there are not enough
examples as yet. Including archival light curves will help greatly improve the
SNe/CV dichotomy. All that is part of on going research, and that should be
kept in mind when judging the current results.

Figure 1: The plots show the performance of the classifier in predicting the
nature of the objects in the CRTS survey. In the plot to the left, black
colour is used to represent the total count of objects in each class while
blue colour is used to represent the number of correctly recognized objects
and red colour is used to indicate the number of contaminants in each class.
The plot to the right shows the same with percent Bayesian confidence and the
counts (in log scale to improve clarity).

As another example, the classifier was used on stellar spectra to classify
them into 98 different classes as given in the Indo-US Stellar
library444http://www.noao.edu/cflib/ (Valdes et al., 2004). The input features
used were the major absorption lines in the spectra and the maximum flux
values from four regions of the continuum in a window of 200Å centered around
wavelengths 3700, 4500, 6300, 8500Å, respectively. These spectra have a
coverage from 3460Å to 9464Å with a few that have missing bands in between. It
was taken with a 0.9m Coud$e^{{}^{\prime}}$ telescope at Kitt Peak National
Observatory in five different grating settings.

The sequence length of the extracted features used by us for this study had 32
most distinctive features in the spectra quantified by their equivalent widths
and the said four flux values. We used 958 examples for training the
classifier and 1104 examples for testing the learning. It was found that 88%
of the classifications were in agreement with the classifications given in the
catalog. The purpose of this work was only to demonstrate that the method
might find some application in spectral analysis, especially for chemical
abundance measures in the stellar atmosphere at a larger scale than what is
possible otherwise.

##  5 Conclusions

We describe a new method to make use of diverse and sparse information from
various sources to classify astronomical data. One practical application we
found is in the case of transients where alerts need to be sent to astronomers
to carry out follow up observations whenever an object that is likely to be of
interest to them is detected. However, it is possible to extend the method to
other applications such as spectroscopic classification where it is difficult
to predefine important absorption/emission lines and we want to cluster them
in a high dimensional space.

## References

  * Djorgovski et al. (2011a) Djorgovski, S. G., Donalek, C., Mahabal, A. A., et al. 2011a, in: Proc. CIDU 2011 Conf., eds. A. Srivasatva, et al., p. 174, NASA Ames Res. Ctr.
  * Djorgovski et al. (2011b) Djorgovski, S. G., Drake, A. J., Mahabal, A. A., et al. 2011b, arXiv:1102.5004
  * Drake et al. (2009) Drake, A. J., Djorgovski, S. G., Mahabal, A., et al. 2009, ApJ, 696, 870
  * Mahabal et al. (2011) Mahabal, A. A., Djorgovski, S. G., Drake, A. J., et al. 2011, BASI, 39, 387
  * Mahabal et al. (2012) Mahabal, A. A., Donalek, C., Djorgovski, S. G., et al. 2012, arXiv:1111.3699, in: Proc. IAU Symp. 285, New Horizons in Time Domain Astronomy, eds. E. Griffin et al., CUP, in press
  * Philip & Joseph (2000) Philip, N. S., Joseph, K. B.,2000, Journal of Intelligent Data Analysis, 4, 463
  * Philip (2009) Philip, N. S., World Conference on Nature and Biologically Inspired Computing,
(NaBIC-2009), IEEE, ISBN: 978-1-4244-5612-3.

  * Philip (2010) Philip, N. S., A Learning Algorithm based on High School Teaching Wisdom,
Paladyn Journal of Behavioral Robotics, 2010, 1(3), 160

  * Valdes et al. (2004) Valdes, F., Gupta, R., Rose, J. A., Singh, H. P., & Bell, D. J. 2004, ApJS, 152, 251
  * von Toussaint (2011) von Toussaint, U. 2011, Reviews of Modern Physics, 83, 943
  * Williams et al. (2009) Williams, R. D., Djorgovski, S. G., Drake, A. J., Graham, M. J., & Mahabal, A. 2009,
Astronomical Data Analysis Software and Systems XVIII, 411, 115

## Acknowledgements

The work at Caltech has been supported in part by the NSF grants AST-0407448,
CNS-0540369, AST-0834235, AST-0909182 and IIS-1118041; the NASA grant
08-AISR08-0085; and by the Ajax and Fishbein Family Foundations. The first
author wishes to thank Prof. Ranjan Gupta for useful discussions and
acknowledges the use of Indo-US Library of Coud$e^{{}^{\prime}}$ Feed Stellar
Spectra.