Rozwiązywanie

Muszę macierzy $A$ i $G$ . Jest niewielkie i z bardzo duża (może być rzędu kilku milionów.) oznacza wysoką osnowę z małe ( ), a każda kolumna może mieć tylko jeden wejście reszta to „s tak, że . jest ogromny, więc bardzo trudno odwrócić i mogę rozwiązać układ liniowy, taki jak $A$ $n\times n$ $n$ $G$ $n\times m$ $m$ $1 \lt m \lt 1000$ $1$ $0$ $G^TG = I$ $A$ $Ax = b$ iteracyjnie przy użyciu metody podprzestrzeni Kryłowa, takiej jak $\mathrm{BiCGStab}(l)$ , ale nie mamwyraźnie $A^{-1}$ .

Chcę rozwiązać układ postaci: $(G^TA^{-1}G)x = b$ , gdzie $x$ i $b$ są wektorami długości $m$ . Jednym ze sposobów jest wykorzystanie algorytmu iteracyjnego w algorytmie iteracyjnym do rozwiązania dla $A^{-1}$ 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.

linear-algebra linear-solver sparse-matrix krylov-method Costis
źródło

Właśnie dodałem do mojej odpowiedzi uwagę na temat wykorzystania struktury 0-1.

Arnold Neumaier

Odpowiedzi:

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

Edycja: Ponieważ jest złożoną symetryczną, podobnie jak powiększona macierz, ale nie ma quasidefinity. Można jednak użyć procedury do obliczenia ; 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. $A$ $Ax$ $A^*x=\overline{A\overline x}$

Edit2: Actually, the (0,1)-structure of $G$ means that you can eliminate $x$ amd the components of $G^Ty$ symbolically, thus ending up with a smaller system to solve. This means messing with the structure of $A$ , and pays only when $A$ 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

A $A$ ? 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

G $G$ 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

(GTA−HBA−1G)x=b $(G^TA^{-H}BA^{-1}G)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

GTA−HBA−1G $G^TA^{-H}BA^{-1}G$ 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

AB−1AH $AB^{-1}A^H$ . I don't think it would work simply with the pseudoinverse of B either.

Costis

Introduce

y:=A−1Gx $y:=A^{-1}Gx$ and

z:=A−HBy $z:=A^{-H}By$ , and proceed in analogy to the case worked out. (POssibly you also need to factor

B $B$ 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=b $Ax=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

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

Let us say that for simplicity that $G$ has $m$ columns and $n$ rows and that exactly the first $m$ rows have ones in them and that be reordering the elements of $x$ I can make it so that $G$ has the $m \times m$ identity matrix at the top and a $n-m \times m$ zero matrix at the bottom. Then I can partition $y=(y_c,y_f)$ into $m$ "constrained" and $n-m$ "free" elements so that $y_c=-b$ . I can also partition $A$ so that $A=\begin{pmatrix} A_{cc} & A_{cf} \\ A_{fc} & A_{ff} \end{pmatrix}$ . From the equation $Ay+Gx=0$ I then get the following:

A c c y c + A c f y f + x = 0, A f c y c + A f f y f = 0

$A_{cc} y_c + A_{cf} y_f + x = 0, \\ A_{fc} y_c + A_{ff} y_f = 0$ and using what we know about

yc $y_c$ we have from the second of these equations

A f f y f = A f c b

$A_{ff} y_f = A_{fc} b$ and consequently

x = A c c b - A c f A - 1 f f A f c b .

$x = A_{cc} b - A_{cf} A_{ff}^{-1} A_{fc} b.$ In other words, the only matrix you have to invert is the subset of

A $A$ whose rows and columns are not mentioned in

G $G$ (the null space of

G $G$ ). This you can easily do: (i) compute

z=Afcb $z=A_{fc} b$ ; (ii) use whatever solver you have to solve

Affh=z $A_{ff} h = z$ ; (iii) compute

x=Accb−Acfh $x = A_{cc} b - A_{cf} h$ .

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

Wolfgang Bangerth
źródło

But we know $G$ , $G^T$ and $A$ , so

$G^T A^{-1} G x = b$

$G G^T A^{-1} G x = G b$

Since $G^T G = I$ , then $G^T = G^{-1}$ , so $G G^T=I$ :

$A^{-1} G x = G b$

$A A^{-1} G x = A G b$

$G x = A G b$

$G^T G x = G^T A G b$

$x = G^T A G b$

Unless I've missed something, you don't need any iteration, or any solver to calculate x given $G$ , $A$ and $b$ .

Phil H
źródło

GT $G^T$ being a left inverse of

G $G$ does not imply that it is also a right inverse. Consider

$G=e_1$ , where

$G^T=e_1^T$ is a left inverse, but

$GG^T=e_1e_1^T\neq I$ .

Jack Poulson

$G\in \mathbb{C}^n\otimes \mathbb{C}^m$ , hence

$G^TG=I_{m\times m}$ , but

$GG^T\ne I_{n\times n}$ . Rather it's a projector on a subspace.

Deathbreath