que - transpuesta de una matriz en python
numéricamente estable inversa de una matriz de 2x2 (5)
No invierta la matriz. Casi siempre, lo que está utilizando lo contrario para lograr se puede hacer más rápido y con mayor precisión sin invertir la matriz. La inversión de la matriz es intrínsecamente inestable, y mezclar eso con números en coma flotante es buscar problemas.
Diciendo C = B . inv(A) C = B . inv(A) es lo mismo que decir que quieres resolver AC = B para C. Puedes lograr esto dividiendo cada B y C en dos columnas. Al resolver A C1 = B1 y A C2 = B2 se producirá C.
En un solucionador numérico estoy trabajando en C, necesito invertir una matriz de 2x2 y luego se multiplica en el lado derecho por otra matriz:
C = B . inv(A)
He estado usando la siguiente definición de una matriz invertida de 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);
En las primeras iteraciones de mi solucionador, esto parece dar las respuestas correctas, sin embargo, después de unos pocos pasos las cosas comienzan a crecer y finalmente explotan.
Ahora, comparando una implementación usando SciPy, encontré que las mismas matemáticas no explotan. La única diferencia que puedo encontrar es que el código SciPy usa scipy.linalg.inv() , que internamente usa LAPACK internamente para realizar la inversión.
Cuando reemplazo la llamada a inv() con los cálculos anteriores, la versión de Python explota, así que estoy bastante seguro de que este es el problema. Se están introduciendo pequeñas diferencias en los cálculos, lo que me lleva a pensar que se trata de un problema numérico, lo cual no es del todo sorprendente para una operación de inversión.
Estoy usando flotadores de doble precisión (64 bits), con la esperanza de que los problemas numéricos no sean un problema, pero aparentemente ese no es el caso.
Pero: me gustaría resolver esto en mi código C sin necesidad de llamar a una biblioteca como LAPACK, porque la razón para portarlo a C puro es hacerlo funcionar en un sistema de destino. Además, me gustaría entender el problema, no solo llamar a una caja negra. Finalmente, me gustaría ejecutarlo con precisión simple también, si es posible.
Entonces, mi pregunta es, para una matriz tan pequeña, ¿hay una forma numéricamente más estable para calcular la inversa de A?
Gracias.
Editar: Actualmente estoy tratando de averiguar si puedo evitar la inversión resolviendo C
Tu código está bien; sin embargo, corre el riesgo de pérdida de precisión de cualquiera de las cuatro restas.
Considere el uso de técnicas más avanzadas, como la utilizada en matfunc.py . Ese código realiza la inversión usando una descomposición QR implementada con reflexiones de Householder . El resultado se mejora aún más con el refinamiento iterativo .
Use el método de Jacobi, que es un método iterativo que implica "invertir" solo la diagonal principal de A, que es muy directa y menos propensa a la inestabilidad numérica que la inversión de toda la matriz.
La computación del determinante no es estable. Una mejor manera es usar Gauss-Jordan con pivote parcial, que aquí se puede resolver de manera explícita fácilmente.
Resolviendo un sistema 2x2
Vamos a resolver el sistema (use c, f = 1, 0 luego c, f = 0, 1 para obtener el inverso)
a * x + b * y = c
d * x + e * y = f
En pseudo código, esto dice
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
Esto es más estable que determinante + comatrix ( beta es el determinante * alguna constante, calculada de manera estable). Puede calcular el equivalente de pivote completo (es decir, potencialmente intercambiando xey, de modo que la primera división por a sea tal que a sea el número más grande en magnitud entre a, b, d, e), y esto puede ser más estable en algunas circunstancias, pero el método anterior ha funcionado bien para mí.
Esto es equivalente a realizar la descomposición de LU (almacenar gamma, beta, a, b, c si desea almacenar esta descomposición de LU).
La computación de la descomposición QR también se puede hacer explícitamente (y también es muy estable siempre que lo haga correctamente), pero es más lenta (e implica tomar raíces cuadradas). La decisión es tuya.
Mejorando la precisión
Si necesita una mayor precisión (el método anterior es estable, pero hay algún error de redondeo, proporcional a la proporción de los valores propios), puede "resolver la corrección".
De hecho, supongamos que resuelve A * x = b para x con el método anterior. Ahora calcula A * x , y encuentra que no es igual a b , que hay un error leve:
A * x - b = db
Ahora, si resuelve para dx en A * dx = db , tiene
A * (x - dx) = b + db - db - ddb = b - ddb
donde ddb es el error inducido por la resolución numérica de A * dx = db , que es típicamente mucho más pequeño que db (dado que db es mucho más pequeño que b ).
Puede repetir el procedimiento anterior, pero típicamente se requiere un paso para restaurar la precisión total de la máquina.
Estoy de acuerdo con Jean-Vicotr en que probablemente debas usar el método Jacobbiano. Aquí está mi ejemplo:
#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