Recordemos la definición de valor y vector propio de una matriz cuadrada \(\mathbf{A}\):

Definición de valor y vector propio de una matriz. Dada una matriz cuadrada \(\mathbf{A}\), \(n\times n\), diremos que \(\mathbf{v}\neq 0\) es un vector propio de valor propio \(\lambda\) de la matriz \(\mathbf{A}\) si \[\mathbf{A}\cdot \mathbf{v}=\lambda\mathbf{v}.\]

Es decir, el “efecto” que tiene la matriz \(\mathbf{A}\) sobre el vector \(\mathbf{v}\) es alargándolo o reduciéndolo un factor \(\lambda\).

El cálculo de valores y vectores propios tiene multitud de aplicaciones en machine learning, entre otras,

• El cálculo de los vectores y valores propios es clave en la diagonalización de una matriz \(\mathbf{A}\), es decir, los vectores propios son representativos de dicha matriz y realizan una tarea similar a los codificadores automáticos en las redes neuronales profundas.

• El análisis de componentes principales (ACP) es una herramienta para encontrar patrones en datos multidimensionales como puede ser en análisis de imágenes. Los profesionales del aprendizaje automático a veces usan el ACP para preprocesar datos para sus redes neuronales. Al centrar, rotar y escalar los datos, el ACP prioriza la dimensionalidad y puede mejorar la velocidad de convergencia de la red neuronal y la calidad general de los resultados.

Recordemos que para hallar los valores propios de una matriz \(\mathbf{A}\), hemos de hallar los ceros de la denominada ecuación característica de \(\mathbf{A}\): \[ p_A(\lambda)=\mathrm{det}(\mathbf{A}-\lambda\cdot\mathbf{I}_n)=0. \] La función \(p_A(\lambda)\) es un polinomio de grado \(n\) en \(\lambda\). Por tanto, hallar valores propios es equivalente a hallar ceros del polinomio característico \(p_A(\lambda)\).

Una vez hallado un valor propio \(\lambda_1\) de la matriz \(A\), para hallar los vectores propios \(\mathbf{v}_1\) de valor propio \(\lambda_1\), hemos de resolver la ecuación: \[ \mathbf{A}\cdot\mathbf{v}_1 =\lambda_1\mathbf{v}_1,\ \Rightarrow (\mathbf{A}-\lambda_1\mathbf{I}_n)\cdot\mathbf{v}_1=\mathbf{0}. \] La ecuación anterior es lineal, homogénea e indeterminada, es decir, tiene muchas soluciones.

Repasar el capítulo de Álgebra lineal numérica para ver algún ejemplo de cálculo así como las propiedades de los valores y vectores propios de una matriz cuadrada \(\mathbf{A}\).

El método de la potencia es un método iterativo que calcula el valor propio de módulo máximo con un vector propio asociado de la matriz cuadrada.

Usando variaciones del método, se pueden hallar los demás valores propios junto con los vectores propios asociados.

Para poder aplicar el método de la potencia, necesitamos suponer que la matriz \(\mathbf{A}\), de dimensiones \(n\times n\) es diagonalizable, es decir, que \(\mathbf{A}\) tiene \(n\) valores propios \(\lambda_1,\lambda_2,\ldots,\lambda_n\) que supondremos ordenados donde \(\lambda_1\) tiene el módulo máximo: \[ |\lambda_1|>|\lambda_2|\geq |\lambda_3|\geq\cdots\geq |\lambda_n|. \] Fijémonos que los valores propios \(\lambda_i\), \(i=2,\ldots,n\) pueden ser repetidos.

Como la matriz \(\mathbf{A}\) es diagonalizable, sean \(\mathbf{v}^{(1)},\mathbf{v}^{(2)},\ldots,\mathbf{v}^{(n)}\) unos vectores propios asociados a cada valor propio \(\lambda_i\), es decir, \(\mathbf{A}\mathbf{v}^{(i)}=\lambda_i\mathbf{v}^{(i)}\), \(i=1,\ldots,n\).

Veamos cómo funciona el método.

Sea \(\mathbf{x}\) un vector cualquiera. Como los vectores propios \(\mathbf{v}^{(1)},\mathbf{v}^{(2)},\ldots,\mathbf{v}^{(n)}\) son linealmente independientes, forman una base del espacio vectorial \(\mathbb{R}^n\), y por tanto, existen unos valores \(\alpha_1,\ldots,\alpha_n\) tal que: \[ \mathbf{x}=\alpha_1\mathbf{v}^{(1)}+\cdots +\alpha_n\mathbf{v}^{(n)}=\sum_{i=1}^n\alpha_i\mathbf{v}^{(i)}. \]

Si calculamos los vectores \(\mathbf{A}\mathbf{x},\mathbf{A}^2\mathbf{x},\ldots,\mathbf{A}^k\mathbf{x},\ldots,\) obtenemos: \[ \begin{align*} \mathbf{A}\mathbf{x} & = \mathbf{A}(\alpha_1\mathbf{v}^{(1)}+\cdots +\alpha_n\mathbf{v}^{(n)})=\alpha_1\mathbf{A}\mathbf{v}^{(1)}+\cdots +\alpha_n\mathbf{A}\mathbf{v}^{(n)}\\ &=\alpha_1\lambda_1\mathbf{v}^{(1)}+\cdots +\alpha_n\lambda_n\mathbf{v}^{(n)}=\sum_{i=1}^n\alpha_i\lambda_i\mathbf{v}^{(i)},\\ & \vdots\\ \mathbf{A}^k\mathbf{x} & = \mathbf{A}^k(\alpha_1\mathbf{v}^{(1)}+\cdots +\alpha_n\mathbf{v}^{(n)})=\alpha_1\mathbf{A}^k\mathbf{v}^{(1)}+\cdots +\alpha_n\mathbf{A}^k\mathbf{v}^{(n)}\\ &=\alpha_1\lambda_1^k\mathbf{v}^{(1)}+\cdots +\alpha_n\lambda_n^k\mathbf{v}^{(n)}=\sum_{i=1}^n\alpha_i\lambda_i^k\mathbf{v}^{(i)},\\ & \vdots \end{align*} \]

Tenemos por tanto que, en general para un entero \(k\) positivo, \(\displaystyle \mathbf{A}^k\mathbf{x}=\sum_{i=1}^n\alpha_i\lambda_i^k\mathbf{v}^{(i)}\). La expresión se puede escribir de la forma siguiente: \[ \mathbf{A}^k\mathbf{x}=\sum_{i=1}^n\alpha_i\lambda_i^k\mathbf{v}^{(i)}=\lambda_1^k\sum_{i=1}^n\alpha_i\left(\frac{\lambda_i}{\lambda_1}\right)^k\mathbf{v}^{(i)}. \] Como \(|\lambda_1|>|\lambda_i|\), para \(i=2,\ldots,n\), tenemos que \(\left|\frac{\lambda_i}{\lambda_1}\right|<1\) y, por tanto, \(\displaystyle\lim_{k\to\infty}\left(\frac{\lambda_i}{\lambda_1}\right)^k=0\).

Usando la expresión anterior, podemos afirmar que si \(\alpha_1\neq 0\), \[ \lim_{k\to\infty}\mathbf{A}^k\mathbf{x}=\lim_{k\to\infty}\lambda_1^k\alpha_1\mathbf{v}^{(1)}=\begin{cases} \mathbf{0}, & \mbox{si }|\lambda_1|<1,\\ \infty, & \mbox{si }|\lambda_1|>1,\\ \alpha_1\mathbf{v}^{(1)}, & \mbox{si }\lambda_1 =1,\\ \mbox{diverge}, & \mbox{si }|\lambda_1|=1\mbox{ y } \lambda_1\neq 1.\end{cases} \]

Una vez vistos los preliminares, expliquemos cómo podemos hallar el valor propio \(\lambda_1\) junto con un vector propio \(\mathbf{v}^{(1)}\) hallando una sucesión de valores reales \((\mu^{(n)})_{n=0,1,\ldots}\) y una sucesión de vectores \((\mathbf{x}^{(n)})_{n=0,1,\ldots}\) que convergen a \(\lambda_1\) y \(\mathbf{v}^{(1)}\) respectivamente.

