2
$\begingroup$

Updated: the goal is to solve $(A^\top A+B^\top B + D)x=c$. Maybe it is not necessary to compute $(A^\top A+B^\top B + D)^{-1}$.

Denote $e=(1,1,\ldots,1)^\top\in\mathbb{R}^n$ and $$A=\begin{bmatrix} e & & & \\ & e & & \\ & & \ddots & \\ & & & e \end{bmatrix}$$ and $$B=\begin{bmatrix}\mathrm{diag}(e) & \mathrm{diag}(e) & \cdots & \mathrm{diag}(e)\end{bmatrix}$$ where $e$ appears $n$ times in $A$ and $n$ times in $B$ (i.e. $A,B\in\mathbb{R}^{n\times n^2}$).

For example, if $n=3$ we have:

$$A=\begin{bmatrix} 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 \end{bmatrix}$$

$$B=\begin{bmatrix} 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0\\ 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0\\ 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 \end{bmatrix}$$

Let $D$ be a $n^2\times n^2$ diagonal matrix with positive elements.

I'm looking for an efficient way to solve the linear equation $$(A^\top A+B^\top B + D)x=c$$ where $x\in\mathbb{R}^{n^2}$ is the variable.

Thank you in advance for any suggestions.

$\endgroup$
  • $\begingroup$ The matrix $A^TA+B^TB+D$ is block diagonal, the blocks having a form $\tilde{D}+ee^T$, where $\tilde{D}$ is diagonal (a block from $D$ plus an identity). To invert them, use the Woodbury formula. $\endgroup$ – Algebraic Pavel Jun 20 '14 at 23:33
  • $\begingroup$ @AlgebraicPavel: thanks for your comment, but I don't think that matrix is block diagonal. Indeed, we have: $\endgroup$ – Khue Jun 24 '14 at 11:13
  • $\begingroup$ $A^\top A=\begin{bmatrix} 1&1&1&0&0&0&0&0&0\\ 1&1&1&0&0&0&0&0&0\\ 1&1&1&0&0&0&0&0&0\\ 0&0&0&1&1&1&0&0&0\\ 0&0&0&1&1&1&0&0&0\\ 0&0&0&1&1&1&0&0&0\\ 0&0&0&0&0&0&1&1&1\\ 0&0&0&0&0&0&1&1&1\\ 0&0&0&0&0&0&1&1&1\\ \end{bmatrix}$ while $B^\top B = \begin{bmatrix} 1&0&0&1&0&0&1&0&0\\ 0&1&0&0&1&0&0&1&0\\ 0&0&1&0&0&1&0&0&1\\ 1&0&0&1&0&0&1&0&0\\ 0&1&0&0&1&0&0&1&0\\ 0&0&1&0&0&1&0&0&1\\ 1&0&0&1&0&0&1&0&0\\ 0&1&0&0&1&0&0&1&0\\ 0&0&1&0&0&1&0&0&1\\ \end{bmatrix}$ $\endgroup$ – Khue Jun 24 '14 at 11:21
  • $\begingroup$ True, I was looking on that from "other" side :) You can still use the Woodbury formula though. Let me think of it a bit. $\endgroup$ – Algebraic Pavel Jun 24 '14 at 11:34
  • $\begingroup$ @AlgebraicPavel: I refreshed the page to see your updated anwser but it disappeared. Did you find some mistakes? $\endgroup$ – Khue Jun 24 '14 at 12:30
0
$\begingroup$

Let $C=A^TA+B^TB+D$. You can write it as $$ C=D+EE^T, \quad E=[A^T,B^T]. $$ The inverse can be then using the Woodbury formula written as $$ \begin{split} C^{-1}&=D^{-1}-D^{-1}E(I+E^TD^{1}E)^{-1}E^TD^{-1}\\ &=D^{-1}-D^{-1}[A^T,B^T]\left\{I+\begin{bmatrix}A\\B\end{bmatrix}D^{-1}[A^T,B^T]\right\}^{-1}\begin{bmatrix}A\\B\end{bmatrix}D^{-1} \end{split}\tag{$*$} $$ Assume that $D=\mathrm{diag}(D_1,\ldots,D_n)$, where each diagonal block $D_i$ is $n\times n$. Since $A=\mathrm{diag}(e^T)$ and $B=\mathrm{diag}(I_n)$, the "small" matrix inside the curly brackets can be written as $$ I+\begin{bmatrix}A\\B\end{bmatrix}D^{-1}[A^T,B^T]= \left[\begin{array}{ccc|c} 1+e^TD_1^{-1}e & & & e^TD_1^{-1} \\ & \ddots & & \vdots \\ & & 1+e^TD_n^{-1}e & e^TD_n^{-1} \\ \hline D_1^{-1}e & & D_n^{-1}e & I+D_1^{-1}+\cdots+D_n^{-1} \end{array}\right]. $$ It can hence be written in the form $$\tag{✿} I+\begin{bmatrix}A\\B\end{bmatrix}D^{-1}[A^T,B^T]=\begin{bmatrix}\Phi&\Theta^T\\\Theta&\Sigma\end{bmatrix}, $$ where the matrices $\Phi=\mathrm{diag}(1+e^TD_i^{-1}e)_{i=1}^n$ and $\Sigma=I+D_1^{-1}+\cdots+D_n^{-1}$ are $n\times n$ diagonal matrices. The $n\times n$ matrix $\Theta$ contains in the $i$th column the entries of the diagonal of $D_i^{-1}$. Now you can use some block inversion formulas, e.g., both $\Phi$ and the Schur complement $\Sigma-\Theta\Phi^{-1}\Theta^T$ are invertible (the Schur complement is normally a dense matrix). You can plug the computed inverse back to ($*$) and compute the inverse of $C^{-1}$.

