returning determinant vector - Matlab

Question

I have a matrix with 3 dimentions Y(i,j,w). I want to get a determinant vector d(w), in which each number would be the determinant of the matrix Y(:,:,w). Is there an elegant syntax for it, or I just have to use a loop? thanks

Accepted Answer

Well, first of all, you virtually NEVER truly want to compute a determinant, you just think you do. In fact, it is almost never a good thing, because determinants are so poorly scaled. Too often they are used to infer the singularity status of a matrix, which is a terrible thing to do in terms of numerical analysis.

Having stated my mini-rant against determinants in general...

OPTION 1:

Convert your 3-d array into a cell array of square matrices, with each plane of the array as one cell. mat2cell will do the trick easily and efficiently.

Next, use cellfun on the cell array. cellfun can apply a function (@det) to every cell, and then will return a vector of determinants. Is this incredibly efficient? Its probably not a huge gain over applying det in a loop, as long as you pre-allocate the vector in advance when you write a loop.

OPTION 2:

If the matrices are small, thus say 2x2 or 3x3 matrices, then expand out the multiplications for the determinant as explicit vector multiplies. I think this is not clear as I am writing it, so for a 2x2 case, where Y is 2x2xn:

```
d = Y(1,1,:).*Y(2,2,:) - Y(1,2,:).*Y(2,1,:);
```

Surely you see that this forms a vector of 2x2 determinants for every plane of the matrix Y. The 3x3 case is simple enough to write also, as six 3-way products of terms. I've not carefully checked the 3x3 case below, but it should be close.

```
d = Y(1,1,:).*Y(2,2,:).*Y(3,3,:) + ...
    Y(2,1,:).*Y(3,2,:).*Y(1,3,:) + ...
    Y(3,1,:).*Y(1,2,:).*Y(2,3,:) - ...
    Y(3,1,:).*Y(2,2,:).*Y(1,3,:) - ...
    Y(2,1,:).*Y(1,2,:).*Y(3,3,:) - ...
    Y(1,1,:).*Y(3,2,:).*Y(2,3,:);
```

As you can see, OPTION 2 will be pretty fast, and it is vectorized.

Edit: as a response to Chris, there is a SIGNIFICANT difference in the time required. Consider the time required for a set of 1e5 matrices.

```
p = 2;
n = 1e5;
Y = rand(p,p,n);

tic,
d0 = squeeze(Y(1,1,:).*Y(2,2,:) - Y(2,1,:).*Y(1,2,:));
toc

Elapsed time is 0.002141 seconds.

tic,
X = squeeze(mat2cell(Y,p,p,ones(1,n)));
d1= cellfun(@det,X);
toc

Elapsed time is 12.041883 seconds.
```

The two calls return the same values to within floating point trash.

```
std(d0-d1)
ans =
   3.8312e-17
```

A loop would be no better, in fact, surely worse. So were I to write a piece of code that would be faced with the task of generating determinants for MANY such matrices in an array, I would special case the code for 2x2 and 3x3 matrices. I might even write it out for 4x4 matrices. Yes, it is a mess to write out, but there is a big difference in the time required.

One reason is that MATLAB's det uses a call to LU, factorizing the matrix. This is better in theory than the multiplies for even medium large matrices, but for a 2x2 or a 3x3, the extra overhead is a killer. (I won't guess where the break even point falls, but one could test that easily enough.)

returning determinant vector - Matlab

Solution

Problem