# Miscellaneous Ops¶

The pyro.ops module implements tensor utilities that are mostly independent of the rest of Pyro.

## Utilities for HMC¶

class DualAveraging(prox_center=0, t0=10, kappa=0.75, gamma=0.05)[source]

Bases: object

Dual Averaging is a scheme to solve convex optimization problems. It belongs to a class of subgradient methods which uses subgradients to update parameters (in primal space) of a model. Under some conditions, the averages of generated parameters during the scheme are guaranteed to converge to an optimal value. However, a counter-intuitive aspect of traditional subgradient methods is “new subgradients enter the model with decreasing weights” (see $$$$). Dual Averaging scheme solves that phenomenon by updating parameters using weights equally for subgradients (which lie in a dual space), hence we have the name “dual averaging”.

This class implements a dual averaging scheme which is adapted for Markov chain Monte Carlo (MCMC) algorithms. To be more precise, we will replace subgradients by some statistics calculated during an MCMC trajectory. In addition, introducing some free parameters such as t0 and kappa is helpful and still guarantees the convergence of the scheme.

References

 Primal-dual subgradient methods for convex problems, Yurii Nesterov

 The No-U-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo, Matthew D. Hoffman, Andrew Gelman

Parameters: prox_center (float) – A “prox-center” parameter introduced in $$$$ which pulls the primal sequence towards it. t0 (float) – A free parameter introduced in $$$$ that stabilizes the initial steps of the scheme. kappa (float) – A free parameter introduced in $$$$ that controls the weights of steps of the scheme. For a small kappa, the scheme will quickly forget states from early steps. This should be a number in $$(0.5, 1]$$. gamma (float) – A free parameter which controls the speed of the convergence of the scheme.
reset()[source]
step(g)[source]

Updates states of the scheme given a new statistic/subgradient g.

Parameters: g (float) – A statistic calculated during an MCMC trajectory or subgradient.
get_state()[source]

Returns the latest $$x_t$$ and average of $$\left\{x_i\right\}_{i=1}^t$$ in primal space.

velocity_verlet(z, r, potential_fn, inverse_mass_matrix, step_size, num_steps=1, z_grads=None)[source]

Second order symplectic integrator that uses the velocity verlet algorithm.

Parameters: Return tuple (z_next, r_next, z_grads, potential_energy): z (dict) – dictionary of sample site names and their current values (type Tensor). r (dict) – dictionary of sample site names and corresponding momenta (type Tensor). potential_fn (callable) – function that returns potential energy given z for each sample site. The negative gradient of the function with respect to z determines the rate of change of the corresponding sites’ momenta r. inverse_mass_matrix (torch.Tensor) – a tensor $$M^{-1}$$ which is used to calculate kinetic energy: $$E_{kinetic} = \frac{1}{2}z^T M^{-1} z$$. Here $$M$$ can be a 1D tensor (diagonal matrix) or a 2D tensor (dense matrix). step_size (float) – step size for each time step iteration. num_steps (int) – number of discrete time steps over which to integrate. z_grads (torch.Tensor) – optional gradients of potential energy at current z. next position and momenta, together with the potential energy and its gradient w.r.t. z_next.
potential_grad(potential_fn, z)[source]

Gradient of potential_fn w.r.t. parameters z.

Parameters: potential_fn – python callable that takes in a dictionary of parameters and returns the potential energy. z (dict) – dictionary of parameter values keyed by site name. tuple of (z_grads, potential_energy), where z_grads is a dictionary with the same keys as z containing gradients and potential_energy is a torch scalar.
class WelfordCovariance(diagonal=True)[source]

Bases: object

Implements Welford’s online scheme for estimating (co)variance (see $$$$). Useful for adapting diagonal and dense mass structures for HMC.

References

 The Art of Computer Programming, Donald E. Knuth

reset()[source]
update(sample)[source]
get_covariance(regularize=True)[source]

## Newton Optimizers¶

newton_step(loss, x, trust_radius=None)[source]

Performs a Newton update step to minimize loss on a batch of variables, optionally constraining to a trust region .

This is especially usful because the final solution of newton iteration is differentiable wrt the inputs, even when all but the final x is detached, due to this method’s quadratic convergence . loss must be twice-differentiable as a function of x. If loss is 2+d-times differentiable, then the return value of this function is d-times differentiable.

When loss is interpreted as a negative log probability density, then the return values mode,cov of this function can be used to construct a Laplace approximation MultivariateNormal(mode,cov).

Warning

