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

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

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

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

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

Ответ 1

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

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

Ответ 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-матрицы украдены со связанной страницы.)

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

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

Ответ 3

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

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

Ответ 4

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

Ответ 5

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

Ответ 6

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

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

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

-Adam

Ответ 7

Обычная старая гауссова элиминация будет работать хорошо.

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

Ответ 8

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