0
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Nox Feb 20 at 16:39
4
$\begingroup$

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)).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.