# Finding Matrix inverse with LU and repeted left division calls

Hello I am in a basic numerical methods class and our teacher has given us an algorithm which can compute the inverse of a matrix other than using MATLAB's built in library function.

A = [1 2 4; 1 3 9; 1 4 16];
b = [7; 13; 21];
[L, U, P] = lu(A);
n = 3;
X = zeros(n,n);
for i=1:n
e = zeros(n,1);
e(i) = 1;
X(:,i) = U\(L\(P*e))
end


I am a bit new to MATLAB, so I Am a bit unsure as to what this code is actually doing to find the inverse. I naively set n = 3 because it was the only value i could find to get the algorithm to work.

Now I know this algorithm finds the LU factorization of A then uses the \ operator repeatedly.

what I do not understand is what x = zeros(n,n) what is the point of this line as well as the inner for loop

e = zeros(n,1); //creates a 3 x 1 matrix of zeroes?
e(i) = 1;       //No idea what this line is doing
X(:,i) = U\(L\(P*e))


As for the last line of the for loop, This is the repeated "\" operation but I have no idea what the X(:,i) stands for as well as well using "\" by P*e

Could someone help me understand this algorithm better? Because I am only understanding fragments of it at the moment.

• This looks like a homework question to me. I would advise you to look at any introduction to Matlab your teacher may have given you, or really any one available and then take a basic book on linear algebra and read it. IMHO deciphering the workings of the code this way is more productive (in terms of learning) than asking - at least for such simple codes. – Nox Feb 20 '19 at 16:39

The inverse of $$A$$ is the solution to the linear system of equations

$$AX=I$$.

If we write $$X$$ in terms of its columns as

$$X=\left[ x^{(1)}, x^{(2)}, \ldots, x^{(n)} \right]$$,

then we can break up $$AX=I$$ into $$n$$ separate systems of equations

$$Ax^{(1)}=\left[ \begin{array}{c} 1 \\ 0 \\ \vdots \\ 0 \end{array} \right]$$,

$$Ax^{(2)}=\left[ \begin{array}{c} 0 \\ 1 \\ \vdots \\ 0 \end{array} \right]$$, $$\vdots$$

$$Ax^{(n)}=\left[ \begin{array}{c} 0 \\ 0 \\ \vdots \\ 0 \\ 1 \end{array} \right]$$.

A common notation for a vector of all zeros except for a 1 in the $$i$$th position is $$e^{(i)}$$. With that notation, these systems are $$Ax^{(i)}=e^{(i)}$$, for $$i=1, 2, \ldots, n$$. In MATLAB you can create a vector $$e^{(i)}$$ with e=zeros(n,1); e(i)=1.

The code you were given solves these $$n$$ linear systems of equations using the LU factorization of $$A$$.

The LU factorization (with pivoting) is

$$PA=LU$$.

So, to solve $$Ax=b$$, multiply by $$P$$ on both sides of the equation to get

$$PAx=Pb$$

or

$$LUx=Pb$$.

Next, if we solve

$$Lw=Pb$$

and

$$Ux=w$$

then

$$PAx=LUx=Lw=Pb$$ and $$Ax=b$$. The MATLAB \ function is used in code to solve the triangular systems. In particular, we get $$w$$ by w=L\(P*b) and then get $$x$$ by x=U\(L\(P*b)).