杂项操作

pyro.ops 模块实现了大部分独立于 Pyro 其他部分的张量工具。

HMC的实用工具

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

基础类:object

双重平均法是一种解决凸优化问题的方案。它属于一类使用次梯度来更新模型参数(在原始空间中)的次梯度方法。在某些条件下,该方案中生成的参数的平均值保证会收敛到一个最优值。然而,传统次梯度方法的一个反直觉方面是“新的次梯度以递减的权重进入模型”(参见\([1]\))。双重平均法通过使用相等的权重更新次梯度(位于对偶空间中)的参数来解决这一现象,因此我们称之为“双重平均”。

该类实现了一种适用于马尔可夫链蒙特卡罗(MCMC)算法的双重平均方案。更准确地说,我们将用MCMC轨迹中计算的一些统计量来替代次梯度。此外,引入一些自由参数如t0kappa是有帮助的,并且仍然保证方案的收敛性。

参考文献

[1] 凸问题的原始对偶次梯度方法, 尤里·涅斯捷罗夫

[2] 无掉头采样器:在哈密顿蒙特卡洛中自适应设置路径长度, Matthew D. Hoffman, Andrew Gelman

Parameters
  • prox_center (float) – 在\([1]\)中引入的“prox-center”参数,它将原始序列拉向它。

  • t0 (float) – 在\([2]\)中引入的一个自由参数,用于稳定方案的初始步骤。

  • kappa (float) – 在\([2]\)中引入的一个自由参数,用于控制方案步骤的权重。 对于较小的kappa,方案将快速忘记早期步骤的状态。这应该是一个在\((0.5, 1]\)范围内的数字。

  • gamma (float) – 一个自由参数,用于控制方案的收敛速度。

reset()[source]
step(g)[source]

根据新的统计量/次梯度 g 更新方案的状态。

Parameters

g (float) – 在MCMC轨迹或次梯度中计算的统计量。

get_state()[source]

返回原始空间中的最新\(x_t\)\(\left\{x_i\right\}_{i=1}^t\)的平均值。

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

使用速度Verlet算法的二阶辛积分器。

Parameters
  • z (dict) – 样本站点名称及其当前值的字典 (类型 Tensor).

  • r (dict) – 样本站点名称和对应动量的字典 (类型 Tensor).

  • potential_fn (callable) – 返回给定z的势能的函数 对于每个采样点。函数的负梯度相对于 z 决定了相应站点的变化率 动量 r

  • kinetic_grad (callable) – 一个计算动能相对于动量变量的梯度的函数。

  • step_size (float) – 每次时间步迭代的步长。

  • num_steps (int) – 进行积分的离散时间步数。

  • z_grads (torch.Tensor) – 当前 z 处的势能梯度(可选)。

Return tuple (z_next, r_next, z_grads, potential_energy)

下一个位置和动量, 以及势能及其相对于 z_next 的梯度。

potential_grad(potential_fn, z)[源代码]

potential_fn 关于参数 z 的梯度。

Parameters
  • potential_fn – 一个可调用的python函数,它接收一个参数字典并返回势能。

  • z (dict) – 以站点名称为键的参数值字典。

Returns

元组 (z_grads, potential_energy),其中 z_grads 是一个字典, 其键与 z 相同,包含梯度,而 potential_energy 是一个 torch 标量。

register_exception_handler(name: str, handler: Callable[[Exception], bool], warn_on_overwrite: bool = True) None[source]

注册一个异常处理器,用于在评估势函数时处理(主要是数值)错误。

Parameters
  • name – 处理程序的名称(必须唯一)。

  • handler – 一个可调用的映射,将异常映射到布尔值。在任何handler中评估为true的异常将在势能的计算中被处理。

  • warn_on_overwrite – 如果为True,当覆盖已注册的处理程序时会发出警告。

class WelfordCovariance(diagonal=True)[source]

基础类:object

实现了Welford的在线方案用于估计(协)方差(参见\([1]\))。 对于适应HMC的对角线和密集质量结构非常有用。

参考文献

[1] 计算机程序设计艺术, Donald E. Knuth

reset()[source]
update(sample)[source]
get_covariance(regularize=True)[source]
class WelfordArrowheadCovariance(head_size=0)[source]

基础类:object

类似于 WelfordCovariance,但推广到了箭头结构。

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

获取箭头形式的协方差:(top, bottom_diag),其中 top = cov[:head_size]bottom_diag = cov.diag()[head_size:]

牛顿优化器

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

执行牛顿更新步骤以最小化一批变量的损失,可选择约束到信任区域 [1]。

这尤其有用,因为牛顿迭代的最终解相对于输入是可微的,即使除了最终的x之外的所有内容都被分离,由于该方法的二次收敛性[2]。loss必须作为x的函数是二次可微的。如果loss2+d次可微的,那么该函数的返回值是d次可微的。

loss被解释为负对数概率密度时,可以使用此函数的返回值mode,cov来构建拉普拉斯近似MultivariateNormal(mode,cov)

警告

在优化循环中使用此函数时,请注意分离其结果。如果在优化过程中忘记分离此函数的结果,那么反向传播将贯穿整个迭代过程,更糟糕的是,每一步都会计算两个额外的导数。

在循环中的示例用法:

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
[1] Yuan, Ya-xiang. Iciam. Vol. 99. 2000.

“优化中信任区域算法的回顾。” ftp://ftp.cc.ac.cn/pub/yyx/papers/p995.pdf

[2] Christianson, Bruce. Optimization Methods and Software 3.4 (1994)

“反向累积和吸引固定点。” http://uhra.herts.ac.uk/bitstream/handle/2299/4338/903839.pdf

Parameters
  • loss (torch.Tensor) – 一个关于 x 的标量函数,需要被最小化。

  • x (torch.Tensor) – 一个形状为 (N, D) 的因变量,其中 N 是批量大小,D 是一个较小的数字。

  • trust_radius (float) – 一个可选的信任区域 trust_radius。该函数的更新值 mode 将在输入 xtrust_radius 范围内。

