我需要为我的工作生成许多随机均值不变正交矩阵。平均不变矩阵具有性质A*1_n=1_n，其中1_n是标量1的大小为n的向量，基本上是np.ones(n)。我使用Python，特别是Numpy来创建矩阵，但我想确保我的方法既正确又高效。此外，我想介绍我对我尝试过的3种不同正交化方法的发现，并希望能解释为什么一种方法比其他方法更快。关于我的发现，我在帖子的最后问了四个问题。
通常，为了创建均值不变的随机正交矩阵a，您需要创建一个随机平方矩阵M1，将其第一列替换为一列1，并使矩阵正交。然后，用另一个矩阵M2再做一次，最终的均值不变随机正交矩阵是A=M1*(M2.T)。这个过程的瓶颈是正交化。正交化有三种主要方式，即使用投影的Gram–Schmidt过程、使用反射的Householder变换和Givens旋转。
使用Numpy创建nxn随机矩阵非常简单：M1 = np.random.normal(size=(n,n))。然后，我将M1的第一列替换为1_n。
据我所知，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执行。
对于10000x10000随机均值不变正交矩阵，此过程在我的计算机上大约需要1小时(8核@3.7GHz，16GB RAM，512GB SSD)。
我发现，我可以用以下方法来正交Numpy中的矩阵，而不是Gram-Schmidt过程：q1, r1 = np.linalg.qr(M1)其中q1是正交矩阵，r1是上三角矩阵(我们不需要保留r1)。我对M2做同样的操作，得到q2。然后，A=q1*(q2.T)。在同一台计算机上，对于相同的10000x10000矩阵，此过程大约需要70秒。我相信linalg.qr()库使用了Householder转换，但我希望有人能证实这一点。
最后，我试图改变产生初始随机矩阵M1和M2的方式。而不是M1 = np.random.normal(size=(n,n))I使用了狄利克雷分布：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()慢得多 - 我知道狄利克雷过程产生了一个几乎正交的矩阵。是因为我们正在创建10000个维度，所以很可能随机得到一个与所有其他维度正交的向量吗？CCD_ 12不关心矩阵与正交性的接近程度
- 有没有更快的方法生成随机正交矩阵？我可以对代码进行任何优化以使其更快/更高效吗
编辑：同一10000x10000随机矩阵上的cupy的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