12
\$\begingroup\$

An orthogonal matrix is a square matrix with real entries whose columns and rows are orthogonal unit vectors (i.e., orthonormal vectors).

This means that M^T M = I, where I is the identity matrix and ^T signifies matrix transposition.

Note that this is orthogonal not "special orthogonal" so the determinant of M can be 1 or -1.

The aim of this challenge is not machine precision so if M^T M = I to within 4 decimal places that will be fine.

The task is to write code that takes a positive integer n > 1 and outputs a random orthogonal n by n matrix. The matrix should be randomly and uniformly chosen from all n by n orthogonal matrices. In this context, "uniform" is defined in terms of the Haar measure, which essentially requires that the distribution does not change if multiplied by any freely chosen orthogonal matrix. This means the values of the matrix will be floating point values in the range -1 to 1.

The input and output can be any form you find convenient.

Please show an explicit example of your code running.

You may not use any existing library function which creates orthogonal matrices. This rule is a little subtle so I will explain more. This rule bans the use of any existing function which takes in some (or no) input and outputs a matrix of size at least n by n which is guaranteed to be orthogonal. As an extreme example, if you want the n by n identity matrix, you will have to create it yourself.

You can use any standard random number generator library for choosing random numbers of your choosing.

Your code should complete within at most a few seconds for n < 50.

\$\endgroup\$
5
  • \$\begingroup\$ So using built-in identity matrix is forbidden? \$\endgroup\$ Nov 10, 2016 at 0:33
  • \$\begingroup\$ @JHM You can't use it to create an n by n identity matrix at least. \$\endgroup\$
    – user9206
    Nov 10, 2016 at 5:55
  • \$\begingroup\$ What about diag? It creates an diagonal matrix which is indeed orthogonal but not always orthonormal. \$\endgroup\$
    – Karl Napf
    Nov 10, 2016 at 7:26
  • \$\begingroup\$ This seems to be an example of "do X without Y", which - so the consensus - should be avoided. \$\endgroup\$
    – flawr
    Nov 10, 2016 at 23:01
  • 1
    \$\begingroup\$ Diagonal matrices aren't orthogonal matrices so diag should be ok. \$\endgroup\$
    – Angs
    Nov 11, 2016 at 5:43

4 Answers 4

8
\$\begingroup\$

Haskell, 169 150 148 141 132 131 bytes

import Numeric.LinearAlgebra
z=(unitary.flatten<$>).randn 1
r 1=asRow<$>z 1
r n=do;m<-r$n-1;(<>diagBlock[m,1]).haussholder 2<$>z n

Recursively extend an orthogonal matrix of size n-1 by adding 1 to lower right corner and apply a random Householder reflection. randn gives a matrix with random values from a gaussian distribution and z d gives a uniformly distributed unit vector in d dimensions.

haussholder tau v returns the matrix I - tau*v*vᵀ which isn't orthogonal when v isn't a unit vector.

Usage:

*Main> m <- r 5
*Main> disp 5 m
5x5
-0.24045  -0.17761   0.01603  -0.83299  -0.46531
-0.94274   0.12031   0.00566   0.29741  -0.09098
-0.02069   0.30417  -0.93612  -0.13759   0.10865
 0.02155  -0.83065  -0.35109   0.32365  -0.28556
-0.22919  -0.41411   0.01141  -0.30659   0.82575
*Main> (<1e-14) . maxElement . abs $ tr m <> m - ident 5
True
\$\endgroup\$
5
  • \$\begingroup\$ Making the 1×1 matrix takes too much space for my taste, a special case just for getting zero from a gaussian random variable :/ (Without it, there's a infinitesimal chance to get a zero column) \$\endgroup\$
    – Angs
    Nov 10, 2016 at 11:40
  • \$\begingroup\$ I like your spirit of making it totally correct, but i think you can drop that requirement. In my code there is also a chance that 2 rows are linearly dependent and no one cares. \$\endgroup\$
    – Karl Napf
    Nov 10, 2016 at 15:22
  • \$\begingroup\$ @KarlNapf well, I figured out a way to lose two bytes from that part anyway, so problem partly solved :) \$\endgroup\$
    – Angs
    Nov 10, 2016 at 15:35
  • \$\begingroup\$ Ah okay, deleting my comments... \$\endgroup\$
    – Karl Napf
    Nov 11, 2016 at 11:11
  • \$\begingroup\$ Always happy when a Haskell answer wins! \$\endgroup\$
    – user9206
    Nov 16, 2016 at 14:23
5
\$\begingroup\$

Python 2+NumPy, 163 bytes

Thanks to xnor for pointing out to use normal distributed random values instead of uniform ones.

from numpy import*
n=input()
Q=random.randn(n,n)
for i in range(n):
 for j in range(i):u=Q[:,j];Q[:,i]-=u*dot(u,Q[:,i])/dot(u,u)
Q/=(Q**2).sum(axis=0)**0.5
print Q

Uses the Gram Schmidt Orthogonalization on a matrix with gaussian random values to have all directions.

The demonstration code is followed by

print dot(Q.transpose(),Q)

n=3:

[[-0.2555327   0.89398324  0.36809917]
 [-0.55727299  0.17492767 -0.81169398]
 [ 0.79003155  0.41254608 -0.45349298]]
