From f6d429b2ce6665033d69449c91fe896f19413a21 Mon Sep 17 00:00:00 2001 From: Morten Hjorth-Jensen Date: Thu, 9 Sep 2021 09:42:41 +0200 Subject: [PATCH] added svd algo --- doc/pub/week36/html/week36-bs.html | 125 ++++++++++--------- doc/pub/week36/html/week36-reveal.html | 57 ++++++++- doc/pub/week36/html/week36-solarized.html | 60 ++++++++- doc/pub/week36/html/week36.html | 60 ++++++++- doc/pub/week36/ipynb/ipynb-week36-src.tar.gz | Bin 192 -> 191 bytes doc/pub/week36/ipynb/week36.ipynb | 77 +++++++++++- doc/src/week36/week36.do.txt | 54 +++++++- 7 files changed, 343 insertions(+), 90 deletions(-) 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 c07b5e271a7db591519421c96b786d43f0383c4e..ce546c0f8d71fe5f6b6f817a458fd4f7e4ee0a7f 100644 GIT binary patch literal 191 zcmV;w06_mAiwFSax;bG01MSbv3c@f92k@Qu6nTQt+IF4=J$MjAe1VQ-uFkcy?a;lw z`v5&DUWN$$UH*iGkQ_Fv^=^~+yN_l=2uV4IA(yEblPsosM5zGEDNPv@j0+&NPZ9x{ z@1&R3d1?Do>hg@*p?){lj}^s-J<}`j%s+7~m4ofR^QFs# literal 192 zcmV;x06+g9iwFRVwmD${1MSbv3c@f92k@Qu6nTQtu3J3|?%+WX@dY}TxjNU*wnO*! z?gR9sco`z}cli?%LUPE~n_U*Uy9*XW2;oTrW3G5QC5h*HLTL&#p_DVmV#_Wl|}V#ZWt@e4}0cU;F*8oSV;@peeWu*Kxv1$)-~J^YsMql_9}-$ uqa9md@Y+d(AXE>cD5R5GiA&fTeKMl4QTXd;JkRsIuRQ<-`4^@D2mk