Returns

一对 (mode, cov),其中 mode 是一个与原始值 x 形状相同的更新张量,而 cov 是一个估计的协方差 DxD 矩阵,其形状为 cov.shape == x.shape[:-1] + (D,D)

Return type

tuple

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

执行牛顿更新步骤以最小化一批一维变量的损失,可选择正则化以约束到信任区域。

详情请参见 newton_step()

Parameters
  • loss (torch.Tensor) – 一个关于 x 的标量函数,需要被最小化。

  • x (torch.Tensor) – 一个最右侧大小为1的因变量。

  • trust_radius (float) – 一个可选的信任区域 trust_radius。此函数的更新值 mode 将在输入 xtrust_radius 范围内。

Returns

一对 (mode, cov) 其中 mode 是一个与原始值 x 形状相同的更新张量,而 cov 是一个估计的协方差 1x1 矩阵,其形状为 cov.shape == x.shape[:-1] + (1,1)

Return type

tuple

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

执行牛顿更新步骤以最小化一批二维变量的损失,可选择正则化以约束到信任区域。

详情请参见 newton_step()

Parameters
  • loss (torch.Tensor) – 一个关于 x 的标量函数,需要被最小化。

  • x (torch.Tensor) – 一个最右侧大小为2的因变量。

  • trust_radius (float) – 一个可选的信任区域 trust_radius。此函数的更新值 mode 将在输入 xtrust_radius 范围内。

Returns

一对 (mode, cov) 其中 mode 是一个与原始值 x 形状相同的更新张量,而 cov 是一个估计的协方差 2x2 矩阵,其形状为 cov.shape == x.shape[:-1] + (2,2)

Return type

tuple

newton_step_3d(loss, x, trust_radius=None)[源代码]

执行牛顿更新步骤以最小化一批三维变量的损失,可选择正则化以约束到信任区域。

详情请参见 newton_step()

Parameters
  • loss (torch.Tensor) – 一个关于 x 的标量函数,需要被最小化。

  • x (torch.Tensor) – 一个最右侧大小为2的因变量。

  • trust_radius (float) – 一个可选的信任区域 trust_radius。此函数的更新值 mode 将在输入 xtrust_radius 范围内。

Returns

一对 (mode, cov) 其中 mode 是一个与原始值 x 形状相同的更新张量,而 cov 是一个估计的协方差 3x3 矩阵,其形状为 cov.shape == x.shape[:-1] + (3,3)

Return type

tuple

特殊函数

safe_log(x)[source]

类似于 torch.log(),但通过将 log(0) 的梯度限制在最多 1 / finfo.eps 来避免无限梯度。

log_beta(x, y, tol=0.0)[source]

计算对数Beta函数。

tol >= 0.02 时,这使用了一个移位的斯特林近似来计算对数贝塔函数。该近似方法调整了斯特林近似来计算对数伽玛函数:

lgamma(z) ≈ (z - 1/2) * log(z) - z + log(2 * pi) / 2

近似对数Beta函数:

log_beta(x, y) ≈ ((x-1/2) * log(x) + (y-1/2) * log(y)
                  - (x+y-1/2) * log(x+y) + log(2*pi)/2)

近似值通过使用递归迭代地移动对数伽马近似值,从而在接近零时进一步提高准确性:

lgamma(x) = lgamma(x + 1) - log(x)

如果此递归应用了n次,则绝对误差被限制在 error < 0.082 / n < tol,因此我们根据用户提供的tol选择n

Parameters
  • x (torch.Tensor) – 一个正张量。

  • y (torch.Tensor) – 一个正张量。

  • tol (float) – 最大绝对误差的界限。默认为0.1。对于非常小的tol,此函数会直接调用log_beta()

Return type

torch.Tensor

log_binomial(n, k, tol=0.0)[source]

计算对数二项式系数。

tol >= 0.02 时,这使用了一个移位的斯特林近似来计算对数 Beta 函数,通过 log_beta()

Parameters
Return type

torch.Tensor

log_I1(orders: int, value: torch.Tensor, terms=250)[source]

计算第一类修正贝塞尔函数的前n个对数 .. math

\log(I_v(z)) = v*\log(z/2) + \log(\sum_{k=0}^\inf \exp\left[2*k*\log(z/2) - \sum_kk^k log(kk)
- \lgamma(v + k + 1)\right])
Parameters
  • orders – 对数修正贝塞尔函数的阶数。

  • value – 用于计算修正贝塞尔函数的值

  • terms – 求和的截断

Returns

0 到 orders 修改的贝塞尔函数

get_quad_rule(num_quad, prototype_tensor)[source]

获取具有指定数量积分点的高斯-埃尔米特积分规则的积分点及相应的对数权重。

示例用法:

quad_points, log_weights = get_quad_rule(32, prototype_tensor)
# transform to N(0, 4.0) Normal distribution
quad_points *= 4.0
# compute variance integral in log-space using logsumexp and exponentiate
variance = torch.logsumexp(quad_points.pow(2.0).log() + log_weights, axis=0).exp()
assert (variance - 16.0).abs().item() < 1.0e-6
Parameters
  • num_quad (int) – 积分点的数量。

  • prototype_tensor (torch.Tensor) – 用于确定返回张量的dtypedevice

Returns

元组形式为torch.Tensor`s的`(quad_points, log_weights)

sparse_multinomial_likelihood(total_count, nonzero_logits, nonzero_value)[source]

以下是等价的:

# Version 1. dense
log_prob = Multinomial(logits=logits).log_prob(value).sum()

# Version 2. sparse
nnz = value.nonzero(as_tuple=True)
log_prob = sparse_multinomial_likelihood(
    value.sum(-1),
    (logits - logits.logsumexp(-1))[nnz],
    value[nnz],
)

张量工具

