Multiplication matricielle dans Clojure vs Numpy
Je travaille sur une application dans Clojure qui doit multiplier les grandes matrices et rencontre de gros problèmes de performances par rapport à une version Numpy identique. Numpy semble être capable de multiplier une matrice 1,000,000x23 par sa transposition en moins d'une seconde, tandis que le code Clojure équivalent prend plus de six minutes. (Je peux imprimer la matrice résultante de Numpy, donc c'est définitivement évaluer tout).
Est-ce que je fais quelque chose de terriblement mal dans ce code Clojure? Être il y a un truc de Numpy que je peux essayer d'imiter?
Voici le python:
import numpy as np
def test_my_mult(n):
A = np.random.rand(n*23).reshape(n,23)
At = A.T
t0 = time.time()
res = np.dot(A.T, A)
print time.time() - t0
print np.shape(res)
return res
# Example (returns a 23x23 matrix):
# >>> results = test_my_mult(1000000)
#
# 0.906938076019
# (23, 23)
Et le clojure:
(defn feature-vec [n]
(map (partial cons 1)
(for [x (range n)]
(take 22 (repeatedly rand)))))
(defn dot-product [x y]
(reduce + (map * x y)))
(defn transpose
"returns the transposition of a `coll` of vectors"
[coll]
(apply map vector coll))
(defn matrix-mult
[mat1 mat2]
(let [row-mult (fn [mat row]
(map (partial dot-product row)
(transpose mat)))]
(map (partial row-mult mat2)
mat1)))
(defn test-my-mult
[n afn]
(let [xs (feature-vec n)
xst (transpose xs)]
(time (dorun (afn xst xs)))))
;; Example (yields a 23x23 matrix):
;; (test-my-mult 1000 i/mmult) => "Elapsed time: 32.626 msecs"
;; (test-my-mult 10000 i/mmult) => "Elapsed time: 628.841 msecs"
;; (test-my-mult 1000 matrix-mult) => "Elapsed time: 14.748 msecs"
;; (test-my-mult 10000 matrix-mult) => "Elapsed time: 434.128 msecs"
;; (test-my-mult 1000000 matrix-mult) => "Elapsed time: 375751.999 msecs"
;; Test from wikipedia
;; (def A [[14 9 3] [2 11 15] [0 12 17] [5 2 3]])
;; (def B [[12 25] [9 10] [8 5]])
;; user> (matrix-mult A B)
;; ((273 455) (243 235) (244 205) (102 160))
Mise à jour: j'ai implémenté le même benchmark en utilisant la bibliothèque JBLAS et j'ai trouvé des améliorations de vitesse massives et massives. Merci à tous pour leur participation! Il est temps d'envelopper cette ventouse dans Clojure. Voici le nouveau code:
(import '[org.jblas FloatMatrix])
(defn feature-vec [n]
(FloatMatrix.
(into-array (for [x (range n)]
(float-array (cons 1 (take 22 (repeatedly rand))))))))
(defn test-mult [n]
(let [xs (feature-vec n)
xst (.transpose xs)]
(time (let [result (.mmul xst xs)]
[(.rows result)
(.columns result)]))))
;; user> (test-mult 10000)
;; "Elapsed time: 6.99 msecs"
;; [23 23]
;; user> (test-mult 100000)
;; "Elapsed time: 43.88 msecs"
;; [23 23]
;; user> (test-mult 1000000)
;; "Elapsed time: 383.439 msecs"
;; [23 23]
(defn matrix-stream [rows cols]
(repeatedly #(FloatMatrix/randn rows cols)))
(defn square-benchmark
"Times the multiplication of a square matrix."
[n]
(let [[a b c] (matrix-stream n n)]
(time (.mmuli a b c))
nil))
;; forma.matrix.jblas> (square-benchmark 10)
;; "Elapsed time: 0.113 msecs"
;; nil
;; forma.matrix.jblas> (square-benchmark 100)
;; "Elapsed time: 0.548 msecs"
;; nil
;; forma.matrix.jblas> (square-benchmark 1000)
;; "Elapsed time: 107.555 msecs"
;; nil
;; forma.matrix.jblas> (square-benchmark 2000)
;; "Elapsed time: 793.022 msecs"
;; nil
9 réponses
La version Python est compilée en boucle en C tandis que la version Clojure construit une nouvelle séquence intermédiaire pour chacun des appels à mapper dans ce code. Il est probable que la différence de performance que vous voyez provient de la différence de structures de données.
Pour obtenir mieux que cela, vous pourriez jouer avec une bibliothèque comme Incanter, ou écrire votre propre version comme expliqué dans ce DONC, la question. voir aussi celui-ci, Neandertal ou nd4j. Si vous voulez vraiment rester avec des séquences pour garder les propriétés d'évaluation paresseuses, etc. ensuite, vous pouvez obtenir un véritable coup de pouce en regardant dans transitoires pour les calculs de matrice interne
EDIT: vous avez oublié d'ajouter la première étape de réglage de clojure, activez "avertir sur la réflexion"
Numpy est lié à des routines Blas / Lapack qui ont été optimisées pendant des décennies au niveau de l'architecture de la machine tandis que le Clojure est une implémentation de la multiplication de la manière la plus simple et la plus naïve.
Chaque fois que vous avez des opérations matricielles/vectorielles non triviales à effectuer, vous devriez probablement créer un lien vers BLAS / LAPACK.
La seule fois où cela ne sera pas plus rapide est pour les petites matrices de langues où la surcharge de traduction de la représentation de données entre l'exécution de la langue et le LAPACK dépassent le temps passé à faire le calcul.
Je viens de mettre en scène une petite fusillade entre Incanter 1.3 et jBLAS 1.2.1. Voici le code:
(ns ml-class.experiments.mmult
[:use [incanter core]]
[:import [org.jblas DoubleMatrix]])
(defn -main [m]
(let [n 23 m (Integer/parseInt m)
ai (matrix (vec (double-array (* m n) (repeatedly rand))) n)
ab (DoubleMatrix/rand m n)
ti (copy (trans ai))
tb (.transpose ab)]
(dotimes [i 20]
(print "Incanter: ") (time (mmult ti ai))
(print " jBLAS: ") (time (.mmul tb ab)))))
Dans mon test, Incanter est toujours plus lent que jBLAS d'environ 45% en multiplication matricielle simple. Cependant, la fonction Incanter trans ne crée pas une nouvelle copie d'une matrice, et donc (.mmul (.transpose ab) ab) dans jBLAS prend deux fois plus de mémoire et est seulement 15% plus rapide que {[4] } dans Incanter.
Étant donné Incanter riche ensemble de fonctionnalités (en particulier, il est bibliothèque de traçage), Je ne pense pas que je vais passer à jBLAS de sitôt. Pourtant, j'aimerais voir une autre fusillade entre Jblas et Parallel Colt, et peut-être vaut-il la peine d'envisager de remplacer Parallel Colt par jBLAS dans Incanter? :-)
EDIT: Voici les nombres absolus (en msec.) Je suis arrivé sur mon PC (plutôt lent):
Incanter: 665.362452
jBLAS: 459.311598
numpy: 353.777885
Pour chaque bibliothèque, j'ai choisi le meilleur temps sur 20 courses, taille de la matrice 23x400000.
PS. Les résultats de Haskell hmatrix sont proches de numpy, mais je ne sais pas comment le comparer correctement.
Le code Numpy utilise des bibliothèques intégrées, écrites dans Fortran au cours des dernières décennies et optimisées par les auteurs, votre fournisseur de CPU et votre distributeur OS (ainsi que les personnes Numpy) pour des performances maximales. Vous venez de faire l'approche complètement directe et évidente de la multiplication matricielle. Il n'est pas surprenant, vraiment, que le rendement diffère.
Mais si vous insistez pour le faire dans Clojure, envisagez de rechercher de meilleurs algorithmes , en utilisant des boucles directes par opposition à des fonctions d'ordre supérieur comme reduce, ou trouver une bibliothèque d'algèbre matricielle appropriée pour Java (je doute qu'il y en ait de bonnes dans Clojure, mais je ne sais pas vraiment) écrites par un mathématicien compétent.
Enfin, regardez comment écrire correctement Clojure rapide. Utilisez des astuces de type, exécutez un profileur sur votre code (surprise! vous dot produit fonction est en utilisant le plus de temps), et laisser tomber les fonctionnalités de haut niveau à l'intérieur de vos boucles serrées.
Comme @littleidea et d'autres l'ont souligné, votre version numpy utilise LAPACK / BLAS / ATLAS qui sera beaucoup plus rapide que tout ce que vous faites dans clojure puisqu'elle a été finement ajustée pendant des années. :)
Cela dit, le plus gros problème avec le code Clojure est qu'il utilise des Doubles, comme dans les doubles en boîte. J'appelle cela le problème "double paresseux" et je l'ai rencontré au travail un certain nombre de fois. À partir de maintenant, même avec 1.3, les collections de clojure ne sont pas amicales primitives. (Vous pouvez en créer un vecteur de primitives mais cela ne vous aidera pas depuis tout le seq. les fonctions finiront par les boxe! Je devrais aussi dire que les améliorations primitives dans 1.3 sont assez agréables et finissent par aider.. nous ne sommes tout simplement pas 100% là WRT support primitif dans les collections.)
Lorsque vous faites n'importe quel type de calcul matriciel dans clojure, vous devez vraiment utiliser des tableaux java ou, mieux encore, des bibliothèques matricielles. Incanter utilise parrelcolt mais vous devez faire attention aux fonctions incanter que vous utilisez... depuis beaucoup de choses d'entre eux rendent les matrices seqable qui finit par boxe les doubles vous donnant des performances similaires à ce que vous voyez actuellement. (BTW, j'ai mes propres wrappers parrelcolt mis en place que je pourrais libérer Si vous pensez qu'ils seraient utiles.)
Pour utiliser les bibliothèques BLAS, vous avez quelques options en java-terre. Avec toutes ces options, vous devez payer une taxe JNA... toutes vos données doivent être copiées avant de pouvoir être traitées. Cette taxe est logique lorsque vous faites CPU lié opérations comme les décompositions matricielles et dont le temps de traitement prend plus de temps que le temps nécessaire pour copier les données. Pour des opérations plus simples avec de petites matrices, rester dans java-land sera probablement plus rapide. Vous avez juste besoin de faire quelques tests comme vous l'avez fait ci-dessus pour voir ce qui fonctionne le mieux pour vous.
Voici vos options pour utiliser BLAS à partir de java:
Http://code.google.com/p/netlib-java/
- API pour CI-DESSUS: http://code.google.com/p/matrix-toolkits-java/
Je dois souligner que parrelcolt utilise le projet netlib-java. Ce qui signifie, Je crois, si vous le configurez correctement, il utilisera BLAS. Cependant, je n'ai pas vérifiée cette. Pour une explication sur les différences entre jblas et netlib-java, voir ce fil que j'ai commencé à ce sujet sur la liste de diffusion de jblas:
Http://groups.google.com/group/jblas-users/browse_thread/thread/c9b3867572331aa5
Je devrait également souligner la bibliothèque de paquets Universal Java Matrix:
Http://sourceforge.net/projects/ujmp/
Il enveloppe toutes les bibliothèques que j'ai mentionnées, puis certaines! Je n'ai pas trop regardé L'API pour savoir à quel point leur abstraction est fuyante. Il semble comme un beau projet. J'ai fini par utiliser mes propres wrappers parrelcolt clojure car ils étaient assez rapides et j'ai vraiment aimé l'API colt. (Colt utilise des objets de fonction, ce qui signifie que j'ai pu il suffit de passer les fonctions clojure avec peu de problèmes!)
Si vous voulez faire des nombres dans Clojure, je vous recommande fortement d'utiliser Incanter plutôt que d'essayer de rouler vos propres fonctions matricielles et autres.
Incanter utilise Colt parallèle sous le capot, ce qui est assez rapide.
Modifier:
Début 2013, si vous voulez faire des numerics dans Clojure, je vous recommande fortement de vérifier core.la matrice
Numpy est hautement optimisé pour l'algèbre linéaire. Certainement pour les grandes matrices, où la plupart du traitement est dans le code C natif.
Afin de faire correspondre cette performance (en supposant que C'est possible en Java), vous devrez supprimer la plupart des abstractions de Clojure: n'utilisez pas map avec des fonctions anonymes lors de l'itération sur de grandes matrices, ajoutez des astuces de type pour activer l'utilisation de tableaux Java bruts, etc.
Probablement la meilleure option est juste d'utiliser une bibliothèque Java prête à l'emploi optimisée pour calcul numérique (http://math.nist.gov/javanumerics / ou similaire).
Je n'ai pas de réponses spécifiques pour vous; juste quelques suggestions.
- Utilisez un profileur pour déterminer où le temps est passé
- définissez warn-on-reflection et utilisez des astuces de type si nécessaire
- vous devrez peut-être abandonner certaines constructions de haut niveau et aller avec loop-recur pour sqeeze cette dernière once de performance
IME, le code Clojure devrait fonctionner assez près de Java (2 ou 3X). Mais vous devez travailler dessus.
N'utilisez map () que si cela a du sens. Ce qui signifie: Si vous avez un problème spécifique comme multiplier deux matrices, n'essayez pas de le mapper (), il suffit de multiplier les matrices.
J'ai tendance à utiliser map() seulement quand cela a un sens linguistique (c'est-à-dire si le programme est vraiment plus lisible que sans lui). Multiplier les matrices est une boucle si évidente que le mappage n'a aucun sens.
Le vôtre.
Pedro Fortuny.