如何在python numpy中创建随机正交矩阵
Posted
技术标签:
【中文标题】如何在python numpy中创建随机正交矩阵【英文标题】：How to create random orthonormal matrix in python numpy 【发布时间】：2016-11-20 10:46:28 【问题描述】：有没有可以调用的方法在 python 中创建随机正交矩阵？可能使用numpy？或者有没有办法使用多个 numpy 方法创建一个正交矩阵？谢谢。
【问题讨论】：
【参考方案1】：scipy 0.18 版具有scipy.stats.ortho_group 和scipy.stats.special_ortho_group。添加它的拉取请求是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.]])
【讨论】：
感谢您的回复。我注意到给出的答案都是关于方阵的。我还能用这个方法得到一个d x k矩阵，其中k 构造一个 d x d 矩阵，然后删除多余的列【参考方案2】：您可以通过对具有元素 i.i.d 的n x n 矩阵执行QR 分解来获得随机的n x n 正交矩阵Q，（均匀分布在n x n 正交矩阵的流形上）。均值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]]
编辑：（在@gg g 发表评论后重新访问此答案。）定理 2.3 建议了上述关于提供均匀分布（在所谓的 Stiefel 流形上）正交矩阵的高斯矩阵的 QR 分解的声明this reference 的 .18-19。请注意，结果的陈述表明“类似 QR”的分解，但是，三角矩阵 R 具有正元素。
显然，scipy(numpy)函数的qr函数不保证R的正对角元素，而对应的Q实际上不一致分散式。这已在this 专着中观察到，秒。 4.6（讨论指的是 MATLAB，但我猜 MATLAB 和 scipy 都使用相同的 LAPACK 例程）。建议将qr 提供的矩阵Q 通过后乘以随机酉对角矩阵来修改。
下面我在上面的参考中重现了实验，绘制了qr提供的“直接”Q矩阵的特征值相位的经验分布（直方图），以及“修改”版本，其中它可以看出，修改后的版本确实具有均匀的特征值相位，正如从均匀分布的正交矩阵所预期的那样。
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');
【讨论】：
你怎么知道得到的矩阵均匀分布在流形上？并根据歧管上的什么度量？ @gg 感谢您的评论！实际上，qr 的“直接”应用并没有提供均匀分布的正交矩阵。请参阅编辑后的答案以获取解决方法。
很好的答案！参考文献中肯，非常感谢！
不需要随机角度，只需使用R对角线的符号即可：Q_modified = Q @ np.diag(np.sign(np.diag(R))【参考方案3】：
这是从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 的测试，https://***.com/a/38426572/901925
【讨论】：
【参考方案4】：创建任何形状 (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)。
【讨论】：
我认为您应该将rand 更改为randn 以获得矩阵的均匀分布。否则，我认为这是最好的答案。【参考方案5】：
如果您想要一个具有正交列向量的非方阵，您可以使用上述任何一种方法创建一个方阵并删除一些列。
【讨论】：
【参考方案6】：from scipy.stats import special_ortho_group
num_dim=3
x = special_ortho_group.rvs(num_dim)
Documentation
【讨论】：
【参考方案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中创建随机正交矩阵的主要内容，如果未能解决你的问题，请参考以下文章
在 Python、NumPy 和 R 中创建相同的随机数序列