时间戳

首页 > 解决方案 > 结合python和c++，或者cython，优化一个函数；最大似然示例；c++的小知识

c++ - 结合python和c++，或者cython，优化一个函数；最大似然示例；c++的小知识

问题描述

我知道 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++ 知识，例如，请参阅扩展 Python 。我正在为了解 Python 但完全不了解 C++ 以及将 Python 与 C++ 结合的方法的人寻找资源。编辑：也许使用我的示例或有关结合 Python 和 C++ 的可能收益的信息来说明如何做到这一点的示例会很有用。

标签： c++python-3.xscipycython

解决方案

如建议的那样，我尝试了 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

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

对我的实施的任何评论将不胜感激。

推荐阅读