• Paso 1. Supongamos que elegimos un vector \(\mathbf{x}^{(0)}\neq \mathbf{0}\) tal que el coeficiente \(\alpha_1\) en la descomposición de \(\mathbf{x}^{(0)}\) usando los vectores propios no es cero: \[ \mathbf{x}^{(0)}=\alpha_1\mathbf{v}^{(1)}+\cdots +\alpha_n\mathbf{v}^{(n)}=\sum_{i=1}^n\alpha_i\mathbf{v}^{(i)},\quad \alpha_1\neq 0. \] Supongamos también que \(\|\mathbf{x}^{(0)}\|_\infty =1\). Si no fuera así, consideramos: \(\frac{\mathbf{x}^{(0)}}{\|\mathbf{x}^{(0)}\|_\infty}\) como vector inicial \(\mathbf{x}^{(0)}\) y ya se verificará. Recordemos que \(\|\mathbf{x}^{(0)}\|_\infty=\max_{i=1,\ldots,n}|x_i^{(0)}|=1\). Por tanto, existirá una componente del vector \(\mathbf{x}^{(0)}\) a la que llamaremos \(m_0\) tal que \(|x^{(0)}_{m_0}|=1\). Definimos \(\displaystyle\mathbf{y}^{(1)}=\mathbf{A}\mathbf{x}^{(0)}=\sum_{i=1}^n\alpha_i\lambda_i\mathbf{v}^{(i)}\). Definimos \(\mu^{(1)}=\frac{y^{(1)}_{m_0}}{x^{(0)}_{m_0}}\), es decir, el cociente entre las componentes \(m_0\) de los vectores \(\mathbf{y}^{(1)}\) y \(\mathbf{x}^{(0)}\).

• Paso 1. (continuación) Podemos escribir \(\mu^{(1)}\) como: \[ \mu^{(1)}=\frac{y^{(1)}_{m_0}}{x^{(0)}_{m_0}}=\frac{\alpha_1\lambda_1{v}^{(1)}_{m_0}+\sum_{i=2}^n\alpha_i\lambda_i{v}^{(i)}_{m_0}}{\alpha_1\mathbf{v}^{(1)}_{m_0}+\sum_{i=2}^n\alpha_i {v}^{(i)}_{m_0}}. \] Definimos \(\mathbf{x}^{(1)}=\frac{\mathbf{y}^{(1)}}{\|\mathbf{y}^{(1)}\|_\infty}\).

• Paso 2. Como \(\|\mathbf{x}^{(1)}\|_\infty =1\), existe una componente de dicho vector \(m_1\) tal que \(|x^{(1)}_{m_1}|=1\). Definimos \(\displaystyle\mathbf{y}^{(2)}=\mathbf{A}\mathbf{x}^{(1)}=\frac{1}{\|\mathbf{y}^{(1)}\|_\infty}\mathbf{A}\mathbf{y}^{(1)}=\frac{1}{\|\mathbf{y}^{(1)}\|_\infty}\mathbf{A}^2\mathbf{x}^{(0)}=\frac{1}{\|\mathbf{y}^{(1)}\|_\infty}\sum_{i=1}^n\alpha_i\lambda_i^2\mathbf{v}^{(i)}\). Definimos \(\mu^{(2)}=\frac{y^{(2)}_{m_1}}{x^{(1)}_{m_1}}\), es decir, el cociente entre las componentes \(m_1\) de los vectores \(\mathbf{y}^{(2)}\) y \(\mathbf{x}^{(1)}\).

• Paso 2. (continuación) Podemos escribir \(\mu^{(2)}\) como: \[ \begin{align*} \mu^{(2)} & =\frac{y^{(2)}_{m_1}}{x^{(1)}_{m_1}}=\frac{(\alpha_1\lambda_1^2{v}^{(1)}_{m_1}+\sum_{i=2}^n\alpha_i\lambda_i^2 {v}^{(i)}_{m_1})/\|\mathbf{y}^{(1)}\|_\infty}{(\alpha_1\lambda_1 {v}^{(1)}_{m_1}+\sum_{i=2}^n\alpha_i\lambda_i {v}^{(i)}_{m_1})/\|\mathbf{y}^{(1)}\|_\infty}\\ & =\lambda_1\left( \frac{\alpha_1 {v}^{(1)}_{m_1}+\sum_{i=2}^n\alpha_i(\lambda_i/\lambda_1)^2 {v}^{(i)}_{m_1}}{\alpha_1 {v}^{(1)}_{m_1}+\sum_{i=2}^n\alpha_i(\lambda_i/\lambda_1) {v}^{(i)}_{m_1}}\right) \end{align*} \] Definimos \(\mathbf{x}^{(2)}=\frac{\mathbf{y}^{(2)}}{\|\mathbf{y}^{(2)}\|_\infty}\).
• \(\vdots\)

• Paso \(k\). En este paso, tenemos definidos
  • dos secuencias de vectores \(\mathbf{x}^{(j)}\), para \(j=0,\ldots,k-1\) y \(\mathbf{y}^{(j)}\), \(j=1,\ldots,k-1\) cumpliendo:
    • \(\mathbf{y}^{(j)}=\mathbf{A}\mathbf{x}^{(j-1)}\),
    • \(\mathbf{x}^{(j)}=\frac{1}{\prod_{l=1}^j \|\mathbf{y}^{(l)}\|_\infty}\mathbf{A}^j\mathbf{x}^{(0)}\) y \(\|\mathbf{x}^{(j)}\|_\infty=1\).
  • una secuencia de valores reales \(\mu^{(j)}\), \(j=1,\ldots,k-1\) y una secuencia de valores enteros \(m_j\), \(j=0,\ldots,k-2\) tal que: \[ \mu^{(j)}=\frac{y^{(j)}_{m_{j-1}}}{x^{(j-1)}_{m_{j-1}}}=\lambda_1\left( \frac{\alpha_1 {v}^{(1)}_{m_{j-1}}+\sum_{i=2}^n\alpha_i(\lambda_i/\lambda_1)^j {v}^{(i)}_{m_{j-1}}}{\alpha_1 {v}^{(1)}_{m_{j-1}}+\sum_{i=2}^n\alpha_i(\lambda_i/\lambda_1)^{j-1} {v}^{(i)}_{m_{j-1}}}\right) \]

• Paso \(k\) (continuación) Como \(\|\mathbf{x}^{(k-1)}\|=1\), existe una componente de dicho vector \(m_{k-1}\) tal que \(|x^{(k-1)}_{m_{k-1}}|=1\). Definimos \[ \mathbf{y}^{(k)}=\mathbf{A}\mathbf{x}^{(k-1)}=\frac{1}{\prod_{l=1}^{k-1}\|\mathbf{y}^{(l)}\|_\infty}\mathbf{A}^k\mathbf{x}^{(0)}=\frac{1}{\prod_{l=1}^{k-1}\|\mathbf{y}^{(l)}\|_\infty}\sum_{i=1}^n\alpha_i\lambda_i^k\mathbf{v}^{(i)}. \] Definimos \(\mu^{(k)}=\frac{y^{(k)}_{m_{k-1}}}{x^{(k-1)}_{m_{k-1}}}\), es decir, el cociente entre las componentes \(m_{k-1}\) de los vectores \(\mathbf{y}^{(k)}\) y \(\mathbf{x}^{(k-1)}\).

• Paso \(k\) (continuación). Podemos escribir \(\mu^{(k)}\) como: \[ \begin{align*} \mu^{(k)} & =\frac{y^{(k)}_{m_{k-1}}}{x^{(k-1)}_{m_{k-1}}}=\frac{(\alpha_1\lambda_1^k{v}^{(1)}_{m_{k-1}}+\sum_{i=2}^n\alpha_i\lambda_i^k {v}^{(i)}_{m_{k-1}})/\prod_{l=1}^{k-1}\|\mathbf{y}^{(l)}\|}{(\alpha_1\lambda_1^{k-1} {v}^{(1)}_{m_{k-1}}+\sum_{i=2}^n\alpha_i\lambda_i^{k-1} {v}^{(i)}_{m_{k-1}})/\prod_{l=1}^{k-1}\|\mathbf{y}^{(l)}\|}\\ & =\lambda_1\left( \frac{\alpha_1 {v}^{(1)}_{m_{k-1}}+\sum_{i=2}^n\alpha_i(\lambda_i/\lambda_1)^k {v}^{(i)}_{m_{k-1}}}{\alpha_1 {v}^{(1)}_{m_{k-1}}+\sum_{i=2}^n\alpha_i(\lambda_i/\lambda_1)^{k-1} {v}^{(i)}_{m_{k-1}}}\right). \end{align*} \] Por último definimos: \(\mathbf{x}^{(k)}=\frac{\mathbf{y}^{(k)}}{\|\mathbf{y}^{(k)}\|_\infty}\).
• \(\vdots\)

