柯西矩阵 ( Wikipedia article ) 是由两个向量(数字数组)确定的矩阵。给定两个向量x
和y
,它们生成的柯西矩阵C
按条目定义为
C[i][j] := 1/(x[i] - y[j])
给定两个 Numpy 数组 x
和 y
,生成柯西矩阵的有效方法是什么?
最佳答案
这是我发现的最有效的方法,使用数组广播来利用矢量化。
1.0 / (x.reshape((-1,1)) - y)
<小时/>
编辑:@HYRY 和 @shx2 建议您可以使用 x[:, 而不是创建副本的
,它返回同一数组的 View 。 @HYRY 还建议 x.reshape((-1,1))
np.newaxis]1.0/np.subtract.outer(x,y)
,这对我来说稍微慢一些,但可能更明确。
示例:
>>> x = numpy.array([1,2,3,4]) #x
>>> y = numpy.array([5,6,7]) #y
>>>
>>> #transpose x, to nx1
... x = x.reshape((-1,1))
>>> x
array([[1],
[2],
[3],
[4]])
>>>
>>> #array of differences x[i] - y[j]
... #an nx1 array minus a 1xm array is an nxm array
... diff_matrix = x-y
>>> diff_matrix
array([[-4, -5, -6],
[-3, -4, -5],
[-2, -3, -4],
[-1, -2, -3]])
>>>
>>> #apply the multiplicative inverse to each entry
... cauchym = 1.0/diff_matrix
>>> cauchym
array([[-0.25 , -0.2 , -0.16666667],
[-0.33333333, -0.25 , -0.2 ],
[-0.5 , -0.33333333, -0.25 ],
[-1. , -0.5 , -0.33333333]])
<小时/>
我尝试了一些其他方法,所有这些方法都明显慢了。
这是一种幼稚的方法,需要列表理解:
cauchym = numpy.array([[ 1.0/(x_i-y_j) for y_j in y] for x_i in x])
这会将矩阵生成为一维数组(节省嵌套 Python 列表的成本),然后将其重新整形为矩阵。它还将除法移动到单个 Numpy 运算:
cauchym = 1.0/numpy.array([(x_i-y_j) for x_i in x for y_j in y]).reshape([len(x),len(y)])
使用numpy.repeat
和numpy.tile
(分别水平和垂直平铺数组)。这种方式会产生不必要的副本:
lenx = len(x)
leny = len(y)
xm = numpy.repeat(x,leny) #the i'th row is s_i
ym = numpy.tile(y,lenx)
cauchym = (1.0/(xm-ym)).reshape([lenx,leny]);
关于python - 从两个 Numpy 数组高效生成柯西矩阵,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/21427687/