Aproximación de valores y vectores propios

Cálculo de los discos de Gerschgorin

En primer lugar vamos a considerar una matriz cuadrada de orden n, con coeficientes reales, y calcular los radios de los discos de Gerschgorin, tanto por filas como por columnas, para dibujar los discos correspondientes.

\(A = {{\(-9.\), 1., 2., 1.}, {1, 2, \(-1\), 0}, {1, 1, 0, 1}, {1, 2, 1, 
          \(-3\)}}; 

    n = \(Dimensions[A]\)[[1]]; 

    rf = Table[
        ∑\+\(j = 1\)\%n Abs[
              A[[i, j]]] - 
          Abs[A[[i, i]]], {i, n}]\)rc = 
      Table[(∑_(i = 1)^n) Abs[
              A[[i, j]]] - 
          Abs[A[[j, j]]], {j, n}]

gf = 
      Table[Graphics[
          Circle[{A[[i, i]], 0}, 
            rf[[i]]], 
          
              
          AspectRatio -> Automatic], {i, 4}];

Show[gf, Axes -> Automatic];

gc = 
      Table[Graphics[
          Circle[{A[[i, i]], 0}, 
            rc[[i]]], 
          
              
          AspectRatio -> Automatic], {i, 4}];

Show[gc, Axes -> Automatic];

En la primera gráfica tenemos los discos de Gerschgorin cuyos radios se han obtenido sumando los valores absolutos de los elementos de A, excluido el de la diagonal, por filas. Para la segunda, los radios se han obtenido sumando por columnas. De ambas gráficas deducimos que existe un valor propio dominante cuyo valor es próximo a -9.

Cálculo del valor propio dominante

El método de las potencias

En muchos problemas de valores propios no es necesario calcular todos todos los valores propios, sino sólo el dominante, es decir el valor propio mayor en valor absoluto. Supuesto que exista un valor propio dominante y que sea posible tener una base de (ℝ^n) formada por vectores propios, existe un método iterativo, llamado método de las potencias, que permite aproximarlo; si fuera necesario aproximar todos los valores propios de una matriz, el método de las potencias se puede combinar con otro, denominado método de deflación (consistente en aplicar el método de las potencias a distintas matrices), para calcular los restantes valores propios.

Como todo método iterativo, parte de una aproximación inicial y obtiene una nueva, que reemplaza a la inicial para repetir el proceso de forma sucesiva. La aproximación inicial no puede ser en este caso un vector nulo.

Para la matriz dada anteriormente, consideremos el vector x, inicial, cuyas componentes son todas iguales a 1.

x = Table[1., {i, n}];

A continuación obtenemos 15 aproximaciones del valor propio dominante, que se almacenan en una lista, aprox:

\(aprox = {}; 

    For[i = 1, i ≤ 15, \(i++\), 
\ \ \ \ \ \ \ 
      y = A.x; 
\ \ \ \ \ \ \ 
      AppendTo[aprox, 
        y[[1]]\/x[[1
            ]]]; 
\ \ \ \ \ \ \ x = y]\)aprox

En este ejemplo, se podía haber usado como aproximantes los cocientes entre las segundas componentes de los sucesivos vectores, o entre las terceras. Cámbielas y procese de nuevo el vector inicial y las aproximaciones del valor propio dominante.
Recordemos que Mathematica, cuando la matriz A no es de orden muy elevado, puede calcular los valores propios directamente con la sentencia Eigenvalues y los vectores propios con Eigenvectors. Por tanto, en este caso podemos comprobar qué tal es la aproximación del valor propio dominante obtenida anteriormente.

N[Eigenvalues[A]]

En el cociente entre las componentes elegidas puede aparecer una división por cero; por tanto, debemos asegurar que esto no va a ocurrir. Para ello lo mejor es usar el cociente en el que el denominador tenga valor máximo.  Otro problema es el crecimiento hacia infinito (cuando el valor propio dominante es mayor que uno en valor absoluto) o el decrecimiento hacia cero (cuando el valor propio dominante es menor que uno) de las sucesivas iteraciones. Para observar este hecho vamos a repetir los cálculos usados en el ejemplo anterior, sólo que escribiendo los vectores x en lugar de los cocientes.

Volvemos a asignar a x el vector inicial considerado anteriormente.

x = Table[1., {i, n}];

Generamos diez iteraciones.

\(aprox = {}; 

    For[i = 1, i ≤ 10, \(i++\), 
\ \ \ \ \ \ \ 
      y = A.x; 
\ \ \ \ \ \ \ AppendTo[aprox, y]; 
      
\ \ \ \ \ \ \ x = y]\)TableForm[aprox, TableAlignments -> Left]