Usando la ecuación de la sucesión \(\mu^{(k)}\), podemos afirmar que la sucesión de números reales \((\mu^{(k)})_{k\geq 1}\) tiende hacia el valor propio \(\lambda_1\): \[ \lim_{k\to\infty}\mu^{(k)}=\lim_{k\to\infty}\lambda_1\left( \frac{\alpha_1 {v}^{(1)}_{m_{k-1}}+\sum_{i=2}^n\alpha_i(\lambda_i/\lambda_1)^k {v}^{(i)}_{m_{k-1}}}{\alpha_1 {v}^{(1)}_{m_{k-1}}+\sum_{i=2}^n\alpha_i(\lambda_i/\lambda_1)^{k-1} {v}^{(i)}_{m_{k-1}}}\right) =\lambda_1, \] ya que \(\left|\frac{\lambda_i}{\lambda_1}\right|<1\) y por tanto \(\displaystyle\lim_{k\to\infty}\left(\frac{\lambda_i}{\lambda_1}\right)^k=0\), para \(k=2,\ldots,n\).

Además como \(\mathbf{x}^{(k)}=\frac{\mathbf{y}^{(k)}}{\|\mathbf{y}^{(k)}\|}\) y \(\mathbf{x}^{(k)}=\frac{1}{\prod_{l=1}^k\|\mathbf{y}^{(l)}\|}\mathbf{A}^k\mathbf{x}^{(0)}\), tenemos que: \[ \begin{align*} \mathbf{y}^{(k)} & =\mathbf{A}\mathbf{x}^{(k-1)}=\frac{1}{\prod_{l=1}^{k-1}\|\mathbf{y}^{(l)}\|}\mathbf{A}^k\mathbf{x}^{(0)},\ \Rightarrow \\ \mathbf{x}^{(k)} & =\frac{\mathbf{y}^{(k)}}{\|\mathbf{y}^{(k)}\|}=\frac{\mathbf{A}^k\mathbf{x}^{(0)}/\prod_{l=1}^{k-1}\|\mathbf{y}^{(l)}\|}{\left\|\mathbf{A}^k\mathbf{x}^{(0)}/\prod_{l=1}^{k-1}\|\mathbf{y}^{(l)}\|\right\|}=\frac{\mathbf{A}^k\mathbf{x}^{(0)}}{\left\|\mathbf{A}^k\mathbf{x}^{(0)}\right\|}\\ & = \frac{\alpha_1\lambda_1^k\mathbf{v}^{(1)}+\sum_{i=2}^n \alpha_i\lambda_i^k\mathbf{v}^{(i)}}{\left\|\alpha_1\lambda_1^k\mathbf{v}^{(1)}+\sum_{i=2}^n \alpha_i\lambda_i^k\mathbf{v}^{(i)}\right\|}\\ & =\left(\frac{\lambda_1}{|\lambda_1|}\right)^k \frac{\alpha_1 \mathbf{v}^{(1)}+\sum_{i=2}^n \alpha_i(\lambda_i/\lambda_1)^k\mathbf{v}^{(i)}}{\left\|\alpha_1 \mathbf{v}^{(1)}+\sum_{i=2}^n \alpha_i(\lambda_i/\lambda_1)^k\mathbf{v}^{(i)}\right\|}. \end{align*} \]

Entonces, \[ \lim_{k\to\infty}\mathbf{x}^{(k)}=\pm\mathbf{v}^{(1)}, \] es decir la sucesión de vectores \(\mathbf{x}^{(k)}\) tiende hacia un vector propio de valor propio \(\lambda_1\).

En resumen, hemos obtenido una sucesión de números reales \((\mu^{(k)})_{k\geq 1}\) que tiende al valor propio de módulo máximo \(\lambda_1\) y una sucesión de vectores \((\mathbf{x}^{(k)})_{k\geq 0}\) de módulo \(1\) que tiende a un vector propio de valor propio \(\lambda_1\).

Observación 1.

Podemos escribir que para \(k\) grande: \[ \mu^{(k)}\approx \lambda_1+O\left(\left|\frac{\lambda_i}{\lambda_1}\right|\right)^k, \] ya que: \[ \begin{align*} \mu^{(k)} & \approx \lambda_1\left( \frac{\alpha_1 {v}^{(1)}_{m_{k-1}}+\sum_{i=2}^n\alpha_i(\lambda_i/\lambda_1)^k {v}^{(i)}_{m_{k-1}}}{\alpha_1 {v}^{(1)}_{m_{k-1}}+\sum_{i=2}^n\alpha_i(\lambda_i/\lambda_1)^{k-1} {v}^{(i)}_{m_{k-1}}}\right)\\ & \approx\lambda_1\left( \frac{\alpha_1 {v}^{(1)}_{m_{k-1}}+\sum_{i=2}^n\alpha_i(\lambda_i/\lambda_1)^k {v}^{(i)}_{m_{k-1}}}{\alpha_1 {v}^{(1)}_{m_{k-1}}}\right)\\ & =\lambda_1 \left(1+\frac{\sum_{i=2}^n\alpha_i(\lambda_i/\lambda_1)^k {v}^{(i)}_{m_{k-1}}}{\alpha_1 {v}^{(1)}_{m_{k-1}}}\right)\approx \lambda_1+O\left(\left|\frac{\lambda_i}{\lambda_1}\right|\right)^k. \end{align*} \]

Observación 1 (continuación).

De la misma manera podemos escribir que para \(k\) grande: \[ \mathbf{x}^{(k)}\approx \pm \frac{\mathbf{v}^{(1)}}{\|\mathbf{v}^{(1)}\|}+O\left(\left|\frac{\lambda_i}{\lambda_1}\right|\right)^k. \] ya que: \[ \begin{align*} \mathbf{x}^{(k)} &\approx\pm \frac{\alpha_1 \mathbf{v}^{(1)}+\sum_{i=2}^n \alpha_i(\lambda_i/\lambda_1)^k\mathbf{v}^{(i)}}{\left\|\alpha_1 \mathbf{v}^{(1)}+\sum_{i=2}^n \alpha_i(\lambda_i/\lambda_1)^k\mathbf{v}^{(i)}\right\|}\\ & \approx \pm \frac{\alpha_1 \mathbf{v}^{(1)}+\sum_{i=2}^n \alpha_i(\lambda_i/\lambda_1)^k\mathbf{v}^{(i)}}{\left\|\alpha_1 \mathbf{v}^{(1)}\right\|}=\pm \frac{\mathbf{v}^{(1)}}{\|\mathbf{v}^{(1)}\|}+O\left(\left|\frac{\lambda_i}{\lambda_1}\right|\right)^k. \end{align*} \]

Observación 1 (continuación).

El resumen, el orden de convergencia de las sucesiones \(\mu^{(k)}\) y \(\mathbf{x}^{(k)}\) vale \(O\left(\frac{\max_{i=2,\ldots,n}|\lambda_i|}{|\lambda_1|}\right)^k\).

Observación 2.

Aunque la matriz \(\mathbf{A}\) tenga un valor propio dominante \(\lambda_1\) con multiplicidad \(r>1\), el método de la potencia seguirá convergiendo a \(\lambda_1\) y los vectores \(\mathbf{x}^{(k)}\) seguirán convergiendo a un vector propio de valor propio \(\lambda_1\).

Observación 3.

