删除 gsl 矩阵的列

Question

在以下代码中，我从 beta 分布 中采样一维 double-valued 矩阵 (mu)，然后我使用此 mu 矩阵采样二维二进制值矩阵 (Z)，使用 伯努利分布。但是，有时我最终会得到一些仅由 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 memoryviews 允许您以“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]