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



Inverting products of matrices
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!


