численно устойчивая обратная матрица 2x2

Question

В числовом решателе, над которым я работаю в C, мне нужно инвертировать матрицу 2x2, а затем умножить ее справа на другую матрицу:

C = B . inv(A)

Я использовал следующее определение перевернутой матрицы 2x2:

a = A[0][0];
b = A[0][1];
c = A[1][0];
d = A[1][1];
invA[0][0] = d/(a*d-b*c);
invA[0][1] = -b/(a*d-b*c);
invA[1][0] = -c/(a*d-b*c);
invA[1][1] = a/(a*d-b*c);

В первые несколько итераций моего решателя это, кажется, дает правильные ответы, однако после нескольких шагов все начинает расти и в конечном итоге взрывается.

Теперь, сравнивая с реализацией, использующей SciPy, я обнаружил, что та же математика не взрывается. Единственное отличие, которое я могу найти, состоит в том, что в коде SciPy используется scipy.linalg.inv(), который внутренне использует LAPACK для выполнения инверсии.

Когда я заменяю вызов inv() приведенными выше вычислениями, версия Python действительно взрывается, так что я почти уверен, что это проблема. Наблюдаются небольшие различия в вычислениях, что заставляет меня полагать, что это численная проблема, что неудивительно для операции инверсии.

Я использую поплавки двойной точности (64-разрядные), надеясь, что числовые проблемы не будут проблемой, но, очевидно, это не так.

Но: я хотел бы решить эту проблему в своем C-коде, не обращаясь к библиотеке, подобной LAPACK, потому что вся причина переноса его на чистый C - запустить его на целевой системе. Более того, я хотел бы понять проблему, а не просто вызвать черный ящик. В конце концов, я бы хотел, чтобы он работал с одинарной точностью, если это возможно.

Итак, у меня такой вопрос: для такой маленькой матрицы есть численно более устойчивый способ вычисления обратного значения A?

Спасибо.

Изменить: В настоящее время пытается выяснить, могу ли я просто избежать инверсии , решив для C.

Answer 1

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

Сказать C = B . inv(A) - это то же самое, что сказать, что вы хотите решить AC = B для C. Вы можете сделать это, разделивкаждый B и C в два столбца.Решая A C1 = B1 и A C2 = B2, мы получим C.

Answer 2

Ваш код в порядке;однако, это рискует потерю точности из-за любого из четырех вычитаний.

Рассмотрите возможность использования более продвинутых методов, таких как matfunc.py .Этот код выполняет инверсию, используя QR-декомпозицию , реализованную с отражениями домохозяина .Результат дополнительно улучшен с помощью итеративного уточнения .

Answer 3

Вычисление определителя не является стабильным. Лучший способ - использовать Gauss-Jordan с частичным поворотом, который вы можете легко легко здесь проработать.

Решение системы 2x2

Давайте решим систему (используйте c, f = 1, 0, затем c, f = 0, 1, чтобы получить обратное)

a * x + b * y = c
d * x + e * y = f

В псевдокоде это выглядит как

if a == 0 and d == 0 then "singular"

if abs(a) >= abs(d):
    alpha <- d / a
    beta <- e - b * alpha
    if beta == 0 then "singular"
    gamma <- f - c * alpha
    y <- gamma / beta
    x <- (c - b * y) / a
else
    swap((a, b, c), (d, e, f))
    restart

Это стабильнее, чем определитель + коматрица (beta - определитель * некоторая константа, вычисленная устойчивым образом). Вы можете получить полный эквивалент поворота (т. Е. Потенциально поменять местами x и y, чтобы первое деление на a было таким, чтобы a было наибольшим числом по величине среди a, b, d, e), и это может быть более стабильным в некоторых обстоятельствах, но вышеописанный метод работал хорошо для меня.

Это эквивалентно выполнению декомпозиции LU (сохранить гамму, бета, a, b, c, если вы хотите сохранить эту декомпозицию LU).

Вычисление QR-разложения также может быть выполнено явно (и также очень стабильно при условии, что вы делаете это правильно), но это медленнее (и включает в себя получение квадратных корней). Выбор за вами.

Повышение точности

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

Действительно, предположим, что вы решили A * x = b для x с помощью вышеуказанного метода. Теперь вы вычисляете A * x и обнаруживаете, что оно не совсем равно b, что есть небольшая ошибка:

A * x - b = db

Теперь, если вы решите для dx в A * dx = db, у вас будет

A * (x - dx) = b + db - db - ddb = b - ddb

, где ddb - ошибка, вызванная численным решением A * dx = db, которое обычно намного меньше, чем db (поскольку db намного меньше, чем b).

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

Answer 4

Я согласен с Жан-Викотром в том, что вам, вероятно, следует использовать метод Якобба. Вот мой пример:

#Helper functions:
def check_zeros(A,I,row, col=0):
"""
returns recursively the next non zero matrix row A[i]
"""
if A[row, col] != 0:
    return row
else:
    if row+1 == len(A):
        return "The Determinant is Zero"
    return check_zeros(A,I,row+1, col)

def swap_rows(M,I,row,index):
"""
swaps two rows in a matrix
"""
swap = M[row].copy()
M[row], M[index] = M[index], swap
swap = I[row].copy()
I[row], I[index] = I[index], swap

# Your Matrix M
M = np.array([[0,1,5,2],[0,4,9,23],[5,4,3,5],[2,3,1,5]], dtype=float)
I = np.identity(len(M))

M_copy = M.copy()
rows = len(M)

for i in range(rows):
index =check_zeros(M,I,i,i)
while index>i:
    swap_rows(M, I, i, index)
    print "swaped"
    index =check_zeros(M,I,i,i) 

I[i]=I[i]/M[i,i]
M[i]=M[i]/M[i,i]   

for j in range(rows):
    if j !=i:
        I[j] = I[j] - I[i]*M[j,i]
        M[j] = M[j] - M[i]*M[j,i]
print M
print I  #The Inverse Matrix

Answer 5

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