Take care to detach the result of this function when used in an optimization loop. If you forget to detach the result of this function during optimization, then backprop will propagate through the entire iteration process, and worse will compute two extra derivatives for each step.

Example use inside a loop:

x = torch.zeros(1000, 2)  # arbitrary initial value
for step in range(100):
x = x.detach()          # block gradients through previous steps
x.requires_grad = True  # ensure loss is differentiable wrt x
loss = my_loss_function(x)
x = newton_step(loss, x, trust_radius=1.0)
# the final x is still differentiable

 Yuan, Ya-xiang. Iciam. Vol. 99. 2000.
“A review of trust region algorithms for optimization.” ftp://ftp.cc.ac.cn/pub/yyx/papers/p995.pdf
 Christianson, Bruce. Optimization Methods and Software 3.4 (1994)
“Reverse accumulation and attractive fixed points.” http://uhra.herts.ac.uk/bitstream/handle/2299/4338/903839.pdf
Parameters: loss (torch.Tensor) – A scalar function of x to be minimized. x (torch.Tensor) – A dependent variable of shape (N, D) where N is the batch size and D is a small number. trust_radius (float) – An optional trust region trust_radius. The updated value mode of this function will be within trust_radius of the input x. A pair (mode, cov) where mode is an updated tensor of the same shape as the original value x, and cov is an esitmate of the covariance DxD matrix with cov.shape == x.shape[:-1] + (D,D). tuple
newton_step_1d(loss, x, trust_radius=None)[source]

Performs a Newton update step to minimize loss on a batch of 1-dimensional variables, optionally regularizing to constrain to a trust region.

See newton_step() for details.

Parameters: loss (torch.Tensor) – A scalar function of x to be minimized. x (torch.Tensor) – A dependent variable with rightmost size of 1. trust_radius (float) – An optional trust region trust_radius. The updated value mode of this function will be within trust_radius of the input x. A pair (mode, cov) where mode is an updated tensor of the same shape as the original value x, and cov is an esitmate of the covariance 1x1 matrix with cov.shape == x.shape[:-1] + (1,1). tuple
newton_step_2d(loss, x, trust_radius=None)[source]

Performs a Newton update step to minimize loss on a batch of 2-dimensional variables, optionally regularizing to constrain to a trust region.

See newton_step() for details.

Parameters: loss (torch.Tensor) – A scalar function of x to be minimized. x (torch.Tensor) – A dependent variable with rightmost size of 2. trust_radius (float) – An optional trust region trust_radius. The updated value mode of this function will be within trust_radius of the input x. A pair (mode, cov) where mode is an updated tensor of the same shape as the original value x, and cov is an esitmate of the covariance 2x2 matrix with cov.shape == x.shape[:-1] + (2,2). tuple
newton_step_3d(loss, x, trust_radius=None)[source]

Performs a Newton update step to minimize loss on a batch of 3-dimensional variables, optionally regularizing to constrain to a trust region.

See newton_step() for details.

Parameters: loss (torch.Tensor) – A scalar function of x to be minimized. x (torch.Tensor) – A dependent variable with rightmost size of 2. trust_radius (float) – An optional trust region trust_radius. The updated value mode of this function will be within trust_radius of the input x. A pair (mode, cov) where mode is an updated tensor of the same shape as the original value x, and cov is an esitmate of the covariance 3x3 matrix with cov.shape == x.shape[:-1] + (3,3). tuple

## Tensor Indexing¶

vindex(tensor, args)[source]

See also the convenience wrapper Vindex.

This is useful for writing indexing code that is compatible with batching and enumeration, especially for selecting mixture components with discrete random variables.

For example suppose x is a parameter with x.dim() == 3 and we wish to generalize the expression x[i, :, j] from integer i,j to tensors i,j with batch dims and enum dims (but no event dims). Then we can write the generalize version using Vindex

xij = Vindex(x)[i, :, j]

batch_shape = broadcast_shape(i.shape, j.shape)
event_shape = (x.size(1),)
assert xij.shape == batch_shape + event_shape


To handle the case when x may also contain batch dimensions (e.g. if x was sampled in a plated context as when using vectorized particles), vindex() uses the special convention that Ellipsis denotes batch dimensions (hence ... can appear only on the left, never in the middle or in the right). Suppose x has event dim 3. Then we can write:

old_batch_shape = x.shape[:-3]
old_event_shape = x.shape[-3:]

xij = Vindex(x)[..., i, :, j]   # The ... denotes unknown batch shape.

