从两个 Numpy 数组有效地生成柯西矩阵答案

【问题标题】：Efficiently generating a Cauchy matrix from two Numpy arrays从两个 Numpy 数组有效地生成柯西矩阵
【发布时间】：2014-02-21 01:20:12
【问题描述】：

柯西矩阵 (Wikipedia article) 是由两个向量（数字数组）确定的矩阵。给定两个向量x和y，由它们生成的柯西矩阵C逐项定义为

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

给定两个 Numpy 数组 x 和 y，生成柯西矩阵的有效方法是什么？

【问题讨论】：

标签： python optimization numpy matrix

【解决方案1】：

这是我发现的最有效的方法，使用数组广播来利用矢量化。

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

编辑：@HYRY 和@shx2 建议您可以使用x[:,np.newaxis]，而不是创建副本的x.reshape((-1,1))，它返回相同数组的视图。 @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])

这个将矩阵生成为一维数组（节省嵌套 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]);

【讨论】：

你也可以使用：1.0/np.subtract.outer(x, y),1.0/(x[:,None]-y)
或更明确的1.0/(x[:,np.newaxis]-y)

【解决方案2】：

我创建了一个函数，希望它能帮助你更好地理解。

# Creating a function in order to form a cauchy matrix
def cauchy_matrix(arr1,arr2):
    """
    Enter two arrays in order to get a cauchy matrix.The input array should be a 1-D array.
    arr1 = First 1-D array
    arr2 = Second 1-D array
    It returns the cauchy matrix having shape equal to m*n, where m is size of arr1 and n is size of arr2.
    """
    my_list = []
    try:
        for i in range(len(arr1)):
            for j in range(len(arr2)):
                z = 1/(arr1[i]-arr2[j])
                my_list.append(z)

        return np.array(my_list).reshape(arr1.shape[0],arr2.shape[0])
    except ZeroDivisionError:
        print("Check if both the arrays has '0' as one of it's element. One array can have a zero but both the arrays having '0' is not acceptable!")

【讨论】：