Alexander Godunov Personal website

Companion Website

Computational Physics

From Equations to Simulation

Chapter 13

Optimization and Data Fitting

These resources accompany Chapter 13 of Computational Physics: From Equations to Simulation by Alexander Godunov. They connect the practical codes with numerical libraries, provide simple benchmark problems, and help you choose methods and assess results in optimization and data fitting.

Chapter 13 codes Chapter 13 exercises & projects Chapter 13 updates & errata Chapter 13 overview

Software and library resources

The companion codes expose the numerical steps and provide a starting point for adapting a method to a new problem. For larger calculations, compare with a tested numerical library. The options below solve related tasks; they do not necessarily implement the same algorithm as the companion program. Read the stopping criteria and inspect the returned status as well as the parameter values.

On a narrow screen, scroll horizontally to view all language columns.

Numerical software by task and language
Numerical taskFortranC++MATLABPython
One-dimensional minimization The companion golden-section routine for a unimodal function on a finite interval Boost.Math brent_find_minima for a smooth function with one minimum in the chosen interval fminbnd for a local minimum in a finite interval SciPy minimize_scalar; use method='bounded' with finite bounds
Multivariable optimization NLopt has a Fortran interface for local and global search algorithms NLopt algorithms, including gradient-based and derivative-free methods fminsearch for unconstrained local search without derivatives SciPy minimize, differential_evolution, and dual_annealing; choose according to smoothness and bounds
Bounds and other constraints NLopt supports bounds; select an algorithm that also supports any required equality or inequality constraints NLopt with explicit bounds and an algorithm appropriate to the constraints fmincon (Optimization Toolbox) minimize with method='L-BFGS-B' and bounds for box constraints; use SLSQP or trust-constr for more general constraints
Linear least squares LAPACK DGELS for full-rank systems; DGELSY or DGELSD when rank deficiency is possible The same LAPACK solvers through a compatible C interface such as LAPACKE A\y for an overdetermined design matrix A numpy.linalg.lstsq; inspect the returned rank and singular values
Nonlinear least squares MINPACK LMDIF or LMDER for unconstrained fitting with a numerical or supplied Jacobian Ceres Solver for structured nonlinear least-squares problems, including parameter bounds lsqcurvefit (Optimization Toolbox) scipy.optimize.least_squares accepts a residual vector, with options for bounds and robust losses; curve_fit also returns an approximate parameter covariance

External libraries require installation and linking as appropriate. Their interfaces may require a newer compiler than the companion Fortran 95 programs. MATLAB fminbnd, fminsearch, and backslash are available in base MATLAB; lsqcurvefit and fmincon require Optimization Toolbox.

Benchmarks with known answers

Start with a problem whose answer is known before substituting a physical model or measured data. The examples below are deterministic checks, independent of random-number generators. For a noisy fitting problem, use statistical checks rather than expecting exact recovery of the parameters that generated one sample.

Simple checks for optimization and fitting routines
ProblemKnown answerWhat to check
f(x) = (x − 2)2 + 1 on [−1,4] x = 2, f = 1 The golden-section example. Vary the interval-width tolerance and compare the location error and function-call count.
f(x,y) = x2 + 4y2 + 2xy Unique minimum at (0,0), f = 0 The gradient-descent example. Its gradient is (2x + 2y, 2x + 8y); the Hessian has eigenvalues 5 − sqrt(13) and 5 + sqrt(13). A constant step converges from every starting point in exact arithmetic when 0 < eta < 2/(5 + sqrt(13)) ≈ 0.232408. The companion code uses backtracking instead.
f(x) = x4 − 14x2 + 24x + 5 on [−4,4] Global minimum: x = −3, f = −112. Other local minimum: x = 2, f = 13. Exercise 1. The derivative factors as 4(x − 1)(x − 2)(x + 3). The stationary point at x = 1 is a local maximum. Compare the basins reached from different starting points.
Rosenbrock function: f(x,y) = (1 − x)2 + 100(y − x2)2 Unique global minimum at (1,1), f = 0 Exercise 2. Compare distance from the known minimum and objective evaluations from several starting points. A narrow curved valley can slow local search.
Fit y = a0 + a1*x to x = [0,1,2,3], y = [1,3,5,7] a0 = 1, a1 = 2; zero residuals in exact arithmetic Check the coefficients and residuals. A zero residual estimate of noise for these constructed data is not evidence of zero uncertainty in real measurements.
Fit y = A*exp(−t/tau) to data generated without noise at t = [0,1,2,3,4], with A = 2, tau = 1.5 A = 2, tau = 1.5; zero residuals in exact arithmetic Check the parameter order, exponential sign, and fit in the original y values. Use several sensible initial guesses.

Choosing an optimization method

Choosing an optimization approach
ProblemUseful starting approachMain check
One variable, with one minimum in a known intervalGolden-section search or a bounded scalar minimizerCheck unimodality and the allowed domain. If the minimum over a closed interval is required, also evaluate its endpoints.
Smooth objective with reliable derivativesGradient descent with step control; Newton or quasi-Newton methods for faster local refinementScale the variables and check derivatives. A small gradient can also occur at a maximum or saddle point.
Local search when derivatives are inconvenientA derivative-free local method such as Nelder–MeadRepeat from different initial points and inspect sensitivity to parameter scaling.
Several basins in a bounded continuous search spaceSimulated annealing, particle swarm, or differential evolutionCompare several seeds at comparable evaluation budgets; refine promising candidates locally.
Parameters with physical bounds or other constraintsA solver that enforces the required bounds and constraintsCheck feasibility and the solver’s optimality measure. At an active bound, the unconstrained gradient need not vanish.
Routes, permutations, or other discrete candidatesAnt colony optimization or a genetic algorithm with a suitable representationEvery proposed route or candidate must remain valid. Include the return leg when evaluating a closed tour.

