Introduction
In this article, we show you how to invert a matrix without going through the burdens of inversions, i.e. without doing 
operations, which is required by Gauss Jordan. Herein, a Newton method is adopted for matrix inversion.
The approach is via Newton’s method and instead of computing explicitly 
, we solve 
whose solution is obviously .
Approach
Given matrix 
, define
(1) 
, since (2) 
are (3) ![\begin{equation*} \pmb{X}_k = \pmb{X}_{k-1} - [\pmb{F}'(\pmb{X}_{k-1} )]^{-1}\pmb{F}(\pmb{X}_{k-1}) \end{equation*}](https://bazziahmad.com/wp-content/ql-cache/quicklatex.com-488c11b45094f30e017a2fd686c3333e_l3.png "Rendered by QuickLaTeX.com")
hence (4) 
is what we trying to estimate by iterations, and hence replace it by the most recent estimate of , namely 
(5) 
(6) 
(7) 
using equation (6) then plugging equation (4) at 
(8) ![\begin{equation*} \begin{split} \pmb{R}_{k+1} = \pmb{I} - \pmb{A}\pmb{X}_{k+1} &= \pmb{I} - \pmb{A}\big[ \pmb{X}_{k} +\pmb{X}_{k} \big( \pmb{I} - \pmb{A} \pmb{X}_{k} \big) \big] \\ &= \pmb{I} - \pmb{A}\pmb{X}_k - \pmb{A} \pmb{X}_{k} \big( \pmb{I} - \pmb{A} \pmb{X}_{k} \big) \\ &= \pmb{I} - \pmb{A}\pmb{X}_k - \pmb{A}\pmb{X}_k + \pmb{A}\pmb{X}_k \pmb{A}\pmb{X}_k \\ &= \pmb{I} - 2\pmb{A}\pmb{X}_k + \pmb{A}\pmb{X}_k \pmb{A}\pmb{X}_k \end{split} \end{equation*}](https://bazziahmad.com/wp-content/ql-cache/quicklatex.com-eba80ec788d8c10131908ffe6f88ac1a_l3.png "Rendered by QuickLaTeX.com")
(9) 
(10) ![\begin{equation*} \begin{split} \pmb{E}_{k+1} = \pmb{A}^{-1} - \pmb{X}_{k+1} &= \pmb{A}^{-1} - \big[ \pmb{X}_{k} +\pmb{X}_{k} \big( \pmb{I} - \pmb{A} \pmb{X}_{k} \big) \big] \\ &= \pmb{A}^{-1} - \pmb{X}_{k} -\pmb{X}_{k} \big( \pmb{I} - \pmb{A} \pmb{X}_{k} \big) \\ &= \pmb{A}^{-1} - \pmb{X}_{k} - \pmb{X}_{k} + \pmb{X}_{k} \pmb{A} \pmb{X}_{k} \\ &= \pmb{A}^{-1} - 2 \pmb{X}_{k} + \pmb{X}_{k} \pmb{A} \pmb{X}_{k} \end{split} \end{equation*}](https://bazziahmad.com/wp-content/ql-cache/quicklatex.com-9a8393cf2f96e32719caaf7089915889_l3.png "Rendered by QuickLaTeX.com")
MATLAB implementation
So, Let’s start with building a function called NewtonInverse, which computes the inverse based on equation (5). This is simply attained as such
function[X] = NewtonInverse(A,N,Niter)
%% NewtonInverse: computes inverse of A
%input
%A - matrix
%N - matrix dimension
%Niter - maximum number of iterations
I = eye(N);
%the true inverse to compute the error.
A_trueinverse = inv(A);
%initial guess
X(:,:,1) = A'/(norm(A,1) * norm(A,inf));
%iterations
for n = 2:Niter
X(:,:,n) = X(:,:,n-1) + X(:,:,n-1)*(I - A*X(:,:,n-1));
error(n-1) = norm(X(:,:,n) - A_trueinverse,2);
endFor sake of completeness, we will see how the iterations converge and compare it with an inverse based off MATLAB’s LU decomposition (which is not iterative herein). So we shall add the following block
%Using LU decompostion
[L1,U] = lu(A);
A_inverse_LU = inv(U)*inv(L1);
error_LU = norm(A_inverse_LU - A_trueinverse,2);
figure
plot(log10(error),'k','Linewidth',2)
hold on
plot(log10(error_LU*ones(Niter,1)),'--*r','Linewidth',2)
xlabel('n iteration')
ylabel('error (log scale)')
legend('Using Newton','Matlabs LU decomposition')
grid on
grid minorNow our function looks like this
function[X] = NewtonInverse(A,N,Niter)
%% NewtonInverse: computes inverse of A
%input
%A - matrix
%N - matrix dimension
%Niter - maximum number of iterations
I = eye(N);
%the true inverse to compute the error.
A_trueinverse = inv(A);
%initial guess
X(:,:,1) = A'/(norm(A,1) * norm(A,inf));
%iterations
for n = 2:Niter
X(:,:,n) = X(:,:,n-1) + X(:,:,n-1)*(I - A*X(:,:,n-1));
error(n-1) = norm(X(:,:,n) - A_trueinverse,2);
end
%Using LU decompostion
[L1,U] = lu(A);
A_inverse_LU = inv(U)*inv(L1);
error_LU = norm(A_inverse_LU - A_trueinverse,2);
figure
plot(log10(error),'k','Linewidth',2)
hold on
plot(log10(error_LU*ones(Niter,1)),'--*r','Linewidth',2)
xlabel('n iteration')
ylabel('error (log scale)')
legend('Using Newton','Matlabs LU decomposition')
grid on
grid minor
Our main function will then test the above by simply generating a random matrix and inputting it to the above function as
A = 50*randn(3,3);
N = 3;
Niter = 30;
[X] = NewtonInverse(A,N,Niter);
Summary
We have presented a method of matrix inversion through Newton by iterating over an appropriate cost function. MATLAB code is detailed as well.
Finally, Subscribe to my channel to support us ! We are very close to 100K subscribers <3 w/ love. The algorithmic chef – Ahmad Bazzi. Click here to subscribe. Check my other articles here.
Recent Comments