一个很好的例行程序，但要注意:
你好,
虽然我发现这对很多矩阵都是一个好方法，但当雅可比矩阵的特征值接近0时，我发现它失败了几个数量级。我比较了这个求解器的中心差分近似的雅可比矩阵和分析表达式的特征值在这个极限，并发现中心差分工作得比这个函数更好。
尽管如此，在大多数情况下，我发现它非常有用，比中心差分更准确。
欢呼,
伊恩

你好，约翰，

在花了一些时间在代码中并将其与正确工作的派生函数进行比较之后。我认为问题出在矩阵的转置上。更具体地说，应该使用。'操作符(这是厄米操作符)而不是。'操作符(这是转置操作符)。
针对上述问题，我提出以下三点建议:
在雅可比矩阵中。m文件:
1.第189行:change from: rombcoefs = rromb\(qromb'*rhs);To rombcoefs = rromb\(qromb.'*rhs);
2.第190行:change from: der_rombs = rombcoefs(1，:)';To der_rombs = rombcoefs(1，:).';
3.第196行:change from: errest = s'*12.7062047361747*sqrt(cov1(1));'*12.7062047361747*sqrt(cov1(1));

请自行检查以上更改。再次感谢你抽出时间。

你好，约翰，我也想感谢你这个很酷很有价值的工具。

当我在使用它的时候，我注意到一个问题，我猜这个问题与以前liuxbin xbin发现的问题非常相似。它发生在试图计算一个有复数输入的向量函数的雅可比矩阵时，其中得到的结果是预期结果的共轭。复制的小代码如下:

Myfun = @(x) [exp(x(1))，exp(x(2))]
值= [1+1i,1+1i]
Res1 = myfun(value)
[res2, err] = jacobianest(@(x) myfun(x)，value)

如果你能看一眼就太酷了。先谢谢你

John，作为一项愿望，我想知道您是否考虑扩展代码以处理返回gpuArray输出的函数。例如，这将会失败，

> > c = gpuArray (2);
>> fun=@(x) c*x;
> >雅可比矩阵(有趣,1)

当然，你可以坚持让用户定义函数，以便结果总是被拉回主机:

>> fun=@(x) gather(c*x);

但是在高维的例子中，这会带来很大的开销。

伟大的工作!对我的研究非常有价值!谢谢楼主分享!

@Tim -尝试使用任何类型的fminsearch，有15个未知数会带来不切实际的期望。15个未知数对于fminsearch来说处理得太好了。Fminsearchcon也好不到哪里去，因为它只是fminsearch的一个覆盖。这是相同的基本算法。

更糟糕的是，目标函数似乎有一个奇点。这表明您的模型参数化过度。至少有一个未知数看起来是不必要的至少在局部是这样。(一般来说，奇点是指参数的某种线性组合允许将一组参数改为另一组参数，而不会对目标函数造成任何损失。)这可能意味着您的模型过于复杂，或者您的数据不足以适应这种复杂性的模型。

你得到了一个非常好的契合，这并没有太大的关系。换一组参数也能得到同样好的结果。反过来，这意味着试图计算任何类型的置信区间都是错误的。这些置信区间并不能告诉你什么，除了我上面所解释的。

这是一个很棒的包装!非常感谢!
我使用雅可比函数来计算雅可比矩阵和我的拟合参数的置信区间，这是基于John D’errico在这里的评论:https://groups.google.com/forum/#!topic/comp.soft-sys.matlab/JJOqLCNJjF0
我的问题是雅可比矩阵的一列有非常小的值(1e-29)，当我计算与inv(J'*J)成比例的西格玛(参见上面的链接)时，我有一个消息:警告:矩阵接近奇异或糟糕
按比例缩小的。结果可能不准确。RCOND =
5.516957 e - 117。

然而，我有一个很好的拟合，拟合参数有很多意义，但他们的误差是如此之大，这是完全不现实的。我知道上面的方法假设我有一个高斯误差系统，但老实说，我不知道我怎么能确定它。我有一个15个参数的非线性方程，我正在使用fminsearchcon找到最小二乘法的最佳拟合。
当高斯假设不成立时，是否有一种方法可以使用现有的Matlab函数来获得置信区间?
欢迎提出任何建议。

非常感谢这个不错的工具和文档。节省了很多时间，这正是我要找的…甚至更好!

干得好，先生

我使用黑森矩阵来计算极大似然估计参数的置信区间估计。而黑森逆矩阵的对角线为负。我该如何处理这个问题?

这段代码很棒。但应该谨慎，因为它似乎有时会给出错误的结果。其他评论者也提到了这一点。我遇到的情况是在计算二维朗格曼函数的梯度时。与解析梯度相比，范围[0,10]*[0,10]中的梯度似乎是不正确的。

dervestsuite的C/ c++端口:
https://github.com/ofrei/cppDERIVEST

这个软件真的很棒!
值得注意的是，有一种情况下DERIVESTsuite似乎给出了错误的答案:在“向后”模式下计算的2阶导数和4阶导数报告的符号相反。

>>衍生est(@exp, 1， ' derivatives ativeorder '， 2， 'Style'， 'backward')
ans =
-2.7183
>>衍生est(@exp, 1， ' derivatives ativeorder '， 4， 'Style'， 'backward')
ans =
-2.7183

@Yaguang杨
它肯定比fminunc好，但我要提醒你不要假设“代码从来没有失败过”。那里仍然有皱纹，就像保罗(2015年11月)报道的这个

>>雅可比矩阵(@(x) x+ 10,1e -12)

ans =

0

非常可靠的工具。我基于John的代码构建了我的修改牛顿代码，并在matlab优化工具箱中将修改后的牛顿与fminunc进行了比较。对于随机生成的问题(一个工程师的问题会有不同的大小和参数)的数百万次试验，代码从未失败过;但是fminunc失败了很多次。包括约翰的DERIVESTsuite在内的改进牛顿的问题是速度比fminunc慢得多。我要问约翰的问题是:你能不能为你的matlab代码构建一个墨西哥文件，这样人们就可以用墨西哥文件而不是m文件来解决大问题了?

@Abhinav,
即使你计算的是正x处的导数，雅可比矩阵所使用的有限差分也没有理由先验地避免在负x处计算函数。有限差分总是涉及到在所期望的点的左右一段距离的采样。但是，如果您正在处理标量x(听起来像是在示例中)，则可以使用衍生方法。derived提供了限制有限差分步长的输入选项，覆盖了自适应选择。

@matt
非常感谢您的回复。我不打算计算x=0处的导数。问题是“雅可比矩阵”仍然试图在负值处求函数的值。例如，当我试图计算x=2.5处的导数时，雅可比矩阵是求函数在252.50和-247.50处的值。

@Abhinav

对于一个函数在一点上是可微的，它必须在该点周围的开放区域上定义。这是导数定义的基本部分。所以f(x)只在x>=0时定义在x=0时是不可微的。

嗨,约翰,
谢谢你的代码!
我想求出一个有正支持度的函数的导数，也就是说，这个函数只对正参数有定义。金宝app我还没能想出对变量施加约束的方法。你能建议一种方法吗?

最好的
AG)。