[[  1.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  0.00000000e+00   1.00000000e+00  -5.55111512e-17]
 [  0.00000000e+00  -5.55111512e-17   1.00000000e+00]]

n=5:

[[-0.63470728  0.41984536  0.41569193  0.25708079  0.42659843]
 [-0.36418389  0.06244462 -0.82734663 -0.24066123  0.3479231 ]
 [ 0.07863783  0.7048799   0.08914089 -0.64230492 -0.27651168]
 [ 0.67691426  0.33798442 -0.05984083  0.17555011  0.62702062]
 [-0.01095148 -0.45688226  0.36217501 -0.65773717  0.47681205]]
[[  1.00000000e+00   1.73472348e-16   5.37764278e-17   4.68375339e-17
   -2.23779328e-16]
 [  1.73472348e-16   1.00000000e+00   1.38777878e-16   3.33066907e-16
   -6.38378239e-16]
 [  5.37764278e-17   1.38777878e-16   1.00000000e+00   1.38777878e-16
    1.11022302e-16]
 [  4.68375339e-17   3.33066907e-16   1.38777878e-16   1.00000000e+00
    5.55111512e-16]
 [ -2.23779328e-16  -6.38378239e-16   1.11022302e-16   5.55111512e-16
    1.00000000e+00]]

It completes within a blink for n=50 and a few seconds for n=500.

\$\endgroup\$
4
  • \$\begingroup\$ I don't think this is uniform. Your starting distribution with a cube, which has more stuff towards the diagonals. Random gaussians would work because they generate a spherically symmetric distribution. \$\endgroup\$
    – xnor
    Nov 10, 2016 at 7:57
  • \$\begingroup\$ @xnor fixed. Luckily, this costed exactly 1 byte. \$\endgroup\$
    – Karl Napf
    Nov 10, 2016 at 8:07
  • \$\begingroup\$ @xnor Even more lucky, this saved the bytes for -0.5 \$\endgroup\$
    – Karl Napf
    Nov 10, 2016 at 8:10
  • \$\begingroup\$ Almost, you need the mean of the normal to be 0, but that's no longer than n. \$\endgroup\$
    – xnor
    Nov 10, 2016 at 18:08
4
+200
\$\begingroup\$

APL (Dyalog Unicode), 63 61 bytes

⎕CY'dfns'
n←⊢÷.5*⍨a←+.×⍨
{{⍵,(n⊢-⊣×a)/⌽⍺,⍵}/⊂∘n¨↓NormRand⍵ ⍵}

Try it online for n=5

-1 byte thanks to @user41805 (improved train for inner dfn)

Uses the modified Gram-Schmidt process on a set of n random unit vectors --- fast.

How?

⎕CY'dfns'   ⍝ Copy the dfns namespace, which includes the NormRand dfn
n←⊢÷.5*⍨+.×⍨  ⍝ n: normalize the vector argument, from https://codegolf.stackexchange.com/a/149160/68261
{{⍵,(n⊢-⊣×+.×⍨)/⌽⍺,⍵}/⊂∘n¨↓NormRand⍵ ⍵} ⍝ Main anonymous function
                           NormRand ⍵ ⍵  ⍝ Generate an n×n array of gaussian variables
                       ⊂∘n¨↓             ⍝ Treat as an array of singular arrays of n-vectors; normalize each
 {                   }/                  ⍝ Reduce in a sneaky manner to get array of orthogonal (e_1, ... e_n)
                                         ⍝ Start with e_arr=⍬ (empty)
                                         ⍝ For each vector v_i:
                ⌽⍺,⍵                      ⍝ Start this step with e_arr←(e_{i-1}, ..., e_1, e_0, v_i)
    (         )/                          ⍝ Reduce this to repeatedly modify from the right:
     n⊢-⊣×+.×⍨                             ⍝ Normalize(v_i - proj_{e_j}(v_i))
  ⍵,                                       ⍝ Append this to e_arr
\$\endgroup\$
2
  • \$\begingroup\$ You should be able to trainify the innermost dfn like spoiler \$\endgroup\$
    – user41805
    Jul 12, 2020 at 9:31
  • \$\begingroup\$ @user41805 Nice! That saves 1 byte over my train and lets me save another since it repeats +.×⍨ \$\endgroup\$ Jul 13, 2020 at 7:29
-1
\$\begingroup\$

Mathematica, 69 bytes, probably non-competing

#&@@QRDecomposition@Array[RandomVariate@NormalDistribution[]&,{#,#}]&

QRDecomposition returns a pair of matrices, the first of which is guaranteed to be orthogonal (and the second of which is not orthogonal, but upper triangular). One could argue that this technically obeys the letter of the restriction in the post: it doesn't output an orthogonal matrix, but a pair of matrices....

Mathematica, 63 bytes, definitely non-competing

Orthogonalize@Array[RandomVariate@NormalDistribution[]&,{#,#}]&

Orthogonalize is unambiguously forbidden by the OP. Still, Mathematica is pretty cool eh?

\$\endgroup\$
2
  • \$\begingroup\$ You may not use any existing library function which creates orthogonal **matrices**. \$\endgroup\$
    – Karl Napf
    Nov 10, 2016 at 8:38
  • \$\begingroup\$ Uses library functions which generate orthoganal matrices. \$\endgroup\$
    – emanresu A
    Jun 27, 2022 at 4:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.