API¶

Implementation of direct approaches to sample from multivariate Gaussian distributions.

See also

pygauss.direct_sampling.sampler_band(mu, A, b, mode='precision', seed=None, size=1)[source]¶

Algorithm dedicated to sample from a multivariate real-valued Gaussian distribution \(\mathcal{N}(\boldsymbol{\mu},\mathbf{A})\) or \(\mathcal{N}(\boldsymbol{\mu},\mathbf{A}^{-1})\) when \(\mathbf{A}\) is a band matrix.

Parameters:	mu (1-D array_like, of length d) – Mean of the d-dimensional Gaussian distribution. A (2-D array_like, of shape (d, d)) – Covariance or precision matrix of the distribution. It must be symmetric and positive-definite for proper sampling. b (int) – Bandwidth of A. mode (string, optional) – Indicates if A refers to the precision or covariance matrix of the Gaussian distribution. seed (int, optional) – Random seed to reproduce experimental results. size (int, optional) – Given a size of for instance T, T independent and identically distributed (i.i.d.) samples are returned.
Returns:	theta – The drawn samples, of shape (d,size), if that was provided. If not, the shape is (d,1).
Return type:	ndarray, of shape (d,size)
Raises:	`ValueError` – If mode is not included in [‘covariance’,’precision’].

Examples

>>> d = 2
>>> mu = np.zeros(d)
>>> A = np.eye(2)
>>> b = 0
>>> mode = "covariance"
>>> size = 1
>>> theta = sampler_band(mu,A,b,mode=mode,seed=2022,size=size)

pygauss.direct_sampling.sampler_circulant(mu, a, M, N, mode='precision', seed=None, size=1)[source]¶

Algorithm dedicated to sample from a multivariate real-valued Gaussian distribution \(\mathcal{N}(\boldsymbol{\mu},\mathbf{A})\) or \(\mathcal{N}(\boldsymbol{\mu},\mathbf{A}^{-1})\) when \(\mathbf{A}\) is a block-circulant matrix with circulant blocks.

Parameters:	mu (1-D array_like, of length d) – Mean of the d-dimensional Gaussian distribution. a (2-D array_like, of shape (N, M)) – Vector built by stacking the first columns associated to the \(M\) blocks of size \(N\) of the matrix \(\mathbf{A}\). M (int) – Number of different blocks. N (int) – Dimension of each block. mode (string, optional) – Indicates if A refers to the precision or covariance matrix of the Gaussian distribution. seed (int, optional) – Random seed to reproduce experimental results. size (int, optional) – Given a size of for instance T, T independent and identically distributed (i.i.d.) samples are returned.
Returns:	theta – The drawn samples, of shape (d,size), if that was provided. If not, the shape is (d,1).
Return type:	ndarray, of shape (d,size)
Raises:	`ValueError` – If mode is not included in [‘covariance’,’precision’].

Examples

>>> d = 2
>>> mu = np.zeros(d)
>>> a = np.matrix([1,0]).T
>>> M = 1
>>> N = 2
>>> mode = "covariance"
>>> size = 1
>>> theta = sampler_circulant(mu,a,M,N,mode=mode,seed=2022,size=size)

pygauss.direct_sampling.sampler_factorization(mu, A, mode='precision', method='Cholesky', seed=None, size=1)[source]¶

Algorithm dedicated to sample from a multivariate real-valued Gaussian distribution \(\mathcal{N}(\boldsymbol{\mu},\mathbf{A})\) or \(\mathcal{N}(\boldsymbol{\mu},\mathbf{A}^{-1})\) based on matrix factorization (e.g., Cholesky or square root).

Parameters:	mu (1-D array_like, of length d) – Mean of the d-dimensional Gaussian distribution. A (2-D array_like, of shape (d, d)) – Covariance or precision matrix of the distribution. It must be symmetric and positive-definite for proper sampling. mode (string, optional) – Indicates if A refers to the precision or covariance matrix of the Gaussian distribution. method (string, optional) – Factorization method. Choose either ‘Cholesky’ or ‘square-root’. seed (int, optional) – Random seed to reproduce experimental results. size (int, optional) – Given a size of for instance T, T independent and identically distributed (i.i.d.) samples are returned.
Returns:	theta – The drawn samples, of shape (d,size), if that was provided. If not, the shape is (d,1).
Return type:	ndarray, of shape (d,size)
Raises:	`ValueError` – If A is not positive definite and symmetric. If mode is not included in [‘covariance’,’precision’]. If method is not included in [‘Cholesky’,’square-root’].

