El método de la potencia visto anteriormente halla, dada una matriz \(\mathbf{A}\), \(n\times n\) el valor propio de módulo máximo junto con el vector propio correspondiente.

Para hallar los demás valores y vectores propios, podemos usar la técnica de la deflación.

En esta sección, cambiaremos la perspectiva de cálculo de valores propios. En lugar de ir calculándolos de uno en uno, vamos a calcularlos todos de una vez. Son los denominados métodos de ortogonalización.

Fijarse que nos concentraremos en el cálculo de todos los valores propios. Los vectores propios pueden hallarse resolviendo un sistema lineal o aplicar el método de la potencia visto anteriormente.

Concretamente, vamos a ver un método que, dada una matriz \(\mathbf{A}\), la transforma en una matriz tridiagonal simétrica o en una matriz Hessenberg superior usando una transformación de similaridad ortogonal.

Es decir, dada una matriz \(\mathbf{A}\), vamos a hallar una matriz ortogonal \(\mathbf{H}\) tal que \(\mathbf{H}^\top\mathbf{A}\mathbf{H}=\mathbf{B}\), donde la matriz \(\mathbf{B}\) será tridiagonal simétrica en el caso en que la matriz original \(\mathbf{A}\) sea simétrica o Hessenberg superior en caso contrario.

Recordemos que en el capítulo de Álgebra lineal numérica vimos que las matrices \(\mathbf{A}\) y \(\mathbf{B}\) tienen los mismos valores propios.

Aunque hemos dicho que no nos concentraremos en el cálculo de los vectores propios, existe la relación siguiente entre dichos vectores propios: si \(\mathbf{w}\) es un vector propio de la matriz \(\mathbf{B}\) de valor propio \(\lambda\), entonces \(\mathbf{v}=\mathbf{H}\mathbf{w}\) es un vector propio de la matriz original \(\mathbf{A}\) del mismo valor propio.

Veámoslo: como \(\mathbf{v}\) es un vector propio de la matriz \(\mathbf{B}\) de valor propio \(\lambda\), se cumple \(\mathbf{B}\mathbf{w}=\lambda\mathbf{w}\), entonces, como \(\mathbf{B}=\mathbf{H}^\top\mathbf{A}\mathbf{H}\), multiplicando por \(\mathbf{H}\) ambos miembros nos queda \(\mathbf{H}\mathbf{B}=\mathbf{H}\mathbf{H}^\top\mathbf{A}\mathbf{H}=\mathbf{I}\mathbf{A}\mathbf{H}=\mathbf{A}\mathbf{H}\). Veamos que el vector \(\mathbf{v}=\mathbf{H}\mathbf{w}\) es un vector propio de la matriz \(\mathbf{A}\) de valor propio \(\lambda\): \[ \mathbf{A}\mathbf{v}=\mathbf{A}\mathbf{H}\mathbf{w}=\mathbf{H}\mathbf{B}\mathbf{w}=\mathbf{H}\lambda\mathbf{w}=\lambda\mathbf{H}\mathbf{w}=\lambda\mathbf{v}, \] como queríamos ver.

Sintentizando lo dicho, hallar los valores propios de una matriz es equivalente a hallar los valores propios de matrices triagonales simétricas (caso de matrices simétricas) o Hessenberg superiores.

Entonces, veamos en primer lugar cómo hallar valores propios de matrices tridiagonales simétricas y Hessenberg superiores.

Sea \(\mathbf{A}\) una matriz tridiagonal simétrica:

\[\mathbf{A}=\begin{bmatrix}a_{1}&b_{2}&0&0&\ldots&0 \\b_{2}&a_{2}&b_{3}&0&\ldots&0 \\0&b_{3}&a_{3}&b_{4}&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&b_{n-1}&a_{n-1}&b_{n} \\0&\ldots&\ldots&0&b_{n}&a_{n} \\\end{bmatrix},\] de la que queremos hallar sus valores propios.

Vamos a construir una sucesión de matrices \(\mathbf{A}^{(1)}=\mathbf{A},\mathbf{A}^{(2)},\ldots,\mathbf{A}^{(n)},\ldots\), todas semejantes a la matriz original \(\mathbf{A}\) y tridiagonales simétricas tal que los elementos de las matrices \(\mathbf{A}^{(n)}\) fuera de la diagonal tienen en valor absoluto valores más pequeños a medida que aumenta \(n\). Es decir, \[ \lim_{n\to\infty}\max_{i\neq j}|a_{ij}^{(n)}|=0. \] Concretamente, se cumplira que para cada valor de \(i\), existirá una matriz ortogonal Hessenberg superior \(\mathbf{Q}^{(i)}\) tal que \(\mathbf{A}^{(i+1)}=(\mathbf{Q}^{(i)})^\top\mathbf{A}^{(i)}\mathbf{Q}^{(i)}\).

Recordemos que una matriz Hessenberg superior era de la forma:

\[\mathbf{Q}=\begin{bmatrix}q_{11}&q_{12}&q_{13}&q_{14}&\ldots&q_{1n} \\q_{21}&q_{22}&q_{23}&q_{24}&\ldots&q_{2n} \\0&q_{32}&q_{33}&q_{34}&\ldots&q_{3n} \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&q_{n-1,n-2}&q_{n-1,n-1}&q_{n-1,n} \\0&\ldots&\ldots&0&q_{n,n-1}&q_{nn} \\\end{bmatrix},\]

De esta manera, los valores propios de las matrices de la sucesión serán los mismos y como la sucesión de matrices \((\mathbf{A}^{(n)})_{n\geq 0}\) tiende a una matriz diagonal, podemos aproximar los valores propios de la matriz original por los elementos diagonales de la matriz \(\mathbf{A}^{(n)}\).

Los valores propios de una matriz diagonal \(\mathbf{D}\),

\[\mathbf{D}=\begin{bmatrix}\lambda_1&0&\ldots&0 \\0&\lambda_2&\ldots&0 \\\vdots&\vdots&\ddots&\vdots \\0&\ldots&0&\lambda_n \\\end{bmatrix}.\] son precisamente los valores \(\lambda_i\), \(i=1,\ldots,n\) ya que si hacemos \(\mathrm{det}(\mathbf{D}-\lambda_i\mathbf{I})=0\) al tener dicha matriz la fila \(i\)-ésima nula.

Para hallar la sucesión de matrices \((\mathbf{A}^{(n)})_{n\geq 0}\) se siguen los pasos siguientes:

• Paso \(1\). Sea \(\mathbf{A}^{(1)}=\mathbf{A}\). Descomponemos la matriz \(\mathbf{A}^{(1)}\) en el producto de una matriz ortogonal Hessenberg superior \(\mathbf{Q}^{(1)}\) y una matriz triangular superior \(\mathbf{S}^{(1)}\): \(\mathbf{A}^{(1)}=\mathbf{Q}^{(1)}\mathbf{S}^{(1)}\).

