在 SciPy 中生成一个又高又瘦的随机正交矩阵

Question

我需要在SciPy中生成一个又高又细的随机列正交矩阵；也就是说，行数 n 比 p 的列数大很多数量级（比如说 n = 1e5 和 p = 100。我知道 scipy.stats.ortho_group 生成正方形正交矩阵。但是，在我的例子中，生成 n-by-n 随机正交矩阵然后保留第一个 p 列是根本不可行的...是否有更省时且 space- 有效的方法？

Answer 1

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

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

这里q就是你要的矩阵

Answer 2

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

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

Answer 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)

在 SciPy 中生成一个又高又瘦的随机正交矩阵

Generating a tall-and-thin random orthonormal matrix in SciPy

python

numpy

matrix

orthogonal

scipy