Code golf une matrice orthogonale aléatoire

9

Une matrice orthogonale est une matrice carrée avec des entrées réelles dont les colonnes et les lignes sont des vecteurs unitaires orthogonaux (c.-à-d. Des vecteurs orthonormaux).

Cela signifie que M ^ TM = I, où I est la matrice d'identité et ^ T signifie la transposition matricielle.

Notez que ce n'est pas orthogonal "spécial orthogonal", donc le déterminant de M peut être 1 ou -1.

Le but de ce défi n'est pas la précision de la machine, donc si M ^ TM = I à 4 décimales près, ça ira.

La tâche consiste à écrire du code qui prend un entier positif n > 1et génère une matrice orthogonale aléatoire n par n . La matrice doit être choisie de manière aléatoire et uniforme parmi toutes les matrices orthogonales n par n. Dans ce contexte, "uniforme" est défini en termes de la mesure de Haar, qui exige essentiellement que la distribution ne change pas si elle est multipliée par une matrice orthogonale librement choisie. Cela signifie que les valeurs de la matrice seront des valeurs à virgule flottante comprises entre -1 et 1.

L'entrée et la sortie peuvent prendre n'importe quelle forme que vous jugez pratique.

Veuillez montrer un exemple explicite de l'exécution de votre code.

Vous ne pouvez utiliser aucune fonction de bibliothèque existante qui crée des matrices orthogonales. Cette règle est un peu subtile donc je vais vous expliquer plus. Cette règle interdit l'utilisation de toute fonction existante qui accepte certaines entrées (ou aucune) et génère une matrice de taille au moins n par n qui est garantie d'être orthogonale. À titre d'exemple extrême, si vous voulez la matrice d'identité n par n, vous devrez la créer vous-même.

Vous pouvez utiliser n'importe quelle bibliothèque de générateur de nombres aléatoires standard pour choisir des nombres aléatoires de votre choix.

Votre code devrait se terminer en quelques secondes au maximum n < 50.


la source
Il est donc interdit d'utiliser la matrice d'identité intégrée?
JungHwan Min
@JHM Vous ne pouvez pas l'utiliser pour créer au moins une matrice d'identité n par n.
Et alors diag? Il crée une matrice diagonale qui est en effet orthogonale mais pas toujours orthonormale.
Karl Napf
Cela semble être un exemple de "faire X sans Y", qui - donc le consensus - devrait être évité.
flawr
1
Les matrices diagonales ne sont pas des matrices orthogonales, donc ça diagdevrait aller.
Angs

Réponses:

7

Haskell, 169 150 148 141 132 131 131 octets

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

Étendez récursivement une matrice orthogonale de taille n-1en ajoutant 1 au coin inférieur droit et appliquez une réflexion aléatoire du chef de ménage. randndonne une matrice avec des valeurs aléatoires d'une distribution gaussienne et z ddonne un vecteur unitaire uniformément distribué en ddimensions.

haussholder tau vrenvoie la matrice I - tau*v*vᵀqui n'est pas orthogonale lorsqu'elle vn'est pas un vecteur unitaire.

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
Angs
la source
Faire la 1×1matrice prend trop de place à mon goût, un cas spécial juste pour obtenir zéro à partir d'une variable aléatoire gaussienne: / (Sans cela, il y a une chance infinitésimale d'obtenir une colonne zéro)
Angs
J'aime votre esprit de le rendre totalement correct, mais je pense que vous pouvez laisser tomber cette exigence. Dans mon code, il y a aussi une chance que 2 lignes soient linéairement dépendantes et que personne ne s'en soucie.
Karl Napf
@KarlNapf bien, j'ai quand même trouvé un moyen de perdre deux octets de cette partie, donc le problème est en partie résolu :)
Angs
Ah d'accord, en supprimant mes commentaires ...
Karl Napf
Toujours heureux quand une réponse Haskell gagne!
4

Python 2 + NumPy, 163 octets

Merci à xnor d'avoir indiqué d'utiliser des valeurs aléatoires distribuées normales au lieu de valeurs 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

Utilise l' orthogonalisation de Gram Schmidt sur une matrice avec des valeurs aléatoires gaussiennes pour avoir toutes les directions.

Le code de démonstration est suivi de

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

Il se termine en un clin d'œil pour n = 50 et quelques secondes pour n = 500.

Karl Napf
la source
Je ne pense pas que ce soit uniforme. Votre distribution de départ avec un cube, qui a plus de choses vers les diagonales. Les gaussiens aléatoires fonctionneraient car ils génèrent une distribution à symétrie sphérique.
xnor
@xnor fixe. Heureusement, cela a coûté exactement 1 octet.
Karl Napf
@xnor Encore plus de chance, cela a sauvé les octets de-0.5
Karl Napf
Presque, vous avez besoin que la moyenne de la normale soit 0, mais ce n'est pas plus long que n.
xnor
-1

Mathematica, 69 octets, probablement non concurrentiel

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

QRDecompositionrenvoie une paire de matrices, dont la première est garantie orthogonale (et la seconde n'est pas orthogonale, mais triangulaire supérieure). On pourrait soutenir que cela obéit techniquement à la lettre de la restriction dans le message: il ne produit pas une matrice orthogonale, mais une paire de matrices ....

Mathematica, 63 octets, définitivement non concurrentiel

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

Orthogonalizeest sans ambiguïté interdit par le PO. Pourtant, Mathematica est plutôt cool hein?

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