Efficiently generating a Cauchy matrix from two Numpy arrays

Question

A Cauchy matrix (Wikipedia article) is a matrix determined by two vectors (arrays of numbers). Given two vectors x and y, the Cauchy matrix C generated by them is defined entry-wise as

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

Given two Numpy arrays x and y, what is an efficient way to generate a Cauchy matrix?

Answer 1

This is the most efficient way I found, using array broadcasting to take advantage of vectorization.

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

---

Edit: @HYRY and @shx2 have suggested that, instead of x.reshape((-1,1)), which makes a copy, you can use x[:,np.newaxis], which returns a view of the same array. @HYRY also suggests 1.0/np.subtract.outer(x,y), which is slightly slower for me but maybe more explicit.

---

Example:

>>> 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]])

---

I tried a few other methods, all of which were significantly slower.

This is the naive approach, which costs list comprehension:

cauchym = numpy.array([[ 1.0/(x_i-y_j) for y_j in y] for x_i in x])

This one generates the matrix as a 1-dimensional array (saving the cost of nested Python lists) and reshapes it to a matrix afterward. It also moves the division to a single Numpy operation:

cauchym = 1.0/numpy.array([(x_i-y_j) for x_i in x for y_j in y]).reshape([len(x),len(y)])

Using numpy.repeat and numpy.tile (which respectively tile the array horizontally and vertically). This way makes unnecessary copies:

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]);