Python中随机正交矩阵的高效生成

问题描述

我需要为我的工作生成大量随机均值不变的正交矩阵。均值不变矩阵具有 A*1_n=1_n 属性，其中 1_n 是标量 1 的大小为 n 的向量，基本上为 np.ones(n)。我使用 Python，特别是 Numpy 来创建我的矩阵，但我想确保我的方法既正确又最有效。此外，我想展示我尝试过的 3 种不同正交化方法的发现，并希望得到关于为什么一种方法比其他方法更快的解释。我在帖子的最后问了四个关于我的发现的问题。

一般来说，为了创建一个均值不变的随机正交矩阵 A，你需要创建一个随机方阵 M1，用一列 1 替换它的第一列并对矩阵进行正交化。然后，您再次使用另一个矩阵 M2 执行此操作，最终的均值不变随机正交矩阵为 A = M1*(M2.T)。这个过程的瓶颈是正交化。正交化有 3 种主要方式，即使用投影的 Gram–Schmidt process、使用反射的 Householder transformation 和 Givens rotation。

使用 Numpy 创建 nxn 随机矩阵非常简单： M1 = np.random.normal(size=(n,n))。然后，我用 1_n 替换 M1 的第一列。

据我所知，Gram–Schmidt 过程在任何非常流行的库中都不存在，所以我发现这段代码运行良好：

def gram_schmidt(M1):
    """Orthogonalize a set of vectors stored as the columns of matrix M1."""
    # Get the number of vectors.
    n = M1.shape[1]
    for j in range(n):
        # To orthogonalize the vector in column j with respect to the
        # prevIoUs vectors,subtract from it its projection onto
        # each of the prevIoUs vectors.
        for k in range(j):
            M1[:,j] -= np.dot(M1[:,k],M1[:,j]) * M1[:,k]
        M1[:,j] = M1[:,j] / np.linalg.norm(M1[:,j])
    return M1

显然，上面的代码必须对 M1 和 M2 都执行。

对于 10,000x10,000 随机均值不变正交矩阵，此过程在我的计算机（8 核 @3.7GHz、16GB RAM、512GB SSD）上大约需要 1 小时。

我发现我可以使用以下方法对 Numpy 中的矩阵进行正交化，而不是 Gram-Schmidt 过程： q1,r1 = np.linalg.qr(M1) 其中 q1 是正交矩阵，r1 是上三角矩阵（我们不需要保留 r1）。我对 M2 做同样的事情并得到 q2。那么，A=q1*(q2.T)。对于相同的 10,000 矩阵，此过程在同一台计算机上大约需要 70 秒。我相信 linalg.qr() 库使用了 Householder 转换，但我希望有人确认。

最后，我尝试改变生成初始随机矩阵 M1 和 M2 的方式。代替 M1 = np.random.normal(size=(n,n)) 我使用了狄利克雷分布： M1 = np.random.default_rng().dirichlet(np.ones(n),size=(n))。然后，我像以前一样使用了 linalg.qr()，并在与 M1 = np.random.normal(size=(n,n)) 大致相同的时间内得到了 10000x10000 矩阵。

我的问题是：

Numpy 的 np.linalg.qr() 方法是否实际使用了 Householder 转换？或者可能是 Givens 轮换？
为什么 Gram-Schmidt 过程比 np.linalg.qr() 慢这么多？
我知道狄利克雷过程产生了一个几乎正交的矩阵。是不是因为我们创建了 10,000 个维度，所以很可能随机得到一个与所有其他维度正交的向量？ np.linalg.qr() 不关心矩阵与正交性的接近程度。
有没有更快的方法来生成随机正交矩阵？是否可以对代码进行任何优化以使其更快/更高效？

编辑：同一个 10,000 随机矩阵上的丘比 cp.linalg.qr() 在我的 2080ti 上只需要 16 秒，而不是 cpu 上的 numpy 需要 70 秒（8 核 @3.7GHz，多线程，16GB RAM和 512GB SSD）。

解决方法

这是制造此类矩阵的另一种方法。我不知道分布是什么样的，但它可能比你描述的方法更快。

首先，这里的家庭反射器是一个矩阵形式

H = I - 2*h*h'/(h'*h)

其中 h 是一个向量。

注意：

H 是正交的，对称的

H 可以应用于 O(dim) 中的向量

如果 x 和 y 是具有相同范数的任意两个向量，那么我们可以找到这样一个矩阵来将 x 映射到 y（h 是 x-y）

制造“随机”正交矩阵的一种方法是计算“随机”Householder 矩阵的乘积

如果 Q 是一个正交矩阵，而 u 是全 1 的向量，并且

q = Q*u

那么 q 和 u 具有相同的范数，所以如果 H 是将 q 映射到 u 的家庭反射器，

R = H*Q is orthogonal
R*u = H*Q*u = H*q = u

dirichlet matrix matrix numpy orthogonal