Вызывающие точечные произведения и операции линейной алгебры в Китоне?

Я пытаюсь использовать точечные продукты, инверсию матрицы и другие основные операции линейной алгебры, доступные в numpy от Cython. Такие функции, как numpy.linalg.inv (инверсия), numpy.dot (точечный продукт), X.t (транспонирование матрицы/массива). Там большие накладные расходы для вызова numpy.* из функций Cython, а остальная часть функции написана в Cython, поэтому я хотел бы избежать этого.

Если я предполагаю, что у пользователей установлен numpy, есть ли способ сделать что-то вроде:

#include "numpy/npy_math.h"

как extern, и вызывать эти функции? Или, альтернативно, напрямую вызывайте BLAS (или что это такое, что numpy вызывает эти основные операции)?

Чтобы привести пример, представьте, что у вас есть функция в Cython, которая делает много вещей и, в конце концов, должна сделать вычисление, включающее точечные продукты и обратные матрицы:

cdef myfunc(...):
  # ... do many things faster than Python could
  # ...
  # compute one value using dot products and inv
  # without using 
  #   import numpy as np 
  #   np.*
  val = gammaln(sum(v)) - sum(gammaln(v)) + dot((v - 1).T, log(x).T)

как это можно сделать? Если там уже есть библиотека, которая реализует их в Cython, я также могу ее использовать, но ничего не нашел. Даже если эти процедуры менее оптимизированы, чем BLAS напрямую, не имея накладных расходов на вызов numpy модуля Python от Cython, все равно будет делать все быстрее.

Примеры функций, которые я бы назвал:

dot product (np.dot)
инверсия матрицы (np.linalg.inv)
матричное умножение
принимает транспонирование (эквивалентно X.t в numpy)
Функция gammaln (например, эквивалент scipy.gammaln, который должен быть доступен в C)

Я понимаю, как он говорит в списке рассылки numpy (https://groups.google.com/forum/?fromgroups=#!topic/cython-users/XZjMVSIQnTE), что, если вы вызываете эти функции на больших матрицах, нет смысла делать это из Cython, так как вызов его из numpy приведет к большей части времени, затраченного на оптимизированный код C, который вызывает numpy. Однако в моем случае у меня много вызовов к этим операциям линейной алгебры на малых матрицах - в этом случае накладные расходы, вводимые путем повторного перехода от Cython обратно к numpy и обратно в Cython, намного перевесят время, затраченное на вычисление операции из BLAS. Поэтому я хотел бы сохранить все на уровне C/Cython для этих простых операций и не проходить через python.

Я бы предпочел не проходить через GSL, поскольку это добавляет еще одну зависимость, и поскольку неясно, активно ли GSL. Поскольку я предполагаю, что у пользователей кода уже установлен scipy/numpy, я могу с уверенностью предположить, что у них есть все связанные с ним C-коды, которые идут вместе с этими библиотеками, поэтому я просто хочу иметь возможность использовать этот код и называть его от Китона.

edit. Я нашел библиотеку, которая обертывает BLAS в Cython (https://github.com/tokyo/tokyo), который близок, но не то, что я ищу для. Я хотел бы вызвать функции numpy/scipy C напрямую (я предполагаю, что пользователь установил эти настройки.)

Ответ 1

Вызов BLAS в комплекте с Scipy "довольно" прямолинейный, здесь один пример вызова DGEMM для вычисления умножения матрицы: https://gist.github.com/pv/5437087 Обратите внимание, что BLAS и LAPACK ожидают, что все массивы (по модулю параметров lda/b/c), следовательно order="F" и double[::1,:], которые необходимы для правильного функционирования.

Вычислительные инверсии можно аналогичным образом выполнить, применив функцию LAPACK dgesv к единице. Для подписи см. здесь. Все это требует перехода к довольно низкоуровневому кодированию, вам нужно самостоятельно распределять временные рабочие массивы и т.д. --- однако они могут быть инкапсулированы в ваши собственные удобные функции или просто повторно использовать код из tokyo, заменив lib_* с указателями функций, полученными из Scipy указанным выше способом.

Если вы используете синтаксис вида Cython memoryview (double[::1,:]), то вы транспонируете то же самое x.T, как обычно. В качестве альтернативы вы можете вычислить транспонирование, написав собственную функцию, которая меняет элементы массива по диагонали. Numpy фактически не содержит эту операцию, x.T изменяет только шаги массива и не перемещает данные вокруг.

Вероятно, можно будет переписать модуль tokyo, чтобы использовать BLAS/LAPACK, экспортированный Scipy, и связать его в scipy.linalg, чтобы вы могли просто сделать from scipy.linalg.blas cimport dgemm. Pull request принимаются, если кто-то хочет перейти к нему.

Как вы можете видеть, все это сводится к передаче указателей функций. Как упоминалось выше, Cython фактически предоставляет собственный протокол для обмена указателями на функции. Например, рассмотрим from scipy.spatial import qhull; print(qhull.__pyx_capi__) - эти функции можно было бы получить через from scipy.spatial.qhull cimport XXXX в Cython (они частные, но не делайте этого).

Однако в настоящее время scipy.special не предлагает этот C-API. Однако на самом деле было бы довольно просто предоставить его, учитывая, что интерфейсный модуль в scipy.special написан в Cython.

Я не думаю, что на данный момент есть разумный и переносимый способ доступа к функции, выполняющей тяжелый подъем для gamln (хотя вы могли бы snoop вокруг объекта UFunc, но это не разумное решение:), поэтому на данный момент, вероятно, лучше всего просто захватить соответствующую часть исходного кода из scipy.special и связать его с вашим проектом или использовать, например, GSL.

Ответ 1

Ответ 2