Examples

>>> d = 2
>>> mu = np.zeros(d)
>>> A = np.eye(d)
>>> mode = "covariance"
>>> method = "Cholesky"
>>> size = 1
>>> theta = sampler_factorization(mu,A,mode=mode,method=method,seed=2022,size=size)

pygauss.direct_sampling.sampler_squareRootApprox(mu, A, lam_l, lam_u, tol, K=100, mode='precision', seed=None, size=1, info=False)[source]¶

Algorithm dedicated to sample from a multivariate real-valued Gaussian distribution \(\mathcal{N}(\boldsymbol{\mu},\mathbf{A})\) or \(\mathcal{N}(\boldsymbol{\mu},\mathbf{A}^{-1})\) based on matrix square root approximation using Chebychev polynomials.

Parameters:	mu (1-D array_like, of length d) – Mean of the d-dimensional Gaussian distribution. A (function) – Linear operator returning the matrix-vector product \(\mathbf{Ax}\) where \(\mathbf{x}) \in \mathbb{R}^d\). lam_l (float) – Lower bound on the eigenvalues of \(\mathbf{A}\). lam_u (float) – Upper bound on the eigenvalues of \(\mathbf{A}\). tol (float) – Tolerance threshold used to optimize the polynomial order \(K\). This threshold stands for the Euclidean distance between the vector computed using order \(K\) and the one computed using order \(L\)\(\leq\)\(K\). K (int, optional) – Polynomial order of the approximation. mode (string, optional) – Indicates if A refers to the precision or covariance matrix of the Gaussian distribution. seed (int, optional) – Random seed to reproduce experimental results. size (int, optional) – Given a size of for instance T, T independent and identically distributed (i.i.d.) samples are returned. info (boolean, optional) – If info is True, returns the order \(K\) used in the polynomial approximation.
Returns:	theta – The drawn samples, of shape (d,size), if that was provided. If not, the shape is (d,1).
Return type:	ndarray, of shape (d,size)
Raises:	`ValueError` – If mode is not included in [‘covariance’,’precision’].

Examples

>>> d = 2
>>> mu = np.zeros(d)
>>> def A(x):
    return np.eye(d).dot(x)
>>> lam_l = 0
>>> lam_u = 1
>>> tol = 1e-4
>>> mode = "covariance"
>>> size = 1
>>> theta = sampler_squareRootApprox(mu,A,lam_l=lam_l,lam_u=lam_u,tol=tol,
mode=mode,seed=2022,size=size)

pygauss.direct_sampling.sampler_CG(mu, A, K, init, tol=0.0001, mode='precision', seed=None, size=1, info=False)[source]¶

Algorithm dedicated to sample from a multivariate real-valued Gaussian distribution \(\mathcal{N}(\boldsymbol{\mu},\mathbf{A})\) or \(\mathcal{N}(\boldsymbol{\mu},\mathbf{A}^{-1})\) based on the conjugate gradient algorithm.

Parameters:	mu (1-D array_like, of length d) – Mean of the d-dimensional Gaussian distribution. A (function) – Linear operator returning the matrix-vector product \(\mathbf{Ax}\) where \(\mathbf{x} \in \mathbb{R}^d\). K (int, optional) – Number of conjugate gradient iterations. init (1-D array_like, of length d) – Vector used to initialize the CG sampler. tol (float, optional) – Tolerance threshold used to stop the conjugate gradient sampler. mode (string, optional) – Indicates if A refers to the precision or covariance matrix of the Gaussian distribution. seed (int, optional) – Random seed to reproduce experimental results. size (int, optional) – Given a size of for instance T, T independent and identically distributed (i.i.d.) samples are returned. info (boolean, optional) – If info is True, returns the number of iterations \(K\).
Returns:	theta – The drawn samples, of shape (d,size), if that was provided. If not, the shape is (d,1).
Return type:	ndarray, of shape (d,size)
Raises:	`ValueError` – If mode is not included in [‘covariance’,’precision’].

