Оптимизация умножения матриц

Если все числа целые, можно попробовать реализовать умножение матриц на шаблонах.

Это делается так: "переменные" засовываются в enum'ы и non-type parameters, c помощью частичной специализации можно устроить некое подобие pattern-matching'а. Циклы делаются с помощью рекурсии.

Ну или кодогенерация.

Элементы A не целые. Кодогенерация естественно. Вопрос был не в этом.

нет

Это можно доказать?

можно преобразовать A в треугольник (или как там оно по научному) до компиляции.

Вопрос в том, можно ли что-то выкинуть, например? В общем случае, когда про вектор совсем ничего не известно — нет.

можно преобразовать A в треугольник

Э... и как это по научному, когда треугольная матрица ведет себя близко к полной?

а что тут доказывать? yi = sum aij xj, xj любые

тогда и икс надо на соответствующую матрицу домножить, свели задачу к исходной

хоть сформулируй четко чего ты хочешь. Что значит «ведет себя»?

Вопрос в том, можно ли что-то выкинуть, например?

Именно. Что типа разложения матрицы в ряд, которое бы выделило главные компоненты (ох, какую, право, антинаучную херню я несу).

В общем случае, когда про вектор совсем ничего не известно — нет.

Целые числа от 0 до 2^14.

которое бы выделило главные компоненты

Таких нѣтъ. Представь, что у тебя вектор из нулей, а одна компонента 100500. Какой бы матрица не была, эта компонента вектора всегда даст "главный" вклад.

Есть ли какой нибудь способ... небольшая (до 30х30) не разрежённая матрица

Нет, ибо: факторизация полезна только для разреженных матриц, а алгоритм Штрассена имеет преимущество при размерах > 64x64.

Если умножать только один раз, то дешевле просто взять и умножить. А вот если много — то действительно лучше найти собственный вектор и упростить матрицу.

А как это для умножения-то поможет, что-то до меня не доходит. Если матрицу привести в собственный базис, то и вектор придётся туда-обратно переводить.

тогда и икс надо на соответствующую матрицу домножить

нафига ? A и полученная из неё треугольная R, эквивалентны. RX считается мягко говоря быстрее чем AX.

А, хрен его знает: вдруг упростится матрица? А то и вообще к диагональному (ну или хотя бы треугольному) виду приведется?

// шучу. Тупость сморозил.

А то и вообще к диагональному виду приведется?

Если все собственные вектора линейно независимы, то в базисе из них матрица будет именно диагональной. Но вектор придётся приводить в этот базис и обратно двумя умножениями на матрицы. Не весело это.

Ты не понял: если умножений на одну и ту же матрицу 100500 штук, то если ты ее приведешь к диагональному виду, ускоришь минимум в одну-полторы степени размера матрицы!

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

Ну, если вообще все вычисления в этот базис удастся перевести, то да.

Кстати, не любую матрицу можно привести к диагональному виду. ЕМНИП, самый простой вид, к которому можно свести любую матрицу, - это нормальная жорданова форма

30*30*14 = 12килобайт. На нормальном проце в кеш врезет. Главное X транспонируй перед умножением, всё летать будет.

Профилирование делал? Проблемы именно в этом умножении?

Алгоритм Штрассена

Большая ли матрица X? Или может быть 30х10000?

В два «треугольника». LU-разложение называется.

A и полученная из неё треугольная R, эквивалентны. RX считается мягко говоря быстрее чем AX.

Ересь. Если Rx=Ax для любых x, то R=A.

30*30*14 = 12килобайт. На нормальном проце в кеш врезет. Главное X транспонируй перед умножением, всё летать будет.

+1. Я бы тоже сосредоточился на эффективном использовании кэшей и, если нужна экстремальная скорость, на SSE. А вот и годная статейка про кэши и умножение матриц: http://www.akkadia.org/drepper/cpumemory.pdf

Да, походу общего рецепта нет. А жаль. Для больших матриц много чего придумали, для маленьких придется делать по-честному.

30*30*14 = 12килобайт.

Недавно, помню, тут был разговор на тему, а зачем 64 бита, если на самом деле приложению доступно только 48. Рассуждать можно долго, но всё равно должно быть кратно 8-ми: 30*30*16=14КБ. Ещё в кэш должен влезть сам вектор.

В два «треугольника». LU-разложение называется.

количество операций останется прежним.

A и полученная из неё треугольная R, эквивалентны. RX считается мягко говоря быстрее чем AX.

Ересь. Если Rx=Ax для любых x, то R=A.

только если det(R)!=0

SSE, наверное. И всякие операции скалярных произведений на потоки разбивать. Только не слушай активистов со Штрассеном, на 30*30 разницы не заметишь. А числа - примитивные типы или какой-то библиотеки? Если библиотека, то в ней уже должны быть оптимизированные операции.

только если det(R)!=0

Пример в студию.

не будет примера. Я не заметил слов «для любых x» :)

Это все понятно. Я знаю про MKL, OpenBLAS (бывший GotoBLAS) и Atlas. Вопрос был в том можно ли так исхитриться что бы считать не все, а только часть для известной ошибки результата.

можно попробовать следующий прикол:

есть такая штука, как разложение Ивасавы. То есть мы разлагаем матрицу на произведение A=ODL, где O - это ортоногальная матрица, D - диагональная и L - треугольная.

тогда мы имеем O^{-1} Y = DLX. правую часть посчитать очевдно быстрее можно, ибо DL - треугольная матрица.

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