Numerical Jacobian in Matlab

278 次查看（过去 30 天）
Daniel Wells
Daniel Wells 2012-2-5
回答： men8th 2024-11-20，15:07
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.)

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

采纳的回答

Anton Semechko
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
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) )
Daniel Wells
Daniel Wells 2012-2-5
I do not necessarily insist on using lsqnonlin to find a jacobian, I am only trying to find a (different, see above) way to do so in Matlab. I was referred to this method from the mathworks support site, and it seems that whoever answered that was wrong. You are right that 'MaxFunEvals' needs to be 0 as well, which she did not mention. I hadn't realized this, and it is indeed the answer to what I was trying to figure out.
Thanks.

请先登录，再进行评论。

更多回答（5 个）

Matthieu Heitz
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.
  1 个评论
Victor
Victor 2024-9-9
编辑：Victor 2024-9-9
it seems that this solution does not work properly.
Let's try the point x = [1; 2] and the Jacobian is not correct.
The code produces J = [2 2; 24 24] instead of [2 4; 24 12]

请先登录，再进行评论。

John D'Errico
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.
http://www.mathworks.com/matlabcentral/fileexchange/13490-adaptive-robust-numerical-differentiation
  1 个评论
Daniel Wells
Daniel Wells 2012-2-5
I really do like your jacobianest a lot, by the way. I use it often, and have never had a problem with it, and don't expect the problem to be there in this case. One just has to be careful with these things.

请先登录，再进行评论。

Daniel Wells
Daniel Wells 2012-2-5
I am using jacobianest already, and was looking for a second opinion on the estimation for the jacobian of a particularly nasty function that I am using, since my code is behaving poorly, which makes me have to question every function I am using already.
  2 个评论
John D'Errico
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
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
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
  1 个评论
Daniel Wells
Daniel Wells 2012-2-5
There is something wrong, because I do not want to Jacobian evaluated at [.7308, .7308], but rather at [1, 1]. It is giving me the Jacobian at an x not of my choosing.

请先登录，再进行评论。

men8th
men8th 2024-11-20，15:07
Ran in:
Matlab has a built-in, but not very well documented function which will compute Jacobians numerically
help numjac
NUMJAC Numerically compute the Jacobian dF/dY of function F(T,Y). [DFDY,FAC] = NUMJAC(F,T,Y,FTY,THRESH,FAC,VECTORIZED) numerically computes the Jacobian of function F(T,Y), returning it as full matrix DFDY. F is a function handle. For a scalar T and a column vector Y, F(T,Y) must return a column vector. T is the independent variable and Y contains the dependent variables. Vector FTY is F evaluated at (T,Y). Column vector THRESH provides a threshold of significance for Y, i.e. the exact value of a component Y(i) with abs(Y(i)) < THRESH(i) is not important. All components of THRESH must be positive. THRESH can also be specified as a cell array, where THRESH{1} provides a threshold of significance for Y, and THRESH{2} provides a typical value of Y. Column FAC is working storage. On the first call, set FAC to []. Do not alter the returned value between calls. VECTORIZED tells NUMJAC whether multiple values of F can be obtained with a single function evaluation. In particular, VECTORIZED=1 indicates that F(t,[y1 y2 ...]) returns [F(t,y1) F(t,y2) ...], and VECTORIZED=2 that F([x1 x2 ...],[y1 y2 ...]) returns [F(x1,y1) F(x2,y2) ...]. When solving ODEs, use ODESET to set the ODE solver 'Vectorized' property to 'on' if the ODE function has been coded so that F(t,[y1 y2 ...]) returns [F(t,y1) F(t,y2) ...]. When solving BVPs, use BVPSET to set the BVP solver 'Vectorized' property to 'on' if the ODE function has been coded so that F([x1 x2 ...],[y1 y2 ...]) returns [F(x1,y1) F(x2,y2) ...]. Vectorizing the function F may speed up the computation of DFDY. [DFDY,FAC,G] = NUMJAC(F,T,Y,FTY,THRESH,FAC,VECTORIZED,S,G) numerically computes a sparse Jacobian matrix DFDY. S is a non-empty sparse matrix of 0's and 1's. A value of 0 for S(i,j) means that component i of the function F(T,Y) does not depend on component j of vector Y (hence DFDY(i,j) = 0). Column vector G is working storage. On the first call, set G to []. Do not alter the returned value between calls. [DFDY,FAC,G,NFEVALS,NFCALLS] = NUMJAC(...) returns the number of values F(T,Y) computed while forming dFdy (NFEVALS) and the number of calls to the function F (NFCALLS). If F is not vectorized, NFCALLS equals NFEVALS. Although NUMJAC was developed specifically for the approximation of partial derivatives when integrating a system of ODE's, it can be used for other applications. In particular, when the length of the vector returned by F(T,Y) is different from the length of Y, DFDY is rectangular. See also COLGROUP, ODE15S, ODE23S, ODE23T, ODE23TB, ODESET. Other uses of numjac customreg/numjac idnlfun/numjac

类别

在 Help Center 和 File Exchange 中查找有关 Stability Analysis 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by