Умножение матриц в Clojure против Numpy

Question

Я работаю над приложением в Clojure, которое должно умножать большие матрицы, и сталкиваюсь с некоторыми большими проблемами производительности по сравнению с идентичной версией Numpy. Кажется, что Numpy может умножить матрицу 1,000,000x23 на транспонирование менее чем за секунду, в то время как эквивалентный код clojure занимает более шести минут. (Я могу распечатать полученную матрицу из Numpy, так что она определенно оценивает все).

Я что-то не так делаю в этом коде Clojure? Есть ли какой-нибудь трюк с Numpy, который я могу попытаться повторить?

Вот этот питон:

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)

И замыкание:

(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))

ОБНОВЛЕНИЕ: я реализовал тот же эталонный тест с использованием библиотеки JBLAS и обнаружил значительные, огромные улучшения скорости. Спасибо всем за их вклад! Время обернуть эту присоску в Clojure. Вот новый код:

(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

Answer 1

Версия Python компилируется в цикл на C, а версия Clojure создает новую промежуточную последовательность для каждого из вызовов map в этом коде.Вероятно, разница в производительности объясняется различием структур данных.

Чтобы стать лучше, вы можете поиграть с такой библиотекой, как Incanter или написать свою собственную версию, как объяснено в в этом вопросе .см. также этот , неандерталец или nd4j .Если вы действительно хотите остаться с последовательностями, чтобы сохранить ленивые свойства оценки и т. Д., То вы можете получить реальный импульс, взглянув на переходные процессы для внутренних вычислений матрицы

РЕДАКТИРОВАТЬ: забыли добавитьпервый шаг в настройке clojure, включите «Предупреждать об отражении»

Answer 2

Numpy связывается с процедурами BLAS / Lapack, которые десятилетиями были оптимизированы на уровне архитектуры машин, в то время как Clojure реализует умножение самым простым и наивным способом.

Каждый раз, когда у вас нетТривиальные матричные / векторные операции для выполнения, вы, вероятно, должны ссылаться на BLAS / LAPACK.

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

Answer 3

Я только что устроил небольшую перестрелку между Инкантер 1.3 и jBLAS 1.2.1. Вот код:

(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)))))

В моем тесте Incanter постоянно медленнее, чем jBLAS, примерно на 45% при умножении в простой матрице. Однако функция Incanter trans не создает новую копию матрицы, и поэтому (.mmul (.transpose ab) ab) в jBLAS занимает в два раза больше памяти и всего на 1011 * 15% быстрее, чем (mmult (trans ai) ai) в Incanter.

Учитывая богатый набор функций Incanter (особенно это графическая библиотека), я не думаю, что переключусь на jBLAS в ближайшее время. Тем не менее, я хотел бы увидеть очередную перестрелку между jBLAS и Parallel Colt, и, может быть, стоит подумать о замене Parallel Colt на jBLAS в Incanter? : -)

---

РЕДАКТИРОВАНИЕ: Вот абсолютные числа (в мсек.), Которые я получил на своем (довольно медленном) ПК:

Incanter: 665.362452
   jBLAS: 459.311598
   numpy: 353.777885

Для каждой библиотеки я выбрал лучшее время из 20 прогонов, размер матрицы 23x400000.

PS. Результаты Haskell Hmatrix близки к значительным, но я не уверен, как правильно его сравнить.

Answer 4

Код Numpy использует встроенные библиотеки, написанные на Фортране за последние несколько десятилетий и оптимизированные авторами, вашим поставщиком ЦП и вашим дистрибьютором ОС (а также сотрудниками Numpy) для максимальной производительности.Вы только что сделали совершенно прямой, очевидный подход к умножению матриц.Неудивительно, что производительность отличается.

Но если вы настаиваете на том, чтобы сделать это в Clojure, подумайте о поиске лучших алгоритмов , используя прямые циклы в отличие от функций более высокого порядкакак reduce, или найти подходящую библиотеку матричной алгебры для Java (я сомневаюсь, что в Clojure есть хорошие библиотеки, но я не знаю точно), написанную компетентным математиком.

Наконец, посмотрите, какправильно пиши быстро Clojure.Используйте подсказки типа, запустите профилировщик для своего кода (удивительно! Вы работаете с точечной функцией больше всего времени) и бросьте высокоуровневые функции в свои узкие петли.

Answer 5

Как @littleidea и другие пользователи указали, что ваша версия numpy использует LAPACK / BLAS / ATLAS, которая будет намного быстрее, чем все, что вы делаете в clojure, поскольку она была настроена годами.:)

