How do I test if a matrix is unitary?

71 views (last 30 days)
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

Accepted Answer

Roger Stafford
Roger Stafford on 9 May 2016
Again, e^(i*H) is not the same as exp(i*H). Check matlab's "mpower" operator.
  2 Comments
Bryan Acer
Bryan Acer on 9 May 2016
Thank you, this was the error that caused the failure to recognize U as a unitary matrix. Replacing "U = exp(i*H)" with "U = 2.718 ^ (i * H)" and then checking U as unitary with the conjugate transpose and inverse does prove that it is unitary. Although this is a crude fix. Can you explain why exp(i*H) is not equal to 2.718^(i*H)? Matlab states, "exp(X) is the exponential of the elements of X, e to the X"
Roger Stafford
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.

Sign in to comment.

More Answers (3)

Azzi Abdelmalek
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
  1 Comment
Bryan Acer
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.

Sign in to comment.


Pawel Tokarczuk
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
Kevin on 9 May 2016
Maybe you have already thought about this but decided not to use this method.
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).

Categories

Find more on Eigenvalues & Eigenvectors in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!