Developer Forum

I had a strange result, using the matrix inversion lema. I have a
simple example. I get the same results with octave. I think I am doing
something wrong. But, what? Has to be simple, but, I cannot see it.

tomdean

// matrix inversion lema
// given a partitioned matrix, X=[A,B;C,D] with A,D invertible,
// define
// E = D - C*inv(A)*B
// F = A - B*inv(D)*C
// then Y is the inverse of X, where Y is defined by
function [Y] = lema(W)
Y11= inv(A) + inv(A)*B*inv(W)*C*inv(A);
Y12=-inv(A)*B*inv(W);
Y21=-inv(W)*C*inv(A);
Y22=inv(W);
Y = [Y11, Y12; Y21, Y22];
endfunction;

//rand("seed",1);
A=[12,43,27;32,98,33;11,1,2];
B=[1,2,3;4,5,6;7,8,9];
C=B+11;
D=[44,87,1;8,4,9;11,43,88];
X = [A,B;C,D];

inv(A)
inv(D)

E = D - C*inv(A)*B;
F = A - B*inv(D)*C;

inv(E)
inv(F)

// assume E invertible - it is
Y=lema(E);
// clean(X*Y) - I should be zero
//disp(clean(X*Y));
disp(clean(X*Y-eye(Y))== zeros(Y));

// assume F invertible - it is
Y=lema(F);
//disp(clean(X*Y));
disp(clean(X*Y-eye(Y))== zeros(Y));

tomdean

Thomas D. Dean a rit :
Your computation is obviously wrong:
assume A and D invertible, and C=0
then E=D, F=A are both invertible, so if your computation were to make
sense Y should be equal to both lema(E) and lema(F),
so Y22 should be equal to both inv(A) and inv(D), so A=D
this is clearly absurd.

What you want to do is probably some variant of
block-gaussian elimination; try to understand the math behind the
computation and you will be able to correct the formulas.

By the way, yout function lema is not correctly defined:
your mixture of parameters (W) and global variables (A,B,C) will
fail on all practical uses.

> Your computation is obviously wrong:
> assume A and D invertible, and C=0
> then E=D, F=A are both invertible, so if your computation were to make
> sense Y should be equal to both lema(E) and lema(F),
> so Y22 should be equal to both inv(A) and inv(D), so A=D
> this is clearly absurd.
I think you mistake the actions in function calls. Please read a book
on programming and you will be able to correct your mistake.

As a brief attempt to teach you something about programming, assuming
your trivial, meaningless reduction, on the first call to lema, Y22 is
equal to inv(D), NOT inv(A). Then, on the second call to lems, Y22 is
equal to inv(A), not inv(D).

Y22 is not a constant, but, a value calculated from function arguments.

Turns out my problem was a typo.