• Paso \(2\). Definimos \(\mathbf{A}^{(2)}=\mathbf{S}^{(1)}\mathbf{Q}^{(1)}\), el producto de la matriz triangular superior por la matriz ortogonal Hessenberg superior. Como \(\mathbf{A}^{(1)}=\mathbf{Q}^{(1)}\mathbf{S}^{(1)}\), multiplicando por \((\mathbf{Q}^{(1)})^\top\), nos queda que \((\mathbf{Q}^{(1)})^\top\mathbf{A}^{(1)}=\mathbf{S}^{(1)}\). Por tanto, \(\mathbf{A}^{(2)}=(\mathbf{Q}^{(1)})^\top\mathbf{A}^{(1)}\mathbf{Q}^{(1)}\), es decir, las matrices \(\mathbf{A}^{(1)}\) y \(\mathbf{A}^{(2)}\) son semejantes y tienen los mismos valores propios.

• La matriz \(\mathbf{A}^{(2)}\)
  • será simétrica ya que \[ \begin{align*} (\mathbf{A}^{(2)})^\top & =((\mathbf{Q}^{(1)})^\top\mathbf{A}^{(1)}\mathbf{Q}^{(1)})^\top =(\mathbf{Q}^{(1)})^\top (\mathbf{A}^{(1)})^\top ((\mathbf{Q}^{(1)})^\top)^\top \\ & = (\mathbf{Q}^{(1)})^\top\mathbf{A}^{(1)}\mathbf{Q}^{(1)}=\mathbf{A}^{(2)}. \end{align*} \]
  • será Hessenberg superior ya que como \(\mathbf{A}^{(2)}\) es el producto de una matriz triangular superior por una matriz Hessenberg superior, \(\mathbf{A}^{(2)}=\mathbf{S}^{(1)}\mathbf{Q}^{(1)}\), la matriz \(\mathbf{A}^{(2)}\) será Hessenberg superior pero como es simétrica, resulta que la matriz \(\mathbf{A}^{(2)}\) será tridiagonal simétrica.
• A continuación descomponemos la matriz \(\mathbf{A}^{(2)}\) en el producto de una matriz ortogonal Hessenberg superior \(\mathbf{Q}^{(2)}\) y una matriz triangular superior \(\mathbf{S}^{(2)}\): \(\mathbf{A}^{(2)}=\mathbf{Q}^{(2)}\mathbf{S}^{(2)}\).

• Paso \(3\). Definimos \(\mathbf{A}^{(3)}=\mathbf{S}^{(2)}\mathbf{Q}^{(2)}\), etc.
• \(\vdots\)
• Paso \(i\). Definimos \(\mathbf{A}^{(i)}=\mathbf{S}^{(i-1)}\mathbf{Q}^{(i-1)}\). Usando un razonamiento inductivo y usando lo visto en el paso \(2\), la matriz \(\mathbf{A}^{(i)}\) será tridiagonal simétrica. Descomponemos la matriz \(\mathbf{A}^{(i)}\) en el producto de una matriz ortogonal Hessenberg superior \(\mathbf{Q}^{(i)}\) y una matriz triangular superior \(\mathbf{S}^{(i)}\): \(\mathbf{A}^{(i)}=\mathbf{Q}^{(i)}\mathbf{S}^{(i)}\).
• Paso \(i+1\). Definimos \(\mathbf{A}^{(i+1)}=\mathbf{S}^{(i)}\mathbf{Q}^{(i)}\). Como \(\mathbf{A}^{(i)}=\mathbf{Q}^{(i)}\mathbf{S}^{(i)}\), multiplicando por \((\mathbf{Q}^{(i)})^\top\), nos queda que \((\mathbf{Q}^{(i)})^\top\mathbf{A}^{(i)}=\mathbf{S}^{(i)}\). Por tanto, \(\mathbf{A}^{(i+1)}=(\mathbf{Q}^{(i)})^\top\mathbf{A}^{(i)}\mathbf{Q}^{(i)}\), es decir, las matrices \(\mathbf{A}^{(i+1)}\) y \(\mathbf{A}^{(i)}\) son semejantes y tienen los mismos valores propios. Además, usando el mismo razonamiento anterior, la matriz \(\mathbf{A}^{(i+1)}\) sería tridiagonal simétrica, etc.

La descomposición de la matriz tridiagonal simétrica del paso \(i\)-ésimo en una matriz ortogonal Hessenberg superior \(\mathbf{Q}^{(i)}\) y una matriz \(\mathbf{S}^{(i)}\) triangular superior está basada en las llamadas matrices de rotación:

Definición de matriz de rotación.

Una matriz \(\mathbf{R}_{ij}=(r_{kl}^{(ij)})_{k,l=1,\ldots,n}\) de dimensiones \(n\times n\) es una matriz de rotación si existen dos valores \(i\) y \(j\) entre \(1\) y \(n\) tal que la matriz \(\mathbf{R}_{ij}\) y la matriz identidad \(\mathbf{I}\) coinciden excepto en las filas y columnas \(i\) y \(j\). Para las filas y columnas \(i\) y \(j\), existe un valor \(\alpha\in [0,2\pi)\) con: \[ r_{kl}=\begin{cases}\cos(\alpha), & \mathrm{si\ }k=l=i\mathrm{\ o\ }k=l=j,\\ \sin(\alpha), & \mathrm{si\ }k=i,l=j,\\ -\sin(\alpha), & \mathrm{si\ }k=j,l=i,\\ 1, &\mathrm{si\ }k=l\neq i,j,\\ 0, &\mathrm{en\ caso\ contrario.} \end{cases} \]

\[\mathbf{R}_{ij}=\begin{bmatrix}1&\ldots&0&\ldots&0&\ldots&0 \\\vdots&\ddots&\vdots&\vdots&\vdots&\vdots&\vdots \\0&\ldots&\cos(\alpha)&\ldots&\sin(\alpha)&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots \\0&\ldots&-\sin(\alpha)&\ldots&\cos(\alpha)&\ldots&0 \\\vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots \\0&\ldots&0&\ldots&0&\ldots&1 \\\end{bmatrix}.\]

Las matrices de rotación son matrices ortogonales tal como afirma la proposición siguiente:

Proposición.

Las matrices de rotación son ortogonales.

Supongamos que estamos en el paso \(i\)-ésimo y tenemos la matriz \(\mathbf{A}^{(i)}\) tridiagonal simétrica:

\[\mathbf{A}^{(i)}=\begin{bmatrix}a_{1}&b_{2}&0&0&\ldots&0 \\b_{2}&a_{2}&b_{3}&0&\ldots&0 \\0&b_{3}&a_{3}&b_{4}&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&b_{n-1}&a_{n-1}&b_{n} \\0&\ldots&\ldots&0&b_{n}&a_{n} \\\end{bmatrix}.\]

