Las matrices y el cálculo numérico

Introducción

Las matrices tienen un relevante papel en el cálculo (y análisis) numérico, pues muchos algoritmos se fundamental en ellas, especialmente los algoritmos de resolución numérica de sistemas de ecuaciones lineales que podemos escribir de la forma $Ax=b$, donde en general, $A$ es una matriz de una aplicación lineal de un endomorfismo $f:V \rightarrow V$ sobre un cuerpo $\mathbb{K}$ —que puede ser $\mathbb{R}$ o bien $\mathbb{C}$— de tamaño $m \times n$; $x$ es el vector columna de las incógnitas (de tamaño $n\times 1$), y $b$ es el vector columna de los términos independientes (de tamaño $m\times 1$). Existen dos clases de métodos numéricos para resolver sistemas de ecuaciones lineales: los métodos directos, entre los que son muy importantes los derivados del método de Gauss —el que se estudia en el bachillerato, junto con el método de Cramer o de los determinantes— y a los que se les da el nombre de métodos gaussianos, como los basados en la factorización $A=LU$ (m. de Doolittle, m. de Crout, y m. de Cholesky); los que se basan en la factorización $A=QR$ (m. de Gram-Schmidt, m. de Householder), y los iterativos (por contraposición a los métodos directos) como son el m. de Gauss-Seidel, el m. de Jacobi, y el m. de sobrerrelajación. La razón de ser de estos métodos numéricos es la de procurar la eficiencia —claramente, por ejemplo, el método de Cramer no es eficiente— así como la estabilidad del mismo en la resolución de un problema, especialmente cuando la matriz $A$ es de gran tamaño.

Definiciones preliminares

Veamos primero algunos tipos de matrices para afianzar los conceptos a la hora de ponernos a estudiar métodos numéricos matriciales, a la vez que iremos haciendo diversos comentarios relacionados con la aplicación de las matrices a los métodos numéricos de resolución de sistemas de ecuaciones lineales, cálculo de determinantes, cálculo de la matriz inversa asociada a una matriz regular, y cálculo de valores.

Los elementos de una matriz $A\in \mathcal{M}_{m,n}$ (espacio vectorial de las matrices en un cuerpo, que, en general, será $\mathbb{C}$, de $m$ filas y $n$ columnas) se denotan por $m_{ij}$; para $i=1,\ldots,m$ y $j=1,\ldots,n$

Dada una matriz $A\in \mathcal{M}_{m,n}$ compleja, designamos por $A^{\top}$ a su matriz traspuesta; por $\bar{A}$ su matriz conjugada, y por $A^{\ast}$ (en algunos libros se denota por $A^\dagger$) a la matriz adjunta de $A$, que es la traspuesta de la matriz conjugada, esto es $A^{\ast}\equiv(\bar{A})^{\top}=\overline{A^{\top}}$.

Propiedades muy útiles en los cálculos con matrices rectangulares

1. $(A+B)^\top=A^\top+B^\top$
2. $\overline{A+B}=\bar{A}+\bar{B}$
3. $(AB)^\top=B^\top\,A^\top$
4. $\overline{AB}=\bar{A}\,\bar{B}$

Hablemos de las matrices cuadradas de orden $n$

