трансляция просмотров нерегулярно NumPy

Question

Предполагая, что я хочу иметь массив с размером (n,m), где n очень большой, но с большим количеством дубликатов, т.е.0:n1 идентичны, n1:n2 идентичны и т. Д. (С n2%n1!=0, то есть не с регулярными интервалами).Есть ли способ сохранить только один набор значений для каждого из дубликатов при просмотре всего массива?

пример:

unique_values = np.array([[1,1,1], [2,2,2] ,[3,3,3]]) #these are the values i want to store in memory
index_mapping = np.array([0,0,1,1,1,2,2]) # a mapping between index of array above, with array below

unique_values_view = np.array([[1,1,1],[1,1,1],[2,2,2],[2,2,2],[2,2,2], [3,3,3],[3,3,3]]) #this is how I want the view to look like for broadcasting reasons

Я планирую умножить массив (просмотр)по некоторому другому массиву размером (1,m) и возьмем скалярное произведение этого произведения:

other_array1 = np.arange(unique_values.shape[1]).reshape(1,-1) # (1,m)
other_array2 = 2*np.ones((unique_values.shape[1],1)) # (m,1)
output = np.dot(unique_values_view * other_array1, other_array2).squeeze()

Выход - это одномерный массив длины n.

Answer 1

Ваше выражение допускает две существенные оптимизации:

• выполнять индексацию в последнюю очередь
• , умножить other_array1 сначала на other_array2, а затем на unique_values

Давайте применим следующие оптимизации:

>>> output_pp = (unique_values @ (other_array1.ravel() * other_array2.ravel()))[index_mapping]

# check for correctness
>>> (output == output_pp).all()
True

# and compare it to @Yakym Pirozhenko's approach
>>> from timeit import timeit
>>> print("yp:", timeit("np.dot(unique_values * other_array1, other_array2).squeeze()[index_mapping]", globals=globals()))
yp: 3.9105667411349714
>>> print("pp:", timeit("(unique_values @ (other_array1.ravel() * other_array2.ravel()))[index_mapping]", globals=globals()))
pp: 2.2684884609188884

Эти оптимизации легко обнаружить, если мы наблюдаем две вещи:

(1), если A является mxn -матрицей иb является n -вектором, тогда

A * b == A @ diag(b)
A.T * b[:, None] == diag(b) @ A.T

(2), если A является mxn -матрицей и I является k -вектором целых чисел из range(m) затем

A[I] == onehot(I) @ A

onehot можно определить как

def onehot(I, m, dtype=int):
    out = np.zeros((I.size, m), dtype=dtype)
    out[np.arange(I.size), I] = 1
    return out

Используя эти факты и сокращая uv, im, oa1 и oa2, мы можем написать

uv[im] * oa1 @ oa2 == onehot(im) @ uv @ diag(oa1) @ oa2

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

onehot(im) @ (uv @ (diag(oa1) @ oa2))

Используя (1) и (2) в обратном направлении, мы получаемоптимизированное выражение с начала этого поста.

Answer 2

Исходя из вашего примера, вы можете просто вычислить отображение индекса до конца:

output2 = np.dot(unique_values * other_array1, other_array2).squeeze()[index_mapping]

assert (output == output2).all()