Block matrix pseudoinverse

In mathematics, a block matrix pseudoinverse is a formula for the pseudoinverse of a partitioned matrix. This is useful for decomposing or approximating many algorithms updating parameters in signal processing, which are based on the least squares method.

Derivation

Consider a column-wise partitioned matrix:

[math]\displaystyle{ \begin{bmatrix}\mathbf A & \mathbf B\end{bmatrix},\quad \mathbf A \in \reals^{m \times n},\quad \mathbf B \in \reals^{m \times p},\quad m \geq n + p. }[/math]

If the above matrix is full column rank, the Moore–Penrose inverse matrices of it and its transpose are

[math]\displaystyle{ \begin{align} \begin{bmatrix}\mathbf A & \mathbf B\end{bmatrix}^+ &= \left( \begin{bmatrix}\mathbf A & \mathbf B\end{bmatrix}^\textsf{T} \begin{bmatrix}\mathbf A & \mathbf B\end{bmatrix} \right)^{-1} \begin{bmatrix}\mathbf A & \mathbf B\end{bmatrix}^\textsf{T}, \\ \begin{bmatrix} \mathbf A^\textsf{T} \\ \mathbf B^\textsf{T} \end{bmatrix}^+ &= \begin{bmatrix}\mathbf A & \mathbf B\end{bmatrix} \left( \begin{bmatrix}\mathbf A & \mathbf B\end{bmatrix}^\textsf{T} \begin{bmatrix}\mathbf A & \mathbf B\end{bmatrix} \right)^{-1}. \end{align} }[/math]

This computation of the pseudoinverse requires (n + p)-square matrix inversion and does not take advantage of the block form.

To reduce computational costs to n- and p-square matrix inversions and to introduce parallelism, treating the blocks separately, one derives

[math]\displaystyle{ \begin{align} \begin{bmatrix}\mathbf A & \mathbf B\end{bmatrix}^+ &= \begin{bmatrix} \mathbf P_B^\perp \mathbf A\left(\mathbf A^\textsf{T} \mathbf P_B^\perp \mathbf A\right)^{-1} \\ \mathbf P_A^\perp \mathbf B\left(\mathbf B^\textsf{T} \mathbf P_A^\perp \mathbf B\right)^{-1} \end{bmatrix} = \begin{bmatrix} \left(\mathbf P_B^\perp\mathbf A\right)^+ \\ \left(\mathbf P_A^\perp\mathbf B\right)^+ \end{bmatrix}, \\ \begin{bmatrix} \mathbf A^\textsf{T} \\ \mathbf B^\textsf{T} \end{bmatrix}^+ &= \begin{bmatrix} \mathbf P_B^\perp \mathbf A\left(\mathbf A^\textsf{T} \mathbf P_B^\perp \mathbf A\right)^{-1},\quad \mathbf P_A^\perp \mathbf B\left(\mathbf B^\textsf{T} \mathbf P_A^\perp \mathbf B\right)^{-1} \end{bmatrix} = \begin{bmatrix} \left(\mathbf A^\textsf{T} \mathbf P_B^\perp\right)^+ & \left(\mathbf B^\textsf{T} \mathbf P_A^\perp\right)^+ \end{bmatrix}, \end{align} }[/math]

where orthogonal projection matrices are defined by

[math]\displaystyle{ \begin{align} \mathbf P_A^\perp &= \mathbf I - \mathbf A \left(\mathbf A^\textsf{T} \mathbf A\right)^{-1} \mathbf A^\textsf{T}, \\ \mathbf P_B^\perp &= \mathbf I - \mathbf B \left(\mathbf B^\textsf{T} \mathbf B\right)^{-1} \mathbf B^\textsf{T}. \end{align} }[/math]

The above formulas are not necessarily valid if [math]\displaystyle{ \begin{bmatrix}\mathbf A & \mathbf B\end{bmatrix} }[/math] does not have full rank – for example, if [math]\displaystyle{ \mathbf A \neq 0 }[/math], then

[math]\displaystyle{ \begin{bmatrix}\mathbf A & \mathbf A\end{bmatrix}^+ = \frac{1}{2}\begin{bmatrix} \mathbf A^+ \\ \mathbf A^+ \end{bmatrix} \neq \begin{bmatrix} \left(\mathbf P_A^\perp\mathbf A\right)^+ \\ \left(\mathbf P_A^\perp\mathbf A\right)^+ \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]\displaystyle{ \mathbf x = \begin{bmatrix} \mathbf x_1 \\ \mathbf x_2 \\ \end{bmatrix} }[/math] solves an over-determined system:

