如何在Python中生成稀疏正交矩阵？含随机生成方法问询

嘿，这个问题确实戳中了稀疏矩阵和正交矩阵的一个矛盾点——正交性的约束很容易让稀疏矩阵变稠密，常规的QR分解确实不太能保留稀疏性。我来分享几个在Python里可行的方案，从结构化到随机版本都有：

一、构造结构化的稀疏正交矩阵 #

最直接且靠谱的就是置换矩阵，它本身就是正交矩阵（转置等于逆），而且极端稀疏（n阶矩阵只有n个非零元素）。用scipy.sparse可以轻松构造：

import numpy as np
from scipy.sparse import csr_matrix

n = 10  # 你需要的矩阵维度
# 生成一个随机置换序列
perm = np.random.permutation(n)
# 定义非零元素的行、列位置和值
row_indices = np.arange(n)
col_indices = perm
non_zero_vals = np.ones(n)
# 创建稀疏正交矩阵
sparse_perm_ortho = csr_matrix((non_zero_vals, (row_indices, col_indices)), shape=(n, n))

# 验证正交性：转置乘原矩阵应该等于单位矩阵
print((sparse_perm_ortho.T @ sparse_perm_ortho).toarray())  # 输出单位矩阵

如果你需要更灵活的稀疏结构，可以试试块对角稀疏正交矩阵：在对角线上放置多个小的正交子矩阵（比如2x2的旋转矩阵），其余位置全为0，整体依然保持正交性，稀疏度可以通过调整块的大小和数量来控制：

from scipy.sparse import block_diag

n = 10
blocks = []
remaining_dim = n

while remaining_dim > 0:
    if remaining_dim >= 2:
        # 生成随机的2x2旋转正交矩阵
        theta = np.random.uniform(0, 2 * np.pi)
        rot_block = np.array([
            [np.cos(theta), -np.sin(theta)],
            [np.sin(theta), np.cos(theta)]
        ])
        blocks.append(rot_block)
        remaining_dim -= 2
    else:
        # 剩下1维的话直接加单位矩阵块
        blocks.append(np.array([[1]]))
        remaining_dim -= 1

# 构造块对角稀疏正交矩阵
block_sparse_ortho = block_diag(blocks, format='csr')
# 验证正交性
print(np.allclose(block_sparse_ortho.T @ block_sparse_ortho, np.eye(n)))  # 应该返回True

二、生成随机稀疏正交矩阵（非结构化） #

如果需要更一般的非结构化随机稀疏正交矩阵，事情会复杂一些，因为直接对随机稀疏矩阵做QR分解会快速丢失稀疏性。这里有两个实用的思路：

1. 基于稀疏Householder变换的迭代构造 #

Householder变换可以用来构造正交矩阵，如果我们只对稀疏的行/列进行变换，就能尽量保留整体的稀疏性。我们可以从单位矩阵开始，多次应用随机的稀疏Householder变换：

import numpy as np
from scipy.sparse import csr_matrix, eye

def create_random_sparse_householder(n, sparsity=0.1):
    # 生成稀疏的Householder向量
    v = csr_matrix((n, 1))
    # 随机选择非零元素的位置
    non_zero_pos = np.random.choice(n, size=int(n * sparsity), replace=False)
    v[non_zero_pos] = np.random.randn(len(non_zero_pos))
    # 归一化向量
    v_norm = np.linalg.norm(v.toarray())
    v = v / v_norm
    # 构造Householder矩阵：H = I - 2vv^T
    householder_mat = eye(n) - 2 * v @ v.T
    return householder_mat

n = 10
# 从单位矩阵开始，应用几次稀疏Householder变换
random_sparse_ortho = eye(n)
# 迭代次数别太多，否则稀疏性会下降
for _ in range(3):
    h_mat = create_random_sparse_householder(n, sparsity=0.2)
    random_sparse_ortho = random_sparse_ortho @ h_mat

# 验证正交性
print(np.allclose(random_sparse_ortho.T @ random_sparse_ortho, np.eye(n)))
# 查看稀疏度
print(f"矩阵稀疏度: {1 - random_sparse_ortho.nnz/(n*n):.2%}")

这种方法能较好地平衡稀疏性和正交性，不过迭代次数不宜过多，否则稀疏性会逐渐降低。

2. 稠密随机正交矩阵剪枝再正交化 #

另一种折中方案是先生成稠密的随机正交矩阵，然后将绝对值较小的元素置为0（剪枝），再对剪枝后的稀疏矩阵重新做正交化。这种方法适合对正交性精度要求不是极高，但需要强随机性的场景：

n = 10
# 生成稠密的随机正交矩阵
dense_random_ortho = np.linalg.qr(np.random.randn(n, n))[0]
# 设置阈值，把小元素置为0
prune_threshold = 0.3
sparse_pruned = np.where(np.abs(dense_random_ortho) < prune_threshold, 0, dense_random_ortho)
# 转成稀疏矩阵
sparse_pruned_csr = csr_matrix(sparse_pruned)
# 对稀疏矩阵做QR分解重新正交化
from scipy.sparse.linalg import qr
final_sparse_ortho, _ = qr(sparse_pruned_csr)

# 验证正交性（注意精度可能略有损失）
print(np.allclose(final_sparse_ortho.T @ final_sparse_ortho, np.eye(n), atol=1e-3))
# 查看稀疏度
print(f"矩阵稀疏度: {1 - final_sparse_ortho.nnz/(n*n):.2%}")

总结一下 #

• 如果可以接受结构化的稀疏正交矩阵，置换矩阵和块对角正交矩阵是最简单、最可靠的选择；
• 如果需要非结构化的随机版本，稀疏Householder变换迭代是更能保留稀疏性的方法；
• 稠密正交矩阵剪枝再正交化则是一种兼顾随机性和稀疏性的折中方案，适合对精度要求不极端的场景。

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