How to invert a square matrix using Gaussian elimination

Hi there,

I have attached below the code I have written which is supposed to invert a square n x n matrix by Gaussian elimination. This m-file is as part of a larger project where it is implemented into a Newton-Raphson iterative method. The formula for this is Xn = Xn - (Df^-1)*(f). When i use this Gaussian code to invert the matrix, the Newton-Raphson method does not converge so any help on solving the problems in my Gaussia code attached below would be brilliant, thank you!

A = [1 1;2 5];
A_inv = Gauss_inverse(A)

A_inv = 2x2

0 1 1 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

A*A_inv

ans = 2x2

1 1 5 2

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

function A_inv = Gauss_inverse(A)
% Find dimensions of input matrix
[p,q] = size(A);
% If input matrix is not square, stop function
if p ~= q
    error('Input must be a square n x n matrix')
end
% Target identity matrix to be transformed into the output
% inverse matrix
A_inv = eye(p);
%The following code actually performs the matrix inversion by working
% on each element of the input
for j = 1 : p
    % Pivot row swapping for improved stability
    [~,max_row] = max(abs(A(j:p, j))); % Find the row index of maximum element in the column
    max_row = max_row + j - 1; % Adjust index to global reference
    if A(max_row,j) ~= 0 % Skip if already zero
        % Swap rows in both A and A_inv
        temp = A(j,:);
        A(j,:) = A(max_row,:);
        A(max_row,:) = temp;
        
        temp = A_inv(j,:);
        A_inv(j,:) = A_inv(max_row,:);
        A_inv(max_row,:) = temp;
    end
end
if A(j,j) ~= 0
    % Scale the pivot row to have 1 on diagonal
    t = 1/A(j,j);
    A(j,:) = t * A(j,:);
    A_inv(j,:) = t * A_inv(j,:);
end

% Eliminate non-zero elements below the pivot
for i = j+1 : p
    if A(i,j) ~= 0
        t = -A(i,j);
        A(i,:) = A(i,:) + t * A(j,:);
        A_inv(i,:) = A_inv(i,:) + t * A_inv(j,:);
    end
end
end

Answer 1

[Torsten]

The newton update is

x_(n+1) = x_n - Df^(-1)*f, 

not

x_(n+1) = x_n + Df^(-1)*f.

But you should not explicitly compute the inverse, but solve the linear system

Df*dx = -f 

and set

x_(n+1) = x_n + dx

And as you can see from the simple example from above: your code to find the inverse of a matrix does not work as expected.

Answer 2

[John D'Errico]

I assume you are forced to use your own code for gaussian elimination here. A bad idea if you are not. NEVER write such code yourself, when high quality code written by professionals in the art is available.

If you are forced to write your own Newton method, An excellent suggestion is to do as Torsten suggests, to instead solve a system of equations. Best of course is to use backslash.

Ok, if you are forced to use an inverse there, and one you have written yourself, then first, verify that your code works! Never just throw something together and not test the pieces!

Once you have verified that piece of your code, then you need to recognize that Newton's method does not always converge. There is no assurance of convergence for all such problems. Sorry, but there is not. So, IF you are seeing non-convergence, you would want to learn to use the debugger. Are the Newton steps going in the correct direction? Does your code do something that makes sense?