如何在SciPy中高效生成高瘦型随机列正交矩阵

Great question—when you're dealing with a matrix where rows n are way larger than columns p (like n=1e5, p=100), generating a full n×n orthogonal matrix and slicing the first p columns is totally impractical (that would require ~80GB of memory for n=1e5!). Instead, there's a much more memory-efficient and faster approach using QR decomposition, which gives you exactly the tall-skinny column-orthogonal matrix you need.

The QR Decomposition Method #

The key insight here is that if you generate a random n×p Gaussian matrix (columns are independent standard normals), performing a QR decomposition on it will produce an n×p matrix Q with orthonormal columns. This Q is uniformly distributed over the Stiefel manifold (the set of all n×p column-orthogonal matrices), which matches the statistical property you'd get from slicing a full orthogonal matrix—without the massive overhead.

Here's the code to implement this:

import numpy as np
from scipy.linalg import qr

# Define dimensions
n = 100000
p = 100

# Generate a random n×p Gaussian matrix (columns are iid standard normals)
random_matrix = np.random.randn(n, p)

# Perform "economic" QR decomposition to get n×p orthonormal Q matrix
Q, R = qr(random_matrix, mode='economic')

# Verify orthogonality (Q^T Q should be the identity matrix, up to numerical precision)
print(np.allclose(Q.T @ Q, np.eye(p)))  # Output: True

Why This Works #

• Memory Efficiency: The Gaussian matrix and resulting Q are both n×p, so for n=1e5 and p=100, that's only ~80MB of memory (using 64-bit floats)—a tiny fraction of the 80GB needed for a full n×n matrix.
• Speed: The QR decomposition runs in O(n p²) time. With p=100, this is computationally trivial even for large n, since the quadratic term is only over the small column dimension.
• Statistical Correctness: The resulting Q has the same uniform distribution over column-orthogonal matrices as slicing the first p columns of a full n×n orthogonal matrix from scipy.stats.ortho_group.

Alternative: Gram-Schmidt Orthogonalization #

While QR decomposition is the recommended approach (it's numerically more stable), you could also implement modified Gram-Schmidt orthogonalization directly on the random Gaussian matrix. However, this requires more manual code and is less efficient than using the optimized QR implementation in SciPy/NumPy, so it's generally not necessary.

内容的提问来源于stack exchange，提问作者p-value