基于稀疏协方差矩阵的正态分布采样优化方案问询

当处理1Mx1M这类大规模稀疏协方差矩阵时，直接转成全矩阵采样会耗尽内存，必须利用稀疏性或矩阵结构来优化。以下是两种实用的Python实现方案：

方案一：稀疏Cholesky分解 + 标准高斯采样 #

高斯分布采样的核心逻辑是：若目标分布为$X \sim \mathcal{N}(0, S)$，则$X = LZ$，其中$L$是$S$的Cholesky分解结果（满足$S = LL^T$），$Z$是标准正态分布$\mathcal{N}(0, I)$的采样向量。利用稀疏矩阵的Cholesky分解，可避免生成全量矩阵：

import numpy as np
from scipy import sparse
from scipy.sparse.linalg import cholesky

# 构建1M规模的三对角稀疏协方差矩阵
n = 10**6
main_diag = np.full(n, 1.0)
off_diag = 0.1 * np.ones(n-1)
S = sparse.diags([main_diag, off_diag, off_diag], [0, -1, 1], format='csc')

# 对稀疏矩阵做Cholesky分解，得到下三角稀疏矩阵L
L = cholesky(S, lower=True)

# 采样标准正态分布向量
Z = np.random.normal(size=n)

# 稀疏矩阵与稠密向量相乘，得到目标采样结果
X = L.dot(Z)

该方案的优势在于：稀疏Cholesky分解会保留原矩阵的稀疏结构，内存占用仅为O(n)量级（远低于全矩阵的O(n²)），矩阵乘法操作也基于稀疏优化，效率极高。

方案二：针对三对角矩阵的结构化优化算法 #

如果你的协方差矩阵是三对角这类具有强结构化的稀疏矩阵，可以直接利用矩阵特性手动实现Cholesky分解和采样递推，进一步压缩内存并提升速度：

import numpy as np

n = 10**6
# 定义三对角协方差矩阵的主对角线和次对角线元素
d = np.full(n, 1.0)  # S[i,i]
e = 0.1 * np.ones(n-1)  # S[i,i+1] = S[i+1,i]

# 手动计算三对角矩阵的Cholesky分解，得到L的主对角线l和次对角线m
l = np.zeros(n)
m = np.zeros(n-1)
l[0] = np.sqrt(d[0])
for i in range(1, n):
    m[i-1] = e[i-1] / l[i-1]
    l[i] = np.sqrt(d[i] - m[i-1]**2)

# 采样标准正态分布向量
Z = np.random.normal(size=n)

# 递推计算采样结果X = L@Z
X = np.zeros(n)
X[0] = l[0] * Z[0]
for i in range(1, n):
    X[i] = l[i] * Z[i] + m[i-1] * X[i-1]

额外优化建议 #

• 若要进一步加速递推过程，可使用numba对循环进行JIT编译，1M规模的计算速度可提升数倍。
• 确保协方差矩阵是正定的：对于三对角矩阵，需满足主对角线元素大于相邻次对角线元素的绝对值之和（本例中1 > 0.1+0.1，满足正定条件），否则Cholesky分解会失败。

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