用Python构建协方差矩阵

Question

问题 我想实现一个算法，从我的主管未发表的论文，作为一部分，我需要建立一个协方差矩阵C使用一些规则在文中给出。我来自Matlab，想借此机会最终学习Python，因此我的问题是：如何在Python（包括numpy和scipy）中以最高效（快速）的方式来学习Python？

子问题1:

• 选项1：我使用2作为循环，循环覆盖所有行和所有列。我想这是最糟糕的事情了。
• 选项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)])

但是，这将构建几个临时数组（在本例中是三个）。如果我们用10^6的观测值处理10000维多项式，上面的方法会使用太多的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

---

总的来说，我只谈到了一些一般的模式。不过，希望这能让你找到正确的方向。