Newton_minimization_2D.html

 <!doctype html>
<html class="no-js" lang="en">
<head>
	<meta charset="utf-8">
	<meta name="viewport" content="width=device-width, initial-scale=1">
	<link rel="stylesheet" href="https://interactivecomputergraphics.github.io/physics-simulation/examples/style.css">
	<script type="text/x-mathjax-config">
	  MathJax.Hub.Config({
		extensions: ["tex2jax.js"],
		jax: ["input/TeX", "output/HTML-CSS"],
		tex2jax: {
		  inlineMath: [ ['$','$'], ["\\(","\\)"] ],
		  displayMath: [ ['$$','$$'], ["\\[","\\]"] ],
		  processEscapes: true
		},
		"HTML-CSS": { fonts: ["TeX"] }
	  });
	</script>
	<script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.4/MathJax.js"></script>
	<script src="https://cdn.plot.ly/plotly-2.5.1.min.js"></script>
	<title>Newton's method for 2D minimization</title>
</head>
<body>

<header class="page-header">
	<h1>Newton's method for 2D minimization</h1>
</header>

<main>
	<!-- Simulation panel: canvas + controls -->
	<div class="card sim-panel">
		<div class="sim-canvas-wrap">
			<div id="plotOutput" style="border:0px solid #000000;border-radius: 0px;background-color:#EEEEEE"></div>
				<p style="margin:6px 0 0 4px;color:#555;font-size:0.9em">Click anywhere on the plot to set the starting point x&#x2080;.</p>
			</div>
	</div>
	<div class="controls-panel" style="width: 100%;align:center;margin-left:auto;margin-right:auto">
		<h3>Controls</h3>
		<div class="controls-grid" style="width: 400px;align:left">
			<label>Newton steps</label>
			<p style="display:flex;align-items:center;gap:6px">
				<button onclick="plot.changeSteps(-1)" style="width:30px;height:30px;font-size:16px;cursor:pointer">&#8722;</button>
				<input type="text" id="textInput" value="1" readonly style="width:36px;text-align:center">
				<button onclick="plot.changeSteps(+1)" style="width:30px;height:30px;font-size:16px;cursor:pointer">&#43;</button>
			</p>
			
			<label>Starting point x&#x2080;</label>
			<span id="x0Display" style="font-family:monospace"></span>
			
			<label for="fct">Function</label>
			<select onchange="plot.reset()" id="fct" size="1">
				<option selected="selected">Elliptic quadratic</option>
				<option>Rosenbrock function</option>
				<option>Himmelblau function</option>
			</select>
		</div>
	</div>

	<!-- Theory section -->
	<div class="card theory">
			<h2>Newton's method for 2D minimization</h2>

			<p>In two dimensions, Newton's method minimizes $f(\mathbf{x})$, $\mathbf{x} = (x_1, x_2)^\top$, by
			iteratively building a local quadratic (second-order Taylor) model and jumping to its minimizer.</p>

			<p>The quadratic model at $\mathbf{x}_n$ is:</p>
			$$q(\mathbf{x}) = f(\mathbf{x}_n) + \nabla f(\mathbf{x}_n)^\top (\mathbf{x}-\mathbf{x}_n)
			+ \frac{1}{2}(\mathbf{x}-\mathbf{x}_n)^\top \mathbf{H}(\mathbf{x}_n)(\mathbf{x}-\mathbf{x}_n)$$
			<p>where $\nabla f$ is the gradient and $\mathbf{H}$ is the Hessian matrix. Setting $\nabla q = 0$ gives the update:</p>
			$$\mathbf{x}_{n+1} = \mathbf{x}_n - \mathbf{H}(\mathbf{x}_n)^{-1} \nabla f(\mathbf{x}_n)$$

			<p>The <strong>ellipses</strong> drawn at each iterate are level curves of the local quadratic model
			$q(\mathbf{x}) = \text{const}$ centred at $\mathbf{x}_n$. Their shape is entirely determined by the
			Hessian: a well-conditioned Hessian gives nearly circular contours and fast convergence;
			an ill-conditioned one gives elongated ellipses and slower convergence.</p>

			<h3>Example functions</h3>

			<p><strong>Elliptic quadratic</strong> — a simple convex bowl with different curvatures along each axis.
			The Hessian is constant, so Newton's method converges in a single step from any starting point.
			$$f(x_1, x_2) = 3(x_1-1)^2 + 8(x_2-0.5)^2, \qquad \mathbf{x}^* = (1,\; 0.5)$$</p>

			<p><strong>Rosenbrock function</strong> — the classic "banana valley" test problem.
			The minimum lies inside a narrow curved valley, producing a highly ill-conditioned Hessian
			and very elongated ellipses. Pure Newton's method still converges quickly once inside the valley.
			$$f(x_1, x_2) = 100(x_2 - x_1^2)^2 + (1 - x_1)^2, \qquad \mathbf{x}^* = (1,\; 1)$$</p>

			<p><strong>Himmelblau's function</strong> — a non-convex function with four minima.
			The starting point determines which minimum Newton's method converges to; near saddle points
			the Hessian is indefinite and the step may point away from any minimum.
			$$f(x_1, x_2) = (x_1^2 + x_2 - 11)^2 + (x_1 + x_2^2 - 7)^2$$
			$$\mathbf{x}^* \in \{(3,2),\;(-2.805,\,3.131),\;(-3.779,\,-3.283),\;(3.584,\,-1.848)\}$$</p>
	</div>	
	
