MATLAB: Incosistent matrix multiplication and probelm with covarince matrix

floating-point errorsindefinite covarince matrixMATLAB

I have a covariance matrix (S) which was computed using
   S=J*C*J';
where C is a diagonal matrix and J is another matrix. Both C and J are functions of some parameters. My problem is that the matrix S should be symmetric as is clear from the above expression. But MatLab does not return a symmetric covariance matrix. I checked it by computing Q:
Q=S-S';
Now, all the elements of Q should be zeros. But they are not as you can see in attached Q.mat file. I understand that this may be due to small some floating point arithmetic but I need it to be exactly symmetric because otherwise, my eigenvalues become complex which is physical implausibility. I have tried a fix for this
 ST=triu(S); % get upper-triangular elemnets of S
 S=ST+ST'-diag(diag(S));
It makes my matrix symmetric. But the next problem is that since it is a covariance matrix it should be positive-definite (at-least semi-positive definite) but it gives me very small negative eigenvalues which again might be due to floating point arithmetic. So I fix this again by forcefully making the negative eigenvalues equal to zero as follows
 [V,D]=eig(S);
 D(D<0)=0;
 S=V*D*V';
It removes the negative eigenvalue problem but it again makes the matrix not exactly-symmetric. So it seems that I am trapped in this cycle. Any suggestion for possible solutions?

Best Answer

The really best solution would be to get rid of redundant rows of J (because why would you want the covariance of redundant variables?), so that S will be strictly positive definite. Using the attached file,
Jr=licols(J')';
S=Jr*C*Jr';
S=(S+S')/2;
Now, even small floating point errors shouldn't produce negative eigenvalues since they are strictly bounded away from zero.
Related Question