Is there a convenient form for this large covariance matrix? Consider the following bivariate vector autoregression: $$X_t=\mu +X_{t-1}A+\varepsilon_t,\ \varepsilon_t \overset{iid}{\sim} MVN(0, V),\ X_t=(X_{1,t},X_{2,t})',$$ 
where the assumptions on the coefficient matrix $A$ are such that the process $\{X_t\}$ is stationary, i.e. there are no unit roots. 
The goal is to find a nice (from a computational perspective) expression for the inverse covariance matrix for the vector $y=(X_{1,1},X_{1,2},\dots , X_{1,T},  X_{2,1} \dots X_{2,T})'$.
In the univariate case, i.e. when $\{X_t\}$ is an AR(1) process such as:
$$
X_t=\mu + aX_{t-1}+\varepsilon_t,\ \varepsilon_t \sim N(0, v)
$$ 
then Chen and Deo (2009) gives the following convenient expression: $\Omega^{-1}=B'B,$ where
$$
B=
\ \left( \begin{array}{ccc}
\sqrt{1-a^2} & 0 & 0 & \dots & 0\\
-a & 1 & 0 & \dots & 0\\
0 & -a & 1  &\dots & 0\\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & -a & 1
\end{array} \right)\ 
$$
and $v\Omega$ is the $T\times T$ covariance matrix for $(X_1,X_2,\dots,X_T)'$.
If it helps, and I think it does, we may assume that $A$ is diagonal so that the bivariate autoregression is really two separate AR(1)-processes with possibly correlated innovations.
The problem is research relate and I need to to speed up numerical computations of bilinear forms such as $w'\Sigma^{-1} u$, for some vectors $w$ and $u$, where $\Sigma$ is the covariance matrix of $y$.
I have a feeling this may be a standard problem solved in some old paper -- I was indeed expecting to find something in a textbook even -- but I can't find any good references. All help, even if just hints and not full answers, is greatly appreciated.
Update
In the case with diagonal $A$, all the elements of $\Sigma$ are easily found. We have $$\mathrm{Cov}(X_{i,t},X_{j,t+h})=A_{j,j}^h \frac{\mathrm{Cov}(\varepsilon_{i,t},\varepsilon_{j,t})}{1-A_{i,i}A_{j,j}}=A_{j,j}^h \frac{V_{i,j}}{1-A_{i,i}A_{j,j}}.$$ Note that the fraction only takes 3 different values, one for $i=j=1$, one for $i=j=2$ and one for $i\neq j$. These values corresponds to different blocks of $\Sigma$. Each of these blocks' inverse can be found directly using Chen and Deo's expression (I think).
Another observation is that the power of the $A_{j,j}$ term is just incremented by $\pm 1$ as we move up/down or left/right within these blocks of the matrix. 
 A: It looks like that one should assume that the system has already evolved for a sufficiently long time and reached its stationary point. Then, one considers $T-1$ steps of subsequent evolution and the covariance in this interval.
