在Python中构建协方差矩阵

Question

问题 我想通过我的主管从未发表的论文中实现算法,作为其中的一部分,我需要使用本文中给出的一些规则来构建协方差矩阵C. 我来自Matlab并希望借此机会最终学习Python,因此我的问题是:我如何以最有效(快速)的方式在Python(包括numpy,scipy)中做到这一点？

子问题1:

• 选项1:我使用2 for循环,循环遍历所有行和所有列.我认为这是最糟糕的事情.
• 选项2:使用列表推导,我构建了一个欧几里德对列表,然后迭代该列表.这就是我现在正在做的事情.

有没有更好的方法？

子问题2

• 选项1:我迭代矩阵中的所有元素.
• 选项2:我只在下三角部分(无对角线)上迭代,然后添加转置(因为协方差矩阵是对称的),然后添加对角线.

我相信子问题1是不费脑筋但我不知道子问题2.我可能也应该说我正在处理的矩阵可能是2*10 ^ 4 x 2*10 ^ 4.

谢谢!

编辑 我不想给出实际的协方差矩阵,但由于人们想要一个例子,假设我们想要构造一个称为"布朗桥"的随机过程的协方差矩阵.它的结构由下式给出:

cov(Xs,Xt)= min {s,t} - st

因为我们说s,t∈{1,...,100}.你会如何建造它？

Answer 1

首先,对于将​​来可能会遇到这个问题的其他人:如果你确实有数据,并且想要估计协方差矩阵,正如几个人所指出的那样,使用np.cov或类似的东西.

从模式构建阵列

但是,您的问题是如何在给定一些预定义规则的情况下构建大型矩阵.为了澄清评论中的一些混淆:你的问题似乎不是关于估计协方差矩阵,而是关于指定一个.换句话说,你问的是如何在给定一些预定义规则的情况下构建一个大型数组.

哪种方式最有效取决于你正在做的细节.在这种情况下,大多数性能技巧将涉及在您正在执行的计算中利用对称性.(例如,一行是否相同？)

如果不确切知道自己在做什么,就很难说清楚.因此,我将重点关注如何做这类事情.(注意:我刚注意到你的编辑.我将在稍后的例子中包含一个布朗桥的例子......)

常量(或简单)行/列

最基本的情况是输出数组中的常量行或列.使用切片语法可以轻松创建数组并为列或行指定值:

import numpy as np

num_vars = 10**4
cov = np.zeros((num_vars, num_vars), dtype=float)

要设置整个列/行:

# Third column will be all 9's
cov[:,2] = 9

# Second row will be all 1's (will overwrite the 9 in col3)
cov[1,:] = 1

您还可以将数组分配给列/行:

# 5th row will have random values
cov[4,:] = np.random.random(num_vars)

# 6th row will have a simple geometric sequence
cov[5,:] = np.arange(num_vars)**2

堆叠阵列

在许多情况下,(但可能不是这种情况)您需要从现有数组构建输出.您可以使用vstack/ hstack/ column_stack/ tile和许多其他类似的功能.

一个很好的例子是,如果我们为多项式的线性反演建立一个矩阵:

import numpy as np

num = 10
x = np.random.random(num) # Observation locations

# "Green's functions" for a second-order polynomial
# at our observed locations
A = np.column_stack([x**i for i in range(3)])

但是,这将构建几个临时数组(在本例中为三个).如果我们使用10000维多项式进行10 ^ 6次观测,则上述方法将使用太多RAM.因此,您可能会迭代列而不是:

ndim = 2
A = np.zeros((x.size, ndim + 1), dtype=float)
for j in range(ndim + 1):
    A[:,j] = x**j

在大多数情况下,不要担心临时数组.colum_stack除非您使用相对较大的数组,否则基于该示例是正确的方法.

最通用的方法

没有任何更多的信息,我们不能利用任何形式的对称性.最通用的方法是迭代.通常你会想要避免这种方法,但有时它是不可避免的(特别是如果计算取决于先前的值).

速度方面,这与嵌套for循环相同,但使用np.ndindex而不是多个for循环更容易(特别是对于> 2D数组):

import numpy as np

num_vars = 10**4
cov = np.zeros((num_vars, num_vars), dtype=float)
for i, j in np.ndindex(cov.shape):
    # Logic presumably in some function...
    cov[i, j] = calculate_value(i, j)

矢量基于索引的计算

如果情况很多,您可以对基于索引的计算进行矢量化.换句话说,直接在输出索引的数组上操作.

假设我们的代码看起来像:

import numpy as np

cov = np.zeros((10, 10)), dtype=float)
for i, j in np.ndindex(cov.shape):
    cov[i,j] = i*j - i

我们可以用以下内容代替:

i, j = np.mgrid[:10, :10]
cov = i*j - i

再举一个例子,让我们建立一个100 x 100"倒锥"的值:

# The complex numbers in "mgrid" give the number of increments
# mgrid[min:max:num*1j, min:max:num*1j] is similar to
# meshgrid(linspace(min, max, num), linspace(min, max, num))
y, x = np.mgrid[-5:5:100j, -5:5:100j]

# Our "inverted cone" is just the distance from 0
r = np.hypot(x, y)

布朗桥

这是一个很容易被矢量化的例子.如果我正确地阅读你的例子,你会想要类似的东西:

import numpy as np

st = np.mgrid[1:101, 1:101]
s, t = st
cov = st.min(axis=0) - s * t

总的来说,我只涉及一些一般模式.但是,希望这会让你指向正确的方向.