• Paso \(1\). En primer lugar, vamos a considerar la matriz de rotación \(\mathbf{R}_{12}\) y vamos a hallar el valor \(\alpha\) tal que la componente \((2,1)\) del producto de matrices \(\mathbf{R}_{12}\mathbf{A}^{(i)}\) sea \(0\): \[ \mathbf{R}_{12}\mathbf{A}^{(i)}=\begin{bmatrix}\cos(\alpha)&\sin(\alpha)&0&\ldots&0 \\-\sin(\alpha)&\cos(\alpha)&0&\ldots&0 \\\vdots&\vdots&\ddots&\vdots&\vdots \\0&0&0&\ldots&1 \\\end{bmatrix}\begin{bmatrix}a_{1}&b_{2}&0&0&\ldots&0 \\b_{2}&a_{2}&b_{3}&0&\ldots&0 \\0&b_{3}&a_{3}&b_{4}&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&b_{n-1}&a_{n-1}&b_{n} \\0&\ldots&\ldots&0&b_{n}&a_{n} \\\end{bmatrix}. \] El elemento \((2,1)\) del producto anterior vale: \(-a_1\sin(\alpha_1)+b_2\cos(\alpha_1)\). Hallamos el valor \(\alpha_1\) tal que: \[ \begin{align*} & -a_1\sin(\alpha_1)+b_2\cos(\alpha_1)=0,\ \Rightarrow a_1\sin(\alpha_1)=b_2\cos(\alpha_1),\\ & \Rightarrow \tan(\alpha_1)=\frac{b_2}{a_1}. \end{align*} \]

Por tanto, \[ \cos(\alpha_1)=\sqrt{\frac{1}{1+\tan^2(\alpha_1)}}=\frac{a_1}{\sqrt{(a_1)^2+(b_2)^2}},\quad\sin(\alpha_1)=\frac{b_2}{\sqrt{(a_1)^2+(b_2)^2}}. \] El producto anterior será de la forma: \[ \mathbf{R}_{12}\mathbf{A}^{(i)}=\begin{bmatrix}*&*&*&*&\ldots&* \\0&a_2^{(2)}&*&*&\ldots&* \\0&b_3^{(2)}&*&*&\ldots&* \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&*&*&* \\0&\ldots&\ldots&0&*&* \\\end{bmatrix}, \] donde el producto anterior será una matriz Hessenberg superior al ser las matrices \(\mathbf{R}_{12}\) y \(\mathbf{A}^{(i)}\) matrices Hessenberg superiores.

• Paso \(2\). A continuación, consideramos la matriz \(\mathbf{R}_{23}\) y hallamos el valor \(\alpha_2\) tal que la componente \((3,2)\) producto de las matrices \(\mathbf{R}_{23}(\mathbf{R}_{12}\mathbf{A}^{(i)})\) sea \(0\).

Razonando de forma similar, el valor \(\alpha_2\) verificará: \[ \cos(\alpha_2)=\sqrt{\frac{1}{1+\tan^2(\alpha_2)}}=\frac{a_2^{(2)}}{\sqrt{(a_2^{(2)})^2+(b_3^{(2)})^2}},\quad\sin(\alpha_2)=\frac{b_3^{(2)}}{\sqrt{(a_2^{(2)})^2+(b_3^{(2)})^2}}. \] La matriz \(\mathbf{R}_{23}(\mathbf{R}_{12}\mathbf{A}^{(i)})\) será Hesenberg superior al ser el producto de dos matrices Hessenberg superiores: \[ \mathbf{R}_{23}(\mathbf{R}_{12}\mathbf{A}^{(i)})=\begin{bmatrix}*&*&*&*&\ldots&* \\0&*&*&*&\ldots&* \\0&0&*&*&\ldots&* \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&*&*&* \\0&\ldots&\ldots&0&*&* \\\end{bmatrix}. \]

De esta forma, en \(n-1\) pasos podemos hacer \(0\) las componentes \((2,1), (3,2),\ldots,(n,n-1)\) obteniendo el producto \(\mathbf{R}_{n-1,n}\cdots \mathbf{R}_{12}\mathbf{A}^{(i)}=\mathbf{S}^{(i)},\) donde la matriz \(\mathbf{S}^{(i)}\) sería triangular superior al ser Hessenberg superior y tener nulas las componentes \((2,1), (3,2),\ldots,(n,n-1)\).

Por último, al ser las matrices \(\mathbf{R}_{j,j+1}\), \(j=1,\ldots,n-1\) ortogonales, tenemos que: \[ \mathbf{A}^{(i)}=\mathbf{R}_{12}^\top\cdots\mathbf{R}_{n-1,n}^\top\mathbf{S}^{(i)}. \] La matriz \(\mathbf{Q}^{(i)}\) será, pues, \(\mathbf{Q}^{(i)}=\mathbf{R}_{12}^\top\cdots\mathbf{R}_{n-1,n}^\top\) y la matriz \(\mathbf{S}^{(i)}\), la matriz \(\mathbf{S}^{(i)}\) de la descomposición anterior.

Observación.

La matriz \(\mathbf{Q}^{(i)}=\mathbf{R}_{12}^\top\cdots\mathbf{R}_{n-1,n}^\top\) es ortogonal ya que: \[ \begin{align*} (\mathbf{Q}^{(i)})^\top\mathbf{Q}^{(i)} & =(\mathbf{R}_{12}^\top\cdots\mathbf{R}_{n-1,n}^\top)^\top\mathbf{R}_{12}^\top\cdots\mathbf{R}_{n-1,n}^\top\\ & =\mathbf{R}_{n-1,n}\cdots\mathbf{R}_{12}\mathbf{R}_{12}^\top\cdots\mathbf{R}_{n-1,n}^\top = \mathbf{I}, \end{align*} \] al ser las matrices \(\mathbf{R}_{j,j+1}\) ortogonales, \(j=1,\ldots,n-1\).

• INPUT matriz del sistema \(\mathbf{A}=\)\(\begin{bmatrix}a_{1}&b_{2}&0&0&\ldots&0 \\b_{2}&a_{2}&b_{3}&0&\ldots&0 \\0&b_{3}&a_{3}&b_{4}&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&b_{n-1}&a_{n-1}&b_{n} \\0&\ldots&\ldots&0&b_{n}&a_{n} \\\end{bmatrix}\).
• Set \(\mathbf{S}=\mathbf{A}\) (en la matriz \(\mathbf{S}\) iremos guardando los productos \(\mathbf{R}_{k,k+1}\cdots\mathbf{R}_{12}\mathbf{A}\))
• Set \(\mathbf{Q}=\mathbf{I}\) (en la matriz \(\mathbf{Q}\) iremos guardando los productos \(\mathbf{R}_{1,2}^\top\cdots\mathbf{R}_{k,k+1}^\top\))

