The residual is zero when we reach the minimum of the quadratic equation @@ -223,9 +208,6 @@ variance, then the matrix \( \hat{A} \), which is called the Hessian, is given by the second-derivative of the function we want to minimize. This quantity is always positive definite. -
-More details will be added here soon. -
@@ -252,7 +234,7 @@ More details will be added here soon.
+
+We denote the initial guess for \( \hat{x} \) as \( \hat{x}_0 \). +We can assume without loss of generality that +$$ +\begin{equation*} +\hat{x}_0=0, +\end{equation*} +$$ - -
#include <cmath>
-#include <iostream>
-#include <fstream>
-#include <iomanip>
-#include "vectormatrixclass.h"
-using namespace std;
-// Main function begins here
-int main(int argc, char * argv[]){
- int dim = 2;
- Vector x(dim),xsd(dim), b(dim),x0(dim);
- Matrix A(dim,dim);
-
- // Set our initial guess
- x0(0) = x0(1) = 0;
- // Set the matrix
- A(0,0) = 3; A(1,0) = 2; A(0,1) = 2; A(1,1) = 6;
- b(0) = 2; b(1) = -8;
- cout << "The Matrix A that we are using: " << endl;
- A.Print();
- cout << endl;
- x = ConjugateGradient(A,b,x0);
- xsd = SteepestDescent(A,b,x0);
- cout << "The approximate solution using Conjugate Gradient is: " << endl;
- x.Print();
- cout << endl;
- cout << "The approximate solution using Steepest Descent is: " << endl;
- xsd.Print();
- cout << endl;
-}
--
@@ -274,7 +237,7 @@ MathJax.Hub.Config({
-
+One can show that the solution \( \hat{x} \) is also the unique minimizer of the quadratic form +$$ +\begin{equation*} + f(\hat{x}) = \frac{1}{2}\hat{x}^T\hat{A}\hat{x} - \hat{x}^T \hat{x} , \quad \hat{x}\in\mathbf{R}^n. +\end{equation*} +$$ + +This suggests taking the first basis vector \( \hat{p}_1 \) +to be the gradient of \( f \) at \( \hat{x}=\hat{x}_0 \), +which equals +$$ +\begin{equation*} +\hat{A}\hat{x}_0-\hat{b}, +\end{equation*} +$$ + +and +\( \hat{x}_0=0 \) it is equal \( -\hat{b} \). - -
Vector SteepestDescent(Matrix A, Vector b, Vector x0){
- int IterMax, i;
- int dim = x0.Dimension();
- const double tolerance = 1.0e-14;
- Vector x(dim),f(dim),z(dim);
- double c,alpha,d;
- IterMax = 30;
- x = x0;
- f = A*x-b;
- i = 0;
- while (i <= IterMax){
- z = A*f;
- c = dot(f,f);
- alpha = c/dot(f,z);
- x = x - alpha*f;
- f = A*x-b;
- if(sqrt(dot(f,f)) < tolerance) break;
- i++;
- }
- return x;
-}
-
- + -
-We will use linear regression as a case study for the gradient descent -methods. Linear regression is a great test case for the gradient -descent methods discussed in the lectures since it has several -desirable properties such as: - -
+Let \( \hat{r}_k \) be the residual at the \( k \)-th step: $$ -y_i = 5x_i^2 + 0.1\xi_i, \ i=1,\cdots,100 +\begin{equation*} +\hat{r}_k=\hat{b}-\hat{A}\hat{x}_k. +\end{equation*} $$ -with \( x_i \in [0,1] \) chosen randomly with a uniform distribution. Additionally \( \xi_i \) represents stochastic noise chosen according to a normal distribution \( \cal {N}(0,1) \). -The linear regression model is given by +Note that \( \hat{r}_k \) is the negative gradient of \( f \) at +\( \hat{x}=\hat{x}_k \), +so the gradient descent method would be to move in the direction \( \hat{r}_k \). +This gives the following expression $$ -h_\beta(x) = \hat{y} = \beta_0 + \beta_1 x, +\begin{equation*} +\hat{p}_{k+1}=\hat{r}_k-\frac{\hat{p}_k^T \hat{A}\hat{r}_k}{\hat{p}_k^T\hat{A}\hat{p}_k} \hat{p}_k. +\end{equation*} $$ +
@@ -262,7 +241,7 @@ $$
- + -
-Let \( \mathbf{y} = (y_1,\cdots,y_n)^T \), \( \mathbf{\hat{y}} = (\hat{y}_1,\cdots,\hat{y}_n)^T \) and \( \beta = (\beta_0, \beta_1)^T \) - -
-It is convenient to write \( \mathbf{\hat{y}} = X\beta \) where \( X \in \mathbb{R}^{100 \times 2} \) is the design matrix given by +
+We can also compute the residual iteratively as $$ -X \equiv \begin{bmatrix} -1 & x_1 \\ -\vdots & \vdots \\ -1 & x_{100} & \\ -\end{bmatrix}. +\begin{equation*} +\hat{r}_{k+1}=\hat{b}-\hat{A}\hat{x}_{k+1}, + \end{equation*} $$ -The loss function is given by +which equals $$ -C(\beta) = ||X\beta-\mathbf{y}||^2 = ||X\beta||^2 - 2 \mathbf{y}^T X\beta + ||\mathbf{y}||^2 = \sum_{i=1}^{100} (\beta_0 + \beta_1 x_i)^2 - 2 y_i (\beta_0 + \beta_1 x_i) + y_i^2 +\begin{equation*} +\hat{b}-\hat{A}(\hat{x}_k+\alpha_k\hat{p}_k), + \end{equation*} $$ -and we want to find \( \beta \) such that \( C(\beta) \) is minimized. +or +$$ +\begin{equation*} +(\hat{b}-\hat{A}\hat{x}_k)-\alpha_k\hat{A}\hat{p}_k, + \end{equation*} +$$ + +which gives + +$$ +\begin{equation*} +\hat{r}_{k+1}=\hat{r}_k-\hat{A}\hat{p}_{k}, + \end{equation*} +$$ +
@@ -254,7 +253,7 @@ and we want to find \( \beta \) such that \( C(\beta) \) is minimized.
-Computing \( \partial C(\beta) / \partial \beta_0 \) and \( \partial C(\beta) / \partial \beta_1 \) we can show that the gradient can be written as -$$ -\nabla_{\beta} C(\beta) = (\partial C(\beta) / \partial \beta_0, \partial C(\beta) / \partial \beta_1)^T = 2\begin{bmatrix} \sum_{i=1}^{100} \left(\beta_0+\beta_1x_i-y_i\right) \\ -\sum_{i=1}^{100}\left( x_i (\beta_0+\beta_1x_i)-y_ix_i\right) \\ -\end{bmatrix} = 2X^T(X\beta - \mathbf{y}), -$$ - -where \( X \) is the design matrix defined above. +
@@ -244,7 +219,7 @@ where \( X \) is the design matrix defined above.
+
+ + +
#include <cmath>
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+#include "vectormatrixclass.h"
+using namespace std;
+// Main function begins here
+int main(int argc, char * argv[]){
+ int dim = 2;
+ Vector x(dim),xsd(dim), b(dim),x0(dim);
+ Matrix A(dim,dim);
+
+ // Set our initial guess
+ x0(0) = x0(1) = 0;
+ // Set the matrix
+ A(0,0) = 3; A(1,0) = 2; A(0,1) = 2; A(1,1) = 6;
+ b(0) = 2; b(1) = -8;
+ cout << "The Matrix A that we are using: " << endl;
+ A.Print();
+ cout << endl;
+ xsd = SteepestDescent(A,b,x0);
+ cout << "The approximate solution using Steepest Descent is: " << endl;
+ xsd.Print();
+ cout << endl;
+}
++
@@ -243,7 +255,7 @@ This result implies that \( C(\beta) \) is a convex function since the matrix \(
-We can now write a program that minimizes \( C(\beta) \) using the gradient descent method with a constant learning rate \( \gamma \) according to -$$ -\beta_{k+1} = \beta_k - \gamma \nabla_\beta C(\beta_k), \ k=0,1,\cdots -$$ - -
-We can use the expression we computed for the gradient and let use a -\( \beta_0 \) be chosen randomly and let \( \gamma = 0.001 \). Stop iterating -when \( ||\nabla_\beta C(\beta_k) || \leq \epsilon = 10^{-8} \). - -
-And finally we can compare our solution for \( \beta \) with the analytic result given by -\( \beta= (X^TX)^{-1} X^T \mathbf{y} \). +
- -
import numpy as np
-
-"""
-The following setup is just a suggestion, feel free to write it the way you like.
-"""
-
-#Setup problem described in the exercise
-N = 100 #Nr of datapoints
-M = 2 #Nr of features
-x = np.random.rand(N) #Uniformly generated x-values in [0,1]
-y = 5*x**2 + 0.1*np.random.randn(N)
-X = np.c_[np.ones(N),x] #Construct design matrix
-
-#Compute beta according to normal equations to compare with GD solution
-Xt_X_inv = np.linalg.inv(np.dot(X.T,X))
-Xt_y = np.dot(X.transpose(),y)
-beta_NE = np.dot(Xt_X_inv,Xt_y)
-print(beta_NE)
+
+Vector SteepestDescent(Matrix A, Vector b, Vector x0){
+ int IterMax, i;
+ int dim = x0.Dimension();
+ const double tolerance = 1.0e-14;
+ Vector x(dim),f(dim),z(dim);
+ double c,alpha,d;
+ IterMax = 30;
+ x = x0;
+ f = A*x-b;
+ i = 0;
+ while (i <= IterMax){
+ z = A*f;
+ c = dot(f,f);
+ alpha = c/dot(f,z);
+ x = x - alpha*f;
+ f = A*x-b;
+ if(sqrt(dot(f,f)) < tolerance) break;
+ i++;
+ }
+ return x;
+}
+
+
@@ -270,7 +251,7 @@ beta_NE = np.35
-Another simple example is here
-
+We will use linear regression as a case study for the gradient descent
+methods. Linear regression is a great test case for the gradient
+descent methods discussed in the lectures since it has several
+desirable properties such as:
-
-
@@ -278,7 +247,7 @@ plt.show()
+Let \( \mathbf{y} = (y_1,\cdots,y_n)^T \), \( \mathbf{\hat{y}} = (\hat{y}_1,\cdots,\hat{y}_n)^T \) and \( \beta = (\beta_0, \beta_1)^T \)
-
-
+It is convenient to write \( \mathbf{\hat{y}} = X\beta \) where \( X \in \mathbb{R}^{100 \times 2} \) is the design matrix given by
+$$
+X \equiv \begin{bmatrix}
+1 & x_1 \\
+\vdots & \vdots \\
+1 & x_{100} & \\
+\end{bmatrix}.
+$$
-x = 2*np.random.rand(100,1)
-y = 4+3*x+np.random.randn(100,1)
+The loss function is given by
+$$
+C(\beta) = ||X\beta-\mathbf{y}||^2 = ||X\beta||^2 - 2 \mathbf{y}^T X\beta + ||\mathbf{y}||^2 = \sum_{i=1}^{100} (\beta_0 + \beta_1 x_i)^2 - 2 y_i (\beta_0 + \beta_1 x_i) + y_i^2
+$$
+
+and we want to find \( \beta \) such that \( C(\beta) \) is minimized.
-xb = np.c_[np.ones((100,1)), x]
-beta_linreg = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
-print(beta_linreg)
-sgdreg = SGDRegressor(n_iter = 50, penalty=None, eta0=0.1)
-sgdreg.fit(x,y.ravel())
-print(sgdreg.intercept_, sgdreg.coef_)
-
@@ -253,7 +239,7 @@ sgdreg.fit(x,y.
-We have also discussed Ridge regression where the loss function contains a regularized given by the \( L_2 \) norm of \( \beta \),
+Computing \( \partial C(\beta) / \partial \beta_0 \) and \( \partial C(\beta) / \partial \beta_1 \) we can show that the gradient can be written as
$$
-C_{\text{ridge}}(\beta) = ||X\beta -\mathbf{y}||^2 + \lambda ||\beta||^2, \ \lambda \geq 0.
-$$
-
-
-In order to minimize \( C_{\text{ridge}}(\beta) \) using GD we only have adjust the gradient as follows
-$$
-\nabla_\beta C_{\text{ridge}}(\beta) = 2\begin{bmatrix} \sum_{i=1}^{100} \left(\beta_0+\beta_1x_i-y_i\right) \\
+\nabla_{\beta} C(\beta) = (\partial C(\beta) / \partial \beta_0, \partial C(\beta) / \partial \beta_1)^T = 2\begin{bmatrix} \sum_{i=1}^{100} \left(\beta_0+\beta_1x_i-y_i\right) \\
\sum_{i=1}^{100}\left( x_i (\beta_0+\beta_1x_i)-y_ix_i\right) \\
-\end{bmatrix} + 2\lambda\begin{bmatrix} \beta_0 \\ \beta_1\end{bmatrix} = 2 (X^T(X\beta - \mathbf{y})+\lambda \beta).
+\end{bmatrix} = 2X^T(X\beta - \mathbf{y}),
$$
-
-We can now extend our program to minimize \( C_{\text{ridge}}(\beta) \) using gradient descent and compare with the analytical solution given by
-$$
-\beta_{\text{ridge}} = \left(X^T X + \lambda I_{2 \times 2} \right)^{-1} X^T \mathbf{y},
-$$
+where \( X \) is the design matrix defined above.
-for \( \lambda = {0,1,10,50,100} \) (\( \lambda = 0 \) corresponds to ordinary least squares).
-We can then compute \( ||\beta_{\text{ridge}}|| \) for each \( \lambda \).
-
-
-
-
-
@@ -286,7 +229,7 @@ beta_ridge = np
-Stochastic gradient descent (SGD) and variants thereof address some of
-the shortcomings of the Gradient descent method discussed above.
-
-
-The underlying idea of SGD comes from the observation that the cost
-function, which we want to minimize, can almost always be written as a
-sum over \( n \) data points \( \{\mathbf{x}_i\}_{i=1}^n \),
+
@@ -247,7 +228,7 @@ $$
-This in turn means that the gradient can be
-computed as a sum over \( i \)-gradients
+We can now write a program that minimizes \( C(\beta) \) using the gradient descent method with a constant learning rate \( \gamma \) according to
$$
-\nabla_\beta C(\mathbf{\beta}) = \sum_i^n \nabla_\beta c_i(\mathbf{x}_i,
-\mathbf{\beta}).
+\beta_{k+1} = \beta_k - \gamma \nabla_\beta C(\beta_k), \ k=0,1,\cdots
$$
-Stochasticity/randomness is introduced by only taking the
-gradient on a subset of the data called minibatches. If there are \( n \)
-data points and the size of each minibatch is \( M \), there will be \( n/M \)
-minibatches. We denote these minibatches by \( B_k \) where
-\( k=1,\cdots,n/M \).
+We can use the expression we computed for the gradient and let use a
+\( \beta_0 \) be chosen randomly and let \( \gamma = 0.001 \). Stop iterating
+when \( ||\nabla_\beta C(\beta_k) || \leq \epsilon = 10^{-8} \).
+
+And finally we can compare our solution for \( \beta \) with the analytic result given by
+\( \beta= (X^TX)^{-1} X^T \mathbf{y} \).
+
+
+
+
@@ -249,7 +255,7 @@ minibatches. We denote these minibatches by \( B_k \) where
-The idea is now to approximate the gradient by replacing the sum over
-all data points with a sum over the data points in one the minibatches
-picked at random in each gradient descent step
-$$
-\nabla_{\beta}
-C(\mathbf{\beta}) = \sum_{i=1}^n \nabla_\beta c_i(\mathbf{x}_i,
-\mathbf{\beta}) \rightarrow \sum_{i \in B_k}^n \nabla_\beta
-c_i(\mathbf{x}_i, \mathbf{\beta}).
-$$
+Another simple example is here
+
+
+
@@ -252,8 +262,6 @@ $$
-Thus a gradient descent step now looks like
-$$
-\beta_{j+1} = \beta_j - \gamma_j \sum_{i \in B_k}^n \nabla_\beta c_i(\mathbf{x}_i,
-\mathbf{\beta})
-$$
-
-where \( k \) is picked at random with equal
-probability from \( [1,n/M] \). An iteration over the number of
-minibathces (n/M) is commonly referred to as an epoch. Thus it is
-typical to choose a number of epochs and for each epoch iterate over
-the number of minibatches, as exemplified in the code below.
+
+
@@ -246,9 +236,6 @@ the number of minibatches, as exemplified in the code below.
+We have also discussed Ridge regression where the loss function contains a regularized given by the \( L_2 \) norm of \( \beta \),
+$$
+C_{\text{ridge}}(\beta) = ||X\beta -\mathbf{y}||^2 + \lambda ||\beta||^2, \ \lambda \geq 0.
+$$
+
+
+In order to minimize \( C_{\text{ridge}}(\beta) \) using GD we only have adjust the gradient as follows
+$$
+\nabla_\beta C_{\text{ridge}}(\beta) = 2\begin{bmatrix} \sum_{i=1}^{100} \left(\beta_0+\beta_1x_i-y_i\right) \\
+\sum_{i=1}^{100}\left( x_i (\beta_0+\beta_1x_i)-y_ix_i\right) \\
+\end{bmatrix} + 2\lambda\begin{bmatrix} \beta_0 \\ \beta_1\end{bmatrix} = 2 (X^T(X\beta - \mathbf{y})+\lambda \beta).
+$$
+
+
+We can now extend our program to minimize \( C_{\text{ridge}}(\beta) \) using gradient descent and compare with the analytical solution given by
+$$
+\beta_{\text{ridge}} = \left(X^T X + \lambda I_{2 \times 2} \right)^{-1} X^T \mathbf{y},
+$$
+
+for \( \lambda = {0,1,10,50,100} \) (\( \lambda = 0 \) corresponds to ordinary least squares).
+We can then compute \( ||\beta_{\text{ridge}}|| \) for each \( \lambda \).
-
-Taking the gradient only on a subset of the data has two important
-benefits. First, it introduces randomness which decreases the chance
-that our opmization scheme gets stuck in a local minima. Second, if
-the size of the minibatches are small relative to the number of
-datapoints (\( M < n \)), the computation of the gradient is much
-cheaper since we sum over the datapoints in the \( k-th \) minibatch and not
-all \( n \) datapoints.
-
@@ -258,10 +268,6 @@ all \( n \) datapoints.
-A natural question is when do we stop the search for a new minimum?
-One possibility is to compute the full gradient after a given number
-of epochs and check if the norm of the gradient is smaller than some
-threshold and stop if true. However, the condition that the gradient
-is zero is valid also for local minima, so this would only tell us
-that we are close to a local/global minimum. However, we could also
-evaluate the cost function at this point, store the result and
-continue the search. If the test kicks in at a later stage we can
-compare the values of the cost function and keep the \( \beta \) that
-gave the lowest value.
+Stochastic gradient descent (SGD) and variants thereof address some of
+the shortcomings of the Gradient descent method discussed above.
+
+
+The underlying idea of SGD comes from the observation that the cost
+function, which we want to minimize, can almost always be written as a
+sum over \( n \) data points \( \{\mathbf{x}_i\}_{i=1}^n \),
+$$
+C(\mathbf{\beta}) = \sum_{i=1}^n c_i(\mathbf{x}_i,
+\mathbf{\beta}).
+$$
@@ -242,11 +228,6 @@ gave the lowest value.
-Another approach is to let the step length \( \gamma_j \) depend on the
-number of epochs in such a way that it becomes very small after a
-reasonable time such that we do not move at all.
+This in turn means that the gradient can be
+computed as a sum over \( i \)-gradients
+$$
+\nabla_\beta C(\mathbf{\beta}) = \sum_i^n \nabla_\beta c_i(\mathbf{x}_i,
+\mathbf{\beta}).
+$$
-As an example, let \( e = 0,1,2,3,\cdots \) denote the current epoch and let \( t_0, t_1 > 0 \) be two fixed numbers. Furthermore, let \( t = e \cdot m + i \) where \( m \) is the number of minibatches and \( i=0,\cdots,m-1 \). Then the function $$\gamma_j(t; t_0, t_1) = \frac{t_0}{t+t_1} $$ goes to zero as the number of epochs gets large. I.e. we start with a step length \( \gamma_j (0; t_0, t_1) = t_0/t_1 \) which decays in time \( t \).
+Stochasticity/randomness is introduced by only taking the
+gradient on a subset of the data called minibatches. If there are \( n \)
+data points and the size of each minibatch is \( M \), there will be \( n/M \)
+minibatches. We denote these minibatches by \( B_k \) where
+\( k=1,\cdots,n/M \).
-
-In this way we can fix the number of epochs, compute \( \beta \) and
-evaluate the cost function at the end. Repeating the computation will
-give a different result since the scheme is random by design. Then we
-pick the final \( \beta \) that gives the lowest value of the cost
-function.
-
-
-
-
-
@@ -272,12 +229,6 @@ j = 0
The residual is zero when we reach the minimum of the quadratic equation
@@ -680,14 +680,150 @@ symmetric. If we search for a minimum of the quantum mechanical
variance, then the matrix \( \hat{A} \), which is called the Hessian, is
given by the second-derivative of the function we want to minimize.
This quantity is always positive definite.
-
-
-More details will be added here soon.
+We denote the initial guess for \( \hat{x} \) as \( \hat{x}_0 \).
+We can assume without loss of generality that
+
+One can show that the solution \( \hat{x} \) is also the unique minimizer of the quadratic form
+
+Let \( \hat{r}_k \) be the residual at the \( k \)-th step:
+
+We can also compute the residual iteratively as
+
@@ -713,11 +849,7 @@ More details will be added here soon.
cout << "The Matrix A that we are using: " << endl;
A.Print();
cout << endl;
- x = ConjugateGradient(A,b,x0);
xsd = SteepestDescent(A,b,x0);
- cout << "The approximate solution using Conjugate Gradient is: " << endl;
- x.Print();
- cout << endl;
cout << "The approximate solution using Steepest Descent is: " << endl;
xsd.Print();
cout << endl;
@@ -729,7 +861,7 @@ More details will be added here soon.
@@ -763,7 +895,7 @@ More details will be added here soon.
We will use linear regression as a case study for the gradient descent
@@ -803,7 +935,7 @@ $$
Let \( \mathbf{y} = (y_1,\cdots,y_n)^T \), \( \mathbf{\hat{y}} = (\hat{y}_1,\cdots,\hat{y}_n)^T \) and \( \beta = (\beta_0, \beta_1)^T \)
@@ -832,7 +964,7 @@ and we want to find \( \beta \) such that \( C(\beta) \) is minimized.
Computing \( \partial C(\beta) / \partial \beta_0 \) and \( \partial C(\beta) / \partial \beta_1 \) we can show that the gradient can be written as
@@ -849,7 +981,7 @@ where \( X \) is the design matrix defined above.
We can now write a program that minimizes \( C(\beta) \) using the gradient descent method with a constant learning rate \( \gamma \) according to
@@ -909,7 +1041,7 @@ beta_NE = np.dot(Xt_X_inv,Xt_y)
Another simple example is here
@@ -959,7 +1091,7 @@ plt.show()
@@ -984,7 +1116,7 @@ sgdreg.fit(x,y.ravel())
We have also discussed Ridge regression where the loss function contains a regularized given by the \( L_2 \) norm of \( \beta \),
@@ -1048,7 +1180,7 @@ beta_ridge = np.dot(Z,np.dot(X.T,y))
Stochastic gradient descent (SGD) and variants thereof address some of
@@ -1068,7 +1200,7 @@ $$
This in turn means that the gradient can be
@@ -1090,7 +1222,7 @@ minibatches. We denote these minibatches by \( B_k \) where
Thus a gradient descent step now looks like
@@ -1137,7 +1269,7 @@ the number of minibatches, as exemplified in the code below.
@@ -1169,7 +1301,7 @@ all \( n \) datapoints.
A natural question is when do we stop the search for a new minimum?
@@ -1186,7 +1318,7 @@ gave the lowest value.
Another approach is to let the step length \( \gamma_j \) depend on the
@@ -1236,355 +1368,6 @@ j = 0
-The success of the CG method for finding solutions of non-linear problems is based
-on the theory of conjugate gradients for linear systems of equations. It belongs
-to the class of iterative methods for solving problems from linear algebra of the type
-
-When we have found the exact solution, \( \hat{r}=0 \).
-
-The residual is zero when we reach the minimum of the quadratic equation
-
-We seek the minimum of the energy or the variance as function of various variational parameters.
-In our case we have thus a function \( f \) whose minimum we are seeking.
-In Newton's method we set \( \nabla f = 0 \) and we can thus compute the next iteration point
-
-In the CG method we define so-called conjugate directions and two vectors
-\( \hat{s} \) and \( \hat{t} \)
-are said to be
-conjugate if
-
-An example is given by the eigenvectors of the matrix
-
-Assume now that we have a symmetric positive-definite matrix \( \hat{A} \) of size
-\( n\times n \). At each iteration \( i+1 \) we obtain the conjugate direction of a vector
-
-The coefficients are given by
-
-If we choose the conjugate vectors \( \hat{p}_k \) carefully,
-then we may not need all of them to obtain a good approximation to the solution
-\( \hat{x} \).
-We want to regard the conjugate gradient method as an iterative method.
-This will us to solve systems where \( n \) is so large that the direct
-method would take too much time.
-
-
-We denote the initial guess for \( \hat{x} \) as \( \hat{x}_0 \).
-We can assume without loss of generality that
-
-One can show that the solution \( \hat{x} \) is also the unique minimizer of the quadratic form
-
-Let \( \hat{r}_k \) be the residual at the \( k \)-th step:
-
-We can also compute the residual iteratively as
-
The residual is zero when we reach the minimum of the quadratic equation
@@ -657,12 +648,131 @@ given by the second-derivative of the function we want to minimize.
This quantity is always positive definite.
-More details will be added here soon.
+
+We denote the initial guess for \( \hat{x} \) as \( \hat{x}_0 \).
+We can assume without loss of generality that
+$$
+\begin{equation*}
+\hat{x}_0=0,
+\end{equation*}
+$$
+
+or consider the system
+$$
+\begin{equation*}
+\hat{A}\hat{z} = \hat{b}-\hat{A}\hat{x}_0,
+\end{equation*}
+$$
+
+instead.
+One can show that the solution \( \hat{x} \) is also the unique minimizer of the quadratic form
+$$
+\begin{equation*}
+ f(\hat{x}) = \frac{1}{2}\hat{x}^T\hat{A}\hat{x} - \hat{x}^T \hat{x} , \quad \hat{x}\in\mathbf{R}^n.
+\end{equation*}
+$$
+
+This suggests taking the first basis vector \( \hat{p}_1 \)
+to be the gradient of \( f \) at \( \hat{x}=\hat{x}_0 \),
+which equals
+$$
+\begin{equation*}
+\hat{A}\hat{x}_0-\hat{b},
+\end{equation*}
+$$
+
+and
+\( \hat{x}_0=0 \) it is equal \( -\hat{b} \).
+
+
+
+
+Let \( \hat{r}_k \) be the residual at the \( k \)-th step:
+$$
+\begin{equation*}
+\hat{r}_k=\hat{b}-\hat{A}\hat{x}_k.
+\end{equation*}
+$$
+
+Note that \( \hat{r}_k \) is the negative gradient of \( f \) at
+\( \hat{x}=\hat{x}_k \),
+so the gradient descent method would be to move in the direction \( \hat{r}_k \).
+This gives the following expression
+$$
+\begin{equation*}
+\hat{p}_{k+1}=\hat{r}_k-\frac{\hat{p}_k^T \hat{A}\hat{r}_k}{\hat{p}_k^T\hat{A}\hat{p}_k} \hat{p}_k.
+\end{equation*}
+$$
+
+
+We can also compute the residual iteratively as
+$$
+\begin{equation*}
+\hat{r}_{k+1}=\hat{b}-\hat{A}\hat{x}_{k+1},
+ \end{equation*}
+$$
+
+which equals
+$$
+\begin{equation*}
+\hat{b}-\hat{A}(\hat{x}_k+\alpha_k\hat{p}_k),
+ \end{equation*}
+$$
+
+or
+$$
+\begin{equation*}
+(\hat{b}-\hat{A}\hat{x}_k)-\alpha_k\hat{A}\hat{p}_k,
+ \end{equation*}
+$$
+
+which gives
+
+$$
+\begin{equation*}
+\hat{r}_{k+1}=\hat{r}_k-\hat{A}\hat{p}_{k},
+ \end{equation*}
+$$
+
+
+
@@ -689,11 +799,7 @@ More details will be added here soon.
cout << "The Matrix A that we are using: " << endl;
A.Print();
cout << endl;
- x = ConjugateGradient(A,b,x0);
xsd = SteepestDescent(A,b,x0);
- cout << "The approximate solution using Conjugate Gradient is: " << endl;
- x.Print();
- cout << endl;
cout << "The approximate solution using Steepest Descent is: " << endl;
xsd.Print();
cout << endl;
@@ -706,7 +812,7 @@ More details will be added here soon.
@@ -742,7 +848,7 @@ More details will be added here soon.
-
We will use linear regression as a case study for the gradient descent
@@ -775,7 +881,7 @@ $$
-
Let \( \mathbf{y} = (y_1,\cdots,y_n)^T \), \( \mathbf{\hat{y}} = (\hat{y}_1,\cdots,\hat{y}_n)^T \) and \( \beta = (\beta_0, \beta_1)^T \)
@@ -800,7 +906,7 @@ and we want to find \( \beta \) such that \( C(\beta) \) is minimized.
Computing \( \partial C(\beta) / \partial \beta_0 \) and \( \partial C(\beta) / \partial \beta_1 \) we can show that the gradient can be written as
@@ -815,7 +921,7 @@ where \( X \) is the design matrix defined above.
We can now write a program that minimizes \( C(\beta) \) using the gradient descent method with a constant learning rate \( \gamma \) according to
@@ -870,7 +976,7 @@ beta_NE = np.dot(Xt_X_inv,Xt_y)
Another simple example is here
@@ -919,7 +1025,7 @@ plt.show()
@@ -943,7 +1049,7 @@ sgdreg.fit(x,y.ravel())
-
We have also discussed Ridge regression where the loss function contains a regularized given by the \( L_2 \) norm of \( \beta \),
@@ -1000,7 +1106,7 @@ beta_ridge = np.dot(Z,np.dot(X.T,y))
Stochastic gradient descent (SGD) and variants thereof address some of
@@ -1018,7 +1124,7 @@ $$
This in turn means that the gradient can be
@@ -1038,7 +1144,7 @@ minibatches. We denote these minibatches by \( B_k \) where
Thus a gradient descent step now looks like
@@ -1081,7 +1187,7 @@ the number of minibatches, as exemplified in the code below.
@@ -1113,7 +1219,7 @@ all \( n \) datapoints.
A natural question is when do we stop the search for a new minimum?
@@ -1130,7 +1236,7 @@ gave the lowest value.
Another approach is to let the step length \( \gamma_j \) depend on the
@@ -1175,324 +1281,6 @@ j = 0
print("gamma_j after %d epochs: %g" % (n_epochs,gamma_j))
-
-The success of the CG method for finding solutions of non-linear problems is based
-on the theory of conjugate gradients for linear systems of equations. It belongs
-to the class of iterative methods for solving problems from linear algebra of the type
-$$
-\begin{equation*}
- \hat{A}\hat{x} = \hat{b}.
-\end{equation*}
-$$
-
-In the iterative process we end up with a problem like
-
-$$
-\begin{equation*}
- \hat{r}= \hat{b}-\hat{A}\hat{x},
-\end{equation*}
-$$
-
-where \( \hat{r} \) is the so-called residual or error in the iterative process.
-
-
-When we have found the exact solution, \( \hat{r}=0 \).
-
-
-
-
-The residual is zero when we reach the minimum of the quadratic equation
-$$
-\begin{equation*}
- P(\hat{x})=\frac{1}{2}\hat{x}^T\hat{A}\hat{x} - \hat{x}^T\hat{b},
-\end{equation*}
-$$
-
-with the constraint that the matrix \( \hat{A} \) is positive definite and symmetric.
-If we search for a minimum of the quantum mechanical variance, then the matrix
-\( \hat{A} \), which is called the Hessian, is given by the second-derivative of the function we want to minimize. This quantity is always positive definite. In our case this corresponds normally to the second derivative of the energy.
-
-
-We seek the minimum of the energy or the variance as function of various variational parameters.
-In our case we have thus a function \( f \) whose minimum we are seeking.
-In Newton's method we set \( \nabla f = 0 \) and we can thus compute the next iteration point
-$$
-\begin{equation*}
-\hat{x}-\hat{x}_i=\hat{A}^{-1}\nabla f(\hat{x}_i).
-\end{equation*}
-$$
-
-Subtracting this equation from that of \( \hat{x}_{i+1} \) we have
-$$
-\begin{equation*}
-\hat{x}_{i+1}-\hat{x}_i=\hat{A}^{-1}(\nabla f(\hat{x}_{i+1})-\nabla f(\hat{x}_i)).
-\end{equation*}
-$$
-
-
-In the CG method we define so-called conjugate directions and two vectors
-\( \hat{s} \) and \( \hat{t} \)
-are said to be
-conjugate if
-$$
-\begin{equation*}
-\hat{s}^T\hat{A}\hat{t}= 0.
-\end{equation*}
-$$
-
-The philosophy of the CG method is to perform searches in various conjugate directions
-of our vectors \( \hat{x}_i \) obeying the above criterion, namely
-$$
-\begin{equation*}
-\hat{x}_i^T\hat{A}\hat{x}_j= 0.
-\end{equation*}
-$$
-
-Two vectors are conjugate if they are orthogonal with respect to
-this inner product. Being conjugate is a symmetric relation: if \( \hat{s} \) is conjugate to \( \hat{t} \), then \( \hat{t} \) is conjugate to \( \hat{s} \).
-
-
-An example is given by the eigenvectors of the matrix
-$$
-\begin{equation*}
-\hat{v}_i^T\hat{A}\hat{v}_j= \lambda\hat{v}_i^T\hat{v}_j,
-\end{equation*}
-$$
-
-which is zero unless \( i=j \).
-
-
-Assume now that we have a symmetric positive-definite matrix \( \hat{A} \) of size
-\( n\times n \). At each iteration \( i+1 \) we obtain the conjugate direction of a vector
-$$
-\begin{equation*}
-\hat{x}_{i+1}=\hat{x}_{i}+\alpha_i\hat{p}_{i}.
-\end{equation*}
-$$
-
-We assume that \( \hat{p}_{i} \) is a sequence of \( n \) mutually conjugate directions.
-Then the \( \hat{p}_{i} \) form a basis of \( R^n \) and we can expand the solution
-$ \hat{A}\hat{x} = \hat{b}$ in this basis, namely
-
-$$
-\begin{equation*}
- \hat{x} = \sum^{n}_{i=1} \alpha_i \hat{p}_i.
-\end{equation*}
-$$
-
-
-The coefficients are given by
-$$
-\begin{equation*}
- \mathbf{A}\mathbf{x} = \sum^{n}_{i=1} \alpha_i \mathbf{A} \mathbf{p}_i = \mathbf{b}.
-\end{equation*}
-$$
-
-Multiplying with \( \hat{p}_k^T \) from the left gives
-
-$$
-\begin{equation*}
- \hat{p}_k^T \hat{A}\hat{x} = \sum^{n}_{i=1} \alpha_i\hat{p}_k^T \hat{A}\hat{p}_i= \hat{p}_k^T \hat{b},
-\end{equation*}
-$$
-
-and we can define the coefficients \( \alpha_k \) as
-
-$$
-\begin{equation*}
- \alpha_k = \frac{\hat{p}_k^T \hat{b}}{\hat{p}_k^T \hat{A} \hat{p}_k}
-\end{equation*}
-$$
-
-
-
-
-If we choose the conjugate vectors \( \hat{p}_k \) carefully,
-then we may not need all of them to obtain a good approximation to the solution
-\( \hat{x} \).
-We want to regard the conjugate gradient method as an iterative method.
-This will us to solve systems where \( n \) is so large that the direct
-method would take too much time.
-
-
-We denote the initial guess for \( \hat{x} \) as \( \hat{x}_0 \).
-We can assume without loss of generality that
-$$
-\begin{equation*}
-\hat{x}_0=0,
-\end{equation*}
-$$
-
-or consider the system
-$$
-\begin{equation*}
-\hat{A}\hat{z} = \hat{b}-\hat{A}\hat{x}_0,
-\end{equation*}
-$$
-
-instead.
-
-
-One can show that the solution \( \hat{x} \) is also the unique minimizer of the quadratic form
-$$
-\begin{equation*}
- f(\hat{x}) = \frac{1}{2}\hat{x}^T\hat{A}\hat{x} - \hat{x}^T \hat{x} , \quad \hat{x}\in\mathbf{R}^n.
-\end{equation*}
-$$
-
-This suggests taking the first basis vector \( \hat{p}_1 \)
-to be the gradient of \( f \) at \( \hat{x}=\hat{x}_0 \),
-which equals
-$$
-\begin{equation*}
-\hat{A}\hat{x}_0-\hat{b},
-\end{equation*}
-$$
-
-and
-\( \hat{x}_0=0 \) it is equal \( -\hat{b} \).
-The other vectors in the basis will be conjugate to the gradient,
-hence the name conjugate gradient method.
-
-
-Let \( \hat{r}_k \) be the residual at the \( k \)-th step:
-$$
-\begin{equation*}
-\hat{r}_k=\hat{b}-\hat{A}\hat{x}_k.
-\end{equation*}
-$$
-
-Note that \( \hat{r}_k \) is the negative gradient of \( f \) at
-\( \hat{x}=\hat{x}_k \),
-so the gradient descent method would be to move in the direction \( \hat{r}_k \).
-Here, we insist that the directions \( \hat{p}_k \) are conjugate to each other,
-so we take the direction closest to the gradient \( \hat{r}_k \)
-under the conjugacy constraint.
-This gives the following expression
-$$
-\begin{equation*}
-\hat{p}_{k+1}=\hat{r}_k-\frac{\hat{p}_k^T \hat{A}\hat{r}_k}{\hat{p}_k^T\hat{A}\hat{p}_k} \hat{p}_k.
-\end{equation*}
-$$
-
-
-We can also compute the residual iteratively as
-$$
-\begin{equation*}
-\hat{r}_{k+1}=\hat{b}-\hat{A}\hat{x}_{k+1},
- \end{equation*}
-$$
-
-which equals
-$$
-\begin{equation*}
-\hat{b}-\hat{A}(\hat{x}_k+\alpha_k\hat{p}_k),
- \end{equation*}
-$$
-
-or
-$$
-\begin{equation*}
-(\hat{b}-\hat{A}\hat{x}_k)-\alpha_k\hat{A}\hat{p}_k,
- \end{equation*}
-$$
-
-which gives
-
-$$
-\begin{equation*}
-\hat{r}_{k+1}=\hat{r}_k-\hat{A}\hat{p}_{k},
- \end{equation*}
-$$
-
diff --git a/doc/pub/Splines/html/Splines.html b/doc/pub/Splines/html/Splines.html
index 6feb11062..30682d1e4 100644
--- a/doc/pub/Splines/html/Splines.html
+++ b/doc/pub/Splines/html/Splines.html
@@ -90,48 +90,39 @@ div { text-align: justify; text-justify: inter-word; }
('More on convex functions', 2, None, '___sec15'),
('Some simple problems', 2, None, '___sec16'),
('Standard steepest descent', 2, None, '___sec17'),
- ('Conjugate gradient method', 2, None, '___sec18'),
+ ('Gradient method', 2, None, '___sec18'),
+ ('Steepest descent method', 2, None, '___sec19'),
+ ('Steepest descent method', 2, None, '___sec20'),
+ ('Gradient descent method', 2, None, '___sec21'),
+ ('Final expressions', 2, None, '___sec22'),
+ ('The Steepest descent algorithm', 2, None, '___sec23'),
('Simple codes for steepest descent and conjugate gradient '
'using a $2\\times 2$ matrix, in c++, Python code to come',
2,
None,
- '___sec19'),
+ '___sec24'),
('The routine for the steepest descent method',
2,
None,
- '___sec20'),
- ('Revisiting our first homework', 2, None, '___sec21'),
- ('Gradient descent example', 2, None, '___sec22'),
- ('The derivative of the cost/loss function', 2, None, '___sec23'),
- ('The Hessian matrix', 2, None, '___sec24'),
- ('Simple program', 2, None, '___sec25'),
- ('Gradient Descent Example', 2, None, '___sec26'),
+ '___sec25'),
+ ('Revisiting our first homework', 2, None, '___sec26'),
+ ('Gradient descent example', 2, None, '___sec27'),
+ ('The derivative of the cost/loss function', 2, None, '___sec28'),
+ ('The Hessian matrix', 2, None, '___sec29'),
+ ('Simple program', 2, None, '___sec30'),
+ ('Gradient Descent Example', 2, None, '___sec31'),
('And a corresponding example using _scikit-learn_',
2,
None,
- '___sec27'),
- ('Gradient descent and Ridge', 2, None, '___sec28'),
- ('Stochastic Gradient Descent', 2, None, '___sec29'),
- ('Computation of gradients', 2, None, '___sec30'),
- ('SGD example', 2, None, '___sec31'),
- ('The gradient step', 2, None, '___sec32'),
- ('Simple example code', 2, None, '___sec33'),
- ('When do we stop?', 2, None, '___sec34'),
- ('Slightly different approach', 2, None, '___sec35'),
- ('Conjugate gradient (CG) method', 2, None, '___sec36'),
- ('Conjugate gradient method', 2, None, '___sec37'),
- ("Conjugate gradient method, Newton's method first",
- 2,
- None,
- '___sec38'),
- ('Conjugate gradient method', 2, None, '___sec39'),
- ('Conjugate gradient method', 2, None, '___sec40'),
- ('Conjugate gradient method', 2, None, '___sec41'),
- ('Conjugate gradient method', 2, None, '___sec42'),
- ('Conjugate gradient method and iterations', 2, None, '___sec43'),
- ('Conjugate gradient method', 2, None, '___sec44'),
- ('Conjugate gradient method', 2, None, '___sec45'),
- ('Conjugate gradient method', 2, None, '___sec46')]}
+ '___sec32'),
+ ('Gradient descent and Ridge', 2, None, '___sec33'),
+ ('Stochastic Gradient Descent', 2, None, '___sec34'),
+ ('Computation of gradients', 2, None, '___sec35'),
+ ('SGD example', 2, None, '___sec36'),
+ ('The gradient step', 2, None, '___sec37'),
+ ('Simple example code', 2, None, '___sec38'),
+ ('When do we stop?', 2, None, '___sec39'),
+ ('Slightly different approach', 2, None, '___sec40')]}
end of tocinfo -->
The residual is zero when we reach the minimum of the quadratic equation
@@ -662,12 +653,131 @@ given by the second-derivative of the function we want to minimize.
This quantity is always positive definite.
-More details will be added here soon.
+
+We denote the initial guess for \( \hat{x} \) as \( \hat{x}_0 \).
+We can assume without loss of generality that
+$$
+\begin{equation*}
+\hat{x}_0=0,
+\end{equation*}
+$$
+
+or consider the system
+$$
+\begin{equation*}
+\hat{A}\hat{z} = \hat{b}-\hat{A}\hat{x}_0,
+\end{equation*}
+$$
+
+instead.
+One can show that the solution \( \hat{x} \) is also the unique minimizer of the quadratic form
+$$
+\begin{equation*}
+ f(\hat{x}) = \frac{1}{2}\hat{x}^T\hat{A}\hat{x} - \hat{x}^T \hat{x} , \quad \hat{x}\in\mathbf{R}^n.
+\end{equation*}
+$$
+
+This suggests taking the first basis vector \( \hat{p}_1 \)
+to be the gradient of \( f \) at \( \hat{x}=\hat{x}_0 \),
+which equals
+$$
+\begin{equation*}
+\hat{A}\hat{x}_0-\hat{b},
+\end{equation*}
+$$
+
+and
+\( \hat{x}_0=0 \) it is equal \( -\hat{b} \).
+
+
+
+
+Let \( \hat{r}_k \) be the residual at the \( k \)-th step:
+$$
+\begin{equation*}
+\hat{r}_k=\hat{b}-\hat{A}\hat{x}_k.
+\end{equation*}
+$$
+
+Note that \( \hat{r}_k \) is the negative gradient of \( f \) at
+\( \hat{x}=\hat{x}_k \),
+so the gradient descent method would be to move in the direction \( \hat{r}_k \).
+This gives the following expression
+$$
+\begin{equation*}
+\hat{p}_{k+1}=\hat{r}_k-\frac{\hat{p}_k^T \hat{A}\hat{r}_k}{\hat{p}_k^T\hat{A}\hat{p}_k} \hat{p}_k.
+\end{equation*}
+$$
+
+
+We can also compute the residual iteratively as
+$$
+\begin{equation*}
+\hat{r}_{k+1}=\hat{b}-\hat{A}\hat{x}_{k+1},
+ \end{equation*}
+$$
+
+which equals
+$$
+\begin{equation*}
+\hat{b}-\hat{A}(\hat{x}_k+\alpha_k\hat{p}_k),
+ \end{equation*}
+$$
+
+or
+$$
+\begin{equation*}
+(\hat{b}-\hat{A}\hat{x}_k)-\alpha_k\hat{A}\hat{p}_k,
+ \end{equation*}
+$$
+
+which gives
+
+$$
+\begin{equation*}
+\hat{r}_{k+1}=\hat{r}_k-\hat{A}\hat{p}_{k},
+ \end{equation*}
+$$
+
+
+
@@ -694,11 +804,7 @@ More details will be added here soon.
cout << "The Matrix A that we are using: " << endl;
A.Print();
cout << endl;
- x = ConjugateGradient(A,b,x0);
xsd = SteepestDescent(A,b,x0);
- cout << "The approximate solution using Conjugate Gradient is: " << endl;
- x.Print();
- cout << endl;
cout << "The approximate solution using Steepest Descent is: " << endl;
xsd.Print();
cout << endl;
@@ -711,7 +817,7 @@ More details will be added here soon.
@@ -747,7 +853,7 @@ More details will be added here soon.
-
We will use linear regression as a case study for the gradient descent
@@ -780,7 +886,7 @@ $$
-
Let \( \mathbf{y} = (y_1,\cdots,y_n)^T \), \( \mathbf{\hat{y}} = (\hat{y}_1,\cdots,\hat{y}_n)^T \) and \( \beta = (\beta_0, \beta_1)^T \)
@@ -805,7 +911,7 @@ and we want to find \( \beta \) such that \( C(\beta) \) is minimized.
Computing \( \partial C(\beta) / \partial \beta_0 \) and \( \partial C(\beta) / \partial \beta_1 \) we can show that the gradient can be written as
@@ -820,7 +926,7 @@ where \( X \) is the design matrix defined above.
We can now write a program that minimizes \( C(\beta) \) using the gradient descent method with a constant learning rate \( \gamma \) according to
@@ -875,7 +981,7 @@ beta_NE = np.
Another simple example is here
@@ -924,7 +1030,7 @@ plt.show()
@@ -948,7 +1054,7 @@ sgdreg.fit(x,y.
-
We have also discussed Ridge regression where the loss function contains a regularized given by the \( L_2 \) norm of \( \beta \),
@@ -1005,7 +1111,7 @@ beta_ridge = np
Stochastic gradient descent (SGD) and variants thereof address some of
@@ -1023,7 +1129,7 @@ $$
This in turn means that the gradient can be
@@ -1043,7 +1149,7 @@ minibatches. We denote these minibatches by \( B_k \) where
Thus a gradient descent step now looks like
@@ -1086,7 +1192,7 @@ the number of minibatches, as exemplified in the code below.
@@ -1118,7 +1224,7 @@ all \( n \) datapoints.
A natural question is when do we stop the search for a new minimum?
@@ -1135,7 +1241,7 @@ gave the lowest value.
Another approach is to let the step length \( \gamma_j \) depend on the
@@ -1180,324 +1286,6 @@ j = 0
print("gamma_j after %d epochs: %g" % (n_epochs,gamma_j))
-
-The success of the CG method for finding solutions of non-linear problems is based
-on the theory of conjugate gradients for linear systems of equations. It belongs
-to the class of iterative methods for solving problems from linear algebra of the type
-$$
-\begin{equation*}
- \hat{A}\hat{x} = \hat{b}.
-\end{equation*}
-$$
-
-In the iterative process we end up with a problem like
-
-$$
-\begin{equation*}
- \hat{r}= \hat{b}-\hat{A}\hat{x},
-\end{equation*}
-$$
-
-where \( \hat{r} \) is the so-called residual or error in the iterative process.
-
-
-When we have found the exact solution, \( \hat{r}=0 \).
-
-
-
-
-The residual is zero when we reach the minimum of the quadratic equation
-$$
-\begin{equation*}
- P(\hat{x})=\frac{1}{2}\hat{x}^T\hat{A}\hat{x} - \hat{x}^T\hat{b},
-\end{equation*}
-$$
-
-with the constraint that the matrix \( \hat{A} \) is positive definite and symmetric.
-If we search for a minimum of the quantum mechanical variance, then the matrix
-\( \hat{A} \), which is called the Hessian, is given by the second-derivative of the function we want to minimize. This quantity is always positive definite. In our case this corresponds normally to the second derivative of the energy.
-
-
-We seek the minimum of the energy or the variance as function of various variational parameters.
-In our case we have thus a function \( f \) whose minimum we are seeking.
-In Newton's method we set \( \nabla f = 0 \) and we can thus compute the next iteration point
-$$
-\begin{equation*}
-\hat{x}-\hat{x}_i=\hat{A}^{-1}\nabla f(\hat{x}_i).
-\end{equation*}
-$$
-
-Subtracting this equation from that of \( \hat{x}_{i+1} \) we have
-$$
-\begin{equation*}
-\hat{x}_{i+1}-\hat{x}_i=\hat{A}^{-1}(\nabla f(\hat{x}_{i+1})-\nabla f(\hat{x}_i)).
-\end{equation*}
-$$
-
-
-In the CG method we define so-called conjugate directions and two vectors
-\( \hat{s} \) and \( \hat{t} \)
-are said to be
-conjugate if
-$$
-\begin{equation*}
-\hat{s}^T\hat{A}\hat{t}= 0.
-\end{equation*}
-$$
-
-The philosophy of the CG method is to perform searches in various conjugate directions
-of our vectors \( \hat{x}_i \) obeying the above criterion, namely
-$$
-\begin{equation*}
-\hat{x}_i^T\hat{A}\hat{x}_j= 0.
-\end{equation*}
-$$
-
-Two vectors are conjugate if they are orthogonal with respect to
-this inner product. Being conjugate is a symmetric relation: if \( \hat{s} \) is conjugate to \( \hat{t} \), then \( \hat{t} \) is conjugate to \( \hat{s} \).
-
-
-An example is given by the eigenvectors of the matrix
-$$
-\begin{equation*}
-\hat{v}_i^T\hat{A}\hat{v}_j= \lambda\hat{v}_i^T\hat{v}_j,
-\end{equation*}
-$$
-
-which is zero unless \( i=j \).
-
-
-Assume now that we have a symmetric positive-definite matrix \( \hat{A} \) of size
-\( n\times n \). At each iteration \( i+1 \) we obtain the conjugate direction of a vector
-$$
-\begin{equation*}
-\hat{x}_{i+1}=\hat{x}_{i}+\alpha_i\hat{p}_{i}.
-\end{equation*}
-$$
-
-We assume that \( \hat{p}_{i} \) is a sequence of \( n \) mutually conjugate directions.
-Then the \( \hat{p}_{i} \) form a basis of \( R^n \) and we can expand the solution
-$ \hat{A}\hat{x} = \hat{b}$ in this basis, namely
-
-$$
-\begin{equation*}
- \hat{x} = \sum^{n}_{i=1} \alpha_i \hat{p}_i.
-\end{equation*}
-$$
-
-
-The coefficients are given by
-$$
-\begin{equation*}
- \mathbf{A}\mathbf{x} = \sum^{n}_{i=1} \alpha_i \mathbf{A} \mathbf{p}_i = \mathbf{b}.
-\end{equation*}
-$$
-
-Multiplying with \( \hat{p}_k^T \) from the left gives
-
-$$
-\begin{equation*}
- \hat{p}_k^T \hat{A}\hat{x} = \sum^{n}_{i=1} \alpha_i\hat{p}_k^T \hat{A}\hat{p}_i= \hat{p}_k^T \hat{b},
-\end{equation*}
-$$
-
-and we can define the coefficients \( \alpha_k \) as
-
-$$
-\begin{equation*}
- \alpha_k = \frac{\hat{p}_k^T \hat{b}}{\hat{p}_k^T \hat{A} \hat{p}_k}
-\end{equation*}
-$$
-
-
-
-
-If we choose the conjugate vectors \( \hat{p}_k \) carefully,
-then we may not need all of them to obtain a good approximation to the solution
-\( \hat{x} \).
-We want to regard the conjugate gradient method as an iterative method.
-This will us to solve systems where \( n \) is so large that the direct
-method would take too much time.
-
-
-We denote the initial guess for \( \hat{x} \) as \( \hat{x}_0 \).
-We can assume without loss of generality that
-$$
-\begin{equation*}
-\hat{x}_0=0,
-\end{equation*}
-$$
-
-or consider the system
-$$
-\begin{equation*}
-\hat{A}\hat{z} = \hat{b}-\hat{A}\hat{x}_0,
-\end{equation*}
-$$
-
-instead.
-
-
-One can show that the solution \( \hat{x} \) is also the unique minimizer of the quadratic form
-$$
-\begin{equation*}
- f(\hat{x}) = \frac{1}{2}\hat{x}^T\hat{A}\hat{x} - \hat{x}^T \hat{x} , \quad \hat{x}\in\mathbf{R}^n.
-\end{equation*}
-$$
-
-This suggests taking the first basis vector \( \hat{p}_1 \)
-to be the gradient of \( f \) at \( \hat{x}=\hat{x}_0 \),
-which equals
-$$
-\begin{equation*}
-\hat{A}\hat{x}_0-\hat{b},
-\end{equation*}
-$$
-
-and
-\( \hat{x}_0=0 \) it is equal \( -\hat{b} \).
-The other vectors in the basis will be conjugate to the gradient,
-hence the name conjugate gradient method.
-
-
-Let \( \hat{r}_k \) be the residual at the \( k \)-th step:
-$$
-\begin{equation*}
-\hat{r}_k=\hat{b}-\hat{A}\hat{x}_k.
-\end{equation*}
-$$
-
-Note that \( \hat{r}_k \) is the negative gradient of \( f \) at
-\( \hat{x}=\hat{x}_k \),
-so the gradient descent method would be to move in the direction \( \hat{r}_k \).
-Here, we insist that the directions \( \hat{p}_k \) are conjugate to each other,
-so we take the direction closest to the gradient \( \hat{r}_k \)
-under the conjugacy constraint.
-This gives the following expression
-$$
-\begin{equation*}
-\hat{p}_{k+1}=\hat{r}_k-\frac{\hat{p}_k^T \hat{A}\hat{r}_k}{\hat{p}_k^T\hat{A}\hat{p}_k} \hat{p}_k.
-\end{equation*}
-$$
-
-
-We can also compute the residual iteratively as
-$$
-\begin{equation*}
-\hat{r}_{k+1}=\hat{b}-\hat{A}\hat{x}_{k+1},
- \end{equation*}
-$$
-
-which equals
-$$
-\begin{equation*}
-\hat{b}-\hat{A}(\hat{x}_k+\alpha_k\hat{p}_k),
- \end{equation*}
-$$
-
-or
-$$
-\begin{equation*}
-(\hat{b}-\hat{A}\hat{x}_k)-\alpha_k\hat{A}\hat{p}_k,
- \end{equation*}
-$$
-
-which gives
-
-$$
-\begin{equation*}
-\hat{r}_{k+1}=\hat{r}_k-\hat{A}\hat{p}_{k},
- \end{equation*}
-$$
-
diff --git a/doc/pub/Splines/ipynb/Splines.ipynb b/doc/pub/Splines/ipynb/Splines.ipynb
index a46ac7740..b6216c22c 100644
--- a/doc/pub/Splines/ipynb/Splines.ipynb
+++ b/doc/pub/Splines/ipynb/Splines.ipynb
@@ -585,7 +585,7 @@
"\n",
"When we have found the exact solution, $\\hat{r}=0$.\n",
"\n",
- "## Conjugate gradient method\n",
+ "## Gradient method\n",
"\n",
"The residual is zero when we reach the minimum of the quadratic equation"
]
@@ -609,7 +609,189 @@
"given by the second-derivative of the function we want to minimize.\n",
"This quantity is always positive definite. \n",
"\n",
- "More details will be added here soon.\n",
+ "\n",
+ "## Steepest descent method\n",
+ "\n",
+ "We denote the initial guess for $\\hat{x}$ as $\\hat{x}_0$. \n",
+ "We can assume without loss of generality that"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\hat{x}_0=0,\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "or consider the system"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\hat{A}\\hat{z} = \\hat{b}-\\hat{A}\\hat{x}_0,\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "instead.\n",
+ "\n",
+ "\n",
+ "## Steepest descent method\n",
+ "One can show that the solution $\\hat{x}$ is also the unique minimizer of the quadratic form"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "f(\\hat{x}) = \\frac{1}{2}\\hat{x}^T\\hat{A}\\hat{x} - \\hat{x}^T \\hat{x} , \\quad \\hat{x}\\in\\mathbf{R}^n.\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This suggests taking the first basis vector $\\hat{p}_1$ \n",
+ "to be the gradient of $f$ at $\\hat{x}=\\hat{x}_0$, \n",
+ "which equals"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\hat{A}\\hat{x}_0-\\hat{b},\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "and \n",
+ "$\\hat{x}_0=0$ it is equal $-\\hat{b}$.\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "## Gradient descent method\n",
+ "Let $\\hat{r}_k$ be the residual at the $k$-th step:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\hat{r}_k=\\hat{b}-\\hat{A}\\hat{x}_k.\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Note that $\\hat{r}_k$ is the negative gradient of $f$ at \n",
+ "$\\hat{x}=\\hat{x}_k$, \n",
+ "so the gradient descent method would be to move in the direction $\\hat{r}_k$. \n",
+ "This gives the following expression"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\hat{p}_{k+1}=\\hat{r}_k-\\frac{\\hat{p}_k^T \\hat{A}\\hat{r}_k}{\\hat{p}_k^T\\hat{A}\\hat{p}_k} \\hat{p}_k.\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Final expressions\n",
+ "We can also compute the residual iteratively as"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\hat{r}_{k+1}=\\hat{b}-\\hat{A}\\hat{x}_{k+1},\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "which equals"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\hat{b}-\\hat{A}(\\hat{x}_k+\\alpha_k\\hat{p}_k),\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "or"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "(\\hat{b}-\\hat{A}\\hat{x}_k)-\\alpha_k\\hat{A}\\hat{p}_k,\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "which gives"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\hat{r}_{k+1}=\\hat{r}_k-\\hat{A}\\hat{p}_{k},\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## The Steepest descent algorithm\n",
+ "\n",
"\n",
"## Simple codes for steepest descent and conjugate gradient using a $2\\times 2$ matrix, in c++, Python code to come"
]
@@ -638,11 +820,7 @@
" cout << \"The Matrix A that we are using: \" << endl;\n",
" A.Print();\n",
" cout << endl;\n",
- " x = ConjugateGradient(A,b,x0);\n",
" xsd = SteepestDescent(A,b,x0);\n",
- " cout << \"The approximate solution using Conjugate Gradient is: \" << endl;\n",
- " x.Print();\n",
- " cout << endl;\n",
" cout << \"The approximate solution using Steepest Descent is: \" << endl;\n",
" xsd.Print();\n",
" cout << endl;\n",
@@ -1289,454 +1467,6 @@
"\n",
"print(\"gamma_j after %d epochs: %g\" % (n_epochs,gamma_j))"
]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Conjugate gradient (CG) method\n",
- "The success of the CG method for finding solutions of non-linear problems is based\n",
- "on the theory of conjugate gradients for linear systems of equations. It belongs\n",
- "to the class of iterative methods for solving problems from linear algebra of the type"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\hat{A}\\hat{x} = \\hat{b}.\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "In the iterative process we end up with a problem like"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\hat{r}= \\hat{b}-\\hat{A}\\hat{x},\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "where $\\hat{r}$ is the so-called residual or error in the iterative process.\n",
- "\n",
- "When we have found the exact solution, $\\hat{r}=0$.\n",
- "\n",
- "\n",
- "\n",
- "\n",
- "## Conjugate gradient method\n",
- "\n",
- "The residual is zero when we reach the minimum of the quadratic equation"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "P(\\hat{x})=\\frac{1}{2}\\hat{x}^T\\hat{A}\\hat{x} - \\hat{x}^T\\hat{b},\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "with the constraint that the matrix $\\hat{A}$ is positive definite and symmetric.\n",
- "If we search for a minimum of the quantum mechanical variance, then the matrix \n",
- "$\\hat{A}$, which is called the Hessian, is given by the second-derivative of the function we want to minimize. This quantity is always positive definite. In our case this corresponds normally to the second derivative of the energy.\n",
- "\n",
- "\n",
- "\n",
- "\n",
- "\n",
- "\n",
- "\n",
- "## Conjugate gradient method, Newton's method first\n",
- "We seek the minimum of the energy or the variance as function of various variational parameters. \n",
- "In our case we have thus a function $f$ whose minimum we are seeking.\n",
- "In Newton's method we set $\\nabla f = 0$ and we can thus compute the next iteration point"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\hat{x}-\\hat{x}_i=\\hat{A}^{-1}\\nabla f(\\hat{x}_i).\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Subtracting this equation from that of $\\hat{x}_{i+1}$ we have"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\hat{x}_{i+1}-\\hat{x}_i=\\hat{A}^{-1}(\\nabla f(\\hat{x}_{i+1})-\\nabla f(\\hat{x}_i)).\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Conjugate gradient method\n",
- "In the CG method we define so-called conjugate directions and two vectors \n",
- "$\\hat{s}$ and $\\hat{t}$\n",
- "are said to be\n",
- "conjugate if"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\hat{s}^T\\hat{A}\\hat{t}= 0.\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "The philosophy of the CG method is to perform searches in various conjugate directions\n",
- "of our vectors $\\hat{x}_i$ obeying the above criterion, namely"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\hat{x}_i^T\\hat{A}\\hat{x}_j= 0.\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Two vectors are conjugate if they are orthogonal with respect to \n",
- "this inner product. Being conjugate is a symmetric relation: if $\\hat{s}$ is conjugate to $\\hat{t}$, then $\\hat{t}$ is conjugate to $\\hat{s}$.\n",
- "\n",
- "\n",
- "\n",
- "## Conjugate gradient method\n",
- "An example is given by the eigenvectors of the matrix"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\hat{v}_i^T\\hat{A}\\hat{v}_j= \\lambda\\hat{v}_i^T\\hat{v}_j,\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "which is zero unless $i=j$.\n",
- "\n",
- "\n",
- "\n",
- "\n",
- "## Conjugate gradient method\n",
- "Assume now that we have a symmetric positive-definite matrix $\\hat{A}$ of size\n",
- "$n\\times n$. At each iteration $i+1$ we obtain the conjugate direction of a vector"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\hat{x}_{i+1}=\\hat{x}_{i}+\\alpha_i\\hat{p}_{i}.\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "We assume that $\\hat{p}_{i}$ is a sequence of $n$ mutually conjugate directions. \n",
- "Then the $\\hat{p}_{i}$ form a basis of $R^n$ and we can expand the solution \n",
- "$ \\hat{A}\\hat{x} = \\hat{b}$ in this basis, namely"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\hat{x} = \\sum^{n}_{i=1} \\alpha_i \\hat{p}_i.\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Conjugate gradient method\n",
- "The coefficients are given by"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\mathbf{A}\\mathbf{x} = \\sum^{n}_{i=1} \\alpha_i \\mathbf{A} \\mathbf{p}_i = \\mathbf{b}.\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Multiplying with $\\hat{p}_k^T$ from the left gives"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\hat{p}_k^T \\hat{A}\\hat{x} = \\sum^{n}_{i=1} \\alpha_i\\hat{p}_k^T \\hat{A}\\hat{p}_i= \\hat{p}_k^T \\hat{b},\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "and we can define the coefficients $\\alpha_k$ as"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\alpha_k = \\frac{\\hat{p}_k^T \\hat{b}}{\\hat{p}_k^T \\hat{A} \\hat{p}_k}\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Conjugate gradient method and iterations\n",
- "\n",
- "If we choose the conjugate vectors $\\hat{p}_k$ carefully, \n",
- "then we may not need all of them to obtain a good approximation to the solution \n",
- "$\\hat{x}$. \n",
- "We want to regard the conjugate gradient method as an iterative method. \n",
- "This will us to solve systems where $n$ is so large that the direct \n",
- "method would take too much time.\n",
- "\n",
- "We denote the initial guess for $\\hat{x}$ as $\\hat{x}_0$. \n",
- "We can assume without loss of generality that"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\hat{x}_0=0,\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "or consider the system"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\hat{A}\\hat{z} = \\hat{b}-\\hat{A}\\hat{x}_0,\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "instead.\n",
- "\n",
- "\n",
- "\n",
- "\n",
- "## Conjugate gradient method\n",
- "One can show that the solution $\\hat{x}$ is also the unique minimizer of the quadratic form"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "f(\\hat{x}) = \\frac{1}{2}\\hat{x}^T\\hat{A}\\hat{x} - \\hat{x}^T \\hat{x} , \\quad \\hat{x}\\in\\mathbf{R}^n.\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "This suggests taking the first basis vector $\\hat{p}_1$ \n",
- "to be the gradient of $f$ at $\\hat{x}=\\hat{x}_0$, \n",
- "which equals"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\hat{A}\\hat{x}_0-\\hat{b},\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "and \n",
- "$\\hat{x}_0=0$ it is equal $-\\hat{b}$.\n",
- "The other vectors in the basis will be conjugate to the gradient, \n",
- "hence the name conjugate gradient method.\n",
- "\n",
- "\n",
- "\n",
- "\n",
- "## Conjugate gradient method\n",
- "Let $\\hat{r}_k$ be the residual at the $k$-th step:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\hat{r}_k=\\hat{b}-\\hat{A}\\hat{x}_k.\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Note that $\\hat{r}_k$ is the negative gradient of $f$ at \n",
- "$\\hat{x}=\\hat{x}_k$, \n",
- "so the gradient descent method would be to move in the direction $\\hat{r}_k$. \n",
- "Here, we insist that the directions $\\hat{p}_k$ are conjugate to each other, \n",
- "so we take the direction closest to the gradient $\\hat{r}_k$ \n",
- "under the conjugacy constraint. \n",
- "This gives the following expression"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\hat{p}_{k+1}=\\hat{r}_k-\\frac{\\hat{p}_k^T \\hat{A}\\hat{r}_k}{\\hat{p}_k^T\\hat{A}\\hat{p}_k} \\hat{p}_k.\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Conjugate gradient method\n",
- "We can also compute the residual iteratively as"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\hat{r}_{k+1}=\\hat{b}-\\hat{A}\\hat{x}_{k+1},\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "which equals"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\hat{b}-\\hat{A}(\\hat{x}_k+\\alpha_k\\hat{p}_k),\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "or"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "(\\hat{b}-\\hat{A}\\hat{x}_k)-\\alpha_k\\hat{A}\\hat{p}_k,\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "which gives"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\hat{r}_{k+1}=\\hat{r}_k-\\hat{A}\\hat{p}_{k},\n",
- "$$"
- ]
}
],
"metadata": {},
diff --git a/doc/pub/Splines/ipynb/ipynb-Splines-src.tar.gz b/doc/pub/Splines/ipynb/ipynb-Splines-src.tar.gz
index 9e4fe066b..b478c41fa 100644
Binary files a/doc/pub/Splines/ipynb/ipynb-Splines-src.tar.gz and b/doc/pub/Splines/ipynb/ipynb-Splines-src.tar.gz differ
diff --git a/doc/pub/Splines/pdf/Splines-minted.pdf b/doc/pub/Splines/pdf/Splines-minted.pdf
index b40638855..c83605159 100644
Binary files a/doc/pub/Splines/pdf/Splines-minted.pdf and b/doc/pub/Splines/pdf/Splines-minted.pdf differ
diff --git a/doc/src/Splines/Splines.do.txt b/doc/src/Splines/Splines.do.txt
index 51367b610..acf8e10b2 100644
--- a/doc/src/Splines/Splines.do.txt
+++ b/doc/src/Splines/Splines.do.txt
@@ -405,7 +405,7 @@ where $\hat{r}$ is the so-called residual or error in the iterative process.
When we have found the exact solution, $\hat{r}=0$.
!split
-===== Conjugate gradient method =====
+===== Gradient method =====
The residual is zero when we reach the minimum of the quadratic equation
!bt
@@ -420,7 +420,104 @@ variance, then the matrix $\hat{A}$, which is called the Hessian, is
given by the second-derivative of the function we want to minimize.
This quantity is always positive definite.
-More details will be added here soon.
+
+!split
+===== Steepest descent method =====
+
+We denote the initial guess for $\hat{x}$ as $\hat{x}_0$.
+We can assume without loss of generality that
+!bt
+\begin{equation*}
+\hat{x}_0=0,
+\end{equation*}
+!et
+or consider the system
+!bt
+\begin{equation*}
+\hat{A}\hat{z} = \hat{b}-\hat{A}\hat{x}_0,
+\end{equation*}
+!et
+instead.
+
+
+!split
+===== Steepest descent method =====
+!bblock
+One can show that the solution $\hat{x}$ is also the unique minimizer of the quadratic form
+!bt
+\begin{equation*}
+ f(\hat{x}) = \frac{1}{2}\hat{x}^T\hat{A}\hat{x} - \hat{x}^T \hat{x} , \quad \hat{x}\in\mathbf{R}^n.
+\end{equation*}
+!et
+This suggests taking the first basis vector $\hat{p}_1$
+to be the gradient of $f$ at $\hat{x}=\hat{x}_0$,
+which equals
+!bt
+\begin{equation*}
+\hat{A}\hat{x}_0-\hat{b},
+\end{equation*}
+!et
+and
+$\hat{x}_0=0$ it is equal $-\hat{b}$.
+
+!eblock
+
+
+!split
+===== Gradient descent method =====
+!bblock
+Let $\hat{r}_k$ be the residual at the $k$-th step:
+!bt
+\begin{equation*}
+\hat{r}_k=\hat{b}-\hat{A}\hat{x}_k.
+\end{equation*}
+!et
+Note that $\hat{r}_k$ is the negative gradient of $f$ at
+$\hat{x}=\hat{x}_k$,
+so the gradient descent method would be to move in the direction $\hat{r}_k$.
+This gives the following expression
+!bt
+\begin{equation*}
+\hat{p}_{k+1}=\hat{r}_k-\frac{\hat{p}_k^T \hat{A}\hat{r}_k}{\hat{p}_k^T\hat{A}\hat{p}_k} \hat{p}_k.
+\end{equation*}
+!et
+!eblock
+
+!split
+===== Final expressions =====
+!bblock
+We can also compute the residual iteratively as
+!bt
+\begin{equation*}
+\hat{r}_{k+1}=\hat{b}-\hat{A}\hat{x}_{k+1},
+ \end{equation*}
+!et
+which equals
+!bt
+\begin{equation*}
+\hat{b}-\hat{A}(\hat{x}_k+\alpha_k\hat{p}_k),
+ \end{equation*}
+!et
+or
+!bt
+\begin{equation*}
+(\hat{b}-\hat{A}\hat{x}_k)-\alpha_k\hat{A}\hat{p}_k,
+ \end{equation*}
+!et
+which gives
+
+!bt
+\begin{equation*}
+\hat{r}_{k+1}=\hat{r}_k-\hat{A}\hat{p}_{k},
+ \end{equation*}
+!et
+!eblock
+
+
+
+!split
+===== The Steepest descent algorithm =====
+
!split
===== Simple codes for steepest descent and conjugate gradient using a $2\times 2$ matrix, in c++, Python code to come =====
@@ -446,11 +543,7 @@ int main(int argc, char * argv[]){
cout << "The Matrix A that we are using: " << endl;
A.Print();
cout << endl;
- x = ConjugateGradient(A,b,x0);
xsd = SteepestDescent(A,b,x0);
- cout << "The approximate solution using Conjugate Gradient is: " << endl;
- x.Print();
- cout << endl;
cout << "The approximate solution using Steepest Descent is: " << endl;
xsd.Print();
cout << endl;
@@ -895,255 +988,6 @@ print("gamma_j after %d epochs: %g" % (n_epochs,gamma_j))
-!split
-===== Conjugate gradient (CG) method =====
-!bblock
-The success of the CG method for finding solutions of non-linear problems is based
-on the theory of conjugate gradients for linear systems of equations. It belongs
-to the class of iterative methods for solving problems from linear algebra of the type
-!bt
-\begin{equation*}
- \hat{A}\hat{x} = \hat{b}.
-\end{equation*}
-!et
-In the iterative process we end up with a problem like
-
-!bt
-\begin{equation*}
- \hat{r}= \hat{b}-\hat{A}\hat{x},
-\end{equation*}
-!et
-where $\hat{r}$ is the so-called residual or error in the iterative process.
-
-When we have found the exact solution, $\hat{r}=0$.
-!eblock
-
-
-!split
-===== Conjugate gradient method =====
-!bblock
-
-The residual is zero when we reach the minimum of the quadratic equation
-!bt
-\begin{equation*}
- P(\hat{x})=\frac{1}{2}\hat{x}^T\hat{A}\hat{x} - \hat{x}^T\hat{b},
-\end{equation*}
-!et
-with the constraint that the matrix $\hat{A}$ is positive definite and symmetric.
-If we search for a minimum of the quantum mechanical variance, then the matrix
-$\hat{A}$, which is called the Hessian, is given by the second-derivative of the function we want to minimize. This quantity is always positive definite. In our case this corresponds normally to the second derivative of the energy.
-!eblock
-
-
-
-
-
-!split
-===== Conjugate gradient method, Newton's method first =====
-!bblock
-We seek the minimum of the energy or the variance as function of various variational parameters.
-In our case we have thus a function $f$ whose minimum we are seeking.
-In Newton's method we set $\nabla f = 0$ and we can thus compute the next iteration point
-!bt
-\begin{equation*}
-\hat{x}-\hat{x}_i=\hat{A}^{-1}\nabla f(\hat{x}_i).
-\end{equation*}
-!et
-Subtracting this equation from that of $\hat{x}_{i+1}$ we have
-!bt
-\begin{equation*}
-\hat{x}_{i+1}-\hat{x}_i=\hat{A}^{-1}(\nabla f(\hat{x}_{i+1})-\nabla f(\hat{x}_i)).
-\end{equation*}
-!et
-!eblock
-
-
-!split
-===== Conjugate gradient method =====
-!bblock
-In the CG method we define so-called conjugate directions and two vectors
-$\hat{s}$ and $\hat{t}$
-are said to be
-conjugate if
-!bt
-\begin{equation*}
-\hat{s}^T\hat{A}\hat{t}= 0.
-\end{equation*}
-!et
-The philosophy of the CG method is to perform searches in various conjugate directions
-of our vectors $\hat{x}_i$ obeying the above criterion, namely
-!bt
-\begin{equation*}
-\hat{x}_i^T\hat{A}\hat{x}_j= 0.
-\end{equation*}
-!et
-Two vectors are conjugate if they are orthogonal with respect to
-this inner product. Being conjugate is a symmetric relation: if $\hat{s}$ is conjugate to $\hat{t}$, then $\hat{t}$ is conjugate to $\hat{s}$.
-!eblock
-
-!split
-===== Conjugate gradient method =====
-!bblock
-An example is given by the eigenvectors of the matrix
-!bt
-\begin{equation*}
-\hat{v}_i^T\hat{A}\hat{v}_j= \lambda\hat{v}_i^T\hat{v}_j,
-\end{equation*}
-!et
-which is zero unless $i=j$.
-!eblock
-
-
-!split
-===== Conjugate gradient method =====
-!bblock
-Assume now that we have a symmetric positive-definite matrix $\hat{A}$ of size
-$n\times n$. At each iteration $i+1$ we obtain the conjugate direction of a vector
-!bt
-\begin{equation*}
-\hat{x}_{i+1}=\hat{x}_{i}+\alpha_i\hat{p}_{i}.
-\end{equation*}
-!et
-We assume that $\hat{p}_{i}$ is a sequence of $n$ mutually conjugate directions.
-Then the $\hat{p}_{i}$ form a basis of $R^n$ and we can expand the solution
-$ \hat{A}\hat{x} = \hat{b}$ in this basis, namely
-
-!bt
-\begin{equation*}
- \hat{x} = \sum^{n}_{i=1} \alpha_i \hat{p}_i.
-\end{equation*}
-!et
-!eblock
-
-!split
-===== Conjugate gradient method =====
-!bblock
-The coefficients are given by
-!bt
-\begin{equation*}
- \mathbf{A}\mathbf{x} = \sum^{n}_{i=1} \alpha_i \mathbf{A} \mathbf{p}_i = \mathbf{b}.
-\end{equation*}
-!et
-Multiplying with $\hat{p}_k^T$ from the left gives
-
-!bt
-\begin{equation*}
- \hat{p}_k^T \hat{A}\hat{x} = \sum^{n}_{i=1} \alpha_i\hat{p}_k^T \hat{A}\hat{p}_i= \hat{p}_k^T \hat{b},
-\end{equation*}
-!et
-and we can define the coefficients $\alpha_k$ as
-
-!bt
-\begin{equation*}
- \alpha_k = \frac{\hat{p}_k^T \hat{b}}{\hat{p}_k^T \hat{A} \hat{p}_k}
-\end{equation*}
-!et
-!eblock
-
-!split
-===== Conjugate gradient method and iterations =====
-!bblock
-
-If we choose the conjugate vectors $\hat{p}_k$ carefully,
-then we may not need all of them to obtain a good approximation to the solution
-$\hat{x}$.
-We want to regard the conjugate gradient method as an iterative method.
-This will us to solve systems where $n$ is so large that the direct
-method would take too much time.
-
-We denote the initial guess for $\hat{x}$ as $\hat{x}_0$.
-We can assume without loss of generality that
-!bt
-\begin{equation*}
-\hat{x}_0=0,
-\end{equation*}
-!et
-or consider the system
-!bt
-\begin{equation*}
-\hat{A}\hat{z} = \hat{b}-\hat{A}\hat{x}_0,
-\end{equation*}
-!et
-instead.
-!eblock
-
-
-!split
-===== Conjugate gradient method =====
-!bblock
-One can show that the solution $\hat{x}$ is also the unique minimizer of the quadratic form
-!bt
-\begin{equation*}
- f(\hat{x}) = \frac{1}{2}\hat{x}^T\hat{A}\hat{x} - \hat{x}^T \hat{x} , \quad \hat{x}\in\mathbf{R}^n.
-\end{equation*}
-!et
-This suggests taking the first basis vector $\hat{p}_1$
-to be the gradient of $f$ at $\hat{x}=\hat{x}_0$,
-which equals
-!bt
-\begin{equation*}
-\hat{A}\hat{x}_0-\hat{b},
-\end{equation*}
-!et
-and
-$\hat{x}_0=0$ it is equal $-\hat{b}$.
-The other vectors in the basis will be conjugate to the gradient,
-hence the name conjugate gradient method.
-!eblock
-
-
-!split
-===== Conjugate gradient method =====
-!bblock
-Let $\hat{r}_k$ be the residual at the $k$-th step:
-!bt
-\begin{equation*}
-\hat{r}_k=\hat{b}-\hat{A}\hat{x}_k.
-\end{equation*}
-!et
-Note that $\hat{r}_k$ is the negative gradient of $f$ at
-$\hat{x}=\hat{x}_k$,
-so the gradient descent method would be to move in the direction $\hat{r}_k$.
-Here, we insist that the directions $\hat{p}_k$ are conjugate to each other,
-so we take the direction closest to the gradient $\hat{r}_k$
-under the conjugacy constraint.
-This gives the following expression
-!bt
-\begin{equation*}
-\hat{p}_{k+1}=\hat{r}_k-\frac{\hat{p}_k^T \hat{A}\hat{r}_k}{\hat{p}_k^T\hat{A}\hat{p}_k} \hat{p}_k.
-\end{equation*}
-!et
-!eblock
-
-!split
-===== Conjugate gradient method =====
-!bblock
-We can also compute the residual iteratively as
-!bt
-\begin{equation*}
-\hat{r}_{k+1}=\hat{b}-\hat{A}\hat{x}_{k+1},
- \end{equation*}
-!et
-which equals
-!bt
-\begin{equation*}
-\hat{b}-\hat{A}(\hat{x}_k+\alpha_k\hat{p}_k),
- \end{equation*}
-!et
-or
-!bt
-\begin{equation*}
-(\hat{b}-\hat{A}\hat{x}_k)-\alpha_k\hat{A}\hat{p}_k,
- \end{equation*}
-!et
-which gives
-
-!bt
-\begin{equation*}
-\hat{r}_{k+1}=\hat{r}_k-\hat{A}\hat{p}_{k},
- \end{equation*}
-!et
-!eblock
Gradient Descent Example
+Revisiting our first homework
# Importing various packages
-from random import random, seed
-import numpy as np
-import matplotlib.pyplot as plt
-from mpl_toolkits.mplot3d import Axes3D
-from matplotlib import cm
-from matplotlib.ticker import LinearLocator, FormatStrFormatter
-import sys
+
+
-x = 2*np.random.rand(100,1)
-y = 4+3*x+np.random.randn(100,1)
+We revisit the example from homework set 1 where we had
+$$
+y_i = 5x_i^2 + 0.1\xi_i, \ i=1,\cdots,100
+$$
-xb = np.c_[np.ones((100,1)), x]
-beta_linreg = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
-print(beta_linreg)
-beta = np.random.randn(2,1)
+with \( x_i \in [0,1] \) chosen randomly with a uniform distribution. Additionally \( \xi_i \) represents stochastic noise chosen according to a normal distribution \( \cal {N}(0,1) \).
+The linear regression model is given by
+$$
+h_\beta(x) = \hat{y} = \beta_0 + \beta_1 x,
+$$
-eta = 0.1
-Niterations = 1000
-m = 100
+such that
+$$
+\hat{y}_i = \beta_0 + \beta_1 x_i.
+$$
-for iter in range(Niterations):
- gradients = 2.0/m*xb.T.dot(xb.dot(beta)-y)
- beta -= eta*gradients
-
-print(beta)
-xnew = np.array([[0],[2]])
-xbnew = np.c_[np.ones((2,1)), xnew]
-ypredict = xbnew.dot(beta)
-ypredict2 = xbnew.dot(beta_linreg)
-plt.plot(xnew, ypredict, "r-")
-plt.plot(xnew, ypredict2, "b-")
-plt.plot(x, y ,'ro')
-plt.axis([0,2.0,0, 15.0])
-plt.xlabel(r'$x$')
-plt.ylabel(r'$y$')
-plt.title(r'Gradient descent example')
-plt.show()
-And a corresponding example using scikit-learn
+Gradient descent example
# Importing various packages
-from random import random, seed
-import numpy as np
-import matplotlib.pyplot as plt
-from sklearn.linear_model import SGDRegressor
+
Gradient descent and Ridge
+The derivative of the cost/loss function
import numpy as np
-
-"""
-The following setup is just a suggestion, feel free to write it the way you like.
-"""
-
-#Setup problem described in the exercise
-N = 100 #Nr of datapoints
-M = 2 #Nr of features
-x = np.random.rand(N)
-y = 5*x**2 + 0.1*np.random.randn(N)
-
-
-#Compute analytic beta for Ridge regression
-X = np.c_[np.ones(N),x]
-XT_X = np.dot(X.T,X)
-
-l = 0.1 #Ridge parameter lambda
-Id = np.eye(XT_X.shape[0])
-
-Z = np.linalg.inv(XT_X+l*Id)
-beta_ridge = np.dot(Z,np.dot(X.T,y))
-
-print(beta_ridge)
-print(np.linalg.norm(beta_ridge)) #||beta||
-
Stochastic Gradient Descent
-
-The Hessian matrix
+The Hessian matrix of \( C(\beta) \) is given by
$$
-C(\mathbf{\beta}) = \sum_{i=1}^n c_i(\mathbf{x}_i,
-\mathbf{\beta}).
+\hat{H} \equiv \begin{bmatrix}
+\frac{\partial^2 C(\beta)}{\partial \beta_0^2} & \frac{\partial^2 C(\beta)}{\partial \beta_0 \partial \beta_1} \\
+\frac{\partial^2 C(\beta)}{\partial \beta_0 \partial \beta_1} & \frac{\partial^2 C(\beta)}{\partial \beta_1^2} & \\
+\end{bmatrix} = 2X^T X.
$$
+This result implies that \( C(\beta) \) is a convex function since the matrix \( X^T X \) always is positive semi-definite.
+
Computation of gradients
+Simple program
import numpy as np
+
+"""
+The following setup is just a suggestion, feel free to write it the way you like.
+"""
+
+#Setup problem described in the exercise
+N = 100 #Nr of datapoints
+M = 2 #Nr of features
+x = np.random.rand(N) #Uniformly generated x-values in [0,1]
+y = 5*x**2 + 0.1*np.random.randn(N)
+X = np.c_[np.ones(N),x] #Construct design matrix
+
+#Compute beta according to normal equations to compare with GD solution
+Xt_X_inv = np.linalg.inv(np.dot(X.T,X))
+Xt_y = np.dot(X.transpose(),y)
+beta_NE = np.dot(Xt_X_inv,Xt_y)
+print(beta_NE)
+
SGD example
-As an example, suppose we have \( 10 \) data points \( (\mathbf{x}_1,\cdots, \mathbf{x}_{10}) \)
-and we choose to have \( M=5 \) minibathces,
-then each minibatch contains two data points. In particular we have
-\( B_1 = (\mathbf{x}_1,\mathbf{x}_2), \cdots, B_5 =
-(\mathbf{x}_9,\mathbf{x}_{10}) \). Note that if you choose \( M=1 \) you
-have only a single batch with all data points and on the other extreme,
-you may choose \( M=n \) resulting in a minibatch for each datapoint, i.e
-\( B_k = \mathbf{x}_k \).
+Gradient Descent Example
# Importing various packages
+from random import random, seed
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+from matplotlib import cm
+from matplotlib.ticker import LinearLocator, FormatStrFormatter
+import sys
+
+x = 2*np.random.rand(100,1)
+y = 4+3*x+np.random.randn(100,1)
+
+xb = np.c_[np.ones((100,1)), x]
+beta_linreg = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
+print(beta_linreg)
+beta = np.random.randn(2,1)
+
+eta = 0.1
+Niterations = 1000
+m = 100
+
+for iter in range(Niterations):
+ gradients = 2.0/m*xb.T.dot(xb.dot(beta)-y)
+ beta -= eta*gradients
+
+print(beta)
+xnew = np.array([[0],[2]])
+xbnew = np.c_[np.ones((2,1)), xnew]
+ypredict = xbnew.dot(beta)
+ypredict2 = xbnew.dot(beta_linreg)
+plt.plot(xnew, ypredict, "r-")
+plt.plot(xnew, ypredict2, "b-")
+plt.plot(x, y ,'ro')
+plt.axis([0,2.0,0, 15.0])
+plt.xlabel(r'$x$')
+plt.ylabel(r'$y$')
+plt.title(r'Gradient descent example')
+plt.show()
+
The gradient step
+And a corresponding example using scikit-learn
# Importing various packages
+from random import random, seed
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.linear_model import SGDRegressor
+x = 2*np.random.rand(100,1)
+y = 4+3*x+np.random.randn(100,1)
+
+xb = np.c_[np.ones((100,1)), x]
+beta_linreg = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
+print(beta_linreg)
+sgdreg = SGDRegressor(n_iter = 50, penalty=None, eta0=0.1)
+sgdreg.fit(x,y.ravel())
+print(sgdreg.intercept_, sgdreg.coef_)
+
Simple example code
+Gradient descent and Ridge
+
+import numpy as np
+
import numpy as np
-n = 100 #100 datapoints
-M = 5 #size of each minibatch
-m = int(n/M) #number of minibatches
-n_epochs = 10 #number of epochs
+"""
+The following setup is just a suggestion, feel free to write it the way you like.
+"""
-j = 0
-for epoch in range(1,n_epochs+1):
- for i in range(m):
- k = np.random.randint(m) #Pick the k-th minibatch at random
- #Compute the gradient using the data in minibatch Bk
- #Compute new suggestion for
- j += 1
+#Setup problem described in the exercise
+N = 100 #Nr of datapoints
+M = 2 #Nr of features
+x = np.random.rand(N)
+y = 5*x**2 + 0.1*np.random.randn(N)
+
+
+#Compute analytic beta for Ridge regression
+X = np.c_[np.ones(N),x]
+XT_X = np.dot(X.T,X)
+
+l = 0.1 #Ridge parameter lambda
+Id = np.eye(XT_X.shape[0])
+
+Z = np.linalg.inv(XT_X+l*Id)
+beta_ridge = np.dot(Z,np.dot(X.T,y))
+
+print(beta_ridge)
+print(np.linalg.norm(beta_ridge)) #||beta||
When do we stop?
+Stochastic Gradient Descent
Slightly different approach
+Computation of gradients
import numpy as np
-
-def step_length(t,t0,t1):
- return t0/(t+t1)
-
-n = 100 #100 datapoints
-M = 5 #size of each minibatch
-m = int(n/M) #number of minibatches
-n_epochs = 500 #number of epochs
-t0 = 1.0
-t1 = 10
-
-gamma_j = t0/t1
-j = 0
-for epoch in range(1,n_epochs+1):
- for i in range(m):
- k = np.random.randint(m) #Pick the k-th minibatch at random
- #Compute the gradient using the data in minibatch Bk
- #Compute new suggestion for beta
- t = epoch*m+i
- gamma_j = step_length(t,t0,t1)
- j += 1
-
-print("gamma_j after %d epochs: %g" % (n_epochs,gamma_j))
-
Conjugate gradient method
+Gradient method
Simple codes for steepest descent and conjugate gradient using a \( 2\times 2 \) matrix, in c++, Python code to come
+Steepest descent method
+
+
+$$
+\begin{equation*}
+\hat{x}_0=0,
+\end{equation*}
+$$
+
+
+or consider the system
+
+$$
+\begin{equation*}
+\hat{A}\hat{z} = \hat{b}-\hat{A}\hat{x}_0,
+\end{equation*}
+$$
+
+
+instead.
+Steepest descent method
+
+$$
+\begin{equation*}
+ f(\hat{x}) = \frac{1}{2}\hat{x}^T\hat{A}\hat{x} - \hat{x}^T \hat{x} , \quad \hat{x}\in\mathbf{R}^n.
+\end{equation*}
+$$
+
+
+This suggests taking the first basis vector \( \hat{p}_1 \)
+to be the gradient of \( f \) at \( \hat{x}=\hat{x}_0 \),
+which equals
+
+$$
+\begin{equation*}
+\hat{A}\hat{x}_0-\hat{b},
+\end{equation*}
+$$
+
+
+and
+\( \hat{x}_0=0 \) it is equal \( -\hat{b} \).
+
+
+Gradient descent method
+
+$$
+\begin{equation*}
+\hat{r}_k=\hat{b}-\hat{A}\hat{x}_k.
+\end{equation*}
+$$
+
+
+Note that \( \hat{r}_k \) is the negative gradient of \( f \) at
+\( \hat{x}=\hat{x}_k \),
+so the gradient descent method would be to move in the direction \( \hat{r}_k \).
+This gives the following expression
+
+$$
+\begin{equation*}
+\hat{p}_{k+1}=\hat{r}_k-\frac{\hat{p}_k^T \hat{A}\hat{r}_k}{\hat{p}_k^T\hat{A}\hat{p}_k} \hat{p}_k.
+\end{equation*}
+$$
+
+Final expressions
+
+$$
+\begin{equation*}
+\hat{r}_{k+1}=\hat{b}-\hat{A}\hat{x}_{k+1},
+ \end{equation*}
+$$
+
+
+which equals
+
+$$
+\begin{equation*}
+\hat{b}-\hat{A}(\hat{x}_k+\alpha_k\hat{p}_k),
+ \end{equation*}
+$$
+
+
+or
+
+$$
+\begin{equation*}
+(\hat{b}-\hat{A}\hat{x}_k)-\alpha_k\hat{A}\hat{p}_k,
+ \end{equation*}
+$$
+
+
+which gives
+
+
+$$
+\begin{equation*}
+\hat{r}_{k+1}=\hat{r}_k-\hat{A}\hat{p}_{k},
+ \end{equation*}
+$$
+
+The Steepest descent algorithm
+Simple codes for steepest descent and conjugate gradient using a \( 2\times 2 \) matrix, in c++, Python code to come
The routine for the steepest descent method
+The routine for the steepest descent method
Revisiting our first homework
+Revisiting our first homework
Gradient descent example
+Gradient descent example
The derivative of the cost/loss function
+The derivative of the cost/loss function
The Hessian matrix
+The Hessian matrix
The Hessian matrix of \( C(\beta) \) is given by
$$
@@ -865,7 +997,7 @@ This result implies that \( C(\beta) \) is a convex function since the matrix \(
Simple program
+Simple program
Gradient Descent Example
+Gradient Descent Example
And a corresponding example using scikit-learn
+And a corresponding example using scikit-learn
Gradient descent and Ridge
+Gradient descent and Ridge
Stochastic Gradient Descent
+Stochastic Gradient Descent
Computation of gradients
+Computation of gradients
SGD example
+SGD example
As an example, suppose we have \( 10 \) data points \( (\mathbf{x}_1,\cdots, \mathbf{x}_{10}) \)
and we choose to have \( M=5 \) minibathces,
then each minibatch contains two data points. In particular we have
@@ -1116,7 +1248,7 @@ $$
The gradient step
+The gradient step
Simple example code
+Simple example code
When do we stop?
+When do we stop?
Slightly different approach
+Slightly different approach
Conjugate gradient (CG) method
-
-$$
-\begin{equation*}
- \hat{A}\hat{x} = \hat{b}.
-\end{equation*}
-$$
-
-
-In the iterative process we end up with a problem like
-
-
-$$
-\begin{equation*}
- \hat{r}= \hat{b}-\hat{A}\hat{x},
-\end{equation*}
-$$
-
-
-where \( \hat{r} \) is the so-called residual or error in the iterative process.
-
-Conjugate gradient method
-
-$$
-\begin{equation*}
- P(\hat{x})=\frac{1}{2}\hat{x}^T\hat{A}\hat{x} - \hat{x}^T\hat{b},
-\end{equation*}
-$$
-
-
-with the constraint that the matrix \( \hat{A} \) is positive definite and symmetric.
-If we search for a minimum of the quantum mechanical variance, then the matrix
-\( \hat{A} \), which is called the Hessian, is given by the second-derivative of the function we want to minimize. This quantity is always positive definite. In our case this corresponds normally to the second derivative of the energy.
-Conjugate gradient method, Newton's method first
-
-$$
-\begin{equation*}
-\hat{x}-\hat{x}_i=\hat{A}^{-1}\nabla f(\hat{x}_i).
-\end{equation*}
-$$
-
-
-Subtracting this equation from that of \( \hat{x}_{i+1} \) we have
-
-$$
-\begin{equation*}
-\hat{x}_{i+1}-\hat{x}_i=\hat{A}^{-1}(\nabla f(\hat{x}_{i+1})-\nabla f(\hat{x}_i)).
-\end{equation*}
-$$
-
-Conjugate gradient method
-
-$$
-\begin{equation*}
-\hat{s}^T\hat{A}\hat{t}= 0.
-\end{equation*}
-$$
-
-
-The philosophy of the CG method is to perform searches in various conjugate directions
-of our vectors \( \hat{x}_i \) obeying the above criterion, namely
-
-$$
-\begin{equation*}
-\hat{x}_i^T\hat{A}\hat{x}_j= 0.
-\end{equation*}
-$$
-
-
-Two vectors are conjugate if they are orthogonal with respect to
-this inner product. Being conjugate is a symmetric relation: if \( \hat{s} \) is conjugate to \( \hat{t} \), then \( \hat{t} \) is conjugate to \( \hat{s} \).
-Conjugate gradient method
-
-$$
-\begin{equation*}
-\hat{v}_i^T\hat{A}\hat{v}_j= \lambda\hat{v}_i^T\hat{v}_j,
-\end{equation*}
-$$
-
-
-which is zero unless \( i=j \).
-Conjugate gradient method
-
-$$
-\begin{equation*}
-\hat{x}_{i+1}=\hat{x}_{i}+\alpha_i\hat{p}_{i}.
-\end{equation*}
-$$
-
-
-We assume that \( \hat{p}_{i} \) is a sequence of \( n \) mutually conjugate directions.
-Then the \( \hat{p}_{i} \) form a basis of \( R^n \) and we can expand the solution
-$ \hat{A}\hat{x} = \hat{b}$ in this basis, namely
-
-
-$$
-\begin{equation*}
- \hat{x} = \sum^{n}_{i=1} \alpha_i \hat{p}_i.
-\end{equation*}
-$$
-
-Conjugate gradient method
-
-$$
-\begin{equation*}
- \mathbf{A}\mathbf{x} = \sum^{n}_{i=1} \alpha_i \mathbf{A} \mathbf{p}_i = \mathbf{b}.
-\end{equation*}
-$$
-
-
-Multiplying with \( \hat{p}_k^T \) from the left gives
-
-
-$$
-\begin{equation*}
- \hat{p}_k^T \hat{A}\hat{x} = \sum^{n}_{i=1} \alpha_i\hat{p}_k^T \hat{A}\hat{p}_i= \hat{p}_k^T \hat{b},
-\end{equation*}
-$$
-
-
-and we can define the coefficients \( \alpha_k \) as
-
-
-$$
-\begin{equation*}
- \alpha_k = \frac{\hat{p}_k^T \hat{b}}{\hat{p}_k^T \hat{A} \hat{p}_k}
-\end{equation*}
-$$
-
-Conjugate gradient method and iterations
-
-$$
-\begin{equation*}
-\hat{x}_0=0,
-\end{equation*}
-$$
-
-
-or consider the system
-
-$$
-\begin{equation*}
-\hat{A}\hat{z} = \hat{b}-\hat{A}\hat{x}_0,
-\end{equation*}
-$$
-
-
-instead.
-Conjugate gradient method
-
-$$
-\begin{equation*}
- f(\hat{x}) = \frac{1}{2}\hat{x}^T\hat{A}\hat{x} - \hat{x}^T \hat{x} , \quad \hat{x}\in\mathbf{R}^n.
-\end{equation*}
-$$
-
-
-This suggests taking the first basis vector \( \hat{p}_1 \)
-to be the gradient of \( f \) at \( \hat{x}=\hat{x}_0 \),
-which equals
-
-$$
-\begin{equation*}
-\hat{A}\hat{x}_0-\hat{b},
-\end{equation*}
-$$
-
-
-and
-\( \hat{x}_0=0 \) it is equal \( -\hat{b} \).
-The other vectors in the basis will be conjugate to the gradient,
-hence the name conjugate gradient method.
-Conjugate gradient method
-
-$$
-\begin{equation*}
-\hat{r}_k=\hat{b}-\hat{A}\hat{x}_k.
-\end{equation*}
-$$
-
-
-Note that \( \hat{r}_k \) is the negative gradient of \( f \) at
-\( \hat{x}=\hat{x}_k \),
-so the gradient descent method would be to move in the direction \( \hat{r}_k \).
-Here, we insist that the directions \( \hat{p}_k \) are conjugate to each other,
-so we take the direction closest to the gradient \( \hat{r}_k \)
-under the conjugacy constraint.
-This gives the following expression
-
-$$
-\begin{equation*}
-\hat{p}_{k+1}=\hat{r}_k-\frac{\hat{p}_k^T \hat{A}\hat{r}_k}{\hat{p}_k^T\hat{A}\hat{p}_k} \hat{p}_k.
-\end{equation*}
-$$
-
-Conjugate gradient method
-
-$$
-\begin{equation*}
-\hat{r}_{k+1}=\hat{b}-\hat{A}\hat{x}_{k+1},
- \end{equation*}
-$$
-
-
-which equals
-
-$$
-\begin{equation*}
-\hat{b}-\hat{A}(\hat{x}_k+\alpha_k\hat{p}_k),
- \end{equation*}
-$$
-
-
-or
-
-$$
-\begin{equation*}
-(\hat{b}-\hat{A}\hat{x}_k)-\alpha_k\hat{A}\hat{p}_k,
- \end{equation*}
-$$
-
-
-which gives
-
-
-$$
-\begin{equation*}
-\hat{r}_{k+1}=\hat{r}_k-\hat{A}\hat{p}_{k},
- \end{equation*}
-$$
-
-
-Conjugate gradient method
+Gradient method
+
+Steepest descent method
+
+
-Simple codes for steepest descent and conjugate gradient using a \( 2\times 2 \) matrix, in c++, Python code to come
+Steepest descent method
+
+
+Gradient descent method
+
+
+Final expressions
+
+
+The Steepest descent algorithm
+
+
+
+Simple codes for steepest descent and conjugate gradient using a \( 2\times 2 \) matrix, in c++, Python code to come
-The routine for the steepest descent method
+The routine for the steepest descent method
Revisiting our first homework
+Revisiting our first homework
Gradient descent example
+Gradient descent example
-The derivative of the cost/loss function
+The derivative of the cost/loss function
-The Hessian matrix
+The Hessian matrix
The Hessian matrix of \( C(\beta) \) is given by
$$
\hat{H} \equiv \begin{bmatrix}
@@ -829,7 +935,7 @@ This result implies that \( C(\beta) \) is a convex function since the matrix \(
-Simple program
+Simple program
-Gradient Descent Example
+Gradient Descent Example
-And a corresponding example using scikit-learn
+And a corresponding example using scikit-learn
Gradient descent and Ridge
+Gradient descent and Ridge
-Stochastic Gradient Descent
+Stochastic Gradient Descent
-Computation of gradients
+Computation of gradients
-SGD example
+SGD example
As an example, suppose we have \( 10 \) data points \( (\mathbf{x}_1,\cdots, \mathbf{x}_{10}) \)
and we choose to have \( M=5 \) minibathces,
then each minibatch contains two data points. In particular we have
@@ -1062,7 +1168,7 @@ $$
-The gradient step
+The gradient step
-Simple example code
+Simple example code
-When do we stop?
+When do we stop?
-Slightly different approach
+Slightly different approach
-
-Conjugate gradient (CG) method
-
-
-Conjugate gradient method
-
-
-Conjugate gradient method, Newton's method first
-
-
-Conjugate gradient method
-
-
-Conjugate gradient method
-
-
-Conjugate gradient method
-
-
-Conjugate gradient method
-
-
-Conjugate gradient method and iterations
-
-
-Conjugate gradient method
-
-
-Conjugate gradient method
-
-
-Conjugate gradient method
-
-Conjugate gradient method
+Gradient method
+
+Steepest descent method
+
+
-Simple codes for steepest descent and conjugate gradient using a \( 2\times 2 \) matrix, in c++, Python code to come
+Steepest descent method
+
+
+Gradient descent method
+
+
+Final expressions
+
+
+The Steepest descent algorithm
+
+
+
+Simple codes for steepest descent and conjugate gradient using a \( 2\times 2 \) matrix, in c++, Python code to come
-The routine for the steepest descent method
+The routine for the steepest descent method
Revisiting our first homework
+Revisiting our first homework
Gradient descent example
+Gradient descent example
-The derivative of the cost/loss function
+The derivative of the cost/loss function
-The Hessian matrix
+The Hessian matrix
The Hessian matrix of \( C(\beta) \) is given by
$$
\hat{H} \equiv \begin{bmatrix}
@@ -834,7 +940,7 @@ This result implies that \( C(\beta) \) is a convex function since the matrix \(
-Simple program
+Simple program
-Gradient Descent Example
+Gradient Descent Example
-And a corresponding example using scikit-learn
+And a corresponding example using scikit-learn
Gradient descent and Ridge
+Gradient descent and Ridge
-Stochastic Gradient Descent
+Stochastic Gradient Descent
-Computation of gradients
+Computation of gradients
-SGD example
+SGD example
As an example, suppose we have \( 10 \) data points \( (\mathbf{x}_1,\cdots, \mathbf{x}_{10}) \)
and we choose to have \( M=5 \) minibathces,
then each minibatch contains two data points. In particular we have
@@ -1067,7 +1173,7 @@ $$
-The gradient step
+The gradient step
-Simple example code
+Simple example code
-When do we stop?
+When do we stop?
-Slightly different approach
+Slightly different approach
-
-Conjugate gradient (CG) method
-
-
-Conjugate gradient method
-
-
-Conjugate gradient method, Newton's method first
-
-
-Conjugate gradient method
-
-
-Conjugate gradient method
-
-
-Conjugate gradient method
-
-
-Conjugate gradient method
-
-
-Conjugate gradient method and iterations
-
-
-Conjugate gradient method
-
-
-Conjugate gradient method
-
-
-Conjugate gradient method
-