Puede haber matrices \(\mathbf{A}\) en las que cuando hallamos la sucesión \(\mathbf{x}^{(k)}\) usando el método de la potencia, puede haber alternancias de signo, es decir, que para valores de \(k\) relativamente grandes, las aproximaciones del vector propio \(\mathbf{v}^{(1)}\) de valor propio dominante \(\lambda_1\) aparecen de la forma siguiente: \[ \ldots, \mathbf{x}^{(k)}\approx\frac{\mathbf{v}^{(1)}}{\|\mathbf{v}^{(1)}\|_\infty},\mathbf{x}^{(k+1)}\approx -\frac{\mathbf{v}^{(1)}}{\|\mathbf{v}^{(1)}\|_\infty},\mathbf{x}^{(k+2)}\approx \frac{\mathbf{v}^{(1)}}{\|\mathbf{v}^{(1)}\|_\infty},\mathbf{x}^{(k+3)}\approx -\frac{\mathbf{v}^{(1)}}{\|\mathbf{v}^{(1)}\|_\infty}\ldots \] Si nos encontramos en situaciones como la descrita, no podemos usar como criterio de parada la condición \(\|\mathbf{x}^{(k)}-\mathbf{x}^{(k-1)}\|_\infty <\epsilon\), donde \(\epsilon\) es la tolerancia permitida ya que \(\|\mathbf{x}^{(k)}-\mathbf{x}^{(k-1)}\|_\infty\approx 2\left\|\frac{\mathbf{v}^{(1)}}{\|\mathbf{v}^{(1)}\|_\infty}\right\|_\infty =2\) y el algoritmo nunca pararía.

Observación 3. (continuación)

Una manera de arreglar esta situación es cambiar el criterio de parada de la manera siguiente: \[ \min\{\|\mathbf{x}^{(k)}-\mathbf{x}^{(k-1)}\|_\infty,\|\mathbf{x}^{(k)}+\mathbf{x}^{(k-1)}\|_\infty\}<\epsilon. \]

• INPUT matriz del sistema \(\mathbf{A}=(a_{ij})_{i=1,\ldots,n,j=1,\ldots,n}\), valor inicial \(\mathbf{X0}=\mathbf{x}\) con \(\|\mathbf{x}\|_\infty =1\), error absoluto o tolerancia TOL, número máximo de iteraciones Nmax.
• Set \(k=1\) (inicializamos el contador de la iteraciones)
• Hallamos el entero más pequeño \(m\) tal que \(|x_m|=\|x\|_\infty\).
• Set \(\mathbf{x}=\frac{\mathbf{x}}{x_m}\).

• While \(k\leq Nmax\) (mientras no se llegue al número máximo de iteraciones, hacer lo siguiente)
  • Set \(\mathbf{y}=\mathbf{A}\mathbf{x}\)
  • Set \(\mu =y_m\).
  • Hallamos el entero más pequeño \(m\) tal que \(|y_m|=\|y\|_\infty\).
  • If \(y_m=0\) then print(la matriz A tiene 0 como valor propio. Escoger otro vector inicial).
  • Set \(\mathbf{x}_1=\frac{\mathbf{y}}{y_m}\)
  • Set error=\(\min\left\{\left\|\mathbf{x}-\mathbf{x}_1\right\|_\infty,\left\|\mathbf{x}+\mathbf{x}_1\right\|_\infty\right\}\).
  • If error<TOL print(\(\mu,\mathbf{x}\)). (El método ha acabado. Tenemos el valor propio y un vector propio con un error menor que la tolerancia permitida.)

• • Set \(\mathbf{x}=\mathbf{x}_1\).
  • Set \(k=k+1\).
• Print Hemos alcanzado el máximo número de iteraciones. El método de la potencia no converge en este caso.

Si aplicamos el método de Aitken (ver el capítulo de ceros) a la sucesión \((\mu^{(n)})_{n\geq 1}\) obtenemos una sucesión que converge más rápidamente al valor propio de módulo máximo \(\lambda_1\).

Concretamente, consideramos la sucesión \((\tilde{\mu}^{(n)})_{n\geq 1}\) definida de la forma siguiente: \[ \tilde{\mu}^{(n)}=\mu^{(n)}-\frac{(\mu^{(n+1)}-\mu^{(n)})^2}{\mu^{(n+2)}-2\mu^{(n+1)}+\mu^{(n)}}. \]

Para aplicar la aceleración anterior, modificamos el pseudódigo del método de la potencia de la manera siguiente:

• INPUT matriz del sistema \(\mathbf{A}=(a_{ij})_{i=1,\ldots,n,j=1,\ldots,n}\), valor inicial \(\mathbf{X0}=\mathbf{x}\) con \(\|\mathbf{x}\|_\infty =1\), error absoluto o tolerancia TOL, número máximo de iteraciones Nmax.
• Set \(k=1\) (inicializamos el contador de la iteraciones)
• Set \(\mu_0=0\).
• Set \(\mu_1=0\). (inicializamos los dos primeros valores de la sucesión \((\mu^{(n)})_{n\geq 1}\) para poder aplicar la aceleración de Aitken.
• Hallamos el entero más pequeño \(m\) tal que \(|x_m|=\|\mathbf{x}\|_\infty\).

• While \(k\leq Nmax\) (mientras no se llegue al número máximo de iteraciones, hacer lo siguiente)
  • Set \(\mathbf{y}=\mathbf{A}\mathbf{x}\)
  • Set \(\mu =y_m\).
  • Hallamos el entero más pequeño \(m\) tal que \(|y_m|=\|\mathbf{y}\|_\infty\).
  • If \(y_m=0\) then print(la matriz A tiene 0 como valor propio. Escoger otro vector inicial).
  • Set \(\hat{\mu}=\mu_0-\frac{(\mu_1-\mu_0)^2}{\mu-2\mu_1+\mu_0}\). (Aplicamos aceleración de Aitken)
  • Set \(\mathbf{x}_1=\frac{\mathbf{y}}{y_m}\).
  • Set error=\(\min\left\{\left\|\mathbf{x}-\mathbf{x}_1\right\|_\infty,\left\|\mathbf{x}+\mathbf{x}_1\right\|_\infty\right\}\).

• • If error<TOL print(\(\mu,\mathbf{x}\)). (El método ha acabado. Tenemos el valor propio y un vector propio con un error menor que la tolerancia permitida.)
  • Set \(\mathbf{x}=\mathbf{x}_1\).
  • Set \(k=k+1\).
  • Set \(\mu_0=\mu_1\).
  • Set \(\mu_1=\mu\). (Actualizamos los valores \(\mu^{(n)}\))
• Print Hemos alcanzado el máximo número de iteraciones. El método de la potencia no converge en este caso.

Recordemos que si la matriz \(\mathbf{A}\) es simétrica existe una base de vectores propios ortogonal, es decir, existen vectores propios \(\mathbf{v}^{(1)},\ldots,\mathbf{v}^{(n)}\) de valores propios \(\lambda_1,\ldots,\lambda_n\), respectivamente, con \((\mathbf{v}^{(i)})^\top \mathbf{v}^{(j)}=0\) si \(i\neq j\).

Igual que en las secciones anteriores, supondremos que \(\lambda_1\) es el valor propio dominante: \(|\lambda_1|>|\lambda_2|\geq \cdots \geq |\lambda_n|\).

En este caso, vamos a modificar el método de la potencia para conseguir una mayor velocidad de convergencia de la sucesión \(\mu^{(n)}\) que converge hacia \(\lambda_1\).

Concretamente, vamos a definir dos sucesiones \(\mu^{(k)}\) y \(\mathbf{x}^{(k)}\) que convergen hacia \(\lambda_1\) y un vector propio de valor propio \(\lambda_1\), la primera en orden \(O\left(\left|\frac{\lambda_i}{\lambda_1}\right|\right)^{2k}\) y la segunda en orden \(O\left(\left|\frac{\lambda_i}{\lambda_1}\right|\right)^{k}\), \(i=2,\ldots,n\).

Otra modificación a tener en cuenta es que usaremos la norma euclídea \(\|\cdot\|_2\) en lugar de la norma infinito \(\|\cdot\|_\infty\).

El algoritmo del método de la potencia para matrices simétricas es el siguiente:

• Paso \(1\):
  • Consideramos \(\mathbf{x}^{(0)}\) un vector inicial con \(\|\mathbf{x}^{(0)}\|_2=1\). Si \(\mathbf{x}^{(0)}\) no tuviese normal euclídea igual a 1, consideraríamos \(\frac{\mathbf{x}^{(0)}}{\|\mathbf{x}^{(0)}\|_2}\) y este último ya tendría norma euclídea igual a \(1\).
  • Definimos \(\mathbf{y}^{(1)}=\mathbf{A}\mathbf{x}^{(0)}\).
  • Definimos \(\mu^{(1)}=(\mathbf{x}^{(0)})^\top \mathbf{y}^{(1)}\).
  • Por último, definimos \(\mathbf{x}^{(1)}=\frac{\mathbf{y}^{(1)}}{\|\mathbf{y}^{(1)}\|_2}\).
• \(\vdots\)

• Paso \(k\). Suponemos definidos \(\mathbf{x}^{(0)},\ldots,\mathbf{x}^{(k-1)}\), \(\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(k-1)}\), \(\mu^{(1)},\ldots,\mu^{(k-1)}\).
  • Definimos \(\mathbf{y}^{(k)}=\mathbf{A}\mathbf{x}^{(k-1)}\).
  • Definimos \(\mu^{(k)}=(\mathbf{x}^{(k-1)})^\top \mathbf{y}^{(k)}\).
  • Definimos \(\mathbf{x}^{(k)}=\frac{\mathbf{y}^{(k)}}{\|\mathbf{y}^{(k)}\|_2}\).

