DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS
DANIELE AGOSTINI AND CARLOS AMÉNDOLA

arXiv:1801.02373v1 [math.AG] 8 Jan 2018

Abstract. We introduce a discrete analogue of the classical multivariate Gaussian distribution.
It is parametrized by the Riemann theta function on the integer lattice. Our real-valued discrete
Gaussian is characterized by the property of maximizing entropy, just as its continuous counterpart.
We capitalize on the theta function representation to derive fundamental properties such as its
characteristic function. Throughout, we exhibit strong connections to the study of abelian varieties
in algebraic geometry.

1. Introduction
Based on the principle of maximum entropy, it is of great interest to find, for a given class of
probability distributions, the ones that maximize entropy. For a random variable X, with density
f (x), its differential entropy as defined by Claude Shannon is
(1.1)

H[X] = E[− log(f (X))]

and it is a measure of the information contained in the distribution. It is well known (see e.g.
[6]) that the maximum entropy distribution for densities supported over a finite set or a bounded
interval in R is given by the corresponding uniform distribution. Since there is no such distribution
supported over the whole real line R, some constraints are needed to make sure there exists a
maximum entropy distribution. If one fixes the mean and the variance, the distribution with
maximum entropy is the Gaussian.
This fact motivates the definition of a discrete Gaussian distribution as the maximum entropy
probability distribution supported over Z with fixed mean and variance.
In their compilation of maximum entropy distributions, Lisman and Van Zuylen in [12] say that
it “cannot be presented in plain terms, because there is no summation procedure for the series
involved”. Kemp in [11] proposes a somewhat explicit form for the discrete Gaussian, and later
Szablowski in [19] observed that the normalization constant can be parametrized in terms of the
Jacobi theta functions.
In this paper, we simplify and extend part of their work to higher dimensions, working over Zg .
Central to our approach is the Riemann theta function [4, Section 8.5],[15, Section II.1]. This is the
holomorphic function

X 1
g
t
t
n τn + n z
(1.2)
θ : C × Hg → C,
θ(z, τ ) =
e
2
g
n∈Z

where e(x) = e2πix and Hg is the Siegel upper half-space, consisting of symmetric g × g complex
matrices whose imaginary part is positive definite.
Using this function, we can define a family of complex probability distributions on Zg , that we
call complex discrete Gaussian distributions.
Definition 1.1 (Complex discrete Gaussian distribution). Fix (z, τ ) ∈ Cg ×Hg such that θ(z, τ ) 6= 0.
We define the discrete Gaussian distribution with parameters (z, τ ) as the complex-valued probability
distribution on Zg given by


1
1 t
t
(1.3)
pθ (n; z, τ ) =
e
n τn + n z
for all n ∈ Zg .
θ(z, τ )
2
Observe that when z and
 τ are purely imaginary (that is, when their real part is zero) the
quantities e 12 nt τ n + nt z are positive real numbers, so that pθ (n; z, τ ) is a standard real-valued
probability density function on the integer lattice Zg . One of our main results, proven in Section 2,
is that this probability distribution is the unique one that maximizes entropy among all those with
same mean and covariance matrix, so that it justifies the name discrete Gaussian.
1

DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS

2

As a consequence, discretizing the kernel of the continuous Gaussian gives indeed a discrete
Gaussian in our sense (maximizing entropy for its own mean and variance). This discretization is
sometimes taken as a definition of discrete Gaussian in the Computer Science literature [1, 16]. Our
observation creates a bridge to papers on the topic and in particular to applications of discrete
Gaussians. In fact, they play a fundamental role in lattice-based cryptography, with one key problem
being how to sample from these distributions efficiently [1, 9, 13]. We believe our approach can
be used to shed light on this problem (e.g. through numerical evaluations of the Riemann theta
function).
The Riemann theta function is a core object in the study of complex abelian varieties: these are
projective varieties that have at the same time the structure of an algebraic group, and they are
ubiquitous in algebraic geometry. We introduce them, their parameter spaces and related concepts
in Section 3. We exhibit several connections between the statistical properties of discrete Gaussians
and the geometry of abelian varieties. For instance, we show in Proposition 3.9 that independent
discrete Gaussians correspond to products of abelian varieties.
In Section 4, we continue exploiting the representation via the Riemann theta function to study
the properties of discrete Gaussians. For example, in Proposition 4.1 we give explicit formulas for
the characteristic function and for the higher moments, answering a question of Szablowski [19]. In
this section we also provide numerical examples where we illustrate simple modeling with discrete
Gaussians and maximum likelihood estimation for them.
Most interestingly, we show in Proposition 3.1 that discrete Gaussians can be interpreted, up
to translation, as points on abelian varieties, and this is our guiding principle. Moreover, this
principle extends in Proposition 3.3 to families of both abelian varieties and discrete Gaussians,
thus obtaining a natural relation between the set of all discrete Gaussian distributions on Zg and
the universal family of abelian varieties over the moduli space Ag [4, Section 8.2].
As a consequence of this principle, every statistical function of the discrete Gaussian that is
invariant under translation, gives a well defined function on abelian varieties. In Section 5, we use
this to define statistical maps of abelian varieties into projective space. These realize abelian varieties
as the moment varieties of the discrete Gaussian distribution. We go on to study the geometry of
these maps and we draw statistical consequences, particularly about moment identifiability. We
show in Theorems 5.9 and 5.14 that complex discrete Gaussians on Z or Z2 with the same parameter
τ are completely determined by the moments up to order three. Moreover, for discrete Gaussians
on Z, we show in Proposition 5.10 how to use the cubic equation of the corresponding elliptic curve
to obtain explicit formulas for the higher moments in terms of the first three.
In summary, in this paper we establish a connection between algebra and statistics that echoes
and complements similar results in Algebraic Statistics. On the one hand, finite discrete exponential
families have been shown to correspond to toric varieties [8]. On the other hand, some continuous
exponential families have been linked to what are called exponential varieties [14]. The present
paper now relates some infinite discrete exponential families to abelian varieties.
2. The discrete Gaussian distribution and the Riemann theta function
We start by extending the definition of the univariate discrete Gaussian in Z to a multivariate
discrete Gaussian in Zg . The multivariate Gaussian distribution given by
1
1
t −1
(2.1)
f (x) = p
e− 2 (x−µ) Σ (x−µ) ,
det(2πΣ)
is a maximum entropy distribution supported over Rg where the mean vector µ ∈ Rg and the
symmetric positive definite g × g covariance matrix Σ are fixed. It is thus natural to attempt the
following definition.
Definition 2.1 (Discrete Gaussian distribution). The g-dimensional discrete Gaussian distribution
is the maximum entropy probability distribution supported over Zg with fixed mean vector µ ∈ Rg
and covariance matrix Σ  0.
A priori, such a distribution may not exist nor be unique. Our main result in this section is that
existence and uniqueness do hold.
To obtain this distribution we will use the classical Riemann theta function; we recall now its
definition. For any positive integer g we denote by Hg the Siegel upper half-space, which consists of

DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS

3

symmetric g × g complex matrices whose imaginary part is positive definite:
Hg = {τ ∈ Cg×g | τ t = τ, Im τ  0}.
Then the Riemann theta function is defined as
(2.2)

g

θ : C × Hg → C,

X

θ(z, τ ) =

n∈Zg


e

1 t
n τ n + nt z
2



where e(x) = e2πix .
Observe that the series is convergent, thanks to the condition Im τ  0: in particular this shows
that the function is holomorphic.
We denote the zero locus of this function by Θg = {(z, τ ) ∈ Cg × Hg | θ(z, τ ) = 0}; this is
sometimes called the universal theta divisor. Now, we define the probability distributions we are
going to work with.
Definition 2.2 (Complex Discrete Gaussian distribution). For (z, τ ) ∈ Cg × Hg \ Θg , we define the
discrete Gaussian distribution with parameters (z, τ ) as the complex-valued probability distribution
on Zg given by


1
1 t
t
(2.3)
pθ (n; z, τ ) =
e
n τn + n z
for all n ∈ Zg .
θ(z, τ )
2
We will denote by X(z,τ ) a random variable with values on Zg and distribution given by pθ (n; z, τ ).
We also give a name to the set of all complex discrete Gaussian distributions. We will see in
Remark 3.4 that this set has a natural structure of complex manifold.
Definition 2.3 (The set Gg ). For a fixed g, we denote by Gg the set of all complex discrete Gaussian
distributions on Zg .
Observe that when z and
 τ are purely imaginary (that is, when their real part is zero) the
quantities e 12 nt τ n + nt z are positive real numbers, so that pθ (n; z, τ ) is a standard real-valued
probability density function on Zg . We will show that this probability distribution is the unique
one that maximizes entropy among all those with same mean and covariance matrix. In particular,
the two definitions of discrete Gaussian agree in this case.
Theorem 2.4. Fix a vector µ ∈ Rg and a positive definite symmetric matrix Σ ∈ Sym2 Rg . Then
there exists a unique distribution supported on Zg with mean vector µ and covariance matrix Σ that
maximizes entropy. This distribution is given by
(2.4)

1
pθ (n; z, τ ) =
e
θ(z, τ )



1 t
n τ n + nt z
2



for some unique z ∈ Cg and τ ∈ Hg purely imaginary.
It is natural to ask how one can compute effectively such z, τ from given µ, Σ. See Remark 2.9
and Subsection 4.1. For a explicit expression of the maximized entropy see Proposition 4.4.
Remark 2.5. When τ has real part zero, it is a positive definite matrix multiplied by the imaginary
unit i. Then the function n 7→ 2πi 12 nt τ n is a negative definite quadratic form. Thus pθ is
a log-concave density (looking at its piecewise linear extension in Rg ), just as its continuous
counterpart.
Remark 2.6. Set g = 1 and
 take two parameters z ∈ C and τ ∈ H1 purely imaginary. Then if we
set q = e (τ ) , λ = e 12 τ + z we can rewrite


n(n−1)
1 2
n τ + nz = λn q 2 .
(2.5)
e
2
This shows that our discrete Gaussian distribution coincides with the one defined by Kemp in [11].
We see now that the density we are proposing is very special in the statistical sense. Indeed, it
belongs to an exponential family.

DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS

4

Exponential families are of fundamental importance in Statistics [5] and we briefly recall their
definition here. Many common examples of distributions belong to an exponential family. An
exponential family on a space X has density
p(η, x) = h(x)ehη,T (x)i−A(η)

(2.6)

for x ∈ X

where η ∈ Rn is a parameter, T : X → Rn is a measurable function, and h·, ·i is the standard scalar
product on Rn . In this case, h is known as the base measure, T as the sufficient statistics and A as
the log-partition function.
Proposition 2.7. Fix z, τ purely imaginary. The density function pθ (n, z, τ ) over Zg in Theorem
2.4 belongs to an exponential family of distributions.
Proof. Indeed, we can rewrite pθ (n, z, τ ) in the exponential family form:
t

