diff --git a/doc/pub/Regression/html/._Regression-bs000.html b/doc/pub/Regression/html/._Regression-bs000.html
index 833e9dec8..e6ce300e7 100644
--- a/doc/pub/Regression/html/._Regression-bs000.html
+++ b/doc/pub/Regression/html/._Regression-bs000.html
@@ -225,7 +225,34 @@ Automatically generated HTML file from DocOnce source
('Elaborating a little bit more', 2, None, '___sec98'),
('The bias', 2, None, '___sec99'),
('The variance', 2, None, '___sec100'),
- ('Summing up', 2, None, '___sec101')]}
+ ('Summing up', 2, None, '___sec101'),
+ ('The one-dimensional Ising model, project 2',
+ 2,
+ None,
+ '___sec102'),
+ ('Example: The one-dimensional Ising model',
+ 2,
+ None,
+ '___sec103'),
+ ('Reformulating the problem to suit regression',
+ 2,
+ None,
+ '___sec104'),
+ ('Linear regression', 2, None, '___sec105'),
+ ('Ordinary least squares', 2, None, '___sec106'),
+ ('Singular Value decomposition', 2, None, '___sec107'),
+ ('Fitting with scikit-learn', 2, None, '___sec108'),
+ ('Ridge regression', 2, None, '___sec109'),
+ ('LASSO regression', 2, None, '___sec110'),
+ ('Performance of the different models', 2, None, '___sec111'),
+ ('Performance as function of the regularization parameter',
+ 2,
+ None,
+ '___sec112'),
+ ('Finding the optimal value of $\\lambda$',
+ 2,
+ None,
+ '___sec113')]}
end of tocinfo -->
@@ -365,6 +392,18 @@ MathJax.Hub.Config({
The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) [2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University
-
Sep 25, 2018 Oct 11, 2018
@@ -423,7 +462,7 @@ MathJax.Hub.Config({
9 10 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 10 11 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 11 12 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 12 13 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 13 14 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 14 15 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 15 16 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 16 17 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 17 18 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 18 19 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 19 20 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 20 21 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 21 22 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 22 23 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 23 24 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 24 25 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 25 26 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 26 27 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 27 28 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 28 29 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) . show()
29 30 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) . show()
30 31 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 31 32 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) cost function.
32 33 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 33 34 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 34 35 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 35 36 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) . show()
36 37 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 37 38 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 38 39 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 39 40 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) R , we can perform similar studies.
40 41 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) . show()
41 42 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 42 43 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 43 44 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 44 45 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 45 46 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 46 47 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 47 48 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 48 49 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) . show()
49 50 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) . show()
50 51 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 51 52 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) . show()
52 53 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) . show()
53 54 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 54 55 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 55 56 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 56 57 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 57 58 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 58 59 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 59 60 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) spread of \( p \) around its mean.
60 61 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 61 62 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 62 63 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 63 64 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 64 65 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 65 66 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 66 67 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 67 68 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 68 69 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 69 70 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 70 71 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 71 72 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 72 73 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 73 74 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 74 75 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 75 76 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 76 77 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 77 78 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 78 79 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 79 80 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 80 81 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 81 82 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 82 83 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 83 84 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 84 85 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 85 86 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 86 87 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 87 88 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 88 89 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 89 90 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) = jackknife(x, stat)
90 91 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 91 92 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 92 93 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 93 94 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 94 95 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 95 96 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) . show()
96 97 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 97 98 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 98 99 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 99 100 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 100 101 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 101 102 ... 103 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) . printOutput()
101 102 103 ... 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 101 102 103 104 ... 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 101 102 103 104 105 ... 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 101 102 103 104 105 106 ... 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 101 102 103 104 105 106 107 ... 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 101 102 103 104 105 106 107 108 ... 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 101 102 103 104 105 106 107 108 109 ... 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 101 102 103 104 105 106 107 108 109 110 ... 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) 101 102 103 104 105 106 107 108 109 110 111 ... 115 » The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \)
diff --git a/doc/pub/Regression/html/._Regression-bs103.html b/doc/pub/Regression/html/._Regression-bs103.html
index 0aa98b52c..32bf8052f 100644
--- a/doc/pub/Regression/html/._Regression-bs103.html
+++ b/doc/pub/Regression/html/._Regression-bs103.html
@@ -194,46 +194,65 @@ Automatically generated HTML file from DocOnce source
'___sec77'),
('Resampling methods: Jackknife', 2, None, '___sec78'),
('Resampling methods: Jackknife estimator', 2, None, '___sec79'),
- ('Resampling methods: Bootstrap', 2, None, '___sec80'),
- ('Resampling methods: Bootstrap background', 2, None, '___sec81'),
+ ('Jackknife code example', 2, None, '___sec80'),
+ ('Resampling methods: Bootstrap', 2, None, '___sec81'),
+ ('Resampling methods: Bootstrap background', 2, None, '___sec82'),
('Resampling methods: More Bootstrap background',
2,
None,
- '___sec82'),
- ('Resampling methods: Bootstrap approach', 2, None, '___sec83'),
- ('Resampling methods: Bootstrap steps', 2, None, '___sec84'),
- ('Resampling methods: Blocking', 2, None, '___sec85'),
- ('Blocking Transformations', 2, None, '___sec86'),
- ('Blocking Transformations', 2, None, '___sec87'),
- ('Blocking Transformations, getting there', 2, None, '___sec88'),
+ '___sec83'),
+ ('Resampling methods: Bootstrap approach', 2, None, '___sec84'),
+ ('Resampling methods: Bootstrap steps', 2, None, '___sec85'),
+ ('Code example for the Bootstrap method', 2, None, '___sec86'),
+ ('Resampling methods: Blocking', 2, None, '___sec87'),
+ ('Blocking Transformations', 2, None, '___sec88'),
+ ('Blocking Transformations', 2, None, '___sec89'),
+ ('Blocking Transformations, getting there', 2, None, '___sec90'),
('Blocking Transformations, final expressions',
2,
None,
- '___sec89'),
+ '___sec91'),
('"Code examples for Blocking, Jackknife and '
'bootstrap":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/ResamplingAnalysisScripts"',
2,
None,
- '___sec90'),
- ('The bias-variance tradeoff', 2, None, '___sec91'),
- ('Training and testing data', 2, None, '___sec92'),
- ('Procedure to find a predictor', 2, None, '___sec93'),
- ('What we want', 2, None, '___sec94'),
- ('The expected generalization error', 2, None, '___sec95'),
- ('Elaborating a little bit more', 2, None, '___sec96'),
- ('The bias', 2, None, '___sec97'),
- ('The variance', 2, None, '___sec98'),
- ('Summing up', 2, None, '___sec99'),
- ('Logistic Regression', 2, None, '___sec100'),
- ('Basics', 2, None, '___sec101'),
- ('Linear classifier', 2, None, '___sec102'),
- ('Some selected properties', 2, None, '___sec103'),
- ('The cross-entropy as a cost function for logistic regression',
+ '___sec92'),
+ ('The bias-variance tradeoff', 2, None, '___sec93'),
+ ('Training and testing data', 2, None, '___sec94'),
+ ('Procedure to find a predictor', 2, None, '___sec95'),
+ ('What we want', 2, None, '___sec96'),
+ ('The expected generalization error', 2, None, '___sec97'),
+ ('Elaborating a little bit more', 2, None, '___sec98'),
+ ('The bias', 2, None, '___sec99'),
+ ('The variance', 2, None, '___sec100'),
+ ('Summing up', 2, None, '___sec101'),
+ ('The one-dimensional Ising model, project 2',
+ 2,
+ None,
+ '___sec102'),
+ ('Example: The one-dimensional Ising model',
+ 2,
+ None,
+ '___sec103'),
+ ('Reformulating the problem to suit regression',
2,
None,
'___sec104'),
- ('Maximum likelihood', 2, None, '___sec105'),
- ('Minimizing the cross entropy', 2, None, '___sec106')]}
+ ('Linear regression', 2, None, '___sec105'),
+ ('Ordinary least squares', 2, None, '___sec106'),
+ ('Singular Value decomposition', 2, None, '___sec107'),
+ ('Fitting with scikit-learn', 2, None, '___sec108'),
+ ('Ridge regression', 2, None, '___sec109'),
+ ('LASSO regression', 2, None, '___sec110'),
+ ('Performance of the different models', 2, None, '___sec111'),
+ ('Performance as function of the regularization parameter',
+ 2,
+ None,
+ '___sec112'),
+ ('Finding the optimal value of $\\lambda$',
+ 2,
+ None,
+ '___sec113')]}
end of tocinfo -->
@@ -351,33 +370,40 @@ MathJax.Hub.Config({
Resampling methods: Jackknife and Bootstrap Resampling methods: Jackknife Resampling methods: Jackknife estimator Resampling methods: Bootstrap Resampling methods: Bootstrap background Resampling methods: More Bootstrap background Resampling methods: Bootstrap approach Resampling methods: Bootstrap steps Resampling methods: Blocking Blocking Transformations Blocking Transformations Blocking Transformations, getting there Blocking Transformations, final expressions "Code examples for Blocking, Jackknife and bootstrap":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/ResamplingAnalysisScripts" The bias-variance tradeoff Training and testing data Procedure to find a predictor What we want The expected generalization error Elaborating a little bit more The bias The variance Summing up Logistic Regression Basics Linear classifier Some selected properties The cross-entropy as a cost function for logistic regression Maximum likelihood Minimizing the cross entropy 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 Resampling methods: Blocking Blocking Transformations Blocking Transformations Blocking Transformations, getting there Blocking Transformations, final expressions "Code examples for Blocking, Jackknife and bootstrap":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/ResamplingAnalysisScripts" The bias-variance tradeoff Training and testing data Procedure to find a predictor What we want The expected generalization error Elaborating a little bit more The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) Linear classifier
+The one-dimensional Ising model, project 2
-Let us start by considering a slightly simpler classifier: a linear classifier that categorizes examples using a weighted linear-combination of the features and an additive offset
+The one-dimensional Ising model with nearest neighbor interaction, no external field and a constant coupling constant \( J \) is given by
+
$$
-\begin{equation}
-s_i = \boldsymbol{x}_i^T\boldsymbol{w} + b_0 \equiv \mathbf{x}_i^T\mathbf{w},
+\begin{align}
+ H = -J \sum_{k}^L s_k s_{k + 1},
\tag{30}
-\end{equation}
+\end{align}
$$
-where we use the short-hand notation
-\( \mathbf{x}_i = (1,\boldsymbol{x}_i) \) and \( \mathbf{w}_i = (b_0,\boldsymbol{w}_i) \).
+where \( s_i \in \{-1, 1\} \) and \( s_{N + 1} = s_1 \). The number of spins in the system is determined by \( L \). For the low temperature limit there is no phase transition.
+
+We will look at a system of \( L = 40 \) spins with a coupling constant of \( J = 1 \). To get enough training data we will generate 10000 states with their respective energies.
+
+
+
+
+
= 40
+n = int (1e4 )
+
+spins = np. random. choice([-1 , 1 ], size= (n, L))
+J = 1.0
+
+energies = np. zeros(n)
+
+for i in range (n):
+ energies[i] = - J * np. dot(spins[i], np. roll(spins[i], 1 ))
+
@@ -428,6 +471,12 @@ where we use the short-hand notation
107 108 109 110 111 112 113 ... 115 » Resampling methods: Jackknife and Bootstrap Resampling methods: Jackknife Resampling methods: Jackknife estimator Resampling methods: Bootstrap Resampling methods: Bootstrap background Resampling methods: More Bootstrap background Resampling methods: Bootstrap approach Resampling methods: Bootstrap steps Resampling methods: Blocking Blocking Transformations Blocking Transformations Blocking Transformations, getting there Blocking Transformations, final expressions "Code examples for Blocking, Jackknife and bootstrap":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/ResamplingAnalysisScripts" The bias-variance tradeoff Training and testing data Procedure to find a predictor What we want The expected generalization error Elaborating a little bit more The bias The variance Summing up Logistic Regression Basics Linear classifier Some selected properties The cross-entropy as a cost function for logistic regression Maximum likelihood Minimizing the cross entropy 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 Resampling methods: Blocking Blocking Transformations Blocking Transformations Blocking Transformations, getting there Blocking Transformations, final expressions "Code examples for Blocking, Jackknife and bootstrap":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/ResamplingAnalysisScripts" The bias-variance tradeoff Training and testing data Procedure to find a predictor What we want The expected generalization error Elaborating a little bit more The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) Some selected properties
+Example: The one-dimensional Ising model
-This function takes values on the entire real axis. In the case of
-logistic regression, however, the labels \( y_i \) are discrete
-variables. One simple way to get a discrete output is to have sign
-functions that map the output of a linear regressor to \( \{0,1\} \),
-\( f(s_i)=sign(s_i)=1 \) if \( s_i\ge 0 \) and 0 if otherwise. Indeed,
-this is commonly known as the "perceptron" in the machine learning
-literature. This model is extremely simple, and it is favorable in
-many cases (e.g. noisy data) to have a ``soft" classifier that outputs
-the probability of a given category. For example, given
-\( \mathbf{x}_i \), the classifier outputs the probability of being in
-category \( m \). One such function is the logistic (or sigmoid) function:
-
-$$
-\begin{equation}
-f(s) = \frac{1}{1+\mathrm e^{-s}}.
-\tag{31}
-\end{equation}
-$$
-
-Note that \( 1-f(s)= f(-s) \), which will be useful shortly.
+Here we use linear (ordinary least squares), ridge and LASSO
+regression to predict the energy in the nearest neighbor
+one-dimensional Ising model on a ring, i.e., the endpoints wrap
+around. We will use the linear regression models to fit a value for
+the coupling constant to achieve this.
+
+
+
import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+import seaborn as sns
+import scipy.linalg as scl
+from sklearn.model_selection import train_test_split
+import sklearn.linear_model as skl
+import tqdm
+sns. set(color_codes= True )
+cmap_args= dict (vmin=-1. , vmax=1. , cmap= 'seismic' )
+
@@ -437,6 +461,13 @@ Note that \( 1-f(s)= f(-s) \), which will be useful shortly.
107 108 109 110 111 112 113 114 ... 115 » Resampling methods: Jackknife and Bootstrap Resampling methods: Jackknife Resampling methods: Jackknife estimator Resampling methods: Bootstrap Resampling methods: Bootstrap background Resampling methods: More Bootstrap background Resampling methods: Bootstrap approach Resampling methods: Bootstrap steps Resampling methods: Blocking Blocking Transformations Blocking Transformations Blocking Transformations, getting there Blocking Transformations, final expressions "Code examples for Blocking, Jackknife and bootstrap":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/ResamplingAnalysisScripts" The bias-variance tradeoff Training and testing data Procedure to find a predictor What we want The expected generalization error Elaborating a little bit more The bias The variance Summing up Logistic Regression Basics Linear classifier Some selected properties The cross-entropy as a cost function for logistic regression Maximum likelihood Minimizing the cross entropy 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 Resampling methods: Blocking Blocking Transformations Blocking Transformations Blocking Transformations, getting there Blocking Transformations, final expressions "Code examples for Blocking, Jackknife and bootstrap":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/ResamplingAnalysisScripts" The bias-variance tradeoff Training and testing data Procedure to find a predictor What we want The expected generalization error Elaborating a little bit more The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) The cross-entropy as a cost function for logistic regression
+Reformulating the problem to suit regression
-The perceptron is an example of a ``hard classification": each datapoint is deterministically assigned to a category (i.e \( y_i=0 \) or \( y_i=1 \)). In many cases, it is favorable to have a "soft" classifier that outputs the probability of a given category rather than a single value. For example, given \( \mathbf{x}_i \), the classifier outputs the probability of being in category \( m \).
-Logistic regression is the most canonical example of a soft classifier. In logistic regression, the probability that a data point \( \boldsymbol{x}_i \) belongs to a category \( y_i=\{0,1\} \) is is given by
-$$
-\begin{eqnarray}
-P(y_i=1|\boldsymbol{x}_i,\boldsymbol{\theta)} &=& \frac{1}{1+\mathrm{e}^{-\mathbf{x}^T_i\mathbf{w}}},\nonumber\\
-P(y_i=0|\boldsymbol{x}_i,\boldsymbol{\theta)} &=& 1 - P(y_i=1|\boldsymbol{x}_i,\boldsymbol{\theta)},
-\end{eqnarray}
-$$
+A more general form for the one-dimensional Ising model is
-where \( \boldsymbol{\theta}=\mathbf{w} \) are the weights we wish to learn from the data.
+$$
+\begin{align}
+ H = - \sum_j^L \sum_k^L s_j s_k J_{jk}.
+\tag{31}
+\end{align}
+$$
-Notice that in terms of the logistic function, we can write
+Here we allow for interactions beyond the nearest neighbors and a more
+adaptive coupling matrix. This latter expression can be formulated as
+a matrix-product on the form
$$
-P(y_i=1) =f(\mathbf{x}_i^T\mathbf{w})=1-P(y_i=0).
+\begin{align}
+ H = X J,
+\tag{32}
+\end{align}
$$
+
+where \( X_{jk} = s_j s_k \) and \( J \) is the matrix consisting of the
+elements \( -J_{jk} \). This form of writing the energy fits perfectly
+with the form utilized in linear regression, viz.
+$$
+\begin{align}
+ y = X\omega + \epsilon,
+\tag{33}
+\end{align}
+$$
+
+
+
+
+
= np. zeros((n, L ** 2 ))
+for i in range (n):
+ X[i] = np. outer(spins[i], spins[i]). ravel()
+y = energies
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.96 )
+
+X_train_own = np. concatenate(
+ (np. ones(len (X_train))[:, np. newaxis], X_train),
+ axis=1
+)
+
+X_test_own = np. concatenate(
+ (np. ones(len (X_test))[:, np. newaxis], X_test),
+ axis=1
+)
+
@@ -432,6 +491,12 @@ $$
107 108 109 110 111 112 113 114 115 » Resampling methods: Jackknife and Bootstrap Resampling methods: Jackknife Resampling methods: Jackknife estimator Resampling methods: Bootstrap Resampling methods: Bootstrap background Resampling methods: More Bootstrap background Resampling methods: Bootstrap approach Resampling methods: Bootstrap steps Resampling methods: Blocking Blocking Transformations Blocking Transformations Blocking Transformations, getting there Blocking Transformations, final expressions "Code examples for Blocking, Jackknife and bootstrap":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/ResamplingAnalysisScripts" The bias-variance tradeoff Training and testing data Procedure to find a predictor What we want The expected generalization error Elaborating a little bit more The bias The variance Summing up Logistic Regression Basics Linear classifier Some selected properties The cross-entropy as a cost function for logistic regression Maximum likelihood Minimizing the cross entropy 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 Resampling methods: Blocking Blocking Transformations Blocking Transformations Blocking Transformations, getting there Blocking Transformations, final expressions "Code examples for Blocking, Jackknife and bootstrap":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/ResamplingAnalysisScripts" The bias-variance tradeoff Training and testing data Procedure to find a predictor What we want The expected generalization error Elaborating a little bit more The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \)
Maximum likelihood
+Linear regression
-We now define the cost function for logistic regression using Maximum
-Likelihood Estimation (MLE). Recall, that in MLE we choose parameters
-to maximize the probability of seeing the observed data. Consider a
-dataset \( \mathcal{D}=\{(y_i,\boldsymbol{x}_i)\} \) with binary labels
-\( y_i\in\{0,1\} \) where the data points are drawn independently. The
-likelihood of the seeing the data under our model is just:
+The problem at hand is to try to fit the equation
$$
\begin{align}
-P(\mathcal{D}|\mathbf{w})& = \prod_{i=1}^n \left[f(\mathbf{x}_i^T\mathbf{w})\right]^{y_i}\left[1-f(\mathbf{x}_i^T\mathbf{w})\right]^{1-y_i}\nonumber \\
-\tag{32}
+ y = f(x) + \epsilon,
+\tag{34}
\end{align}
$$
-from which we can readily compute the log-likelihood:
+
+where \( f(x) \) is some unknown function of the data \( x \) and \( \epsilon \)
+is normally distributed with mean zero noise with standard deviation
+\( \sigma_{\epsilon} \). Our job is to try to find a predictor which
+estimates the function \( f(x) \). In linear regression we assume that we
+can formulate the problem as
+
$$
-\begin{equation}
-l(\mathbf{w}) = \sum_{i=1}^n y_i\log f(\mathbf{x}_i^T\mathbf{w}) + (1-y_i)\log\left[1-f(\mathbf{x}_i^T\mathbf{w})\right].
-\tag{33}
-\end{equation}
+\begin{align}
+ y = X\omega + \epsilon,
+\tag{35}
+\end{align}
$$
+
+where \( X \) and \( \omega \) are now matrices. Our job at hand is now to
+find a cost function \( C \), which we wish to minimize in order to find
+the best estimate of \( \omega \).
+
@@ -435,6 +467,12 @@ $$
107 108 109 110 111 112 113 114 115 » Resampling methods: Jackknife and Bootstrap Resampling methods: Jackknife Resampling methods: Jackknife estimator Resampling methods: Bootstrap Resampling methods: Bootstrap background Resampling methods: More Bootstrap background Resampling methods: Bootstrap approach Resampling methods: Bootstrap steps Resampling methods: Blocking Blocking Transformations Blocking Transformations Blocking Transformations, getting there Blocking Transformations, final expressions "Code examples for Blocking, Jackknife and bootstrap":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/ResamplingAnalysisScripts" The bias-variance tradeoff Training and testing data Procedure to find a predictor What we want The expected generalization error Elaborating a little bit more The bias The variance Summing up Logistic Regression Basics Linear classifier Some selected properties The cross-entropy as a cost function for logistic regression Maximum likelihood Minimizing the cross entropy 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 Resampling methods: Blocking Blocking Transformations Blocking Transformations Blocking Transformations, getting there Blocking Transformations, final expressions "Code examples for Blocking, Jackknife and bootstrap":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/ResamplingAnalysisScripts" The bias-variance tradeoff Training and testing data Procedure to find a predictor What we want The expected generalization error Elaborating a little bit more The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) Ordinary least squares
+
+
+In the ordinary least squares method we choose the cost function
$$
-\hat{\mathbf{w}} = \sum_{i=1}^n y_i\log f(\mathbf{x}_i^T\mathbf{w}) + (1-y_i)\log\left[1-f(\mathbf{x}_i^T\mathbf{w})\right].
+\begin{align}
+ C(X, \omega) = ||X\omega - y||^2
+ = (X\omega - y)^T(X\omega - y)
+\tag{36}
+\end{align}
$$
-Since the cost (error) function is just the negative log-likelihood, for logistic regression we have that
+We then find the extremal point of \( C \) by taking the derivative with respect to \( \omega \) and setting it to zero, i.e.,
+
$$
-\begin{eqnarray}
-\mathcal{C}(\mathbf{w}) &=& - l(\mathbf{w}) \\
-&=& \sum_{i=1}^n -y_i\log f(\mathbf{x}_i^T\mathbf{w}) - (1-y_i)\log\left[1-f(\mathbf{x}_i^T\mathbf{w})\right].\nonumber
-\end{eqnarray}
+\begin{align}
+ \dfrac{\mathrm{d}C}{\mathrm{d}\omega}
+ = 0.
+\tag{37}
+\end{align}
$$
-This equation is known in statistics as the \emph{cross entropy}. Finally, we note that just as in linear regression,
-in practice we usually supplement the cross-entropy with additional regularization terms, usually \( L_1 \) and \( L_2 \) regularization as we did for Ridge and Lasso regression.
+This yields the expression for \( \omega \) to be
+$$
+\begin{align}
+ \omega = \frac{X^T y}{X^T X},
+\tag{38}
+\end{align}
+$$
+
+which immediately imposes some requirements on \( X \) as there must exist
+an inverse of \( X^T X \). If the expression we are modelling contains an
+intercept, i.e., a constant expression we must make sure that the
+first column of \( X \) consists of \( 1 \).
+
+
+
+
+
def get_ols_weights_naive (x: np. ndarray, y: np. ndarray) -> np. ndarray:
+ return scl. inv(x. T @ x) @ (x. T @ y)
+omega = get_ols_weights_naive(X_train_own, y_train)
+
@@ -425,6 +479,12 @@ in practice we usually supplement the cross-entropy with additional regularizati
107 108 109 110 111 112 113 114 115 » Resampling methods: Jackknife and Bootstrap Resampling methods: Jackknife Resampling methods: Jackknife estimator Resampling methods: Bootstrap Resampling methods: Bootstrap background Resampling methods: More Bootstrap background Resampling methods: Bootstrap approach Resampling methods: Bootstrap steps Resampling methods: Blocking Blocking Transformations Blocking Transformations Blocking Transformations, getting there Blocking Transformations, final expressions "Code examples for Blocking, Jackknife and bootstrap":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/ResamplingAnalysisScripts" The bias-variance tradeoff Training and testing data Procedure to find a predictor What we want The expected generalization error Elaborating a little bit more The bias The variance Summing up Logistic Regression Basics Linear classifier Some selected properties The cross-entropy as a cost function for logistic regression Maximum likelihood Minimizing the cross entropy 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 Resampling methods: Blocking Blocking Transformations Blocking Transformations Blocking Transformations, getting there Blocking Transformations, final expressions "Code examples for Blocking, Jackknife and bootstrap":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/ResamplingAnalysisScripts" The bias-variance tradeoff Training and testing data Procedure to find a predictor What we want The expected generalization error Elaborating a little bit more The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) Minimizing the cross entropy
-
-
-The cross entropy is a convex function of the weights \( \mathbf{w} \) and,
-therefore, any local minimizer is a global minimizer. Minimizing this
-cost function leads to the following equation
+
Singular Value decomposition
+Doing the inversion directly turns out to be a bad idea as the matrix
+\( X^TX \) is singular. An alternative approach is to use the singular
+value decomposition . Using the definition of the Moore-Penrose
+pseudoinverse we can write the equation for \( \omega \) as
$$
-\begin{equation}
-\boldsymbol{0}=\boldsymbol{\nabla} \mathcal{C}(\mathbf{w}) = \sum_{i=1}^n\left[f(\mathbf{x}_i^T\mathbf{w})-y_i\right]\mathbf{x}_i,
-\tag{34}
-\end{equation}
+\begin{align}
+ \omega = X^{+}y,
+\tag{39}
+\end{align}
+$$
+
+where the pseudoinverse of \( X \) is given by
+$$
+\begin{align}
+ X^{+} = \frac{X^T}{X^T X}.
+\tag{40}
+\end{align}
$$
-where we made use of the logistic function identity \( \partial_z f(z) =
-f(z)[1-f(z)] \). This equation defines a transcendental equation for
-\( \mathbf{w} \), the solution of which, unlike linear regression, cannot
-be written in a closed form.
-Here we need gradient descent methods!
+Using singular value decomposition we have that \( X = U\Sigma V^T \),
+where \( X^{+} = V\Sigma^{+} U^T \). This reduces the equation for
+\( \omega \) to
+$$
+\begin{align}
+ \omega = V\Sigma^{+} U^T y.
+\tag{41}
+\end{align}
+$$
+
+Note that solving this equation by actually doing the pseudoinverse
+(which is what we will do) is not a good idea as this operation scales
+as \( \mathcal{O}(n^3) \), where \( n \) is the number of elements in a
+general matrix. Instead, doing \( QR \)-factorization and solving the
+linear system as an equation would reduce this down to
+\( \mathcal{O}(n^2) \) operations.
+
+
+
+
+
def get_ols_weights (x: np. ndarray, y: np. ndarray) -> np. ndarray:
+ u, s, v = scl. svd(x)
+ return v. T @ scl. pinv(scl. diagsvd(s, u. shape[0 ], v. shape[0 ])) @ u. T @ y
+
+Before passing in the data to the function we append a column with ones to the training data.
+
+
+
+
+
= get_ols_weights(X_train_own,y_train)
+
diff --git a/doc/pub/Regression/html/._Regression-bs109.html b/doc/pub/Regression/html/._Regression-bs109.html
new file mode 100644
index 000000000..804d9cf08
--- /dev/null
+++ b/doc/pub/Regression/html/._Regression-bs109.html
@@ -0,0 +1,515 @@
+
+
+
+Data Analysis and Machine Learning: Linear Regression and more Advanced Regression Analysis
+
+
+
+
+
+
+
+
+
+
Fitting with scikit-learn
+
+
+Next we fit a LinearRegression-model from Scikit-learn for comparison.
+
+
+
+
+
= skl. LinearRegression(). fit(X_train, y_train)
+
+
+Extracting the \( J \)-matrix from both our own method and the Scikit-learn model where we make sure to remove the intercept.
+
+
+
+
+
= omega[1 :]. reshape(L, L)
+J_sk = clf. coef_. reshape(L, L)
+
+
+A way of looking at the coefficients in \( J \) is to plot the matrices as images.
+
+
+
+
+
= plt. figure(figsize= (20 , 14 ))
+im = plt. imshow(J_own, ** cmap_args)
+plt. title("Home-made OLS" , fontsize=18 )
+plt. xticks(fontsize=18 )
+plt. yticks(fontsize=18 )
+cb = fig. colorbar(im)
+cb. ax. set_yticklabels(cb. ax. get_yticklabels(), fontsize=18 )
+
+fig = plt. figure(figsize= (20 , 14 ))
+im = plt. imshow(J_sk, ** cmap_args)
+plt. title("LinearRegression from Scikit-learn" , fontsize=18 )
+plt. xticks(fontsize=18 )
+plt. yticks(fontsize=18 )
+cb = fig. colorbar(im)
+cb. ax. set_yticklabels(cb. ax. get_yticklabels(), fontsize=18 )
+plt. show()
+
+
+We can see that our model for the least squares method performes close
+to the benchmark from Scikit-learn. It is interesting to note that OLS
+considers both \( J_{j, j + 1} = -0.5 \) and \( J_{j, j - 1} = -0.5 \) as
+valid matrix elements for \( J \).
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/doc/pub/Regression/html/._Regression-bs110.html b/doc/pub/Regression/html/._Regression-bs110.html
new file mode 100644
index 000000000..1994449a9
--- /dev/null
+++ b/doc/pub/Regression/html/._Regression-bs110.html
@@ -0,0 +1,528 @@
+
+
+
+Data Analysis and Machine Learning: Linear Regression and more Advanced Regression Analysis
+
+
+
+
+
+
+
+
+
+
Ridge regression
+
+
+Having explored the ordinary least squares we move on to ridge
+regression. In ridge regression we include a regularizer . This
+involves a new cost function which leads to a new estimate for the
+weights \( \omega \). This results in a penalized regression problem. The
+cost function is given by
+
+$$
+\begin{align}
+ C(X, \omega; \lambda) = ||X\omega - y||^2 + \lambda ||\omega||^2
+ = (X\omega - y)^T(X\omega - y) + \lambda \omega^T\omega.
+\tag{42}
+\end{align}
+$$
+
+Finding the extremum of this function yields the weights
+
+$$
+\begin{align}
+ \omega(\lambda) = \frac{X^Ty}{X^TX + \lambda} \to \frac{\omega_{\text{LS}}}{1 + \lambda},
+\tag{43}
+\end{align}
+$$
+
+
+where \( \omega_{\text{LS}} \) is the weights from ordinary least
+squares. The last assumption assumes that \( X \) is orthogonal, which it
+is not. We will therefore resort to solving the equation as it stands
+on the left hand side.
+
+
+
+
+
def get_ridge_weights (x: np. ndarray, y: np. ndarray, _lambda: float ) -> np. ndarray:
+ return x. T @ y @ scl. inv(
+ x. T @ x + np. eye(x. shape[1 ], x. shape[1 ]) * _lambda
+ )
+lambda = 0.1
+omega_ridge = get_ridge_weights(X_train_own, y_train, np. array([_lambda]))
+clf_ridge = skl. Ridge(alpha= _lambda). fit(X_train, y_train)
+J_ridge_own = omega_ridge[1 :]. reshape(L, L)
+J_ridge_sk = clf_ridge. coef_. reshape(L, L)
+fig = plt. figure(figsize= (20 , 14 ))
+im = plt. imshow(J_ridge_own, ** cmap_args)
+plt. title("Home-made ridge regression" , fontsize=18 )
+plt. xticks(fontsize=18 )
+plt. yticks(fontsize=18 )
+cb = fig. colorbar(im)
+cb. ax. set_yticklabels(cb. ax. get_yticklabels(), fontsize=18 )
+
+fig = plt. figure(figsize= (20 , 14 ))
+im = plt. imshow(J_ridge_sk, ** cmap_args)
+plt. title("Ridge from Scikit-learn" , fontsize=18 )
+plt. xticks(fontsize=18 )
+plt. yticks(fontsize=18 )
+cb = fig. colorbar(im)
+cb. ax. set_yticklabels(cb. ax. get_yticklabels(), fontsize=18 )
+
+plt. show()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/doc/pub/Regression/html/._Regression-bs111.html b/doc/pub/Regression/html/._Regression-bs111.html
new file mode 100644
index 000000000..0e1b20c58
--- /dev/null
+++ b/doc/pub/Regression/html/._Regression-bs111.html
@@ -0,0 +1,500 @@
+
+
+
+Data Analysis and Machine Learning: Linear Regression and more Advanced Regression Analysis
+
+
+
+
+
+
+
+
+
+
LASSO regression
+
+
+In the Least Absolute Shrinkage and Selection Operator (LASSO)-method we get a third cost function.
+$$
+\begin{align}
+ C(X, \omega; \lambda) =
+ ||X\omega - y||^2 + \lambda ||\omega||
+ = (X\omega - y)^T(X\omega - y) + \lambda \sqrt{\omega^T\omega}.
+\tag{44}
+\end{align}
+$$
+
+Finding the extremal point of this cost function is not so straight-forward as in least squares and ridge. We will therefore rely solely on the function ``Lasso`` from Scikit-learn.
+
+
+
+
+
= skl. Lasso(alpha= _lambda). fit(X_train, y_train)
+J_lasso_sk = clf_lasso. coef_. reshape(L, L)
+fig = plt. figure(figsize= (20 , 14 ))
+im = plt. imshow(J_lasso_sk, ** cmap_args)
+plt. title("Lasso from Scikit-learn" , fontsize=18 )
+plt. xticks(fontsize=18 )
+plt. yticks(fontsize=18 )
+cb = fig. colorbar(im)
+cb. ax. set_yticklabels(cb. ax. get_yticklabels(), fontsize=18 )
+
+plt. show()
+
+
+It is quite striking how LASSO breaks the symmetry of the coupling
+constant as opposed to ridge and OLS. We get a sparse solution with
+\( J_{j, j + 1} = -1 \).
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/doc/pub/Regression/html/._Regression-bs112.html b/doc/pub/Regression/html/._Regression-bs112.html
new file mode 100644
index 000000000..4c095c632
--- /dev/null
+++ b/doc/pub/Regression/html/._Regression-bs112.html
@@ -0,0 +1,495 @@
+
+
+
+Data Analysis and Machine Learning: Linear Regression and more Advanced Regression Analysis
+
+
+
+
+
+
+
+
+
+
Performance of the different models
+
+
+In order to judge which model performs best at varying values of \( \lambda \) (for ridge and LASSO) we compute \( R^2 \) which is given by
+
+$$
+\begin{align}
+ R^2 = 1 - \frac{(y - \hat{y})^2}{(y - \bar{y})^2},
+\tag{45}
+\end{align}
+$$
+
+where \( y \) is a vector with the true values of the energy, \( \hat{y} \) is the predicted values of \( y \) from the models and \( \bar{y} \) is the mean of \( \hat{y} \).
+
+
+
+
+
def r_squared (y, y_hat):
+ return 1 - np. sum((y - y_hat) ** 2 ) / np. sum((y - np. mean(y_hat)) ** 2 )
+
+
+This is the same metric used by Scikit-learn for their regression models when scoring.
+
+
+
+
= clf. predict(X_test)
+r_test = r_squared(y_test, y_hat)
+sk_r_test = clf. score(X_test, y_test)
+
+assert abs (r_test - sk_r_test) < 1e-2
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/doc/pub/Regression/html/._Regression-bs113.html b/doc/pub/Regression/html/._Regression-bs113.html
new file mode 100644
index 000000000..f33282cad
--- /dev/null
+++ b/doc/pub/Regression/html/._Regression-bs113.html
@@ -0,0 +1,542 @@
+
+
+
+Data Analysis and Machine Learning: Linear Regression and more Advanced Regression Analysis
+
+
+
+
+
+
+
+
+
+
Performance as function of the regularization parameter
+
+
+We see how the different models perform for a different set of values for \( \lambda \).
+
+
+
+
+
= np. logspace(-4 , 5 , 10 )
+
+train_errors = {
+ "ols_own" : np. zeros(lambdas. size),
+ "ols_sk" : np. zeros(lambdas. size),
+ "ridge_own" : np. zeros(lambdas. size),
+ "ridge_sk" : np. zeros(lambdas. size),
+ "lasso_sk" : np. zeros(lambdas. size)
+}
+
+test_errors = {
+ "ols_own" : np. zeros(lambdas. size),
+ "ols_sk" : np. zeros(lambdas. size),
+ "ridge_own" : np. zeros(lambdas. size),
+ "ridge_sk" : np. zeros(lambdas. size),
+ "lasso_sk" : np. zeros(lambdas. size)
+}
+
+plot_counter = 1
+
+fig = plt. figure(figsize= (32 , 54 ))
+
+for i, _lambda in enumerate (tqdm. tqdm(lambdas)):
+ omega = get_ols_weights(X_train_own, y_train)
+ y_hat_train = X_train_own @ omega
+ y_hat_test = X_test_own @ omega
+
+ train_errors["ols_own" ][i] = r_squared(y_train, y_hat_train)
+ test_errors["ols_own" ][i] = r_squared(y_test, y_hat_test)
+
+ plt. subplot(10 , 5 , plot_counter)
+ plt. imshow(omega[1 :]. reshape(L, L), ** cmap_args)
+ plt. title("Home made OLS" )
+ plot_counter += 1
+
+ omega = get_ridge_weights(X_train_own, y_train, _lambda)
+ y_hat_train = X_train_own @ omega
+ y_hat_test = X_test_own @ omega
+
+ train_errors["ridge_own" ][i] = r_squared(y_train, y_hat_train)
+ test_errors["ridge_own" ][i] = r_squared(y_test, y_hat_test)
+
+ plt. subplot(10 , 5 , plot_counter)
+ plt. imshow(omega[1 :]. reshape(L, L), ** cmap_args)
+ plt. title(r"Home made ridge, $\lambda = %.4f $" % _lambda)
+ plot_counter += 1
+
+ for key, method in zip (
+ ["ols_sk" , "ridge_sk" , "lasso_sk" ],
+ [skl. LinearRegression(), skl. Ridge(alpha= _lambda), skl. Lasso(alpha= _lambda)]
+ ):
+ method = method. fit(X_train, y_train)
+
+ train_errors[key][i] = method. score(X_train, y_train)
+ test_errors[key][i] = method. score(X_test, y_test)
+
+ omega = method. coef_. reshape(L, L)
+
+ plt. subplot(10 , 5 , plot_counter)
+ plt. imshow(omega, ** cmap_args)
+ plt. title(r" %s , $\lambda = %.4f $" % (key, _lambda))
+ plot_counter += 1
+
+plt. show()
+
+
+We can see that LASSO quite fast reaches a good solution for low
+values of \( \lambda \), but will "wither" when we increase \( \lambda \) too
+much. Ridge is more stable over a larger range of values for
+\( \lambda \), but eventually also fades away.
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/doc/pub/Regression/html/._Regression-bs114.html b/doc/pub/Regression/html/._Regression-bs114.html
new file mode 100644
index 000000000..9294abbed
--- /dev/null
+++ b/doc/pub/Regression/html/._Regression-bs114.html
@@ -0,0 +1,514 @@
+
+
+
+Data Analysis and Machine Learning: Linear Regression and more Advanced Regression Analysis
+
+
+
+
+
+
+
+
+
+
Finding the optimal value of \( \lambda \)
+
+
+To determine which value of \( \lambda \) is best we plot the accuracy of
+the models when predicting the training and the testing set. We expect
+the accuracy of the training set to be quite good, but if the accuracy
+of the testing set is much lower this tells us that we might be
+subject to an overfit model. The ideal scenario is an accuracy on the
+testing set that is close to the accuracy of the training set.
+
+
+
+
+
= plt. figure(figsize= (20 , 14 ))
+
+colors = {
+ "ols_own" : "b" ,
+ "ridge_own" : "g" ,
+ "ols_sk" : "r" ,
+ "ridge_sk" : "y" ,
+ "lasso_sk" : "c"
+}
+
+for key in train_errors:
+ plt. semilogx(
+ lambdas,
+ train_errors[key],
+ colors[key],
+ label= "Train {0}" . format(key),
+ linewidth=4.0
+ )
+
+for key in test_errors:
+ plt. semilogx(
+ lambdas,
+ test_errors[key],
+ colors[key] + "--" ,
+ label= "Test {0}" . format(key),
+ linewidth=4.0
+ )
+#plt.semilogx(lambdas, train_errors["ols_own"], label="Train (OLS own)")
+#plt.semilogx(lambdas, test_errors["ols_own"], label="Test (OLS own)")
+
+plt. legend(loc= "best" , fontsize=18 )
+plt. xlabel(r"$\lambda$" , fontsize=18 )
+plt. ylabel(r"$R^2$" , fontsize=18 )
+plt. tick_params(labelsize=18 )
+plt. show()
+
+
+From the above figure we can see that LASSO with \( \lambda = 10^{-2} \)
+achieve a very good accuracy on the test set. This by far surpases the
+other models for all values of \( \lambda \).
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/doc/pub/Regression/html/Regression-bs.html b/doc/pub/Regression/html/Regression-bs.html
index 833e9dec8..e6ce300e7 100644
--- a/doc/pub/Regression/html/Regression-bs.html
+++ b/doc/pub/Regression/html/Regression-bs.html
@@ -225,7 +225,34 @@ Automatically generated HTML file from DocOnce source
('Elaborating a little bit more', 2, None, '___sec98'),
('The bias', 2, None, '___sec99'),
('The variance', 2, None, '___sec100'),
- ('Summing up', 2, None, '___sec101')]}
+ ('Summing up', 2, None, '___sec101'),
+ ('The one-dimensional Ising model, project 2',
+ 2,
+ None,
+ '___sec102'),
+ ('Example: The one-dimensional Ising model',
+ 2,
+ None,
+ '___sec103'),
+ ('Reformulating the problem to suit regression',
+ 2,
+ None,
+ '___sec104'),
+ ('Linear regression', 2, None, '___sec105'),
+ ('Ordinary least squares', 2, None, '___sec106'),
+ ('Singular Value decomposition', 2, None, '___sec107'),
+ ('Fitting with scikit-learn', 2, None, '___sec108'),
+ ('Ridge regression', 2, None, '___sec109'),
+ ('LASSO regression', 2, None, '___sec110'),
+ ('Performance of the different models', 2, None, '___sec111'),
+ ('Performance as function of the regularization parameter',
+ 2,
+ None,
+ '___sec112'),
+ ('Finding the optimal value of $\\lambda$',
+ 2,
+ None,
+ '___sec113')]}
end of tocinfo -->
@@ -365,6 +392,18 @@ MathJax.Hub.Config({
The bias The variance Summing up The one-dimensional Ising model, project 2 Example: The one-dimensional Ising model Reformulating the problem to suit regression Linear regression Ordinary least squares Singular Value decomposition Fitting with scikit-learn Ridge regression LASSO regression Performance of the different models Performance as function of the regularization parameter Finding the optimal value of \( \lambda \) [2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University
-
Sep 25, 2018 Oct 11, 2018
@@ -423,7 +462,7 @@ MathJax.Hub.Config({
9 10 ... 103 115 » [2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University
Sep 25, 2018 Oct 11, 2018
@@ -3541,6 +3541,624 @@ flexible statistical methods have higher variance.
+
+The one-dimensional Ising model, project 2
+
+
+The one-dimensional Ising model with nearest neighbor interaction, no external field and a constant coupling constant \( J \) is given by
+
+
+We will look at a system of \( L = 40 \) spins with a coupling constant of \( J = 1 \). To get enough training data we will generate 10000 states with their respective energies.
+
+
+
+
+
40
+n = int (1e4 )
+
+spins = np.random.choice([-1 , 1 ], size=(n, L))
+J = 1.0
+
+energies = np.zeros(n)
+
+for i in range (n):
+ energies[i] = - J * np.dot(spins[i], np.roll(spins[i], 1 ))
+
+
+
+
+Example: The one-dimensional Ising model
+
+
+Here we use linear (ordinary least squares), ridge and LASSO
+regression to predict the energy in the nearest neighbor
+one-dimensional Ising model on a ring, i.e., the endpoints wrap
+around. We will use the linear regression models to fit a value for
+the coupling constant to achieve this.
+
+
+
+
import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+import seaborn as sns
+import scipy.linalg as scl
+from sklearn.model_selection import train_test_split
+import sklearn.linear_model as skl
+import tqdm
+sns.set(color_codes=True )
+cmap_args=dict (vmin=-1. , vmax=1. , cmap='seismic' )
+
+
+
+
+Reformulating the problem to suit regression
+
+
+A more general form for the one-dimensional Ising model is
+
+
+Here we allow for interactions beyond the nearest neighbors and a more
+adaptive coupling matrix. This latter expression can be formulated as
+a matrix-product on the form
+
+where \( X_{jk} = s_j s_k \) and \( J \) is the matrix consisting of the
+elements \( -J_{jk} \). This form of writing the energy fits perfectly
+with the form utilized in linear regression, viz.
+
+
+
+
2 ))
+for i in range (n):
+ X[i] = np.outer(spins[i], spins[i]).ravel()
+y = energies
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.96 )
+
+X_train_own = np.concatenate(
+ (np.ones(len (X_train))[:, np.newaxis], X_train),
+ axis=1
+)
+
+X_test_own = np.concatenate(
+ (np.ones(len (X_test))[:, np.newaxis], X_test),
+ axis=1
+)
+
+
+
+
+Linear regression
+
+
+The problem at hand is to try to fit the equation
+
+where \( f(x) \) is some unknown function of the data \( x \) and \( \epsilon \)
+is normally distributed with mean zero noise with standard deviation
+\( \sigma_{\epsilon} \). Our job is to try to find a predictor which
+estimates the function \( f(x) \). In linear regression we assume that we
+can formulate the problem as
+
+
+where \( X \) and \( \omega \) are now matrices. Our job at hand is now to
+find a cost function \( C \), which we wish to minimize in order to find
+the best estimate of \( \omega \).
+
+
+
+
+Ordinary least squares
+
+
+In the ordinary least squares method we choose the cost function
+
+which immediately imposes some requirements on \( X \) as there must exist
+an inverse of \( X^T X \). If the expression we are modelling contains an
+intercept, i.e., a constant expression we must make sure that the
+first column of \( X \) consists of \( 1 \).
+
+
+
+
+
def get_ols_weights_naive (x: np.ndarray, y: np.ndarray) -> np.ndarray:
+ return scl.inv(x.T @ x) @ (x.T @ y)
+omega = get_ols_weights_naive(X_train_own, y_train)
+
+
+
+
+Singular Value decomposition
+Doing the inversion directly turns out to be a bad idea as the matrix
+\( X^TX \) is singular. An alternative approach is to use the singular
+value decomposition . Using the definition of the Moore-Penrose
+pseudoinverse we can write the equation for \( \omega \) as
+
+
+Using singular value decomposition we have that \( X = U\Sigma V^T \),
+where \( X^{+} = V\Sigma^{+} U^T \). This reduces the equation for
+\( \omega \) to
+
+Note that solving this equation by actually doing the pseudoinverse
+(which is what we will do) is not a good idea as this operation scales
+as \( \mathcal{O}(n^3) \), where \( n \) is the number of elements in a
+general matrix. Instead, doing \( QR \)-factorization and solving the
+linear system as an equation would reduce this down to
+\( \mathcal{O}(n^2) \) operations.
+
+
+
+
+
def get_ols_weights (x: np.ndarray, y: np.ndarray) -> np.ndarray:
+ u, s, v = scl.svd(x)
+ return v.T @ scl.pinv(scl.diagsvd(s, u.shape[0 ], v.shape[0 ])) @ u.T @ y
+
+Before passing in the data to the function we append a column with ones to the training data.
+
+
+
+
+
+
+
+
+Fitting with scikit-learn
+
+
+Next we fit a LinearRegression-model from Scikit-learn for comparison.
+
+
+
+
+
+Extracting the \( J \)-matrix from both our own method and the Scikit-learn model where we make sure to remove the intercept.
+
+
+
+
+
1 :].reshape(L, L)
+J_sk = clf.coef_.reshape(L, L)
+
+A way of looking at the coefficients in \( J \) is to plot the matrices as images.
+
+
+
+
+
20 , 14 ))
+im = plt.imshow(J_own, **cmap_args)
+plt.title("Home-made OLS" , fontsize=18 )
+plt.xticks(fontsize=18 )
+plt.yticks(fontsize=18 )
+cb = fig.colorbar(im)
+cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18 )
+
+fig = plt.figure(figsize=(20 , 14 ))
+im = plt.imshow(J_sk, **cmap_args)
+plt.title("LinearRegression from Scikit-learn" , fontsize=18 )
+plt.xticks(fontsize=18 )
+plt.yticks(fontsize=18 )
+cb = fig.colorbar(im)
+cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18 )
+plt.show()
+
+We can see that our model for the least squares method performes close
+to the benchmark from Scikit-learn. It is interesting to note that OLS
+considers both \( J_{j, j + 1} = -0.5 \) and \( J_{j, j - 1} = -0.5 \) as
+valid matrix elements for \( J \).
+
+
+
+
+Ridge regression
+
+
+Having explored the ordinary least squares we move on to ridge
+regression. In ridge regression we include a regularizer . This
+involves a new cost function which leads to a new estimate for the
+weights \( \omega \). This results in a penalized regression problem. The
+cost function is given by
+
+
+where \( \omega_{\text{LS}} \) is the weights from ordinary least
+squares. The last assumption assumes that \( X \) is orthogonal, which it
+is not. We will therefore resort to solving the equation as it stands
+on the left hand side.
+
+
+
+
+
def get_ridge_weights (x: np.ndarray, y: np.ndarray, _lambda: float ) -> np.ndarray:
+ return x.T @ y @ scl.inv(
+ x.T @ x + np.eye(x.shape[1 ], x.shape[1 ]) * _lambda
+ )
+lambda = 0.1
+omega_ridge = get_ridge_weights(X_train_own, y_train, np.array([_lambda]))
+clf_ridge = skl.Ridge(alpha=_lambda).fit(X_train, y_train)
+J_ridge_own = omega_ridge[1 :].reshape(L, L)
+J_ridge_sk = clf_ridge.coef_.reshape(L, L)
+fig = plt.figure(figsize=(20 , 14 ))
+im = plt.imshow(J_ridge_own, **cmap_args)
+plt.title("Home-made ridge regression" , fontsize=18 )
+plt.xticks(fontsize=18 )
+plt.yticks(fontsize=18 )
+cb = fig.colorbar(im)
+cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18 )
+
+fig = plt.figure(figsize=(20 , 14 ))
+im = plt.imshow(J_ridge_sk, **cmap_args)
+plt.title("Ridge from Scikit-learn" , fontsize=18 )
+plt.xticks(fontsize=18 )
+plt.yticks(fontsize=18 )
+cb = fig.colorbar(im)
+cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18 )
+
+plt.show()
+
+
+
+
+LASSO regression
+
+
+In the Least Absolute Shrinkage and Selection Operator (LASSO)-method we get a third cost function.
+
+
+
+
20 , 14 ))
+im = plt.imshow(J_lasso_sk, **cmap_args)
+plt.title("Lasso from Scikit-learn" , fontsize=18 )
+plt.xticks(fontsize=18 )
+plt.yticks(fontsize=18 )
+cb = fig.colorbar(im)
+cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18 )
+
+plt.show()
+
+It is quite striking how LASSO breaks the symmetry of the coupling
+constant as opposed to ridge and OLS. We get a sparse solution with
+\( J_{j, j + 1} = -1 \).
+
+
+
+
+Performance of the different models
+
+
+In order to judge which model performs best at varying values of \( \lambda \) (for ridge and LASSO) we compute \( R^2 \) which is given by
+
+
+
+
+
def r_squared (y, y_hat):
+ return 1 - np.sum((y - y_hat) ** 2 ) / np.sum((y - np.mean(y_hat)) ** 2 )
+
+This is the same metric used by Scikit-learn for their regression models when scoring.
+
+
+
+
assert abs (r_test - sk_r_test) < 1e-2
+
+
+
+
+Performance as function of the regularization parameter
+
+
+We see how the different models perform for a different set of values for \( \lambda \).
+
+
+
+
+
4 , 5 , 10 )
+
+train_errors = {
+ "ols_own" : np.zeros(lambdas.size),
+ "ols_sk" : np.zeros(lambdas.size),
+ "ridge_own" : np.zeros(lambdas.size),
+ "ridge_sk" : np.zeros(lambdas.size),
+ "lasso_sk" : np.zeros(lambdas.size)
+}
+
+test_errors = {
+ "ols_own" : np.zeros(lambdas.size),
+ "ols_sk" : np.zeros(lambdas.size),
+ "ridge_own" : np.zeros(lambdas.size),
+ "ridge_sk" : np.zeros(lambdas.size),
+ "lasso_sk" : np.zeros(lambdas.size)
+}
+
+plot_counter = 1
+
+fig = plt.figure(figsize=(32 , 54 ))
+
+for i, _lambda in enumerate (tqdm.tqdm(lambdas)):
+ omega = get_ols_weights(X_train_own, y_train)
+ y_hat_train = X_train_own @ omega
+ y_hat_test = X_test_own @ omega
+
+ train_errors["ols_own" ][i] = r_squared(y_train, y_hat_train)
+ test_errors["ols_own" ][i] = r_squared(y_test, y_hat_test)
+
+ plt.subplot(10 , 5 , plot_counter)
+ plt.imshow(omega[1 :].reshape(L, L), **cmap_args)
+ plt.title("Home made OLS" )
+ plot_counter += 1
+
+ omega = get_ridge_weights(X_train_own, y_train, _lambda)
+ y_hat_train = X_train_own @ omega
+ y_hat_test = X_test_own @ omega
+
+ train_errors["ridge_own" ][i] = r_squared(y_train, y_hat_train)
+ test_errors["ridge_own" ][i] = r_squared(y_test, y_hat_test)
+
+ plt.subplot(10 , 5 , plot_counter)
+ plt.imshow(omega[1 :].reshape(L, L), **cmap_args)
+ plt.title(r"Home made ridge, $\lambda = %.4f$" % _lambda)
+ plot_counter += 1
+
+ for key, method in zip (
+ ["ols_sk" , "ridge_sk" , "lasso_sk" ],
+ [skl.LinearRegression(), skl.Ridge(alpha=_lambda), skl.Lasso(alpha=_lambda)]
+ ):
+ method = method.fit(X_train, y_train)
+
+ train_errors[key][i] = method.score(X_train, y_train)
+ test_errors[key][i] = method.score(X_test, y_test)
+
+ omega = method.coef_.reshape(L, L)
+
+ plt.subplot(10 , 5 , plot_counter)
+ plt.imshow(omega, **cmap_args)
+ plt.title(r"%s, $\lambda = %.4f$" % (key, _lambda))
+ plot_counter += 1
+
+plt.show()
+
+We can see that LASSO quite fast reaches a good solution for low
+values of \( \lambda \), but will "wither" when we increase \( \lambda \) too
+much. Ridge is more stable over a larger range of values for
+\( \lambda \), but eventually also fades away.
+
+
+
+
+Finding the optimal value of \( \lambda \)
+
+
+To determine which value of \( \lambda \) is best we plot the accuracy of
+the models when predicting the training and the testing set. We expect
+the accuracy of the training set to be quite good, but if the accuracy
+of the testing set is much lower this tells us that we might be
+subject to an overfit model. The ideal scenario is an accuracy on the
+testing set that is close to the accuracy of the training set.
+
+
+
+
+
20 , 14 ))
+
+colors = {
+ "ols_own" : "b" ,
+ "ridge_own" : "g" ,
+ "ols_sk" : "r" ,
+ "ridge_sk" : "y" ,
+ "lasso_sk" : "c"
+}
+
+for key in train_errors:
+ plt.semilogx(
+ lambdas,
+ train_errors[key],
+ colors[key],
+ label="Train {0}" .format(key),
+ linewidth=4.0
+ )
+
+for key in test_errors:
+ plt.semilogx(
+ lambdas,
+ test_errors[key],
+ colors[key] + "--" ,
+ label="Test {0}" .format(key),
+ linewidth=4.0
+ )
+#plt.semilogx(lambdas, train_errors["ols_own"], label="Train (OLS own)")
+#plt.semilogx(lambdas, test_errors["ols_own"], label="Test (OLS own)")
+
+plt.legend(loc="best" , fontsize=18 )
+plt.xlabel(r"$\lambda$" , fontsize=18 )
+plt.ylabel(r"$R^2$" , fontsize=18 )
+plt.tick_params(labelsize=18 )
+plt.show()
+
+From the above figure we can see that LASSO with \( \lambda = 10^{-2} \)
+achieve a very good accuracy on the test set. This by far surpases the
+other models for all values of \( \lambda \).
+
+
+
diff --git a/doc/pub/Regression/html/Regression-solarized.html b/doc/pub/Regression/html/Regression-solarized.html
index df9b24537..95341ec36 100644
--- a/doc/pub/Regression/html/Regression-solarized.html
+++ b/doc/pub/Regression/html/Regression-solarized.html
@@ -245,7 +245,34 @@ div { text-align: justify; text-justify: inter-word; }
('Elaborating a little bit more', 2, None, '___sec98'),
('The bias', 2, None, '___sec99'),
('The variance', 2, None, '___sec100'),
- ('Summing up', 2, None, '___sec101')]}
+ ('Summing up', 2, None, '___sec101'),
+ ('The one-dimensional Ising model, project 2',
+ 2,
+ None,
+ '___sec102'),
+ ('Example: The one-dimensional Ising model',
+ 2,
+ None,
+ '___sec103'),
+ ('Reformulating the problem to suit regression',
+ 2,
+ None,
+ '___sec104'),
+ ('Linear regression', 2, None, '___sec105'),
+ ('Ordinary least squares', 2, None, '___sec106'),
+ ('Singular Value decomposition', 2, None, '___sec107'),
+ ('Fitting with scikit-learn', 2, None, '___sec108'),
+ ('Ridge regression', 2, None, '___sec109'),
+ ('LASSO regression', 2, None, '___sec110'),
+ ('Performance of the different models', 2, None, '___sec111'),
+ ('Performance as function of the regularization parameter',
+ 2,
+ None,
+ '___sec112'),
+ ('Finding the optimal value of $\\lambda$',
+ 2,
+ None,
+ '___sec113')]}
end of tocinfo -->
@@ -287,7 +314,7 @@ MathJax.Hub.Config({
[2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University
-
Sep 25, 2018 Oct 11, 2018
+
The one-dimensional Ising model, project 2
+
+
+The one-dimensional Ising model with nearest neighbor interaction, no external field and a constant coupling constant \( J \) is given by
+
+$$
+\begin{align}
+ H = -J \sum_{k}^L s_k s_{k + 1},
+\label{_auto19}
+\end{align}
+$$
+
+where \( s_i \in \{-1, 1\} \) and \( s_{N + 1} = s_1 \). The number of spins in the system is determined by \( L \). For the low temperature limit there is no phase transition.
+
+
+We will look at a system of \( L = 40 \) spins with a coupling constant of \( J = 1 \). To get enough training data we will generate 10000 states with their respective energies.
+
+
+
+
+
40
+n = int (1e4 )
+
+spins = np.random.choice([-1 , 1 ], size=(n, L))
+J = 1.0
+
+energies = np.zeros(n)
+
+for i in range (n):
+ energies[i] = - J * np.dot(spins[i], np.roll(spins[i], 1 ))
+
+
Example: The one-dimensional Ising model
+
+
+Here we use linear (ordinary least squares), ridge and LASSO
+regression to predict the energy in the nearest neighbor
+one-dimensional Ising model on a ring, i.e., the endpoints wrap
+around. We will use the linear regression models to fit a value for
+the coupling constant to achieve this.
+
+
+
+
import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+import seaborn as sns
+import scipy.linalg as scl
+from sklearn.model_selection import train_test_split
+import sklearn.linear_model as skl
+import tqdm
+sns.set(color_codes=True )
+cmap_args=dict (vmin=-1. , vmax=1. , cmap='seismic' )
+
+
Reformulating the problem to suit regression
+
+
+A more general form for the one-dimensional Ising model is
+
+$$
+\begin{align}
+ H = - \sum_j^L \sum_k^L s_j s_k J_{jk}.
+\label{_auto20}
+\end{align}
+$$
+
+
+Here we allow for interactions beyond the nearest neighbors and a more
+adaptive coupling matrix. This latter expression can be formulated as
+a matrix-product on the form
+$$
+\begin{align}
+ H = X J,
+\label{_auto21}
+\end{align}
+$$
+
+
+where \( X_{jk} = s_j s_k \) and \( J \) is the matrix consisting of the
+elements \( -J_{jk} \). This form of writing the energy fits perfectly
+with the form utilized in linear regression, viz.
+$$
+\begin{align}
+ y = X\omega + \epsilon,
+\label{_auto22}
+\end{align}
+$$
+
+
+
+
+
2 ))
+for i in range (n):
+ X[i] = np.outer(spins[i], spins[i]).ravel()
+y = energies
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.96 )
+
+X_train_own = np.concatenate(
+ (np.ones(len (X_train))[:, np.newaxis], X_train),
+ axis=1
+)
+
+X_test_own = np.concatenate(
+ (np.ones(len (X_test))[:, np.newaxis], X_test),
+ axis=1
+)
+
+
Linear regression
+
+
+The problem at hand is to try to fit the equation
+$$
+\begin{align}
+ y = f(x) + \epsilon,
+\label{_auto23}
+\end{align}
+$$
+
+
+where \( f(x) \) is some unknown function of the data \( x \) and \( \epsilon \)
+is normally distributed with mean zero noise with standard deviation
+\( \sigma_{\epsilon} \). Our job is to try to find a predictor which
+estimates the function \( f(x) \). In linear regression we assume that we
+can formulate the problem as
+
+$$
+\begin{align}
+ y = X\omega + \epsilon,
+\label{_auto24}
+\end{align}
+$$
+
+
+where \( X \) and \( \omega \) are now matrices. Our job at hand is now to
+find a cost function \( C \), which we wish to minimize in order to find
+the best estimate of \( \omega \).
+
+
+
Ordinary least squares
+
+
+In the ordinary least squares method we choose the cost function
+$$
+\begin{align}
+ C(X, \omega) = ||X\omega - y||^2
+ = (X\omega - y)^T(X\omega - y)
+\label{_auto25}
+\end{align}
+$$
+
+We then find the extremal point of \( C \) by taking the derivative with respect to \( \omega \) and setting it to zero, i.e.,
+
+$$
+\begin{align}
+ \dfrac{\mathrm{d}C}{\mathrm{d}\omega}
+ = 0.
+\label{_auto26}
+\end{align}
+$$
+
+This yields the expression for \( \omega \) to be
+$$
+\begin{align}
+ \omega = \frac{X^T y}{X^T X},
+\label{_auto27}
+\end{align}
+$$
+
+
+which immediately imposes some requirements on \( X \) as there must exist
+an inverse of \( X^T X \). If the expression we are modelling contains an
+intercept, i.e., a constant expression we must make sure that the
+first column of \( X \) consists of \( 1 \).
+
+
+
+
+
def get_ols_weights_naive (x: np.ndarray, y: np.ndarray) -> np.ndarray:
+ return scl.inv(x.T @ x) @ (x.T @ y)
+omega = get_ols_weights_naive(X_train_own, y_train)
+
+
Singular Value decomposition
+Doing the inversion directly turns out to be a bad idea as the matrix
+\( X^TX \) is singular. An alternative approach is to use the singular
+value decomposition . Using the definition of the Moore-Penrose
+pseudoinverse we can write the equation for \( \omega \) as
+
+$$
+\begin{align}
+ \omega = X^{+}y,
+\label{_auto28}
+\end{align}
+$$
+
+where the pseudoinverse of \( X \) is given by
+$$
+\begin{align}
+ X^{+} = \frac{X^T}{X^T X}.
+\label{_auto29}
+\end{align}
+$$
+
+
+Using singular value decomposition we have that \( X = U\Sigma V^T \),
+where \( X^{+} = V\Sigma^{+} U^T \). This reduces the equation for
+\( \omega \) to
+$$
+\begin{align}
+ \omega = V\Sigma^{+} U^T y.
+\label{_auto30}
+\end{align}
+$$
+
+
+Note that solving this equation by actually doing the pseudoinverse
+(which is what we will do) is not a good idea as this operation scales
+as \( \mathcal{O}(n^3) \), where \( n \) is the number of elements in a
+general matrix. Instead, doing \( QR \)-factorization and solving the
+linear system as an equation would reduce this down to
+\( \mathcal{O}(n^2) \) operations.
+
+
+
+
+
def get_ols_weights (x: np.ndarray, y: np.ndarray) -> np.ndarray:
+ u, s, v = scl.svd(x)
+ return v.T @ scl.pinv(scl.diagsvd(s, u.shape[0 ], v.shape[0 ])) @ u.T @ y
+
+Before passing in the data to the function we append a column with ones to the training data.
+
+
+
+
+
+
Fitting with scikit-learn
+
+
+Next we fit a LinearRegression-model from Scikit-learn for comparison.
+
+
+
+
+
+Extracting the \( J \)-matrix from both our own method and the Scikit-learn model where we make sure to remove the intercept.
+
+
+
+
+
1 :].reshape(L, L)
+J_sk = clf.coef_.reshape(L, L)
+
+A way of looking at the coefficients in \( J \) is to plot the matrices as images.
+
+
+
+
+
20 , 14 ))
+im = plt.imshow(J_own, **cmap_args)
+plt.title("Home-made OLS" , fontsize=18 )
+plt.xticks(fontsize=18 )
+plt.yticks(fontsize=18 )
+cb = fig.colorbar(im)
+cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18 )
+
+fig = plt.figure(figsize=(20 , 14 ))
+im = plt.imshow(J_sk, **cmap_args)
+plt.title("LinearRegression from Scikit-learn" , fontsize=18 )
+plt.xticks(fontsize=18 )
+plt.yticks(fontsize=18 )
+cb = fig.colorbar(im)
+cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18 )
+plt.show()
+
+We can see that our model for the least squares method performes close
+to the benchmark from Scikit-learn. It is interesting to note that OLS
+considers both \( J_{j, j + 1} = -0.5 \) and \( J_{j, j - 1} = -0.5 \) as
+valid matrix elements for \( J \).
+
+
+
Ridge regression
+
+
+Having explored the ordinary least squares we move on to ridge
+regression. In ridge regression we include a regularizer . This
+involves a new cost function which leads to a new estimate for the
+weights \( \omega \). This results in a penalized regression problem. The
+cost function is given by
+
+$$
+\begin{align}
+ C(X, \omega; \lambda) = ||X\omega - y||^2 + \lambda ||\omega||^2
+ = (X\omega - y)^T(X\omega - y) + \lambda \omega^T\omega.
+\label{_auto31}
+\end{align}
+$$
+
+Finding the extremum of this function yields the weights
+
+$$
+\begin{align}
+ \omega(\lambda) = \frac{X^Ty}{X^TX + \lambda} \to \frac{\omega_{\text{LS}}}{1 + \lambda},
+\label{_auto32}
+\end{align}
+$$
+
+
+where \( \omega_{\text{LS}} \) is the weights from ordinary least
+squares. The last assumption assumes that \( X \) is orthogonal, which it
+is not. We will therefore resort to solving the equation as it stands
+on the left hand side.
+
+
+
+
+
def get_ridge_weights (x: np.ndarray, y: np.ndarray, _lambda: float ) -> np.ndarray:
+ return x.T @ y @ scl.inv(
+ x.T @ x + np.eye(x.shape[1 ], x.shape[1 ]) * _lambda
+ )
+lambda = 0.1
+omega_ridge = get_ridge_weights(X_train_own, y_train, np.array([_lambda]))
+clf_ridge = skl.Ridge(alpha=_lambda).fit(X_train, y_train)
+J_ridge_own = omega_ridge[1 :].reshape(L, L)
+J_ridge_sk = clf_ridge.coef_.reshape(L, L)
+fig = plt.figure(figsize=(20 , 14 ))
+im = plt.imshow(J_ridge_own, **cmap_args)
+plt.title("Home-made ridge regression" , fontsize=18 )
+plt.xticks(fontsize=18 )
+plt.yticks(fontsize=18 )
+cb = fig.colorbar(im)
+cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18 )
+
+fig = plt.figure(figsize=(20 , 14 ))
+im = plt.imshow(J_ridge_sk, **cmap_args)
+plt.title("Ridge from Scikit-learn" , fontsize=18 )
+plt.xticks(fontsize=18 )
+plt.yticks(fontsize=18 )
+cb = fig.colorbar(im)
+cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18 )
+
+plt.show()
+
+
LASSO regression
+
+
+In the Least Absolute Shrinkage and Selection Operator (LASSO)-method we get a third cost function.
+$$
+\begin{align}
+ C(X, \omega; \lambda) =
+ ||X\omega - y||^2 + \lambda ||\omega||
+ = (X\omega - y)^T(X\omega - y) + \lambda \sqrt{\omega^T\omega}.
+\label{_auto33}
+\end{align}
+$$
+
+Finding the extremal point of this cost function is not so straight-forward as in least squares and ridge. We will therefore rely solely on the function ``Lasso`` from Scikit-learn.
+
+
+
+
+
20 , 14 ))
+im = plt.imshow(J_lasso_sk, **cmap_args)
+plt.title("Lasso from Scikit-learn" , fontsize=18 )
+plt.xticks(fontsize=18 )
+plt.yticks(fontsize=18 )
+cb = fig.colorbar(im)
+cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18 )
+
+plt.show()
+
+It is quite striking how LASSO breaks the symmetry of the coupling
+constant as opposed to ridge and OLS. We get a sparse solution with
+\( J_{j, j + 1} = -1 \).
+
+
+
Performance of the different models
+
+
+In order to judge which model performs best at varying values of \( \lambda \) (for ridge and LASSO) we compute \( R^2 \) which is given by
+
+$$
+\begin{align}
+ R^2 = 1 - \frac{(y - \hat{y})^2}{(y - \bar{y})^2},
+\label{_auto34}
+\end{align}
+$$
+
+where \( y \) is a vector with the true values of the energy, \( \hat{y} \) is the predicted values of \( y \) from the models and \( \bar{y} \) is the mean of \( \hat{y} \).
+
+
+
+
+
def r_squared (y, y_hat):
+ return 1 - np.sum((y - y_hat) ** 2 ) / np.sum((y - np.mean(y_hat)) ** 2 )
+
+This is the same metric used by Scikit-learn for their regression models when scoring.
+
+
+
+
assert abs (r_test - sk_r_test) < 1e-2
+
+
Performance as function of the regularization parameter
+
+
+We see how the different models perform for a different set of values for \( \lambda \).
+
+
+
+
+
4 , 5 , 10 )
+
+train_errors = {
+ "ols_own" : np.zeros(lambdas.size),
+ "ols_sk" : np.zeros(lambdas.size),
+ "ridge_own" : np.zeros(lambdas.size),
+ "ridge_sk" : np.zeros(lambdas.size),
+ "lasso_sk" : np.zeros(lambdas.size)
+}
+
+test_errors = {
+ "ols_own" : np.zeros(lambdas.size),
+ "ols_sk" : np.zeros(lambdas.size),
+ "ridge_own" : np.zeros(lambdas.size),
+ "ridge_sk" : np.zeros(lambdas.size),
+ "lasso_sk" : np.zeros(lambdas.size)
+}
+
+plot_counter = 1
+
+fig = plt.figure(figsize=(32 , 54 ))
+
+for i, _lambda in enumerate (tqdm.tqdm(lambdas)):
+ omega = get_ols_weights(X_train_own, y_train)
+ y_hat_train = X_train_own @ omega
+ y_hat_test = X_test_own @ omega
+
+ train_errors["ols_own" ][i] = r_squared(y_train, y_hat_train)
+ test_errors["ols_own" ][i] = r_squared(y_test, y_hat_test)
+
+ plt.subplot(10 , 5 , plot_counter)
+ plt.imshow(omega[1 :].reshape(L, L), **cmap_args)
+ plt.title("Home made OLS" )
+ plot_counter += 1
+
+ omega = get_ridge_weights(X_train_own, y_train, _lambda)
+ y_hat_train = X_train_own @ omega
+ y_hat_test = X_test_own @ omega
+
+ train_errors["ridge_own" ][i] = r_squared(y_train, y_hat_train)
+ test_errors["ridge_own" ][i] = r_squared(y_test, y_hat_test)
+
+ plt.subplot(10 , 5 , plot_counter)
+ plt.imshow(omega[1 :].reshape(L, L), **cmap_args)
+ plt.title(r"Home made ridge, $\lambda = %.4f$" % _lambda)
+ plot_counter += 1
+
+ for key, method in zip (
+ ["ols_sk" , "ridge_sk" , "lasso_sk" ],
+ [skl.LinearRegression(), skl.Ridge(alpha=_lambda), skl.Lasso(alpha=_lambda)]
+ ):
+ method = method.fit(X_train, y_train)
+
+ train_errors[key][i] = method.score(X_train, y_train)
+ test_errors[key][i] = method.score(X_test, y_test)
+
+ omega = method.coef_.reshape(L, L)
+
+ plt.subplot(10 , 5 , plot_counter)
+ plt.imshow(omega, **cmap_args)
+ plt.title(r"%s, $\lambda = %.4f$" % (key, _lambda))
+ plot_counter += 1
+
+plt.show()
+
+We can see that LASSO quite fast reaches a good solution for low
+values of \( \lambda \), but will "wither" when we increase \( \lambda \) too
+much. Ridge is more stable over a larger range of values for
+\( \lambda \), but eventually also fades away.
+
+
+
Finding the optimal value of \( \lambda \)
+
+
+To determine which value of \( \lambda \) is best we plot the accuracy of
+the models when predicting the training and the testing set. We expect
+the accuracy of the training set to be quite good, but if the accuracy
+of the testing set is much lower this tells us that we might be
+subject to an overfit model. The ideal scenario is an accuracy on the
+testing set that is close to the accuracy of the training set.
+
+
+
+
+
20 , 14 ))
+
+colors = {
+ "ols_own" : "b" ,
+ "ridge_own" : "g" ,
+ "ols_sk" : "r" ,
+ "ridge_sk" : "y" ,
+ "lasso_sk" : "c"
+}
+
+for key in train_errors:
+ plt.semilogx(
+ lambdas,
+ train_errors[key],
+ colors[key],
+ label="Train {0}" .format(key),
+ linewidth=4.0
+ )
+
+for key in test_errors:
+ plt.semilogx(
+ lambdas,
+ test_errors[key],
+ colors[key] + "--" ,
+ label="Test {0}" .format(key),
+ linewidth=4.0
+ )
+#plt.semilogx(lambdas, train_errors["ols_own"], label="Train (OLS own)")
+#plt.semilogx(lambdas, test_errors["ols_own"], label="Test (OLS own)")
+
+plt.legend(loc="best" , fontsize=18 )
+plt.xlabel(r"$\lambda$" , fontsize=18 )
+plt.ylabel(r"$R^2$" , fontsize=18 )
+plt.tick_params(labelsize=18 )
+plt.show()
+
+From the above figure we can see that LASSO with \( \lambda = 10^{-2} \)
+achieve a very good accuracy on the test set. This by far surpases the
+other models for all values of \( \lambda \).
+
diff --git a/doc/pub/Regression/html/Regression.html b/doc/pub/Regression/html/Regression.html
index 3681d2511..af11e2365 100644
--- a/doc/pub/Regression/html/Regression.html
+++ b/doc/pub/Regression/html/Regression.html
@@ -250,7 +250,34 @@ div { text-align: justify; text-justify: inter-word; }
('Elaborating a little bit more', 2, None, '___sec98'),
('The bias', 2, None, '___sec99'),
('The variance', 2, None, '___sec100'),
- ('Summing up', 2, None, '___sec101')]}
+ ('Summing up', 2, None, '___sec101'),
+ ('The one-dimensional Ising model, project 2',
+ 2,
+ None,
+ '___sec102'),
+ ('Example: The one-dimensional Ising model',
+ 2,
+ None,
+ '___sec103'),
+ ('Reformulating the problem to suit regression',
+ 2,
+ None,
+ '___sec104'),
+ ('Linear regression', 2, None, '___sec105'),
+ ('Ordinary least squares', 2, None, '___sec106'),
+ ('Singular Value decomposition', 2, None, '___sec107'),
+ ('Fitting with scikit-learn', 2, None, '___sec108'),
+ ('Ridge regression', 2, None, '___sec109'),
+ ('LASSO regression', 2, None, '___sec110'),
+ ('Performance of the different models', 2, None, '___sec111'),
+ ('Performance as function of the regularization parameter',
+ 2,
+ None,
+ '___sec112'),
+ ('Finding the optimal value of $\\lambda$',
+ 2,
+ None,
+ '___sec113')]}
end of tocinfo -->
@@ -292,7 +319,7 @@ MathJax.Hub.Config({
[2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University
-
Sep 25, 2018 Oct 11, 2018
+
The one-dimensional Ising model, project 2
+
+
+The one-dimensional Ising model with nearest neighbor interaction, no external field and a constant coupling constant \( J \) is given by
+
+$$
+\begin{align}
+ H = -J \sum_{k}^L s_k s_{k + 1},
+\label{_auto19}
+\end{align}
+$$
+
+where \( s_i \in \{-1, 1\} \) and \( s_{N + 1} = s_1 \). The number of spins in the system is determined by \( L \). For the low temperature limit there is no phase transition.
+
+
+We will look at a system of \( L = 40 \) spins with a coupling constant of \( J = 1 \). To get enough training data we will generate 10000 states with their respective energies.
+
+
+
+
+
= 40
+n = int (1e4 )
+
+spins = np. random. choice([-1 , 1 ], size= (n, L))
+J = 1.0
+
+energies = np. zeros(n)
+
+for i in range (n):
+ energies[i] = - J * np. dot(spins[i], np. roll(spins[i], 1 ))
+
+
Example: The one-dimensional Ising model
+
+
+Here we use linear (ordinary least squares), ridge and LASSO
+regression to predict the energy in the nearest neighbor
+one-dimensional Ising model on a ring, i.e., the endpoints wrap
+around. We will use the linear regression models to fit a value for
+the coupling constant to achieve this.
+
+
+
+
import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+import seaborn as sns
+import scipy.linalg as scl
+from sklearn.model_selection import train_test_split
+import sklearn.linear_model as skl
+import tqdm
+sns. set(color_codes= True )
+cmap_args= dict (vmin=-1. , vmax=1. , cmap= 'seismic' )
+
+
Reformulating the problem to suit regression
+
+
+A more general form for the one-dimensional Ising model is
+
+$$
+\begin{align}
+ H = - \sum_j^L \sum_k^L s_j s_k J_{jk}.
+\label{_auto20}
+\end{align}
+$$
+
+
+Here we allow for interactions beyond the nearest neighbors and a more
+adaptive coupling matrix. This latter expression can be formulated as
+a matrix-product on the form
+$$
+\begin{align}
+ H = X J,
+\label{_auto21}
+\end{align}
+$$
+
+
+where \( X_{jk} = s_j s_k \) and \( J \) is the matrix consisting of the
+elements \( -J_{jk} \). This form of writing the energy fits perfectly
+with the form utilized in linear regression, viz.
+$$
+\begin{align}
+ y = X\omega + \epsilon,
+\label{_auto22}
+\end{align}
+$$
+
+
+
+
+
= np. zeros((n, L ** 2 ))
+for i in range (n):
+ X[i] = np. outer(spins[i], spins[i]). ravel()
+y = energies
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.96 )
+
+X_train_own = np. concatenate(
+ (np. ones(len (X_train))[:, np. newaxis], X_train),
+ axis=1
+)
+
+X_test_own = np. concatenate(
+ (np. ones(len (X_test))[:, np. newaxis], X_test),
+ axis=1
+)
+
+
Linear regression
+
+
+The problem at hand is to try to fit the equation
+$$
+\begin{align}
+ y = f(x) + \epsilon,
+\label{_auto23}
+\end{align}
+$$
+
+
+where \( f(x) \) is some unknown function of the data \( x \) and \( \epsilon \)
+is normally distributed with mean zero noise with standard deviation
+\( \sigma_{\epsilon} \). Our job is to try to find a predictor which
+estimates the function \( f(x) \). In linear regression we assume that we
+can formulate the problem as
+
+$$
+\begin{align}
+ y = X\omega + \epsilon,
+\label{_auto24}
+\end{align}
+$$
+
+
+where \( X \) and \( \omega \) are now matrices. Our job at hand is now to
+find a cost function \( C \), which we wish to minimize in order to find
+the best estimate of \( \omega \).
+
+
+
Ordinary least squares
+
+
+In the ordinary least squares method we choose the cost function
+$$
+\begin{align}
+ C(X, \omega) = ||X\omega - y||^2
+ = (X\omega - y)^T(X\omega - y)
+\label{_auto25}
+\end{align}
+$$
+
+We then find the extremal point of \( C \) by taking the derivative with respect to \( \omega \) and setting it to zero, i.e.,
+
+$$
+\begin{align}
+ \dfrac{\mathrm{d}C}{\mathrm{d}\omega}
+ = 0.
+\label{_auto26}
+\end{align}
+$$
+
+This yields the expression for \( \omega \) to be
+$$
+\begin{align}
+ \omega = \frac{X^T y}{X^T X},
+\label{_auto27}
+\end{align}
+$$
+
+
+which immediately imposes some requirements on \( X \) as there must exist
+an inverse of \( X^T X \). If the expression we are modelling contains an
+intercept, i.e., a constant expression we must make sure that the
+first column of \( X \) consists of \( 1 \).
+
+
+
+
+
def get_ols_weights_naive (x: np. ndarray, y: np. ndarray) -> np. ndarray:
+ return scl. inv(x. T @ x) @ (x. T @ y)
+omega = get_ols_weights_naive(X_train_own, y_train)
+
+
Singular Value decomposition
+Doing the inversion directly turns out to be a bad idea as the matrix
+\( X^TX \) is singular. An alternative approach is to use the singular
+value decomposition . Using the definition of the Moore-Penrose
+pseudoinverse we can write the equation for \( \omega \) as
+
+$$
+\begin{align}
+ \omega = X^{+}y,
+\label{_auto28}
+\end{align}
+$$
+
+where the pseudoinverse of \( X \) is given by
+$$
+\begin{align}
+ X^{+} = \frac{X^T}{X^T X}.
+\label{_auto29}
+\end{align}
+$$
+
+
+Using singular value decomposition we have that \( X = U\Sigma V^T \),
+where \( X^{+} = V\Sigma^{+} U^T \). This reduces the equation for
+\( \omega \) to
+$$
+\begin{align}
+ \omega = V\Sigma^{+} U^T y.
+\label{_auto30}
+\end{align}
+$$
+
+
+Note that solving this equation by actually doing the pseudoinverse
+(which is what we will do) is not a good idea as this operation scales
+as \( \mathcal{O}(n^3) \), where \( n \) is the number of elements in a
+general matrix. Instead, doing \( QR \)-factorization and solving the
+linear system as an equation would reduce this down to
+\( \mathcal{O}(n^2) \) operations.
+
+
+
+
+
def get_ols_weights (x: np. ndarray, y: np. ndarray) -> np. ndarray:
+ u, s, v = scl. svd(x)
+ return v. T @ scl. pinv(scl. diagsvd(s, u. shape[0 ], v. shape[0 ])) @ u. T @ y
+
+Before passing in the data to the function we append a column with ones to the training data.
+
+
+
+
+
= get_ols_weights(X_train_own,y_train)
+
+
Fitting with scikit-learn
+
+
+Next we fit a LinearRegression-model from Scikit-learn for comparison.
+
+
+
+
+
= skl. LinearRegression(). fit(X_train, y_train)
+
+Extracting the \( J \)-matrix from both our own method and the Scikit-learn model where we make sure to remove the intercept.
+
+
+
+
+
= omega[1 :]. reshape(L, L)
+J_sk = clf. coef_. reshape(L, L)
+
+A way of looking at the coefficients in \( J \) is to plot the matrices as images.
+
+
+
+
+
= plt. figure(figsize= (20 , 14 ))
+im = plt. imshow(J_own, ** cmap_args)
+plt. title("Home-made OLS" , fontsize=18 )
+plt. xticks(fontsize=18 )
+plt. yticks(fontsize=18 )
+cb = fig. colorbar(im)
+cb. ax. set_yticklabels(cb. ax. get_yticklabels(), fontsize=18 )
+
+fig = plt. figure(figsize= (20 , 14 ))
+im = plt. imshow(J_sk, ** cmap_args)
+plt. title("LinearRegression from Scikit-learn" , fontsize=18 )
+plt. xticks(fontsize=18 )
+plt. yticks(fontsize=18 )
+cb = fig. colorbar(im)
+cb. ax. set_yticklabels(cb. ax. get_yticklabels(), fontsize=18 )
+plt. show()
+
+We can see that our model for the least squares method performes close
+to the benchmark from Scikit-learn. It is interesting to note that OLS
+considers both \( J_{j, j + 1} = -0.5 \) and \( J_{j, j - 1} = -0.5 \) as
+valid matrix elements for \( J \).
+
+
+
Ridge regression
+
+
+Having explored the ordinary least squares we move on to ridge
+regression. In ridge regression we include a regularizer . This
+involves a new cost function which leads to a new estimate for the
+weights \( \omega \). This results in a penalized regression problem. The
+cost function is given by
+
+$$
+\begin{align}
+ C(X, \omega; \lambda) = ||X\omega - y||^2 + \lambda ||\omega||^2
+ = (X\omega - y)^T(X\omega - y) + \lambda \omega^T\omega.
+\label{_auto31}
+\end{align}
+$$
+
+Finding the extremum of this function yields the weights
+
+$$
+\begin{align}
+ \omega(\lambda) = \frac{X^Ty}{X^TX + \lambda} \to \frac{\omega_{\text{LS}}}{1 + \lambda},
+\label{_auto32}
+\end{align}
+$$
+
+
+where \( \omega_{\text{LS}} \) is the weights from ordinary least
+squares. The last assumption assumes that \( X \) is orthogonal, which it
+is not. We will therefore resort to solving the equation as it stands
+on the left hand side.
+
+
+
+
+
def get_ridge_weights (x: np. ndarray, y: np. ndarray, _lambda: float ) -> np. ndarray:
+ return x. T @ y @ scl. inv(
+ x. T @ x + np. eye(x. shape[1 ], x. shape[1 ]) * _lambda
+ )
+lambda = 0.1
+omega_ridge = get_ridge_weights(X_train_own, y_train, np. array([_lambda]))
+clf_ridge = skl. Ridge(alpha= _lambda). fit(X_train, y_train)
+J_ridge_own = omega_ridge[1 :]. reshape(L, L)
+J_ridge_sk = clf_ridge. coef_. reshape(L, L)
+fig = plt. figure(figsize= (20 , 14 ))
+im = plt. imshow(J_ridge_own, ** cmap_args)
+plt. title("Home-made ridge regression" , fontsize=18 )
+plt. xticks(fontsize=18 )
+plt. yticks(fontsize=18 )
+cb = fig. colorbar(im)
+cb. ax. set_yticklabels(cb. ax. get_yticklabels(), fontsize=18 )
+
+fig = plt. figure(figsize= (20 , 14 ))
+im = plt. imshow(J_ridge_sk, ** cmap_args)
+plt. title("Ridge from Scikit-learn" , fontsize=18 )
+plt. xticks(fontsize=18 )
+plt. yticks(fontsize=18 )
+cb = fig. colorbar(im)
+cb. ax. set_yticklabels(cb. ax. get_yticklabels(), fontsize=18 )
+
+plt. show()
+
+
LASSO regression
+
+
+In the Least Absolute Shrinkage and Selection Operator (LASSO)-method we get a third cost function.
+$$
+\begin{align}
+ C(X, \omega; \lambda) =
+ ||X\omega - y||^2 + \lambda ||\omega||
+ = (X\omega - y)^T(X\omega - y) + \lambda \sqrt{\omega^T\omega}.
+\label{_auto33}
+\end{align}
+$$
+
+Finding the extremal point of this cost function is not so straight-forward as in least squares and ridge. We will therefore rely solely on the function ``Lasso`` from Scikit-learn.
+
+
+
+
+
= skl. Lasso(alpha= _lambda). fit(X_train, y_train)
+J_lasso_sk = clf_lasso. coef_. reshape(L, L)
+fig = plt. figure(figsize= (20 , 14 ))
+im = plt. imshow(J_lasso_sk, ** cmap_args)
+plt. title("Lasso from Scikit-learn" , fontsize=18 )
+plt. xticks(fontsize=18 )
+plt. yticks(fontsize=18 )
+cb = fig. colorbar(im)
+cb. ax. set_yticklabels(cb. ax. get_yticklabels(), fontsize=18 )
+
+plt. show()
+
+It is quite striking how LASSO breaks the symmetry of the coupling
+constant as opposed to ridge and OLS. We get a sparse solution with
+\( J_{j, j + 1} = -1 \).
+
+
+
Performance of the different models
+
+
+In order to judge which model performs best at varying values of \( \lambda \) (for ridge and LASSO) we compute \( R^2 \) which is given by
+
+$$
+\begin{align}
+ R^2 = 1 - \frac{(y - \hat{y})^2}{(y - \bar{y})^2},
+\label{_auto34}
+\end{align}
+$$
+
+where \( y \) is a vector with the true values of the energy, \( \hat{y} \) is the predicted values of \( y \) from the models and \( \bar{y} \) is the mean of \( \hat{y} \).
+
+
+
+
+
def r_squared (y, y_hat):
+ return 1 - np. sum((y - y_hat) ** 2 ) / np. sum((y - np. mean(y_hat)) ** 2 )
+
+This is the same metric used by Scikit-learn for their regression models when scoring.
+
+
+
+
= clf. predict(X_test)
+r_test = r_squared(y_test, y_hat)
+sk_r_test = clf. score(X_test, y_test)
+
+assert abs (r_test - sk_r_test) < 1e-2
+
+
Performance as function of the regularization parameter
+
+
+We see how the different models perform for a different set of values for \( \lambda \).
+
+
+
+
+
= np. logspace(-4 , 5 , 10 )
+
+train_errors = {
+ "ols_own" : np. zeros(lambdas. size),
+ "ols_sk" : np. zeros(lambdas. size),
+ "ridge_own" : np. zeros(lambdas. size),
+ "ridge_sk" : np. zeros(lambdas. size),
+ "lasso_sk" : np. zeros(lambdas. size)
+}
+
+test_errors = {
+ "ols_own" : np. zeros(lambdas. size),
+ "ols_sk" : np. zeros(lambdas. size),
+ "ridge_own" : np. zeros(lambdas. size),
+ "ridge_sk" : np. zeros(lambdas. size),
+ "lasso_sk" : np. zeros(lambdas. size)
+}
+
+plot_counter = 1
+
+fig = plt. figure(figsize= (32 , 54 ))
+
+for i, _lambda in enumerate (tqdm. tqdm(lambdas)):
+ omega = get_ols_weights(X_train_own, y_train)
+ y_hat_train = X_train_own @ omega
+ y_hat_test = X_test_own @ omega
+
+ train_errors["ols_own" ][i] = r_squared(y_train, y_hat_train)
+ test_errors["ols_own" ][i] = r_squared(y_test, y_hat_test)
+
+ plt. subplot(10 , 5 , plot_counter)
+ plt. imshow(omega[1 :]. reshape(L, L), ** cmap_args)
+ plt. title("Home made OLS" )
+ plot_counter += 1
+
+ omega = get_ridge_weights(X_train_own, y_train, _lambda)
+ y_hat_train = X_train_own @ omega
+ y_hat_test = X_test_own @ omega
+
+ train_errors["ridge_own" ][i] = r_squared(y_train, y_hat_train)
+ test_errors["ridge_own" ][i] = r_squared(y_test, y_hat_test)
+
+ plt. subplot(10 , 5 , plot_counter)
+ plt. imshow(omega[1 :]. reshape(L, L), ** cmap_args)
+ plt. title(r"Home made ridge, $\lambda = %.4f $" % _lambda)
+ plot_counter += 1
+
+ for key, method in zip (
+ ["ols_sk" , "ridge_sk" , "lasso_sk" ],
+ [skl. LinearRegression(), skl. Ridge(alpha= _lambda), skl. Lasso(alpha= _lambda)]
+ ):
+ method = method. fit(X_train, y_train)
+
+ train_errors[key][i] = method. score(X_train, y_train)
+ test_errors[key][i] = method. score(X_test, y_test)
+
+ omega = method. coef_. reshape(L, L)
+
+ plt. subplot(10 , 5 , plot_counter)
+ plt. imshow(omega, ** cmap_args)
+ plt. title(r" %s , $\lambda = %.4f $" % (key, _lambda))
+ plot_counter += 1
+
+plt. show()
+
+We can see that LASSO quite fast reaches a good solution for low
+values of \( \lambda \), but will "wither" when we increase \( \lambda \) too
+much. Ridge is more stable over a larger range of values for
+\( \lambda \), but eventually also fades away.
+
+
+
Finding the optimal value of \( \lambda \)
+
+
+To determine which value of \( \lambda \) is best we plot the accuracy of
+the models when predicting the training and the testing set. We expect
+the accuracy of the training set to be quite good, but if the accuracy
+of the testing set is much lower this tells us that we might be
+subject to an overfit model. The ideal scenario is an accuracy on the
+testing set that is close to the accuracy of the training set.
+
+
+
+
+
= plt. figure(figsize= (20 , 14 ))
+
+colors = {
+ "ols_own" : "b" ,
+ "ridge_own" : "g" ,
+ "ols_sk" : "r" ,
+ "ridge_sk" : "y" ,
+ "lasso_sk" : "c"
+}
+
+for key in train_errors:
+ plt. semilogx(
+ lambdas,
+ train_errors[key],
+ colors[key],
+ label= "Train {0}" . format(key),
+ linewidth=4.0
+ )
+
+for key in test_errors:
+ plt. semilogx(
+ lambdas,
+ test_errors[key],
+ colors[key] + "--" ,
+ label= "Test {0}" . format(key),
+ linewidth=4.0
+ )
+#plt.semilogx(lambdas, train_errors["ols_own"], label="Train (OLS own)")
+#plt.semilogx(lambdas, test_errors["ols_own"], label="Test (OLS own)")
+
+plt. legend(loc= "best" , fontsize=18 )
+plt. xlabel(r"$\lambda$" , fontsize=18 )
+plt. ylabel(r"$R^2$" , fontsize=18 )
+plt. tick_params(labelsize=18 )
+plt. show()
+
+From the above figure we can see that LASSO with \( \lambda = 10^{-2} \)
+achieve a very good accuracy on the test set. This by far surpases the
+other models for all values of \( \lambda \).
+
diff --git a/doc/pub/Regression/ipynb/Regression.ipynb b/doc/pub/Regression/ipynb/Regression.ipynb
index c3f5a08af..907fa601a 100644
--- a/doc/pub/Regression/ipynb/Regression.ipynb
+++ b/doc/pub/Regression/ipynb/Regression.ipynb
@@ -10,7 +10,7 @@
" \n",
"**Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University\n",
"\n",
- "Date: **Sep 25, 2018**\n",
+ "Date: **Oct 11, 2018**\n",
"\n",
"Copyright 1999-2018, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license\n",
"\n",
@@ -4505,7 +4505,903 @@
"estimate for our model should not vary too much between training\n",
"sets. However, if a method has high variance then small changes in\n",
"the training data can result in large changes in the model. In general, more\n",
- "flexible statistical methods have higher variance."
+ "flexible statistical methods have higher variance.\n",
+ "\n",
+ "\n",
+ "## The one-dimensional Ising model, project 2\n",
+ "\n",
+ "The one-dimensional Ising model with nearest neighbor interaction, no external field and a constant coupling constant $J$ is given by"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "
\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ " H = -J \\sum_{k}^L s_k s_{k + 1},\n",
+ "\\label{_auto19} \\tag{30}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where $s_i \\in \\{-1, 1\\}$ and $s_{N + 1} = s_1$. The number of spins in the system is determined by $L$. For the low temperature limit there is no phase transition.\n",
+ "\n",
+ "We will look at a system of $L = 40$ spins with a coupling constant of $J = 1$. To get enough training data we will generate 10000 states with their respective energies."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "L = 40\n",
+ "n = int(1e4)\n",
+ "\n",
+ "spins = np.random.choice([-1, 1], size=(n, L))\n",
+ "J = 1.0\n",
+ "\n",
+ "energies = np.zeros(n)\n",
+ "\n",
+ "for i in range(n):\n",
+ " energies[i] = - J * np.dot(spins[i], np.roll(spins[i], 1))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example: The one-dimensional Ising model\n",
+ "\n",
+ "\n",
+ "Here we use linear (ordinary least squares), ridge and LASSO\n",
+ "regression to predict the energy in the nearest neighbor\n",
+ "one-dimensional Ising model on a ring, i.e., the endpoints wrap\n",
+ "around. We will use the linear regression models to fit a value for\n",
+ "the coupling constant to achieve this."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
+ "import seaborn as sns\n",
+ "import scipy.linalg as scl\n",
+ "from sklearn.model_selection import train_test_split\n",
+ "import sklearn.linear_model as skl\n",
+ "import tqdm\n",
+ "sns.set(color_codes=True)\n",
+ "cmap_args=dict(vmin=-1., vmax=1., cmap='seismic')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Reformulating the problem to suit regression\n",
+ "\n",
+ "A more general form for the one-dimensional Ising model is"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "
\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ " H = - \\sum_j^L \\sum_k^L s_j s_k J_{jk}.\n",
+ "\\label{_auto20} \\tag{31}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Here we allow for interactions beyond the nearest neighbors and a more\n",
+ "adaptive coupling matrix. This latter expression can be formulated as\n",
+ "a matrix-product on the form"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "
\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ " H = X J,\n",
+ "\\label{_auto21} \\tag{32}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where $X_{jk} = s_j s_k$ and $J$ is the matrix consisting of the\n",
+ "elements $-J_{jk}$. This form of writing the energy fits perfectly\n",
+ "with the form utilized in linear regression, viz."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "
\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ " y = X\\omega + \\epsilon,\n",
+ "\\label{_auto22} \\tag{33}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "X = np.zeros((n, L ** 2))\n",
+ "for i in range(n):\n",
+ " X[i] = np.outer(spins[i], spins[i]).ravel()\n",
+ "y = energies\n",
+ "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.96)\n",
+ "\n",
+ "X_train_own = np.concatenate(\n",
+ " (np.ones(len(X_train))[:, np.newaxis], X_train),\n",
+ " axis=1\n",
+ ")\n",
+ "\n",
+ "X_test_own = np.concatenate(\n",
+ " (np.ones(len(X_test))[:, np.newaxis], X_test),\n",
+ " axis=1\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Linear regression\n",
+ "\n",
+ "The problem at hand is to try to fit the equation"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "
\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ " y = f(x) + \\epsilon,\n",
+ "\\label{_auto23} \\tag{34}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where $f(x)$ is some unknown function of the data $x$ and $\\epsilon$\n",
+ "is normally distributed with mean zero noise with standard deviation\n",
+ "$\\sigma_{\\epsilon}$. Our job is to try to find a predictor which\n",
+ "estimates the function $f(x)$. In linear regression we assume that we\n",
+ "can formulate the problem as"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "
\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ " y = X\\omega + \\epsilon,\n",
+ "\\label{_auto24} \\tag{35}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where $X$ and $\\omega$ are now matrices. Our job at hand is now to\n",
+ "find a **cost function** $C$, which we wish to minimize in order to find\n",
+ "the best estimate of $\\omega$.\n",
+ "\n",
+ "## Ordinary least squares\n",
+ "\n",
+ "In the ordinary least squares method we choose the cost function"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "
\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ " C(X, \\omega) = ||X\\omega - y||^2\n",
+ " = (X\\omega - y)^T(X\\omega - y)\n",
+ "\\label{_auto25} \\tag{36}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We then find the extremal point of $C$ by taking the derivative with respect to $\\omega$ and setting it to zero, i.e.,"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "
\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ " \\dfrac{\\mathrm{d}C}{\\mathrm{d}\\omega}\n",
+ " = 0.\n",
+ "\\label{_auto26} \\tag{37}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This yields the expression for $\\omega$ to be"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "
\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ " \\omega = \\frac{X^T y}{X^T X},\n",
+ "\\label{_auto27} \\tag{38}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "which immediately imposes some requirements on $X$ as there must exist\n",
+ "an inverse of $X^T X$. If the expression we are modelling contains an\n",
+ "intercept, i.e., a constant expression we must make sure that the\n",
+ "first column of $X$ consists of $1$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "def get_ols_weights_naive(x: np.ndarray, y: np.ndarray) -> np.ndarray:\n",
+ " return scl.inv(x.T @ x) @ (x.T @ y)\n",
+ "omega = get_ols_weights_naive(X_train_own, y_train)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Singular Value decomposition\n",
+ "Doing the inversion directly turns out to be a bad idea as the matrix\n",
+ "$X^TX$ is singular. An alternative approach is to use the **singular\n",
+ "value decomposition**. Using the definition of the Moore-Penrose\n",
+ "pseudoinverse we can write the equation for $\\omega$ as"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "
\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ " \\omega = X^{+}y,\n",
+ "\\label{_auto28} \\tag{39}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where the pseudoinverse of $X$ is given by"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "
\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ " X^{+} = \\frac{X^T}{X^T X}.\n",
+ "\\label{_auto29} \\tag{40}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Using singular value decomposition we have that $X = U\\Sigma V^T$,\n",
+ "where $X^{+} = V\\Sigma^{+} U^T$. This reduces the equation for\n",
+ "$\\omega$ to"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "
\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ " \\omega = V\\Sigma^{+} U^T y.\n",
+ "\\label{_auto30} \\tag{41}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Note that solving this equation by actually doing the pseudoinverse\n",
+ "(which is what we will do) is not a good idea as this operation scales\n",
+ "as $\\mathcal{O}(n^3)$, where $n$ is the number of elements in a\n",
+ "general matrix. Instead, doing $QR$-factorization and solving the\n",
+ "linear system as an equation would reduce this down to\n",
+ "$\\mathcal{O}(n^2)$ operations."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "def get_ols_weights(x: np.ndarray, y: np.ndarray) -> np.ndarray:\n",
+ " u, s, v = scl.svd(x)\n",
+ " return v.T @ scl.pinv(scl.diagsvd(s, u.shape[0], v.shape[0])) @ u.T @ y"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Before passing in the data to the function we append a column with ones to the training data."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "omega = get_ols_weights(X_train_own,y_train)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Fitting with scikit-learn\n",
+ "\n",
+ "Next we fit a `LinearRegression`-model from Scikit-learn for comparison."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 22,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "clf = skl.LinearRegression().fit(X_train, y_train)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Extracting the $J$-matrix from both our own method and the Scikit-learn model where we make sure to remove the intercept."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 23,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "J_own = omega[1:].reshape(L, L)\n",
+ "J_sk = clf.coef_.reshape(L, L)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "A way of looking at the coefficients in $J$ is to plot the matrices as images."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 24,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "fig = plt.figure(figsize=(20, 14))\n",
+ "im = plt.imshow(J_own, **cmap_args)\n",
+ "plt.title(\"Home-made OLS\", fontsize=18)\n",
+ "plt.xticks(fontsize=18)\n",
+ "plt.yticks(fontsize=18)\n",
+ "cb = fig.colorbar(im)\n",
+ "cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)\n",
+ "\n",
+ "fig = plt.figure(figsize=(20, 14))\n",
+ "im = plt.imshow(J_sk, **cmap_args)\n",
+ "plt.title(\"LinearRegression from Scikit-learn\", fontsize=18)\n",
+ "plt.xticks(fontsize=18)\n",
+ "plt.yticks(fontsize=18)\n",
+ "cb = fig.colorbar(im)\n",
+ "cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We can see that our model for the least squares method performes close\n",
+ "to the benchmark from Scikit-learn. It is interesting to note that OLS\n",
+ "considers both $J_{j, j + 1} = -0.5$ and $J_{j, j - 1} = -0.5$ as\n",
+ "valid matrix elements for $J$.\n",
+ "\n",
+ "## Ridge regression\n",
+ "\n",
+ "Having explored the ordinary least squares we move on to ridge\n",
+ "regression. In ridge regression we include a **regularizer**. This\n",
+ "involves a new cost function which leads to a new estimate for the\n",
+ "weights $\\omega$. This results in a penalized regression problem. The\n",
+ "cost function is given by"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "
\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ " C(X, \\omega; \\lambda) = ||X\\omega - y||^2 + \\lambda ||\\omega||^2\n",
+ " = (X\\omega - y)^T(X\\omega - y) + \\lambda \\omega^T\\omega.\n",
+ "\\label{_auto31} \\tag{42}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Finding the extremum of this function yields the weights"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "
\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ " \\omega(\\lambda) = \\frac{X^Ty}{X^TX + \\lambda} \\to \\frac{\\omega_{\\text{LS}}}{1 + \\lambda},\n",
+ "\\label{_auto32} \\tag{43}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where $\\omega_{\\text{LS}}$ is the weights from ordinary least\n",
+ "squares. The last assumption assumes that $X$ is orthogonal, which it\n",
+ "is not. We will therefore resort to solving the equation as it stands\n",
+ "on the left hand side."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 25,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "def get_ridge_weights(x: np.ndarray, y: np.ndarray, _lambda: float) -> np.ndarray:\n",
+ " return x.T @ y @ scl.inv(\n",
+ " x.T @ x + np.eye(x.shape[1], x.shape[1]) * _lambda\n",
+ " )\n",
+ "lambda = 0.1\n",
+ "omega_ridge = get_ridge_weights(X_train_own, y_train, np.array([_lambda]))\n",
+ "clf_ridge = skl.Ridge(alpha=_lambda).fit(X_train, y_train)\n",
+ "J_ridge_own = omega_ridge[1:].reshape(L, L)\n",
+ "J_ridge_sk = clf_ridge.coef_.reshape(L, L)\n",
+ "fig = plt.figure(figsize=(20, 14))\n",
+ "im = plt.imshow(J_ridge_own, **cmap_args)\n",
+ "plt.title(\"Home-made ridge regression\", fontsize=18)\n",
+ "plt.xticks(fontsize=18)\n",
+ "plt.yticks(fontsize=18)\n",
+ "cb = fig.colorbar(im)\n",
+ "cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)\n",
+ "\n",
+ "fig = plt.figure(figsize=(20, 14))\n",
+ "im = plt.imshow(J_ridge_sk, **cmap_args)\n",
+ "plt.title(\"Ridge from Scikit-learn\", fontsize=18)\n",
+ "plt.xticks(fontsize=18)\n",
+ "plt.yticks(fontsize=18)\n",
+ "cb = fig.colorbar(im)\n",
+ "cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)\n",
+ "\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## LASSO regression\n",
+ "\n",
+ "In the **Least Absolute Shrinkage and Selection Operator** (LASSO)-method we get a third cost function."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "
\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ " C(X, \\omega; \\lambda) =\n",
+ " ||X\\omega - y||^2 + \\lambda ||\\omega||\n",
+ " = (X\\omega - y)^T(X\\omega - y) + \\lambda \\sqrt{\\omega^T\\omega}.\n",
+ "\\label{_auto33} \\tag{44}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Finding the extremal point of this cost function is not so straight-forward as in least squares and ridge. We will therefore rely solely on the function ``Lasso`` from Scikit-learn."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "clf_lasso = skl.Lasso(alpha=_lambda).fit(X_train, y_train)\n",
+ "J_lasso_sk = clf_lasso.coef_.reshape(L, L)\n",
+ "fig = plt.figure(figsize=(20, 14))\n",
+ "im = plt.imshow(J_lasso_sk, **cmap_args)\n",
+ "plt.title(\"Lasso from Scikit-learn\", fontsize=18)\n",
+ "plt.xticks(fontsize=18)\n",
+ "plt.yticks(fontsize=18)\n",
+ "cb = fig.colorbar(im)\n",
+ "cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)\n",
+ "\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "It is quite striking how LASSO breaks the symmetry of the coupling\n",
+ "constant as opposed to ridge and OLS. We get a sparse solution with\n",
+ "$J_{j, j + 1} = -1$.\n",
+ "\n",
+ "\n",
+ "## Performance of the different models\n",
+ "\n",
+ "In order to judge which model performs best at varying values of $\\lambda$ (for ridge and LASSO) we compute $R^2$ which is given by"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "
\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ " R^2 = 1 - \\frac{(y - \\hat{y})^2}{(y - \\bar{y})^2},\n",
+ "\\label{_auto34} \\tag{45}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where $y$ is a vector with the true values of the energy, $\\hat{y}$ is the predicted values of $y$ from the models and $\\bar{y}$ is the mean of $\\hat{y}$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 27,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "def r_squared(y, y_hat):\n",
+ " return 1 - np.sum((y - y_hat) ** 2) / np.sum((y - np.mean(y_hat)) ** 2)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This is the same metric used by Scikit-learn for their regression models when scoring."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 28,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "y_hat = clf.predict(X_test)\n",
+ "r_test = r_squared(y_test, y_hat)\n",
+ "sk_r_test = clf.score(X_test, y_test)\n",
+ "\n",
+ "assert abs(r_test - sk_r_test) < 1e-2"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Performance as function of the regularization parameter\n",
+ "\n",
+ "We see how the different models perform for a different set of values for $\\lambda$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 29,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "lambdas = np.logspace(-4, 5, 10)\n",
+ "\n",
+ "train_errors = {\n",
+ " \"ols_own\": np.zeros(lambdas.size),\n",
+ " \"ols_sk\": np.zeros(lambdas.size),\n",
+ " \"ridge_own\": np.zeros(lambdas.size),\n",
+ " \"ridge_sk\": np.zeros(lambdas.size),\n",
+ " \"lasso_sk\": np.zeros(lambdas.size)\n",
+ "}\n",
+ "\n",
+ "test_errors = {\n",
+ " \"ols_own\": np.zeros(lambdas.size),\n",
+ " \"ols_sk\": np.zeros(lambdas.size),\n",
+ " \"ridge_own\": np.zeros(lambdas.size),\n",
+ " \"ridge_sk\": np.zeros(lambdas.size),\n",
+ " \"lasso_sk\": np.zeros(lambdas.size)\n",
+ "}\n",
+ "\n",
+ "plot_counter = 1\n",
+ "\n",
+ "fig = plt.figure(figsize=(32, 54))\n",
+ "\n",
+ "for i, _lambda in enumerate(tqdm.tqdm(lambdas)):\n",
+ " omega = get_ols_weights(X_train_own, y_train)\n",
+ " y_hat_train = X_train_own @ omega\n",
+ " y_hat_test = X_test_own @ omega\n",
+ "\n",
+ " train_errors[\"ols_own\"][i] = r_squared(y_train, y_hat_train)\n",
+ " test_errors[\"ols_own\"][i] = r_squared(y_test, y_hat_test)\n",
+ "\n",
+ " plt.subplot(10, 5, plot_counter)\n",
+ " plt.imshow(omega[1:].reshape(L, L), **cmap_args)\n",
+ " plt.title(\"Home made OLS\")\n",
+ " plot_counter += 1\n",
+ "\n",
+ " omega = get_ridge_weights(X_train_own, y_train, _lambda)\n",
+ " y_hat_train = X_train_own @ omega\n",
+ " y_hat_test = X_test_own @ omega\n",
+ "\n",
+ " train_errors[\"ridge_own\"][i] = r_squared(y_train, y_hat_train)\n",
+ " test_errors[\"ridge_own\"][i] = r_squared(y_test, y_hat_test)\n",
+ "\n",
+ " plt.subplot(10, 5, plot_counter)\n",
+ " plt.imshow(omega[1:].reshape(L, L), **cmap_args)\n",
+ " plt.title(r\"Home made ridge, $\\lambda = %.4f$\" % _lambda)\n",
+ " plot_counter += 1\n",
+ "\n",
+ " for key, method in zip(\n",
+ " [\"ols_sk\", \"ridge_sk\", \"lasso_sk\"],\n",
+ " [skl.LinearRegression(), skl.Ridge(alpha=_lambda), skl.Lasso(alpha=_lambda)]\n",
+ " ):\n",
+ " method = method.fit(X_train, y_train)\n",
+ "\n",
+ " train_errors[key][i] = method.score(X_train, y_train)\n",
+ " test_errors[key][i] = method.score(X_test, y_test)\n",
+ "\n",
+ " omega = method.coef_.reshape(L, L)\n",
+ "\n",
+ " plt.subplot(10, 5, plot_counter)\n",
+ " plt.imshow(omega, **cmap_args)\n",
+ " plt.title(r\"%s, $\\lambda = %.4f$\" % (key, _lambda))\n",
+ " plot_counter += 1\n",
+ "\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We can see that LASSO quite fast reaches a good solution for low\n",
+ "values of $\\lambda$, but will \"wither\" when we increase $\\lambda$ too\n",
+ "much. Ridge is more stable over a larger range of values for\n",
+ "$\\lambda$, but eventually also fades away.\n",
+ "\n",
+ "## Finding the optimal value of $\\lambda$\n",
+ "\n",
+ "To determine which value of $\\lambda$ is best we plot the accuracy of\n",
+ "the models when predicting the training and the testing set. We expect\n",
+ "the accuracy of the training set to be quite good, but if the accuracy\n",
+ "of the testing set is much lower this tells us that we might be\n",
+ "subject to an overfit model. The ideal scenario is an accuracy on the\n",
+ "testing set that is close to the accuracy of the training set."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 30,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "fig = plt.figure(figsize=(20, 14))\n",
+ "\n",
+ "colors = {\n",
+ " \"ols_own\": \"b\",\n",
+ " \"ridge_own\": \"g\",\n",
+ " \"ols_sk\": \"r\",\n",
+ " \"ridge_sk\": \"y\",\n",
+ " \"lasso_sk\": \"c\"\n",
+ "}\n",
+ "\n",
+ "for key in train_errors:\n",
+ " plt.semilogx(\n",
+ " lambdas,\n",
+ " train_errors[key],\n",
+ " colors[key],\n",
+ " label=\"Train {0}\".format(key),\n",
+ " linewidth=4.0\n",
+ " )\n",
+ "\n",
+ "for key in test_errors:\n",
+ " plt.semilogx(\n",
+ " lambdas,\n",
+ " test_errors[key],\n",
+ " colors[key] + \"--\",\n",
+ " label=\"Test {0}\".format(key),\n",
+ " linewidth=4.0\n",
+ " )\n",
+ "#plt.semilogx(lambdas, train_errors[\"ols_own\"], label=\"Train (OLS own)\")\n",
+ "#plt.semilogx(lambdas, test_errors[\"ols_own\"], label=\"Test (OLS own)\")\n",
+ "\n",
+ "plt.legend(loc=\"best\", fontsize=18)\n",
+ "plt.xlabel(r\"$\\lambda$\", fontsize=18)\n",
+ "plt.ylabel(r\"$R^2$\", fontsize=18)\n",
+ "plt.tick_params(labelsize=18)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "From the above figure we can see that LASSO with $\\lambda = 10^{-2}$\n",
+ "achieve a very good accuracy on the test set. This by far surpases the\n",
+ "other models for all values of $\\lambda$."
]
}
],
diff --git a/doc/pub/Regression/ipynb/ipynb-Regression-src.tar.gz b/doc/pub/Regression/ipynb/ipynb-Regression-src.tar.gz
index 424c6b90f..0a0242504 100644
Binary files a/doc/pub/Regression/ipynb/ipynb-Regression-src.tar.gz and b/doc/pub/Regression/ipynb/ipynb-Regression-src.tar.gz differ
diff --git a/doc/pub/Regression/pdf/Regression-beamer-handouts2x3.pdf b/doc/pub/Regression/pdf/Regression-beamer-handouts2x3.pdf
index c84b672da..b80d0630e 100644
Binary files a/doc/pub/Regression/pdf/Regression-beamer-handouts2x3.pdf and b/doc/pub/Regression/pdf/Regression-beamer-handouts2x3.pdf differ
diff --git a/doc/pub/Regression/pdf/Regression-beamer.pdf b/doc/pub/Regression/pdf/Regression-beamer.pdf
index b40b486fc..62dc14085 100644
Binary files a/doc/pub/Regression/pdf/Regression-beamer.pdf and b/doc/pub/Regression/pdf/Regression-beamer.pdf differ
diff --git a/doc/pub/Regression/pdf/Regression-minted.pdf b/doc/pub/Regression/pdf/Regression-minted.pdf
index 1819d5d39..c3db8b023 100644
Binary files a/doc/pub/Regression/pdf/Regression-minted.pdf and b/doc/pub/Regression/pdf/Regression-minted.pdf differ
diff --git a/doc/src/Regression/Regression.do.txt b/doc/src/Regression/Regression.do.txt
index 085e76228..2b8d3974a 100644
--- a/doc/src/Regression/Regression.do.txt
+++ b/doc/src/Regression/Regression.do.txt
@@ -2880,3 +2880,509 @@ sets. However, if a method has high variance then small changes in
the training data can result in large changes in the model. In general, more
flexible statistical methods have higher variance.
+
+!split
+===== The one-dimensional Ising model, project 2 =====
+
+The one-dimensional Ising model with nearest neighbor interaction, no external field and a constant coupling constant $J$ is given by
+
+!bt
+\begin{align}
+ H = -J \sum_{k}^L s_k s_{k + 1},
+\end{align}
+!et
+where $s_i \in \{-1, 1\}$ and $s_{N + 1} = s_1$. The number of spins in the system is determined by $L$. For the low temperature limit there is no phase transition.
+
+We will look at a system of $L = 40$ spins with a coupling constant of $J = 1$. To get enough training data we will generate 10000 states with their respective energies.
+
+
+!bc pycod
+L = 40
+n = int(1e4)
+
+spins = np.random.choice([-1, 1], size=(n, L))
+J = 1.0
+
+energies = np.zeros(n)
+
+for i in range(n):
+ energies[i] = - J * np.dot(spins[i], np.roll(spins[i], 1))
+!ec
+
+
+
+!split
+===== Example: The one-dimensional Ising model =====
+
+
+Here we use linear (ordinary least squares), ridge and LASSO
+regression to predict the energy in the nearest neighbor
+one-dimensional Ising model on a ring, i.e., the endpoints wrap
+around. We will use the linear regression models to fit a value for
+the coupling constant to achieve this.
+!bc pycod
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+import seaborn as sns
+import scipy.linalg as scl
+from sklearn.model_selection import train_test_split
+import sklearn.linear_model as skl
+import tqdm
+sns.set(color_codes=True)
+cmap_args=dict(vmin=-1., vmax=1., cmap='seismic')
+!ec
+
+!split
+===== Reformulating the problem to suit regression =====
+
+A more general form for the one-dimensional Ising model is
+
+!bt
+\begin{align}
+ H = - \sum_j^L \sum_k^L s_j s_k J_{jk}.
+\end{align}
+!et
+
+Here we allow for interactions beyond the nearest neighbors and a more
+adaptive coupling matrix. This latter expression can be formulated as
+a matrix-product on the form
+!bt
+\begin{align}
+ H = X J,
+\end{align}
+!et
+
+where $X_{jk} = s_j s_k$ and $J$ is the matrix consisting of the
+elements $-J_{jk}$. This form of writing the energy fits perfectly
+with the form utilized in linear regression, viz.
+!bt
+\begin{align}
+ y = X\omega + \epsilon,
+\end{align}
+!et
+
+
+
+!bc pycod
+X = np.zeros((n, L ** 2))
+for i in range(n):
+ X[i] = np.outer(spins[i], spins[i]).ravel()
+y = energies
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.96)
+
+X_train_own = np.concatenate(
+ (np.ones(len(X_train))[:, np.newaxis], X_train),
+ axis=1
+)
+
+X_test_own = np.concatenate(
+ (np.ones(len(X_test))[:, np.newaxis], X_test),
+ axis=1
+)
+!ec
+
+!split
+===== Linear regression =====
+
+The problem at hand is to try to fit the equation
+!bt
+\begin{align}
+ y = f(x) + \epsilon,
+\end{align}
+!et
+
+where $f(x)$ is some unknown function of the data $x$ and $\epsilon$
+is normally distributed with mean zero noise with standard deviation
+$\sigma_{\epsilon}$. Our job is to try to find a predictor which
+estimates the function $f(x)$. In linear regression we assume that we
+can formulate the problem as
+
+!bt
+\begin{align}
+ y = X\omega + \epsilon,
+\end{align}
+!et
+
+where $X$ and $\omega$ are now matrices. Our job at hand is now to
+find a _cost function_ $C$, which we wish to minimize in order to find
+the best estimate of $\omega$.
+
+!split
+===== Ordinary least squares =====
+
+In the ordinary least squares method we choose the cost function
+!bt
+\begin{align}
+ C(X, \omega) = ||X\omega - y||^2
+ = (X\omega - y)^T(X\omega - y)
+\end{align}
+!et
+We then find the extremal point of $C$ by taking the derivative with respect to $\omega$ and setting it to zero, i.e.,
+
+!bt
+\begin{align}
+ \dfrac{\mathrm{d}C}{\mathrm{d}\omega}
+ = 0.
+\end{align}
+!et
+This yields the expression for $\omega$ to be
+!bt
+\begin{align}
+ \omega = \frac{X^T y}{X^T X},
+\end{align}
+!et
+
+which immediately imposes some requirements on $X$ as there must exist
+an inverse of $X^T X$. If the expression we are modelling contains an
+intercept, i.e., a constant expression we must make sure that the
+first column of $X$ consists of $1$.
+
+
+!bc pycod
+def get_ols_weights_naive(x: np.ndarray, y: np.ndarray) -> np.ndarray:
+ return scl.inv(x.T @ x) @ (x.T @ y)
+omega = get_ols_weights_naive(X_train_own, y_train)
+!ec
+
+
+!split
+===== Singular Value decomposition =====
+Doing the inversion directly turns out to be a bad idea as the matrix
+$X^TX$ is singular. An alternative approach is to use the _singular
+value decomposition_. Using the definition of the Moore-Penrose
+pseudoinverse we can write the equation for $\omega$ as
+
+!bt
+\begin{align}
+ \omega = X^{+}y,
+\end{align}
+!et
+where the pseudoinverse of $X$ is given by
+!bt
+\begin{align}
+ X^{+} = \frac{X^T}{X^T X}.
+\end{align}
+!et
+
+Using singular value decomposition we have that $X = U\Sigma V^T$,
+where $X^{+} = V\Sigma^{+} U^T$. This reduces the equation for
+$\omega$ to
+!bt
+\begin{align}
+ \omega = V\Sigma^{+} U^T y.
+\end{align}
+!et
+
+Note that solving this equation by actually doing the pseudoinverse
+(which is what we will do) is not a good idea as this operation scales
+as $\mathcal{O}(n^3)$, where $n$ is the number of elements in a
+general matrix. Instead, doing $QR$-factorization and solving the
+linear system as an equation would reduce this down to
+$\mathcal{O}(n^2)$ operations.
+
+
+!bc pycod
+def get_ols_weights(x: np.ndarray, y: np.ndarray) -> np.ndarray:
+ u, s, v = scl.svd(x)
+ return v.T @ scl.pinv(scl.diagsvd(s, u.shape[0], v.shape[0])) @ u.T @ y
+!ec
+Before passing in the data to the function we append a column with ones to the training data.
+
+
+!bc pycod
+omega = get_ols_weights(X_train_own,y_train)
+!ec
+
+!split
+===== Fitting with scikit-learn =====
+
+Next we fit a `LinearRegression`-model from Scikit-learn for comparison.
+
+!bc pycod
+clf = skl.LinearRegression().fit(X_train, y_train)
+!ec
+
+Extracting the $J$-matrix from both our own method and the Scikit-learn model where we make sure to remove the intercept.
+
+
+!bc pycod
+J_own = omega[1:].reshape(L, L)
+J_sk = clf.coef_.reshape(L, L)
+!ec
+
+A way of looking at the coefficients in $J$ is to plot the matrices as images.
+
+
+!bc pycod
+fig = plt.figure(figsize=(20, 14))
+im = plt.imshow(J_own, **cmap_args)
+plt.title("Home-made OLS", fontsize=18)
+plt.xticks(fontsize=18)
+plt.yticks(fontsize=18)
+cb = fig.colorbar(im)
+cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)
+
+fig = plt.figure(figsize=(20, 14))
+im = plt.imshow(J_sk, **cmap_args)
+plt.title("LinearRegression from Scikit-learn", fontsize=18)
+plt.xticks(fontsize=18)
+plt.yticks(fontsize=18)
+cb = fig.colorbar(im)
+cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)
+plt.show()
+!ec
+
+We can see that our model for the least squares method performes close
+to the benchmark from Scikit-learn. It is interesting to note that OLS
+considers both $J_{j, j + 1} = -0.5$ and $J_{j, j - 1} = -0.5$ as
+valid matrix elements for $J$.
+
+!split
+===== Ridge regression =====
+
+Having explored the ordinary least squares we move on to ridge
+regression. In ridge regression we include a _regularizer_. This
+involves a new cost function which leads to a new estimate for the
+weights $\omega$. This results in a penalized regression problem. The
+cost function is given by
+
+!bt
+\begin{align}
+ C(X, \omega; \lambda) = ||X\omega - y||^2 + \lambda ||\omega||^2
+ = (X\omega - y)^T(X\omega - y) + \lambda \omega^T\omega.
+\end{align}
+!et
+Finding the extremum of this function yields the weights
+
+!bt
+\begin{align}
+ \omega(\lambda) = \frac{X^Ty}{X^TX + \lambda} \to \frac{\omega_{\text{LS}}}{1 + \lambda},
+\end{align}
+!et
+
+where $\omega_{\text{LS}}$ is the weights from ordinary least
+squares. The last assumption assumes that $X$ is orthogonal, which it
+is not. We will therefore resort to solving the equation as it stands
+on the left hand side.
+
+
+!bc pycod
+def get_ridge_weights(x: np.ndarray, y: np.ndarray, _lambda: float) -> np.ndarray:
+ return x.T @ y @ scl.inv(
+ x.T @ x + np.eye(x.shape[1], x.shape[1]) * _lambda
+ )
+lambda = 0.1
+omega_ridge = get_ridge_weights(X_train_own, y_train, np.array([_lambda]))
+clf_ridge = skl.Ridge(alpha=_lambda).fit(X_train, y_train)
+J_ridge_own = omega_ridge[1:].reshape(L, L)
+J_ridge_sk = clf_ridge.coef_.reshape(L, L)
+fig = plt.figure(figsize=(20, 14))
+im = plt.imshow(J_ridge_own, **cmap_args)
+plt.title("Home-made ridge regression", fontsize=18)
+plt.xticks(fontsize=18)
+plt.yticks(fontsize=18)
+cb = fig.colorbar(im)
+cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)
+
+fig = plt.figure(figsize=(20, 14))
+im = plt.imshow(J_ridge_sk, **cmap_args)
+plt.title("Ridge from Scikit-learn", fontsize=18)
+plt.xticks(fontsize=18)
+plt.yticks(fontsize=18)
+cb = fig.colorbar(im)
+cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)
+
+plt.show()
+!ec
+
+!split
+===== LASSO regression =====
+
+In the _Least Absolute Shrinkage and Selection Operator_ (LASSO)-method we get a third cost function.
+!bt
+\begin{align}
+ C(X, \omega; \lambda) =
+ ||X\omega - y||^2 + \lambda ||\omega||
+ = (X\omega - y)^T(X\omega - y) + \lambda \sqrt{\omega^T\omega}.
+\end{align}
+!et
+Finding the extremal point of this cost function is not so straight-forward as in least squares and ridge. We will therefore rely solely on the function ``Lasso`` from Scikit-learn.
+
+!bc pycod
+clf_lasso = skl.Lasso(alpha=_lambda).fit(X_train, y_train)
+J_lasso_sk = clf_lasso.coef_.reshape(L, L)
+fig = plt.figure(figsize=(20, 14))
+im = plt.imshow(J_lasso_sk, **cmap_args)
+plt.title("Lasso from Scikit-learn", fontsize=18)
+plt.xticks(fontsize=18)
+plt.yticks(fontsize=18)
+cb = fig.colorbar(im)
+cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)
+
+plt.show()
+!ec
+
+It is quite striking how LASSO breaks the symmetry of the coupling
+constant as opposed to ridge and OLS. We get a sparse solution with
+$J_{j, j + 1} = -1$.
+
+
+!split
+===== Performance of the different models =====
+
+In order to judge which model performs best at varying values of $\lambda$ (for ridge and LASSO) we compute $R^2$ which is given by
+
+!bt
+\begin{align}
+ R^2 = 1 - \frac{(y - \hat{y})^2}{(y - \bar{y})^2},
+\end{align}
+!et
+where $y$ is a vector with the true values of the energy, $\hat{y}$ is the predicted values of $y$ from the models and $\bar{y}$ is the mean of $\hat{y}$.
+
+
+!bc pycod
+def r_squared(y, y_hat):
+ return 1 - np.sum((y - y_hat) ** 2) / np.sum((y - np.mean(y_hat)) ** 2)
+!ec
+
+This is the same metric used by Scikit-learn for their regression models when scoring.
+!bc pycod
+y_hat = clf.predict(X_test)
+r_test = r_squared(y_test, y_hat)
+sk_r_test = clf.score(X_test, y_test)
+
+assert abs(r_test - sk_r_test) < 1e-2
+!ec
+
+
+!split
+===== Performance as function of the regularization parameter =====
+
+We see how the different models perform for a different set of values for $\lambda$.
+
+
+!bc pycod
+lambdas = np.logspace(-4, 5, 10)
+
+train_errors = {
+ "ols_own": np.zeros(lambdas.size),
+ "ols_sk": np.zeros(lambdas.size),
+ "ridge_own": np.zeros(lambdas.size),
+ "ridge_sk": np.zeros(lambdas.size),
+ "lasso_sk": np.zeros(lambdas.size)
+}
+
+test_errors = {
+ "ols_own": np.zeros(lambdas.size),
+ "ols_sk": np.zeros(lambdas.size),
+ "ridge_own": np.zeros(lambdas.size),
+ "ridge_sk": np.zeros(lambdas.size),
+ "lasso_sk": np.zeros(lambdas.size)
+}
+
+plot_counter = 1
+
+fig = plt.figure(figsize=(32, 54))
+
+for i, _lambda in enumerate(tqdm.tqdm(lambdas)):
+ omega = get_ols_weights(X_train_own, y_train)
+ y_hat_train = X_train_own @ omega
+ y_hat_test = X_test_own @ omega
+
+ train_errors["ols_own"][i] = r_squared(y_train, y_hat_train)
+ test_errors["ols_own"][i] = r_squared(y_test, y_hat_test)
+
+ plt.subplot(10, 5, plot_counter)
+ plt.imshow(omega[1:].reshape(L, L), **cmap_args)
+ plt.title("Home made OLS")
+ plot_counter += 1
+
+ omega = get_ridge_weights(X_train_own, y_train, _lambda)
+ y_hat_train = X_train_own @ omega
+ y_hat_test = X_test_own @ omega
+
+ train_errors["ridge_own"][i] = r_squared(y_train, y_hat_train)
+ test_errors["ridge_own"][i] = r_squared(y_test, y_hat_test)
+
+ plt.subplot(10, 5, plot_counter)
+ plt.imshow(omega[1:].reshape(L, L), **cmap_args)
+ plt.title(r"Home made ridge, $\lambda = %.4f$" % _lambda)
+ plot_counter += 1
+
+ for key, method in zip(
+ ["ols_sk", "ridge_sk", "lasso_sk"],
+ [skl.LinearRegression(), skl.Ridge(alpha=_lambda), skl.Lasso(alpha=_lambda)]
+ ):
+ method = method.fit(X_train, y_train)
+
+ train_errors[key][i] = method.score(X_train, y_train)
+ test_errors[key][i] = method.score(X_test, y_test)
+
+ omega = method.coef_.reshape(L, L)
+
+ plt.subplot(10, 5, plot_counter)
+ plt.imshow(omega, **cmap_args)
+ plt.title(r"%s, $\lambda = %.4f$" % (key, _lambda))
+ plot_counter += 1
+
+plt.show()
+!ec
+
+We can see that LASSO quite fast reaches a good solution for low
+values of $\lambda$, but will "wither" when we increase $\lambda$ too
+much. Ridge is more stable over a larger range of values for
+$\lambda$, but eventually also fades away.
+
+!split
+===== Finding the optimal value of $\lambda$ =====
+
+To determine which value of $\lambda$ is best we plot the accuracy of
+the models when predicting the training and the testing set. We expect
+the accuracy of the training set to be quite good, but if the accuracy
+of the testing set is much lower this tells us that we might be
+subject to an overfit model. The ideal scenario is an accuracy on the
+testing set that is close to the accuracy of the training set.
+
+
+!bc pycod
+fig = plt.figure(figsize=(20, 14))
+
+colors = {
+ "ols_own": "b",
+ "ridge_own": "g",
+ "ols_sk": "r",
+ "ridge_sk": "y",
+ "lasso_sk": "c"
+}
+
+for key in train_errors:
+ plt.semilogx(
+ lambdas,
+ train_errors[key],
+ colors[key],
+ label="Train {0}".format(key),
+ linewidth=4.0
+ )
+
+for key in test_errors:
+ plt.semilogx(
+ lambdas,
+ test_errors[key],
+ colors[key] + "--",
+ label="Test {0}".format(key),
+ linewidth=4.0
+ )
+#plt.semilogx(lambdas, train_errors["ols_own"], label="Train (OLS own)")
+#plt.semilogx(lambdas, test_errors["ols_own"], label="Test (OLS own)")
+
+plt.legend(loc="best", fontsize=18)
+plt.xlabel(r"$\lambda$", fontsize=18)
+plt.ylabel(r"$R^2$", fontsize=18)
+plt.tick_params(labelsize=18)
+plt.show()
+!ec
+
+From the above figure we can see that LASSO with $\lambda = 10^{-2}$
+achieve a very good accuracy on the test set. This by far surpases the
+other models for all values of $\lambda$.