sistema - transpuesta de una matriz en python numpy

llamando a productos de punto y operaciones de álgebra lineal en Cython? (3)

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

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

#include "numpy/npy_math.h"

como un extern , y llamar a estas funciones? O, alternativamente, llame a BLAS directamente (o lo que sea que numpy llama 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 punto e inversos de matriz:

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 estos 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 numpy Python de Cython aún hará que todo sea más rápido.

Ejemplo de funciones a las que me gustaría llamar:

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

Me doy cuenta, como dice en la lista de correo numpy ( https://groups.google.com/forum/?fromgroups=#!topic/cython-users/XZjMVSIQnTE ) que si llamas 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 dedicado 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 regreso a Cython superará con creces el tiempo invertido en calcular la operación de 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 no está claro si GSL se mantiene activamente. Como supongo que los usuarios del código ya tienen instalado scipy / numpy, puedo suponer con seguridad que tienen todos los códigos C asociados que van con estas bibliotecas, así que solo quiero poder acceder a ese código y llamarlo de Cython.

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

Como acabo de encontrar el mismo problema y escribí algunas funciones adicionales, las incluiré aquí en caso de que alguien las encuentre útiles. Codigo alguna multiplicación de matriz, y también llamo funciones LAPACK para inversión de matriz, descomposición determinante y cholesky. Pero deberías considerar hacer cosas de álgebra lineal fuera de cualquier bucle, si es que tienes alguno, como hago here . Y, por cierto, la función determinante aquí no está funcionando si tienes sugerencias. Además, tenga en cuenta que no realizo 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) 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

Llamar a BLAS incluido 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 Fortran-contiguas (módulo los parámetros lda / b / c), por lo tanto order="F" y double[::1,:] que son necesarios para el correcto funcionamiento.

Las inversas informáticas se pueden hacer de manera similar aplicando la función LAPACK dgesv en la matriz de identidad. Para la firma, mira here . Todo esto requiere bajar a una codificación de bajo nivel, necesita asignar arreglos temporales de trabajo, etc., etc., sin embargo, estos pueden ser encapsulados 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 en la forma anterior.

Si usas la sintaxis de vista de memoria de Cython ( double[::1,:] ), la transposición es la misma xT que 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 sea posible reescribir el módulo de tokyo para usar el BLAS / LAPACK exportado por Scipy y scipy.linalg en scipy.linalg , para que pueda hacer simplemente from scipy.linalg.blas cimport dgemm . Se aceptan solicitudes de extracción si alguien quiere llegar a ella.

Como puede ver, todo se reduce a pasar indicadores de función a su alrededor. Como se mencionó anteriormente, Cython de hecho proporciona su propio protocolo para intercambiar punteros de función. Por ejemplo, considere from scipy.spatial import qhull; print(qhull.__pyx_capi__) from scipy.spatial import qhull; print(qhull.__pyx_capi__) --- se puede acceder a esas funciones a través from scipy.spatial.qhull cimport XXXX en Cython (aunque son privadas, no lo haga).

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

No creo que haya en este momento una forma sana y portátil para acceder a la función que hace el trabajo pesado para gamln (aunque se podía husmear en el objeto UFunc, pero esa no es una solución sensata :), así que por el 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 el 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 debería encargarse del pedido de Fortran vs C. No hay muchas características nuevas de GSL, pero puede suponer de forma segura que se mantiene activamente. CythonGSL es más completo en comparación con tokyo; por ejemplo, presenta productos de matriz simétrica que están ausentes en numpy.