在parfor循环中运行派生est有什么问题吗?我已经成功地使用了一段时间，今天我将它嵌入到parfor循环中，我得到了一个错误，我无法单独重现。

一个名为MaxPol的通用框架最近已经发布，并在这里提供
//www.tatmou.com/matlabcentral/fileexchange/63294-maxpol-smoothing-and-differentiation-package

MaxPol提供了一个框架来设计各种数值微分内核，其属性如下:
(1)无侧高程伪影的截止(低通)设计(用于抗噪声情况)
(2)任意的微分顺序
(3)任意多项式精度
(4)导数矩阵设计
(5)带转向力矩的二维导数核
(6)信号与图像处理中的直观例子

使用此代码进行MLE计算，我发现它快速、无缝且构思良好。毫无疑问……John D'Errico生产的代码简单、创新、防弹。我感谢他的贡献。
马克

这是一个很好的工具箱，但当我使用jacobianest/hessianest函数时，可能会出现一些严重的问题。例如，当我在脚本文件fun1.m中定义非常简单的函数fun1时:

函数res = fun1 (x)
如果x < = 20
错误('x必须大于20');
结束
res = sin (x);
结束
那么fun1可以在点21处计算，值为0.8367。然而，当使用jacobianest函数时，如下代码
雅可比矩阵(@sin, 21);
我发现以下错误通知:
使用fun1出错(第3行)
20 !

雅可比矩阵误差(第112行)
Fun (swapelement(x0,i,x0_i - delta(j)));

我不知道为什么。期待您的回复。
致以最亲切的问候
Pingjun邓

关闭所有
清除所有
clc

n = 100;

m1 = 1, m2 = 2;
s1 = 1; s2 = 1;

x1 = normrnd (m1, s1, n, 1);
x2 = normrnd (m2、s2, n, 1);

h1 = 0.9 *分钟(std (x1)、差(x1) / 1.34) *长度(x1) ^ (-0.2);
h2 = 0.9 *分钟(std (x2)、差(x2) / 1.34) *长度(x2) ^ (-0.2);

F1inv = @(在)ksdensity (x1,,“功能”,“icdf”,“宽度”,h1); %“支持”,“积极的”)金宝app;
F2 = @(在)ksdensity (x2,“功能”,“提供”,“宽度”,h2); %“支持”,“积极的”);金宝app

装置= @ (t) 1-F2 (F1inv (1 - t));

t = 0:0.01:1;
情节(t,装置(t));
title('我想要这条曲线在任意点(0

%=============================
我想求rochat的导数at=0.8 0.6 0.4 0.2 0.1
myderivatives 1= derivatives (@(x) rochat(x)， 0.8)
myderivatives 1= derivatives (@(x) rochat(x)， 0.6)
myderivatives 1= derivatives (@(x) rochat(x)， 0.4)
myderivatives 1= derivatives (@(x) rochat(x)， 0.2)
myderivatives 1= derivatives (@(x) rochat(x)， 0.1)

为什么前三个是0?

神奇的工具。伟大的发现错别字的梯度和雅可比矩阵，是馈送给求解器和优化器。

@Peng, derivatives计算一个变量的函数的导数，即另一个变量的函数。黑森矩阵是一个n维变量函数的二阶导数的方阵。

你好,
我在使用代码时遇到了一些问题。我写:
d =派生(@ (x) exp ([1 3] * x - 0.1) + exp ([1 3] * x - 0.1) + exp ([1 0] * x - 0.1), (0, 0))
，但它返回一个错误:
使用*错误
矩阵内部维数必须一致。

错误GeneralFunction > @ (x) exp ([1,3] * x - 0.1) + exp ([1,3] * x - 0.1) + exp ((1,0) * x - 0.1)

派生函数中的错误(第362行)
F_plusdel = fun(x0i+h*delta);
但是我们在MATLAB中使用命令hessian:
h =黑森(@ (x) exp ([1 3] * x - 0.1) + exp ([1 3] * x - 0.1) + exp ([1 0] * x - 0.1), (0, 0))
它返回正确答案。
你能帮我看看导数有什么问题吗?

嗨,约翰,

谢谢分享这个伟大的包!:)我有一个小问题:

我可以在雅可比矩阵中选择方法(向前等)和精度吗?
只有(?)雅可比矩阵允许向量化代码，但没有提供选项。

谢谢马特·J。

我现在有其他的问题，虽然是hessdiag。M，它进入这个for循环

对于ind = 1:nx
(高清(印第安纳州),犯错(印第安纳州),finaldelta(印第安纳州)]=推导(…
@(xi) fun(swapelement(x0,ind,xi))，…
x0(印第安纳州),“导数”,2、“矢量化”,“不”);
结束

它传递了x0(ind)所以它传递了x0向量的第一个元素去求导。然后在derived中:

