0

对于上下文，我在 MATLAB 中有一个小项目，我尝试复制一个涉及牛顿算法优化的算法。虽然我的问题主要与 MATLAB 相关，但也许是我缺乏深厚的背景知识使我无法找到解决方案，因此如果需要，请随时将我重定向到适当的 StackExchange 站点。

我需要计算梯度向量和 Hessian 矩阵以进行优化的函数是：

function [zi] = Zi(lambda,j)
    zi = m(j)*exp(-(lambda*v_tilde(j,:).'));
end

function [z] = Z(lambda)
    res = arrayfun(@(x) Zi(lambda,x),1:length(omega));
    z = sum(res);

end

function [f] = F(lambda)
    f = log(Z(lambda));
end

其中omega和v_tilde是 n 个 d 维向量的矩阵，lambda是函数的 d 维参数。（现在，m(j)只是选择器（1 或 0），但算法允许对它们进行细化，因此不应删除它们。我使用 Derivest Suite 以数值方式计算梯度和 Hessian，虽然逻辑上对于高维很慢, 算法作为一个整体有效。我使用 sym 包实现了相同的解决方案，这样我就可以提前计算一些修正 n 和 d 的梯度和 Hessian，这样就可以在需要时快速评估它们。这将是符号版本：

V_TILDE = sym('v_tilde',[d,n]) 
syms n k

lambda = sym('lambda',[d,1]); 
F = log(M*exp(-(transpose(V_TILDE)*lambda)));

matlabFunction( grad_F, 'File', sprintf('Grad_%d_dim_%d_n.m',d,n_max),  'vars',{a,lambda,V_TILDE});
matlabFunction( hesse_F, 'File', sprintf('Hesse_%d_dim_%d_n.m',d,n_max), 'vars',{a,lambda,V_TILDE});

由于 n 是固定的，因此无需遍历 omega。this的梯度和Hessian可以通过sym的对应函数得到，然后存储为matlabFunctions。

但是，当我针对某些值测试这两种实现时，它们不匹配，令人惊讶的是，hessian 矩阵的值匹配而梯度的值不匹配（数值计算正确），并且牛顿算法迭代直到这些值只是NaN。这些是 d=2 和 n=8 的一些示例值：

Omega:
12.6987   91.3376
95.7507   96.4889
15.7613   97.0593
95.7167   48.5376
70.6046    3.1833
27.6923    4.6171
9.7132   82.3458
95.0222    3.4446

v:
61.2324
52.2271

gNum =          HNum =    1.0e+03 *

  8.3624                 1.4066   -0.5653
 -1.1496                -0.5653    1.6826

gSym =          HSym =    1.0e+03 *

 -52.8700                1.4066   -0.5653
 -53.3768               -0.5653    1.6826

如您所见，HNum 和 HSym 的值匹配，但梯度不匹配。

我很乐意提供更多上下文信息、代码片段或任何有帮助的东西。先感谢您！

编辑：根据要求，这是一个最小的测试。问题基本上是 gNum 和 gSym 的值不匹配（上面更长的解释）：

omega = [[12.6987,   91.3376];[95.7507, 96.4889];[15.7613, 97.0593];
[95.7167, 48.5376];[70.6046, 3.1833];[27.6923, 4.6171];[9.7132,  82.3458];
[95.0222, 3.4446]];

v = [61.2324;52.2271];

gradStr = sprintf('Grad_%d_dim_%d_n',length(omega(1,:)),length(omega));
hesseStr = sprintf('Hesse_%d_dim_%d_n',length(omega(1,:)),length(omega));
g = str2func(gradStr);
H = str2func(hesseStr);
selector = ones(1,length(omega)); %this will change, when n_max>n

vtilde = zeros(length(omega),length(omega(1,:)));
for i = 1:length(omega)
    vtilde(i,:) = omega(i,:)-v;
end

lambda = zeros(1,length(omega(1,:))); % start of the optimization

[gNum,~,~] = gradest(@F,lambda)
[HNum,~] = hessian(@F,lambda)
gSym = g(selector,lambda.',omega.')
HSym = H(selector,lambda.',omega.')

注意：DerivestSuite 是一个小型库（~6 个源文件），可以在https://de.mathworks.com/matlabcentral/fileexchange/13490-adaptive-robust-numerical-differentiation下获得

4

0 回答 0