Efficient Cython for large matrix creation/manipulation
Question: I have been trying to speed up a section of code that creates and manipulates
a very large matrix of data (approx. 15,000 x 15,000; type double). For now, I
don't think the size of the matrix is that important because I do not see
speedup even for a small 10 x 10 matrix (in fact, the compiled cython code is
slower than pure python for small matrices, whereas the time is nearly
identical between cython and python for the large matrices). Please be patient
with me, as I have only been coding python for a week (newly converted from
Matlab) and I am only a humble chemical engineer.

The goal of the code is to take a 1D array (length L) as input, for a example:

    
    
    [ 16.66  16.85  16.93  16.98  17.08  17.03  17.09  16.76  16.67  16.72]
    

And produce a matrix (height L, width L-1) as output:

    
    
    [[ 16.66  16.85  16.93  16.98  17.08  17.03  17.09  16.76  16.67]
     [ 16.85  16.93  16.98  17.08  17.03  17.09  16.76  16.67  16.72]
     [ 16.93  16.98  17.08  17.03  17.09  16.76  16.67  16.72   0.  ]
     [ 16.98  17.08  17.03  17.09  16.76  16.67  16.72   0.     0.  ]
     [ 17.08  17.03  17.09  16.76  16.67  16.72   0.     0.     0.  ]
     [ 17.03  17.09  16.76  16.67  16.72   0.     0.     0.     0.  ]
     [ 17.09  16.76  16.67  16.72   0.     0.     0.     0.     0.  ]
     [ 16.76  16.67  16.72   0.     0.     0.     0.     0.     0.  ]
     [ 16.67  16.72   0.     0.     0.     0.     0.     0.     0.  ]
     [ 16.72   0.     0.     0.     0.     0.     0.     0.     0.  ]]
    

I hope it is clear from the example above and from the code below what I am
trying to achieve. The algorithm needs to scale to very large matrices, which
it currently does without error, it's just slow!

Here is my cython code:

    
    
    from scipy.sparse import spdiags
    import numpy as np
    cimport numpy as np
    cimport cython
    
    @cython.boundscheck(False)
    @cython.wraparound(False)
    def sfmat(np.ndarray[double, ndim=1] data):
        cdef int h = data.shape[0]   
        cdef np.ndarray[double, ndim=2] m = np.zeros([h, h-1])
        m = np.flipud(spdiags(np.tril(np.tile(data,[h-1,1]).T,0),range(1-h,1), h, h-1).todense())
        return m
    

I also tried more verbose code, which may be more clear to read:

    
    
    from scipy.sparse import spdiags
    import numpy as np
    cimport numpy as np
    cimport cython
    
    DTYPE = np.float
    ctypedef np.float_t DTYPE_t
    
    @cython.boundscheck(False)
    @cython.wraparound(False)
    def sfmat(np.ndarray[DTYPE_t, ndim=1] data):
        assert data.dtype == DTYPE
        cdef int h = data.shape[0]   
        cdef np.ndarray[DTYPE_t, ndim=2] m = np.zeros([h, h-1], dtype=DTYPE)
        cdef np.ndarray[DTYPE_t, ndim=2] s1 = np.zeros([h, h-1], dtype=DTYPE)
        cdef np.ndarray[DTYPE_t, ndim=2] s2 = np.zeros([h, h-1], dtype=DTYPE)
        cdef np.ndarray[DTYPE_t, ndim=2] s3 = np.zeros([h, h-1], dtype=DTYPE)
    
        s1 = np.tile(data,[h-1,1]).T
        s2 = np.tril(s1,0)
        s3 = spdiags(s2,range(1-h,1), h, h-1).todense()
        m = np.flipud(s3)
        return m
    

Any help with the cython implemenation would be very much appreciated. If
there is any other way to speed up this algorithm, that would help too. Thank
you for any help!

Because I am new to this, here are more details, which may or may not be
preventing me from speeding this up. I am running 64 bit Windows 7 Pro, and
compiling the cython code successfuly using the Windows SDK C/C++ compiler. (I
followed the directions on github
[here](https://github.com/cython/cython/wiki/64BitCythonExtensionsOnWindows)
with success). The simple "hello world" cython examples compile fine and run
fine in 64 bit mode, and the code above also compiles and runs with no errors.
For manipulation of the entire 15,000 x 15,000 matrix, 64bit architecture is
required, or at least I believe so, because running the code after compiling
for 32bit resulted in a memory error. For this question, please assume that
breaking up the matrix into smaller chunks is not possible. Please let me know
if there is any other information required to answer this question.

Cheers, scientistR

**UPDATE**

I thought that avoiding for loops would be the best approach, however, spdiags
is the main bottleneck. Thus, a new algorithm works better (4x improvement on
my computer):

    
    
    import numpy as np
    cimport numpy as np
    cimport cython
    
    @cython.boundscheck(False)
    @cython.wraparound(False)
    def sfmat(np.ndarray[double, ndim=1] data):
         cdef int i
         cdef np.ndarray[double, ndim=2] m = np.zeros([data.shape[0], data.shape[0]-1])
         for i in range(data.shape[0]-1):
             m[:,i] = np.roll(data,-i);
         return m
    

But Cython does not offer any improvement over pure Python. Please help. As
the commentators have pointed out, there may not be a way to improve this,
besides a more optimizes algorithm, but I am hopeful. Thanks! Also, is there a
faster algorithm, cython or python?


Answer: This might be a bit of an old question, but no question should be left
unanswered `:)`. I was able to speed up your Cython code by about 8x by using
straightforward for-loops (which are actually fast in Cython), for array size
of 7000.. Note by the way that your implementation using `np.roll` does not
produce the array that you want (!), but I used that function to compare
timings with.

_Edited code to use Typed Memoryviews and`np.empty` instead of `np.zeros`_

    
    
    def sfmat(double[:] data):
         cdef int n = data.shape[0]
         cdef np.ndarray[double, ndim=2] out = np.empty((n, n-1))
         cdef double [:, :] out_v = out  # "typed memoryview"
    
         cdef int i, j
         for i in range(n-1):
            out_v[0, i] = data[i]
    
         for i in range(1, n):
            for j in range(n-i):
                out_v[i, j] = data[i+j]
            for j in range(n-i, n-1):
                out_v[i, j] = 0.
         return out
    

Unfortunately, the Cython effort is only ~1.2x faster than running the
following code in a regular Python session:

    
    
    def sfmat(data):
        n = len(data)
        out = np.empty((n, n-1))
        out[0, :] = data[:n-1]
        for i in xrange(1, n):
            out[i, :n-i] = data[i:]
            out[i, n-i:] = 0
        return out
    

However as already discussed in the comments, blowing up your original fairly
small matrix in this way is probably not the most efficient way to tackle your
actual, overall problem anyway. If all you wanted to do initially was avoid
using for-loops, in Cython there's simply no need to do that!

