Computation of a Moore-Penrose Pseudoinverse
============================================

In this example it is the goal to compute derivatives of the Moore-Penrose pseudoinverse.
We compare different possibilities:

    1) A naive approach where A^T A is explicitly computed (numerically unstable)
    2) A QR approach where at first a QR decomposition of A is formed and the inverse
       is computed by a forward and then back substitution of R. T
    3) A direct approach where an analytic formula for the derivatives of the
       Moore-Penrose Formula is derived. Then usage of the QR decomposition is
       used to make the computation of the analytical formula numerically stable.

More explicitly, we compute

.. math::

    A^\dagger = (A^T A)^{-1} A^T

where :math:`A \in \mathbb R^{M \times N}` with :math:`M \geq N` with possibly
bad condition number but full column rank.

.. literalinclude:: moore_penrose_pseudoinverse.py


Running the code yields the output shown below. One can see that computing 
the Moore-Penrose pseudoinverse with the naive approach yields a solution
(including first directional derivatives)
that only poorly fulfills the defining equations (c.f. e.g. http://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_pseudoinverse).
The other two methods yield a solution that results in a similarly small residual
of the defining equations.
Note that the result is a (D,P,M,N) array, i.e., the first block of numbers is the 
residual of the nominal solution and the second block of numbers is the residual
in the first derivatives::
    
    naive approach: A Apinv A - A = 0 
    [[[[  1.38013414e-06   1.38012688e-06]
       [  1.38012190e-06   1.38018107e-06]
       [  1.38013414e-06   1.38012688e-06]
       [  1.38012190e-06   1.38018107e-06]
       [  1.38013414e-06   1.38012688e-06]]]
    
    
     [[[ -3.14117999e-06  -3.14127715e-06]
       [  2.91085793e-05   2.91090665e-05]
       [ -3.14117999e-06  -3.14127715e-06]
       [  2.91085793e-05   2.91090665e-05]
       [ -3.14117999e-06  -3.14127715e-06]]]]
    naive approach: Apinv A Apinv - Apinv = 0 
    [[[[ 0.22144346 -0.33215987  0.22144346 -0.33215987  0.22144346]
       [-0.22144008  0.3321555  -0.22144008  0.3321555  -0.22144008]]]
    
    
     [[[-1.0878914   0.01930878 -1.0878914   0.01930878 -1.0878914 ]
       [ 1.08787702 -0.01930766  1.08787702 -0.01930766  1.08787702]]]]
    naive approach: (Apinv A)^T - Apinv A = 0 
    [[[[  0.00000000e+00   5.64176865e-08]
       [ -5.64176865e-08   0.00000000e+00]]]
    
    
     [[[  0.00000000e+00   6.45011636e+00]
       [ -6.45011636e+00   0.00000000e+00]]]]
    naive approach: (A Apinv)^T - A Apinv = 0 
    [[[[  0.00000000e+00   6.20026697e-10   0.00000000e+00   6.20026697e-10
          0.00000000e+00]
       [ -6.20026697e-10   0.00000000e+00  -6.20026697e-10   0.00000000e+00
         -6.20026697e-10]
       [  0.00000000e+00   6.20026697e-10   0.00000000e+00   6.20026697e-10
          0.00000000e+00]
       [ -6.20026697e-10   0.00000000e+00  -6.20026697e-10   0.00000000e+00
         -6.20026697e-10]
       [  0.00000000e+00   6.20026697e-10   0.00000000e+00   6.20026697e-10
          0.00000000e+00]]]
    
    
     [[[  0.00000000e+00   6.45163836e-06   0.00000000e+00   6.45163836e-06
          0.00000000e+00]
       [ -6.45163836e-06   0.00000000e+00  -6.45163836e-06   0.00000000e+00
         -6.45163836e-06]
       [  0.00000000e+00   6.45163836e-06   0.00000000e+00   6.45163836e-06
          0.00000000e+00]
       [ -6.45163836e-06   0.00000000e+00  -6.45163836e-06   0.00000000e+00
         -6.45163836e-06]
       [  0.00000000e+00   6.45163836e-06   0.00000000e+00   6.45163836e-06
          0.00000000e+00]]]]
    QR approach: A Apinv A - A = 0 
    [[[[ -3.18345350e-12  -3.18389759e-12]
       [  1.54691815e-11   1.54698476e-11]
       [ -3.18345350e-12  -3.18389759e-12]
       [  1.54691815e-11   1.54698476e-11]
       [ -3.18345350e-12  -3.18389759e-12]]]
    
    
     [[[  6.38786801e-11   6.38786801e-11]
       [ -4.68758365e-11  -4.68767247e-11]
       [  6.38786801e-11   6.38786801e-11]
       [ -4.68758365e-11  -4.68767247e-11]
       [  6.38786801e-11   6.38786801e-11]]]]
    QR approach: Apinv A Apinv - Apinv = 0 
    [[[[  1.94382301e-06  -3.84821760e-06   1.94382301e-06  -3.84821760e-06
          1.94382301e-06]
       [ -1.94371387e-06   3.84805026e-06  -1.94371387e-06   3.84805026e-06
         -1.94371387e-06]]]
    
    
     [[[ -6.20307401e-06   1.85714744e-05  -6.20307401e-06   1.85714744e-05
         -6.20307401e-06]
       [  6.20371429e-06  -1.85723184e-05   6.20371429e-06  -1.85723184e-05
          6.20371429e-06]]]]
    QR approach: (Apinv A)^T - Apinv A = 0 
    [[[[  0.00000000e+00   3.73022496e-06]
       [ -3.73022496e-06   0.00000000e+00]]]
    
    
     [[[  0.00000000e+00  -2.96084671e-05]
       [  2.96084671e-05   0.00000000e+00]]]]
    QR approach: (A Apinv)^T - A Apinv = 0 
    [[[[  0.00000000e+00   8.84625706e-13   0.00000000e+00   8.84625706e-13
          0.00000000e+00]
       [ -8.84625706e-13   0.00000000e+00  -8.84625706e-13   0.00000000e+00
         -8.84625706e-13]
       [  0.00000000e+00   8.84625706e-13   0.00000000e+00   8.84625706e-13
          0.00000000e+00]
       [ -8.84625706e-13   0.00000000e+00  -8.84625706e-13   0.00000000e+00
         -8.84625706e-13]
       [  0.00000000e+00   8.84625706e-13   0.00000000e+00   8.84625706e-13
          0.00000000e+00]]]
    
    
     [[[  0.00000000e+00   9.56390522e-12   0.00000000e+00   9.56390522e-12
          0.00000000e+00]
       [ -9.56390522e-12   0.00000000e+00  -9.56390522e-12   0.00000000e+00
         -9.56390522e-12]
       [  0.00000000e+00   9.56390522e-12   0.00000000e+00   9.56390522e-12
          0.00000000e+00]
       [ -9.56390522e-12   0.00000000e+00  -9.56390522e-12   0.00000000e+00
         -9.56390522e-12]
       [  0.00000000e+00   9.56390522e-12   0.00000000e+00   9.56390522e-12
          0.00000000e+00]]]]
    analytical approach: A Apinv A - A = 0 
    [[[[ -3.18345350e-12  -3.18389759e-12]
       [  1.54691815e-11   1.54698476e-11]
       [ -3.18345350e-12  -3.18389759e-12]
       [  1.54691815e-11   1.54698476e-11]
       [ -3.18345350e-12  -3.18389759e-12]]]
    
    
     [[[  2.81423151e-09   2.70797518e-09]
       [ -4.25869606e-09  -4.09931600e-09]
       [  2.81423151e-09   2.70797518e-09]
       [ -4.25869606e-09  -4.09931600e-09]
       [  2.81423151e-09   2.70797518e-09]]]]
    analytical approach: Apinv A Apinv - Apinv = 0 
    [[[[  1.94382301e-06  -3.84821760e-06   1.94382301e-06  -3.84821760e-06
          1.94382301e-06]
       [ -1.94371387e-06   3.84805026e-06  -1.94371387e-06   3.84805026e-06
         -1.94371387e-06]]]
    
    
     [[[  8.85959818e-01  -1.32856906e+00   8.85959818e-01  -1.32856906e+00
          8.85959818e-01]
       [ -8.85947414e-01   1.32855046e+00  -8.85947414e-01   1.32855046e+00
         -8.85947414e-01]]]]
    analytical approach: (Apinv A)^T - Apinv A = 0 
    [[[[  0.00000000e+00   3.73022496e-06]
       [ -3.73022496e-06   0.00000000e+00]]]
    
    
     [[[  0.00000000e+00  -1.39551660e-03]
       [  1.39551660e-03   0.00000000e+00]]]]
    analytical approach: (A Apinv)^T - A Apinv = 0 
    [[[[  0.00000000e+00   8.84625706e-13   0.00000000e+00   8.84625706e-13
          0.00000000e+00]
       [ -8.84625706e-13   0.00000000e+00  -8.84625706e-13   0.00000000e+00
         -8.84625706e-13]
       [  0.00000000e+00   8.84625706e-13   0.00000000e+00   8.84625706e-13
          0.00000000e+00]
       [ -8.84625706e-13   0.00000000e+00  -8.84625706e-13   0.00000000e+00
         -8.84625706e-13]
       [  0.00000000e+00   8.84625706e-13   0.00000000e+00   8.84625706e-13
          0.00000000e+00]]]
    
    
     [[[  0.00000000e+00  -1.37629996e-09   0.00000000e+00  -1.37629996e-09
          0.00000000e+00]
       [  1.37629996e-09   0.00000000e+00   1.37629996e-09   0.00000000e+00
          1.37629996e-09]
       [  0.00000000e+00  -1.37629996e-09   0.00000000e+00  -1.37629996e-09
          0.00000000e+00]
       [  1.37629996e-09   0.00000000e+00   1.37629996e-09   0.00000000e+00
          1.37629996e-09]
       [  0.00000000e+00  -1.37629996e-09   0.00000000e+00  -1.37629996e-09
          0.00000000e+00]]]]