as_complex(x)[source]

类似于 torch.view_as_complex(),但在步幅不是两的倍数的情况下会复制数据。

block_diag_embed(mat)[source]

接受一个形状为 (…, B, M, N) 的张量,并返回一个形状为 (…, B x M, B x N) 的块对角张量。

Parameters

mat (torch.Tensor) – 一个具有3个或更多维度的输入张量

Returns torch.Tensor

一个具有维度 m.dim() - 1 的块对角张量

block_diagonal(mat, block_size)[源代码]

接受一个形状为 (…, B x M, B x N) 的块对角张量,并返回一个形状为 (…, B, M, N) 的张量。

Parameters
  • mat (torch.Tensor) – 一个具有2个或更多维度的输入张量

  • block_size (int) – 块的数量 B。

Returns torch.Tensor

一个维度为 mat.dim() + 1 的张量

periodic_repeat(tensor, size, dim)[source]

重复一个period大小的张量直到达到给定的size。例如:

>>> x = torch.tensor([[1, 2, 3], [4, 5, 6]])
>>> periodic_repeat(x, 4, 0)
tensor([[1, 2, 3],
        [4, 5, 6],
        [1, 2, 3],
        [4, 5, 6]])
>>> periodic_repeat(x, 4, 1)
tensor([[1, 2, 3, 1],
        [4, 5, 6, 4]])

这对于计算时间序列模型中的静态季节性非常有用。

Parameters
  • 张量 (torch.Tensor) – 一个差异的张量。

  • size (int) – 沿着维度 dim 的结果的期望大小。

  • dim (int) – 沿着该张量维度进行重复。

periodic_cumsum(tensor, period, dim)[source]

沿给定维度计算周期性累积和。例如,如果 dim=0:

for t in range(period):
    assert result[t] == tensor[t]
for t in range(period, len(tensor)):
    assert result[t] == tensor[t] + result[t - period]

这对于计算时间序列模型中的漂移季节性非常有用。

Parameters
  • tensor (torch.Tensor) – 一个差异的张量。

  • period (int) – 重复的周期。

  • dim (int) – 沿此张量维度进行累积。

periodic_features(duration, max_period=None, min_period=None, **options)[source]

max_periodmin_period创建周期性(sin,cos)特征。

这在时间序列模型中非常有用,其中长期不均匀的季节性可以通过回归来处理。当仅指定max_period时,这会在所有长度尺度上生成周期性特征。当还指定min_period时,这会在较大的长度尺度上生成周期性特征,但会忽略高频特征。这在将长期季节性的回归与其他技术(如periodic_repeat()periodic_cumsum())结合用于短时间尺度时非常有用。例如,要将年度季节性回归结合到一周的尺度,可以设置max_period=365.25min_period=7

Parameters
  • duration (int) – 离散时间步数。

  • max_period (float) – 可选的最大周期,默认为 duration

  • min_period (float) – 可选的最小周期(不包括),默认为 2 = 奈奎斯特截止频率。

  • **options – 张量构造选项,例如 dtypedevice

Returns

一个形状为(duration, 2 * ceil(max_period / min_period) - 2)的张量,其特征已归一化到[-1,1]范围内。

Return type

Tensor

next_fast_len(size)[source]

返回下一个最大的数字 n >= size,其质因数全部为2、3或5。这些大小对于快速傅里叶变换是高效的。等同于 scipy.fftpack.next_fast_len()

Parameters

size (int) – 一个正数。

Returns

可能更大的数字。

Rtype int

convolve(signal, kernel, mode='full')[source]

使用FFT计算信号与核的一维卷积。 两个参数应具有相同的最右侧维度,但除此之外可以是任意可广播的。

Parameters
  • signal (torch.Tensor) – 一个要进行卷积的信号。

  • kernel (torch.Tensor) – 一个卷积核。

  • mode (str) – 其中之一:'full', 'valid', 'same'。

Returns

一个具有广播形状的张量。假设 m = signal.size(-1)n = kernel.size(-1),结果的最右边大小将是: m + n - 1 如果模式是‘full’; max(m, n) - min(m, n) + 1 如果模式是‘valid’;或者 max(m, n) 如果模式是‘same’。

Rtype torch.Tensor

repeated_matmul(M, n)[source]

将一批矩阵 M 作为输入,并返回执行 n 次矩阵乘法 \(M\), \(M^2\), …, \(M^n\) 的堆叠结果。 并行成本在 n 中是对数的。

Parameters
  • M (torch.Tensor) – 一批形状为 (…, N, N) 的方形张量。

  • n (int) – 最大乘积的阶数 \(M^n\)

Returns torch.Tensor

一批形状为 (n, …, N, N) 的方形张量

dct(x, dim=- 1)[source]

类型II的离散余弦变换,缩放为正交归一化。

这是idct_ii()的逆操作,并且等同于 scipy.fftpack.dct() 使用 norm="ortho"

Parameters
  • x (Tensor) – 输入信号。

  • dim (int) – 计算DCT的维度。

Return type

张量

idct(x, dim=- 1)[source]

类型II的逆离散余弦变换,缩放为正交归一化。

这是dct_ii()的逆函数,并且等同于 scipy.fftpack.idct() 使用 norm="ortho"

Parameters
  • x (Tensor) – 输入信号。

  • dim (int) – 计算DCT的维度。

Return type

张量

haar_transform(x)[source]

离散哈尔变换。

沿最后一个维度执行Haar变换。 这是inverse_haar_transform()的逆操作。

Parameters

x (Tensor) – 输入信号。

Return type

张量

inverse_haar_transform(x)[source]

沿最后一个维度执行逆哈尔变换。 这是haar_transform()的逆操作。

Parameters

x (Tensor) – 输入信号。

Return type

张量

safe_cholesky(x)[source]
cholesky_solve(x, y)[source]
matmul(x, y)[source]
matvecmul(x, y)[source]
triangular_solve(x, y, upper=False, transpose=False)[源代码]
precision_to_scale_tril(P)[source]
safe_normalize(x, *, p=2)[source]

