Существует ли численно оптимальный порядок умножения матриц? - PullRequest
Новая
       45

Существует ли численно оптимальный порядок умножения матриц?

8 голосов
/ 08 июля 2019

TL; DR: вопрос о умножении ТОЧНОСТЬ

Я должен умножить матрицы A (100x8000), B (8000x27) и C (27x1).

Поскольку матрицы B и C являются постоянными, а A - переменными, я предпочитаю рассчитывать их как: ABC = np.dot(A, np.dot(B, C)). Однако мне интересно, что это может быть численно хуже (с точки зрения точность ), чем np.dot(np.dot(a, B), C).

Что может быть важно: матрицы A и B содержат 8000 выборок (соответственно) 100 и 27 коррелированных признаков.

Существует ли численно оптимальный (с точки зрения точность ) порядок умножения? Если да - как я могу это определить?

Особый случай

Можно предположить, что матрицы A и B неотрицательны. Кроме того: 

C = np.linalg.solve(cov(B, k), X)

, где X - матрица 27x1 из 27 (возможно, коррелированных) случайных величин неизвестного распределения, cov = lambda X, k: np.dot(X.T, X) + k * np.eye(X.shape[1]), а k - неотрицательная константа, минимизирующая выражение: 

sum((X[i, 0] - np.dot(np.dot(B[:, [i]].T, drop(B, i)),
                      np.linalg.solve(cov(drop(B, i), k),
                                      np.delete(X, i, axis=0))) **2
    for i in range(27))

Функция drop() определяется как lambda X, i: np.delete(X, i, axis=1).

Еще более особый случай

Можно предположить, что np.cov(B.T, B) является ковариационной матрицей X, которая следует многомерному распределению Гаусса.

Ответы [ 4 ]

4 голосов
/ 11 июля 2019

На данный момент лучшая идея (для определенного набора матриц), которую я имею, - провести следующий численный эксперимент:

  1. Рассчитать справочную матрицу как среднее для продуктов, рассчитанных с высокой точностью (например, `np.float128).
  2. Расчет тестовых продуктов с более низкой точностью (np.float64, np.float32, даже np.float16),
  3. Анализ ошибок, рассчитанных как разница между тестируемыми продуктами и эталонной матрицей. Ожидается, что ошибки уменьшатся, поскольку точность выше.
2 голосов
/ 16 июля 2019

Предложите использовать формулу, которая приводит к наименьшей ошибке накопления с плавающей запятой.

  • (AB)C: каждый элемент (AB) является точечным произведением каждой строки в A с каждымстолбец в B.Для каждого столбца в B найдите отношение элемента max к элементу min, а затем найдите максимальное из этих соотношений на столбец.Назовите это B_ratio.
  • A(BC): пусть T = (BC), в результате получается матрица 8000x1.Каждый элемент A(BC) является точечным произведением каждой строки в A с каждым [одним] столбцом в T.Для каждого [одного] столбца в T найдите отношение элемента max к элементу min.Назовите это T_ratio.

При условии равномерного распределения значений в A, в зависимости от того, какое из приведенных выше решений имеет меньшее абсолютное отношение , оно является более численно устойчивым.то есть, если fabs(B_ratio) < fabs(T_ratio), тогда ожидайте, что (AB)C будет лучше.

Обоснование: Ошибка накапливается при добавлении больших и малых чисел - младшие биты меньших чисел "теряются в шуме".При сохранении коэффициентов в пределах меньшего абсолютного спреда вероятность потери вкладов отдельных мелких членов уменьшается.

Подробности с плавающей запятой

При добавлении чисел с плавающей запятой IEEE 754 z = x + y один изМогут возникнуть следующие случаи:

  1. z.exponent == max(x.exponent, y.exponent) && z.mantissa != max(x,y).mantissa, что означает, что мантиссы суммируются без выполнения.Ошибка не внесена.
  2. z.exponent == max(x.exponent, y.exponent)+1, что означает мантиссы, суммированные с проведением.Наименьший бит точности теряется, в результате чего вводится ошибка 2^-(z.exponent - MANTISSA_BITS).
  3. z == max(x, y), означающая, что x или y были больше в 2 ^ MANTISSA_BITS, , что приводит к полной потередругой вход (он даже ниже порога шума).

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

Дополнительное чтение

Что должен знать каждый специалист по вычислительной технике об арифметике с плавающей точкой

0 голосов
/ 08 июля 2019

Хм, я не думаю, что есть разница в скорости для разных порядков умножения.Расчет количества должен быть таким же.Также я не вижу возможного улучшения кеширования, как в других матричных вычислениях (например, для порядка итераций в массиве).

Единственным моментом будет сохранение np.dot(B,C) в дополнительной матрице, если онидействительно не изменится, и если вам нужен результат в многократных вычислениях.

0 голосов
/ 08 июля 2019

Не будет ли умножение трех матриц всегда медленнее, чем умножение только двух?

У вас действительно есть только два варианта: (AB)C и A(BC). Поскольку BC = const, вы можете иметь постоянную T = BC формы 8000x1 и затем умножать AT без необходимости пересчитывать T.

...