Discrete Gaussian Samplers over Lattices#

This file implements oracles which return samples from a lattice following a discrete Gaussian distribution. That is, if \(\sigma\) is big enough relative to the provided basis, then vectors are returned with a probability proportional to \(\exp(-|x - c|_2^2 / (2\sigma^2))\). More precisely lattice vectors in \(x \in \Lambda\) are returned with probability:

\[\frac{\exp(-|x - c|_2^2 / (2\sigma^2))}{\sum_{v \in \Lambda} \exp(-|v|_2^2 / (2\sigma^2))}.\]

AUTHORS:

Martin Albrecht (2014-06-28): initial version
Gareth Ma (2023-09-22): implement non-spherical sampling

EXAMPLES:

Sage

sage: D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ^10, 3.0)
sage: D(), D(), D()  # random
((3, 0, -5, 0, -1, -3, 3, 3, -7, 2), (4, 0, 1, -2, -4, -4, 4, 0, 1, -4), (-3, 0, 4, 5, 0, 1, 3, 2, 0, -1))
sage: a = D()
sage: a.parent()
Ambient free module of rank 10 over the principal ideal domain Integer Ring

Python

>>> from sage.all import *
>>> D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ**Integer(10), RealNumber('3.0'))
>>> D(), D(), D()  # random
((3, 0, -5, 0, -1, -3, 3, 3, -7, 2), (4, 0, 1, -2, -4, -4, 4, 0, 1, -4), (-3, 0, 4, 5, 0, 1, 3, 2, 0, -1))
>>> a = D()
>>> a.parent()
Ambient free module of rank 10 over the principal ideal domain Integer Ring

class sage.stats.distributions.discrete_gaussian_lattice.DiscreteGaussianDistributionLatticeSampler(B, sigma=1, c=0, r=None, precision=None, sigma_basis=False)[source]#

Bases: SageObject

GPV sampler for Discrete Gaussians over Lattices.

EXAMPLES:

Sage

sage: D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ^10, 3.0); D
Discrete Gaussian sampler with Gaussian parameter σ = 3.00000000000000, c=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0) over lattice with basis

[1 0 0 0 0 0 0 0 0 0]
[0 1 0 0 0 0 0 0 0 0]
[0 0 1 0 0 0 0 0 0 0]
[0 0 0 1 0 0 0 0 0 0]
[0 0 0 0 1 0 0 0 0 0]
[0 0 0 0 0 1 0 0 0 0]
[0 0 0 0 0 0 1 0 0 0]
[0 0 0 0 0 0 0 1 0 0]
[0 0 0 0 0 0 0 0 1 0]
[0 0 0 0 0 0 0 0 0 1]

Python

>>> from sage.all import *
>>> D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ**Integer(10), RealNumber('3.0')); D
Discrete Gaussian sampler with Gaussian parameter σ = 3.00000000000000, c=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0) over lattice with basis
<BLANKLINE>
[1 0 0 0 0 0 0 0 0 0]
[0 1 0 0 0 0 0 0 0 0]
[0 0 1 0 0 0 0 0 0 0]
[0 0 0 1 0 0 0 0 0 0]
[0 0 0 0 1 0 0 0 0 0]
[0 0 0 0 0 1 0 0 0 0]
[0 0 0 0 0 0 1 0 0 0]
[0 0 0 0 0 0 0 1 0 0]
[0 0 0 0 0 0 0 0 1 0]
[0 0 0 0 0 0 0 0 0 1]

We plot a histogram:

Sage

sage: import warnings
sage: warnings.simplefilter('ignore', UserWarning)
sage: D = distributions.DiscreteGaussianDistributionLatticeSampler(identity_matrix(2), 3.0)
sage: S = [D() for _ in range(2^12)]
sage: l = [vector(v.list() + [S.count(v)]) for v in set(S)]
sage: list_plot3d(l, point_list=True, interpolation='nn')                       # needs sage.plot
Graphics3d Object

Python

>>> from sage.all import *
>>> import warnings
>>> warnings.simplefilter('ignore', UserWarning)
>>> D = distributions.DiscreteGaussianDistributionLatticeSampler(identity_matrix(Integer(2)), RealNumber('3.0'))
>>> S = [D() for _ in range(Integer(2)**Integer(12))]
>>> l = [vector(v.list() + [S.count(v)]) for v in set(S)]
>>> list_plot3d(l, point_list=True, interpolation='nn')                       # needs sage.plot
Graphics3d Object

