Exo 1

Soient λ une valeur propre A et x = (x₁, ..., x_n) un vecteur propre associé.

Pour $i$ compris entre 1 et $n$, on a

$(\lambda -a_{ii})x_{i} = \sum_{j\neq i} a_{ij} x_{j}$

Choisissons un indice i pour lequel le module de $x_i$ est maximal.

Puisque x est un vecteur propre, $|x_i|$ est non nul et il est possible de former le quotient $\frac {|x_j|}{|x_i|} \leq 1$.

On a :

$|a_{ii}-\lambda |$ $=\left|\sum _{j\neq i}a_{ij}{\frac {x_{j}}{x_{i}}}\right|$ $\leq \sum _{j\neq i}\left|a_{ij}{\frac {x_{j}}{x_{i}}}\right|$ $\leq \sum _{j\neq i}|a_{ij}|$

---

Exo 2

On veut $\min |\lambda_i - \lambda | \leq \epsilon$

J'ai fait au tableau par l'absurde :

$ (A - \lambda) u = (A - \lambda) \sum u_i v_i = A (\sum u_i v_i) - \lambda (\sum u_i v_i)$

$= \sum u_i A v_i - \sum u_i (\lambda v_i) = \sum u_i (\lambda_i v_i) - \sum u_i (\lambda v_i)$

$= \sum u_i (\lambda_i - \lambda ) v_i$

D'après Pythagore

$\| (A - \lambda) u \|^2 = \sum u_i^2 |\lambda_i - \lambda |^2$

et si $\min |\lambda_i - \lambda | > \epsilon$ alors

$ \| (A - \lambda) u \|^2 > \epsilon^2 \| u \|^2$

import numpy as np

import matplotlib.pyplot as plt

---

Méthode de la puissance:

est un algorithme permettant d'estimer

• la valeur propre de plus grand module
• un vecteur propre associé

d'une matrice $A$. L'idée est que, si on définit une suite $(x_n)$ de vecteurs de $\mathbb R^n$ par la relation $$x_{n+1}=Ax_n,$$ alors les vecteurs $x_n$ vont "se tourner" dans la direction des vecteurs propres associés à la plus grande des valeurs propres.

---

Théorème : Soit $A\in\mathcal M_n(\mathbb R)$ une matrice diagonalisable

• $\lambda_1,\cdots,\lambda_n$ ses valeurs propres,
• on suppose que $|\lambda_1| > |\lambda_2| \geq \dots\geq|\lambda_n|)$.
• $v_1$ un vecteur propre associé à $\lambda_1$.

Alors, si $x_0$ est un vecteur non orthogonal à $v_1$, et si $(x_n)$ est définie par la relation de récurrence $x_{n+1}=Ax_n$, alors

• $\frac{x_n}{\|x_n\|}$ converge vers un multiple de $v_1$
• $\frac{\langle x_{n+1},x_n\rangle}{\|x_n\|^2_2}$ converge vers $\lambda_1$.

---

1. Pour éviter que l'algorithme ne provoque un dépassement de capacité, on calcule souvent directement une version normalisée de $x_n$.
2. Remarquons aussi que l'hypothèse sur $x_0$ n'est pas discriminante, si on choisit le vecteur au hasard, il a une probabilité nulle d'être orthogonal à $v_1$.
3. Lorsque les multiplicités algébriques et géométriques associées à la valeur propre $\lambda _1$ sont égales, le taux de convergence de l'algorithme se comporte en $O\left(|\lambda _2/\lambda _1|^n\right)$ où $\lambda _1,\lambda _2$ sont la plus grande et la seconde plus grande valeurs propres (en valeur absolue).
4. L'intérêt de la méthode des puissances itérées est de ne pas nécessiter d'hypothèses très particulières sur la matrice, comme être symétrique. En revanche, elle demande une configuration particulière des valeurs propres qui n'est pas toujours vérifiée.

