MatLab - algorithm for finding inverse of matrix

Question

I am trying to write an algorithm in MatLab which takes as its input a lower triangular matrix. The output should be the inverse of this matrix (which also should be in lower triangular form). I have almost managed to solve this, but one part of my algorithm still leaves me scratching my head. So far I have:

function AI = inverse(A)
n = length(A);
I = eye(n);
AI = zeros(n);
for k = 1:n
    AI(k,k) = (I(k,k) - A(k,1:(k-1))*AI(1:(k-1),k))/A(k,k);
    for i = k+1:n
        AI(i,k) = (I(i,k) - (??????????????))/A(i,i);
    end
end

I have marked with question marks the part I am unsure of. I have tried to find a pattern for this part of the code by writing out the procedure on paper, but I just can't seem to find a proper way to solve this part.

If anyone can help me out, I would be very grateful!

Answer 1

Here is my code to get the inverse of a lower triangular matrix by using row transformation:

function AI = inverse(A)
    len = length(A);
    I  = eye(len);
    M  = [A I];
    for row = 1:len
        M(row,:) = M(row,:)/M(row,row);
        for idx = 1:row-1
            M(row,:) = M(row,:) - M(idx,:)*M(row,idx);
        end
    end
    AI = M(:,len+1:end);
end

Answer 2

You can see how it's done on Octave's source. This seems to be implemented in different places depending on the class of the matrix. For Float type Diagonal Matrix it's on liboctave/array/fDiagMatrix.cc, for Complex Diagonal matrix it's on liboctave/array/CDiagMatrix.cc, etc...

One of the advantages of free (as in freedom) software is that you are free to study how things are implemented ;)

Answer 3

Thanks for all the input! I was actually able to find a very nice and easy way to solve this problem today, given that the input is a lower triangular matrix:

function AI = inverse(A)
n = length(A);
I = eye(n);
AI = zeros(n);
for k = 1:n
    for i = 1:n
        AI(k,i) = (I(k,i) - A(k,1:(k-1))*AI(1:(k-1),i))/A(k,k);
    end
end