我们需要乐趣吗?
if (rem(par. derivatives ativeorder,2) == 0) || ~strncmpi(par.Style，'central'，7)
如果strcmpi (par.Vectorized,“是的”)
F_x0 = fun(x0);
其他的
%未向量化，因此循环
F_x0 = 0 (size(x0));
对于j = 1:numel(x0)
F_x0 (j) = fun(x0(j));
结束
结束
其他的
F_x0 = [];
结束

失败:
f_x0(j) = fun(x0(j))
在赋值A(I) = B中，B和I中的元素数量必须相同。

因为我们传递了x0(ind)它是一个单独的数而不是整个向量。什么好主意吗?我不想在这个doe周围混日子，并试图解决它，因为我认为我在第一时间做错了什么。谢谢

@Huck,

这不需要对derivative est进行任何修改。您可以使用这里描述的技术将任何带有额外形参的函数转换为单参数函数，

//www.tatmou.com/help/optim/ug/passing-extra-parameters.html

你能修改这个工具箱，使约束函数也能接受一组参数吗?例如，如果我试图计算约束的雅可比因约束函数是这样的:

C = constraints(x, auxdata);

其中auxdata是我的约束函数需要运行的参数单元格数组

谢谢!

它显然比优化工具箱中的有限差分要健壮得多。然而，我也看到了Paul(2015年11月23日)报道的失败案例。一旦这些皱纹被消除，我会把我的评分提高到5颗星。

嗨,约翰,
感谢分享此代码。我发现了一些数值问题，当估计雅可比矩阵时，如果起点接近零。

例如:
雅可比矩阵(@(x) x+ 10,1e -12)
返回0，尽管它应该是1。起始点为1e-11就可以了。Zero也能很好地工作，因为它在代码中是作为特殊情况处理的。问题是将非常小的delta与目标函数中相对较大的+10组合在一起，因为如果起始点为~= 0 (jacobianest.m中的第100行)，则delta会按起始点缩放。一种变通方法是不按起点缩放。

约翰,

正是因为我知道差分公式是非常病态的，所以我正在测试一个奇特的有限差分格式什么时候会失效。我循环不同的coeff值来得到这个感觉，每次我碰到提到的coeff值，我看到导数在向量v的每一点都很好，除了v==4附近。换句话说，你对我信念的假设是错误的，问题仍然存在:当噪声均匀随机分布在区间[-1e-10,1e-10]时，v==4有什么特别之处?
如果您使用这里可用的派生版本多次运行我的代码，您应该明白我的意思。

不，不，不。没有错误，只是您认为衍生工具是用来解决您提出的问题的。

一个包含随机噪声的函数根据定义是处处不可微的。导数不是设计来微分一个有噪声的函数的，也不应该期望它这样做。当导数函数试图对微小的噪声(本质上是最不重要的位)不敏感时，你就不能使用一个试图微分一个实际上有导数的函数的工具了。

虽然你的问题中的基本函数有一个导数，但一旦你给它加上噪声，结果是不可微的。

注意，微分是一个不适定问题，从某种意义上说，它倾向于放大函数中的任何噪声。我不确定这是否重要，但是你尝试用更高阶的方法来产生导数。这将导致更多的噪音放大。

您可能希望使用一种工具，从本质上消除噪音，区分结果。有很多这样的工具，这取决于你的应用程序和你最终需要什么。这可以是任何工具，从Savitsky-Golay工具到最小二乘样条或局部多项式模型，或其他工具。

约翰,

我用噪声测试有限差分，如下所示:
多项式系数= 1平台以及;
v = 0:0.1:2 *π;
(d ~) =派生(@ (x) sin (x) +多项式系数*(1 + 2 *兰德(大小(x))), v,“导数”,2,‘MethodOrder’,4);
情节(v d ' - o ')

每次我运行这段代码(即在区间[-coeff coeff]中使用rand生成不同的噪声)，在v==4附近的误差要大得多。这似乎不寻常，我认为这是一个bug。我没查过衍生案，但我想我会举报的。

感谢您的投稿!!

罗德里戈。

你好约翰,

伟大的工作。很有帮助。

我认为hessian计算可以显著加快，如果被测函数的计算很慢(就像我的情况一样):
在黑森。M你叫hessiandiag。M。它需要计算特定样本点的函数值。
年级。M，也被称为hessian。M，计算与hessiandiag.m相同样本点上的函数值。所以没有必要做两次。

@Simon Tardivel:我可以证实你的观察和你的解决方案。谢谢你！

你好约翰,

非常感谢您提供的代码。它写得很好，使用简单，足够快。奇怪的是，这些函数在官方内置函数中仍然缺失。

但是，我可能发现了一个错误/意外行为。不幸的是，我不能分享代码(加上它需要数千行代码…)，但这里有一个问题和一个工作修复:

情况:我正在计算的hessdiag应该产生[a b c]，其中a,b,c都不等于0。我知道它们都不等于0，这是物理定律的要求。然而，b相当小，通常比a小几个数量级。a和c具有可比性，符号相反。

问题:hessdiag输出[a,0,0]，估计误差为[10^-7*a, 0,0]。

修正:我去衍生和改变默认RombergTerms从2到3。B和c又出现了。B估计为10^-3*a, c约为-a，这是它们的正确值。

我不知道根本的问题是什么。我刚刚注意到hessian有默认的Romberg terms = 3，但没有将这个参数传递给假设Romberg terms = 2的梯度计算。事实上，Romberg terms = 3似乎只用于上三角形(实际上是3+1)。也许这就是问题所在?我对数值微分方法一无所知，所以这可能完全偏离了主题……

最好的
西蒙。

你好,

我想知道你能否就以下问题提出你的看法。
对于我的函数，雅可比矩阵不成立我想原因可能是我试图在函数域边界附近求雅可比矩阵的值。这是一个合理的猜测吗?如果是这样，我该如何解决这个问题?我的直觉是MaxStep太大了，代码试图在没有定义的地方计算函数。但是我可以选择一个更小的MaxStep值吗?

谢谢!

你知道你输入的函数实际上不是有效的MATLAB语法吗?

f=@(x)(1/2*x'(Ax) + x'b)

