randomOrthogonalMatrix¶

Generates a pseudorandom orthogonal matrix or a correlation matrix.

Synopsis¶

randomOrthogonalMatrix (n)

Required Arguments¶

int n (Input): The order of the matrix to be generated.

Return Value¶

n by n random orthogonal matrix.

Optional Arguments¶

eigenvalues, float (Input): A vector of length n containing the eigenvalues of the correlation matrix to be generated. The elements of eigenvalues must be positive, they must sum to n, and they cannot all be equal.
aMatrix, float (Input): n by n random orthogonal matrix. A random correlation matrix is generated using the orthogonal matrix input in aMatrix. The option eigenvalues must also be supplied if aMatrix is used.

Description¶

Function randomOrthogonalMatrix generates a pseudorandom orthogonal matrix from the invariant Haar measure. For each column, a random vector from a uniform distribution on a hypersphere is selected and then is projected onto the orthogonal complement of the columns already formed. The method is described by Heiberger (1978). (See also Tanner and Thisted 1982.)

If the optional argument eigenvalues is used, a correlation matrix is formed by applying a sequence of planar rotations to the matrix \(A^TDA\), where D = diag(eigenvalues[0], …, eigenvalues [n-1]), so as to yield ones along the diagonal. The planar rotations are applied in such an order that in the two by two matrix that determines the rotation, one diagonal element is less than 1.0 and one is greater than 1.0. This method is discussed by Bendel and Mickey (1978) and by Lin and Bendel (1985).

The distribution of the correlation matrices produced by this method is not known. Bendel and Mickey (1978) and Johnson and Welch (1980) discuss the distribution.

For larger matrices, rounding can become severe.

Example¶

In this example, randomOrthogonalMatrix is used to generate a 4 by 4 pseudorandom correlation matrix with eigenvalues in the ratio 1:2:3:4.

from numpy import *
from pyimsl.stat.randomOrthogonalMatrix import randomOrthogonalMatrix
from pyimsl.stat.randomSeedSet import randomSeedSet
from pyimsl.stat.writeMatrix import writeMatrix

n_random = 4
ev = [1., 2., 3., 4.]

for i in range(0, n_random):
    ev[i] = 4. * ev[i] / 10.

randomSeedSet(123457)

a = randomOrthogonalMatrix(n_random)

writeMatrix("Random orthogonal matrix", a)

cor = randomOrthogonalMatrix(n_random,
                             eigenvalues=ev, aMatrix=a)

writeMatrix("Random correlation matrix", cor)

Output¶

 
              Random orthogonal matrix
             1            2            3            4
1      -0.8804      -0.2417       0.4065      -0.0351
2       0.3088      -0.3002       0.5520       0.7141
3      -0.3500       0.5256      -0.3874       0.6717
4      -0.0841      -0.7584      -0.6165       0.1941
 
              Random correlation matrix
             1            2            3            4
1        1.000       -0.236       -0.326       -0.110
2       -0.236        1.000        0.191       -0.017
3       -0.326        0.191        1.000       -0.435
4       -0.110       -0.017       -0.435        1.000