¿Llamando productos de puntos y operaciones de álgebra lineal en Cython?

Estoy tratando de usar productos de puntos, inversión de matrices y otras operaciones básicas de álgebra lineal que están disponibles en muchos de Cython. Funciones como numpy.linalg.inv (inversión), numpy.dot (producto de puntos), Xt (transposición de matriz / matriz). Hay una gran sobrecarga para llamar numpy.* Desde las funciones de Cython y el rest de la función está escrita en Cython, así que me gustaría evitar esto.

Si asumo que los usuarios tienen un numpy instalado, ¿hay alguna manera de hacer algo como:

 #include "numpy/npy_math.h" 

Como extern , y llamar a estas funciones? ¿O también puede llamar a BLAS directamente (o lo que sea que el número llame para estas operaciones centrales)?

Para dar un ejemplo, imagine que tiene una función en Cython que hace muchas cosas y al final necesita hacer un cálculo que involucre productos de puntos y matrices inversas:

 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) 

¿Cómo se puede hacer esto? Si ya hay una biblioteca que implementa esto en Cython, también puedo usar eso, pero no he encontrado nada. Incluso si esos procedimientos están menos optimizados que BLAS directamente, no tener la sobrecarga de llamar numpy módulo de Python desde Cython hará que las cosas en general sean más rápidas.

Ejemplo de funciones que me gustaría llamar:

  • producto np.dot ( np.dot )
  • inversión de matriz ( np.linalg.inv )
  • multiplicación de matrices
  • teniendo transposición (equivalente a xT en numpy)
  • Función gammaln (como scipy.gammaln equivalente, que debería estar disponible en C)

Me doy cuenta de lo que dice en la lista de correo numpy ( https://groups.google.com/forum/?fromgroups=#!topic/cython-users/XZjMVSIQnTE ) que si llama a estas funciones en matrices grandes, no tiene sentido hacerlo desde Cython, ya que llamarlo desde numpy solo resultará en la mayor parte del tiempo empleado en el código C optimizado que numpy llama. Sin embargo, en mi caso, tengo muchas llamadas a estas operaciones de álgebra lineal en matrices pequeñas ; en ese caso, la sobrecarga introducida al pasar repetidamente de Cython a numpy y de nuevo a Cython superará con creces el tiempo empleado en calcular la operación desde BLAS. Por lo tanto, me gustaría mantener todo en el nivel C / Cython para estas operaciones simples y no pasar por Python.

Preferiría no pasar por GSL, ya que eso agrega otra dependencia y ya que no está claro si GSL se mantiene activamente. Ya que asumo que los usuarios del código ya tienen instalado scipy / numpy, puedo asumir con seguridad que tienen todo el código C asociado que acompaña a estas bibliotecas, así que solo quiero poder acceder a ese código y llamarlo de Cython.

edición : encontré una biblioteca que contiene BLAS en Cython ( https://github.com/tokyo/tokyo ) que está cerca pero no es lo que estoy buscando. Me gustaría llamar a las funciones numpy / scipy C directamente (supongo que el usuario las tiene instaladas).

Llamar a BLAS con Scipy es “bastante” sencillo, aquí hay un ejemplo para llamar a DGEMM para calcular la multiplicación de matrices: https://gist.github.com/pv/5437087 Tenga en cuenta que BLAS y LAPACK esperan que todas las matrices sean contiguas a Fortran (módulo los parámetros lda / b / c), por lo tanto, order="F" y double[::1,:] que se requieren para un correcto funcionamiento.

La computación inversa se puede hacer de manera similar aplicando la función LAPACK dgesv en la matriz de identidad. Para la firma, ver aquí . Todo esto requiere descender a una encoding de bajo nivel, debe asignar arreglos de trabajo temporales, etc., etc. Sin embargo, estos pueden encapsularse en sus propias funciones de conveniencia, o simplemente reutilizar el código de tokyo reemplazando las funciones lib_* . con punteros de función obtenidos de Scipy de la manera anterior.

Si usa la syntax de la vista de memoria de Cython ( double[::1,:] ), la transposición es la misma xT como de costumbre. Alternativamente, puede calcular la transposición escribiendo una función propia que intercambie elementos de la matriz a través de la diagonal. Numpy en realidad no contiene esta operación, xT solo cambia los pasos de la matriz y no mueve los datos.

Probablemente sería posible volver a escribir el módulo de tokyo para usar el BLAS / LAPACK exportado por Scipy y scipy.linalg en scipy.linalg , para que pueda hacerlo from scipy.linalg.blas cimport dgemm . Se aceptan solicitudes de extracción si alguien quiere hacerlo.


Como puede ver, todo se reduce a pasar punteros de función alrededor. Como se mencionó anteriormente, Cython de hecho proporciona su propio protocolo para intercambiar punteros de función. Para un ejemplo, considere from scipy.spatial import qhull; print(qhull.__pyx_capi__) from scipy.spatial import qhull; print(qhull.__pyx_capi__) — se puede acceder a estas funciones from scipy.spatial.qhull cimport XXXX en Cython (aunque son privadas, así que no lo hagas).

Sin embargo, en el presente, scipy.special no ofrece este C-API. Sin embargo, sería bastante sencillo proporcionarlo, dado que el módulo de interfaz en scipy.special está escrito en Cython.

No creo que en este momento exista una forma sensata y portátil de acceder a la función haciendo el trabajo pesado para el gamln , (aunque podría husmear alrededor del objeto UFunc, pero eso no es una solución sensata :), por lo que en este momento es Probablemente sea mejor tomar la parte relevante del código fuente de scipy.special y agruparlo con su proyecto, o usar, por ejemplo, GSL.

Quizás la forma más fácil si acepta usar GSL sería usar esta interfaz GSL-> cython https://github.com/twiecki/CythonGSL y llamar a BLAS desde allí (vea el ejemplo https://github.com/twiecki /CythonGSL/blob/master/examples/blas2.pyx ). También debe hacerse cargo de los pedidos Fortran vs C. No hay muchas características nuevas de GSL, pero puede asumir con seguridad que se mantiene activamente. El CythonGSL es más completo en comparación con tokio; por ejemplo, presenta productos de matriz simétrica que están ausentes en el número.

Como acabo de encontrar el mismo problema y escribí algunas funciones adicionales, las incluiré aquí en caso de que alguien más las encuentre útiles. Codifico alguna multiplicación de matrices, y también llamo a las funciones LAPACK para la inversión matricial, determinante y descomposición de Cholesky. Pero deberías considerar intentar hacer cosas de álgebra lineal fuera de cualquier bucle, si tienes alguna, como lo hago aquí . Y, por cierto, la función determinante aquí no funciona del todo si tienes sugerencias. Además, tenga en cuenta que no hago ninguna comprobación para ver si las entradas son 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) nxn array to invert B : memoryview (numpy array) nxn array to use within the function, function will modify this matrix in place to become the inverse of A work : memoryview (numpy array) nxn 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) nxn array to compute determinant of work : memoryview (numpy array) nxn 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) nxn matrix to compute cholesky decomposition B : memoryview (numpy array) nxn 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 (nxm) and B (mxl) Parameters ---------- A : memoryview (numpy array) nxm left matrix B : memoryview (numpy array) mxr right matrix out : memoryview (numpy array) nxr 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