Para evitar que las componentes de los vectores x crezcan hacia infinito o decrezcan hacia cero debemos normalizar el vector x en cada iteración, es decir multiplicarlos por una cantidad de forma que sus componentes tengan valores similares en cada iteración. Por ejemplo, si el vector x se divide por la componente mayor, ésta pasa a valer 1 y las demás son menores que 1 (en valor absoluto); las componentes del siguiente vector y quedan divididas por la misma cantidad. Por tanto el cociente no varía. Por otra parte si usamos para el cociente entre componentes aquélla que en x vale 1, no es necesario efectuar la división y además evitamos el problema que habíamos comentado sobre una posible división por cero.

Inicializamos x, una vez más.

x = Table[1., {i, n}];

Buscamos la componente de valor absoluto máximo del vector y. Para ello usamos las órdenes Max[lista], que proporciona el máximo de la lista especificada, y Position[lista, valor], que da como resultado la posición de valor en la lista. Este resultado será una lista, pues el máximo puede alcanzarse en diferentes posiciones. Por tanto, nos quedaremos con la componente de subíndices (1, 1), es decir, con la primera coincidencia.

\(y = A.x; 

    For[i = 1, i ≤ 15, \(i++\), 
\ \ \ \ \ \ \ 
      m = Max[Abs[y]]; 
\ \ \ \ \ \ \ \ 
      c = \(componentemaxima = 
          \(Position[Abs[y], m]\)[[1, 1
            ]]\); 
\ \ \ \ \ \ \ \ 
      x = y\/y[[c]]; 
      
\ \ \ \ \ \ \ \ y = A.x]\)\(Print["Valor Propio dominante:  ", 
      y[[c]]]; 

    Print["Vector propio asociado:", y]; \)

Los vectores propios pueden obtenerse con la sentencia Eigenvectors. Los da todos, siendo el primero de ellos el correspondiente al valor propio dominante. Recuerde que, dado un vector propio, cualquier otro vector proporcional a él es también vector propio asociado al mismo valor propio. 
El obtenido anteriormente es unas diez veces mayor que el dado por la sentencia Eigenvectors (y de signo opuesto).

Eigenvectors[A]

Aunque el método de las potencias lo hemos usado para obtener el valor propio dominante, (cuando existe), puede extenderse al cálculo de los demás valores propios, realizando una transformación adecuada en la matriz A una vez que se ha calculado el valor propio dominante y repitiendo el proceso. Posteriormente trataremos esta situación.

Cocientes de Rayleigh

El método de los cocientes de Rayleigh es otro procedimiento para calcular el valor propio dominante de una matriz. En las mismas condiciones que el método de las potencias, si la matriz a la que calculamos una aproximación del valor propio dominante es simétrica, el método de los cocientes de Rayleigh produce una sucesión convergente a éste con mayor rapidez que el cociente entre componenetes. 

Vamos a aplicarlo a un caso concreto.

\(M = {{4., 1., 0., 0., 0., 0.}, {1., 4., 1., 0., 0., 0.}, 
        
\ \ {0., 1., 4., 1., 0., 0.}, {0., 0., 1., 4., 1., 
          0.}, 
\ \ {0., 0., 0., 1., 4., 1.}, {0., 0., 0., 
          0., 1., 4.}}; 
MatrixForm[M]\)

Partimos de la siguiente aproximación inicial:

x = {1., 1., 1., 1., 1., 1.};

Generemos 10 términos de la sucesión de aproximaciones del valor propio dominante y para la última imprimimos el vector propio asociado:

\(aprox = {}; 
y = M.x; 

    For[i = 1, i ≤ 10, \(i++\), 
\ \ \ \ \ \ \ 
      x = y; 
\ \ \ \ \ \ \ y = M.x; 
      
\ \ \ \ \ \ \ AppendTo[aprox, x.y\/x.x]]; 

    Print["Aproximaciones del valor propio dominante:"]; 
aprox\)Print["Aproximación del vector propio asociado: ", 
      (y/(Sqrt[y.y]))]

Comparemos la aproximación del valor propio dominante con la proporcionada por Eigenvalues.

Eigenvalues[M]

Con sólo 10 iteraciones hemos obtenido una aproximación que coincide con la proporcionada por la orden Eigenvalues hasta la tercera cifra decimal.

Puede comprobarse que el método de las potencias requiere 33 iteraciones para obtener igual número de cifras coincidentes.

El método de deflación de Wielandt

Conocido un valor propio de una matriz, el método de deflación consiste en construir otra matriz cuyos valores propios sean los mismos que los de la inicial excepto el considerado, que pasa a ser cero. Así, calculado el valor propio dominante, podemos pasar a aplicar de nuevo el método de las potencias para calcular el siguiente (en módulo) y el proceso puede continuar hasta calcularlos todos, supuesto que tengan módulo diferente.

