Я пытаюсь реализовать некоторые основные операции линейной алгебры, и одна из этих операций - инверсия треугольной (верхней и/или нижней) матрицы. Есть ли простой и стабильный алгоритм для этого?
Спасибо.
Да, используйте обратную замену. Стандартный алгоритм инвертирования матрицы состоит в том, чтобы найти ее LU-декомпозицию (разложение на нижнетреугольную и верхнетреугольную матрицу), использовать обратную подстановку на треугольных кусках, а затем объединить результаты, чтобы получить обратную исходную матрицу.
Не инвертируйте его, если сможете. Это одна из основных заповедей числовой линейной алгебры.
Гораздо быстрее и численно более стабильно хранить матрицу L в памяти и вычислять
inv(L)bОбратите внимание, что обычный алгоритм его инвертирования требует решения систем
inv(L)[1 0 0 ...],
inv(L)[0 1 0 ....],
inv(L)[0 0 1 ....]Учитывая нижнюю треугольную матрицу L, backsubstitution позволяет вам решить систему L x = b быстро для любой правой части b.
Чтобы инвертировать L, вы можете решить эту систему для правых частей e1 = (1,0,..., 0), e2 = (0,1,..., 0),..., en = (0,0,..., 1) и объединить полученные векторы решения в единую (обязательно нижнетреугольную) матрицу.
Если вы заинтересованы в решении с закрытой формой, диагональные элементы обратного являются инверсиями исходных диагональных элементов, а формула для остальных элементов обратного становится все сложнее, когда вы двигаетесь из диагонали.
Если вы говорите о действиях с одиночной точностью, посмотрите исходный код для подпрограмм LAPACK STRTRI и STRTI2.
Ничего себе, что практически половина содержания курса численного анализа. Стандартные алгоритмы сделают это, и есть куча консервированного кода здесь. Конечным источником этой и большинства других обычных задач численного анализа является Численные рецепты.