# Newton - Raphson method : maxima of function in 2 variables

I am computing the maximum of a function (with two-variables) using Newton-Raphson method. The function is : $e^{-(x \ - x_0)^2 - (y \ - y_0)^2}$, whose maxima exists at $(x_0,y_0)$.

The Jacobian will be a function of the second derivatives (including mixed) of this function. It can be written as : $J = e^{-(x \ - x_o)^2 - (y \ - y_o)^2} \ \begin{Bmatrix} -2 \ + \ 4 \times (x \ - \ x_o)^2 && 4 \times (x \ - \ x_o) \ (y \ - \ y_o) \\ 4 \times (x \ - \ x_o) \ (y \ - \ y_o) && -2 \ + \ 4 \times (y \ - \ y_o)^2 \end{Bmatrix}$

The above expression analytically calculates the second derivatives to evaluate the Jacobian. The Newton-Raphson can now be formulated as :

$$\vec{x}^{\, n+1} \ = \ \vec{x}^{\, n} \ - \ J^{-1} \ f(\vec{x}^{\,n})$$

where $\vec{x}$ includes $x$ and $y$.

An alternative way to evaluate the Jacobian (second derivatives) includes using Finite Differences like : $J_{i,j} (\vec{x}) \ = \ \frac{f_i(\vec{x} \ + \ h_j \vec{e_j}) \ - \ f_i(\vec{x})}{h_j}$, where $h_j$ is approximately taken to be $x_j \ \times \sqrt{\epsilon}$, with $\epsilon$ being the computer precision (generally taken as $10^{-16}$). This can be roughly called as Secant Method (though not exactly, as it means).

Practically, analytical way of calculating the Jacobian should yield faster convergence (in a few iterations) in comparison to the aternative way described above. However, I seem to get the converse result (for the guess values ($x_0$ and $y_0$ respectively) taken to be 3.6 and 4.6 in the code below !). In other words, I seem to get the converged solution (for the case of Analytical Jacobi) after 26 iterations, whereas 15 for the other case.

Here is the code, where I have implemented this.

#include<iostream>
#include<cmath>
#include<fstream>
#include<iomanip>

using namespace std;

double eps = 1e-16;

void user_input(double x0, double y0)
{
cout << "Enter the value of x0 - to find the maximum of the exponential function\n" << endl;
cin >> x0;
cout << "Enter the value of y0 - to find the maximum of the exponential function\n" << endl;
cin >> y0;

}

//evaluates the exponential at x0 and y0
double exp_fn(double x, double y, double x0, double y0)
{
double val;
val = exp(-(x - x0)*(x - x0) - (y - y0)*(y - y0));
//cout << val << endl;

return val;
}

//x- gradient of the given function
double f_x(double x, double y, double x0, double y0)
{

double value;
value = -2*exp_fn(x,y,x0,y0)*(x - x0);

return (value);
}

//y- gradient of the given function
double f_y(double x, double y, double x0, double y0)
{

double value;
value = -2*exp_fn(x,y,x0,y0)*(y - y0);

return (value);
}

//first component in Jacobian
double f_xx(double x, double y, double x0, double y0)
{
double value;
value = exp_fn(x,y,x0,y0)*( -2 + 4*(x - x0)*(x - x0));

return value;
}

//second component in Jacobian
double f_xy(double x, double y, double x0, double y0) //returns the same value as f_yx
{
double value;
value = exp_fn(x,y,x0,y0)*(x - x0)*(y - y0)*4;

return value;
}

//last component in Jacobian
double f_yy(double x, double y, double x0, double y0)
{
double value;
value = exp_fn(x,y,x0,y0)*( -2 + 4*(y - y0)*(y - y0));

return value;
}

//equivalent to f_xx, but using finite difference, not analytical
double f_xx_sec(double x, double y, double x0, double y0)
{
double value, h_x;
h_x = x*sqrt(eps);

value = (f_x(x + h_x, y, x0, y0) - f_x(x, y, x0, y0))/h_x;

return value;
}

//equivalent to f_xy, but using finite difference, not analytical
double f_xy_sec(double x, double y, double x0, double y0)
{
double value, h_y;
h_y = y*sqrt(eps);

value = (f_x(x, y + h_y, x0, y0) - f_x(x, y, x0, y0))/h_y;

return value;
}

//equivalent to f_yx, but using finite difference, not analytical
double f_yx_sec(double x, double y, double x0, double y0)
{
double value, h_x;
h_x = x*sqrt(eps);

value = (f_y(x + h_x, y, x0, y0) - f_y(x, y, x0, y0))/h_x;

return value;
}

//equivalent to f_yy, but using finite difference, not analytical
double f_yy_sec(double x, double y, double x0, double y0)
{
double value, h_y;
h_y = y*sqrt(eps);

value = (f_y(x, y + h_y, x0, y0) - f_y(x, y, x0, y0))/h_y;

return value;
}

//determinant of the matrix
double det(double a, double b, double c, double d)
{
double det_val;
det_val = (a*d - b*c);
return det_val;

}

