calling dot products and linear algebra operations in Cython?

Question

I'm trying to use dot products, matrix inversion and other basic linear algebra operations that are available in numpy from Cython. Functions like numpy.linalg.inv (inversion), numpy.dot (dot product), X.t (transpose of matrix/array). There's a large overhead to calling numpy.* from Cython functions and the rest of the function is written in Cython, so I'd like to avoid this.

If I assume users have numpy installed, is there a way to do something like:

#include "numpy/npy_math.h"

as an extern, and call these functions? Or alternatively call BLAS directly (or whatever it is that numpy calls for these core operations)?

To give an example, imagine you have a function in Cython that does many things and in the end needs to make a computation involving dot products and matrix inverses:

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)

how can this be done? If there's a library that implements these in Cython already, I can also use that, but have not found anything. Even if those procedures are less optimized than BLAS directly, not having the overhead of calling numpy Python module from Cython will still make things overall faster.

Example functions I'd like to call:

• dot product (np.dot)
• matrix inversion (np.linalg.inv)
• matrix multiplication
• taking transpose (equivalent of x.T in numpy)
• gammaln function (like scipy.gammaln equivalent, which should be available in C)

I realize as it says on numpy mailing list (https://groups.google.com/forum/?fromgroups=#!topic/cython-users/XZjMVSIQnTE) that if you call these functions on large matrices, there is no point in doing it from Cython, since calling it from numpy will just result in the majority of the time spent in the optimized C code that numpy calls. However, in my case, I have many calls to these linear algebra operations on small matrices -- in that case, the overhead introduced by repeatedly going from Cython back to numpy and back to Cython will far outweigh the time spent actually computing the operation from BLAS. Therefore, I'd like to keep everything at the C/Cython level for these simple operations and not go through python.

I'd prefer not to go through GSL, since that adds another dependency and since it's unclear if GSL is actively maintained. Since I'm assuming users of the code already have scipy/numpy installed, I can safely assume that they have all the associated C code that goes along with these libraries, so I just want to be able to tap into that code and call it from Cython.

edit: I found a library that wraps BLAS in Cython (https://github.com/tokyo/tokyo) which is close but not what I'm looking for. I'd like to call the numpy/scipy C functions directly (I'm assuming the user has these installed.)

Answer 1

Calling BLAS bundled with Scipy is "fairly" straightforward, here's one example for calling DGEMM to compute matrix multiplication: https://gist.github.com/pv/5437087 Note that BLAS and LAPACK expect all arrays to be Fortran-contiguous (modulo the lda/b/c parameters), hence order="F" and double[::1,:] which are required for correct functioning.

Computing inverses can be similarly done by applying the LAPACK function dgesv on the identity matrix. For the signature, see here. All this requires dropping down to rather low-level coding, you need to allocate temporary work arrays yourself etc etc. --- however these can be encapsulated into your own convenience functions, or just reuse the code from tokyo by replacing the lib_* functions with function pointers obtained from Scipy in the above way.

If you use Cython's memoryview syntax (double[::1,:]) you transpose is the same x.T as usual. Alternatively, you can compute the transpose by writing a function of your own that swaps elements of the array across the diagonal. Numpy doesn't actually contain this operation, x.T only changes the strides of the array and doesn't move the data around.

It would probably be possible to rewrite the tokyo module to use the BLAS/LAPACK exported by Scipy and bundle it in scipy.linalg, so that you could just do from scipy.linalg.blas cimport dgemm. Pull requests are accepted if someone wants to get down to it.

---

As you can see, it all boils down to passing function pointers around. As alluded to above, Cython does in fact provide its own protocol for exchanging function pointers. For an example, consider from scipy.spatial import qhull; print(qhull.__pyx_capi__) --- those functions could be accessed via from scipy.spatial.qhull cimport XXXX in Cython (they're private though, so don't do that).

However, at the present, scipy.special does not offer this C-API. It would however in fact be quite simple to provide it, given that the interface module in scipy.special is written in Cython.

I don't think there is at the moment any sane and portable way to access the function doing the heavy lifting for gamln, (although you could snoop around the UFunc object, but that's not a sane solution :), so at the moment it's probably best to just grab the relevant part of source code from scipy.special and bundle it with your project, or use e.g. GSL.

Answer 2

Perhaps the easiest way if you do accept using the GSL would be to use this GSL->cython interface https://github.com/twiecki/CythonGSL and call BLAS from there (see the example https://github.com/twiecki/CythonGSL/blob/master/examples/blas2.pyx). It should also take care of the Fortran vs C ordering. There aren't many new GSL features, but you can safely assume it is actively maintained. The CythonGSL is more complete compared to tokyo; e.g., it features symmetric-matrix products that are absent in numpy.

Answer 3

As I've just encountered the same problem, and wrote some additional functions, I'll include them here in case someone else finds them useful. I code up some matrix multiplication, and also call LAPACK functions for matrix inversion, determinant and cholesky decomposition. But you should consider trying to do linear algebra stuff outside any loops, if you have any, like I do here. And by the way, the determinant function here isn't quite working if you have suggestions. Also, please note that I don't do any checking to see if inputs are conformable.

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