How to efficiently and randomly produce an orthonormal matrix?

The following code can produce an orthonormal matrix for me:

import numpy as np
from scipy.linalg import orth

def get_matrix(M, N):
    Phi = np.random.randn(N, N).astype(np.float32)
    return orth(Phi)[:M]

Phi = get_matrix(10, 10)
print(np.matmul(Phi, np.transpose(Phi, [1,0])))  # should be very close to identity matrix

I want to use it in each network forwarding process, during training, but the above implementation with numpy and scipy may be too slow for me.

How to implement the above function with high efficiency by using PyTorch?

Hi Bin!

You can use pytorch’s svd() to do what I am pretty sure is the same
thing as scipy.linalg.orth() is doing:

>>> import torch
>>> torch.__version__
'1.11.0'
>>> _ = torch.manual_seed (2022)
>>> gaus = torch.randn (5, 5)
>>> svd = torch.linalg.svd (gaus)
>>> orth = svd[0] @ svd[2]
>>> orth
tensor([[-0.5558, -0.6665, -0.2910, -0.1637, -0.3679],
        [-0.1533,  0.2271, -0.8099,  0.4482,  0.2612],
        [ 0.5413, -0.3854,  0.0067,  0.6283, -0.4044],
        [ 0.1234, -0.5894,  0.0896,  0.0411,  0.7923],
        [ 0.5995, -0.0906, -0.5014, -0.6131, -0.0722]])
>>> orth @ orth.T
tensor([[ 1.0000e+00,  5.4137e-08,  1.2895e-07,  1.8028e-07,  7.5459e-08],
        [ 5.4137e-08,  1.0000e+00, -1.2912e-07,  1.9165e-07,  2.2009e-07],
        [ 1.2895e-07, -1.2912e-07,  1.0000e+00,  2.5904e-08, -3.9519e-08],
        [ 1.8028e-07,  1.9165e-07,  2.5904e-08,  1.0000e+00, -1.5845e-08],
        [ 7.5459e-08,  2.2009e-07, -3.9519e-08, -1.5845e-08,  1.0000e+00]])

I believe* (but it would be worth double checking) that the generated
matrix, orth, is randomly distributed over the space of orthonormal
matrices (with respect to the standard invariant measure).

If you don’t want a square matrix, you can drop some rows or columns,
as in your get_matrix() function.

*) See, for example, this answer on stackoverflow:

Best.

K. Frank

1 Like

Thank you very much. My question is addressed by your answer.