Как удобно использовать операции с массивами numpy fortran contiguos?

Question

Некоторые numpy функции, такие как np.matmul(a, b), имеют удобное поведение для стеков матриц.

В руководстве указано:

Если один из аргументов равен ND, N> 2, он обрабатывается как стек матриц, находящихся в последних двух индексах и передаваемых соответственно.

Таким образом, для a.shape = (10 , 2, 4) и b.shape(10, 4, 2) утверждение a @ b имеет смысл и будет иметь форму (10, 2, 2)

Тем не менее, я пришел из мира линейной алгебры, где я привык к непрерывному расположению массива на Фортране.

Тот же a, представленный в виде непрерывного массива Фортрана, будет иметь форму (4, 2, 10) и аналогично b.shape = (2, 4, 10).

Для выполнения a @ b, как и раньше, мне придется вызывать (a.T @ b.T).T .

Еще хуже, предположим, что вы наивно создали один и тот же непрерывный массив Фортрана a с учетом поведения matmul, так что он имеет форму (10, 4, 2). Затем a.strides = (8, 80, 320) с наименьшим шагом в индексе «стека», который на самом деле должен иметь самый высокий шаг.

Это действительно путь к go или я что-то упустил?

Answer 1

Хотя numpy может обрабатывать все виды макетов, многие детали разработаны с учетом макета "C". Хорошие примеры - это то, как вложенные списки преобразуются в массивы, и как numpy операции пакетируют избыточные измерения, как в случае с матмулами.

Правильно, что в результате numpy, как правило, не зависит от макет массива (FORTRAN, C, несмежный); Скорость, тем не менее, определенно имеет место и в значительной степени так:

rng = np.random.default_rng()
a = rng.random((100,111,200))
b = rng.random((111,77,200))
af = np.array(a,order="F")
bf = np.array(b,order="F")

np.allclose((b.T@a.T).T,(bf.T@af.T).T)
# True
timeit(lambda:(b.T@a.T).T,number=10)
# 5.972857117187232
timeit(lambda:(bf.T@af.T).T,number=10)
# 0.1994628761895001

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

timeit(lambda:(np.array(b.T,order="C")@np.array(a.T,order="C")).T,number=10)
# 0.3931349152699113

Мой совет: если вам нужна скорость и удобство, вероятно, лучше всего go с макетом "C", чтобы привыкнуть к нему, уйдет не так много времени и избавит вас от многих потенциальных головных болей.

Answer 2

Матричное умножение

numpy работает независимо от внутренней структуры массива. Например, вот два C -упорядоченных массива:

>>> import numpy as np
>>> a = np.random.rand(10, 2, 4)
>>> b = np.random.rand(10, 4, 2)
>>> print('a', a.shape, a.strides)
>>> print('b', b.shape, b.strides)
a (10, 2, 4) (64, 32, 8)
b (10, 4, 2) (64, 16, 8)

Вот эквивалентные массивы в порядке Фортрана:

>>> af = np.asfortranarray(a)
>>> bf = np.asfortranarray(b)
>>> print('af', af.shape, af.strides)
>>> print('bf', bf.shape, bf.strides)
af (10, 2, 4) (8, 80, 160)
bf (10, 4, 2) (8, 80, 320)

Numpy обрабатывает эквивалентные массивы как эквивалентные, независимо от того, их внутренней компоновки:

>>> np.allclose(a, af) and np.allclose(b, bf)
True

Результаты умножения матриц не зависят от внутренней компоновки:

>>> np.allclose(a @ b, af @ bf)
True

, и вы даже можете смешивать макеты, если вы будете sh:

>>> np.allclose(a @ bf, af @ b)
True

Короче говоря, наиболее удобный способ использования массивов, упорядоченных по Фортрану, в numpy - не беспокоиться о внутренней компоновке массива: форма - это все, что имеет значение.

Если ваш формы массива отличаются от ожидаемых numpy matmul API , лучше всего изменить форму массивов, например, используя a.transpose(2, 0, 1) @ b.transpose(2, 0, 1) или аналогичный, в зависимости от того, что подходит для вашего варианта использования, но не волнуйтесь: для C или смежных массивов Fortran эта операция только корректирует метаданные вокруг представления массива, но не приводит к копированию или переупорядочению основного буфера данных.