Para matrices cuadradas, $m=n$, designamos por $I$ a la matriz identidad, donde sus elementos son $\delta_{ij}$ (delta de Kronecker), es decir, $\delta_{ij}=\left\{\begin{matrix}1&\text{si}&i=j\\0&\text{si}&i\neq j \end{matrix}\right.$

Dada una matriz cuadrada $A$ de tamaño $m \times m$, llamamos matrices principales $(A)_k$ de dicha matriz a las submatrices $k\times k$ formadas por las primeras $k$ filas y las primeras $k$ columnas. Para designar a los determinantes de esas matrices principales (menores principales) escribiremos $\Delta_k \equiv \text{det}(A)_k$

Es importante distinguir entre los diversos tipos que puede presentar una matriz cuadrada $A$:

• $A$ es singular si $\text{det}\,A=0$
• $A$ es regular si $\text{det}\,A\neq 0$ (si una matriz es regular, entonces tiene asociada una matriz inversa $A^{-1}$ que es única, de manera que $AA^{-1}=A^{-1}A=I$), y lo recíproco también es cierto. Algunas proposiciones importantes en el cálculo numérico en relación a las matrices invesibles son las siguientes:
  • Toda matriz inversible admite una factorización $A=QR$, donde $Q$ es una matriz ortogonal ($Q=Q^\top$) y $R$ es una matriz triangular superior. Para obtener la matriz $Q$ deberemos obtener una base ortogonal. Dos métodos importantes para ortogonalizar son los de Gram-Schmidt y el de Householder.
• $A$ es hermitiana si $A^{\ast}=A$. En el caso de que el cuerpo sea el de los números reales, esto se traduce en $A^{\top}=A$ y decimos que $A$ es una matriz simétrica.
• $A$ es unitaria si $A^{-1}=A^{\ast}$. En el caso de que el cuerpo sea el de los números reales, esto se traduce en $A^{-1}=A^{\top}$ y decimos que $A$ es una matriz ortogonal.
• $A$ es una matriz banda $(p,q)$ si $a_{ij}=0$ para $i\ge j+p$ ó $j \ge i+q$ y puede ser:
  • $A$ Hessenberg superior si $p=2$ y $q=n$.
  • $A$ Hessenberg inferior si $p=n$ y $q=2$.
  • $A$ triangular superior si $p=1$ y $q=n$.
  • $A$ triangular inferior si $p=n$ y $q=1$.
  • $A$ es: diagonal, tridiagonal, cuatridiagonal, pentadiagonal, ... si $p=q=1,2,3,5,\ldots$, respectivamente.
• $A$ es definida positiva si es hermitiana ($A^{\ast}=A$) y $x^{\ast}Ax\gt 0$ para todo vector $x\neq 0$.Algunas propiedades muy útilies de las matrices definidas positivas son las siguientes:
  • Recordemos el importante criterio de Sylvester que proporciona una condición necesaria y suficiente para que una matriz sea definida positiva (véase [5]), que es válido únicamente para matrices simétricas (en el caso de que $\mathbb{K}=\mathbb{R}$) o hermitianas (en el caso de que $\mathbb{K}=\mathbb{C}$). Así, para una matriz simétrica (definida sobre el cuerpo de los números reales), $A^\top=A$, decimos que es definida positiva si y sólo si todos sus menores principales —los menores principales de una matriz cuadrada son los determinantes de las submatrices principales de orden $k\le n$, esto es, los menores de orden inferior o igual a $n$ cuyos elementos son simétricos con respecto de la diagonal principal de la matriz $A$— son positivos: $\Delta_k\gt 0$, $\forall\, k=1,\ldots,n$.
  • Si una matriz $A$ es definida positiva, entonces $A$ es regular, y por tanto inversible.
  • Si $A$ es una matriz definida positiva, entonces puede resolverse un sistema de ecuaciones lineales $Ax=b$ por el método de Gauss sin intercambiar ni ecuaciones ni incógnitas; una consecuencia de ello es que $A$ admite una descomposición como el producto de una matriz triangular inferior $L$ y una matriz triangular superior $U$, de la forma $A=LU$.
  • Una matriz definida positiva y simétrica admite una factorización $A=LL^\top$ (factorización de Cholesky), donde $L$ es una matriz triangular inversible. Lo recíproco también es cierto.
• $A$ es estrictamente diagonal dominante si $|a_{ii}|\gt \displaystyle \sum_{j=i+1}^{n}\,|a_{ij}|$, para $i=1,\dots,n$. Algunas de las propiedades de estas matrices son las siguientes:
  1. Si $A$ es una matriz estrictamente diagonal dominante, entonces $A$ es regular, y por tanto inversible.
  2. Si $A$ es una matriz estrictamente diagonal dominante, entonces puede resolverse un sistema lineal $Ax=b$ por Gauss, sin intercambiar ni ecuaciones ni incóngitas, con lo cual $A$ admite una descomposición $A=LU$
• Si los determinantes principales (menores principales) $\Delta_k$ de $A$ no se anulan, entonces $A$ admite una factorización $A=LU$
• Diremos que dos matrices $A$ y $B$ son semejantes si existe una matriz inversible $C$ —a la que denominaremos matriz de transformación de semejanza de $A$ y $B$— tal que se verifique $B=C^{-1}AC$. En el caso de que $A$ sea semejante a una matriz diagonal, diremos que $A$ es diagonalizable; y si se puede encontrar una por cajas de Jordan (por cajas) que sea semejante a $A$, diremos que se puede descomponer por Jordan.
  • Lema de Shur. Toda matriz $A$ es semejante a una matriz triangular superior mediante una transformación de semejanza con una matriz unitaria.

Subrayemos también estas propiedades referidas al determinante y a la matriz inversa

1. $\text{det}\,A^\top=\text{det}\,A$
2. $\text{det}\,AB=\text{det}\,A\cdot\text{det}\,B$
3. En general, $\text{det}\,A+B\neq\text{det}\,A+\text{det}\,B$
4. $(AB)^{-1}=B^{-1}\,A^{-1}$

Acerca de los valores y vectores propios de una matriz cuadrada $A$ de orden $n$

El cálculo de los vectores propios de una matriz es muy importante en los métodos numéricos matriciales —los valors propios son los ceros del polinomio característico de una matriz $A$, y los vectores propios son los vectores asociados a los valores propios&mdash. Para calcularlos, se utilizan métodos como el de Krylov y el método de Givens, entre otros que derivan de los métodos para la resolución de sistemas de ecuaciones lineales.

Dada una matriz $A$ de orden $n$, decimos que $v\in V\equiv \mathbb{R}^n$ (espacio vectorial con el que se construye la aplicación lineal) es un vector propio de $A$ de valor propio $\lambda$ cuando se cumple que $A\,v=\lambda\,v$, condición que puede escribirse de la forma $(A-\lambda\,I)\,v=0_V$, y tal cosa es cierta si y sólo si se cumple que su valor propio $\lambda$ satiface la ecuación característica $p_{A}(\lambda)\equiv \text{det}\,(A-\lambda\,I)=0$. El primer miembro de esa igualdad es un polinomio al cual denominamos polinomio característico de $A$. Un concepto importante en el análisis numérico con métodos matriciales es el de radio espectral de $A$, que se suele designar como $\rho(A):=\text{máx}\,\{|\lambda_j|\}$, con $j=1, \ldots,r$ valores propios de $A$, es decir, se define como el módulo máximo de todos los valores propios de $A$. Se demuestra que el radio espectral es menor o igual que el valor de cualquier norma natural de $A$, $\rho(A)\le \left\|A\right\|_{.}$. Es sabido que los métodos iterativos para el cálculo de los valores y vectores propios, así como para resolver sistemas de ecuaciones lineales, y —tal como se demuestra a partir de resultados que todavía no se han expuesto— convergen si $\rho(A)\lt 1$, lo cual equivale a $\displaystyle \lim_{n\,\rightarrow\, \infty}\,A^n=0_{matriz}$, así como a $\displaystyle \lim_{n\,\rightarrow\, \infty}\,\left\|A^n\right\|_{.}=0$. Acerca de los vectores y valores pripios de $A$ podemos destacar las siguientes propiedades:

1. Los vectores propios asociados a valores propios distintos son linealmente independientes, luego si los valores propios de una matriz $A$ son todos distintos ésta es diagonalizable ($A$ es semejante a una matriz diagonal).
2. Dos matrices semejantes $A$ y $B$ tienen los mismos valores propios, y se cumple que $v$ es un vector propio de valor propio $\lambda$ de $A$ si y sólo si $S^{-1}v$ es un vector propio de valor propio $\lambda$ de $B$, donde $S$ es la matriz de la transformación de semejanza $B=S^{-1}AS$.
3. Si $A$ es semejante a una matriz diagonal $D$ por una transformación de semejanza $S$, esto es $D=S^{-1}AS$, entonces los elementos de la diagonal de $D$ son los valores propios de $A$, y los vectores asociados a las columnas de la matriz de la transformación de semejanza $S$ constituyen una base de vectores propios de $A$ (vectores propios por la derecha), y los vectores asociados a las filas de la matriz $P^{-1}$ constituyen una base de vectores propios por la izquierda de $A$.
4. Las matrices $A$ y $A^\top$ tienen los mismos valores propios.
5. Los vectores propios por la derecha (los vectores propios de $A$) son ortogonales a los vectores propios por la izquierda (los vectores propios de $A^\top$) para valores propios diferentes.
6. Una matriz $A$ es regular (inversible) si y sólo si todos sus valores propios son distintos de cero; y, en ese caso, los valores propios de $A^{-1}$ son los recíprocos de los de $A$; si $v$ es un vector propio de valor propio $\lambda$ de $A$ si y sólo si $v$ es un vector propio de valor propio $\dfrac{1}{\lambda}$ de $A^{-1}$.
7. Sea $p$ es un polinomio no nulo, entonces $\lambda$ es un valor propio de $A$ si y sólo si $p(\lambda)$ es un valor propio de $p(A)$. Y $v$ es un vector propio de valor propio $\lambda$ de $A$ si y sólo si es también vector propio de valor propio de $p(\lambda)$ de $p(A)$.
8. Sea una matriz $A$ por bloques de matrices cuadradas, de la forma $A=\begin{pmatrix}A_{11}&A_{12}\\0&A_{22}\end{pmatrix}$. Entonces, el conjunto de valores propios de $A$ es la unión de los conjuntos de valores propios de las submatrices cuadradas de la diagonal por bloques, $A_{11}$ y $A_{22}$. Luego, los valores propios de una matriz triangular son los elementos de la diagonal principal.
9. Sean $\lambda_1,\ldots,\lambda_n$ los valores propios de una matriz $A$ de orden $n$, que aprecen escritos tantas veces repetidos como indique su multiplicidad). Entonces:
   • $\displaystyle \sum_{i=1}^{n}\,\lambda_i=\text{traza}\,A\equiv \sum_{i=1}^{n}\,a_{ii}$
   • $\displaystyle \prod_{i=1}^{n}\,\lambda_i=\text{det}\,A$
