Binomial iterated expectation in Cython

Question

I am brand new to Cython. How to convert the Python function called Values below to Cython? With factors=2 and i=60 this takes 2.8 secs on my big Linux box. The goal is sub 1 sec with factors=2 and i=360.

Here's the code. Thanks!

import numpy as np
import itertools

class Numeraire:
    def __init__(self, rate):
        self.rate = rate
    def __call__(self, timenext, time, state):
        return np.exp(-self.rate*(timenext - time))

def Values(values, i1, i0=0, numeraire=Numeraire(0.)):
    factors=len(values.shape)
    norm=0.5**factors
    for i in np.arange(i1-1, i0-1, -1):
        for j in itertools.product(np.arange(i+1), repeat=factors):
            value = 0.
            for k in itertools.product(np.arange(2), repeat=factors):
                value += values[tuple(np.array(j) + np.array(k))]
            values[j] = value*norm*numeraire(i+1, i, j)
    return values

factors = 2
i = 60
values = np.ones([i+1]*factors)
Values(values, i, numeraire=Numeraire(0.05/12))
print values[(0,)*factors], np.exp(-0.05/12*i)

You should probably try to use standard NumPy vectorized operations before bringing in Cython. In particular, you'll pretty quickly be bottlenecked by `numeraire` unless you start computing that in a vectorized manner or move it into Cython as well. — user2357112, May 21 '15 at 21:08
I think you'd get a lot more help if you explained what this code does, including the effects on a small array. With your sample, `values` is smallest in the upper left corner, and increases in 'circles' to 1 on the right and bottom. Why? — hpaulj, May 22 '15 at 21:40
Because it's attempting to implement a binomial lattice, used in option pricing. Here's a reference: http://www.investopedia.com/features/eso/eso3.asp — Steve Schulist, May 26 '15 at 21:41

score 2 · Answer 1 · edited Aug 05 '20 at 07:38

Before using Cython, you should optimize your code with Numpy. Here, vectorizing the third and second inner for loops, yields a x40 speed-up,

In [1]: import numpy as np
   ...: import itertools
   ...: 
   ...:  # define Numaire and Values functions from the question above
   ...:
   ...: def Values2(values, i1, i0=0, numeraire=Numeraire(0.)):
   ...:     factors=len(values.shape)
   ...:     norm=0.5**factors
   ...:     k = np.array(list(itertools.product(np.arange(2), repeat=factors)))
   ...:     for i in np.arange(i1-1, i0-1, -1):
   ...:         j = np.array(list(itertools.product(np.arange(i+1), repeat=factors)))
   ...:         mask_all = j[:,:,np.newaxis] + k.T[np.newaxis, :, :]
   ...:         mask_x, mask_y = np.swapaxes(mask_all, 2, 1).reshape(-1, 2).T
   ...:      
   ...:         values_tmp = values[mask_x, mask_y].reshape((j.shape[0], k.shape[0]))
   ...:         values_tmp = values_tmp.sum(axis=1)
   ...:         values[j[:,0], j[:,1]] = values_tmp*norm*numeraire(i+1, i, j)
   ...:     return values
   ...:
   ...: factors = 2
   ...: i = 60
   ...: values = lambda : np.ones([i+1]*factors)
   ...: print values()[(0,)*factors], np.exp(-0.05/12*i)
   ...:
   ...: res = Values(values(), i, numeraire=Numeraire(0.05/12))
   ...: res2 = Values2(values(), i, numeraire=Numeraire(0.05/12))
   ...: np.testing.assert_allclose(res, res2)
   ...:
   ...: %timeit  Values(values(), i, numeraire=Numeraire(0.05/12))
   ...: %timeit  Values2(values(), i, numeraire=Numeraire(0.05/12))
   ...:
1.0 0.778800783071
1 loops, best of 3: 1.26 s per loop
10 loops, best of 3: 31.8 ms per loop

The next step would be to replace the line,

j = np.array(list(itertools.product(np.arange(i+1), repeat=factors)

with it's Numpy equivalent, taken from this answer (not very pretty),

def itertools_product_numpy(some_list, some_length):
    return some_list[np.rollaxis(
          np.indices((len(some_list),) * some_length), 0, some_length + 1)
    .reshape(-1, some_length)]

k = itertools_product_numpy(np.arange(i+1), factors)

this result in an overall x160 speed up and the code runs in 1.2 second on my laptop for i=360 and factors = 2.

In this last version, I don't think that you will get much speed up, if you port it to Cython, since there is just one loop remaining and it has only ~360 iterations. Rather, some fine-tuned Python/Numpy optimizations should be performed to get a further speed increase.

Alternatively, you can try applying Cython to your original implementation. However because it is based on itertools.product, which is slow when called repeatedly in a loop, Cython will not help there.

Wow, thank you for your wonderful answer! On my machine your solution runs in 1.4 sec! Sorry to be greedy but, can your function be modified to allow the factors = 1 case, for example? — Steve Schulist, May 26 '15 at 21:34
@SteveSchulist Aww, yes sorry, I have not realized to that `values` has `factors` columns, this was indeed done for the case `factors=2` (i.e. 2D array) and would need to be generalized to ND arrays (or at least 1D). Unfortunately, I don't think I could look into this week though... — rth, May 26 '15 at 22:54

Steve Schulist · Answer 2 · 2015-06-02T17:13:30.367

Here's my latest answer (no Cython!), which runs in 125 msec for the factor=2, i=360 case.

import numpy as np
import itertools

slices = (slice(None, -1, None), slice(1, None, None))

def Expectation(values, numeraire, i, i0=0):
    def Values(values, i):
        factors = values.ndim
        expect = np.zeros((i,)*factors)
        for j in itertools.product(slices, repeat=factors):
            expect += values[j]
        return expect*0.5**factors*numeraire(i, i-1)
    return reduce(Values, range(i, i0, -1), values)

class Numeraire:
    def __init__(self, factors, rate=0):
        self.factors = factors
        self.rate = rate
    def __call__(self, timenext, time):
        return np.full((time+1,)*factors, np.exp(-self.rate*(timenext - time)))

factors = 2
i = 360
values, numeraire = np.ones((i+1,)*factors), Numeraire(factors, 0.05/12)
%timeit Expectation(values, numeraire, i)
Expectation(values, numeraire, i)[(0,)*factors], np.exp(-0.05/12*i)

Binomial iterated expectation in Cython

2 Answers2