Data fitting: checks that matter

  • Use a suitable model. A model is linear in its parameters even when its basis functions are nonlinear in the independent variable. Polynomial fitting is a linear least-squares problem.
  • Match weights to uncertainties. For independent errors with known standard deviations sigma_i, minimize the sum of ((y_i − model_i)/sigma_i)2. Divide each row of the design matrix and the corresponding observation by sigma_i before a linear least-squares solve.
  • Inspect residuals. Plot y_i − model_i against the independent variable and, when useful, against observation order. Curvature, changing spread, or serial patterns can reveal assumptions that need attention.
  • Check parameter sensitivity. Try different initial guesses for nonlinear models. Strong parameter correlations or a nearly rank-deficient design matrix can make parameter estimates unstable even when the fitted curve looks good.
  • Interpret uncertainty correctly. Standard errors are not confidence intervals. For a correct linear model with independent Gaussian errors and known measurement standard deviations, reduced chi-squared has mean one and standard deviation sqrt(2/nu), where nu = m − k. For nu = 18, this is about 0.33, not a strict acceptance band. The result is approximate for nonlinear fits.
  • Compare predictive performance. Increasing polynomial degree cannot increase the minimum training sum of squared residuals for nested polynomial models in exact arithmetic. Use held-out data or cross-validation, physical plausibility, and parameter stability to judge whether the extra terms help.

Prefer QR or singular-value decomposition for general linear least-squares calculations rather than explicitly forming and inverting the normal-equation matrix. Centering and scaling the independent variable can improve polynomial conditioning.

Where parameter uncertainties come from

For m independent observations and k fitted parameters, let J be the Jacobian of the fitted residual vector, evaluated at the solution. In a linear model these formulas give the usual least-squares covariance; for a nonlinear model they give a local approximation. Assume m > k, full column rank, and an interior solution.

Covariance scaling for least-squares fits
Error informationJacobian and residual sumParameter covariance C
Known absolute standard deviations sigma_iUse residuals r_i = (y_i − model_i)/sigma_i and their weighted Jacobian J.C ≈ (JTJ)−1. Do not rescale by the fitted reduced chi-squared.
Unknown common standard deviationUse unweighted residuals and Jacobian. Let S = sum(r_i2) and s2 = S/(m − k).C ≈ s2(JTJ)−1.

The parameter standard errors are sqrt(C_jj). Compute covariance using a factorization rather than explicitly inverting the normal-equation matrix. Rank deficiency, active bounds, strong nonlinearity, or robust losses need additional treatment; the formulas above should not be applied mechanically in those cases.

In Python, curve_fit returns popt, pcov. Supply measured standard deviations as sigma and set absolute_sigma=True when their absolute values are known. The default absolute_sigma=False rescales the covariance using the fitted residual variance; with sigma=None, it assumes an unknown common noise scale. In contrast, least_squares returns residuals and a Jacobian, but no parameter covariance. For ordinary squared loss, its cost is half the sum of squared supplied residuals.

When measurement standard deviations are unavailable and independent errors with a common variance are a reasonable assumption, estimate that variance from the residual sum of squares divided by m − k. Count every parameter estimated from the data: fitting both slope and intercept uses k = 2; fitting only a slope with a known, fixed intercept uses k = 1. Report each parameter with its units, identify the uncertainty as a standard error or a confidence interval, and state how the noise variance was obtained. For a confidence interval, also give its confidence level and construction method.

Common pitfalls

  • Stopping is not proof of convergence. Record the objective value, stopping reason, feasibility, and objective-evaluation count. A completed search budget or a tiny step can indicate stagnation; neither establishes global optimality.
  • Iterations measure different amounts of work. Compare population methods and single-point searches using objective-evaluation budgets, together with achieved accuracy.
  • Transforming data changes the fit. Taking logarithms changes the fitting objective and error model. Fitting log-transformed observations is generally different from fitting an exponential in the original measurement scale.

Random seeds and bootstrap uncertainty

Keep a fixed seed while testing, debugging, and comparing changes. The stochastic companion programs provide a commented alternative for seeding from the clock. For a reliability study, run several recorded seeds and summarize the spread of results. The same seed need not produce the same random sequence across languages, generators, or Fortran compilers.

In a pairs bootstrap, sample observation pairs (x_i,y_i) together, with replacement, and refit every resample. Sampling x and y independently destroys their observed relationship. The standard deviation of the bootstrap estimates gives a bootstrap standard error; percentile bounds give a different summary and need not be symmetric about the original estimate. Record failed or degenerate fits and increase the number of resamples to check the stability of the reported uncertainty.

Choose the resampling scheme to match the experiment. Pairs resampling is natural for independent sampled pairs; fixed-design measurements may call for residual or parametric resampling under suitable error assumptions. Correlated measurements require a method that preserves the relevant dependence. Bootstrap resampling cannot repair systematic measurement errors or an unsuitable physical model. The SciPy bootstrap documentation provides a library comparison: for pairs use paired=True, and for percentile intervals select method='percentile' explicitly. Its default interval method differs from a percentile bootstrap.

Selected references and documentation

← Back to Chapter 13