REFERENCES:

[GPV2008]

__init__(B, sigma=1, c=0, r=None, precision=None, sigma_basis=False)[source]#

Construct a discrete Gaussian sampler over the lattice \(\Lambda(B)\) with parameter sigma and center \(c\).

INPUT:

B – a (row) basis for the lattice, one of the following:
- an integer matrix,
- an object with a .matrix() method, e.g. ZZ^n, or
- an object where matrix(B) succeeds, e.g. a list of vectors
sigma – Gaussian parameter, one of the following:
- a real number \(\sigma > 0\) (spherical),
- a positive definite matrix \(\Sigma\) (non-spherical), or
- any matrix-like S, equivalent to \(\Sigma = SS^T\), when sigma_basis is set
c – (default: 0) center \(c\), any vector in \(\ZZ^n\) is supported, but \(c \in \Lambda(B)\) is faster
r – (default: None) rounding parameter \(r\) as defined in [Pei2010]; ignored for spherical Gaussian parameter; if not provided, set to be the maximal possible such that \(\Sigma - rBB^T\) is positive definite
precision – bit precision \(\geq 53\).
sigma_basis – (default: False) When set, sigma is treated as
a (row) basis, i.e. the covariance matrix is computed by \(\Sigma = SS^T\)

Todo

Rename class methods like .f and hide most of them (at least behind something like .data).

EXAMPLES:

Sage

sage: n = 2; sigma = 3.0
sage: D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ^n, sigma)
sage: f = D.f
sage: nf = D._normalisation_factor_zz(); nf                                 # needs sage.symbolic
56.5486677646...

sage: from collections import defaultdict
sage: counter = defaultdict(Integer); m = 0
sage: def add_samples(i):
....:     global counter, m
....:     for _ in range(i):
....:         counter[D()] += 1
....:         m += 1

sage: v = vector(ZZ, n, (-3, -3))
sage: v.set_immutable()
sage: while v not in counter: add_samples(1000)
sage: while abs(m*f(v)*1.0/nf/counter[v] - 1.0) >= 0.1:                     # needs sage.symbolic
....:     add_samples(1000)

sage: counter = defaultdict(Integer); m = 0
sage: v = vector(ZZ, n, (0, 0))
sage: v.set_immutable()
sage: while v not in counter:
....:     add_samples(1000)
sage: while abs(m*f(v)*1.0/nf/counter[v] - 1.0) >= 0.1:                     # needs sage.symbolic
....:     add_samples(1000)

Python

>>> from sage.all import *
>>> n = Integer(2); sigma = RealNumber('3.0')
>>> D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ**n, sigma)
>>> f = D.f
>>> nf = D._normalisation_factor_zz(); nf                                 # needs sage.symbolic
56.5486677646...

>>> from collections import defaultdict
>>> counter = defaultdict(Integer); m = Integer(0)
>>> def add_samples(i):
...     global counter, m
...     for _ in range(i):
...         counter[D()] += Integer(1)
...         m += Integer(1)

>>> v = vector(ZZ, n, (-Integer(3), -Integer(3)))
>>> v.set_immutable()
>>> while v not in counter: add_samples(Integer(1000))
>>> while abs(m*f(v)*RealNumber('1.0')/nf/counter[v] - RealNumber('1.0')) >= RealNumber('0.1'):                     # needs sage.symbolic
...     add_samples(Integer(1000))

>>> counter = defaultdict(Integer); m = Integer(0)
>>> v = vector(ZZ, n, (Integer(0), Integer(0)))
>>> v.set_immutable()
>>> while v not in counter:
...     add_samples(Integer(1000))
>>> while abs(m*f(v)*RealNumber('1.0')/nf/counter[v] - RealNumber('1.0')) >= RealNumber('0.1'):                     # needs sage.symbolic
...     add_samples(Integer(1000))

Spherical covariance are automatically handled:

Sage

sage: distributions.DiscreteGaussianDistributionLatticeSampler(ZZ^3, sigma=Matrix(3, 3, 2))
Discrete Gaussian sampler with Gaussian parameter σ = 2.00000000000000, c=(0, 0, 0) over lattice with basis

