nabbla (nabbla1) wrote,
nabbla
nabbla1

Categories:

Решение системы линейных уравнений через LDLT-разложение

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

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


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

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



Начало работы точно такое же: мы проводим LDLT - разложение, как описано здесь:



подставляем это в уравнение:



DLTx - это некоторый вектор N×1, обозначим его c. Получаем:


То есть, получилась новая система линейных уравнений относительно c:



Решается она замечательно, "прямой подстановкой" сверху вниз:






Как видно, и здесь мы можем делать вычисления "на месте" - как только мы находим c1, нам больше не нужен b1, и так далее.

Поскольку нам нужно найти лишь N значений, решение получается очень "дешёвым": для нахождения cj мы совершаем j-1 умножений и столько же сложений, что даёт в общей сложности (N-1)N/2 умножений и столько же сложений.

Теперь выпишем, что же такое вектор c:



Умножаем обе части слева на D-1:


Как видно, полезно уже при осуществлении LDLT-разложения находить не D, а сразу D-1 - там по ходу дела эти обратные величины оказываются задействованы всё равно! В любом случае, тут добавляется ещё N операций - либо умножений, либо делений, в зависимости от реализации. Как всегда, операцию удаётся выполнить "на месте", приходя к новому вектору c'.

Остаётся решить уравнение


или


Такое уравнение ничуть не сложнее предыдущего, только в этот раз надо идти снизу вверх:





Те же самые (N-1)N/2 умножений и столько же сложений.

Как видно, вместо 2/3N3 операций при обращении матрицы, решение системы линейных уравнений требует лишь порядка N3/6 операций - в 4 раза меньше.

Один из способов обращения матрицы состоит в следующем: решим одно за другим такие системы уравнений:





и ещё 2 таких же. Иными словами, найти обратную матрицу столбец за столбцом. Если сделать это "в лоб", то нам потребуется по N2/2 умножений на каждый столбец при первой подстановке и ещё столько же при второй, то есть суммарно: N3 операций.

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

Выходит, что первый столбец мы считаем "честно", затратив (N-1)N/2 умножений, второй - не досчитав 1 значения, затратив (N-2)(N-1)/2 умножений, ну а последний вообще работы не требует - там уже лежит правильное значение. В итоге количество умножений на второй подстановке (из c' в x) падает до (N-1)N(N+1)/6, то есть с порядка N3/2 до N3/6, фактически втрое!

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

Поэтому имеем общую сложность N3/6 (в LDLT разложении) + N3/2 (первая подстановка) + N3/6 (вторая подстановка) = 5/6N3 умножений. Это больше, чем в методе, описанном нами ранее.

Но здесь возможен кое-какой финт ушами, о котором расскажем в следующей части.
Tags: математика, программки, работа
Subscribe

  • Нахождение двух самых отдалённых точек

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

  • Слишком общительный счётчик

    Вчера я чуть поторопился отсинтезировать проект,параметры не поменял: RomWidth = 8 вместо 7, RamWidth = 9 вместо 8, и ещё EnableByteAccess=1, чтобы…

  • Балансируем конвейер QuatCore

    В пятницу у нас всё замечательно сработало на симуляции, первые 16 миллисекунд полёт нормальный. А вот прошить весь проект на ПЛИС и попробовать "в…

  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your IP address will be recorded 

  • 0 comments