How to create random orthonormal matrix in python numpy
Version 0.18 of scipy has 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
For example,
In [24]: from scipy.stats import ortho_group # Requires version 0.18 of scipyIn [25]: m = ortho_group.rvs(dim=3)In [26]: mOut[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.]])
You can obtain a random n x n orthogonal matrix Q, (uniformly distributed over the manifold of n x n orthogonal matrices) by performing a QR factorization of an n x n matrix with elements i.i.d. Gaussian random variables of mean 0 and variance 1. Here is an example:
import numpy as npfrom scipy.linalg import qrn = 3H = 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]]EDIT: (Revisiting this answer after the comment by @g g.) The claim above on the QR decomposition of a Gaussian matrix providing a uniformly distributed (over the, so called, Stiefel manifold) orthogonal matrix is suggested by Theorems 2.3.18-19 of this reference. Note that the statement of the result suggests a "QR-like" decomposition, however, with the triangular matrix R having positive elements.
Apparently, the qr function of scipy (numpy) function does not guarantee positive diagonal elements for R and the corresponding Q is actually not uniformly distributed. This has been observed in this monograph, Sec. 4.6 (the discussion refers to MATLAB, but I guess both MATLAB and scipy use the same LAPACK routines). It is suggested there that the matrix Q provided by qr is modified by post multiplying it with a random unitary diagonal matrix.
Below I reproduce the experiment in the above reference, plotting the empirical distribution (histogram) of phases of eigenvalues of the "direct" Q matrix provided by qr, as well as the "modified" version, where it is seen that the modified version does indeed have a uniform eigenvalue phase, as would be expected from a uniformly distributed orthogonal matrix.
from scipy.linalg import qr, eigvalsfrom seaborn import distplotn = 50repeats = 10000angles = []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');
This is the rvs method pulled from the https://github.com/scipy/scipy/pull/5622/files, with minimal change - just enough to run as a stand alone numpy function.
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 HIt matches Warren's test, https://stackoverflow.com/a/38426572/901925