new_batch_shape = broadcast_shape(old_batch_shape, i.shape, j.shape)
new_event_shape = (x.size(1),)
assert xij.shape = new_batch_shape + new_event_shape


Note that this special handling of Ellipsis differs from the NEP .

Formally, this function assumes:

1. Each arg is either Ellipsis, slice(None), an integer, or a batched torch.LongTensor (i.e. with empty event shape). This function does not support Nontrivial slices or torch.BoolTensor masks. Ellipsis can only appear on the left as args.
2. If args is not Ellipsis then tensor is not batched, and its event dim is equal to len(args).
3. If args is Ellipsis then tensor is batched and its event dim is equal to len(args[1:]). Dims of tensor to the left of the event dims are considered batch dims and will be broadcasted with dims of tensor args.

Note that if none of the args is a tensor with .dim() > 0, then this function behaves like standard indexing:

if not any(isinstance(a, torch.Tensor) and a.dim() for a in args):
assert Vindex(x)[args] == x[args]


References

introduces vindex as a helper for vectorized indexing. The Pyro implementation is similar to the proposed notation x.vindex[] except for slightly different handling of Ellipsis.
Parameters: tensor (torch.Tensor) – A tensor to be indexed. args (tuple) – An index, as args to __getitem__. A nonstandard interpetation of tensor[args]. torch.Tensor
class Vindex(tensor)[source]

Bases: object

Convenience wrapper around vindex().

The following are equivalent:

Vindex(x)[..., i, j, :]
vindex(x, (Ellipsis, i, j, slice(None)))

Parameters: tensor (torch.Tensor) – A tensor to be indexed. An object with a special __getitem__() method.

## Tensor Contraction¶

contract_expression(equation, *shapes, **kwargs)[source]

Wrapper around opt_einsum.contract_expression() that optionally uses Pyro’s cheap optimizer and optionally caches contraction paths.

Parameters: cache_path (bool) – whether to cache the contraction path. Defaults to True.
contract(equation, *operands, **kwargs)[source]

Wrapper around opt_einsum.contract() that optionally uses Pyro’s cheap optimizer and optionally caches contraction paths.

Parameters: cache_path (bool) – whether to cache the contraction path. Defaults to True.
einsum(equation, *operands, **kwargs)[source]

Generalized plated sum-product algorithm via tensor variable elimination.

This generalizes contract() in two ways:

1. Multiple outputs are allowed, and intermediate results can be shared.
2. Inputs and outputs can be plated along symbols given in plates; reductions along plates are product reductions.

The best way to understand this function is to try the examples below, which show how einsum() calls can be implemented as multiple calls to contract() (which is generally more expensive).

To illustrate multiple outputs, note that the following are equivalent:

z1, z2, z3 = einsum('ab,bc->a,b,c', x, y)  # multiple outputs

z1 = contract('ab,bc->a', x, y)
z2 = contract('ab,bc->b', x, y)
z3 = contract('ab,bc->c', x, y)


To illustrate plated inputs, note that the following are equivalent:

assert len(x) == 3 and len(y) == 3
z = einsum('ab,ai,bi->b', w, x, y, plates='i')

z = contract('ab,a,a,a,b,b,b->b', w, *x, *y)


When a sum dimension a always appears with a plate dimension i, then a corresponds to a distinct symbol for each slice of a. Thus the following are equivalent:

assert len(x) == 3 and len(y) == 3
z = einsum('ai,ai->', x, y, plates='i')

z = contract('a,b,c,a,b,c->', *x, *y)


When such a sum dimension appears in the output, it must be accompanied by all of its plate dimensions, e.g. the following are equivalent:

assert len(x) == 3 and len(y) == 3
z = einsum('abi,abi->bi', x, y, plates='i')

z = torch.stack([z0, z1, z2])


Note that each plate slice through the output is multilinear in all plate slices through all inptus, thus e.g. batch matrix multiply would be implemented without plates, so the following are all equivalent:

xy = einsum('abc,acd->abd', x, y, plates='')
xy = torch.stack([xa.mm(ya) for xa, ya in zip(x, y)])
xy = torch.bmm(x, y)


Among all valid equations, some computations are polynomial in the sizes of the input tensors and other computations are exponential in the sizes of the input tensors. This function raises NotImplementedError whenever the computation is exponential.