10. El determinante de una matriz triangular es igual al producto de los elementos de la diagonal principal (consecuencia de propiedades anteriores).
11. Sea $A$ una matriz hermítica, esto es, una matriz tal que $A^\ast\equiv (\bar{A})^\top= A$. Entonces, todos sus valores propios son números reales.
12. Toda matriz hermítica, $A=A^\ast$, es semejante a una matriz diagonal $D$ mediante una transformación de semejanza cuya matriz de la transformación es una matriz unitaria $U$ (es decir $D=U^{\ast}AU$) —recordemos $U$ es unitaria si $U^{\ast}=U$—, entonces las columnas de $U$ forman una base de vectores propios ortonormales de $A$. Y en el caso de que el cuerpo sobre el que esté definida $A$ sea el de los números reales, $A$ es diagonalizable mediante una transformación de semejanza cuya matriz de la transformación es una matriz ortogonal $S=S^\top$, estando sus columnas asociadas a los respectivos vectores propios ortonormales de $A$.

Sobre las normas matriciales

Es sabido que podemos definir un espacio normado a partir de un espacio vectorial real o complejo. También podemos hablar de normas matriciales. Existe una definición genérica para la norma matricial de una matriz si bien en casos relativamente sencillos, y para matrices cuadradas, podemos trabajar con normas matriciales que se definen como herencia de las normas vectoriales y constituyen herramientas importante a la hora de analizar el comportamiento (convergencia, estabilidad,...) de algoritmos que hagan uso de matrices.

