Is there faster way to apply `det` function along the third dimension?

I am trying to calculate the det of many 4*4 matrix. I store the data in a matrix with shape 4*4*n. (n =4000000). I am using the for loop to get the result as below. Is there any way to accelerate the progress? Like bsxfun or arrayfun, or some simillar ideas?

v = zeros(length(m),1);
for i=1:size(m,3)
    v(i)=det(m(:,:,i));
end

[Walter Roberson]

function dets = det4(x)
   in1 = reshape(x,[],size(x,3));
   dets = in1(1,:).*in1(6,:).*in1(11,:).*in1(16,:)-in1(1,:).*in1(6,:).*in1(15,:).*in1(12,:)-in1(1,:).*in1(10,:).*in1(7,:).*in1(16,:)+in1(1,:).*in1(10,:).*in1(15,:).*in1(8,:)+in1(1,:).*in1(14,:).*in1(7,:).*in1(12,:)-in1(1,:).*in1(14,:).*in1(11,:).*in1(8,:)-in1(5,:).*in1(2,:).*in1(11,:).*in1(16,:)+in1(5,:).*in1(2,:).*in1(15,:).*in1(12,:)+in1(5,:).*in1(10,:).*in1(3,:).*in1(16,:)-in1(5,:).*in1(10,:).*in1(15,:).*in1(4,:)-in1(5,:).*in1(14,:).*in1(3,:).*in1(12,:)+in1(5,:).*in1(14,:).*in1(11,:).*in1(4,:)+in1(9,:).*in1(2,:).*in1(7,:).*in1(16,:)-in1(9,:).*in1(2,:).*in1(15,:).*in1(8,:)-in1(9,:).*in1(6,:).*in1(3,:).*in1(16,:)+in1(9,:).*in1(6,:).*in1(15,:).*in1(4,:)+in1(9,:).*in1(14,:).*in1(3,:).*in1(8,:)-in1(9,:).*in1(14,:).*in1(7,:).*in1(4,:)-in1(13,:).*in1(2,:).*in1(7,:).*in1(12,:)+in1(13,:).*in1(2,:).*in1(11,:).*in1(8,:)+in1(13,:).*in1(6,:).*in1(3,:).*in1(12,:)-in1(13,:).*in1(6,:).*in1(11,:).*in1(4,:)-in1(13,:).*in1(10,:).*in1(3,:).*in1(8,:)+in1(13,:).*in1(10,:).*in1(7,:).*in1(4,:)
end

This will vectorize the calculation.

[wei zhang]

Thanks a lot at first! But how do you get the formula? Did you have an automatic way to form this formula? When the size of x is bigger? Maybe 5, 6, or 10. It is error-prone to calculate the formula manually.

[Bruno Luong]

Formula is here

https://en.wikipedia.org/wiki/Leibniz_formula_for_determinants

As I noted earlier, this formula is only useful in theory or small size matrix, it is bad for numerical precision, and the complexity goes up like n!

For n=10, you need roughly to write down exlicitly a code with roughly

>> factorial(10)

ans =

     3628800

operations. Good luck.

[Walter Roberson]

I got the formula by using

syms x [4 4]
d = det(x)
F = matlabFunction(d, 'vars', x(:), 'filename', 'det4.m')

and then editing the results very slightly to change the input variable name from in1 to x and putting the reshape in to move the 3rd dimension to be row elements.

[Walter Roberson]

Bruno's criticism of the precision problems and the high length of the formula for increasing n, are valid criticisms. Generally speaking, it is often the case that making code faster comes at the price of making it less accurate towards the margins. Theoretical definitions that suppose infinite precision get substituted for more nuanced checks that deal with floating-point realities, and time gets saved by not making the checks to figure out what compensation is needed for each case.

If you have some time you should look at how hypot() (finding the length of a hypotenuse) has to be implemented in practice in order to maintain accuracy. Consider sqrt(A^2 + B^2) under the circumstance that A or B is smaller than sqrt(realmin) and so squaring it might underflow to 0...

[wei zhang]

@Bruno Luong Thank you for introducing the formula in wiki.

@Walter Roberson.

1. I am glad to know the matlabfunction command. It is much faster than combination of subs and syms.

2. Thank you for introducing the cause of error. It seems the error is increasing fast with the size of matrix increasing. May cause terrible results.Wow.

Answer 1

[Bruno Luong]

Similar question addressed here

I would note that Walter's solution that use recursive formal determinant formula might be fast but might be sensitive to umerical errors, that is how I was tough (for medium/large size matrix).