Determinant of a unitary matrix

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.

[Christine Tobler]

I can't think of a way to use the fact that the input is unitary in DET. The usual tricks for computing the determinant would be to factorize into triagular matrices (as DET does with LU), and there's nothing particularly useful about a unitary matrix there. Another way would be to split the matrix into blocks and use Schur-complement, but since the blocks of a unitary matrix aren't unitary, I don't think this can lead far.

Well, at least you can compute the absolute value returned by DET to check how accurate the result is.

[Paul]

What is the criteria for "better"? Accuracy? Fewer computations?

[Bruno Luong]

Mainly faster with similar accuracy.

But any idea is welcome; for instant if there is some algorithm that is slow in matlab but requires less flops theoretically I'll take it.

[Christine Tobler]

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]

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.

Answer 1

[John D'Errico]

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.

[Bruno Luong]

Alright if John and Christine couldn't find alternative method the it should not be possible. ;-)

BTW I end up using Laplace recursive method alternatively for dimension n < 4. I hope the precision is OK as we deal with unitary matrix, and there is no issue from substraction of 2 large quantities. Recursive formula is vectorizable.

My tests show that the cpu time is about the same for n=4, and faster when n smaller and vice versa.

Also currently for input real matrix A, if U is obtained by

[U,~] = qr(A);

then U has determinant = (-1)^(n-1).

However I can't guess the determinant of U if A is complex.

[Christine Tobler]

That's quite interesting about the QR: it would make sense since U here is computed by Householder transformations, and I think it would be possible to compute the determinant more efficiently given Householder transformations as opposed to just an arbitrary unitary matrix.

I tried a bit more. So as you said, for QR on real the determinant is always determined by just the matrix size:

>> bins = []; for n=1:10, for ii=1:1e4, rng(ii); [U, ~] = qr(randn(n)); bins(n, ii) = det(U); end; end; [(1:10)', mean(bins, 2)]
ans =
     1     1
     2    -1
     3     1
     4    -1
     5     1
     6    -1
     7     1
     8    -1
     9     1
    10    -1

For ORTH, that's not the case, but we still see a bias towards the same sign from QR (this makes sense since ORTH uses SVD, and that uses Householder transformations as a first step):

>> bins = []; for n=1:10, for ii=1:1e4, rng(ii); U = orth(randn(n)); bins(n, ii) = det(U); end; end; [(1:10)', mean(bins, 2)]
ans =
    1.0000   -0.0012
    2.0000   -1.0000
    3.0000    1.0000
    4.0000   -0.7414
    5.0000    0.7882
    6.0000   -0.8346
    7.0000    0.7786
    8.0000   -0.6912
    9.0000    0.6648
   10.0000   -0.6244

For QR on complex, there's no clear trend (here's a plot of the different determinants from 1e4 random complex matrices fed into QR:

>> bins = []; for n=1:10, for ii=1:1e4, rng(ii); [U, ~] = qr(randn(n)+1i*randn(n)); bins(n, ii) = det(U); end; end; [(1:10)', mean(bins, 2)]
ans =
   1.0000 + 0.0000i  -0.6401 + 0.0026i
   2.0000 + 0.0000i   0.5001 + 0.0058i
   3.0000 + 0.0000i  -0.4136 + 0.0049i
   4.0000 + 0.0000i   0.3527 - 0.0102i
   5.0000 + 0.0000i  -0.3066 + 0.0011i
   6.0000 + 0.0000i   0.2798 - 0.0065i
   7.0000 + 0.0000i  -0.2438 + 0.0096i
   8.0000 + 0.0000i   0.2314 - 0.0009i
   9.0000 + 0.0000i  -0.2262 + 0.0041i
  10.0000 + 0.0000i   0.2035 + 0.0046i

Apart from the scalar case, there's a bit of bias on the real part, but nothing you could actually use directly.

That being said, if your complex matrices are constructed by doing Householder transformations on a matrix, I think you could likely design a way of scaling those Householder transformations so the determinant is always 1 or -1, too.

Is that where your unitary matrices are coming from?

[Bruno Luong]

"Is that where your unitary matrices are coming from?"

Actually I can use this file exhange to perform qr factorization, and U matrix always has determinant of (-1)^n, even for complex entry.

This package however becomes slow for large n.

[Paul]

I'm finding the following to be more accurate, though slower, for unitary matrices developed as in your example code:

[U,~]=qr(rand(5)+1i*rand(5));
d2=exp(1i*sum(angle(eig(U))));

By "more accurate" I mean that:

abs(1-abs(d2)) < abs(1-abs(det(U)))

[Bruno Luong]

Well we all know

abs(det(U)) = 1;

to "improve" accuracy, we just need to renormalize it

d = det(U);
detu = d/abs(d);

Essentially your formula does just that in a complicated way IMO.

[Christine Tobler]

So the background here is that we're using a LAPACK function, which uses different scaling conventions in the real and complex cases:

• In the real case Householder transformations have a factor tau that is always represented as nonnegative. This implies that the diagonal values of R may be positive or negative.
• In the complex case Householder transformations have a factor tau that always has nonnegative real part but is usually complex. This is chosen so that the diagonal values of R are always real and nonnegative.

So the choice for the real case is useful to you here, and the one for the complex case isn't. I'd assume you're using the same choices of scalings in the real and complex case in the MultipleQR package.

For performance, possibly you could write a mex file that calls LAPACK's QR factorization and returns the Householder vectors and scalings, and then modify those scalings and call the second LAPACK function that constructs Q from those Householder vectors.

Caveat: Changing the scalings while computing these Householder vectors would definitely work, but doing so after the fact might result in a wrong result.

[Bruno Luong]

Neat idea this lapack calling.

However this goes a little overboard for the time I can dedicate to create a file exchange for the moment. But yes I definitively keep the idea in mind. Thanks.

[Paul]

I don't think it's a normalization. Just a different way to compute prod(eig(U)). It doesn't alwas yield a result == 1 as would be the case with normalization.

[Bruno Luong]

@Paul

exp(1i*sum(angle(eig(U)))) = prod(exp(1i*angle(eig(U)))

Since for all lambda in C

exp(1i*angle(lambda)) = lambda / abs(lambda)

So if we take lambda = eig(U), and recall det(U):=prod(lambda)

prod(exp(1i*angle(eig(U))) = 
prod(lambda / abs(lambda)) = 
prod(lambda) / abs(prod(lambda)) =
det(U) / abs(det(U))

Your formulation is just a normalization on top of compute det from eig

lambda =eig(U)
d = prod(lambda)
d = d/ abs(d) % this is unit complex, not necessary 1 

[Paul]

I should have said "... doesn't always yield a result such that abs(result) == 1 ..."

I still don't see how taking advantage of the fact that all eigenvalues of U (theoretically) have a unit norm to compute prod(eig(U)) is the same as normalizing after the fact. In other words

prod(eig(U))=
prod(abs(eig(U)).*exp(1i*angle(eig(U))))=
prod(exp(1i*angle(eig(U))))=
exp(1i*sum(angle(eig(U)))))

If prod(eig(U)) isn't normalized than neither is exp(1i*sum(angle(eig(U)))).

[Bruno Luong]

anglething = sum(angle(eig(U)) is REAL (because angle is REAL) so

exp(1i*anglething) MUST have modulus 1 (or norm 1 if you like). Please check this out

Your formula does nothing but normalizes the product of eigen values.

I show by A+B it's the same thing, I can't help you if you don't believe the math.

[Paul]

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]

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.

Answer 2

[Bruno Luong]

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.