例如，当你在纸上写x'b时，它可能是可以的，但MATLAB需要一个*运算符来理解乘法(点积)。

接下来，你这样写:

f=@(x)(1/2*x'(Ax) + x'b)
获得(f0, x0),

定义函数f，然后用f0调用导数函数。它怎么知道你的意图不是你所写的呢?

我猜你做事不小心。使用有效的语法。

我有一个二次函数
1/2x'Ax + x'b
但它似乎不起作用，它一般起作用吗?还是我说错了?

我试过:
f=@(x)(1/2*x'(Ax) + x'b)
派生(f0,x0)，但总是得到以下错误:
输入参数不够。

派生函数中的错误(第362行)
F_plusdel = fun(x0i+h*delta);

有什么办法能让它正常工作吗?
谢谢!

你可能认为多元函数的梯度是导数，但它不是。导数至少是在一个变量的函数上定义的，除非他们在我学这门课之后改变了微积分。

导数计算函数y(x)的导数，dy/dx，一个单变量的函数。

事实上，我确实在包中提供了一个梯度计算器，所以您可能已经尝试过了!看看帮助。看看演示。

最后，我不会那么仔细地查看代码是否有意义。但快速扫一眼就能看到这样的线条:

P = max(P,0.0000001);

这将使你的函数不可微，并且根本没有任何目的，因为你首先要计算的是一个对数似然函数!当然，它看起来并不相关，因为在代码中没有在后面使用p。如果你想把它转储到位桶中，为什么要计算它呢?

下面这几行是有用的，它们确实会使函数不可微。

Y =重塑(Y,max(alt)，max(ncs));
X =重塑(exp(v(b))，max(alt)，max(ncs));

所以我认为你需要考虑你在做什么，如果有一个优化器失败的原因，如果它假设它要优化的函数是可微的，当你的函数不是的时候。

你好,
我目前正在研究用MATLAB优化一些计量经济学模型的LogLikelihood函数。问题是内部优化器似乎在最大化这些函数方面存在问题。因此，我想看看如果我提供一个好的梯度数值估计，情况是否会变得更好。于是，我找到了这个工具箱。然而，我有一些问题使它运行。下面是我试图优化的函数的一个非常简单的例子:

