Numerical instability?

Question

I am working in a program that concerns the optimization of some objective function obj over the scalar beta. The true global minimum beta0 is set at beta0=1.

In the mwe below you can see that obj is constructed as the sum of the 100-R (here I use R=3) smallest eigenvalues of the 100x100 symmetric matrix u'*u. While around the true global minimum obj "looks good" when I plot the objective function evaluated at much larger values of beta the objective function becomes very unstable (here or running the mwe you can see that multiple local minima (and maxima) appear, associated with values of obj(beta) smaller than the true global minimum).

My guess is that there is some sort of "numerical instability" going on, but I am unable to find the source.

%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

This evaluates the objective function at obj(beta=1e6). I have noticed that some of the eigenvalues from eig(u'*u) are negative (see object ev), when by construction the matrix u'*u is positive semidefinite

I am guessing this may have to do with floating point arithmetic issues and may (partly) be the answer to the instability of my function, but I am not sure.

Finally, this is what the objective function obj evaluated at a wide range of values for betalooks like:

% 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 gives this figure, which shows how unstable it becomes for extreme values for beta.

Answer 1

The key here is noticing that computing eigenvalues can be a hard problem. Actually the condition number for this problem is K = norm(A) * norm(inv(A)) (don't compute it this way, use cond(). This means the the an (relative) perturbation in the inpute (i.e. the matrix entries) gets amplified by the condition number when computing the output. I modified your code a little bit to compute and plot the condition number in each step. It turns out that for a large part of the range you are interested in it is greater than 10^17, which is abysmal. (Note that the double floating point numbers are accurate to not quite 16 significant (decimal) digits. This means even the representation error of double floating point numbers will here produce errors that make every digit "insignificant".) This already explains the bad behaviour. You should note that usually we can compute the largest eigenvalues quite accurately, the errors in the smaller (in magnitude) ones usually increase.

If the condition number was better (closer to 1) I would have suggested computing the singular values, as they happen to be the eigenvalues (due to the symmetry). The svd is numerically more stable, but with this really bad condition even this will not help. In the following modification of the final snippet I added a graph that plots the condition number.

The only case where anything is salvageable is for R=0, then we actually want to compute the sum of all eigenvalues, which happens to be the trace of our matrix, which can easily be computed by just summing the diagonal entries.

To summarize: This problem seems to have an inherent bad condition, so it doesn't really matter how you compute it. If you have a completely different formulation for the same problem that might help.

% 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');