Векторизация индексной операции для матрицы scipy.sparse - PullRequest
6 голосов
/ 08 марта 2010

Следующий код выполняется слишком медленно, хотя кажется, что все векторизовано.

from numpy import *
from scipy.sparse import *

n = 100000;
i = xrange(n); j = xrange(n);
data = ones(n);

A=csr_matrix((data,(i,j)));

x = A[i,j]

Кажется, проблема в том, что операция индексирования реализована как функция python, и при вызове A[i,j] получается следующий результат профилирования

         500033 function calls in 8.718 CPU seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
   100000    7.933    0.000    8.156    0.000 csr.py:265(_get_single_element)
        1    0.271    0.271    8.705    8.705 csr.py:177(__getitem__)
(...)

А именно, функция python _get_single_element вызывается 100000 раз, что действительно неэффективно. Почему это не реализовано в чистом C? Кто-нибудь знает способ обойти это ограничение и ускорить приведенный выше код? Должен ли я использовать другой тип разреженной матрицы?

Ответы [ 2 ]

7 голосов
/ 09 марта 2010

Вы можете использовать A.diagonal(), чтобы получить диагональ намного быстрее (0,0009 секунды против 3,8 секунды на моем аппарате). Однако, если вы хотите выполнить произвольную индексацию, тогда это более сложный вопрос, потому что вы используете не столько слайсы, сколько список индексов. Функция _get_single_element вызывается 100000 раз, потому что она просто перебирает итераторы (i и j), которые вы ей передали. Срез будет A [30: 60,10] или что-то подобное.

Кроме того, я бы использовал csr_matrix(eye(n,n)) для создания той же матрицы, которую вы создали с помощью итераторов, просто для простоты.

Обновление:

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

  • Почему это не реализовано в чистом C?

Ответ прост: никто не дошел до этого. Из того, что я видел, в области скудных матричных модулей предстоит еще много работы. Одна часть, которая реализована в C, - это преобразования между различными форматами разреженных матриц.

  • Кто-нибудь знает способ обойти это ограничение и ускорить приведенный выше код?

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

Тем не менее, было все еще быстрее использовать lil_matrix (хотя медленнее инициализировать матрицу), но мне пришлось делать доступ с пониманием списка, потому что матрицы lil не настроены для типа индексации, который вы делаете. Использование понимания списка для csr_matrix все еще оставляет способ матрицы lil далеко впереди. В конечном счете, матрица lil быстрее для доступа к случайным элементам, потому что она не сжата.

Использование lil_matrix в его первоначальном виде занимает примерно пятую часть времени кода, который вы перечислили выше. Если я возьму некоторые связанные проверки и вызову метод _get1 () lil_matrix напрямую, я смогу сократить время еще на 7% по сравнению с первоначальным временем. Для ясности это ускорение с 3,4-3,8 секунды до 0,261 секунды.

Наконец, я попытался создать свою собственную функцию, которая напрямую обращается к данным матрицы lil и избегает повторных вызовов функций. Время для этого было около 0,136 секунды. Это не использовало преимущества сортировки данных, что является еще одной потенциальной оптимизацией (в частности, если вы обращаетесь ко множеству элементов, которые находятся в одних и тех же строках).

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

  • Должен ли я использовать другой тип разреженной матрицы?

Хорошо, я предлагаю матрицу lil, если вы хотите получить доступ ко всем элементам, но все зависит от того, что вам нужно сделать. Вам также нужно умножить матрицы, например? Просто помните, что переключение между матрицами может, по крайней мере, иногда (в определенных обстоятельствах) быть довольно быстрым, поэтому не исключайте возможность перехода на другой формат матрицы для выполнения различных операций.

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

0 голосов
/ 09 марта 2010

Я думаю, что _get_single_element вызывается только тогда, когда используется dtype по умолчанию для объекта. Вы пытались предоставить dtype, например csr_matrix((data, (i,j)), dtype=int32)

...