diff --git a/doc/pub/week36/html/week36-bs.html b/doc/pub/week36/html/week36-bs.html index 9def9c7a3..92eecfe0a 100644 --- a/doc/pub/week36/html/week36-bs.html +++ b/doc/pub/week36/html/week36-bs.html @@ -56,6 +56,10 @@ Automatically generated HTML file from DocOnce source 2, None, 'and-finally-boldsymbol-x-boldsymbol-x-t'), + ('Code for SVD and Inversion of Matrices', + 2, + None, + 'code-for-svd-and-inversion-of-matrices'), ('Ridge and LASSO Regression', 2, None, @@ -269,65 +273,66 @@ MathJax.Hub.Config({
  • Linear Regression and the SVD
  • What does it mean?
  • And finally \( \boldsymbol{X}\boldsymbol{X}^T \)
    • -
  • Ridge and LASSO Regression
    • -
  • From OLS to Ridge and Lasso
    • -
  • Deriving the Ridge Regression Equations
    • -
  • SVD analysis
    • -
  • Interpreting the Ridge results
    • -
  • More interpretations
    • -
  • Deriving the Lasso Regression Equations
    • -
  • Simple example to illustrate Ordinary Least Squares, Ridge and Lasso Regression
    • -
  • Ridge Regression
    • -
  • Lasso Regression
    • -
  • Yet another Example
    • -
  • The OLS case
    • -
  • The Ridge case
    • -
  • Writing the Cost Function
    • -
  • Lasso case
    • -
  • The first Case
    • -
  • Simple code for solving the above problem
    • -
  • With Lasso Regression
    • -
  • Another Example, now with a polynomial fit
    • -
  • Linking the regression analysis with a statistical interpretation
    • -
  • Assumptions made
    • -
  • Expectation value and variance
    • -
  • Expectation value and variance for \( \boldsymbol{\beta} \)
    • -
  • Deriving OLS from a probability distribution
    • -
  • Independent and Identically Distrubuted (iid)
    • -
  • Maximum Likelihood Estimation (MLE)
    • -
  • A new Cost Function
    • -
  • More basic Statistics and Bayes' theorem
    • -
  • Marginal Probability
    • -
  • Conditional Probability
    • -
  • Bayes' Theorem
    • -
  • Frliday September 10
    • -
  • Why resampling methods
    • -
  • Resampling methods
    • -
  • Resampling approaches can be computationally expensive
    • -
  • Why resampling methods ?
    • -
  • Statistical analysis
    • -
  • Resampling methods
    • -
  • Resampling methods: Jackknife and Bootstrap
    • -
  • Resampling methods: Jackknife
    • -
  • Jackknife code example
    • -
  • Resampling methods: Bootstrap
    • -
  • Resampling methods: Bootstrap background
    • -
  • Resampling methods: More Bootstrap background
    • -
  • Resampling methods: Bootstrap approach
    • -
  • Resampling methods: Bootstrap steps
    • -
  • Code example for the Bootstrap method
    • -
  • Various steps in cross-validation
    • -
  • How to set up the cross-validation for Ridge and/or Lasso
    • -
  • Cross-validation in brief
    • -
  • Code Example for Cross-validation and \( k \)-fold Cross-validation
    • -
  • The bias-variance tradeoff
    • -
  • Example code for Bias-Variance tradeoff
    • -
  • Understanding what happens
    • -
  • Summing up
    • -
  • Another Example from Scikit-Learn's Repository
    • -
  • More examples on bootstrap and cross-validation and errors
    • -
  • The same example but now with cross-validation
    • -
  • Cross-validation with Ridge
    • +
  • Code for SVD and Inversion of Matrices
    • +
  • Ridge and LASSO Regression
    • +
  • From OLS to Ridge and Lasso
    • +
  • Deriving the Ridge Regression Equations
    • +
  • SVD analysis
    • +
  • Interpreting the Ridge results
    • +
  • More interpretations
    • +
  • Deriving the Lasso Regression Equations
    • +
  • Simple example to illustrate Ordinary Least Squares, Ridge and Lasso Regression
    • +
  • Ridge Regression
    • +
  • Lasso Regression
    • +
  • Yet another Example
    • +
  • The OLS case
    • +
  • The Ridge case
    • +
  • Writing the Cost Function
    • +
  • Lasso case
    • +
  • The first Case
    • +
  • Simple code for solving the above problem
    • +
  • With Lasso Regression
    • +
  • Another Example, now with a polynomial fit
    • +
  • Linking the regression analysis with a statistical interpretation
    • +
  • Assumptions made
    • +
  • Expectation value and variance
    • +
  • Expectation value and variance for \( \boldsymbol{\beta} \)
    • +
  • Deriving OLS from a probability distribution
    • +
  • Independent and Identically Distrubuted (iid)
    • +
  • Maximum Likelihood Estimation (MLE)
    • +
  • A new Cost Function
    • +
  • More basic Statistics and Bayes' theorem
    • +
  • Marginal Probability
    • +
  • Conditional Probability
    • +
  • Bayes' Theorem
    • +
  • Frliday September 10
    • +
  • Why resampling methods
    • +
  • Resampling methods
    • +
  • Resampling approaches can be computationally expensive
    • +
  • Why resampling methods ?
    • +
  • Statistical analysis
    • +
  • Resampling methods
    • +
  • Resampling methods: Jackknife and Bootstrap
    • +
  • Resampling methods: Jackknife
    • +
  • Jackknife code example
    • +
  • Resampling methods: Bootstrap
    • +
  • Resampling methods: Bootstrap background
    • +
  • Resampling methods: More Bootstrap background
    • +
  • Resampling methods: Bootstrap approach
    • +
  • Resampling methods: Bootstrap steps
    • +
  • Code example for the Bootstrap method
    • +
  • Various steps in cross-validation
    • +
  • How to set up the cross-validation for Ridge and/or Lasso
    • +
  • Cross-validation in brief
    • +
  • Code Example for Cross-validation and \( k \)-fold Cross-validation
    • +
  • The bias-variance tradeoff
    • +
  • Example code for Bias-Variance tradeoff
    • +
  • Understanding what happens
    • +
  • Summing up
    • +
  • Another Example from Scikit-Learn's Repository
    • +
  • More examples on bootstrap and cross-validation and errors
    • +
  • The same example but now with cross-validation
    • +
  • Cross-validation with Ridge
    • @@ -386,7 +391,7 @@ MathJax.Hub.Config({
  • 9
  • 10
  • ...
    • -
  • 66
    • +
  • 67
  • »
    • diff --git a/doc/pub/week36/html/week36-reveal.html b/doc/pub/week36/html/week36-reveal.html index 2eff9beef..83b18fd30 100644 --- a/doc/pub/week36/html/week36-reveal.html +++ b/doc/pub/week36/html/week36-reveal.html @@ -328,6 +328,55 @@ values and the column vectors of \( \boldsymbol{V} \). +
    +

    Code for SVD and Inversion of Matrices

    + +

    +How do we use the SVD to invert a matrix \( \boldsymbol{X}^\boldsymbol{X} \) which is singular or near singular? +The simple answer is to use the linear algebra function for pseudoinvers, that is +

    + + +

    Ainv = np.linlag.pinv(A)
+
    +

    +Let us first look at a matrix which does not causes problems and write our own function where we just use the SVD. + +

    + + +

    import numpy as np
+# SVD inversion
+def SVDinv(A):
+    ''' Takes as input a numpy matrix A and returns inv(A) based on singular value decomposition (SVD).
+    SVD is numerically more stable than the inversion algorithms provided by
+    numpy and scipy.linalg at the cost of being slower.
+    '''
+    U, s, VT = np.linalg.svd(A)
+    print('test U')
+    print( (np.transpose(U) @ U - U @np.transpose(U)))
+    print('test VT')
+    print( (np.transpose(VT) @ VT - VT @np.transpose(VT)))
+
+
+    D = np.zeros((len(U),len(VT)))
+    D = np.diag(s)
+    UT = np.transpose(U); V = np.transpose(VT); invD = np.linalg.inv(D)
+    return np.matmul(V,np.matmul(invD,UT))
+
+
+#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]])
+print(X)
+A = np.transpose(X) @ X
+# Brute force inversion
+B = np.linalg.inv(A)
+C = SVDinv(A)
+print(np.abs(B-C))
+
    +
    + +

    Ridge and LASSO Regression

    @@ -977,14 +1026,14 @@ lambdas = np.logspace(-4, for i in range(nlambdas): lmb = lambdas[i] Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y -# print(Ridgebeta) + print(Ridgebeta) # and then make the prediction ypredictRidge = X @ Ridgebeta MSERidgePredict[i] = MSE(y,ypredictRidge) -# print(MSEPredict[i]) RegLasso = linear_model.Lasso(lmb) RegLasso.fit(X,y) ypredictLasso = RegLasso.predict(X) + print(RegLasso_coef_) MSELassoPredict[i] = MSE(y,ypredictLasso) # Now plot the results plt.figure() @@ -1040,13 +1089,9 @@ OLSbeta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train print(OLSbeta) # and then make the prediction ytildeOLS = X_train @ OLSbeta -print("Training R2 for OLS") -print(R2(y_train,ytildeOLS)) print("Training MSE for OLS") print(MSE(y_train,ytildeOLS)) ypredictOLS = X_test @ OLSbeta -print("Test R2 for OLS") -print(R2(y_test,ypredictOLS)) print("Test MSE OLS") print(MSE(y_test,ypredictOLS)) diff --git a/doc/pub/week36/html/week36-solarized.html b/doc/pub/week36/html/week36-solarized.html index fc95d417d..3b9c421d7 100644 --- a/doc/pub/week36/html/week36-solarized.html +++ b/doc/pub/week36/html/week36-solarized.html @@ -76,6 +76,10 @@ div { text-align: justify; text-justify: inter-word; } 2, None, 'and-finally-boldsymbol-x-boldsymbol-x-t'), + ('Code for SVD and Inversion of Matrices', + 2, + None, + 'code-for-svd-and-inversion-of-matrices'), ('Ridge and LASSO Regression', 2, None, @@ -438,6 +442,54 @@ values and the column vectors of \( \boldsymbol{V} \).











    +

    Code for SVD and Inversion of Matrices

    + +

    +How do we use the SVD to invert a matrix \( \boldsymbol{X}^\boldsymbol{X} \) which is singular or near singular? +The simple answer is to use the linear algebra function for pseudoinvers, that is +

    + + +

    Ainv = np.linlag.pinv(A)
+
    +

    +Let us first look at a matrix which does not causes problems and write our own function where we just use the SVD. + +

    + + +

    import numpy as np
+# SVD inversion
+def SVDinv(A):
+    ''' Takes as input a numpy matrix A and returns inv(A) based on singular value decomposition (SVD).
+    SVD is numerically more stable than the inversion algorithms provided by
+    numpy and scipy.linalg at the cost of being slower.
+    '''
+    U, s, VT = np.linalg.svd(A)
+    print('test U')
+    print( (np.transpose(U) @ U - U @np.transpose(U)))
+    print('test VT')
+    print( (np.transpose(VT) @ VT - VT @np.transpose(VT)))
+
+
+    D = np.zeros((len(U),len(VT)))
+    D = np.diag(s)
+    UT = np.transpose(U); V = np.transpose(VT); invD = np.linalg.inv(D)
+    return np.matmul(V,np.matmul(invD,UT))
+
+
+#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]])
+print(X)
+A = np.transpose(X) @ X
+# Brute force inversion
+B = np.linalg.inv(A)
+C = SVDinv(A)
+print(np.abs(B-C))
+
    +

    +









    +

    Ridge and LASSO Regression

    @@ -995,14 +1047,14 @@ lambdas = np.logspace(-4, for i in range(nlambdas): lmb = lambdas[i] Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y -# print(Ridgebeta) + print(Ridgebeta) # and then make the prediction ypredictRidge = X @ Ridgebeta MSERidgePredict[i] = MSE(y,ypredictRidge) -# print(MSEPredict[i]) RegLasso = linear_model.Lasso(lmb) RegLasso.fit(X,y) ypredictLasso = RegLasso.predict(X) + print(RegLasso_coef_) MSELassoPredict[i] = MSE(y,ypredictLasso) # Now plot the results plt.figure() @@ -1057,13 +1109,9 @@ OLSbeta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train print(OLSbeta) # and then make the prediction ytildeOLS = X_train @ OLSbeta -print("Training R2 for OLS") -print(R2(y_train,ytildeOLS)) print("Training MSE for OLS") print(MSE(y_train,ytildeOLS)) ypredictOLS = X_test @ OLSbeta -print("Test R2 for OLS") -print(R2(y_test,ypredictOLS)) print("Test MSE OLS") print(MSE(y_test,ypredictOLS)) diff --git a/doc/pub/week36/html/week36.html b/doc/pub/week36/html/week36.html index ba6797338..69f561b9b 100644 --- a/doc/pub/week36/html/week36.html +++ b/doc/pub/week36/html/week36.html @@ -81,6 +81,10 @@ div { text-align: justify; text-justify: inter-word; } 2, None, 'and-finally-boldsymbol-x-boldsymbol-x-t'), + ('Code for SVD and Inversion of Matrices', + 2, + None, + 'code-for-svd-and-inversion-of-matrices'), ('Ridge and LASSO Regression', 2, None, @@ -443,6 +447,54 @@ values and the column vectors of \( \boldsymbol{V} \).











    +

    Code for SVD and Inversion of Matrices

    + +

    +How do we use the SVD to invert a matrix \( \boldsymbol{X}^\boldsymbol{X} \) which is singular or near singular? +The simple answer is to use the linear algebra function for pseudoinvers, that is +

    + + +

    Ainv = np.linlag.pinv(A)
+
    +

    +Let us first look at a matrix which does not causes problems and write our own function where we just use the SVD. + +

    + + +

    import numpy as np
+# SVD inversion
+def SVDinv(A):
+    ''' Takes as input a numpy matrix A and returns inv(A) based on singular value decomposition (SVD).
+    SVD is numerically more stable than the inversion algorithms provided by
+    numpy and scipy.linalg at the cost of being slower.
+    '''
+    U, s, VT = np.linalg.svd(A)
+    print('test U')
+    print( (np.transpose(U) @ U - U @np.transpose(U)))
+    print('test VT')
+    print( (np.transpose(VT) @ VT - VT @np.transpose(VT)))
+
+
+    D = np.zeros((len(U),len(VT)))
+    D = np.diag(s)
+    UT = np.transpose(U); V = np.transpose(VT); invD = np.linalg.inv(D)
+    return np.matmul(V,np.matmul(invD,UT))
+
+
+#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]])
+print(X)
+A = np.transpose(X) @ X
+# Brute force inversion
+B = np.linalg.inv(A)
+C = SVDinv(A)
+print(np.abs(B-C))
+
    +

    +









    +

    Ridge and LASSO Regression

    @@ -1000,14 +1052,14 @@ lambdas = np.for i in range(nlambdas): lmb = lambdas[i] Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y -# print(Ridgebeta) + print(Ridgebeta) # and then make the prediction ypredictRidge = X @ Ridgebeta MSERidgePredict[i] = MSE(y,ypredictRidge) -# print(MSEPredict[i]) RegLasso = linear_model.Lasso(lmb) RegLasso.fit(X,y) ypredictLasso = RegLasso.predict(X) + print(RegLasso_coef_) MSELassoPredict[i] = MSE(y,ypredictLasso) # Now plot the results plt.figure() @@ -1062,13 +1114,9 @@ OLSbeta = np.print(OLSbeta) # and then make the prediction ytildeOLS = X_train @ OLSbeta -print("Training R2 for OLS") -print(R2(y_train,ytildeOLS)) print("Training MSE for OLS") print(MSE(y_train,ytildeOLS)) ypredictOLS = X_test @ OLSbeta -print("Test R2 for OLS") -print(R2(y_test,ypredictOLS)) print("Test MSE OLS") print(MSE(y_test,ypredictOLS)) diff --git a/doc/pub/week36/ipynb/ipynb-week36-src.tar.gz b/doc/pub/week36/ipynb/ipynb-week36-src.tar.gz index c07b5e271..ce546c0f8 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 9395e12c4..a17bc5486 100644 --- a/doc/pub/week36/ipynb/week36.ipynb +++ b/doc/pub/week36/ipynb/week36.ipynb @@ -254,6 +254,75 @@ "values and the column vectors of $\\boldsymbol{V}$.\n", "\n", "\n", + "\n", + "## Code for SVD and Inversion of Matrices\n", + "\n", + "How do we use the SVD to invert a matrix $\\boldsymbol{X}^\\boldsymbol{X}$ which is singular or near singular?\n", + "The simple answer is to use the linear algebra function for pseudoinvers, that is" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "Ainv = np.linlag.pinv(A)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let us first look at a matrix which does not causes problems and write our own function where we just use 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", + " ''' Takes as input a numpy matrix A and returns inv(A) based on singular value decomposition (SVD).\n", + " SVD is numerically more stable than the inversion algorithms provided by\n", + " numpy and scipy.linalg at the cost of being slower.\n", + " '''\n", + " U, s, VT = np.linalg.svd(A)\n", + " print('test U')\n", + " print( (np.transpose(U) @ U - U @np.transpose(U)))\n", + " print('test VT')\n", + " print( (np.transpose(VT) @ VT - VT @np.transpose(VT)))\n", + "\n", + "\n", + " D = np.zeros((len(U),len(VT)))\n", + " D = np.diag(s)\n", + " UT = np.transpose(U); V = np.transpose(VT); invD = np.linalg.inv(D)\n", + " return np.matmul(V,np.matmul(invD,UT))\n", + "\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", + "print(X)\n", + "A = np.transpose(X) @ X\n", + "# Brute force inversion\n", + "B = np.linalg.inv(A)\n", + "C = SVDinv(A)\n", + "print(np.abs(B-C))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ "## 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", @@ -1222,14 +1291,14 @@ "for i in range(nlambdas):\n", " lmb = lambdas[i]\n", " Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y\n", - "# print(Ridgebeta)\n", + " print(Ridgebeta)\n", " # and then make the prediction\n", " ypredictRidge = X @ Ridgebeta\n", " MSERidgePredict[i] = MSE(y,ypredictRidge)\n", - "# print(MSEPredict[i])\n", " RegLasso = linear_model.Lasso(lmb)\n", " RegLasso.fit(X,y)\n", " ypredictLasso = RegLasso.predict(X)\n", + " print(RegLasso_coef_)\n", " MSELassoPredict[i] = MSE(y,ypredictLasso)\n", "# Now plot the results\n", "plt.figure()\n", @@ -1293,13 +1362,9 @@ "print(OLSbeta)\n", "# and then make the prediction\n", "ytildeOLS = X_train @ OLSbeta\n", - "print(\"Training R2 for OLS\")\n", - "print(R2(y_train,ytildeOLS))\n", "print(\"Training MSE for OLS\")\n", "print(MSE(y_train,ytildeOLS))\n", "ypredictOLS = X_test @ OLSbeta\n", - "print(\"Test R2 for OLS\")\n", - "print(R2(y_test,ypredictOLS))\n", "print(\"Test MSE OLS\")\n", "print(MSE(y_test,ypredictOLS))\n", "\n", diff --git a/doc/src/week36/week36.do.txt b/doc/src/week36/week36.do.txt index 201bbac0f..c3d0d4b3f 100644 --- a/doc/src/week36/week36.do.txt +++ b/doc/src/week36/week36.do.txt @@ -144,6 +144,52 @@ of our data (the columns of $\bm{X}$, the quantity of interest for us are the no values and the column vectors of $\bm{V}$. + +!split +===== Code for SVD and Inversion of Matrices ===== + +How do we use the SVD to invert a matrix $\bm{X}^\bm{X}$ which is singular or near singular? +The simple answer is to use the linear algebra function for pseudoinvers, that is +!bc pycod +Ainv = np.linlag.pinv(A) +!ec + +Let us first look at a matrix which does not causes problems and write our own function where we just use the SVD. + +!bc pycod +import numpy as np +# SVD inversion +def SVDinv(A): + ''' Takes as input a numpy matrix A and returns inv(A) based on singular value decomposition (SVD). + SVD is numerically more stable than the inversion algorithms provided by + numpy and scipy.linalg at the cost of being slower. + ''' + U, s, VT = np.linalg.svd(A) + print('test U') + print( (np.transpose(U) @ U - U @np.transpose(U))) + print('test VT') + print( (np.transpose(VT) @ VT - VT @np.transpose(VT))) + + + D = np.zeros((len(U),len(VT))) + D = np.diag(s) + UT = np.transpose(U); V = np.transpose(VT); invD = np.linalg.inv(D) + return np.matmul(V,np.matmul(invD,UT)) + + +#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]]) +print(X) +A = np.transpose(X) @ X +# Brute force inversion +B = np.linalg.inv(A) +C = SVDinv(A) +print(np.abs(B-C)) +!ec + + + + !split ===== Ridge and LASSO Regression ===== @@ -698,14 +744,14 @@ lambdas = np.logspace(-4, 4, nlambdas) for i in range(nlambdas): lmb = lambdas[i] Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y -# print(Ridgebeta) + print(Ridgebeta) # and then make the prediction ypredictRidge = X @ Ridgebeta MSERidgePredict[i] = MSE(y,ypredictRidge) -# print(MSEPredict[i]) RegLasso = linear_model.Lasso(lmb) RegLasso.fit(X,y) ypredictLasso = RegLasso.predict(X) + print(RegLasso_coef_) MSELassoPredict[i] = MSE(y,ypredictLasso) # Now plot the results plt.figure() @@ -758,13 +804,9 @@ OLSbeta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train print(OLSbeta) # and then make the prediction ytildeOLS = X_train @ OLSbeta -print("Training R2 for OLS") -print(R2(y_train,ytildeOLS)) print("Training MSE for OLS") print(MSE(y_train,ytildeOLS)) ypredictOLS = X_test @ OLSbeta -print("Test R2 for OLS") -print(R2(y_test,ypredictOLS)) print("Test MSE OLS") print(MSE(y_test,ypredictOLS))