(2.7)

1

tτ n

e2πi(n z+ 2 n
pθ (n, z, τ ) =
θ(z, τ )

)

1
t
= ehi(z,τ ),2π(n, 2 nn )i−log(θ(z,τ )) ,

so that we have identity base measure h(n) = 1, sufficient statistics T (n) = 2π(n, 12 nnt ) and
log-partition function A(z, τ ) = log(θ(z, τ ))

Corollary 2.8. There is a bijection between the set of purely imaginary canonical parameters
(z, τ ) ∈ Cg × Hg and the moments (µ, Σ) ∈ Rg × Sym2 Rg with Σ  0. It is given up to scaling by
1
the gradient of the log-partition function ∇A(z, τ ) =
∇θ(z, τ ).
θ(z, τ )
Proof. From classical theory of exponential families [5, 14], there is a bijection between the space of
canonical parameters {η|A(η) < ∞} and the space of sufficient statistics conv(T (x)|x ∈ X ), and it
is given by the gradient of the log-partition function ∇A(η).

Remark 2.9. Furthermore, the inverse map (µ, Σ) → (z, τ ) is known to give the maximum
likelihood estimate (MLE). That is, if we want to estimate (z, τ ) from a sample of an unknown
1
discrete Gaussian, solving the equation ∇A(z, τ ) =
∇θ(z, τ ) is equivalent to solving the
θ(z, τ )
likelihood equations for the sample mean µ and covariance Σ. See Subsection 4.1 for an example.
Finally, exponential families are precisely the right candidates to be maximum entropy distributions
[6]. We present the argument for this in our case.
Corollary 2.10. Fix z, τ purely imaginary and let µ, Σ be the moment vector and the covariance
matrix of the corresponding discrete Gaussian. Among all densities q(n) of support Zg with fixed
mean vector µ and covariance matrix Σ, the distributon pθ (n) = p(n, z, τ ) maximizes entropy.
Proof. We observe that q having the first two matching moments to pθ means that when we sum
q(n) against the linear nt z and against the quadratic form nt τ n, it is equivalent to doing it with
pθ (n). That is,
X
X
X
X
q(n)nt z =
p(n)nt z ,
q(n)nt τ n =
p(n)nt τ n.
n∈Zg

n∈Zg

n∈Zg

n∈Zg

Thus,
X

q(n) log(pθ (n)) =

n∈Zg

X

p(n) log(pθ (n)).

n∈Zg

With this in mind, we proceed:
H(q) = −

X

q(n) log(q(n))

n∈Zg

X
q(n)
−
q(n) log(pθ (n))
pθ (n)
n∈Zg
n∈Zg
X
= −KL(q|pθ ) −
q(n) log(pθ (n))

=−

X

≤−

X

q(n) log

n∈Zg

n∈Zg

= H(pθ )

p(n) log(pθ (n))

DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS

5

where KL(q|pθ ) is the Kullback-Leibler divergence (always non-negative by the classical Jensen’s
inequality). Further, equality holds if and only if q = pθ almost everywhere. Since we are working
with the counting measure over the integer lattice, we have uniqueness.

Combining Corollary 2.8 and Corollary 2.10 we obtain the main result stated in Theorem 2.4.
3. Abelian varieties and discrete Gaussians
The Riemann theta function is a central object in the study of complex abelian varieties. These
are projective varieties that have at the same time the structure of an algebraic group, and they are
fundamental objects in algebraic geometry. Especially important are principally polarized abelian
varieties: these are pairs (A, Θ), where A is an abelian variety and Θ is an ample divisor on A such
that h0 (A, Θ) = 1.
The theta function can be used to build a universal family of these varieties. This is a very
classical construction that we recall briefly here. For more details one can look into [4, Chapter 8].
We should however point out a slight difference between our construction and the classical one:
for some sources, for example [4], an isomorphism of principally polarized abelian varieties (A, Θ)
and (A0 , Θ0 ) is considered to be an isomorphism of groups F : A → A0 such that F −1 (Θ0 ) and Θ
differ by a translation. However, for our purposes we will need to fix the theta divisors, so that, we
will define an isomorphism between (A, Θ) and (A0 , Θ0 ) to be an isomorphism of varieties F : A → A0
such that F −1 (Θ0 ) = Θ.
3.1. Parameter spaces of abelian varieties. Fix an integer g ≥ 1, and for any τ ∈ Hg consider
the subgroup Λτ = {m + τ n | m, n ∈ Zg } ⊆ Cg . Since τ ∈ Hg one can see that Λτ is a lattice, so
that Aτ = Cg /Λτ is a complex torus. Moreover, the theta function is quasiperiodic with respect to
this lattice, meaning that for all n, m ∈ Zg and z ∈ Cg , τ ∈ Hg we have:


1 t
(3.1)
θ(z + m + τ n, τ ) = e − n τ n − nz θ(z, τ ).
2
In particular, the theta divisor Θτ := {z ∈ Cg | θ(z, τ ) = 0} is invariant under Λτ , so that it descends
to a divisor on Aτ that we denote by Θτ again.
Riemann proved that (Aτ , Θτ ) is a principally polarized abelian variety, and moreover he also
showed that every principally polarized abelian variety of dimension g is isomorphic to one of these.
More precisely, for any principally polarized abelian variety (A, Θ), there is a certain τ ∈ Hg and an
isomorphism F : A → Aτ such that F −1 (Θτ ) = Θ. Hence, with this interpretation, the space Hg
becomes a parameter space for all principally polarized abelian varieties.
One can go further and construct a moduli space of these varieties. The space Hg has a natural
action of the symplectic group





0 Ig
0 Ig
t
(3.2)
Sp(2g, Z) = M ∈ M (2g, Z) M
M=
−Ig 0
−Ig 0
defined as follows:

(3.3)

M=

α β
γ δ


∈ Sp(2g, Z), τ ∈ Hg

M τ := (ατ + β)(γτ + δ)−1 .

For M and τ as above, we can define an invertible C-linear map
(3.4)
feM,τ : Cg → Cg ,
z 7→ (γτ + δ)−t z
which in turn induces an isomorphism of abelian varieties fM,τ : Aτ → AM τ . However, this is in
−1
general not an isomorphism of polarized abelian varieties, since the pullback fM,τ
(ΘM τ ) could differ
from Θτ by a translation. To fix this, one defines (see [4, Lemma 8.4.1])
1
1
(3.5)
c1 = diag(γδ t ),
c2 = diag(αβ t ),
cM,τ = (M τ )c1 + c2
2
2
and then considers the affine transformation
(3.6)
FeM,τ : Cg → Cg ,
FeM,τ (z) = feM,τ (z) + cM,τ
which in turn induces a map FM,τ : Aτ → AM τ . One can see that this is again an isomorphism and
−1
moreover FM,τ
(ΘM τ ) = Θτ . Hence, this shows that the two polarized abelian varieties (Aτ , Θτ ) and
(AM τ , ΘM τ ) are isomorphic. Moreover, it turns out that two polarized abelian varieties (Aτ , Θτ ) and

DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS

6

(Aτ 0 , Θτ 0 ) are isomorphic if and only if τ 0 = M τ for a certain M ∈ Sp(2g, Z), and the isomorphism
is then given by FM,τ as before.
Hence, the quotient Ag = Hg /Sp(2g, Z) is a natural parameter space for isomorphism classes of
principally polarized abelian varieties: it is usually called the moduli space of abelian varieties of
dimension g, and it is an algebraic variety itself [4, Remark 8.10.4].
3.2. Universal families of abelian varieties. Out of the previous discussion we can actually
construct also universal families of abelian varieties. The group Zg ⊕ Zg acts on Cg × Hg by
(3.7)

(m, n) · (z, τ ) = (z + m + τ n, τ ).

The quotient Xg = (Cg × Hg )/(Zg ⊕ Zg ) is a complex manifold equipped with a map Xg → Hg such
that the fiber over τ is precisely the abelian variety Aτ . Moreover, by quasiperiodicity (3.1), the
universal theta divisor Θg = {(z, τ ) ∈ Cg × Hg | θ(z, τ ) = 0} passes to the quotient Xg and defines
another divisor in Xg whose restriction to Aτ is precisely the theta divisor Θτ . Hence, the map
Xg → Hg , together with Θg , can be considered as the universal family of abelian varieties over Hg .
We can do something similar also with the moduli space. Indeed, the action of Sp(2g, Z) on Hg ,
extends (see [4, Lemma 8.8.1]) to an action of the group Gg = (Zg ⊕ Zg ) o Sp(2g, Z) on Cg × Hg ,
given by
(3.8)

((m, n), M ) · (z, τ ) = (FeM,τ (z) + m + (M τ )n, M τ ).

The quotient Ug = (Cg × Hg )/Gg by this action can be also seen as a quotient Xg /Sp(2g, Z). In
particular, it has a natural map Ug → Ag . Moreover, the Theta Transformation Formula

[4, Theorem

8.6.1] tells us explicitly how the theta function changes under Gg : for each M = αγ βδ ∈ Sp(2g, Z)
and z ∈ Cg , τ ∈ Hg we have


1 t
−1
t
−1
(3.9)
θ(FM (z, τ ), M τ ) = C(M, τ )e
z (γτ + δ) γz − c1 (γτ + δ) z θ(z, τ )
2
for a certain explicit C(M, τ ) 6= 0. In particular, this tells us that the universal theta divisor Θg
is invariant under the action of Gg , so that it passes to the quotient Ug . In this setting, the map
Ug → Ag , together with Θg is sometimes called the universal family over the moduli space Ag .

3.3. Discrete Gaussian distributions and abelian varieties. In the following, we study further
statistical properties of discrete Gaussians, in light of the connection to abelian varieties.
The key observation is that the quasiperiodicity property has a transparent interpretation in
terms of the discrete Gaussian. Recall that if (z, τ ) ∈ Cg × Hg \ Θg , we denote by X(z,τ ) a random
variable with discrete Gaussian distribution of parameters (z, τ ). We will use the standard notation
X ∼ Y to denote that two random variables have the same distribution.
Proposition 3.1 (Quasiperiodicity). Let (z, τ ) ∈ Cg × Hg \ Θg . Then for every m, n ∈ Zg we have
that
(3.10)

X(z+m+τ n,τ ) ∼ X(z,τ ) − n.

Proof. First we observe that by the quasiperiodicity of the theta function (3.1) we have that
(z + m + τ n, τ ) ∈
/ Θg , so that it makes sense to speak of the variable X(z+m+τ n,τ ) . Now, we fix more
generally (z, τ ) ∈ Cg × Hg and m, n ∈ Zg and for all h ∈ Zg we can compute






