Reputation: 555
I am trying to write a method that checks if a matrix is orthogonal and return TRUE if it is or FALSE if it isn't My problem is that my isequal() is not working how I want it to. Basically I can do the check in two ways based on the two formulas:
ONE way is check to see if the transpose of matrix R is equal to the inverse of matrix R. If they are equal then it is orthogonal. (R'=inv(R))
ANOTHER way is to check and see if matrix R times the transpose of matrix R equals the Identity matrix of R. (R'R=I) If yes then the matrix is orthogonal. I have most been using isequal() but it keeps yielding false. Can someone look at my code and tell me why this would be so?
I use Z=orth(randn(3,3)) to generate random orthogonal matrix and i call my method isortho(Z)
function R = isortho(r)
%isortho(R), which returns true if R is orthogonal matrix, otherwise returns false.
if ismatrix(r) && size(r,1)==size(r,2) %checks if input is square matrix
'------'
trans=transpose(r)
inverted=inv(r)
isequal(trans,inverted)
trans==inverted
isequal(transpose(r),inv(r)) %METHOD ONE
i=size(r,1);
I=eye(i) %creating Identity matrix based on size of r
r*transpose(r)
r*transpose(r)==I %METHOD TWO
%check if transpose of r is times inverse of r equals Identity matrix of r
if (r*transpose(r)==I)
R= 'True';
else
R= 'False';
end
end
end
this is my output:
>> isortho(Z)
ans =
------
trans =
-0.2579 -0.7291 -0.6339
0.8740 0.1035 -0.4747
0.4117 -0.6765 0.6106
inverted =
-0.2579 -0.7291 -0.6339
0.8740 0.1035 -0.4747
0.4117 -0.6765 0.6106
ans = ////isequal(trans,inverted) which yielded 0 false
0
ans = ////trans==inverted
0 1 0
1 0 0
0 1 1
ans = ////isequal(transpose(r),inv(r))
0
I =
1 0 0
0 1 0
0 0 1
ans =
1.0000 0 0.0000
0 1.0000 0.0000
0.0000 0.0000 1.0000
ans =
1 1 0
1 1 0
0 0 1
ans =
False
>>
could someone help me fix this or tell my why the isequal() is failing when matrix inverted and trans appear to be the same?
Upvotes: 2
Views: 3255
Reputation: 8459
As stated in the comments, you are running into computer precision issues. For more detail see Why is 24.0000 not equal to 24.0000 in MATLAB? and http://matlabgeeks.com/tips-tutorials/floating-point-comparisons-in-matlab/. This is not a Matlab specific thing, it's a computer thing, and you just have to deal with it.
In your case, you are trying to see whether two things are equal, but the two things are the result of a lot of floating point operations. So they will virtually never be exactly the same, but should always be very close. So, set a tolerance, say 1e-12, and say that the two things are equal if some measure of their difference is below that tolerance, e.g.:
norm(r.'-inv(r))<tol
Which finds the 2-norm of the difference between the two matrices, and then if it is less that tol
, this will evaluate to 1, or true.
If I set tol=1e-12
, then everything works well. If I set tol=1e-15
, everything works well. But if I set tol=1e-16
, then everything stops working! This is because the amount of computer precsion error is larger than 1e-16, so the answer to norm(r.'-inv(r))
cannot be accurate to that tolerance. The smallest amount Matlab can distinguish between on my computer is roughly 2.2x10^(-16), so you have to ensure that you tolerance is set well above this value. Setting tol
too large will, of course, mean you say some non-orthogonal matrices are orthogonal, but I would not expect tol=1e-14
to give you any significant issues.
Upvotes: 4