численно стабильной обратной матрицы 2х2

Question

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

C = B . inv(A) 

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

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?

Спасибо.

Редактировать: В настоящее время пытается выяснить, могу ли я только avoid the inversion путем решения для C.

Answer 1

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

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

Answer 2

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

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

Answer 3

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

Answer 4

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

Решая систему 2х2

Решим систему (использование в, е = 1, 0, то с, е = 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 

Это стабильнее, чем детерминант + comatrix (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 5

Я согласен с Jean-Vicotr в том, что вы, вероятно, должны использовать метод Jacobbian. Вот мой пример:

#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