c++ - 结合python和c++，或者cython，优化一个函数；最大似然示例；对c++的了解很少

Question

我懂 Python，但我不懂 C++。我正在尝试最大化需要很长时间才能评估的功能。我相信一个好的工作流程是用 C++ 编写评估函数的函数，并将此函数与 scipy.optim.minimize 一起使用以找到最佳值。例如，假设我正在最大化可能性。

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)

可能有更有效的方法用纯 Python 编写 llf 函数，并且肯定有加速优化例程的方法(例如，通过选择适合问题的特定优化器，或者通过提供衍生品)，但这不是这个问题的重点。我选择这个特定示例是因为我有一个循环(我正在使用所有数据，包括某些列缺少值的行)来评估可能性，这在纯 python 中需要花费大量时间，尤其是当我的样本量增加时.

如何用 C++ 编写似然函数并将其与 Python 最小化例程相结合？请记住，我没有使用 C++ 的经验，但愿意学习。然而，许多可用于此的资源似乎假定 C++ 知识，请参阅 Extending Python例如。我正在为了解 Python 但完全不了解 C++ 以及将 Python 与 C++ 结合的方法的人寻找资源。 编辑:也许使用我的示例或有关结合 Python 和 C++ 可能获得的 yield 的信息的示例会很有用。

Answer 1

按照建议，我尝试了 Cython 解决方案。由于我以前从未使用过 Cython，因此我将完成用于实现 Cython 解决方案的步骤。

首先，我安装了 Cython。然后我编写了一个名为 faSTLlf.pyx 的文件，其中包含以下 Cython 代码:

#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)

然后我创建了一个 setup.py 文件，其中包含以下内容:

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

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

接下来，我在终端中使用以下命令编译了 Cython 代码。

$ python3 setup.py build_ext --inplace

最后，我比较了旧的纯 Python 实现(稍微修改为使用数组而不是数据帧)和 Cython 实现之间的结果。

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)

我得到了以下结果:

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

这使得最大似然优化变得更加可行，尤其是当我的问题变得更加复杂时。唯一的问题是我需要找到在 Cython 中编写 llf 函数所需的数学和统计函数，或者我需要编写自己的函数，就像我对上面的普通 pdf 所做的那样。

如对我的实现提出任何意见，我们将不胜感激。