Examples

>>> d = 2
>>> mu = np.zeros(d)
>>> def A(x):
    return np.eye(d).dot(x)
>>> K = 2
>>> init = mu
>>> theta = sampler_CG(mu,A,K,init)

class pygauss.direct_sampling.sampler_PO(mu1, mu2, K, init, tol=0.0001, seed=None, size=1)[source]¶

Algorithm dedicated to sample from a multivariate real-valued Gaussian distribution \(\mathcal{N}(\boldsymbol{\mu},\mathbf{Q}^{-1})\) where \(\mathbf{Q}\) is a symmetric and positive definite precision matrix. We assume here that \(\mathbf{Q} = \mathbf{G}_1^T\mathbf{\Lambda}_1^{-1}\mathbf{G}_1 + \mathbf{G}_2^T\mathbf{\Lambda}_2^{-1}\mathbf{G}_2\). The mean vector is assumed to have the form \(\boldsymbol{\mu} = \mathbf{G}_1^T\mathbf{\Lambda}_1^{-1}\boldsymbol{\mu}_1 + \mathbf{G}_2^T\mathbf{\Lambda}_2^{-1}\boldsymbol{\mu}_2\). Sampling from the corresponding multivariate Gaussian distribution is done with the perturbation-optimization sampler.

__init__(mu1, mu2, K, init, tol=0.0001, seed=None, size=1)[source]¶

Parameters:

mu1 (1-D array_like, of length d) –
mu2 (1-D array_like, of length d) –
K (int, optional) – Number of conjugate gradient iterations to solve the linear system \(\mathbf{Q}\boldsymbol{\theta} = \boldsymbol{\eta}\).
init (1-D array_like, of length d) – Vector used to initialize the CG algorithm.
tol (float, optional) – Tolerance threshold used to stop the conjugate gradient algorithm.
seed (int, optional) – Random seed to reproduce experimental results.
size (int, optional) – Given a size of for instance T, T independent and identically distributed (i.i.d.) samples are returned.

circu_diag_band(Lamb1, g, M, N, Q2, b2)[source]¶

We assume here that \(\mathbf{G}_1\) is a circulant matrix, \(\mathbf{\Lambda}_1\) is diagonal, \(\mathbf{G}_2\) is the identity matrix and \(\mathbf{Q}_2 = \mathbf{\Lambda}_2^{-1}\) is a band matrix.

Parameters:	Lamb1 (1-D array_like, of length d) – Diagonal elements of \(\mathbf{\Lambda}_1\). g (2-D array_like, of shape (N, M)) – Vector built by stacking the first columns associated to the \(M\) blocks of size \(N\) of the matrix \(\mathbf{G}_1\). M (int) – Number of different blocks in \(\mathbf{G}_1\). N (int) – Dimension of each block in \(\mathbf{G}_1\). Q2 (2-D array_like, of shape (d, d)) – Precision matrix \(\mathbf{Q}_2\). b2 (int) – Bandwidth of \(\mathbf{Q}_2\).
Returns:	theta – The drawn samples, of shape (d,size), if that was provided. If not, the shape is (d,1).
Return type:	ndarray, of shape (d,size)

Examples

>>> d = 15
>>> mu1 = np.zeros(d)
>>> mu2 = np.zeros(d)
>>> K = 15
>>> init = np.zeros(d)
>>> Lamb1 = np.random.normal(2,0.1,d)
>>> g =  np.reshape(np.random.normal(2,0.1,d),(d,1))
>>> M = 1
>>> N = d
>>> Q2 = np.diag(np.random.normal(2,0.1,d))
>>> b2 = 0
>>> size = 10000
>>> S = sampler_PO(mu1,mu2,K,init,size=10000)
>>> theta = S.circu_diag_band(Lamb1,g,M,N,Q21,b2)