1
$\begingroup$

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

$\endgroup$
1
  • $\begingroup$ By $A'$, do you mean the transpose? Also, why do you prefer direct methods? If $A$ is sparse, and if you're having memory issues, then Krylov methods could be advantageous. $\endgroup$
    – artificial_moonlet
    Jun 10, 2015 at 19:48

1 Answer 1

Reset to default
0
$\begingroup$

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.
$\endgroup$
4
  • $\begingroup$ 1. Yes, $A'$ standards for the transpose of $A$. 2. $A$ is full-row rank, and it only has explicit formula for its right inverse, say, $A^+$. Then $A^+=A'(AA')^{-1}$, see en.wikipedia.org/wiki/Moore–Penrose_pseudoinverse . $\endgroup$
    – Bayes
    Jun 11, 2015 at 4:02
  • $\begingroup$ Besides, I also have a new idea, but it also does work, seelink. 1. [R1, p1, s1] = chol(Q, 'vector'); 2. R1A(s1, :) = R1'\A(:,s1)'; % this step most of the time fails for some Q, with my Mac frozen. 3. X1 = factorize(R1A'*R1A)\B; % factorize is a function in SuiteSparse package. $\endgroup$
    – Bayes
    Jun 11, 2015 at 4:11
  • $\begingroup$ You are correct: the left inverses do not exist. I will either modify or delete my answer. $\endgroup$
    – artificial_moonlet
    Jun 14, 2015 at 0:40
  • $\begingroup$ I think it is better to always write the name of the sources you link to instead of Here or There. Someday the links won't be valid, the names will be searchable. $\endgroup$
    – Royi
    Apr 1, 2020 at 13:51

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .