在 Python 中构建协方差矩阵
Building a covariance matrix in Python
问题
我想从我的主管未发表的论文中实现一个算法，作为其中的一部分，我需要使用论文中给出的一些规则构造一个协方差矩阵 C。我来自 Matlab，想借此机会最终学习 Python，因此我的问题是：如何在 Python（包括 numpy、scipy)?
子问题1：
- 选项 1：我使用 2 个 for 循环，遍历所有行和所有列。我认为这是最糟糕的事情。
- 选项 2：使用列表推导式，我构建了一个欧氏对列表，然后遍历该列表。这就是我现在正在做的。
有没有更好的方法？
子问题2
- 选项 1：我遍历矩阵中的所有元素。
- 选项 2：我只迭代下三角部分（没有对角线），然后添加转置（因为协方差矩阵是对称的），然后添加对角线。
我相当确信子问题 1 很简单，但我不知道子问题 2。我可能还应该说我正在处理的矩阵可能是 2*10^4 x 2* 10^4.
谢谢！
编辑
我不想给出实际的协方差矩阵，但由于人们想要一个例子，假设我们想要构建一个名为 'Brownian bridge' 的随机过程的协方差矩阵。它的结构如下：
cov(Xs, Xt) = min{s,t} − st
假设 s,t ∈ {1,...,100}。你会如何建造它？
首先，对于将来可能会遇到此问题的其他人：如果您确实有数据并且想要估计协方差矩阵，正如几个人指出的那样，请使用 np.cov 或类似的东西。
从模式构建数组
但是，您的问题是关于如何在给定一些预定义规则的情况下构建大型矩阵。为了消除评论中的一些混乱：你的问题似乎不是关于估计协方差矩阵，而是关于指定一个。换句话说，你问的是如何根据一些预定义的规则构建一个大数组。
哪种方式最有效将取决于您正在做的细节。在这种情况下，大多数性能技巧将涉及在您正在执行的计算中利用对称性。 （例如，一行是否相同？）
如果不确切知道自己在做什么，很难说出任何具体的内容。因此，我将重点介绍一般情况下如何做这类事情。 （注意：我刚刚注意到您的编辑。稍后我将包含一个布朗桥的示例...）
恒定（或简单）Row/Column
最基本的情况是输出数组中的常量行或列。使用切片语法很容易创建数组并将值分配给列或行：
import numpy as np
num_vars = 10**4
cov = np.zeros((num_vars, num_vars), dtype=float)
设置整个column/row:
# 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
您还可以将数组分配给 columns/rows:
# 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 "inverted cone" 的值：
# 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
总的来说，我只谈到了一些一般模式。但是，希望这能让您指明正确的方向。
问题 我想从我的主管未发表的论文中实现一个算法，作为其中的一部分，我需要使用论文中给出的一些规则构造一个协方差矩阵 C。我来自 Matlab，想借此机会最终学习 Python，因此我的问题是：如何在 Python（包括 numpy、scipy)?
子问题1：
- 选项 1：我使用 2 个 for 循环，遍历所有行和所有列。我认为这是最糟糕的事情。
- 选项 2：使用列表推导式，我构建了一个欧氏对列表，然后遍历该列表。这就是我现在正在做的。
有没有更好的方法？
子问题2
- 选项 1：我遍历矩阵中的所有元素。
- 选项 2：我只迭代下三角部分（没有对角线），然后添加转置（因为协方差矩阵是对称的），然后添加对角线。
我相当确信子问题 1 很简单，但我不知道子问题 2。我可能还应该说我正在处理的矩阵可能是 2*10^4 x 2* 10^4.
谢谢！
编辑 我不想给出实际的协方差矩阵，但由于人们想要一个例子，假设我们想要构建一个名为 'Brownian bridge' 的随机过程的协方差矩阵。它的结构如下：
cov(Xs, Xt) = min{s,t} − st
假设 s,t ∈ {1,...,100}。你会如何建造它？
首先，对于将来可能会遇到此问题的其他人：如果您确实有数据并且想要估计协方差矩阵，正如几个人指出的那样，请使用 np.cov 或类似的东西。
从模式构建数组
但是，您的问题是关于如何在给定一些预定义规则的情况下构建大型矩阵。为了消除评论中的一些混乱：你的问题似乎不是关于估计协方差矩阵，而是关于指定一个。换句话说，你问的是如何根据一些预定义的规则构建一个大数组。
哪种方式最有效将取决于您正在做的细节。在这种情况下，大多数性能技巧将涉及在您正在执行的计算中利用对称性。 （例如，一行是否相同？）
如果不确切知道自己在做什么，很难说出任何具体的内容。因此，我将重点介绍一般情况下如何做这类事情。 （注意：我刚刚注意到您的编辑。稍后我将包含一个布朗桥的示例...）
恒定（或简单）Row/Column
最基本的情况是输出数组中的常量行或列。使用切片语法很容易创建数组并将值分配给列或行：
import numpy as np
num_vars = 10**4
cov = np.zeros((num_vars, num_vars), dtype=float)
设置整个column/row:
# 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
您还可以将数组分配给 columns/rows:
# 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 "inverted cone" 的值：
# 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
总的来说，我只谈到了一些一般模式。但是，希望这能让您指明正确的方向。