Векторизация индексной операции для матрицы scipy.sparse

Question

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

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?Кто-нибудь знает способ обойти это ограничение и ускорить выполнение приведенного выше кода?Должен ли я использовать другой тип разреженной матрицы?

Answer 1

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

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

Обновлять:

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

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

Ответ прост:никто до этого не дошел.Судя по тому, что я видел, предстоит еще много работы в области разреженных матричных модулей Scipy.Одна часть, реализованная в 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 секунды.При этом не использовалась сортировка данных, что является еще одной потенциальной оптимизацией (особенно, если вы получаете доступ к множеству элементов, находящихся в одних и тех же строках).

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

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

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

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

Answer 2

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