如何在 python numpy 中创建随机正交矩阵

Question

在 python 中是否可以调用一种方法来创建随机正交矩阵？可能使用 numpy？或者有没有办法使用多种 numpy 方法创建正交矩阵？谢谢

Answer 1

scipy 的 0.18 版有 scipy.stats.ortho_group and scipy.stats.special_ortho_group. The pull request where it was added is https://github.com/scipy/scipy/pull/5622

例如，

In [24]: from scipy.stats import ortho_group  # Requires version 0.18 of scipy

In [25]: m = ortho_group.rvs(dim=3)

In [26]: m
Out[26]: 
array([[-0.23939017,  0.58743526, -0.77305379],
       [ 0.81921268, -0.30515101, -0.48556508],
       [-0.52113619, -0.74953498, -0.40818426]])

In [27]: np.set_printoptions(suppress=True)

In [28]: m.dot(m.T)
Out[28]: 
array([[ 1.,  0., -0.],
       [ 0.,  1.,  0.],
       [-0.,  0.,  1.]])

Answer 2

这是从 https://github.com/scipy/scipy/pull/5622/files 中提取的 rvs 方法，变化很小 - 足以运行作为一个独立的 numpy 函数。

import numpy as np    

def rvs(dim=3):
     random_state = np.random
     H = np.eye(dim)
     D = np.ones((dim,))
     for n in range(1, dim):
         x = random_state.normal(size=(dim-n+1,))
         D[n-1] = np.sign(x[0])
         x[0] -= D[n-1]*np.sqrt((x*x).sum())
         # Householder transformation
         Hx = (np.eye(dim-n+1) - 2.*np.outer(x, x)/(x*x).sum())
         mat = np.eye(dim)
         mat[n-1:, n-1:] = Hx
         H = np.dot(H, mat)
         # Fix the last sign such that the determinant is 1
     D[-1] = (-1)**(1-(dim % 2))*D.prod()
     # Equivalent to np.dot(np.diag(D), H) but faster, apparently
     H = (D*H.T).T
     return H

它符合 Warren 的测试，

Answer 3

您可以通过对 QR 进行 QR 分解来获得随机 n x n 正交矩阵 Q，（均匀分布在 n x n 正交矩阵的流形上） =13=] 矩阵，元素为 i.i.d。均值 0 和方差 1 的高斯随机变量。这是一个例子：

import numpy as np
from scipy.linalg import qr

n = 3
H = np.random.randn(n, n)
Q, R = qr(H)

print (Q.dot(Q.T))

[[  1.00000000e+00  -2.77555756e-17   2.49800181e-16]
 [ -2.77555756e-17   1.00000000e+00  -1.38777878e-17]
 [  2.49800181e-16  -1.38777878e-17   1.00000000e+00]]

编辑：（在@g g 发表评论后重新审视这个答案。）上面关于提供均匀分布（在所谓的 Stiefel 流形上）正交矩阵的高斯矩阵的 QR 分解的声明是由this reference 的定理 2.3.18-19。请注意，结果的陈述表明 "QR-like" 分解，但是 三角矩阵 R 具有正元素 .

显然，scipy (numpy) 函数 的 qr 函数不保证 R 的正对角线元素和相应的 Q 实际上不是均匀分布。这已在 this 专着第 1 节中观察到。 4.6（讨论指的是MATLAB，但我猜MATLAB和scipy都使用相同的LAPACK例程）。那里建议将 qr 提供的矩阵 Q 修改为 post 将其与随机酉对角矩阵相乘。

下面我重现了上述参考文献中的实验，绘制了 qr 提供的 "direct" Q 矩阵的特征值相位的经验分布（直方图），以及"modified" 版本，其中可以看出修改后的版本确实具有统一的特征值相位，正如均匀分布的正交矩阵所预期的那样。

from scipy.linalg import qr, eigvals
from seaborn import distplot

n = 50
repeats = 10000

angles = []
angles_modified = []
for rp in range(repeats):
    H = np.random.randn(n, n)
    Q, R = qr(H)
    angles.append(np.angle(eigvals(Q)))
    Q_modified = Q @ np.diag(np.exp(1j * np.pi * 2 * np.random.rand(n)))
    angles_modified.append(np.angle(eigvals(Q_modified))) 

fig, ax = plt.subplots(1,2, figsize = (10,3))
distplot(np.asarray(angles).flatten(),kde = False, hist_kws=dict(edgecolor="k", linewidth=2), ax= ax[0])
ax[0].set(xlabel='phase', title='direct')
distplot(np.asarray(angles_modified).flatten(),kde = False, hist_kws=dict(edgecolor="k", linewidth=2), ax= ax[1])
ax[1].set(xlabel='phase', title='modified');

Answer 4

如果您想要 none 具有正交列向量的方阵，您可以使用上述任何一种方法创建一个方阵并删除一些列。

Answer 5

创建任意形状 (n x m) 正交矩阵的简单方法：

import numpy as np

n, m = 3, 5

H = np.random.rand(n, m)
u, s, vh = np.linalg.svd(H, full_matrices=False)
mat = u @ vh

print(mat @ mat.T) # -> eye(n)

注意，如果n > m，会得到mat.T @ mat = eye(m)。

Answer 6

from scipy.stats import special_ortho_group
num_dim=3
x = special_ortho_group.rvs(num_dim)

Documentation

Answer 7

Numpy 也有 qr 分解。 https://numpy.org/doc/stable/reference/generated/numpy.linalg.qr.html

import numpy as np

a = np.random.rand(3, 3)
q, r = np.linalg.qr(a)

q @ q.T
# array([[ 1.00000000e+00,  8.83206468e-17,  2.69154044e-16],
#        [ 8.83206468e-17,  1.00000000e+00, -1.30466244e-16],
#        [ 2.69154044e-16, -1.30466244e-16,  1.00000000e+00]])

如何在 python numpy 中创建随机正交矩阵

How to create random orthonormal matrix in python numpy

python

numpy

linear-algebra

orthogonal