int main(int argv, char** argc)
{

int user = 0, iter = 1; //user = 0 determines the default x0 and y0 values to be used.
//Else, for user to input these values, equate user = 1

int secant_met = 1; //1 causes secant method of Jacobian computation. 0 uses the conventional way

double x0 = 4, y0 = 7, x_old =3.6, y_old =6.7, error = 1., magnitude; //x0, y0 - guess  values

double x_new, y_new, mult_x, mult_y, lambda = 1.0; // lambda refers to the relaxation factor

ofstream fout;
fout.open("convergence_data.dat");
fout << "Error" << "\t" << "X_updated" << "\t" << "Y_updated" << "\t" << "Iteration" << endl;

if (user == 1)
user_input(x0,y0);

while(error > 1e-8)
{

if (secant_met ==1) //using the secant way of evaluating Jacobian - Finite differences
{
mult_x = f_x(x_old, y_old, x0, y0)*f_yy_sec(x_old, y_old, x0, y0) - f_y(x_old, y_old, x0, y0)*f_xy_sec(x_old, y_old, x0, y0);
mult_y = f_y(x_old, y_old, x0, y0)*f_xx_sec(x_old, y_old, x0, y0) - f_x(x_old, y_old, x0, y0)*f_yx_sec(x_old, y_old, x0, y0);

x_new = x_old - lambda*mult_x/det(f_xx_sec(x_old, y_old, x0, y0), f_xy_sec(x_old, y_old, x0, y0), f_yx_sec(x_old, y_old, x0, y0), f_yy_sec(x_old, y_old, x0, y0)); //Generalized newton's method
y_new = y_old - lambda*mult_y/det(f_xx_sec(x_old, y_old, x0, y0), f_xy_sec(x_old, y_old, x0, y0), f_yx_sec(x_old, y_old, x0, y0), f_yy_sec(x_old, y_old, x0, y0));
}

else // implementing conventional way of computing Jacobian
{
mult_x = f_x(x_old, y_old, x0, y0)*f_yy(x_old, y_old, x0, y0) - f_y(x_old, y_old, x0, y0)*f_xy(x_old, y_old, x0, y0);
mult_y = f_y(x_old, y_old, x0, y0)*f_xx(x_old, y_old, x0, y0) - f_x(x_old, y_old, x0, y0)*f_xy(x_old, y_old, x0, y0);

cout << setprecision(7) << f_xx(x_old, y_old, x0, y0) << "\t" << setprecision(7) << f_xy(x_old, y_old, x0, y0) << "\t" << setprecision(7)<< f_yy(x_old, y_old, x0, y0) << endl;

x_new = x_old - lambda*mult_x/det(f_xx(x_old, y_old, x0, y0), f_xy(x_old, y_old, x0, y0), f_xy(x_old, y_old, x0, y0), f_yy(x_old, y_old, x0, y0)); //Generalized newton's method
y_new = y_old - lambda*mult_y/det(f_xx(x_old, y_old, x0, y0), f_xy(x_old, y_old, x0, y0), f_xy(x_old, y_old, x0, y0), f_yy(x_old, y_old, x0, y0));
}

magnitude = fabs( (x_new - x_old)*(x_new - x_old) + (y_new - y_old)*(y_new - y_old));
error = sqrt(magnitude); //calculating the error for convergence

x_old = x_new;
y_old = y_new;

cout << fixed;
cout << "Error = " <<  error << "\t" << "Iteration" << "\t" <<  iter << endl << "X updated" << "\t" << x_new << "\t" << "Y_updated" << "\t" << y_new << endl;
fout << fixed;
fout << error << "\t" << x_new << "\t" << y_new << "\t" << iter << endl;

iter++;
}

fout.close();
return 0;
}


Can someone give me an insight as to where I am going wrong in the code, wherein I have implemented both these methods ? Any kind of help would be deeply appreciated.

EDIT:

This is the plot of convergence error vs iterations which I seem to get when using both the methods, as described above. • Is this true for all starting positions, or just for your particular choice? And is it also true if you reduce the tolerance to something smaller than 1e-8? – Wolfgang Bangerth Nov 28 '17 at 15:34
• @WolfgangBangerth I tried changing $x_0$ and $y_0$ to other values, and faced the same problem. For this case (with $x_0$ and $y_0$ of 4 and 7, my N-R solver diverges for guess values below 3.6 and 6.7 respectively (as a guess close to the true solution is necessary). At these guess values, I face the problem as mentioned in the post (it remains the same even if I reduce the tolerance to, say 1e-5). However, for guess values exceeding 3.6 and 6.7, both methods converge for the same number of iterations (4-5 I think). – Sthavishtha Bhopalam Nov 28 '17 at 19:40
• I'm asking this because (i) the convergence rates for Newton are in some sense averages over all starting points [that's not entirely correct, but you can occasionally get lucky with a secant method and be better than Newton], and (ii) they are asymptotic, i.e. only cover the case where you make the tolerance very small. You cannot consider a single starting point and a large tolerance and draw any inferences. In particular, if the iteration diverges because you start from too far away, then there is no inference to be drawn. – Wolfgang Bangerth Nov 28 '17 at 23:24
• I do agree with the point that one cannot infer anything at a specific point, when divergence occurs at guess values below these (or both the methods exhibit same behavior beyond this value). However, in my case, even if i make the tolerance error very small, I don't think that the results will change, as you can see from the figure (updated in the post). Surprisingly, the N-R solver using analytical expression of Jacobian takes a larger number of iterations to yield error close to 1 (it was giving 0.998 - 0.995 around) when compared to the secant method. – Sthavishtha Bhopalam Nov 29 '17 at 16:09
• As N-R method has a second order convergence and the other method has less than second order (due to the approximation induced of the Jacobian terms), I thought that N-R will require less number of iterations than the other method. Thanks for clearing this confusion, that the convergence rate cannot be related to the number of iterations, but rather the slope of the curve. – Sthavishtha Bhopalam Nov 30 '17 at 8:21