In this worksheet we will focus on solving the linear system given by \begin{equation} Ax = b, \end{equation} where $A$ is a non-singular $n\times n$ matrix.
In tis worksheet you will find the $LU$ factorization of the Matrix $A$, such that \begin{equation} L\cdot U = A, \end{equation} where $L$ is lower triangular and $U$ is upper triangular.
If the LU factors are known, then solving the linear system can be performed in two triangular solves. \begin{equation} Lp = b, \end{equation} and \begin{equation} Ux = p. \end{equation}
To warm up we will compute the LU factorization of a $2\times 2$ matrix.
We define the matrix
In [13]:
A = [1. 3.; 2. 5.]
Out[13]:
We check that the matrix is invertible by computing its determinant, which should be non-zero.
In [14]:
det(A)
Out[14]:
We make space for the LU factors
In [15]:
L = zeros(Float64, 2, 2);
U = zeros(Float64, 2, 2);
We have that
\begin{equation} A = \left [ \begin{array}{cc} a_{1,1} & a_{1,2} \\ a_{2,1} & a_{2,2} \end{array} \right ], \qquad L = \left [ \begin{array}{cc} l_{1,1} & 0 \\ l_{2,1} & l_{2,2} \end{array} \right ], \qquad U = \left [ \begin{array}{cc} u_{1,1} & u_{1,2} \\ 0 & u_{2,2} \end{array} \right ]. \end{equation}Then we have that \begin{equation} A = \left [ \begin{array}{cc} a_{1,1} & a_{1,2} \\ a_{2,1} & a_{2,2} \end{array} \right ] = \left [ \begin{array}{cc} l_{1,1}u_{1,1} & l_{1,1}u_{1,2} \\ l_{2,1}u_{1,1} & l_{2,1}u_{1,2} + u_{2,2}l_{2,2} \end{array} \right ]. \end{equation}
We need to find a combination of $l_{1,1}u_{1,1} = a_{1,1}$.
In [16]:
L[:,1] = A[:,1]/A[1,1]
U[1,:] = A[1,:]
println("U = \n",U)
println("L = \n",L)
We have already factorized the first column and row of $A$. So we can continue with the next row and column.
In this case we have that
\begin{equation}
l_{2,1}u_{1,2} + u_{2,2}l_{2,2} = a_{2,2},
\end{equation}
or, given that $l_{2,1}$ and $u_{1,2}$ are known we have that
\begin{equation}
u_{2,2}l_{2,2} = a_{2,2} - l_{2,1}u_{1,2}.
\end{equation}
By convention we have that $l_{2,2} = 1$, then we have that $u_{2,2} = a_{2,2} - l_{2,1}u_{1,2}$.
Which can be written in julia as
In [17]:
L[2,2] = 1;
U[2,2] = A[2,2] - L[2,1]*U[1,2];
And we can test that the factorization is accurate by
In [18]:
vecnorm(L*U - A)
Out[18]:
Now we will generalize the LU factorization for any matrix, for simplicity this will performed recursively.
Let $A$ be a $n\times n$ matrix given by \begin{equation} A = \left [ \begin{array}{cc} a_{1,1} & A_{1,2} \\ A_{2,1} & A_{2,2} \end{array} \right ] \end{equation} where
We have \begin{equation} L = \left [ \begin{array}{cc} l_{1,1} & 0 \\ L_{2,1} & L_{2,2} \end{array} \right ], \qquad U = \left [ \begin{array}{cc} u_{1,1} & U_{1,2} \\ 0 & U_{2,2} \end{array} \right ]. \end{equation} Where
In [27]:
m = 4
A = rand(m,m) + m*eye(m)
L = zeros(A)
U = zeros(A)
Out[27]:
Given that we will modify A, we will save a copy to make the comparison afterwards.
In [28]:
Acopy = copy(A)
Out[28]:
As before we have that $l_{1,1} = 1$ so $u_{1,1} = a_{1,1}$.
But we have that
\begin{equation}
A = \left [ \begin{array}{cc}
a_{1,1} & A_{1,2} \\
A_{2,1} & A_{2,2}
\end{array}
\right ] = \left [ \begin{array}{cc}
l_{1,1}u_{1,1} & l_{1,1}U_{1,2} \\
L_{2,1}u_{1,1} & L_{2,1}U_{1,2} + L_{2,2}U_{2,2}
\end{array}
\right ].
\end{equation}
which means that
$L_{2,1} = A_{2,1}/A_{1,1}$ and $U_{1,2} = A_{1,2}$.
We can then save the first row of $U$ and the first columns of $L$ using
In [29]:
U[1,:] = A[1,:];
L[:,1] = A[:,1]/A[1,1];
println("U = \n",U)
println("L = \n",L)
Then we have that
\begin{equation}
L_{2,2}U_{2,2} = A_{2,2} - L_{2,1}U_{1,2},
\end{equation}
where $L_{2,1}$ and $U_{1,2}$ are known.
This means that $L_{2,2}$ and $U_{2,2}$ are the LU factors of the $(n-1\times n-1)$ matrix $A_{2,2} - L_{2,1}U_{1,2}$.
We have "eliminated" the first rown and column of A, and we have reduced to problem to finding the LU factors of a smaller matrix.
You can write this in Julia as
In [30]:
A[1:end,1:end] = A[1:end,1:end] - L[1:end,1]*U[1,1:end]
Out[30]:
You can see that the first column and row of $A$ have been eliminated.
Now we can apply the same procedure to the modified matrix.
We eliminate the second row and second columns.
In [32]:
U[2,2:end] = A[2,2:end]
L[2:end,2] = A[2:end,2]/A[2,2]
println("U = \n",U)
println("L = \n",L)
And we compute the smaller matrix
In [33]:
A[2:end,2:end] = A[2:end,2:end] - L[2:end,2]*U[2,2:end]
println("A = \n",A)
And we continue
In [35]:
U[3,3:end] = A[3,3:end]
L[3:end,3] = A[3:end,3]/A[3,3]
println("U = \n",U)
println("L = \n",L)
We eliminate another colum and row
In [37]:
A[3:end,3:end] = A[3:end,3:end] - L[3:end,3]*U[3,3:end]
println("A = ",A)
Out[37]:
In [ ]:
We need to compute the LU factors of just a number,
In [38]:
L[end,end] = 1
U[end,end] = A[end,end]
Finally, we check that our LU factorization is accurate
In [39]:
println(maximum(abs(L*U-Acopy)))
For homework 7, you will need to manipulate a very simple sparse matrix.
A sparse matrix, is a matrix which is sparsely populated, in other words, it has a few non-zero entries per row. When you have a matrix like that, storing all the zeros can become extremely expensive. In order to overcome this problem, we can represent sparse matrices without storing the zeros explicitely.
On typical example is the identity matrix, that has a special contructor
In [41]:
Isp = speye(m)
Out[41]:
You can observe that the sparse matrix contains only the non zero elements.
This fact is extremely valuable not only for storage (a sparse matrix uses much less memory than its dense counterpart), but for any operation. If a $n\times n$ sparse matrix has $\mathcal{O}(1)$ non-zeros per row, or it has only $\mathcal{O}(n)$ non-zeros entries, then a matrix vector multiplication can be performed in $\mathcal{O}(n)$ instead of $\mathcal{O}(n^2)$.
You can run the snippet below, and see how the run-times scale with respect to $n$, as $n$ increases
In [56]:
n = 10000
r = rand(n)
Msp = spdiagm(r, 0);
M = diagm(r,0);
b = rand(n,1)
@time xsp = Msp*b;
@time x =M*b;
norm(x-xsp)
Out[56]:
In julia Sparse matrices behave exaxtly as normal matrices. You can index them, slice them and reorder them.
For example you can swap rows by writting
In [64]:
i = 2;
j = 3;
Isp = speye(4)
# we use full to convert back the sparse amtrix to a dense matrix
# so it "prints" better in your console
println("original matrix \n ",full(Isp))
buffer = Isp[i,:]
Isp[i,:] = Isp[j,:]
Isp[j,:] = buffer
println("permuted matrix \n",full(Isp))
In [ ]: