c - 删除 gsl 矩阵的列

Question

在下面的代码中，我采样了一维 double-valued来自 beta 分布的矩阵( mu )，然后我使用这个 mu使用伯努利分布对二维二进制值矩阵 ( Z ) 进行采样。然而有时我最终会得到一些由 zero 填充的列。值(value)观。如何有效地编写一个函数，丢弃被零占用的列，并丢弃矩阵 mu 中的相应值不会在功能中造成任何冲突 func其中 gsl 矩阵 Z , mu是首先定义的？

由于我需要经常调用删除代码中零值列的函数，我想知道定义动态 gsl 矩阵并分配特定大小但能够 resize 的最佳方法是什么？一遍又一遍的矩阵没有出现任何错误？这是到目前为止我的代码:

import numpy as np
cimport numpy as np
cimport cython
from cython.view cimport array as cvarray
from libc.stdlib cimport malloc, free
from libc.math cimport log, exp
from cython_gsl cimport *
import ctypes
from libc.string cimport memset

cdef extern from "stdlib.h":         
    long rand()                    
    void srand(int seed)                    

cdef extern from "time.h":
    ctypedef long time_t
    time_t time(time_t*)

cdef gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937)
srand(time(NULL))
@cython.cdivision(True)     
@cython.boundscheck(False)
@cython.wraparound(False)     
cdef void initPar( double* alpha, int* K, int* N, gsl_matrix* mu, gsl_matrix *Z ):
     cdef:
          int i, j
          int seed = rand()
     #generate random seed
     gsl_rng_set(r, seed)
     cdef double prev_mu
     for i from 0 <= i < K[0]:
         if i >= 1:
            prev_mu *= gsl_ran_beta(r, alpha[0], 1)                
         else:
            prev_mu = gsl_ran_beta(r, alpha[0], 1)

         gsl_matrix_set(mu, i, 0, prev_mu)        

     cdef double mi
     for i from 0 <= i < K[0]:
         mi= gsl_matrix_get(mu, i, 0)
         for j from 0 <= j < N[0]: 
             gsl_matrix_set(Z, j, i, gsl_ran_bernoulli(r, mi) )

     return

def func(np.ndarray[ndim=2, dtype=np.float64_t] inputs, int LFeatures = 5, double alpha):
    cdef:
        int *K_init= &LFeatures
        int N    = inputs.shape[0]
        int D    = inputs.shape[1]
        gsl_matrix *mu = gsl_matrix_alloc(K_init[0], 1)
        gsl_matrix *Z  = gsl_matrix_alloc(N, K_init[0])
        int i, j
        gsl_matrix *X    = gsl_matrix_alloc(N, D)
    for i from 0 <= i < N:
         for j from 0 <= j < D:
             gsl_matrix_set(X, i, j, inputs[i,j])

    gsl_matrix_set_zero(mu)
    gsl_matrix_set_zero(Z)
    initPar(&alpha, K_init, &N, mu, Z )

谢谢!

Answer 1

这个答案并不是您直接想要的，因为它只是 Numpy。然而，我确实认为使用 GSL 会使您的代码变得更加复杂，但没有任何实际好处。

您声称使用 GSL 是因为您相信它具有“C 速度”，而 Numpy 则不然 - 这根本不是事实。 Numpy 最终是用 C 编写的，因此如果您可以一次对整个数组执行操作，那么速度会非常快。 Numpy 数组的迭代速度较慢，但​​ Cython+类型化内存 View 允许您以“C 速度”执行此操作。

使用 GSL 的充分理由是 1) 与使用 GSL 的现有 C 库交互，或者 2) 如果 GSL 提供了一个您无法从 Numpy 获得的有用函数(在这种情况下，您可以只使用该函数+numpy )。我认为这两者都不适用于此。

<小时/>

您可以轻松地将初始化函数简化为简短的纯 Python+Numpy，该函数将在 Numpy 内优化的 C 循环中执行:

def initPar(alpha, N, K):
    mu = np.cumprod(np.random.beta(alpha,1,size=(K,)))

    # Bernoulli distribution is a special case of the binomial distribution
    # use array broadcasting to get the shape of Z
    Z = np.random.binomial(np.ones((N,1),dtype=np.int32),
                           mu[np.newaxis,:])
    return mu, Z

删除所有全为零的列(以及 mu 的相应元素)的问题在 numpy 中又是直接且快速的:

mask = np.any(Z,axis=0) # true where any element of the set - i.e. the columns to keep
Z = Z[:,mask]
mu = mu[mask]