Determinant of a unitary matrix

10 visualizaciones (últimos 30 días)
Bruno Luong
Bruno Luong el 23 de Nov. de 2020
Editada: Bruno Luong el 25 de Nov. de 2020
Is there any other (better) way to compute the determinant of the unitay matrix beside det (that calls lu factorization)?
>> [U,~]=qr(rand(5)+1i*rand(5))
U =
-0.4354 - 0.1474i -0.2285 - 0.0527i -0.0673 - 0.1461i 0.5989 + 0.0097i 0.3444 - 0.4800i
-0.0104 - 0.3044i -0.1395 - 0.1222i -0.6371 + 0.1020i -0.4880 - 0.2927i 0.3406 - 0.1294i
-0.1929 - 0.4992i -0.0791 - 0.2610i -0.2843 + 0.1059i 0.2578 + 0.0370i -0.6394 + 0.2658i
-0.5246 - 0.3650i 0.4425 + 0.2340i 0.2840 - 0.3511i -0.3396 - 0.1282i -0.0556 - 0.0476i
-0.0303 - 0.0159i -0.6434 - 0.4143i 0.4108 - 0.3052i -0.3370 - 0.0652i -0.1474 - 0.1081i
>> det(U)
ans =
-0.8370 - 0.5472i
It seems not but I could miss some obscure algebra properties.
  5 comentarios
Christine Tobler
Christine Tobler el 24 de Nov. de 2020
BTW, I'd be interested in why you need to know the determinant of this unitary matrix. If it's computed through QR, do you also need the determinant of the R factor?
Bruno Luong
Bruno Luong el 24 de Nov. de 2020
It's a subtask when I want to generate a random matrix in SU(n).
So I generate matrix in U(n), compute its determinant and divide one of the vector by the determinant.
I just submit the FEX earlier today.

Iniciar sesión para comentar.

Respuesta aceptada

John D'Errico
John D'Errico el 23 de Nov. de 2020
Editada: John D'Errico el 23 de Nov. de 2020
Gosh. I wonder, if there were really much better ways to compute the determinant, they might have used it? ;-)
There are no special properties you can use, at least none I can think of. You could compute the product of the eigenvalues, but eig should generally be slower than lu.
[U,~]=qr(rand(5)+1i*rand(5));
det(U)
ans = -0.9789 - 0.2042i
If the matrix is real, then the determinant would be 1. But for the complex case, all you can know is the magnitude of the determinant should be 1. I think that is all you get from the matrix being unitary.
abs(det(U))
ans = 1.0000
timeit(@() det(U))
ans = 8.4200e-06
prod(eig(U))
ans = -0.9789 - 0.2042i
timeit(@() prod(eig(U)))
ans = 1.3420e-05
And of course, you could use more foolish ways, like decomposing it as an expansion by minors. That would be an exponentially bad idea. Actually, "factorially" might be a better word, as I recall.
But I don't think you can do much better than the lu scheme.
  13 comentarios
Paul
Paul el 25 de Nov. de 2020
It's self evident that the sum of the angles is real and that exp(1i*anglething) should have norm 1. But it doesn't always have norm 1 because of numerical inaccuracies. For example:
>> anglething=-3.817809607026706e+00;
>> d=exp(1i*anglething);
>> abs(d)
ans =
9.999999999999999e-01
To force abs(d) == 1 then normalize:
>> d=d/abs(d);
>> abs(d)
ans =
1
That last operation ( /abs(d)) is what I would call "normalization," which was not included in the original formula that started this subthread. I guess we just have different ideas of what normalization means.
Bruno Luong
Bruno Luong el 25 de Nov. de 2020
Editada: Bruno Luong el 25 de Nov. de 2020
I simply claim your method is "equivalent" to a normalization, but perform in a complicate way. I never claim it gives the exact result as the normalization.
Thanks for teaching me the fact that if one compute numerically 2 thing differently it leads to some small numerical error, and sin(x)^2+cos(c)^2 is not always == 1 numerically. (x == anglething). That not's really a new knowledge? or is it?
The angle(...) takes atan2 of imaginary and real part of lambda, then exp(1i*..) takes the cos() and sin() then for the complex number. That's a complicated way to normalization it to me, and you are free to tell it's not. Granted, in between you also use the fact that exp of the sum is a product of exp.
Last example
>> U=randn(5)+1i*randi(5) % I purposingly use non unitary matrix here
U =
0.7269 + 4.0000i -1.1471 + 4.0000i 0.3252 + 4.0000i -0.2414 + 4.0000i -0.1649 + 4.0000i
-0.3034 + 4.0000i -1.0689 + 4.0000i -0.7549 + 4.0000i 0.3192 + 4.0000i 0.6277 + 4.0000i
0.2939 + 4.0000i -0.8095 + 4.0000i 1.3703 + 4.0000i 0.3129 + 4.0000i 1.0933 + 4.0000i
-0.7873 + 4.0000i -2.9443 + 4.0000i -1.7115 + 4.0000i -0.8649 + 4.0000i 1.1093 + 4.0000i
0.8884 + 4.0000i 1.4384 + 4.0000i -0.1022 + 4.0000i -0.0301 + 4.0000i -0.8637 + 4.0000i
>> exp(1i*sum(angle(eig(U))))
ans =
0.0771 + 0.9970i
>> d=det(U); d=d/abs(d)
d =
0.0771 + 0.9970i
>> d=prod(eig(U)); d=d/abs(d)
d =
0.0771 + 0.9970i
Those three methods give the exact same value numerically? I wouldn't dare to claim that, but close enough are to me they are equivalent and I call them all normalization of determinant.
OK I get your formula. It's totally revolutionary (well not really, sorry because if I have to select one of the three for my implementation, I'll pick the second method).
I have to move on.

Iniciar sesión para comentar.

Más respuestas (1)

Bruno Luong
Bruno Luong el 24 de Nov. de 2020
Editada: Bruno Luong el 25 de Nov. de 2020
Here is one way of computing determinant of a complex matrix U, based on Q-less qr
n = size(U,1);
d = 1;
for k=1:n
i = k:n;
u = U(k,i);
a = norm(u)*u(1)/abs(u(1));
u(1) = u(1)+a;
j = k+1:n;
v = U(j,i)*u';
v = v*(2/sum(u.*conj(u)));
U(j,j) = U(j,j) - v*u(2:end);
d = d*a;
end
% here is d is det(U), original U is detroyed
With such matlab implementation it expected to be slower than det(U). But can be a base for C-implementation.

Categorías

Más información sobre Sparse Matrices en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by