【发布时间】：2021-03-09 01:09:35
【问题描述】：
我需要为我的工作生成大量随机均值不变正交矩阵。均值不变矩阵具有属性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 process、使用反射的Householder transformation和Givens rotation。
使用 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 随机均值不变正交矩阵，此过程在我的计算机（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,000x10,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 转换？或者可能是吉文斯轮换？ - 为什么 Gram-Schmidt 过程比
np.linalg.qr()慢这么多？ - 我知道 Dirichlet 过程会产生一个几乎正交的矩阵。是不是因为我们正在创建一万维，所以很可能随机得到一个与其他所有向量正交的向量？
np.linalg.qr()不关心矩阵与正交性的接近程度。 - 是否有更快的方法来生成随机正交矩阵？我可以对我的代码进行任何优化以使其更快/更高效吗？
编辑：cupy 在同一个 10,000x10,000 随机矩阵上的 cp.linalg.qr() 在我的 2080ti 上只需要 16 秒，而不是 CPU 上的 numpy 需要 70 秒（8 核 @3.7GHz，多线程，16GB RAM 和 512GB固态硬盘）。
【问题讨论】：
标签： numpy matrix orthogonal dirichlet