函数logli = logl(b,v,ncs,alt,y)
Y =重塑(Y,max(alt)，max(ncs));
X =重塑(exp(v(b))，max(alt)，max(ncs));
P = x(y == 1)./sum(x)';
P = max(P,0.0000001);
Logli = -sum(log(x(y == 1)./sum(x)'));
结束

在哪里
V是一个匿名函数
V = @(b) b(1)*x1 + b(2)*x2 + b(3)*x3
B是具有起始值的参数向量
B = [0,0,0]
NCS和Alt是标识符变量
Y是自变量

如果我想求一阶导数

(der, er) =派生(@ (b) logl (b, v, nc, alt, y), b)

但是得到错误消息:
使用派生函数时出错(第409行)
Fun没有返回正确的大小结果(Fun必须向量化)

如果我指定:
(der, er) =派生(@ (b) logl (b, v, nc, alt, y), b,“矢量化”,“不”)
它似乎只经过了b的一维
试图访问b(2);索引越界，因为数字(b)=1。

这些问题的数据结构为:

NCS Alt y x1 x2 x3
11 1 x11 x21 x31
12 0 x12 x22 x32
13 0 x13 x23 x33
2 10 x11 x21 x31
22 1 x12 x22 x32
23 0 x13 x23 x33
.．.

我很乐意听到任何关于如何在这种情况下使用导数函数的建议。

干杯
塞巴斯蒂安。

简单的答案是使用hessian函数。例如，我要计算二阶导数

d ^ 2 f /

函数的:

Y = sin(x + Y²)

在点(x,y) =(2,3)处。我将使用这些工具来验证，然后再使用符号工具箱来验证答案。

趣味= @(xy) sin(xy(1) + xy(2).^2);
格式long g
3[2],黑森(有趣)
ans =
0.999990206550529 - 5.99994123927655
5.99994123927655 - 36.0084988318012

你想要的第二个偏导是(1,2)元素，也就是黑森矩阵的非对角线元素。

现在使用符号方法验证该值。

符号x y
Fs = sin(x + y.²);

Dfs_dxdy = diff(diff(fs,x)，y)
dfs_dxdy =
-2*y*sin(y²+ x)

潜艇(dfs_dxdy {x, y}, 3 [2])
ans =
6 * sin (11)

vpa (ans)
ans =
5.9999412393042207423093893941531

这很好，尽管它计算了我们不需要的附加项。

我也可以用导数函数来做，尽管看起来我必须告诉导数函数这个函数没有正确地向量化。

Fun = @(x,y) sin(x + y.^2);

作为y的函数，在x = 2处，我们有

dfdx = @ (y)得到(@ (x)有趣的(x, y), 2,“矢量化”,“不”);

获得(dfdx 3“矢量化”,“不”)
ans =
5.99994123930711

好文件，我有个问题。我有一个目标函数，我对它求导wrt输入matlab。我想对它求导两次，一次对输入求导另一次对参数求导。但是如果我使用GRADEST来获得wrt输入的梯度，现在我有一个数字的答案，这是一个梯度。我怎么能第二次求梯度当我没有某种形式的一阶导数的时候。

伟大的文件，非常感谢。

然而，我目前面临一些问题与我的函数(糟糕的缩放参数)，因此与二阶导数矩阵(Hessian)。

这里是交易:我正在最小化代价函数(x 2 = (y - ymodel)^2)，其中我的模型是3非线性参数1，使用Nelder Mead算法。

使用您的包中的|hessian|函数使用最佳拟合参数，产生了一些奇怪的结果(因为Parameter_1的范围从[0.1:0.9]，Parameter_2从[10:30]和Parameter_3从[150:250])，也就是说，我需要重新缩放它们。我的问题是，实现这一点的最佳方法是什么，您的代码是否能够内生地扩展结果?

非常感谢你，

负载比fmincon的输出更好

很棒的提交，节省了我很多工作，谢谢!

Retrobytes——并非所有函数都能使用复杂参数，因此derived避免了这种方法。超对偶数似乎也有同样的缺陷，至少不需要投入大量的精力。

非常有用的工具，约翰，谢谢!

你们有人想过用超对偶数的方法来解决这个问题吗?我最近看到了这个:
http://www.stanford.edu/~jfike/SNUpresentation_rev4.pdf
和:
http://www.stanford.edu/~jfike/hyperdual.html
我想知道你们有没有人试过用复数代替复数?

绝对是很棒的工作!很容易上手，做得很好。这可能是自切片面包以来人类发生的最好的事情。还免费放上去?约翰，你是个好人!

是伟大的

非常感谢这个有用的工具!

我有一个问题，程序会从最初的点评估乐趣有多远。我有一个函数，它只在初始点周围的某个范围内存在(但在我的意义上足够大，例如从.1开始，需要>0)。但是，在计算车费金额后，该程序将在评估超出范围的点时崩溃。你能告诉我为什么会这样吗?

(顺便说一下，我的函数有一些随机性，这有关系吗?)

非常感谢!

太好了!它的一些前辈也做精确的导数，但这个是特别容易使用梯度和黑森。

约翰H,

看起来不错。我试着用黑森函数求最大似然的费雪信息矩阵。然而，黑森程序选择的参数是不允许在我选择的模型。我如何限制黑森算法为允许的参数值

你的问题在于你调用代码的方式。

H=hessian(@(p) @reg2002， [0.03 3 -5])

下面是错误信息的一部分:

“? ?从function_handle转换为double时发生以下错误:
使用==> double时出错
不可能从function_handle转换为double。”

reg2002是一个函数m文件吗?您试图将变量p传递到这个文件中?如果是这样，那么你可以这样做

H=hessian('reg2002'， [0.03 3 -5])

或者，你可以这样做

H=hessian(@reg2002， [0.03 3 -5])

在这两种情况下，它都假设reg2002是一个m文件，并传入一个包含两个元素的变量作为vector。它会返回一个2x2的黑森矩阵。使用@reg2002将其转换为函数句柄。

请查看您所看到的错误消息。Matlab看到:

@ @reg2002 (p)

并尝试将函数句柄作为函数的结果返回。然后我的代码查看它，并尝试使用函数句柄作为一个数字。当然失败了。

最后，你可能还做过这个

H=hessian(@(p) reg2002(p)， [0.03 3 -5])

在这里，您将创建一个匿名函数，一个变量p的函数。它是代码可以正确使用的另一个函数句柄。

你好，约翰，非常感谢你给我这么棒的套房!

我在使用m文件中的函数时遇到了一个问题。我的代码是:

H=hessian(@(p) @reg2002， [0.03 3 -5])

?？?从function_handle转换为double时发生以下错误:
使用==> double时出错
不可能从function_handle转换为double。

在==>中衍生330的错误
F_x0 (j) = fun(x0(j));

错误==> hessdiag在59
(高清(印第安纳州),犯错(印第安纳州),finaldelta(印第安纳州)]=推导(…

错误在==> hessian在74
[hess,err] = hessdiag(fun,x0);

我不知道哪里不对。我犯了什么错误吗?非常感谢。

艾玛

Felipe -请注意，提交的#15235和#26807实际上出现在我的工具(#13490)之后。提交的文件在FEX上被分配顺序编号。我看不出有什么理由承认以后会提交有关外汇交易所的报告。一般来说，由于向前确认似乎需要某种类型的预知，这是相当困难的。事实上，维基百科文章中的大部分引用都是在这个工具编写并发布在FEX上几年后出现的。

我会根据你的评论修改标题，因为这个短语可能会让一些人感到困惑。

标题中的“自动”可能会让观众感到困惑，因为它实现了“自动区分”，参见<http://en.wikipedia.org/wiki/Automatic_differentiation>和提交# 26807。

您可能需要使用确认字段，以便相关的提交#15235反向链接到您的。

还没试过，这就是我不给它评分的原因。

嗨,约翰,

非常感谢你给我这么棒的套房!

关于hessian函数的信息明确指出，它不适合经常用于代价高昂的目标函数上。我还是试了，因为这就是我需要的。大多数情况下，我会得到一些合理的信息(很难检查，但我以后会担心这个问题)，但有时我会得到NaN。知道是什么问题吗?

也许需要更多的细节:我想微分的函数是一个似然函数。它是时间的函数，我需要每个时间步的最大值处的黑森。(最大值是与最后一个时间段相关联的一个点。)该函数有5个参数。

任何帮助或想法，感谢!

非常感谢，
二

如果你正在寻找所有的二阶偏微分，那么是的，在黑森计算上的循环似乎是解决方案。

嗨,约翰,
非常感谢您的快速回复。
雅可比矩阵的确是我所需要的，我应该先研究一下整个包装。

我在看hessian函数，它不允许我使用向量函数(也许可以理解，否则输出将是一个张量)。

继续看矩阵向量的例子
@ * x (x);
我是否应该遍历黑森的行并调用黑森m次来得到所有的黑森?
i = 1: m
f = @(x) A(i，:)*x;//返回标量
[h(i) herr(i)] = hessian(f, x0);
结束

谢谢,
——一位

嗨,艾登,

诀窍在于，导数函数本身只取单个变量的函数。假设我有一个含三个变量的标量函数。因此，它接受一组三个变量，返回一个数字。如果你想计算这个函数的偏导数，你得到的是一个数。但这是偏导数，不是导数。三个变量的标量函数的导数由梯度向量给出，由三个偏导数组成。所以我提供了工具等级最高的。将它用于多变量的标量函数。

现在，假设我有一个n个变量的函数，产生m个输出结果?(顺便说一下，这是你有的。)现在你必须计算一个雅可比矩阵，一个m × n的偏导数矩阵。再说一遍，derivative不是工具，gradest也不是。这是一个不同的问题需要不同的工具，所以我写了雅可比矩阵。

>> A =魔法(3)
一个=
8 1 6
3 5 7
4 9 2

>> [J,Jerr] = jacobianest(fun，[1;2;3])
J =
8 1 6
3 5 7
4 9 2

Jerr =
7.7667e-14 9.7084e-15 5.3057e-14
4.3706e-14 7.09e-14
3.8833e-14 9.4059e-14 1.9417e-14

在这种情况下，您真的只想知道结果的特定元素的导数，相对于几个输入中的特定单个变量，您需要定义一个函数，其中OTHER变量保持固定。因此……

>> x1 = 1;
>> x2 = 2;
>> fun2 = @(x3) [0 1 0]*fun([x1;x2;x3]);
查理——> > [c15,德尔·]=派生(fun2 3“矢量化”,“不”)
c15 =
7
derr =
4.785 e-14

因此，fun2修复了前两个变量，然后提取出第二个输出结果。只要我告诉它这个函数不是向量化的，导数函数就可以起作用。它没有设置为获取输入向量并为每个输入生成一个结果。Fun2只接受一个标量输入，也只产生一个标量输出。

约翰

嗨,约翰,
谢谢你可爱的工具。
但是我不知道怎么用它来微分多变量函数。例如,
f = @(x) A*x

A是一个矩阵，x是一个向量。

你的代码似乎是这样做的:

%遍历x0的元素，将其简化为
%标量问题。抱歉，矢量化不是
%完成，但这只是一个循环。
Der = 0 (nx0);
Errest = der;
Finaldelta = der;
对于I = 1:n
.．.
结束

而在我看来，A*x并没有真正化简为单个标量问题。有什么建议吗?

对不起，约翰，解决了这个问题，是我的错…再次感谢这些有用的代码。
亚历山德罗

嗨,约翰,

我试着运行Rosenbrock示例，包括在您的例程hessian的帮助下，即我运行

Rosen = @(x) (1-x(1))。^2 + 105*(x(2)-x(1).^2).^2;
[h,err] = hessian(rosen，[1 1])

我得到以下错误:

“? ?使用==> subsindex时出错
函数subsindex没有为类function_handle的值定义。"

你知道怎么解吗?
谢谢你！
亚历山德罗

我相信fmincon已经改进了，所以返回的hessian现在是一个更好的近似值。当然，如果您提供了粗麻绳，则不需要更新。

嗨Errico,

在另一篇文章中，你说matlab不能很好地计算黑森。这就是黑森在fmincon中不能更新的原因吗?其次，提供一个用户提供的Hessian(例如，通过派生est例程)是否能够提高fmincon在所有情况下更新Hessian的能力?谢谢!

嗨Errico

我试着扩展你的代码来推导双线性矩阵不等式。该程序仍然是一个粗略的形式，特别是我没有实现估计的误差部分，但我能够正确地产生BMI函数的导数。如果您有兴趣将项目寄给您，请让我知道。总有一天我会尝试实现误差项。

Kimonas

嗨,安德鲁,

重点是，导数函数使用了一系列嵌套的近似，对步骤使用了几何序列的点。例如，假设我们选择计算tan(x)在x = *2/3处的导数。

[D,errest] =导数(@tan,pi*2/3)

D =
4

呃=
7.0302 e-13

Tan(x)在x = /2处有一个奇点，所以如果步长选择得足够大，近似值实际上会跨越这个奇点。这将导致相应的估计是无用的，但由于产生了多个估计，其中一些将是好的，并且我可以找出序列中哪些估计是好的，哪些是坏的。

我选择了预测误差最小的估计值，以及一些附加的标准。

近似的嵌套序列的一个好处是它允许我处理有缩放问题的函数。我几乎总能找到可行的解决方案。

约翰

在derivest_demo。m，我找到了下面的声明:
“尽管一个中心规则可能会把一些样本放在错误的地方，但它仍然可能成功。”

它如何把样本放在“错误的地方”?

嗨Prasanna,
理论上，我可以用导数来提供更高的导数。在实际应用中，高阶数值微分是比较困难的。微分是一个噪声放大器。它吸收你身体系统中的任何微小噪音，并增强它。当你达到四阶导数的时候，最低有效位的噪声就会开始导致数值问题。因此，我将其限制在这里。如果再进一步，我就不再能够相信一些测试用例中的结果了。因为我必须在某个地方设置极限，我选择4作为它的最大阶导数。

您可以使用提供的误差估计来衡量我认为我能够做得有多好。举个例子，试试一个简单的三角函数。在x = /4处求导sin(x)结果总是正负根号(2)/2，取决于导数的顺序。

>> fun = @罪恶;
>> [D,DE] = derivatives est(fun,pi/4，' derivatives ativeorder'，1)
D =
0.70711
德=
7.0617 e15汽油

对于一阶导数，我几乎可以给出完全的精度。然而，二阶导数并不完全是可验证可计算的，误差估计反映了这一事实。

>> [D,DE] = derivatives est(fun,pi/4，' derivatives ativeorder'，2)
D =
-0.70711
德=
1.666 e-12

>> [D,DE] = derivatives est(fun,pi/4，' derivatives ativeorder'，3)
D =
-0.70711
德=
6.6146 e-11

>> [D,DE] = derivatives est(fun,pi/4，' derivatives ativeorder'，4)
D =
0.70711
德=
1.9122 e-05

当我求四阶导数的时候，导数函数告诉我只相信结果中的前四五位数字。如果我在这里允许更高的阶数，结果将迅速成为一个虚拟随机数生成器。

现在，只要小心，我们可以修改代码，使之适用于五阶导数。只需查找par.DerivativeOrder，并查看所做的工作，将代码扩展为5阶导数。快速检查代码会发现只有一两个地方需要认真理解所使用的方法。(老实说，这样的代码修改还是留给我吧。)

你能嵌套派生调用吗?是的，理论上你可以这样做来产生更高阶的导数，误差估计会告诉你它是否认为它产生的结果是垃圾。这个测试计算sin(x)在/4处的五阶导数。

>>格式长g
>> d2 = @(x) derivatives est(fun,pi/4 + x，'derivative '，2);

>> [d5,d5err] = derivatives est(d2,0，' derivatives '，3)
d5 =
0.707106775613828
d5err =
4.17633428512851 e-08

> >√6 (2)/ 2
ans =
0.707106781186548

请注意，derived认为它能够很好地做到这一点，并且它刚好告诉您不要再相信结果中的数字。事实上，在这种情况下，嵌套的结果可能比直接计算5阶导数更准确。所以二阶导数的三阶导数比直接计算五阶导数更准确(但有点慢!)

HTH,
约翰

你好约翰,
这段代码很棒。我有一个问题，我有一个函数，我没有解析表达式，因此不能向量化它。但是我可以对它进行快速的数值计算我知道它有平滑的一阶，二阶，三阶和四阶导数通过对它应用导数测试函数。我还需要一个导数来完成我的计算，我想知道我是否应该用导数测试来计算四阶导数，并取它的一阶导数来得到五阶导数。我想知道是否有什么特殊的原因使得这段代码只能达到四阶导数?

非常感谢，你帮了大忙。这是一组很好的函数，非常受欢迎。

毫不奇怪，在fminsearch中调用一个循环来对几个点进行评分并不是世界上最快的事情，但这是一种权衡，而且它似乎确实完成了工作!

我已经提供了相关的函数gradest来计算梯度，压缩与衍生。例如，我将为参数k1,k2,k3,k4,k5选择任意值，然后定义一个函数句柄。

K1 = 1/2;
K2 = 1/3;
K3 = 1/4;
K4 = 1/5;
K5 = 1/6;

乐趣= @(xyz) k1*xyz(1)*(1-…
exp (k2 * (xyz (2) ^ k3) * (xyz (3) ^ k4) * (xyz (1) ^ k5)))

现在，叫gradest…

Gradest (fun，[1 2 3])
ans =
0.21997 0.018836 0.010046

对于多个点，围绕gradest调用进行循环。

约翰

我可能忽略了一些明显的东西，但是有没有一种聪明的方法来使用导数来计算一组偏导数?

我有一个这样的函数
F = k1 * x *(1 - exp(-k2*(y^k3)*(z^k4)*(x^k5))
并且希望能够计算偏导数df/dx, df/dy和df/dz对于特定的x y z集合。

我这样做是为了计算最小二乘拟合的误差权重——我在x、y和z中有不确定性，我想把它们映射到f上，以确定最小化的权重因子(所以我想要每个点的梯度)。我更倾向于用数字来确定导数因为这样我可以更灵活地尝试不同形式的函数，等等。

多谢。

较大的相对报告误差可能意味着计算出的导数与其中的不确定性相比较小。这也可能意味着在计算中存在一些问题，使得导数不确定。仅从这一描述很难知道发生了什么事。

一种选择是用导数计算问题中的导数，而不是hessian工具。虽然这样做的效率较低，但对于克服我在derived本身中提供的问题，还有更多的选择。这样你就能确定问题是什么，并有希望克服它。

你好,

我有一些问题，有时当我计算黑森时，我得到err包含的条款比报告的黑森中的条款更大。造成这种情况的原因是什么?有没有什么通用的方法来重写一些东西来避免这个问题?

非常感谢您的快速回复
我会得到第14或第7版

嗨,V,

您使用的一定是旧的MATLAB版本。下面是我在R2010a中尝试您的示例时发生的情况:

>> myfun1 = inline('exp(x)');
[d,犯错]=派生(myfun1 1“矢量化”,“不”)

d =

2.7183

呃=

2.1641 e-14

所以我更仔细地看了你的评论。我看到警告信息泄露了答案。

文件:C:\MATLAB6p5\work\ derivative。m行:463列:1
无与伦比的“结束”。

您正在运行版本6.5(即版本13)。派生程序是在版本14或版本7中编写的。因此，它使用了一些会使您的MATLAB版本非常不满意的东西。

对不起,
约翰

嗨,约翰,

我试过了

Myfun1 = inline('exp(x)');
[d,犯错]=派生(myfun1 1“矢量化”,“不”)

我得到了下面的误差，
谢谢你能提供的任何帮助

文件:C:\MATLAB6p5\work\ derivative。m行:463列:1
无与伦比的“结束”。
(输入“warning off MATLAB:m_warning_end_without_block”来屏蔽此警告。)
文件:C:\MATLAB6p5\work\ derivative。m行:523列:1
无与伦比的“结束”。
(输入“warning off MATLAB:m_warning_end_without_block”来屏蔽此警告。)
文件:C:\MATLAB6p5\work\ derivative。m行:538列:1
无与伦比的“结束”。
(输入“warning off MATLAB:m_warning_end_without_block”来屏蔽此警告。)
文件:C:\MATLAB6p5\work\ derivative。m行:573列:1
无与伦比的“结束”。
(输入“warning off MATLAB:m_warning_end_without_block”来屏蔽此警告。)
文件:C:\MATLAB6p5\work\ derivative。m行:659列:1
无与伦比的“结束”。
(输入“warning off MATLAB:m_warning_end_without_block”来屏蔽此警告。)
文件:C:\MATLAB6p5\work\ derivative。m行:759列:1
无与伦比的“结束”。
(输入“warning off MATLAB:m_warning_end_without_block”来屏蔽此警告。)
?？?使用==>阶乘出错
N必须是正整数。

C:\MATLAB6p5\work\derivest.m(fdamat)
第564行==> c = 1./阶乘(1:2:(2*nterms));

C:\MATLAB6p5\work\derivest.m
第275行==> fdarule = [1 0]/fdamat(par.StepRatio,1,2);

另一个非常优秀的代码，我希望TMW每年为他们广泛扩展他们的软件而给你买一份漂亮的圣诞礼物。

你好,

黑森不工作，无论是在任何例子中，也当我试图使用它。它给出的信息是:
函数subsindex没有为类function_handle的值定义。
所有其他功能似乎都很棒。
最好的

亚当

嗨,约翰,

感谢您及时详细的回复!我假定派生程序处理所有情况。现在我知道雅可比矩阵是一个独立的函数。再次感谢!还有很棒的代码!

米歇尔

嗨,米歇尔,

问题是，当您定义myfun时，它不是一个标量值函数。它返回一个向量参数。因此我们可以这样做:

myfun =内联(“[exp (x); exp (2 * x)] ");

myfun (2)
ans =
7.3891
54.598

查看标量输入返回两个数字。

(顺便说一下，我强烈建议使用函数句柄和匿名函数，而不是内联函数。内联函数的计算相对来说非常慢。)

导数代码计算标量值函数的导数。将函数更改为返回标量值，如下面的myfun1和myfun2所做的那样，派生函数将按设计工作。

Myfun1 = inline('exp(x)');
[d,犯错]=派生(myfun1 1“矢量化”,“不”)

d =
2.7183
呃=
9.3127 e15汽油

Myfun2 = inline('exp(2*x)');
[d,犯错]=派生(myfun2 1“矢量化”,“不”)
d =
14.778
呃=
1.8401 e-13

那么如何定义向量值函数的导数呢?这在数学上被定义为雅可比矩阵。事实上，在最导数套件中，我已经提供了雅可比矩阵函数来计算你想要的东西。它作用于向量值函数，返回雅可比矩阵。这里，由于这些函数的参数是一个标量，雅可比矩阵就是一个长度为2的向量。

[d,err] = jacobianest(myfun,1)
d =
2.7183
14.778
呃=
3.6841 e-14
1.4736 e-13

约翰

非常感谢您提供如此有用的套件!

我在测试代码时输入了以下两个命令:

> > myfun =内联(“[exp (x); exp (2 * x)]”)
> > (d,犯错)=派生(myfun 1“矢量化”,“不”)

但得到的错误信息是

?？?在赋值A(I) = B中，B和中的元素数量
我一定也一样。

错误在==> 的导数在362
F_plusdel (j) = fun(x0i+h*delta(j));

我把它贴在这里，因为我认为你可能会比我更快地发现问题。你能告诉我吗?

提前感谢你的帮助!

Benn,

我猜您已经编写了自己的阶乘函数，可能是作为家庭作业，因为错误消息显示阶乘的递归实现。这将导致派生失败，因为您的版本没有正确地向量化。通过向量化，我的意思是当使用数字列表调用阶乘时，Mathworks版本的阶乘计算每个数字的阶乘。

你可以通过在命令行输入以下命令来了解调用哪个阶乘版本:

哪个阶乘-全部

MATLAB将显示将使用阶乘函数的哪个版本。您应该看到您编写的版本列在最顶部。

起初，我认为您可能只是在使用旧版本的MATLAB，但这是不可能的，因为Mathworks永远不会写出具有递归定义的阶乘函数。(递归阶乘定义是教学生递归编码风格的好方法，但在实践中，这是一种糟糕的编码风格，因为仅仅为了将一组连续数字相乘，就会产生无关的、浪费的函数调用开销。)

解决方案是从搜索路径中删除factorial的坏版本(家庭作业)。无论如何，与MATLAB已经提供的阶乘版本相比，它都是低效的，而且您的版本没有正确地向量化。至少，把它移动到搜索路径的另一个地方，这样matlab就会使用正确的阶乘版本。

约翰

的帮助!只是下载了工具箱，并尝试在派生文档中运行最基本的示例:

(d, e) =派生(@ (x) exp (x) 1)

得到了一堆这样的错误

错误==> !在8
X = n*阶乘(n-1);

错误==> !在8
X = n*阶乘(n-1);

错误在==>衍生>fdamat在564
C = 1./factorial(1:2:(2*nterms));

错误在==>派生在283
fdarule = [0 1 0]/fdamat(par.StepRatio,1,3);

工作真的很好!很有帮助!谢谢……

这比fminsearch的工作更好，当被用作阻尼牛顿方法寻找根具有极高的灵敏度。Fminsearch速度太慢，而且经常超出根目录，即使容忍度很高，而且fzero根本不起作用。推导很快。
这是我第二次使用John D’erico的优秀代码。
谢谢!

我刚刚上传了一个新的版本，修复了复杂域分化的问题。例如，当前版本现在得到了正确的结果:

exp (1 + 1)
ans =
1.4687 + 2.2874i

获得(@ (z) exp (z), 1 + i)
ans =
1.4687 + 2.2874i

当我用这个函数求复数处的导数时，问题出现了。例如:
获得(@ (z) exp (z), 1 + i)
ans =
1.46869393991590 - 2.28735528717885i
然而，exp(z)在z(复数)处的导数应为exp(z)
exp (1 + 1)
ans =
1.46869393991589 + 2.28735528717884i
它是上面的共轭。

摇滚

摇滚

工作得很好。我遇到的情况是，使用hessians是有问题的，因为要计算的函数具有瞬时复杂的性质。hessians可以很好地用于从不使用任何复杂算术的函数，但如果内部函数使用任何复杂算术，则要注意使用复杂步骤，即使输入和输出都是实数。下载了这个包，修改了函数名，它就可以开箱即用了——速度慢一些，但“防弹”有时比“快速”要好。好的文档等等。我很感激hessiancs使用了与你使用的类似的形式-当我需要更多的速度和有时间修复我的损失函数时，我会回到它。

斯科特

谁能告诉我什么与光速有冲突?如果可能的话，我会试着把它修好。它可能是我选择的函数名。否则，这只是基本的Matlab代码。

非常有用!我能够使用它来检查我的符号衍生的黑森矩阵在一个优化例程…optimset(' derivative vecheck '，'on')用数字检查梯度，但显然不是黑森。

以防其他人有问题，我发现这与matlab Lightspeed工具箱不兼容，但这无疑是Lightspeed工具箱的错误。移除光速修复了这个问题。

这个网站很有帮助:)

我有时会得到这样的错误:
?？?对矩形空矩阵的不正确赋值。

错误在==>派生在457
[errest(i)，ind] = min(误差);

-为什么会这样?

最好的!