final update for lecture Thursday

doc/src/week36/week36.do.txt

@@ -34,7 +34,7 @@ We used the SVD to analyse the matrix to invert in ordinary lineat regression
 Since the matrices here have dimension $p\times p$, with $p$ corresponding to the singular values, we defined last week the matrix
 !bt
 \[
-\bm{\Sigma}^T\bm{\Sigma} = \begin{bmatrix} \tilde{\bm{\Sigma}} & \bm{0}\\ \end{bmatrix}\begin{bmatrix} \tilde{\bm{\Sigma}} \\ \bm{0}\\ \end{bmatrix},
+\bm{\Sigma}^T\bm{\Sigma} = \begin{bmatrix} \tilde{\bm{\Sigma}} & \bm{0}\\ \end{bmatrix}\begin{bmatrix} \tilde{\bm{\Sigma}} \\ \bm{0}\end{bmatrix},
 \]
 !et
 where the tilde-matrix $\tilde{\bm{\Sigma}}$ is a matrix of dimension $p\times p$ containing only the singular values $\sigma_i$, that is
@@ -178,15 +178,64 @@ def SVDinv(A):
 
 
 #X = np.array([ [1.0, -1.0, 2.0], [1.0, 0.0, 1.0], [1.0, 2.0, -1.0], [1.0, 1.0, 0.0] ])
-X = np.array( [ [1,2],[2,3]])
+# Non-singular square matrix
+X = np.array( [ [1,2,3],[2,4,5],[3,5,6]])
 print(X)
 A = np.transpose(X) @ X
 # Brute force inversion
-B = np.linalg.inv(A)
+B = np.linalg.inv(A)  # here we could use np.linalg.pinv(A)
 C = SVDinv(A)
 print(np.abs(B-C))
+
 !ec
 
+!split
+===== Inverse of Rectangular Matrix =====
+
+Although our matrix to invert $\bm{X}^T\bm{X}$ is a square matrix, our matrix may be singular. 
+
+The pseudoinverse is the generalization of the matrix inverse for square matrices to
+rectangular matrices where the number of rows and columns are not equal.
+
+It is also called the the Moore-Penrose Inverse after two independent discoverers of the method or the Generalized Inverse.
+It is used for the calculation of the inverse for singular or near singular matrices and for rectangular matrices.
+
+Using the SVD we can obtain the pseudoinverse of a matrix $\bm{A}$ (labeled here as $\bm{A}_{\mathrm{PI}}$
+!bt
+\[
+\bm{A}_{\mathrm{PI}}= \bm{V}\bm{D}_{\mathrm{PI}}\bm{U}^T,
+\]
+!et
+where $\bm{D}_{\mathrm{PI}}$ can be calculated by creating a diagonal matrix from $\bm{Sigma}$ where we only keep the singular values (the non-zero values). The following code computes the pseudoinvers of the matrix based on the SVD.
+
+
+!bc pycod
+import numpy as np
+# SVD inversion
+def SVDinv(A):
+    U, s, VT = np.linalg.svd(A)
+    # reciprocals of singular values of s
+    d = 1.0 / s
+    # create m x n D matrix
+    D = np.zeros(A.shape)
+    # populate D with n x n diagonal matrix
+    D[:A.shape[1], :A.shape[1]] = np.diag(d)
+    UT = np.transpose(U)
+    V = np.transpose(VT)
+    return np.matmul(V,np.matmul(D.T,UT))
+
+
+A = np.array([ [0.3, 0.4], [0.5, 0.6], [0.7, 0.8],[0.9, 1.0]])
+print(A)
+# Brute force inversion of super-collinear matrix
+B = np.linalg.pinv(A)
+print(B)
+# Compare our own algorithm with pinv
+C = SVDinv(A)
+print(np.abs(C-B))
+
+!ec
+As you can see from this example, our own decomposition based on the SVD agrees the pseudoinverse algorithm provided by _Numpy_.