I don't know if the issue is from the Hessian matrix that you are reading in or from the rounding that you are manually performing in your program.
Why do you think you need to round the Sigma matrix to 0.001? If any element is less than 0.0005 in absolute value, it will be truncated to 0. Is that what you want?
Here are some IML helper functions that you can use to diagnose and fix the issue:
/* test symmetry */
start TestSymmetry(M);
diff = max(abs(M - M`));
name = strip(parentname("M"));
labl = "Max Abs Diff: " + name + " - T(" + name + ")";
print diff[L=labl];
finish;
/* make a matrix symmetric */
start MakeSymmetric(M);
return (M + M`) / 2;
finish;
The first module prints the maximum absolute difference of any unsymmetric terms. The second module returns a symmetric matrix that is close to the input matrix.
I suggest you test the input Hessian (J) for symmetry by using RUN TestSymmetry(J); If the module prints 0, then the input is not a problem/
My suspicion is that the problem in in the computation of Sigma. If J is an ill-conditioned matrix, then computing the inverse is an unstable operation and Sigma might not be symmetric. I suggest you try
S = inv(J)*K*inv(J);
run TestSymmetry(S);
Sigma= MakeSymmetric(S); /* ROUND here, if required */
I didn't use the ROUND function because I don't understand why it is needed. If required, you can add Sigma = round(Sigma, 0.001); after the call to the MakeSymmetric function.
... View more