# Removing for loop in computing the finite difference Jacobian

22 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.

Show 1 older comment
Matt J on 9 Jan 2018
Correct.
Hassan on 9 Jan 2018
Alright, suppose I need to compute J^T*v, where v is any vector with dimension equal to the rows of J. Could you please modify the code so that J^T*v is computed instead of J only? In this case, the output will be a vector of dimension same as the columns of J. How is the vector w=J^T*v be assembled using the for loop? or any other MATLAB function?
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.