Combine python and c++, or cython, to optimize a function; maximum likelihood example; little knowledge of c++

Question

I know Python, but I don't know C++. I'm trying to maximize a function that takes a long time to evaluate. I believe a good workflow would be to write the function that evaluates the function in C++ and use this function with scipy.optim.minimize to find the optimum. As an example, suppose I am maximizing a likelihood.

import pandas as pd
import numpy as np
from scipy.optimize import minimize
from scipy.stats import norm

# simulating data
means = np.array([10, 20, 30])
cov = np.diag([1, 4, 10])

N = 1000

df = pd.DataFrame(np.random.multivariate_normal(mean=means, cov=cov, size=N),
    columns=['a', 'b', 'c'])
df[np.random.choice([True, False], size=(N, 3), p=[0.3, 0.7])] = np.nan

# a function to print parameters used in likelihood function
def print_params(params):
    print('Means: {}'.format(params[:3]))
    print('Variances: {}'.format(np.exp(params[3:])**2))

# defining likelihood
def llf(params):
    logll = 0
    for i in df.index:
        for j,col in enumerate(['a', 'b', 'c']):
            if not np.isnan(df.loc[i, col]):
                m = params[j]
                sd = np.exp(params[j+3])
                logll += np.log(norm.pdf(df.loc[i, col], loc=m, scale=sd))

    print_params(params)
    return -logll


opt = minimize(llf, x0=np.array([0, 0, 0, 1, 1, 1]), options={'maxiter':30})
print_params(opt.x)

There may be more efficient ways to write the llf function in pure Python, and there are definitely ways to speed up the optimization routine (e.g., by choosing a specific optimizer suited to the problem, or by supply derivatives), but this is not the focus of this question. I chose this particular example because I have a loop (I'm using all of the data, including rows where some columns have missing values) to evaluate the likelihood, which takes a lot of time in pure python, especially if my sample size grows.

How can I write the likelihood function in C++ and combine that with the Python minimize routine? Keep in mind, I have no experience with C++ but am willing to learn. However, a lot of the resources available for this seem to presume C++ knowledge, see Extending Python for example. I'm looking for resources for someone who knows Python, but is completely ignorant of C++ and methods for combining Python with C++. EDIT: Perhaps an example of how to do this using my example or information about the likely gains from combining Python and C++ would be useful.

Answer 1

As suggested, I've tried a Cython solution. As I've never used Cython before, I'll be complete in the steps I used to implement the Cython solution.

First, I installed Cython. Then I wrote a file called fastllf.pyx which contained the following Cython code:

#cython: boundscheck=False, wraparound=False, nonecheck=False

from libc.math cimport exp, sqrt, pi, log, isnan

cdef double SQ_PI = sqrt(2*pi)


cdef double norm_pdf(double x, double loc, double scale):
    return (1/(SQ_PI*scale))*exp(-(0.5)*((x - loc)**2)/(scale**2))

cdef double llf_c(double[:, :] X, double[:] params):

    cdef double logll = 0
    cdef int N = X.shape[0]
    cdef int K = X.shape[1]
    cdef int i, j
    cdef double m, sd

    for i in range(N):
        for j in range(K):
            if not isnan(X[i, j]):
                m = params[j]
                sd = exp(params[j+K])

                logll += log(norm_pdf(X[i, j], m, sd))
    return -logll

def llf(double[:, :] X, double[:] params):
    return llf_c(X, params)

Then I created a setup.py file which included the following:

from distutils.core import setup
from Cython.Build import cythonize

setup(name="fastllf", ext_modules=cythonize('fastllf.pyx'))

Next I compiled the Cython code using the following command in the terminal.

$ python3 setup.py build_ext --inplace

Finally, I compared the results between my old, pure Python implementation (slightly modified to use arrays instead of dataframes) and the Cython implementation.

import numpy as np
from scipy.stats import norm
import time
from fastllf import llf as cython_llf

# simulating data
means = np.array([10, 20, 30])
cov = np.diag([1, 4, 10])

N = 100000
np.random.seed(10)

X = np.random.multivariate_normal(mean=means, cov=cov, size=N)
X[np.random.choice([True, False], size=(N, 3), p=[0.3, 0.7])] = np.nan

def norm_pdf(x, loc, scale):
    return (1/(np.sqrt(2*np.pi)*scale))*np.exp(-(0.5)*((x-loc)**2)/(scale**2))

def llf(X, params):

    logll = 0
    N = X.shape[0]
    K = X.shape[1]

    for i in range(N):
        for j in range(K):
            if not np.isnan(X[i, j]):
                m = params[j]
                sd = np.exp(params[j+K])

                logll += np.log(norm_pdf(X[i, j], loc=m, scale=sd))    
    return -logll    

def timeit(fun, *args):
    start = time.time()
    rslt = fun(*args)
    end = time.time()
    print(rslt)
    print(end - start)

params = np.array([1.,1,1,1,1,1])
timeit(llf, X, params)
timeit(cython_llf, X, params)

And I obtained the following results:

Python Value: 6570173.7597125955
Python Time:  1.9558300971984863 seconds
Cython Value: 6570173.7597125955
Cython Time:  0.016242027282714844 seconds

This makes optimization by maximum likelihood much more feasible, especially when my problem gets more complicated. The only problem is that I need to find the math and statistical functions I need to write an llf function in Cython or I need to write my own, as I did for the normal pdf above.

Any comments on my implementation would be appreciated.