Para demostrar que la sucesión \(\mu^{(k)}\) tiene orden de convergencia \(O\left(\left|\frac{\lambda_i}{\lambda_1}\right|\right)^{2k}\), necesitamos usar el resultado siguiente:

Proposición.

Para todo valor de \(k\geq 0\), el vector \(\mathbf{x}^{(k)}\) puede escribirse como: \[\mathbf{x}^{(k)}=\frac{\mathbf{A}^k\mathbf{x}^{(0)}}{\|\mathbf{A}^k\mathbf{x}^{(0)}\|_2}=\frac{\alpha_1\lambda_1^k\mathbf{v}^{(1)}+\sum_{i=2}^n\alpha_i\lambda_i^k\mathbf{v}^{(i)}}{\|\alpha_1\lambda_1^k\mathbf{v}^{(1)}+\sum_{i=2}^n\alpha_i\lambda_i^k\mathbf{v}^{(i)}\|_2}.\]

La sucesión \(\mu^{(k)}\), usando la proposición anterior, puede escribirse como: \[ \begin{align*} \mu^{(k)} & =(\mathbf{x}^{(k-1)})^\top \mathbf{y}^{(k)}=(\mathbf{x}^{(k-1)})^\top \mathbf{A}\mathbf{x}^{(k-1)}\\ & =\frac{(\alpha_1\lambda_1^{k-1}\mathbf{v}^{(1)}+\sum_{i=2}^n\alpha_i\lambda_i^{k-1}\mathbf{v}^{(i)})^\top \mathbf{A}(\alpha_1\lambda_1^{k-1}\mathbf{v}^{(1)}+\sum_{i=2}^n\alpha_i\lambda_i^{k-1}\mathbf{v}^{(i)})}{\|\alpha_1\lambda_1^{k-1}\mathbf{v}^{(1)}+\sum_{i=2}^n\alpha_i\lambda_i^{k-1}\mathbf{v}^{(i)}\|_2^2}\\ & =\frac{(\alpha_1\lambda_1^{k-1}\mathbf{v}^{(1)}+\sum_{i=2}^n\alpha_i\lambda_i^{k-1}\mathbf{v}^{(i)})^\top (\alpha_1\lambda_1^{k}\mathbf{v}^{(1)}+\sum_{i=2}^n\alpha_i\lambda_i^{k}\mathbf{v}^{(i)})}{\|\alpha_1\lambda_1^{k-1}\mathbf{v}^{(1)}+\sum_{i=2}^n\alpha_i\lambda_i^{k-1}\mathbf{v}^{(i)}\|_2^2}. \end{align*} \]

Usando la ortogonalidad de los vectores propios \(\mathbf{v}^{(i)}\), podemos escribir \(\mu^{(k)}\) como: \[ \begin{align*} \mu^{(k)} & =\frac{\alpha_1^2\lambda_1^{2k-1}+\sum_{i=2}^n \alpha_i^2\lambda_i^{2k-1}}{\alpha_1^2\lambda_1^{2k-2}+\sum_{i=2}^n \alpha_i^2\lambda_i^{2k-2}}=\frac{\alpha_1^2\lambda_1^{2k-1}\left(1+\sum_{i=2}^n(\alpha_i/\alpha_1)^2 (\lambda_i/\lambda_1)^{2k-1}\right)}{\alpha_1^2\lambda_1^{2k-2}\left(1+\sum_{i=2}^n(\alpha_i/\alpha_1)^2 (\lambda_i/\lambda_1)^{2k-2}\right)}\\ & = \lambda_1 \frac{\left(1+\sum_{i=2}^n(\alpha_i/\alpha_1)^2 (\lambda_i/\lambda_1)^{2k-1}\right)}{\left(1+\sum_{i=2}^n(\alpha_i/\alpha_1)^2 (\lambda_i/\lambda_1)^{2k-2}\right)}\\ & \approx \lambda_1 \left(1+\sum_{i=2}^n(\alpha_i/\alpha_1)^2 (\lambda_i/\lambda_1)^{2k-1}\right)=\lambda_1+O\left(\left|\frac{\lambda_i}{\lambda_1}\right|\right)^{2k}, \end{align*} \] tal como queríamos ver.

Los vectores \(\mathbf{x}^{(k)}\), usando la proposición anterior, se pueden escribir como: \[ \begin{align*} \mathbf{x}^{(k)} & =\frac{\alpha_1\lambda_1^k\mathbf{v}^{(1)}+\sum_{i=2}^n\alpha_i\lambda_i^k\mathbf{v}^{(i)}}{\|\alpha_1\lambda_1^k\mathbf{v}^{(1)}+\sum_{i=2}^n\alpha_i\lambda_i^k\mathbf{v}^{(i)}\|_2}\\ & =\frac{\alpha_1\lambda_1^k\left(\mathbf{v}^{(1)}+\sum_{i=2}^n(\alpha_i/\alpha_1)(\lambda_i/\lambda_1)^k\mathbf{v}^{(i)}\right)}{\|\alpha_1\lambda_1^k\left(\mathbf{v}^{(1)}+\sum_{i=2}^n(\alpha_i/\alpha_1)(\lambda_i/\lambda_1)^k\mathbf{v}^{(i)}\right)\|_2}\\ & \approx \frac{\alpha_1\lambda_1^k\left(\mathbf{v}^{(1)}+\sum_{i=2}^n(\alpha_i/\alpha_1)(\lambda_i/\lambda_1)^k\mathbf{v}^{(i)}\right)}{\|\alpha_1\lambda_1^k\left(\mathbf{v}^{(1)}\right)\|_2}=\pm \frac{\mathbf{v}^{(1)}}{\|\mathbf{v}^{(1)}\|_2}+O\left(\left|\frac{\lambda_i}{\lambda_1}\right|\right)^{k}, \end{align*} \] es decir, la sucesión \(\mathbf{x}^{(k)}\) converge hacia un vector propio de valor propio \(\lambda_1\) con orden de convergencia \(O\left(\left|\frac{\lambda_i}{\lambda_1}\right|\right)^{k}\).

• INPUT matriz del sistema \(\mathbf{A}=(a_{ij})_{i=1,\ldots,n,j=1,\ldots,n}\), valor inicial \(\mathbf{X0}=\mathbf{x}\neq \mathbf{0}\) con \(\|\mathbf{x}\|_2 =1\), error absoluto o tolerancia TOL, número máximo de iteraciones Nmax.
• Set \(k=1\) (inicializamos el contador de la iteraciones)
• Set \(\mathbf{x}=\frac{\mathbf{x}}{\|\mathbf{x}\|_2}\).

• While \(k\leq Nmax\) (mientras no se llegue al número máximo de iteraciones, hacer lo siguiente)
  • Set \(\mathbf{y}=\mathbf{A}\mathbf{x}\)
  • Set \(\mu =\mathbf{x}^\top \mathbf{y}\).
  • If \(\|\mathbf{y}\|_2=0\) then print(A tiene 0 como valor propio. Escoge otro vector propio inicial).
  • Set \(\mathbf{x}_1=\frac{\mathbf{y}}{\|\mathbf{y}\|_2}\)
  • Set error=\(\min\left\{\left\|\mathbf{x}-\mathbf{x}_1\right\|_2,\left\|\mathbf{x}+\mathbf{x}_1\right\|_2\right\}\).
  • If error<TOL print(\(\mu,\mathbf{x}\)). (El método ha acabado. Tenemos el valor propio y un vector propio con un error menor que la tolerancia permitida.)

