Not sure best way to code orthogonal diagonalization

Just having an issue with highlighted part, will attach rest of code below.

function []=symmetric(A)

if size(A,1)~=size(A,2)

disp('Matrix is not symmetric')

return

else

if all(abs(A-A')<.000001)

[P,D]=eig(A);

if all(abs(inv(P)-P')<.00001)

disp('P:')

disp(P)

disp('D:')

disp(D)

disp('P is orthogonal matrix')

else

disp('P:')

disp(P)

disp('D:')

disp(D)

disp('P is orthogonal matrix')

end

else

disp('Matrix is not symmetric')

end

end

end

[John D'Errico]

Note that this test for symmetry is generally inadequate:

if all(abs(A-A')<.000001)

A simple counter-example is:

A = [0 -1e-7;1e-7,0];

which is as clearly a non-symmetric matrix as you could imagine.

Next, this is not the test you were required to perform:

if all(abs(inv(P)-P')<.00001)

[David Goodmanson]

Hi Thomas,

The A matrices you have been given have integer elements so with the 1e-6 test it's perfectly clear whether or not they are symmetric. But in general you have to be more wary.

To expand on the point John is making, suppose that [P D] = eig(A). If P is a real orthogonal matrix, and with P' being the transpose of P, then by definition P'*P = the identity matrix. Now suppose A passes the 1e-6 test but is not symmetric:

A = [1   3e-07
    8e-07    2]

Here the eigenvalues of A are very nearly equal to 1 and 2 and are well separated. Then

format short g
[P D] = eig(A)
P'*P
ans  =
            1       -5e-07
       -5e-07            1

so the P matrix misses being orthogonal by about the same amount that A misses being symmetrical. So you might say, ok, the 1e-6 test makes sense in that regard. But now suppose the two eigenvalues are nearly identical:

A = [1   3e-07
    8e-07    1]
format short g
[P D] = eig(A)
P'*P
ans =
            1      0.45455
      0.45455            1

In this case, even though the 1e-6 test is still in effect, P is not even close to being an orthogonal matrix. Hence the inadequacy of the test on its own.

[Thomas Stowell]

Im not entirely sure of a good way to test for symmetry and orthogonality simultaneously. I believe I understand why my code is inadequate but I'm new to matlab and am not so good with code. Do you have suggestions?

[Thomas Stowell]

Is there any suggested code or is this thread primarily to point out how bad other code is?

[David Goodmanson]

Hi Thomas,

That's a good point. Because of your other question I had a feeling that the stuff here was not really helping.

Since the entries of A in this case are integers you can do a check for exact equality, A - A' == 0. The other tests are going to take tolerance checks, such as all(A*P-P*D < tolerance). As far as the acceptable tolerance, it looks like the instructor has provided that with closetozeroroundoff (although I have no idea how it works). For orthogonality, you can have all(inv(P) -P' < tolerance) as you are doing. Probably better, especially for large matrices, is not doing the inverse. For an orthogonal matrix P*P' = eye(size(P)) so you can check all(P*P'-eye(size(P))< tolerance).

Answer 1

[David Hill]

For homework, we like to help you figure it out yourself. The code below might help and you could likely find a better way.

function a=symmetric(A)
   if issymmetric(A)
       [P,D]=eig(A);
       if closetozeroroundoff(P*D,A*P)&&closetozeroroundoff(inv(P),P')
           a='orthogonal diagonal';
       else
           a='symmetric';
       end
   else
       a='nonsymmetric';
   end
end
function yn=closetozeroroundoff(a,b)
yn=0;
    if norm(a-b)<1e-14
        yn=1;
    end
end