page-wise matrix determinant or eigenvalues

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!

Answer 1

[Heiko Weichelt]

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.

[Henry Brinkerhoff]

Great, thank you!

[David Ober]

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".

Answer 2

[Henry Brinkerhoff]

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.

[John D'Errico]

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]

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!

Answer 3

[John D'Errico]

@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.

Answer 4

[lala]

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