0

In the following code, I followed a procedure to create a random positive definite matrix P.

At first, I created a singular value decomposition [U,S,V] of a random array A and I am trying to verify that in fact U'*U==U*U'=I (where I is the identity matrix) which is known from theory. The problem is that if I check the contents of the matrices on my own I can verify that, but Matlab produces a logical matrix which does not verifies that since zero is represented as -0.000 or 0.0000 so the result is 1 only if signs match. Why is that?

But the bigger problem arises in a few lines below this, where the positive definite (all its eigenvalues are positive) matrix P is produced and I simply want to check that P=P'. By clicking at x and y I can see that their contents are exactly the same but Matlab can't verify this either. I don't understand why this is happening since in this case we are talking about the same variable P here which is simply transposed.

Here is the code:

%Dimension of the problem
n = 100;

%Random matrix A
A = rand(n, n)*10 - 5;

%Singular value decomposition
[U, S, V] = svd(A);

%Verification that U*U'=U'*U=I
U'*U == U*U'

%Minimum eigenvalue of S
l_min = min(diag(S));

%Maximum eigenvalue of S
l_max = max(diag(S));

%The rest of the eigenvalues are distributed randomly between the minimum
%and the maximum value
z = l_min + (l_max - l_min)*rand(n - 2, 1);

%The vector of the eigenvalues
eig_p = [l_min; l_max; z];

%The Lambda diagonal matrix
Lambda = diag(eig_p);

%The P matrix 
P = U*Lambda*U';

%Verification that P is positive definite
all(eig(P) > 0)

%Verification that P=P'
x=P;
y=P';
x==y
mgus
  • 808
  • 4
  • 17
  • 39

1 Answers1

2

Your problem is not that 0 is represented as -0.000 or 0.0000 (you can verify this by running (-0) == 0), but rather that the values are actually just really small (relative to the other numbers in the matrix). Operations like what you are doing will always have small errors from the way that floating point numbers are represented and manipulated.

Instead of verifying that P=P', try verifying that abs(P-P') < 0.000000001, or whatever threshold is appropriate for you given the application.

zelanix
  • 3,326
  • 1
  • 25
  • 35