• For \(k=1,\ldots,n-1\)
  • Set \(c=\frac{s_{k,k}}{\sqrt{(s_{k,k})^2+(s_{k+1,k})^2}}\).
  • Set \(s=\frac{s_{k+1,k}}{\sqrt{(s_{k,k})^2+(s_{k+1,k})^2}}\).
  • Set \(\mathbf{R}_{k,k+1}=(r_{ij})_{i,j=1,\ldots,n}=\mathbf{I}\) (definimos la matriz \(\mathbf{R}_{k,k+1}\) inicialmente como la identidad)

• (Cambiamos las componentes \((k,k),(k+1,k),(k,k+1)\) y \((k+1,k+1)\))
  • Set \(r_{k,k}=c\).
  • Set \(r_{k,k+1}=s\).
  • Set \(r_{k+1,k}=-s\).
  • Set \(r_{k+1,k+1}=c\).
  • Set \(\mathbf{S}=R_{k,k+1}\mathbf{S}\). (Cambiamos la matriz \(\mathbf{S}\))
  • Set \(\mathbf{Q}=\mathbf{Q}\mathbf{R}_{k,k+1}^\top\) (Cambiamos la matriz \(\mathbf{Q}\))
• Print \(\mathbf{Q}\), \(\mathbf{S}\).

• INPUT matriz del sistema \(\mathbf{A}=\mathbf{A}^{(1)}=\)\(\begin{bmatrix}a_1^{(1)}&b_{2}^{(1)}&0&0&\ldots&0 \\b_{2}^{(1)}&a_{2}^{(1)}&b_{3}^{(1)}&0&\ldots&0 \\0&b_{3}^{(1)}&a_{3}^{(1)}&b_{4}^{(1)}&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&b_{n-1}^{(1)}&a_{n-1}^{(1)}&b_{n}^{(1)} \\0&\ldots&\ldots&0&b_{n}^{(1)}&a_{n}^{(1)} \\\end{bmatrix}\), tolerancia TOL, número máximo de iteraciones Nmax.
• Set \(k=1\).

• While \(k\leq Nmax\)
  • If \(\max_{i=2,\ldots,n}|b_i^{(k)}|\leq TOL\)
    • Print valores propios \(a_1^{(k)},\ldots,a_n^{(k)}\)
    • STOP
  • Else
    • Descomponemos \(\mathbf{A}^{(k)}=\mathbf{Q}^{(k)}\mathbf{S}^{(k)}\).
    • Set \(\mathbf{A}^{(k+1)}=\mathbf{S}^{(k)}\mathbf{Q}^{(k)}\). (definimos \(\mathbf{A}^{(k+1)}\))
    • Set \(k=k+1\) (aumentamos \(k\))
• Print: el algoritmo no converge en \(Nmax\) pasos.

Sea \(\mathbf{A}\) una matriz Hessenberg superior:

\[ \mathbf{A}=\begin{bmatrix}a_{11}&a_{12}&a_{13}&a_{14}&\ldots&a_{1n} \\a_{21}&a_{22}&a_{23}&a_{24}&\ldots&a_{2n} \\0&a_{32}&a_{33}&a_{34}&\ldots&a_{3n} \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&a_{n-1,n-2}&a_{n-1,n-1}&a_{n-1,n} \\0&\ldots&\ldots&0&a_{n,n-1}&a_{nn} \\\end{bmatrix}. \] Vamos a aplicar el algoritmo anterior a este tipo de matrices.

Es decir, vamos a construir una sucesión de matrices \(\mathbf{A}^{(1)}=\mathbf{A},\mathbf{A}^{(2)},\ldots,\mathbf{A}^{(n)},\ldots,\) todas semejantes a la matriz original \(\mathbf{A}\) y Hesenberg superiores tal que los elementos de los matrices \(\mathbf{A}^{(n)}\) por debajo de la diagonal tienen en valor absoluto valores cada vez más pequeños a medida que aumenta \(n\). Es decir, \[ \lim_{n\to\infty}\max_{i> j}|a_{ij}^{(n)}|=0. \] Entonces los valores propios de la matriz original \(\mathbf{A}\) serán aproximadamente los valores de la diagonal de las matrices \(\mathbf{A}^{(n)}\).

Para hallar la sucesión de matrices \((\mathbf{A}^{(n)})_{n\geq 0}\) se siguen los pasos siguientes tal como hicimos en el caso de matrices tridiagonales simétricas:

• Paso \(1\). Sea \(\mathbf{A}^{(1)}=\mathbf{A}\). Descomponemos la matriz \(\mathbf{A}^{(1)}\) en el producto de una matriz ortogonal Hessenberg superior \(\mathbf{Q}^{(1)}\) y una matriz triangular superior \(\mathbf{S}^{(1)}\): \(\mathbf{A}^{(1)}=\mathbf{Q}^{(1)}\mathbf{S}^{(1)}\).

• Paso \(2\). Definimos \(\mathbf{A}^{(2)}=\mathbf{S}^{(1)}\mathbf{Q}^{(1)}\), el producto de la matriz triangular superior por la matriz ortogonal Hessenberg superior. Como \(\mathbf{A}^{(2)}=(\mathbf{Q}^{(1)})^\top\mathbf{A}^{(1)}\mathbf{Q}^{(1)}\) (ver el caso de matrices tridiagonales simétricas), las matrices \(\mathbf{A}^{(1)}\) y \(\mathbf{A}^{(2)}\) son semejantes y tienen los mismos valores propios. La matriz \(\mathbf{A}^{(2)}\) será Hessenberg superior ya que como \(\mathbf{A}^{(2)}\) es el producto de una matriz triangular superior por una matriz Hessenberg superior, \(\mathbf{A}^{(2)}=\mathbf{S}^{(1)}\mathbf{Q}^{(1)}\), la matriz \(\mathbf{A}^{(2)}\) será Hessenberg superior.

• A continuación descomponemos la matriz \(\mathbf{A}^{(2)}\) en el producto de una matriz ortogonal Hessenberg superior \(\mathbf{Q}^{(2)}\) y una matriz triangular superior \(\mathbf{S}^{(2)}\): \(\mathbf{A}^{(2)}=\mathbf{Q}^{(2)}\mathbf{S}^{(2)}\).

• Paso \(3\). Definimos \(\mathbf{A}^{(3)}=\mathbf{S}^{(2)}\mathbf{Q}^{(2)}\), etc.

• \(\vdots\)