[math]\displaystyle{ \begin{bmatrix} \mathbf A, & \mathbf B \end{bmatrix} \begin{bmatrix} \mathbf x_1 \\ \mathbf x_2 \\ \end{bmatrix} = \mathbf d,\quad \mathbf d \in \reals^{m \times 1}. }[/math]

Using the block matrix pseudoinverse, we have

[math]\displaystyle{ \mathbf x = \begin{bmatrix} \mathbf A, & \mathbf B \end{bmatrix}^+\,\mathbf d = \begin{bmatrix} \left(\mathbf P_B^\perp \mathbf A\right)^+ \\ \left(\mathbf P_A^\perp \mathbf B\right)^+ \end{bmatrix}\mathbf d. }[/math]

Therefore, we have a decomposed solution:

[math]\displaystyle{ \mathbf x_1 = \left(\mathbf P_B^\perp \mathbf A\right)^+\,\mathbf d,\quad \mathbf x_2 = \left(\mathbf P_A^\perp \mathbf B\right)^+\,\mathbf d. }[/math]

Row-wise partitioning in under-determined least squares

Suppose a solution [math]\displaystyle{ \mathbf x }[/math] solves an under-determined system:

[math]\displaystyle{ \begin{bmatrix} \mathbf A^\textsf{T} \\ \mathbf B^\textsf{T} \end{bmatrix}\mathbf x = \begin{bmatrix} \mathbf e \\ \mathbf f \end{bmatrix},\quad \mathbf e \in \reals^{n \times 1},\quad \mathbf f \in \reals^{p \times 1}. }[/math]

The minimum-norm solution is given by

[math]\displaystyle{ \mathbf x = \begin{bmatrix} \mathbf A^\textsf{T} \\ \mathbf B^\textsf{T} \end{bmatrix}^+\, \begin{bmatrix} \mathbf e \\ \mathbf f \end{bmatrix}. }[/math]

Using the block matrix pseudoinverse, we have

[math]\displaystyle{ \mathbf x = \begin{bmatrix} \left(\mathbf A^\textsf{T}\mathbf P_B^\perp\right)^+ & \left(\mathbf B^\textsf{T}\mathbf P_A^\perp\right)^+ \end{bmatrix} \begin{bmatrix} \mathbf e \\ \mathbf f \end{bmatrix} = \left(\mathbf A^\textsf{T}\mathbf P_B^\perp\right)^+\,\mathbf e + \left(\mathbf B^\textsf{T}\mathbf P_A^\perp\right)^+\,\mathbf f. }[/math]

Comments on matrix inversion

Instead of [math]\displaystyle{ \mathbf \left(\begin{bmatrix}\mathbf A & \mathbf B\end{bmatrix}^\textsf{T} \begin{bmatrix}\mathbf A & \mathbf B\end{bmatrix}\right)^{-1} }[/math], we need to calculate directly or indirectly

[math]\displaystyle{ \left(\mathbf A^\textsf{T} \mathbf A\right)^{-1},\quad \left(\mathbf B^\textsf{T} \mathbf B\right)^{-1},\quad \left(\mathbf A^\textsf{T} \mathbf P_B^\perp \mathbf A\right)^{-1},\quad \left(\mathbf B^\textsf{T} \mathbf P_A^\perp \mathbf B\right)^{-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]\displaystyle{ \left(\mathbf A^\textsf{T} \mathbf A\right)^{-1} }[/math] and [math]\displaystyle{ \left(\mathbf B^\textsf{T} \mathbf B\right)^{-1} }[/math] in parallel. Then, we finish to compute [math]\displaystyle{ \left(\mathbf A^\textsf{T} \mathbf P_B^\perp \mathbf A\right)^{-1} }[/math] and [math]\displaystyle{ \left(\mathbf B^\textsf{T} \mathbf P_A^\perp \mathbf B\right)^{-1} }[/math] also in parallel.

See also

• Invertible matrix § Blockwise inversion

External links

• The Matrix Reference Manual by Mike Brookes
• Linear Algebra Glossary by John Burkardt
• The Matrix Cookbook by Kaare Brandt Petersen
• Lecture 8: Least-norm solutions of undetermined equations by Stephen P. Boyd