Подтвердить что ты не робот

Матричное умножение в Clojure vs Numpy

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

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

Здесь 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)

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

UPDATE: я реализовал те же тесты, используя библиотеку 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
4b9b3361

Ответ 1

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

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

EDIT: забыл добавить первый шаг в настройке clojure, включить "предупреждать об отражении"

Ответ 2

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

В любой момент, когда вы выполняете нетривиальные операции с матрицей/вектором, вы, вероятно, должны ссылаться на BLAS/LAPACK.

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

Ответ 3

Я только что организовал небольшую перестрелку между Incanter 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 занимает вдвое больше памяти и только 15% быстрее, чем (mmult (trans ai) ai) в Incanter.

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


EDIT: Здесь приведены абсолютные числа (в мс.). Я попал на свой (довольно медленный) ПК:

Incanter: 665.362452
   jBLAS: 459.311598
   numpy: 353.777885

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

PS. Результаты Haskell hmatrix близки к numpy, но я не уверен, как правильно их исправить.

Ответ 4

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

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

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

Ответ 5

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

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

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

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

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

http://jblas.org/

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

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

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

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

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

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

Ответ 6

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

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

EDIT:

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

Ответ 7

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

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

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

Ответ 8

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

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

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

Ответ 9

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

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

С уважением.

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