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: 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

>>> 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: 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]

>>> 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: 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

>>> 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:

__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 – boolean (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: 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)

>>> 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: 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]

>>> 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: 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...

>>> 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: 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)

>>> 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: 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

>>> 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: 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]

>>> 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: 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

>>> 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: 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

>>> 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: 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

sage: import numpy
sage: M = matrix(ZZ, [[1,2],[0,1]])
sage: D = distributions.DiscreteGaussianDistributionLatticeSampler(M, 20.0)
sage: L = [D() for _ in range(2^12)]                                               # long time
sage: div = numpy.mean([abs(x) for x,y in L]) / numpy.mean([abs(y) for x,y, in L]) # long time
sage: 0.9 < div < 1.1                                                              # long time
True

>>> 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

>>> import numpy
>>> M = matrix(ZZ, [[Integer(1),Integer(2)],[Integer(0),Integer(1)]])
>>> D = distributions.DiscreteGaussianDistributionLatticeSampler(M, RealNumber('20.0'))
>>> L = [D() for _ in range(Integer(2)**Integer(12))]                                               # long time
>>> div = numpy.mean([abs(x) for x,y in L]) / numpy.mean([abs(y) for x,y, in L]) # long time
>>> RealNumber('0.9') < div < RealNumber('1.1')                                                              # 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: 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

>>> 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: 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)

>>> 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 – integer \(>= 53\) nor None

  • sigma – if precision is None then the precision of sigma is used

EXAMPLES:

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

>>> 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: 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

>>> 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: 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]

>>> 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: D = distributions.DiscreteGaussianDistributionLatticeSampler(ZZ^3, 3.0, c=(1,0,0))
sage: D.sigma()
3.00000000000000

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