Numpy точка слишком умная о симметричных умножениях

Question

Кто-нибудь знает о документации для этого поведения?

import numpy as np
A  = np.random.uniform(0,1,(10,5))
w  = np.ones(5)
Aw = A*w
Sym1 = Aw.dot(Aw.T)
Sym2 = (A*w).dot((A*w).T)
diff = Sym1 - Sym2

diff.max () близка к точности, отличной от нуля, например, 4.4e-16.

Это (расхождение с 0) обычно хорошо... в мире с конечной точностью мы не должны удивляться.

Более того, я предполагаю, что Numpy умно относится к симметричным продуктам, чтобы сохранить флопы и обеспечить симметричный выход...

Но я имею дело с хаотическими системами, и это небольшое несоответствие быстро становится заметным при отладке. Поэтому я хотел бы точно знать, что происходит.

Answer 1

Это поведение является результатом изменения, внесенного для NumPy 1.11.0, в запросе № 6932. Из заметок о выпуске для 1.11.0:

Ранее операции gemm BLAS использовались для всех продуктов матрицы. Теперь, если матричный продукт находится между матрицей и ее транспонированием, он будет использовать операции syrk BLAS для повышения производительности. Эта оптимизация была расширена до @, numpy.dot, numpy.inner и numpy.matmul.

В изменениях для этого PR можно найти этот комментарий:

/*
 * Use syrk if we have a case of a matrix times its transpose.
 * Otherwise, use gemm for all other cases.
 */

Поэтому NumPy делает явную проверку для случая, когда матрица умножает ее транспонирование, и в этом случае вызывает другую базовую функцию BLAS. Как отмечает @hpaulj в комментарии, такая проверка обходится NumPy дешево, поскольку транспонированный 2d-массив - это просто представление исходного массива с инвертированной формой и шагами, поэтому достаточно проверить несколько фрагментов метаданных в массивах (вместо того, чтобы сравнивать фактические данные массива).

Вот немного более простой случай, который показывает несоответствие. Обратите внимание, что с помощью .copy на одном из аргументов dot достаточно, чтобы победить специальный корпус NumPy.

import numpy as np
random = np.random.RandomState(12345)
A = random.uniform(size=(10, 5))
Sym1 = A.dot(A.T)
Sym2 = A.dot(A.T.copy())
print(abs(Sym1 - Sym2).max())

Я предполагаю, что одно из преимуществ этого специального корпуса, помимо очевидного потенциала для ускорения, заключается в том, что вы гарантированы (я надеюсь, но на практике это будет зависеть от реализации BLAS), чтобы получить идеально симметричный результат, когда syrk используется, а не матрица, которая просто симметрична с точностью до числовой ошибки. В качестве (по общему признанию, не очень хорошего) теста для этого я попытался:

import numpy as np
random = np.random.RandomState(12345)
A = random.uniform(size=(100, 50))
Sym1 = A.dot(A.T)
Sym2 = A.dot(A.T.copy())
print("Sym1 symmetric: ", (Sym1 == Sym1.T).all())
print("Sym2 symmetric: ", (Sym2 == Sym2.T).all())

Результаты на моей машине:

Sym1 symmetric:  True
Sym2 symmetric:  False

Answer 2

Я подозреваю, что это связано с продвижением промежуточных регистров с плавающей запятой с точностью до 80 бит. Отчасти подтверждение этой гипотезы заключается в том, что если мы используем меньше чисел с плавающей точкой, мы последовательно получаем 0 в наших результатах, ала

A  = np.random.uniform(0,1,(4,2))
w  = np.ones(2)
Aw = A*w
Sym1 = Aw.dot(Aw.T)
Sym2 = (A*w).dot((A*w).T)
diff = Sym1 - Sym2
# diff is all 0's (ymmv)