En el proceso de cálculo se requiere hallar sucesivamente valores propios, por lo que precisamos una rutina, que llamaremos pot,  para aplicar el método de las potencias a las sucesivas matrices que aparecen. El módulo proporcionará como salida una lista con dos posiciones: la primera será la aproximación del valor propio dominante y la segunda su vector propio asociado. Las entradas son la matriz, el vector inicial y el número de iteraciones que se realizan para obtener la aproximación.

pot[matriz_, vectorinicial_, iterac_] := 
      
Module[{A = matriz, x = vectorinicial, y, m}, 
        
                y = A.x; 
        
                
        For[i = 1, i ≤ iterac, i++, 
          
                       
          x = y; 
                       m = Max[Abs[x]]; 
          
                       
          c = Position[Abs[x], m][[1, 1
              ]]; 
          
                       
          x = (x/x)[[c]]; 
          
                       
          y = A.x]; 
        
                {
          y[[c]], y}]

Para la matriz M y el vector x anteriormente definidos, y 30 iteraciones se considera la siguiente orden:

pot[M, x, 30]

El método de Wielandt es un caso particular del procedimiento general, y vamos a ilustrarlo haciendo todos los cálculos correspondientes a una etapa.

Volvemos a trabajar con la matriz A, que redefinimos, aunque sigue activa.

A = {{-9., 1., 2., 1.}, {1, 2, -1, 0}, {1, 1, 0, 1}, {1, 2, 
          1, -3}};

Le aplicamos el método de las potencias.

lista = pot[A, {1, 1, 1, -8}, 100]

Definimos variables para las aproximaciones del valor propio dominante y el vector propio asociado:

valp1 = lista[[1]]vecp1 = lista[[2]]

Definimos una variable para el índice correspondiente a la primera coordenada no nula de ese vector propio:

ind = 1;

Formamos el vector que nos permite modificar la matriz de partida:

v = 
      A[[ind]]/(valp1 
          vecp1[[ind]])

Construimos una nueva matriz, B, a partir de A.

B = 
      A - Table[{vecp1[[i]]}, {i, 
                Length[vecp1]}].{
              A[[ind]]}/
          vecp1[[ind]]

Esta matriz y la inicial tienen los mismos valores propios, excepto el dominante, que ha pasado a ser cero. La fila asociada al índice ind  (la primera, en este caso) debería ser nula, pero los errores de redondeo afectan al cálculo, aunque obtenemos valores muy pequeños.

Eigenvalues[A]Eigenvalues[B]

Aplicamos el método de las potencias a esta matriz B.

lista = pot[B, {1, 1, 1, -8}, 100]

La primera posición de la lista anterior es el mayor valor propio de B y segundo (en módulo) de A.

valp2 = lista[[1]]

La segunda posición contiene un vector propio asociado de B.

w = lista[[2]]

Finalizamos calculando un vector propio de A asociado al valor propio calculado, para lo que hacemos uso del vector propio de B que acabamos de mostrar.

vecp2 = (valp2 - valp1) w + valp1 v.w vecp1

Como

Eigenvectors[A][[2]]

y

-((vecp2 0.0581705)/1.22366)

el vector propio calculado es proporcional al segundo de los vectores dado por la orden Eigenvectors[A]. Análogamente podríamos continuar hasta obtener todos los valores y vectores propios.

Desplazamientos de una matriz

Si la matriz A es desplazada por la matriz c I, para obtener una matriz B = A - c I, entonces los valores propios de B son los valores propios de A restándoles c. Comprobémoslo usando la orden directa Eigenvalues para calcular los valores propios de A1 y de A1d = A1 - 8  I. A los de A1d hay que sumarle 8 para obtener los de A1.

\(A1 = Table[Random[], {i, 5}, {j, 5}]; 
MatrixForm[A1]\)A1d = A1 - 8 IdentityMatrix[5];

Eigenvalues[A1]Eigenvalues[A1d] + 8

Los vectores propios de ambas matrices son los mismos o proporcionales caso de ser reales. Si las componentes son complejas pueden ser conjudados unos respecto de otros.

Eigenvectors[A1]Eigenvectors[A1d]

Podemos desplazar los valores propios de la matriz M, dada más arriba, de forma que la dominancia entre del primero con respecto al segundo aumente. Como los valores propios de M son 5.80194, 5.24698, 4.44504, 3.55496, 2.75302 y 2.19806, el desplazamiento óptimo es restar a M la matriz identidad multiplicada por la semisuma entre el valor propio segundo y el último, es decir 3.72252, con lo cual los valores propios de la nueva matriz disminuyen en dicha cantidad y serán 2.07942, 1.52446, 0.72252, -0.16756, -0.9695 y -1.52446. Ahora el cociente entre el segundo valor propio y el primero tiene valor mínimo y, por tanto, la convergencia es más rápida.

Md = M - 3.72252 IdentityMatrix[6];

