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

Question

I have matrices A and G. A is sparse and is nxn with n very large (can be on the order of several million.) G is an nxm tall matrix with m rather small (1 < m < 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 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.

Answer 1

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.

Answer 2

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

Answer 3

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