diff --git a/doc/pub/Regression/ipynb/Regression.ipynb b/doc/pub/Regression/ipynb/Regression.ipynb index 52985f172..4bc47575b 100644 --- a/doc/pub/Regression/ipynb/Regression.ipynb +++ b/doc/pub/Regression/ipynb/Regression.ipynb @@ -381,9 +381,7 @@ { "cell_type": "code", "execution_count": 1, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", @@ -888,9 +886,7 @@ { "cell_type": "code", "execution_count": 2, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# matrix inversion to find beta\n", @@ -909,9 +905,7 @@ { "cell_type": "code", "execution_count": 3, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "fit = np.linalg.lstsq(X, Energies, rcond =None)[0]\n", @@ -928,9 +922,7 @@ { "cell_type": "code", "execution_count": 4, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "Masses['Eapprox'] = ytilde\n", @@ -960,9 +952,7 @@ { "cell_type": "code", "execution_count": 5, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "def R2(y_data, y_model):\n", @@ -979,9 +969,7 @@ { "cell_type": "code", "execution_count": 6, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "print(R2(Energies,ytilde))" @@ -997,9 +985,7 @@ { "cell_type": "code", "execution_count": 7, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "def MSE(y_data,y_model):\n", @@ -1019,9 +1005,7 @@ { "cell_type": "code", "execution_count": 8, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "def RelativeError(y_data,y_model):\n", @@ -1398,9 +1382,7 @@ { "cell_type": "code", "execution_count": 9, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# Common imports\n", @@ -1519,9 +1501,7 @@ { "cell_type": "code", "execution_count": 10, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "import os\n", @@ -1638,9 +1618,7 @@ { "cell_type": "code", "execution_count": 11, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", @@ -1660,9 +1638,7 @@ { "cell_type": "code", "execution_count": 12, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "from sklearn.datasets import load_boston\n", @@ -1684,9 +1660,7 @@ { "cell_type": "code", "execution_count": 13, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "boston = pd.DataFrame(boston_dataset.data, columns=boston_dataset.feature_names)\n", @@ -1704,9 +1678,7 @@ { "cell_type": "code", "execution_count": 14, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# check for missing values in all the columns\n", @@ -1723,9 +1695,7 @@ { "cell_type": "code", "execution_count": 15, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# set the size of the figure\n", @@ -1746,9 +1716,7 @@ { "cell_type": "code", "execution_count": 16, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# compute the pair wise correlation for all columns \n", @@ -1768,9 +1736,7 @@ { "cell_type": "code", "execution_count": 17, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(20, 5))\n", @@ -1798,9 +1764,7 @@ { "cell_type": "code", "execution_count": 18, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "X = pd.DataFrame(np.c_[boston['LSTAT'], boston['RM']], columns = ['LSTAT','RM'])\n", @@ -1817,9 +1781,7 @@ { "cell_type": "code", "execution_count": 19, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "from sklearn.model_selection import train_test_split\n", @@ -1843,9 +1805,7 @@ { "cell_type": "code", "execution_count": 20, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "from sklearn.linear_model import LinearRegression\n", @@ -1884,9 +1844,7 @@ { "cell_type": "code", "execution_count": 21, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# plotting the y_test vs y_pred\n", @@ -1971,9 +1929,7 @@ { "cell_type": "code", "execution_count": 22, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# Common imports\n", @@ -2839,11 +2795,33 @@ }, { "cell_type": "code", - "execution_count": 23, - "metadata": { - "collapsed": false - }, - "outputs": [], + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[ 1. -1. 2.]\n", + " [ 1. 0. 1.]\n", + " [ 1. 2. -1.]\n", + " [ 1. 1. 0.]]\n", + "[[ 4. 2. 2.]\n", + " [ 2. 6. -4.]\n", + " [ 2. -4. 6.]]\n", + "[[-9.57425734e-17 8.16496581e-01 -5.77350269e-01]\n", + " [-7.07106781e-01 4.08248290e-01 5.77350269e-01]\n", + " [ 7.07106781e-01 4.08248290e-01 5.77350269e-01]]\n", + "[1.00000000e+01 6.00000000e+00 9.10898112e-32]\n", + "[[ 3.33066907e-17 -7.07106781e-01 7.07106781e-01]\n", + " [ 8.16496581e-01 4.08248290e-01 4.08248290e-01]\n", + " [ 5.77350269e-01 -5.77350269e-01 -5.77350269e-01]]\n", + "[[-3.65939208e+30 3.65939208e+30 3.65939208e+30]\n", + " [ 3.65939208e+30 -3.65939208e+30 -3.65939208e+30]\n", + " [ 3.65939208e+30 -3.65939208e+30 -3.65939208e+30]]\n" + ] + } + ], "source": [ "import numpy as np\n", "# SVD inversion\n", @@ -3198,11 +3176,20 @@ }, { "cell_type": "code", - "execution_count": 24, - "metadata": { - "collapsed": false - }, - "outputs": [], + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.06921280734901214\n", + "4.282897881229527\n", + "[[ 1.16689268 3.45320489]\n", + " [ 3.45320489 11.11096825]]\n" + ] + } + ], "source": [ "# Importing various packages\n", "import numpy as np\n", @@ -3231,11 +3218,20 @@ }, { "cell_type": "code", - "execution_count": 25, - "metadata": { - "collapsed": false - }, - "outputs": [], + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.07629059232706191\n", + "1.3689323216991471\n", + "[[1. 0.63967543]\n", + " [0.63967543 1. ]]\n" + ] + } + ], "source": [ "import numpy as np\n", "n = 100\n", @@ -3277,11 +3273,40 @@ }, { "cell_type": "code", - "execution_count": 26, - "metadata": { - "collapsed": false - }, - "outputs": [], + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[-0.55783318 -2.36855737]\n", + " [-0.48038212 -1.12123656]\n", + " [-0.90619932 -1.43972046]\n", + " [ 0.0509597 0.25649197]\n", + " [-0.00589746 1.0155279 ]\n", + " [ 0.42557278 1.62052588]\n", + " [ 0.76914274 -1.09850401]\n", + " [-0.13134882 0.66162568]\n", + " [ 0.85273632 2.51549691]\n", + " [-0.01675064 -0.04164992]]\n", + " 0 1\n", + "0 -0.557833 -2.368557\n", + "1 -0.480382 -1.121237\n", + "2 -0.906199 -1.439720\n", + "3 0.050960 0.256492\n", + "4 -0.005897 1.015528\n", + "5 0.425573 1.620526\n", + "6 0.769143 -1.098504\n", + "7 -0.131349 0.661626\n", + "8 0.852736 2.515497\n", + "9 -0.016751 -0.041650\n", + " 0 1\n", + "0 1.000000 0.657956\n", + "1 0.657956 1.000000\n" + ] + } + ], "source": [ "import numpy as np\n", "import pandas as pd\n", @@ -3309,11 +3334,49 @@ }, { "cell_type": "code", - "execution_count": 27, - "metadata": { - "collapsed": false - }, - "outputs": [], + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 0 1 2 3 4 5 6 7 \\\n", + "0 0.0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 \n", + "1 0.0 0.098274 0.088019 0.102816 0.095130 0.088038 0.095734 0.089350 \n", + "2 0.0 0.088019 0.080153 0.091357 0.085370 0.079873 0.085008 0.079956 \n", + "3 0.0 0.102816 0.091357 0.113039 0.104455 0.096480 0.108457 0.101346 \n", + "4 0.0 0.095130 0.085370 0.104455 0.097137 0.090339 0.100399 0.094305 \n", + "5 0.0 0.088038 0.079873 0.096480 0.090339 0.084640 0.092879 0.087723 \n", + "6 0.0 0.095734 0.085008 0.108457 0.100399 0.092879 0.106219 0.099527 \n", + "7 0.0 0.089350 0.079956 0.101346 0.094305 0.087723 0.099527 0.093665 \n", + "8 0.0 0.083518 0.075349 0.094820 0.088708 0.082989 0.093363 0.088259 \n", + "9 0.0 0.078178 0.071145 0.088809 0.083553 0.078629 0.087665 0.083254 \n", + "10 0.0 0.087401 0.077865 0.101048 0.093860 0.087128 0.100476 0.094478 \n", + "11 0.0 0.081937 0.073490 0.094954 0.088611 0.082657 0.094714 0.089415 \n", + "12 0.0 0.076958 0.069504 0.089383 0.083807 0.078560 0.089433 0.084768 \n", + "13 0.0 0.072416 0.065869 0.084280 0.079402 0.074802 0.084581 0.080492 \n", + "14 0.0 0.068264 0.062550 0.079595 0.075355 0.071347 0.080112 0.076548 \n", + "\n", + " 8 9 10 11 12 13 14 \n", + "0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 \n", + "1 0.083518 0.078178 0.087401 0.081937 0.076958 0.072416 0.068264 \n", + "2 0.075349 0.071145 0.077865 0.073490 0.069504 0.065869 0.062550 \n", + "3 0.094820 0.088809 0.101048 0.094954 0.089383 0.084280 0.079595 \n", + "4 0.088708 0.083553 0.093860 0.088611 0.083807 0.079402 0.075355 \n", + "5 0.082989 0.078629 0.087128 0.082657 0.078560 0.074802 0.071347 \n", + "6 0.093363 0.087665 0.100476 0.094714 0.089433 0.084581 0.080112 \n", + "7 0.088259 0.083254 0.094478 0.089415 0.084768 0.080492 0.076548 \n", + "8 0.083545 0.079175 0.088938 0.084512 0.080443 0.076693 0.073228 \n", + "9 0.079175 0.075391 0.083801 0.079956 0.076416 0.073148 0.070125 \n", + "10 0.088938 0.083801 0.096178 0.090983 0.086211 0.081816 0.077756 \n", + "11 0.084512 0.079956 0.090983 0.086383 0.082149 0.078242 0.074627 \n", + "12 0.080443 0.076416 0.086211 0.082149 0.078403 0.074940 0.071730 \n", + "13 0.076693 0.073148 0.081816 0.078242 0.074940 0.071881 0.069039 \n", + "14 0.073228 0.070125 0.077756 0.074627 0.071730 0.069039 0.066534 \n" + ] + } + ], "source": [ "# Common imports\n", "import numpy as np\n", @@ -3956,9 +4019,7 @@ { "cell_type": "code", "execution_count": 28, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# Importing various packages\n", @@ -3981,9 +4042,7 @@ { "cell_type": "code", "execution_count": 29, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", @@ -5146,9 +5205,7 @@ { "cell_type": "code", "execution_count": 30, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "from numpy import *\n", @@ -5289,9 +5346,7 @@ { "cell_type": "code", "execution_count": 31, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "from numpy import *\n", @@ -5437,9 +5492,7 @@ { "cell_type": "code", "execution_count": 32, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", @@ -5668,9 +5721,7 @@ { "cell_type": "code", "execution_count": 33, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", @@ -5739,9 +5790,7 @@ { "cell_type": "code", "execution_count": 34, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", @@ -5840,9 +5889,7 @@ { "cell_type": "code", "execution_count": 35, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "\"\"\"\n", @@ -5928,9 +5975,7 @@ { "cell_type": "code", "execution_count": 36, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# Common imports\n", @@ -6024,9 +6069,7 @@ { "cell_type": "code", "execution_count": 37, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# Common imports\n", @@ -6108,9 +6151,7 @@ { "cell_type": "code", "execution_count": 38, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", @@ -6191,9 +6232,7 @@ { "cell_type": "code", "execution_count": 39, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", @@ -6306,9 +6345,7 @@ { "cell_type": "code", "execution_count": 40, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "X = np.zeros((n, L ** 2))\n", @@ -6372,9 +6409,7 @@ { "cell_type": "code", "execution_count": 41, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "X_train_own = np.concatenate(\n", @@ -6390,9 +6425,7 @@ { "cell_type": "code", "execution_count": 42, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "def ols_inv(x: np.ndarray, y: np.ndarray) -> np.ndarray:\n", @@ -6477,9 +6510,7 @@ { "cell_type": "code", "execution_count": 43, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "def ols_svd(x: np.ndarray, y: np.ndarray) -> np.ndarray:\n", @@ -6490,9 +6521,7 @@ { "cell_type": "code", "execution_count": 44, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "beta = ols_svd(X_train_own,y_train)" @@ -6508,9 +6537,7 @@ { "cell_type": "code", "execution_count": 45, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "J = beta[1:].reshape(L, L)" @@ -6526,9 +6553,7 @@ { "cell_type": "code", "execution_count": 46, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(20, 14))\n", @@ -6594,9 +6619,7 @@ { "cell_type": "code", "execution_count": 47, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", @@ -6702,9 +6725,7 @@ { "cell_type": "code", "execution_count": 48, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "X = np.zeros((n, L ** 2))\n", @@ -6734,9 +6755,7 @@ { "cell_type": "code", "execution_count": 49, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "clf = skl.LinearRegression().fit(X_train, y_train)" @@ -6752,9 +6771,7 @@ { "cell_type": "code", "execution_count": 50, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "J_sk = clf.coef_.reshape(L, L)" @@ -6770,9 +6787,7 @@ { "cell_type": "code", "execution_count": 51, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(20, 14))\n", @@ -6828,9 +6843,7 @@ { "cell_type": "code", "execution_count": 52, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "_lambda = 0.1\n", @@ -6881,9 +6894,7 @@ { "cell_type": "code", "execution_count": 53, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "clf_lasso = skl.Lasso(alpha=_lambda).fit(X_train, y_train)\n", @@ -6917,9 +6928,7 @@ { "cell_type": "code", "execution_count": 54, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "lambdas = np.logspace(-4, 5, 10)\n", @@ -6982,9 +6991,7 @@ { "cell_type": "code", "execution_count": 55, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(20, 14))\n", @@ -7029,7 +7036,25 @@ ] } ], - "metadata": {}, + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.3" + } + }, "nbformat": 4, "nbformat_minor": 2 } diff --git a/doc/src/Regression/Regression.do.txt b/doc/src/Regression/Regression.do.txt index 33b034e21..1c0cbc1b9 100644 --- a/doc/src/Regression/Regression.do.txt +++ b/doc/src/Regression/Regression.do.txt @@ -1312,7 +1312,7 @@ This may however not the be case in general and a standard matrix inversion algorithm based on say LU, QR or Cholesky decomposition may lead to singularities. We will see examples of this below. -There is however a way to partially circumvent this problem and also gain some insight about the ordinary least squares approach. +There is however a way to partially circumvent this problem and also gain some insights about the ordinary least squares approach, and later shrinkage methods like Ridge and Lasso regressions. This is given by the _Singular Value Decomposition_ algorithm, perhaps the most powerful linear algebra algorithm. Let us look at a @@ -1519,6 +1519,56 @@ If $p > n$, then only the first $n$ columns of $\bm{V}$ are computed and $\bm{\S The $n=p$ case is obvious, we retain the full SVD. In general the economy-size SVD leads to less FLOPS and still conserving the desired accuracy. +!split +===== Codes for the SVD ===== + +!bc pycod +import numpy as np +# SVD inversion +def SVDinv(A): + ''' Takes as input a numpy matrix A and returns inv(A) based on singular value decomposition (SVD). + SVD is numerically more stable than the inversion algorithms provided by + numpy and scipy.linalg at the cost of being slower. + ''' + U, s, VT = np.linalg.svd(A) +# print('test U') +# print( (np.transpose(U) @ U - U @np.transpose(U))) +# print('test VT') +# print( (np.transpose(VT) @ VT - VT @np.transpose(VT))) + print(U) + print(s) + print(VT) + + D = np.zeros((len(U),len(VT))) + for i in range(0,len(VT)): + D[i,i]=s[i] + UT = np.transpose(U); V = np.transpose(VT); invD = np.linalg.inv(D) + return np.matmul(V,np.matmul(invD,UT)) + + +X = np.array([ [1.0, -1.0, 2.0], [1.0, 0.0, 1.0], [1.0, 2.0, -1.0], [1.0, 1.0, 0.0] ]) +print(X) +A = np.transpose(X) @ X +print(A) +# Brute force inversion of super-collinear matrix +#B = np.linalg.inv(A) +#print(B) +C = SVDinv(A) +print(C) + +!ec + +The matrix $\bm{X}$ has columns that are linearly dependent. The first +column is the row-wise sum of the other two columns. The rank of a +matrix (the column rank) is the dimension of space spanned by the +column vectors. The rank of the matrix is the number of linearly +independent columns, in this case just $2$. We see this from the +singular values when running the above code. Running the standard +inversion algorithm for matrix inversion with $\bm{X}^T\bm{X}$ results +in the program terminating due to a singular matrix. + + + !split ===== Mathematical Properties ===== @@ -1753,54 +1803,6 @@ We will come back to more interpreations after we have gone through some of the For more discussions of Ridge and Lasso regression, "Wessel van Wieringen's":"https://arxiv.org/abs/1509.09169" article is highly recommended. Similarly, "Mehta et al's article":"https://arxiv.org/abs/1803.08823" is also recommended. -!split -===== Codes for the SVD ===== - -!bc pycod -import numpy as np -# SVD inversion -def SVDinv(A): - ''' Takes as input a numpy matrix A and returns inv(A) based on singular value decomposition (SVD). - SVD is numerically more stable than the inversion algorithms provided by - numpy and scipy.linalg at the cost of being slower. - ''' - U, s, VT = np.linalg.svd(A) -# print('test U') -# print( (np.transpose(U) @ U - U @np.transpose(U))) -# print('test VT') -# print( (np.transpose(VT) @ VT - VT @np.transpose(VT))) - print(U) - print(s) - print(VT) - - D = np.zeros((len(U),len(VT))) - for i in range(0,len(VT)): - D[i,i]=s[i] - UT = np.transpose(U); V = np.transpose(VT); invD = np.linalg.inv(D) - return np.matmul(V,np.matmul(invD,UT)) - - -X = np.array([ [1.0, -1.0, 2.0], [1.0, 0.0, 1.0], [1.0, 2.0, -1.0], [1.0, 1.0, 0.0] ]) -print(X) -A = np.transpose(X) @ X -print(A) -# Brute force inversion of super-collinear matrix -#B = np.linalg.inv(A) -#print(B) -C = SVDinv(A) -print(C) - -!ec - -The matrix $\bm{X}$ has columns that are linearly dependent. The first -column is the row-wise sum of the other two columns. The rank of a -matrix (the column rank) is the dimension of space spanned by the -column vectors. The rank of the matrix is the number of linearly -independent columns, in this case just $2$. We see this from the -singular values when running the above code. Running the standard -inversion algorithm for matrix inversion with $\bm{X}^T\bm{X}$ results -in the program terminating due to a singular matrix. - !split ===== A better understanding of regularization =====