• Paso \(i\). Definimos \(\mathbf{A}^{(i)}=\mathbf{S}^{(i-1)}\mathbf{Q}^{(i-1)}\). Usando un razonamiento inductivo y usando lo visto en el paso \(2\), la matriz \(\mathbf{A}^{(i)}\) será Hessenberg superior. Descomponemos la matriz \(\mathbf{A}^{(i)}\) en el producto de una matriz ortogonal Hessenberg superior \(\mathbf{Q}^{(i)}\) y una matriz triangular superior \(\mathbf{S}^{(i)}\): \(\mathbf{A}^{(i)}=\mathbf{Q}^{(i)}\mathbf{S}^{(i)}\).
• Paso \(i+1\). Definimos \(\mathbf{A}^{(i+1)}=\mathbf{S}^{(i)}\mathbf{Q}^{(i)}\). Como \(\mathbf{A}^{(i+1)}=(\mathbf{Q}^{(i)})^\top\mathbf{A}^{(i)}\mathbf{Q}^{(i)}\) (ver el caso de matrices tridiagonales simétricas), , las matrices \(\mathbf{A}^{(i+1)}\) y \(\mathbf{A}^{(i)}\) son semejantes y tienen los mismos valores propios. Además, usando el mismo razonamiento anterior, la matriz \(\mathbf{A}^{(i+1)}\) sería Hessenberg superior, etc.

La descomposición de una matriz Hesenberg superior \(\mathbf{A}\) en una matriz ortogonal \(\mathbf{Q}\) y una matriz triangular superior \(\mathbf{S}\), \(\mathbf{A}=\mathbf{Q}\mathbf{S}\) se hace siguiendo los mismos pasos que en el caso de matrices tridiagonales simétricas.

Los pseudocódigos para la descomposición anterior y para el cálculos de la sucesión \((\mathbf{A}^{(n)})_{n\geq 0}\) son los mismos que los pseudocódigos de las matrices tridiagonales simétricas.

Recapitulemos lo visto hasta el momento:

Si la matriz a hallar los valores propios es tridiagonal simétrica o Hessenberg superior, podemos hallar una sucesión de matrices \((\mathbf{A}^{(n)})_{n\geq 0}\) de matrices tridiagonales simétricas en el primer caso o Hessenberg superior en el segundo tal que los valores propios de la matriz original \(\mathbf{A}\) se aproximan como los valores diagonales de dicha sucesión de matrices a medida que \(n\) se hace grande.

Ahora bien, ¿qué ocurre si nuestra matriz \(\mathbf{A}\) no es de ninguno de los dos tipos mencionados antes?

En este caso, como ya comentamos en la introducción del capítulo, vamos a hallar una matriz ortogonal \(\mathbf{H}\) tal que \(\mathbf{H}^\top\mathbf{A}\mathbf{H}=\mathbf{B}\), donde \(\mathbf{B}\) será tridiagonal simétrica si la matriz \(\mathbf{A}\) es simétrica, o Hessenberg superior en caso contrario.

Como las matrices \(\mathbf{A}\) y \(\mathbf{B}\) tienen los mismos valores propios, para hallar los valores propios de \(\mathbf{B}\) usamos los métodos descritos anteriormente y ya tendremos una aproximación de los valores propios de la matriz original \(\mathbf{A}\).

La transformación de la matriz \(\mathbf{A}\) en la matriz \(\mathbf{B}\) se denomina método de Householder y la matriz de transformación \(\mathbf{H}\), matriz de Householder.

Empecemos, pues, a definir el concepto de matriz de Householder.

Definición de matriz de Householder.

Sea \(\mathbf{w}\in\mathbb{R}^n\) un vector de \(\mathbb{R}^n\) de norma euclídea igual a \(1\), es decir, \(\mathbf{w}^\top \mathbf{w}=1\).

La matriz \(n\times n\) \[ \mathbf{P}(\mathbf{w})=\mathbf{I}-2\mathbf{w}\mathbf{w}^\top, \] se denomina matriz de Householder.

Las matrices de Householder son simétricas y ortogonales:

Teorema. Sea \(\mathbf{w}\in\mathbb{R}^n\) un vector de \(\mathbb{R}^n\) con \(\|\mathbf{w}\|_2=1\), o \(\mathbf{w}^\top\mathbf{w}=1\). Entonces la matriz \(\mathbf{P}(\mathbf{w})\) es simétrica y ortogonal.

Como ya hemos comentado, el objetivo es transformar nuestra matriz original \(\mathbf{A}\) en una matriz Hessenberg superior \(\mathbf{B}\).

Realizaremos dicha transformación en \(n-1\) pasos. Sea \(\mathbf{A}^{(1)}=\mathbf{A}\).

• Paso \(1\). Hallaremos una matriz de Householder \(\mathbf{P}^{(1)}\) tal que la matriz \(\mathbf{A}^{(2)}=(\mathbf{P}^{(1)})^\top\mathbf{A}\mathbf{P}^{(1)}=\mathbf{P}^{(1)}\mathbf{A}\mathbf{P}^{(1)}\) sea una matriz con la primera columna nula a partir de la tercera componente: \(a_{i1}^{(2)}=0\), \(i=3,\ldots,n\):

\[ \mathbf{A}^{(2)}=\begin{bmatrix}a_{11}^{(2)}&a_{12}^{(2)}&a_{13}^{(2)}&\ldots&a_{1n}^{(2)} \\a_{21}^{(2)}&a_{22}^{(2)}&a_{23}^{(2)}&\ldots&a_{2n}^{(2)} \\0&a_{32}^{(2)}&a_{33}^{(2)}&\ldots&a_{3n}^{(2)} \\\vdots&\vdots&\vdots&\ddots&\vdots \\0&a_{n2}^{(2)}&\ldots&a_{n,n-1}^{(2)}&a_{nn}^{(2)} \\\end{bmatrix}. \]

• Paso \(1\) (continuación). Las matrices \(\mathbf{A}^{(2)}\) y \(\mathbf{A}^{(1)}\) son semejantes con matriz de transformación ortogonal y simétrica \(\mathbf{P}^{(1)}\). Por tanto, tienen los mismos valores propios. Además, si la matriz \(\mathbf{A}=\mathbf{A}^{(1)}\) es simétrica, también lo es \(\mathbf{A}^{(2)}\) ya que \((\mathbf{A}^{(2)})^\top =(\mathbf{P}^{(1)}\mathbf{A}\mathbf{P}^{(1)})^\top = \mathbf{P}^{(1)}\mathbf{A}^\top\mathbf{P}^{(1)}=\mathbf{P}^{(1)}\mathbf{A}\mathbf{P}^{(1)}=\mathbf{A}^{(2)}\).