Hence, let us first find the stationary point. Suppose that we started with an initial value $\tilde{X}_{1}$. Iterating the autoregression gives
\begin{equation}
\tilde{X}_{t} = (I-A)^{-1}(I-A^{t-1})\mu + A^{t-1}\tilde{X}_{1} +\sum_{s=0}^{t-2} A^{s}\varepsilon_{t-s}
\end{equation}
Assuming that this is a stationary process, large powers of $A$ can be ignored, and we obtain the following result:
\begin{equation}
\begin{split}
\tilde{X}_{\infty} &= (I-A)^{-1}\mu + \sum_{s=0}^{\infty} A^{s}\varepsilon_{1-s}\\
&\equiv X_{1}.
\end{split}
\end{equation}
This will be our starting point for further evolution. From the fact that $\varepsilon_{t}$ are i.i.d. $\rm{MVN}(0,V)$, we can deduce that
\begin{equation}
\rm{E}(X_{1}) = (I-A)^{-1}\mu,
\end{equation}
\begin{equation}
\rm{Cov}(X_{1}) = \sum_{s=0}^{\infty} A^{s}V A^{\prime\,s} \equiv V_{1}.
\end{equation}
The covariance $V_{1}$ satisfies the following discrete Lyapunov equation:
\begin{equation}
AV_{1}A^{\prime} - V_{1} +V = 0.
\end{equation}
Looking back at the equation at the top, it is obvious that subsequent evolution from $X_{1}$ is simply given by
\begin{equation}
X_{t} = (I-A)^{-1}(I-A^{t-1})\mu + A^{t-1}X_{1} +\sum_{s=0}^{t-2} A^{s}\varepsilon_{t-s}
\end{equation}
As we are ultimately concerned with covariance, considering the deviation makes life easier:
\begin{equation}
\delta X_{t}\equiv X_{t} - \rm{E}(X_{t}),
\end{equation}
whose evolution is as follows:
\begin{equation}
\delta X_{t} = A^{t-1}\delta X_{1} +\sum_{s=0}^{t-2} A^{s}\varepsilon_{t-s}.
\end{equation}
Notice that $\delta X_{1}$, $\epsilon_{2}$, $\epsilon_{3}$, ... are independent (vector) variables of zero mean and covariance $V_{1}$, $V$, $V$, ....
That is, the covariance matrix of $\boldsymbol{\varepsilon} \equiv (\delta X_{1},\varepsilon_{2},\ldots,\varepsilon_{T})^{\prime}$ is
\begin{equation}
\textbf{V} \equiv \left(\begin{array}{cccccc} V_{1}&0&0&\cdots&0&0\\ 0&V&0&\cdots&0&0\\ 0&0&V&\cdots&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\0&0&0&\cdots&V&0\\0&0&0&\cdots&0&V\end{array}\right).
\end{equation}
We need to turn this into the covariance of $\delta\textbf{X} \equiv (\delta X_{1},\delta X_{2},\ldots,\delta X_{T})^{\prime}$. To do this, we note that
\begin{equation}
\delta\textbf{X} = \textbf{M} \boldsymbol{\varepsilon},
\end{equation}
where
\begin{equation}
\textbf{M} \equiv \left(\begin{array}{cccccc} I&0&0&\cdots&0&0\\ A&I&0&\cdots&0&0\\ A^2&A&I&\cdots&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\A^{T-2}&A^{T-3}&A^{T-4}&\cdots&I&0\\A^{T-1}&A^{T-2}&A^{T-3}&\cdots&A&I\end{array}\right).
\end{equation}
This follows from the expression for $\delta X_{t}$ given above. Therefore, the covariance matrix of the transformed variable $\delta\textbf{X}$, and hence that of $\textbf{X} \equiv (X_{1},X_{2},\ldots,X_{T})^{\prime}$ itself, is given by
\begin{equation}
\textbf{W} = \textbf{M}\textbf{V}\textbf{M}^{\prime}.
\end{equation}
Its inverse is
\begin{equation}
\textbf{W}^{-1} \equiv \textbf{B}^{\prime}\textbf{V}^{-1}\textbf{B},
\end{equation}
where
\begin{equation}
\textbf{B}\equiv\textbf{M}^{-1} = \left(\begin{array}{cccccc} I&0&0&\cdots&0&0\\ -A&I&0&\cdots&0&0\\ 0&-A&I&\cdots&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\0&0&0&\cdots&I&0\\0&0&0&\cdots&-A&I\end{array}\right).
\end{equation}
As a final remark, notice that in the single-variable case, $V_{1}$ is simply given by
\begin{equation}
V_{1} = \frac{V}{1-A^{2}},
\end{equation}
which leads to the formula OP quoted.
Update: To make it easier to work in a rearranged basis, one can explicitly carry out the matrix multiplication in $\textbf{W}^{-1} \equiv \textbf{B}^{\prime}\textbf{V}^{-1}\textbf{B}$. The result is
\begin{equation}
\textbf{W}^{-1} = \left(\begin{array}{cccccc} V_{1}^{-1}+P&-Q^{\prime}&0&\cdots&0&0\\ -Q&V^{-1}+P&-Q^{\prime}&\cdots&0&0\\ 0&-Q&V^{-1}+P&\cdots&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\0&0&0&\cdots&V^{-1}+P&-Q^{\prime}\\0&0&0&\cdots&-Q&V^{-1}\end{array}\right),
\end{equation}
where $P\equiv A^{\prime}V^{-1}A$ and $Q\equiv V^{-1}A$. One can now rearrange the rows and columns of this matrix to transform to a basis that is differently ordered.