安全地将向量投影到球体上,相对于p-范数。这通过将零映射到向量[1, 0, 0, ..., 0]来避免在零处的奇点。

Parameters
  • x (torch.Tensor) – 一个向量

  • p (float) – 范数指数,默认为2,即欧几里得范数。

Returns

归一化版本 x / ||x||_p

Return type

张量

张量索引

index(tensor, args)[source]

使用嵌套元组进行索引。

另请参阅便捷包装器 Index

这对于编写与多种解释兼容的索引代码非常有用,例如标量评估、向量化评估或重塑。

例如,假设 x 是一个参数,且 x.dim() == 2,我们希望 推广表达式 x[..., t],其中 t 可以是以下任意一种:

  • 一个标量 t=1x[..., 1];

  • 一个切片 t=slice(None) 等同于 x[..., :];或者

  • 一个重塑操作 t=(Ellipsis, None) 等同于 x.unsqueeze(-1)

虽然简单的索引适用于前两个例子,但第三个例子会导致一个嵌套的元组 (Ellipsis, (Ellipsis, None))。这个辅助函数会将该嵌套元组展平并合并连续的 Ellipsis

Parameters
  • 张量 (torch.Tensor) – 一个需要被索引的张量。

  • args (tuple) – 一个索引,作为 __getitem__ 的参数。

Returns

tensor[args] 的扁平化解释。

Return type

torch.Tensor

class Index(tensor)[source]

基础类:object

围绕 index() 的便捷包装器。

以下是等价的:

Index(x)[..., i, j, :]
index(x, (Ellipsis, i, j, slice(None)))
Parameters

tensor (torch.Tensor) – 一个要被索引的张量。

Returns

一个具有特殊__getitem__()方法的对象。

vindex(tensor, args)[source]

具有广播语义的矢量化高级索引。

另请参阅便捷包装器 Vindex

这对于编写与批处理和枚举兼容的索引代码非常有用,特别是对于使用离散随机变量选择混合组件的情况。

例如,假设 x 是一个参数,且 x.dim() == 3,我们希望将表达式 x[i, :, j] 从整数 i,j 推广到具有批次维度和枚举维度(但没有事件维度)的张量 i,j。然后我们可以使用 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

处理当x可能还包含批次维度的情况时(例如,如果x在使用向量化粒子时在板状上下文中采样),vindex()使用特殊约定,即Ellipsis表示批次维度(因此...只能出现在左侧,不能出现在中间或右侧)。假设x具有事件维度3。那么我们可以写成:

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

请注意,这种对Ellipsis的特殊处理与NEP [1]不同。

正式地,此函数假设:

  1. 每个参数要么是Ellipsis,要么是slice(None),一个整数,或者一个批处理的torch.LongTensor(即具有空事件形状)。此函数不支持非平凡的切片或torch.BoolTensor掩码。Ellipsis只能作为args[0]出现在左侧。

  2. 如果 args[0] is not Ellipsis 那么 tensor 不是 批处理的,并且其事件维度等于 len(args)

  3. 如果 args[0] is Ellipsis 那么 tensor 是批处理的,并且 其事件维度等于 len(args[1:])tensor 的事件维度左侧的维度 被视为批处理维度,并将与张量参数的维度进行广播。

请注意,如果没有任何参数是具有.dim() > 0的张量,那么此函数的行为类似于标准索引:

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

参考文献

[1] https://www.numpy.org/neps/nep-0021-advanced-indexing.html

介绍了vindex作为向量化索引的辅助工具。 Pyro的实现与提议的符号x.vindex[]类似,除了对Ellipsis的处理略有不同。

Parameters
  • tensor (torch.Tensor) – 一个要被索引的张量。

  • args (tuple) – 一个索引,作为 __getitem__ 的参数。

Returns

tensor[args]的非标准解释。

Return type

torch.Tensor

class Vindex(tensor)[source]

基础类:object

围绕vindex()的便捷包装器。

以下是等价的:

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

tensor (torch.Tensor) – 一个要被索引的张量。

Returns

一个具有特殊__getitem__()方法的对象。

张量收缩

contract(equation, *operands, **kwargs)[source]

围绕opt_einsum.contract()的包装器,可选择使用Pyro的廉价优化器,并可选择缓存收缩路径。

Parameters

cache_path (bool) – 是否缓存收缩路径。 默认为 True。

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

围绕opt_einsum.contract_expression()的包装器,可选择使用Pyro的廉价优化器,并可选择缓存收缩路径。

Parameters

cache_path (bool) – 是否缓存收缩路径。 默认为 True。

einsum(equation, *operands, **kwargs)[source]

通过张量变量消除的广义板求和积算法。

这以两种方式概括了contract()

  1. 允许多个输出,并且可以共享中间结果。

  2. 输入和输出可以沿着plates中给出的符号进行排列;沿着plates的减少是乘积减少。

理解此函数的最佳方法是尝试下面的示例,这些示例展示了如何将einsum()调用实现为对contract()的多次调用(这通常更昂贵)。

为了说明多个输出,请注意以下是等价的:

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)

为了说明plated输入,请注意以下是等价的:

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)

当求和维度 a 总是与板维度 i 一起出现时,a 对应于 a 的每个切片的独特符号。因此,以下是等价的:

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)

当这样的求和维度出现在输出中时,它必须伴随其所有的板维度,例如以下是等价的:

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

z0 = contract('ab,ac,ad,ab,ac,ad->b', *x, *y)
z1 = contract('ab,ac,ad,ab,ac,ad->c', *x, *y)
z2 = contract('ab,ac,ad,ab,ac,ad->d', *x, *y)
z = torch.stack([z0, z1, z2])

请注意,输出中的每个板块切片在所有输入的所有板块切片中都是多线性的,因此例如批量矩阵乘法将不会使用plates来实现,因此以下内容都是等效的:

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)