Parameters: equation (str) – An einsum equation, optionally with multiple outputs. operands (torch.Tensor) – A collection of tensors. plates (str) – An optional string of plate symbols. backend (str) – An optional einsum backend, defaults to ‘torch’. cache (dict) – An optional shared_intermediates() cache. modulo_total (bool) – Optionally allow einsum to arbitrarily scale each result plate, which can significantly reduce computation. This is safe to set whenever each result plate denotes a nonnormalized probability distribution whose total is not of interest. a tuple of tensors of requested shape, one entry per output. tuple ValueError – if tensor sizes mismatch or an output requests a plated dim without that dim’s plates. NotImplementedError – if contraction would have cost exponential in the size of any input tensor.
ubersum(equation, *operands, **kwargs)[source]

Deprecated, use einsum() instead.

## Gaussian Contraction¶

class Gaussian(log_normalizer, info_vec, precision)[source]

Bases: object

Non-normalized Gaussian distribution.

This represents an arbitrary semidefinite quadratic function, which can be interpreted as a rank-deficient scaled Gaussian distribution. The precision matrix may have zero eigenvalues, thus it may be impossible to work directly with the covariance matrix.

Parameters: log_normalizer (torch.Tensor) – a normalization constant, which is mainly used to keep track of normalization terms during contractions. info_vec (torch.Tensor) – information vector, which is a scaled version of the mean info_vec = precision @ mean. We use this represention to make gaussian contraction fast and stable. precision (torch.Tensor) – precision matrix of this gaussian.
dim()[source]
batch_shape[source]
expand(batch_shape)[source]
reshape(batch_shape)[source]
__getitem__(index)[source]

Index into the batch_shape of a Gaussian.

static cat(parts, dim=0)[source]

Concatenate a list of Gaussians along a given batch dimension.

event_pad(left=0, right=0)[source]

Pad along event dimension.

event_permute(perm)[source]

Permute along event dimension.

__add__(other)[source]

Adds two Gaussians in log-density space.

log_density(value)[source]

Evaluate the log density of this Gaussian at a point value:

-0.5 * value.T @ precision @ value + value.T @ info_vec + log_normalizer


This is mainly used for testing.

condition(value)[source]

Condition this Gaussian on a trailing subset of its state. This should satisfy:

g.condition(y).dim() == g.dim() - y.size(-1)


Note that since this is a non-normalized Gaussian, we include the density of y in the result. Thus condition() is similar to a functools.partial binding of arguments:

left = x[..., :n]
right = x[..., n:]
g.log_density(x) == g.condition(right).log_density(left)

marginalize(left=0, right=0)[source]

Marginalizing out variables on either side of the event dimension:

g.marginalize(left=n).event_logsumexp() = g.logsumexp()
g.marginalize(right=n).event_logsumexp() = g.logsumexp()


and for data x:

g.condition(x).event_logsumexp()
= g.marginalize(left=g.dim() - x.size(-1)).log_density(x)
event_logsumexp()[source]

Integrates out all latent state (i.e. operating on event dimensions).

class AffineNormal(matrix, loc, scale)[source]

Bases: object

Represents a conditional diagonal normal distribution over a random variable Y whose mean is an affine function of a random variable X. The likelihood of X is thus:

AffineNormal(matrix, loc, scale).condition(y).log_density(x)


which is equivalent to:

Normal(x @ matrix + loc, scale).to_event(1).log_prob(y)

Parameters: matrix (torch.Tensor) – A transformation from X to Y. Should have rightmost shape (x_dim, y_dim). loc (torch.Tensor) – A constant offset for Y’s mean. Should have rightmost shape (y_dim,). scale (torch.Tensor) – Standard deviation for Y. Should have rightmost shape (y_dim,).
condition(value)[source]

Condition on a Y value.

Parameters: Return Gaussian: value (torch.Tensor) – A value of Y. A gaussian likelihood over X.
mvn_to_gaussian(mvn)[source]

Convert a MultivaiateNormal distribution to a Gaussian.

Parameters: mvn (MultivariateNormal) – A multivariate normal distribution. An equivalent Gaussian object. Gaussian
matrix_and_mvn_to_gaussian(matrix, mvn)[source]

Convert a noisy affine function to a Gaussian. The noisy affine function is defined as:

y = x @ matrix + mvn.sample()

Parameters: matrix (Tensor) – A matrix with rightmost shape (x_dim, y_dim). mvn (MultivariateNormal) – A multivariate normal distribution. A Gaussian with broadcasted batch shape and .dim() == x_dim + y_dim. Gaussian
gaussian_tensordot(x, y, dims=0)[source]