1 t
1 t
1 t
t
t
t
t
t
t
e
h τ h + h (z + m + τ n) = e
h τh + h z + h m + h τn = e
h τh + h z + h τn
2
2
2


1
1 t
1 t
1 t
t
t
t
t
(h + n) τ (h + n) + (h + n) z − n τ h − h τ n − n τ n − n z + h τ n
=e
2
2
2
2

 
 

1 t
1
1 t
1 t
t
t
t
t
= e − n τn − n z e
(h + n) τ (h + n) + (h + n) z e − h τ n − n τ h + h τ n .
2
2
2
2

DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS

7

However, since τ is symmetric we see that the last factor in this expression is 1. Notice that this
shows in particular that the theta function is quasiperiodic as in (3.1). Now we can compute

e 21 ht τ h + ht (z + m + τ n)
P(X(z+m+τ n,τ ) = h) =
θ(z + m + τ n, τ )


1 t
e − 2 n τ n − nt z e 12 (h + n)t τ (h + n) + (h + n)t z

=
e − 12 nt τ n − nt z θ(z, τ )

e 12 (h + n)t τ (h + n) + (h + n)t z
= P(X(z,τ ) = h + n)
=
θ(z, τ )
which is precisely what we want.



Remark 3.2. This result tells us that if we fix a parameter τ ∈ Hg , then discrete Gaussian
distributions X(z,τ ) correspond, up to translation, to points in the open subset Aτ \ Θτ ⊆ Aτ .
Hence, we see a direct connection between discrete Gaussians and abelian varieties. More precisely,
the next Proposition 3.3 shows that we can relate the set Gg of all discrete Gaussian distributions
to the universal family Ug .
Proposition 3.3 (Equivalence of discrete Gaussians). Let (z, τ ), (z 0 , τ 0 ) ∈ Cg × Hg \ Θg . Then
X(z,τ ) ∼ X(z 0 ,τ 0 ) if and only if
(3.11)

z = z0 +

1
diag(B) + a,
2

τ = τ0 + B

where B ∈ Sym2 Zg is a symmetric matrix with integer coefficients, diag(B) is the diagonal of B
and a ∈ Zg .
Proof. The two discrete Gaussians X(z,τ ) and X(z 0 ,τ 0 ) have the same distribution if and only if


e 12 nt τ 0 n + nt z 0
e 21 nt τ n + nt z
(3.12)
=
for all n ∈ Zg .
θ(z, τ )
θ(z 0 , τ 0 )
We can rewrite this as


1 t
θ(z, τ )
0
t
0
(3.13)
e
n (τ − τ )n + n (z − z ) =
2
θ(z 0 , τ 0 )

for all n ∈ Zg .

In particular, this shows that the left hand side is independent of n, so that plugging n = 0 gives
that


1 t
0
t
0
(3.14)
e
n (τ − τ )n + n (z − z ) = 1
for all n ∈ Zg .
2
Conversely, if this last condition holds, then it is clear that θ(z, τ ) = θ(z 0 , τ 0 ), so that (3.13) holds
as well. This proves that X(z,τ ) ∼ X(z 0 ,τ 0 ) is equivalent to (3.14). We can rewrite (3.14) as
(3.15)

1 t
n (τ − τ 0 )n + nt (z − z 0 ) ∈ Z
2

Now, let’s set τ − τ 0 = B and z − z 0 =

1
2

for all n ∈ Zg .

diag(B) + a: the condition becomes

1 t
1
n Bn + nt diag(B) + nt a ∈ Z
for all n ∈ Zg .
2
2
We need to prove that this holds if and only if a, B have integer coefficients. Suppose first that a, B
have integer coefficients: then (3.16) holds if and only if
(3.16)

1 2
1
n(n + 1)
n Bii + nBii =
Bii ∈ Z
for all n ∈ Z, i ∈ {1, . . . , g}
2
2
2
but this is true because Bii ∈ Z and n(n + 1) is even.
Conversely, suppose that (3.16) holds and let’s prove that a, B have integer coefficients. Let Ei
denote the i-th vector of the canonical basis. Taking n = Ei and n = −Ei in (3.16) we get that
Bii ∈ Z and ai ∈ Z for all i. To conclude, we just need to check that Bij ∈ Z for all i 6= j and this
follows again from (3.16) taking n = Ei + Ej .

(3.17)

DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS

8

Remark 3.4. We can interpret this via universal families of abelian varieties as follows: the group
Ng = Zg × Sym2 Zg embeds in the group Gg = (Zg ⊕ Zg ) o Sp(2g, Z) via the map



Ig B
(3.18)
Ng ,→ Gg
(a, B) 7→ (a, 0),
.
0 Ig
In particular, the action (3.8) of Gg on Cg × Hg restricts to an action of Ng as follows:


1
(3.19)
(a, B) · (z, τ ) = z + a + diag(B), τ + B .
2
Then, Proposition 3.3 says precisely that the set Gg of all discrete Gaussian distributions is naturally
identified with the quotient ((Cg × Hg ) \ Θg )/Ng . In particular, since the action of Ng is free and
properly discontinuous, we see that Gg has a natural structure of a complex manifold.
Moreover, since Ng is a subgroup of Gg , we have a natural map Gg −→ Ug , so that the set Gg of
discrete Gaussians is an intermediate space between Cg × Hg and the universal family Ug .
3.4. Discrete Gaussians and affine transformations. As in the case of continuous Gaussians,
the class of discrete Gaussians is invariant under affine automorphisms. The case of translations is
covered by Proposition 3.1, so we just need to consider linear automorphisms.
Proposition 3.5. Let X(z,τ ) be a discrete Gaussian on Zg . Then, for any A ∈ GL(g, Z) we have
that
AX(z,τ ) ∼ X(A−1 z,A−t τ A−1 ) .

(3.20)

Proof. We just write down everything explicitly: if n ∈ Zg , then we see that




1 −1 t
1 t −1 −1
−1
−t
t
−t
e
(3.21)
(A n) τ (A n) + (A n)z = e
n (A τ A )n + n (A z)
2
2
In particular, summing over Zg we get that θ(z, τ ) = θ(A−1 z, A−t τ A−1 ). Putting these two together,
it is straightforward get what we want.

Remark 3.6. We can also interpret this result as in Remark 3.4. Indeed, the group GL(g, Z)
embeds in the group Gg = (Zg ⊕ Zg ) o Sp(2g, Z) via the map

 −t

A
0
(3.22)
GL(g, Z) ,→ Gg ,
A 7→ (0, 0),
.
0
A
In particular, the action (3.8) of Gg on Cg × Hg restricts to an action of Ng as follows:

(3.23)
A · (z, τ ) = A−1 z, A−t τ A ,
which is exactly what appears in Proposition 3.5.
Remark 3.7 (Parity). In particular, it follows immediately from Proposition 3.5 that X(−z,τ ) ∼
−X(z,τ ) . This gives a parity property for discrete Gaussians that reflects the analogous property for
the Riemann Theta function.
However, in general it is not true that an arbitrary linear transformation of a discrete Gaussian
is again a discrete Gaussian. Still, we can control the difference, and thanks to Proposition 3.5, it is
enough to do it for projections.
Proposition 3.8. Let X = (X1, X2 ) be a discrete Gaussian random variable on Zg1 × Zg2 with
τ11 τ12
t τ
parameters z = (z1 , z2 ) and τ = τ12
.
22
Then X1 has distribution
(3.24)

P(X1 = n1 ) =

