我需要为我的工作生成大量随机均值不变正交矩阵。均值不变矩阵具有属性A*1_n=1_n，其中 1_n 是标量 1 的大小为 n 的向量，基本为np.ones(n)。我使用 Python，特别是 Numpy 来创建我的矩阵，但我想确保我的方法既正确又最有效。此外，我想介绍我尝试过的 3 种独立正交化方法的发现，并希望能解释为什么一种方法比其他方法更快。我在帖子的末尾就我的发现提出了四个问题。
一般来说，为了创建一个均值不变的随机正交矩阵 A，您需要创建一个随机方阵 M1，将其第一列替换为一列并正交化矩阵。然后，使用另一个矩阵 M2 再次执行此操作，最终均值不变随机正交矩阵为 A = M1*(M2.T)。这个过程的瓶颈是正交化。正交化的主要方法有 3 种，即使用投影的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 都执行。
对于 10,000x10,000 随机均值不变正交矩阵，此过程在我的计算机上大约需要1 小时（8 核 @3.7GHz，16GB RAM，512GB SSD）。
我发现代替 Gram-Schmidt 过程，我可以用以下方法对 Numpy 中的矩阵进行正交化：
q1, r1 = np.linalg.qr(M1)
其中 q1 是正交化矩阵，r1 是上三角矩阵（我们不需要保留 r1）。我对 M2 做同样的事情并得到 q2。那么，A=q1*(q2.T)。对于相同的 10,000x10,000 矩阵，此过程在同一台计算机上大约需要70 秒。我相信linalg.qr()图书馆使用 Householder 转换，但我希望有人确认它。
最后，我尝试改变初始随机矩阵 M1 和 M2 的生成方式。而不是
M1 = np.random.normal(size=(n,n))我使用 Dirichlet 分布：
M1 = np.random.default_rng().dirichlet(np.ones(n),size=(n)). 然后，我linalg.qr()之前使用过类似的东西，我得到了 10000x10000 的矩阵M1 = np.random.normal(size=(n,n))。
我的问题是：
- Numpy 的
np.linalg.qr()方法是否真的使用了 Householder 变换？或者可能是吉文斯轮换？ - 为什么 Gram-Schmidt 过程比 慢得多
np.linalg.qr()？ - 我知道狄利克雷过程会产生一个几乎正交的矩阵。是不是因为我们正在创建一万维，所以很可能随机得到一个与其他所有向量正交的向量？
np.linalg.qr()不关心矩阵与正交性的接近程度。 - 有没有更快的方法来生成随机正交矩阵？我可以对我的代码进行任何优化以使其更快/更高效吗？
编辑：cp.linalg.qr()在我的 2080ti 上，cupy 在同一个 10,000x10,000 随机矩阵上只需要 16 秒，而不是 CPU 上的 70 秒（8 核 @3.7GHz，多线程，16GB RAM 和 512GB SSD）。