# How do I test if a matrix is unitary?

96 views (last 30 days)
Bryan Acer on 9 May 2016
Answered: Pawel Tokarczuk on 9 May 2016
My current test for a unitary matrix is shown in the code, I'd also like to know if U = e^(iH) [i is the complex number] is coded correctly. Thanks!
U = exp(i*H)
Uinverse = inv(U)
UConjTran = U'
if UConjTran == Uinverse
disp('U is unitary')
else
disp('U is NOT unitary')
end

Roger Stafford on 9 May 2016
Again, e^(i*H) is not the same as exp(i*H). Check matlab's "mpower" operator.
##### 2 CommentsShow 1 older commentHide 1 older comment
Roger Stafford on 9 May 2016
To get an accurate value you should have used exp(1) instead of 2.718.
As I mentioned, you should read up on matlab's 'mpower' operator. If you have A^B where either A or B is a matrix, it has a different meaning than A.^B. For example, if
A = [1 2;3 4];
then
A^2 = A*A = [7 10;15 22]
whereas
A.^2 = [1 4;9 16];
In the case of your e^(i*H), the eigenvalues of H must be real, since H is hermitian, so the matrix e^(i*H) will have eigenvalues that all lie in the unit circle of the complex plane and that will yield a unitary matrix.

Azzi Abdelmalek on 9 May 2016
Your code is correct, but when the inverse is calculated, there maybe some errors (even if it's equal to 10^(-6)) caused by the precision allowed by any computer (software and hardware). It's recommended to test using a tolerance.
U = exp(i*H);
Uinverse = inv(U);
UConjTran = U';
tol=10^(-6);
er=abs(UConjTran-Uinverse)
if sum(er(:))<tol
disp('U is unitary')
else
disp('U is NOT unitary')
end
Bryan Acer on 9 May 2016
Great! Thank you, U = exp(i*H) was also not providing the correct values as that is not equal to the same expression when eulers number "e" is plugged in explicitly as 2.718. Roger Stafford addresses that above.

Pawel Tokarczuk on 9 May 2016
Your notation suggests that what you need is the matrix exponential:
U = expm(i*H);
This is not the same thing (in general) as the element-wise exponential:
V = exp(i*H);
Try it and be convinced. Anyway, the test for a unitary matrix is:
U*U' = U'*U = I, to some floating-point tolerance, where I is the unit matrix.
Finally, bear in mind that the evolution operator U takes on a more complicated (time-ordered) form when Hamiltonians H evaluated at different times do not commute.
See: Merzbacher ("Quantum Mechanics"), Constantinescu and Magyari ("Problems in Quantum Mechanics"), etc.

Kevin on 9 May 2016
Another way to test if a matrix is unitary is to check if (U * U') == square identity matrix (with some threshold). This way, you can avoid matrix inversion.
Another way that I can think of is to use SVD (Matlab function svd).