Инвертировать матрицу 4x4 - необходимо самое устойчивое численное решение

Question

Я хочу инвертировать матрицу 4x4. Мои номера хранятся в формате с фиксированной запятой (точнее 1.15.16).

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

Примечание. В арифметике с фиксированной запятой я всегда отбрасываю некоторые младшие биты непосредственных результатов.

Итак - Какой самый числовой стабильный способ инвертировать матрицу? Я не особо беспокоюсь о производительности, но просто перейти к плавающей точке будет означать замедление моей целевой архитектуры.

Answer 1

Мета-ответ: действительно ли это общая матрица 4х4? Если ваша матрица имеет специальную форму, то есть прямые формулы для инвертирования, которые будут быстрыми и уменьшат количество операций.

Например, если это стандартное однородное преобразование координат из графики, например:

[ux vx wx tx]
[uy vy wy ty]
[uz vz wz tz]
[ 0  0  0  1]

(при условии композиции вращения, масштаба, матрицы перевода)

тогда есть легко выводимая прямая формула , которая равна

[ux uy uz -dot(u,t)]
[vx vy vz -dot(v,t)]
[wx wy wz -dot(w,t)]
[ 0  0  0     1    ]

(матрицы ASCII, украденные со связанной страницы.)

Вы, вероятно, не можете победить это из-за потери точности в фиксированной точке.

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

Answer 2

Я хотел бы ответить на второй вопрос, поднятый Джейсоном С. Вы уверены, что вам нужно инвертировать матрицу? Это почти никогда не нужно. Мало того, это часто плохая идея. Если вам нужно решить Ax = b, численно более стабильно решать систему напрямую, чем умножать b на A в обратном порядке.

Даже если вам придется снова и снова решать Ax = b для многих значений b, инвертировать A. все равно не очень хорошая идея. сохраните факторы, чтобы не повторять эту работу каждый раз, но вы все равно решали бы систему каждый раз, используя факторизацию.

Answer 3

Я думаю, что ответ на этот вопрос зависит от точной формы матрицы. Стандартный метод декомпозиции (LU, QR, Cholesky и т. Д.) С поворотом (необходим) достаточно хорош для фиксированной точки, особенно для небольшой матрицы 4x4. См. Книгу «Численные рецепты» Press et al. для описания этих методов.

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

Answer 4

Позвольте мне задать другой вопрос: вам определенно нужно инвертировать матрицу (назовите ее M) или вам нужно использовать матрицу, обратную для решения других уравнений? (например, Mx = b для известных M, b) Часто существуют другие способы сделать это без явного вычисления обратного. Или, если матрица M является функцией времени и она изменяется медленно, вы можете рассчитать полное обратное значение один раз, и есть итерационные способы его обновления.

Answer 5

Вы можете рассмотреть возможность удвоения до 1,31, прежде чем выполнять свой обычный алгоритм. Это удвоит количество умножений, но вы делаете инвертирование матрицы, и все, что вы делаете, будет довольно сильно привязано к множителю в вашем процессоре.

Для тех, кто заинтересован в нахождении уравнений для инвертирования 4x4, вы можете использовать символический математический пакет, чтобы решить их для себя. TI-89 сделает это даже, хотя это займет несколько минут.

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

-Adam

Answer 6

Чтобы свести к минимуму ошибки усечения и другие ошибки, используйте «поворот» - см. Главу по инвертированию матриц в числовых рецептах. У них есть лучшее объяснение, которое я нашел до сих пор.

Answer 7

Если матрица представляет аффинное преобразование (часто это имеет место с матрицами 4x4, если вы не вводите компонент масштабирования), инверсия - это просто транспонирование верхней части вращения 3x3 с отрицанием последнего столбца. Очевидно, что если вам требуется обобщенное решение, то поиск исключения Гаусса, вероятно, самый простой.

Answer 8

Простое старое гауссовское исключение будет хорошо работать.

Это зависит от того, какие библиотеки / классы / структуры вы используете. Вы можете взглянуть на GSL .