Computing determinants of a 3D array

Let’s say I have an NxNxL array. L is typically 10^4-10^5, and N is typically 10^0-10^1.

My goal is to calculate a vector of length L, where the i-th cell contains the determinant of (:,:,i).

I currently use for-loop, as det(A) accepts 3D arrays with the last dimension being a singleton, so this code works:

   for i = 1:L

       Vec(i) = det(Mat(:,:,i));

   end

However, it seems weird that i can't implement this in a vectorial fasion. Can anybody think of any way to get rid of the for loop? Note that the format of an 3D array is pretty stiff, i.e. i can't change the input to cell array and use cellfun.

Thanks!

Yanir.

[Matt J]

Hopefully, the reason you're asking for this is not for the purpose of solving many linear systems.

[Yanir Hainick]

Hi Matt,

The context is correct, however my purpose is different. Here's the bigger picture: I attempt to solve a differential equation, subjected to boundary conditions.

Those boundary conditions may be formulated in the form A(b)*p = 0, where:

A is a matrix which depends on a parameter b, and represents the equations of boundary conditions; p is a vector of the coefficients of the differential equation solutions.

I actually search for the b's, for whom det(A(b)) = 0, and p is not the trivial one.

Yanir.

Answer 1

[Sean de Wolski]

What's wrong with the for-loop?

n=100;
pages = 1e4;
X = rand(n,n,pages);
D = zeros(pages,1);

tic;
for ii = 1:pages
  D(ii) = det(X(:,:,ii));
end
toc;
%Elapsed time is 1.805616 seconds.

Please see my comments here: FOR loops are actually fast enough...

[James Tursa]

My 2 cents to all:

The performance advantage of vectorized approaches to this is that one wishes to avoid the data copy involved with the X(:,:,ii) slices and the overhead of the loop. That being said, this data copy & loop overhead is in all likelihood swamped by the numerical calculations (and possible data copy) involved in the determinant calculation itself. So even if one were to get a vectorized one-liner to this that did not involve explicit slices, it probably wouldn't run significantly faster (if at all) than straight forward loops (the small size explicit code e.g. that Matt shows excepted).

Answer 2

[Matt J]

For small N, it would be an advantage to vectorize the determinant formula explicitly (example for for N=2 below). For larger N, maybe you could do the same thing recursively.

    %fake data
    N=2;
    L=1e5;
    Mat=rand(N,N,L);

    tic
    Vec=zeros(1,L);
    for ii=1:L

        Vec(ii)=det(Mat(:,:,ii));
    end
    toc;
    %Elapsed time is 0.194555 seconds.

    Mat=reshape(Mat,[],L);

    tic;
    Vec =Mat(1,:).*Mat(4,:) - Mat(2,:).*Mat(3,:);
    toc
    %Elapsed time is 0.000378 seconds.

[Pi Ting]

 The line
   Vec =Mat(1,:).*Mat(3,:) - Mat(2,:).*Mat(4,:);
should be
  Vec =Mat(1,:).*Mat(4,:) - Mat(2,:).*Mat(3,:)?

[Al in St. Louis]

I had to use Pi Ting's expression to get the correct answers. This is exactly what I need to do!

Answer 3

[Yanir Hainick]

Hello everybody,

Thanks for all the elaborate answers!

It seems as if there isn't a better way (within the frame of work of Matlab) to computing those L determinants, but with a for loop (arrayfun, which was the best candidate, is indeed slower...).

Computing and vectorizing the determinant explicitly will work for N=2 and even N=3, but it becomes 'exponentially' cumbersome with N ( try N=5 :(, and N=8 is rather common)

although the for loop is fast for every day purposes, the profiler shows that this is ~50% of my total run time. Since i know how to deal with the other 50%, tackling the determinant issue should prove worthy of my efforts.

Again - forgive me for my ignorance, but is there a way to transfer the array to a c++ code from within Matlab? (assume the c++ code knows how to compute an NxN determinant, and implements a for loop in the same way). will it do any good? Regards,

Yanir.