cython 中没有 python 个对象的高效矩阵运算
Efficient matrix operations in cython with no python objects
我正在编写一个函数，我想使用 Cython 尽可能多地将其转换为 C。为此，我需要使用线性代数运算。这是我的功能。 编辑：我学到的教训是尝试处理循环之外的线性代数，这在很大程度上是我能够做到的。否则，求助于包装 LAPACK/BLAS 或编写我自己的函数。
import numpy as np
from scipy.stats import multivariate_normal as mv
import itertools
def llf(data, rho, mu, sigma, A, V, n):
'''evaluate likelihood by guass-hermite quadrature
Parameters
----------
data : array
N x J matrix, columns are measurements
rho : array
length L vector of weights for mixture of normals
mu : array
L x K vector of means of mixture of normals
sigma : array
K x L*K matrix of variance matrices for mixture of normals
A : array
J x (K + 1) matrix of loadings
V : array
J x J variance matrix of measurement errors
n : int
number of sample points for quadrature
'''
N = data.shape[0]
L, K = mu.shape
# getting weights and sample points for approximating integral
v, w = np.polynomial.hermite.hermgauss(n)
totllf = 0
for i in range(N):
M_i = data[i, :]
totllf_i = 0
for l in range(L):
rho_l = rho[l]
sigma_l = sigma[:, K*l:K*(l+1)]
mu_l = mu[l, :]
chol_l = np.linalg.cholesky(sigma_l)
for ix in itertools.product(*(list(range(n)) for k in range(K))):
wt = np.prod(w[list(ix)])
pt = np.sqrt(2)*chol_l.dot(v[list(ix)]) + mu_l
totllf_i += wt*rho_l*mv.pdf(M_i, A[:, 0] + A[:, 1:].dot(pt), V)
totllf += np.log(totllf_i)
return totllf
为了完成这个，我需要有矩阵乘法、转置、行列式、矩阵求逆和 cholesky 分解的函数。我看过 some posts 使用 BLAS 函数，但我真的不清楚如何使用这些函数。
编辑 2018 年 4 月 29 日
按照建议，我采用了内存视图方法并在循环之前初始化了所有内容。我的新函数写成
def llf_c(double[:, ::1] data, double[::1] rho, double[:, ::1] mu,
double[:, ::1] sigma, double[:, ::1] A, double[:, ::1] V, int n):
'''evaluate likelihood by guass-hermite quadrature
Parameters
----------
data : array
N x J matrix, columns are measurements
rho : array
length L vector of weights for mixture of normals
mu : array
L x K vector of means of mixture of normals
sigma : array
K x L*K matrix of variance matrices for mixture of normals
A : array
J x (K + 1) matrix of loadings
V : array
J x J variance matrix of measurement errors
n : int
number of sample points for quadrature
'''
cdef Py_ssize_t N = data.shape[0], J = data.shape[1], L = mu.shape[0], K = mu.shape[1]
# initializing indexing variables
cdef Py_ssize_t i, l, j, k
# getting weights and sample points for approximating integral
v_a, w_a = np.polynomial.hermite.hermgauss(n)
cdef double[::1] v = v_a
cdef double[::1] w = w_a
cdef double[::1] v_ix = np.zeros(K, dtype=np.float)
# initializing memory views for cholesky decomposition of sigma matrices
sigma_chol_a = np.zeros((K, L*K), dtype=np.float)
for l in range(L):
sigma_chol_a[:, K*l:K*(l+1)] = np.linalg.cholesky(sigma[:, K*l:K*(l+1)])
cdef double[:, ::1] sigma_chol = sigma_chol_a
# intializing V inverse and determinant
cdef double[:, ::1] V_inv = np.linalg.inv(V)
cdef double V_det = np.linalg.det(V)
# initializing memoryviews for work matrices
cdef double[::1] work_K = np.zeros(K, dtype=np.float)
cdef double[::1] work_J = np.zeros(J, dtype=np.float)
# initializing memoryview for quadrature points
cdef double[::1] pt = np.zeros(K, dtype=np.float)
# initializing memorview for means for multivariate normal
cdef double[::1] loc = np.zeros(J, dtype=np.float)
# initializing values for loop
cdef double[::1] totllf = np.zeros(N, dtype=np.float)
cdef double wt, pdf_init = 1./sqrt(((2*pi)**J)*V_det)
cdef int[:, ::1] ix_vals = np.vstack(itertools.product(*(list(range(n)) for k in range(K)))).astype(np.int32)
cdef Py_ssize_t ix_len = ix_vals.shape[0]
for ix_row in range(ix_len):
ix = ix_vals[ix_row]
# weights and points for quadrature
wt = 1.
for k in range(K):
wt *= w[ix[k]]
v_ix[k] = v[ix[k]]
for l in range(L):
# change of variables
dotmv_c(sigma_chol[:, K*l:K*(l+1)], v_ix, work_K)
for k in range(K):
pt[k] = sqrt(2)*work_K[k]
addvv_c(pt, mu[l, :], pt)
for i in range(N):
# generating demeaned vector for multivariate normal pdf
dotmv_c(A[:, 1:], pt, work_J)
addvv_c(A[:, 0], work_J, work_J)
for j in range(J):
loc[j] = data[i, j] - work_J[j]
# performing matrix products in exponential
# print(wt, rho[l], np.asarray(work_J))
dotvm_c(loc, V_inv, work_J)
totllf[i] += wt*rho[l]*pdf_init*exp(-0.5*dotvv_c(work_J, loc))
return np.log(np.asarray(totllf)).sum()
dotvm_c、dotmv_c和addvv_c是执行向量与矩阵的矩阵乘法、矩阵与向量的矩阵乘法以及两个向量的元素相加的函数。我也在 Cython 中编写了这些，但为简洁起见不包括在内。我不再像在使用 numpy 循环之前处理所有其他线性代数那样包装任何 LAPACK 函数。我还有几个问题。为什么我的循环中仍然有黄色？ （见下面的简介）。我认为现在一切都应该在 C 中。另外，如果您对新的实施有任何其他建议，请告诉我。
例如，在第 221 行，我在编译时收到此消息："Index should be typed for more efficient access." 但我认为我输入了索引 k。另外，由于 addvv_c 显示为黄色，我将在下面向您展示它的定义。
cpdef void addvv_c(double[:] a, double[:] b, double[:] out):
'''add two vectors elementwise
'''
cdef Py_ssize_t i, n = a.shape[0]
for i in range(n):
out[i] = a[i] + b[i]
关于你优化的Cython/BLAS函数的几点小意见：
ipiv_a = np.zeros(n).astype(np.int32)
cdef int[::1] ipiv = ipiv_a
可以有两个简单的改进：它不必通过临时变量，你可以直接创建一个类型 np.int32 的数组，而不是创建一个不同类型的数组然后转换：
cdef int[::1] ipiv = np.zeros(n,dtype=np.int32)
同样，在这两个函数中，您可以通过执行
以更少的步骤初始化 B
cdef double[:, ::1] B = A.copy()
而不是创建零数组并复制
第二个（更重要的）变化是 to use C arrays 用于临时变量，例如 Fortran 工作区。我仍然将 return 值之类的东西保留为 numpy 数组，因为引用计数和将它们发送回 Python 的能力 非常 有用。
cdef double* work = <double*>malloc(n*n*sizeof(double))
try:
# rest of function
finally:
free(work)
您需要从 libc.stdlib 导入 malloc 和 free。 try: ... finally: 确保内存被正确释放。不要太过头了 - 例如，如果在何处释放 C 数组并不明显，那么只需使用 numpy。
要查看的最后一个选项是没有 return 值但修改输入：
cdef void inv_c(double[:,::1] A, double[:,::1] B):
# check that A and B are the right size, then just write into B
# ...
这样做的好处是，如果您需要在具有相同大小输入的循环中调用它，那么您只需为整个循环进行一次分配。您可以将其扩展为也包含工作数组，尽管这可能会稍微复杂一些。
我正在编写一个函数，我想使用 Cython 尽可能多地将其转换为 C。为此，我需要使用线性代数运算。这是我的功能。 编辑：我学到的教训是尝试处理循环之外的线性代数，这在很大程度上是我能够做到的。否则，求助于包装 LAPACK/BLAS 或编写我自己的函数。
import numpy as np
from scipy.stats import multivariate_normal as mv
import itertools
def llf(data, rho, mu, sigma, A, V, n):
'''evaluate likelihood by guass-hermite quadrature
Parameters
----------
data : array
N x J matrix, columns are measurements
rho : array
length L vector of weights for mixture of normals
mu : array
L x K vector of means of mixture of normals
sigma : array
K x L*K matrix of variance matrices for mixture of normals
A : array
J x (K + 1) matrix of loadings
V : array
J x J variance matrix of measurement errors
n : int
number of sample points for quadrature
'''
N = data.shape[0]
L, K = mu.shape
# getting weights and sample points for approximating integral
v, w = np.polynomial.hermite.hermgauss(n)
totllf = 0
for i in range(N):
M_i = data[i, :]
totllf_i = 0
for l in range(L):
rho_l = rho[l]
sigma_l = sigma[:, K*l:K*(l+1)]
mu_l = mu[l, :]
chol_l = np.linalg.cholesky(sigma_l)
for ix in itertools.product(*(list(range(n)) for k in range(K))):
wt = np.prod(w[list(ix)])
pt = np.sqrt(2)*chol_l.dot(v[list(ix)]) + mu_l
totllf_i += wt*rho_l*mv.pdf(M_i, A[:, 0] + A[:, 1:].dot(pt), V)
totllf += np.log(totllf_i)
return totllf
为了完成这个，我需要有矩阵乘法、转置、行列式、矩阵求逆和 cholesky 分解的函数。我看过 some posts 使用 BLAS 函数，但我真的不清楚如何使用这些函数。
编辑 2018 年 4 月 29 日
按照建议，我采用了内存视图方法并在循环之前初始化了所有内容。我的新函数写成
def llf_c(double[:, ::1] data, double[::1] rho, double[:, ::1] mu,
double[:, ::1] sigma, double[:, ::1] A, double[:, ::1] V, int n):
'''evaluate likelihood by guass-hermite quadrature
Parameters
----------
data : array
N x J matrix, columns are measurements
rho : array
length L vector of weights for mixture of normals
mu : array
L x K vector of means of mixture of normals
sigma : array
K x L*K matrix of variance matrices for mixture of normals
A : array
J x (K + 1) matrix of loadings
V : array
J x J variance matrix of measurement errors
n : int
number of sample points for quadrature
'''
cdef Py_ssize_t N = data.shape[0], J = data.shape[1], L = mu.shape[0], K = mu.shape[1]
# initializing indexing variables
cdef Py_ssize_t i, l, j, k
# getting weights and sample points for approximating integral
v_a, w_a = np.polynomial.hermite.hermgauss(n)
cdef double[::1] v = v_a
cdef double[::1] w = w_a
cdef double[::1] v_ix = np.zeros(K, dtype=np.float)
# initializing memory views for cholesky decomposition of sigma matrices
sigma_chol_a = np.zeros((K, L*K), dtype=np.float)
for l in range(L):
sigma_chol_a[:, K*l:K*(l+1)] = np.linalg.cholesky(sigma[:, K*l:K*(l+1)])
cdef double[:, ::1] sigma_chol = sigma_chol_a
# intializing V inverse and determinant
cdef double[:, ::1] V_inv = np.linalg.inv(V)
cdef double V_det = np.linalg.det(V)
# initializing memoryviews for work matrices
cdef double[::1] work_K = np.zeros(K, dtype=np.float)
cdef double[::1] work_J = np.zeros(J, dtype=np.float)
# initializing memoryview for quadrature points
cdef double[::1] pt = np.zeros(K, dtype=np.float)
# initializing memorview for means for multivariate normal
cdef double[::1] loc = np.zeros(J, dtype=np.float)
# initializing values for loop
cdef double[::1] totllf = np.zeros(N, dtype=np.float)
cdef double wt, pdf_init = 1./sqrt(((2*pi)**J)*V_det)
cdef int[:, ::1] ix_vals = np.vstack(itertools.product(*(list(range(n)) for k in range(K)))).astype(np.int32)
cdef Py_ssize_t ix_len = ix_vals.shape[0]
for ix_row in range(ix_len):
ix = ix_vals[ix_row]
# weights and points for quadrature
wt = 1.
for k in range(K):
wt *= w[ix[k]]
v_ix[k] = v[ix[k]]
for l in range(L):
# change of variables
dotmv_c(sigma_chol[:, K*l:K*(l+1)], v_ix, work_K)
for k in range(K):
pt[k] = sqrt(2)*work_K[k]
addvv_c(pt, mu[l, :], pt)
for i in range(N):
# generating demeaned vector for multivariate normal pdf
dotmv_c(A[:, 1:], pt, work_J)
addvv_c(A[:, 0], work_J, work_J)
for j in range(J):
loc[j] = data[i, j] - work_J[j]
# performing matrix products in exponential
# print(wt, rho[l], np.asarray(work_J))
dotvm_c(loc, V_inv, work_J)
totllf[i] += wt*rho[l]*pdf_init*exp(-0.5*dotvv_c(work_J, loc))
return np.log(np.asarray(totllf)).sum()
dotvm_c、dotmv_c和addvv_c是执行向量与矩阵的矩阵乘法、矩阵与向量的矩阵乘法以及两个向量的元素相加的函数。我也在 Cython 中编写了这些，但为简洁起见不包括在内。我不再像在使用 numpy 循环之前处理所有其他线性代数那样包装任何 LAPACK 函数。我还有几个问题。为什么我的循环中仍然有黄色？ （见下面的简介）。我认为现在一切都应该在 C 中。另外，如果您对新的实施有任何其他建议，请告诉我。
例如，在第 221 行，我在编译时收到此消息："Index should be typed for more efficient access." 但我认为我输入了索引 k。另外，由于 addvv_c 显示为黄色，我将在下面向您展示它的定义。
cpdef void addvv_c(double[:] a, double[:] b, double[:] out):
'''add two vectors elementwise
'''
cdef Py_ssize_t i, n = a.shape[0]
for i in range(n):
out[i] = a[i] + b[i]
关于你优化的Cython/BLAS函数的几点小意见：
ipiv_a = np.zeros(n).astype(np.int32)
cdef int[::1] ipiv = ipiv_a
可以有两个简单的改进：它不必通过临时变量，你可以直接创建一个类型 np.int32 的数组，而不是创建一个不同类型的数组然后转换：
cdef int[::1] ipiv = np.zeros(n,dtype=np.int32)
同样，在这两个函数中，您可以通过执行
以更少的步骤初始化B
cdef double[:, ::1] B = A.copy()
而不是创建零数组并复制
第二个（更重要的）变化是 to use C arrays 用于临时变量，例如 Fortran 工作区。我仍然将 return 值之类的东西保留为 numpy 数组，因为引用计数和将它们发送回 Python 的能力 非常 有用。
cdef double* work = <double*>malloc(n*n*sizeof(double))
try:
# rest of function
finally:
free(work)
您需要从 libc.stdlib 导入 malloc 和 free。 try: ... finally: 确保内存被正确释放。不要太过头了 - 例如，如果在何处释放 C 数组并不明显，那么只需使用 numpy。
要查看的最后一个选项是没有 return 值但修改输入：
cdef void inv_c(double[:,::1] A, double[:,::1] B):
# check that A and B are the right size, then just write into B
# ...
这样做的好处是，如果您需要在具有相同大小输入的循环中调用它，那么您只需为整个循环进行一次分配。您可以将其扩展为也包含工作数组，尽管这可能会稍微复杂一些。