</main>

<script id="simulation_code" type="text/javascript">
	// ---------------------------------------------------------------------------
	// Small 2×2 matrix / linear-algebra helpers
	// ---------------------------------------------------------------------------

	// Solve H*d = g for d  (Cramer's rule, H is 2×2)
	function solve2x2(H, g)
	{
		const a = H[0][0], b = H[0][1], c = H[1][0], d = H[1][1];
		const det = a*d - b*c;
		if (Math.abs(det) < 1e-14) return [0, 0];
		return [(d*g[0] - b*g[1]) / det,
		        (a*g[1] - c*g[0]) / det];
	}

	// Eigenvalues of a symmetric 2×2 matrix
	function eigenvalues2x2sym(H)
	{
		const a = H[0][0], b = H[0][1], d = H[1][1];
		const tr   = a + d;
		const disc = Math.sqrt(Math.max(0, (a-d)*(a-d) + 4*b*b));
		return [(tr + disc) / 2, (tr - disc) / 2];
	}

	// Unit eigenvectors of a symmetric 2×2 matrix; returns [v1, v2]
	function eigenvectors2x2sym(H)
	{
		const a = H[0][0], b = H[0][1], d = H[1][1];
		const disc = Math.sqrt(Math.max(0, (a-d)*(a-d) + 4*b*b));
		const lam1 = ((a + d) + disc) / 2;
		let v1;
		if (Math.abs(b) > 1e-12)
			v1 = [lam1 - d, b];
		else
			v1 = (Math.abs(a) >= Math.abs(d)) ? [1, 0] : [0, 1];
		const len = Math.sqrt(v1[0]*v1[0] + v1[1]*v1[1]);
		v1 = [v1[0]/len, v1[1]/len];
		const v2 = [-v1[1], v1[0]];
		return [v1, v2];
	}

	// Parametric ellipse centred at xn1 that passes exactly through xn.
	// d = xn - xn1.  We use the ABSOLUTE eigenvalues of H so the ellipse is
	// valid even when H is indefinite or negative definite — the shape still
	// reflects the local curvature magnitude.
	function quadraticLevelEllipse(xn1, H, d)
	{
		const lambdas = eigenvalues2x2sym(H);
		const vecs    = eigenvectors2x2sym(H);

		// Absolute eigenvalues: valid for any symmetric H
		const absL1 = Math.abs(lambdas[0]);
		const absL2 = Math.abs(lambdas[1]);
		if (absL1 <= 1e-12 || absL2 <= 1e-12) return null;

		// Project d onto the eigenvectors
		const p1 = d[0]*vecs[0][0] + d[1]*vecs[0][1];
		const p2 = d[0]*vecs[1][0] + d[1]*vecs[1][1];

		// c = 0.5*(|λ1|*p1² + |λ2|*p2²): the level at which the ellipse passes through xn
		const c = 0.5 * (absL1*p1*p1 + absL2*p2*p2);
		if (c <= 1e-14) return null;   // d ≈ 0: iterates have converged

		const r1 = Math.sqrt(2 * c / absL1);
		const r2 = Math.sqrt(2 * c / absL2);
		if (!isFinite(r1) || !isFinite(r2)) return null;

		const ex = [], ey = [];
		const N = 120;
		for (let k = 0; k <= N; k++)
		{
			const theta = 2 * Math.PI * k / N;
			const u = r1 * Math.cos(theta);
			const v = r2 * Math.sin(theta);
			ex.push(xn1[0] + u*vecs[0][0] + v*vecs[1][0]);
			ey.push(xn1[1] + u*vecs[0][1] + v*vecs[1][1]);
		}
		return { x: ex, y: ey };
	}

	// ---------------------------------------------------------------------------
	// 2-D test functions — each returns { f, grad:[gx,gy], H:[[…],[…]] }
	// ---------------------------------------------------------------------------

	// f = 3(x-1)^2 + 8(y-0.5)^2  — elongated quadratic bowl, minimum at (1, 0.5)
	// The condition number 8/3 ≈ 2.7 produces visibly non-circular ellipses.
	function ellipticQuadratic(x, y)
	{
		const dx = x - 1.0, dy = y - 0.5;
		return {
			f:    3*dx*dx + 8*dy*dy,
			grad: [6*dx, 16*dy],
			H:    [[6, 0], [0, 16]]
		};
	}

	// Rosenbrock: f = 100(y-x^2)^2 + (1-x)^2  — banana valley, minimum at (1,1)
	// Ill-conditioned Hessian gives very elongated ellipses; convergence takes ~6 steps.
	function rosenbrock(x, y)
	{
		const A = y - x*x, B = 1 - x;
		return {
			f:    100*A*A + B*B,
			grad: [-400*x*A - 2*B,  200*A],
			H:    [[-400*(A - 2*x*x) + 2, -400*x],
			       [-400*x,                 200]]
		};
	}

	// Himmelblau: f = (x^2+y-11)^2 + (x+y^2-7)^2 — four minima, non-convex landscape
	function himmelblau(x, y)
	{
		const A = x*x + y - 11, B = x + y*y - 7;
		return {
			f:    A*A + B*B,
			grad: [4*x*A + 2*B,  2*A + 4*y*B],
			H:    [[12*x*x + 4*y - 44 + 2,  4*x + 4*y],
			       [4*x + 4*y,   2 + 12*y*y + 4*x - 28]]
		};
	}

	// ---------------------------------------------------------------------------
	class Plot
	{
		constructor()
		{
			this.x0 = null;
			this.currentFct = null;
			this.reset();
		}

		defaultX0(fct)
		{
			if (fct === 'Elliptic quadratic')  return [-1.5, -1.5];
			if (fct === 'Rosenbrock function') return [-1.2,  0.8];
			if (fct === 'Himmelblau function') return [ 0.0,  0.0];
			return [0.0, 0.0];
		}

		updateX0Display()
		{
			const el = document.getElementById('x0Display');
			if (el) el.textContent = `(${this.x0[0].toFixed(3)}, ${this.x0[1].toFixed(3)})`;
		}

		reset()
		{
			this.num_newton_steps = parseInt(document.getElementById('textInput').value);
			const newFct = document.getElementById('fct').value;
			// Reset starting point only when the function changes
			if (newFct !== this.currentFct)
			{
				this.currentFct = newFct;
				this.x0 = this.defaultX0(newFct);
			}
			this.fct = newFct;
			this.plotFunctions();
		}

		changeSteps(delta)
		{
			const el = document.getElementById('textInput');
			const next = Math.max(1, Math.min(20, parseInt(el.value) + delta));
			el.value = next;
			this.reset();
		}

		// Newton step in 2D:  x_{n+1} = x_n - H^{-1} grad f
		newton_step_2d(fctFn, xn)
		{
			const { grad, H } = fctFn(xn[0], xn[1]);
			const d = solve2x2(H, grad);
			return [xn[0] - d[0], xn[1] - d[1]];
		}

		addMinima(minima, data)
		{
			data.push({
				x: minima.map(m => m[0]),
				y: minima.map(m => m[1]),
				mode: 'markers',
				marker: { color: 'red', size: 12, symbol: 'star',
				          line: { color: 'white', width: 1.5 } },
				name: 'minimum',
				showlegend: true
			});
		}

		computeData(fctFn, x0, xRange, yRange, data)
		{
			const NX = 200, NY = 200;
			const [xlo, xhi] = xRange;
			const [ylo, yhi] = yRange;

			// Build contour grid (z[j][i] = f at (xGrid[i], yGrid[j]))
			const xGrid = [], yGrid = [], zGrid = [];
			for (let i = 0; i < NX; i++)
				xGrid.push(xlo + (xhi - xlo) * i / (NX - 1));
			for (let j = 0; j < NY; j++)
			{
				const y = ylo + (yhi - ylo) * j / (NY - 1);
				yGrid.push(y);
				const row = [];
				for (let i = 0; i < NX; i++)
					row.push(fctFn(xGrid[i], y).f);
				zGrid.push(row);
			}

			// Height-map background: smooth colour gradient (blue=low, red=high)
			// with contour lines overlaid for orientation
			data.push({
				type: 'contour',
				x: xGrid, y: yGrid, z: zGrid,
				colorscale: 'Jet',
				contours: { coloring: 'heatmap', showlines: true,
				            line: { color: 'rgba(0,0,0,0.25)', width: 0.8 } },
				showscale: true,
				colorbar: { title: 'f(x₁,x₂)' },
				ncontours: 24,
				name: 'f(x₁,x₂)',
				showlegend: false
			});

			// Compute Newton iterates
			const iterates = [x0.slice()];
			let cur = x0.slice();
			for (let i = 0; i < this.num_newton_steps; i++)
			{
				cur = this.newton_step_2d(fctFn, cur);
				iterates.push(cur.slice());
			}

			// Ellipses: level curve of the local curvature model at each iterate.
			// Centre = x_{n+1}; the ellipse passes exactly through x_n.
			// Absolute eigenvalues of H are used so the ellipse is always valid,
			// even when H is indefinite or negative definite.
			for (let i = 0; i < iterates.length - 1; i++)
			{
				const xn    = iterates[i];
				const xn1   = iterates[i + 1];
				const { H } = fctFn(xn[0], xn[1]);
				const d     = [xn[0] - xn1[0], xn[1] - xn1[1]];

				const ellipse = quadraticLevelEllipse(xn1, H, d);
				if (ellipse)
				{
					data.push({
						x: ellipse.x, y: ellipse.y,
						mode: 'lines',
						line: { color: 'rgb(255,255,0)', width: 2.5 },
						name: 'quadratic model',
						legendgroup: 'quadratic model',
						showlegend: i === 0
					});
				}
			}

			// Dotted path connecting iterates
			data.push({
				x: iterates.map(p => p[0]),
				y: iterates.map(p => p[1]),
				mode: 'lines',
				line: { color: 'white', width: 1.5, dash: 'dot' },
				name: 'Newton path',
				showlegend: true
			});

			// Iterate dots with subscript labels
			data.push({
				x: iterates.map(p => p[0]),
				y: iterates.map(p => p[1]),
				mode: 'markers+text',
				marker: { color: 'white', size: 9, symbol: 'circle' },
				text: iterates.map((_, i) => `x<sub>${i}</sub>`),
				textposition: 'top right',
				textfont: { size: 13, color: 'white' },
				name: 'iterates',
				showlegend: true
			});
		}

		plotFunctions()
		{
			const data = [];
			const layout = {
				title: "Newton's method — 2D minimization",
				autosize: true,
				height: 600,
				margin: { l: 50, r: 10, t: 50, b: 50 },
				xaxis: { title: 'x₁' },
				yaxis: { title: 'x₂', scaleanchor: 'x', scaleratio: 1 },
				legend: {
					x: 0.01, y: 0.99,
					xanchor: 'left', yanchor: 'top',
					bgcolor: 'rgba(50,50,50,0.85)',
					bordercolor: '#aaa', borderwidth: 1,
					font: { color: 'white' }
				}
			};

			if (this.fct === 'Elliptic quadratic')
			{
				layout.xaxis.range = [-6, 6];
				layout.yaxis.range = [-4, 4];
				this.computeData(ellipticQuadratic, this.x0, [-6, 6], [-4, 4], data);
				this.addMinima([[1.0, 0.5]], data);
			}
			else if (this.fct === 'Rosenbrock function')
			{
				layout.xaxis.range = [-6, 6];
				layout.yaxis.range = [-4, 4];
				this.computeData(rosenbrock, this.x0, [-6, 6], [-4, 4], data);
				this.addMinima([[1.0, 1.0]], data);
			}
			else if (this.fct === 'Himmelblau function')
			{
				layout.xaxis.range = [-6, 6];
				layout.yaxis.range = [-4, 4];
				this.computeData(himmelblau, this.x0, [-6, 6], [-4, 4], data);
				// four minima; only (3, 2) lies inside [-3,3]×[-3,3]
				this.addMinima([[3.0, 2.0], [-2.805, 3.131], [-3.779, -3.283], [3.584, -1.848]], data);
			}

			this.updateX0Display();
			Plotly.newPlot('plotOutput', data, layout).then(gd =>
			{
				gd.on('plotly_click', e =>
				{
					if (e.points && e.points.length > 0)
					{
						this.x0 = [e.points[0].x, e.points[0].y];
						this.plotFunctions();
					}
				});
			});
		}
	}

	plot = new Plot();
	plot.reset();
</script>

</body>
</html>