performance - matriz - multiplicacion de matrices 3x3
¿Producto matricial de funcionalidad pura razonablemente eficiente en Haskell? (3)
Tan eficiente como, por ejemplo, Java. Para ser concretos, supongamos que estoy hablando de un diseño simple de triple lazo, precisión simple, columna contigua (float [], no float [] []) y matrices de tamaño 1000x1000, y una CPU de un solo núcleo. (Si obtiene 0.5-2 operaciones de coma flotante por ciclo, probablemente esté en el estadio)
Así que algo así como
public class MatrixProd {
static float[] matProd(float[] a, int ra, int ca, float[] b, int rb, int cb) {
if (ca != rb) {
throw new IllegalArgumentException("Matrices not fitting");
}
float[] c = new float[ra*cb];
for(int i = 0; i < ra; ++i) {
for(int j = 0; j < cb; ++j) {
float sum = 0;
for(int k = 0; k < ca; ++k) {
sum += a[i*ca+k]*b[k*cb+j];
}
c[i*cb+j] = sum;
}
}
return c;
}
static float[] mkMat(int rs, int cs, float x, float d) {
float[] arr = new float[rs*cs];
for(int i = 0; i < rs; ++i) {
for(int j = 0; j < cs; ++j) {
arr[i*cs+j] = x;
x += d;
}
}
return arr;
}
public static void main(String[] args) {
int sz = 100;
float strt = -32, del = 0.0625f;
if (args.length > 0) {
sz = Integer.parseInt(args[0]);
}
if (args.length > 1) {
strt = Float.parseFloat(args[1]);
}
if (args.length > 2) {
del = Float.parseFloat(args[2]);
}
float[] a = mkMat(sz,sz,strt,del);
float[] b = mkMat(sz,sz,strt-16,del);
System.out.println(a[sz*sz-1]);
System.out.println(b[sz*sz-1]);
long t0 = System.currentTimeMillis();
float[] c = matProd(a,sz,sz,b,sz,sz);
System.out.println(c[sz*sz-1]);
long t1 = System.currentTimeMillis();
double dur = (t1-t0)*1e-3;
System.out.println(dur);
}
}
¿Supongo? (No había leído las especificaciones correctamente antes de la codificación, por lo que el diseño es principal, pero dado que el patrón de acceso es el mismo, eso no hace una diferencia ya que mezclar diseños, así que supongo que está bien).
No he dedicado tiempo a pensar en un algoritmo inteligente o trucos de optimización de bajo nivel (de todos modos, no lograría mucho en Java). Acabo de escribir el bucle simple, porque
No quiero que esto suene como un desafío, pero tenga en cuenta que Java puede satisfacer todo lo anterior fácilmente
Y eso es lo que Java da fácilmente , así que lo tomaré.
(Si obtiene 0.5-2 operaciones de coma flotante por ciclo, probablemente esté en el estadio)
No estoy cerca, me temo, ni en Java ni en Haskell. Demasiadas fallas de caché para alcanzar ese rendimiento con el triple lazo simple.
Haciendo lo mismo en Haskell, una vez más sin pensar en ser inteligente, un sencillo bucle triple directo:
{-# LANGUAGE BangPatterns #-}
module MatProd where
import Data.Array.ST
import Data.Array.Unboxed
matProd :: UArray Int Float -> Int -> Int -> UArray Int Float -> Int -> Int -> UArray Int Float
matProd a ra ca b rb cb =
let (al,ah) = bounds a
(bl,bh) = bounds b
{-# INLINE getA #-}
getA i j = a!(i*ca + j)
{-# INLINE getB #-}
getB i j = b!(i*cb + j)
{-# INLINE idx #-}
idx i j = i*cb + j
in if al /= 0 || ah+1 /= ra*ca || bl /= 0 || bh+1 /= rb*cb || ca /= rb
then error $ "Matrices not fitting: " ++ show (ra,ca,al,ah,rb,cb,bl,bh)
else runSTUArray $ do
arr <- newArray (0,ra*cb-1) 0
let outer i j
| ra <= i = return arr
| cb <= j = outer (i+1) 0
| otherwise = do
!x <- inner i j 0 0
writeArray arr (idx i j) x
outer i (j+1)
inner i j k !y
| ca <= k = return y
| otherwise = inner i j (k+1) (y + getA i k * getB k j)
outer 0 0
mkMat :: Int -> Int -> Float -> Float -> UArray Int Float
mkMat rs cs x d = runSTUArray $ do
let !r = rs - 1
!c = cs - 1
{-# INLINE idx #-}
idx i j = cs*i + j
arr <- newArray (0,rs*cs-1) 0
let outer i j y
| r < i = return arr
| c < j = outer (i+1) 0 y
| otherwise = do
writeArray arr (idx i j) y
outer i (j+1) (y + d)
outer 0 0 x
y el módulo que llama
module Main (main) where
import System.Environment (getArgs)
import Data.Array.Unboxed
import System.CPUTime
import Text.Printf
import MatProd
main :: IO ()
main = do
args <- getArgs
let (sz, strt, del) = case args of
(a:b:c:_) -> (read a, read b, read c)
(a:b:_) -> (read a, read b, 0.0625)
(a:_) -> (read a, -32, 0.0625)
_ -> (100, -32, 0.0625)
a = mkMat sz sz strt del
b = mkMat sz sz (strt - 16) del
print (a!(sz*sz-1))
print (b!(sz*sz-1))
t0 <- getCPUTime
let c = matProd a sz sz b sz sz
print $ c!(sz*sz-1)
t1 <- getCPUTime
printf "%.6f/n" (fromInteger (t1-t0)*1e-12 :: Double)
Así que estamos haciendo casi exactamente las mismas cosas en ambos idiomas. Compila el Haskell con -O2 , el Java con javac
$ java MatrixProd 1000 "-13.7" 0.013
12915.623
12899.999
8.3592897E10
8.193
$ ./vmmult 1000 "-13.7" 0.013
12915.623
12899.999
8.35929e10
8.558699
Y los tiempos resultantes son bastante cercanos.
Y si compilamos el código Java en native, con gcj -O3 -Wall -Wextra --main=MatrixProd -fno-bounds-check -fno-store-check -o jmatProd MatrixProd.java ,
$ ./jmatProd 1000 "-13.7" 0.013
12915.623
12899.999
8.3592896512E10
8.215
todavía no hay una gran diferencia.
Como bonificación especial, el mismo algoritmo en C (gcc-O3):
$ ./cmatProd 1000 "-13.7" 0.013
12915.623047
12899.999023
8.35929e+10
8.079759
Así que esto no revela ninguna diferencia fundamental entre Java directo y Haskell directo cuando se trata de tareas computacionalmente intensivas usando números de coma flotante (cuando se trata de números enteros en números grandes o medianos, el uso de GMP por GHC hace que Haskell supere a BigInteger de Java por un gran margen para muchas tareas, pero eso es, por supuesto, un problema de biblioteca, no de lenguaje), y ambos están cerca de C con este algoritmo.
Con toda justicia, sin embargo, eso se debe a que el patrón de acceso causa una falta de caché cada dos nanosegundos, por lo que en los tres idiomas este cálculo está ligado a la memoria.
Si mejoramos el patrón de acceso multiplicando una matriz principal de fila con una matriz columna-mayor, todos se vuelven más rápidos, la C compilada por gcc termina 1.18s, Java toma 1.23s y la Haskell compilada por ghc toma alrededor de 5.8s, que se puede reducir a 3 segundos utilizando el servidor de fondo llvm.
Aquí, el control de rango de la biblioteca de arreglo realmente duele. Usando el acceso a la matriz no verificada (como debería, después de buscar errores, ya que las comprobaciones ya están hechas en el código que controla los bucles), el backend nativo de GHC termina en 2.4s, pasando por el backend de llvm permite que el cálculo termine en 1.55s, que es decente, aunque significativamente más lento que C y Java. Utilizando las primitivas de GHC.Prim lugar de la biblioteca de matriz , el backend llvm produce código que se ejecuta en 1.16s (nuevamente, sin verificación de límites en cada acceso, pero que solo se producen índices válidos durante el cálculo, en este caso puede ser fácilmente probado antes, así que aquí, no se sacrifica la seguridad de la memoria¹; la comprobación de cada acceso aumenta el tiempo a 1.96s, aún significativamente mejor que la comprobación de límites de la biblioteca de arreglos ).
En pocas palabras: GHC necesita una ramificación (mucho) más rápida para la verificación de límites, y hay margen de mejora en el optimizador, pero en principio, "el enfoque de Haskell (pureza codificada en el sistema de tipos) es compatible con eficiencia, memoria-seguridad y simplicidad ", todavía no estamos allí. Por el momento, uno tiene que decidir cuánto de qué punto está dispuesto a sacrificar.
¹ Sí, ese es un caso especial, en general, omite los límites: el control sacrifica la seguridad de la memoria, o al menos es más difícil probar que no lo hace.
Conozco a Haskell un poco, y me pregunto si es posible escribir algo así como un producto de matriz matricial en Haskell que sea todo lo siguiente:
- Funcionalidad pura : sin
IOo mónadas deStateen su firma de tipo (no me importa lo que sucede en el cuerpo de la función. Es decir, no me importa si el cuerpo de la función usa mónadas, siempre que la función sea pura ) . Es posible que desee utilizar este producto matriz-matriz en una función pura. - Memoria segura: sin malloc o punteros . Sé que es posible "escribir C" en Haskell, pero se pierde la seguridad de la memoria. En realidad, escribir este código en C e interconectarlo con Haskell también pierde seguridad en la memoria.
- Tan eficiente como, por ejemplo, Java. Para ser concretos, supongamos que estoy hablando de un diseño simple de triple lazo, precisión simple, columna contigua (
float[], nofloat[][]) y matrices de tamaño 1000x1000, y una CPU de un solo núcleo. (Si obtiene 0.5-2 operaciones de punto flotante por ciclo, probablemente esté en el estadio).
(No quiero que esto suene como un desafío, pero tenga en cuenta que Java puede satisfacer todo lo anterior más fácilmente).
eso ya lo se
- La implementación de triple circuito no es la más eficiente . Es bastante caché-ajeno. Es mejor usar una implementación BLAS bien escrita en este caso particular. Sin embargo, uno no siempre puede contar con que una biblioteca C esté disponible para lo que se intenta hacer. Me pregunto si se puede escribir un código razonablemente eficiente en Haskell normal.
- Algunas personas escribieron trabajos de investigación completos que demuestran # 3 . Sin embargo, no soy un investigador de ciencias de la computación. Me pregunto si es posible mantener simples las cosas simples en Haskell.
- The Gentle Introduction to Haskell tiene una implementación de producto matricial . Sin embargo, no satisfaría los requisitos anteriores.
Dirigiendo comentarios:
Tengo tres razones: primero, el requisito de "no malloc o punteros" aún no está bien definido (lo desafío a que escriba cualquier parte del código Haskell que no utilice punteros);
Vi muchos programas Haskell que no usaban Ptr . ¿Quizás se refiere al hecho de que en el nivel de instrucción de la máquina, se usarán punteros? Eso no es lo que quise decir. Me refería al nivel de abstracción del código fuente de Haskell.
segundo, el ataque a la investigación de CS está fuera de lugar (y además no puedo imaginar nada más simple que usar un código que alguien más ya haya escrito para ti); tercero, hay muchos paquetes de matriz en Hackage (y el trabajo de preparación para hacer esta pregunta debe incluir revisar y rechazar cada uno).
Parece que sus # 2 y # 3 son iguales ("use las bibliotecas existentes"). Estoy interesado en el producto de la matriz como una prueba simple de lo que Haskell puede hacer por sí mismo, y si le permite "mantener simples las cosas simples". Fácilmente podría haber surgido un problema numérico que no tiene ninguna biblioteca preparada, pero luego tendría que explicar el problema, mientras que todos ya saben qué es un producto de matriz.
¿Cómo puede Java posiblemente satisfacer 1.? Cualquier método de Java es esencialmente
:: IORef Arg -> ... -> IORef This -> IO Ret
Esto va a la raíz de mi pregunta, en realidad (+1). Si bien Java no pretende rastrear la pureza, Haskell sí. En Java, si la función es pura o no está indicada en los comentarios. Puedo afirmar que el producto de la matriz es puro, aunque hago una mutación en el cuerpo de la función. La pregunta es si el enfoque de Haskell (pureza codificada en el sistema de tipos) es compatible con la eficiencia, la seguridad de la memoria y la simplicidad.
Al igual que Java, Haskell no es el mejor lenguaje para escribir código numérico.
La generación de código numérico pesado de Haskell es ... promedio. No ha contado con los años de investigación que tienen los gustos de Intel y GCC.
Lo que Haskell le ofrece, en su lugar, es una forma de interconectar limpiamente su código "rápido" con el resto de su aplicación. Recuerde que el 3% del código es responsable del 97% del tiempo de ejecución de su aplicación. 1
Con Haskell, tiene una manera de llamar a estas funciones altamente optimizadas de una manera que interactúa extremadamente bien con el resto de su código: a través de la muy agradable interfaz C Foreign Function. De hecho, si lo desea, podría escribir su código numérico en el lenguaje ensamblador de su arquitectura y obtener aún más rendimiento. La inmersión en C para obtener partes pesadas de tu aplicación no es un error, es una característica.
Pero yo divago.
Al tener estas funciones altamente optimizadas aisladas, y con una interfaz similar al resto de su código Haskell, puede realizar optimizaciones de alto nivel con las poderosas reglas de reescritura de Haskell, que le permiten escribir reglas como reverse . reverse == id reverse . reverse == id que reduce automágicamente expresiones complejas en tiempo de compilación 2 . Esto conduce a bibliotecas extremadamente rápidas, puramente funcionales y fáciles de usar, como Data.Text 3 y Data.Vector [4].
Al combinar altos y bajos niveles de optimización, terminamos con una implementación mucho más optimizada, con cada mitad ("C / asm" y "Haskell") relativamente fácil de leer. La optimización de bajo nivel se realiza en su lengua nativa (C o ensamblaje), la optimización de alto nivel obtiene un DSL especial (reglas de reescritura de Haskell), y el resto del código no lo reconoce por completo.
En conclusión, sí, Haskell puede ser más rápido que Java. Pero hace trampa al pasar por C para los FLOPS crudos. Esto es mucho más difícil de hacer en Java (además de tener una sobrecarga mucho mayor para el FFI de Java), por lo que se evita. En Haskell, es natural. Si su aplicación pasa una cantidad exorbitante de tiempo haciendo cálculos numéricos, entonces tal vez en lugar de mirar a Haskell o Java, mire a Fortran para sus necesidades. Si su aplicación pasa una gran parte de su tiempo en una pequeña parte del código sensible al rendimiento, entonces el Haskell FFI es su mejor opción. Si su aplicación no pierde tiempo en el código numérico ... entonces use lo que quiera. =)
Haskell (ni Java, para el caso) no es Fortran.
1 Estos números fueron inventados, pero entiendes mi punto.
2 http://www.cse.unsw.edu.au/~dons/papers/CLS07.html
3 http://hackage.haskell.org/package/text
[4] http://hackage.haskell.org/package/vector
Ahora que eso está fuera del camino, para responder a su pregunta real:
No, actualmente no es inteligente escribir tus multiplicaciones de matriz en Haskell. Por el momento, REPA es la manera canónica de hacer esto [5]. La implementación rompe parcialmente la seguridad de la memoria (usan unsliceSlice), pero la "seguridad de la memoria rota" está aislada de esa función, en realidad es muy segura (pero el compilador no la verifica fácilmente) y es fácil de eliminar si las cosas van mal (reemplazar) inseguroSlice "con" slice ").
¡Pero este es Haskell! Muy rara vez se toman las características de rendimiento de una función de forma aislada. Eso puede ser algo malo (en el caso de fugas de espacio), o algo muy, muy bueno.
Aunque el algoritmo de multiplicación de la matriz utilizado es ingenuo, tendrá un peor rendimiento en un punto de referencia sin procesar. Pero rara vez nuestro código se ve como puntos de referencia.
¿Qué pasa si fueras un científico con millones de puntos de datos y quieres multiplicar matrices enormes? [7]
Para esas personas, tenemos mmultP [6]. Esto realiza la multiplicación de la matriz, pero es paralela a los datos y está sujeta al paralelismo de datos anidados de REPA. También tenga en cuenta que el código esencialmente no se modifica desde la versión secuencial.
Para aquellas personas que no multiplican grandes matrices, y en su lugar multiplican muchas matrices pequeñas, tiende a haber otro código que interactúa con dichas matrices. Posiblemente cortarlo en vectores de columna y encontrar sus productos de puntos, tal vez encontrar sus valores propios, tal vez algo completamente distinto. A diferencia de C, Haskell sabe que, aunque le gusta resolver problemas de forma aislada, la solución más eficiente no suele encontrarse allí.
Al igual que ByteString, Text y Vector, las matrices REPA están sujetas a fusión. 2 Por cierto, debería leer 2 , es un documento muy bien escrito. Esto, combinado con la incorporación agresiva del código relevante y la naturaleza altamente paralela de REPA, nos permite expresar estos conceptos matemáticos de alto nivel con optimizaciones de alto nivel muy avanzadas entre bastidores.
Aunque actualmente no se conoce un método para escribir una multiplicación eficiente de matrices en lenguajes funcionales puros, podemos acercarnos un poco (sin vectorización automática, algunas desreferencias excesivas para llegar a los datos reales, etc.), pero nada cerca de lo que IFORT o GCC puede hacer. Pero los programas no existen en una isla, y hacer que la isla en su conjunto funcione bien es mucho, mucho más fácil en Haskell que en Java.
[7] Acutally, la mejor manera de hacerlo es mediante el uso de la GPU. Hay algunas DSL de GPU disponibles para Haskell que hacen esto posible de forma nativa. ¡Están realmente limpios!
Hay dos ángulos para atacar este problema.
La investigación, en esta línea, está en curso. Ahora, hay muchos programadores de Haskell que son más inteligentes que yo; un hecho que constantemente recuerdo y me siento humilde. Uno de ellos puede venir y corregirme, pero no conozco ninguna forma simple de componer primitivos Haskell seguros en una rutina de multiplicación de matrices de primera línea. Esos documentos de los que hablas suenan como un buen comienzo.
Sin embargo, no soy un investigador de ciencias de la computación. Me pregunto si es posible mantener simples las cosas simples en Haskell.
Si citan esos documentos, tal vez podríamos ayudarlos a descifrarlos.
La ingeniería del software, en esta línea, es bien comprendida, directa e incluso fácil. Un codificador inteligente de Haskell usaría un envoltorio delgado alrededor de BLAS, o buscaría una envoltura de este tipo en Hackage.
Descifrar la investigación de vanguardia es un proceso continuo que traslada el conocimiento de los investigadores a los ingenieros. Fue un investigador en informática, CAR Hoare, quien descubrió el quicksort y publicó un artículo al respecto. Hoy en día, es un graduado de ciencias de la computación raro que no puede implementar personalmente quicksort de memoria (al menos, los que se graduaron recientemente).
Un poco de historia
Casi esta misma pregunta se ha formulado en la historia varias veces antes.
¿Es posible escribir aritmética de matriz en Fortran que sea tan rápido como el ensamblaje?
¿Es posible escribir la aritmética de la matriz en C que es tan rápido como Fortran?
¿Es posible escribir una aritmética de matriz en Java que sea tan rápida como C?
¿Es posible escribir una aritmética de matriz en Haskell que sea tan rápida como Java?
Hasta ahora, la respuesta siempre ha sido "todavía no", seguido de "lo suficientemente cerca". Los avances que lo hacen posible provienen de las mejoras en la escritura del código, las mejoras en los compiladores y las mejoras en el lenguaje de programación en sí.
Como ejemplo específico, C no fue capaz de superar Fortran en muchas aplicaciones del mundo real hasta que los compiladores C99 se generalizaron en la última década. En Fortran, se supone que las diferentes matrices tienen un almacenamiento distinto unas de otras, mientras que en C no suele ser el caso. Por lo tanto, se permitió a los compiladores de Fortran realizar optimizaciones que los compiladores de C no podían. Bueno, no hasta que salga C99 y puedas agregar el calificador restrict a tu código.
Los compiladores Fortran esperaron. Finalmente, los procesadores se volvieron lo suficientemente complejos como para que la escritura de un buen ensamblado se volviera más difícil, y los compiladores se volvieron lo suficientemente sofisticados como para que el Fortran fuera rápido.
Luego, los programadores de C esperaron hasta la década de 2000 para escribir código que coincidiera con Fortran. Hasta ese momento, usaron bibliotecas escritas en Fortran o ensamblador (o ambas), o aguantaron la velocidad reducida.
Los programadores Java, asimismo, tuvieron que esperar compiladores JIT y tuvieron que esperar a que aparecieran optimizaciones específicas. Los compiladores JIT fueron originalmente un concepto de investigación esotérica hasta que se convirtieron en parte de la vida cotidiana. La optimización de comprobación de límites también fue necesaria para evitar una prueba y una bifurcación para cada acceso a la matriz.
De vuelta a Haskell
Entonces, está claro que los programadores Haskell están "esperando", al igual que los programadores Java, C y Fortran antes que ellos. ¿Qué estamos esperando?
Tal vez solo estamos esperando que alguien escriba el código y nos muestre cómo se hace.
Quizás estamos esperando que los compiladores mejoren.
Tal vez estamos esperando una actualización del lenguaje Haskell en sí.
Y tal vez estamos esperando alguna combinación de lo anterior.
Acerca de la pureza
La pureza y las mónadas se fusionan mucho en Haskell. La razón de esto es porque en Haskell, las funciones impuras siempre usan la mónada IO . Por ejemplo, la mónada de State es 100% pura. Por lo tanto, cuando dices "firma pura" y "tipo no usa la mónada de State ", en realidad son requisitos completamente independientes y separados.
Sin embargo, también puede usar la mónada IO en la implementación de funciones puras, y de hecho, es bastante fácil:
addSix :: Int -> Int
addSix n = unsafePerformIO $ return (n + 6)
De acuerdo, sí, es una función estúpida, pero es pura. Incluso es obviamente puro. La prueba de pureza es doble:
¿Da el mismo resultado para las mismas entradas? Sí.
¿Produce algún efecto secundario semánticamente significativo? No.
La razón por la que nos gusta la pureza es porque las funciones puras son más fáciles de componer y manipular que las funciones impuras. Cómo se implementan no importa tanto. No sé si eres consciente de esto, pero Integer y ByteString son básicamente envoltorios alrededor de funciones C impuras, aunque la interfaz sea pura. (Hay trabajo en una nueva implementación de Integer , no sé qué tan lejos está).
Respuesta final
La pregunta es si el enfoque de Haskell (pureza codificada en el sistema de tipos) es compatible con la eficiencia, la seguridad de la memoria y la simplicidad.
La respuesta a esa parte es "sí", ya que podemos tomar funciones simples de BLAS y ponerlas en un envoltorio puro y seguro. El tipo de envoltura codifica la seguridad de la función, a pesar de que el compilador Haskell no puede probar que la implementación de la función es pura. Nuestro uso de unsafePerformIO en su implementación es un reconocimiento de que hemos probado la pureza de la función, y también es una concesión que no pudimos encontrar una manera de expresar esa prueba en el sistema de tipos de Haskell.
Pero la respuesta es "todavía no", ya que no sé cómo implementar la función completamente en Haskell como tal.
La investigación en esta área está en curso. La gente está buscando sistemas de prueba como Coq y nuevos lenguajes como Agda , así como desarrollos en GHC. Para ver qué tipo de sistema necesitaríamos para demostrar que las rutinas BLAS de alto rendimiento se pueden usar de forma segura. Estas herramientas también se pueden usar con otros lenguajes como Java. Por ejemplo, podría escribir una prueba en Coq de que su implementación de Java es pura.
Me disculpo por la respuesta "sí y no", pero ninguna otra respuesta reconocería tanto las contribuciones de los ingenieros (a quienes les importa "sí") como a los investigadores (a quienes les importa "todavía no").
PD Por favor cite los papeles.