En álgebra lineal numérica , el algoritmo de valor propio de Jacobi es un método iterativo para el cálculo de los valores propios y los vectores propios de una matriz simétrica real (un proceso conocido como diagonalización ). Lleva el nombre de Carl Gustav Jacob Jacobi , quien propuso por primera vez el método en 1846, [1] pero solo se volvió ampliamente utilizado en la década de 1950 con la llegada de las computadoras. [2]
Descripción
Dejar ser una matriz simétrica, y ser una matriz de rotación de Givens . Luego:
es simétrico y similar a.
Además, tiene entradas:
dónde y .
Desde es ortogonal, y tienen la misma norma de Frobenius (la suma de la raíz cuadrada de los cuadrados de todos los componentes), sin embargo, podemos elegir tal que , en ese caso tiene una mayor suma de cuadrados en la diagonal:
Establezca esto igual a 0 y reorganice:
Si
Para optimizar este efecto, S ij debería ser el elemento fuera de la diagonal con el mayor valor absoluto, llamado pivote .
El método del valor propio de Jacobi realiza rotaciones repetidas veces hasta que la matriz se vuelve casi diagonal. A continuación, los elementos de la diagonal son aproximaciones de los valores propios (reales) de S .
Convergencia
Si es un elemento pivote, entonces, por definición por . Dejar denotar la suma de cuadrados de todas las entradas fuera de la diagonal de . Desde tiene exactamente elementos fuera de la diagonal, tenemos o . Ahora. Esto implica o , es decir, la secuencia de rotaciones de Jacobi converge al menos linealmente por un factor a una matriz diagonal.
Un numero de Las rotaciones de Jacobi se denominan barrido; dejardenotar el resultado. La estimación anterior rinde
- ,
es decir, la secuencia de barridos converge al menos linealmente con un factor ≈ .
Sin embargo, el siguiente resultado de Schönhage [3] produce una convergencia cuadrática local. Con este fin, supongamos que S tenga m valores propios distintos con multiplicidades y sea d > 0 la distancia más pequeña de dos valores propios diferentes. Llamemos a un número de
Jacobi rota un barrido de Schönhage. Si denota el resultado entonces
- .
Por tanto, la convergencia se vuelve cuadrática tan pronto como
Costo
Cada rotación de Jacobi se puede realizar en O ( n ) pasos cuando se conoce el elemento pivote p . Sin embargo, la búsqueda de p requiere la inspección de todos los N ≈ 1/2 n 2 elementos fuera de la diagonal. También podemos reducir esto a la complejidad O ( n ) si introducimos una matriz de índice adicional con la propiedad que es el índice del elemento más grande en la fila i , ( i = 1, ..., n - 1) del S actual . Entonces los índices del pivote ( k , l ) deben ser uno de los pares. También la actualización de la matriz de índice puede hacerse en O ( n ) en el caso promedio complejidad: En primer lugar, la entrada máxima en la filas actualizado k y l se puede encontrar en O ( n ) pasos. En las otras filas i , solo cambian las entradas de las columnas k y l . Bucle sobre estas filas, sino es ni k ni l , basta con comparar el antiguo máximo en a las nuevas entradas y actualización si necesario. Sidebe ser igual a k o l y la entrada correspondiente disminuyeron durante la actualización, el máximo sobre la fila i tiene que ser encontrado desde cero en O ( n ) la complejidad. Sin embargo, esto sucederá en promedio solo una vez por rotación. Por lo tanto, cada rotación tiene O ( n ) y un barrido O ( n 3 ) complejidad de caso promedio, lo que equivale a una multiplicación de matrices. Además, eldebe inicializarse antes de que comience el proceso, lo que se puede realizar en n 2 pasos.
Normalmente, el método de Jacobi converge dentro de la precisión numérica después de un pequeño número de barridos. Tenga en cuenta que múltiples valores propios reducen el número de iteraciones desde.
Algoritmo
El siguiente algoritmo es una descripción del método Jacobi en notación matemática. Calcula un vector e que contiene los autovalores y una matriz E que contiene los autovectores correspondientes, es decir es un valor propio y la columna un vector propio ortonormal para , yo = 1, ..., n .
procedimiento jacobi ( S ∈ R n × n ; fuera e ∈ R n ; fuera E ∈ R n × n ) var i , k , l , m , estado ∈ N s , c , t , p , y , d , r ∈ R ind ∈ N n cambiado ∈ L n función maxind ( k ∈ N ) ∈ N ! índice de grande fuera de la diagonal elemento en la fila k m : = k 1 para i : = k 2 a n hacer si │ S ki │> │ S km │ entonces m : = i ENDIF endfor retorno m EndFunc actualización del procedimiento ( k ∈ N ; t ∈ R )! actualizar e k y su estado y : = e k ; e k : = y + t si cambió k y ( y = e k ) luego cambió k : = falso; estado : = estado −1 elsif (no cambió k ) y ( y ≠ e k ) luego cambió k : = verdadero; estado : = estado +1 endif endproc procedimiento rotar ( k , l , i , j ∈ N )! realizar la rotación de S ij , S kl ┌ ┐ ┌ ┐┌ ┐ │ S kl │ │ c - s ││ S kl │ │ │: = │ ││ │ │ S ij │ │ s c ││ S ij │ └ ┘ └ ┘└ ┘ endproc ! init e, E y arrays ind, cambiado E : = I ; Estado : = n para k : = 1 a n hacer ind k : = maxind ( k ); e k : = S kk ; cambiado k : = final verdadero para el estado while ≠ 0 hacer ! siguiente rotación m : = 1! encontrar índice (k, l) de pivote P para k : = 2 a n -1 hacer si │ S k ind k │> │ S m ind m │ entonces m : = k ENDIF endfor k : = m ; l : = ind m ; p : = S kl ! calcular c = cos φ, s = sen φ y : = ( e l - e k ) / 2; d : = │ y │ + √ ( p 2 + y 2 ) r : = √ ( p 2 + d 2 ); c : = d / r ; s : = p / r ; t : = p 2 / d si y <0 entonces s : = - s ; t : = - t endif S kl : = 0.0; actualizar ( k , - t ); actualizar ( l , t ) ! filas girar y columnas K y L para i : = 1 a k -1 hacer rotate ( i , k , i , l ) endfor para i : = k 1 a l -1 hago rotate ( k , i , i , l ) endfor para i : = l 1 a n hago rotate ( k , i , l , i ) endfor ! vectores propios rotar para i : = 1 a n hacer ┌ ┐ ┌ ┐┌ ┐ │ E ik │ │ c - s ││ E ik │ │ │: = │ ││ │ │ E il │ │ s c ││ E il │ └ ┘ └ ┘└ ┘ endfor ! filas k, l han cambiado, actualizar filas ind k , ind l ind k : = maxind ( k ); ind l : = maxind ( l ) bucle endproc
Notas
1. La matriz lógica modificada mantiene el estado de cada valor propio. Si el valor numérico de o cambia durante una iteración, el componente correspondiente de cambiado se establece en verdadero , de lo contrario, en falso . El estado entero cuenta el número de componentes de cambiado que tienen el valor verdadero . La iteración se detiene tan pronto como state = 0. Esto significa que ninguna de las aproximacionesha cambiado recientemente su valor y, por lo tanto, no es muy probable que esto suceda si continúa la iteración. Aquí se supone que las operaciones de coma flotante se redondean de manera óptima al número de coma flotante más cercano.
2. El triángulo superior de la matriz S se destruye mientras que el triángulo inferior y la diagonal permanecen sin cambios. Por lo tanto, es posible restaurar S si es necesario de acuerdo con
para k : = 1 a n -1 hacer ! restaurar matriz S para l : = k 1 a n do S kl : = S lk endfor endfor
3. Los valores propios no están necesariamente en orden descendente. Esto se puede lograr mediante un algoritmo de clasificación simple.
para k : = 1 a n -1 hacer m : = k para l : = k 1 a n hacer si e l > e m entonces m : = l ENDIF endfor si k ≠ m luego de intercambio e m , e k de intercambio E m , E k endif endfor
4. El algoritmo se escribe utilizando notación matricial (matrices basadas en 1 en lugar de basadas en 0).
5. Al implementar el algoritmo, la parte especificada mediante notación matricial debe realizarse simultáneamente.
6. Esta implementación no tiene en cuenta correctamente el caso en el que una dimensión es un subespacio independiente. Por ejemplo, si se le da una matriz diagonal, la implementación anterior nunca terminará, ya que ninguno de los valores propios cambiará. Por lo tanto, en implementaciones reales, se debe agregar lógica adicional para dar cuenta de este caso.
Ejemplo
Dejar
Entonces jacobi produce los siguientes autovalores y autovectores después de 3 barridos (19 iteraciones):
Aplicaciones de matrices simétricas reales
Cuando se conocen los autovalores (y autovectores) de una matriz simétrica, los siguientes valores se calculan fácilmente.
- Valores singulares
- Los valores singulares de una matriz (cuadrada) A son las raíces cuadradas de los valores propios (no negativos) de . En el caso de una matriz simétrica S tenemos de , por lo tanto, los valores singulares de S son los valores absolutos de los valores propios de S
- 2-norma y radio espectral
- La norma 2 de una matriz A es la norma basada en la norma vectorial euclidiana, es decir, el valor más grande cuando x recorre todos los vectores con . Es el mayor valor singular de A . En el caso de una matriz simétrica, es el valor absoluto más grande de sus vectores propios y, por lo tanto, es igual a su radio espectral.
- Número de condición
- El número de condición de una matriz A no singular se define como . En el caso de una matriz simétrica, es el valor absoluto del cociente del valor propio más grande y más pequeño. Las matrices con números de condición grandes pueden causar resultados numéricamente inestables: una pequeña perturbación puede resultar en grandes errores. Las matrices de Hilbert son las matrices mal condicionadas más famosas. Por ejemplo, la matriz de Hilbert de cuarto orden tiene una condición de 15514, mientras que para el orden 8 es 2,7 × 10 8 .
- Rango
- Una matriz A tiene rango r si tiene r columnas que son linealmente independientes mientras que las columnas restantes son linealmente dependientes de estas. De manera equivalente, r es la dimensión de la gama de A . Además, es el número de valores singulares distintos de cero.
- En el caso de una matriz simétrica, r es el número de valores propios distintos de cero. Desafortunadamente, debido a errores de redondeo, las aproximaciones numéricas de valores propios cero pueden no ser cero (también puede suceder que una aproximación numérica sea cero mientras que el valor verdadero no lo es). Por lo tanto, solo se puede calcular el rango numérico tomando una decisión sobre cuáles de los valores propios están lo suficientemente cerca de cero.
- Pseudo-inverso
- La pseudo inversa de una matriz A es la matriz única para los cuales AX y XA son simétricos y para los cuales AXA = A, XAX = X se cumple. Si A no es singular, entonces ' .
- Cuando se llama al procedimiento jacobi (S, e, E), entonces la relación se mantiene donde Diag ( e ) denota la matriz diagonal con el vector e en la diagonal. Dejar denotar el vector donde es reemplazado por Si y por 0 si es (numéricamente cercano a) cero. Dado que la matriz E es ortogonal, se deduce que la pseudoinversa de S está dada por .
- Solución de mínimos cuadrados
- Si la matriz A no tiene rango completo, puede que no haya una solución del sistema lineal Ax = b . Sin embargo, se puede buscar un vector x para el cual es mínimo. La solucion es . En el caso de una matriz simétrica S como antes, uno tiene .
- Matriz exponencial
- De uno encuentra donde exp e es el vector donde es reemplazado por . De la misma manera, f ( S ) se puede calcular de manera obvia para cualquier función (analítica) f .
- Ecuaciones diferenciales lineales
- La ecuación diferencial x ' = Ax , x (0) = a tiene la solución x ( t ) = exp ( t A ) a . Para una matriz simétrica S , se deduce que . Si es la expansión de a por los autovectores de S , entonces .
- Dejar ser el espacio vectorial generado por los autovectores de S que corresponden a un autovalor negativo y análogamente para los valores propios positivos. Si luego es decir, el punto de equilibrio 0 es atractivo ax ( t ). Si luego , es decir, 0 es repulsivo ax ( t ). y son llamados estables y inestables colectores para S . Si a tiene componentes en ambos colectores, entonces un componente es atraído y otro es repelido. Por tanto, x ( t ) se aproxima como .
Generalizaciones
El Método Jacobi se ha generalizado a matrices hermitianas complejas , matrices complejas y reales no simétricas generales, así como matrices de bloques.
Dado que los valores singulares de una matriz real son las raíces cuadradas de los valores propios de la matriz simétrica también se puede utilizar para el cálculo de estos valores. Para este caso, el método se modifica de tal manera que S no debe calcularse explícitamente, lo que reduce el peligro de errores de redondeo . Tenga en cuenta que con .
El Método Jacobi también es adecuado para el paralelismo.
Referencias
- ↑ Jacobi, CGJ (1846). "Über ein leichtes Verfahren, die in der Theorie der Säkularstörungen vorkommenden Gleichungen numerisch aufzulösen" . Diario de Crelle (en alemán). 1846 (30): 51–94. doi : 10.1515 / crll.1846.30.51 .
- ^ Golub, GH; van der Vorst, HA (2000). "Cálculo de valores propios en el siglo XX" . Revista de Matemática Computacional y Aplicada . 123 (1–2): 35–65. doi : 10.1016 / S0377-0427 (00) 00413-1 .
- ^ Schönhage, A. (1964). "Zur quadratischen Konvergenz des Jacobi-Verfahrens". Numerische Mathematik (en alemán). 6 (1): 410–412. doi : 10.1007 / BF01386091 . Señor 0174171 .
Otras lecturas
- Presione, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007), "Sección 11.1. Transformaciones de Jacobi de una matriz simétrica" , Recetas numéricas: El arte de la informática científica (3ª ed.), Nueva York: Cambridge University Press, ISBN 978-0-521-88068-8
- Rutishauser, H. (1966). "Álgebra lineal serie manual: el método de Jacobi para matrices simétricas reales". Numerische Mathematik . 9 (1): 1–10. doi : 10.1007 / BF02165223 . Señor 1553948 .
- Sameh, AH (1971). "Sobre los algoritmos de Jacobi y Jacobi para una computadora paralela" . Matemáticas de la Computación . 25 (115): 579–590. doi : 10.1090 / s0025-5718-1971-0297131-6 . JSTOR 2005221 . Señor 0297131 .
- Shroff, Gautam M. (1991). "Un algoritmo paralelo para los autovalores y autovectores de una matriz compleja general". Numerische Mathematik . 58 (1): 779–805. CiteSeerX 10.1.1.134.3566 . doi : 10.1007 / BF01385654 . Señor 1098865 .
- Veselić, K. (1979). "En una clase de procedimientos similares a Jacobi para diagonalizar matrices reales arbitrarias". Numerische Mathematik . 33 (2): 157-172. doi : 10.1007 / BF01399551 . Señor 0549446 .
- Veselić, K .; Wenzel, HJ (1979). "Un método similar a Jacobi cuadráticamente convergente para matrices reales con valores propios complejos". Numerische Mathematik . 33 (4): 425–435. doi : 10.1007 / BF01399324 . Señor 0553351 .
enlaces externos
- Implementación de Matlab del algoritmo de Jacobi que evita funciones trigonométricas
- Implementación de C ++ 11