• Paso \(2\). Hallaremos una matriz de Householder \(\mathbf{P}^{(2)}\) tal que la matriz \(\mathbf{A}^{(3)}=(\mathbf{P}^{(2)})^\top\mathbf{A}^{(2)}\mathbf{P}^{(2)}=\mathbf{P}^{(2)}\mathbf{A}^{(2)}\mathbf{P}^{(2)}\) sea una matriz de la forma siguiente, con la primera columna nula a partir de la tercera componente y con la segunda columna nula a partir de la cuarta componente: \(a_{i1}^{(3)}=0\), \(i=3,\ldots,n\), \(a_{i2}^{(3)}=0\), \(i=4,\ldots,n\).

\[ \mathbf{A}^{(3)}=\begin{bmatrix}a_{11}^{(3)}&a_{12}^{(3)}&a_{13}^{(3)}&\ldots&a_{1n}^{(3)} \\a_{21}^{(3)}&a_{22}^{(3)}&a_{23}^{(3)}&\ldots&a_{2n}^{(3)} \\0&a_{32}^{(3)}&a_{33}^{(3)}&\ldots&a_{3n}^{(3)} \\0&0&a_{43}^{(3)}&\ldots&a_{4n}^{(3)} \\\vdots&\vdots&\vdots&\ddots&\vdots \\0&0&a_{n3}^{(3)}&\ldots&a_{nn}^{(3)} \\\end{bmatrix}. \]

• Paso \(2\) (continuación). Las matrices \(\mathbf{A}^{(3)}\) y \(\mathbf{A}^{(2)}\) son semejantes con matriz de transformación ortogonal y simétrica \(\mathbf{P}^{(2)}\). Por tanto, tienen los mismos valores propios. Además, si la matriz \(\mathbf{A}=\mathbf{A}^{(2)}\) es simétrica, también lo es \(\mathbf{A}^{(3)}\) ya que \((\mathbf{A}^{(3)})^\top =(\mathbf{P}^{(2)}\mathbf{A}^{(2)}\mathbf{P}^{(2)})^\top = \mathbf{P}^{(2)}(\mathbf{A}^{(2)})^\top\mathbf{P}^{(2)}=\mathbf{P}^{(2)}\mathbf{A}^{(2)}\mathbf{P}^{(2)}=\mathbf{A}^{(3)}\).

• Paso \(k\)-ésimo. Suponemos que la matriz \(\mathbf{A}^{(k)}\) es de la forma:

\[ \mathbf{A}^{(i)}=\begin{bmatrix}a_{11}^{(k)}&a_{12}^{(k)}&\ldots&a_{1,k-1}^{(k)}&a_{1,k}^{(k)}&\ldots&a_{1,n}^{(k)} \\a_{21}^{(k)}&a_{22}^{(k)}&\ldots&a_{2,k-1}^{(k)}&a_{2,k}^{(k)}&\ldots&a_{2,n}^{(k)} \\0&a_{32}^{(k)}&\ldots&a_{3,k-1}^{(k)}&a_{3,k}^{(k)}&\ldots&a_{3,n}^{(k)} \\\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots \\0&0&\ldots&a_{k-1,k-1}^{(k)}&a_{k-1,k}^{(k)}&\ldots&a_{k-1,n}^{(k)} \\0&0&\ldots&a_{k,k-1}^{(k)}&a_{k,k}^{(k)}&\ldots&a_{k,n}^{(k)} \\0&0&\ldots&0&a_{k+1,k}^{(k)}&\ldots&a_{k+1,n}^{(k)} \\\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots \\0&0&\ldots&0&a_{n,k}^{(k)}&\ldots&a_{n,n}^{(k)} \\\end{bmatrix}. \]

• Paso \(k\)-ésimo. Es decir, \(a_{jl}^{(k)}=0\) para \(l=1,\ldots,k-1\), \(j=l+2,\ldots,n\). Dicho en otras palabras, la matriz \(\mathbf{A}^{(k)}\) es una matriz Hessenberg superior hasta la columna \(k-1\)-ésima.

Hallaremos una matriz de Householder \(\mathbf{P}^{(k)}\) tal que la matriz \(\mathbf{A}^{(k+1)}=(\mathbf{P}^{(k)})^\top\mathbf{A}^{(k)}\mathbf{P}^{(k)}=\mathbf{P}^{(k)}\mathbf{A}^{(k)}\mathbf{P}^{(k)}\) sea una matriz de Householder hasta la columna \(k\)-ésima:

\[ \mathbf{A}^{(k+1)}=\begin{bmatrix}a_{11}^{(k+1)}&a_{12}^{(k+1)}&\ldots&a_{1,k}^{(k+1)}&a_{1,k+1}^{(k+1)}&\ldots&a_{1,n}^{(k+1)} \\a_{21}^{(k+1)}&a_{22}^{(k+1)}&\ldots&a_{2,k}^{(k+1)}&a_{2,k+1}^{(k+1)}&\ldots&a_{2,n}^{(k+1)} \\0&a_{32}^{(k+1)}&\ldots&a_{3,k}^{(k+1)}&a_{3,k+1}^{(k+1)}&\ldots&a_{3,n}^{(k+1)} \\\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots \\0&0&\ldots&a_{k,k}^{(k+1)}&a_{k,k+1}^{(k+1)}&\ldots&a_{k,n}^{(k+1)} \\0&0&\ldots&a_{k+1,k}^{(k+1)}&a_{k+1,k+1}^{(k+1)}&\ldots&a_{k+1,n}^{(k+1)} \\0&0&\ldots&0&a_{k+2,k+1}^{(k+1)}&\ldots&a_{k+2,n}^{(k+1)} \\\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots \\0&0&\ldots&0&a_{n,k+1}^{(k+1)}&\ldots&a_{n,n}^{(k+1)} \\\end{bmatrix}. \]

• Paso \(n-1\)-ésimo. La matriz \(\mathbf{B}\) será \(\mathbf{B}=\mathbf{A}^{(n-1)}\) ya que sería Hessenberg superior en el caso general y tridiagonal simétrica en el caso en que la matriz original sea simétrica. La última afirmación se deduce del hecho de que las matrices que se van obteniendo \(\mathbf{A}^{(1)},\mathbf{A}^{(2)},\ldots,\mathbf{A}^{(n-1)}\) son Hessenberg superiores y simétricas y toda matriz Hessenberg superior simétrica es trivialmente tridiagonal simétrica.

De hecho podemos escribir la matriz \(\mathbf{A}^{(n-1)}\) de la forma siguiente: \[ \mathbf{B}=\mathbf{A}^{(n-1)}=\mathbf{P}^{(n-1)}\cdots\mathbf{P}^{(1)}\mathbf{A}\mathbf{P}^{(1)}\cdots\mathbf{P}^{(n-1)}. \]

Veamos cómo calcular las matrices de transformación o de Householder en los pasos anteriores:

• Cálculo de \(\mathbf{P}^{(1)}\). Hemos de hallar un vector \(\mathbf{w}\) tal que \(\mathbf{P}^{(1)}=\mathbf{I}-2\mathbf{w}(\mathbf{w})^\top\). La condición que debe verificar \(\mathbf{P}^{(1)}\) es la siguiente: \[ \mathbf{P}^{(1)}\mathbf{A}^{(1)}\mathbf{P}^{(1)}=\mathbf{A}^{(2)}, \] con \(a_{i1}^{(2)}=0\), para \(i=3,\ldots,n\).

Elegimos el vector \(\mathbf{w}\) de la forma \(\mathbf{w}=(0,\hat{\mathbf{w}})^\top\), es decir, elegimos el vector \(\mathbf{w}\) con la primera componente igual a \(0\) y \(\hat{\mathbf{w}}=\begin{bmatrix}w_2\\\vdots\\ w_n\end{bmatrix}\).

Además, imponemos que se verifique que: \(\mathbf{P}^{(1)}\mathbf{a}^{(1)}=\begin{bmatrix}a_{11}^{(1)}\\ \alpha\\ 0\\\vdots\\ 0\end{bmatrix}\), donde \(\mathbf{a}^{(1)}\) representa la primera columna de la matriz \(\mathbf{A}^{(1)}\) y \(\alpha\) es un valor que calcularemos más tarde.

Escribimos la matriz a hallar \(\mathbf{P}^{(1)}\) de la forma siguiente:

\[ \mathbf{P}^{(1)}=\begin{bmatrix}1&0&\ldots&0 \\0&*&\ldots&* \\\vdots&\vdots&\vdots&\vdots \\0&*&\ldots&* \\\end{bmatrix}. \] Llamamos \(\hat{\mathbf{P}}\) a la submatriz anterior de dimensiones \((n-1)\times (n-1)\) donde aparecen los asteriscos: \(\mathbf{P}^{(1)}=\begin{bmatrix}1&\mathbf{0}\\\mathbf{0}&\hat{\mathbf{P}}\end{bmatrix}\).

Como \(\mathbf{P}^{(1)}=\mathbf{I}-2\mathbf{w}\mathbf{w}^\top\), la matriz \(\hat{\mathbf{P}}\) también será de Householder, (y, por tanto, simétrica y ortogonal) con vector asociado \(\hat{\mathbf{w}}\): \(\hat{\mathbf{P}}=\mathbf{I}-2\hat{\mathbf{w}}\hat{\mathbf{w}}^\top\).

Como \(\mathbf{P}^{(1)}\mathbf{a}^{(1)}=\begin{bmatrix}a_{11}^{(1)}\\ \alpha\\ 0\\\vdots\\ 0\end{bmatrix}\), tendremos que: \[ \begin{bmatrix}1&\mathbf{0}\\\mathbf{0}&\hat{\mathbf{P}}\end{bmatrix}\begin{bmatrix}a_{11}^{(1)}\\ a_{21}^{(1)}\\\vdots\\ a_{n1}^{(1)}\end{bmatrix}=\begin{bmatrix}a_{11}^{(1)}\\ \alpha\\ 0\\\vdots\\ 0\end{bmatrix},\ \Rightarrow \hat{\mathbf{P}}\begin{bmatrix}a_{21}^{(1)}\\ \\\vdots\\ a_{n1}^{(1)}\end{bmatrix}=\begin{bmatrix}\alpha\\ 0\\\vdots\\ 0\end{bmatrix}. \]

Llamamos \(\hat{\mathbf{a}^{(1)}}\) a la subcolumna de la matriz \(\mathbf{A}^{(1)}\) formada por las \(n-1\) últimas componentes. La condición anterior será, pues: \[ \hat{\mathbf{P}}\hat{\mathbf{a}^{(1)}}=\begin{bmatrix}\alpha\\ 0\\\vdots\\ 0\end{bmatrix},\ \Rightarrow (\mathbf{I}-2\hat{\mathbf{w}}\hat{\mathbf{w}}^\top)\hat{\mathbf{a}^{(1)}}=\begin{bmatrix}\alpha\\ 0\\\vdots\\ 0\end{bmatrix},\ \Rightarrow \hat{\mathbf{a}^{(1)}}-2((\hat{\mathbf{w}})^\top \hat{\mathbf{a}^{(1)}})\hat{\mathbf{w}}=\begin{bmatrix}\alpha\\ 0\\\vdots\\ 0\end{bmatrix}. \] Sea \(\beta =(\hat{\mathbf{w}})^\top \hat{\mathbf{a}^{(1)}}\). La condición anterior, escrita en componentes, es la siguiente: \[ a_{21}^{(1)}-2\beta w_2=\alpha,\ a_{j1}^{(1)}-2\beta w_j=0,\ j=3,\ldots,n. \]

Entonces: \[ 2\beta w_2=a_{21}^{(1)}-\alpha,\ 2\beta w_j=a_{j1}^{(1)},\ j=3,\ldots,n. \] Como \(\hat{\mathbf{w}}^\top \hat{\mathbf{w}}=1\), tendremos que \(\displaystyle\sum_{j=2}^n w_j^2=1\). Entonces: \[ 4\beta^2 \left(w_2^2+\sum_{j=3}^n w_j^2\right)=4\beta^2\sum_{j=2}^n w_j^2=4\beta^2 =(a_{21}^{(1)}-\alpha)^2+\sum_{j=3}^n (a_{j1}^{(1)})^2. \] Ahora bien, recordemos que \(\hat{\mathbf{P}}\hat{\mathbf{a}^{(1)}}=\begin{bmatrix}\alpha\\ 0\\\vdots\\ 0\end{bmatrix}\), es decir \(\hat{\mathbf{a}^{(1)}}=\hat{\mathbf{P}}\begin{bmatrix}\alpha\\ 0\\\vdots\\ 0\end{bmatrix}.\)

Como \(\hat{\mathbf{P}}\) es ortogonal: \[ \begin{align*} \|\hat{\mathbf{a}^{(1)}}\|_2^2 & =(\hat{\mathbf{a}^{(1)}})^\top \hat{\mathbf{a}^{(1)}}=\left(\hat{\mathbf{P}}\begin{bmatrix}\alpha\\ 0\\\vdots\\ 0\end{bmatrix}\right)^\top \hat{\mathbf{P}}\begin{bmatrix}\alpha\\ 0\\\vdots\\ 0\end{bmatrix}=(\alpha,0,\ldots,0)\hat{\mathbf{P}}^\top\hat{\mathbf{P}}\begin{bmatrix}\alpha\\ 0\\\vdots\\ 0\end{bmatrix}\\ & =(\alpha,0,\ldots,0)\begin{bmatrix}\alpha\\ 0\\\vdots\\ 0\end{bmatrix}=\alpha^2. \end{align*} \]

