Tak (function)

2013-04-01T00:58:39Z

86.69.126.183:

{{refimprove|date=December 2010}}
In mathematics, '''block matrix pseudoinverse''' is a formula of [[pseudoinverse]] of a [[partitioned matrix]]. This is useful for decomposing or approximating many algorithms updating parameters in [[signal processing]], which are based on [[least squares]] method.

== Derivation ==
Consider a column-wise partitioned matrix:
:<math> [\mathbf A, \mathbf B], \qquad \mathbf A \in \reals^{n\times m}, \qquad \mathbf B \in \reals^{p\times m}, \qquad m \geq n+p.</math>

If the above matrix is full rank, the [[pseudoinverse]] matrices of it and its transpose are as follows.
:<math>
\begin{bmatrix}
\mathbf A, & \mathbf B
\end{bmatrix}
^{+} = ([\mathbf A, \mathbf B]^T [\mathbf A, \mathbf B])^{-1} [\mathbf A, \mathbf B]^T,
</math>
:<math>
\begin{bmatrix}
\mathbf A^T \\ \mathbf B^T
\end{bmatrix}
^{+} = [\mathbf A, \mathbf B] ([\mathbf A, \mathbf B]^T [\mathbf A, \mathbf B])^{-1}.
</math>
The pseudoinverse requires (''n'' + ''p'')-square matrix inversion.

To reduce complexity and introduce parallelism, we derive the following decomposed formula. From a block matrix inverse<math> \mathbf ([\mathbf A, \mathbf B]^T [\mathbf A, \mathbf B])^{-1}</math>, we can have{{citation needed|date=December 2010}}{{Original research|date=December 2010}}
:<math>
\begin{bmatrix}
\mathbf A, & \mathbf B
\end{bmatrix}
^{+} = \left[\mathbf P_B^\perp \mathbf A( \mathbf A^T \mathbf P_B^\perp \mathbf A)^{-1}, \quad \mathbf P_A^\perp \mathbf B(\mathbf B^T \mathbf P_A^\perp \mathbf B)^{-1}\right]^T,
</math>
:<math>
\begin{bmatrix}
\mathbf A^T \\ \mathbf B^T
\end{bmatrix}
^{+} = \left[\mathbf P_B^\perp \mathbf A( \mathbf A^T \mathbf P_B^\perp \mathbf A)^{-1}, \quad \mathbf P_A^\perp \mathbf B(\mathbf B^T \mathbf P_A^\perp \mathbf B)^{-1}\right],
</math>
where [[orthogonal projection]] matrices are defined by
::<math>
\begin{align}
\mathbf P_A^\perp & = \mathbf I - \mathbf A (\mathbf A^T \mathbf A)^{-1} \mathbf A^T, \\ \mathbf P_B^\perp & = \mathbf I - \mathbf B (\mathbf B^T \mathbf B)^{-1} \mathbf B^T.
\end{align}
</math>

Interestingly, from the [[idempotence]] of projection matrix, we can verify that the pseudoinverse of block matrix consists of pseudoinverse of projected matrices:
:<math>
\begin{bmatrix}
\mathbf A, & \mathbf B
\end{bmatrix}
^{+}
=
\begin{bmatrix}
(\mathbf P_B^{\perp}\mathbf A)^{+}
\\
(\mathbf P_A^{\perp}\mathbf B)^{+}
\end{bmatrix}, </math>
:<math>
\begin{bmatrix}
\mathbf A^T \\ \mathbf B^T
\end{bmatrix}
^{+}
= [(\mathbf A^T \mathbf P_B^{\perp})^{+},
\quad (\mathbf B^T \mathbf P_A^{\perp})^{+} ]. </math>

Thus, we decomposed the block matrix pseudoinverse into two submatrix pseudoinverses, which cost ''n''- and ''p''-square matrix inversions, respectively.

Note that the above formulae are not necessarily valid if <math>[\mathbf A, \mathbf B]</math> does not have full rank – for example, if <math>\mathbf A \neq 0</math>, then
:<math>
\begin{bmatrix}
\mathbf A, & \mathbf A
\end{bmatrix}
^{+}
= \frac{1}{2}
\begin{bmatrix}
\mathbf A^{+} \\ \mathbf A^{+}
\end{bmatrix}
\neq
\begin{bmatrix}
(\mathbf P_A^{\perp}\mathbf A)^{+}
\\
(\mathbf P_A^{\perp}\mathbf A)^{+}
\end{bmatrix}
= 0
</math>

== Application to least squares problems ==

Given the same matrices as above, we consider the following least squares problems, which
appear as multiple objective optimizations or constrained problems in signal processing.
Eventually, we can implement a parallel algorithm for least squares based on the following results.

=== Column-wise partitioning in over-determined least squares ===

