从两个Numpy阵列有效生成柯西矩阵

Question

柯西矩阵（Wikipedia article）是由两个向量（数组数组）确定的矩阵。给定两个矢量x和y，由它们所产生的矩阵柯西C定义条目明智作为

C[i][j] := 1/(x[i] - y[j]) 

给定两个numpy的阵列x和y，什么是生成矩阵柯西一种有效的方法是什么？

Answer 1

这是我发现的最有效的方式，使用阵列广播来利用矢量化。

1.0/(x.reshape((-1,1)) - y) 

---

编辑：@HYRY和@ shx2建议，与其x.reshape((-1,1))，这使得副本，您可以使用x[:,np.newaxis]，返回同一阵列的视图。 @HYRY也暗示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]) 

这一个生成矩阵作为1维阵列（保存嵌套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]);