Definición genérica de una norma matricial de una matriz cuadrada

En general, siempre podemos definir una norma matricial $\left\|\,\right\|$ consistente con una norma vectorial dada (que notamos de la misma manera) mediante $$\left\| A \right\|:= \text{máx}\, \left\{ \dfrac{\left\|Ax\right\|}{\left\|x\right\|}: x \neq 0_V\right\}$$ Nota: en algunos libros (p.e. en [1]), se da esta otra condición: $$\left\| A \right\|:= \text{máx}\, \left\{ \left\| Ax \right\|: \left\| x \right\|=1 \right\}$$

Revisemos, a continuación, lo que sabemos sobre un espacio vectorial normado, para a continuación ir adentrándonos en la noción de normas matriciales.

Espacio vectorial normado

Consideremos un espacio vectorial $V$ sobre un cuerpo $\mathbb{K}$ (ya sea éste $\mathbb{R}$ o bien $\mathbb{C}^n$). Se define una norma en $V\equiv \mathbb{R}^n$ (o en $V\equiv \mathbb{C}^n$) como una aplicación $\left\|\,\right\|: V \rightarrow \mathbb{R}^{+}\cup\{0\};\,x\mapsto \left||x\right\|$, que cumple las siguiente propiedades:

1. Sea $x\in V$, entonces $\left||x\right\|=0_{\mathbb{R}}\Leftrightarrow x=0_V$
2. $\left||kx\right\|=|k|\,\left||x\right\|\;\forall\,\text{escalar}\,c\,\text{y}\,\forall\, x\in V$
3. $\left||x+y\right\|\le \left||x\right\|+\left||y\right\|\;\forall\, x,y\in V$

