(EIGVALEIGVEC)=

# 2.2 Eigenvalores y eigenvectores

```{admonition} Notas para contenedor de docker:

Comando de docker para ejecución de la nota de forma local:

nota: cambiar `<ruta a mi directorio>` por la ruta de directorio que se desea mapear a `/datos` dentro del contenedor de docker.

`docker run --rm -v <ruta a mi directorio>:/datos --name jupyterlab_optimizacion -p 8888:8888 -d palmoreck/jupyterlab_optimizacion:2.1.4`

password para jupyterlab: `qwerty`

Detener el contenedor de docker:

`docker stop jupyterlab_optimizacion`

Documentación de la imagen de docker `palmoreck/jupyterlab_optimizacion:2.1.4` en [liga](https://github.com/palmoreck/dockerfiles/tree/master/jupyterlab/optimizacion).

```

---

Nota generada a partir de [liga](https://www.dropbox.com/s/s4ch0ww1687pl76/3.2.2.Factorizaciones_matriciales_SVD_Cholesky_QR.pdf?dl=0).

```{admonition} Al final de esta nota el y la lectora:
:class: tip

* Aprenderá las definiciones más relevantes en el tema de eigenvalores y eigenvectores para su uso en el desarrollo de algoritmos en el análisis numérico en la resolución de problemas del álgebra lineal numérica. En específico las definiciones de: diagonalizable o *non defective* y similitud son muy importantes.

* Comprenderá el significado geométrico de calcular los eigenvalores y eigenvectores de una matriz simétrica para una forma cuadrática que define a una elipse.

* Aprenderá cuáles problemas en el cálculo de eigenvalores y eigenvectores de una matriz son bien y mal condicionados.

```

En esta nota **asumimos** que $A \in \mathbb{R}^{n \times n}$.

## Eigenvalor (valor propio o característico)

```{admonition} Definición

El número $\lambda$ (real o complejo) se denomina *eigenvalor* de A si existe $v \in \mathbb{C}^n - \{0\}$ tal que $Av = \lambda v$. El vector $v$ se nombra eigenvector (vector propio o característico) de $A$ correspondiente al eigenvalor $\lambda$.
```

```{admonition} Observación
:class: tip 

Observa que si $Av=\lambda v$ y $v \in \mathbb{C}^n-\{0\}$ entonces la matriz $A-\lambda I_n$ es singular por lo que su determinante es cero.

```

```{admonition} Comentarios

* Una matriz con componentes reales puede tener eigenvalores y eigenvectores con valores en $\mathbb{C}$ o $\mathbb{C}^n$ respectivamente.
* El conjunto de eigenvalores se le nombra **espectro de una matriz** y se denota como:

$$\lambda(A) = \{ \lambda | \det(A-\lambda I_n) = 0\}.$$

* El polinomio 

$$p(z) = \det(A-zI_n) = (-1)^nz^n + a_{n-1}z^{n-1}+ \dots + a_1z + a_0$$

se le nombra **polinomio característico asociado a $A$** y sus raíces o ceros son los eigenvalores de $A$.

* La multiplicación de $A$ por un eigenvector es un reescalamiento y posible cambio de dirección del eigenvector.
* Si consideramos que nuestros espacios vectoriales se definen sobre $\mathbb{C}$ entonces siempre podemos asegurar que $A$ tiene un eigenvalor con eigenvector asociado. En este caso $A$ tiene $n$ eigenvalores y pueden o no repetirse.

* Se puede probar que el determinante de $A$: $\det(A) = \displaystyle \prod_{i=1}^n \lambda_i$ y la traza de $A$: $tr(A) = \displaystyle \sum_{i=1}^n \lambda_i$.
```

### Ejemplo


```python
import numpy as np
```


```python
np.set_printoptions(precision=3, suppress=True)
```


```python
A=np.array([[10,-18],[6,-11]])
```


```python
print(A)
```




    array([[ 10, -18],
           [  6, -11]])



**En *NumPy* con el módulo [numpy.linalg.eig](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eig.html) podemos obtener eigenvalores y eigenvectores** 


```python
evalue, evector = np.linalg.eig(A)
```


```python
print('eigenvalores:')
print(evalue)
print('eigenvectores:')
print(evector)
```

    eigenvalores:
    [ 1. -2.]
    eigenvectores:
    [[0.894 0.832]
     [0.447 0.555]]


```{margin}

$Av_1 = \lambda_1 v_1$.
```


```python
print('matriz * eigenvector:')
print(A@evector[:,0])
print('eigenvalor * eigenvector:')
print(evalue[0]*evector[:,0])
```

    matriz * eigenvector:
    [0.894 0.447]
    eigenvalor * eigenvector:
    [0.894 0.447]


```{margin}

$Av_2 = \lambda_2 v_2$.
```


```python
print('matriz * eigenvector:')
print(A@evector[:,1])
print('eigenvalor * eigenvector:')
print(evalue[1]*evector[:,1])
```

    matriz * eigenvector:
    [-1.664 -1.109]
    eigenvalor * eigenvector:
    [-1.664 -1.109]


### Ejemplo 

Si $v$ es un eigenvector entonces $cv$ es eigenvector donde: $c$ es una constante distinta de cero.



```python
const = -2
const_evector = const*evector[:,0]
print(const_evector)
```

    [-1.789 -0.894]


```{margin}

$cv$ es un eigenvector con eigenvalor asociado $c\lambda$ pues $A(cv) = \lambda(cv)$ se satisface si $Av = \lambda v$ y $c \neq 0$.

```


```python
print('matriz * (constante * eigenvector):')
print(A@const_evector)
print('eigenvalor * (constante * eigenvector):')
print(evalue[0]*const_evector)
```

    matriz * (constante * eigenvector):
    [-1.789 -0.894]
    eigenvalor * (constante * eigenvector):
    [-1.789 -0.894]


### Ejemplo 

Una matriz con entradas reales puede tener eigenvalores y eigenvectores complejos:


```python
A=np.array([[3,-5],[1,-1]])
```


```python
print(A)
```

    [[ 3 -5]
     [ 1 -1]]



```python
evalue, evector = np.linalg.eig(A)
```

```{margin}

Para $A \in \mathbb{R}^{n \times n}$ se tiene: $\lambda \in \mathbb{C}$ si y sólo si $\bar{\lambda} \in \mathbb{C}$ con $\bar{\lambda}$ el conjugado de $\lambda$.

```


```python
print('eigenvalores:')
print(evalue)
print('eigenvectores:')
print(evector)
```

    eigenvalores:
    [1.+1.j 1.-1.j]
    eigenvectores:
    [[0.913+0.j    0.913-0.j   ]
     [0.365-0.183j 0.365+0.183j]]


### Ejemplo

Los eigenvalores de una matriz diagonal son iguales a su diagonal y sus eigenvectores son los vectores canónicos $e_1, e_2, \dots e_n$.


```python
A = np.diag([2, 2, 2, 2])
```


```python
print(A)
```

    [[2 0 0 0]
     [0 2 0 0]
     [0 0 2 0]
     [0 0 0 2]]



```python
evalue, evector = np.linalg.eig(A)
```


```python
print('eigenvalores:')
print(evalue)
print('eigenvectores:')
print(evector)
```

    eigenvalores:
    [2. 2. 2. 2.]
    eigenvectores:
    [[1. 0. 0. 0.]
     [0. 1. 0. 0.]
     [0. 0. 1. 0.]
     [0. 0. 0. 1.]]


```{admonition} Definición

La **multiplicidad algebraica** de un eigenvalor es su multiplicidad considerado como raíz/cero del polinomio característico $p(z)$. Si no se repite entonces tal eigenvalor se le nombra de multiplicidad **simple**.

La **multiplicidad geométrica** de un eigenvalor es el número máximo de eigenvectores linealmente independientes asociados a éste.

```

```{admonition} Observación
:class: tip

En el ejemplo anterior cada eigenvalor tiene una multiplicidad simple y la multiplicidad geométrica de cada eigenvalor es $1$.

```

### Ejemplo

Los eigenvalores de una matriz triangular son iguales a su diagonal.


```python
A=np.array([[10,0, -1],
            [6,10, 10],
            [3, 4, 11.0]])
A = np.triu(A)
```


```python
print(A)
```




    array([[10.,  0., -1.],
           [ 0., 10., 10.],
           [ 0.,  0., 11.]])




```python
evalue, evector = np.linalg.eig(A)
```

```{margin}

Observa que el eigenvalor igual a $10$ está repetido dos veces (multiplicidad algebraica igual a $2$) y se tienen dos eigenvectores linealmente independientes asociados a éste (multiplicidad geométrica igual a $2$). 
```


```python
print('eigenvalores:')
print(evalue)
print('eigenvectores:')
print(evector)
```

    eigenvalores:
    [10. 10. 11.]
    eigenvectores:
    [[ 1.     0.    -0.099]
     [ 0.     1.     0.99 ]
     [ 0.     0.     0.099]]


**Otro ejemplo:**


```python
A=np.array([[10,18, -1],
            [6,10, 10],
            [3, 4, 11.0]])
A = np.triu(A)
```


```python
print(A)
```




    array([[10., 18., -1.],
           [ 0., 10., 10.],
           [ 0.,  0., 11.]])




```python
evalue, evector = np.linalg.eig(A)
```

```{margin}

Observa que en este ejemplo el eigenvalor $10$ está repetido dos veces (multiplicidad algebraica es igual a $2$) y sus eigenvectores asociados son linealmente dependientes (multiplicidad geométrica es igual a $1$).
```


```python
print('eigenvalores:')
print(evalue)
print('eigenvectores:')
print(evector)
```

    eigenvalores:
    [10. 10. 11.]
    eigenvectores:
    [[ 1.    -1.     0.998]
     [ 0.     0.     0.056]
     [ 0.     0.     0.006]]


### Ejemplo

Un eigenvalor puede estar repetido y tener un sólo eigenvector linealmente independiente: 


```python
A = np.array([[2, 1, 0],
              [0, 2, 1],
              [0, 0, 2]])
```


```python
evalue, evector = np.linalg.eig(A)
```

```{margin}

Observa que en este ejemplo el eigenvalor $2$ está repetido tres veces (multiplicidad algebraica es igual a $3$) y sus eigenvectores asociados son linealmente dependientes (multiplicidad geométrica es igual a $1$).
```


```python
print('eigenvalores:')
print(evalue)
print('eigenvectores:')
print(evector)
```

    eigenvalores:
    [2. 2. 2.]
    eigenvectores:
    [[ 1. -1.  1.]
     [ 0.  0. -0.]
     [ 0.  0.  0.]]


```{admonition} Definición

Si $(\lambda, v)$ es una pareja de eigenvalor-eigenvector de $A$ tales que $Av = \lambda v$ entonces $v$ se le nombra eigenvector derecho. Si $(\lambda, v)$ es una pareja de eigenvalor-eigenvector de $A^T$ tales que $A^Tv = \lambda v$ (que es equivalente a $v^TA=\lambda v^T$) entonces $v$ se le nombra eigenvector izquierdo.
```

```{admonition} Observaciones
:class: tip

* En todos los ejemplos anteriores se calcularon eigenvectores derechos.

* Los eigenvectores izquierdos y derechos para una matriz simétrica son iguales.

```

(DIAGONALIZABLE)=

## $A$ diagonalizable

```{admonition} Definición

Si $A$ tiene $n$ eigenvectores linealmente independientes entonces $A$ se nombra diagonalizable o *non defective*. En este caso si $x_1, x_2, \dots, x_n$ son eigenvectores de $A$ con $Ax_i = \lambda_i x_i$ para $i=1,\dots,n$ entonces la igualdad anterior se escribe en ecuación matricial como:

$$AX = X \Lambda$$

o bien:

$$A = X \Lambda X^{-1}$$

donde: $X$ tiene por columnas los eigenvectores de $A$ y $\Lambda$ tiene en su diagonal los eigenvalores de $A$.

A la descomposición anterior $A = X \Lambda X^{-1}$ para $A$ diagonalizable o *non defective* se le nombra ***eigen decomposition***.
```

```{admonition} Observación
:class: tip

* Si $A = X \Lambda X^{-1}$ entonces $X^{-1}A = \Lambda X^{-1}$ y los renglones de $X^{-1}$ (o equivalentemente las columnas de $X^{-T}$) son eigenvectores izquierdos.

* Si $A = X \Lambda X^{-1}$ y $b = Ax = (X \Lambda X^{-1}) x$ entonces:

$$\tilde{b} = X^{-1}b = X^{-1}  (Ax) = X^{-1} (X \Lambda X^{-1}) x = \Lambda X^{-1}x = \Lambda \tilde{x}.$$

Lo anterior indica que el producto matricial $Ax$ para $A$ diagonalizable es equivalente a multiplicar una matriz diagonal por un vector denotado como $\tilde{x}$ que contiene los coeficientes de la combinación lineal de las columnas de $X$ para el vector $x$ . El resultado de tal multiplicación es un vector denotado como $\tilde{b}$ que también contiene los coeficientes de la combinación lineal de las columnas de $X$ para el vector $b$. En resúmen, si $A$ es diagonalizable o *non defective* la multiplicación $Ax$ es equivalente a la multiplicación por una matriz diagonal $\Lambda \tilde{x}$ (salvo un cambio de bases, ver [Change of basis](https://en.wikipedia.org/wiki/Change_of_basis)).

* Si una matriz $A$ tiene eigenvalores distintos entonces es diagonalizable y más general: si $A$ tiene una multiplicidad geométrica igual a su multiplicidad algebraica de cada eigenvalor entonces es diagonalizable. 

```

### Ejemplo

La matriz:

$$A = \left[
\begin{array}{ccc}
1 & -4 & -4\\
8 & -11 & -8\\
-8 & 8 & 5
\end{array}
\right]
$$

es diagonalizable.


```python
A = np.array([[1, -4, -4],
              [8, -11, -8],
              [-8, 8, 5.0]])
```


```python
print(A)
```




    array([[  1.,  -4.,  -4.],
           [  8., -11.,  -8.],
           [ -8.,   8.,   5.]])




```python
evalue, evector = np.linalg.eig(A)
```


```python
print('eigenvalores:')
print(evalue)
```

    eigenvalores:
    [ 1. -3. -3.]


```{margin}

Se verifica que los eigenvectores de este ejemplo es un conjunto linealmente independiente por lo que $A=X\Lambda X^{-1}$.

```


```python
print('eigenvectores:')
print(evector)
```

    eigenvectores:
    [[ 0.333 -0.717 -0.241]
     [ 0.667 -0.02  -0.796]
     [-0.667 -0.697  0.555]]



```python
X = evector
Lambda = np.diag(evalue)
```

```{margin}

Observa que si $Z$ es desconocida y $X^T Z^T = \Lambda$ entonces $Z^T = X^{-T} \Lambda$ y por tanto $XZ =X\Lambda X^{-1}$.
```


```python
X@np.linalg.solve(X.T, Lambda).T
```




    array([[  1.,  -4.,  -4.],
           [  8., -11.,  -8.],
           [ -8.,   8.,   5.]])




```python
print(A)
```




    array([[  1.,  -4.,  -4.],
           [  8., -11.,  -8.],
           [ -8.,   8.,   5.]])



$A$ es diagonalizable pues: $X^{-1} A X = \Lambda$

```{margin}

Observa que si $Z$ es desconocida y $XZ = A$ entonces $Z = X^{-1}A$ y por tanto $ZX = X^{-1} A X$.
```


```python
print(np.linalg.solve(X, A)@X)
```

    [[ 1.  0. -0.]
     [-0. -3. -0.]
     [ 0. -0. -3.]]



```python
print(Lambda)
```




    array([[ 1.,  0.,  0.],
           [ 0., -3.,  0.],
           [ 0.,  0., -3.]])



```{admonition} Observación
:class: tip

Observa que **no necesariamente** $X$ en la *eigen decomposition* es una matriz ortogonal.

```

```{margin}

Aquí se toma $X[1:3,1]$ como la primera columna de $X$ y se satisface $X[1:3,1]^TX[1:3,1] = 1$ en este ejemplo pero en general esto no se cumple.

```


```python
X[:,0].dot(X[:,0])
```




    1.0



```{margin}

$X[1:3,1]^TX[1:3,2] \neq 0$ por lo que la primera y segunda columna de $X$ no son ortogonales.

```


```python
X[:,0].dot(X[:,1])
```




    0.2117366284008178



**Eigenvectores derechos:**

```{margin}

`x_1` es la primer columna de $X$: $X[1:3, 1]$ y `lambda_1` el eigenvalor asociado.
```


```python
x_1 = X[:,0]
lambda_1 = Lambda[0,0]
```


```python
print(A@x_1)
```




    array([ 0.333,  0.667, -0.667])



```{margin}

$Ax_1 = \lambda_1 x_1$.
```


```python
print(lambda_1*x_1)
```




    array([ 0.333,  0.667, -0.667])



```{margin}

`x_2` es la segunda columna de $X$: $X[1:3, 2]$ y `lambda_2` el eigenvalor asociado.
```


```python
x_2 = X[:,1]
lambda_2 = Lambda[1,1]
```


```python
print(A@x_2)
```




    array([2.151, 0.061, 2.09 ])



```{margin}

$Ax_2 = \lambda_2 x_2$.
```


```python
print(lambda_2*x_2)
```




    array([2.151, 0.061, 2.09 ])



**Eigenvectores izquierdos:**

```{admonition} Observación
:class: tip

Para los eigenvectores izquierdos se deben tomar los renglones de $X^{-1}$ (o equivalentemente las columnas de $X^{-T}$) sin embargo no se utiliza el método [inv](https://numpy.org/doc/stable/reference/generated/numpy.linalg.inv.html) de *NumPy* pues es más costoso computacionalmente y amplifica los errores por redondeo. En su lugar se utiliza el método [solve](https://numpy.org/doc/stable/reference/generated/numpy.linalg.solve.html) y se resuelve el sistema: $X^{-T} z  = e_i$ para $e_i$ $i$-ésimo vector canónico.

```


```python
e1 = np.zeros((X.shape[0],1))
```


```python
e1[0] = 1
```


```python
print(e1)
```




    array([[1.],
           [0.],
           [0.]])



```{margin}

`x_inv_1` es el primer renglón de $X^{-1}$: $X^{-1}[1, 1:3]$.
```


```python
x_inv_1 = np.linalg.solve(X.T, e1)
```


```python
print(x_inv_1)
```




    array([[ 3.],
           [-3.],
           [-3.]])




```python
print(A.T@x_inv_1)
```




    array([[ 3.],
           [-3.],
           [-3.]])



```{margin}

$A^TX^{-T}[1:3,1] = \lambda_1 X^{-T}[1:3,1]$, `lambda_1` el eigenvalor asociado a `x_inv_1`.
```


```python
print(lambda_1*x_inv_1)
```




    array([[ 3.],
           [-3.],
           [-3.]])




```python
e2 = np.zeros((X.shape[0],1))
```


```python
e2[1] = 1
```

```{margin}

`x_inv_2` es el segundo renglón de $X^{-1}$: $X^{-1}[2, 1:3]$.
```


```python
x_inv_2 = np.linalg.solve(X.T, e2)
```


```python
print(x_inv_2)
```




    array([[-0.851],
           [-0.13 ],
           [-0.556]])




```python
print(A.T@x_inv_2)
```




    array([[2.552],
           [0.391],
           [1.667]])



```{margin}

$A^TX^{-T}[1:3,2] = \lambda_2 X^{-T}[1:3,2]$, `lambda_2` el eigenvalor asociado a `x_inv_2`.
```


```python
print(lambda_2*x_inv_2)
```




    array([[2.552],
           [0.391],
           [1.667]])



```{admonition} Ejercicio
:class: tip

¿Es la siguiente matriz diagonalizable?

$$A = \left [
\begin{array}{ccc}
-1 & -1 & -2\\
8 & -11 & -8\\
-10 & 11 & 7
\end{array}
\right]
$$

si es así encuentra su *eigen decomposition* y diagonaliza a $A$.
```

(DESCESP)=

### Resultado: $A$ simétrica

Si A es simétrica entonces tiene eigenvalores reales. Aún más: $A$ tiene eigenvectores reales linealmente independientes, forman un conjunto ortonormal y se escribe como un producto de tres matrices nombrado **descomposición espectral o *symmetric eigen decomposition***: 

$$A = Q \Lambda Q^T$$ 

donde: $Q$ es una matriz ortogonal cuyas columnas son eigenvectores de $A$ y $\Lambda$ es una matriz diagonal con eigenvalores de $A$.


```{admonition} Comentarios

* Por lo anterior una matriz simétrica es **ortogonalmente diagonalizable**, ver {ref}`A diagonalizable <DIAGONALIZABLE>`.

* Los eigenvalores de $A$ simétrica se pueden ordenar:

$$\lambda_n(A) \leq \lambda_{n-1}(A) \leq \dots \leq \lambda_1(A)$$

con:

$\lambda_{max}(A) = \lambda_1(A)$, $\lambda_{min}(A) = \lambda_n(A)$.

* Se prueba para $A$ simétrica:

$$\lambda_{max}(A) = \displaystyle \max_{x \neq 0} \frac{x^TAx}{x^Tx}$$

$$\lambda_{min}(A) = \displaystyle \min_{x \neq 0} \frac{x^TAx}{x^Tx}.$$

por lo tanto:

$$\lambda_{min}(A) \leq \frac{x^TAx}{x^Tx} \leq \lambda_{max}(A) \forall x \neq 0.$$

* $||A||_2 = \displaystyle \max\{|\lambda_1(A)|, |\lambda_n(A)|\}$.

* $||A||_F = \left( \displaystyle \sum_{i=1}^n \lambda_i ^2 \right)^{1/2}$.

* Los valores singulares de $A$ son el conjunto $\{|\lambda_1(A)|, \dots, |\lambda_{n-1}(A)|, |\lambda_n(A)|\}$.
```

### Ejemplo

Matriz simétrica y descomposición espectral de la misma:


```python
A=np.array([[5,4,2],[4,5,2],[2,2,2]])
```


```python
print(A)
```

    [[5 4 2]
     [4 5 2]
     [2 2 2]]



```python
evalue, evector = np.linalg.eigh(A)
```

```{margin}

Como $A$ es simétrica sus eigenvalores son reales y sus eigenvectores forman un conjunto linealmente independiente. Por lo anterior $A$ tiene descomposción espectral.
```


```python
print('eigenvalores:')
print(evalue)
print('eigenvectores:')
print(evector)
```

    eigenvalores:
    [ 1.  1. 10.]
    eigenvectores:
    [[ 0.482  0.569  0.667]
     [-0.727 -0.166  0.667]
     [ 0.49  -0.806  0.333]]


```{margin}

$A = Q \Lambda Q^T$
```


```python
print('descomposición espectral:')
Lambda = np.diag(evalue)
Q = evector
print('QLambdaQ^T:')
print(Q@Lambda@Q.T)
print('A:')
print(A)
```

    descomposición espectral:
    QLambdaQ^T:
    [[5. 4. 2.]
     [4. 5. 2.]
     [2. 2. 2.]]
    A:
    [[5 4 2]
     [4 5 2]
     [2 2 2]]


A es diagonalizable pues: $Q^T A Q = \Lambda$


```python
print(Q.T@A@Q)
```

    [[ 1.  0.  0.]
     [-0.  1. -0.]
     [ 0.  0. 10.]]



```python
print(Lambda)
```




    array([[ 1.,  0.,  0.],
           [ 0.,  1.,  0.],
           [ 0.,  0., 10.]])



Ver [numpy.linalg.eigh](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.linalg.eigh.html).

## Condición del problema del cálculo de eigenvalores y eigenvectores

La condición del problema del cálculo de eigenvalores y eigenvectores de una matriz, es la sensibilidad de los mismos ante perturbaciones en la matriz, ver {ref}`Condición de un problema y estabilidad de un algoritmo <CPEA>`. Diferentes eigenvalores o eigenvectores de una matriz no necesariamente son igualmente sensibles a perturbaciones en la matriz.

```{admonition} Observación
:class: tip

La condición del problema del cálculo de eigenvalores y eigenvectores de una matriz **no** es igual a la condición del problema de resolver un sistema de ecuaciones lineales, ver {ref}`Número de condición de una matriz <NCM>`.

```

Se prueba que la condición de un eigenvalor **simple** de una matriz $A$ está dado por $\frac{1}{|y^Tx|}$ con $x$ eigenvector derecho, $y$ eigenvector izquierdo de $A$ ambos asociados al eigenvalor simple y normalizados esto es: $x^Tx = y^Ty=1$. 

```{admonition} Comentarios

* Para los casos en que: $\lambda$ eigenvalor de $A$ sea simple, $A$ sea diagonalizable, existen eigenvectores izquierdos y derechos asociados a un eigenvalor de $A$ tales que $y^Tx \neq 0$. En tales casos, el análisis del condicionamiento del problema del cálculo de eigenvalores y eigenvectores es más sencillo de realizar que para matrices no diagonalizables o eigenvalores con multiplicidad algebraica mayor a $1$. En particular, los eigenvalores de una matriz simétrica están muy bien condicionados: las perturbaciones en $A$ únicamente perturban a los eigenvalores en una magnitud medida con la norma de las perturbaciones y no depende de otros factores, por ejemplo del número de condición de $A$.

* La sensibilidad de un eigenvector depende de la sensibilidad de su eigenvalor asociado y de la distancia de tal eigenvalor de otros eigenvalores.

* Los eigenvalores que son "cercanos" o aquellos de multiplicidad mayor a $1$ pueden ser mal condicionados y por lo tanto difíciles de calcularse de forma exacta y precisa en especial si la matriz es defectuosa (no diagonalizable). Puede mejorarse el número de condición si se escala el problema por una matriz diagonal y similar a $A$, ver {ref}`similitud <SIMILITUD>`.
```

(SIMILITUD)=

## Similitud

```{admonition} Definición

Si existe $X \in \mathbb{R}^{n \times n}$ tal que $B = XAX^{-1}$ con $A, B \in \mathbb{R}^{n \times n}$ entonces $A$ y $B$ se nombran similares.

```

```{admonition} Observación
:class: tip

Las matrices que son similares tienen el mismo espectro, de hecho: $Ax = \lambda x$ si y sólo si $By = \lambda y$ para $y=Xx$. Lo anterior quiere decir que los eigenvalores de una matriz son **invariantes** ante cambios de bases o representación en coordenadas distintas.

```

### Ejemplo

Dada la matriz 

$$A=
\left [
\begin{array}{cccc}
-1 & -1 & -1 & -1\\
0 & -5 & -16 & -22\\
0 & 3 & 10 & 14\\
4 & 8 & 12 & 14
\end{array}
\right ]
$$

Definir matrices $B_1, B_2$ similares a $A$ a partir de las matrices:

$$
\begin{array}{l}
X_1 = 
\left [
\begin{array}{cccc}
2 & -1 & 0 & 0\\
-1 & 2 & -1 & 0\\
0 & -1 & 2 & -1\\
0 & 0 & -1 & 1
\end{array}
\right ],
X_2 = \left [
\begin{array}{cccc}
2 & -1 & 1 & 0\\
-1 & 2 & 0 & 0\\
0 & -1 & 0 & 0\\
0 & 0 & 0 & 1
\end{array}
\right ]
\end{array}
$$

y verificar que los eigenvalores de $A$ son los mismos que los de $B_1, B_2$, esto es, tienen el mismo espectro.


```python
A = np.array([[-1, -1 , -1, -1],
              [0, -5, -16, -22],
              [0, 3, 10, 14],
              [4, 8, 12, 14.0]])
```


```python
X1 = np.array([[2, -1, 0, 0],
              [-1, 2, -1, 0],
              [0, -1, 2, -1],
              [0, 0, -1, 1.0]])
```

$B_1 = X_1^{-1}AX_1$:

```{margin}

Calculamos $B1$ explícitamente para revisar qué forma tiene pero no es necesario.

```


```python
B1 = np.linalg.solve(X1, A)@X1

```


```python
print(B1)
```

    [[ 1.  2. -0.  0.]
     [ 3.  4. -0.  0.]
     [ 0.  0.  5.  6.]
     [ 0.  0.  7.  8.]]



```python
X2 = np.array([[2, -1, 1, 0],
               [-1, 2, 0, 0],
               [0, -1, 0, 0],
               [0, 0, 0, 1.0]])
```

$B_2 = X_2^{-1}AX_2$:

```{margin}

Calculamos $B2$ explícitamente para revisar qué forma tiene pero no es necesario.

```


```python
B2 = np.linalg.solve(X2, A)@X2

```


```python
print(B2)
```

    [[  1.   2.   0.  -6.]
     [  3.   4.   0. -14.]
     [ -0.   0.  -1.  -3.]
     [  0.   0.   4.  14.]]


**$B1$ y $B2$ son similares a $A$ por tanto tienen los mismos eigenvalores:**


```python
evalue, evector = np.linalg.eig(A)
```

```{margin}

`evalue` son los eigenvalores de $A$.

```


```python
print(evalue)
```

    [13.152  5.372 -0.152 -0.372]



```python
evalue_B1, evector_B1 = np.linalg.eig(B1)
```

```{margin}

`evalue_B1` son los eigenvalores de $B_1$, obsérvese que son iguales a los de $A$ salvo el orden.

```


```python
print(evalue_B1)
```

    [-0.372  5.372 -0.152 13.152]



```python
evalue_B2, evector_B2 = np.linalg.eig(B2)
```

```{margin}

`evalue_B2` son los eigenvalores de $B_2$, obsérvese que son iguales a los de $A$ salvo el orden.

```


```python
print(evalue_B2)
```

    [ 5.372 -0.372 -0.152 13.152]


Los eigenvectores **no son los mismos** pero pueden obtenerse vía multiplicación de matrices:

```{margin}

Elegimos un eigenvalor de $A$.

```


```python
print(evalue[1])
```




    5.372281323269013



```{margin}

Y elegimos el mismo eigenvalor en el *array* `evalue_B1` para $B_1$ que para este ejemplo corresponde al índice $1$ (el mismo que en `evalue` pero podría haber sido otro).

```


```python
print(evalue_B1[1])
```




    5.3722813232690125



```{margin}

Su correspondiente eigenvector en el índice $1$ del *array* `evector_B1`.

```


```python
print(evector_B1[:,1])
```




    array([-0.416, -0.909,  0.   ,  0.   ])



**$X^{-1}x$ es eigenvector de $B_1$ para $x$ eigenvector de $A$**:

```{margin}

`evector[:,1]` es el eigenvector de $A$ correspondiente al eigenvalor `evalue[1]`. En esta celda se hace el producto $X_1^{-1}x$ y `evector[:,1]` representa a $x$.

```


```python
X1_inv_evector = np.linalg.solve(X1, evector[:,1])
```


```python
print(X1_inv_evector)
```




    array([ 0.249,  0.543, -0.   , -0.   ])




```python
print(B1@(X1_inv_evector))
```




    array([ 1.335,  2.919, -0.   , -0.   ])



```{margin}

Se verifica que $B1(X_1^{-1}x) = \lambda (X_1^{-1}x)$ con $\lambda$ igual al valor `evalue_B1[1]`.

```


```python
print(evalue_B1[1]*(X1_inv_evector))
```




    array([ 1.335,  2.919, -0.   , -0.   ])



```{admonition} Observación
:class: tip

Obsérvese que son los mismos eigenvectores salvo una constante distinta de cero.

```


```python
print(evector_B1)
```

    [[-0.825 -0.416  0.     0.   ]
     [ 0.566 -0.909 -0.     0.   ]
     [ 0.     0.    -0.759 -0.593]
     [ 0.     0.     0.651 -0.805]]


```{margin}

El valor `1.33532534` es la primera entrada de `X1_inv_evector` que es $X_1^{-1}x$ y $x$ eigenvector de $A$. El valor `2.91920903` es la segunda entrada de `X_1_inv_evector`. Las entradas restantes son cercanas a cero.
```


```python
print(1.33532534e+00/evector_B1[0,1])
```




    -3.2101207266138467




```python
print(2.91920903e+00/evector_B1[1,1])
```




    -3.2101207350977647



La constante es aproximadamente $-3.21$:

```{margin}

`evector_B1` fue calculado con la función `eig` pero en la siguiente celda se observa que no es necesario si se tiene un eigenvector de $A$.
```


```python
print(evector_B1[:,1]*(-3.21))
```




    array([ 1.335,  2.919, -0.   , -0.   ])



```{margin}

Recuerda que `X_1_inv_evector` es $X_1^{-1}x$ con $x$ eigenvector de $A$ que en este caso se utilizó `evector[:,1]`.

```


```python
print(B1@(X1_inv_evector))
```




    array([ 1.335,  2.919, -0.   , -0.   ])



```{margin}

Se comprueba que $X_1^{-1}x$ es eigenvector de $B$ si $x$ es eigenvector de $A$.
```


```python
print(evalue_B1[1]*(X1_inv_evector))
```




    array([ 1.335,  2.919, -0.   , -0.   ])



Como $A$ tiene eigenvalores distintos entonces es diagonalizable, esto es existen $X_3, \Lambda$ tales que $X_3^{-1} A X_3 = \Lambda$.


```python
X_3 = evector
Lambda = np.diag(evalue)
```


```python
print(A)
```




    array([[ -1.,  -1.,  -1.,  -1.],
           [  0.,  -5., -16., -22.],
           [  0.,   3.,  10.,  14.],
           [  4.,   8.,  12.,  14.]])




```python
print(np.linalg.solve(X_3, A)@X_3)
```

    [[13.152 -0.     0.    -0.   ]
     [ 0.     5.372 -0.     0.   ]
     [-0.     0.    -0.152  0.   ]
     [-0.     0.    -0.    -0.372]]



```python
print(Lambda)
```




    array([[13.152,  0.   ,  0.   ,  0.   ],
           [ 0.   ,  5.372,  0.   ,  0.   ],
           [ 0.   ,  0.   , -0.152,  0.   ],
           [ 0.   ,  0.   ,  0.   , -0.372]])



```{admonition} Comentario

**$X_1$ diagonaliza a $A$ por bloques, $X_2$ triangulariza a $A$ por bloques y $X_3$ diagonaliza a $A$.** Las tres matrices representan al mismo operador lineal (que es una transformación lineal del espacio vectorial sobre sí mismo) pero en **coordenadas diferentes**. Un aspecto muy **importante** en el álgebra lineal es representar a tal operador lineal en unas coordenadas lo más simple posible. En el ejemplo la matriz $X_3$, que en sus columnas están los eigenvectores de $A$, ayuda a representarlo de forma muy simple.

```

```{admonition} Observación
:class: tip

$X_3$ es una matriz que diagonaliza a $A$ y tiene en sus columnas a eigenvectores de $A$, si el objetivo es diagonalizar a una matriz **no es necesario** resolver un problema de cálculo de eigenvalores-eigenvectores pues cualquier matriz $X$ no singular puede hacer el trabajo. Una opción es considerar una factorización para $A$ simétrica del tipo $LDL^T$ (que tiene un costo computacional bajo para calcularse), la matriz $L$ no es ortogonal y la matriz $D$ tiene los pivotes que se calculan en la eliminación Gaussiana, ver {ref}` Operaciones y transformaciones básicas del Álgebra Lineal Numérica <OTBALN>`.

```

```{admonition} Ejercicio
:class: tip

Considera 

$$A=
\left [
\begin{array}{cccc}
-2 & -1 & -5 & 2\\
-9 & 0 & -8 & -2\\
2 & 3 & 11 & 5\\
3 & -5 & -13 & -7
\end{array}
\right ]
$$

Define $X_1$ tal que $X_1^{-1}AX_1$ sea diagonal.
```

### Ejemplo


```python
import sympy
import matplotlib.pyplot as plt
```

Considérese la siguiente ecuación cuadrática:

$$\frac{\tilde{x}^2}{16} + \frac{\tilde{y}^2}{9} = 1.$$

Con Geometría Analítica sabemos que tal ecuación representa una elipse inclinada. Además si



```python
D = sympy.Matrix([[sympy.Rational(1,16), 0],
                  [0, sympy.Rational(1,9)]])
```


```python
print(D)
```




$\displaystyle \left[\begin{matrix}\frac{1}{16} & 0\\0 & \frac{1}{9}\end{matrix}\right]$



Entonces el producto 

$$\left [ \begin{array}{c}
\tilde{x}\\
\tilde{y}
\end{array}
\right ] ^TD 
\left [
\begin{array}{c}
\tilde{x}\\
\tilde{y}
\end{array}
\right ]
$$

es:


```python
x_tilde, y_tilde = sympy.symbols("x_tilde, y_tilde")
x_y_tilde = sympy.Matrix([x_tilde, y_tilde])
```


```python
print((x_y_tilde.T*D*x_y_tilde)[0])
```




$\displaystyle \frac{x_{tilde}^{2}}{16} + \frac{y_{tilde}^{2}}{9}$



Rotemos al [eje mayor de la elipse](https://en.wikipedia.org/wiki/Semi-major_and_semi-minor_axes) un ángulo de $\theta = \frac{\pi}{3}$ con una {ref}`transformación de rotación <TROT>` que genera la ecuación matricial:


$$\begin{array}{l}
\left[
\begin{array}{c}
x\\
y
\end{array}
\right ]
=
\left [
\begin{array}{cc}
\cos(\theta) & -\sin(\theta)\\
\sin(\theta) & \cos(\theta)
\end{array}
\right ]
\left[
\begin{array}{c}
\tilde{x}\\
\tilde{y}
\end{array}
\right ]
=
\left [
\begin{array}{cc}
\frac{1}{2} & -\frac{\sqrt{3}}{2}\\
\frac{\sqrt{3}}{2} & \frac{1}{2}
\end{array}
\right ]
\left[
\begin{array}{c}
\tilde{x}\\
\tilde{y}
\end{array}
\right ]
=
Q\left[
\begin{array}{c}
\tilde{x}\\
\tilde{y}
\end{array}
\right ]
\end{array}
$$

donde: $Q$ es la matriz de rotación en sentido contrario a las manecillas del reloj por el ángulo $\theta$.

Esto es:

$$
\begin{eqnarray}
x =\frac{\tilde{x}}{2} - \frac{\tilde{y}\sqrt{3}}{2} \nonumber \\
y =\frac{\tilde{x}\sqrt{3}}{2} + \frac{\tilde{y}}{2} \nonumber
\end{eqnarray}
$$

Despejando $\tilde{x},\tilde{y}$ y sustituyendo en $\frac{\tilde{x}^2}{16} + \frac{\tilde{y}^2}{9} = 1$ resulta en la ecuación:



```{margin}

Esta es una ecuación de una elipse inclinada.

```


```python
theta = sympy.pi/3
Q = sympy.Matrix([[sympy.cos(theta), -sympy.sin(theta)],
                  [sympy.sin(theta), sympy.cos(theta)]])
x,y = sympy.symbols("x, y")
x_tilde = (Q.T*sympy.Matrix([x,y]))[0]
y_tilde = (Q.T*sympy.Matrix([x,y]))[1]

print((x_tilde**2/16 + y_tilde**2/9).expand()*576)

```




$\displaystyle 57 x^{2} - 14 \sqrt{3} x y + 43 y^{2}$



O equivalentemente el producto 

$$\left [ \begin{array}{c}
x\\
y
\end{array}
\right ]^T A 
\left [
\begin{array}{c}
x\\
y
\end{array}
\right ]
$$ 



```python
x_y = sympy.Matrix([x,y])
A = Q*D*Q.T
print(((x_y.T*A*x_y)[0]).expand()*576)
```




$\displaystyle 57 x^{2} - 14 \sqrt{3} x y + 43 y^{2}$



con $A$ matriz dada por $A=QDQ^T$:

```{margin}

Observa que $A$ es **simétrica**.

```


```python
print(A)
```




$\displaystyle \left[\begin{matrix}\frac{19}{192} & - \frac{7 \sqrt{3}}{576}\\- \frac{7 \sqrt{3}}{576} & \frac{43}{576}\end{matrix}\right]$



En este ejemplo la matriz $Q$ de rotación es la matriz que diagonaliza ortogonalmente a $A$ pues: $Q^TAQ = D.$


Para realizar la **gráfica** de la elipse con *NumPy* observar que:

```{margin}

Estas ecuaciones nos indican que la misma elipse se puede representar en diferentes coordenadas. El cambio de coordenadas se realiza con la matriz $Q$.

```

$$
\begin{eqnarray}
1&=&57x^2 - 14 \sqrt{3}xy + 43 y^2 \nonumber \\
&=& \left [ \begin{array}{c}
x\\
y
\end{array}
\right ]^T A 
\left [
\begin{array}{c}
x\\
y
\end{array}
\right ] \nonumber \\
&=& \left [ \begin{array}{c}
x\\
y
\end{array}
\right ]^T QDQ^T \left [
\begin{array}{c}
x\\
y
\end{array}
\right ] \nonumber \\
&=& \left(Q^T \left [ \begin{array}{c}
x\\
y
\end{array}
\right ]\right)^TD\left(Q^T \left [ \begin{array}{c}
x\\
y
\end{array}
\right ]\right) \nonumber \\
&=& \left [ \begin{array}{c}
\tilde{x}\\
\tilde{y}
\end{array}
\right ] ^TD 
\left [
\begin{array}{c}
\tilde{x}\\
\tilde{y}
\end{array}
\right ] \nonumber \\
&=& \frac{\tilde{x}^2}{16} + \frac{\tilde{y}^2}{9} \nonumber
\end{eqnarray}
$$

**Gráfica para eigenvalores ordenados de forma decreciente en la diagonal de la matriz $D$.**

```{margin}

Usamos [eig](https://numpy.org/doc/stable/reference/generated/numpy.linalg.eig.html) para el cálculo numérico de eigenvalores, eigenvectores de $A$.

```


```python
A_np = np.array(A.evalf(),dtype = float)
evalue, evector = np.linalg.eig(A_np)
Q_np = evector
D_np = np.diag([evalue[0], evalue[1]])
```


```python
print(Q_np)
```




    array([[ 0.866,  0.5  ],
           [-0.5  ,  0.866]])




```python
Q = sympy.Matrix([Q[1,:],-Q[0,:]])
print(Q)
```




$\displaystyle \left[\begin{matrix}\frac{\sqrt{3}}{2} & \frac{1}{2}\\- \frac{1}{2} & \frac{\sqrt{3}}{2}\end{matrix}\right]$




```python
print(D_np)
```




    array([[0.111, 0.   ],
           [0.   , 0.062]])




```python
print(sympy.Matrix([[D[1,1], 0],
              [0, D[0,0]]]))
```




$\displaystyle \left[\begin{matrix}\frac{1}{9} & 0\\0 & \frac{1}{16}\end{matrix}\right]$




```python
d1_inv=1/4
d2_inv=1/3
density=1e-2
x=np.arange(-1/d1_inv,1/d1_inv,density)
y1=1/d2_inv*np.sqrt(1-(d1_inv*x)**2)
y2=-1/d2_inv*np.sqrt(1-(d1_inv*x)**2)
#transform
x_y1_hat = np.column_stack((x,y1))
x_y2_hat = np.column_stack((x,y2))
apply_Q = lambda vec : np.transpose(Q@np.transpose(vec))
Q_to_vector_1 = apply_Q(x_y1_hat)
Q_to_vector_2 = apply_Q(x_y2_hat)
fig = plt.figure(figsize=(10, 5))
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
#first plot
ax1.plot(Q_to_vector_1[:,0],Q_to_vector_1[:,1],'g',
         Q_to_vector_2[:,0],Q_to_vector_2[:,1],'g')
ax1.set_title("57x^2-14(3^(1/2))xy+43y^2=1")
ax1.set_xlabel("Ejes coordenados típicos")
ax1.axhline(color='r')
ax1.axvline(color='r')
#second plot
evector_1 = 1/d1_inv*Q_np[:,0]
evector_2 = 1/d2_inv*Q_np[:,1]
Evector_1 = np.row_stack((np.zeros(2), evector_1))
Evector_2 = np.row_stack((np.zeros(2), evector_2))
ax2.plot(Q_to_vector_1[:,0],Q_to_vector_1[:,1],
         color='g', label = "Elipse")
ax2.plot(Q_to_vector_2[:,0],Q_to_vector_2[:,1],
         color='g', label = "_nolegend_")
ax2.plot(Evector_1[:,0], Evector_1[:,1],
         color='b', label = "Eigenvector Q[:,0], define al semieje mayor principal de la elipse")
ax2.plot(-Evector_1[:,0], -Evector_1[:,1],
         color='b', label = "_nolegend_")
ax2.plot(Evector_2[:,0], Evector_2[:,1],
         color='m', label = "Eigenvector Q[:,1],  define al semieje menor principal de la elipse")
ax2.plot(-Evector_2[:,0], -Evector_2[:,1],
         color='m', label = "_nolegend_")
ax2.set_title("x_tilde^2(1/16) + y_tilde^2(1/9)=1")
ax2.set_xlabel("Ejes coordenados rotados")
ax2.legend(bbox_to_anchor=(1, 1))
fig.suptitle("Puntos en el plano que cumplen x_y_tilde^TDx_y_tilde = 1 o  x_y^TAx_y =1")
ax2.grid()
plt.show()
```

```{margin}

Recuerda que $A = Q D Q^T$, $A$ es similar a $D$ matriz diagonal y $Q$ es ortogonal.

```

En la gráfica anterior se representa la rotación de los ejes coordenados definidos por los vectores canónicos $e_1, e_2$ y los rotados definidos por los eigenvectores de $A$. Los eigenvectores de $A$ están en las columnas de $Q$. La primera columna de $Q$ define al eje mayor principal de la elipse y la segunda columna al eje menor principal. La longitud de los semiejes están dados respectivamente por la raíz cuadrada de los recíprocos de los eigenvalores de $A$ que en este caso son: $\frac{1}{9}, \frac{1}{16}$, esto es: $4$ y $16$. Ver por ejemplo: [Principal_axis_theorem](https://en.wikipedia.org/wiki/Principal_axis_theorem), [Diagonalizable_matrix](https://en.wikipedia.org/wiki/Diagonalizable_matrix).

```{admonition} Ejercicio
:class: tip

Rotar los ejes coordenados $45^o$ la ecuación de la elipse:

$$13x^2+10xy+13y^2=72$$

para representar tal ecuación alineando los ejes mayor y menor de la elipse a sus eigenvectores. Encontrar las matrices $Q, D$ tales que $A=QDQ^T$ con $Q$ ortogonal y $D$ diagonal.

```

## Algunos algoritmos para calcular eigenvalores y eigenvectores

Dependiendo de las siguientes preguntas es el tipo de algoritmo que se utiliza:

* ¿Se requiere el cómputo de todos los eigenvalores o de sólo algunos?

* ¿Se requiere el cómputo de únicamente los eigenvalores o también de los eigenvectores?

* ¿$A$ tiene entradas reales o complejas?

* ¿$A$ es de dimensión pequeña y es densa o grande y rala?

* ¿$A$ tiene una estructura especial o es una matriz general?

Para la última pregunta a continuación se tiene una tabla que resume las estructuras en las matrices que son relevantes para problemas del cálculo de eigenvalores-eigenvectores:

|Estructura|Definición|
|:---:|:---:|
|Simétrica|$A=A^T$|
|Ortogonal|$A^TA=AA^T=I_n$|
|Normal|$A^TA = AA^T$|

Ver {ref}`Ejemplos de matrices normales <EJMN>`.

(EJMN)=

### Una opción (inestable numéricamente respecto al redondeo): encontrar raíces del polinomio característico... 

Por definición, los eigenvalores de $A \in \mathbb{R}^{n \times n}$ son las raíces o ceros del polinomio característico $p(z)$ por lo que un método es calcularlas vía tal polinomio. Sin embargo, **no es un buen método** calcular tales raíces o ceros pues para una $n > 4$ [Abel](https://en.wikipedia.org/wiki/Abel%E2%80%93Ruffini_theorem) probó de forma teórica que las raíces en general no son posibles expresarlas por una fórmula cerrada que involucren los coeficientes, operaciones aritméticas y raíces $\sqrt[n]{\cdot}$ . Por lo anterior para calcular eigenvalores de matrices con dimensión $n>4$ requiere de un **método iterativo**.

```{margin}

Como ejemplo de este enunciado considérese:

$$A=\left[
\begin{array}{cc}
1 & \epsilon\\
\epsilon & 1\\
\end{array}
\right]
$$

cuyos eigenvalores son $1 + \epsilon$, $1 - \epsilon$ con $\epsilon$ menor que $\epsilon_{maq}$. Usando aritmética en el SPF se prueba que las raíces del polinomio característico es $1$ de multiplicidad $2$.
```

Además de lo anterior, en ciertas bases de polinomios, por ejemplo $\{1, x, x^2, \dots, x^n\}$, los coeficientes de los polinomios  numéricamente no están bien determinados por los errores por redondeo y las raíces de los polinomios son muy sensibles a perturbaciones de los coeficientes, esto es, es un problema mal condicionado, ver {ref}`condición de un problema y estabilidad de un algoritmo <CPEA>`. Ver [Wilkinson's polynomial](https://en.wikipedia.org/wiki/Wilkinson%27s_polynomial) para un ejemplo.

### Alternativas

Revisaremos en la nota {ref}`Algoritmos y aplicaciones de eigenvalores, eigenvectores de una matriz <AAEVALEVEC>` algunos algoritmos como:

* Método de la potencia y método de la potencia inversa o iteración inversa.

* Iteración por el cociente de Rayleigh.

* Algoritmo QR.

* Método de rotaciones de Jacobi.

---

## Ejemplos de matrices normales

```{sidebar} Descomposición espectral para matrices normales

Las matrices normales generalizan al caso de entradas en $\mathbb{C}$ la diagonalización ortogonal al ser **unitariamente diagonalizables**. $A \in \mathbb{C}^{n \times n}$ es normal si y sólo si $A = U \Lambda U^H$ con $U$ matriz unitaria (generalización de una matriz ortogonal a entradas $\mathbb{C}$), $U^H$ la conjugada transpuesta de $U$ y $\Lambda$ matriz diagonal. Para $A \in \mathbb{R}^{n \times n}$ lo anterior se escribe como: $A$ es simétrica si y sólo si es ortogonalmente diagonalizable: $A = Q \Lambda Q^T$ (ver {ref}`descomposición espectral <DESCESP>`).

```

$$\begin{array}{l}
\left[
\begin{array}{cc}
1 &-2 \\
2 &1
\end{array}
\right],
\left[
\begin{array}{ccc}
1 &2 & 0\\
0 & 1 & 2\\
2 & 0 & 1
\end{array}
\right]
\end{array}
$$

Otro ejemplo: 

$$A = 
\left[
\begin{array}{ccc}
1 &1 & 0\\
0 & 1 & 1\\
1 & 0 & 1
\end{array}
\right]
$$


```python
A = np.array([[1, 1, 0],
              [0, 1, 1],
              [1, 0, 1.0]])
```


```python
print(A.T@A)
```




    array([[2., 1., 1.],
           [1., 2., 1.],
           [1., 1., 2.]])



```{margin}

Como $A$ es normal entonces se cumple que $AA^T=A^TA$.

```


```python
print(A@A.T)
```




    array([[2., 1., 1.],
           [1., 2., 1.],
           [1., 1., 2.]])




```python
evalue, evector = np.linalg.eig(A)
```


```python
print('eigenvalores:')
print(evalue)
```

    eigenvalores:
    [0.5+0.866j 0.5-0.866j 2. +0.j   ]


```{margin}

Se verifica que los eigenvectores de este ejemplo forman un conjunto linealmente independiente pues $A$ es normal.

```


```python
print('eigenvectores:')
print(evector)
```

    eigenvectores:
    [[-0.289+0.5j -0.289-0.5j -0.577+0.j ]
     [-0.289-0.5j -0.289+0.5j -0.577+0.j ]
     [ 0.577+0.j   0.577-0.j  -0.577+0.j ]]


```{margin}

Para una matriz normal $A$ se cumple que es unitariamente diagonalizable y $A = Q \Lambda Q^H$ donde: $Q^H$ es la conjugada transpuesta de $Q$.

```


```python
print('descomposición espectral:')
Lambda = np.diag(evalue)
Q = evector

```

    descomposición espectral:



```python
print('QLambdaQ^H:')
print(Q@Lambda@Q.conjugate().T)

```

    QLambdaQ^H:
    [[ 1.+0.j  1.+0.j -0.+0.j]
     [ 0.+0.j  1.+0.j  1.+0.j]
     [ 1.+0.j -0.+0.j  1.+0.j]]



```python
print(A)
```




    array([[1., 1., 0.],
           [0., 1., 1.],
           [1., 0., 1.]])



```{margin}

Observa que $Q^HQ=QQ^H = I_3$ donde: $Q^H$ es la conjugada transpuesta de $Q$.
```


```python
print(Q.conjugate().T@Q)
```

    [[1.+0.j 0.-0.j 0.+0.j]
     [0.+0.j 1.+0.j 0.-0.j]
     [0.-0.j 0.+0.j 1.+0.j]]


```{admonition} Observación
:class: tip

El problema del cálculo de eigenvalores para matrices normales es bien condicionado.

```

**Preguntas de comprehensión:**

1)¿Qué son los eigenvalores de una matriz y qué nombre recibe el conjunto de eigenvalores de una matriz?

2)¿Cuántos eigenvalores como máximo puede tener una matriz?

3)¿Qué característica geométrica tiene multiplicar una matriz por su eigenvector?

4)¿A qué se le nombra matriz diagonalizable o *non defective*?

5)¿Cuál es el número de condición del problema de cálculo de eigenvalores con multiplicidad simple para una matriz simétrica?

