6 Algebra
First, a method for solving a linear system by a direct algorithm is studied, then an iterative method will be used to determine the eigenvector associated to the largest eigenvalue of a matrix. Finally, the groups generated by a set of permutations will be studied.
Exercise 6.1 - LU decomposition
The LU decomposition consists in decomposing a matrix \(A\) of size \(n \times n\) into the form \(A=LU\), where \(L\) is a lower triangular matrix with 1 on the diagonal and \(U\) an upper triangular matrix. Once the decomposition \(A=LU\) of a matrix is known, it is then very easy to solve the linear problem \(A\boldsymbol{x} = \boldsymbol{b}\) by solving first \(L\boldsymbol{y} = \boldsymbol{b}\) then \(U\boldsymbol{x} = \boldsymbol{y}\). The advantage of the LU decomposition over the Gauss algorithm, for example, is that once the LU decomposition is known, it is possible to solve the linear system quickly with different right-hand sides.
Since \(l_{ik}=0\) if \(k>i\), we have:
\[ a_{ij} = \sum_{k=1}^n l_{ik}u_{kj} = l_{ii}u_{ij} + \sum_{k=1}^{i-1} l_{ik}u_{kj} \,, \]
and therefore as \(l_{ii}=1\):
\[ u_{ij} = a_{ij} - \sum_{k=1}^{i-1} l_{ik}u_{kj} \,. \]
On the contrary, since \(u_{kj}=0\) if \(k>j\), then:
\[ a_{ij} = \sum_{k=1}^n l_{ik}u_{kj} = l_{ij}u_{jj} + \sum_{k=1}^{j-1} l_{ik}u_{kj} \,, \]
and therefore if \(u_{jj}\neq0\):
\[ l_{ij} = \frac{1}{u_{jj}} \left( a_{ij} - \sum_{k=1}^{j-1} l_{ik}u_{kj} \right) \,. \]
Thus, if the first \((i-1)\) rows of \(U\) and the first \((i-1)\) columns of \(L\) are known, it is possible to determine the \(i\)-th row of \(U\) by:
\[ u_{ij} = a_{ij} - \sum_{k=1}^{i-1} l_{ik}u_{kj} \,, \quad j \geq i \,, \]
then, the \(i\)-th column of \(L\) by:
\[ l_{ji} = \frac{1}{u_{ii}} \left( a_{ji} - \sum_{k=1}^{i-1} l_{jk}u_{ki} \right) \,, \quad j > i\ \,. \]
This algorithm for the LU decomposition of a matrix \(A\) requires that \(u_{ii}\) is never zero. This is indeed the case when the matrix \(A\) and all its principal submatrices are invertible.
a. Write a function LU(A) that returns the result of the LU decomposition of a matrix.
b. Given the LU decomposition of a matrix \(A\), write a function solve(L,U,b) that solves the linear system \(Ax = \boldsymbol{b}\).
c. Write a function LU_inplace(A) that does not create new arrays for \(L\) and \(U\) but modifies \(A\) so that its lower triangular part (without the diagonal) matches \(L\) and its upper triangular part (with the diagonal) matches \(U\). Also make a version solve_inplace that takes as input the output of LU_inplace and returns the solution \(\boldsymbol{x}\), without using any new arrays.
d. Using the LU decomposition of the matrix \(A\), write a function det(A) that returns its determinant.
Exercise 6.2 - Power iteration method
The goal of this exercise is to determine the eigenvector of a matrix associated to the largest eigenvalue (in modulus), assuming that this one is unique. Given a real matrix \(A\) of size \(n \times n\) and a vector \(\boldsymbol{x}_0\in\mathbb{R}^n\), the sequence of vectors \((\boldsymbol{x}_k)_ {k\in\mathbb{N}}\) is defined by:
\[ \boldsymbol{x}_{k+1} = \frac{A\boldsymbol{x}_k}{\Vert A\boldsymbol{x}_k\Vert} \,. \]
Assuming, for example, that the matrix \(A\) is diagonalizable, it is then possible to show that the sequence \((\boldsymbol{x}_k)_ {k\in\mathbb{N}}\) converges up to a sign to the eigenvector associated to the largest eigenvalue of \(A\) in absolute value. The convergence takes place almost surely for all choices of \(\boldsymbol{x}_0\).
a. Define a function power(A, x0, k) that returns \(\boldsymbol{x}_k\). With this function, determine the largest eigenvector of the matrix:
\[A = \begin{pmatrix}0.5 & 0.5\\ 0.2 & 0.8 \end{pmatrix} \,.\]
b. Determine the eigenvalue associated with the previous eigenvector.
c. Write a function to automatically compute the largest eigenvalue (in modulus) and the associated eigenvector of a square matrix with a given precision. We will choose the initial vector \(\boldsymbol{x}_0\) randomly.
d. Assuming the matrix \(A\) is diagonalizable with a single eigenvalue of largest modulus, show that the sequence \(\boldsymbol{x}_k\) converges up to one sign to the eigenvector associated with this largest eigenvalue for almost all choices of \(\boldsymbol{x}_0\).
e. Look at the NumPy documentation to find the functions to compute the eigenvectors and eigenvalues of a matrix.
f. Compare the speed of the previous code and the NumPy functions for matrices of size \(n\times n\) with \(n=10,100,1 000\).
Exercise 6.3 - Exponential of matrices
The goal of this exercise is to develop an algorithm to calculate the exponential of a real square matrix. If \(A\) is a real square matrix, its exponential is defined by the series:
\[ \exp(A) = \sum_{k=0}^\infty \frac{A^k}{k!} \,, \]
by analogy with the exponential on real numbers. Here, \(A^k\) represents the matrix product of \(A\) with itself \(k\) times.
a. Define the NumPy arrays corresponding to the matrices \(A_1\) and \(A_2\) defined by:
\[ \begin{align*} A_{1} &= \begin{pmatrix} 1 & 0.8 & 0.6\\ 0.8 & 0.2 & 0.8\\ 0 & 1.2 & 0.9 \end{pmatrix} , & A_{2} &= \begin{pmatrix} 2 & 3 & 2\\ 1 & 2 & 3\\ 4 & 3 & 5.2 \end{pmatrix} . \end{align*} \]
b. Define a function matrix_power(A,n=20) returning an approximation of \(\exp(A)\) obtained by keeping only the first \(n+1\) terms of the series, i.e., the terms from \(k=0\) to \(k=n\).
c. Test on the matrices \(A_1\) and \(A_2\) and compare with the results of the function expm of the module linalg of SciPy.
d. Without using the norm function of NumPy or SciPy, define a function computing the Frobenius norm \(\Vert A\Vert_{F}\) of a matrix \(A\) of size \(m \times m\) defined by:
\[ \Vert A\Vert_{F}^{2}=\operatorname{tr}(AA^{t})=\sum_{i=1}^{m}\sum_{j=1}^{m}|a_{ij}|^{2} \,. \]
Compute the Frobenius norms of the matrices \(A_1\) and \(A_2\).
e. For matrices \(A_1\) and \(A_2\), plot the error in the Frobenius norm between the result of matrix_power and the result of expm as a function of the number of terms \(n\) kept. Put a logarithmic scale on the y-axis.
From a theoretical point of view, it is possible to show that the error is bounded by:
\[ \left\Vert \exp(A)-\sum_{k=0}^{n}\frac{A^{k}}{k!}\right\Vert_{F} \leq\frac{e^{\Vert A\Vert_{F}}}{(n+1)!}\Vert A\Vert_{F}^{n+1} \,. \]
f. Plot this bound as a function of \(n\) for different values of the \(\Vert A\Vert_{F}\) ranging from \(2\) to \(20\), with also a logarithmic scale on the y-axis. Roughly deduce the number of terms to keep so that the bound is lower than the machine precision of \(10^{-15}\) if \(\Vert A\Vert_{F}=20\). Compare the theoretical bound with what was observed in the previous question.
A basic idea to improve the convergence of the series when the norm of the matrix is large is to perform a rescaling using the relation:
\[ \exp{A} = \bigl(\exp(A/p)\bigr)^p \,, \]
for \(p \geq 1\) a well-chosen integer such that \(\Vert A/p \Vert_F\) is small, for example, less than one.
g. Using the previous property, write a function matrix_power_opt(A,n=20) based on this property.
h. Redo the same graph as at point e but with this new function and comment.
Exercise 6.4 - Groups of permutations
The goal is to study the groups of permutations by developing algorithms to characterize them. A group of permutations \(G \subset S_n\) is defined as being generated by a number of permutations: \(G = \langle g_1, g_2,\dots,g_k \rangle\), with \(g_i\in S_n\) a permutation of the set \(\{1,2,\dots,n\}\). A permutation:
\[g = \begin{pmatrix}1 & 2 & 3 & 4\\ 4 & 1 & 2 & 3 \end{pmatrix} ,\]
can be represented in Python by the tuple g = (0, 4, 1, 2, 3). The zero is added in order to conform with Python’s zero-based indexing and thus simplify the implementations a bit. This simply means that vertex 0 goes on vertex 0. Note that this exercise lends itself particularly well to an object-oriented implementation, and in this case the questions can be adapted accordingly.
a. Write a function product(g1,g2) that computes the product of two permutations.
b. Write a function inverse(g) that computes the inverse of a permutation.
c. Write a function that returns the decomposition of a permutation into cycles represented by a list of tuples.
d. Write a function that returns the permutation corresponding to a list of cycles, i.e., that does the inverse of the previous function.
e. ! In Python, a group \(G = \langle g_1,g_2,\dots,g_k \rangle\) generated by a family of permutations, can be represented by a list of permutations G = [g1,g2,...,gk]. Write a function orbit(G,x) that returns the orbit of a point \(x\in\{1,2,\dots,n\}\):
\[ O_x = Gx = \{gx, g \in G\} \,. \]
f. !! Define a function stabilizer(G,x) that returns the stabilizer of a point \(x\in\{1,2,\dots,n\}\):
\[ G_x = \{g \in G : g x = x \} \,, \]
in the form of a set of generators.