De matrix $$\mathbf{A}$$ is reëel en telt $$M$$ rijen en $$N$$ kolommen ($$\mathbf{A}$$ is dus een $$M \times N$$ matrix). Een (niet noodzakelijk vierkante) matrix $$\mathbf{A^+}$$, die voldoet aan

$$\mathbf{AA^+A} = \mathbf{A}$$

wordt een pseudo-inverse van $$\mathbf{A}$$ genoemd. Het bestaan van minstens 1 dergelijke matrix is verzekerd.

In deze oefening willen we deze matrix $$\mathbf{A^+}$$ berekenen, en wel op een iteratieve manier. Hiertoe construeren we een rij van matrices $$\mathbf{A_i}$$, gedefinieerd via:

$$\mathbf{A_{i+1}} = \mathbf{2A_i} - \mathbf{A_i A A_i}$$

Als startwaarde voor de iteratie kiezen we voor de matrix $$\mathbf{A_0}$$:

$$\mathbf{A_0} = (\mathbf{A^T A} + \delta\mathbf{I_N})^{-1}\mathbf{A^T}$$

Hierin stelt $$\delta$$ een positief reëel getal voor, $$\mathbf{I_N}$$ de eenheidsmatrix met $$N$$ rijen en $$N$$ kolommen, en $$\mathbf{A^T}$$ de getransponeerde van de originele matrix $$\mathbf{A}$$.

Schrijf een functie pseudo_inverse() met drie argumenten, namelijk:

Het laatste argument geeft aan wanneer de iteratie moet eindigen. Dit is het geval indien

$$||\mathbf{AA^+A} - \mathbf{A}|| < eps$$

Hierin stelt $$||\mathbf{X}||$$ de norm van de matrix $$||\mathbf{X}||$$ voor, gedefinieerd als de vierkantswortel uit de som van de kwadraten van alle elementen van $$\mathbf{X}$$. Deze norm kan je berekenen via de functie np.linalg.norm() (met als enig argument de 2-dimensionale NumPy-rij die de matrix in kwestie voorstelt. Het resultaat van de functie is een 2-dimensionale NumPy-tabel (type array), die de berekende pseudo-inverse van $$\mathbf{A}$$ voorstelt.

Voorbeeld

C =  np.array([[ 0.,  1.,  2.,  3.], 
[-1., -2., 5., 7.], 
[3., 2., 9., 1.], 
[-5., -3., 2., 0.]]) 

print(pseudo_inverse(C, 0.1, 0.01))
#[[-0.38415626  0.15852859  0.04299677 -0.20566935]
#[ 0.62303578 -0.26681444 -0.00175955  0.05194953]
#[-0.02616147 -0.00376664  0.10486554  0.06368012]
#[ 0.14200467  0.09188593 -0.06927242 -0.05998626]]