6

What Is a Random Orthogonal Matrix?

 3 years ago
source link: https://nhigham.com/2020/04/22/what-is-a-random-orthogonal-matrix/
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.

What Is a Random Orthogonal Matrix?

Various explicit parametrized formulas are available for constructing orthogonal matrices. To construct a random orthogonal matrix we can take such a formula and assign random values to the parameters. For example, a Householder matrix H = I - 2uu^T/(u^Tu) is orthogonal and symmetric and we can choose the nonzero vector u randomly. Such an example is rather special, though, as it is a rank-1 perturbation of the identity matrix.

What is usually meant by a random orthogonal matrix is a matrix distributed according to the Haar measure over the group of orthogonal matrices. The Haar measure provides a uniform distribution over the orthogonal matrices. Indeed it is invariant under multiplication on the left and the right by orthogonal matrices: if Q is from the Haar distribution then so is UQV for any orthogonal (possibly non-random) U and V. A random Householder matrix is not Haar distributed.

A matrix from the Haar distribution can be generated as the orthogonal factor in the QR factorization of a random matrix with elements from the standard normal distribution (mean 0, variance 1). In MATLAB this is done by the code

[Q,R] = qr(randn(n));
Q = Q*diag(sign(diag(R)));

The statement [Q,R] = qr(randn(n)), which returns the orthogonal factor Q, is not enough on its own to give a Haar distributed matrix, because the QR factorization is not unique. The second line adjusts the signs so that Q is from the unique factorization in which the triangular factor R has nonnegative diagonal elements. This construction requires 2n^3 flops.

A more efficient construction is possible, as suggested by Stewart (1980). Let x_k be an (n-k+1)-vector of elements from the standard normal distribution and let H_k be the Householder matrix that reduces x_k to r_{kk}e_1, where e_1 is the first unit vector. Then Q = DH'_1H'_2\dots H'_{n-1} is Haar distributed, where H'_k = \mathrm{diag}(I_{k-1}, H_k), D = \mathrm{diag}(\mathrm{sign}(r_{kk})), and r_{nn} = x_n. This construction expresses Q as the product of n-1 Householder matrices of growing effective dimension, and the product can be formed from right to left in 4n^3/3 flops. The MATLAB statement Q = gallery('qmult',n) carries out this construction.

A similar construction can be made using Givens rotations (Anderson et al., 1987).

Orthogonal matrices from the Haar distribution can also be formed as Q = A (A^TA)^{-1/2}, where the elements of A are from the standard normal distribution. This Q is the orthogonal factor in the polar decomposition A = QH (where H is symmetric positive semidefinite).

Random orthogonal matrices arise in a variety of applications, including Monte Carlo simulation, random matrix theory, machine learning, and the construction of test matrices with known eigenvalues or singular values.

All these ideas extend to random unitary (complex) matrices. In MATLAB, Haar distributed unitary matrices can be constructed by the code

[Q,R] = qr(complex(randn(n),randn(n)));
Q = Q*diag(sign(diag(R)));

This code exploits the fact that the R factor computed by MATLAB has real diagonal entries. If the diagonal of R were complex then this code would need to be modified to use the complex sign function given by \mathrm{sign}(z) = z/|z|.

References

This is a minimal set of references, which contain further useful references within.

Related Blog Posts


About Joyk


Aggregate valuable and interesting links.
Joyk means Joy of geeK