?

Log in

No account? Create an account

О реализации "индийского" обращения матрицы
QuatCore
nabbla1
Когда мы рассматривали "классический" способ обращения симметричной положительно определённой матрицы через LDLT-разложение, обращение нижней унитреугольной матрицы и умножение трёх матриц, мы особенно подчёркивали, в каком порядке надо проводить вычисления, чтобы все они были выполнены "на месте", не требуя дополнительной памяти вообще!

Но особую индийскую магию мы пока что оставили без деталей реализации, продемонстрировав лишь сам принцип.

Как оказывается, в этом методе всё-таки нужна вспомогательная память, длиной в одну строку (N элементов), может хватить и N-1 элемента.

И ещё не исключено, что индийцы пошли тропой Рамануджана, заново переоткрыв в 2011 году метод, который давным-давно был реализован на Фортране, просто никто туда "не лез"... Впрочем, доказывая мой исходный тезис про Write-only :)

Read more...Collapse )

И ещё замечание. Мы всё время говорили, что исходная матрица должна быть положительно определённой. Это не совсем так: для методов, основанных на LDLT-разложении матрица может быть любой (просто симметричной), но только в случае положительной определённости метод обладает хорошей устойчивостью. В противном случае рекомендуется делать pivoting, то бишь менять местами строки и столбцы, так, чтобы поместить на диагональ большие элементы.

LDLT-разложение: особая индийская магия!
ook
nabbla1
Мы изложили "стандартный" подход к обращению симметричной положительно определённой матрицы через LDLT-разложение с последующим обращением матриц L и D, и получением итогового результата как результат умножения 3 матриц.

В итоге количество умножений получилось порядка 2/3N3, что конечно получше, чем N3 в алгоритме обращения матрицы "общего вида", но не настолько лучше, как хотелось бы. Казалось, что когда у тебя вдвое меньше "уникальных" элементов как на входе, так и на выходе, то и количество операций должно упасть почти вдвое. Но всё-таки "окольные пути" взяли своё.

Но не так давно (в 2011 году с правками в 2013) два индийских математика опубликовали статью, в которой заявили, что их метод позволяет обратить такую матрицу всего за 1/2N3 умножений, что на 1/6 лучше существующих методов. Что ещё интереснее, они проверяли работу "с фиксированной точкой" и получили, что такой метод даёт чуть меньше ошибок вычислений (так я на эту статью и вышел - искал подсказки, как это всё получше выполнить без FPU). Я, конечно, не удивлюсь, если всё это уже было у Гаусса, но пока не слышал о таком :)

Сейчас попытаюсь объяснить, что же там происходит...
Read more...Collapse )

Решение системы линейных уравнений через LDLT-разложение
QuatCore
nabbla1
Мы продемонстрировали, как с помощью LDLT-разложения можно обратить матрицу. Это всё равно довольно затратный процесс, требующий порядка 2/3N3 умножений. В нашей задаче это было необходимо, поскольку обратная матрица - это ковариационная матрица шума измерений, её мы должны "скормить" фильтру Калмана вместе с вычисленными параметрами, что позволит ему правильно объединить эти данные с данными от других датчиков.

Но если всё, что нам нужно - это решить систему линейных уравнений


(где A - симметричная положительно определённая размером N×N, в задачах оптимизации, когда мы пытаемся свести к минимуму некоторую "невязку", они встречаются сплошь и рядом,
x - неизвестный вектор N×1,
b - известный вектор N×1)

то обращать матрицу не надо, есть решение гораздо более эффективное. Настолько эффективное, что после кое-какой доработки индийским напильником оно позволяет обращать матрицу даже быстрее, чем способ, рассмотренный ранее...

Read more...Collapse )

Обращение матрицы через LDLT-разложение: итоги
QuatCore
nabbla1
Мы описали, как это делается:
раз
два
три.

Теперь подведём итоги по количеству необходимых операций. Поначалу кажется, что это удаление гланд армейским способом - мы выполняем операцию ну очень окольными путями. Что ж, сейчас сравним.

Read more...Collapse )

Обращение симм. пол. опр. матрицы - последний этап
QuatCore
nabbla1
У нас была матрица A - симметричная N×N, положительно определённая. В памяти хранится только N(N+1)/2 её коэффициентов.

Чтобы обратить её, мы сначала сделали LDLT-разложение "на месте" (хотя можно и на новое место, если исходная матрица ещё пригодится).

Затем обратили матрицы L и D, также "на месте".

