Почему один код (matmul) быстрее другого (Python)

Question

Проблема касается работы матрицы c. В следующем коде c1 по существу равен c2. Однако первый способ вычислений намного быстрее, чем второй. На самом деле, сначала я подумал, что для первого способа необходимо выделить матрицу ab, которая в два раза больше матрицы a, следовательно, может быть медленнее. Оказывается, все наоборот. Почему?

import time
import numpy as np
a = np.random.rand(20000,100)+np.random.rand(20000,100)*1j

tic = time.time()
b = np.vstack((a.real,a.imag))
c1 = b.T @ b
t1 = time.time()-tic

tic = time.time()
c2 = a.real.T @ a.real+a.imag.T@a.imag
t2 = time.time()-tic

print('t1=%f. t2=%f.'%(t1,t2))

Пример результата:

t1=0.037965. t2=4.375873.

Answer 1

a.real и a.imag находятся на месте, а np.vstack создает новую копию. Оператор @ (matmul()) работает с a.real, а a.imag занимает больше времени. Чтобы сделать это быстрее, вы можете создать копию каждого и затем передать ее @ или использовать np.dot(a.real.T, a.real) и np.dot(a.imag.T, a.imag) (я не уверен в каждой реализации в BLAS).

Для больших матриц первый метод в следующем коде все еще должен быть немного быстрее:

a = np.random.rand(20000,100)+np.random.rand(20000,100)*1j

tic = time.time()
b = np.vstack((a.real,a.imag))
c1 = b.T @ b
t1 = time.time()-tic

tic = time.time()
b = a.real.copy()
c = a.imag.copy()
c2 = b.T @ b + c.T @ c
t2 = time.time()-tic

print('t1=%f. t2=%f.'%(t1,t2))  

Вывод:

t1=0.031620. t2=0.021769.

РЕДАКТИРОВАТЬ: Погружение глубже:

Давайте кратко рассмотрим различные макеты памяти:

a = np.random.rand(20000,100)+np.random.rand(20000,100)*1j
print('a.flags\n', a.flags)
print('a.real flags\n', a.real.flags)


a.flags
C_CONTIGUOUS : True
F_CONTIGUOUS : False
OWNDATA : True
...
a.real flags
C_CONTIGUOUS : False
F_CONTIGUOUS : False
OWNDATA : False

Итак, a равно C_CONTIGUOUS, а a.real нет. Я не уверен, как @ реализует вычисления, но я думаю, что это отличается от трюков и шагов кэша и развернутых циклов. Я оставлю это экспертам для объяснения. Теперь array.copy() по умолчанию равно C_CONTIGUOUS ( БУДЬТЕ ОСТОРОЖНЫ: np.copy() не равно C_CONTIGUOUS по умолчанию. ), и поэтому второй подход выше так же быстр, как первый (где b также C_CONTIGUOUS).

@ George C: Я думаю, что первый немного быстрее, потому что np.vstack создает новый объект C_CONTIGUOUS, который может использовать уловки кеша в одном месте, в то время как во втором подходе вывод a.real.T @ a.real и a.imag.T@a.imag находятся в разных местах памяти и требуют дополнительных усилий для расчета. Вот ссылка на дополнительное объяснение.
Отказ от ответственности: Любой эксперт, который знает детали реализации NumPy, может редактировать этот пост.