Resolviendo

22

Tengo matrices y G . A es escaso y es n × n con n muy grande (puede ser del orden de varios millones). G es una matriz alta de n × m con m bastante pequeño ( 1 < m < 1000 ) y cada columna solo puede tener una sola 1 entrada con el resto siendo 0 's, de tal manera que G T G = I . A es enorme, por lo que es muy difícil invertir, y puedo resolver un sistema lineal como AAGAn×nnGn×mm1<m<100010GTG=IA usando iterativamente un método de subespacio de Krylov como B i C G S t a b ( l ) , pero no tengo A - 1 explícitamente.Ax=bBiCGStab(l)A1

Quiero resolver un sistema de la forma: , donde x y b son vectores de longitud m . Una forma de hacerlo es usar un algoritmo iterativo dentro de un algoritmo iterativo para resolver A - 1 para cada iteración del algoritmo iterativo externo. Sin embargo, esto sería extremadamente costoso computacionalmente. Me preguntaba si hay una manera computacionalmente más fácil de resolver este problema.(GTA1G)x=bxbmA1

Costis
fuente
Acabo de agregar a mi respuesta un comentario sobre la explotación de la estructura 0-1.
Arnold Neumaier

Respuestas:

19

Introduzca el vector y resuelva el sistema acoplado grande A y + G x = 0 , G T y = - b para ( y , x ) simultáneamente, utilizando un método iterativo. Si A es simétrico (como parece probable aunque no lo diga explícitamente), entonces el sistema es simétrico (pero indefinido, aunque cuasidefinito si Ay:=A1GxAy+Gx=0GTy=b(y,x)AAes positivo definido), lo que podría ayudarlo a elegir un método apropiado. (palabras clave relevantes: matriz KKT, matriz cuasidefinita).

Editar: como es simétrico complejo, también lo es la matriz aumentada, pero no hay cuasidefinidad. Sin embargo, puede utilizar la rutina A x para calcular A x = ¯ A ¯ x ; por lo tanto, podría adaptar un método como QMR ftp://ftp.math.ucla.edu/pub/camreport/cam92-19.pdf (diseñado para sistemas reales, pero puede reescribirlo fácilmente para sistemas complejos, utilizando el adjunto en lugar de la transposición) para resolver su problema.AAxAx=Ax¯¯

Edit2: En realidad, la estructura (0,1) de significa que puede eliminar x y los componentes de G T y simbólicamente, terminando así con un sistema más pequeño para resolver. Esto significa jugar con la estructura de A , y paga solo cuando A se da explícitamente en formato disperso en lugar de como un operador lineal.GxGTyAA

Arnold Neumaier
fuente
¡Gracias! A es complejo simétrico. ¿Hay alguna razón para esperar que la condición de la matriz aumentada sea peor que la de la matriz original ? Si m es pequeño, la matriz aumentada tiene un tamaño marginalmente mayor que A, por lo que sospecho que resolver este sistema aumentado de forma iterativa no debería ser mucho más difícil que resolver un sistema con A. A
Costis
El número de condición de los dos sistemas generalmente no está relacionado; depende mucho de lo que es - Agregué a mi respuesta información sobre cómo explotar la simetría compleja. G
Arnold Neumaier
¡Hola chicos! Gracias por todas las respuestas; Este lugar es genial! Una extensión de la pregunta original: supongamos ahora que tengo , donde G y A tienen el mismo significado que en la pregunta original pero B es una matriz nxn deficiente en rango ( mismo tamaño que A) y todo G T A - H B A - 1 G es rango completo. ¿Cómo resolvería el nuevo sistema, ya que ahora no existe el inverso de B, por lo que no puede tener A B - 1 A H(GTAHBA1G)x=bGTAHBA1GAB1AH. No creo que funcione simplemente con el pseudoinverso de B tampoco.
Costis
1
Introduzca y z : = A - H B y , y proceda en analogía con el caso resuelto. (Posiblemente también necesite factorizar B en matrices de rango completo e introducir un vector intermedio adicional).y:=A1Gxz:=AHByB
Arnold Neumaier
Hola arnold Gracias, esto realmente funciona! Lo probé con algunos ejemplos de prueba muy pequeños, y funciona muy bien. Sin embargo, mi solucionador iterativo está teniendo grandes problemas para invertir la matriz aumentada. Si bien solo se necesitan unas 80 iteraciones (unos segundos) para resolver un sistema de la forma con la matriz A original, el sistema con la matriz aumentada (que es 2n + mx 2n + mo 2n-mx 2n- m usando el enfoque de @ wolfgang-bangerth) toma decenas de miles de iteraciones (varias horas) para resolver un RHS. ¿Hay alguna estrategia para acelerar la convergencia? tal vez un preacondicionador? Ax=b
Costis
7

Después de la respuesta de Arnold, hay algo que puede hacer para simplificar el problema. Específicamente, reescriba el sistema como . Luego, tenga en cuenta que, a partir de la afirmación de que G es alto y angosto y que cada fila tiene solo un 1 y ceros, de lo contrario, la afirmación G T y = - b significa que un subconjunto de los elementos de y tiene un valor fijo, es decir, los elementos de - b .Ay+Gx=0,GTy=bGGTy=byb

Digamos que para simplificar que tiene m columnas yn filas y que exactamente las primeras m filas tienen unas y que reordenando los elementos de x , puedo hacer que G tenga la matriz de identidad m × m en la parte superior y una matriz n - m × m cero en la parte inferior. Entonces puedo dividir y = ( y c , y f ) en m elementos "restringidos" y n - m "libres" para queGmnmxGm×mnm×my=(yc,yf)mnm . También puedo dividir A para que A = ( A c c A c f A f c A f f ) . De la ecuación A y + G x = 0 , obtengo lo siguiente: A c c y c + A c f y f + x = 0 ,yc=bAA=(AccAcfAfcAff)Ay+Gx=0 y usando lo que sabemos sobre y c que tenemos de la segunda de estas ecuaciones A f f y f = A f c b y consecuentemente x = A c c b - A c f A - 1 f f A f c b . En otras palabras, la única matriz que tiene que invertir es el subconjunto de A

Accyc+Acfyf+x=0,Afcyc+Affyf=0
yc
Affyf=Afcb
x=AccbAcfAff1Afcb.
Acuyas filas y columnas no se mencionan en (el espacio nulo de G ). Esto se puede hacer fácilmente: (i) calcular z = A f c b ; (ii) use cualquier solucionador que tenga para resolver A f f h = z ; (iii) calcule x = A c c b - A c f h .GGz=AfcbAffh=zx=AccbAcfh

En otras palabras, dada la estructura de , resolviendo el sistema lineal que tiene en realidad no es más difícil que la solución de un solo sistema lineal con una .GA

Wolfgang Bangerth
fuente
0

GGTA

GTA1Gx=b

GGTA1Gx=Gb

GTG=IGT=G1GGT=I

A1Gx=Gb

AA1Gx=AGb

Gx=AGb

GTGx=GTAGb

x=GTAGb

GAb

Phil H
fuente
3
GTGG=e1GT=e1TGGT=e1e1TI
1
GCnCmGTG=Im×mGGTIn×n