在所有有效的方程中,有些计算是输入张量大小的多项式,而其他计算是输入张量大小的指数。当计算是指数时,此函数会引发NotImplementedError

Parameters
  • equation (str) – 一个einsum方程,可选地包含多个输出。

  • 操作数 (torch.Tensor) – 一组张量。

  • plates (str) – 一个可选的板符号字符串。

  • backend (str) – 一个可选的einsum后端,默认为‘torch’。

  • cache (dict) – 一个可选的 shared_intermediates() 缓存。

  • modulo_total (bool) – 可选地允许einsum任意缩放每个结果板,这可以显著减少计算量。只要每个结果板表示一个非归一化的概率分布,且其总和不感兴趣,就可以安全地设置此选项。

Returns

一个请求形状的张量元组,每个输出一个条目。

Return type

tuple

Raises
  • ValueError – 如果张量大小不匹配或输出请求了一个没有该维度板的维度。

  • NotImplementedError – 如果收缩在任何输入张量的大小上会导致指数级成本。

ubersum(equation, *operands, **kwargs)[source]

已弃用,请改用 einsum()

高斯收缩

class Gaussian(log_normalizer: torch.Tensor, info_vec: torch.Tensor, precision: torch.Tensor)[source]

基础类:object

非归一化的高斯分布。

这表示一个任意的半定二次函数,可以解释为一个秩不足的缩放高斯分布。精度矩阵可能具有零特征值,因此可能无法直接使用协方差矩阵。

Parameters
  • log_normalizer (torch.Tensor) – 一个归一化常数,主要用于在收缩过程中跟踪归一化项。

  • info_vec (torch.Tensor) – 信息向量,它是均值的缩放版本 info_vec = precision @ mean。我们使用这种表示法来使高斯收缩快速且稳定。

  • precision (torch.Tensor) – 该高斯分布的精度矩阵。

dim()[source]
property batch_shape
expand(batch_shape) pyro.ops.gaussian.Gaussian[source]
reshape(batch_shape) pyro.ops.gaussian.Gaussian[source]
__getitem__(index) pyro.ops.gaussian.Gaussian[来源]

索引到高斯分布的批次形状中。

static cat(parts, dim=0) pyro.ops.gaussian.Gaussian[source]

沿给定的批次维度连接高斯分布列表。

event_pad(left=0, right=0) pyro.ops.gaussian.Gaussian[source]

沿事件维度填充。

event_permute(perm) pyro.ops.gaussian.Gaussian[source]

沿着事件维度进行置换。

__add__(other: Union[pyro.ops.gaussian.Gaussian, int, float, torch.Tensor]) pyro.ops.gaussian.Gaussian[source]

在log-density空间中添加两个高斯分布。

log_density(value: torch.Tensor) torch.Tensor[source]

评估此高斯分布在某点的对数密度:

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

这主要用于测试。

rsample(sample_shape=torch.Size([]), noise: Optional[torch.Tensor] = None) torch.Tensor[source]

重新参数化的采样器。

condition(value: torch.Tensor) pyro.ops.gaussian.Gaussian[source]

在这个高斯分布的尾部子集上进行条件化。 这应该满足:

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

请注意,由于这是一个非归一化的高斯分布,我们在结果中包含了y的密度。因此,condition()类似于functools.partial的参数绑定:

left = x[..., :n]
right = x[..., n:]
g.log_density(x) == g.condition(right).log_density(left)
left_condition(value: torch.Tensor) pyro.ops.gaussian.Gaussian[source]

将这个高斯条件设置在其状态的前导子集上。 这应该满足:

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

请注意,由于这是一个非归一化的高斯分布,我们在结果中包含了y的密度。因此,condition()类似于functools.partial的参数绑定:

left = x[..., :n]
right = x[..., n:]
g.log_density(x) == g.left_condition(left).log_density(right)
marginalize(left=0, right=0) pyro.ops.gaussian.Gaussian[source]

在事件维度的任一侧边缘化变量:

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

对于数据 x

g.condition(x).event_logsumexp()

= g.marginalize(left=g.dim() - x.size(-1)).log_density(x)

event_logsumexp() torch.Tensor[来源]

整合所有潜在状态(即在事件维度上操作)。

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

基础类:object

表示一个随机变量 Y 的条件对角正态分布,其均值是随机变量 X 的仿射函数。 因此,X 的似然函数为:

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

这相当于:

Normal(x @ matrix + loc, scale).to_event(1).log_prob(y)
Parameters
  • matrix (torch.Tensor) – 从 XY 的转换。 应具有最右侧的形状 (x_dim, y_dim)

  • loc (torch.Tensor) – Y 的均值的常数偏移量。 应具有最右侧形状 (y_dim,)

  • scale (torch.Tensor) – Y 的标准差。 应具有最右侧形状 (y_dim,)

property batch_shape
condition(value)[source]
left_condition(value)[source]

如果 value.size(-1) == x_dim,这将返回一个具有 event_dim=1 的正态分布。应用此方法后,绘制样本的成本是 O(y_dim) 而不是 O(y_dim ** 3)

rsample(sample_shape=torch.Size([]), noise: Optional[torch.Tensor] = None) torch.Tensor[source]

重新参数化的采样器。

to_gaussian()[来源]
expand(batch_shape)[source]
reshape(batch_shape)[source]
__getitem__(index)[source]
event_permute(perm)[source]
__add__(other)[来源]
marginalize(left=0, right=0)[source]
mvn_to_gaussian(mvn)[source]

将多元正态分布转换为高斯分布。

Parameters

mvn (MultivariateNormal) – 一个多元正态分布。

Returns

一个等价的高斯对象。

Return type

Gaussian

matrix_and_gaussian_to_gaussian(matrix: torch.Tensor, y_gaussian: pyro.ops.gaussian.Gaussian) pyro.ops.gaussian.Gaussian[source]

