appel aux produits dot et aux opérations d'algèbre linéaire à Cython?
j'essaie d'utiliser les produits dot, l'inversion de matrice et d'autres opérations de base de l'algèbre linéaire qui sont disponibles en numpy de Cython. Fonctions comme numpy.linalg.inv (inversion), numpy.dot (produit de point), X.t (transposition de matrice/tableau). Il y a un grand overhead pour appeler numpy.* à partir des fonctions Cython et le reste de la fonction est écrite en Cython, donc j'aimerais éviter cela.
si je suppose que les utilisateurs ont numpy installé, est-il une façon de faire quelque chose comme:
#include "numpy/npy_math.h"
comme extern , et appelez ces fonctions? Ou bien appeler directement BLAS (ou quoi que ce soit que numpy appelle pour ces opérations centrales)?
pour donner un exemple, imaginez que vous avez une fonction à Cython qui fait beaucoup de choses et qui a finalement besoin de faire un calcul impliquant des produits de points et des inverses de matrice:
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)
comment faire? Si il y a un bibliothèque qui met en œuvre ces À Cython déjà, je peux également utiliser cela, mais n'ai rien trouvé. Même si ces procédures sont moins optimisées que BLAS directement, ne pas avoir la charge d'appeler numpy module Python de Cython rendra les choses globalement plus rapides.
exemples de fonctions que j'aimerais appeler:
- produit scalaire (
np.dot) - inversion de matrice (
np.linalg.inv) - la multiplication de matrice
- prendre "transpose" (équivalent de
x.Tdans numpy) - fonction gammaln (comme
scipy.gammalnéquivalent, qui devrait être disponible en C)
je me rends compte comme il est dit sur la liste de diffusion numpy ( https://groups.google.com/forum/?fromgroups=#!topic / cython-users / XZjMVSIQnTE ) que si vous appelez ces fonctions sur de grandes matrices, il n'y a pas de point dans le faire à partir de Cython, car l'appeler à partir de numpy ne fera que passer la majorité du temps dans le code C optimisé que numpy appelle. Cependant, dans mon cas, j'ai beaucoup d'appels à ces opérations d'algèbre linéaire sur de petites matrices -- dans ce cas, le overhead introduit par aller à plusieurs reprises de Cython retour à num PY et retour à Cython sera de loin plus que le temps passé réellement le calcul de l'opération de BLAS. Par conséquent, je voudrais garder tout à la Niveau C / Cython pour ces opérations simples et ne pas passer par python.
je préférerais ne pas passer par GSL, car cela ajoute une autre dépendance et parce qu'il n'est pas clair si GSL est maintenu activement. Puisque je suppose que les utilisateurs du code ont déjà scipy / num PY installé, je peux en toute sécurité supposer qu'ils ont tout le code C associé qui va avec ces bibliothèques, donc je veux juste être en mesure de puiser dans ce code et L'appeler de Cython.
edit : j'ai trouvé une bibliothèque qui enveloppe BLAS à Cython ( https://github.com/tokyo/tokyo ) qui est proche mais pas ce que je cherche. J'aimerais appeler les fonctions C de numpy/scipy directement (je suppose que l'utilisateur les a installées.)
3 réponses
appeler BLAS avec Scipy est "assez" simple, voici un exemple pour appeler DGEMM pour calculer la multiplication matricielle: https://gist.github.com/pv/5437087 notez que BLAS et LAPACK s'attendent à ce que tous les tableaux soient contigus (modulo les paramètres lda/b/c), donc order="F" et double[::1,:] qui sont nécessaires pour un fonctionnement correct.
les inverses de calcul peuvent être faites de la même façon en appliquant la fonction LAPACK dgesv sur la matrice d'identité. Pour la signature, voir ici . Tout cela nécessite de descendre à un niveau de codage assez bas, vous devez allouer des tableaux de travail temporaires vous-même etc etc --- cependant, ceux-ci peuvent être encapsulés dans vos propres fonctions de commodité, ou tout simplement réutiliser le code de tokyo en remplaçant les lib_* fonctions avec des pointeurs de fonction obtenus de Scipy de la manière ci-dessus.
si vous utilisez Cython's memoryview la syntaxe ( double[::1,:] ) que vous transposez est la même x.T que d'habitude. Alternativement, vous pouvez calculer la transposition en écrivant une fonction de votre propre qui échange des éléments du tableau à travers la diagonale. Num PY ne contient pas réellement cette opération, x.T ne fait que changer les pas du tableau et ne déplace pas les données.
il serait probablement possible de réécrire le module tokyo pour utiliser le BLAS / LAPACK exporté par Scipy et le regrouper en scipy.linalg , pour que tu puisses juste faire from scipy.linalg.blas cimport dgemm . Pull request sont acceptés si quelqu'un veut descendre.
comme vous pouvez le voir, tout se résume à passer des indicateurs de fonction autour. Comme mentionné ci-dessus, Cython fournit en fait son propre protocole pour échanger des pointeurs de fonction. Pour un exemple, considérons from scipy.spatial import qhull; print(qhull.__pyx_capi__) - - - ces fonctions peuvent être accessibles via from scipy.spatial.qhull cimport XXXX à Cython (elles sont privées bien que, donc ne faites pas cela).
cependant, à l'heure actuelle, scipy.special n'offre pas ce C-API. Il serait en fait assez simple de le fournir, étant donné que le module d'interface en scipy.spécial est écrit en Cython.
Je ne pense pas qu'il y ait pour le moment de manière saine et portable pour accéder à la fonction de levage lourd pour gamln , (bien que vous pourriez fouiner autour de L'objet UFunc, mais ce n'est pas une solution saine :), donc à l'heure actuelle, il est probablement préférable de simplement saisir la partie pertinente du code source de scipy.spécial et le regrouper avec votre projet, ou utiliser par exemple GSL.
peut-être que la façon la plus simple si vous acceptez d'utiliser la GSL serait d'utiliser cette interface GSL - >cython https://github.com/twiecki/CythonGSL et appelez BLAS de là (voir l'exemple https://github.com/twiecki/CythonGSL/blob/master/examples/blas2.pyx ). Il doit également prendre soin de la commande Fortran vs C. Il n'y a pas beaucoup de nouvelles fonctionnalités GSL, mais vous pouvez en toute sécurité supposer qu'il est activement maintenu. Le CythonGSL est plus complet comparé à tokyo, par exemple, il présente des produits à matrice symétrique qui sont absents de num PY.
comme je viens de rencontrer le même problème, et a écrit quelques fonctions supplémentaires, je vais les inclure ici au cas où quelqu'un d'autre les trouve utiles. Je code vers le haut d'une certaine multiplication de matrice, et aussi appellent des fonctions de LAPACK pour l'inversion de matrice, le déterminant et la décomposition cholesky. Mais vous devriez envisager d'essayer de faire des trucs d'algèbre linéaire en dehors des boucles, si vous en avez, comme je fais ici . Et d'ailleurs, la fonction déterminante ici ne fonctionne pas vraiment si vous avez des suggestions. Aussi, s'il vous plaît noter que je ne fais aucune vérification pour voir si les entrées sont conformes.
from scipy.linalg.cython_lapack cimport dgetri, dgetrf, dpotrf
cpdef void double[:, ::1] inv_c(double[:, ::1] A, double[:, ::1] B,
double[:, ::1] work, double[::1] ipiv):
'''invert float type square matrix A
Parameters
----------
A : memoryview (numpy array)
n x n array to invert
B : memoryview (numpy array)
n x n array to use within the function, function
will modify this matrix in place to become the inverse of A
work : memoryview (numpy array)
n x n array to use within the function
ipiv : memoryview (numpy array)
length n array to use within function
'''
cdef int n = A.shape[0], info, lwork
B[...] = A
dgetrf(&n, &n, &B[0, 0], &n, &ipiv[0], &info)
dgetri(&n, &B[0,0], &n, &ipiv[0], &work[0,0], &lwork, &info)
cpdef double det_c(double[:, ::1] A, double[:, ::1] work, double[::1] ipiv):
'''obtain determinant of float type square matrix A
Notes
-----
As is, this function is not yet computing the sign of the determinant
correctly, help!
Parameters
----------
A : memoryview (numpy array)
n x n array to compute determinant of
work : memoryview (numpy array)
n x n array to use within function
ipiv : memoryview (numpy array)
length n vector use within function
Returns
-------
detval : float
determinant of matrix A
'''
cdef int n = A.shape[0], info
work[...] = A
dgetrf(&n, &n, &work[0,0], &n, &ipiv[0], &info)
cdef double detval = 1.
cdef int j
for j in range(n):
if j != ipiv[j]:
detval = -detval*work[j, j]
else:
detval = detval*work[j, j]
return detval
cdef void chol_c(double[:, ::1] A, double[:, ::1] B):
'''cholesky factorization of real symmetric positive definite float matrix A
Parameters
----------
A : memoryview (numpy array)
n x n matrix to compute cholesky decomposition
B : memoryview (numpy array)
n x n matrix to use within function, will be modified
in place to become cholesky decomposition of A. works
similar to np.linalg.cholesky
'''
cdef int n = A.shape[0], info
cdef char uplo = 'U'
B[...] = A
dpotrf(&uplo, &n, &B[0,0], &n, &info)
cdef int i, j
for i in range(n):
for j in range(n):
if j > i:
B[i, j] = 0
cpdef void dotmm_c(double[:, :] A, double[:, :] B, double[:, :] out):
'''matrix multiply matrices A (n x m) and B (m x l)
Parameters
----------
A : memoryview (numpy array)
n x m left matrix
B : memoryview (numpy array)
m x r right matrix
out : memoryview (numpy array)
n x r output matrix
'''
cdef Py_ssize_t i, j, k
cdef double s
cdef Py_ssize_t n = A.shape[0], m = A.shape[1]
cdef Py_ssize_t l = B.shape[0], r = B.shape[1]
for i in range(n):
for j in range(r):
s = 0
for k in range(m):
s += A[i, k]*B[k, j]
out[i, j] = s