[1 0 0]
[0 1 0]
[0 0 1]

Python

>>> from sage.all import *
>>> distributions.DiscreteGaussianDistributionLatticeSampler(ZZ**Integer(3), sigma=Matrix(Integer(3), Integer(3), Integer(2)))
Discrete Gaussian sampler with Gaussian parameter σ = 2.00000000000000, c=(0, 0, 0) over lattice with basis
<BLANKLINE>
[1 0 0]
[0 1 0]
[0 0 1]

The sampler supports non-spherical covariance in the form of a Gram matrix:

Sage

sage: n = 3
sage: Sigma = Matrix(ZZ, [[5, -2, 4], [-2, 10, -5], [4, -5, 5]])
sage: c = vector(ZZ, [7, 2, 5])
sage: D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ^n, Sigma, c)
sage: f = D.f
sage: nf = D._normalisation_factor_zz(); nf # This has not been properly implemented
78.6804...

Python

>>> from sage.all import *
>>> n = Integer(3)
>>> Sigma = Matrix(ZZ, [[Integer(5), -Integer(2), Integer(4)], [-Integer(2), Integer(10), -Integer(5)], [Integer(4), -Integer(5), Integer(5)]])
>>> c = vector(ZZ, [Integer(7), Integer(2), Integer(5)])
>>> D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ**n, Sigma, c)
>>> f = D.f
>>> nf = D._normalisation_factor_zz(); nf # This has not been properly implemented
78.6804...

We can compute the expected number of samples before sampling a vector:

Sage

sage: v = vector(ZZ, n, (11, 4, 8))
sage: v.set_immutable()
sage: 1 / (f(v) / nf)
2553.9461...

sage: counter = defaultdict(Integer); m = 0
sage: while v not in counter:
....:     add_samples(1000)
sage: sum(counter.values())  # random
3000
sage: while abs(m*f(v)*1.0/nf/counter[v] - 1.0) >= 0.1:                     # needs sage.symbolic
....:     add_samples(1000)

Python

>>> from sage.all import *
>>> v = vector(ZZ, n, (Integer(11), Integer(4), Integer(8)))
>>> v.set_immutable()
>>> Integer(1) / (f(v) / nf)
2553.9461...

>>> counter = defaultdict(Integer); m = Integer(0)
>>> while v not in counter:
...     add_samples(Integer(1000))
>>> sum(counter.values())  # random
3000
>>> while abs(m*f(v)*RealNumber('1.0')/nf/counter[v] - RealNumber('1.0')) >= RealNumber('0.1'):                     # needs sage.symbolic
...     add_samples(Integer(1000))

If the covariance provided is not positive definite, an error is thrown:

Sage

sage: Sigma = Matrix(ZZ, [[0, 1], [1, 0]])
sage: distributions.DiscreteGaussianDistributionLatticeSampler(ZZ^2, Sigma)
Traceback (most recent call last):
...
RuntimeError: Sigma(=[0.000000000000000  1.00000000000000]
[ 1.00000000000000 0.000000000000000]) is not positive definite

Python

>>> from sage.all import *
>>> Sigma = Matrix(ZZ, [[Integer(0), Integer(1)], [Integer(1), Integer(0)]])
>>> distributions.DiscreteGaussianDistributionLatticeSampler(ZZ**Integer(2), Sigma)
Traceback (most recent call last):
...
RuntimeError: Sigma(=[0.000000000000000  1.00000000000000]
[ 1.00000000000000 0.000000000000000]) is not positive definite

The sampler supports passing a basis for the covariance:

Sage

sage: n = 3
sage: S = Matrix(ZZ, [[2, 0, 0], [-1, 3, 0], [2, -1, 1]])
sage: D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ^n, S, sigma_basis=True)
sage: D.sigma()
[ 4.00000000000000 -2.00000000000000  4.00000000000000]
[-2.00000000000000  10.0000000000000 -5.00000000000000]
[ 4.00000000000000 -5.00000000000000  6.00000000000000]

Python

>>> from sage.all import *
>>> n = Integer(3)
>>> S = Matrix(ZZ, [[Integer(2), Integer(0), Integer(0)], [-Integer(1), Integer(3), Integer(0)], [Integer(2), -Integer(1), Integer(1)]])
>>> D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ**n, S, sigma_basis=True)
>>> D.sigma()
[ 4.00000000000000 -2.00000000000000  4.00000000000000]
[-2.00000000000000  10.0000000000000 -5.00000000000000]
[ 4.00000000000000 -5.00000000000000  6.00000000000000]