构建一个条件高斯分布 p(y|x),其中 y - x @ matrix ~ y_gaussian

Parameters
  • matrix (torch.Tensor) – 一个右作用的变换矩阵。

  • y_gaussian (Gaussian) – 一个关于噪声的分布,噪声为 y - x@matrix

Return type

Gaussian

matrix_and_mvn_to_gaussian(matrix, mvn)[source]

将噪声仿射函数转换为高斯函数。噪声仿射函数定义为:

y = x @ matrix + mvn.sample()
Parameters
  • matrix (Tensor) – 一个具有最右侧形状 (x_dim, y_dim) 的矩阵。

  • mvn (MultivariateNormal) – 一个多元正态分布。

Returns

一个具有广播批量形状的高斯分布,且.dim() == x_dim + y_dim

Return type

Gaussian

gaussian_tensordot(x: pyro.ops.gaussian.Gaussian, y: pyro.ops.gaussian.Gaussian, dims: int = 0) pyro.ops.gaussian.Gaussian[source]

计算两个高斯函数的积分:

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

其中 x 是变量 (a,b) 上的高斯分布,y 是变量 (b,c) 上的高斯分布,(a,b,c) 可以是零个或多个变量的集合,dims 是 b 的大小。

Parameters
  • x – 一个高斯实例

  • y – 一个高斯实例

  • dims – 要收缩的变量数量

sequential_gaussian_tensordot(gaussian: pyro.ops.gaussian.Gaussian) pyro.ops.gaussian.Gaussian[源代码]

对高斯x进行积分,其最右边的批次维度是时间,计算如下:

x[..., 0] @ x[..., 1] @ ... @ x[..., T-1]
Parameters

gaussian (Gaussian) – 一个批处理的高斯分布,其最右边的维度是时间。

Returns

高斯分布沿其时间维度的马尔可夫乘积。

Return type

Gaussian

sequential_gaussian_filter_sample(init: pyro.ops.gaussian.Gaussian, trans: pyro.ops.gaussian.Gaussian, sample_shape: Tuple[int, ...] = (), noise: Optional[torch.Tensor] = None) torch.Tensor[source]

通过并行扫描前向滤波后向采样,从高斯的马尔可夫积中绘制一个重新参数化的样本。

Parameters
  • init (Gaussian) – 表示初始状态的高斯分布。

  • trans (Gaussian) – 一个表示一系列状态转换的高斯分布,时间作为最右边的批次维度。这个高斯分布的事件维度必须是init的两倍:trans.dim() == 2 * init.dim()

  • sample_shape (tuple) – 一个可选的额外形状,用于绘制样本。

  • 噪声 (torch.Tensor) – 一个可选的标准白噪声张量,形状为 sample_shape + batch_shape + (duration, state_dim),其中 duration = 1 + trans.batch_shape[-1] 是要采样的时间点数, state_dim = init.dim() 是状态维度。 这对于计算均值(传递零)、变化温度(传递缩放噪声)和对立采样(传递 cat([z,-z]))非常有用。

Returns

一个重新参数化的样本,形状为 sample_shape + batch_shape + (duration, state_dim)

Return type

torch.Tensor

统计工具

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

计算样本链的R-hat。要求 input.size(sample_dim) >= 2input.size(chain_dim) >= 2

Parameters
  • 输入 (torch.Tensor) – 输入张量。

  • chain_dim (int) – 链的维度。

  • sample_dim (int) – 样本维度。

Returns torch.Tensor

input 的 R-hat。

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

计算样本链上的R-hat。要求input.size(sample_dim) >= 4

Parameters
  • 输入 (torch.Tensor) – 输入张量。

  • chain_dim (int) – 链的维度。

  • sample_dim (int) – 样本维度。

Returns torch.Tensor

分割 input 的 R-hat。

autocorrelation(input, dim=0)[source]

计算样本在维度 dim 上的自相关。

参考:https://en.wikipedia.org/wiki/Autocorrelation#Efficient_computation

Parameters
  • 输入 (torch.Tensor) – 输入张量。

  • dim (int) – 计算自相关的维度。

Returns torch.Tensor

input的自相关。

autocovariance(input, dim=0)[source]

计算样本在维度 dim 上的自协方差。

Parameters
  • 输入 (torch.Tensor) – 输入张量。

  • dim (int) – 计算自相关的维度。

Returns torch.Tensor

input的自相关。

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

计算输入的有效样本大小。

参考:

[1] Introduction to Markov Chain Monte Carlo,

查尔斯·J·盖尔

[2] Stan Reference Manual version 2.18,

斯坦开发团队

Parameters
  • 输入 (torch.Tensor) – 输入张量。

  • chain_dim (int) – 链的维度。

  • sample_dim (int) – 样本维度。

Returns torch.Tensor

input 的有效样本量。

resample(input, num_samples, dim=0, replacement=False)[source]

input的维度dim中抽取num_samples个样本。

Parameters
  • 输入 (torch.Tensor) – 输入张量。

  • num_samples (int) – 从 input 中抽取的样本数量。

  • dim (int) – 从 input 中提取的维度。

Returns torch.Tensor

input中随机抽取的样本。

quantile(input, probs, dim=0)[source]

计算inputprobs处的分位数。如果probs是标量,输出将在dim处被压缩。

Parameters
  • 输入 (torch.Tensor) – 输入张量。

  • probs (list) – 分位数位置。

  • dim (int) – 从 input 中获取分位数的维度。

Returns torch.Tensor

inputprobs 的分位数。

weighed_quantile(input: torch.Tensor, probs: Union[List[float], Tuple[float, ...], torch.Tensor], log_weights: torch.Tensor, dim: int = 0) torch.Tensor[源代码]

计算在probs处的加权input样本的分位数。