Note that you actually don't need the inverse explicitly. Assuming that the matrix (✿) has a Cholesky factorisation $LL^T$, all you need is to solve the system $$ LX=\begin{bmatrix}A\\B\end{bmatrix} $$ with multiple right-hand sides. Then $$ C^{-1}= D^{-1}-D^{-1}[A^T,B^T](LL^T)^{-1}\begin{bmatrix}A\\B\end{bmatrix}D^{-1} =D^{-1}-D^{-1}X^TXD^{-1}. $$

To solve a system with $C$, you don't need to compute $X$ but still you need to factorise the matrix (✿). You can use ($*$) to solve a system $Cx=b$ as follows:

  • Compute $y=D^{-1}b$,
  • Compute $z=E^Ty$,
  • Solve the system with (✿) by $Lu=z$ and $L^Tv=u$,
  • Compute $s=D^{-1}Ev$,
  • Set $x=y-s$.

Multiplying with $E$ and $E^T$ is easy as it involves only some summations and copies.

NOTE: I'm afraid that there might be no nicely looking solution except in cases where $D$ is a multiple of identity (or at least the matrices $D_i$ are multiples of identity). Otherwise, $\Theta$ is a general $n\times n$ matrix (with positive entries) so the inverse of the Schur complement $\Sigma-\Theta\Phi^{-1}\Theta^T$ cannot be written in some simple form. However, in the mentioned special cases the matrix $\Theta\Phi^{-1}\Theta^T$ would have rank one so one could use the Sherman-Morrison formula to invert the Schur complement.

NOTE 2: The matrix (✿) has a special form so this can be somewhat exploited in the Cholesky factorisation. The matrix can be factorised as $$ \begin{bmatrix}\Phi&\Theta^T\\\Theta&\Sigma\end{bmatrix} = \begin{bmatrix}I&0\\\Theta\Phi^{-1}&I\end{bmatrix} \begin{bmatrix}\Phi&0\\0&\Pi\end{bmatrix} \begin{bmatrix}I&\Phi^{-1}\Theta^T\\0&I\end{bmatrix}, $$ where $\Pi=\Sigma-\Theta\Phi^{-1}\Theta^T$ is the already mentioned Schur complement. So if $\tilde{L}\tilde{L}^T=\Pi$ is the Cholesky factorisation of $\Pi$, then $$ \begin{bmatrix}\Phi&\Theta^T\\\Theta&\Sigma\end{bmatrix} = \underbrace{\begin{bmatrix}\Phi^{1/2}&0\\\Theta\Phi^{-1/2}&\tilde{L}\end{bmatrix}}_{L} \underbrace{\begin{bmatrix}\Phi^{1/2}&\Phi^{-1/2}\Theta^T\\0&\tilde{L}^T\end{bmatrix}}_{L^T} $$ is the Cholesky factorisation of (✿). Of course, some rearrangements can be made to avoid computing the square roots.

$\endgroup$
  • $\begingroup$ Thanks Algebraic Pavel for your effort. I'm checking carefully your post. There's a small mistake in the application of the Woodbury formula, it should be $E^TD^{-1}E$ instead of $E^TDE$, but I hope this does not affect too much the other arguments. I'm still checking... $\endgroup$ – Khue Jun 25 '14 at 13:22
  • $\begingroup$ @Khue Thanks, it's fixed. $\endgroup$ – Algebraic Pavel Jun 25 '14 at 13:36
  • $\begingroup$ Sorry it took me a while to check your answer. Everything seems to be fine except a typo in the Cholesky factorization of (✿): you forgot the transpose sign in the last matrix (it should read $\Phi^{-1}\Theta^T$ instead of $\Phi^{-1}\Theta$). I will try to implement to see if this solution gives reasonable performance (as I see the most expensive steps are solving $Lu=z,L^Tv=u$ and the Cholesky factorization of $\Pi$). Thanks again for your help. $\endgroup$ – Khue Jun 27 '14 at 14:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.