Suppose a solution
<math> \mathbf x = \begin{bmatrix}
\mathbf x_1 \\
\mathbf x_2 \\
\end{bmatrix}
</math> solves an over-determined system:
:<math>
\begin{bmatrix}
\mathbf A, & \mathbf B
\end{bmatrix}
\begin{bmatrix}
\mathbf x_1 \\
\mathbf x_2 \\
\end{bmatrix}
=
\mathbf d
,
\qquad \mathbf d \in \reals^{m\times 1}.</math>

Using the block matrix pseudoinverse, we have
:<math>
\mathbf x
=
\begin{bmatrix}
\mathbf A, & \mathbf B
\end{bmatrix}
^{+}\,
\mathbf d
=
\begin{bmatrix}
(\mathbf P_B^{\perp} \mathbf A)^{+}\\
(\mathbf P_A^{\perp} \mathbf B)^{+}
\end{bmatrix}
\mathbf d
.
</math>
Therefore, we have a decomposed solution:
:<math>
\mathbf x_1
=
(\mathbf P_B^{\perp} \mathbf A)^{+}\,
\mathbf d
,
\qquad
\mathbf x_2
=
(\mathbf P_A^{\perp} \mathbf B)^{+}
\,
\mathbf d
.
</math>

=== Row-wise partitioning in under-determined least squares ===

Suppose a solution <math> \mathbf x </math> solves an under-determined system:
:<math>
\begin{bmatrix}
\mathbf A^T \\ \mathbf B^T
\end{bmatrix}
\mathbf x
=
\begin{bmatrix}
\mathbf e \\ \mathbf f
\end{bmatrix},
\qquad \mathbf e \in \reals^{n\times 1},
\qquad \mathbf f \in \reals^{p\times 1}.</math>

The minimum-norm solution is given by
:<math>
\mathbf x
=
\begin{bmatrix}
\mathbf A^T \\ \mathbf B^T
\end{bmatrix}
^{+}\,
\begin{bmatrix}
\mathbf e \\ \mathbf f
\end{bmatrix}.
</math>

Using the block matrix pseudoinverse, we have
:<math>
\mathbf x
=
[(\mathbf A^T\mathbf P_B^{\perp})^{+},
\quad (\mathbf B^T\mathbf P_A^{\perp})^{+} ]
\begin{bmatrix}
\mathbf e \\ \mathbf f
\end{bmatrix}
=
(\mathbf A^T\mathbf P_B^{\perp})^{+}\,\mathbf e
+
(\mathbf B^T\mathbf P_A^{\perp} )^{+}\,\mathbf f
.
</math>

== Comments on matrix inversion ==

Instead of <math> \mathbf ([\mathbf A, \mathbf B]^T [\mathbf A, \mathbf B])^{-1}</math>,
we need to calculate directly or indirectly{{citation needed|date=December 2010}}{{original research?|date=December 2010}}

:<math>
\quad (\mathbf A^T \mathbf A)^{-1},
\quad (\mathbf B^T \mathbf B)^{-1},
\quad (\mathbf A^T \mathbf P_B^{\perp} \mathbf A)^{-1},
\quad (\mathbf B^T \mathbf P_A^{\perp} \mathbf B)^{-1}
.
</math>

In a dense and small system, we can use [[singular value decomposition]], [[QR decomposition]], or [[Cholesky decomposition]] to replace the matrix inversions with numerical routines. In a large system, we may employ [[iterative methods]] such as Krylov subspace methods.

Considering [[parallel algorithms]], we can compute <math>(\mathbf A^T \mathbf A)^{-1}</math> and
<math>(\mathbf B^T \mathbf B)^{-1}</math> in parallel. Then, we finish to compute <math>(\mathbf A^T \mathbf P_B^{\perp} \mathbf A)^{-1}</math> and <math>(\mathbf B^T \mathbf P_A^{\perp} \mathbf B)^{-1}</math> also in parallel.

