Numerical Jacobian in Matlab
조회 수: 152 (최근 30일)
이전 댓글 표시
I am trying to estimate the Jacobian of a multivariable function in Matlab. I came across a method here:
<http://www.mathworks.com/support/solutions/en/data/1-CMM53T/index.html?product=OP&solution=1-CMM53T>
that recommends using the function lsqnonlin. As an example, I define the function
f = @(x)[x(1)^2 + x(2)^2, x(1)^3.*x(2)^3];
and the point
A = [1 1];
The Jacobian of f at this point should pretty clearly be [2, 2;3, 3];
However, when I use the method described in the above link, I find that
[x, ~, ~, ~, ~, ~, jac] = lsqnonlin(f, A,[],[], optimset('MaxIter', 0))
x = [.7308 .7308]
and jac = [1.4615, 1.4615; .6252, .6252];
Clearly, something is a foot. Any ideas what is going on? Any ideas of a better way, using a native MATLAB function, to evaluate a jacobian numerically? (I don't want to do it symbolically.)
댓글 수: 0
채택된 답변
Anton Semechko
2012년 2월 5일
lsqnonlin is not designed for the purpose of evaluating Jacobians. I thought your main concern was that it does not evaluate the Jacobian accurately. But given the above example, it does.
댓글 수: 2
Anton Semechko
2012년 2월 5일
but if you insist on using lsqnonlin for the sole purpose of evaluating Jacobians then instead of MaxIter=0 you have to set MaxFunEvals=0 (i.e. optimset('MaxFunEvals', 0) )
추가 답변 (4개)
Matthieu Heitz
2017년 2월 2일
Hi !
If I may provide another answer that uses finite differences.
% Jacobian functor
J = @(x,h,F)(F(repmat(x,size(x'))+diag(h))-F(repmat(x,size(x'))))./h';
% Your function
f = @(x)[x(1)^2 + x(2)^2; x(1)^3.*x(2)^3];
% Point at which to estimate it
x = [1;1];
% Step to take on each dimension (has to be small enough for precision)
h = 1e-5*ones(size(x));
% Compute the jacobian
J(x,h,f)
This J function works with any function f. Note that x and the result of f must be in a vertical vector for it to work.
댓글 수: 0
John D'Errico
2012년 2월 5일
Whats the problem?
f = @(x)[x(1)^2 + x(2)^2, x(1)^3.*x(2)^3];
jacobianest(f,[1 1])
ans =
2 2
3 3
It is on the file exchange.
Daniel Wells
2012년 2월 5일
댓글 수: 2
John D'Errico
2012년 2월 5일
You can always use derivest to estimate individual components of the Jacobian. It will allow you much more control over the estimation procedure, with many more options. You can vary the order of the approximation, apply one sided or two sided estimators, change limits on the steps taken, etc. If these estimates are consistent with the Jacobian, taking into account the error estimates it provides, then you can assume it is doing its job reasonably well. Jacobianest is not as sophisticated as is derivest, with far fewer options.
Of course, one can always come up with a nasty enough function that any adaptive scheme will go wrong. For proof of that, just look at the number of people who post questions asking why quad/quadgk etc. has failed on their function. And those people are not even actively trying to make those tools fail.
Hannah Trachtman
2019년 6월 19일
Hi John,
I realize this is many years later, but I am trying to do exactly this and having trouble. I'm working from your example code below:
fun = @(x,y) x.^2 + y.^2;
xy = [2 3];
gradvec = [derivest(@(x) fun(x,xy(2)),xy(1),'d',1), ...
derivest(@(y) fun(xy(1),y),xy(2),'d',1)]
But my function is a large vector-valued function stored in an m-file (since I was originally using jacobianest). It looks something like this:
function [calcMoments] = calculator(parameters)
alpha = parameters(1)
beta = parameters(2)
calcMoments = alpha^2 + beta^2
end
How do I modify your code to apply to this kind of function? I'm a relatively inexperienced Matlab user so I might be missing something trivial. Can you point me in the right direction?
Thank you!
Hannah
Anton Semechko
2012년 2월 5일
Actually there is non problem what so ever. The Jacobian of your function is:
J=@(x) [2*x(1) 2*x(2);3*x(2).*(x(1).*x(2)).^2 3*x(1).*(x(1).*x(2)).^2];
Evaluating J at x = [.7308 .7308]:
jac=J([.7308 .7308])
jac =
1.4616 1.4616
0.6253 0.6253
you get the same answer as that returned by lsqnonlin
참고 항목
카테고리
Help Center 및 File Exchange에서 Function Creation에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!