pot[Md, x, 15]

Observamos que en este caso con 15 iteraciones el método de las potencias da las 6 cifras coincidentes con las obtenidas con la orden directa Eigenvalues[Md], mientras que anteriormente con 30 iteraciones sólo conseguimos la mitad de cifras coincidentes

Eigenvalues[Md]

Evidentemente el desplazamiento óptimo, como el anterior, no podemos realizarlo en la práctica, pero desplazamientos que mejoren sensiblemente la convergencia es posible realizarlos en muchos problemas, ya que se puede tener alguna información sobre los valores propios. Por ejemplo, la obtenida aplicando el teorema de Gerschgorin.

Método de potencias inversas

Si A tiene inversa no puede tener por valor propio al cero y los valores propios de la inversa de A son los inversos de los valores propios de A. En tal caso, el método de las potencias aplicado a la inversa de A nos da el inverso del valor propio de A más cercano al cero (supuesto que sea menor en módulo que los restantes). Veámoslo aplicando el método de las potencias a la inversa de la matriz A dada a continuación.

\(A = {{17, 1, 3, 0}, {1, 2.5, 3.5, 4.5}, 
        
\ \ {0, 1, 0, 1}, {1, \(-1\), \(-1\), \(-2\)}}; 

    
pot[Inverse[A], {1., 1., 1., 1.}, 20]\)

Eigenvalues[A]

El inverso de -2.75016 es -0.363615 y el valor propio de menor módulo de la matriz A obtenido con la orden directa es -0.363616. 
Las iteraciones podíamos haberlas planteado resolviendo un sistema lineal en cada paso. Veámoslo

\(x0 = {1., 1, 1., 1}; 
x1 = LinearSolve[A, x0]; 

    x2 = LinearSolve[A, x1]; 
x3 = LinearSolve[A, x2]; 

    x4 = LinearSolve[A, x3]; 

    x4[[3]]\/x3[[3
        ]]\)

Podemos generar un bucle con vectores usando índices del tipo x[i], y de esa forma realizar el número de iteraciones que se desee.

Los vectores propios de una matriz y su inversa son los mismos o proporcionales

Eigenvectors[A]

Eigenvectors[Inverse[A]]

Ejercicios

1.- Se considera la matriz

		A=(({{8.1, 1., 0., 0.}, {1., -5., 1., 0.}, {0., 1., 0., 1.}, {0., 0., 1., -1.5}})).

a) Calcule cinco aproximaciones, partiendo de {1.,1.,1.,1.}, con el método de las potencias sin normalizar para obtener una aproximación del valor propio dominante y de un vector propio asociado.

b) Modifique la parte relativa a la normalización de forma que se divida en cada paso por la primera componente del vector encontrado y obtenga la aproximación del valor propio dominante como cociente entre primeras componentes.

c) Aplique cinco iteraciones para hallar la aproximación del valor propio dominante usando el método de los coeficientes de Rayleigh.

d) Obtenga los valores y vectores propios de A con las órdenes directas correspondientes y observe las diferencias entre el valor propio dominante y las aproximaciones obtenidas en los apartados a) y c).

e) Aplique el método de deflación para aproximar el segundo valor propio de A.

2.- Repita los mismos apartados del ejercicio anterior para la matriz

		(({{-34.5, 9., 7.}, {1., 0., 1.}, {0.75, 0.2, 4.5}})).

3.- Un alumno pretendió comprobar la buena marcha del método de las potencias. Para ello construyó una matriz cuyos valores y vectores propios conocía de antemano. Concretamente, consideró que los valores propios serían 7, 5 y -1 y que los vectores propios serían las filas de la matriz

		c=(({{0, 1, 1}, {1, 0, 1}, {-2, 3, 0}})).
		
Definió, pues, la matrix A como A=P (({{7, 0, 0}, {0, 5, 0}, {0, 0, -1}}))(P^(-1)), siendo P = (c^T). Obtuvo

		A=(({{17, 12, -12}, {-24, -17, 24}, {-6, -4, 11}})).

Aplicó el método de las potencias a la matriz A obteniendo las aproximaciones como cociente entre las primeras componenetes de dos iteraciones sucesivas; pero por más que cambiaba de vector inicial siempre terminaba dando una aproximación del segundo valor propio, 5, en lugar del primero, 7. Compruébelo. ¿A qué se debe?, ¿qué debió hacer?

4.- A partir de la matriz A del ejercicio 1 y del valor propio dominante obtenido en ese ejercicio, vp1, calcule una matriz B, B=A-vp1 I, siendo I la matriz identidad de orden 4. Aplique el método de las potencias a B y deshaga el desplazamiento. Observará que ha aproximado el  valor propio dominante de entre los de signo opuesto al primero.

5.- Aproxime el valor propio de móludo menor de la matriz A del ejercicio 1.