删除gsl矩阵的列

Question

在下面的代码中，我从贝塔分布double-valued 中抽取了一个一维mu矩阵(mu)，然后用这个mu矩阵来采样一个二维二进制值矩阵(Z)，使用the Bernoulli分布。但是，有时我会得到一些列，这些列是由zero值填充的。如何有效地编写一个函数，该函数丢弃由零占据的列，同时在矩阵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+typed内存视图允许您以"C速度“执行此操作。

使用GSL的好理由是：( 1)与使用GSL或2)的现有C库交互，如果GSL提供一个您无法从Numpy获得的有用函数(在这种情况下，您只需使用该function+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]