Computes the integral over two gaussians:

(x @ y)(a,c) = log(integral(exp(x(a,b) + y(b,c)), b)),

where x is a gaussian over variables (a,b), y is a gaussian over variables (b,c), (a,b,c) can each be sets of zero or more variables, and dims is the size of b.

Parameters: x – a Gaussian instance y – a Gaussian instance dims – number of variables to contract

## Statistical Utilities¶

gelman_rubin(input, chain_dim=0, sample_dim=1)[source]

Computes R-hat over chains of samples. It is required that input.size(sample_dim) >= 2 and input.size(chain_dim) >= 2.

Parameters: Returns torch.Tensor: input (torch.Tensor) – the input tensor. chain_dim (int) – the chain dimension. sample_dim (int) – the sample dimension. R-hat of input.
split_gelman_rubin(input, chain_dim=0, sample_dim=1)[source]

Computes R-hat over chains of samples. It is required that input.size(sample_dim) >= 4.

Parameters: Returns torch.Tensor: input (torch.Tensor) – the input tensor. chain_dim (int) – the chain dimension. sample_dim (int) – the sample dimension. split R-hat of input.
autocorrelation(input, dim=0)[source]

Computes the autocorrelation of samples at dimension dim.

Parameters: Returns torch.Tensor: input (torch.Tensor) – the input tensor. dim (int) – the dimension to calculate autocorrelation. autocorrelation of input.
autocovariance(input, dim=0)[source]

Computes the autocovariance of samples at dimension dim.

Parameters: Returns torch.Tensor: input (torch.Tensor) – the input tensor. dim (int) – the dimension to calculate autocorrelation. autocorrelation of input.
effective_sample_size(input, chain_dim=0, sample_dim=1)[source]

Computes effective sample size of input.

Reference:

 Introduction to Markov Chain Monte Carlo,
Charles J. Geyer
 Stan Reference Manual version 2.18,
Stan Development Team
Parameters: Returns torch.Tensor: input (torch.Tensor) – the input tensor. chain_dim (int) – the chain dimension. sample_dim (int) – the sample dimension. effective sample size of input.
resample(input, num_samples, dim=0, replacement=False)[source]

Draws num_samples samples from input at dimension dim.

Parameters: Returns torch.Tensor: input (torch.Tensor) – the input tensor. num_samples (int) – the number of samples to draw from input. dim (int) – dimension to draw from input. samples drawn randomly from input.
quantile(input, probs, dim=0)[source]

Computes quantiles of input at probs. If probs is a scalar, the output will be squeezed at dim.

Parameters: Returns torch.Tensor: input (torch.Tensor) – the input tensor. probs (list) – quantile positions. dim (int) – dimension to take quantiles from input. quantiles of input at probs.
pi(input, prob, dim=0)[source]

Computes percentile interval which assigns equal probability mass to each tail of the interval.

Parameters: Returns torch.Tensor: input (torch.Tensor) – the input tensor. prob (float) – the probability mass of samples within the interval. dim (int) – dimension to calculate percentile interval from input. quantiles of input at probs.
hpdi(input, prob, dim=0)[source]

Computes “highest posterior density interval” which is the narrowest interval with probability mass prob.

Parameters: Returns torch.Tensor: input (torch.Tensor) – the input tensor. prob (float) – the probability mass of samples within the interval. dim (int) – dimension to calculate percentile interval from input. quantiles of input at probs.
waic(input, log_weights=None, pointwise=False, dim=0)[source]

Computes “Widely Applicable/Watanabe-Akaike Information Criterion” (WAIC) and its corresponding effective number of parameters.

Reference:

 WAIC and cross-validation in Stan, Aki Vehtari, Andrew Gelman

Parameters: input (torch.Tensor) – the input tensor, which is log likelihood of a model. log_weights (torch.Tensor) – weights of samples along dim. dim (int) – the sample dimension of input. tuple of WAIC and effective number of parameters.
fit_generalized_pareto(X)[source]

Given a dataset X assumed to be drawn from the Generalized Pareto Distribution, estimate the distributional parameters k, sigma using a variant of the technique described in reference , as described in reference .

References  ‘A new and efficient estimation method for the generalized Pareto distribution.’ Zhang, J. and Stephens, M.A. (2009).  ‘Pareto Smoothed Importance Sampling.’ Aki Vehtari, Andrew Gelman, Jonah Gabry

Parameters: torch.Tensor – the input data X tuple of floats (k, sigma) corresponding to the fit parameters