En álgebra lineal numérica , el algoritmo de valores propios de Jacobi es un método iterativo para el cálculo de los valores propios y vectores propios de una matriz simétrica real (un proceso conocido como diagonalización ). Recibe su nombre en honor a Carl Gustav Jacob Jacobi , quien propuso por primera vez el método en 1846, [1] pero recién comenzó a usarse ampliamente en la década de 1950 con la llegada de las computadoras. [2]
Descripción
Sea una matriz simétrica y una matriz de rotación de Givens . Entonces:
es simétrico y similar a .
Además, tiene entradas:
donde y .
Dado que es ortogonal y tienen la misma norma de Frobenius (la raíz cuadrada de la suma de los cuadrados de todos los componentes), sin embargo, podemos elegir tal que , en cuyo caso tiene una suma mayor de cuadrados en la diagonal:
Establezca esto igual a 0 y reorganice:
si
Para optimizar este efecto, S ij debe ser el elemento fuera de la diagonal con el valor absoluto más grande, llamado pivote .
El método de valores propios de Jacobi realiza rotaciones repetidas veces hasta que la matriz se vuelve casi diagonal. Entonces, los elementos en la diagonal son aproximaciones de los valores propios (reales) de S .
Convergencia
Si es un elemento pivote, entonces por definición para . Sea la suma de los cuadrados de todas las entradas fuera de la diagonal de . Como 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 número de rotaciones de Jacobi se denomina barrido; sea α el resultado. La estimación anterior arroja
- ;
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. Para ello, supongamos que S tiene m valores propios distintos con multiplicidades y que d > 0 es la distancia más pequeña entre dos valores propios diferentes. Llamemos a un número
Jacobi hace rotaciones con un barrido de Schönhage. Si denota el resultado, entonces
- .
Por lo tanto, la convergencia se vuelve cuadrática tan pronto como
Costo
Cada rotación de Givens 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/2n 2 elementos fuera de la diagonal. También podemos reducir esto a una complejidad O( n ) si introducimos una matriz de índice adicional con la propiedad de 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 se puede hacer en una complejidad de caso promedio O( n ): Primero, la entrada máxima en las filas actualizadas k y l se puede encontrar en O( n ) pasos. En las otras filas i , solo cambian las entradas en las columnas k y l . Al recorrer estas filas, si no es ni k ni l , basta con comparar el antiguo máximo en con las nuevas entradas y actualizar si es necesario. Si debe ser igual a k o l y la entrada correspondiente disminuyó durante la actualización, el máximo sobre la fila i debe encontrarse desde cero en una complejidad O( n ). Sin embargo, esto ocurrirá en promedio solo una vez por rotación. Por lo tanto, cada rotación tiene O( n ) y una complejidad de caso promedio de barrido O( n 3 ), lo que equivale a una multiplicación de matrices. Además, se debe inicializar antes de que comience el proceso, lo que se puede hacer en n 2 pasos.
Por lo general, el método de Jacobi converge dentro de la precisión numérica después de una pequeña cantidad de barridos. Tenga en cuenta que los valores propios múltiples reducen la cantidad de iteraciones ya que .
Algoritmo
El siguiente algoritmo es una descripción del método de Jacobi en notación matemática. Calcula un vector e que contiene los valores propios y una matriz E que contiene los vectores propios correspondientes; es decir, es un valor propio y la columna un vector propio ortonormal para , i = 1, ..., n .
procedimiento jacobi( S ∈ R n × n ; salida e ∈ R n ; salida 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 del elemento más grande fuera de la diagonal en la fila k m := k +1 para i := k +2 hasta n hacer si │ S ki │ > │ S km │ entonces m := i finsi finpara devolver m finfunc procedimiento update( k ∈ N ; t ∈ R ) ! actualiza e k y su estado y := e k ; e k := y + t si cambió k y ( y = e k ) entonces cambió k := falso; estado := estado −1 elsif (no cambió k ) y ( y ≠ e k ) entonces cambió k := verdadero; estado := estado +1 fin si finproc procedimiento rotate( k , l , i , j ∈ N ) ! realiza la rotación de S ij , S kl ┌ ┐ ┌ ┐┌ ┐ │ S kl │ │ c − s ││ S kl │ │ │ := │ ││ │ │Sij││Sij│ └ ┘ └ ┘└ ┘ proceso final ! inite e, E, y matrices ind, cambiaron E := I ; estado := n para k := 1 a n hacer ind k := maxind( k ); e k := S kk ; cambiaron k := verdadero fin para mientras estado ≠0 hacer ! siguiente rotación m := 1 ! encontrar el índice (k,l) del pivote p para k := 2 a n −1 hacer si │ S k ind k │ > │ S m ind m │ entonces m := k fin si fin para 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 fin si S kl := 0.0; actualizar( k ,− t ); actualizar( l , t ) ! rotar filas y columnas k y l para i := 1 a k −1 hacer rotate( i , k , i , l ) fin para i := k +1 a l −1 hacer rotate( k , i , i , l ) fin para i := l +1 a n hacer rotate( k , i , l , i ) fin para ! rotar vectores propios para i := 1 a n hacer ┌ ┐ ┌ ┐┌ ┐ │ E ik │ │ c − s ││ E ik │ │ │ := │ ││ │ │ E il │ │ s c ││ E il │ └ ┘ └ ┘└ ┘ ¡ endfor ! actualiza todos los ind i potencialmente modificados para i := 1 a n hace ind i := maxind( i ) endfor bucle endproc
Notas
1. La matriz lógica changed contiene el estado de cada valor propio. Si el valor numérico de o cambia durante una iteración, el componente correspondiente de changed se establece en true , de lo contrario, en false . El entero state cuenta la cantidad de componentes de changed que tienen el valor true . La iteración se detiene tan pronto como state = 0. Esto significa que ninguna de las aproximaciones ha cambiado recientemente su valor y, por lo tanto, no es muy probable que esto suceda si la iteración continúa. Aquí se supone que las operaciones de punto flotante se redondean de manera óptima al número de punto 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 inalterados. Por lo tanto, es posible restaurar S si es necesario según
para k := 1 a n −1 hacer ! restaurar la matriz S para l := k +1 a n hacer S kl := S lk fin para fin para
3. Los valores propios no están necesariamente en orden descendente. Esto se puede lograr mediante un algoritmo de ordenación simple.
para k := 1 a n −1 do m := k para l := k +1 a n hacer si e l > e m entonces m := l finsi finpara si k ≠ m entonces intercambiar e m , e k intercambiar E m , E k finsi finpara
4. El algoritmo está escrito 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 da una matriz diagonal, la implementación anterior nunca terminará, ya que ninguno de los valores propios cambiará. Por lo tanto, en las implementaciones reales, se debe agregar lógica adicional para tener en cuenta este caso.
Ejemplo
Dejar
Luego, Jacobi produce los siguientes valores propios y vectores propios después de 3 barridos (19 iteraciones):
Aplicaciones de matrices simétricas reales
Cuando se conocen los valores propios (y vectores propios) de una matriz simétrica, se calculan fácilmente los siguientes valores.
- Valores singulares
- Los valores singulares de una matriz (cuadrada) son las raíces cuadradas de los valores propios (no negativos) de . En el caso de una matriz simétrica tenemos de , por lo tanto, los valores singulares de son los valores absolutos de los valores propios de
- 2-norma y radio espectral
- La 2-norma de una matriz A es la norma basada en la norma vectorial euclidiana; es decir, el valor más grande cuando x pasa por todos los vectores con . Es el valor singular más grande de . En el caso de una matriz simétrica, es el valor absoluto más grande de sus vectores propios y, por lo tanto, igual a su radio espectral.
- Número de condición
- El número de condición de una matriz no singular se define como . En el caso de una matriz simétrica, es el valor absoluto del cociente entre el valor propio mayor y el menor. Las matrices con números de condición grandes pueden causar resultados numéricamente inestables: una pequeña perturbación puede dar como resultado 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 tiene rango si tiene columnas que son linealmente independientes mientras que las columnas restantes dependen linealmente de ellas. De manera equivalente, es la dimensión del rango de . 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. Lamentablemente, 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 sea). 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-inversa
- La pseudoinversa de una matriz es la única matriz para la cual y son simétricas y para la cual se cumple. Si no es singular, entonces .
- Cuando se llama al procedimiento jacobi (S, e, E), entonces se cumple la relación donde Diag( e ) denota la matriz diagonal con el vector e en la diagonal. Sea el vector donde se reemplaza 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 no tiene rango completo, puede que no haya una solución del sistema lineal . Sin embargo, se puede buscar un vector x para el cual sea mínimo. La solución es . En el caso de una matriz simétrica S como antes, se tiene .
- Matriz exponencial
- De uno se obtiene donde exp es el vector donde se reemplaza por . De la misma manera, se puede calcular de manera obvia para cualquier función (analítica) .
- Ecuaciones diferenciales lineales
- La ecuación diferencial tiene como solución . Para una matriz simétrica , se deduce que . Si es la expansión de por los vectores propios de , entonces .
- Sea el espacio vectorial generado por los vectores propios de los cuales corresponden a un valor propio negativo y análogamente para los valores propios positivos. Si entonces ; es decir, el punto de equilibrio 0 es atractivo para . Si entonces ; es decir, 0 es repulsivo para . y se denominan variedades estables e inestables para . Si tiene componentes en ambas variedades, entonces un componente es atraído y un componente es repelido. Por lo tanto, se aproxima como .
Implementación de Julia
El siguiente código es una implementación sencilla de la descripción matemática del algoritmo de valor propio de Jacobi en el lenguaje de programación Julia .
Usando Álgebra Lineal , Prueba función find_pivot ( Sprime ) n = tamaño ( Sprime , 1 ) pivot_i = pivot_j = 0 pivot = 0.0 para j = 1 : n para i = 1 : ( j - 1 ) si abs ( Sprime [ i , j ]) > pivot pivot_i = i pivot_j = j pivot = abs ( Sprime [ i , j ]) fin fin fin retorno ( pivot_i , pivot_j , pivot ) fin # en la práctica no se debe instanciar explícitamente la función de matriz de rotación de Givens givens_rotation_matrix ( n , i , j , θ ) G = Matrix { Float64 }( I ,( n , n )) G [ i , i ] = G [ j , j ] = cos ( θ ) G [ i , j ] = sin ( θ ) G [ j , i ] = - sin ( θ ) return G end # S es una matriz simétrica n por n n = 4 sqrtS = randn ( n , n ); S = sqrtS * sqrtS ' ; # el elemento fuera de la diagonal más grande permitido de U' * S * U # donde U son los vectores propios tol = 1e-14 Sprime = copiar ( S ) U = Matriz { Float64 }( I ,( n , n )) mientras sea verdadero ( pivot_i , pivot_j , pivot ) = find_pivot ( Sprime ) si pivote < tol break fin θ = atan ( 2 * Sprime [ pivote_i , pivot_j ] / ( Sprime [ pivote_j , pivot_j ] - Sprime [ pivote_i , pivot_i ] )) / 2 G = matriz_de_rotación_datos ( n , pivote_i , pivote_j , θ ) # actualizar Sprime y U Sprime .= G '* Sprime * G U .= U * G fin # Sprime es ahora (casi) una matriz diagonal # extraer valores propios λ = diag ( Sprime ) # ordena los valores propios (y los vectores propios correspondientes U) aumentando los valores i = sortperm ( λ ) λ = λ [ i ] U = U [ : , i ] # S debe ser igual a U * diagm(λ) * U' @test S ≈ U * diagm ( λ ) * U '
Generalizaciones
El método de Jacobi se ha generalizado a matrices hermíticas complejas , matrices reales y complejas no simétricas generales, así como a 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. S2CID 199546177.
- ^ 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. S2CID 118301078.
Lectura adicional
- Press, 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 computación científica (3.ª ed.), Nueva York: Cambridge University Press, ISBN 978-0-85-0-312-0 978-0-521-88068-8, archivado desde el original el 11-08-2011 , consultado el 13-08-2011
- Rutishauser, H. (1966). "Serie de manuales de álgebra lineal: el método de Jacobi para matrices simétricas reales". Matemática numérica . 9 (1): 1–10. doi :10.1007/BF02165223. SEÑOR 1553948. S2CID 120520713.
- Sameh, AH (1971). "Sobre algoritmos de Jacobi y similares a 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. MR 0297131.
- Shroff, Gautam M. (1991). "Un algoritmo paralelo para los valores propios y vectores propios de una matriz compleja general". Numerische Mathematik . 58 (1): 779–805. CiteSeerX 10.1.1.134.3566 . doi :10.1007/BF01385654. MR 1098865. S2CID 13904356.
- Veselić, K. (1979). "Sobre una clase de procedimientos de tipo Jacobi para diagonalizar matrices reales arbitrarias". Numerische Mathematik . 33 (2): 157–172. doi :10.1007/BF01399551. MR 0549446. S2CID 119919630.
- Veselić, K.; Wenzel, HJ (1979). "Un método de tipo Jacobi cuadráticamente convergente para matrices reales con valores propios complejos". Numerische Mathematik . 33 (4): 425–435. doi :10.1007/BF01399324. MR 0553351. S2CID 119554420.
- Yousef Saad: "Revisitando el método de rotación del subespacio de Jacobi (en bloques) para el problema de valores propios simétricos", Numerical Algorithms, vol.92 (2023), pp.917-944. https://doi.org/10.1007/s11075-022-01377-w .
Enlaces externos
- Implementación en Matlab del algoritmo de Jacobi que evita las funciones trigonométricas
- Implementación de C++11