Computing determinants of a 3D array

12 次查看（过去 30 天）
Yanir Hainick
Yanir Hainick 2012-11-26
编辑： Matt J 2020-4-8
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.
  2 个评论
Matt J
Matt J 2012-11-26
Hopefully, the reason you're asking for this is not for the purpose of solving many linear systems.
Yanir Hainick
Yanir Hainick 2012-11-29
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Sean de Wolski
Sean de Wolski 2012-11-26
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...
  1 个评论
James Tursa
James Tursa 2012-11-27
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).

请先登录，再进行评论。

更多回答（2 个）

Matt J
Matt J 2012-11-26
编辑：Matt J 2020-4-8
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.
  2 个评论
Pi Ting
Pi Ting 2017-11-8
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
Al in St. Louis 2020-4-8
I had to use Pi Ting's expression to get the correct answers. This is exactly what I need to do!

请先登录，再进行评论。

Yanir Hainick
Yanir Hainick 2012-11-29
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.

类别

在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by