6)¿Verdadero o Falso? 

a.Una matriz es diagonalizable entonces tiene eigenvalores distintos.

b.Una matriz con eigenvalores distintos es diagonalizable.

c.Si $A=XDX^{-1}$ con $X$ matriz invertible entonces en la diagonal de $D$ y en las columnas de $X$ encontramos eigenvalores y eigenvectores derechos de $A$ respectivamente.

7)Describe la descomposición espectral de una matriz simétrica.

8)¿Qué característica tienen las matrices similares?

**Referencias:**

1. M. T. Heath, Scientific Computing. An Introductory Survey, McGraw-Hill, 2002.

2.  G. H. Golub, C. F. Van Loan, Matrix Computations, John Hopkins University Press, 2013.

3. L. Trefethen, D. Bau, Numerical linear algebra, SIAM, 1997.

4. C. Meyer, Matrix Analysis and Applied Linear Algebra, SIAM, 2000.

47b283e28bec7f1e4f382ff308cabd38f9415830
ipynb
Jupyter Notebook
libro_optimizacion/temas/II.computo_matricial/2.2/Eigenvalores_y_eigenvectores.ipynb
ltejadal/analisis-numerico-computo-cientifico
9695a26779d69acf4d8899e39667567335fcb342
Apache-2.0
libro_optimizacion/temas/II.computo_matricial/2.2/Eigenvalores_y_eigenvectores.ipynb
ltejadal/analisis-numerico-computo-cientifico
9695a26779d69acf4d8899e39667567335fcb342
Apache-2.0
libro_optimizacion/temas/II.computo_matricial/2.2/Eigenvalores_y_eigenvectores.ipynb
ltejadal/analisis-numerico-computo-cientifico
9695a26779d69acf4d8899e39667567335fcb342
Apache-2.0
Qwen/Qwen-72B
1. YES
2. YES