Inverting products of matrices

Question

I need to compute a large number of inverses of the following form:

$(A \Lambda_k A^\top)^{-1}$

where $A \in \mathbb{R}^{m \times n}$, $n > m$ and $\Lambda_k = \text{diag}(\lambda_1, ..., \lambda_n)$ with $\lambda_i > 0$. Is there an efficient way to do this?

In the end, I want to sample from Gaussians with means $(A \Lambda_k A^\top)^{-1} \mu_k$ and covariances $(A\Lambda_kA)^{-1}$.

I actually can choose $A$ freely, as long as the rows of $A$ are orthogonal to the rows of a matrix $B \in \mathbb{R}^{(n - m) \times n}$, i.e. $A \cdot B^\top = 0$. Oh, and $\text{rank}(A) + \text{rank}(B) = n$.

Answer 1

Now that you've provided some more information, I think I can make some useful suggestions.

First, a quick review of linear transformations of multivariate normal random vectors. If $z$ is an MVN vector with mean $\mu$ and covariance matrix $C$, and $M$ is a matrix of the appropriate size and $y=Mz$, then $y$ is MVN with $E[y]=M\mu$ and $Cov(y)=MCM^{T}$.

This is a key fact that can be very useful in algorithms for generating MVN random numbers with desired mean and covariance. If you want to generate an MVN vector $x$ with mean $\mu$ and covariance $C$, then you can do this by

1. Compute the Cholesky factorization of $C$, $C=R^{T}R$.
2. Let $z$ be an N(0,I) random vector.
3. Let $x=R^{T}z+\mu$.
4. Then $E[x]=R^{T}E[z]+\mu=\mu$ and $Cov(x)=R^{T}IR=C$.

You could simply apply this algorithm to your problem by computing the Cholesky factorization of the covariance matrix:

$(A\Lambda_{k} A^{T})=R^{T}R$.

Then

$(A\Lambda_{k} A^{T})^{-1}=R^{-1}R^{-T}$.

Then you could generate the desired random vector $x$ with

$x=R^{-1}z+R^{-1}R^{-T}\mu_{k}$.

Computing the matrix $A\Lambda_{k} A^{T}$ takes $O(m^2n)$ time. Computing the Cholesky factorization takes $O(m^3)$ time, with a somewhat larger constant factor. Computing the inverse of $R$ can be done quickly by backsolving, but still takes $O(m^3)$ time. However, if $n$ is much bigger than $m$, you'll end up spending most of your time computing $A\Lambda_{k} A^{T}$.

Because of the orthogonality constraints on the rows of $A$, you can assume that $A$ and $A\Lambda_{k} A^{T}$ will be fully dense. Thus there doesn't appear to be anything you can gain here by using sparse matrix techniques.

My first version of this answer suggested using the $QR$ factorization of $A$. This would leave you with

$A\Lambda A^{T}=QR\Lambda_{k} R^{T}Q$

where $Q$ is an $m$ by $m$ orthogonal matrix and $R$ is an $m$ by $n$ upper triangular matrix. Unfortunately, you'd then have to go on and compute a Cholesky factorization of $R\Lambda_{k} R^{T}$, which is just as much work as computing the Cholesky factorization of $A\Lambda_{k} A^{T}$. So I don't think the QR factorization is worth while after all. However, if it happened that $m=n$, then the QR factorization approach would be very helpful!

Answer 2

You haven't really told us much about the problem.

Are you working in conventional single or double precision floating point arithmetic, or are you working in some obscure field?

Are the elements $\lambda_{k}$ positive, or could some of them be negative?

Are your matrices $A$ sparse, or are they fully dense?

Why are you computing these inverses- do you actually need some or all of the elements of the inverse matrix, or are you just going to be multiplying the inverse times one or more vectors once you get it? If that's the case, then you should probably not be computing the inverse but rather computing a matrix factorization and using the factorization to solve systems of equations that are equivalent to multiplication by the inverse.

Systems of equations of the form $ADA^{T}x=b$ where $D$ is a diagonal matrix with positive elements arise frequently in interior point methods for linear programming (as well as interior point methods for SOCP and SDP.) In that application it turns out that sparse Cholesky factorization with matrix reordering to improve the sparsity of the Cholesky factors is the way to go.