Parameters
  • 输入 (torch.Tensor) – 输入张量。

  • probs (list) – 分位数位置。

  • log_weights (torch.Tensor) – 样本权重张量。

  • dim (int) – 从 input 中获取分位数的维度。

Returns torch.Tensor

inputprobs 的分位数。

示例:

>>> from pyro.ops.stats import weighed_quantile
>>> import torch
>>> input = torch.Tensor([[10, 50, 40], [20, 30, 0]])
>>> probs = torch.Tensor([0.2, 0.8])
>>> log_weights = torch.Tensor([0.4, 0.5, 0.1]).log()
>>> result = weighed_quantile(input, probs, log_weights, -1)
>>> torch.testing.assert_close(result, torch.Tensor([[40.4, 47.6], [9.0, 26.4]]))
pi(input, prob, dim=0)[source]

计算百分位数区间,该区间为区间的每个尾部分配相等的概率质量。

Parameters
  • 输入 (torch.Tensor) – 输入张量。

  • prob (float) – 区间内样本的概率质量。

  • dim (int) – 从 input 计算百分位区间的维度。

Returns torch.Tensor

inputprobs 的分位数。

hpdi(input, prob, dim=0)[source]

计算“最高后验密度区间”,这是概率质量prob的最窄区间。

Parameters
  • 输入 (torch.Tensor) – 输入张量。

  • prob (float) – 区间内样本的概率质量。

  • dim (int) – 从 input 计算百分位区间的维度。

Returns torch.Tensor

inputprobs 的分位数。

waic(input, log_weights=None, pointwise=False, dim=0)[source]

计算“广泛适用/渡边-赤池信息准则”(WAIC)及其对应的有效参数数量。

参考:

[1] WAIC 和 Stan 中的交叉验证, Aki Vehtari, Andrew Gelman

Parameters
  • 输入 (torch.Tensor) – 输入张量,它是模型的对数似然。

  • log_weights (torch.Tensor) – 沿着 dim 维度的样本权重。

  • dim (int) – input的样本维度。

Returns tuple

WAIC 和有效参数数量的元组。

fit_generalized_pareto(X)[source]

给定一个假设从广义帕累托分布中抽取的数据集X,使用参考文献[1]中描述的技术变体来估计分布参数k和sigma,如参考文献[2]中所述。

参考文献 [1] ‘一种新的高效广义帕累托分布估计方法。’ 张杰和斯蒂芬斯,M.A. (2009). [2] ‘帕累托平滑重要性采样。’ 阿基·维塔里,安德鲁·格尔曼,乔纳·加布里

Parameters

torch.Tensor – 输入数据 X

Returns tuple

浮点数元组 (k, sigma) 对应于拟合参数

crps_empirical(pred, truth)[source]

计算一组样本 pred 和真实数据 truth 之间的负连续排名概率评分 CRPS* [1]。这使用了一个 n log(n) 时间算法来计算一个量,该量在样本数量 n 上原本具有二次复杂度:

CRPS* = E|pred - truth| - 1/2 E|pred - pred'|
      = (pred - truth).abs().mean(0)
      - (pred - pred.unsqueeze(1)).abs().mean([0, 1]) / 2

请注意,对于单个样本,这简化为绝对误差。

参考文献

[1] Tilmann Gneiting, Adrian E. Raftery (2007)

严格适当的评分规则、预测和估计 https://www.stat.washington.edu/raftery/Research/PDF/Gneiting2007jasa.pdf

Parameters
  • pred (torch.Tensor) – 一组在最右侧维度上批处理的样本预测。 这应该具有形状 (num_samples,) + truth.shape

  • 真实值 (torch.Tensor) – 一个包含真实观测值的张量。

Returns

一个形状为 truth.shape 的张量。

Return type

torch.Tensor

energy_score_empirical(pred: torch.Tensor, truth: torch.Tensor) torch.Tensor[source]

计算一组多元样本 pred 和真实数据向量 truth 之间的负能量评分 ES*(参见[1]中的公式22)。运行时间与样本数量 n 的平方成正比。在单变量样本的情况下,输出与 CRPS 一致:

ES* = E|pred - truth| - 1/2 E|pred - pred'|

请注意,对于单个样本,这简化为样本predtruth之间差异的欧几里得范数。

这是一个严格的适当分数,因此对于根据分布 \(P\) 分布的 pred 和根据分布 \(Q\) 分布的 truth,我们有 \(ES^{*}(P,Q) \ge ES^{*}(Q,Q)\),当且仅当 \(P=Q\) 时等号成立,即如果 \(P\)\(Q\) 具有相同的多元分布(仅当 \(P\)\(Q\) 具有相同的边缘分布时,等号成立是不够的)。

参考文献

[1] Tilmann Gneiting, Adrian E. Raftery (2007)

严格适当的评分规则、预测和估计 https://www.stat.washington.edu/raftery/Research/PDF/Gneiting2007jasa.pdf

Parameters
  • pred (torch.Tensor) – 一组样本预测,批处理在第二个最左边的维度上。 最左边的维度是多变量样本的维度。

  • 真实值 (torch.Tensor) – 一个与pred形状相同的真实观测值的张量,除了最左边的第二个维度可以有任何值或被省略。

Returns

一个形状为 truth.shape 的张量。

Return type

torch.Tensor

流式统计

class StreamingStats[source]

基础类:abc.ABC

用于张量树的可流式统计的抽象基类。

派生类必须实现 update(), merge(), 和 get().

abstract update(sample) None[source]

从单个样本更新状态。

这会改变 self 并且不返回任何内容。更新应该是独立于顺序的,即样本应该是可交换的。

Parameters

sample – 一个样本值,它是一个嵌套的字典,包含 torch.Tensor 叶子节点。它可以有任意的嵌套和 形状,但假设在调用 .update() 时形状保持不变。

abstract merge(other) pyro.ops.streaming.StreamingStats[source]

选择两个聚合统计量,例如来自不同的MCMC链。

这是一个纯函数:它返回一个新的StreamingStats对象,并且不会修改selfother

