Removing for loop in computing the finite difference Jacobian

2 次查看（过去 30 天）
Hassan
Hassan 2018-1-9
评论： Matt J 2018-1-9
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

请先登录，再回答此问题。

回答（1 个）

Matt J
Matt J 2018-1-9
编辑：Matt J 2018-1-9
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 个评论
显示 2更早的评论
Hassan
Hassan 2018-1-9
编辑：Hassan 2018-1-9
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 2018-1-9
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by