Codifique golf una matriz ortogonal aleatoria

9

Una matriz ortogonal es una matriz cuadrada con entradas reales cuyas columnas y filas son vectores unitarios ortogonales (es decir, vectores ortonormales).

Esto significa que M ^ TM = I, donde I es la matriz de identidad y ^ T significa la transposición de la matriz.

Tenga en cuenta que esto es ortogonal no "ortogonal especial", por lo que el determinante de M puede ser 1 o -1.

El objetivo de este desafío no es la precisión de la máquina, por lo que si M ^ TM = I dentro de 4 decimales estará bien.

La tarea es escribir código que tome un número entero positivo n > 1y genere una matriz ortogonal aleatoria n por n . La matriz debe elegirse aleatoria y uniformemente de todas las matrices n por n ortogonales. En este contexto, "uniforme" se define en términos de la medida de Haar, que esencialmente requiere que la distribución no cambie si se multiplica por cualquier matriz ortogonal elegida libremente. Esto significa que los valores de la matriz serán valores de coma flotante en el rango de -1 a 1.

La entrada y la salida pueden ser de cualquier forma que le resulte conveniente.

Muestre un ejemplo explícito de su código en ejecución.

No puede usar ninguna función de biblioteca existente que cree matrices ortogonales. Esta regla es un poco sutil, así que explicaré más. Esta regla prohíbe el uso de cualquier función existente que tome alguna (o ninguna) entrada y genere una matriz de tamaño al menos n por n que se garantice que sea ortogonal. Como ejemplo extremo, si desea la matriz de identidad n por n, tendrá que crearla usted mismo.

Puede usar cualquier biblioteca generadora de números aleatorios estándar para elegir números aleatorios de su elección.

Su código debe completarse en unos segundos como máximo n < 50.


fuente
Entonces, ¿está prohibido usar una matriz de identidad incorporada?
JungHwan Min
@JHM No puede usarlo para crear una matriz de identidad n por n al menos.
¿Qué hay de diag? Crea una matriz diagonal que es de hecho ortogonal pero no siempre ortonormal.
Karl Napf
Este parece ser un ejemplo de "hacer X sin Y", que, por lo tanto, el consenso debe evitarse.
falla
1
Las matrices diagonales no son matrices ortogonales, por lo que diagdeberían estar bien.
Angs

Respuestas:

7

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

Extienda recursivamente una matriz de tamaño ortogonal n-1agregando 1 a la esquina inferior derecha y aplique una reflexión aleatoria de Householder. randnda una matriz con valores aleatorios de una distribución gaussiana y z dda un vector unitario distribuido uniformemente en ddimensiones.

haussholder tau vdevuelve la matriz I - tau*v*vᵀque no es ortogonal cuando vno es un vector unitario.

Uso:

*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
Angs
fuente
Hacer la 1×1matriz toma demasiado espacio para mi gusto, un caso especial solo para obtener cero de una variable aleatoria gaussiana: / (Sin ella, existe una posibilidad infinitesimal de obtener una columna cero)
Angs
Me gusta tu espíritu de hacerlo totalmente correcto, pero creo que puedes dejar ese requisito. En mi código también existe la posibilidad de que 2 filas sean linealmente dependientes y a nadie le importe.
Karl Napf
@KarlNapf bueno, descubrí una forma de perder dos bytes de esa parte de todos modos, así que el problema se resolvió en parte :)
Angs
Ah, está bien, borrando mis comentarios ...
Karl Napf
¡Siempre feliz cuando una respuesta de Haskell gana!
4

Python 2 + NumPy, 163 bytes

Gracias a xnor por señalar el uso de valores aleatorios distribuidos normales en lugar de valores uniformes.

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

Utiliza la ortogonalización de Gram Schmidt en una matriz con valores aleatorios gaussianos para tener todas las direcciones.

El código de demostración es seguido por

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]]

Se completa en un abrir y cerrar por n = 50 y unos segundos por n = 500.

Karl Napf
fuente
No creo que esto sea uniforme. Su distribución inicial con un cubo, que tiene más cosas hacia las diagonales. Los gaussianos aleatorios funcionarían porque generan una distribución esféricamente simétrica.
xnor
@xnor arreglado. Afortunadamente, esto costó exactamente 1 byte.
Karl Napf
@xnor Aún más afortunado, esto guardó los bytes para-0.5
Karl Napf
Casi, necesita que la media de lo normal sea 0, pero eso no es más que n.
xnor
-1

Mathematica, 69 bytes, probablemente no compitiendo

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

QRDecompositiondevuelve un par de matrices, la primera de las cuales se garantiza que es ortogonal (y la segunda no es ortogonal, sino triangular superior). Se podría argumentar que esto técnicamente obedece la letra de la restricción en la publicación: no genera una matriz ortogonal, sino un par de matrices ...

Mathematica, 63 bytes, definitivamente no compite

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

Orthogonalizeestá inequívocamente prohibido por el OP. Aún así, Mathematica es genial, ¿eh?

Greg Martin
fuente
You may not use any existing library function which creates orthogonal **matrices**.
Karl Napf