Тем не менее, самая большая проблема с кодом Clojure заключается в том, что он использует Double, как в штучной упаковке.Я называю это проблемой «ленивый двойник», и я сталкивался с ней на работе несколько раз.На данный момент, даже с 1.3, коллекции clojure не являются примитивно дружественными.(Вы можете создать вектор примитивов, но это вам не поможет, так как все функции seq. В конечном итоге их упаковывают! Я должен также сказать, что примитивные улучшения в 1.3 довольно хороши и в конечном итоге помогают ... мы простоне существует 100% поддержки примитивов WRT в коллекциях.)

При выполнении любой математической математики в clojure вам действительно необходимо использовать Java-массивы или, еще лучше, матричные библиотеки.Incanter действительно использует parrelcolt, но вы должны быть осторожны с тем, какие функции incanter вы используете ... так как многие из них делают матрицы пригодными для секвенции, что приводит к тому, что они удваивают двойные значения, давая вам производительность, аналогичную той, что вы видите в настоящее время.(Кстати, у меня есть свои настроенные упаковщики parrelcolt, которые я мог бы выпустить, если вы думаете, что они будут полезны.)

Чтобы использовать библиотеки BLAS, у вас есть несколько вариантов в java-land.Со всеми этими опциями вы должны заплатить налог JNA ... все ваши данные должны быть скопированы, прежде чем они могут быть обработаны.Этот налог имеет смысл, когда вы выполняете связанные с процессором операции, такие как матричные декомпозиции, и время обработки которых занимает больше времени, чем время, необходимое для копирования данных.Для более простых операций с небольшими матрицами пребывание в java-земле, вероятно, будет быстрее.Вам просто нужно сделать несколько тестов, как вы сделали выше, чтобы увидеть, что лучше всего подходит для вас.

Вот ваши варианты использования BLAS из Java:

http://jblas.org/

http://code.google.com/p/netlib-java/

• API для выше: http://code.google.com/p/matrix-toolkits-java/

Я должен отметить, что parrelcolt использует проект netlib-java.Это значит, я верю, что если вы настроите его правильно, он будет использовать BLAS.Однако я не проверил это.Объяснение различий между jblas и netlib-java смотрите в этой ветке, которую я начал об этом в списке рассылки jblas:

http://groups.google.com/group/jblas-users/browse_thread/thread/c9b3867572331aa5

Я также должен указать на библиотеку Universal Java Matrix Package:

http://sourceforge.net/projects/ujmp/

Он охватывает все библиотеки, которые я упомянул, а затем и некоторые!Я не слишком много смотрел на API, хотя, чтобы знать, насколько утечки их абстракция.Это похоже на хороший проект.Я закончил тем, что использовал свои собственные оболочки parrelcolt clojure, поскольку они были достаточно быстрыми, и мне действительно очень понравился API-интерфейс colt.(Colt использует функциональные объекты, что означает, что я смог просто передать закрывающие функции без особых проблем!)

Answer 6

Если вы хотите использовать числовые значения в Clojure, я настоятельно рекомендую использовать Incanter вместо того, чтобы пытаться свернуть свои собственные матричные функции и тому подобное.

Incanter использует Parallel Colt под капотом, что довольно быстро.

EDIT:

По состоянию на начало 2013 года, если вы хотите делать цифры в Clojure, я настоятельно рекомендую проверить core.matrix

Answer 7

Numpy высоко оптимизирован для линейной алгебры. Конечно, для больших матриц, где большая часть обработки выполняется в собственном коде C.

Чтобы соответствовать этой производительности (предполагая, что это возможно в Java), вам придется убрать большинство абстракций Clojure: не используйте map с анонимными функциями при переборе больших матриц, добавляйте подсказки типов, чтобы разрешить использование сырой Java массивы и т. д.

Вероятно, лучший вариант - просто использовать готовую библиотеку Java, оптимизированную для численных вычислений (http://math.nist.gov/javanumerics/ или аналогичную).

Answer 8

Используйте map (), только если это имеет смысл. Это означает: если у вас есть конкретная проблема, такая как умножение двух матриц, не пытайтесь отобразить ее (), просто умножьте матрицы.

Я склонен использовать map () только тогда, когда это имеет лингвистический смысл (т. Е. Если программа действительно более читабельна, чем без нее). Умножение матриц является настолько очевидным циклом, что отображать его бессмысленно.

Ваш.

Педро Фортуни.

Answer 9

У меня нет конкретных ответов для вас; только несколько предложений.

1. Используйте профилировщик, чтобы выяснить, где тратится время
2. установить предупреждение об отражении и использовать подсказки типа при необходимости
3. Возможно, вам придется отказаться от некоторых высокоуровневых конструкций и пойти с loop-recur, чтобы подвести последнюю унцию производительности

IME, код Clojure должен работать довольно близко к Java (2 или 3X). Но над этим нужно работать.