== Block matrix inversion ==
Let a block matrix be
:<math>\begin{bmatrix}
A & B \\
C & D
\end{bmatrix}
.</math>
We can get an inverse formula by combining the previous results in.<ref>{{Cite web |url=http://ccrma.stanford.edu/~jos/lattice/Block_matrix_decompositions.html|title=Block matrix decompositions |title=Block matrix decompositions |last=Smith |first=Julius O |archiveurl=http://web.archive.org/web/20080915181439/http://ccrma.stanford.edu/~jos/lattice/Block_matrix_decompositions.html |archivedate=2008-09-15}}</ref>
:<math>\begin{bmatrix}
A & B \\
C & D
\end{bmatrix}^{-1}
=
\begin{bmatrix}
(A - BD^{-1}C)^{-1} & -A^{-1}B(D - CA^{-1}B)^{-1} \\
-D^{-1}C(A - BD^{-1}C)^{-1} & (D - CA^{-1}B)^{-1}
\end{bmatrix}
=
\begin{bmatrix}
S^{-1}_D & -A^{-1}BS^{-1}_A \\
-D^{-1}CS^{-1}_D & S^{-1}_A
\end{bmatrix}
,</math>
where <math>S_A</math> and <math>S_D</math>, respectively, [[Schur complements]] of <math>A</math>
and <math>D</math>, are defined by <math>S_A = D - C A^{-1}B</math>, and <math>S_D =
A - BD^{-1}C</math>. This relation is derived by using Block Triangular
Decomposition. It is called ''simple block matrix inversion.''<ref>[http://ieeexplore.ieee.org/xpls/abs_all.jsp?isnumber=30419&arnumber=1399280&count=249&index=181 S. Jo, S. W. Kim and T. J. Park, "Equally constrained affine projection algorithm," ''in Conference Record of the Thirty-Eighth Asilomar Conference on Signals, Systems and Computers,'' vol. 1, pp. 955–959, Nov. 7–10, 2004.]</ref>

Now we can obtain the inverse of the symmetric block matrix:
:<math>
\begin{bmatrix}
\mathbf A^T \mathbf A & \mathbf A^T \mathbf B \\
\mathbf B^T \mathbf A & \mathbf B^T \mathbf B
\end{bmatrix}^{-1}
=
\begin{bmatrix}
(\mathbf A^T \mathbf A-\mathbf A^T \mathbf B(\mathbf B^T \mathbf B)^{-1}\mathbf B^T \mathbf A)^{-1}
& -(\mathbf A^T \mathbf A)^{-1}\mathbf A^T \mathbf B(\mathbf B^T \mathbf B-\mathbf B^T \mathbf A(\mathbf A^T \mathbf A)^{-1}\mathbf A^T \mathbf B)^{-1}
\\
-(\mathbf B^T \mathbf B)^{-1}\mathbf B^T \mathbf A(\mathbf A^T \mathbf A-\mathbf A^T \mathbf B(\mathbf B^T \mathbf B)^{-1}\mathbf B^T \mathbf A)^{-1}
& (\mathbf B^T \mathbf B-\mathbf B^T \mathbf A(\mathbf A^T \mathbf A)^{-1}\mathbf A^T \mathbf B)^{-1}
\end{bmatrix}
</math>
:::<math>
=
\begin{bmatrix}
(\mathbf A^T \mathbf P_B^\perp \mathbf A)^{-1}
& -(\mathbf A^T \mathbf A)^{-1}\mathbf A^T \mathbf B(\mathbf B^T \mathbf P_A^\perp \mathbf B)^{-1}
\\
-(\mathbf B^T \mathbf B)^{-1}\mathbf B^T \mathbf A(\mathbf A^T \mathbf P_B^\perp \mathbf A)^{-1}
& (\mathbf B^T \mathbf P_A^{\perp} \mathbf B)^{-1}
\end{bmatrix}
</math>
Since the block matrix is symmetric, we also have
:<math>
\begin{bmatrix}
\mathbf A^T \mathbf A & \mathbf A^T \mathbf B \\
\mathbf B^T \mathbf A & \mathbf B^T \mathbf B
\end{bmatrix}^{-1}
=
\begin{bmatrix}
(\mathbf A^T \mathbf P_B^{\perp} \mathbf A)^{-1}
&
-(\mathbf A^T \mathbf P_B^{\perp} \mathbf A)^{-1}
\mathbf A^T \mathbf B(\mathbf B^T \mathbf B)^{-1}
\\
-(\mathbf B^T \mathbf P_A^{\perp} \mathbf B)^{-1}
\mathbf B^T \mathbf A (\mathbf A^T \mathbf A)^{-1}
& (\mathbf B^T \mathbf P_A^{\perp} \mathbf B)^{-1}
\end{bmatrix}.
</math>

Then, we can see how the Schur complements are connected to the projection matrices of the symmetric, partitioned matrix.

== See also ==
*[[Invertible matrix#Blockwise inversion]]

==References ==
{{Reflist}}

== External links ==
* [http://www.ee.ic.ac.uk/hp/staff/dmb/matrix/intro.html The Matrix Reference Manual] by [http://www.ee.ic.ac.uk/hp/staff/dmb/dmb.html Mike Brookes]
* [http://www.csit.fsu.edu/~burkardt/papers/linear_glossary.html Linear Algebra Glossary] by [http://www.csit.fsu.edu/~burkardt/ John Burkardt]
* [http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3274/pdf/imm3274.pdf The Matrix Cookbook] by [http://www2.imm.dtu.dk/pubdb/views/publication_details.php?id=3274/ Kaare Brandt Petersen]
* [http://see.stanford.edu/materials/lsoeldsee263/08-min-norm.pdf Lecture 8: Least-norm solutions of undetermined equations] by [http://www.stanford.edu/~boyd/ Stephen P. Boyd]

{{Numerical linear algebra}}

{{DEFAULTSORT:Block Matrix Pseudoinverse}}
[[Category:Numerical linear algebra]]
[[Category:Matrix theory]]

formulasearchengine - User contributions [en]

Tak (function)