Inverse of $AQ^{-1}A'$

Question

Suppose that $A$ is a $m\times n$ full row rank sparse matrix, and $Q$ is an $n\times n$ symmetric positive definite sparse matrix with $m<n$. Besides, $m$ is about $10^5$, and $n$ is about $10^6$. There is no other special structure for $A$ and $Q$ (i.e., not circulant, not Toeplitz, etc.). I need to solve the linear system $AQ^{-1}A'X=B$, where $B$ is an $m\times m$ matrix. Is there any way to find the explicit formula for $(AQ^{-1}A')^{-1}$?

BTW, calculation of cholesky decomposition for $Q$ is not a problem on my laptop, but computing $Q^{-1}A'=L'\backslash(L\backslash A')$ breaks down with out of memory errors in MATLAB. If I compute $Q^{-1}A'$ with parfor loop in MATLAB, it is extremely slow. Direct solvers for these linear systems are preferred.

Thanks, Pulong

Answer 1

Having tried several variants myself, I'm convinced that direct methods are not feasible for the typical laptop computer and systems of this size.

You could try the following:

1. Compute $Y = Q^{-1}A'$ via CG (conjugate gradients). Matlab has a built-in version, pcg. To deal with the multiple right-hand sides $A'$, you can either use a for-loop (which is slow), or see what they do in A CG Method for Multiple Right Hand Sides and Multiple Shifts in Lattice QCD Calculations.
2. Then compute $(AY)^{-1}B$ via the backslash operator or an LU factorization. $AY$ will be full. In fact, this is the inherent issue in your problem: you want the inverse of a full matrix $AQ^{-1}A'$. Therefore, you cannot make use of sparse routines at this point.

It's possible step 2 won't work-- for instance, my laptop cannot even store a full array of dimension $10000 \times 10000$. In which case, you may want to run this problem on a high-performance computer, if your work or campus has one.

Update:

You can reformulate your problem as the computation of a bilinear form of a matrix function, $W^Tf(Q)W$, where $W = A^T$, and $f(z) = z^{-1}$. This form relates to a Riemann-Stieltjes integral, which can be approximated by certain quadrature rules. For how to express $z^{-1}$ as a Riemann-Stieltjes integral, see my post How Can One Write $ {z}^{−1} $ as a Stieltjes Function?. I only recently stumbled on this topic myself via the following sources:

1. Chapter 7 in Matrices, Moments and Quadrature with Applications by Golub and Meurant.
2. Matrices, moments and quadrature II; How to compute the norm of the error in iterative methods also by Golub and Meaurant.