数值不稳定？

Question

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

在下面的mwe中你可以看到obj是100x100对称矩阵的100-R（这里我使用R=3）最小特征值的总和u'*u.在真正的全局最小值 obj "looks good" 附近时，当我绘制 objective 函数时，在更大的 beta 值下评估 objective 函数变得非常不稳定（here 或运行 mwe 您可以看到出现了多个局部最小值（和最大值），与小于真实全局最小值的 obj(beta) 值相关。

我的猜测是有某种 "numerical instability" 正在发生，但我无法找到来源。


%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) 处的 objective 函数。我注意到 eig(u'*u) 的一些特征值是负的（见对象 ev），当通过构造矩阵 u'*u 是半正定的

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

最后，objective 函数 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 浮点数精确到不到 16 位有效（十进制）数字。这意味着即使 double 浮点数的表示错误也会在这里产生错误使每个数字"insignificant"。）这已经解释了不良行为。您应该注意，通常我们可以非常准确地计算最大的特征值，较小（幅度）的误差通常会增加。

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

唯一可以挽救的情况是R=0，那么我们实际上想要计算 all 个特征值的总和，这恰好是我们矩阵的 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');

数值不稳定？

Numerical instability?

algorithm

precision

matlab

convergence