Matlab code has a systematic error

Question

I have the following code to calculate the number of real eigenvectors of a normally distributed matrix:

x=zeros(1,100); %make an 1x100 array of 0's
n = 1000;

for i=1:100
   Nr=0;
   a=normrnd(0,1,n); % make a nxn matrix N(0,1)
   C=eig(a);  %make c the eigen values of a
   for k=1:n
      if (isreal(C(k))==1) %loop through each value of c andcheck if it is real
          Nr=Nr+1; %if it is real increment Nr
      end
   end
   x(i)=Nr/sqrt(n);
end
Estimation_Of_C=mean (x)
Estimation_Of_Error=std (x)

According to http://math.mit.edu/~edelman/homepage/papers/howmany.pdf this should be equal to $\sqrt{\frac{2}{\pi}}$ however for whichever value of $n$ I use (bigger and smaller than 1000), I more often than not get Estimation_ofC> 0.8 so I appear to have a small systematic error.

Am I simply wanting the code to be too precise or do I have an actual error in my code?

Answer 1

The value $\sqrt \frac{2}{\pi}$ is in the limit when $n\rightarrow \infty$. If you use the asymptotic formula instead $$\frac{N_r}{\sqrt{n}} = \sqrt{\frac{2}{\pi}}+\frac{1}{2}\frac{1}{\sqrt{n}}+O(\frac{1}{n})$$ in your case $\sqrt{\frac{2}{\pi}}+\frac{1}{2}\frac{1}{\sqrt{1000}} \approx 0.8136959$

Answer 2

I would imagine that the issue is with your use of isreal. eig is an iterative solver, which means that it is possible that some eigenvalues are being computed with a very small but non-zero imaginary part due to truncation error. For instance

isreal(4+1e-12i)

will return false.

Instead, I would condition a check on the complex part of each eigenvalue using imag. Rather than saying isreal(C(k)) == 1, I would do something like:

epsilon = 1e-10;
if abs(imag(C(k)) < epsilon
    ...
end

This might improve your results.

Answer 3

You need to reset Nr to zero at the start of each iteration. Though you can make this moot by avoiding the inner loop completely. Also note that isreal( complex(1) ) will return false and not true. See help isreal for an explanation. Also, the numerical issue is valid, so using something the like following is probably the proper way to go about it (though using angle off the real axis might be better):

x=zeros(1,100); %make an 1x100 array of 0's
n = 1000;
for i=1:100
   a=normrnd(0,1,n); % make a nxn matrix N(0,1)
   C=eig(a);  %make c the eigen values of a
   Nr = length( find( abs(imag(C)) < 10*eps ) );
   x(i)=Nr/sqrt(n);
end
Estimation_Of_C=mean (x)
Estimation_Of_Error=std (x)

Estimation_Of_C =

         0.813337814195307


Estimation_Of_Error =

         0.131508727674423

Here is the histogram for 1000 iterations: