Diagonalización de matrices, 
Semejanza.

Introducción

Si A es una matriz real cuadrada de dimensión n, el problema de la diagonalización de A consiste en:
- hallar una matriz regular P;
- y hallar una matriz diagonal D tal  que:

                                                       A=P D (P^(-1)).
NOTA:
Cuando dos matrices cuadradas A y D cualesquiera verifican una relación del tipo anterior, 
diremos que ambas matrices son semejantes.

Supuesto que A sea diagonalizable, veamos quienes han de ser estas matrices P y D. 
Escribamos P fijándonos en sus columnas, es decir P =(({{\(\(
\)\(x\_\(\(1\)\(\ \)\) | \(x\_2 | \(\(...\)\(\ \)\(\(|\)\(\ \ \)\(x\_\(\(n\)\(\ \)\)\)\)\)\)\)\)}, { }})) siendo (Subscript[x, i]) la i-ésima columna de P,
y llamemos \(TraditionalForm\`λ\_1, \ λ\_2, \ \(\(...\) \(\(,\)\(\ \)\(λ\_n\)\)\)\) a los números que aparecen en la diagonal de la matriz D.
Entonces, reescribiendo la relación entre A, P y D de la siguiente forma 

                                     A=P D (P^(-1))   =>    AP = P D ,
                                                     
y fijándonos en la i-ésima columna de las dos igualdades obtenemos:

                                                       A (Subscript[x, i ])=(Subscript[λ, i]) (Subscript[x, i]),     

para cada i de 1 a n. 
Cuando se encuentra un número real λ  y un vector no nulo x que verifican la relación A x  =  λ  x, diremos que:
 - λ es un valor propio o autovalor de la matriz A;
 - x es un vector propio o autovector de A asociado al valor propio λ. 

Por lo tanto, encontrar la matriz D equivale a encontrar n valores propios (no necesariamente distintos), que serán los que conformen la diagonal de D, y encontrar la matriz  P equivale a encontrar n vectores propios linealmente independientes (para que P sea regular) asociados a sus respectivos valores propios que conformarán las columnas de P. Así, las columnas de P formarían una base de (ℝ^n) compuesta únicamente por vectores propios de la matriz A. 
Es fundamental, pues, hallar los valores propios de A y los vectores propios asociados. Como los primeros hacen que el sistema A x  =  λ  x
tenga solución x distinta de la solución cero, la matriz de coeficientes A - λ I (donde I denota la matriz identidad de orden n) debe tener determinante no nulo.
Este determinante det(A-λ I) es un polinomio en λ de grado n y se denomina polinomio característico de A.
Por lo tanto los valores propios de A serán los ceros del polinomio característico de A al que se suele notar como 

                                                (Subscript[p, A])(λ) = det(A - λ I).

Por otro lado, el conjunto de vectores propios de A asociados a un mismo valor propio λ forman un subespacio vectorial de (R^n) que se llama subespacio propio asociado al valor propio λ.
Para concluir si una matriz A es o no diagonalizable bastará pues averiguar si hay "suficientes" valores propios para construir D y si hay "suficientes" vectores propios linealmente independientes asociados;  esta información nos la dará la dimensión de los subespacios propios y queda recogida en el siguiente resultado.

Ejemplo 1, matriz diagonalizable.

Veamos cómo se pueden usar las órdenes de Mathematica para estudiar si una matriz es diagonalizable.

A={{2, 2, -2}, {6, 9, -17}, {6, -3, -5}};
A//MatrixForm

i=IdentityMatrix[3];

valorp=Eigenvalues[A]

Como todos los valores propios son distintos, todos los subespacios propios son de dimensión uno y por lo tanto la matriz es diagonalizable. Hallamos una base del correspondiente al primer valor propio, le vamos añadiendo los vectores que forman las base de los restantes subespacios propios y así disponemos de las columnas de la matriz P de paso.

Aunque podríamos usar la orden Eigenvectors[A]  para determinar los vectores propios de A, vamos a determinarlos secuencialmente,
hallando las soluciones de los sistemas lineales homogéneos correspondientes a cada uno de ellos mediante la orden NullSpace[  ]

vectorp=NullSpace[A-valorp[[1]]*i];
MatrixForm[Transpose[vectorp]]

vectorp=AppendTo[vectorp,NullSpace[A-valorp[[2]]*i][[1]]];
MatrixForm[Transpose[vectorp]]

vectorp=AppendTo[vectorp,NullSpace[A-valorp[[3]]*i][[1]]];
MatrixForm[Transpose[vectorp]]

La matriz diagonal se forma fácilmente con el comando DiagonalMatrix[ lista ]

diag=DiagonalMatrix[valorp];
diag//MatrixForm

Y la matriz de paso P la construímos colocando los vectores propios por columnas:

P=Transpose[vectorp];
P//MatrixForm

Es inmediato comprobar la exactitud de los resultados encontrados: 
basta estudiar la anulación de la matriz A-P.diag.Inverse[P].

A==P.diag.Inverse[P]

Ejemplo 2, matriz no diagonalizable.

Consideramos, ahora, la matriz A, a la que le calculamos los valores propios.

A={{0,1,0},{0,0,1},{2,-5,4}};
valorp=Eigenvalues[A]

El valor propio 1 es doble y el valor propio 2 es simple. Para que sea diagonalizable, el subespacio propio asociado al valor propio 1 tendrá que tener dimensión 2.

vectorp=NullSpace[A-valorp[[1]]*i]

Como es un subespacio de dimensión uno y el valor propio 1 es doble, la matriz A no es diagonalizable.

Ejemplo 3, matriz simétrica.

Veamos ahora cómo diagonalizar una matriz simétrica de manera que la matriz de paso resulte ortogonal.

A={{-1,1,0},{1,-1,0},{0,0,-2}};
A//MatrixForm

valorp = Eigenvalues[A]

En este caso, el valor propio λ = -2 es doble y (como sabemos que A es diagonalizable, por ser simétrica) el subespacio propio tiene dimensión 2.

vectorp=NullSpace[A-(-2)*IdentityMatrix[3]]

Para que P resulte ortogonal necesitamos que los dos vectores base de este espacio propio sean ortonormales. En este caso, ya son perpendiculares, por lo que únicamente hay que normalizar (si no lo fuesen podemos usar el método de Gram Schmidt):

\(\(\(vectorp[\([1]\)] = 
        vectorp[\([1]\)]\/\@\(vectorp[\([1]\)] . vectorp[\([1]\)]\);\)\(
\)
    \)\)
vectorp[[2]] = 
        vectorp[[2]]/(Sqrt[vectorp[[2]] . vectorp[[2]]]);

Por último, añadimos el vector propio asociado al valor propio  λ = 0, (es decir, una base del núcleo de A). Como sabemos, vectores propios asociados a valores propios diferentes siempre son linealmente independientes y además, en el caso de matrices simétricas, son perpendiculares, por lo que bastará normalizar este tercer vector propio para tener construida la base ortonormal:

vectorp=AppendTo[vectorp,NullSpace[A][[1]]];
vectorp[[3]]=vectorp[[3]]/Sqrt[vectorp[[3]].vectorp[[3]]]

Construímos la matriz diagonal respetando el mismo orden de valores propios {-2, -2, 0}:

diag=DiagonalMatrix[valorp];
diag//MatrixForm

Y la matriz de paso P la construímos colocando los vectores propios por columnas:

P=Transpose[vectorp];
P//MatrixForm

Por último, comprobamos tanto que P resulta ortogonal, como la relación de Semejanza:

P.Transpose[P]==IdentityMatrix[3]
A==P.diag.Transpose[P]

Aplicación: potencias de una matriz.

En ocasiones es necesario calcular potencias "elevadas" de una matriz dada o incluso "raices". En estos casos es útil disponer de una matriz semejante a la dada a la que sea fácil calcularle esas potencias o ríces para disminuir el coste computacional. La simplificación se basa en que:

                                       A=P D (P^(-1))     =>     (A^m) =  P (D^m) (P^(-1 )),
                                       
Si  la matriz D es diagonal, sus potencias se calculan simplemente haciendo la misma potencia (incluso fraccionaria) de cada uno de los elementos de la diagonal, lo que conlleva menos cómputo que hacerlo para una matriz arbitraria.

Para la sigiente matriz
                                                        A=(({{1, 2, ..., 5}, {2, 3, ..., 6}, {:, :, \, :}, {5, 6, ..., 10}})), 
    
calculemos por un lado (A^8000) y por otro, calculemos una matriz B tal que (B^3) = A. A la matriz B se le llama "una raíz cúbica de la matriz A".

A = Table[i + j - 1, {i, 5}, {j, 5}];
MatrixForm[A]

Como la matriz es simétrica, sabemos que es diagonalizable. La diagonalizamos.

i=IdentityMatrix[5];

valorp=Eigenvalues[A]

Dado que el conjunto de valores propios sigue siendo sencillo, repetimos el proceso anterior. Hágase también usando la orden Eigenvectors y obsérvense los resultados.

vectorp=NullSpace[A-valorp[[1]]*i];

vectorp=AppendTo[vectorp,NullSpace[A-N[valorp[[4]]]*i][[1]]];

vectorp=AppendTo[vectorp,NullSpace[A-N[valorp[[5]]]*i][[1]]];

Obsérvese que hemos usado el valor numerico de los valores propios. ¿Qué ocurre si no lo hacemos?

Por último construimos D y P.

diag=DiagonalMatrix[valorp];
MatrixForm[diag]

P=Transpose[vectorp];
MatrixForm[P]

Comprobamos la semejanza y calculamos las potencias pedidas.
Puesto que el penúltimo valor propio es negativo, para evitar una raíz cúbica compleja lo hacemos positivo y luego le volvemos a colocar el signo negativo en su sitio.

A==P.diag.Inverse[P]

valorp8000 = N[valorp]^8000;
valorpcubic = Sign[valorp] N[Abs[valorp]]^(1/3);
diag8000 = DiagonalMatrix[valorp8000];
diagcubic = DiagonalMatrix[valorpcubic];

B  =  P . diagcubic . Inverse[P];
MatrixForm[B]

\(Print[\*""\!\(A\^8000\) = ""\ ]\)
MatrixForm[P . diag8000 . Inverse[P] ]

Verificamos la construcción de B.
También comprobamos (en esta matriz 5 x 5) la diferencia de tiempo si usamos la orden MatrixPower.

MatrixPower[B, 3] == A

MatrixPower[A, 8000];

Ejercicios

1.- Diagonaliza la siguiente matriz:                         A=(({{5, 2, 1}, {2, 7, 0}, {1, 0, -3}})).

2.- Verifica, con los datos obtenidos en el ejercicio anterior, la relación:

                                                                                              (A^5) =  P (D^5) (P^(-1 )).

3.- Calcula los valores propios de la matriz  A=(({{2, 5, 0}, {0, 3, 5}, {0, 0, 3}})) y las dimensiones de los subespacios propios asociados a cada uno.
Determina si es o no diagonalizable.

4.- Comprueba que la matriz simétrica A=(({{-1, 1, 0}, {1, -1, 0}, {0, 0, -2}})) es diagonalizable calculando los valores propios y los subespacios propios.