Numerically Computing the Gradient of a Function in MATLAB

Question

I have a question on using Matlab's gradient function. Here is sample code:

npts=100;
x1 = linspace(-10,10,npts);
x2 = linspace(-10,10,npts);
x3 = linspace(-10,10,npts);
f1 = x1.^2 + 2*x2.^2 + 2*x3.^3 + 2*x1.*x2 + 2*x2.*x3;
f2 = @(x1,x2,x3) x1^2 + 2*x2^2 + 2*x3^3 + 2*x1*x2 + 2*x2*x3;

I want to numerically compute the partial derivatives of function $f_2$:

$$ dx(i)=\frac{\partial f_2}{\partial x} \bigg |_{(x1(i),x2(i),x3(i))}\\ dy(i)=\frac{\partial f_2}{\partial y} \bigg |_{(x1(i),x2(i),x3(i))} \\ dz(i)=\frac{\partial f_2}{\partial z} \bigg |_{(x1(i),x2(i),x3(i))} $$

where $i$ goes from 1 to npts.

I also defined f1, but even after reading documentation, I could not figure out how to find gradient at given points. Any help is appreciated.

Answer 1

You can use Finite Differences method for that.

Have a look at How does one compute a single finite differences in Matlab efficiently?.
You'll find in my answer an efficient code to do so.

Answer 2

A way better approach (numerically more stable, no issue of choosing the perturbation hyperparameter, accurate up to machine precision) is to use algorithmic/automatic differentiation. For this you need the Matlab Deep Learning Toolbox. Then you can use dlgradient to compute the gradient. Below you find the source code attached corresponding to your example.

Most importantly, you can examine the error and observe that the deviation of the automatic approach from the analytical solution is indeed machine precision, while for the finite difference approach (I choose second order central differences) the error is orders of magnitude higher. For 100 points and a range of $[-10, 10]$ this errors are somewhat tolerable, but if you play a bit with Rand_Max and n_points you observe that the errors become larger and larger.

Error of algorithmic / automatic diff. is: 1.4755528111219851e-14
Error of finite difference diff. is: 1.9999999999348703e-01 for perturbation 1.0000000000000001e-01
Error of finite difference diff. is: 1.9999999632850161e-03 for perturbation 1.0000000000000000e-02
Error of finite difference diff. is: 1.9999905867860374e-05 for perturbation 1.0000000000000000e-03
Error of finite difference diff. is: 1.9664569947425062e-07 for perturbation 1.0000000000000000e-04
Error of finite difference diff. is: 1.0537897883625319e-07 for perturbation 1.0000000000000001e-05
Error of finite difference diff. is: 1.5469326944467290e-06 for perturbation 9.9999999999999995e-07
Error of finite difference diff. is: 1.3322061696937969e-05 for perturbation 9.9999999999999995e-08
Error of finite difference diff. is: 1.7059535957436630e-04 for perturbation 1.0000000000000000e-08
Error of finite difference diff. is: 4.9702408787320664e-04 for perturbation 1.0000000000000001e-09

Source Code:

f2.m

function y = f2(x)

  x1 = x(:, 1);
  x2 = x(:, 2);
  x3 = x(:, 3);
  
  y = x1.^2 + 2*x2.^2 + 2*x3.^3 + 2*x1.*x2 + 2*x2.*x3;

f2_grad_analytic.m:

function grad = f2_grad_analytic(x)
  x1   = x(:, 1);
  x2   = x(:, 2);
  x3   = x(:, 3);
  
  grad(:, 1) = 2*x1 + 2*x2;
  grad(:, 2) = 4*x2 + 2*x1 + 2 * x3;
  grad(:, 3) = 6*x3.^2 + 2*x2;

f2_grad_AD.m:

function grad = f2_grad_AD(x)
  x1 = x(:, 1);
  x2 = x(:, 2);
  x3 = x(:, 3);
  
  y = x1.^2 + 2*x2.^2 + 2*x3.^3 + 2*x1.*x2 + 2*x2.*x3;
  grad = dlgradient(y, x);

CalcNumericalGradient.m:

function NumericalGrad = CalcNumericalGradient(InputPoints, eps)

% (Central, second order accurate FD)
NumericalGrad = zeros(size(InputPoints) );
for i = 1:size(InputPoints, 2)
  perturb             = zeros(size(InputPoints));
  perturb(:, i)       = eps;
  NumericalGrad(:, i) = (f2(InputPoints + perturb) - f2(InputPoints - perturb)) / (2 * eps);
end

main.m:

clear;
close all;
clc;

n_points = 100;
Rand_Max = 20;
x_test_FD = rand(n_points, 3) * Rand_Max - Rand_Max/2;

% Calculate analytical solution
grad_analytic = f2_grad_analytic(x_test_FD);

grad_AD = zeros(n_points, 3);
for i = 1:n_points
  x_test_dl    = dlarray(x_test_FD(i,:) );
  grad_AD(i,:) = dlfeval(@f2_grad_AD, x_test_dl);
end
Err_AD = norm(grad_AD - grad_analytic);
fprintf("Error of algorithmic / automatic diff. is: %.16e\n", Err_AD);

eps_range = [1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8, 1e-9];
for i = 1:length(eps_range)
  eps = eps_range(i);
  grad_FD = CalcNumericalGradient(x_test_FD, eps);
  Err_FD = norm(grad_FD - grad_analytic);
  fprintf("Error of finite difference diff. is: %.16e for perturbation %.16e\n", Err_FD, eps);
end