-The method of steepest descent The basic idea of gradient descent is
-that a function \( F(\mathbf{x}) \),
-\( \mathbf{x} \equiv (x_1,\cdots,x_n) \), decreases fastest if one goes from \( \bf {x} \) in the
-direction of the negative gradient \( -\nabla F(\mathbf{x}) \).
+In our discussion on Logistic Regression we defined we studied first the
+case of
+two classes, with \( y_i \) either
+\( 0 \) or \( 1 \). Furthermore we assumed also that we have only two
+parameters \( \beta \) in our fitting of the Sigmoid function, that is we
+defined probabilities
-
-It can be shown that if
$$
-\mathbf{x}_{k+1} = \mathbf{x}_k - \gamma_k \nabla F(\mathbf{x}_k),
+\begin{align*}
+p(y_i=1|x_i,\hat{\beta}) &= \frac{\exp{(\beta_0+\beta_1x_i)}}{1+\exp{(\beta_0+\beta_1x_i)}},\nonumber\\
+p(y_i=0|x_i,\hat{\beta}) &= 1 - p(y_i=1|x_i,\hat{\beta}),
+\end{align*}
$$
-with \( \gamma_k > 0 \).
-
-
-For \( \gamma_k \) small enough, then \( F(\mathbf{x}_{k+1}) \leq
-F(\mathbf{x}_k) \). This means that for a sufficiently small \( \gamma_k \)
-we are always moving towards smaller function values, i.e a minimum.
+where \( \hat{\beta} \) are the weights we wish to extract from data, in our case \( \beta_0 \) and \( \beta_1 \).
@@ -213,7 +243,7 @@ we are always moving towards smaller function values, i.e a minimum.
-The previous observation is the basis of the method of steepest
-descent, which is also referred to as just gradient descent (GD). One
-starts with an initial guess \( \mathbf{x}_0 \) for a minimum of \( F \) and
-computes new approximations according to
+Our compact equations used a definition of a vector \( \hat{y} \) with \( n \)
+elements \( y_i \), an \( n\times p \) matrix \( \hat{X} \) which contains the
+\( x_i \) values and a vector \( \hat{p} \) of fitted probabilities
+\( p(y_i\vert x_i,\hat{\beta}) \). We rewrote in a more compact form
+the first derivative of the cost function as
$$
-\mathbf{x}_{k+1} = \mathbf{x}_k - \gamma_k \nabla F(\mathbf{x}_k), \ \ k \geq 0.
+\frac{\partial \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}} = -\hat{X}^T\left(\hat{y}-\hat{p}\right).
$$
-The parameter \( \gamma_k \) is often referred to as the step length or
-the learning rate within the context of Machine Learning.
+If we in addition define a diagonal matrix \( \hat{W} \) with elements
+\( p(y_i\vert x_i,\hat{\beta})(1-p(y_i\vert x_i,\hat{\beta}) \), we can obtain a compact expression of the second derivative as
+
+$$
+\frac{\partial^2 \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}\partial \hat{\beta}^T} = \hat{X}^T\hat{W}\hat{X}.
+$$
+
+This defines what we call the Hessian.
@@ -209,7 +248,7 @@ the learning rate within the context of Machine Learning.
-Ideally the sequence \( \{\mathbf{x}_k \}_{k=0} \) converges to a global
-minimum of the function \( F \). In general we do not know if we are in a
-global or local minimum. In the special case when \( F \) is a convex
-function, all local minima are also global minima, so in this case
-gradient descent can converge to the global solution. The advantage of
-this scheme is that it is conceptually simple and straightforward to
-implement. However the method in this form has some severe
-limitations:
+If we can set up these equations, Newton-Raphson's iterative method is the nomrally the method of choice. It requires however that we setting the matrices that define the first and second derivatives.
-In machine learing we are often faced with non-convex high dimensional
-cost functions with many local minima. Since GD is deterministic we
-will get stuck in a local minimum, if the method converges, unless we
-have a very good intial guess. This also implies that the scheme is
-sensitive to the chosen initial condition.
+Our iterative scheme is then given by
+
+$$
+\hat{\beta}^{\mathrm{new}} = \hat{\beta}^{\mathrm{old}}-\left(\frac{\partial^2 \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}\partial \hat{\beta}^T}\right)^{-1}\times \left(\frac{\partial \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}}\right)_{\hat{\beta}^{\mathrm{old}}},
+$$
+
+or in matrix form as
+
+$$
+\hat{\beta}^{\mathrm{new}} = \hat{\beta}^{\mathrm{old}}-\left(\hat{X}^T\hat{W}\hat{X} \right)^{-1}\times \left(-\hat{X}^T(\hat{y}-\hat{p}) \right)_{\hat{\beta}^{\mathrm{old}}}.
+$$
+
+The right-hand side is computed with the old values of \( \beta \).
-Note that the gradient is a function of \( \mathbf{x} =
-(x_1,\cdots,x_n) \) which makes it expensive to compute numerically.
+If we can compute these matrices, in particular the Hessian, the above is often the easiest method to implement.
@@ -217,7 +249,7 @@ Note that the gradient is a function of \( \mathbf{x} =
-The gradient descent method
-is sensitive to the choice of learning rate \( \gamma_k \). This is due
-to the fact that we are only guaranteed that \( F(\mathbf{x}_{k+1}) \leq
-F(\mathbf{x}_k) \) for sufficiently small \( \gamma_k \). The problem is to
-determine an optimal learning rate. If the learning rate is chosen too
-small the method will take a long time to converge and if it is too
-large we can experience erratic behavior.
+Let us quicly remind ourselves how we derive the above method.
-Many of these shortcomings can be alleviated by introducing
-randomness. One such method is that of Stochastic Gradient Descent
-(SGD), see below.
+Perhaps the most celebrated of all one-dimensional root-finding
+routines is Newton's method, also called the Newton-Raphson
+method. This method is distinguished from the previously discussed
+methods by the fact that it requires the evaluation of both the
+function \( f \) and its derivative \( f' \) at arbitrary points. In this
+sense, it is taylored to cases with e.g., transcendental equations.
+If you can only calculate the derivative
+numerically and/or your function is not of the smooth type, we
+discourage the use of this method.
@@ -211,7 +243,7 @@ randomness. One such method is that of Stochastic Gradient Descent
-Ideally we want our cost/loss function to be convex(concave).
+The Newton-Raphson formula consists geometrically of extending the
+tangent line at a current point until it crosses zero, then setting
+the next guess to the abscissa of that zero-crossing. The mathematics
+behind this method is rather simple. Employing a Taylor expansion for
+\( x \) sufficiently close to the solution \( s \), we have
+
+$$
+ f(s)=0=f(x)+(s-x)f'(x)+\frac{(s-x)^2}{2}f''(x) +\dots.
+ \tag{1}
+$$
-First we give the definition of a convex set: A set \( C \) in
-\( \mathbb{R}^n \) is said to be convex if, for all \( x \) and \( y \) in \( C \) and
-all \( t \in (0,1) \) , the point \( (1 − t)x + ty \) also belongs to
-C. Geometrically this means that every point on the line segment
-connecting \( x \) and \( y \) is in \( C \) as discussed below.
+For small enough values of the function and for well-behaved
+functions, the terms beyond linear are unimportant, hence we obtain
+
+$$
+ f(x)+(s-x)f'(x)\approx 0,
+$$
+
+yielding
+$$
+ s\approx x-\frac{f(x)}{f'(x)}.
+$$
-The convex subsets of \( \mathbb{R} \) are the intervals of
-\( \mathbb{R} \). Examples of convex sets of \( \mathbb{R}^2 \) are the
-regular polygons (triangles, rectangles, pentagons, etc...).
+Having in mind an iterative procedure, it is natural to start iterating with
+$$
+ x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}.
+$$
-Convex function: Let \( X \subset \mathbb{R}^n \) be a convex set. Assume that the function \( f: X \rightarrow \mathbb{R} \) is continuous, then \( f \) is said to be convex if $$f(tx_1 + (1-t)x_2) \leq tf(x_1) + (1-t)f(x_2) $$ for all \( x_1, x_2 \in X \) and for all \( t \in [0,1] \). If \( \leq \) is replaced with a strict inequaltiy in the definition, we demand \( x_1 \neq x_2 \) and \( t\in(0,1) \) then \( f \) is said to be strictly convex. For a single variable function, convexity means that if you draw a straight line connecting \( f(x_1) \) and \( f(x_2) \), the value of the function on the interval \( [x_1,x_2] \) is always below the line as illustrated below.
+The above is Newton-Raphson's method. It has a simple geometric
+interpretation, namely \( x_{n+1} \) is the point where the tangent from
+\( (x_n,f(x_n)) \) crosses the $x-$axis. Close to the solution,
+Newton-Raphson converges fast to the desired result. However, if we
+are far from a root, where the higher-order terms in the series are
+important, the Newton-Raphson formula can give grossly inaccurate
+results. For instance, the initial guess for the root might be so far
+from the true root as to let the search interval include a local
+maximum or minimum of the function. If an iteration places a trial
+guess near such a local extremum, so that the first derivative nearly
+vanishes, then Newton-Raphson may fail totally
-In the following we state first and second-order conditions which
-ensures convexity of a function \( f \). We write \( D_f \) to denote the
-domain of \( f \), i.e the subset of \( R^n \) where \( f \) is defined. For more
-details and proofs we refer to: S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press.
+Newton's method can be generalized to systems of several non-linear equations
+and variables. Consider the case with two equations
+$$
+ \begin{array}{cc} f_1(x_1,x_2) &=0\\
+ f_2(x_1,x_2) &=0\end{array},
+$$
+
+which we Taylor expand to obtain
+
+$$
+ \begin{array}{cc} 0=f_1(x_1+h_1,x_2+h_2)=&f_1(x_1,x_2)+h_1
+ \partial f_1/\partial x_1+h_2
+ \partial f_1/\partial x_2+\dots\\
+ 0=f_2(x_1+h_1,x_2+h_2)=&f_2(x_1,x_2)+h_1
+ \partial f_2/\partial x_1+h_2
+ \partial f_2/\partial x_2+\dots
+ \end{array}.
+$$
+
+Defining the Jacobian matrix \( {\bf \hat{J}} \) we have
+$$
+ {\bf \hat{J}}=\left( \begin{array}{cc}
+ \partial f_1/\partial x_1 & \partial f_1/\partial x_2 \\
+ \partial f_2/\partial x_1 &\partial f_2/\partial x_2
+ \end{array} \right),
+$$
+
+we can rephrase Newton's method as
+$$
+\left(\begin{array}{c} x_1^{n+1} \\ x_2^{n+1} \end{array} \right)=
+\left(\begin{array}{c} x_1^{n} \\ x_2^{n} \end{array} \right)+
+\left(\begin{array}{c} h_1^{n} \\ h_2^{n} \end{array} \right),
+$$
+
+where we have defined
+$$
+ \left(\begin{array}{c} h_1^{n} \\ h_2^{n} \end{array} \right)=
+ -{\bf \hat{J}}^{-1}
+ \left(\begin{array}{c} f_1(x_1^{n},x_2^{n}) \\ f_2(x_1^{n},x_2^{n}) \end{array} \right).
+$$
+
+We need thus to compute the inverse of the Jacobian matrix and it
+is to understand that difficulties may
+arise in case \( {\bf \hat{J}} \) is nearly singular.
-
-
-
-Suppose \( f \) is differentiable (i.e \( \nabla f(x) \) is well defined for
-all \( x \) in the domain of \( f \)). Then \( f \) is convex if and only if \( D_f \)
-is a convex set and $$f(y) \geq f(x) + \nabla f(x)^T (y-x) $$ holds
-for all \( x,y \in D_f \). This condition means that for a convex function
-the first order Taylor expansion (right hand side above) at any point
-a global under estimator of the function. To convince yourself you can
-make a drawing of \( f(x) = x^2+1 \) and draw the tangent line to \( f(x) \) and
-note that it is always below the graph.
-
-
-
-
-
-
-
-
-Assume that \( f \) is twice
-differentiable, i.e the Hessian matrix exists at each point in
-\( D_f \). Then \( f \) is convex if and only if \( D_f \) is a convex set and its
-Hessian is positive semi-definite for all \( x\in D_f \). For a
-single-variable function this reduces to \( f''(x) \geq 0 \). Geometrically this means that \( f \) has nonnegative curvature
-everywhere.
-
-
-
-
-
-This condition is particularly useful since it gives us an procedure for determining if the function under consideration is convex, apart from using the definition.
+It is rather straightforward to extend the above scheme to systems of
+more than two non-linear equations. In our case, the Jacobian matrix is given by the Hessian that represents the second derivative of cost function.
@@ -239,7 +282,7 @@ This condition is particularly useful since it gives us an procedure for determi
-The next result is of great importance to us and the reason why we are
-going on about convex functions. In machine learning we frequently
-have to minimize a loss/cost function in order to find the best
-parameters for the model we are considering.
+The method of steepest descent The basic idea of gradient descent is
+that a function \( F(\mathbf{x}) \),
+\( \mathbf{x} \equiv (x_1,\cdots,x_n) \), decreases fastest if one goes from \( \bf {x} \) in the
+direction of the negative gradient \( -\nabla F(\mathbf{x}) \).
-Ideally we want the
-global minimum (for high-dimensional models it is hard to know
-if we have local or global minimum). However, if the cost/loss function
-is convex the following result provides invaluable information:
+It can be shown that if
+$$
+\mathbf{x}_{k+1} = \mathbf{x}_k - \gamma_k \nabla F(\mathbf{x}_k),
+$$
+
+with \( \gamma_k > 0 \).
-
-
-
-Consider the problem of finding \( x \in \mathbb{R}^n \) such that \( f(x) \)
-is minimal, where \( f \) is convex and differentiable. Then, any point
-\( x^* \) that satisfies \( \nabla f(x^*) = 0 \) is a global minimum.
-
-
-
-
-
-This result means that if we know that the cost/loss function is convex and we are able to find a minimum, we are guaranteed that it is a global minimum.
+For \( \gamma_k \) small enough, then \( F(\mathbf{x}_{k+1}) \leq
+F(\mathbf{x}_k) \). This means that for a sufficiently small \( \gamma_k \)
+we are always moving towards smaller function values, i.e a minimum.
@@ -227,7 +252,7 @@ This result means that if we know that the cost/loss function is convex and we a
Show that \( f(x)=x^2 \) is convex for \( x \in \mathbb{R} \) using the definition of convexity. Hint: If you re-write the definition, \( f \) is convex if the following holds for all \( x,y \in D_f \) and any \( \lambda \in [0,1] \) $\lambda f(x)+(1-\lambda)f(y)-f(\lambda x + (1-\lambda) y ) \geq 0$.
-
Using the second order condition show that the following functions are convex on the specified domain.
+
+The previous observation is the basis of the method of steepest
+descent, which is also referred to as just gradient descent (GD). One
+starts with an initial guess \( \mathbf{x}_0 \) for a minimum of \( F \) and
+computes new approximations according to
-
-
\( f(x) = e^x \) is convex for \( x \in \mathbb{R} \).
-
\( g(x) = -\ln(x) \) is convex for \( x \in (0,\infty) \).
Let \( f(x) = x^2 \) and \( g(x) = e^x \). Show that \( f(g(x)) \) and \( g(f(x)) \) is convex for \( x \in \mathbb{R} \). Also show that if \( f(x) \) is any convex function than \( h(x) = e^{f(x)} \) is convex.
-
A norm is any function that satisfy the following properties
-
-
-
\( f(\alpha x) = |\alpha| f(x) \) for all \( \alpha \in \mathbb{R} \).
-
\( f(x+y) \leq f(x) + f(y) \)
-
\( f(x) \leq 0 \) for all \( x \in \mathbb{R}^n \) with equality if and only if \( x = 0 \)
-
-
-
-
-Using the definition of convexity, try to show that a function satisfying the properties above is convex (the third condition is not needed to show this).
+
+The parameter \( \gamma_k \) is often referred to as the step length or
+the learning rate within the context of Machine Learning.
@@ -224,7 +248,7 @@ Using the definition of convexity, try to show that a function satisfying the pr
-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:
+Ideally the sequence \( \{\mathbf{x}_k \}_{k=0} \) converges to a global
+minimum of the function \( F \). In general we do not know if we are in a
+global or local minimum. In the special case when \( F \) is a convex
+function, all local minima are also global minima, so in this case
+gradient descent can converge to the global solution. The advantage of
+this scheme is that it is conceptually simple and straightforward to
+implement. However the method in this form has some severe
+limitations:
-
-
An analytical solution (recall homework set 1).
-
The gradient can be computed analytically.
-
The cost function is convex which guarantees that gradient descent converges for small enough learning rates
-
+
+In machine learing we are often faced with non-convex high dimensional
+cost functions with many local minima. Since GD is deterministic we
+will get stuck in a local minimum, if the method converges, unless we
+have a very good intial guess. This also implies that the scheme is
+sensitive to the chosen initial condition.
-We revisit the example from homework set 1 where we had
-$$
-y_i = 5x_i^2 + 0.1\xi_i, \ i=1,\cdots,100
-$$
-
-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,
-$$
-
-such that
-$$
-\hat{y}_i = \beta_0 + \beta_1 x_i.
-$$
+
+Note that the gradient is a function of \( \mathbf{x} =
+(x_1,\cdots,x_n) \) which makes it expensive to compute numerically.
-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 \)
+The gradient descent method
+is sensitive to the choice of learning rate \( \gamma_k \). This is due
+to the fact that we are only guaranteed that \( F(\mathbf{x}_{k+1}) \leq
+F(\mathbf{x}_k) \) for sufficiently small \( \gamma_k \). The problem is to
+determine an optimal learning rate. If the learning rate is chosen too
+small the method will take a long time to converge and if it is too
+large we can experience erratic behavior.
-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}.
-$$
-
-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.
+Many of these shortcomings can be alleviated by introducing
+randomness. One such method is that of Stochastic Gradient Descent
+(SGD), see below.
@@ -222,7 +248,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}),
-$$
+Ideally we want our cost/loss function to be convex(concave).
-where \( X \) is the design matrix defined above.
+
+First we give the definition of a convex set: A set \( C \) in
+\( \mathbb{R}^n \) is said to be convex if, for all \( x \) and \( y \) in \( C \) and
+all \( t \in (0,1) \) , the point \( (1 − t)x + ty \) also belongs to
+C. Geometrically this means that every point on the line segment
+connecting \( x \) and \( y \) is in \( C \) as discussed below.
+
+
+The convex subsets of \( \mathbb{R} \) are the intervals of
+\( \mathbb{R} \). Examples of convex sets of \( \mathbb{R}^2 \) are the
+regular polygons (triangles, rectangles, pentagons, etc...).
@@ -212,7 +249,7 @@ where \( X \) is the design matrix defined above.
-This result implies that \( C(\beta) \) is a convex function since the matrix \( X^T X \) always is positive semi-definite.
+
+Convex function: Let \( X \subset \mathbb{R}^n \) be a convex set. Assume that the function \( f: X \rightarrow \mathbb{R} \) is continuous, then \( f \) is said to be convex if $$f(tx_1 + (1-t)x_2) \leq tf(x_1) + (1-t)f(x_2) $$ for all \( x_1, x_2 \in X \) and for all \( t \in [0,1] \). If \( \leq \) is replaced with a strict inequaltiy in the definition, we demand \( x_1 \neq x_2 \) and \( t\in(0,1) \) then \( f \) is said to be strictly convex. For a single variable function, convexity means that if you draw a straight line connecting \( f(x_1) \) and \( f(x_2) \), the value of the function on the interval \( [x_1,x_2] \) is always below the line as illustrated below.
@@ -211,7 +237,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
-$$
+In the following we state first and second-order conditions which
+ensures convexity of a function \( f \). We write \( D_f \) to denote the
+domain of \( f \), i.e the subset of \( R^n \) where \( f \) is defined. For more
+details and proofs we refer to: S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press.
-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} \).
+
+
+
+Suppose \( f \) is differentiable (i.e \( \nabla f(x) \) is well defined for
+all \( x \) in the domain of \( f \)). Then \( f \) is convex if and only if \( D_f \)
+is a convex set and $$f(y) \geq f(x) + \nabla f(x)^T (y-x) $$ holds
+for all \( x,y \in D_f \). This condition means that for a convex function
+the first order Taylor expansion (right hand side above) at any point
+a global under estimator of the function. To convince yourself you can
+make a drawing of \( f(x) = x^2+1 \) and draw the tangent line to \( f(x) \) and
+note that it is always below the graph.
+
+
+
-And finally we can compare our solution for \( \beta \) with the analytic result given by
-\( \beta= (X^TX)^{-1} X^T \mathbf{y} \).
+
+
+
+Assume that \( f \) is twice
+differentiable, i.e the Hessian matrix exists at each point in
+\( D_f \). Then \( f \) is convex if and only if \( D_f \) is a convex set and its
+Hessian is positive semi-definite for all \( x\in D_f \). For a
+single-variable function this reduces to \( f''(x) \geq 0 \). Geometrically this means that \( f \) has nonnegative curvature
+everywhere.
+
+
+
+
+This condition is particularly useful since it gives us an procedure for determining if the function under consideration is convex, apart from using the definition.
-
-
importnumpyasnp
-
-"""
-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)
-
-Another simple example is here
+The next result is of great importance to us and the reason why we are
+going on about convex functions. In machine learning we frequently
+have to minimize a loss/cost function in order to find the best
+parameters for the model we are considering.
+
+Ideally we want the
+global minimum (for high-dimensional models it is hard to know
+if we have local or global minimum). However, if the cost/loss function
+is convex the following result provides invaluable information:
-
-
# Importing various packages
-fromrandomimport random, seed
-importnumpyasnp
-importmatplotlib.pyplotasplt
-frommpl_toolkits.mplot3dimport Axes3D
-frommatplotlibimport cm
-frommatplotlib.tickerimport LinearLocator, FormatStrFormatter
-importsys
+
+
+
+
+Consider the problem of finding \( x \in \mathbb{R}^n \) such that \( f(x) \)
+is minimal, where \( f \) is convex and differentiable. Then, any point
+\( x^* \) that satisfies \( \nabla f(x^*) = 0 \) is a global minimum.
+
+This result means that if we know that the cost/loss function is convex and we are able to find a minimum, we are guaranteed that it is a global minimum.
-eta =0.1
-Niterations =1000
-m =100
-
-foriterinrange(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()
-
Show that \( f(x)=x^2 \) is convex for \( x \in \mathbb{R} \) using the definition of convexity. Hint: If you re-write the definition, \( f \) is convex if the following holds for all \( x,y \in D_f \) and any \( \lambda \in [0,1] \) $\lambda f(x)+(1-\lambda)f(y)-f(\lambda x + (1-\lambda) y ) \geq 0$.
+
Using the second order condition show that the following functions are convex on the specified domain.
Let \( f(x) = x^2 \) and \( g(x) = e^x \). Show that \( f(g(x)) \) and \( g(f(x)) \) is convex for \( x \in \mathbb{R} \). Also show that if \( f(x) \) is any convex function than \( h(x) = e^{f(x)} \) is convex.
+
A norm is any function that satisfy the following properties
+
+
+
\( f(\alpha x) = |\alpha| f(x) \) for all \( \alpha \in \mathbb{R} \).
+
\( f(x+y) \leq f(x) + f(y) \)
+
\( f(x) \leq 0 \) for all \( x \in \mathbb{R}^n \) with equality if and only if \( x = 0 \)
+
+
+
+
+Using the definition of convexity, try to show that a function satisfying the properties above is convex (the third condition is not needed to show this).
-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_)
-
-We have also discussed Ridge regression where the loss function contains a regularized given by the \( L_2 \) norm of \( \beta \),
+Before we proceed, we would like to mention the approach called the standard Steepest descent, which again leads to us having to be able to compute a matrix.
+
+
+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
$$
-C_{\text{ridge}}(\beta) = ||X\beta -\mathbf{y}||^2 + \lambda ||\beta||^2, \ \lambda \geq 0.
+\begin{equation*}
+\hat{A}\hat{x} = \hat{b}.
+\end{equation*}
$$
-In order to minimize \( C_{\text{ridge}}(\beta) \) using GD we only have adjust the gradient as follows
+In the iterative process we end up with a problem like
+
$$
-\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).
+\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.
+
-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},
-$$
+When we have found the exact solution, \( \hat{r}=0 \).
-for \( \lambda = {0,1,10,50,100} \) (\( \lambda = 0 \) corresponds to ordinary least squares).
-We can then compute \( ||\beta_{\text{ridge}}|| \) for each \( \lambda \).
-
-
-
-
-
importnumpyasnp
-
-"""
-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 (SGD) and variants thereof address some of
-the shortcomings of the Gradient descent method discussed above.
+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*}
+$$
-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}).
-$$
+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.
+
+
Simple codes for steepest descent and conjugate gradient using a \( 2\times 2 \) matrix, in c++, Python code to come
+
+
+
-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}).
-$$
+
+
#include<cmath>
+#include<iostream>
+#include<fstream>
+#include<iomanip>
+#include"vectormatrixclass.h"
+usingnamespace std;
+// Main function begins here
+intmain(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;
+}
+
-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 \).
+
+
+
@@ -217,7 +274,7 @@ minibatches. We denote these minibatches by \( B_k \) where
-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 \).
-
+
The routine for the steepest descent method
+
+
+
-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}).
-$$
+
+
+
Vector SteepestDescent(Matrix A, Vector b, Vector x0){
+ int IterMax, i;
+ int dim = x0.Dimension();
+ constdouble 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;
+}
+
-Thus a gradient descent step now looks like
+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:
+
+
+
An analytical solution (recall homework set 1).
+
The gradient can be computed analytically.
+
The cost function is convex which guarantees that gradient descent converges for small enough learning rates
+
+
+We revisit the example from homework set 1 where we had
$$
-\beta_{j+1} = \beta_j - \gamma_j \sum_{i \in B_k}^n \nabla_\beta c_i(\mathbf{x}_i,
-\mathbf{\beta})
+y_i = 5x_i^2 + 0.1\xi_i, \ i=1,\cdots,100
$$
-
-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.
+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,
+$$
+
+such that
+$$
+\hat{y}_i = \beta_0 + \beta_1 x_i.
+$$
@@ -216,7 +262,7 @@ the number of minibatches, as exemplified in the code below.
importnumpyasnp
-
-n =100#100 datapoints
-M =5#size of each minibatch
-m =int(n/M) #number of minibatches
-n_epochs =10#number of epochs
-
-j =0
-for epoch inrange(1,n_epochs+1):
- for i inrange(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
-
-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.
+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}.
+$$
+
+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.
-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.
+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.
-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 result implies that \( C(\beta) \) is a convex function since the matrix \( X^T X \) always is positive semi-definite.
-
-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 \).
-
-
-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.
-
-
-
-
-
importnumpyasnp
-
-defstep_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 inrange(1,n_epochs+1):
- for i inrange(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))
-
-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.
+
Simple program
-When we have found the exact solution, \( \hat{r}=0 \).
-
-
+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} \).
+
+
+
+
importnumpyasnp
+
+"""
+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)
+
@@ -229,7 +270,7 @@ 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*}
-$$
+Another simple example is here
+
-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*}
-$$
+
And a corresponding example using scikit-learn
-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
+
Gradient descent and Ridge
+
+
+We have also discussed Ridge regression where the loss function contains a regularized given by the \( L_2 \) norm of \( \beta \),
$$
-\begin{equation*}
-\hat{s}^T\hat{A}\hat{t}= 0.
-\end{equation*}
+C_{\text{ridge}}(\beta) = ||X\beta -\mathbf{y}||^2 + \lambda ||\beta||^2, \ \lambda \geq 0.
$$
-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
+
+In order to minimize \( C_{\text{ridge}}(\beta) \) using GD we only have adjust the gradient as follows
$$
-\begin{equation*}
-\hat{x}_i^T\hat{A}\hat{x}_j= 0.
-\end{equation*}
+\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).
$$
-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} \).
-
-
+
+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 \).
+
+
+
+
+
importnumpyasnp
+
+"""
+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||
+
@@ -225,6 +283,10 @@ this inner product. Being conjugate is a symmetric relation: if \( \hat{s} \) is
-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*}
-$$
+
Stochastic Gradient Descent
-which is zero unless \( i=j \).
-
-
+
+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}).
+$$
@@ -212,6 +243,11 @@ 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
+
Computation of gradients
+
+
+This in turn means that the gradient can be
+computed as a sum over \( i \)-gradients
$$
-\begin{equation*}
-\hat{x}_{i+1}=\hat{x}_{i}+\alpha_i\hat{p}_{i}.
-\end{equation*}
+\nabla_\beta C(\mathbf{\beta}) = \sum_i^n \nabla_\beta c_i(\mathbf{x}_i,
+\mathbf{\beta}).
$$
-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*}
-$$
-
-
-
+
+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 \).
-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
+
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 \).
+
+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
$$
-\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*}
+\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}).
$$
-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.
+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})
+$$
-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.
-
-
-
+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.
-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*}
-$$
+
Simple example code
-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.
-
-
+
+
importnumpyasnp
+n =100#100 datapoints
+M =5#size of each minibatch
+m =int(n/M) #number of minibatches
+n_epochs =10#number of epochs
+
+j =0
+for epoch inrange(1,n_epochs+1):
+ for i inrange(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
+
+
+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.
@@ -220,6 +253,15 @@ 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*}
-$$
-
-
+
When do we stop?
+
+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.
+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.
-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*}
-$$
-
-
+
+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 \).
+
+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.
+
+
importnumpyasnp
+
+defstep_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 inrange(1,n_epochs+1):
+ for i inrange(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))
+
[2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University
-
Sep 21, 2018
+
Sep 27, 2018
@@ -174,7 +174,237 @@ some approximative/numerical method to compute the minimum.
-
Steepest descent
+
Revisiting our Logistic Regression case
+
+
+In our discussion on Logistic Regression we defined we studied first the
+case of
+two classes, with \( y_i \) either
+\( 0 \) or \( 1 \). Furthermore we assumed also that we have only two
+parameters \( \beta \) in our fitting of the Sigmoid function, that is we
+defined probabilities
+
+
+
+where \( \hat{\beta} \) are the weights we wish to extract from data, in our case \( \beta_0 \) and \( \beta_1 \).
+
+
+
+
+
The equations to solve
+
+
+Our compact equations used a definition of a vector \( \hat{y} \) with \( n \)
+elements \( y_i \), an \( n\times p \) matrix \( \hat{X} \) which contains the
+\( x_i \) values and a vector \( \hat{p} \) of fitted probabilities
+\( p(y_i\vert x_i,\hat{\beta}) \). We rewrote in a more compact form
+the first derivative of the cost function as
+
+
+If we in addition define a diagonal matrix \( \hat{W} \) with elements
+\( p(y_i\vert x_i,\hat{\beta})(1-p(y_i\vert x_i,\hat{\beta}) \), we can obtain a compact expression of the second derivative as
+
+
+If we can set up these equations, Newton-Raphson's iterative method is the nomrally the method of choice. It requires however that we setting the matrices that define the first and second derivatives.
+
+
+
+The right-hand side is computed with the old values of \( \beta \).
+
+
+If we can compute these matrices, in particular the Hessian, the above is often the easiest method to implement.
+
+
+
+
+
Brief reminder on Newton-Raphson's method
+
+
+Let us quicly remind ourselves how we derive the above method.
+
+
+Perhaps the most celebrated of all one-dimensional root-finding
+routines is Newton's method, also called the Newton-Raphson
+method. This method is distinguished from the previously discussed
+methods by the fact that it requires the evaluation of both the
+function \( f \) and its derivative \( f' \) at arbitrary points. In this
+sense, it is taylored to cases with e.g., transcendental equations.
+If you can only calculate the derivative
+numerically and/or your function is not of the smooth type, we
+discourage the use of this method.
+
+
+
+
+
The equations
+
+
+The Newton-Raphson formula consists geometrically of extending the
+tangent line at a current point until it crosses zero, then setting
+the next guess to the abscissa of that zero-crossing. The mathematics
+behind this method is rather simple. Employing a Taylor expansion for
+\( x \) sufficiently close to the solution \( s \), we have
+
+
+For small enough values of the function and for well-behaved
+functions, the terms beyond linear are unimportant, hence we obtain
+
+
+$$
+ f(x)+(s-x)f'(x)\approx 0,
+$$
+
+
+yielding
+
+$$
+ s\approx x-\frac{f(x)}{f'(x)}.
+$$
+
+
+
+Having in mind an iterative procedure, it is natural to start iterating with
+
+$$
+ x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}.
+$$
+
+
+
+
+
+
Simple geometric interpretation
+
+
+The above is Newton-Raphson's method. It has a simple geometric
+interpretation, namely \( x_{n+1} \) is the point where the tangent from
+\( (x_n,f(x_n)) \) crosses the $x-$axis. Close to the solution,
+Newton-Raphson converges fast to the desired result. However, if we
+are far from a root, where the higher-order terms in the series are
+important, the Newton-Raphson formula can give grossly inaccurate
+results. For instance, the initial guess for the root might be so far
+from the true root as to let the search interval include a local
+maximum or minimum of the function. If an iteration places a trial
+guess near such a local extremum, so that the first derivative nearly
+vanishes, then Newton-Raphson may fail totally
+
+
+
+
+
Extending to more than one variable
+
+
+Newton's method can be generalized to systems of several non-linear equations
+and variables. Consider the case with two equations
+
+
+We need thus to compute the inverse of the Jacobian matrix and it
+is to understand that difficulties may
+arise in case \( {\bf \hat{J}} \) is nearly singular.
+
+
+It is rather straightforward to extend the above scheme to systems of
+more than two non-linear equations. In our case, the Jacobian matrix is given by the Hessian that represents the second derivative of cost function.
+
+
+
+
+
Steepest descent
The method of steepest descent The basic idea of gradient descent is
@@ -200,7 +430,7 @@ we are always moving towards smaller function values, i.e a minimum.
-
More on Steepest descent
+
More on Steepest descent
The previous observation is the basis of the method of steepest
@@ -221,7 +451,7 @@ the learning rate within the context of Machine Learning.
-
The ideal
+
The ideal
Ideally the sequence \( \{\mathbf{x}_k \}_{k=0} \) converges to a global
@@ -247,7 +477,7 @@ Note that the gradient is a function of \( \mathbf{x} =
-
The sensitiveness of the gradient descent
+
The sensitiveness of the gradient descent
The gradient descent method
@@ -266,7 +496,7 @@ randomness. One such method is that of Stochastic Gradient Descent
-
Convex functions
+
Convex functions
Ideally we want our cost/loss function to be convex(concave).
@@ -286,7 +516,7 @@ regular polygons (triangles, rectangles, pentagons, etc...).
-
Convex function
+
Convex function
Convex function: Let \( X \subset \mathbb{R}^n \) be a convex set. Assume that the function \( f: X \rightarrow \mathbb{R} \) is continuous, then \( f \) is said to be convex if
In the following we state first and second-order conditions which
@@ -338,7 +568,7 @@ This condition is particularly useful since it gives us an procedure for determi
-
More on convex functions
+
More on convex functions
The next result is of great importance to us and the reason why we are
@@ -367,7 +597,7 @@ This result means that if we know that the cost/loss function is convex and we a
-
Some simple problems
+
Some simple problems
Show that \( f(x)=x^2 \) is convex for \( x \in \mathbb{R} \) using the definition of convexity. Hint: If you re-write the definition, \( f \) is convex if the following holds for all \( x,y \in D_f \) and any \( \lambda \in [0,1] \) $\lambda f(x)+(1-\lambda)f(y)-f(\lambda x + (1-\lambda) y ) \geq 0$.
@@ -394,7 +624,146 @@ Using the definition of convexity, try to show that a function satisfying the pr
-
Revisiting our first homework
+
Standard steepest descent
+
+
+Before we proceed, we would like to mention the approach called the standard Steepest descent, which again leads to us having to be able to compute a matrix.
+
+
+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
+
+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.
+
+
+More details will be added here soon.
+
+
+
+
+
Simple codes for steepest descent and conjugate gradient using a \( 2\times 2 \) matrix, in c++, Python code to come
+
+
+
+
+
+
#include<cmath>
+#include<iostream>
+#include<fstream>
+#include<iomanip>
+#include"vectormatrixclass.h"
+usingnamespace std;
+// Main function begins here
+intmain(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;
+}
+
+
+
+
+
+
+
+
The routine for the steepest descent method
+
+
+
+
+
+
Vector SteepestDescent(Matrix A, Vector b, Vector x0){
+ int IterMax, i;
+ int dim = x0.Dimension();
+ constdouble 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;
+}
+
+
+
+
+
+
+
+
Revisiting our first homework
We will use linear regression as a case study for the gradient descent
@@ -434,7 +803,7 @@ $$
-
Gradient descent example
+
Gradient descent example
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 \)
@@ -463,7 +832,7 @@ and we want to find \( \beta \) such that \( C(\beta) \) is minimized.
-
The derivative of the cost/loss function
+
The derivative of the cost/loss function
Computing \( \partial C(\beta) / \partial \beta_0 \) and \( \partial C(\beta) / \partial \beta_1 \) we can show that the gradient can be written as
@@ -480,7 +849,7 @@ where \( X \) is the design matrix defined above.
-
The Hessian matrix
+
The Hessian matrix
The Hessian matrix of \( C(\beta) \) is given by
$$
@@ -496,7 +865,7 @@ This result implies that \( C(\beta) \) is a convex function since the matrix \(
-
Simple program
+
Simple program
We can now write a program that minimizes \( C(\beta) \) using the gradient descent method with a constant learning rate \( \gamma \) according to
@@ -540,7 +909,7 @@ beta_NE = np.dot(Xt_X_inv,Xt_y)
-
Gradient Descent Example
+
Gradient Descent Example
Another simple example is here
@@ -590,7 +959,7 @@ plt.show()
-
And a corresponding example using scikit-learn
+
And a corresponding example using scikit-learn
@@ -615,7 +984,7 @@ sgdreg.fit(x,y.ravel())
-
Gradient descent and Ridge
+
Gradient descent and Ridge
We have also discussed Ridge regression where the loss function contains a regularized given by the \( L_2 \) norm of \( \beta \),
@@ -679,7 +1048,7 @@ beta_ridge = np.dot(Z,np.dot(X.T,y))
-
Stochastic Gradient Descent
+
Stochastic Gradient Descent
Stochastic gradient descent (SGD) and variants thereof address some of
@@ -699,7 +1068,7 @@ $$
-
Computation of gradients
+
Computation of gradients
This in turn means that the gradient can be
@@ -721,7 +1090,7 @@ minibatches. We denote these minibatches by \( B_k \) where
-
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
@@ -747,7 +1116,7 @@ $$
-
The gradient step
+
The gradient step
Thus a gradient descent step now looks like
@@ -768,7 +1137,7 @@ the number of minibatches, as exemplified in the code below.
-
Simple example code
+
Simple example code
@@ -800,7 +1169,7 @@ all \( n \) datapoints.
-
When do we stop?
+
When do we stop?
A natural question is when do we stop the search for a new minimum?
@@ -817,7 +1186,7 @@ gave the lowest value.
-
Slightly different approach
+
Slightly different approach
Another approach is to let the step length \( \gamma_j \) depend on the
@@ -868,7 +1237,7 @@ j = 0
-
Conjugate gradient (CG) method
+
Conjugate gradient (CG) method
@@ -902,7 +1271,7 @@ When we have found the exact solution, \( \hat{r}=0 \).
-
Conjugate gradient method
+
Conjugate gradient method
@@ -923,7 +1292,7 @@ If we search for a minimum of the quantum mechanical variance, then the matrix
-
Conjugate gradient method, Newton's method first
+
Conjugate gradient method, Newton's method first
@@ -951,7 +1320,7 @@ $$
-
Conjugate gradient method
+
Conjugate gradient method
@@ -984,7 +1353,7 @@ this inner product. Being conjugate is a symmetric relation: if \( \hat{s} \) is
-
Conjugate gradient method
+
Conjugate gradient method
@@ -1003,7 +1372,7 @@ which is zero unless \( i=j \).
-
Conjugate gradient method
+
Conjugate gradient method
@@ -1033,7 +1402,7 @@ $$
-
Conjugate gradient method
+
Conjugate gradient method
@@ -1070,7 +1439,7 @@ $$
-
Conjugate gradient method and iterations
+
Conjugate gradient method and iterations
@@ -1107,7 +1476,7 @@ instead.
-
Conjugate gradient method
+
Conjugate gradient method
@@ -1140,7 +1509,7 @@ hence the name conjugate gradient method.
-
[2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University
-
Sep 21, 2018
+
Sep 27, 2018
@@ -167,7 +188,209 @@ some approximative/numerical method to compute the minimum.
-
Steepest descent
+
Revisiting our Logistic Regression case
+
+
+In our discussion on Logistic Regression we defined we studied first the
+case of
+two classes, with \( y_i \) either
+\( 0 \) or \( 1 \). Furthermore we assumed also that we have only two
+parameters \( \beta \) in our fitting of the Sigmoid function, that is we
+defined probabilities
+
+$$
+\begin{align*}
+p(y_i=1|x_i,\hat{\beta}) &= \frac{\exp{(\beta_0+\beta_1x_i)}}{1+\exp{(\beta_0+\beta_1x_i)}},\nonumber\\
+p(y_i=0|x_i,\hat{\beta}) &= 1 - p(y_i=1|x_i,\hat{\beta}),
+\end{align*}
+$$
+
+where \( \hat{\beta} \) are the weights we wish to extract from data, in our case \( \beta_0 \) and \( \beta_1 \).
+
+
+
+
+
The equations to solve
+
+
+Our compact equations used a definition of a vector \( \hat{y} \) with \( n \)
+elements \( y_i \), an \( n\times p \) matrix \( \hat{X} \) which contains the
+\( x_i \) values and a vector \( \hat{p} \) of fitted probabilities
+\( p(y_i\vert x_i,\hat{\beta}) \). We rewrote in a more compact form
+the first derivative of the cost function as
+
+$$
+\frac{\partial \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}} = -\hat{X}^T\left(\hat{y}-\hat{p}\right).
+$$
+
+
+If we in addition define a diagonal matrix \( \hat{W} \) with elements
+\( p(y_i\vert x_i,\hat{\beta})(1-p(y_i\vert x_i,\hat{\beta}) \), we can obtain a compact expression of the second derivative as
+
+$$
+\frac{\partial^2 \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}\partial \hat{\beta}^T} = \hat{X}^T\hat{W}\hat{X}.
+$$
+
+This defines what we call the Hessian.
+
+
+
+
+
Solving using Newton-Raphson's method
+
+
+If we can set up these equations, Newton-Raphson's iterative method is the nomrally the method of choice. It requires however that we setting the matrices that define the first and second derivatives.
+
+
+Our iterative scheme is then given by
+
+$$
+\hat{\beta}^{\mathrm{new}} = \hat{\beta}^{\mathrm{old}}-\left(\frac{\partial^2 \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}\partial \hat{\beta}^T}\right)^{-1}\times \left(\frac{\partial \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}}\right)_{\hat{\beta}^{\mathrm{old}}},
+$$
+
+or in matrix form as
+
+$$
+\hat{\beta}^{\mathrm{new}} = \hat{\beta}^{\mathrm{old}}-\left(\hat{X}^T\hat{W}\hat{X} \right)^{-1}\times \left(-\hat{X}^T(\hat{y}-\hat{p}) \right)_{\hat{\beta}^{\mathrm{old}}}.
+$$
+
+The right-hand side is computed with the old values of \( \beta \).
+
+
+If we can compute these matrices, in particular the Hessian, the above is often the easiest method to implement.
+
+
+
+
+
Brief reminder on Newton-Raphson's method
+
+
+Let us quicly remind ourselves how we derive the above method.
+
+
+Perhaps the most celebrated of all one-dimensional root-finding
+routines is Newton's method, also called the Newton-Raphson
+method. This method is distinguished from the previously discussed
+methods by the fact that it requires the evaluation of both the
+function \( f \) and its derivative \( f' \) at arbitrary points. In this
+sense, it is taylored to cases with e.g., transcendental equations.
+If you can only calculate the derivative
+numerically and/or your function is not of the smooth type, we
+discourage the use of this method.
+
+
+
+
+
The equations
+
+
+The Newton-Raphson formula consists geometrically of extending the
+tangent line at a current point until it crosses zero, then setting
+the next guess to the abscissa of that zero-crossing. The mathematics
+behind this method is rather simple. Employing a Taylor expansion for
+\( x \) sufficiently close to the solution \( s \), we have
+
+$$
+ f(s)=0=f(x)+(s-x)f'(x)+\frac{(s-x)^2}{2}f''(x) +\dots.
+ \label{eq:taylornr}
+$$
+
+
+For small enough values of the function and for well-behaved
+functions, the terms beyond linear are unimportant, hence we obtain
+
+$$
+ f(x)+(s-x)f'(x)\approx 0,
+$$
+
+yielding
+$$
+ s\approx x-\frac{f(x)}{f'(x)}.
+$$
+
+
+Having in mind an iterative procedure, it is natural to start iterating with
+$$
+ x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}.
+$$
+
+
+
+
+
Simple geometric interpretation
+
+
+The above is Newton-Raphson's method. It has a simple geometric
+interpretation, namely \( x_{n+1} \) is the point where the tangent from
+\( (x_n,f(x_n)) \) crosses the $x-$axis. Close to the solution,
+Newton-Raphson converges fast to the desired result. However, if we
+are far from a root, where the higher-order terms in the series are
+important, the Newton-Raphson formula can give grossly inaccurate
+results. For instance, the initial guess for the root might be so far
+from the true root as to let the search interval include a local
+maximum or minimum of the function. If an iteration places a trial
+guess near such a local extremum, so that the first derivative nearly
+vanishes, then Newton-Raphson may fail totally
+
+
+
+
+
Extending to more than one variable
+
+
+Newton's method can be generalized to systems of several non-linear equations
+and variables. Consider the case with two equations
+$$
+ \begin{array}{cc} f_1(x_1,x_2) &=0\\
+ f_2(x_1,x_2) &=0\end{array},
+$$
+
+which we Taylor expand to obtain
+
+$$
+ \begin{array}{cc} 0=f_1(x_1+h_1,x_2+h_2)=&f_1(x_1,x_2)+h_1
+ \partial f_1/\partial x_1+h_2
+ \partial f_1/\partial x_2+\dots\\
+ 0=f_2(x_1+h_1,x_2+h_2)=&f_2(x_1,x_2)+h_1
+ \partial f_2/\partial x_1+h_2
+ \partial f_2/\partial x_2+\dots
+ \end{array}.
+$$
+
+Defining the Jacobian matrix \( {\bf \hat{J}} \) we have
+$$
+ {\bf \hat{J}}=\left( \begin{array}{cc}
+ \partial f_1/\partial x_1 & \partial f_1/\partial x_2 \\
+ \partial f_2/\partial x_1 &\partial f_2/\partial x_2
+ \end{array} \right),
+$$
+
+we can rephrase Newton's method as
+$$
+\left(\begin{array}{c} x_1^{n+1} \\ x_2^{n+1} \end{array} \right)=
+\left(\begin{array}{c} x_1^{n} \\ x_2^{n} \end{array} \right)+
+\left(\begin{array}{c} h_1^{n} \\ h_2^{n} \end{array} \right),
+$$
+
+where we have defined
+$$
+ \left(\begin{array}{c} h_1^{n} \\ h_2^{n} \end{array} \right)=
+ -{\bf \hat{J}}^{-1}
+ \left(\begin{array}{c} f_1(x_1^{n},x_2^{n}) \\ f_2(x_1^{n},x_2^{n}) \end{array} \right).
+$$
+
+We need thus to compute the inverse of the Jacobian matrix and it
+is to understand that difficulties may
+arise in case \( {\bf \hat{J}} \) is nearly singular.
+
+
+It is rather straightforward to extend the above scheme to systems of
+more than two non-linear equations. In our case, the Jacobian matrix is given by the Hessian that represents the second derivative of cost function.
+
+
+
+
+
Steepest descent
The method of steepest descent The basic idea of gradient descent is
@@ -191,7 +414,7 @@ we are always moving towards smaller function values, i.e a minimum.
-
More on Steepest descent
+
More on Steepest descent
The previous observation is the basis of the method of steepest
@@ -210,7 +433,7 @@ the learning rate within the context of Machine Learning.
-
The ideal
+
The ideal
Ideally the sequence \( \{\mathbf{x}_k \}_{k=0} \) converges to a global
@@ -236,7 +459,7 @@ Note that the gradient is a function of \( \mathbf{x} =
-
The sensitiveness of the gradient descent
+
The sensitiveness of the gradient descent
The gradient descent method
@@ -255,7 +478,7 @@ randomness. One such method is that of Stochastic Gradient Descent
-
Convex functions
+
Convex functions
Ideally we want our cost/loss function to be convex(concave).
@@ -275,7 +498,7 @@ regular polygons (triangles, rectangles, pentagons, etc...).
-
Convex function
+
Convex function
Convex function: Let \( X \subset \mathbb{R}^n \) be a convex set. Assume that the function \( f: X \rightarrow \mathbb{R} \) is continuous, then \( f \) is said to be convex if $$f(tx_1 + (1-t)x_2) \leq tf(x_1) + (1-t)f(x_2) $$ for all \( x_1, x_2 \in X \) and for all \( t \in [0,1] \). If \( \leq \) is replaced with a strict inequaltiy in the definition, we demand \( x_1 \neq x_2 \) and \( t\in(0,1) \) then \( f \) is said to be strictly convex. For a single variable function, convexity means that if you draw a straight line connecting \( f(x_1) \) and \( f(x_2) \), the value of the function on the interval \( [x_1,x_2] \) is always below the line as illustrated below.
@@ -283,7 +506,7 @@ regular polygons (triangles, rectangles, pentagons, etc...).
-
Conditions on convex functions
+
Conditions on convex functions
In the following we state first and second-order conditions which
@@ -325,7 +548,7 @@ This condition is particularly useful since it gives us an procedure for determi
-
More on convex functions
+
More on convex functions
The next result is of great importance to us and the reason why we are
@@ -355,7 +578,7 @@ This result means that if we know that the cost/loss function is convex and we a
-
Some simple problems
+
Some simple problems
Show that \( f(x)=x^2 \) is convex for \( x \in \mathbb{R} \) using the definition of convexity. Hint: If you re-write the definition, \( f \) is convex if the following holds for all \( x,y \in D_f \) and any \( \lambda \in [0,1] \) $\lambda f(x)+(1-\lambda)f(y)-f(\lambda x + (1-\lambda) y ) \geq 0$.
@@ -379,10 +602,147 @@ This result means that if we know that the cost/loss function is convex and we a
Using the definition of convexity, try to show that a function satisfying the properties above is convex (the third condition is not needed to show this).
+
+
+
+
Standard steepest descent
+
+
+Before we proceed, we would like to mention the approach called the standard Steepest descent, which again leads to us having to be able to compute a matrix.
+
+
+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 \).
+
+
+
+
+
Conjugate gradient method
+
+
+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.
+
+
+More details will be added here soon.
+
+
+
+
+
Simple codes for steepest descent and conjugate gradient using a \( 2\times 2 \) matrix, in c++, Python code to come
+
+
+
+
+
+
+
#include<cmath>
+#include<iostream>
+#include<fstream>
+#include<iomanip>
+#include"vectormatrixclass.h"
+usingnamespace std;
+// Main function begins here
+intmain(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;
+}
+
+
+
+
+
+
+
+
+
The routine for the steepest descent method
+
+
+
+
+
+
+
Vector SteepestDescent(Matrix A, Vector b, Vector x0){
+ int IterMax, i;
+ int dim = x0.Dimension();
+ constdouble 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;
+}
+
+
+
+
+
-
Revisiting our first homework
+
Revisiting our first homework
We will use linear regression as a case study for the gradient descent
@@ -415,7 +775,7 @@ $$
-
Gradient descent example
+
Gradient descent example
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 \)
@@ -440,7 +800,7 @@ and we want to find \( \beta \) such that \( C(\beta) \) is minimized.
-
The derivative of the cost/loss function
+
The derivative of the cost/loss function
Computing \( \partial C(\beta) / \partial \beta_0 \) and \( \partial C(\beta) / \partial \beta_1 \) we can show that the gradient can be written as
@@ -455,7 +815,7 @@ where \( X \) is the design matrix defined above.
-
The Hessian matrix
+
The Hessian matrix
The Hessian matrix of \( C(\beta) \) is given by
$$
\hat{H} \equiv \begin{bmatrix}
@@ -469,7 +829,7 @@ This result implies that \( C(\beta) \) is a convex function since the matrix \(
-
Simple program
+
Simple program
We can now write a program that minimizes \( C(\beta) \) using the gradient descent method with a constant learning rate \( \gamma \) according to
@@ -510,7 +870,7 @@ beta_NE = np.dot(Xt_X_inv,Xt_y)
-
Gradient Descent Example
+
Gradient Descent Example
Another simple example is here
@@ -559,7 +919,7 @@ plt.show()
-
And a corresponding example using scikit-learn
+
And a corresponding example using scikit-learn
@@ -583,7 +943,7 @@ sgdreg.fit(x,y.ravel())
-
Gradient descent and Ridge
+
Gradient descent and Ridge
We have also discussed Ridge regression where the loss function contains a regularized given by the \( L_2 \) norm of \( \beta \),
@@ -640,7 +1000,7 @@ beta_ridge = np.dot(Z,np.dot(X.T,y))
-
Stochastic Gradient Descent
+
Stochastic Gradient Descent
Stochastic gradient descent (SGD) and variants thereof address some of
@@ -658,7 +1018,7 @@ $$
-
Computation of gradients
+
Computation of gradients
This in turn means that the gradient can be
@@ -678,7 +1038,7 @@ minibatches. We denote these minibatches by \( B_k \) where
-
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
@@ -702,7 +1062,7 @@ $$
-
The gradient step
+
The gradient step
Thus a gradient descent step now looks like
@@ -721,7 +1081,7 @@ the number of minibatches, as exemplified in the code below.
-
Simple example code
+
Simple example code
@@ -753,7 +1113,7 @@ all \( n \) datapoints.
-
When do we stop?
+
When do we stop?
A natural question is when do we stop the search for a new minimum?
@@ -770,7 +1130,7 @@ gave the lowest value.
-
Slightly different approach
+
Slightly different approach
Another approach is to let the step length \( \gamma_j \) depend on the
@@ -818,7 +1178,7 @@ j = 0
-
Conjugate gradient (CG) method
+
Conjugate gradient (CG) method
@@ -849,7 +1209,7 @@ When we have found the exact solution, \( \hat{r}=0 \).
-
Conjugate gradient method
+
Conjugate gradient method
@@ -871,7 +1231,7 @@ If we search for a minimum of the quantum mechanical variance, then the matrix
-
Conjugate gradient method, Newton's method first
+
Conjugate gradient method, Newton's method first
@@ -896,7 +1256,7 @@ $$
-
Conjugate gradient method
+
Conjugate gradient method
@@ -926,7 +1286,7 @@ this inner product. Being conjugate is a symmetric relation: if \( \hat{s} \) is
-
Conjugate gradient method
+
Conjugate gradient method
@@ -944,7 +1304,7 @@ which is zero unless \( i=j \).
-
Conjugate gradient method
+
Conjugate gradient method
@@ -971,7 +1331,7 @@ $$
-
Conjugate gradient method
+
Conjugate gradient method
@@ -1003,7 +1363,7 @@ $$
-
Conjugate gradient method and iterations
+
Conjugate gradient method and iterations
@@ -1039,7 +1399,7 @@ instead.
-
Conjugate gradient method
+
Conjugate gradient method
@@ -1069,7 +1429,7 @@ hence the name conjugate gradient method.
[2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University
-
Sep 21, 2018
+
Sep 27, 2018
@@ -172,7 +193,209 @@ some approximative/numerical method to compute the minimum.
-
Steepest descent
+
Revisiting our Logistic Regression case
+
+
+In our discussion on Logistic Regression we defined we studied first the
+case of
+two classes, with \( y_i \) either
+\( 0 \) or \( 1 \). Furthermore we assumed also that we have only two
+parameters \( \beta \) in our fitting of the Sigmoid function, that is we
+defined probabilities
+
+$$
+\begin{align*}
+p(y_i=1|x_i,\hat{\beta}) &= \frac{\exp{(\beta_0+\beta_1x_i)}}{1+\exp{(\beta_0+\beta_1x_i)}},\nonumber\\
+p(y_i=0|x_i,\hat{\beta}) &= 1 - p(y_i=1|x_i,\hat{\beta}),
+\end{align*}
+$$
+
+where \( \hat{\beta} \) are the weights we wish to extract from data, in our case \( \beta_0 \) and \( \beta_1 \).
+
+
+
+
+
The equations to solve
+
+
+Our compact equations used a definition of a vector \( \hat{y} \) with \( n \)
+elements \( y_i \), an \( n\times p \) matrix \( \hat{X} \) which contains the
+\( x_i \) values and a vector \( \hat{p} \) of fitted probabilities
+\( p(y_i\vert x_i,\hat{\beta}) \). We rewrote in a more compact form
+the first derivative of the cost function as
+
+$$
+\frac{\partial \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}} = -\hat{X}^T\left(\hat{y}-\hat{p}\right).
+$$
+
+
+If we in addition define a diagonal matrix \( \hat{W} \) with elements
+\( p(y_i\vert x_i,\hat{\beta})(1-p(y_i\vert x_i,\hat{\beta}) \), we can obtain a compact expression of the second derivative as
+
+$$
+\frac{\partial^2 \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}\partial \hat{\beta}^T} = \hat{X}^T\hat{W}\hat{X}.
+$$
+
+This defines what we call the Hessian.
+
+
+
+
+
Solving using Newton-Raphson's method
+
+
+If we can set up these equations, Newton-Raphson's iterative method is the nomrally the method of choice. It requires however that we setting the matrices that define the first and second derivatives.
+
+
+Our iterative scheme is then given by
+
+$$
+\hat{\beta}^{\mathrm{new}} = \hat{\beta}^{\mathrm{old}}-\left(\frac{\partial^2 \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}\partial \hat{\beta}^T}\right)^{-1}\times \left(\frac{\partial \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}}\right)_{\hat{\beta}^{\mathrm{old}}},
+$$
+
+or in matrix form as
+
+$$
+\hat{\beta}^{\mathrm{new}} = \hat{\beta}^{\mathrm{old}}-\left(\hat{X}^T\hat{W}\hat{X} \right)^{-1}\times \left(-\hat{X}^T(\hat{y}-\hat{p}) \right)_{\hat{\beta}^{\mathrm{old}}}.
+$$
+
+The right-hand side is computed with the old values of \( \beta \).
+
+
+If we can compute these matrices, in particular the Hessian, the above is often the easiest method to implement.
+
+
+
+
+
Brief reminder on Newton-Raphson's method
+
+
+Let us quicly remind ourselves how we derive the above method.
+
+
+Perhaps the most celebrated of all one-dimensional root-finding
+routines is Newton's method, also called the Newton-Raphson
+method. This method is distinguished from the previously discussed
+methods by the fact that it requires the evaluation of both the
+function \( f \) and its derivative \( f' \) at arbitrary points. In this
+sense, it is taylored to cases with e.g., transcendental equations.
+If you can only calculate the derivative
+numerically and/or your function is not of the smooth type, we
+discourage the use of this method.
+
+
+
+
+
The equations
+
+
+The Newton-Raphson formula consists geometrically of extending the
+tangent line at a current point until it crosses zero, then setting
+the next guess to the abscissa of that zero-crossing. The mathematics
+behind this method is rather simple. Employing a Taylor expansion for
+\( x \) sufficiently close to the solution \( s \), we have
+
+$$
+ f(s)=0=f(x)+(s-x)f'(x)+\frac{(s-x)^2}{2}f''(x) +\dots.
+ \label{eq:taylornr}
+$$
+
+
+For small enough values of the function and for well-behaved
+functions, the terms beyond linear are unimportant, hence we obtain
+
+$$
+ f(x)+(s-x)f'(x)\approx 0,
+$$
+
+yielding
+$$
+ s\approx x-\frac{f(x)}{f'(x)}.
+$$
+
+
+Having in mind an iterative procedure, it is natural to start iterating with
+$$
+ x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}.
+$$
+
+
+
+
+
Simple geometric interpretation
+
+
+The above is Newton-Raphson's method. It has a simple geometric
+interpretation, namely \( x_{n+1} \) is the point where the tangent from
+\( (x_n,f(x_n)) \) crosses the $x-$axis. Close to the solution,
+Newton-Raphson converges fast to the desired result. However, if we
+are far from a root, where the higher-order terms in the series are
+important, the Newton-Raphson formula can give grossly inaccurate
+results. For instance, the initial guess for the root might be so far
+from the true root as to let the search interval include a local
+maximum or minimum of the function. If an iteration places a trial
+guess near such a local extremum, so that the first derivative nearly
+vanishes, then Newton-Raphson may fail totally
+
+
+
+
+
Extending to more than one variable
+
+
+Newton's method can be generalized to systems of several non-linear equations
+and variables. Consider the case with two equations
+$$
+ \begin{array}{cc} f_1(x_1,x_2) &=0\\
+ f_2(x_1,x_2) &=0\end{array},
+$$
+
+which we Taylor expand to obtain
+
+$$
+ \begin{array}{cc} 0=f_1(x_1+h_1,x_2+h_2)=&f_1(x_1,x_2)+h_1
+ \partial f_1/\partial x_1+h_2
+ \partial f_1/\partial x_2+\dots\\
+ 0=f_2(x_1+h_1,x_2+h_2)=&f_2(x_1,x_2)+h_1
+ \partial f_2/\partial x_1+h_2
+ \partial f_2/\partial x_2+\dots
+ \end{array}.
+$$
+
+Defining the Jacobian matrix \( {\bf \hat{J}} \) we have
+$$
+ {\bf \hat{J}}=\left( \begin{array}{cc}
+ \partial f_1/\partial x_1 & \partial f_1/\partial x_2 \\
+ \partial f_2/\partial x_1 &\partial f_2/\partial x_2
+ \end{array} \right),
+$$
+
+we can rephrase Newton's method as
+$$
+\left(\begin{array}{c} x_1^{n+1} \\ x_2^{n+1} \end{array} \right)=
+\left(\begin{array}{c} x_1^{n} \\ x_2^{n} \end{array} \right)+
+\left(\begin{array}{c} h_1^{n} \\ h_2^{n} \end{array} \right),
+$$
+
+where we have defined
+$$
+ \left(\begin{array}{c} h_1^{n} \\ h_2^{n} \end{array} \right)=
+ -{\bf \hat{J}}^{-1}
+ \left(\begin{array}{c} f_1(x_1^{n},x_2^{n}) \\ f_2(x_1^{n},x_2^{n}) \end{array} \right).
+$$
+
+We need thus to compute the inverse of the Jacobian matrix and it
+is to understand that difficulties may
+arise in case \( {\bf \hat{J}} \) is nearly singular.
+
+
+It is rather straightforward to extend the above scheme to systems of
+more than two non-linear equations. In our case, the Jacobian matrix is given by the Hessian that represents the second derivative of cost function.
+
+
+
+
+
Steepest descent
The method of steepest descent The basic idea of gradient descent is
@@ -196,7 +419,7 @@ we are always moving towards smaller function values, i.e a minimum.
-
More on Steepest descent
+
More on Steepest descent
The previous observation is the basis of the method of steepest
@@ -215,7 +438,7 @@ the learning rate within the context of Machine Learning.
-
The ideal
+
The ideal
Ideally the sequence \( \{\mathbf{x}_k \}_{k=0} \) converges to a global
@@ -241,7 +464,7 @@ Note that the gradient is a function of \( \mathbf{x} =
-
The sensitiveness of the gradient descent
+
The sensitiveness of the gradient descent
The gradient descent method
@@ -260,7 +483,7 @@ randomness. One such method is that of Stochastic Gradient Descent
-
Convex functions
+
Convex functions
Ideally we want our cost/loss function to be convex(concave).
@@ -280,7 +503,7 @@ regular polygons (triangles, rectangles, pentagons, etc...).
-
Convex function
+
Convex function
Convex function: Let \( X \subset \mathbb{R}^n \) be a convex set. Assume that the function \( f: X \rightarrow \mathbb{R} \) is continuous, then \( f \) is said to be convex if $$f(tx_1 + (1-t)x_2) \leq tf(x_1) + (1-t)f(x_2) $$ for all \( x_1, x_2 \in X \) and for all \( t \in [0,1] \). If \( \leq \) is replaced with a strict inequaltiy in the definition, we demand \( x_1 \neq x_2 \) and \( t\in(0,1) \) then \( f \) is said to be strictly convex. For a single variable function, convexity means that if you draw a straight line connecting \( f(x_1) \) and \( f(x_2) \), the value of the function on the interval \( [x_1,x_2] \) is always below the line as illustrated below.
@@ -288,7 +511,7 @@ regular polygons (triangles, rectangles, pentagons, etc...).
-
Conditions on convex functions
+
Conditions on convex functions
In the following we state first and second-order conditions which
@@ -330,7 +553,7 @@ This condition is particularly useful since it gives us an procedure for determi
-
More on convex functions
+
More on convex functions
The next result is of great importance to us and the reason why we are
@@ -360,7 +583,7 @@ This result means that if we know that the cost/loss function is convex and we a
-
Some simple problems
+
Some simple problems
Show that \( f(x)=x^2 \) is convex for \( x \in \mathbb{R} \) using the definition of convexity. Hint: If you re-write the definition, \( f \) is convex if the following holds for all \( x,y \in D_f \) and any \( \lambda \in [0,1] \) $\lambda f(x)+(1-\lambda)f(y)-f(\lambda x + (1-\lambda) y ) \geq 0$.
@@ -384,10 +607,147 @@ This result means that if we know that the cost/loss function is convex and we a
Using the definition of convexity, try to show that a function satisfying the properties above is convex (the third condition is not needed to show this).
+
+
+
+
Standard steepest descent
+
+
+Before we proceed, we would like to mention the approach called the standard Steepest descent, which again leads to us having to be able to compute a matrix.
+
+
+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 \).
+
+
+
+
+
Conjugate gradient method
+
+
+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.
+
+
+More details will be added here soon.
+
+
+
+
+
Simple codes for steepest descent and conjugate gradient using a \( 2\times 2 \) matrix, in c++, Python code to come
+
+
+
+
+
+
+
#include<cmath>
+#include<iostream>
+#include<fstream>
+#include<iomanip>
+#include"vectormatrixclass.h"
+usingnamespace std;
+// Main function begins here
+intmain(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;
+}
+
+
+
+
+
+
+
+
+
The routine for the steepest descent method
+
+
+
+
+
+
+
Vector SteepestDescent(Matrix A, Vector b, Vector x0){
+ int IterMax, i;
+ int dim = x0.Dimension();
+ constdouble 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;
+}
+
+
+
+
+
-
Revisiting our first homework
+
Revisiting our first homework
We will use linear regression as a case study for the gradient descent
@@ -420,7 +780,7 @@ $$
-
Gradient descent example
+
Gradient descent example
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 \)
@@ -445,7 +805,7 @@ and we want to find \( \beta \) such that \( C(\beta) \) is minimized.
-
The derivative of the cost/loss function
+
The derivative of the cost/loss function
Computing \( \partial C(\beta) / \partial \beta_0 \) and \( \partial C(\beta) / \partial \beta_1 \) we can show that the gradient can be written as
@@ -460,7 +820,7 @@ where \( X \) is the design matrix defined above.
-
The Hessian matrix
+
The Hessian matrix
The Hessian matrix of \( C(\beta) \) is given by
$$
\hat{H} \equiv \begin{bmatrix}
@@ -474,7 +834,7 @@ This result implies that \( C(\beta) \) is a convex function since the matrix \(
-
Simple program
+
Simple program
We can now write a program that minimizes \( C(\beta) \) using the gradient descent method with a constant learning rate \( \gamma \) according to
@@ -515,7 +875,7 @@ beta_NE = np.
-
Gradient Descent Example
+
Gradient Descent Example
Another simple example is here
@@ -564,7 +924,7 @@ plt.show()
-
And a corresponding example using scikit-learn
+
And a corresponding example using scikit-learn
@@ -588,7 +948,7 @@ sgdreg.fit(x,y.
-
Gradient descent and Ridge
+
Gradient descent and Ridge
We have also discussed Ridge regression where the loss function contains a regularized given by the \( L_2 \) norm of \( \beta \),
@@ -645,7 +1005,7 @@ beta_ridge = np
-
Stochastic Gradient Descent
+
Stochastic Gradient Descent
Stochastic gradient descent (SGD) and variants thereof address some of
@@ -663,7 +1023,7 @@ $$
-
Computation of gradients
+
Computation of gradients
This in turn means that the gradient can be
@@ -683,7 +1043,7 @@ minibatches. We denote these minibatches by \( B_k \) where
-
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
@@ -707,7 +1067,7 @@ $$
-
The gradient step
+
The gradient step
Thus a gradient descent step now looks like
@@ -726,7 +1086,7 @@ the number of minibatches, as exemplified in the code below.
-
Simple example code
+
Simple example code
@@ -758,7 +1118,7 @@ all \( n \) datapoints.
-
When do we stop?
+
When do we stop?
A natural question is when do we stop the search for a new minimum?
@@ -775,7 +1135,7 @@ gave the lowest value.
-
Slightly different approach
+
Slightly different approach
Another approach is to let the step length \( \gamma_j \) depend on the
@@ -823,7 +1183,7 @@ j =0
-
Conjugate gradient (CG) method
+
Conjugate gradient (CG) method
@@ -854,7 +1214,7 @@ When we have found the exact solution, \( \hat{r}=0 \).
-
Conjugate gradient method
+
Conjugate gradient method
@@ -876,7 +1236,7 @@ If we search for a minimum of the quantum mechanical variance, then the matrix
-
Conjugate gradient method, Newton's method first
+
Conjugate gradient method, Newton's method first
@@ -901,7 +1261,7 @@ $$
-
Conjugate gradient method
+
Conjugate gradient method
@@ -931,7 +1291,7 @@ this inner product. Being conjugate is a symmetric relation: if \( \hat{s} \) is
-
Conjugate gradient method
+
Conjugate gradient method
@@ -949,7 +1309,7 @@ which is zero unless \( i=j \).
-
Conjugate gradient method
+
Conjugate gradient method
@@ -976,7 +1336,7 @@ $$
-
Conjugate gradient method
+
Conjugate gradient method
@@ -1008,7 +1368,7 @@ $$
-
Conjugate gradient method and iterations
+
Conjugate gradient method and iterations
@@ -1044,7 +1404,7 @@ instead.
-
Conjugate gradient method
+
Conjugate gradient method
@@ -1074,7 +1434,7 @@ hence the name conjugate gradient method.
-
Conjugate gradient method
+
Conjugate gradient method
@@ -1103,7 +1463,7 @@ $$
-
Conjugate gradient method
+
Conjugate gradient method
diff --git a/doc/pub/Splines/ipynb/Splines.ipynb b/doc/pub/Splines/ipynb/Splines.ipynb
index 688508aaa..a46ac7740 100644
--- a/doc/pub/Splines/ipynb/Splines.ipynb
+++ b/doc/pub/Splines/ipynb/Splines.ipynb
@@ -10,7 +10,7 @@
" \n",
"**Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University\n",
"\n",
- "Date: **Sep 21, 2018**\n",
+ "Date: **Sep 27, 2018**\n",
"\n",
"Copyright 1999-2018, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license\n",
"\n",
@@ -28,6 +28,327 @@
"analytically, however this is not possible in general and we must use\n",
"some approximative/numerical method to compute the minimum.\n",
"\n",
+ "\n",
+ "## Revisiting our Logistic Regression case\n",
+ "\n",
+ "In our discussion on Logistic Regression we defined we studied first the \n",
+ "case of\n",
+ "two classes, with $y_i$ either\n",
+ "$0$ or $1$. Furthermore we assumed also that we have only two\n",
+ "parameters $\\beta$ in our fitting of the Sigmoid function, that is we\n",
+ "defined probabilities"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\begin{align*}\n",
+ "p(y_i=1|x_i,\\hat{\\beta}) &= \\frac{\\exp{(\\beta_0+\\beta_1x_i)}}{1+\\exp{(\\beta_0+\\beta_1x_i)}},\\nonumber\\\\\n",
+ "p(y_i=0|x_i,\\hat{\\beta}) &= 1 - p(y_i=1|x_i,\\hat{\\beta}),\n",
+ "\\end{align*}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where $\\hat{\\beta}$ are the weights we wish to extract from data, in our case $\\beta_0$ and $\\beta_1$. \n",
+ "\n",
+ "## The equations to solve\n",
+ "\n",
+ "Our compact equations used a definition of a vector $\\hat{y}$ with $n$\n",
+ "elements $y_i$, an $n\\times p$ matrix $\\hat{X}$ which contains the\n",
+ "$x_i$ values and a vector $\\hat{p}$ of fitted probabilities\n",
+ "$p(y_i\\vert x_i,\\hat{\\beta})$. We rewrote in a more compact form\n",
+ "the first derivative of the cost function as"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\frac{\\partial \\mathcal{C}(\\hat{\\beta})}{\\partial \\hat{\\beta}} = -\\hat{X}^T\\left(\\hat{y}-\\hat{p}\\right).\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "If we in addition define a diagonal matrix $\\hat{W}$ with elements \n",
+ "$p(y_i\\vert x_i,\\hat{\\beta})(1-p(y_i\\vert x_i,\\hat{\\beta})$, we can obtain a compact expression of the second derivative as"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\frac{\\partial^2 \\mathcal{C}(\\hat{\\beta})}{\\partial \\hat{\\beta}\\partial \\hat{\\beta}^T} = \\hat{X}^T\\hat{W}\\hat{X}.\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This defines what we call the Hessian.\n",
+ "\n",
+ "## Solving using Newton-Raphson's method\n",
+ "\n",
+ "If we can set up these equations, Newton-Raphson's iterative method is the nomrally the method of choice. It requires however that we setting the matrices that define the first and second derivatives. \n",
+ "\n",
+ "Our iterative scheme is then given by"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\hat{\\beta}^{\\mathrm{new}} = \\hat{\\beta}^{\\mathrm{old}}-\\left(\\frac{\\partial^2 \\mathcal{C}(\\hat{\\beta})}{\\partial \\hat{\\beta}\\partial \\hat{\\beta}^T}\\right)^{-1}\\times \\left(\\frac{\\partial \\mathcal{C}(\\hat{\\beta})}{\\partial \\hat{\\beta}}\\right)_{\\hat{\\beta}^{\\mathrm{old}}},\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "or in matrix form as"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\hat{\\beta}^{\\mathrm{new}} = \\hat{\\beta}^{\\mathrm{old}}-\\left(\\hat{X}^T\\hat{W}\\hat{X} \\right)^{-1}\\times \\left(-\\hat{X}^T(\\hat{y}-\\hat{p}) \\right)_{\\hat{\\beta}^{\\mathrm{old}}}.\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The right-hand side is computed with the old values of $\\beta$. \n",
+ "\n",
+ "If we can compute these matrices, in particular the Hessian, the above is often the easiest method to implement. \n",
+ "\n",
+ "\n",
+ "## Brief reminder on Newton-Raphson's method\n",
+ "\n",
+ "Let us quicly remind ourselves how we derive the above method.\n",
+ "\n",
+ "Perhaps the most celebrated of all one-dimensional root-finding\n",
+ "routines is Newton's method, also called the Newton-Raphson\n",
+ "method. This method is distinguished from the previously discussed\n",
+ "methods by the fact that it requires the evaluation of both the\n",
+ "function $f$ and its derivative $f'$ at arbitrary points. In this\n",
+ "sense, it is taylored to cases with e.g., transcendental equations.\n",
+ "If you can only calculate the derivative\n",
+ "numerically and/or your function is not of the smooth type, we\n",
+ "discourage the use of this method.\n",
+ "\n",
+ "## The equations\n",
+ "\n",
+ "The Newton-Raphson formula consists geometrically of extending the\n",
+ "tangent line at a current point until it crosses zero, then setting\n",
+ "the next guess to the abscissa of that zero-crossing. The mathematics\n",
+ "behind this method is rather simple. Employing a Taylor expansion for\n",
+ "$x$ sufficiently close to the solution $s$, we have"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "
\n",
+ "\n",
+ "$$\n",
+ "f(s)=0=f(x)+(s-x)f'(x)+\\frac{(s-x)^2}{2}f''(x) +\\dots.\n",
+ " \\label{eq:taylornr} \\tag{1}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "For small enough values of the function and for well-behaved\n",
+ "functions, the terms beyond linear are unimportant, hence we obtain"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "f(x)+(s-x)f'(x)\\approx 0,\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "yielding"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "s\\approx x-\\frac{f(x)}{f'(x)}.\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Having in mind an iterative procedure, it is natural to start iterating with"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "x_{n+1}=x_n-\\frac{f(x_n)}{f'(x_n)}.\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Simple geometric interpretation\n",
+ "\n",
+ "The above is Newton-Raphson's method. It has a simple geometric\n",
+ "interpretation, namely $x_{n+1}$ is the point where the tangent from\n",
+ "$(x_n,f(x_n))$ crosses the $x-$axis. Close to the solution,\n",
+ "Newton-Raphson converges fast to the desired result. However, if we\n",
+ "are far from a root, where the higher-order terms in the series are\n",
+ "important, the Newton-Raphson formula can give grossly inaccurate\n",
+ "results. For instance, the initial guess for the root might be so far\n",
+ "from the true root as to let the search interval include a local\n",
+ "maximum or minimum of the function. If an iteration places a trial\n",
+ "guess near such a local extremum, so that the first derivative nearly\n",
+ "vanishes, then Newton-Raphson may fail totally\n",
+ "\n",
+ "\n",
+ "## Extending to more than one variable\n",
+ "\n",
+ "Newton's method can be generalized to systems of several non-linear equations\n",
+ "and variables. Consider the case with two equations"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\begin{array}{cc} f_1(x_1,x_2) &=0\\\\\n",
+ " f_2(x_1,x_2) &=0\\end{array},\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "which we Taylor expand to obtain"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\begin{array}{cc} 0=f_1(x_1+h_1,x_2+h_2)=&f_1(x_1,x_2)+h_1\n",
+ " \\partial f_1/\\partial x_1+h_2\n",
+ " \\partial f_1/\\partial x_2+\\dots\\\\\n",
+ " 0=f_2(x_1+h_1,x_2+h_2)=&f_2(x_1,x_2)+h_1\n",
+ " \\partial f_2/\\partial x_1+h_2\n",
+ " \\partial f_2/\\partial x_2+\\dots\n",
+ " \\end{array}.\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Defining the Jacobian matrix ${\\bf \\hat{J}}$ we have"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "{\\bf \\hat{J}}=\\left( \\begin{array}{cc}\n",
+ " \\partial f_1/\\partial x_1 & \\partial f_1/\\partial x_2 \\\\\n",
+ " \\partial f_2/\\partial x_1 &\\partial f_2/\\partial x_2\n",
+ " \\end{array} \\right),\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "we can rephrase Newton's method as"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\left(\\begin{array}{c} x_1^{n+1} \\\\ x_2^{n+1} \\end{array} \\right)=\n",
+ "\\left(\\begin{array}{c} x_1^{n} \\\\ x_2^{n} \\end{array} \\right)+\n",
+ "\\left(\\begin{array}{c} h_1^{n} \\\\ h_2^{n} \\end{array} \\right),\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where we have defined"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\left(\\begin{array}{c} h_1^{n} \\\\ h_2^{n} \\end{array} \\right)=\n",
+ " -{\\bf \\hat{J}}^{-1}\n",
+ " \\left(\\begin{array}{c} f_1(x_1^{n},x_2^{n}) \\\\ f_2(x_1^{n},x_2^{n}) \\end{array} \\right).\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We need thus to compute the inverse of the Jacobian matrix and it\n",
+ "is to understand that difficulties may\n",
+ "arise in case ${\\bf \\hat{J}}$ is nearly singular.\n",
+ "\n",
+ "It is rather straightforward to extend the above scheme to systems of\n",
+ "more than two non-linear equations. In our case, the Jacobian matrix is given by the Hessian that represents the second derivative of cost function. \n",
+ "\n",
+ "\n",
+ "\n",
"## Steepest descent\n",
"\n",
"The method of steepest descent The basic idea of gradient descent is\n",
@@ -218,6 +539,154 @@
"\n",
"Using the definition of convexity, try to show that a function satisfying the properties above is convex (the third condition is not needed to show this).\n",
"\n",
+ "\n",
+ "## Standard steepest descent\n",
+ "\n",
+ "\n",
+ "Before we proceed, we would like to mention the approach called the **standard Steepest descent**, which again leads to us having to be able to compute a matrix.\n",
+ "\n",
+ "[The success of the CG method](https://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf)\n",
+ "for finding solutions of non-linear problems is based on the theory\n",
+ "of conjugate gradients for linear systems of equations. It belongs to\n",
+ "the class of iterative methods for solving problems from linear\n",
+ "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",
+ "## 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\n",
+ "symmetric. If we search for a minimum of the quantum mechanical\n",
+ "variance, then the matrix $\\hat{A}$, which is called the Hessian, is\n",
+ "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",
+ "## Simple codes for steepest descent and conjugate gradient using a $2\\times 2$ matrix, in c++, Python code to come"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ " #include \n",
+ " #include \n",
+ " #include \n",
+ " #include \n",
+ " #include \"vectormatrixclass.h\"\n",
+ " using namespace std;\n",
+ " // Main function begins here\n",
+ " int main(int argc, char * argv[]){\n",
+ " int dim = 2;\n",
+ " Vector x(dim),xsd(dim), b(dim),x0(dim);\n",
+ " Matrix A(dim,dim);\n",
+ " \n",
+ " // Set our initial guess\n",
+ " x0(0) = x0(1) = 0;\n",
+ " // Set the matrix\n",
+ " A(0,0) = 3; A(1,0) = 2; A(0,1) = 2; A(1,1) = 6;\n",
+ " b(0) = 2; b(1) = -8;\n",
+ " 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",
+ " }\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## The routine for the steepest descent method"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ " Vector SteepestDescent(Matrix A, Vector b, Vector x0){\n",
+ " int IterMax, i;\n",
+ " int dim = x0.Dimension();\n",
+ " const double tolerance = 1.0e-14;\n",
+ " Vector x(dim),f(dim),z(dim);\n",
+ " double c,alpha,d;\n",
+ " IterMax = 30;\n",
+ " x = x0;\n",
+ " f = A*x-b;\n",
+ " i = 0;\n",
+ " while (i <= IterMax){\n",
+ " z = A*f;\n",
+ " c = dot(f,f);\n",
+ " alpha = c/dot(f,z);\n",
+ " x = x - alpha*f;\n",
+ " f = A*x-b;\n",
+ " if(sqrt(dot(f,f)) < tolerance) break;\n",
+ " i++;\n",
+ " }\n",
+ " return x;\n",
+ " }\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
"\n",
"## Revisiting our first homework\n",
"\n",
@@ -397,7 +866,9 @@
{
"cell_type": "code",
"execution_count": 1,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"import numpy as np\n",
@@ -431,30 +902,11 @@
},
{
"cell_type": "code",
- "execution_count": 6,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "[[4.11631855]\n",
- " [2.78555876]]\n",
- "[[4.11631855]\n",
- " [2.78555876]]\n"
- ]
- },
- {
- "data": {
- "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEWCAYAAAB1xKBvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VPW9//HXJ4GwibKjogRRwLoruEtKW1ut2tpW2yvFtfZy1duqrXax/Fq30v1W7e1tLW1dk7rUtmq9tlerwaAIGlRENkEhLCIEEERAIMnn98c5gSHMTE6SmTkzk/fz8ZhHZs76PWcm53O+6zF3R0REpDUlcSdAREQKgwKGiIhEooAhIiKRKGCIiEgkChgiIhKJAoaIiESigCE5Z2ZLzez08P33zOwPMaVjnJmtiGPfnYWZ3WNmP4w7HZIZChiyGzO7wMxmmtlmM1sTvr/KzCwb+3P3H7n7Vzu6HTMbZmZuZl0yka646UIr+UgBQ3Yys+uAO4CfA/sCg4ErgFOBshTrlOYsgSISKwUMAcDM9gFuAa5y90fcfZMHXnX3Ce6+LVzuHjP7rZk9aWabgY+Z2dlm9qqZvW9my83sphbbvsjM6sxsnZlNajHvJjOrTPh8kplNN7MNZjbbzMYlzJtqZrea2QtmtsnMnjKzAeHsmvDvBjP7wMxOTnKMPcL0v2dm84DjW8zf38z+Ymb1ZrbEzK5OmHeCmdWGx7jazH6ZMO+0hDQvN7NLw+ndzOwXZrYsXOdOM+sRzhtnZivM7LowJ7fKzC4L500EJgDfDo/l7ym+s0PN7GkzW29mC83sS+H0MjN7zcy+Hn4uDc/ZDxKO5cUwvavM7NdmVpawXQ9zlYvC83yrmR0cHuP7ZvZw8/IJx/E9M1sbFjdOSJbecPlzwrRtCLd3VKplJQ+5u156AZwJNABdWlnuHmAjQa6jBOgOjAOODD8fBawGPhcufxjwAVABdAN+Ge7n9HD+TUBl+H4IsA44K9zWJ8PPA8P5U4G3gJFAj/DzT8J5wwBPl37gJ8A0oB9wIPAGsCKcVwLMAn5AkJsaDrwNnBHOfxG4KHy/F3BS+L4c2ASMB7oC/YFjwnm3AY+H++sN/B34cThvXHgebgnXOwvYAvRNOM8/THMsvYDlwGVAF+BYYC1wWDj/COA94CPAJGAGUBrOGw2cFK43DJgPXJuwbQceA/YGDge2Ac+E52QfYB5wSYvj+GX4/X4U2AyMankcYRrXACcCpcAlwFKgW9y/f72ivZTDkGYDgLXu3tA8IeGueauZVSQs+5i7v+DuTe7+obtPdfc54efXgQcILhwA5wNPuHuNB7mU7wNNKdJwIfCkuz8ZbutpoJbgYtrsbnd/0923Ag8Dx7ThGL8ETHb39e6+HPhVwrzjCQLTLe6+3d3fBn4PXBDO3wEcYmYD3P0Dd58RTv8y8C93f8Ddd7j7Ond/LazzmQh8I9zfJuBHCdtr3uYt4XpPEgTWURGP5Rxgqbvf7e4N7v4q8BfgiwDu/gbwQ+BR4HqCYNcYzpvl7jPC9ZYCv2PX99XsZ+7+vrvPJQisT7n72+6+EfgHwcU/0ffdfZu7Pwf8b3iuW5oI/M7dZ7p7o7vfSxCMTop4zBIzBQxptg4YkFhp7O6nuHufcF7ib2V54opmdqKZVYdFORsJ6j2ai4r2T1ze3TeH20umHPhiGKQ2mNkG4DRgv4Rl3k14v4Xgbj+q3dIC1LXY9/4t9v09gnocgMsJcjYLzOxlMzsnnH4gQa6npYFAT2BWwvb+GU5vti4xQLfxeMqBE1ukdwJB3VOze8PlnnT3Rc0TzWykmT1hZu+a2fsEgWwAu1ud8H5rks+J6Xwv/F6b1RGc62Rpvq5Fmg9MsazkIQUMafYiwd3euRGWbTnE8Z8Iil4OdPd9gDuB5lZVqwguCgCYWU+CYptklgP3u3ufhFcvd/9JO9KUzG5pAYa22PeSFvvu7e5nAbj7IncfDwwCfgo8YmbNxUIHJ9nXWoIL6+EJ29vH3aMGhNaOZznwXIv07uXuVyYs8xvgCeAMMzstYfpvgQXACHffmyAwdqQVXN/wXDQbCryTIs2TW6S5p7s/0IF9Sw4pYAgA7r4BuBn4jZmdb2a9zazEzI4hKC9Ppzew3t0/NLMTCIppmj0CnBNWDJcRlNmn+t1VAp8xszPCitruYaXqAREOoZ6gqGt4mmUeBm4ws77hNr+eMO8lYJOZfSesHC81syPM7HgAM7vQzAa6exOwIVynCagCTjezL5lZFzPrb2bHhMv9HrjNzAaF2xhiZmdEOBYI7ujTHcsTwEgLGhR0DV/Hm9lHwn1dRFBXcSlwNXCvmTUHq97A+8AHZnYocOWem2+zm8PK9rEExWV/TrLM74ErwhypmVkvCxpM9M7A/iUHFDBkJ3f/GfBN4NsEF6zVBOXb3wGmp1n1KuAWM9tEUGn8cMI25wL/SZALWUVQEZu0s1xYr3AuwR1vPcEd6beI8Dt19y3AZOCFsLgjWbn4zQTFJUuAp4D7E9ZvJLjQHRPOXwv8gaCSF4JGAXPN7AOCpscXuPtWd19GUMdyHbAeeA04OlznO8BiYEZY9PMvotdR/BE4LDyWR5Mc7ybgUwR1Iu8QFNX9FOhmZkOB24GLw/qWPxHUBd0Wrn49QVDfRHARfyhimlJ5l+B7fYcggF7h7guSpLkW+Hfg1+HyiwkCmhQIc9cDlESkfSxo9lzp7lFygVLglMMQEZFIFDBERCQSFUmJiEgkymGIiEgkBTWy54ABA3zYsGFxJ0NEpKDMmjVrrbsPbH3J9AoqYAwbNoza2tq4kyEiUlDMrK71pVqnIikREYlEAUNERCJRwBARkUgUMEREJBIFDBERiUQBQ0REIlHAEBGRSBQwREQkEgUMERGJRAFDREQiyXrAMLO7zGyNmb2RZN51ZuZm1vIB9CIikmdykcO4h+DxlrsxswMJHjG5LAdpEBGRDsp6wHD3GoJnHbd0G8Gzo/VADhGRAhBLHYaZnQusdPfZEZadaGa1ZlZbX1+fg9SJiEgyOQ8YZtYT+B7wgyjLu/sUdx/j7mMGDuzwcO4iItJOceQwDgYOAmab2VLgAOAVM9s3hrSIiEhEOX+AkrvPAQY1fw6Dxhh3X5vrtIiISHS5aFb7APAiMMrMVpjZ5dnep4iIZF7WcxjuPr6V+cOynQYREek49fQWEZFIFDBERCQSBQwREYlEAUNERCJRwBARkUgUMEREJBIFDBERiUQBQ0REIlHAEBGRSBQwREQkEgUMERGJRAFDREQiUcAQEZFIFDBERCQSBQwREYlEAUNERCJRwBARkUgUMEREJBIFDBERiUQBQ0REIsl6wDCzu8xsjZm9kTDt52a2wMxeN7O/mVmfbKdDREQ6Jhc5jHuAM1tMexo4wt2PAt4EbshBOkREpAOyHjDcvQZY32LaU+7eEH6cARyQ7XSIiEjH5EMdxleAf6SaaWYTzazWzGrr6+tzmCwREUkUa8Aws0lAA1CVahl3n+LuY9x9zMCBA3OXOBER2U2XuHZsZpcC5wCfcHePKx0iIhJNLAHDzM4Evg181N23xJEGERFpm1w0q30AeBEYZWYrzOxy4NdAb+BpM3vNzO7MdjpERKRjsp7DcPfxSSb/Mdv7FRGRzMqHVlIiIlIAFDBERCQSBQwREYlEAUNERCJRwBARybSqKhg2DEpKgr9VKfsmF5TYOu6JiBSlqiqYOBG2hF3M6uqCzwATJsSXrgxQDkNEJJMmTdoVLJpt2RJML3AKGCIimbRsWdumFxAFDBGRTBo6tG3TC4gChohIJk2eDD177j6tZ89geoFTwBARyaQJE2DKFCgvB7Pg75QpBV/hDWolJSKSeRMmFEWAaEk5DBGJV5H2WShGymGISHyKuM9CMVIOQ0TiU8R9FoqRAoaIxKeI+yzsoQiK3hQwRCQ+hdRnoSMX/Oait7o6cN9V9FZgQUMBQ0TiUyh9Fjp6wS+SojcFDBGJT6H0WejoBb9Iit7USkpE4lUIfRY6esEfOjTIlSSbXkCynsMws7vMbI2ZvZEwrZ+ZPW1mi8K/fbOdDhGRdutoXUuhFL21IhdFUvcAZ7aY9l3gGXcfATwTfhYRyU8dveDnougtB62wzN0zvtE9dmI2DHjC3Y8IPy8Exrn7KjPbD5jq7qNa286YMWO8trY2q2kVEUmqqiqos1i2LMhZTJ6cP0VpLTtAQhDQwqBkZrPcfUxHdxNXpfdgd18Vvn8XGJxqQTObaGa1ZlZbX1+fm9SJSOfT2h36hAmwdCk0NQV/8yVYQM5aYcXeSsqDLE7KbI67T3H3Me4+ZuDAgTlMmUgBKYJOYbEq9H4SOWqFFVfAWB0WRRH+XRNTOkQKX6Ff7DKpvYEz7n4SHQ34OeoAGVfAeBy4JHx/CfBYTOkQKXxxX+zyRUcCZ5z9JDIR8FNUyq+//Fs8ePX0jCU165XeZvYAMA4YAKwGbgQeBR4GhgJ1wJfcfX1r21Klt0gSJSXBhaYls6C8vbMYNix5X4fy8qDOIVvrdlSG9u2VVTRc9x26rHmHdV0GMbnxBm73a8K5man0znrHPXcfn2LWJ7K9b5FOoUg6hXVYR3IJkycnb2WUi34S7Ux3U0MT8/7+FjUPvUPN9K5MW/lx3mlaAUD/xnWM3XcRvzxhKmO/MIjjL0m7qcjU01uk0MV5scsnHQmczS2e4mg2GzHdDR828OpDb1LzyBpqanvw/OoRrPcRwAiGlKziowe+TcUpi6i4YH8OPWs4JV1O2rVyhgIG7l4wr9GjR7uIJFFZ6V5e7m4W/K2sjDtFHdfWY6qsdO/Z0z0ooAtePXvm/7lIke5tU+725371mt/6iWr/ZL9a78WmnbMP6brEvzKixu/56jR/+7ll3tTYlHYXQK1n4BocexBoy0sBQ6STaO/Fv1ADZ2WlNw450JswX99tsN/U/Ud+Eff4Esq9EfOV7O93HzDJH7r2BX/n1XfbvPlMBYyc9PTOFFV6i3QScVZC58jahet4/p7F1Dy1lZoFg3h1yyiaKKWUBm4o+y++33AjZU3bdq2Q0HM7pRS90TPV01sBQ0TyTxG2/Frx8iqm3beEmmcbqHlrf+ZtOwSA7mzlxH0WUnHUBirO2ZuTLh7JXicd0faAmWZ4ELvwQgUMKUD5PB6P5M/3U+A5DG9y3qpeRk3lMmqmGTV1Q1nSEFRi9+Z9Th2wkIrjNlPxuX6MmTCKbnt3230D7QmYac6Z1dUVRrNakZ1a3gE1d1ACBY18kE/fT4G1/GpqaGLuY4upeWgVNTO6UrPyYN5tKgfKGWBrGbvvYq4+cQkV5w/iqPNG0KX78ek32J4WX7nofJiJipBcvVTpXeDKy3evxGx+lZfHnTJxz7/vJ48rsLdv3u4z73rDf352tX9m8Azva+t3nq4DSlf6l8uf9zu//JzP+/vi9C2YUh1jeyr903x/qJWUFByz5D9os7hTJu76ftLYsm6LT739Vb/l49V+eosmriO7vu2Xj6zxe/99mi+ZtrzVJq47tRYUMtisWAFDCk++3cHK7nL9/eRxDmLj8o3+j1tf9htOrvZTe8/2Mj4MYieNflT3Bf61I6f6w9+Y7qtmr27/TrJxvlOc00wFDFV6S+608pAXiVkuv588+y3Uz18bNHF9ehs1Cwbx2taRNFFKF3YwZq8FjP3IOirO7Mmpl42k70F9MrPTHLYEU7NaKUz50gpHksvV9xNzK6jlM98JmrhWN1Lz1hDmbz8YCJq4ntxnAWOP2hg0cb1kFL0G9cpOInJ4DhQwRKRw5fDu2pucRU8vZdoDK4ImrsvKWdpwIAB7s5FTB74ZNHH9fH9Gjx+5ZxPXbMlhLitTAUPNakUk9/r1g3Xr9pyegRF2mxqaeOPR5iauZdS8czCrmw4CDmKg1VOx/2KuPeFtKr44mKPOG0FpWZomrtnMccU54GE7KWCISG5VVcH77+85vaxs934WVVVwzTW7Akv//nDHHXtcUHds2cErDyyk5q9rqantyfP1I9ngI4GRHFi6kk+Wv8nYUxZQ8eUDGHXmQVhJxEc956JfyoQJeR0gWlKRlIjkVqqy+/79Ye3a4H1VFVx2GezYsfsyZWVsu+NOpm85hmlPbKTmtb158b1RbCGoZxhV9jZjD1pBxbgSKi4aRvmpB2Q+nYXQ27xFzmhgXd2SevfhHd2sAoZIoSn0hgNR6i9SXayBOoYyjDqMJo7u8SZjR6ym4vQyxl52CIOPiJh7yFQ681GSupHR0DTLvbSjm1aRlEghyafhO9qrlWEv1sytZ2DdMizV6izjiRtf5tTLRtKn/FDg0FjSmbeSPOPdoCQTm87IRkQkR5JcDNiyJSjrHzYsuCseNiwILPlq8uSgNVCChi7d+MP2i/hIt7cYfMRA6kh9UTYzzh7xJn3K98l5OvN5PKudMjl2VAsKGCL5qKoqeQBIdTFYty64G3bflevIw6DhTc7Cfifz7HHXsbbLYJowllLOxQ1/5PpV13Nwn3X85MypbP/KFXjXrik24kHgbJbqXHXUhAlBE9fy8qAYqry8MDqZZjMHlInu4u19Ad8A5gJvAA8A3dMtr6FBpFNIN8ZQquEkkr1KS2MfdqNhW4O/+uAC/9V5U/38IdN9kK3ZmbxBtsbPHzLd7/jCVH/1wQXesK1h95UrK9MfX/MyhfhY1mxKck6Og0bPxDU7Extp145hCLAE6BF+fhi4NN06ChjSKaQbYyjZBTLKK0cX0W2btvn0373uP/10tZ89aKbvw4ZdyS9d7hcNn+ZTLnrOFzz5VrRB+lKdC7P0AbR5PKY8Hq8qq1oc9wB424sgYCwH+hFUvj8BfCrdOgoY0m6FdOFobdTYlsdSUhItaGRhEMHN9Zv9mV+84jeNq/aP953lPdi8c3eHli32iYc+5/df8bwvfX55+3ZQWZn6fDSfg3QBRbkPd/fiGHzQzK4BJgNbgafcfY/CQTObCEwEGDp06Oi6FE3tRFLKs4HuWtXW9v+Wqj1RkuU62Bx0Q91GXrj7Tab9czM1c/tT+8EodlCG0cQxPRZSMWoNY08vY+ylhzDo8Aw1cU11fGapWzKVlwd/C7UfRYZlamiQOHMYfYFngYFAV+BR4MJ06yiHUaDivrsvtGHV23pnHLVeo+XxRvhe3p2zxv/8zel+9dFT/Zge891odHDvyjY/pfds/+5J1f6/N73k7y3dkOGTkKB//9THk+5c6fkeO1Hoz8MAvgj8MeHzxcBv0q2jgBGK+wLcFvlQLFCIF462fMdR6jVanvMU38uaH/zK7/uPaf7vhz7no8re2jmrB5v9E31n+c0fq/Zn/+sV31y/OdtnYFc6y8r2PJ6uXVt/0FCh3ShkUTEEjBMJWkj1BAy4F/h6unUUMDw/LsBtkQ//tPmQhmyrrEx9J96//56/jxTnZAnlDu77sMHPGTTTf3ZWtb/4+zm+bdO29PvO1g1Mqu+uf//W1y20/5UsKviAERwDNwMLCJrV3g90S7e8AoYX3sUvH+7u23LhKKTcWzJp0t+wrcFf+dN8v/3zU72R5N9LE+avPZSkiWu6/WXzopzq99P8G4qS+yrk7zNDiiJgtPWlgOH5cQFui3wJcFEuHNm4+LXnucwZusBt27TNX7jzdf/xGdV+1sCXfO+EJq7LGZKZ7yXb32+U+plOmmtoi5wFDOBp4OhM7KyjLwUMz58LcFSpyteTFZPkOl0tL8yZPrdtDUDtWT7hGLb+5o/+r5/N8hs/Wu0f6/PKbk1cP1K22P/jI8951VXP+7IZKzMXHLN9AxO130m+/v7zRC4DxnFANXA3sF8mdtrelwKGF2a5bKry9bjSneocpiv6aI+2BqC2LF9Z6U09euy23Af09PFUegkNflyPeX7tsVP9r99+0dfMq099Hjqam8nFDUxiOjP9HXUSOS+SAs4DZgM3EvbOzvVLASNUiOWy+ZQzSpWW0tLMprGtd9/p7qDdfdXs1f7wN6b7146c6ivZP+lyW/bZ1zcu39i+9LZHrm9g8ul3VEByGjDCVkxHAFcAa4EVwEWZSEBbXgoYBSyf6l7S3alm8uLX1otbioDVQImP7Pr2riTxQcpK61jOZy5vYAoxh50HMhUwWh2t1sxeAFYCtxEM53EpMA44wcymROodKJJqBM04ni2Qap/No5FmanTSNgyP7U2ONzYm3YzRxKh+a/j52VOZeddcNmwuo6Q8j87nhAlBz+mmpuBvNnvPF+oIssWitYgCHE74ZL4k8+ZnImpFfSmHUcByeWfY2h1vHqSlYVuDz6qc57d9bqp/fr8XfYDV+xLKk+YamoYOTb5d3Wm3XSEW52YA+dCsFhieiUREfSlgFLhc/LNGvZDm+MLx4cYP/fnfzPYffarazxzwkn+F3/sSyr0R8+UM8d8O+r5Xj/2+N3bbvSLbu3YNGgwkS2cnvfi1W7LfRnPxZJGfv7wIGLl+KWAUkWxd7PKkUvSD1R/40z+d5d8fW+3j+rzi3dmyMynXl/6Xf2jdkge1xPPSv/+ew2IoF9F+rfXpKOJzq4AhhSubxSmtDXedpTvy9W+/54//v5l+/ZhqP7HXHO/Cdgf3Ehp8dM+5/o3jqv1v353h9QvWRg9qmQx+yo2kb+xQ5K2tMhUwYh3evK3GjBnjtbW1cSdDOqqtw3dnYtv9+8PWrRkb4nzVa6uZdu9b1Dyzg5pF+/HGh4fglFDGNk7YewEVR7xHxVl7cfIlI9n7gL13X7mkJLg8tdRy+PFUy0Hq6ckU2vDu2ZLqt5EoA0PA56NMDW+ugCG5F/WC2R6pLo49egTPvW4pQpDyJmfp8yuoub+OmuecaUsPYNGOgwDoxQec0m8hFcduouLcvpxw0Si69+mePo1RA2aq5czg/vujX+yzGaALSbLfRktFek4K/nkY7XmpSKpIZLueIVnxSxv6gTQ1Nvncxxb5b8c/5+PLn/cDSlfuXLyfrfPP7jvDf3FOtb90z1zfvnl729PUv39QmR2lYj7d0+aiyqc+MHFLHAKm5XlRHUarr9iDQFteChhFIo4moWmC1I6tO7z2/nn+y3Or/XNhE9fm2fuVrPJ/O/AF/59/m+pz/vqmN+5obPu+kx1vWVnq1k+JUpW1t+VinycNAfJOJ6rXUcCQwpbrf9YkF+1tJd38x71u8b14f+fk4V2W+qWH1Phdl9X44meWelNjU8f33ZELdiYu9uqz0ekpYMjuOtHdUlttWrXJn/pxrT88cpKvsv28EfMllPt4Kv3wbm/6lYc/5w98/QVf8fI72UlAR4qEMnWx1++jU1PAkF3ivIOMeiHqyAWrjeuuW7zeH/veDL9udLUf3+sNL2WHg3spO3xMz7n+zdHV/ugNM3ztm+uip6EjOppLiNJzXcFA0lDAkF3iKqNuS6/q9ga0COuunLXKH7z6Bb/qiKl+RLc3dy52MXf7KtvPmzDf0mdf3/Lff8jwCYgomwFdxU2dQwdvChQwZJe4WsF0tANaaWnr/wAp1t3Ua7BfNqLGD+6ydOfkvXjfP9X/Zf/h6dU+/+IfeVOPPLqQFnnPdsmiDNwUZCpgqB9GMYirnX0mOqA1S9GRzEtKgoG8W2jCGGhrGTt4ERUnbKXivEEc86WRdOneJVigs/Q9yGafFskPGfgtZ6ofRqvDm0sBaMMw2hkVdcjyKENub9kCkybR8GEDtffN45fnTuVz+81khQ9JunjDwP2o396HR1edyDcfG8eYiw/bFSwAli1Lvp9U0wtVPg0bL9mRR7/lWAOGmfUxs0fMbIGZzTezk+NMT8HKxjMCqqqCO5uSkuBvVdWey0QNVMmWS6Kpbhl9e2zl+EsO47rHxzF33WCeHHwZDaXd9thH2W0/o6RLmp9vZ7mQxnWzILmTT7/lTJRrtfcF3At8NXxfBvRJt7zqMHKkLWWmbWwl1WTmjVaStNx9Jfv7VUdM9QevfsFXzlrV9n209xgKnVpJFbc8qsOIM1jsAywhxcOZkr0UMHIkgxWpa99c54/eMMO/Obrax/Sc66Xs8PFU+gfs/g/Q1CXNcx/aSxdSKRZ50koqtkpvMzsGmALMA44GZgHXuPvmFstNBCYCDB06dHRda6NNyp6qqmDSpKDMc+jQoLgiXXFVBypSV9auYtp9S6h5toGaxfsxd9sIALrxISfts4CxR26g4uzejN17Nt1/dktQmVdSsud2O+NoqiJZUvCDDwJjgAbgxPDzHcCt6dYp+hxGNu6I25OdjZjDaGps8kX/Wup3XVbjlx5S48MTmrj2ZqOfOeAln/zJap/2P7P9w40fRkubmoaKZBxFkMPYF5jh7sPCz2OB77r72anWKepmtdl6ZkF7muSlSEvTnb9jbs8TmPbnVdRM70rNyuGsatoXgP62jop9FzH2hA+pOH8wR58/YvdWS21JWzM1DRXJiILPYYSBahowKnx/E/DzdMsXZQ4jcbjl9txlt5YraW+nvspKbzpwqDdhvqHHYP/F3jf5V/ndbs+h/nX/H/hvxz/ncx9b1L5RXFt7AlpbB9jT8BkiSVHold7BMXAMUAu8DjwK9E23fGwBI1sXm9aKZFq7sEcpbmpDBfaWdVt86u2v+q2fqPZP9qv1Xmzaufg3Sm5L/Rzq9koXKNuy7dbOQ2dqMSWSRFEEjLa+YgkY2bzYtPZQ+tbusqMEgzTp37h8o//j1pf9e6dU+2l7v+ZlfBjEKBr9qO4L/GtHTvWHrn3B33n13ewMQZEqYPbv37bz21raNHyGdHIKGLkS9WLTnlxIa0UyrQWmqMVNCX0gNu+9r99fPslH95zrJTQ4BKO4nthrjn/r+Gp//P/N9HWL17d/X22Vidxba2nTE+ekk1PAyJUoF5v25kJaq7to7/phMFv+0jteddXzfsVhz/lh3RbtnN2dLT6uzyv+g4pqf/qns3zTqk27jiPVxTuf79KzncNQ/YcUOAWMXIlysWnvBamjxV1J1t9e2s1/O+j7flCXup2Te7PRPz3wJf/xGdX+/G/a0MS1UOoBspn2fD5ukYgUMLIh2Z1klAtGR5+o1o6718YdjT77zwv9nydM8jWlg3d7itwAq/cv7P+i3/75qT6rcp43bGtofYNR60Py9U47W62k8jlnJRJRpgKGhjdvlq4fBKTvKZ2DobR3bNnBKw8sZNrf1lJT24Pn14xZh6q3AAAOfElEQVTkPe8LwAGl7/DRA5dQcWojYy8YwqFnDcdKrG070DDZyem8SBHIVD+MVnpWdSKTJu0eLGDnkNssXZq+89zkycmDTQdGDN26fisz71tIzeMbmDa7N9PXH8oWjgBgZNclnDdyDhXjShh74TDKTxmClezf7n0BQSBMFvSKbXTXttJ5EdklE9mUXL2yWiTV0ZY0V14ZPEGurRXXoQ11G/zJm1/y755U7af2nu1d2ebjqfQlDPVGzNd1Hewvnnmjr5q9uv3HmE4+ldXnU9FXPp0XkXZCdRgZ1pGy6nQd8FJcXNbMq/e/fOtFv/bYqX5cj3k7m7h2YbuftNfr/qfhk3xH1+65vVDlw4W6stK9rGz34y4riz9oxH1eRDogUwFDdRjNslGH0ay8nOUPTafm3iXUVDcy7e0hzN9+MADd2crJfRZQcfRGKj6zDydeNJJeg3p1nkeMtjRgAKxbt+f0/v1h7drcp0ekCBTFWFJtfe2Rw8j0nV+mW0mFr0Zs58e92eBnDXzJf3Jmtb9w5+u+bdO25GlItb1i72yW5jyKSPvQ6XMY2RrdtaUId/peXo6leb7umtLBPHjuQ4w9fzBHnTeC0rLS5AsmO6Y0+y1KlqZ1VwH9VkXySaZyGIUbMHJVZJOiWaVj/PzTzzJtVk8Gr5nNHVxLL/a80HvPnljUINZa0VZneKiQiqREMi5TAaMkE4mJRao7+jR3+u2SovlkHUP5zj/GsWjDAEoOHcWrn7iehsFDgpmlYQ6ivDx6sID0aS8vL/5gAXDHHdC16+7TunYNpotIrAo3YPTrl3x6BtrHb1y2kSdvfpnvnjSVG9+9gs303G3+duvG6rMu5d059SzYNpwp8ys47V830+XdFUFupKEh+Juu/0ZVVZCjKCkJ/lZVqW0/BOfr7ruDAGkW/L377uIPlCKFIBMVIbl6jR49OqgU7t8/eaVo167tqvhe/cYaf+T66X710VP92B7z3Gjc2cT15L1e9wcOnuSb99nXmzJZuZ6sIv3KK9M/H0Pt/0WkHeiUld7Dh3vt6tWpK4Wb6xuSNX1NsOzFldTct5SaqUET1wVhE9cebOHkvguoOPr9oInrxaPoOaBn0m10SLr6l8mTgya8qeoyir3SW0QyrnNWenfr5rXbt0dbOKwg9vFf5s3/W0LNAyuoeb6EacuGUdd4AAD7sJHTBi2kYvQWKr4wgOMuGEnZXmVZPIJQlPGJoo5hVFWVvo+IiHR6nTNgmHlbuu2tKR3MkU1zWOMDARhk9VTsv5iKk7ZT8aV9OeJzh6Ru4ppNUVp4RVkmV02LRaSgdc5WUmVtu/sf0LiGMw5ayO8vnsbCfy7h3YYB/HnFyXz9kY9y9JdGxRMsIMgF9GxR1NVysMIoy6QbMDFRsgp2EZG2ykRFSK5eow86yJt69NitIrgJvIGS5JXE+fzMgii91FtbJptPA+ysNG6UFCGKpdLbzEqBWmClu5+Tbtm9Sg/3zzd9h1v5AUNZxru2L/868CscemRXjn/mp9iHW3ct3BmKZjJVtCUBFfFJkSqmIqlrgPlRFz7gpKHMv/nPbKp7j/2b3uHiuh9ywhM3Yn/4/e5t9zvDP3mUYqtcdXAsBlGL+EQ6qVhzGGZ2AHAvMBn4Zms5jKyOVluoWmslpRxGdHq6nhSpYslh3A58G0j532hmE82s1sxq6+vrc5eyQjFhQnDhb2pK3rM8Si5EAql62qsHvggQY8Aws3OANe4+K91y7j7F3ce4+5iBAwfmKHVFZMKEoHiusxXXtYeCq0hacT7T+1Tgs2Z2FtAd2NvMKt39whjTVJwmTFCAiKL5HKkjpEhSseUw3P0Gdz/A3YcBFwDPFmWwyFUfCPW1yIzWivhEOrE4cxjFr2Uzzbq64DNk9kKUq/2ISKcWez+Mtii4VlK5aqGkllAikkaxtJIqbrnqA6G+FiKSA4UVMObMKawy+lw101RzUBHJgcIKGNu3Bx2rmsvo8z1o5KqZppqDikgOFFbASFQIQzbkqg+E+lqISA4UVqV3y+dhZHrIBj2MSESKkCq9IbNl9M1NU+vqCqPYS/0uRCTHCjeHkelhpwupaaqG4RaRNuicOYyysuyV0eeqaWomcgYahltEYlBYPb2PPBKy1XFv6NDkOYxsFHt1tEe2+l2ISAwKK4eRTblompqpnIH6XYhIDBQwmuWiaWqmcgbqdyEiMVDASJTtkUozlTNQvwsRiYECRi5lMmegYbhFJMcUMHJJOQMRKWCF1UqqGOjpdyJSoJTDEBGRSBQwREQkEgUMERGJRAFDREQiUcAQEZFIYgsYZnagmVWb2Twzm2tm18SVFhERaV2czWobgOvc/RUz6w3MMrOn3X1ejGkSEZEUYsthuPsqd38lfL8JmA8MiSs9IiKSXl7UYZjZMOBYYGaSeRPNrNbMauvr63OdNBERCcUeMMxsL+AvwLXu/n7L+e4+xd3HuPuYgQMH5j6BIiICxBwwzKwrQbCocve/xpkWERFJL85WUgb8EZjv7r+MKx0iIhJNnDmMU4GLgI+b2Wvh66wY0yMiImnE1qzW3Z8HLK79i4hI28Re6S0iIoVBAUNERCJRwBARkUgUMEREJBIFDBERiUQBQ0REIlHAEBGRSBQwREQkEgUMERGJRAFDREQiUcAQEZFIFDBERCQSBQwREYlEAUNERCJRwBARkUgUMEREJBIFDBERiUQBQ0REIlHAEBGRSBQwREQkklgDhpmdaWYLzWyxmX03zrSIiEh6sQUMMysF/gf4NHAYMN7MDosrPSIikl6cOYwTgMXu/ra7bwceBM6NMT0iIpJGlxj3PQRYnvB5BXBiy4XMbCIwMfy4zczeyEHaOmoAsDbuRESgdGZOIaQRlM5MK5R0jsrERuIMGJG4+xRgCoCZ1br7mJiT1CqlM7MKIZ2FkEZQOjOtkNKZie3EWSS1Ejgw4fMB4TQREclDcQaMl4ERZnaQmZUBFwCPx5geERFJI7YiKXdvMLOvAf8HlAJ3ufvcVlabkv2UZYTSmVmFkM5CSCMonZnWqdJp7p6J7YiISJFTT28REYlEAUNERCLJm4DR2jAhZtbNzB4K5880s2EJ824Ipy80szNiTOM3zWyemb1uZs+YWXnCvEYzey18ZbVyP0I6LzWz+oT0fDVh3iVmtih8XRJzOm9LSOObZrYhYV5OzqeZ3WVma1L1/7HAr8JjeN3MjkuYl8tz2Vo6J4Tpm2Nm083s6IR5S8Ppr2Wq+WUH0jnOzDYmfLc/SJiXs6GEIqTzWwlpfCP8PfYL5+XkfJrZgWZWHV5z5prZNUmWyezv091jfxFUer8FDAfKgNnAYS2WuQq4M3x/AfBQ+P6wcPluwEHhdkpjSuPHgJ7h+yub0xh+/iCPzuWlwK+TrNsPeDv82zd83zeudLZY/usEDSNyfT4rgOOAN1LMPwv4B2DAScDMXJ/LiOk8pXn/BMPxzEyYtxQYkCfncxzwREd/L9lOZ4tlPwM8m+vzCewHHBe+7w28meR/PaO/z3zJYUQZJuRc4N7w/SPAJ8zMwukPuvs2d18CLA63l/M0unu1u28JP84g6FuSax0ZcuUM4Gl3X+/u7wFPA2fmSTrHAw9kKS0puXsNsD7NIucC93lgBtDHzPYjt+ey1XS6+/QwHRDfbzPK+Uwlp0MJtTGdcf02V7n7K+H7TcB8ghE0EmX095kvASPZMCEtD3znMu7eAGwE+kdcN1dpTHQ5QWRv1t3Mas1shpl9LgvpaxY1neeFWdRHzKy5A2WuzmWb9hUW7R0EPJswOVfnszWpjiOX57KtWv42HXjKzGZZMBRP3E42s9lm9g8zOzyclpfn08x6Elxo/5IwOefn04Ii+mOBmS1mZfT3mfdDgxQiM7sQGAN8NGFyubuvNLPhwLNmNsfd34onhfwdeMDdt5nZfxDk3D4eU1qiuAB4xN0bE6bl0/ksGGb2MYKAcVrC5NPCczkIeNrMFoR32HF4heC7/cDMzgIeBUbElJYoPgO84O6JuZGcnk8z24sgYF3r7u9naz+QPzmMKMOE7FzGzLoA+wDrIq6bqzRiZqcDk4DPuvu25unuvjL8+zYwleBuIBtaTae7r0tI2x+A0VHXzWU6E1xAiyx/Ds9na1IdR94NfWNmRxF83+e6+7rm6Qnncg3wN7JTpBuJu7/v7h+E758EuprZAPLwfIbS/Tazfj7NrCtBsKhy978mWSSzv89sV8xErLzpQlDpchC7KrQOb7HMf7J7pffD4fvD2b3S+22yU+kdJY3HElTMjWgxvS/QLXw/AFhElirsIqZzv4T3nwdm+K6KsCVhevuG7/vFlc5wuUMJKhEtjvMZ7mMYqStpz2b3SsWXcn0uI6ZzKEH93iktpvcCeie8nw6cGWM6923+rgkutMvCcxvp95KrdIbz9yGo5+gVx/kMz8t9wO1plsno7zNrJ7sdB38WQS3/W8CkcNotBHfqAN2BP4c/+peA4QnrTgrXWwh8OsY0/gtYDbwWvh4Pp58CzAl/5HOAy2M+lz8G5obpqQYOTVj3K+E5XgxcFmc6w883AT9psV7OzifB3eMqYAdBOe/lwBXAFeF8I3gQ2FthWsbEdC5bS+cfgPcSfpu14fTh4XmcHf4mJsWczq8l/DZnkBDgkv1e4kpnuMylBA1uEtfL2fkkKFZ04PWE7/WsbP4+NTSIiIhEki91GCIikucUMEREJBIFDBERiUQBQ0REIlHAEBGRSBQwREQkEgUMERGJRAFDpAPC5xF8Mnz/QzP777jTJJItGnxQpGNuBG4JB5o7FvhszOkRyRr19BbpIDN7DtgLGOfBcwlEipKKpEQ6wMyOJHjy2XYFCyl2Chgi7RQ+uayK4KlmH5hZ1p6oJ5IPFDBE2iF80tpfgevcfT5wK0F9hkjRUh2GiIhEohyGiIhEooAhIiKRKGCIiEgkChgiIhKJAoaIiESigCEiIpEoYIiISCT/H1a0aT1GwF0xAAAAAElFTkSuQmCC\n",
- "text/plain": [
- ""
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
+ "execution_count": 2,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
@@ -509,7 +961,9 @@
{
"cell_type": "code",
"execution_count": 3,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"# Importing various packages\n",
@@ -593,7 +1047,9 @@
{
"cell_type": "code",
"execution_count": 4,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"import numpy as np\n",
@@ -739,7 +1195,9 @@
{
"cell_type": "code",
"execution_count": 5,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"import numpy as np \n",
@@ -801,7 +1259,9 @@
{
"cell_type": "code",
"execution_count": 6,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"import numpy as np \n",
@@ -1279,25 +1739,7 @@
]
}
],
- "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.7.0"
- }
- },
+ "metadata": {},
"nbformat": 4,
"nbformat_minor": 2
}
diff --git a/doc/pub/Splines/ipynb/ipynb-Splines-src.tar.gz b/doc/pub/Splines/ipynb/ipynb-Splines-src.tar.gz
index 34845ba2e..9e4fe066b 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 b63252b15..b40638855 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 950f1cad9..51367b610 100644
--- a/doc/src/Splines/Splines.do.txt
+++ b/doc/src/Splines/Splines.do.txt
@@ -15,6 +15,203 @@ the cost function. Ideally we would be able to solve for $\beta$
analytically, however this is not possible in general and we must use
some approximative/numerical method to compute the minimum.
+
+!split
+===== Revisiting our Logistic Regression case =====
+
+In our discussion on Logistic Regression we defined we studied first the
+case of
+two classes, with $y_i$ either
+$0$ or $1$. Furthermore we assumed also that we have only two
+parameters $\beta$ in our fitting of the Sigmoid function, that is we
+defined probabilities
+
+!bt
+\begin{align*}
+p(y_i=1|x_i,\hat{\beta}) &= \frac{\exp{(\beta_0+\beta_1x_i)}}{1+\exp{(\beta_0+\beta_1x_i)}},\nonumber\\
+p(y_i=0|x_i,\hat{\beta}) &= 1 - p(y_i=1|x_i,\hat{\beta}),
+\end{align*}
+!et
+where $\hat{\beta}$ are the weights we wish to extract from data, in our case $\beta_0$ and $\beta_1$.
+
+!split
+===== The equations to solve =====
+
+Our compact equations used a definition of a vector $\hat{y}$ with $n$
+elements $y_i$, an $n\times p$ matrix $\hat{X}$ which contains the
+$x_i$ values and a vector $\hat{p}$ of fitted probabilities
+$p(y_i\vert x_i,\hat{\beta})$. We rewrote in a more compact form
+the first derivative of the cost function as
+
+!bt
+\[
+\frac{\partial \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}} = -\hat{X}^T\left(\hat{y}-\hat{p}\right).
+\]
+!et
+
+If we in addition define a diagonal matrix $\hat{W}$ with elements
+$p(y_i\vert x_i,\hat{\beta})(1-p(y_i\vert x_i,\hat{\beta})$, we can obtain a compact expression of the second derivative as
+
+!bt
+\[
+\frac{\partial^2 \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}\partial \hat{\beta}^T} = \hat{X}^T\hat{W}\hat{X}.
+\]
+!et
+This defines what we call the Hessian.
+
+!split
+===== Solving using Newton-Raphson's method =====
+
+If we can set up these equations, Newton-Raphson's iterative method is the nomrally the method of choice. It requires however that we setting the matrices that define the first and second derivatives.
+
+Our iterative scheme is then given by
+
+!bt
+\[
+\hat{\beta}^{\mathrm{new}} = \hat{\beta}^{\mathrm{old}}-\left(\frac{\partial^2 \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}\partial \hat{\beta}^T}\right)^{-1}\times \left(\frac{\partial \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}}\right)_{\hat{\beta}^{\mathrm{old}}},
+\]
+!et
+or in matrix form as
+
+!bt
+\[
+\hat{\beta}^{\mathrm{new}} = \hat{\beta}^{\mathrm{old}}-\left(\hat{X}^T\hat{W}\hat{X} \right)^{-1}\times \left(-\hat{X}^T(\hat{y}-\hat{p}) \right)_{\hat{\beta}^{\mathrm{old}}}.
+\]
+!et
+The right-hand side is computed with the old values of $\beta$.
+
+If we can compute these matrices, in particular the Hessian, the above is often the easiest method to implement.
+
+
+!split
+===== Brief reminder on Newton-Raphson's method =====
+
+Let us quicly remind ourselves how we derive the above method.
+
+Perhaps the most celebrated of all one-dimensional root-finding
+routines is Newton's method, also called the Newton-Raphson
+method. This method is distinguished from the previously discussed
+methods by the fact that it requires the evaluation of both the
+function $f$ and its derivative $f'$ at arbitrary points. In this
+sense, it is taylored to cases with e.g., transcendental equations.
+If you can only calculate the derivative
+numerically and/or your function is not of the smooth type, we
+discourage the use of this method.
+
+!split
+===== The equations =====
+
+The Newton-Raphson formula consists geometrically of extending the
+tangent line at a current point until it crosses zero, then setting
+the next guess to the abscissa of that zero-crossing. The mathematics
+behind this method is rather simple. Employing a Taylor expansion for
+$x$ sufficiently close to the solution $s$, we have
+
+
+!bt
+\[
+ f(s)=0=f(x)+(s-x)f'(x)+\frac{(s-x)^2}{2}f''(x) +\dots.
+ \label{eq:taylornr}
+\]
+!et
+
+For small enough values of the function and for well-behaved
+functions, the terms beyond linear are unimportant, hence we obtain
+
+
+!bt
+\[
+ f(x)+(s-x)f'(x)\approx 0,
+\]
+!et
+yielding
+!bt
+\[
+ s\approx x-\frac{f(x)}{f'(x)}.
+\]
+!et
+
+Having in mind an iterative procedure, it is natural to start iterating with
+!bt
+\[
+ x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}.
+\]
+!et
+
+!split
+===== Simple geometric interpretation =====
+
+The above is Newton-Raphson's method. It has a simple geometric
+interpretation, namely $x_{n+1}$ is the point where the tangent from
+$(x_n,f(x_n))$ crosses the $x-$axis. Close to the solution,
+Newton-Raphson converges fast to the desired result. However, if we
+are far from a root, where the higher-order terms in the series are
+important, the Newton-Raphson formula can give grossly inaccurate
+results. For instance, the initial guess for the root might be so far
+from the true root as to let the search interval include a local
+maximum or minimum of the function. If an iteration places a trial
+guess near such a local extremum, so that the first derivative nearly
+vanishes, then Newton-Raphson may fail totally
+
+
+!split
+===== Extending to more than one variable =====
+
+Newton's method can be generalized to systems of several non-linear equations
+and variables. Consider the case with two equations
+!bt
+\[
+ \begin{array}{cc} f_1(x_1,x_2) &=0\\
+ f_2(x_1,x_2) &=0\end{array},
+\]
+!et
+which we Taylor expand to obtain
+
+!bt
+\[
+ \begin{array}{cc} 0=f_1(x_1+h_1,x_2+h_2)=&f_1(x_1,x_2)+h_1
+ \partial f_1/\partial x_1+h_2
+ \partial f_1/\partial x_2+\dots\\
+ 0=f_2(x_1+h_1,x_2+h_2)=&f_2(x_1,x_2)+h_1
+ \partial f_2/\partial x_1+h_2
+ \partial f_2/\partial x_2+\dots
+ \end{array}.
+\]
+!et
+Defining the Jacobian matrix ${\bf \hat{J}}$ we have
+!bt
+\[
+ {\bf \hat{J}}=\left( \begin{array}{cc}
+ \partial f_1/\partial x_1 & \partial f_1/\partial x_2 \\
+ \partial f_2/\partial x_1 &\partial f_2/\partial x_2
+ \end{array} \right),
+\]
+!et
+we can rephrase Newton's method as
+!bt
+\[
+\left(\begin{array}{c} x_1^{n+1} \\ x_2^{n+1} \end{array} \right)=
+\left(\begin{array}{c} x_1^{n} \\ x_2^{n} \end{array} \right)+
+\left(\begin{array}{c} h_1^{n} \\ h_2^{n} \end{array} \right),
+\]
+!et
+where we have defined
+!bt
+\[
+ \left(\begin{array}{c} h_1^{n} \\ h_2^{n} \end{array} \right)=
+ -{\bf \hat{J}}^{-1}
+ \left(\begin{array}{c} f_1(x_1^{n},x_2^{n}) \\ f_2(x_1^{n},x_2^{n}) \end{array} \right).
+\]
+!et
+We need thus to compute the inverse of the Jacobian matrix and it
+is to understand that difficulties may
+arise in case ${\bf \hat{J}}$ is nearly singular.
+
+It is rather straightforward to extend the above scheme to systems of
+more than two non-linear equations. In our case, the Jacobian matrix is given by the Hessian that represents the second derivative of cost function.
+
+
+
!split
===== Steepest descent =====
@@ -178,6 +375,118 @@ o A norm is any function that satisfy the following properties
Using the definition of convexity, try to show that a function satisfying the properties above is convex (the third condition is not needed to show this).
+
+!split
+===== Standard steepest descent =====
+
+
+Before we proceed, we would like to mention the approach called the _standard Steepest descent_, which again leads to us having to be able to compute a matrix.
+
+"The success of the CG method":"https://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf"
+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$.
+
+!split
+===== Conjugate gradient method =====
+
+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.
+
+More details will be added here soon.
+
+!split
+===== Simple codes for steepest descent and conjugate gradient using a $2\times 2$ matrix, in c++, Python code to come =====
+!bblock
+!bc cppcod
+#include
+#include
+#include
+#include
+#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;
+}
+!ec
+!eblock
+
+!split
+===== The routine for the steepest descent method =====
+!bblock
+!bc cppcod
+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;
+}
+!ec
+!eblock
+
+
!split
===== Revisiting our first homework =====
@@ -839,3 +1148,4 @@ which gives
+