2015-10-21 43 views
1

我正在寻找一种向量化方法来应用将二维数组返回给二维数组的每一行并产生的函数一个三维数组。NumPy:应用将矩阵返回矩阵的每一行的函数的一般向量化方法

更具体地说,我有一个函数,它需要一个长度为p的向量并返回一个二维数组(m乘n)。以下是我的函数的程式化版本:

import numpy as np 
def test_func(x, m, n): 
    # this function is just an example and does not do anything useful. 
    # but, the dimensions of input and output is what I want to convey. 
    np.random.seed(x.sum()) 
    return np.random.randint(5, size=(m, n)) 

我有一件T由P 2维输入数据:

t = 5 
p = 6 
input_data = np.arange(t*p).reshape(t, p) 
input_data 
Out[403]: 
array([[ 0, 1, 2, 3, 4, 5], 
     [ 6, 7, 8, 9, 10, 11], 
     [12, 13, 14, 15, 16, 17], 
     [18, 19, 20, 21, 22, 23], 
     [24, 25, 26, 27, 28, 29]]) 

我想申请test_func到input_data的每一行。由于test_func返回一个矩阵,我期望创建一个3维(t乘m乘n)的数组。我可以用下面的代码产生我想要的结果:

output_data = np.array([test_func(x, m=3, n=2) for x in input_data]) 
output_data 
Out[405]: 
array([[[0, 4], 
     [0, 4], 
     [3, 3], 
     [1, 0]], 

     [[1, 0], 
     [1, 0], 
     [4, 1], 
     [2, 4]], 

     [[3, 3], 
     [3, 0], 
     [1, 4], 
     [0, 2]], 

     [[2, 4], 
     [2, 1], 
     [3, 2], 
     [3, 1]], 

     [[3, 4], 
     [4, 3], 
     [0, 3], 
     [3, 0]]]) 

但是,这段代码似乎并不是最优的代码。它有一个明确的降低速度,并使用不必要地分配额外内存的中间列表。所以,我喜欢找到一个矢量化的解决方案。我最好的猜测是以下代码,但它不起作用。

output = np.apply_along_axis(test_func, m=3, n=2, axis=1, arr=input_data) 
Traceback (most recent call last): 

    File "<ipython-input-406-5bef44da348f>", line 1, in <module> 
    output = np.apply_along_axis(test_func, m=3, n=2, axis=1, arr=input_data) 

    File "C:\Anaconda\lib\site-packages\numpy\lib\shape_base.py", line 117, in apply_along_axis 
    outarr[tuple(i.tolist())] = res 

ValueError: could not broadcast input array from shape (3,2) into shape (3) 

请您建议一种有效的方法来解决这个问题。

UPDATE

下面是我想申请的实际功能。它执行多维古典缩放。问题的目的不是优化函数的内部运作,而是找到一个向量化函数apply的泛化方法。但是,本着全面披露的精神,我将实际功能放在这里。请注意,此功能仅当p == M *(M-1)/ 2

def mds_classical_scaling(v, m, n):  

    # create a symmetric distance matrix from the elements in vector v 
    D = np.zeros((m, m)) 
    D[np.triu_indices(4, k=1)] = v 
    D = (D + D.T) 

    # Transform the symmetric matrix 
    A = -0.5 * (D**2) 
    # Create centering matrix  
    H = np.eye(m) - np.ones((m, m))/m 
    # Doubly center A and store in B 
    B = H*A*H 

    # B should be positive definite otherwise the function 
    # would not work. 
    mu, V = eig(B) 

    #index of largest eigen values 
    ndx = (-mu).argsort() 

    # calculate the point configuration from largest eigen values 
    # and corresponding eigen vectors 
    Mu1 = diag(mu[ndx][:n]) 
    V1 = V[:, ndx[:n]] 
    X = V1*sqrt(Mu1)  

    return X 

任何性能提升,我从量化得到的是微不足道的比较实际的功能。主要原因是学习:)

+5

你可以使用'np.vectorize'或'np。apply_along_axis'使任意的Python函数以“矢量化”的方式运行,但是这些通用的解决方案与标准的Python for循环相比,性能优势可以忽略不计。为了获得任何有意义的性能改进,您需要具体说明要矢量化的实际功能。 –

+0

谢谢ali_m。你的评论和奥利弗的回答提供了我正在寻找的答案。 – Sina

回答

2

ali_m的评论是专注:对于严重的速度增益,您应该更具体地了解该功能的功能。

话虽这么说,如果你仍然想使用np.apply_along_axis得到一个(可能)小速度提升,再考虑(重读that function's docstring后),您可以轻松地

  1. 换你的函数产生一维数组,
  2. 使用np.apply_along_axis与包装和
  3. 重塑所得阵列:

    def test_func_wrapper(*args, **kwargs): 
        return test_func(*args, **kwargs).ravel() 
    
    output = np.apply_along_axis(test_func_wrapper, m=3, n=2, axis=1, arr=input_data) 
    np.allclose(output.reshape(5,3, -1), output_data) 
    # output: True 
    

请注意,这是一个通用方式来加速这样的循环。如果您使用更具体到实际问题的功能,您可能会获得更好的性能。

相关问题