Оптимизация матричного произведения в numpy

Что такое симметричная неквадратная матрица?

Прошу прощения, описался. Имелось в виду, учесть результат вычисления A.T.dot(A) есть симметричная матрица.

Ничего не знаю про numpy, но точно нужна именно сама матрица A^T * A ? Может в самой задаче нужно не именно это, а результат умножения этой матрицы на что-то? Тогда можно пробовать как-нибудь по-умному вычислять результат.

А вообще, если сходу, можно посчитать только те элементы (x_ij) искомой матрицы, у которых i>=j. А потом x_ji = x_ij. Но это уточнение от Кэпа, разумеется.

Ну и получишь ты сложность ~N^3, в то время как умножение двух матриц имеет сложность ~N^(2+x), где x<1. ЕМНИП есть алгоитмы где x~=0.5.

Оптимизировать имеет смысл там, где оптимизация нужна. А с такими задачами ТС вряд ли столкнулся, во всех остальных случаях надо отдавать предпочтения ясности кода, даже в ущерб его «оптимальности». Кроме того, вызывает некий скептицизм сама задача — оптимизировать скриптовый код для перемножения матрицы, в сравнении с функциями написанными на си, которые уже есть в штатных библиотеках питона.

Да, мне действительно не нужна эта матрица «сама в себе». Она используется в процессе Levenberg–Marquardt algorithm применительно к нейронным сетям.

Итоговое выражение:

inv(J.T.dot(J) + mu * I).dot(J.T.dot(e))

где J - матрица (numpy 2-d array) размером P на N, mu - скаляр, I - единичная матрица N на N, e - вектор длинной P (numpy 1-d array).

Таким образом, резульатом выражения является вектор N на 1.

Если кто-то поможет записать это быстрее, буду благодарен.

Размеры матриц какие? Примерно.

Средствами numpy никак не оптимизируешь. Самая долгая операция тут это операция обращения матрицы. ~N^3 она асимптотически давит перемножение. А ее тут никак не оптимизировать.

Гм.. а почему же тогда исключение J.T.dot(J) из выражения даёт такой большой прирост в скорости?

import numpy as np

N = 768
P = 1024

J = np.random.random((P, N))
e = np.random.random(P)

timeit('np.linalg.inv(J.T.dot(J)).dot(J.T.dot(e))')
# 5 loops, best of 3: 4.18 s per loop

H = J.T.dot(J)
timeit('np.linalg.inv(H).dot(J.T.dot(e))')
# 5 loops, best of 3: 745 ms per loop

Потому что «в асимптотике». А у тебя размеры маленькие.

В общем, размеры примерно такие и есть.

По-прежнему ищуются мудрецы, могущие подсказать, можно ли реализовать это быстрее малой кровью?

Это сложно. Даже переход на cython на большинстве операций ничего не даст - numpy и так уже оптимизирован по самые помидоры. Возможно здесь имеются живые пользователи фортрана, которые поделятся своим опытом - по идее он быстрее.

учтя каким-либо образом информацию о том, что результат - симметричная матрица

Ну в два раза и можно (если сделать оптимизацию настолько же хорошо, как авторы numpy). Легче станет от этого?

Бывший живой пользователь фортрана, в нем кроме coarrays и синтаксиса ничего экстраординарного нет.

Что такое ^T?

Просто транспонирование.

А, ты в ОП не поправил, что А — несимметричная.

Тем, кто нагуглит этот тред: выражением inv(J.T.dot(J) + mu * I).dot(J.T.dot(e)) можно оптимизировать при помощи numpy.linalg.solve.