Rozwiązywanie

22

Muszę macierzy AA i GG . Jest niewielkie i n x n z n bardzo duża (może być rzędu kilku milionów.) G oznacza n x m wysoką osnowę z m małe ( 1 < m < 1000 ), a każda kolumna może mieć tylko jeden 1 wejście reszta to 0 „s tak, że G T G = I . A jest ogromny, więc bardzo trudno odwrócić i mogę rozwiązać układ liniowy, taki jak A.An×nnGn×mm1<m<100010GTG=IAx = bAx=b iteracyjnie przy użyciu metody podprzestrzeni Kryłowa, takiej jak B i C G S t a b ( l )BiCGStab(l) , ale nie mamwyraźnie A - 1A1 .

Chcę rozwiązać układ postaci: ( G T A - 1 G ) x = b(GTA1G)x=b , gdzie xx i bb są wektorami długości mm . Jednym ze sposobów jest wykorzystanie algorytmu iteracyjnego w algorytmie iteracyjnym do rozwiązania dla A - 1A1 dla każdej iteracji zewnętrznego algorytmu iteracyjnego. Byłoby to jednak niezwykle kosztowne obliczeniowo. Zastanawiałem się, czy istnieje obliczeniowo łatwiejszy sposób rozwiązania tego problemu.

Costis
źródło
Właśnie dodałem do mojej odpowiedzi uwagę na temat wykorzystania struktury 0-1.
Arnold Neumaier

Odpowiedzi:

19

Wprowadź wektor y : = - A - 1 G x i rozwiąż duży układ sprzężony A y + G x = 0 , G T y = - b dla ( y , x ) jednocześnie, stosując metodę iteracyjną. Jeśli A jest symetryczny (co wydaje się prawdopodobne, choć nie podajesz tego wprost), to system jest symetryczny (ale nieokreślony, choć quasidefinite, jeśli Ay:=A1GxAy+Gx=0GTy=b(y,x)AAjest pozytywnie określony), co może pomóc ci wybrać odpowiednią metodę. (odpowiednie słowa kluczowe: macierz KKT, macierz quasidefinite).

Edycja: Ponieważ A jest złożoną symetryczną, podobnie jak powiększona macierz, ale nie ma quasidefinity. Można jednak użyć procedury A x do obliczenia A x = ¯ A ¯ x ; dlatego możesz dostosować metodę, taką jak QMR ftp://ftp.math.ucla.edu/pub/camreport/cam92-19.pdf (zaprojektowaną dla prawdziwych systemów, ale możesz łatwo przepisać ją dla złożonych systemów, używając funkcji łączenia w miejsce transpozycji), aby rozwiązać problem.AAxAx=Ax¯¯¯¯¯¯¯¯¯¯

Edit2: Actually, the (0,1)-structure of GG means that you can eliminate xx amd the components of GTyGTy symbolically, thus ending up with a smaller system to solve. This means messing with the structure of AA, and pays only when AA is given explicitly in sparse format rather than as a linear operator.

Arnold Neumaier
źródło
Thank you! A is complex symmetric. Is there reason to expect the condition of the augmented matrix to be worse than that of the original matrix AA? If m is small, the augmented matrix is only marginally larger in size than A, so I would suspect that solving this augmented system iteratively should not be much tougher than solving a system with A?
Costis
The condition number of the two systems is generally quite unrelated; it depends very much on what GG is. - I added to my answer information on how to exploit complex symmetry.
Arnold Neumaier
Hi guys! Thanks for all the replies; this place is great! An extension to the original question: Assume now that I have (GTAHBA1G)x=b(GTAHBA1G)x=b, where G and A have the same meaning as in the original question but B is a rank deficient nxn matrix (same size as A) and the whole GTAHBA1GGTAHBA1G is full rank. How would you go about solving the new system, since now the inverse of B does not exist so you cannot have AB1AHAB1AH. I don't think it would work simply with the pseudoinverse of B either.
Costis
1
Introduce y:=A1Gxy:=A1Gx and z:=AHByz:=AHBy, and proceed in analogy to the case worked out. (POssibly you also need to factor BB into full rank matrices and introduce an additional intermediate vector.)
Arnold Neumaier
Hi Arnold. Thanks, this indeed does work! I tested it with some very small test examples, and it works great. However, my iterative solver is having huge issues inverting the augmented matrix. While it takes only about 80 iterations (a few seconds) to solve a system of the form Ax=bAx=b with the original A matrix, the system with the augmented matrix (which is 2n+m x 2n+m or 2n-m x 2n-m using @wolfgang-bangerth 's approach) takes over tens of thousands of iterations (several hours) to solve for one RHS. Are there any strategies for accelerating the convergence? perhaps a preconditioner?
Costis
7

Following Arnold's reply, there is something you can do to simplify the problem. Specifically, rewrite the system as Ay+Gx=0,GTy=bAy+Gx=0,GTy=b. Then note that from the statement that GG is tall and narrow and each row has only one 1 and zeros otherwise, then the statement GTy=bGTy=b means that a subset of the elements of yy have a fixed value, namely the elements of bb.

Let us say that for simplicity that GG has mm columns and nn rows and that exactly the first mm rows have ones in them and that be reordering the elements of xx I can make it so that GG has the m×mm×m identity matrix at the top and a nm×mnm×m zero matrix at the bottom. Then I can partition y=(yc,yf)y=(yc,yf) into mm "constrained" and nmnm "free" elements so that yc=byc=b. I can also partition AA so that A=(AccAcfAfcAff)A=(AccAfcAcfAff). From the equation Ay+Gx=0Ay+Gx=0 I then get the following: Accyc+Acfyf+x=0,Afcyc+Affyf=0

Accyc+Acfyf+x=0,Afcyc+Affyf=0
and using what we know about ycyc we have from the second of these equations Affyf=Afcb
Affyf=Afcb
and consequently x=AccbAcfA1ffAfcb.
x=AccbAcfA1ffAfcb.
In other words, the only matrix you have to invert is the subset of AA whose rows and columns are not mentioned in GG (the null space of GG). This you can easily do: (i) compute z=Afcbz=Afcb; (ii) use whatever solver you have to solve Affh=zAffh=z; (iii) compute x=AccbAcfhx=AccbAcfh.

In other words, given the structure of GG, solving the linear system you have is really not more difficult than solving a single linear system with AA.

Wolfgang Bangerth
źródło
0

But we know GG, GTGT and AA, so

GTA1Gx=bGTA1Gx=b

GGTA1Gx=GbGGTA1Gx=Gb

Since GTG=IGTG=I, then GT=G1GT=G1, so GGT=IGGT=I:

A1Gx=GbA1Gx=Gb

AA1Gx=AGbAA1Gx=AGb

Gx=AGbGx=AGb

GTGx=GTAGbGTGx=GTAGb

x=GTAGbx=GTAGb

Unless I've missed something, you don't need any iteration, or any solver to calculate x given GG, AA and bb.

Phil H
źródło
3
GTGT being a left inverse of GG does not imply that it is also a right inverse. Consider G=e1, where GT=eT1 is a left inverse, but GGT=e1eT1I.
Jack Poulson
1
GCnCm, hence GTG=Im×m, but GGTIn×n. Rather it's a projector on a subspace.
Deathbreath