algorithm - 数值不稳定？

Question

我在一个程序中工作，该程序涉及在标量 beta 上优化某些目标函数 obj。真正的全局最小值 beta0 设置为 beta0=1。

在下面的 mwe 中，您可以看到 obj 构造为 100x100 对称的 100-R(这里我使用 R=3)最小特征值的总和矩阵 u'*u。当我绘制以更大的 beta 值评估的目标函数时，围绕真正的全局最小值 obj “看起来不错”，目标函数变得非常不稳定( here 或运行mwe 你可以看到多个局部最小值(和最大值)出现，与 obj(beta) 的值相关联，小于真正的全局最小值)。

我的猜测是存在某种“数值不稳定”，但我无法找到来源。

%Matrix dimensions
N=100;
T=100;

%Reproducibility
rng('default');

%True global minimum
beta0=1;

%Generating data

l=1+randn(N,2);
s=randn(T+1,2);
la=1+randn(N,2);
X(1,:,:)=1+(3*l+la)*(3*s(1:T,:)+s(2:T+1,:))';
s=s(1:T,:);
a=(randn(N,T));
Y=beta0*squeeze(X(1,:,:))+l*s'+a;

%Give "beta" a large value
beta=1e6;

%Compute objective function
u=Y-beta*squeeze(X(1,:,:));

ev=sort(eig(u'*u));  % sort eigenvalues

obj=sum(ev(1:100-3))/(N*T); % "obj" is sum of 97 smallest eigenvalues

这会评估 obj(beta=1e6) 处的目标函数。我注意到 eig(u'*u) 中的一些特征值是负的(参见对象 ev)，当构造矩阵 u'*u 是半正定的

我猜这可能与浮点运算问题有关，并且可能(部分)是我函数不稳定的答案，但我不确定。

最后，目标函数 obj 在 beta 的各种值下评估如下:

% Now plot "obj" for a wide range of values of "beta"
clear obj
betaGrid=-5e5:100:5e5;

for i=1:length(betaGrid)
    u=Y-betaGrid(i)*squeeze(X(1,:,:));
    ev=sort(eig(u'*u));
    obj(i)=sum(ev(1:100-3))/(N*T); 
end

plot(betaGrid,obj,"*")
xlabel('\beta') 
ylabel('obj')

这给出了 this figure ，它显示了 beta 的极值变得多么不稳定。

Answer 1

这里的关键是要注意计算特征值可能是一个难题。实际上 condition number因为这个问题是K = norm(A) * norm(inv(A))(不要这样计算，使用cond()。这意味着在计算输出时，输入(即矩阵项)中的(相对)扰动被条件数放大。我稍微修改了您的代码以计算并绘制每个步骤中的条件数。事实证明，对于一个大的您感兴趣的范围的一部分大于 10^17，这很糟糕。(请注意，double float 精确到不完全是 16 位有效(十进制)数字。这意味着即使double float 的表示错误会在这里产生错误，使每个数字都“无关紧要”。)这已经解释了不良行为。您应该注意，通常我们可以非常准确地计算最大的特征值，较小(幅度)的错误通常会增加。

如果条件数更好(接近 1)我会建议计算奇异值，因为它们恰好是特征值(由于对称性)。 svd 在数值上更稳定，但这真的很糟糕条件即使这样也无济于事。在以下修改中最后一个片段我添加了一个图表来绘制条件数。

唯一可以挽救的情况是 R=0，那么我们实际上想要计算所有特征值的总和，这恰好是我们矩阵的 trace，只需将对角线条目。

总结一下:这个问题似乎有一个固有的坏条件，所以你如何计算它并不重要。如果您对同一问题有完全不同的表述，可能会有所帮助。

% Now plot "obj" for a wide range of values of "beta"
clear obj
L = 5e5; % decrease to 5e-1 to see that the condition number is still >1e9 around the optimum
betaGrid=linspace(-L,L,1000);
condition = nan(size(betaGrid));
for i=1:length(betaGrid)
    disp(i/length(betaGrid))
    u=Y-betaGrid(i)*squeeze(X(1,:,:));
    A = u'*u;
    ev=sort(eig(A));
    condition(i) = cond(A);
    obj(i)=sum(ev(1:100-3))/(N*t); % for R=0 use trace(A)/(N*T);
end

subplot(1,2,1);
plot(betaGrid,obj,"*")
xlabel('\beta') 
ylabel('obj')
subplot(1,2,2);
semilogy(betaGrid, condition);
title('condition number');