A = np.array([99,1,0,1,100,1,0,1,98]).reshape(3,3)  # - 97 * np.identity(3)
A

array([[ 99,   1,   0],
       [  1, 100,   1],
       [  0,   1,  98]])

np.linalg.eig(A)[0]

array([100.87938524,  98.65270364,  97.46791111])

---

Calculate the Gershgorin circles

why not it's easy ?

centers  = np.diag(A).astype(complex)
radii = np.sum(np.abs(A - np.diag(centers)), axis = 1)

fig, ax = plt.subplots()
ax.set_title('Gershgorin circles')

for cc, rr in zip(centers, radii):
    # this is dumb but must plot the centers first
    ax.plot([cc.real],[cc.imag])
    ax.add_patch(plt.Circle((cc.real,cc.imag),rr, 
                             color='r', alpha=.5))

---

Coding the algorithm

• isn't memory efficient for just finding the eigenvector
• since we are going to study convergence it's the best
• and it is easy to read.

MAXITER = 1000
EPSILON = 10**-11

# I know how many there will be so don't need a dynamic list
x =  np.zeros((MAXITER + 1, 3))
# initial vector should check not eigenvector
x[0] =  A[0]

for k in range(MAXITER):    
    #iterate
    x[k+1] = A @ x[k]
    # normalise see above
    x[k+1] /= np.linalg.norm(x[k+1])   
    #stopping condition
    if np.linalg.norm(x[k+1] - x[k]) < EPSILON: break
    
vecs = x[:k]

#eigenvalue

vecs[-1]@A@vecs[-1]

100.87938524157182

def pow_meth(A,
            MAXITER = 1000,
            EPSILON = 10**-11):
    
    # I know how many there will be so don't need a dynamic list
    x =  np.zeros((MAXITER + 1,3))
    # initial vector should check not eigenvector
    x[0,0] =  1

    for k in range(MAXITER):    
        #iterate
        x[k+1] = A @ x[k]
        # normalise see above
        x[k+1] /= np.linalg.norm(x[k+1])   
        #stopping condition
        if np.linalg.norm(x[k+1] - x[k]) < EPSILON: break
    else:
        raise ValueError

    return x[:k]

Memory

low memory use version

MAXITER = 1000
EPSILON = 10**-11

# initial vector should check not eigenvector
# you should program it
x = A[0].astype(float)

for _ in range(MAXITER):    
    x, y = A@x, x
    # normalise 
    x /= np.linalg.norm(x) 
    #stopping condition
    if np.linalg.norm(x - y) < EPSILON: break

# eigenvalue, eigenvector
y@A@y, y

(100.8793852415718, array([0.44909879, 0.84402963, 0.29312841]))

and compare it with numpy.

np.linalg.eigvals(A)[0]

100.87938524157187

---

Convergence

import matplotlib.pyplot as plt

vecs = pow_meth(A )
Y = [np.linalg.norm(v) for v in vecs[:-1] - vecs[-1]]

YY = np.log(Y) 
plt.plot(YY)

# calculate the rate of convergence
a,b = 450,200
slope = (YY[b] - YY[a])/(b-a)
plt.title(f'log plot for $A$ slope = {slope}');

the slope should be $\log( \lambda_2/\lambda_1)$.

vv = np.linalg.eig(A)[0]
np.log(vv[1]/vv[0])

-0.022319959154051165

---

vecs = pow_meth(A - 97*np.identity(3))

Y = np.array([np.linalg.norm(v) for v in vecs[:-1] - vecs[-1]]) 

YY = np.log(Y) 
plt.plot(YY)
a,b = 5,20
slope = (YY[b] - YY[a])/(b-a)
plt.title(f'log plot for  $A - 97I_3$ slope = {slope}');

vv = np.linalg.eig(A - 97*np.identity(3))[0]

np.log(vv[1]/vv[0])

-0.8532641787457773

'/home/macbuse'