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

Question

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

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

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

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

Answer 1

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

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

Answer 2

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

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

[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 3

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

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

Answer 4

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

Answer 5

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

Answer 6

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

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

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

-Адам

Answer 7

Простое старое исключение Гаусса сработало бы хорошо.

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

Answer 8

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