Решение системы линейных уравнений через LDLT-разложение
Мы продемонстрировали, как с помощью LDLT-разложения можно обратить матрицу. Это всё равно довольно затратный процесс, требующий порядка 2/3N3 умножений. В нашей задаче это было необходимо, поскольку обратная матрица - это ковариационная матрица шума измерений, её мы должны "скормить" фильтру Калмана вместе с вычисленными параметрами, что позволит ему правильно объединить эти данные с данными от других датчиков.
Но если всё, что нам нужно - это решить систему линейных уравнений
(где A - симметричная положительно определённая размером N×N, в задачах оптимизации, когда мы пытаемся свести к минимуму некоторую "невязку", они встречаются сплошь и рядом,
x - неизвестный вектор N×1,
b - известный вектор N×1)
то обращать матрицу не надо, есть решение гораздо более эффективное. Настолько эффективное, что после кое-какой доработки индийским напильником оно позволяет обращать матрицу даже быстрее, чем способ, рассмотренный ранее...
Начало работы точно такое же: мы проводим LDLT - разложение, как описано здесь:
подставляем это в уравнение:
DLTx - это некоторый вектор N×1, обозначим его c. Получаем:
То есть, получилась новая система линейных уравнений относительно c:
%5Cleft(%5Cbegin{array}{c}c_1%5C%5Cc_2%5C%5Cc_3%5C%5Cc_4%5Cend{array}%5Cright)=%5Cleft(%5Cbegin{array}{c}b_1%5C%5Cb_2%5C%5Cb_3%5C%5Cb_4%5Cend{array}%5Cright))
Решается она замечательно, "прямой подстановкой" сверху вниз:
Как видно, и здесь мы можем делать вычисления "на месте" - как только мы находим c1, нам больше не нужен b1, и так далее.
Поскольку нам нужно найти лишь N значений, решение получается очень "дешёвым": для нахождения cj мы совершаем j-1 умножений и столько же сложений, что даёт в общей сложности (N-1)N/2 умножений и столько же сложений.
Теперь выпишем, что же такое вектор c:
Умножаем обе части слева на D-1:
)
Как видно, полезно уже при осуществлении LDLT-разложения находить не D, а сразу D-1 - там по ходу дела эти обратные величины оказываются задействованы всё равно! В любом случае, тут добавляется ещё N операций - либо умножений, либо делений, в зависимости от реализации. Как всегда, операцию удаётся выполнить "на месте", приходя к новому вектору c'.
Остаётся решить уравнение
или
%5Cleft(%5Cbegin{array}{c}x_1%5C%5Cx_2%5C%5Cx_3%5C%5Cx_4%5Cend{array}%5Cright)=%5Cleft(%5Cbegin{array}{c}c%27_1%5C%5Cc%27_2%5C%5Cc%27_3%5C%5Cc%27_4%5Cend{array}%5Cright))
Такое уравнение ничуть не сложнее предыдущего, только в этот раз надо идти снизу вверх:
Те же самые (N-1)N/2 умножений и столько же сложений.
Как видно, вместо 2/3N3 операций при обращении матрицы, решение системы линейных уравнений требует лишь порядка N3/6 операций - в 4 раза меньше.
Один из способов обращения матрицы состоит в следующем: решим одно за другим такие системы уравнений:
=%5Cleft(%5Cbegin{array}{c}1%5C%5C0%5C%5C0%5C%5C0%5Cend{array}%5Cright),)
=%5Cleft(%5Cbegin{array}{c}0%5C%5C1%5C%5C0%5C%5C0%5Cend{array}%5Cright),)
и ещё 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 умножений. Это больше, чем в методе, описанном нами ранее.
Но здесь возможен кое-какой финт ушами, о котором расскажем в следующей части.
Но если всё, что нам нужно - это решить систему линейных уравнений
(где 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 умножений. Это больше, чем в методе, описанном нами ранее.
Но здесь возможен кое-какой финт ушами, о котором расскажем в следующей части.