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

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

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 умнее о симметричных продуктах, чтобы сохранить флопы и обеспечить симметричный вывод...

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

Ответ 1

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

Ранее операции GAMM 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

Ответ 2

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

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 (ymmv)