El desafío es escribir el código más rápido posible para calcular el hafniano de una matriz .
El hafniano de una matriz simétrica 2n
-por- se define como:2n
A
Aquí S 2n representa el conjunto de todas las permutaciones de los enteros de 1
a 2n
, es decir [1, 2n]
.
El enlace de wikipedia también ofrece una fórmula de aspecto diferente que puede ser de interés (e incluso existen métodos más rápidos si busca más en la web). La misma página wiki habla de matrices de adyacencia, pero su código también debería funcionar para otras matrices. Puede suponer que todos los valores serán enteros, pero no que sean todos positivos.
También hay un algoritmo más rápido, pero parece difícil de entender. y Christian Sievers fue el primero en implementarlo (en Haskell).
En esta pregunta, todas las matrices son cuadradas y simétricas con una dimensión par.
Implementación de referencia (tenga en cuenta que esto está utilizando el método más lento posible).
Aquí hay un ejemplo de código de Python del Sr. Xcoder.
from itertools import permutations
from math import factorial
def hafnian(matrix):
my_sum = 0
n = len(matrix) // 2
for sigma in permutations(range(n*2)):
prod = 1
for j in range(n):
prod *= matrix[sigma[2*j]][sigma[2*j+1]]
my_sum += prod
return my_sum / (factorial(n) * 2 ** n)
print(hafnian([[-1, 1, 1, -1, 0, 0, 1, -1], [1, 0, 1, 0, -1, 0, -1, -1], [1, 1, -1, 1, -1, -1, 0, -1], [-1, 0, 1, -1, -1, 1, -1, 0], [0, -1, -1, -1, -1, 0, 0, -1], [0, 0, -1, 1, 0, 0, 1, 1], [1, -1, 0, -1, 0, 1, 1, 0], [-1, -1, -1, 0, -1, 1, 0, 1]]))
4
M = [[1, 1, 0, 0, 0, 0, 0, 1, 0, 0], [1, 1, -1, 0, -1, 1, 1, 1, 0, -1], [0, -1, -1, -1, 0, -1, -1, 0, -1, 1], [0, 0, -1, 1, -1, 1, -1, 0, 1, -1], [0, -1, 0, -1, -1, -1, -1, 1, -1, 1], [0, 1, -1, 1, -1, 1, -1, -1, 1, -1], [0, 1, -1, -1, -1, -1, 1, 0, 0, 0], [1, 1, 0, 0, 1, -1, 0, 1, 1, -1], [0, 0, -1, 1, -1, 1, 0, 1, 1, 1], [0, -1, 1, -1, 1, -1, 0, -1, 1, 1]]
print(hafnian(M))
-13
M = [[-1, 0, -1, -1, 0, -1, 0, 1, -1, 0, 0, 0], [0, 0, 0, 0, 0, -1, 0, 1, -1, -1, -1, -1], [-1, 0, 0, 1, 0, 0, 0, 1, -1, 1, -1, 0], [-1, 0, 1, -1, 1, -1, -1, -1, 0, -1, -1, -1], [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, -1, 0], [-1, -1, 0, -1, 0, 0, 1, 1, 1, 1, 1, 0], [0, 0, 0, -1, 0, 1, 1, -1, -1, 0, 1, 0], [1, 1, 1, -1, 0, 1, -1, 1, -1, -1, -1, -1], [-1, -1, -1, 0, 0, 1, -1, -1, -1, 1, -1, 0], [0, -1, 1, -1, 1, 1, 0, -1, 1, -1, 1, 1], [0, -1, -1, -1, -1, 1, 1, -1, -1, 1, 0, -1], [0, -1, 0, -1, 0, 0, 0, -1, 0, 1, -1, 1]]
print(hafnian(M))
13
M = [[-1, 1, 0, 1, 0, -1, 0, 0, -1, 1, -1, 1, 0, -1], [1, -1, 1, -1, 1, 1, -1, 0, -1, 1, 1, 0, 0, -1], [0, 1, 1, 1, -1, 1, -1, -1, 0, 0, -1, 0, -1, -1], [1, -1, 1, -1, 1, 0, 1, 1, -1, -1, 0, 0, 1, 1], [0, 1, -1, 1, 0, 1, 0, 1, -1, -1, 1, 1, 0, -1], [-1, 1, 1, 0, 1, 1, -1, 0, 1, -1, -1, -1, 1, -1], [0, -1, -1, 1, 0, -1, -1, -1, 0, 1, -1, 0, 1, -1], [0, 0, -1, 1, 1, 0, -1, 0, 0, -1, 0, 0, 0, 1], [-1, -1, 0, -1, -1, 1, 0, 0, 1, 1, 0, 1, -1, 0], [1, 1, 0, -1, -1, -1, 1, -1, 1, 1, 1, 0, 1, 0], [-1, 1, -1, 0, 1, -1, -1, 0, 0, 1, -1, 0, -1, 0], [1, 0, 0, 0, 1, -1, 0, 0, 1, 0, 0, 1, 1, 1], [0, 0, -1, 1, 0, 1, 1, 0, -1, 1, -1, 1, 1, -1], [-1, -1, -1, 1, -1, -1, -1, 1, 0, 0, 0, 1, -1, -1]]
print(hafnian(M))
83
La tarea
Debería escribir un código que, dado 2n
por una 2n
matriz, genere su hafniano.
Como tendré que probar su código, sería útil si pudiera darme una forma simple de dar una matriz como entrada a su código, por ejemplo, leyendo desde el estándar. Probaré su código en matrices elegidas al azar con elementos seleccionado de {-1, 0, 1}. El propósito de pruebas como esta es reducir la posibilidad de que el Hafnian sea un valor muy grande.
Idealmente, su código podrá leerse en matrices exactamente como las tengo en los ejemplos de esta pregunta directamente desde el estándar. Así es como se vería la entrada, [[1,-1],[-1,-1]]
por ejemplo. Si desea utilizar otro formato de entrada, solicítelo y haré todo lo posible para adaptarlo.
Puntuaciones y lazos
Probaré su código en matrices aleatorias de tamaño creciente y me detendré la primera vez que su código tarde más de 1 minuto en mi computadora. Las matrices de puntuación serán consistentes para todas las presentaciones a fin de garantizar la equidad.
Si dos personas obtienen el mismo puntaje, entonces el ganador es el que tiene el valor más rápido n
. Si están dentro de 1 segundo el uno del otro, entonces es el publicado primero.
Idiomas y bibliotecas
Puede usar cualquier idioma y bibliotecas disponibles que desee, pero ninguna función preexistente para calcular el hafniano. Siempre que sea posible, sería bueno poder ejecutar su código, así que incluya una explicación completa sobre cómo ejecutar / compilar su código en Linux, si es posible.
Mi máquina Los tiempos se ejecutarán en mi máquina de 64 bits. Esta es una instalación estándar de ubuntu con 8 GB de RAM, procesador AMD FX-8350 de ocho núcleos y Radeon HD 4250. Esto también significa que necesito poder ejecutar su código.
Solicite respuestas en más idiomas.
Sería genial obtener respuestas en su lenguaje de programación súper rápido favorito. Para empezar, ¿qué tal fortran , nim y óxido ?
Tabla de clasificación
- 52 por millas usando C ++ . 30 segundos.
- 50 por ngn usando C . 50 segundos
- 46 por Christian Sievers usando Haskell . 40 segundos
- 40 por millas usando Python 2 + pypy . 41 segundos
- 34 por ngn usando Python 3 + pypy . 29 segundos
- 28 por Dennis usando Python 3 . 35 segundos (Pypy es más lento)
Respuestas:
Haskell
Esto implementa una variación del Algoritmo 2 de Andreas Björklund: Contando combinaciones perfectas tan rápido como Ryser .
Compile usando
ghc
opciones de tiempo de compilación-O3 -threaded
y use opciones de tiempo de ejecución+RTS -N
para la paralelización. Toma entrada de stdin.fuente
parallel
yvector
debe ser instalado?Python 3
Pruébalo en línea!
fuente
C ++ (gcc)
Pruébalo en línea! (13s para n = 24)
Basado en la implementación más rápida de Python en mi otra publicación. Edite la segunda línea
#define S 3
en su máquina de 8 núcleos y compile cong++ -pthread -march=native -O2 -ftree-vectorize
.Divide el trabajo por la mitad, por lo que el valor de
S
debería serlog2(#threads)
. Los tipos se pueden cambiar fácilmente entreint
,long
,float
ydouble
modificando el valor de#define TYPE
.fuente
tr -d \ < matrix52.txt > matrix52s.txt
Python 3
Esto calcula haf (A) como una suma memorizada (A [i] [j] * haf (A sin filas y cols i y j)).
fuente
C
Otra implicación del artículo de Andreas Björklund , que es mucho más fácil de entender si también observa el código Haskell de Christian Sievers . Para los primeros niveles de la recursividad, distribuye subprocesos round-robin sobre las CPU disponibles. El último nivel de la recursión, que representa la mitad de las invocaciones, se optimiza a mano.
Compilar con:
gcc -O3 -pthread -march=native
; gracias @Dennis por una aceleración 2xn = 24 en 24s en TIO
Algoritmo:
La matriz, que es simétrica, se almacena en forma triangular inferior izquierda. Los índices de triángulo
i,j
corresponden al índice linealT(max(i,j))+min(i,j)
dondeT
es una macro parai*(i-1)/2
. Los elementos de la matriz son polinomios de gradon
. Un polinomio se representa como una matriz de coeficientes ordenados desde el término constante (p[0]
) hasta el coeficiente de x n (p[n]
). Los valores iniciales de la matriz -1,0,1 se convierten primero en polinomios constantes.Realizamos un paso recursivo con dos argumentos: la media matriz (es decir, el triángulo)
a
de polinomios y un polinomio separadop
(denominado beta en el documento). Reducimos elm
problema de tamaño (inicialmentem=2*n
) recursivamente a dos problemas de tamañom-2
y devolvemos la diferencia de sus hafnianos. Una de ellas es usar lo mismoa
sin sus dos últimas filas, y lo mismop
. Otra es usar el triángulob[i][j] = a[i][j] + shmul(a[m-1][i],a[m-2][j]) + shmul(a[m-1][j],a[m-2][i])
(dondeshmul
está la operación shift-multiply en polinomios; es como un producto polinomial como de costumbre, adicionalmente multiplicado por la variable "x"; las potencias superiores a x ^ n son ignoradas), y el polinomio separadoq = p + shmul(p,a[m-1][m-2])
. Cuando la repetición realiza un tamaño-0a
, volvemos el mayor coeficiente de p:p[n]
.La operación shift-and-multiply se implementa en función
f(x,y,z)
. Se modificaz
en el lugar. Hablando libremente, lo hacez += shmul(x,y)
. Esta parece ser la parte más crítica para el rendimiento.Una vez que la recursión ha finalizado, debemos corregir el signo del resultado multiplicando por (-1) n .
fuente
-march=native
Hará una gran diferencia aquí. Al menos en TIO, casi reduce el tiempo de la pared a la mitad.Pitón
Esta es una implementación de referencia directa del Algoritmo 2 del documento mencionado . Los únicos cambios fueron mantener solo el valor actual de B , eliminando los valores de β solo actualizando g cuando i ∈ X , y multiplicando el polinomio truncado calculando solo los valores hasta el grado n .
Pruébalo en línea!
Aquí hay una versión más rápida con algunas de las optimizaciones fáciles.
Pruébalo en línea!
Para mayor diversión, aquí hay una implementación de referencia en J.
fuente
pmul
, usefor j in range(n-i):
y evite elif
j != k
conj < k
. Copia una submatriz en el caso else que se puede evitar cuando manejamos y eliminamos las dos últimas en lugar de las dos primeras filas y columnas. Y cuando calcula conx={1,2,4}
y más tarde conx={1,2,4,6}
, repite sus cálculos hastai=5
. Reemplacé la estructura de los dos bucles externos con el primer bucle activadoi
y luego asumiendo recursivamentei in X
yi not in X
. - Por cierto, puede ser interesante observar el crecimiento del tiempo necesario en comparación con los otros programas más lentos.Octava
Esto es básicamente una copia de la entrada de Dennis , pero optimizado para Octave. La optimización principal se realiza utilizando la matriz de entrada completa (y su transposición) y la recursión utilizando solo índices de matriz, en lugar de crear matrices reducidas.
La principal ventaja es la copia reducida de matrices. Mientras que Octave no tiene una diferencia entre punteros / referencias y valores y funcionalmente solo pasa por valor, es una historia diferente detrás de escena. Allí, se utiliza la copia en escritura (copia diferida). Eso significa que, para el código
a=1;b=a;b=b+1
, la variableb
solo se copia en una nueva ubicación en la última instrucción, cuando se cambia. Dadomatin
ymatransp
no se cambian, ellos nunca por copiado. La desventaja es que la función pasa más tiempo calculando los índices correctos. Puede que tenga que probar diferentes variaciones entre índices numéricos y lógicos para optimizar esto.Nota importante: la matriz de entrada debe ser
int32
! Guarde la función en un archivo llamadohaf.m
Ejemplo de script de prueba:
He probado esto usando TIO y MATLAB (en realidad nunca he instalado Octave). Me imagino que hacerlo funcionar es tan simple como
sudo apt-get install octave
. El comandooctave
cargará la interfaz gráfica de usuario de Octave. Si es más complicado que esto, eliminaré esta respuesta hasta que haya proporcionado instrucciones de instalación más detalladas.fuente
Recientemente Andreas Bjorklund, Brajesh Gupt y yo publicamos un nuevo algoritmo para hafnianos de matrices complejas: https://arxiv.org/pdf/1805.12498.pdf . Para una matriz n \ times n se escala como n ^ 3 2 ^ {n / 2}.
Si entiendo correctamente el algoritmo original de Andreas de https://arxiv.org/pdf/1107.4466.pdf se escala como n ^ 4 2 ^ {n / 2} o n ^ 3 log (n) 2 ^ {n / 2} si usaste transformadas de Fourier para hacer multiplicaciones polinómicas.
Nuestro algoritmo está específicamente diseñado para matrices complejas, por lo que no será tan rápido como los desarrollados aquí para matrices {-1,0,1}. Sin embargo, me pregunto si uno puede usar algunos de los trucos que usó para mejorar nuestra implementación. Además, si las personas están interesadas, me gustaría ver cómo funciona su implementación cuando se les dan números complejos en lugar de números enteros. Finalmente, cualquier comentario, crítica, mejoras, errores, mejoras son bienvenidos en nuestro repositorio https://github.com/XanaduAI/hafnian/
¡Salud!
fuente