The non-spherical sampler supports offline computation to speed up sampling. This will be useful when changing the center \(c\) is supported. The difference is more significant for larger matrices. For 128x128 we observe a 4x speedup (86s -> 20s):

Sage

sage: D.offline_samples = []
sage: T = 2**12
sage: L = [D() for _ in range(T)] # 560ms
sage: D.add_offline_samples(T)    # 150ms
sage: L = [D() for _ in range(T)] # 370ms

Python

>>> from sage.all import *
>>> D.offline_samples = []
>>> T = Integer(2)**Integer(12)
>>> L = [D() for _ in range(T)] # 560ms
>>> D.add_offline_samples(T)    # 150ms
>>> L = [D() for _ in range(T)] # 370ms

We can also initialise with matrix-like objects:

Sage

sage: qf = matrix(3, [2, 1, 1,  1, 2, 1,  1, 1, 2])
sage: D = distributions.DiscreteGaussianDistributionLatticeSampler(qf, 3.0); D
Discrete Gaussian sampler with Gaussian parameter σ = 3.00000000000000, c=(0, 0, 0) over lattice with basis

[2 1 1]
[1 2 1]
[1 1 2]
sage: D().parent() is D.c().parent()
True

Python

>>> from sage.all import *
>>> qf = matrix(Integer(3), [Integer(2), Integer(1), Integer(1),  Integer(1), Integer(2), Integer(1),  Integer(1), Integer(1), Integer(2)])
>>> D = distributions.DiscreteGaussianDistributionLatticeSampler(qf, RealNumber('3.0')); D
Discrete Gaussian sampler with Gaussian parameter σ = 3.00000000000000, c=(0, 0, 0) over lattice with basis
<BLANKLINE>
[2 1 1]
[1 2 1]
[1 1 2]
>>> D().parent() is D.c().parent()
True

__call__()[source]#

Return a new sample.

EXAMPLES:

Sage

sage: D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ^3, 3.0, c=(1,0,0))
sage: L = [D() for _ in range(2^12)]
sage: mean_L = sum(L) / len(L)
sage: norm(mean_L.n() - D.c()) < 0.25
True

sage: D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ^3, 3.0, c=(1/2,0,0))
sage: L = [D() for _ in range(2^12)]    # long time
sage: mean_L = sum(L) / len(L)          # long time
sage: norm(mean_L.n() - D.c()) < 0.25   # long time
True

Python

>>> from sage.all import *
>>> D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ**Integer(3), RealNumber('3.0'), c=(Integer(1),Integer(0),Integer(0)))
>>> L = [D() for _ in range(Integer(2)**Integer(12))]
>>> mean_L = sum(L) / len(L)
>>> norm(mean_L.n() - D.c()) < RealNumber('0.25')
True

>>> D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ**Integer(3), RealNumber('3.0'), c=(Integer(1)/Integer(2),Integer(0),Integer(0)))
>>> L = [D() for _ in range(Integer(2)**Integer(12))]    # long time
>>> mean_L = sum(L) / len(L)          # long time
>>> norm(mean_L.n() - D.c()) < RealNumber('0.25')   # long time
True

add_offline_samples(cnt=1)[source]#

Precompute samples from \(B^{-1}D_1\) to be used in _call_non_spherical().

EXAMPLES:

Sage

sage: Sigma = Matrix([[5, -2, 4], [-2, 10, -5], [4, -5, 5]])
sage: D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ^3, Sigma)
sage: assert not D.is_spherical
sage: D.add_offline_samples(2^12)
sage: L = [D() for _ in range(2^12)] # Takes less time

Python

>>> from sage.all import *
>>> Sigma = Matrix([[Integer(5), -Integer(2), Integer(4)], [-Integer(2), Integer(10), -Integer(5)], [Integer(4), -Integer(5), Integer(5)]])
>>> D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ**Integer(3), Sigma)
>>> assert not D.is_spherical
>>> D.add_offline_samples(Integer(2)**Integer(12))
>>> L = [D() for _ in range(Integer(2)**Integer(12))] # Takes less time

c()[source]#

Center \(c\).

Samples from this sampler will be centered at \(c\).

