Removing for loop in computing the finite difference Jacobian

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

Answer 1

[Matt J]

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.

[Hassan]

Assuming f: R^n to R^m, where m>=n. In this case, f is general, so you mean vectorizing the code further will not be possible?

[Matt J]

Correct.

[Hassan]

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]

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.