Nota: Recordemos que una norma en un espacio vectorial induce una métrica, que nos permite manejar distancias.

Unas normas vectoriales muy empleadas con las p-normas de Hölder: $\left\|x\right\|_{p}:=\left(\displaystyle \sum_{i=1}^{n}\,|x_i|^p\right)^{\frac{1}{p}}\;\text{donde}\;p\ge 1$, donde $|x|$ indica el módulo de un vector. Así, por ejemplo, se utilizan las siguientes normas particulares para todo conjunto de vectores $\{x_i\}\,,\,i=1\ldots,n$ de $V$ :

1. La norma suma de módulos se define haciendo $p=1$ en la norma genérica de Hölder: $\left\|x\right\|_1=\displaystyle \sum_{i=1}^{n}\,|x_i|$
2. La norma euclídea corresponde a poner $p=2$ en la norma genérica de Hölder: $\left\|x\right\|_1=\displaystyle \left(\sum_{i=1}^{n}\,|x_i|^2\right)^{\frac{1}{2}}$
3. La norma del máximo se define como el límite de la norma genérica de Hölder al hacer tender $p\rightarrow +\infty$, esto es $\left\|x\right\|_{\infty}:=\displaystyle \lim_{p\rightarrow +\infty}\,\left( \sum_{i=1}^{n}\,|x_i|^p\right)^{\frac{1}{p}}=\text{máx}_{i=1,\ldots\,n}\,\{|x_i|\}$

Normas matriciales subordinadas a alguna norma vectorial

Dicho esto, entedemos una norma matricial como una norma en el espacio vectorial de las matrices cuadradas $\mathcal{M}_{n,n}$ que cumpla el que sea multiplicativa, esto es, $\left\|AB\right\| \le \left\|A \right\|\,\left\|B \right\|\;\forall A,B\in \mathcal{M}_{n,n}$, y diremos que dicha norma matrical $\left\|\,\right\|$ es consistente con una cierta norma vectorial (está subordianda a una norma vectorial) si y sólo si $\left\|Ax\right\|\le \left\|A\right\|\,\left\|x\right\|\,\forall \,A\in \mathcal{M}_{n,n}\,\text{y}\,\forall x\in \mathbb{K}^n$

Ya se ha comentado más arriba que, en casos sencillos, pueden determinarse expresiones que proporcionen la norma matricial en función de los coeficientes de una matriz. Así, en particular, podemos hablar de las siguientes normas subordinadas a las normas vectoriales que hemos comentado arriba:

1. Norma matricial subordinada 1. Es la norma matricial subordinada a la norma vectorial suma de módulos (normal subordinada 1): $\left\|A\right\|_{1}:=\text{máx}_{1\le j \le n}\,\left\{\displaystyle \sum_{i=1}^{n}\,|a_{ij}|\right\}$ (esto es, el máximo por columnas de la suma de módulos de los elementos de todas y cada una de las columnas de la matriz)
2. Norma matricial subordinada 2. Es la norma subordinada a la norma vectorial euclídea: $\left\|A\right\|_{2}:=\rho(A^\top\,A)^{\frac{1}{2}}$ —que también puede escribirse como $\left\|A\right\|_{2}:=\sqrt{\rho(A^\top\,A)}$—, donde recordemos que $\rho$ denota el radio espectral (en este caso, de la matriz $A^\top\,A$); esto es, $\rho(A^\top\,A)=\text{máx}\,\{\text{valores propios en valor absoluto de}\,A^\top\,A\}$.
   
   En el caso de que $A$ sea simétrica ($A=A^\top$), esto se reduce a $\rho(A)$.
   
   Cabe también comentar el siguiente teorema: $\rho(A)=\text{ínfimo}\,\left\{\left\|A\right\|:\, \text{para cualquier norma matricial subordinada}\,\left\|.\right\|\right\}$. Como corolario de este teorema, se demuestra un resultado que ya había comentado más rarriba: (i) $\displaystyle \lim_{n\,\rightarrow\, \infty}\,A^n=0_{matriz}$, (ii) $\displaystyle \lim_{n\,\rightarrow\, \infty}\,\left\|A^n\right\|=0$, y (iii) $\rho(A)\lt 1$.
3. Norma matricial subordinada $\infty$. Es la norma matricial subordinada a la norma vectorial infinito: $\left\|A\right\|_{\infty}:=\text{máx}_{1\le i \le n}\,\left\{\displaystyle \sum_{j=1}^{n}\,|a_{ij}|\right\}$ (esto es, el máximo por filas de la suma de módulos de los elementos de todas y cada una de las filas de la matriz)

Otra norma matricial: la norma de Fröbenius

Además de las normas matriciales subordinadas que acabamos de comentar, también manejaremos esta otra norma específica, referida al producto escalar de matrices cuadradas, producto que se define como $A:B\,:=\text{traza}\,(A^\top\,B)$. Dicha norma es la que llamamos norma de Fröbenius y se define así: $$\left\|A\right\|_{F}:=(A:A)^{\frac{1}{2}}=\displaystyle \left( \sum_{i,j}^{n}\,a_{ij}^2 \right)^{\frac{1}{2}}$$ Téngase en cuenta que todas las normas matriciales toman el valor $1$ para la matriz identidad $A\equiv I$, pero la norma de Fröbenius toma el valor $\sqrt{n}$ para dicha matriz.

Acerca de la inestabilidad de los sitemas de eucaciones lineales $Ax=b$

En relación con la norma matricial (de la matriz de los coeficientes del sistema) es relevante también mencionar que existe una cantidad que denominamos número condición, $\text{Cond}(A)$ o $K(A)$ y que definimos de la forma $\text{Cond}(A):=\left\|A\right\|\cdot \left\|A^{-1}\right\|$. Se espera que la matriz se comporte bien si $\text{Cond}(A)\approx 1$, y un comportamiento defectuoso (por mal condicionamiento) si $\text{Cond}(A)$ es significativamente mayor que $1$.
$\diamond$

¿Y los sitemas de ecuaciones no lineales? ¿Cómo pueden resolverse numéricamente?

Solamente comentaré aquí, sin entrar en detalles, que existen métodos iterativos, como el método de Newton, que consisten en encontrar matrices apropiadas para construir una matriz iteradora. En el caso del método de Newton, dicha matriz se construye a partir de la matriz jacobiana (matriz de las derivadas parciales $\dfrac{\partial\,f_i}{\partial\,x_j}$ de las funciones $f_{i}(x_1,\ldots,x_n)=0$ ($i,j=1\ldots,n) de varias variables que corresponden a cada una de las ecuaciones del sistema (véase, p.e., [4] pp. 176-177, y un ejemplo de resolución de un sistema no lineal en las pp. 225-226).

-oOo-

Referencias:
[1] C. Moreno, Introducción al cálculo numérico (UNED, Madrid, 2011).
[2] A. Aubanell, A. Benseny y A. Deslshams, Útiles de cálculo numérico (Labor, Barcelona, 1993).
[3] B.P., Demidovich y I.A. Maron, Cálculo numérico fundamental (Paraninfo, Madrid, 1988).
[4] F. García y A. Nevot, Ejercicios resueltos de cálculo numérico (Paraninfo, Madrid, 1992).
[5] vv.aa., https://ca.wikipedia.org/wiki/Criteri_de_Sylvester (Vikipèdia).