• • Set \(\mathbf{x}=\mathbf{x}_1\).
  • Set \(k=k+1\).
• Print Hemos alcanzado el máximo número de iteraciones. El método de la potencia no converge en este caso.

Dada una matriz \(\mathbf{A}\), el método de la potencia nos proporciona el vector propio junto con el vector propio correspondiente de módulo máximo de dicha matriz.

Sin embargo, el método de la potencia puede modificarse para que halle el valor propio junto con el vector propio correspondiente más cercano a un determinado valor real \(a\).

Antes de ver cómo poder realizar dicha modificación necesitamos el resultado siguiente:

Proposición.

Sea \(\mathbf{A}\) una matriz \(n\times n\) con valores propios \(\lambda_1,\lambda_2,\ldots,\lambda_n\) junto con vectores propios asociados \(\mathbf{v}^{(1)},\mathbf{v}^{(2)},\ldots,\mathbf{v}^{(n)}\). Sea \(a\) un número real. Entonces la matriz \((\mathbf{A}-a\mathbf{I})^{-1}\) tiene como valores propios \(\frac{1}{a-\lambda_1},\frac{1}{a-\lambda_2},\ldots,\frac{1}{a-\lambda_n}\) junto con los mismos vectores propios anteriores asociados a dichos valores propios.

Dado un valor real \(a\), usando la proposición anterior, si aplicamos el método de la potencia a la matriz \((\mathbf{A}-a\mathbf{I})^{-1}\) hallaremos el valor propio de la matriz \(\mathbf{A}\) más cercano a \(a\) junto con el vector propio correspondiente.

Veamos por qué: sean \(\lambda_1,\ldots,\lambda_n\) los valores propios de la matriz \(\mathbf{A}\). La proposición anterior dice que \(\frac{1}{\lambda_1-a},\ldots,\frac{1}{\lambda_n-a}\) son los valores propios de la matriz \((\mathbf{A}-a\mathbf{I})^{-1}\). Por tanto el método de la potencia hallará el valor propio \(\lambda_j\) tal que: \[ \left|\frac{1}{\lambda_j-a}\right|=\frac{1}{|\lambda_j-a|}=\max_{i=1,\ldots,n}\left|\frac{1}{\lambda_i-a}\right|=\frac{1}{\min_{i=1,\ldots,n}|\lambda_i-a|}, \] es decir, \(|\lambda_j-a|=\min_{i=1,\ldots,n}|\lambda_i-a|\), lo que significa que \(\lambda_j\) es el valor propio de la matriz \(\mathbf{A}\) más cercano al valor \(a\).

• INPUT matriz del sistema \(\mathbf{A}=(a_{ij})_{i=1,\ldots,n,j=1,\ldots,n}\), valor inicial \(\mathbf{X0}=\mathbf{x}\) con \(\|\mathbf{x}\|_\infty =1\), valor \(a\), error absoluto o tolerancia TOL, número máximo de iteraciones Nmax.
• Set \(k=1\) (inicializamos el contador de la iteraciones)
• Hallamos el entero más pequeño \(m\) tal que \(|x_m|=\|\mathbf{x}\|_\infty\).

• While \(k\leq Nmax\) (mientras no se llegue al número máximo de iteraciones, hacer lo siguiente)
  • Resolver el sistema \((\mathbf{A}-a\mathbf{I})\mathbf{y}=\mathbf{x}\)
  • If el sistema anterior no tiene solución única, significa que \(a\) es un valor propio de la matriz \(\mathbf{A}\).
  • Set \(\mu =y_m\).
  • Hallamos el entero más pequeño \(m\) tal que \(|y_m|=\|\mathbf{y}\|_\infty\).
  • Set \(\mathbf{x}_1=\frac{\mathbf{y}}{\|\mathbf{y}\|_\infty}\)
  • Set error=\(\min\left\{\left\|\mathbf{x}-\mathbf{x}_1\right\|_\infty,\left\|\mathbf{x}+\mathbf{x}_1\right\|_\infty\right\}\).
  • If error<TOL print(\(a+\frac{1}{\mu},\mathbf{x}\)). (El método ha acabado. Tenemos el valor propio de la matriz \(\mathbf{A}\) más cercano a \(a\) y un vector propio con un error menor que la tolerancia permitida.)

• • Set \(\mathbf{x}=\mathbf{x}_1\).
  • Set \(k=k+1\).
• Print Hemos alcanzado el máximo número de iteraciones. El método de la potencia no converge en este caso.

Observaciones:

• El valor \(\mathbf{y}\) en el paso \(k\)-ésimo del método de la potencia aplicado a la matriz \((\mathbf{A}-a\mathbf{I})^{-1}\) es \(\mathbf{y}=(\mathbf{A}-a\mathbf{I})^{-1}\mathbf{x}\). Sin embargo, hallar una inversa es muy costoso desde el punto de vista numérico. En su lugar, como \((\mathbf{A}-a\mathbf{I})\mathbf{y}=\mathbf{x}\), hemos de resolver un sistema de ecuaciones. Se suele usar el método de Gauss con pivotaje o la descomposición \(LU\). Como la matriz del sistema es la misma en todas las iteraciones, la descomposición \(LU\) o la transformación de la matriz del sistema en una matriz triangular superior sólo se realiza una vez.
• El método de la potencia aplicado a la matriz \((\mathbf{A}-a\mathbf{I})^{-1}\) nos da el valor propio de la forma \(\mu_j =\frac{1}{\lambda_j-a}\). Por tanto, \(\lambda_j=\frac{1}{\mu_j}+a\) que es el valor que nos interesa cuando el error es menor que la tolerancia.

Los métodos de la potencia y el de la potencia inversa nos permiten hallar valores propios de una matriz \(\mathbf{A}\) que cumplan una serie de condiciones: el primero nos calcula el valor propio de módulo máximo y el segundo nos calcula el valor propio más cercano a un número concreto \(a\).

En esta sección vamos a ver una técnica denominada deflación de Wielandt que nos permite hallar los valores propios restantes de una matriz \(\mathbf{A}\) una vez hallado uno de ellos.

Concretamente, dada una matriz \(\mathbf{A}\) de la que hemos hallado un valor propio \(\lambda_1\) junto con un correspondiente vector propio \(\mathbf{v}^{(1)}\), hallaremos otra matriz \(\mathbf{B}\) que tiene los mismos valores propios que la matriz \(\mathbf{A}\) excepto el valor propio \(\lambda_1\). Además, los vectores propios de la matriz original \(\mathbf{A}\) se podrán hallar a partir de los vectores propios de la nueva matriz \(\mathbf{B}\):

Teorema. Deflación de Wielandt.

Sea \(\mathbf{A}\) una matriz con valores propios \(\lambda_1,\lambda_2\ldots,\lambda_n\) junto con vectores propios asociados \(\mathbf{v}^{(1)},\mathbf{v}^{(2)},\ldots,\mathbf{v}^{(n)}\). Supongamos que el valor propio \(\lambda_1\) tiene multiplicidad \(1\). Sea \(\mathbf{x}\) un vector que verifica \(\mathbf{x}^\top \mathbf{v}^{(1)}=1\). Entonces la matriz: \[ \mathbf{B}=\mathbf{A}-\lambda_1\mathbf{v}^{(1)}\mathbf{x}^\top, \] tiene valores propios \(0,\lambda_2,\ldots,\lambda_n\) con vectores asociados \(\mathbf{v}^{(1)},\mathbf{w}^{(2)},\ldots,\mathbf{w}^{(n)}\), donde la relación entre los vectores propios \(\mathbf{v}^{(i)}\) y \(\mathbf{w}^{(i)}\) es la siguiente: \[ \mathbf{v}^{(i)}=(\lambda_i-\lambda_1)\mathbf{w}^{(i)}+\lambda_1 (\mathbf{x}^\top \mathbf{w}^{(i)})\mathbf{v}^{(1)}, \] para \(i=2,3,\ldots,n\).

