Код гольф случайной ортогональной матрицы

9

Ортогональная матрица представляет собой квадратную матрицу с вещественными элементами , чьи столбцы и строки являются ортогональные орты (т.е. ортонормальные векторы).

Это означает, что M ^ TM = I, где I - единичная матрица, а ^ T означает транспонирование матрицы.

Обратите внимание, что это ортогональный, а не «специальный ортогональный», поэтому определитель М может быть 1 или -1.

Целью этой задачи является не точность машины, поэтому, если M ^ TM = I с точностью до 4 десятичных знаков, это будет хорошо.

Задача состоит в том, чтобы написать код, который принимает положительное целое число n > 1и выводит случайную ортогональную матрицу n на n . Матрица должна выбираться случайным образом и равномерно из всех n по n ортогональных матриц. В этом контексте «униформа» определяется в терминах меры Хаара, которая, по сути, требует, чтобы распределение не изменялось при умножении на произвольно выбранную ортогональную матрицу. Это означает, что значения матрицы будут значениями с плавающей запятой в диапазоне от -1 до 1.

Вход и выход могут быть любой формы, которая вам удобна.

Пожалуйста, покажите явный пример выполнения вашего кода.

Вы не можете использовать любую существующую библиотечную функцию, которая создает ортогональные матрицы. Это правило немного тонкое, поэтому я объясню больше. Это правило запрещает использование любой существующей функции, которая принимает некоторые (или не принимает) входные данные и выводит матрицу размером не менее n на n, которая гарантированно будет ортогональной. В качестве крайнего примера, если вы хотите, чтобы n - n единичная матрица, вам придется создать ее самостоятельно.

Вы можете использовать любую стандартную библиотеку генератора случайных чисел для выбора случайных чисел по вашему выбору.

Ваш код должен завершиться в течение максимум нескольких секунд n < 50.


источник
Так что использование встроенной идентификационной матрицы запрещено?
JungHwan Мин
@JHM Вы не можете использовать его для создания по меньшей мере n-n-идентичной матрицы.
Как насчет diag? Это создает диагональную матрицу, которая действительно ортогональна, но не всегда ортонормирована.
Карл Напф
Похоже, это пример «делай X без Y», чего следует избегать.
flawr
1
Диагональные матрицы не являются ортогональными матрицами, поэтому diagдолжно быть в порядке.
Angs

Ответы:

7

Haskell, 169 150 148 141 132 131 байт

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

Рекурсивно расширяйте ортогональную матрицу размера n-1, добавляя 1 в нижний правый угол и применяя случайное отражение Домохозяина. randnдает матрицу со случайными значениями из гауссовского распределения и z dдает равномерно распределенный единичный вектор по dизмерениям.

haussholder tau vвозвращает матрицу, I - tau*v*vᵀкоторая не ортогональна, когда vне является единичным вектором.

Применение:

*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
источник
Создание 1×1матрицы занимает слишком много места на мой вкус, особый случай только для получения нуля из гауссовой случайной величины: / (Без этого есть бесконечно малый шанс получить нулевой столбец)
Angs
Мне нравится твой дух сделать это полностью правильным, но я думаю, что ты можешь отбросить это требование. В моем коде также есть вероятность того, что 2 строки линейно зависимы и никому нет до этого дела.
Карл Напф
@KarlNapf хорошо, я нашел способ потерять два байта из этой части, так что проблема частично решена :)
Angs
Ах, хорошо, удаляя мои комментарии ...
Карл Напф
Всегда счастлив, когда победит ответ на Haskell!
4

Python 2 + NumPy, 163 байта

Спасибо xnor за указание использовать обычные распределенные случайные значения вместо равномерных.

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

Использует ортогонализацию Грама Шмидта на матрице с гауссовыми случайными значениями, чтобы иметь все направления.

Код демонстрации сопровождается

print dot(Q.transpose(),Q)

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

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

Это завершается за мгновение для n = 50 и несколько секунд для n = 500.

Карл Напф
источник
Я не думаю, что это единообразно. Ваш начальный дистрибутив с кубом, который имеет больше материала по диагонали. Случайные гауссианы будут работать, потому что они генерируют сферически симметричное распределение.
xnor
@xnor исправлено. К счастью, это стоило ровно 1 байта.
Карл Напф
@xnor Еще более удачно, что это спасло байты для-0.5
Карл Напф
Почти необходимо, чтобы среднее значение нормали было равно 0, но это не больше, чем n.
xnor
-1

Mathematica, 69 байт, вероятно, не конкурирует

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

QRDecompositionвозвращает пару матриц, первая из которых гарантированно будет ортогональной (а вторая из которых не ортогональна, а верхняя треугольная). Можно утверждать, что это технически подчиняется букве ограничения в посте: он не выводит ортогональную матрицу, но пару матриц ....

Mathematica, 63 байта, определенно неконкурентный

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

Orthogonalizeоднозначно запрещено ОП. Тем не менее, Mathematica довольно круто, а?

Грег Мартин
источник
You may not use any existing library function which creates orthogonal **matrices**.
Карл Напф