page-wise matrix determinant or eigenvalues

42 次查看（过去 30 天）
Henry Brinkerhoff
Henry Brinkerhoff 2023-3-10
编辑： Walter Roberson 2023-10-24
I'm really loving the vectorized improvement in a lot of my code from incorporating Matlab's new page-wise matrix operations, since I'm regularly running code where I need to do, for example, inverses of each of the M*N number of 3x3 matrices in a 3x3xMxN array.
I'd really like to also be able to take a determinant of each matrix in this same fashion, or at least get the eigenvalues so I can quickly compute the determinant as their product. pagesvd.m has been implemented, and it seems straightforward to similarly implement a pageeigs.m or pagedet.m in a vectorized fashion, but unfortunately it doesn't seem like these functions exist.
Has anybody else encountered this issue and found a workaround for a quick, vectorized way (without "for" loops) to implement a page-wise determinant?
Thanks!

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

采纳的回答

Heiko Weichelt
Heiko Weichelt 2023-3-16
R2023a, released March 15th 2023, ships PAGEEIG as new function in the page* family:
https://www.mathworks.com/help/matlab/release-notes.html?rntext=pageeig&startrelease=R2023a&endrelease=R2023a&groupby=release&sortby=descending&searchHighlight=pageeig
Try it out and let us know if it meets your requirements. We're happy to hear about more use-cases and workflows for further page* functions.
  2 个评论
David Ober
David Ober 2023-3-28
I use determinant to solve items in the decomposition of the projection matrix for over 100k 3x3 matricies. So I wrote my own det33 (for the 3x3 that is vectorized) using the Rule of Sarrus. My custom version is faster (0.5X the time) than the following one liner for 100k matricies.
fvdet = prod(pageeig(mmM,'vector'),1);
a Page-wise "diag" would be very helpful to check the signs of the decomposition. Thankfully, since the decomposition is a 3x3, I can write my own vectorized "diag33".

请先登录，再进行评论。

更多回答（3 个）

Henry Brinkerhoff
Henry Brinkerhoff 2023-3-10
To anyone reading this in the future: I found something that works for my specific case, using symmetric and positive-definite matrices:
The product of the diagonal elements of the diagonal matrix in the singular value decomposition (which exists as a page-wise function) is still the determinant of the matrix, if the input matrix is positive definite! So you can use that function to generate the outputs you need to quickly calculate the determinant of each matrix in the array.
  2 个评论
John D'Errico
John D'Errico 2023-3-10
Ran in:
Um, close. In fact, the product of the singular values will be the absolute value of the determinant. And that will apply regardless of whether your matrix is positive definite. For example:
A = randn(4)
A = 4×4
-1.7274 0.2800 -0.3666 0.5376 -0.1166 0.6912 0.7849 0.0446 -0.8638 0.2588 -0.2872 1.0908 -0.0641 -0.0214 -0.4188 -2.0971
A is clearly never going to be positive definite.
det(A)
ans = -0.9566
prod(svd(A))
ans = 0.9566
But, as I said, the product of the singular values is the same, to within a factor of -1
Henry Brinkerhoff
Henry Brinkerhoff 2023-3-10
Right- which is why this only works for positive definite matrices, which have determinant greater than 0- absolute value of a positive number is of course the same number!

请先登录，再进行评论。

John D'Errico
John D'Errico 2023-3-10
Ran in:
@Henry Brinkerhoff seems to have found a semi-viable solution, in the form of pagesvd. It will be valid, within a factor of -1. However, that are other methods around that would work as well, and maybe better.
For example, the determinant can also be gained from the product of the diagonal elements of the U factor in an LU decomposition, times the number of column swaps involved in the permutation matrix P. For example:
A = rand(4);
[L,U,P] = lu(A)
L = 4×4
1.0000 0 0 0 0.0826 1.0000 0 0 0.0628 0.2774 1.0000 0 0.3842 -0.0464 0.9820 1.0000
U = 4×4
0.9228 0.6367 0.8450 0.6547 0 0.8852 0.5847 0.1327 0 0 0.6889 0.1395 0 0 0 0.0161
P = 4×4
0 0 0 1 0 1 0 0 1 0 0 0 0 0 1 0
[det(A),prod(diag(U))]
ans = 1×2
0.0091 0.0091
And see that we can convert the matrix P into an identity matrix by swapping one pair of the columns. So there is one column swap, and therefore one factor of -1. If you don't care about the sign of the determinant, you can ignore the permutation matrix P.
But a nice thing is, if we could convert the matrix A into a block diagonal matrix., then we could compute the determinants of MANY pages quickly.
For example,
A = rand(5,5,1000);
C = mat2cell(A,5,5,ones(1,1000));
[L,U] = lu(blkdiag(C{:}));
D = prod(reshape(diag(U),5,1000),1)
D = 1×1000
0.1952 0.0029 0.0551 -0.1405 -0.0759 0.0420 -0.0197 0.0186 0.0042 -0.0554 0.0081 0.0220 -0.0723 -0.0025 0.2211 0.0389 -0.1037 0.0396 0.0198 -0.0928 -0.1030 0.1175 0.0091 0.0426 0.0561 -0.0259 0.0892 0.0777 -0.0657 -0.0994
det(A(:,:,1))
ans = -0.1952
det(A(:,:,2))
ans = 0.0029
Again, to within an arbitrary factor of -1, the two are the same. Unfortunately, there is no paged version of the LU. The only thing missing is the fact that we need the matrix to be a sparse block diagonal matrix. Then it would be very efficient. We can fix that too.
tic
C = sparse(reshape(A,5,5*1000));
S = mat2cell(C,5,repmat(5,1,1000));
[L,U] = lu(blkdiag(S{:}));D1 = prod(reshape(diag(U),5,1000),1);
toc
Elapsed time is 0.041744 seconds.
tic,D2 = squeeze(prod(pagesvd(A),1));toc
Elapsed time is 0.019279 seconds.
Surprisingly, PAGESVD is still faster than the sparse block diagonal LU. I guess we need a paged LU to be competitive.
lala
lala 2023-10-24
编辑：Walter Roberson 2023-10-24
Ran in:
det with for-loop is fatest.
format long g
[tmr,tf] = testpagemdet(100000)
tmr = 1×3
1.0e+00 * 7e-06 1.3e-05 1.7e-05
tf = logical
1
function [tmr,tf] = testpagemdet(times)
A = rand(3,3,8);
tmr = zeros(1,3);
for i = 1:times
t = tic;
d1 = pagemdet(A);
tmr(1) = toc(t);
t = tic;
d2 = det33(A);
tmr(2) = toc(t);
t = tic;
d3 = real(prod(pageeig(A,'vector'),1));
tmr(3) = toc(t);
end
tol = 1e-6;
tf = all(abs(d1-d2)<tol & abs(d1-d3)<tol);
end
function d = pagemdet(A)
sz = size(A);
d = zeros([1,1,sz(3:end)]);
nm = prod(sz(3:end));
for i = 1:nm
d(i) = det(A(:,:,i));
end
end
function d = det33(A)
d = A(1,1,:).*A(2,2,:).*A(3,3,:)...
+A(1,2,:).*A(2,3,:).*A(3,1,:)...
+A(1,3,:).*A(2,1,:).*A(3,2,:)...
-A(1,3,:).*A(2,2,:).*A(3,1,:)...
-A(1,2,:).*A(2,1,:).*A(3,3,:)...
-A(1,1,:).*A(2,3,:).*A(3,2,:);
end

类别

在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by