Source code for estimator.nd

# -*- coding: utf-8 -*-

from dataclasses import dataclass

from sage.all import binomial, ceil, exp, log, oo, parent, pi, RealField, RR, sqrt


[docs] def stddevf(sigma): """ Gaussian width parameter σ → standard deviation. :param sigma: Gaussian width parameter σ EXAMPLE:: >>> from estimator.nd import stddevf >>> stddevf(64.0) 25.532... >>> stddevf(64) 25.532... >>> stddevf(RealField(256)(64)).prec() 256 """ try: prec = parent(sigma).prec() except AttributeError: prec = 0 if prec > 0: FF = parent(sigma) else: FF = RR return FF(sigma) / FF(sqrt(2 * pi))
[docs] def sigmaf(stddev): """ Standard deviation → Gaussian width parameter σ. :param stddev: standard deviation EXAMPLE:: >>> from estimator.nd import stddevf, sigmaf >>> n = 64.0 >>> sigmaf(stddevf(n)) 64.000... >>> sigmaf(RealField(128)(1.0)) 2.5066282746310005024157652848110452530 >>> sigmaf(1.0) 2.506628274631... >>> sigmaf(1) 2.506628274631... """ RR = parent(stddev) # check that we got ourselves a real number type try: if abs(RR(0.5) - 0.5) > 0.001: RR = RealField(53) # hardcode something except TypeError: RR = RealField(53) # hardcode something return RR(sqrt(2 * pi)) * stddev
[docs] @dataclass class NoiseDistribution: """ All noise distributions are instances of this class. """ # cut-off for Gaussian distributions gaussian_tail_bound = 2 # probability that a coefficient falls within the cut-off gaussian_tail_prob = 1 - 2 * exp(-4 * pi) stddev: float mean: float = 0 n: int = None bounds: tuple = (None, None) density: float = 1.0 # Hamming weight / dimension. tag: str = "" def __lt__(self, other): """ We compare distributions by comparing their standard deviation. EXAMPLE:: >>> from estimator.nd import NoiseDistribution as ND >>> ND.DiscreteGaussian(2.0) < ND.CenteredBinomial(18) True >>> ND.DiscreteGaussian(3.0) < ND.CenteredBinomial(18) False >>> ND.DiscreteGaussian(4.0) < ND.CenteredBinomial(18) False """ try: return self.stddev < other.stddev except AttributeError: return self.stddev < other def __le__(self, other): """ We compare distributions by comparing their standard deviation. EXAMPLE:: >>> from estimator.nd import NoiseDistribution as ND >>> ND.DiscreteGaussian(2.0) <= ND.CenteredBinomial(18) True >>> ND.DiscreteGaussian(3.0) <= ND.CenteredBinomial(18) True >>> ND.DiscreteGaussian(4.0) <= ND.CenteredBinomial(18) False """ try: return self.stddev <= other.stddev except AttributeError: return self.stddev <= other def __str__(self): """ EXAMPLE:: >>> from estimator.nd import NoiseDistribution as ND >>> ND.DiscreteGaussianAlpha(0.01, 7681) D(σ=30.64) """ if self.n: return f"D(σ={float(self.stddev):.2f}, μ={float(self.mean):.2f}, n={int(self.n)})" else: return f"D(σ={float(self.stddev):.2f}, μ={float(self.mean):.2f})" def __repr__(self): if self.mean == 0.0: return f"D(σ={float(self.stddev):.2f})" else: return f"D(σ={float(self.stddev):.2f}, μ={float(self.mean):.2f})" def __hash__(self): """ EXAMPLE:: >>> from estimator.nd import NoiseDistribution as ND >>> hash(ND(3.0, 1.0)) == hash((3.0, 1.0, None)) True """ return hash((self.stddev, self.mean, self.n)) def __len__(self): """ EXAMPLE:: >>> from estimator.nd import NoiseDistribution as ND >>> D = ND.SparseTernary(1024, p=128, m=128) >>> len(D) 1024 >>> int(round(len(D) * float(D.density))) 256 """ if self.n is not None: return self.n else: raise ValueError("Distribution has no length.") @property def is_Gaussian_like(self): return ("Gaussian" in self.tag) or ("CenteredBinomial" in self.tag) @property def is_bounded(self): return (self.bounds[1] - self.bounds[0]) < oo @property def is_sparse(self): """ We consider a distribution "sparse" if its density is < 1/2. """ # NOTE: somewhat arbitrary return self.density < 0.5
[docs] def support_size(self, n=None, fraction=1.0): """ Compute the size of the support covering the probability given as fraction. EXAMPLE:: >>> from estimator.nd import NoiseDistribution as ND >>> D = ND.Uniform(-3,3, 64) >>> D.support_size(fraction=.99) 1207562882759477428726191443614714994252339953407098880 >>> D = ND.SparseTernary(64, 8) >>> D.support_size() 32016101348447354880 """ if not n: if not self.n: raise ValueError(f"Length required to determine support size, but n was {n}.") n = self.n if "SparseTernary" in self.tag: h = self.get_hamming_weight(n) # TODO: this is assuming that the non-zero entries are uniform over {-1,1} # need p and m for more accurate calculation size = 2**h * binomial(n, h) * RR(fraction) elif self.is_bounded: # TODO: this might be suboptimal/inaccurate for binomial distribution a, b = self.bounds size = RR(fraction) * (b - a + 1) ** n else: # Looks like nd is Gaussian # -> we'll treat it as bounded (with failure probability) t = self.gaussian_tail_bound p = self.gaussian_tail_prob if p**n < fraction: raise NotImplementedError( f"TODO(nd.support-size): raise t. {RR(p ** n)}, {n}, {fraction}" ) b = 2 * t * sigmaf(self.stddev) + 1 return (2 * b + 1) ** n return ceil(size)
[docs] def get_hamming_weight(self, n=None): if not n: if not self.n: raise ValueError("Length required to determine hamming weight.") n = self.n return round(n * float(self.density))
[docs] @staticmethod def DiscreteGaussian(stddev, mean=0, n=None): """ A discrete Gaussian distribution with standard deviation ``stddev`` per component. EXAMPLE:: >>> from estimator.nd import NoiseDistribution as ND >>> ND.DiscreteGaussian(3.0, 1.0) D(σ=3.00, μ=1.00) """ b_val = oo if n is None else ceil(log(n, 2) * stddev) return NoiseDistribution( stddev=RR(stddev), mean=RR(mean), n=n, bounds=(-b_val, b_val), density=1 - min(RR(1 / (sqrt(2 * pi) * stddev)), 1.0), tag="DiscreteGaussian", )
[docs] @staticmethod def DiscreteGaussianAlpha(alpha, q, mean=0, n=None): """ A discrete Gaussian distribution with standard deviation α⋅q/√(2π) per component. EXAMPLE:: >>> from estimator.nd import NoiseDistribution as ND >>> ND.DiscreteGaussianAlpha(0.001, 2048) D(σ=0.82) """ stddev = stddevf(alpha * q) return NoiseDistribution.DiscreteGaussian(stddev=RR(stddev), mean=RR(mean), n=n)
[docs] @staticmethod def CenteredBinomial(eta, n=None): """ Sample a_1, …, a_η, b_1, …, b_η and return Σ(a_i - b_i). EXAMPLE:: >>> from estimator.nd import NoiseDistribution as ND >>> ND.CenteredBinomial(8) D(σ=2.00) """ stddev = sqrt(eta / 2.0) return NoiseDistribution( stddev=RR(stddev), density=1 - binomial(2 * eta, eta) * 2 ** (-2 * eta), mean=RR(0), n=n, bounds=(-eta, eta), tag="CenteredBinomial", )
[docs] @staticmethod def Uniform(a, b, n=None): """ Uniform distribution ∈ ``[a,b]``, endpoints inclusive. EXAMPLE:: >>> from estimator.nd import NoiseDistribution as ND >>> ND.Uniform(-3, 3) D(σ=2.00) >>> ND.Uniform(-4, 3) D(σ=2.29, μ=-0.50) """ if b < a: raise ValueError(f"upper limit must be larger than lower limit but got: {b} < {a}") m = b - a + 1 mean = (a + b) / RR(2) stddev = sqrt((m**2 - 1) / RR(12)) if a <= 0 and b >= 0: density = 1.0 - 1.0 / m else: density = 0.0 return NoiseDistribution( n=n, stddev=stddev, mean=mean, bounds=(a, b), density=density, tag="Uniform" )
[docs] @staticmethod def UniformMod(q, n=None): """ Uniform mod ``q``, with balanced representation. EXAMPLE:: >>> from estimator.nd import NoiseDistribution as ND >>> ND.UniformMod(7) D(σ=2.00) >>> ND.UniformMod(8) D(σ=2.29, μ=-0.50) """ a = -(q // 2) b = -a -1 if q % 2 == 0 else -a return NoiseDistribution.Uniform(a, b, n=n)
[docs] @staticmethod def SparseTernary(n, p, m=None): """ Distribution of vectors of length ``n`` with ``p`` entries of 1 and ``m`` entries of -1, rest 0. EXAMPLE:: >>> from estimator.nd import NoiseDistribution as ND >>> ND.SparseTernary(100, p=10) D(σ=0.45) >>> ND.SparseTernary(100, p=10, m=10) D(σ=0.45) >>> ND.SparseTernary(100, p=10, m=8) D(σ=0.42, μ=0.02) """ if m is None: m = p if n == 0: # this might happen in the dual attack return NoiseDistribution( stddev=0, mean=0, density=0, bounds=(-1, 1), tag="SparseTernary", n=0 ) mean = RR(p / n - m / n) stddev = sqrt(p / n * (1 - mean)**2 + m / n * (-1 - mean)**2 + (n - (p + m)) / n * (mean)**2) density = RR((p + m) / n) return NoiseDistribution( stddev=stddev, mean=mean, density=density, bounds=(-1, 1), tag="SparseTernary", n=n )