Solving $(G^TA^{-1}G)x = b$ without inverting $A$

I have matrices $A$ and $G$. $A$ is sparse and is $n\times n$ with $n$ very large (can be on the order of several million.) $G$ is an $n\times m$ tall matrix with $m$ rather small ($1 \lt m \lt 1000$) and each column can only have a single $1$ entry with the rest being $0$'s, such that $G^TG = I$. $A$ is huge, so it is really tough to invert, and I can solve a linear system such as $Ax = b$ iteratively using a Krylov subspace method such as $\mathrm{BiCGStab}(l)$, but I do not have $A^{-1}$ explicitly.

I want to solve a system of the form: $(G^TA^{-1}G)x = b$, where $x$ and $b$ are $m$ length vectors. One way to do it is to use an iterative algorithm within an iterative algorithm to solve for $A^{-1}$ for each iteration of the outer iterative algorithm. This would be extremely computationally expensive, however. I was wondering if there is a computationally easier way to go about solving this problem.

• I just added to my answer a remark on exploiting the 0-1 structure. – Arnold Neumaier May 16 '12 at 19:46

Introduce the vector $y:=-A^{-1}Gx$ and solve the large coupled system $Ay+Gx=0$, $G^Ty=-b$ for $(y,x)$ simultaneously, using an iterative method. If $A$ is symmetric (as seems likely though you don't state it explicitly) then the system is symmetric (but indefinite, though quasidefinite if $A$ is positive definite), which might help you to choose an appropriate method. (relevant keywords: KKT matrix, quasidefinite matrix).

Edit: As $A$ is complex symmetric, so is the augmented matrix, but there is no quasidefiniteness. You can however use the $Ax$ routine to compute $A^*x=\overline{A\overline x}$; therefore you could adapt a method such as QMR ftp://ftp.math.ucla.edu/pub/camreport/cam92-19.pdf (designed for real systems, but you can easily rewrite it for complex systems, using the adjoint in place of the transpose) to solve your problem.

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.

• 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$? 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 May 16 '12 at 9:25
• The condition number of the two systems is generally quite unrelated; it depends very much on what $G$ is. - I added to my answer information on how to exploit complex symmetry. – Arnold Neumaier May 16 '12 at 10:02
• Hi guys! Thanks for all the replies; this place is great! An extension to the original question: Assume now that I have $(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 $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^{-1}A^H$. I don't think it would work simply with the pseudoinverse of B either. – Costis May 16 '12 at 21:58
• Introduce $y:=A^{-1}Gx$ and $z:=A^{-H}By$, and proceed in analogy to the case worked out. (POssibly you also need to factor $B$ into full rank matrices and introduce an additional intermediate vector.) – Arnold Neumaier May 17 '12 at 12:30
• 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$ 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 May 19 '12 at 8:30

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_{cc} y_c + A_{cf} y_f + x = 0, \\ A_{fc} y_c + A_{ff} y_f = 0$$ and using what we know about $y_c$ we have from the second of these equations $$A_{ff} y_f = A_{fc} b$$ and consequently $$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$ whose rows and columns are not mentioned in $G$ (the null space of $G$). This you can easily do: (i) compute $z=A_{fc} b$; (ii) use whatever solver you have to solve $A_{ff} h = z$; (iii) compute $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$.

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

• $G^T$ being a left inverse of $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 May 16 '12 at 15:48
• $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 May 16 '12 at 17:32