Python生成行相关系数≈R的标准正态矩阵性能优化问题

原代码的核心问题是显式构造N×N的协方差矩阵，当N=100000时，这个矩阵包含10^10个元素，需要约800GB内存才能存储，完全超出常规硬件的承受范围；同时np.random.multivariate_normal处理这种超大协方差矩阵时计算复杂度极高，直接导致耗时剧增。

针对你需要的行间两两相关系数为R、列间独立、元素服从标准正态分布的矩阵，我们可以利用复合对称协方差结构的数学分解，彻底避免构造超大协方差矩阵，直接生成符合要求的样本。

数学原理 #

行间两两相关系数为R的协方差矩阵可分解为：
$$\Sigma = (1-R)I + R\mathbf{1}\mathbf{1}^T$$
其中$I$是N阶单位矩阵，$\mathbf{1}$是N维全1向量。

基于该分解，样本可通过以下方式生成：
$$X = \sqrt{1-R} \cdot Z + \sqrt{R} \cdot U$$

• $Z$：N×M的独立标准正态分布矩阵（所有元素互相独立）
• $U$：1×M的独立标准正态分布向量（每一列对应一个值，所有行共享该列的U值）

此方法生成的矩阵满足：

1. 任意两行的相关系数为R
2. 任意两列相互独立
3. 每个元素服从标准正态分布（方差为1）

优化后的代码 #

import numpy as np

def draw_effects(N, M, R):
    '''
    创建一个N行M列的矩阵，元素取自标准正态分布，行间相关系数约为R，列间互不相关
    
    参数:
        N: 行数
        M: 列数
        R: 行间相关系数，取值范围[0,1]
    '''
    # 生成独立标准正态矩阵Z
    Z = np.random.normal(size=(N, M))
    # 生成每列共享的标准正态向量U
    U = np.random.normal(size=(1, M))
    
    # 组合得到目标矩阵
    if R == 0:
        return Z
    elif R == 1:
        return np.tile(U, (N, 1))
    else:
        sqrt_1mr = np.sqrt(1 - R)
        sqrt_r = np.sqrt(R)
        return sqrt_1mr * Z + sqrt_r * U

性能对比 #

• 原代码：当N=100000时，无法构造协方差矩阵（内存不足），更无法完成计算
• 优化后代码：生成100000×3000的矩阵仅需数秒（取决于硬件），内存占用约2.4GB（float64类型），完全在常规服务器内存范围内

验证 #

可通过小样本验证行间相关系数的准确性：

# 测试小样本
N_test = 1000
M_test = 100
R_test = 0.5
mat = draw_effects(N_test, M_test, R_test)

# 计算行间相关系数矩阵的均值（排除对角线）
corr_mat = np.corrcoef(mat)
mean_corr = np.mean(corr_mat[np.triu_indices_from(corr_mat, k=1)])
print(f"实际行间平均相关系数: {mean_corr:.4f}")

输出结果会接近设定的R值，验证方法的正确性。

内容的提问来源于stack exchange，提问作者Matthew Howell