$$
-\boldsymbol{\Sigma}^T\boldsymbol{\Sigma} = \begin{bmatrix} \tilde{\boldsymbol{\Sigma}} & \boldsymbol{0}\\ \end{bmatrix}\begin{bmatrix} \tilde{\boldsymbol{\Sigma}} \\ \boldsymbol{0}\\ \end{bmatrix},
+\boldsymbol{\Sigma}^T\boldsymbol{\Sigma} = \begin{bmatrix} \tilde{\boldsymbol{\Sigma}} & \boldsymbol{0}\\ \end{bmatrix}\begin{bmatrix} \tilde{\boldsymbol{\Sigma}} \\ \boldsymbol{0}\end{bmatrix},
$$
+Although our matrix to invert \( \boldsymbol{X}^T\boldsymbol{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 \( \boldsymbol{A} \) (labeled here as \( \boldsymbol{A}_{\mathrm{PI}} \)
+
+
+
+
+As you can see from this example, our own decomposition based on the SVD agrees the pseudoinverse algorithm provided by Numpy.
+
+Although our matrix to invert \( \boldsymbol{X}^T\boldsymbol{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 \( \boldsymbol{A} \) (labeled here as \( \boldsymbol{A}_{\mathrm{PI}} \)
+$$
+\boldsymbol{A}_{\mathrm{PI}}= \boldsymbol{V}\boldsymbol{D}_{\mathrm{PI}}\boldsymbol{U}^T,
+$$
+
+where \( \boldsymbol{D}_{\mathrm{PI}} \) can be calculated by creating a diagonal matrix from \( \boldsymbol{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.
+
+
+
+
+
+As you can see from this example, our own decomposition based on the SVD agrees the pseudoinverse algorithm provided by Numpy.
+
+
+
diff --git a/doc/pub/week36/html/week36.html b/doc/pub/week36/html/week36.html
index 69f561b9b..b3cf5190d 100644
--- a/doc/pub/week36/html/week36.html
+++ b/doc/pub/week36/html/week36.html
@@ -85,6 +85,10 @@ div { text-align: justify; text-justify: inter-word; }
2,
None,
'code-for-svd-and-inversion-of-matrices'),
+ ('Inverse of Rectangular Matrix',
+ 2,
+ None,
+ 'inverse-of-rectangular-matrix'),
('Ridge and LASSO Regression',
2,
None,
@@ -338,7 +342,7 @@ $$
Since the matrices here have dimension \( p\times p \), with \( p \) corresponding to the singular values, we defined last week the matrix
$$
-\boldsymbol{\Sigma}^T\boldsymbol{\Sigma} = \begin{bmatrix} \tilde{\boldsymbol{\Sigma}} & \boldsymbol{0}\\ \end{bmatrix}\begin{bmatrix} \tilde{\boldsymbol{\Sigma}} \\ \boldsymbol{0}\\ \end{bmatrix},
+\boldsymbol{\Sigma}^T\boldsymbol{\Sigma} = \begin{bmatrix} \tilde{\boldsymbol{\Sigma}} & \boldsymbol{0}\\ \end{bmatrix}\begin{bmatrix} \tilde{\boldsymbol{\Sigma}} \\ \boldsymbol{0}\end{bmatrix},
$$
where the tilde-matrix \( \tilde{\boldsymbol{\Sigma}} \) is a matrix of dimension \( p\times p \) containing only the singular values \( \sigma_i \), that is
@@ -484,17 +488,72 @@ Let us first look at a matrix which does not causes problems and write our own f
#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))
+Although our matrix to invert \( \boldsymbol{X}^T\boldsymbol{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 \( \boldsymbol{A} \) (labeled here as \( \boldsymbol{A}_{\mathrm{PI}} \)
+$$
+\boldsymbol{A}_{\mathrm{PI}}= \boldsymbol{V}\boldsymbol{D}_{\mathrm{PI}}\boldsymbol{U}^T,
+$$
+
+where \( \boldsymbol{D}_{\mathrm{PI}} \) can be calculated by creating a diagonal matrix from \( \boldsymbol{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.
+
+
+
+
+
+As you can see from this example, our own decomposition based on the SVD agrees the pseudoinverse algorithm provided by Numpy.
+
+
+
diff --git a/doc/pub/week36/ipynb/ipynb-week36-src.tar.gz b/doc/pub/week36/ipynb/ipynb-week36-src.tar.gz
index ce546c0f8..b5f94b110 100644
Binary files a/doc/pub/week36/ipynb/ipynb-week36-src.tar.gz and b/doc/pub/week36/ipynb/ipynb-week36-src.tar.gz differ
diff --git a/doc/pub/week36/ipynb/week36.ipynb b/doc/pub/week36/ipynb/week36.ipynb
index a17bc5486..ee2d20f06 100644
--- a/doc/pub/week36/ipynb/week36.ipynb
+++ b/doc/pub/week36/ipynb/week36.ipynb
@@ -61,7 +61,7 @@
"metadata": {},
"source": [
"$$\n",
- "\\boldsymbol{\\Sigma}^T\\boldsymbol{\\Sigma} = \\begin{bmatrix} \\tilde{\\boldsymbol{\\Sigma}} & \\boldsymbol{0}\\\\ \\end{bmatrix}\\begin{bmatrix} \\tilde{\\boldsymbol{\\Sigma}} \\\\ \\boldsymbol{0}\\\\ \\end{bmatrix},\n",
+ "\\boldsymbol{\\Sigma}^T\\boldsymbol{\\Sigma} = \\begin{bmatrix} \\tilde{\\boldsymbol{\\Sigma}} & \\boldsymbol{0}\\\\ \\end{bmatrix}\\begin{bmatrix} \\tilde{\\boldsymbol{\\Sigma}} \\\\ \\boldsymbol{0}\\end{bmatrix},\n",
"$$"
]
},
@@ -310,11 +310,12 @@
"\n",
"\n",
"#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] ])\n",
- "X = np.array( [ [1,2],[2,3]])\n",
+ "# Non-singular square matrix\n",
+ "X = np.array( [ [1,2,3],[2,4,5],[3,5,6]])\n",
"print(X)\n",
"A = np.transpose(X) @ X\n",
"# Brute force inversion\n",
- "B = np.linalg.inv(A)\n",
+ "B = np.linalg.inv(A) # here we could use np.linalg.pinv(A)\n",
"C = SVDinv(A)\n",
"print(np.abs(B-C))"
]
@@ -323,6 +324,77 @@
"cell_type": "markdown",
"metadata": {},
"source": [
+ "## Inverse of Rectangular Matrix\n",
+ "\n",
+ "Although our matrix to invert $\\boldsymbol{X}^T\\boldsymbol{X}$ is a square matrix, our matrix may be singular. \n",
+ "\n",
+ "The pseudoinverse is the generalization of the matrix inverse for square matrices to\n",
+ "rectangular matrices where the number of rows and columns are not equal.\n",
+ "\n",
+ "It is also called the the Moore-Penrose Inverse after two independent discoverers of the method or the Generalized Inverse.\n",
+ "It is used for the calculation of the inverse for singular or near singular matrices and for rectangular matrices.\n",
+ "\n",
+ "Using the SVD we can obtain the pseudoinverse of a matrix $\\boldsymbol{A}$ (labeled here as $\\boldsymbol{A}_{\\mathrm{PI}}$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\boldsymbol{A}_{\\mathrm{PI}}= \\boldsymbol{V}\\boldsymbol{D}_{\\mathrm{PI}}\\boldsymbol{U}^T,\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where $\\boldsymbol{D}_{\\mathrm{PI}}$ can be calculated by creating a diagonal matrix from $\\boldsymbol{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."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false,
+ "editable": true
+ },
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "# SVD inversion\n",
+ "def SVDinv(A):\n",
+ " U, s, VT = np.linalg.svd(A)\n",
+ " # reciprocals of singular values of s\n",
+ " d = 1.0 / s\n",
+ " # create m x n D matrix\n",
+ " D = np.zeros(A.shape)\n",
+ " # populate D with n x n diagonal matrix\n",
+ " D[:A.shape[1], :A.shape[1]] = np.diag(d)\n",
+ " UT = np.transpose(U)\n",
+ " V = np.transpose(VT)\n",
+ " return np.matmul(V,np.matmul(D.T,UT))\n",
+ "\n",
+ "\n",
+ "A = np.array([ [0.3, 0.4], [0.5, 0.6], [0.7, 0.8],[0.9, 1.0]])\n",
+ "print(A)\n",
+ "# Brute force inversion of super-collinear matrix\n",
+ "B = np.linalg.pinv(A)\n",
+ "print(B)\n",
+ "# Compare our own algorithm with pinv\n",
+ "C = SVDinv(A)\n",
+ "print(np.abs(C-B))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "As you can see from this example, our own decomposition based on the SVD agrees the pseudoinverse algorithm provided by **Numpy**.\n",
+ "\n",
+ "\n",
+ "\n",
"## Ridge and LASSO Regression\n",
"\n",
"Let us remind ourselves about the expression for the standard Mean Squared Error (MSE) which we used to define our cost function and the equations for the ordinary least squares (OLS) method, that is \n",
diff --git a/doc/src/week36/week36.do.txt b/doc/src/week36/week36.do.txt
index c3d0d4b3f..96697790e 100644
--- a/doc/src/week36/week36.do.txt
+++ b/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_.
@@ -366,17 +366,74 @@ Let us first look at a matrix which does not causes problems and write our own f
#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))
+Inverse of Rectangular Matrix
+
+
+$$
+\boldsymbol{A}_{\mathrm{PI}}= \boldsymbol{V}\boldsymbol{D}_{\mathrm{PI}}\boldsymbol{U}^T,
+$$
+
+
+where \( \boldsymbol{D}_{\mathrm{PI}} \) can be calculated by creating a diagonal matrix from \( \boldsymbol{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.
+
+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))
+
Ridge and LASSO Regression
diff --git a/doc/pub/week36/html/week36-solarized.html b/doc/pub/week36/html/week36-solarized.html
index 3b9c421d7..2d52143dc 100644
--- a/doc/pub/week36/html/week36-solarized.html
+++ b/doc/pub/week36/html/week36-solarized.html
@@ -80,6 +80,10 @@ div { text-align: justify; text-justify: inter-word; }
2,
None,
'code-for-svd-and-inversion-of-matrices'),
+ ('Inverse of Rectangular Matrix',
+ 2,
+ None,
+ 'inverse-of-rectangular-matrix'),
('Ridge and LASSO Regression',
2,
None,
@@ -333,7 +337,7 @@ $$
Since the matrices here have dimension \( p\times p \), with \( p \) corresponding to the singular values, we defined last week the matrix
$$
-\boldsymbol{\Sigma}^T\boldsymbol{\Sigma} = \begin{bmatrix} \tilde{\boldsymbol{\Sigma}} & \boldsymbol{0}\\ \end{bmatrix}\begin{bmatrix} \tilde{\boldsymbol{\Sigma}} \\ \boldsymbol{0}\\ \end{bmatrix},
+\boldsymbol{\Sigma}^T\boldsymbol{\Sigma} = \begin{bmatrix} \tilde{\boldsymbol{\Sigma}} & \boldsymbol{0}\\ \end{bmatrix}\begin{bmatrix} \tilde{\boldsymbol{\Sigma}} \\ \boldsymbol{0}\end{bmatrix},
$$
where the tilde-matrix \( \tilde{\boldsymbol{\Sigma}} \) is a matrix of dimension \( p\times p \) containing only the singular values \( \sigma_i \), that is
@@ -479,17 +483,72 @@ Let us first look at a matrix which does not causes problems and write our own f
#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))
+Inverse of Rectangular Matrix
+
+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))
+
+
Ridge and LASSO Regression
+Inverse of Rectangular Matrix
+
+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))
+
+
Ridge and LASSO Regression