# Removing for loop in computing the finite difference Jacobian

16 views (last 30 days)
Hassan on 9 Jan 2018
Commented: Matt J on 9 Jan 2018
Hello friends,
I am trying to remove the for loop i.e. to vectorize the following code that computes the numerical Jacobian using finite difference. I think the vectorized version of the code will be faster than the current version especially for computing the Jacobian of a large-scale problem.
In addition, I will like to multiply J^T with v, where v is any vector with dimension equal to the rows of J, how can I include this in my code? Please, can anybody help me out?
if true
function [J]=numjac(f,x)
% computes the Jacobian of a function
tic
n=length(x);
fx=feval(f,x);
eps=1.e-8;
xperturb=x;
for i=1:n
xperturb(i)=xperturb(i)+eps;
J(:,i)=(feval(f,xperturb)-fx)/eps;
xperturb(i)=x(i);
end
toc
end

Matt J on 9 Jan 2018
Edited: Matt J on 9 Jan 2018
I don't think you can vectorize further (for general functions f), but not pre-allocating J will definitely hurt you in large problems,
J=nan(length(fx),length(x)); %pre-allocate
for i=1:n
xperturb(i)=xperturb(i)+eps;
J(:,i)=(feval(f,xperturb)-fx)/eps;
xperturb(i)=x(i);
end
You could also use a PARFOR loop if you have the Parallel Computing Toolbox. Furthermore, the Optimization Toolbox least squares solvers have options to do all of the above for you internally.
Matt J on 9 Jan 2018
The efficient way to do that is always problem-specific. That's why Optimization Toolbox solvers like lsqnonlin give you the option of supplying a customized JacobianMultiplyFcn.