t n ,τ )
θ(z1 , τ11 )θ(z2 + τ12
1 22
P(X(z1 ,τ11 ) = n1 ).
θ(z, τ )

DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS

9

Proof. This is a straightforward computation:
X
P(X1 = n1 ) =
P((X1 , X2 ) = (n1 , n2 ))
n2 ∈Zg2


 
 
X 1
1
τ11 τ12
n1
z
t
t
t
t
(n1 , n2 ) t
+ (n1 , n2 ) 1
=
e
τ
τ
n
z2
θ(z, τ )
2
22
2
12
n2 ∈Zg2


X
1
1 t
1 t
t t
t
t
=
n τ11 n1 + n2 τ12 n1 + n2 τ22 n2 + n1 z1 + n2 z2
e
θ(z, τ )
2 1
2
g
2
n2 ∈Z
 X 


1
1 t
1 t
t
t
t
=
e
n τ11 n1 + n1 z1
e
n τ22 n2 + n2 (z2 + τ12 n1 )
θ(z, τ )
2 1
2 2
n2 ∈Zg2


t n ,τ )
t n ,τ )
θ(z1 , τ11 )θ(z2 + τ12
1 t
θ(z2 + τ12
1 22
1 22
t
e
n1 τ11 n1 + n1 z1 =
P(X(z1 ,τ11 ) = n1 )
=
θ(z, τ )
2
θ(z, τ )

In particular, the components of a discrete Gaussian are not discrete Gaussian themselves, but
we are going to see now that things work well when they are independent.
3.5. Independence for discrete Gaussians. We can characterize independence of joint discrete
Gaussians in a way analogous to continuous Gaussians.
Proposition 3.9. Let X = (X1 , X2 ) be a discrete random variable on Zg1 × Zg2 . Then the following
are equivalent:
(1) X is a discrete Gaussian with X1 , X2 independent. 
(2) X is a discrete Gaussian with parameter τ = τ011 τ022 .
(3) X1 , X2 are independent discrete Gaussians with parameters τ11 and τ22 respectively.
Moreover, if X is a discrete Gaussian with parameters τ and z purely imaginary, these conditions
are equivalent to
(4) X1 , X2 are uncorrelated.
Proof. We first show that
 τ (1)
implies (2). Suppose first that X is a discrete Gaussian with parameters
11 τ12
t
z = (z1 , z2 ) and τ = τ12 τ22 . With the same computations as in Lemma 3.8 we see that


1
1 t
1 t
t
t
t
P((X1 , X2 ) = (n1 , n2 )) =
e
n τ11 n1 + n1 τ12 n2 + n2 τ22 n2 + n1 z1 + n2 z2 ,
θ(z, τ )
2 1
2


t
θ(z2 + τ12 n1 , τ22 )
1 t
t
e
n τ11 n1 + n1 z1 ,
P(X1 = n1 ) =
θ(z, τ )
2 1


θ(z1 + τ12 n2 , τ11 )
1 t
P(X2 = n2 ) =
e
n2 τ22 n2 + nt2 z2 .
θ(z, τ )
2
Hence P((X1 , X2 ) = (n1 , n2 )) = P(X1 = n1 )P(X2 = n2 ) if and only if

t
(3.25)
θ(z, τ ) · e nt1 τ12 n2 = θ(z1 + τ12 n2 , τ11 )θ(z2 + τ12
n1 , τ22 )
Suppose now that X1 , X2 are independent: this means precisely that the relation (3.25) holds for
every n1 , n2 . Choosing n1 = 0 shows that θ(z1 + τ12 n2 , τ11 ) is a constant function of n2 , and choosing
t n , τ ) is a constant function of n . Hence, it must be that e nt τ n
n2 = 0 shows that θ(z2 + τ12
1 22
1
1 12 2
is independent of n1 , n2 , and by choosing one of the two to be zero, we see that e nt1 τ12 n2 = 1.
This is the same as asking that nt1 τ12 n2 ∈ Z for all n1 , n2 , which in turn means that τ12 has integer
coefficients. Then, thanks to Proposition 3.3, we can assume τ12 = 0, which proves (2).
To show that (2) implies (3) observe that if we set τ12 = 0 in (3.25), then the relation reduces to
θ(z, τ ) = θ(z1 , τ11 )θ(z2 , τ22 ), which is true by an explicit computation. Hence, we get that X1 , X2
are independent. To show that they are discrete Gaussians, it is enough to plug τ12 = 0 in Lemma
3.8 and use again that θ(z, τ ) = θ(z1 , τ11 )θ(z2 , τ22 ).
To show that (3) implies (1) it is enough to use again (3.25) and the fact that θ(z, τ ) =
θ(z1 , τ11 )θ(z2 , τ22 ).

DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS

10

To conclude, suppose that X is a discrete Gaussian with parameters z, τ purely imaginary. It
is clear that (1) implies (4). We will show here that (4) implies (1). To do this, let Y1 , Y2 be two
independent discrete Gaussians with purely imaginary parameters which have the same mean vector
and covariance matrix as X1 , X2 respectively. Then we see from (3) that Y = (Y1 , Y2 ) is a discrete
Gaussian distribution with purely imaginary parameters. Moreover, since X1 , X2 are uncorrelated,
Y and X have the same mean and covariance matrix. Then, Theorem 2.4 shows that X, Y have the
same distribution, and we are done.

This result tells us that two joint discrete Gaussians are independent if and only if the matrix
parameter τ is diagonal. This evokes the special property of continuous Gaussians that for the two
components X1 , X2 of a Gaussian vector (X1 , X2 ) one has that X1 , X2 are independent if and only
if their joint covariance Σ12 is zero. Moreover, this is actually true for discrete Gaussians with the
parameters z, τ are purely imaginary.
Remark
 3.10. Moreover, the abelian variety Aτ corresponding to a parameter of the form τ =
τ11 0
0 τ22 is a product of abelian varieties of smaller dimension. More precisely: Aτ = Aτ11 × Aτ22 .
This is suggestive of the relation in Algebraic Statistics between statistical independence of finite
random variables and the Segre embedding of Pn1 × Pn2 .
4. Characteristic function and moments
In this section we use the Riemann theta function to compute various statistical properties of the
discrete Gaussian, such as its characteristic function, its moments, its cumulants and its entropy.
We will use the following notation: if a = (a1 , . . . , ag ) is a multi-index, with ai ∈ N, ai ≥ 0, then
we set |a| = a1 + · · · + ag and a! = a1 ! . . . ag !. Also, if a, b ∈ Ng , we say that a ≤ b if ai ≤ bi for all
i = 1, . . . , g.
a
If v = (v1 , . . . , vg ) ∈ Cg is a vector of complex numbers we set v a = v1a1 . . . vg g . If f (z1 , . . . , zg ) is
a holomorphic function, or a formal power series, we set
(4.1)






∂ a1 f
∂f
∂f a1
∂f ag
∂ ag f
∂f
a
a
.
Dz f =
,...,
,
(Dz f ) =
·...·
· . . . · ag ,
Dz f =
∂z1a1
∂z1
∂zg
∂z1
∂zg
∂zg
Moreover, we denote by Dza (log f ) the higher logarithmic derivatives of f . Recall that they are
computed by taking formal derivatives of log f : for example, if f = f (z) is a function of one variable
we have
f0
f · f 00 − f 02
(4.2)
Dz1 (log f ) = ,
Dz2 (log f ) =
.
f
f2
If X = (X1 , . . . , Xg ) is a random variable (possibly with a complex-valued distribution) we denote
the corresponding means as µi = E[Xi ] and the mean vector as µ = (µ1 , . . . , µg )t . Moreover for
each multi-index a ∈ Ng we define the higher moments as
(4.3)

a

µa [X] = µa1 ,...,ag [X] = E[X1a1 · . . . · Xg g ]

ai ∈ N

and the higher central moments as
(4.4)

ma [X] = µa [X − µ] = E[(X1 − µ1 )a1 · . . . · (Xg − µg )ag ].

The moments are encoded by the characteristic function of X: this is the formal power series
X i|a|
t
(4.5)
φX (v) = E[eiv X ] =
µa [X]v a .
a!
g
a∈N

We also define the cumulants κa [X] through a generating function: more precisely, we formally
define
X κa [X]
t
(4.6)
v a = log E[ev X ].
a!
g
a∈N

We rephrase this by saying that the cumulants correspond to the logarithmic derivatives of the
t
formal power series v 7→ E[ev X ] evaluated at 0:
(4.7)

κa [X] = Dva (log E[ev

tX

])|v=0 .

DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS

11

In the case of the discrete Gaussian distribution we can compute everything explicitly using the
Riemann theta function.
Proposition 4.1. Fix (z, τ ) ∈ Cg × Hg \ Θg . Then the characteristic function of X(z,τ ) is given by

v
θ z + 2π
,τ
1
1 X
1
iv t X(z,τ )
(4.8)
φX(z,τ ) (v) = E[e
]=
=
(Dza θ) · v a .
|a|
θ(z, τ )
θ
a!
(2π)
g
a∈N

In particular, the moments, central moments and cumulants are given by
1 a
1
D θ,
|a|
(2πi) θ z
 
1
1 X a (−1)|b|
ma [X(z,τ ) ] =
(Dz θ)b Dza−b θ,
(2πi)|a| θ 0≤b≤a b
θ|b|

(4.9)

µa [X(z,τ ) ] =

(4.10)
(4.11)

κa [X(z,τ ) ] =

1
Da (log θ).
(2πi)|a| z

Proof. By definition,
(4.12)
(4.13)
(4.14)

 


 

1 t
1 X
1 t
1 t
t
φX(z,τ ) (v) = E[e
]=E e
v X(z,τ )
=
e
n τn + n z e
vn
2π
θ
2
2π
n∈Zg



v

,τ
θ z + 2π
1 X
v 
1 t
t
=
n τn + n z +
=
e
θ
2
2π
θ(z, τ )
n∈Zg
1 X
1
1
=
(Dza θ) · v a
|a|
θ
a!
(2π)
a∈Ng
iv t X(z,τ )

v
where the last equality follows from the Taylor expansion of θ(z + 2π
, τ ) around z. The formula
for µa follows immediately. To get the central moments, we observe that if µ = µ[X(z,τ ) ] then
ma [X(z,τ ) ] = µa [X(z,τ ) − µ] and it is straightforward to compute the characteristic function of
X(z,τ ) − µ:
!
X i|c|
t
(4.15)
µc [X(z,τ ) ]v c
φX(z,τ ) −µ (v) = φµ (−v)φX(z,τ ) (v) = e−iv µ
c!
g
c∈N
!
!
|b|
X i
X i|c|
|b| b b
c
=
(−1) µ v
µc [X(z,τ ) ]v
(4.16)
b!
c!
g
g
c∈N
b∈N
!
X
X i|a|
(4.17)
=
(−1)|b| µb · µc [X(z,τ ) ] v a
b!c!
a∈Ng b+c=a


X i|a|
X a

(4.18)
=
(−1)|b| µb · µc [X(z,τ ) ] v a .
a!
b
g
a∈N

1
c
and µc [X(z,τ ) ] = 1θ (2πi)
|c| Dz θ, we get


 
|b|
X i|a| 1
X
1
a (−1)

φX(z,τ ) −µ (v) =
(Dz θ)b · Dza−b θ v a
|a|
|b|
a!
θ
b
(2πi)
θ
a∈Ng
0≤b≤a

Substituting µ =
(4.19)

0≤b≤a

1 1
θ 2πi Dz θ

which is what we want.
To conclude, we compute the cumulants. By definition we have

 

1 t
1 t
1 X
v t X(z,τ )
t
log E[e
] = log
e
vn e
n τn + z n
θ
2πi
2
n∈Ng

t !


1 X
1 t
1
v
= log
e
n τn + z +
v n = log θ z +
, τ − log θ(z, τ )
θ
2
2πi
2πi
g
n∈N

DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS

12

And taking derivatives with respect to v we get
(4.20)

1
Da (log θ)
(2πi)|a| z

κa [X(z,τ ) ] =


Remark 4.2 (Mean and covariance). In particular, we can write down explicitly the mean vector
and the covariance matrix of X(z,τ ) . Indeed, the first moments always coincide with the first
cumulants, so that
1 ∂ log θ
1 1 ∂θ
=
= κi [X(z,τ ) ].
(4.21)
µi [X(z,τ ) ] =
2πi θ ∂zi
2πi ∂zi
For the covariances we also have that the second central moments coincide with the second cumulants,
so that


1
∂2θ
∂θ ∂θ
1 ∂ 2 log θ
1
(4.22)
mij [X(z,τ ) ] =
θ
−
=
= κij [X(z,τ ) ].
2
2
(2πi) θ
∂zi ∂zj
∂zi ∂zj
(2πi)2 ∂zi ∂zj
We can compare these formulas with those obtained from the theory of exponential families in
Corollary 2.8. In particular, from the expression for the covariances, we get that
1 ∂2θ
∂θ
=
∂τii
4πi ∂zi2

(4.23)

1 ∂2θ
1 ∂θ
=
2 ∂τij
4πi ∂zi ∂zj

and

if i < j

which is precisely the heat equation for the theta function [4, Proposition 8.5.5].
Remark 4.3. As a consequence of Lemma 3.7, we see immediately that the moments µa [X(z,τ ) ], the
central moments ma [X(z,τ ) ] and the cumulants κa [X(z,τ ) ] are even or odd functions of z, depending
on whether |a| is even or odd. In particular, we get immediately that µ[X(0,τ ) ] = 0.
We can also compute the entropy of discrete Gaussian random variables.
Proposition 4.4. Fix (z, τ ) ∈ Cg × Hg \ Θg , and set µ = µ[X(z,τ ) ] and Σ = Cov[X(z,τ ) ]. Then the
entropy of X(z,τ ) is given by


1
1
(4.24)
H[X(z,τ ) ] = log θ −
hz, Dz θi + hτ, Dτ θ + diag(Dτ θ)i
θ
2
= log θ − 2πihz, µi − πihτ, Σ + µµt i

(4.25)

where h·, ·i denotes the standard scalar product.
Proof. We compute:
H[X(z,τ ) ] = −

X
n∈Zg

log

e

1 t
2n τn

+ nt z
θ(z, τ )

!

e

1 t
2n τn

+ nt z
θ(z, τ )

!
=

1 t
2n τn


tz
+
n
=
log θ(z, τ ) − πint τ n − 2πint z
θ(z, τ )
g
n∈Z



!
X
X
1
1 t
1 t
t
t
t
t
log θ(z, τ ) −
πi
(n τ n)e
n τ n + n z + 2πi
(n z)e
n τn + n z
.
θ
2
2
g
g
e

X

n∈Z

n∈Z

Now, differentiating the theta function term by term, we see that


X
∂θ
1 t
t
(4.26)
= 2πi
ni e
n τn + n z
∂zi
2
g
n∈Z

whereas
(4.27)


X
1 t
∂θ
2
t
= πi
ni e
n τn + n z
∂τii
2
g
n∈Z

and



X
1 t
∂θ
t
= πi
2ni nj e
n τn + n z
∂τij
2
g

if i < j.

n∈Z

With this, the first equality in (4.24) follows. For the second, we know from Proposition 4.1 that
1
1 Dτ θ+diag(Dτ θ)
= Σ + µµt . To do this, it is enough to
θ Dz θ = µ, so we just need to check that θ
2
apply Proposition 4.1 together with the heat equation for the theta function (4.23)


DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS

13

Remark 4.5. We should point out that, since θ is a complex function, the logarithm log θ is
defined just up to a constant. However, when z, τ are purely imaginary, we can take the standard
determination of the logarithm on C, which is the one that coincides with the usual logarithm on
the positive real axis. In that case the above formula gives the entropy of the real discrete Gaussian.
4.1. Numerical examples. A nice consequence of the bijection in Corollary 2.8 and the expressions
in Proposition 4.1 is that we can compute numerically the purely imaginary parameters z, τ from
given moments µ, Σ by solving the gradient system
1 1
(4.28)
µ=
Dz θ
2πi θ
1 1
(4.29)
Σ + µµt =
Dτ θ.
πi θ
It provides an advantage with respect to parametrizations like in [11] where the constant is not easy
to compute. In fact, geometers have been interested in computing the Riemann Theta function
and there are several numerical implementations [7]. Thus, in principle, the discrete Gaussian
distribution we propose can be computed effectively for applications.
In fact, in SAGE using the package abelfunctions and with the help of Lynn Chua (who has
maintained and implemented similar functions), we were able to compute with numerical precision
the corresponding z, τ for given choices of µ, Σ. A desirable step following this work would be to
have a fully efficient implementation for all g ≥ 1, for instance as an R package.
For example, let g = 1 and consider the sample with 10 data points (1, 0, 1, −2, 1, 2, 3, −2, 1, −1)
that measure some discrete error. Then the sample mean is µ̂ = 0.4 and the sample standard
deviation is σ̂ = 1.6465. Assuming a univariate discrete Gaussian model and solving numerically
the system (4.28), (4.29) we get ẑ = −0.023i and τ̂ = 0.0587i as the maximum likelihood estimates.
Now suppose we want to know the corresponding purely imaginary parameters z, τ for a “standard”
discrete Gaussian on Z with mean µ = 0 and variance σ 2 = 1. By the parity property of Remark 3.7,
we must have z = 0 (see also Remark 4.3). Solving numerically for τ in (4.29), we get τ ≈ 0.1591549i.
We can generalize this to a standard discrete Gaussian on Zg with mean vector µ = (0, 0, . . . , 0)
and identity covariance matrix Σ = Id. Indeed, the parity property gives again z = (0, 0, . . . , 0).
Moreover, Proposition 3.9 implies that


0.1591549 0
0


..
(4.30)
τ ≈
 i.
.
0
0
0
0 0.1591549
We can also go in the other direction and compute the mean vector and covariance matrix
starting from the parameters z and τ . As an explicit example we take the Fermat quartic C =
{X04 + X14 + X24 = 0} in P2 . This is a smooth curve of genus 3, so that its Jacobian is a principally
polarized abelian variety of dimension g = 3. We can compute the corresponding matrix parameter
τ in MAPLE by using the function periodmatrix. We get that


1+i
−0.5 − 0.5i
−0.5 + 0.5i
(4.31)
τ = −0.5 − 0.5i 0.25 + 1.25i −0.25 − 0.75i .
−0.5 + 0.5i −0.25 − 0.75i −0.25 + 0.75i
Then, we can compute the mean vector and the covariance matrix of the discrete Gaussian on Z3
with parameters z = 0 and τ . Since z = 0, the mean vector is also zero by Remark 4.3. We can
compute the covariance matrix numerically and we get
−3.2885 + 8.1041 · 10−9 i
Σ ≈ 0.1335 · 10−6 + 0.8432 · 10−7 i
5.989 + 0.635 · 10−7 i


0.1335 · 10−6 + 0.8432 · 10−7 i
−6.283 + 5.696i
−6.283 + 5.696i


5.989 + 0.635 · 10−7 i
−6.283 + 5.696i  .
−18.262 + 5.696i

5. Statistical theta functions
Now, let’s fix a parameter τ ∈ Hg and let (Aτ , Θτ ) be the corresponding principally polarized
abelian variety. We have seen before that discrete Gaussians X(z,τ ) correspond, up to translations, to
points in the open abelian variety Aτ \ Θτ . Hence, every function that is invariant under translations
of random variables gives a well-defined function on Aτ \ Θτ . Natural choices of such functions

DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS

14

are the higher central moments and the higher cumulants. We show here that they actually define
meromorphic functions on Aτ .
Proposition 5.1. Fix τ ∈ Hg and a ∈ Ng with |a| > 1. Then the higher central moment ma and
the cumulant κa define meromorphic functions on Aτ with poles only along the theta divisor Θτ of
order at most |a|. Hence they can be seen as global sections in H 0 (Aτ , |a|Θτ ). Moreover, they have
poles precisely of order |a| if and only if (Dz θ)a is not identically zero along Θτ . In particular, this
happens if Θτ is irreducible.
Proof. First, we consider the higher central moments. The explicit expression for ma given in
Proposition 4.1, shows that ma is a meromorphic function on Aτ with poles at most along the theta
divisor Θτ . Moreover, looking again at the explicit form, we see that the summand corresponding
to b has a pole of order at most |b| + 1. The maximum value possible for |b| is |a| and it corresponds
to the unique case b = a. The corresponding summand is
(5.1)

1
(−1)|a|
1
(−1)|a|
(−1)|a| (Dz θ)a
a 0
a
(D
θ)
D
θ
=
(D
θ)
θ
=
,
z
z
z
(2πi)|a| θ|a|+1
(2πi)|a| θ|a|+1
(2πi)|a| θ|a|

which has poles only along Θτ of order at most |a|. This shows immediately that ma has poles only
along Θτ of order at most |a|. For the cumulants κa one can follow a similar reasoning by writing
an explicit expression for Dza log θ, or one can use the fact that the cumulants are determined by
the central moments through some explicit expressions [2, eqn (7)].
Now we ask ourselves when do ma and κa attain poles of order |a|. Consider again the central
moments first. Looking at the expression of Proposition 4.1, we consider the summands with
|b| = |a| − 1: these are of the form a − Ei , where i is an index such that ai > 0. Taking the sum
over these we get
(−1)|a|−1 1 X
(−1)|a|−1 1 X
a−Ei Ei
θ
=
D
(5.2)
a
(D
θ)
ai (Dz θ)a
i
z
z
(2πi)|a| θ|a|
(2πi)|a| θ|a| i
{i|ai >0}

(5.3)

=

(−1)|a|−1 (Dz θ)a X
(−1)|a| (Dz θ)a
a
=
(−|a|).
i
|a| θ |a|
(2πi)|a| θ|a|
(2πi)
i

So we can write
(5.4)

µa =

(−1)|a| (Dz θ)a
(1 − |a|) + fa
(2πi)|a| θ|a|

where fa is a meromorphic function with poles only along Θτ of order strictly less than |a|. Since
we are assuming |a| > 1, we see that 1 − |a| 6= 0, hence ma has a pole of order |a| along Θτ if and
only if (Dz θ)a is not identically zero along Θτ . For the cumulants, the term with the highest pole
along Θ is exactly the same as (5.4): one can prove this by using an explicit expression for Dza log θ
or by the formula relating cumulants with central moments [20].
To conclude we need to show that if the theta divisor Θτ is irreducible, then Dza θ does not vanish
∂θ
identically along it. Since Θ is irreducible and Dza θ is a product of terms of the form ∂z
it is enough
i
to check that each one of these terms does not vanish identically along Θτ . But this follows from [4,
Proposition 4.4.1].

Remark 5.2. We see that the cumulants correspond to the logarithmic derivatives of the theta
function, which are classical meromorphic functions on abelian varieties [15, Section I.6, Method
III]. However, the central moments give, to the best of our knowledge, new meromorphic functions.
5.1. Statistical maps of abelian varieties. By definition, abelian varieties can be embedded
into projective space. In algebraic geometry it is of great interest to produce rational maps from
an abelian variety Aτ to different PN . We are showing that this is possible using statistical data,
namely central moments and cumulants.
Let us fix τ ∈ Hg and consider all the higher central moments [ma ]|a|≤d,|a|6=1 of order d ≥ 2.

Then, according to Proposition 5.1, we get rational maps Aτ 99K PNg,d where Ng,d = g+d
− g − 1.
d
Observe that we do not consider the linear central moments, because they are always zero by
definition. These maps can be seen as similar to the moment varieties for Gaussians and their
mixtures studied in [2],[3]. We proceed to study here the geometry of these maps and derive some
statistical consequences.

DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS

15

First, as in [2], we can replace the central moments ma [X(z,τ ) ] with the cumulants κa [X(z,τ ) ],
which they are often easier to work with. Moreover, we want to work with holomorphic functions
instead than with meromorphic functions, so that we multiply the coordinates everywhere by θd .
So, the maps that we want to consider are the following:
Definition 5.3 (Statistical maps of abelian varieties). With the above notation, we define the
statistical maps


g+d
Ng,d
d
(5.5)
φd,τ : Aτ 99K P
,
φd,τ = [θ κa ]|a|≤d,|a|6=1 ,
Ng,d =
− g − 1.
d
To study the geometry of these maps it will be useful to consider their restriction to the theta
divisor. It turns out that this can be described in terms of the classical Gauss map of the theta
divisor. We recall briefly the definition of the Gauss map here, and for more details one can look at
[4, Section 4.4].
First, if we differentiate the quasiperiodicity relation (3.1) we get:




∂θ
1 t
1 t
∂θ
t
t
(5.6)
(z + m + τ n, τ ) = −2πini e − n τ n − n z θ(z, τ ) + e
(z, τ ).
n τn + n z
∂zi
2
2
∂zi
In particular, when θ(z, τ ) = 0 this gives


1 t
∂θ
∂θ
t
(z + m + τ n, τ ) = e
(z, τ ).
n τn + n z
(5.7)
∂zi
2
∂zi
Hence, on the theta divisor Θτ ⊆ Aτ we have a well defined rational map


∂θ
∂θ
g−1
(5.8)
γ : Θτ 99K P
γ=
,...,
∂z1
∂zg
called the Gauss map. By definition, the Gauss map is not defined at the points of Θτ where all the
partial derivatives vanish. This means that the Gauss map is not defined precisely at the singular
points of Θτ . One can show that this is the same map induced by the complete linear system
H 0 (Θτ , OΘτ (Θτ )).
In the next lemma we show how to describe the restriction of φd,τ to Θτ in terms of the Gauss
map.
Lemma 5.4. The restriction φd,τ |Θτ : Θτ 99K PNg,d corresponds to the composition
(5.9)

γ

vd

Θτ 99K Pg−1 ,→ P(

g−1+d
d

)−1 ,→ PNg,d

where vd isthe d-th Veronese embeddding of Pg−1 , and the last map is the linear embedding in the
last g−1+d
coordinates that sets all the others to zero.
d
Proof. By definition, the coordinates of φd,τ are given by θd κa , for |a| ≤ d, |a| 6= 1. The proof of
Proposition 5.1 shows that we can write them in the form
(−1)|a|
θd
(1
−
|a|)
·
(Dz θ)a + θd ga
(2πi)|a|
θ|a|
where the ga are meromorphic functions with poles only along Θτ and of order strictly smaller than
|a|. In particular, since |a| ≤ d, it follows that the functions θd ga vanish identically along the theta
d
divisor. The same is true for the functions θθ|a| (Dz θ)a , when |a| < d. Hence, when we restrict the
θd κa to the theta divisor Θτ , the only nonzero terms are those corresponding to |a| = d. In this
case, we get
(5.10)

(5.11)

θ d κa =

(θd κa )|Θτ =

(−1)d
(1 − d) · (Dz θ)a|Θτ ,
(2πi)d

d

(−1)
The common nonzero constant (2πi)
d (1−d) does not matter since we are taking a map into projective
space. Hence, we just need to consider the functions {(Dz θ)a|Θτ }|a|=d . However, since the (Dz θ)a
correspond to taking monomials of degree d in the partial derivatives Dzi θ, this is the same as the
d-th Veronese embedding composed with the Gauss map, and we are done.


Using this lemma we can draw consequences on the maps φd,τ . For example, it is immediate to
see where the φd,τ are defined.

DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS

16

Corollary 5.5. The rational map φd,τ : Aτ 99K PNg,d is not defined precisely at the singular points
of Θτ . In particular, it is defined everywhere if and only if Θτ is smooth.
Proof. Since the first coordinate of the map is θd , the map is defined everywhere outside of Θτ .
Moreover, Lemma 5.4 shows that the restriction of the map to Θτ is essentially given by the Gauss
map composed with an embedding. Hence, it is not defined if and only if the Gauss map is not
defined, which happens precisely at the singular points of Θτ .

We usually want maps that are nondegenerate, in the sense that its image cannot be contained
in a hyperplane. We now give a sufficient condition for this not to happen to the maps φd,τ .
Lemma 5.6. Suppose that the theta divisor Θτ is irreducible. Then the image of the rational map
φd,τ : Aτ 99K PNg,d is not contained in a hyperplane.
Proof. To show that the image is not contained in any hyperplane is the same as proving that
the functions θd κa for |a| ≤ d, |a| =
6 1 are linearly independent. we are going to prove linear
independence by induction on d. The base cases d = 0 and d = 1 are immediate. So, assume we
have a linear relation
X
X
(5.12)
λ0 θ d +
λa θd κa + · · · +
λa θd κa = 0,
for certain λa ∈ C.
|a|=2

|a|=d

By induction hypothesis, it is enough to prove that λa = 0 for all |a| = d. Since Θτ is irreducible,
Lemma 5.1 shows that each κa has poles of order exactly |a| along Θτ . Hence, if we restrict the
above relation to Θτ , we get that (θd κa )|Θτ = 0 for all |a| < d, and we are left with
X
(5.13)
λa (θd κa )|Θτ = 0.
|a|=d

It is then enough to show that the functions (θd κa )|Θτ , for |a| = d, are linearly independent.
g−1+d
Equivalently, we need to show that the image of the map Θ 99K P( d )−1 that these functions
τ

induce is not contained in any hyperplane. We know from Lemma 5.4 that this map is a composition
(5.14)

γ

vd

g−1+d
d

Θτ 99K Pg−1 ,→ P(

)−1

where γ is the Gauss map and vd is a Veronese embedding. Since Θτ is irreducible, the Gauss map
is dominant [4, Proposition 4.4.2] and the image of the Veronese map is also not contained in any
hyperplane. Hence, the image of the composite map is not contained in any hyperplane, and we are
done.

Remark 5.7. We can give a geometric interpretation of this result, which will help us study the
maps φd,τ . By construction, φd,τ is defined by a subset of functions in H 0 (Aτ , dΘτ ). It is known
[4, Proposition 4.1.5] that if d ≥ 2, the complete linear system H 0 (Aτ , dΘτ ) induces everywhere
defined maps
(5.15)

g −1

ψd,τ : Aτ −→ Pd

.

If Θτ is irreducible, Lemma 5.6 says that our maps φd,τ can be realized as a composition
(5.16)

ψd,τ

Aτ −→ Pd

g −1

π

99K PNg,d

where the second map is a projection from a linear space. Hence, we can study φd,τ through the
maps ψd,τ , which are well-known for algebraic geometers. For example, the fibers of the map ψ2,τ
correspond exactly to the points z, −z exchanged by the involution. In contrast, the maps ψd,τ are
all closed embeddings for d ≥ 3.
We close this subsection with the following proposition.
Proposition 5.8. Let F : Aτ → Aτ 0 be an isomorphism of polarized abelian varieties of dimension
g. Then for every d there exists a linear isomorphism Q : PNg,d → PNg,d such that
(5.17)

φd,τ 0 ◦ F = Q ◦ φd,τ

DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS

17

Proof. If two polarized abelian varieties (Aτ , Θτ ) and (Aτ 0 , Θτ 0 ) are isomorphic, we know from the
discussion in Section 3.2 that it must be that τ 0 = M τ for a certain M ∈ Sp(2g, Z), and then
the isomorphism is given by FM,τ as in (3.6). Then, we can compute the relation between the
logarithmic derivatives of θ(z, τ ) and θ(FM,τ (z), M τ ), using the Theta Transformation Formula
(3.9), and the conclusion follows.

What this result is saying is that when we take two isomorphic abelian varieties (Aτ , Θτ ) ∼
=
(Aτ 0 , Θτ 0 ), the two maps φd,τ , φd,τ 0 differ only by a linear change of coordinates. Hence, we can
say that φd,τ depends only on the isomorphism class of the polarized abelian variety (Aτ , Θτ ). In
particular, the statistical map φd can be defined for every abstract principally polarized abelian
variety (A, Θ).
In the next two subsections we study in more detail the univariate case g = 1 and the bivariate
case g = 2.
5.2. The case g = 1. We first consider univariate discrete Gaussians, which correspond to elliptic
curves on the abelian varieties side. In this case, we are going to see that actually, the maps φd,τ
coincide with the well-studied maps ψd,τ . This means that elliptic normal curves are parametrized
by cumulants of univariate discrete Gaussians.
More precisely, fix a τ ∈ H1 , and for every d ≥ 2, consider the map
(5.18)

φd,τ : Aτ −→ Pd−1 ,

φd,τ = [θd , θd κ2 , . . . , θd κd ].

For elliptic curves, the theta divisor consists of a single point, so that it is immediately smooth
and irreducible. By Corollary 5.5, the statistical map φd,τ is defined everywhere, and according to
Remark 5.7, the map is a composition
(5.19)

ψd,τ

π

Aτ −→ Pd−1 99K PN1,d

where the second map is a linear projection. But N1,d = d − 1, so that the projection π must be the
identity. Hence, we can identify φd,τ with ψd,τ . This means the following:
• d = 2: the map φ2,τ : Aτ → P1 is surjective of degree 2. The fibers correspond to opposite
points z, −z of Aτ . In particular, it is ramified exactly at the four points of order two in Aτ .
• d ≥ 3: the map φd,τ : Aτ ,→ Pd−1 is a closed embedding. The image is usually called an
elliptic normal curve of degree d. In particular, the map φ3,τ embeds Aτ as a smooth cubic
curve in P2 .
This sheds light on the relation between univariate discrete Gaussians and their moments.
Theorem 2.4 tells us that when the parameters z, τ are purely imaginary, the discrete Gaussians are
uniquely determined by the moments up to order two. In contrast, we are going to show now that
arbitrary (that is, complex) univariate discrete Gaussians of a fixed parameter τ are determined by
the moments up to order three.
Theorem 5.9. Fix τ ∈ H1 and let X(z,τ ) , X(z 0 ,τ ) be two discrete Gaussians with the same parameter
τ . If the first three moments coincide, then they have the same distribution. Furthermore, if we
consider discrete Gaussians with different parameters τ , the first two moments are not enough in
general.
Proof. We prove now the first part. If the first three central moments coincide, then the cumulants
κ2 , κ3 coincide as well. This means that the two points z, z 0 ∈ Aτ have the same image under the
map φ3,τ : Aτ → P2 . But we know from the previous discussion that this map is injective, hence,
it must be that z, z 0 are the same point in Aτ , which is to say that z 0 = z + m + τ n for certain
n, m ∈ Z. Then, Proposition 3.1 shows that X(z 0 ,τ ) ∼ X(z,τ ) − n, but since X(z,τ ) , X(z 0 ,τ ) have the
same mean, it must be that n = 0, and we conclude.

A statistical consequence is that all higher moments of a Gaussian distribution are uniquely
determined by the first three. We can make this explicit using the equation of the cubic curve
φ3,τ (Aτ ) ⊆ P2 .
To compute this equation we set some notations. For every fixed τ ∈ H1 , we consider the three
points 0, 12 , τ2 ∈ C. If we take their images in Aτ , these give three of the 2-torsion points of Aτ . The

DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS

other torsion point is given by the image of
[15, Lemma 4.1]. Now, we define
(5.20)
e1 (τ ) :=

Dz2 (log θ(z, τ )) (0, τ ) ,

e2 (τ ) :=

and the quantities
1
(5.21) A(τ ) = 2 (e1 + e2 + e3 ),
π
Then the following holds:

1
2

+ τ2 , which coincides with the theta divisor Θτ ⊆ Aτ

Dz2 (log θ(z, τ ))

B(τ ) =

18



1
,τ
2


,

e3 (τ ) :=

1
(e1 e2 + e1 e3 + e2 e3 ),
4π 4

Dz2 (log θ(z, τ ))

C(τ ) =

τ
2

,τ ,

1
e1 e2 e3 .
16π 6

Proposition 5.10. Fix τ ∈ H1 . Then the image of the map φ3,τ : Aτ → P2[X0 ,X1 ,X2 ] is the cubic
(5.22)

X0 X22 = −4X13 − A(τ )X0 X12 − B(τ )X02 X1 − C(τ )X03 .

Equivalently, for each (z, τ ) ∈ C × H1 \ Θ1 we have the relation
(5.23)

κ23 = −4κ32 − A(τ )κ22 − B(τ )κ2 − C(τ ).

Proof. For notational clarity we set νa = (2πi)a κa = Dza log θ. Now, we fix τ and we consider the
meromorphic function on Aτ given by ν32 + 4ν23 : we can check, for example by a computer, that this
has a pole of order at most four along Θτ . Let H 0 (Aτ , 4Θτ ) be the space of meromorphic functions
on Aτ with poles at most of order 4 along Θτ . This space has dimension 4 by the classical theorem
of Riemann-Roch, and the functions 1, ν2 , ν3 , ν22 are linearly independent elements in H 0 (Aτ , 4Θτ ),
since they have poles of different order along Θτ . Hence, they form a basis of the space, so that
(5.24)

ν32 = −4ν23 + a(τ ) · ν22 + b(τ ) · ν2 + c(τ ) + d(τ ) · ν3

for certain a(τ ), b(τ ), c(τ ), d(τ ) independent of z. Now we observe that the functions ν32 and
−4ν23 + a(τ )ν22 + b(τ )ν2 + c(τ ) are even, whereas the function d(τ ) · ν3 is odd. Hence, it must be
that d(τ ) = 0. The equation reduces to:
(5.25)

ν32 = −4ν23 + a(τ )ν22 + b(τ )ν2 + c(τ ).

Now, consider the map φ2,τ = [1, ν2 ] : Aτ → P1 : we know from a previous discussion that this is
ramified precisely at the two torsion points of Aτ . Since Dz (ν2 ) = ν3 , this means precisely that
ν3 (0) = ν3 ( 12 ) = ν3 ( τ2 ) = 0. Plugging this into (5.25), we get that e1 = ν2 (0), e2 = ν2 ( 12 ), e3 = ν2 ( τ2 )
are roots of the polynomial −4x3 + a(τ )x2 + b(τ )x + c(τ ). Then, we can write
(5.26)

a(τ ) = −4(e1 + e2 + e3 ),

b(τ ) = 4(e1 e2 + e1 e3 + e2 e3 ),

To conclude, it is enough to replace νa = (2πi)a κa .

c(τ ) = −4e1 e2 e3 .


Remark 5.11. This relation tells us that the second cumulant determines the third cumulant up to
a sign. We can also get the higher cumulants as follows: using the same notation as in the previous
proof, consider the equation
(5.27)

ν32 = −4ν23 + a(τ )ν22 + b(τ )ν2 + c(τ ).

Then, taking derivatives with respect to z we get
b(τ )
.
2
Passing to the cumulants, this gives a formula to compute the fourth cumulant from the second.
Differentiating again, we get a formula for the fifth cumulant and so on. All together, this gives
explicit formulas to compute the higher moments, starting from the first three.
It is worth mentioning that the equation we obtain in Proposition 5.10 is the same as the one
obtained through the Weierstrass ℘-function [15, Section 1.6]. In addition, the relations we obtain
by differentiating can be interpreted through the Korteweg-de Vries equation, as in [15, Section 1.6].
(5.28)

ν4 = −6ν22 + a(τ )ν2 +

As an application of the cubic equation of Proposition 5.10, we are going to show that univariate
discrete Gaussians of arbitrary parameters z, τ are not uniquely determined by their moments up to
order two.
To do this, we consider the map

(5.29)
Φ : (C × H1 ) \ Θ1 → C2
(z, τ ) 7→ (µ, Σ) = µ[X(z,τ ) ], Var[X(z,τ ) ] .



DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS

19

This assigns to each couple of parameters (z, τ ) the mean and the covariance of the corresponding
discrete Gaussian distribution. By construction, this map factors through the space G1 of all
univariate discrete Gaussian distributions. Hence, we get a map
e : G1 → C2
Φ

(5.30)

and we are asking whether this map is injective. We show here that this is not the case.
e is not injective. Hence, mean and covariance do not determine an
Theorem 5.9 says the map Φ
univariate discrete Gaussian in general. We conclude the proof now.
Proof. (Theorem 5.9) We know from Remark 3.4 that G1 is a complex manifold of dimension two.
e is injective, its differential should be everywhere an isomorphism between tangent spaces.
Hence, if Φ
To check this, we may as well work with Φ, for which we can compute the differential explicitly. We
have



1
1
(5.31)
Φ(z, τ ) = µ[X(z,τ ) ], Var[X(z,τ ) ] =
ν1 ,
ν2
(2πi)
(2πi)2
where ν1 , ν2 are the same as in the proof of Proposition 5.10. Up to multiplication by a nonzero
constant, the differential of Φ is given by
 1
 

Dz ν1 Dτ1 ν1
ν2 Dτ1 ν1
(5.32)
DΦ =
=
.
Dz1 ν2 Dτ1 ν2
ν3 Dτ1 ν2
However, we see that Dτ1 ν1 = Dτ1 Dz1 (log θ) = Dz1 Dτ1 (log θ) and the same holds for ν2 . Then, we can
use the heat equation (4.23), and we get that

1 Dz2 θ
1
1
Dτ1 θ
=
=
Dz2 (log θ) + Dz1 (log θ)2 =
(ν2 + ν12 )
θ
4πi θ
4πi
4πi
Thus, up to multiplying by a nonzero constant, we can write


ν2
ν3 + 2ν1 ν2
(5.34)
det DΦ = det
ν3 ν4 + 2ν22 + 2ν1 ν3


ν2
ν3
= det
(5.35)
ν3 ν4 + 2ν22


ν
ν3
= det 2
(5.36)
= ν2 ν4 + 2ν23 − ν32 .
ν3 ν4 + 2ν22
(5.33)

Dτ1 (log θ) =

Using the equations (5.25) and (5.28), we get that:
b(τ )
ν2 + c(τ )
2
where b(τ ), c(τ ) are as in the proof of Proposition 5.10. We want to show that this vanishes on some
points of C × H1 \ Θ1 . It is enough to show that b(τ ) is not identically zero. Indeed, if we fix a
τ such that b(τ ) 6= 0, the function − b(τ2 ) ν2 + c(τ ) is a meromorphic function on Aτ with poles of
order two along the theta divisor Θτ , so that it must have zeroes on Aτ \ Θτ . A SAGE computation
reveals that b(τ ) is indeed not identically zero (by evaluating its explicit expression (5.26)), and we
are done.

(5.37)

ν2 ν4 + 2ν23 − ν32 = −

5.3. The case g = 2. In the case of dimension two the situation becomes more complicated.
First, let’s fix a parameter τ ∈ H2 such that the theta divisor Θτ is irreducible. In this case [4,
Corollary 11.8.2] Θτ is also smooth: more precisely, Θτ is a smooth curve of genus 2 and Aτ is its
Jacobian variety.
Remark 5.12. It will be useful to single out some results about the Gauss map
(5.38)

γ : Θτ −→ P1

in this situation. We have seen before that this map corresponds to the one induced by the line
bundle OΘτ (Θτ ). Since Θτ is a smooth curve, the adjunction formula [10, Formula II.8.20] shows
that OΘτ (Θτ ) is the canonical bundle of Θτ . Since Θτ has genus two, it follows that γ is a double
cover of P1 . We can identify the fibers explicitly: indeed we see from the explicit form (5.8) that
the Gauss map is invariant under the involution z 7→ −z of Θτ . Since this map has degree precisely

DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS

20

two by the previous discussion, it follows that the fibers of γ consist precisely of opposite points
z, −z, with z ∈ Θτ .
Now set d ≥ 2. Since Θτ is smooth, Corollary 5.5 shows that the map
(d + 2)(d + 1)
−3
2
is everywhere defined. Moreover, since Θτ is irreducible, Remark 5.7 shows that this map can be
realized as a composition
(5.39)

φτ,d : Aτ → PN2,d ,

ψd,τ

N2,d =

2 −1

Aτ −→ Pd

(5.40)

π

99K PN2,d

where π is a linear projection. We study this situation when d = 2 and d = 3.
• d = 2. In this case, we see that d2 − 1 = N2,d = 3, so that the projection π must be the
identity. Hence we can identify the map φ2,τ with ψ2,τ . This map is well studied: it induces
a degree 2 cover of Aτ onto a quartic surface Xτ ⊆ P3 with 16 nodes. This surface is called
the Kummer surface of Aτ , and the map φ2,τ : Aτ → Xτ realizes it as the quotient of Aτ
under the involution z 7→ −z. In particular, the 16 nodes correspond to the points of order
two in Aτ .
• d = 3. In this case, we see that N2,d = 7, whereas d2 − 1 = 8. Hence, the map φ3,τ is the
composition of the closed embedding ψ3,τ : Aτ ,→ P8 , with a projection π : P8 99K P7 from a
point P . We want to identify geometrically this point. To do this, take any z ∈ Θτ which
is not a point of order two: then we claim that φ3,τ (z) = φ3,τ (−z). Taking this claim for
granted, we can rephrase it by saying that the two points ψ3,τ (z), ψ3,τ (−z) must lie on a
common line with the point P . Hence, the point P can be characterized as the common
intersection point of all the lines in P8 spanned by the couples ψ3,τ (z), ψ3,τ (−z), for z ∈ Θτ .
To conclude we need to prove the claim. We know that the restriction of φ3,τ to the
theta divisor Θτ is the composition of the Gauss map γ : Θτ → P1 with a closed embedding.
Hence, it is enough to show that γ(z) = γ(−z) for all z ∈ Θτ . However, this follows from
the discussion in Remark 5.12.
What if Θτ is reducible? Then [4, Corollary 11.8.2] shows that the polarized abelian surface
(Aτ , Θτ ) must be a product of elliptic curves. If we want
 to study the geometry of the maps φd,τ ,
τ1 0
we can use Proposition 5.8, and assume that τ = 0 τ2 , so that (Aτ , Θτ ) = (Aτ1 , Θτ1 ) × (Aτ2 , Θτ2 ).
By Proposition 3.9, this corresponds to the case of two independent discrete Gaussians, and in
particular it is straightforward to compute the cumulants κ(a1 ,a2 ) . We have that κ(a1 ,0) ((z1 , z2 ), τ ) =
κa1 (z1 , τ1 ), κ(0,a2 ) ((z1 , z2 ), τ ) = κa2 (z1 , τ1 ), and the mixed cumulants κ(a1 ,a2 ) ((z1 , z2 ), τ ) vanish
whenever one of a1 , a2 is nonzero. Hence, we can write the rational maps φd,τ as the composition of
the two maps
(5.41)

A τ1 × A τ2

φd,τ1 ×φd,τ2

−→

fd

Pd−1 × Pd−1 99K P2d−2 ,

where the second map fd is given by
(5.42)

([X0 , . . . , Xd−1 ], [Y0 , . . . , Yd−1 ]) 7→ [X0 Y0 , Y0 X1 , . . . , Y0 Xd−1 , X0 Y1 , . . . , X0 Yd−1 ]

Remark 5.13. In particular, we observe that the map φd,τ1 × φd,τ2 sends the open subset Aτ \ Θτ
into the open subset {X0 6= 0}×{Y0 6= 0}, and the map fd restricted to this subset is an isomorphism
onto the image.
As in Theorem 5.9, we can use the statistical map φ3,τ to see that discrete Gaussians of dimension
2 with the same parameter τ are determined by the moments up to order three.
Theorem 5.14. Fix τ ∈ H2 and let X(z,τ ) , X(z 0 ,τ ) be two discrete Gaussians with the same parameter
τ . If all the moments up to order three coincide, then they have the same distribution.
Proof. Proceeding in the same way as in dimension 1, it is enough to show that the map φ3,τ : Aτ →
P3 is injective on the open set Aτ \ Θτ .
Suppose first that Θτ is irreducible. If we look at Aτ as embedded in P8 by ψ3,τ , the injectivity
of φ3,τ is equivalent to saying that no two points x, y in Aτ \ Θτ lie on a line passing through the
point P . Suppose that this happens and take two distinct points z, −z ∈ Θτ exchanged by the
involution. Then the line that they span is also passing through P , so that the four distinct points

DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS

21

x, y, z, −z span a two-dimensional plane in P8 . However, Lemma 5.15 below shows that in this case
x, y belong to Θτ as well, which is absurd.

Suppose instead that Θτ is reducible: then as before we can assume τ = τ01 τ02 , so that the map
φ3,τ is a composition of the map φ3,τ1 × φ3,τ2 and of the map f3 , as in (5.41). Since φ3,τ1 × φ3,τ2 is
a product of closed embeddings, it is injective. Moreover, Remark 5.13 shows that the map f3 is
injective when restricted to the image of Aτ \ Θτ , and we are done.

To conclude the previous proof, we need a more technical lemma from Algebraic Geometry. In
particular, we will make use some intersection theory on surfaces, for which we refer to [10, Section
V.1].
Lemma 5.15. Fix τ ∈ H2 such that Θτ is irreducible and consider the embedding ψ3,τ : Aτ ,→ P8 .
Take two opposite points z, −z ∈ Θτ and other two points x, y ∈ Aτ . Suppose that the points
ψ3,τ (z), ψ3,τ (−z), ψ3,τ (x), ψ3,τ (y) span a plane of dimension two in P8 . Then x, y belong to Θτ as
well.
Proof. With these assumptions, the proof of [17, Theorem 5.7] shows that the points z, −z, x, y lie
on a smooth curve C ⊆ Aτ of genus two, such that (C · Θτ ) = 2. Since (C 2 ) = 2, an application
of the Hodge Index Theorem [10, Theorem V.1.9] shows that C and Θτ are actually numerically
equivalent curves on Aτ . Since Aτ is an abelian variety, it follows [4, Theorem 4.11.1] that C is
actually the translate of Θτ by a point c ∈ Aτ , so that C = Θτ + c. If c = 0, we are done. Suppose
that c 6= 0: this means [4, Lemma 11.3.1] that the restriction OΘτ (C) is a line bundle of degree
two which is not the canonical line bundle on the curve Θτ . However, the line bundle OΘτ (C)
corresponds by definition to the divisor Θτ ∩ C = {z, −z}, and we have seen in Remark 5.12 that
this is a canonical divisor on Θτ . This gives a contradiction and we are done.

5.4. Open questions. From these results in dimension one and two, there are natural questions
that arise:
Question 5.16. We have seen in Theorem 5.9 and Theorem 5.14 that two Gaussian distributions
in dimension 1 or 2 with the same parameter τ can be distinguished by the moments up to order
three. Is this true for every dimension g?
Following the proofs of Theorems 5.9 and 5.14, a geometric way to prove this would be to show
that the maps
(5.43)

φ3,τ : Aτ \ Θτ −→ PN3,g

are injective for every τ ∈ Hg . Moreover, it could be that for a fixed parameter τ , two Gaussian
distributions can be distinguished just by the moments up to order two. We would like to identify
when this happens too.
We can generalize this question to all Gaussian distributions at the same time:
Question 5.17. We have seen in Theorem 5.9 that in general mean and covariance are not enough
to determine a discrete Gaussian. Could this be true if we take all moments up to order three? Or
up to an higher order d? Geometrically, this would mean to study the fibers of the map
(5.44)

Φd,g : Cg × Hg \ Θg −→ PNd,g ,

(z, τ ) 7→ φd,τ (z).

Remark 5.18. Taking our correspondence further, mixtures of discrete Gaussians correspond to
secants of abelian varieties. More precisely, if τ is fixed then mixtures of two discrete Gaussians with
parameters (z1 , τ ), (z2 , τ ) form a secant line in the corresponding abelian variety, while mixtures
with (z1 , τ1 ), (z2 , τ2 ) lie in a secant line to the universal family Ug . This should be a very interesting
connection to explore, with natural recurring questions such as identifiability from moments (the
case of continuous Gaussian mixtures is treated in [3]).
Acknowledgements
We would like to give special thanks to Bernd Sturmfels for dreaming and suggesting the
connection between multivariate discrete Gaussians and theta functions. We thank Lynn Chua for
her help with SAGE computations. Thanks to Samuel Grushevsky for his interest and his valuable
comments. Thanks also to Daniel Dadush, Oded Regev and Noah Stephens for offering the CS
perspective of discrete Gaussians, and to Robert Gray for his interest in this work. Daniele Agostini

DISCRETE GAUSSIAN DISTRIBUTIONS VIA THETA FUNCTIONS

22

was supported by the DAAD, the DFG Graduiertenkolleg 1800, the DFG Schwerpunkt 1489 and
the Berlin Mathematical School. He also thanks the Stony Brook Mathematics Department for its
hospitality while writing this work. Carlos Améndola was supported by the Einstein Foundation
Berlin.
References
[1] Aggarwal, D., Dadush, D. and Stephens-Davidowitz N. (2015) Solving the closest vector problem in 2n time - the
discrete Gaussian strikes again!. FOCS, IEEE 56th Annual Symposium, 563-582.
[2] Améndola, C., Faugère, J.C. and Sturmfels, B. (2016) Moment varieties of Gaussian mixtures. Journal of Algebraic
Statistics 7, 14-28.
[3] Améndola, C., Ranestad, K. and Sturmfels, B. (2018) Algebraic identifiability of Gaussian mixtures. International
Mathematics Research Notices.
[4] Birkenhake, C. and Lange, H. (2013). Complex Abelian Varieties, Vol 302. Springer Science & Business Media.
[5] Brown, L.D. (1986) Fundamentals of Statistical Exponential Families with Applications in Statistical Decision
Theory. Lecture Notes-monograph series 9, i-279.
[6] Cover, T.M. and Thomas, J.A. (2006). Elements of Information Theory. Wiley, 409-413.
[7] Deconinck, B., Heil, M., Bobenko, A., Van Hoeij, M. and Schmies M. (2003). Computing Riemann Theta functions.
Mathematics of Computation, 73, 247, 1417-1442.
[8] Drton, M., Sturmfels B. and Sullivant, S. (2009) Lectures on Algebraic Statistics, Oberwolfach Seminars, Vol 39,
Birkhäuser, Basel.
[9] Dwarakanath, N.C. and Galbraith, S.D. (2014). Sampling from discrete Gaussians for lattice-based cryptography
on a constrained device. Applicable Algebra in Engineering, Communication and Computing, 25(3), 159-180.
[10] Hartshorne, R. (1977). Algebraic Geometry. Springer.
[11] Kemp, A. W. (1997). Characterizations of a discrete normal distribution. Journal of Statistical Planning and
Inference, 63(2), 223-229.
[12] Lisman, J.H. and Zuylen, M.C. (1972). Note on the generation of most probable frequency distributions. Statistica
Neerlandica, 26, 19-23.
[13] Micciacio, D. and Walter, M. (2017). Gaussian sampling over the integers. IACR Cryptology ePrint Archive, p.
259.
[14] Michalek, M., Sturmfels, B., Uhler C. and Zwiernik, P. (2016). Exponential varieties. Proceedings of the London
Mathematical Society 112, 27-56.
[15] Mumford, D. (1983). Tata Lectures on Theta I. Birkhäuser.
[16] Regev, O. and Stephens-Davidowitz, N. (2017). An inequality for Gaussians on lattices. SIAM Journal on
Discrete Mathematics, 31(2), 749-757.
[17] Terakawa, H. (1998). The k-very ampleness and k-spannedness on polarized abelian surfaces. Mathematische
Nachrichten, 195, 237-250.
[18] Swierczewski, C. and Deconinck, B. (2013). Computing Riemann Theta Functions in Sage with applications.
Mathematics and Computers in Simulation, 127, 263-272.
[19] Szablowski, P. J. (2001). Discrete normal distribution and its relationship with Jacobi Theta functions. Statistics
& probability letters, 52(3), 289-299.
[20] Willink, R. (2003). Relationships between central moments and cumulants, with formulae for the central moments
of Gamma distributions. Communications in Statistics-Theory and Methods, Vol. 32(4), 701-704.
Humboldt-Universität zu Berlin, Institut für Mathematik, Unter den Linden 6, 10099, Berlin
Germany
E-mail address: daniele.agostini@math.hu-berlin.de
Technische Universität Berlin, Institut für Mathematik, Straße des 17. Juni 136, 10623, Berlin
Germany
E-mail address: amendola@math.tu-berlin.de