Parameters

other – 另一个相同类型的流统计实例。

abstract get() Any[source]

返回聚合统计量。

class StatsOfDict(types: Dict[Hashable, Callable[[], pyro.ops.streaming.StreamingStats]] = {}, default: Callable[[], pyro.ops.streaming.StreamingStats] = <class 'pyro.ops.streaming.CountStats'>)[source]

基础类:pyro.ops.streaming.StreamingStats

统计样本,这些样本是具有恒定键集的字典。

例如,以下是等价的:

# Version 1. Hand encode statistics.
>>> a_stats = CountStats()
>>> b_stats = CountMeanStats()
>>> a_stats.update(torch.tensor(0.))
>>> b_stats.update(torch.tensor([1., 2.]))
>>> summary = {"a": a_stats.get(), "b": b_stats.get()}

# Version 2. Collect samples into dictionaries.
>>> stats = StatsOfDict({"a": CountStats, "b": CountMeanStats})
>>> stats.update({"a": torch.tensor(0.), "b": torch.tensor([1., 2.])})
>>> summary = stats.get()
>>> summary
{'a': {'count': 1}, 'b': {'count': 1, 'mean': tensor([1., 2.])}}
Parameters
  • default – 字典值的统计默认类型。默认为廉价的 CountStats

  • types (dict) – 字典映射键到应记录的统计类型,该统计类型对应于该键的值。

update(sample: Dict[Hashable, Any]) None[source]
merge(other: pyro.ops.streaming.StatsOfDict) pyro.ops.streaming.StatsOfDict[source]
get() Dict[Hashable, Any][source]
Returns

统计数据的字典。此字典的键与更新此对象的样本的键相同。

Return type

dict

class StackStats[source]

基础类:pyro.ops.streaming.StreamingStats

统计收集一系列张量到一个单一的堆叠张量。

update(sample: torch.Tensor) None[source]
merge(other: pyro.ops.streaming.StackStats) pyro.ops.streaming.StackStats[source]
get() Dict[str, Union[int, torch.Tensor]][source]
Returns

一个包含键count: int和(如果已收集任何样本)samples: torch.Tensor的字典。

Return type

dict

class CountStats[source]

基础类:pyro.ops.streaming.StreamingStats

仅统计样本数量的统计跟踪。

例如:

>>> stats = CountStats()
>>> stats.update(torch.randn(3, 3))
>>> stats.get()
{'count': 1}
update(sample) None[source]
merge(other: pyro.ops.streaming.CountStats) pyro.ops.streaming.CountStats[来源]
get() Dict[str, int][source]
Returns

一个带有键count: int的字典。

Return type

dict

class CountMeanStats[source]

基础类:pyro.ops.streaming.StreamingStats

统计跟踪单个torch.Tensor的计数和平均值。

update(sample: torch.Tensor) None[来源]
merge(other: pyro.ops.streaming.CountMeanStats) pyro.ops.streaming.CountMeanStats[来源]
get() Dict[str, Union[int, torch.Tensor]][source]
Returns

一个包含键count: int和(如果已收集任何样本)mean: torch.Tensor的字典。

Return type

dict

class CountMeanVarianceStats[source]

基础类:pyro.ops.streaming.StreamingStats

统计跟踪单个torch.Tensor的计数、均值和(对角线)方差。

update(sample: torch.Tensor) None[source]
merge(other: pyro.ops.streaming.CountMeanVarianceStats) pyro.ops.streaming.CountMeanVarianceStats[source]
get() Dict[str, Union[int, torch.Tensor]][source]
Returns

一个包含键count: int的字典,以及(如果已收集任何样本)mean: torch.Tensorvariance: torch.Tensor

Return type

dict

状态空间模型和GP工具

class MaternKernel(nu=1.5, num_gps=1, length_scale_init=None, kernel_scale_init=None)[源代码]

基础类:pyro.nn.module.PyroModule

提供了将具有Matern核的单变量高斯过程(GPs)表示为状态空间模型的基本构建块。

Parameters
  • nu (float) – Matern核的阶数(0.5、1.5或2.5之一)

  • num_gps (int) – GP的数量

  • length_scale_init (torch.Tensor) – 可选的num_gps维向量,用于初始化长度尺度

  • kernel_scale_init (torch.Tensor) – 可选的 num_gps 维向量,用于初始化核尺度

参考文献

[1] Kalman Filtering and Smoothing Solutions to Temporal Gaussian Process Regression Models,

Jouni Hartikainen 和 Simo Sarkka。

[2] Stochastic Differential Equation Methods for Spatio-Temporal Gaussian Process Regression,

阿尔诺·索林。

transition_matrix(dt)[source]

计算GP潜在空间的(指数化的)转移矩阵。 生成的矩阵布局为(num_gps, old_state, new_state),即此 矩阵从右侧乘以状态。

详情请参见参考文献[1]中的第5节。

Parameters

dt (float) – GP潜在空间演化的时间间隔。

Returns torch.Tensor

一个形状为 (num_gps, state_dim, state_dim) 的转移矩阵的三维张量。

stationary_covariance()[source]

计算稳态协方差。参见参考文献[2]中的公式3.26。

Returns torch.Tensor

一个形状为 (num_gps, state_dim, state_dim) 的协方差矩阵的三维张量。

process_covariance(A)[source]

给定一个使用transition_matrix计算的转移矩阵A,按照参考文献[2]中的公式3.11计算过程协方差。

Returns torch.Tensor

一个形状为 (num_gps, state_dim, state_dim) 的批量协方差矩阵

transition_matrix_and_covariance(dt)[来源]

获取与时间间隔dt对应的转移矩阵和过程协方差。

Parameters

dt (float) – GP潜在空间演化的时间间隔。

Returns tuple

(transition_matrix, process_covariance) 都是形状为 (num_gps, state_dim, state_dim) 的三维张量

training: bool