Normas matriciales y condicionamiento
Normas inducidas y número de condición de matrices
Hasta ahora hemos medido el tamaño de un vector con una norma. Para hablar del condicionamiento de un sistema lineal completo necesitamos también medir el tamaño de una matriz, y resulta que no cualquier manera razonable de hacerlo sirve para ese propósito.
1 De vectores a matrices: la norma de Frobenius
Una matriz \(A\) de tamaño \(m\times n\) se puede acomodar como un vector de \(mn\) entradas, simplemente apilando sus columnas una encima de la otra. Cualquier norma vectorial aplicada a ese vector aplanado es, automáticamente, una función que cumple los tres axiomas de norma sobre el espacio de matrices (positividad, homogeneidad, desigualdad triangular), porque esos axiomas se preservan bajo cualquier reacomodo lineal de las entradas.
La elección más natural es la norma 2 del vector aplanado, que recibe nombre propio.
\[\|A\|_F = \left(\sum_{i=1}^m\sum_{j=1}^n |A_{ij}|^2\right)^{1/2}\]
Es exactamente lo que calcula la función norm de Julia cuando se le pasa una matriz.
Demuestre que la norma de Frobenius es, en efecto, una norma.
2 ¿Por qué exigir submultiplicatividad?
Los tres axiomas de norma no son, sin embargo, todo lo que en análisis numérico le pedimos a una norma matricial. La propiedad adicional que realmente importa es
\[\|AB\| \le \|A\|\cdot\|B\|,\]
llamada submultiplicatividad. Es la propiedad que garantiza que el producto de matrices sea una operación continua: sin ella no se puede acotar de forma consistente el error al multiplicar matrices, ni definir el número de condición como vamos a hacerlo más abajo. (Esto es, de hecho, un caso particular de un resultado general de análisis funcional sobre continuidad de aplicaciones bilineales acotadas; ver por ejemplo Rudin o Kreyszig, se sale del curso.)
Lo importante es que no toda norma vectorial aplicada a la matriz aplanada es submultiplicativa, aunque cumpla los tres axiomas sin problema. Considera \(\|A\|_{\max} := \max_{ij}|A_{ij}|\), que es perfectamente válida como norma vectorial. Con
\[A = B = \begin{bmatrix} 1 & 1 \\ 1 & 1\end{bmatrix}\]
tenemos \(\|A\|_{\max}=\|B\|_{\max}=1\), pero
\[AB = \begin{bmatrix} 2 & 2 \\ 2 & 2\end{bmatrix}, \qquad \|AB\|_{\max}=2 > 1 = \|A\|_{\max}\|B\|_{\max}.\]
Entonces \(\|\cdot\|_{\max}\), aunque es una norma legítima sobre el espacio de matrices, no es una norma matricial en el sentido que necesitamos.
Usando la desigualdad de Cauchy–Schwarz, demuestre que la norma de Frobenius sí es submultiplicativa.
3 Normas inducidas
Por la razón anterior, en la práctica se prefiere construir la norma matricial a partir de una norma vectorial de otra manera, viendo a \(A\) como una transformación lineal en vez de como una lista de números.
Dada una norma vectorial \(\|\cdot\|_p\), la norma matricial que induce sobre cualquier matriz \(A\) de \(m\times n\) es
\[\|A\|_p = \max_{\|x\|_p=1} \|Ax\|_p = \max_{x\neq 0} \frac{\|Ax\|_p}{\|x\|_p}.\]
El máximo existe (no es solo un supremo) porque la esfera unitaria \(\{x:\|x\|_p=1\}\) es compacta en dimensión finita y \(x\mapsto \|Ax\|_p\) es continua, así que por el teorema de Weierstrass el máximo se alcanza en algún punto.
Vale la pena señalar también que en dimensión finita (\(\mathbb{R}^n\), \(\mathcal{M}_{m\times n}\), \(P_{\leq n}\)) todas las normas (matriciales o no) son equivalentes: generan la misma noción de convergencia. Concretamente, si \(\|\cdot\|\) y \(\|\cdot\|'\) son dos normas cualesquiera sobre un espacio de dimensión finita, existen constantes \(\alpha,\beta>0\) tales que \(\alpha\|x\|\le\|x\|'\le\beta\|x\|\) para todo \(x\). Este es otro resultado de análisis funcional que no vamos a demostrar aquí.
Demuestre que \(\|\cdot\|_1\), \(\|\cdot\|_2\) y \(\|\cdot\|_\infty\) son equivalentes en \(\mathbb{R}^n\), encontrando explícitamente las constantes correspondientes.
A diferencia de lo que pasó con \(\|\cdot\|_{\max}\), la submultiplicatividad de una norma inducida no cuesta ningún trabajo extra: viene incluida gratis en la definición. Para cualquier \(x\) con \(\|x\|=1\),
\[\|ABx\| = \|A(Bx)\| \le \|A\|\cdot\|Bx\| \le \|A\|\cdot\|B\|\cdot\|x\| = \|A\|\cdot\|B\|,\]
y tomando el máximo sobre \(x\) se obtiene \(\|AB\|\le\|A\|\cdot\|B\|\) directamente, sin necesidad de ninguna desigualdad auxiliar como Cauchy–Schwarz.
Como el máximo de \(\|Ax\|\) se alcanza en algún \(x_{\max}\) de la esfera unitaria, y toda norma cumple \(\|{-x}\|=\|x\|\), el vector \(-x_{\max}\) también está en la esfera y también alcanza el máximo (porque \(A(-x_{\max})=-Ax_{\max}\) tiene la misma norma). Es decir, el maximizador nunca es único: siempre viene en un par \(\pm x_{\max}\).
4 Interpretación geométrica
Si transformamos la esfera unitaria (en la norma que estemos usando) mediante \(A\), la imagen es una especie de elipse (o hiperelipsoide) que queda encerrada en la esfera de radio \(\|A\|\), y la toca en al menos ese par de puntos antípodas. En ese sentido, \(\|A\|\) mide cuánto puede expandir \(A\) a un vector unitario en el mejor de los casos. El análogo para medir cuánto puede encoger es el mínimo en vez del máximo,
\[\sigma_{\min}(A) := \min_{\|x\|=1}\|Ax\|,\]
que también se alcanza por compacidad. En el capítulo de análisis matricial vamos a ver que \(\|A\|_2\) y \(\sigma_{\min}(A)\) son, respectivamente, el mayor y el menor valor singular de \(A\).
5 Fórmulas explícitas para \(\|A\|_\infty\) y \(\|A\|_1\)
Para la mayoría de las normas inducidas, calcular \(\|A\|_p\) exactamente requiere resolver un problema de optimización. Dos casos son la excepción.
\[\|A\|_\infty = \max_{1\le i\le m} \sum_{j=1}^n |A_{ij}| \qquad\text{(máxima suma por filas)}\] \[\|A\|_1 = \max_{1\le j\le n} \sum_{i=1}^m |A_{ij}| \qquad\text{(máxima suma por columnas)}\]
Demostración para \(\|A\|_1\). Sea \(C_j=\sum_i|A_{ij}|\) y \(C=\max_j C_j\). Para \(x\) con \(\|x\|_1\le1\), intercambiando el orden de las sumas,
\[\|Ax\|_1 = \sum_i\Big|\sum_j A_{ij}x_j\Big| \le \sum_j |x_j|\sum_i|A_{ij}| = \sum_j |x_j|\,C_j \le C\sum_j|x_j| \le C,\]
así que \(\|A\|_1\le C\). Para ver que se alcanza, sea \(j^*\) la columna con la suma máxima y toma \(x=\mathbf{e}_{j^*}\) (el vector canónico), que da \(Ax\) igual a la columna \(j^*\) de \(A\), así que \(\|Ax\|_1=C_{j^*}=C\). Combinando, \(\|A\|_1=C\).
Demuestre la fórmula para la norma \(\infty\).
Nótese que en ambas demostraciones el vector que alcanza el máximo es un vértice de la bola unitaria correspondiente y por eso hay una fórmula cerrada tan simple: basta identificar el vértice correcto. La bola unitaria \(\ell_2\), en cambio, es redonda y no tiene vértices, así que no existe una fórmula así de simple para \(\|A\|_2\) en términos de las entradas de \(A\).
6 El número de condición
Para \(A\) invertible, el número de condición (en una norma inducida dada) es
\[\kappa(A) = \|A\|\cdot\|A^{-1}\|.\]
Como \(1=\|I\|=\|AA^{-1}\|\le\|A\|\|A^{-1}\|\), siempre se tiene \(\kappa(A)\ge1\), y \(\kappa(A)=1\) es el mejor condicionamiento posible. Si \(A\) es singular, convenimos \(\kappa(A)=\infty\); y en la práctica, si \(\kappa(A)>\epsilon_{\text{mach}}^{-1}\), la matriz es efectivamente singular para efectos computacionales.
6.1 Perturbando \(b\)
Sea \(Ax=b\) y considera el sistema perturbado \(A(x+\Delta x) = b+\Delta b\) (mismo \(A\), lado derecho perturbado). Restando el sistema original,
\[A\Delta x = \Delta b \quad\Longrightarrow\quad \Delta x = A^{-1}\Delta b \quad\Longrightarrow\quad \|\Delta x\| \le \|A^{-1}\|\,\|\Delta b\|.\]
Dividiendo por \(\|x\|\) y multiplicando y dividiendo por \(\|A\|\),
\[\frac{\|\Delta x\|}{\|x\|} \le \|A^{-1}\|\|A\|\cdot\frac{\|\Delta b\|}{\|A\|\|x\|} \le \kappa(A)\,\frac{\|\Delta b\|}{\|b\|},\]
usando en el último paso que \(\|b\|=\|Ax\|\le\|A\|\|x\|\).
6.2 Perturbando \(A\)
Ahora considera \((A+\Delta A)(x+\Delta x) = b\), con el mismo \(b\). Expandiendo y usando \(Ax=b\),
\[A\Delta x + \Delta A\,x + \Delta A\,\Delta x = 0 \quad\Longrightarrow\quad \Delta x = -A^{-1}\Delta A\,(x+\Delta x).\]
Tomando normas,
\[\|\Delta x\| \le \|A^{-1}\|\,\|\Delta A\|\,\|x+\Delta x\| \quad\Longrightarrow\quad \frac{\|\Delta x\|}{\|x+\Delta x\|} \le \kappa(A)\,\frac{\|\Delta A\|}{\|A\|}.\]
Cuando \(\Delta A\to0\), también \(\Delta x\to0\) (de la propia ecuación despejada arriba), así que \(\|x+\Delta x\|\to\|x\|\), y en ese límite se obtiene
\[\frac{\|\Delta x\|}{\|x\|} \le \kappa(A)\,\frac{\|\Delta A\|}{\|A\|} \qquad (\|\Delta A\|\to0).\]
Es la razón por la que en el enunciado general aparece la aclaración “en el límite \(\|\Delta A\|\to0\)”: la desigualdad, tal como está, es de primer orden, no exacta para cualquier \(\Delta A\).
Si además se sabe que \(\|A^{-1}\Delta A\|<1\), el ejercicio de la serie de Neumann que ya resolvió (¿\((I-A)^{-1}=\sum_k A^k\) cuando \(\|A\|<1\)?) permite volver exacta la cota anterior, sin tomar ningún límite: reagrupando \(\Delta x=-A^{-1}\Delta A\,x - A^{-1}\Delta A\,\Delta x\) se llega a \((I+A^{-1}\Delta A)\Delta x=-A^{-1}\Delta A\,x\), y usando \(\|(I+E)^{-1}\|\le \frac{1}{1-\|E\|}\) (la misma cota de la serie geométrica) se obtiene
\[\frac{\|\Delta x\|}{\|x\|} \le \frac{\kappa(A)\dfrac{\|\Delta A\|}{\|A\|}}{1-\kappa(A)\dfrac{\|\Delta A\|}{\|A\|}}.\]
7 Residuo y error hacia atrás
Suponga que \(Ax=b\) y que \(\tilde x\) es una solución aproximada (por ejemplo, calculada en punto flotante). El residuo es
\[r = b - A\tilde x.\]
Por construcción, \(A\tilde x = b - r\), es decir: \(\tilde x\) es la solución exacta de un sistema con el mismo \(A\) pero lado derecho perturbado \(\tilde b = b-r\). El error hacia atrás relativo es entonces \(\|r\|/\|b\|\), y combinándolo con el teorema de la sección anterior,
\[\frac{\|\tilde x - x\|}{\|x\|} \le \kappa(A)\,\frac{\|r\|}{\|b\|}.\]
Esta desigualdad separa dos preguntas que conviene no mezclar: qué tan pequeño es el error hacia atrás (una propiedad del algoritmo), y cómo ese error se traduce en error hacia adelante \(\|\tilde x-x\|/\|x\|\) (una propiedad del problema, gobernada por \(\kappa(A)\)).
El residuo es una cantidad local y verificable: calcularlo es solo una multiplicación matriz-vector y una resta, y no depende en absoluto de qué tan sensible sea el problema que se está resolviendo. Si \(\tilde x\) proviene de un algoritmo bien diseñado (eliminación gaussiana con pivoteo parcial, por ejemplo), cada operación de punto flotante individual comete un error relativo del orden de \(\epsilon_{\text{mach}}\), y un algoritmo estable hacia atrás es, por definición, uno en el que esos errores se acumulan de manera controlada, sin amplificarse entre sí, de modo que el resultado final \(\tilde x\) termina siendo la solución exacta de un problema cercano: \(\|\Delta A\|/\|A\|\) (o \(\|\Delta b\|/\|b\|\)) del orden de \(\epsilon_{\text{mach}}\). Esta garantía es una propiedad del algoritmo y no depende de \(\kappa(A)\) en absoluto.
Esto es exactamente el eje de “problema vs. algoritmo” que hemos venido trabajando, llevado a su forma más concreta. El número de condición mide qué tan sensible es el problema: cuánto se amplifica una perturbación de los datos al pasar a la solución. La estabilidad hacia atrás mide qué tan bien se comporta el algoritmo: si logra mantener pequeño el residuo, sin importar qué tan sensible sea el problema. Son propiedades independientes: un algoritmo puede ser perfectamente estable (residuo minúsculo) y aun así producir un error hacia adelante enorme, si se aplica a un problema con \(\kappa(A)\) muy grande. La desigualdad de esta sección es precisamente la que traduce ese error hacia atrás minúsculo en el error hacia adelante, amplificado por \(\kappa(A)\).
Una demostración numérica natural para ilustrar esta distinción: tomar una matriz de Hilbert (mal condicionada), resolver \(Ax=b\) con \ en Julia, y verificar que el residuo \(\|b-A\tilde x\|\) sale minúsculo (estabilidad hacia atrás intacta) mientras que el error hacia adelante \(\|\tilde x-x\|/\|x\|\) sale grande, amplificado por \(\kappa(A)\).