EXAMPLES:

Sage

sage: D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ^3, 3.0, c=(1,0,0)); D
Discrete Gaussian sampler with Gaussian parameter σ = 3.00000000000000, c=(1, 0, 0) over lattice with basis

[1 0 0]
[0 1 0]
[0 0 1]

sage: D.c()
(1, 0, 0)

Python

>>> from sage.all import *
>>> D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ**Integer(3), RealNumber('3.0'), c=(Integer(1),Integer(0),Integer(0))); D
Discrete Gaussian sampler with Gaussian parameter σ = 3.00000000000000, c=(1, 0, 0) over lattice with basis
<BLANKLINE>
[1 0 0]
[0 1 0]
[0 0 1]

>>> D.c()
(1, 0, 0)

static compute_precision(precision, sigma)[source]#

Compute precision to use.

INPUT:

precision – an integer \(>= 53\) nor None.
sigma – if precision is None then the precision of sigma is used.

EXAMPLES:

Sage

sage: DGL = distributions.DiscreteGaussianDistributionLatticeSampler
sage: DGL.compute_precision(100, RR(3))
100
sage: DGL.compute_precision(100, RealField(200)(3))
100
sage: DGL.compute_precision(100, 3)
100
sage: DGL.compute_precision(None, RR(3))
53
sage: DGL.compute_precision(None, RealField(200)(3))
200
sage: DGL.compute_precision(None, 3)
53

Python

>>> from sage.all import *
>>> DGL = distributions.DiscreteGaussianDistributionLatticeSampler
>>> DGL.compute_precision(Integer(100), RR(Integer(3)))
100
>>> DGL.compute_precision(Integer(100), RealField(Integer(200))(Integer(3)))
100
>>> DGL.compute_precision(Integer(100), Integer(3))
100
>>> DGL.compute_precision(None, RR(Integer(3)))
53
>>> DGL.compute_precision(None, RealField(Integer(200))(Integer(3)))
200
>>> DGL.compute_precision(None, Integer(3))
53

f(x)[source]#

Compute the Gaussian \(\rho_{\Lambda, c, \Sigma}\).

EXAMPLES:

Sage

sage: Sigma = Matrix(ZZ, [[5, -2, 4], [-2, 10, -5], [4, -5, 5]])
sage: D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ^3, Sigma)
sage: D.f([1, 0, 1])
0.802518797962478
sage: D.f([1, 0, 3])
0.00562800641440405

Python

>>> from sage.all import *
>>> Sigma = Matrix(ZZ, [[Integer(5), -Integer(2), Integer(4)], [-Integer(2), Integer(10), -Integer(5)], [Integer(4), -Integer(5), Integer(5)]])
>>> D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ**Integer(3), Sigma)
>>> D.f([Integer(1), Integer(0), Integer(1)])
0.802518797962478
>>> D.f([Integer(1), Integer(0), Integer(3)])
0.00562800641440405

set_c(c)[source]#

Modifies center \(c\).

EXAMPLES:

Sage

sage: D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ^3, 3.0, c=(1,0,0))
sage: D.set_c((2, 0, 0))
sage: D
Discrete Gaussian sampler with Gaussian parameter σ = 3.00000000000000, c=(2, 0, 0) over lattice with basis

[1 0 0]
[0 1 0]
[0 0 1]

Python

>>> from sage.all import *
>>> D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ**Integer(3), RealNumber('3.0'), c=(Integer(1),Integer(0),Integer(0)))
>>> D.set_c((Integer(2), Integer(0), Integer(0)))
>>> D
Discrete Gaussian sampler with Gaussian parameter σ = 3.00000000000000, c=(2, 0, 0) over lattice with basis
<BLANKLINE>
[1 0 0]
[0 1 0]
[0 0 1]

sigma()[source]#

Gaussian parameter \(\sigma\).

If \(\sigma\) is a real number, samples from this sampler will have expected norm \(\sqrt{n}\sigma\) where \(n\) is the dimension of the lattice.

EXAMPLES:

Sage

sage: D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ^3, 3.0, c=(1,0,0))
sage: D.sigma()
3.00000000000000

Python

>>> from sage.all import *
>>> D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ**Integer(3), RealNumber('3.0'), c=(Integer(1),Integer(0),Integer(0)))
>>> D.sigma()
3.00000000000000