algoritmo de valores propios de Jacobi

En álgebra lineal numérica , el algoritmo de valores propios de Jacobi es un método iterativo para el cálculo de los valores y vectores propios de una matriz simétrica real (un proceso conocido como diagonalización ). Recibe su nombre de Carl Gustav Jacob Jacobi , quien propuso el método por primera vez en 1846, ^{[ 1 ]} pero no se popularizó hasta la década de 1950 con la llegada de las computadoras. ^{[ 2 ]}

Este algoritmo es intrínsecamente un algoritmo para matrices densas : obtiene poca o ninguna ventaja al aplicarse a una matriz dispersa, y destruirá la dispersión al generar relleno. Del mismo modo, no preservará estructuras como la disposición en bandas de la matriz sobre la que opera.

Dejar${\displaystyle S}$sea ​​una matriz simétrica y${\displaystyle G=G(i,j,\theta )}$Sea una matriz de rotación de Givens . Entonces:

${\displaystyle S'=G^{\top }SG\,}$

es simétrico y similar a${\displaystyle S}$.

Además,${\displaystyle S^{\prime }}$tiene entradas:

${\displaystyle {\begin{aligned}S'_{ii}&=c^{2}\,S_{ii}-2\,sc\,S_{ij}+s^{2}\,S_{jj}\\S'_{jj}&=s^{2}\,S_{ii}+2sc\,S_{ij}+c^{2}\, S_{jj}\\S'_{ij}&=S'_{ji}=(c^{2}-s^{2})\,S_{ij}+sc\,(S_{ii}-S_{jj})\\S'_{ik}&=S'_{ki}=c\,S_{ik}-s\,S_{jk}&k\neq i,j\\S'_{jk}&=S'_{kj}=s\,S_{ik}+c\,S_{jk}&k\neq i,j\\S'_{kl}&=S_{kl}&k,l\neq i,j\end{aligned}}}$

dónde${\displaystyle s=\sin(\theta )}$y${\displaystyle c=\cos(\theta )}$.

Desde${\displaystyle G}$es ortogonal,${\displaystyle S}$y${\displaystyle S^{\prime }}$tienen la misma norma de Frobenius${\displaystyle ||\cdot ||_{F}}$(la raíz cuadrada de la suma de los cuadrados de todos los componentes), sin embargo podemos elegir${\displaystyle \theta }$de tal manera que${\displaystyle S_{ij}^{\prime }=0}$, en cuyo caso${\displaystyle S^{\prime }}$tiene una suma de cuadrados mayor en la diagonal:

${\displaystyle S'_{ij}=\cos(2\theta )S_{ij}+{\tfrac {1}{2}}\sin(2\theta )(S_{ii}-S_{jj})}$

Establezca esto igual a 0 y reorganice:

${\displaystyle \tan(2\theta )={\frac {2S_{ij}}{S_{jj}-S_{ii}}}}$

si${\displaystyle S_{jj}=S_{ii}}$

${\displaystyle \theta ={\frac {\pi }{4}}}$

Para optimizar este efecto, S _ij debe ser el elemento fuera de la diagonal con el mayor valor absoluto , llamado pivote .

El método de valores propios de Jacobi realiza rotaciones repetidamente hasta que la matriz se vuelve casi diagonal. Entonces , los elementos de la diagonal son aproximaciones de los valores propios (reales) de S.

Si ${\displaystyle p=S_{kl}}$ es un elemento pivote , entonces por definición ${\displaystyle |S_{ij}|\leq |p|}$ para ${\displaystyle 1\leq i,j\leq n,i\neq j}$. Dejar${\displaystyle \Gamma (S)^{2}}$denota la suma de los cuadrados de todas las entradas fuera de la diagonal de${\displaystyle S}$. Desde${\displaystyle S}$tiene exactamente ${\displaystyle 2N:=n(n-1)}$ elementos fuera de la diagonal, tenemos ${\displaystyle p^{2}\leq \Gamma (S)^{2}\leq 2Np^{2}}$ o${\displaystyle 2p^{2}\geq \Gamma (S)^{2}/N}$. Ahora${\displaystyle \Gamma (S^{J})^{2}=\Gamma (S)^{2}-2p^{2}}$Esto implica ${\displaystyle \Gamma (S^{J})^{2}\leq (1-1/N)\Gamma (S)^{2}}$ o ${\displaystyle \Gamma (S^{J})\leq (1-1/N)^{1/2}\Gamma (S)}$; es decir, la secuencia de rotaciones de Jacobi converge al menos linealmente por un factor ${\displaystyle (1-1/N)^{1/2}}$ a una matriz diagonal .

Una serie de${\displaystyle N}$Las rotaciones de Jacobi se llaman barrido;${\displaystyle S^{\sigma }}$denotemos el resultado. La estimación anterior arroja

${\displaystyle \Gamma (S^{\sigma })\leq \left(1-{\frac {1}{N}}\right)^{N/2}\Gamma (S)}$;

es decir, la secuencia de barridos converge al menos linealmente con un factor ≈${\displaystyle e^{1/2}}$.

Sin embargo, el siguiente resultado de Schönhage ^{[ 3 ]} produce una convergencia cuadrática local. Para ello, supongamos que S tiene m autovalores distintos. ${\displaystyle \lambda _{1},...,\lambda _{m}}$ con multiplicidades ${\displaystyle \nu _{1},...,\nu _{m}}$ y sea d > 0 la distancia más pequeña entre dos autovalores diferentes. Llamemos a un número de

${\displaystyle N_{S}:={\frac {n(n-1)}{2}}-\sum _{\mu =1}^{m}{\frac {1}{2}}\nu _{\mu }(\nu _{\mu }-1)\leq N}$

Jacobi gira un barrido de Schönhage. Si ${\displaystyle S^{s}}$ denota el resultado entonces

${\displaystyle \Gamma (S^{s})\leq {\sqrt {{\frac {n}{2}}-1}}\left({\frac {\gamma ^{2}}{d-2\gamma }}\right),\quad \gamma$ :=\Gamma (S)} .

Por lo tanto, la convergencia se vuelve cuadrática tan pronto como ${\displaystyle \Gamma (S)<{\frac {d}{2+{\sqrt {{\frac {n}{2}}-1}}}}}$

Cada rotación de Givens se puede realizar en${\displaystyle 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, lo que significa que esta búsqueda domina la complejidad general y empuja la complejidad computacional de un barrido en ^el algoritmo clásico de Jacobi a  ${\displaystyle O(n^{4})}$. Los algoritmos competidores alcanzan${\displaystyle O(n^{3})}$complejidad para una diagonalización completa.

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,${\displaystyle e_{i}}$ es un valor propio y la columna${\displaystyle E_{i}}$ un vector propio ortonormal para${\displaystyle e_{i}}$, 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 nfunción maxind( k ∈ N ) ∈ N ! índice del elemento fuera de la diagonal más grande en la fila k m := k +1 para i := k +2 hasta n hacer si │ S ki │ > │ S km │ entonces m := i fin si fin para devolver m fin funciónprocedimiento update( k ∈ N ; t ∈ R ) ! actualizar e k y su estado y := e k ; e k := y + t si ha cambiado k y ( y = e k ) entonces ha cambiado k := falso; estado := estado −1 si no si (no ha cambiado k ) y ( y ≠ e k ) entonces ha cambiado k := verdadero; estado := estado +1 fin si finprocprocedimiento rotate( k , l , i , j ∈ N ) ! realizar rotación de S ij , S kl ┌ ┐ ┌ ┐┌ ┐ │ S kl │ │ c − s ││ S kl │ │ │ := │ ││ │ │ S ij │ │ s c ││ S ij │ └ ┘ └ ┘└ ┘ proceso final ! inicializar e, E y los arreglos ind, cambiado E := I ; estado := n para k := 1 a n hacer ind k := maxind( k ); e k := S kk ; cambiado k := verdadero fin para mientras estado ≠0 hacer ! siguiente rotación m := 1 ! encontrar í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 = sin φ 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 if y <0 then s := − s ; t := − t endif S kl := 0.0; update( k ,− t ); update( l , t ) ! rotar filas y columnas k y l para i := 1 a k −1 hacer rotate( i , k , i , l ) fin para para i := k +1 a l −1 hacer rotate( k , i , i , l ) fin para 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 │ └ ┘ └ ┘└ ┘ fin del bucle ! actualizar todos los ind i potencialmente cambiados para i := 1 a n hacer ind i := maxind( i ) fin del bucle fin del proceso

para k := 1 a n −1 hacer ! restaurar matriz S para l := k +1 a n hacer S kl := S lk fin para fin para

para k := 1 a n −1 hacer m := k para l := k +1 a n hacer si e l > e m entonces m := l fin si fin para si k ≠ m entonces intercambiar e m , e k intercambiar E m , E k fin si fin para

Cuando se conocen los valores propios (y los vectores propios) de una matriz simétrica, los siguientes valores se calculan fácilmente.

Valores singulares

Los valores singulares de una matriz (cuadrada)${\displaystyle A}$son las raíces cuadradas de los valores propios (no negativos) de${\displaystyle A^{T}A}$En caso de una matriz simétrica${\displaystyle S}$tenemos de${\displaystyle S^{T}S=S^{2}}$, por lo tanto los valores singulares de${\displaystyle S}$son los valores absolutos de los valores propios de${\displaystyle S}$.

norma 2 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${\displaystyle \|Ax\|_{2}}$cuando x recorre todos los vectores con${\displaystyle \|x\|_{2}=1}$Es el valor singular más grande de${\displaystyle A}$. En el caso de una matriz simétrica, es el mayor valor absoluto de sus autovectores 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${\displaystyle A}$se define como ${\displaystyle {\mbox{cond}}(A)=\|A\|_{2}\|A^{-1}\|_{2}}$En el caso de una matriz simétrica, es el valor absoluto del cociente del mayor y el menor valor propio. Las matrices con números de condición grandes pueden causar resultados numéricamente inestables: una pequeña perturbación puede generar grandes errores. Las matrices de Hilbert son las matrices mal condicionadas más conocidas. 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⁸ ^.

Rango

Una matriz${\displaystyle A}$tiene rango${\displaystyle r}$si tiene${\displaystyle r}$columnas que son linealmente independientes mientras que las columnas restantes dependen linealmente de estas. De manera equivalente,${\displaystyle r}$es la dimensión del rango de ${\displaystyle A}$Además, es el número de valores singulares distintos de cero.

En el caso de una matriz simétrica, r representa el número de autovalores distintos de cero. Desafortunadamente, debido a errores de redondeo, las aproximaciones numéricas de autovalores cero pueden no serlo (también puede ocurrir que una aproximación numérica sea cero mientras que el valor real no lo sea). Por lo tanto, solo se puede calcular el rango numérico decidiendo cuáles de los autovalores están lo suficientemente cerca de cero.

Pseudoinversa

La pseudoinversa de una matriz${\displaystyle A}$es la matriz única${\displaystyle X=A^{+}}$para qué${\displaystyle AX}$y${\displaystyle XA}$son simétricas y para las cuales${\displaystyle AXA=A,XAX=X}$se sostiene. Si${\displaystyle A}$es no singular, entonces${\displaystyle A^{+}=A^{-1}}$.

Cuando se llama al procedimiento jacobi (S, e, E), entonces la relación${\displaystyle S=E^{T}{\mbox{Diag}}(e)E}$donde Diag( e ) denota la matriz diagonal con el vector e en la diagonal. Sea${\displaystyle e^{+}}$denotemos el vector donde ${\displaystyle e_{i}}$es reemplazado por${\displaystyle 1/e_{i}}$si ${\displaystyle e_{i}\leq 0}$y por 0 si ${\displaystyle e_{i}}$es (numéricamente cercano a) cero. Dado que la matriz E es ortogonal, se deduce que la pseudoinversa de S viene dada por ${\displaystyle S^{+}=E^{T}{\mbox{Diag}}(e^{+})E}$.

Solución de mínimos cuadrados

Si matriz${\displaystyle A}$ no tiene rango completo, puede que no haya solución del sistema lineal ${\displaystyle Ax=b}$. Sin embargo, se puede buscar un vector x para el cual ${\displaystyle \|Ax-b\|_{2}}$es mínimo. La solución es${\displaystyle x=A^{+}b}$. En el caso de una matriz simétrica S como antes, se tiene${\displaystyle x=S^{+}b=E^{T}{\mbox{Diag}}(e^{+})Eb}$.

exponencial matricial

De${\displaystyle S=E^{T}{\mbox{Diag}}(e)E}$uno encuentra${\displaystyle \exp S=E^{T}{\mbox{Diag}}(\exp e)E}$donde exp ${\displaystyle e}$ es el vector donde${\displaystyle e_{i}}$es reemplazado por${\displaystyle \exp e_{i}}$. Del mismo modo,${\displaystyle f(S)}$ se puede calcular de forma obvia para cualquier función (analítica)${\displaystyle f}$.

Ecuaciones diferenciales lineales

La ecuación diferencial ${\displaystyle x'=Ax,x(0)=a}$tiene la solución${\displaystyle x(t)=\exp(tA)}$Para una matriz simétrica${\displaystyle S}$De ello se deduce que ${\displaystyle x(t)=E^{T}{\mbox{Diag}}(\exp te)Ea}$. Si ${\displaystyle a=\sum _{i=1}^{n}a_{i}E_{i}}$es la expansión de${\displaystyle a}$ por los autovectores de${\displaystyle S}$, entonces${\displaystyle x(t)=\sum _{i=1}^{n}a_{i}\exp(te_{i})E_{i}}$.

Dejar ${\displaystyle W^{s}}$sea ​​el espacio vectorial generado por los autovectores de${\displaystyle S}$que corresponden a un valor propio negativo y${\displaystyle W^{u}}$De forma análoga para los autovalores positivos. Si ${\displaystyle a\in W^{s}}$entonces${\displaystyle {\mbox{lim}}_{t\rightarrow \infty }x(t)=0}$; es decir, el punto de equilibrio 0 es atractivo para${\displaystyle x(t)}$. Si${\displaystyle a\in W^{u}}$entonces${\displaystyle {\mbox{lim}}_{t\rightarrow \infty }x(t)=\infty }$; es decir, 0 es repulsivo a ${\displaystyle x(t)}$.${\displaystyle W^{s}}$y${\displaystyle W^{u}}$se denominan variedades estables e inestables para${\displaystyle S}$. Si${\displaystyle a}$Si tiene componentes en ambas variedades, entonces un componente es atraído y el otro repelido. Por lo tanto${\displaystyle x(t)}$aproches${\displaystyle W^{u}}$como${\displaystyle t\to \infty }$.

El siguiente código es una implementación directa de la descripción matemática del algoritmo de valores propios de Jacobi en el lenguaje de programación Julia .

usando Álgebra Lineal , Pruebafunción find_pivot ( Sprime ) n = tamaño ( Sprime , 1 ) pivot_i = pivot_j = 0 pivot = 0.0para j = 1 : n para i = 1 : ( j - 1 ) si abs ( Sprime [ i , j ]) > pivote pivote_i = i pivote_j = j pivote = abs ( Sprime [ i , j ]) fin fin findevolver ( 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-14Sprime = copiar ( S ) U = Matriz { Float64 }( I ,( n , n ))mientras sea verdadero ( pivot_i , pivot_j , pivot ) = encontrar_pivote ( Sprime )si pivot < tol break finθ = atan ( 2 * Sprime [ pivote_i , pivot_j ] / ( Sprime [ pivote_j , pivot_j ] - Sprime [ pivote_i , pivot_i ] )) / 2G = givens_rotation_matrix ( n , pivot_i , pivot_j , θ )# actualizar Sprime y U Sprime .= G '* Sprime * G U .= U * G fin# Sprime ahora es (casi) una matriz diagonal # extraer valores propios λ = diag ( Sprime )# Ordenar los valores propios (y los vectores propios correspondientes U) por valores crecientes i = sortperm ( λ ) λ = λ [ i ] U = U [ : , i ]# S debería ser igual a U * diagm(λ) * U' @test S ≈ U * diagm ( λ ) * U '

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 por bloques.

Dado que los valores singulares de una matriz real son las raíces cuadradas de los valores propios de la matriz simétrica${\displaystyle S=A^{T}A}$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${\displaystyle JSJ^{T}=JA^{T}AJ^{T}=JA^{T}J^{T}JAJ^{T}=B^{T}B}$ con ${\displaystyle B\,:=JAJ^{T}}$.

• Implementación en Matlab del algoritmo de Jacobi que evita las funciones trigonométricas.
• Implementación en C++11