Removing for loop in computing the finite difference Jacobian

9 views (last 30 days)
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

Answers (1)

Matt J
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.
  4 Comments
Hassan
Hassan on 9 Jan 2018
Edited: 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
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.

Sign in to comment.

Categories

Find more on Performance and Memory in Help Center and File Exchange

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!