Computing the determinant of a block matrix

I am having trouble using a well-known formula for computing the determinant of a block matrix. That is, if

in which A and D are square matrices and exists then

.

The code snippet below should achieve this but it returns two different results. I'm sure it's a very trivial mistake I've made but I've probably been staring at it too long now to find it.

As a follow up, does anyone know if there is a suitably adjusted version of this formula for the case when it is B and C which are the square matrices?

function [detM,detBlock] = BlockDeterminant()

    M = magic(4);
    detM = det(M);
    
    A = M(1:3,1:3);
    B = M(1:3,4);
    C = M(4,1:3);
    D = M(4,4);
    detBlock = det(A)*det(D-C*inv(A)*B);

end

Answer 1

[KSSV]

A = rand(2) ; 
B = rand(2) ;
C = rand(2) ;
D = rand(2) ;
M = [A B; C D] ;
det(M)

ans = -0.0221

det(A)*det(D-C*inv(A)*B)

ans = -0.0221

In your case, D is a scalar.

[Abderrahim. B]

Also your code:

A = rand(1) 

A = 0.9162

B = rand(1) ;
C = rand(1) ;
D = rand(1) ;
M = [A B; C D] ;
det(M)

ans = 0.0193

det(A)*det(D-C*inv(A)*B)

ans = 0.0193

[Bruno Luong]

"In your case, D is a scalar. "

It doesn't matter. The formula is still correct.

[RickyBoy]

Sorry @Bruno Luong, is your comment here about the case where B and C become the square matrices? In my case there will always be one 1x1 matrix, i.e. scalar, which in any case offers simplifications.

[RickyBoy]

Actually, I think I may have answered this question myself... if B and C are the square matrices then M is no longer block diagonal, so the conditions for the formula are not met.

[Bruno Luong]

The formula is correct even for m ~= n, where m is the size of the first diagobal block and n is the second block;

m=10; n=4;
A=rand(m,m); B=rand(m,n); C=rand(n,m); D=rand(n,n);
det([A B; C D])

ans = -0.0163

det(A)*det(D-C*inv(A)*B) 

ans = -0.0163

[RickyBoy]

Sure, but in the ideal case for my purposes, I would have

A: m x 1

B: m x m

C: 1 x 1

D: 1 x m

In this case M = [A,B;C,D] is not block-diagonal in the strict sense so I question whether the formula no longer applies or if there is any related formula for this specification.

[Bruno Luong]

You reverse the columns

M = M(:,end:-1:1)

it will become diagonal block again and the determinant of both sides are multiplied by (-1)^(m+n) so the formulat becomes

m=10; n=3;
A=rand(m,n); B=rand(m,m); C=rand(n,n); D=rand(n,m);

det([A B; C D])

ans = -0.0039

det(B)*det(C-D*inv(B)*A) 

ans = -0.0039

[RickyBoy]

Can we verify if the multiplication factor should be (-1)^(m+n-1) ?

[Bruno Luong]

There is no multiplication factor as they cancel out on both sides. The formula is what I write in the code.

det(M) = det(B)*det(C-D*inv(B)*A) 

[RickyBoy]

OK, I may have communicated badly ... I'll post the code I currently have shortly. Should be clearer what I'm trying to express.

[RickyBoy]

@Bruno Luong Here is the code snippet:

function [detM,detBlock,detBlockAlt] = BlockDeterminant2()
    
    m=4;
    n=2;

    A = rand(m,n);
    B = rand(m,m);
    C = rand(n,n);
    D = rand(n,m);

    detM = det([A,B;C,D]);
    
    detBlock = (-1)^(m+n-1)*det(B)*det(C-D*inv(B)*A);
    detBlockAlt = (-1)^(m+n-1)*det(C)*det(B-A*inv(C)*D);

[Bruno Luong]

detBlock = (-1)^(m*n)*det(B)*det(C-D*inv(B)*A)
detBlockAlt = (-1)^(m*n)*det(C)*det(B-A*inv(C)*D)

[RickyBoy]

Aha!

OK, so for a bit of enlightenment, @Bruno Luong, let me demonstrate how this is so helpful...

In my setup (and it happens to be that . Therefore using detBlockAlt, we can achieve

.

Which is a much simpler form of the determinant than a general formulation will result in.

Thank you so much for your engagement!

[Bruno Luong]

Good. I don't get why contribution is so helpful, but if you say so...

Answer 2

[Walter Roberson]

Just round-off error. You are calculating two different order of operations and that affects the results .

If you use M = sym(magic(4)) you will get exact zero for each