Кстати, матрицы L-1 с D-1 и сами по себе представляют интерес: если A-1, которую мы ищем - это ковариационная матрица, то L-1 совместно с D-1 позволяет генерировать вектор случайных величин, удовлетворяющий этой ковариационной матрице. Хотя LLT-разложение в этой задаче, наверное, чуть лучше подходит...

Наконец, нужно перемножить матрицы, чтобы найти обратную матрицу:


Посмотрим, как можно произвести это умножение "за один этап", и притом по-прежнему "на месте".

Read more...Collapse )

В следующем посте немного подведём итоги и сравним с методом "в лоб".

Обращение нижней унитреугольной матрицы
QuatCore
nabbla1
Нам хочется обратить симметричную положительно определённую матрицу А.

Первым делом мы сделали LDLT-разложение:


Здесь D-диагональная матрица, а L-нижняя треугольная матрица с единицами на главной диагонали, что иногда называют термином унитреугольная матрица.

Выразим обратную матрицу:


Обратить диагональную матрицу D легко, мы просто заменяем каждое число обратной величиной:


С нижней унитреугольной матрицей L самую чуточку посложнее.

Read more...Collapse )

Остаётся последний этап: помножить обращённую матрицу L-1 на D-1 и на себя же транспонированную, и тоже очень желательно "на месте", раз уж начали!

LDLT-разложение
QuatCore
nabbla1
В некотором смысле, "заметки для себя": алгоритмы линейной алгебры обладают нехорошим свойством write-only, особенно если там участвует "треугольная адресация". Пока разбираешь - всё вроде понятно, а чуть-чуть погодя смотришь на это безобразие - и оно превращается в набор вложенных циклов и десяток индексов массива, которые обновляются на каждой итерации не вполне понятным образом. Поэтому, пока не забыл - запишу :)

У нас есть симметричная положительно определённая матрица А:



мы часто будем пропускать верхний треугольник, поскольку он такой же. По сути, мы и в памяти предпочтём не дублировать эти компоненты (а тем более "в лоб" считать их по 2 раза).

Чтобы её обратить, можно воспользоваться довольно забавным трёхступенчатым методом, который начинается с LDLT-разложения:



Как видно, матрица L (от слова Lower) -нижняя треугольная матрица с единицами на гл. диагонали. D-диагональная матрица, причём в случае положительно определённой матрицы A, все её компоненты должны получиться ненулевыми. Данный метод - разновидность разложения Холецкого, но только в "обычном" разложении Холецкого мы записываем A = LLT, и приходится посчитать аж N (для матрицы N×N) квадратных корней. А если корни считать не хочется, то LDLT-разложение - то, что доктор прописал.

Две эти матрицы вполне помещаются в памяти, выделенной под матрицу A:



и это разложение может быть выполнено "на месте", не требуя временного хранилища.

Read more...Collapse )

Дальше расскажем, как обратить нижнюю треугольную матрицу с единицами на главной диагонали "на месте".

Всё это уже было в Симпсонах?
slow loris
nabbla1
Мне что-то в последнее время стали на YouTube предлагать доклады с CPPcon. Когда-то я в институте изучал C++, чуть-чуть немножко довелось на нём поковыряться, а по работе имею дело с "чистым" C (для всяких вычислителей), и немножко с программами рабочего места, которые тут были написаны в Borland C++ Builder.

И что-то я немножко в офиге от того, что вижу: вещи, которые Борланд решил ещё тогда, 20 лет назад, только сейчас очень осторожно и с кучей матов пытаются внедрить в стандарт...

Read more...Collapse )

Немножко ржу
beaver with chainsaw
nabbla1
800 000 человек в Калифорнии остались без света, потому как энергетическая компания PG&E испугалась сильных ветров, которые могут немножко порвать проводки и вызвать лесные пожары (в прошлом году они проиграли кучу исков и должны были из своего кармана исправлять последствия пожаров).

Read more...Collapse )

Хороший алгоритм, да рангом не вышел!
QuatCore
nabbla1
С нахрапу реализовать алгоритм сопровождения на ассемблере (худо-бедно описание здесь), в целых числах не получилось. Решил пока отставить в сторону масштабирование и сразу посмотреть, а что там с итоговой матрицей 6х6 и вектором 6х1?

Пока я всё это дело ковырял в плавающей точке, матрица получалась вот такой (с 300 метров, прямой наводкой):


(матрица симметричная, верхний треугольник не рисуем, чтобы не захламлять рисунок)

Числа прям подходящие для целочисленного исполнения - диагональные не превышают 65536 (а мы знаем, что на диагонали у нас неотрицательные значения), все остальные не превышают 32768 - прямо то что надо.

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

Read more...Collapse )