在 SciPy 中生成高瘦随机正交矩阵答案

【问题标题】：Generating a tall-and-thin random orthonormal matrix in SciPy在 SciPy 中生成高瘦随机正交矩阵
【发布时间】：2018-12-06 01:40:30
【问题描述】：

我需要在 SciPy 中生成一个高而细的随机列正交矩阵；也就是说，n 的行数远远大于p 的列数多个数量级（比如n = 1e5 和p = 100。我知道scipy.stats.ortho_group 生成一个正交正交矩阵。但是，在我的情况下，生成一个n-by-n 随机正交矩阵然后保留第一个p 列是根本不可行的......有没有更节省时间和空间的方法？

【问题讨论】：

对不起，请刷新我的记忆。正交矩阵不应该是正方形的吗？
@VivekKalyanarangan 是的，但是矩形矩阵可以是行正交或列正交。

标签： python numpy matrix scipy orthogonal

【解决方案1】：

可以先生成一个又高又瘦的随机矩阵，然后进行qr分解。

a = np.random.random(size=(100000, 100))
q, _ = np.linalg.qr(a)

这里q是你想要的矩阵。

【讨论】：

啊，好主意。尽管 a 应该是随机均匀或随机正态以使 q 成为真正的随机正交矩阵是否重要？
我认为正态分布可能会更好，但我并不擅长这个领域。您可能会发现这篇论文arxiv.org/pdf/math-ph/0609050.pdf 很有帮助，尤其是第 4 节，其中详细讨论了此类问题。
请注意，这不适用于短而宽的随机矩阵，例如>>> a = np.random.normal(0,1,(20,100)); >>> q, _ = np.linalg.qr(a); >>> q.shape; (20, 20)
@ScottGigante 短而宽的矩阵如何成为“列正交”？列数大于空间维度。

【解决方案2】：

对我来说scipy.linalg.orth 比numpy.linalg.qr 快一点：

a = np.random.random(size=(100000, 100))
q = scipy.linalg.orth(a)

【讨论】：

【解决方案3】：

这是一个基准答案。请注意，我进行了一些转置，以便无论该矩阵是高而薄（给出列正交）还是短而宽（给出行正交），这都会起作用。

def qr_method(n, m):
    X = np.random.normal(0,1,(n,m))
    if n < m:
        X = X.T
    Q, _ = np.linalg.qr(X)
    if n < m:
        Q = Q.T
    return Q

def orth_method(n, m):
    X = np.random.normal(0,1,(n,m))
    if n < m:
        X = X.T
    Q = scipy.linalg.orth(X)
    if n < m:
        Q = Q.T
    return Q

def ortho_group_method(n, m):
    Q = scipy.stats.ortho_group.rvs(max(n, m))[:min(n, m),:]
    if m < n:
        Q = Q.T
    return Q

ortho_group 方法（也就是制作一个方阵然后取一个子集）太慢了，以至于我没有将它与其他方法一起进行基准测试：

%timeit ortho_group_method(500, 20)
2.73 s ± 57.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

另外两个，差别可以忽略不计，通过QR略快。

%timeit qr_method(10000, 200)
168 ms ± 3.78 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit orth_method(10000, 200)
193 ms ± 4.09 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

矩阵的高度有区别吗？对于一个非常高的矩阵，它们几乎是等价的。

%timeit qr_method(100000, 20)
122 ms ± 1.1 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit orth_method(100000, 20)
130 ms ± 6.5 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

对于方阵，QR 要快得多。

%timeit qr_method(500, 500)
47.5 ms ± 202 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit orth_method(500, 500)
137 ms ± 1.32 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

【讨论】：