Una posible elección del vector \(\mathbf{x}\) del Teorema anterior podría ser la siguiente. Sea \(v_i^{(1)}\neq 0\) una componente del vector \(\mathbf{v}^{(1)}\) diferente de \(0\). Entonces, consideramos: \[ \mathbf{x}=\frac{1}{\lambda_1 v_i^{(1)}}\begin{bmatrix}a_{i1}\\ a_{i2}\\\vdots\\ a_{in}\end{bmatrix}, \] donde el vector columna \(\begin{bmatrix}a_{i1}\\ a_{i2}\\\vdots\\ a_{in}\end{bmatrix}\) es la fila \(i\)-ésima de la matriz \(\mathbf{A}\).

Véamoslo:

\[ \mathbf{x}^\top \mathbf{v}^{(1)}=\frac{1}{\lambda_1 v_i^{(1)}}(a_{i1}, a_{i2},\ldots, a_{in})\begin{bmatrix}v_{1}^{(1)}\\ v_{2}^{(1)}\\\vdots\\ v_{n}^{(1)}\end{bmatrix}=\frac{1}{\lambda_1 v_i^{(1)}}\sum_{j=1}^n a_{ij}v_j^{(1)}. \] Como el vector \(\mathbf{v}^{(1)}\) es un vector propio de la matriz \(\mathbf{A}\) de valor propio \(\lambda_1\), tenemos que \(\mathbf{A}\mathbf{v}^{(1)}=\lambda_1\mathbf{v}^{(1)}\). La componente \(i\)-ésima de la igualdad anterior será: \(\displaystyle\sum_{j=1}^n a_{ij}v_j^{(1)}=\lambda_1 v_i^{(1)}\). Sustituyendo a la expresión anterior, tenemos que: \[ \mathbf{x}^\top \mathbf{v}^{(1)}=\frac{1}{\lambda_1 v_i^{(1)}}\sum_{j=1}^n a_{ij}v_j^{(1)}=\frac{1}{\lambda_1 v_i^{(1)}}\lambda_1 v_i^{(1)}=1, \] tal como queríamos ver.

Para poder aplicar el Teorema anterior, dada una matriz \(\mathbf{A}\) imaginemos que hemos calculado un valor propio \(\lambda_1\) junto con un vector propio asociado \(\mathbf{v}^{(1)}\) usando por ejemplo el método de la potencia o de la potencia inversa.

Entonces, para hallar los demás, calculamos la matriz \(\mathbf{B}\) dada por el Teorema anterior y aplicamos el método de la potencia a la matriz \(\mathbf{B}\). Como en dicha matriz, el valor propio \(\lambda_1\) ha sido “sustituido” por el valor propio \(0\), dicho valor propio nunca será hallado ya que tiene módulo mínimo y recordemos que el método de la potencia calcula el valor propio de módulo máximo.

Además, como iremos calculando los vectores propios \(\mathbf{w}^{(i)}\), \(i=2,\ldots,n\) de la matriz \(\mathbf{B}\) la relación entre los vectores propios que da el teorema permite calcular los vectores propios \(\mathbf{v}^{(i)}\) de la matriz \(\mathbf{A}\).

Por último, veamos cómo reducir la dimensión de la matriz \(\mathbf{B}\) en una unidad, simplificando de esta manera el coste computacional para el cálculo de los valores propios restantes de la matriz \(\mathbf{A}\).

Sea \(\lambda\neq 0\) un valor propio distinto de \(0\) de la matriz \(\mathbf{B}\) junto con un vector propio asociado \(\mathbf{w}\). Supongamos que hemos elegido el vector \(\mathbf{x}\) como hemos indicado antes: \(\mathbf{x}=\frac{1}{\lambda_1 v_i^{(1)}}\begin{bmatrix}a_{i1}\\ a_{i2}\\\vdots\\ a_{in}\end{bmatrix}.\)

Para fijar ideas y no perdernos con las notaciones, supongamos sin pérdida de generalidad que \(v_1^{(1)}\neq 0\), es decir, el vector \(\mathbf{x}\) sería: \(\mathbf{x}=\frac{1}{\lambda_1 v_1^{(1)}}\begin{bmatrix}a_{11}\\ a_{12}\\\vdots\\ a_{1n}\end{bmatrix}.\)

Entonces, en este caso \(w_1=0\), es decir, la primera componente del vector propio de valor propio \(\lambda\) de la matriz \(\mathbf{B}\) es nulo. Veámoslo: \[ \mathbf{B}\mathbf{w}=\lambda\mathbf{w},\ \Rightarrow (\mathbf{A}-\lambda_1\mathbf{v}^{(1)}\mathbf{x}^\top)\mathbf{w}=\lambda\mathbf{w}. \]

La matriz \(\mathbf{v}^{(1)}\mathbf{x}^\top\) sería: \[ \begin{align*} \mathbf{v}^{(1)}\mathbf{x}^\top & =\begin{bmatrix} v_1^{(1)}\\ v_2^{(1)}\\\vdots\\ v_n^{(1)}\end{bmatrix}(x_1,x_2,\ldots,x_n)=\begin{bmatrix}v_1^{(1)}x_1 &v_1^{(1)}x_2 & \ldots &v_1^{(1)}x_n\\ v_2^{(1)}x_1 &v_2^{(1)}x_2 & \ldots &v_2^{(1)}x_n\\ \vdots & \vdots & \vdots &\vdots\\ v_n^{(1)}x_1 &v_n^{(1)}x_2 & \ldots &v_n^{(1)}x_n \end{bmatrix}\\ & = \begin{bmatrix}v_1^{(1)}\frac{1}{\lambda_1 v_1^{(1)}}a_{11} &v_1^{(1)}\frac{1}{\lambda_1 v_1^{(1)}}a_{12} & \ldots &v_1^{(1)}\frac{1}{\lambda_1 v_1^{(1)}}a_{1n}\\ v_2^{(1)}x_1 &v_2^{(1)}x_2 & \ldots &v_2^{(1)}x_n\\ \vdots & \vdots & \vdots &\vdots\\ v_n^{(1)}x_1 &v_n^{(1)}x_2 & \ldots &v_n^{(1)}x_n \end{bmatrix} \end{align*} \]

\[ \mathbf{v}^{(1)}\mathbf{x}^\top = \begin{bmatrix}\frac{a_{11}}{\lambda_1} &\frac{a_{12}}{\lambda_1} & \ldots &\frac{a_{1n}}{\lambda_1}\\ v_2^{(1)}x_1 &v_2^{(1)}x_2 & \ldots &v_2^{(1)}x_n\\ \vdots & \vdots & \vdots &\vdots\\ v_n^{(1)}x_1 &v_n^{(1)}x_2 & \ldots &v_n^{(1)}x_n \end{bmatrix} \] Entonces como \((\mathbf{A}-\lambda_1\mathbf{v}^{(1)}\mathbf{x}^\top)\mathbf{w}=\lambda \mathbf{w}\), la primera componente de la igualdad anterior será: \[ \sum_{j=1}^n \left(a_{1j}-\lambda_1 \frac{a_{1j}}{\lambda_1}\right)w_j =\lambda w_1,\ \sum_{j=1}^n \left(a_{1j}-a_{1j}\right)w_j =\lambda w_1,\ \Rightarrow 0=w_1, \] tal como queríamos ver.

Entonces en la igualdad \(\mathbf{B}\mathbf{w}=\lambda\mathbf{w}\), la primera columna de \(\mathbf{B}\) no interviene para nada ya que el vector \(\mathbf{w}\) es de la forma: \(\mathbf{w}=\begin{bmatrix}0\\ w_2\\\vdots\\ w_n\end{bmatrix}\) y si escribimos la igualdad \(\mathbf{B}\mathbf{w}=\lambda\mathbf{w}\) en componentes, tendremos: \[ \begin{bmatrix} b_{11} & b_{12}& \ldots & b_{1n}\\ b_{21} & b_{22}& \ldots & b_{2n}\\ \vdots & \vdots & \vdots & \vdots \\ b_{n1} & b_{n2} & \ldots & b_{nn} \end{bmatrix} \begin{bmatrix} 0\\ w_2\\ \vdots \\ w_n \end{bmatrix} =\lambda\begin{bmatrix} 0\\ w_2\\ \vdots \\ w_n \end{bmatrix}. \]