El valor de \(\alpha^2\) será: \[ \alpha^2 = \|\hat{\mathbf{a}^{(1)}}\|_2^2=\sum_{j=2}^n (a_{j1}^{(1)})^2. \] El valor de \(\beta^2\) verificará: \[ \begin{align*} 4\beta^2 & =(a_{21}^{(1)}-\alpha)^2+\sum_{j=3}^n (a_{j1}^{(1)})^2=\alpha^2-2\alpha a_{21}^{(1)}+(a_{21}^{(1)})^2+\sum_{j=3}^n (a_{j1}^{(1)})^2\\ & =\sum_{j=2}^n (a_{j1}^{(1)})^2-2\alpha a_{21}^{(1)}+\sum_{j=2}^n (a_{j1}^{(1)})^2 =2\sum_{j=2}^n (a_{j1}^{(1)})^2-2\alpha a_{21}^{(1)}, \end{align*} \]

de donde deducimos que: \[ \beta^2 =\frac{1}{2}\left(\sum_{j=2}^n (a_{j1}^{(1)})^2-\alpha a_{21}^{(1)}\right). \] Como \(\beta^2\geq 0\), tenemos que elegir \(\alpha\) como: \[ \alpha =-\mathrm{signo}(a_{21}^{(1)})\sqrt{\sum_{j=2}^n (a_{j1}^{(1)})^2}, \] donde la función signo está definida como \[ \mathrm{signo}(x)=\begin{cases}1, & \mathrm{si\ }x\geq 0,\\ -1,& \mathrm{si\ }x<0.\end{cases} \]

El valor de \(\beta^2\), será, pues: \[ \beta^2 = \frac{1}{2}\left(\sum_{j=2}^n (a_{j1}^{(1)})^2+ |a_{21}^{(1)}|\sqrt{\sum_{j=2}^n (a_{j1}^{(1)})^2}\right). \]

Las componentes \(w_i\), \(i=2,\ldots,n\) del vector \(\mathbf{w}\) (recordemos que elegimos \(w_1=0\)) son: \[ w_2 = \frac{a_{21}^{(1)}-\alpha}{2\beta},\ w_j =\frac{a_{j1}^{(1)}}{2\beta},\ j=3,\ldots,n. \]

Ya tenemos definida la matriz \(\hat{\mathbf{P}}\) y la matriz \(\mathbf{P}\).

Veamos que es la matriz que buscamos:

Proposición.

La matriz \(\mathbf{P}^{(1)}=\begin{bmatrix}1&\mathbf{0}\\\mathbf{0}&\hat{\mathbf{P}}\end{bmatrix}\), con \(\hat{\mathbf{P}}\) la matriz hallada anteriormente verifica que: \[ \mathbf{P}^{(1)}\mathbf{A}^{(1)}\mathbf{P}^{(1)}=\begin{bmatrix}a_{11}^{(2)}&a_{12}^{(2)}&a_{13}^{(2)}&\ldots&a_{1n}^{(2)} \\a_{21}^{(2)}&a_{22}^{(2)}&a_{23}^{(2)}&\ldots&a_{2n}^{(2)} \\0&a_{32}^{(2)}&a_{33}^{(2)}&\ldots&a_{3n}^{(2)} \\\vdots&\vdots&\vdots&\ddots&\vdots \\0&a_{n2}^{(2)}&\ldots&a_{n,n-1}^{(2)}&a_{nn}^{(2)} \\\end{bmatrix}. \]

• Paso \(k\)-ésimo. Elegimos el vector \(\mathbf{w}\) de la forma: \[ \mathbf{w}=\begin{bmatrix}\mathbf{0}\\\hat{\mathbf{w}}\end{bmatrix}=\begin{bmatrix}\mathbf{0}\\ w_{k+1}\\\vdots\\ w_n\end{bmatrix},\ \hat{\mathbf{w}}=\begin{bmatrix}w_{k+1}\\\vdots\\ w_n\end{bmatrix}, \] y la matriz \(\mathbf{P}^{(k)}\) de la forma: \(\mathbf{P}^{(k)}=\mathbf{I}-2\mathbf{w}\mathbf{w}^\top.\)

Los valores \(w_{k+1},\ldots,w_n\) son los siguientes: \[ w_{k+1}=\frac{a_{k+1,k}^{(k)}-\alpha}{2\beta},\ w_j=\frac{a_{jk}^{(k)}}{2\beta},\ j=k+2,\ldots,n, \] con \[ \alpha = -\mathrm{signo}(a_{k+1,k}^{(k)})\sqrt{\sum_{j=k+1}^n (a_{jk}^{(k)})^2},\ \beta =\sqrt{\frac{1}{2}\left(\alpha^2-\alpha a_{k+1,k}^{(k)}\right)}. \]

• INPUT matriz del sistema \(\mathbf{A}=(a_{ij})_{i=1,\ldots,n,j=1,\ldots,n}\)
• Set \(\mathbf{A}^{(1)}=\mathbf{A}\).
• For \(k=1,2,\ldots,n-1\)
  • Set \(\displaystyle s=\sum_{j=k+1}^n (a_{jk}^{(k)})^2\).
  • If \(a_{k+1,k}^{(k)}=0\)
    • Set \(\alpha = -\sqrt{s}\).
  • Else
    • Set \(\alpha =-\frac{\sqrt{s}a_{k+1,k}^{(k)}}{|a_{k+1,k}^{(k)}|}\).

• • Set \(\beta = \sqrt{\frac{1}{2}\left(\alpha^2-\alpha a_{k+1,k}^{(k)}\right)}\).
  • Set \(\hat{w}_{1}=\frac{a_{k+1,k}^{(k)}-\alpha}{2\beta}\). (definimos el vector \(\hat{w}\))
  • For \(j=k+2,\ldots,n\)
    • Set \(\hat{w}_{j-k}=\frac{a_{jk}^{(k)}}{2\beta}\).
  • Set \(\hat{\mathbf{P}}=\mathbf{I}-2\hat{w}\hat{w}^\top\) (definimos la matriz \(\hat{\mathbf{P}}\))

• • Set \(\mathbf{P}^{(k)}=\mathbf{I}_n\) (definimos inicialmente la matriz \(\mathbf{P}^{(k)}\) como la matriz identidad)
  • For \(i=k+1,\ldots,n\)
    • For \(j=k+1,\ldots,n\)
      • Set \(p_{ij}^{(k)}=\hat{p}_{i-k,j-k}\).
  • Set \(\mathbf{A}^{(k+1)}=\mathbf{P}^{(k)}\mathbf{A}^{(k)}\mathbf{P}^{(k)}\) (definimos la matriz \(\mathbf{A}^{(k+1)}\))
• Print \(\mathbf{A}^{(n-1)}\). (damos la matriz tridiagonal simétrica si la matriz original \(\mathbf{A}\) es simétrica o Hessenberg superior en caso contrario.)