Numpy: создать матрицу из декартового произведения двух векторов (или одного с самим собой), применяя функцию ко всем парам

Question

Чтобы пояснить немного, я хочу создать ковариационную матрицу, в которой каждый элемент определяется функцией ядра k(x, y), и я хочу сделать это для одного вектора.Это должно быть что-то вроде:

# This is given
x = [x1, x2, x3, x4, ...]

# This is what I want to compute
result = [[k(x1, x1), k(x1, x2), k(x1, x3), ...],
          [k(x2, x1), k(x2, x2), ...],
          [k(x3, x1), k(x3, x2), ...],
          ...]

, но, конечно, это должно быть сделано в массивах, в идеале без использования Python, из-за производительности.Если бы я не заботился о производительности, я бы просто написал:

result = np.zeros((len(x), len(x)))

for i in range(len(x)):
  for j in range(len(x)):
     result[i, j] = k(x[i], x[j])

Но я чувствую, что должен быть более идиоматичный способ написания этого паттерна.

Answer 1

Другое решение состоит в том, чтобы позволить numpy делать само вещание:

import numpy as np

def k(x,y):
    return x**2+y

def meshgrid1D_view(x):
    shp = (len(x),len(x))
    mesh1 = np.broadcast_to(x,shp)
    mesh2 = np.broadcast_to(x[:,None],shp)    
    return mesh1, mesh2

x = np.random.randint(0,10,(10000))
b=k(a[:,None],a[None,:])

def sol0(x):
    k(x[:,None],x[None,:])

def sol1(x):
    x,y=np.meshgrid(x,x)
    k(x,y)

def sol2(x):
    x,y=meshgrid1D_view(x)
    k(x,y)

%timeit sol0(x)
165 ms ± 1.67 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit sol1(x)
655 ms ± 6.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit sol2(x)
341 ms ± 2.91 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Вы видите, что это более эффективно и меньше кода.

Answer 2

Если k работает с 2D-массивами, вы можете использовать np.meshgrid. Но это может привести к дополнительным затратам памяти. Одним из альтернативных вариантов может быть создание 2D видов сетки, аналогично np.meshgrid, например,

def meshgrid1D_view(x):
    shp = (len(x),len(x))
    mesh1 = np.broadcast_to(x,shp)
    mesh2 = np.broadcast_to(x[:,None],shp)    
    return mesh1, mesh2

Пробный прогон -

In [140]: x
Out[140]: array([3, 5, 6, 8])

In [141]: np.meshgrid(x,x)
Out[141]: 
[array([[3, 5, 6, 8],
        [3, 5, 6, 8],
        [3, 5, 6, 8],
        [3, 5, 6, 8]]), array([[3, 3, 3, 3],
        [5, 5, 5, 5],
        [6, 6, 6, 6],
        [8, 8, 8, 8]])]
In [142]: meshgrid1D(x)
Out[142]: 
(array([[3, 5, 6, 8],
        [3, 5, 6, 8],
        [3, 5, 6, 8],
        [3, 5, 6, 8]]), array([[3, 3, 3, 3],
        [5, 5, 5, 5],
        [6, 6, 6, 6],
        [8, 8, 8, 8]]))

Как это помогает?

Это помогает с эффективностью памяти и, следовательно, производительности. Давайте проверим на больших массивах, чтобы увидеть разницу -

In [143]: x = np.random.randint(0,10,(10000))

In [144]: %timeit np.meshgrid(x,x)
10 loops, best of 3: 171 ms per loop

In [145]: %timeit meshgrid1D(x)
100000 loops, best of 3: 6.91 µs per loop