Fijémonos que podríamos cambiar la primera columna \(\begin{bmatrix}b_{11}\\b_{21}\\\vdots\\ b_{n1}\end{bmatrix}\) de la matrix \(\mathbf{B}\) como deseamos y la igualdad anterior seguiría siendo cierta.

La idea entonces es considerar la matriz \(\mathbf{B}'\) eliminando la primera fila y columna de la matriz anterior \(\mathbf{B}\) de cara a calcular los valores propios y vectores propios: \[ \mathbf{B}'=\begin{bmatrix} b_{12}& \ldots & b_{1n}\\ b_{22}& \ldots & b_{2n}\\ \vdots & \vdots & \vdots \\ b_{n2} & \ldots & b_{nn} \end{bmatrix}. \]

De la igualdad \(\mathbf{B}\mathbf{w}=\lambda\mathbf{w}\), podemos deducir la igualdad \(\mathbf{B}'\mathbf{w}'=\lambda\mathbf{w}'\), con \(\mathbf{w}'=\begin{bmatrix} w_2\\ \vdots \\ w_n \end{bmatrix}\), es decir, \(\mathbf{w}'\) es un vector propio de la matriz \(\mathbf{B}'\) del mismo valor propio \(\lambda\). Entonces \(\mathbf{B}'\) tiene los mismos valores propios que la matriz \(\mathbf{B}\). Por tanto, podemos aplicar el método de la potencia a la matriz \(\mathbf{B}'\) y una vez hallado \(\lambda\) y \(\mathbf{w}'\) hallamos \(\mathbf{w}\) como \(\mathbf{w}=\begin{bmatrix}0\\\mathbf{w}'\end{bmatrix}\).

Si en lugar de ser la primera componente del vector propio \(\mathbf{v}^{(1)}\) la que es diferente de cero, \(v_1^{(1)}\neq 0\), es la \(i\)-ésima, \(v_i^{(1)}\neq 0\), hacemos el mismo proceso anterior pero en lugar de eliminar la primera fila y primera columna para hallar la matriz \(\mathbf{B}'\), eliminamos la fila y la columna \(i\)-ésimas.

Veamos cómo hallar todos los valores propios y vectores propios de una matriz \(\mathbf{A}\).

Primeramente, de cara a simplificar los cálculos, como la matriz \(\mathbf{B}\) vale \(\mathbf{B}=\mathbf{A}-\lambda_1\mathbf{v}^{(1)}\mathbf{x}^\top\) y \(\mathbf{x}=\frac{a_{i\cdot}}{\lambda_1 v_i^{(1)}}\), donde \(a_{i\cdot}\) indica la fila \(i\)-ésima de la matriz \(\mathbf{A}\), podemos escribir la matriz \(\mathbf{B}\) como \(\mathbf{B}=\mathbf{A}-\frac{1}{v_i^{(1)}}\mathbf{v}^{(1)}a_{i\cdot}^\top\)

• Paso \(1\). Hallamos el valor propio de módulo máximo \(\lambda_1\) junto con un correspondiente vector propio \(\mathbf{v}^{(1)}\) usando el método de la potencia.

• Paso \(2\). Consideramos una componente de \(\mathbf{v}^{(1)}\) diferente de cero, \(v_i^{(1)}\neq 0\). Para minimizar los errores de redondeo consideramos \(i\) tal que \(|v_i^{(1)}|=\max_{j=1,\ldots,n}|v_j^{(1)}|\).
  • Definimos \(\mathbf{x}=\frac{1}{v_i^{(1)}}\begin{bmatrix}a_{i1}\\ a_{i2}\\\vdots\\ a_{in}\end{bmatrix}\).
  • Consideramos la matriz \(\mathbf{B}=\mathbf{A}-\frac{1}{v_i^{(1)}}\mathbf{v}^{(1)}\mathbf{a_{i\cdot}}^\top\).
• Paso \(3\). Eliminamos las filas y columnas \(i\)-ésimas de la matriz \(\mathbf{B}\) obteniendo la matriz \(\mathbf{B}'\).

• Paso \(4\). Aplicamos el método de la potencia a la matriz \(\mathbf{B}'\) hallando el siguiente valor propio de módulo máximo \(\lambda_2\) de la matriz \(\mathbf{A}\). Sea \(\mathbf{w}'\) el vector propio asociado de la matriz \(\mathbf{B}'\).
• Paso \(5\). Calculamos el vector propio \(\mathbf{w}\) de la matriz \(\mathbf{B}\): \(\mathbf{w}=\begin{bmatrix} w_1'\\\vdots \\ 0(\mathrm{componente\ i})\\ \vdots\\ w_n'\end{bmatrix}\).
• Paso \(6\). Calculamos el vector propio \(\mathbf{v}^{(2)}\) de valor propio \(\lambda_2\) de la matriz \(\mathbf{A}\): \(\mathbf{v}^{(2)}=(\lambda_2-\lambda_1)\mathbf{w}+\frac{1}{v_i^{(1)}}(\mathbf{a_{i\cdot}}^\top \mathbf{w})\mathbf{v}^{(1)}.\)

Primero veamos los lemas siguientes:

Lema 1.

Sean \(\mathbf{a},\mathbf{b}\in\mathbb{R}^n\), entonces: \((\mathbf{a}^\top \mathbf{b})\mathbf{a}=\mathbf{a}(\mathbf{b}^\top \mathbf{a})\).

La demostración del lema anterior se realiza usando que el producto escalar es conmutativo.

Lema 2.

Sean \(\mathbf{a},\mathbf{b},\mathbf{c}\in\mathbb{R}^n\), entonces: \((\mathbf{a}\mathbf{b}^\top)\mathbf{c}=\mathbf{a}(\mathbf{b}^\top \mathbf{c})=(\mathbf{b}^\top \mathbf{c})\mathbf{a}\).

• INPUT matriz del sistema \(\mathbf{A}=(a_{ij})_{i=1,\ldots,n,j=1,\ldots,n}\), aproximación del valor propio \(\lambda\) con vector propio \(\mathbf{v}\), error absoluto o tolerancia TOL, número máximo de iteraciones Nmax.
• Set \(i\) el menor entero entre \(1\) y \(n\) tal que \(|v_i|=\max{|v_j|}\).
• For \(k=1,\ldots,n\) (definimos la matriz \(\mathbf{B}=\mathbf{A}-\frac{1}{v_i} \mathbf{v}\mathbf{x}^\top\))
  • For \(j=1,\ldots,n\)
    • Set \(b_{kj}=a_{kj}-\frac{v_k a_{ij}}{v_i}\).
• Set \(\mathbf{B}'=\mathbf{B}[-i,-i]\) (definimos la matriz \(\mathbf{B}'\) como la matriz \(\mathbf{B}\) de la que hemos eliminados la fila y columna \(i\)-ésima)

• Aplicamos método de la potencia a la matriz \(\mathbf{B}'\) con aproximación inicial \(\mathbf{x}\). Obtenemos valor propio \(\mu\) y vector propio \(\mathbf{w}'\in\mathbb{R}^{n-1}\).
• If \(i\neq 1\) y \(i\neq n\) (definimos el vector \(\mathbf{w}\))
  • For \(j=1,\ldots,i-1\)
    • Set \(w_j=w_j'\)
  • Set \(w_i=0\).
  • For \(j=i+1,\ldots,n\)
    • Set \(w_j=w_{j-1}'\)

• Else
  • If \(i=1\)
    • Set \(w_1=0\)
    • For \(j=2,\ldots,n\)
      • Set \(w_j=w_{j-1}'\)
  • Else (\(i=n\))
    • For \(j=1,\ldots,n-1\)
      • Set \(w_j=w_j'\)
    • Set \(w_n=0\)

• Set \(\mathbf{v}^{(i)}=(\mu-\lambda)\mathbf{w}+\left(\sum_{j=1}^n \frac{a_{ij}w_j}{v_i}\right)\mathbf{v}\) (calculamos el vector propio \(\mathbf{v}^{(i)}\) de la matriz \(\mathbf{A}\))
• Print \(\mu\), \(\mathbf{v}^{(i)}\) (damos el nuevo valor propio \(\mu\) con el vector propio correspondiente \(\mathbf{v}^{(i)}\))