Menu
Home
People
Places
Arts
History
Plants & Animals
Science
Life & Culture
Technology
Reference.org
Block matrix pseudoinverse
open-in-new
<h2 id="derivation">Derivation</h2> <p>Consider a column-wise partitioned matrix: </p> [ A B ] , A ∈ R m × n , B ∈ R m × p , m ≥ n + p . {\displaystyle {\begin{bmatrix}\mathbf {A} &\mathbf {B} \end{bmatrix}},\quad \mathbf {A} \in \mathbb {R} ^{m\times n},\quad \mathbf {B} \in \mathbb {R} ^{m\times p},\quad m\geq n+p.} <p>If the above matrix is full column rank, the <a href="/facts/Moore%25E2%2580%2593Penrose_inverse/CU73tntq">Moore–Penrose inverse</a> matrices of it and its transpose are </p> [ A B ] + = ( [ A B ] T [ A B ] ) − 1 [ A B ] T , [ A T B T ] + = [ A B ] ( [ A B ] T [ A B ] ) − 1 . {\displaystyle {\begin{aligned}{\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{aligned}}} <p>This computation of the pseudoinverse requires (<i>n</i> + <i>p</i>)-square matrix inversion and does not take advantage of the block form. </p><p>To reduce computational costs to <i>n</i>- and <i>p</i>-square matrix inversions and to introduce parallelism, treating the blocks separately, one derives <a class="footnote-ref" id="fnref:1" href="#fn:1"><sup>1</sup></a> </p> [ A B ] + = [ P B ⊥ A ( A T P B ⊥ A ) − 1 P A ⊥ B ( B T P A ⊥ B ) − 1 ] = [ ( P B ⊥ A ) + ( P A ⊥ B ) + ] , [ A T B T ] + = [ P B ⊥ A ( A T P B ⊥ A ) − 1 , P A ⊥ B ( B T P A ⊥ B ) − 1 ] = [ ( A T P B ⊥ ) + ( B T P A ⊥ ) + ] , {\displaystyle {\begin{aligned}{\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{aligned}}} <p>where <a href="/facts/Orthogonal_projection/sElXGkxD">orthogonal projection</a> matrices are defined by </p> P A ⊥ = I − A ( A T A ) − 1 A T , P B ⊥ = I − B ( B T B ) − 1 B T . {\displaystyle {\begin{aligned}\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{aligned}}} <p>The above formulas are not necessarily valid if [ A B ] {\displaystyle {\begin{bmatrix}\mathbf {A} &\mathbf {B} \end{bmatrix}}} does not have full rank – for example, if A ≠ 0 {\displaystyle \mathbf {A} \neq 0} , then </p> [ A A ] + = 1 2 [ A + A + ] ≠ [ ( P A ⊥ A ) + ( P A ⊥ A ) + ] = 0 {\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} <h2 id="application-to-least-squares-problems">Application to least squares problems</h2> <p>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. </p> <h3>Column-wise partitioning in over-determined least squares</h3> <p>Suppose a solution x = [ x 1 x 2 ] {\displaystyle \mathbf {x} ={\begin{bmatrix}\mathbf {x} _{1}\\\mathbf {x} _{2}\\\end{bmatrix}}} solves an over-determined system: </p> [ A , B ] [ x 1 x 2 ] = d , d ∈ R m × 1 . {\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 \mathbb {R} ^{m\times 1}.} <p>Using the block matrix pseudoinverse, we have </p> x = [ A , B ] + d = [ ( P B ⊥ A ) + ( P A ⊥ B ) + ] d . {\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} .} <p>Therefore, we have a decomposed solution: </p> x 1 = ( P B ⊥ A ) + d , x 2 = ( P A ⊥ B ) + d . {\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} .} <h3>Row-wise partitioning in under-determined least squares</h3> <p>Suppose a solution x {\displaystyle \mathbf {x} } solves an under-determined system: </p> [ A T B T ] x = [ e f ] , e ∈ R n × 1 , f ∈ R p × 1 . {\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 \mathbb {R} ^{n\times 1},\quad \mathbf {f} \in \mathbb {R} ^{p\times 1}.} <p>The minimum-norm solution is given by </p> x = [ A T B T ] + [ e f ] . {\displaystyle \mathbf {x} ={\begin{bmatrix}\mathbf {A} ^{\textsf {T}}\\\mathbf {B} ^{\textsf {T}}\end{bmatrix}}^{+}\,{\begin{bmatrix}\mathbf {e} \\\mathbf {f} \end{bmatrix}}.} <p>Using the block matrix pseudoinverse, we have </p> x = [ ( A T P B ⊥ ) + ( B T P A ⊥ ) + ] [ e f ] = ( A T P B ⊥ ) + e + ( B T P A ⊥ ) + f . {\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} .} <h2 id="comments-on-matrix-inversion">Comments on matrix inversion</h2> <p>Instead of ( [ A B ] T [ A B ] ) − 1 {\displaystyle \mathbf {\left({\begin{bmatrix}\mathbf {A} &\mathbf {B} \end{bmatrix}}^{\textsf {T}}{\begin{bmatrix}\mathbf {A} &\mathbf {B} \end{bmatrix}}\right)} ^{-1}} , we need to calculate directly or indirectly[<i>original research?</i>] </p> ( A T A ) − 1 , ( B T B ) − 1 , ( A T P B ⊥ A ) − 1 , ( B T P A ⊥ B ) − 1 . {\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}.} <p>In a dense and small system, we can use <a href="/facts/Singular_value_decomposition/8t3v6wk3">singular value decomposition</a>, <a href="/facts/QR_decomposition/NswfKLGY">QR decomposition</a>, or <a href="/facts/Cholesky_decomposition/Q3IjkuTs">Cholesky decomposition</a> to replace the matrix inversions with numerical routines. In a large system, we may employ <a href="/facts/Iterative_methods/uKMAJhEh">iterative methods</a> such as Krylov subspace methods. </p><p>Considering <a href="/facts/Parallel_algorithms/nUPzEF4C">parallel algorithms</a>, we can compute ( A T A ) − 1 {\displaystyle \left(\mathbf {A} ^{\textsf {T}}\mathbf {A} \right)^{-1}} and ( B T B ) − 1 {\displaystyle \left(\mathbf {B} ^{\textsf {T}}\mathbf {B} \right)^{-1}} in parallel. Then, we finish to compute ( A T P B ⊥ A ) − 1 {\displaystyle \left(\mathbf {A} ^{\textsf {T}}\mathbf {P} _{B}^{\perp }\mathbf {A} \right)^{-1}} and ( B T P A ⊥ B ) − 1 {\displaystyle \left(\mathbf {B} ^{\textsf {T}}\mathbf {P} _{A}^{\perp }\mathbf {B} \right)^{-1}} also in parallel. </p> <h2 id="see-also">See also</h2> <ul><li><a href="/facts/Invertible_matrix/fPqXk3V8">Invertible matrix § Blockwise inversion</a></li></ul> <h2 id="external-links">External links</h2> <ul><li><a href="http://www.ee.ic.ac.uk/hp/staff/dmb/matrix/intro.html">The Matrix Reference Manual</a> by <a href="http://www.ee.ic.ac.uk/hp/staff/dmb/dmb.html">Mike Brookes</a></li> <li><a href="https://web.archive.org/web/20060414125709/http://www.csit.fsu.edu/~burkardt/papers/linear_glossary.html">Linear Algebra Glossary</a> by <a href="http://www.csit.fsu.edu/~burkardt/">John Burkardt</a></li> <li><a href="http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3274/pdf/imm3274.pdf">The Matrix Cookbook</a> by <a href="http://www2.imm.dtu.dk/pubdb/views/publication_details.php?id=3274/">Kaare Brandt Petersen</a></li> <li><a href="http://see.stanford.edu/materials/lsoeldsee263/08-min-norm.pdf">Lecture 8: Least-norm solutions of undetermined equations</a> by <a href="https://stanford.edu/~boyd/">Stephen P. Boyd</a></li></ul> <h2 id="references">References</h2> <ol> <li id="fn:1"><p>J.K. Baksalary and O.M. Baksalary (2007). "Particular formulae for the Moore–Penrose inverse of a columnwise partitioned matrix". Linear Algebra Appl. 421: 16–23. doi:10.1016/j.laa.2006.03.031. <a href="/wiki/Jerzy_Baksalary" target="_blank">/wiki/Jerzy_Baksalary</a> <a href="#fnref:1" class="footnote-back-ref">↩</a></p></li> </ol>