3.3.4.3 Creación de matrices sparse
El uso de la función sparse no es nunca una buena opción para crear una matriz sparse porque no tiene absolutamente ninguna lógica. El motivo es claro, si la existencia de las matrices sparse es ahorrar memoria y aumentar la velocidad no tiene ningún sentido que tengamos que crear antes la matriz original y luego tener que liberar memoria.
Entonces pensamos que sería muy cómodo contar con las mismas funciones de creación de matrices convencionales para el caso sparse:
- speye
- Crea la matriz unidad en formato sparse
- spones
- Crea una matriz sparse con unos en los elementos que deseemos
- sprand
- Crea una matriz sparse de números aleatorios en posiciones aleatorias o siguiendo un determinado patrón.
Por ejemplo, si queremos una matriz sparse de números aleatorios con una densidad de 0.01 y de dimensiones 100ױ00: >> randsparse=sprand(100,100,0.01);
Para ver cómo se reparten aleatoriamente los números en la matriz utilizaremos la función spy que ya conocemos: >> spy(randsparse)
Y como resultado obtenemos la figura 3.1:
Figure 3.1: Elementos no nulos de una matriz sparse creada con sprand
- spdiags
- Extrae e introduce diagonales en matrices sparse.
Al igual que su función homóloga diag tiene un gran sentido físico. De hecho tiene mucha mayor utilidad en la función de creación spdiags que diag porque las matrices en banda entran en la definición de matrices sparse. La única consideración a tener en cuenta es que las matrices en banda se resuelven rápidamente de manera exacta mediante una factorización LU diseñada para matrices sparse. Utilizaremos los métodos iterativos con matrices que no cumplan ningún patrón evidente o sencillo.
3.3.4.4 Manipulación y operaciones con matrices sparse
- sphcat
- Concatena horizontalmente dos matrices sparse de dimensiones compatibles.
- spvcat
- Concatena verticalmente dos matrices sparse de dimensiones compatibles.
- spabs
- Calcula el valor absoluto de los elementos complejos de una matriz sparse sin cambiar el carácter de la misma.
- spreal
- Calcula la parte real de los elementos complejos de una matriz sparse sin cambiar el carácter de la misma.
- spimag
- Calcula la parte imaginaria de los elementos de una matriz sparse sin cambiar el carácter de la misma.
- spinv
- Calcula la matiz inversa de una matriz sparse.
Sin embargo no es recomendable que utilicemos spinv para resolver sistemas de ecuaciones lineales. Como en el caso de sistemas de ecuaciones en matrices convencionales es recomendable que optemos por el operador para invertir la matriz del sistema.
3.3.5 Matrices tridiagonales (Octave)
Son uno de los tipos de matrices más buscadas; tienen un gran sentido físico (provienen de algunas ecuaciones en diferencias) y son muy adecuadas para el cálculo numérico. Son un caso particular de matrices sparse pero al ser tan sencillas pueden resolverse directamente sin realizar demasiados cálculos. La función de resolución de sistemas de ecuaciones con una matriz tridiagonal es trisolve.
3.3.6 Matrices no regulares.
Una de las características más importantes del operador es que también sirve en el caso de tener un sistema con exceso o defecto de ecuaciones. El operador ya no calcula la matriz inversa sino que busca una solución cuyo error cumple la condición de la norma mínima o busca el subespacio solución. En este caso se dice que en vez de buscar la inversa se calcula una pseudoinversa que resuelve la ecuación. El tratamiento de este tipo de problemas suele omitirse en la mayoría de cursos de álgebra en ingeniería, sin embargo tiene una gran utilidad en cálculo numérico y es una herramienta relativamente importante. Creo que no es en vano si intentamos entender con un poco más de profundidad el tratamiento de estas aplicaciones lineales sin entrar en formalismos.
3.3.6.1 Singular Value Decomposition (SVD)
Se demuestra que cualquier matriz de aplicación lineal acepta una descomposición del tipo:
A=UwV⊤
donde U y V⊤ son matrices ortogonales y w una matriz diagonal que contiene los valores singulares. Las matrices ortogonales, o más concretamente las matrices de una transformación ortogonal, tienen la propiedad de mantener los ángulos en la transformación. Las matrices de giro son un ejemplo de transformación ortogonal. La matriz w, al ser diagonal, no produce ningún tipo de rotación, simplemente expande o contrae determinadas direcciones. La demostración de la existencia de esta descomposición va mucho más allá de los objetivos de este libro. Se utiliza la notación V⊤ para esta matriz debido a que para cumplir la definición de matiz de transformación ortogonal sus columnas deben ser vectores ortogonales; es sólo un formalismo porque se demuestra también que la inversa de una transformación ortogonal es también una transformación ortogonal y es equivalente a su traspuesta. Ayuda escribir la anterior ecuación especificando las dimensiones:
A(NנM)=U(NנM)w(MנM)V(MנM)⊤
Vemos también que las dimensiones son correctas, aplicando que V⊤≡ V-1:
A(NנM)V(MנM)=U(N M)w(MנM)
Ahora, un repaso de álgebra de parte de [14]. Si nuestra matriz A es la matriz de una aplicación lineal
f:S—→ T
podemos definir dos subespacios en S; el Núcleo, subespacio que forman los vectores que tienen al vector nulo como imagen y el subespacio Imagen que es el subespacio f(S)de T. Por ejemplo, sean S y T los siguientes espacios reales:
- S son las funciones polinómicas p(x)=a+bx+cx2+dx3
- T el espacio de funciones de R en R
Consideremos la aplicación lineal de f:S—→ T definida por f(p(x))=p′(x), la derivada según la variable independiente; o sea: f(a+bx+cx2+dx3)=b+2cx+3dx2. La imagen de f es el subespacio de T formado por los polinomios de grado menor o igual que 2. El núcleo de f es el subespacio de T que forman los polinomios constantes.
Si volvemos a la forma de la descomposición propuesta por la SVD,
A(NנM)=U(NנM)w(MנM)V(MנM)⊤
Tomemos la condición de N<M, es decir, un sistema con menos ecuaciones que incógnitas. En este caso la aplicación lineal será f:RM—→RN, si ninguna de las filas es combinación lineal del resto tendremos una imagen de dimensión N y un núcleo de dimensión M-N. Se cumple además que la cantidad de valores singulares no nulos será N y que las M-N columnas asociadas a los valores singulares nulos formarán la base del núcleo. La demostración de lo anterior está también fuera de los objetivos de este libro.
Esta operación es perfectamente válida también para el caso N>M en el que tenemos más ecuaciones que incógnitas.
Para obtener las tres matrices de la descomposición Matlab proporciona la siguiente función
- svd
- Descomposición en valores singulares. Si se pide únicamente un valor de salida retorna un vector con valores singulares. Con tres valores de salida retorna la descomposición matricial completa.
Para ejemplificar el uso de la ecuación vamos a intentar descomponer en valores singulares la ya clásica matriz... octave:3> a=[1,2,3;4,5,6;7,8,9];
octave:4> [U,w,V]=svd(a)
U =
-0.21484 0.88723 0.40825
-0.52059 0.24964 -0.81650
-0.82634 -0.38794 0.40825
w =
16.84810 0.00000 0.00000
0.00000 1.06837 0.00000
0.00000 0.00000 0.00000
V =
-0.479671 -0.776691 -0.408248
-0.572368 -0.075686 0.816497
-0.665064 0.625318 -0.408248
Efectivamente, uno de los valores singulares es nulo debido a que la matriz es singular. Para analizar más en profundidad las matrices, Matlab proporciona además las siguientes funciones:
- rank
- Calcula el rango de la matriz, es decir, el número de valores singulares no nulos o la dimensión de la imagen de la aplicación.
- null
- Calcula la dimensión del núcleo de una aplicación lineal expresada por una matriz, esto es, el número de valores singulares nulos.
3.3.6.2 Problemas con defecto o exceso de ecuaciones.
Si lo único que nos interesa es resolver ax=b... ¿Qué calcula exactamente el operador cuando N<M? Lo que hace es calcular la pseudoinversa a partir de la descomposición en valores singulares. La fórmula de la pseudoinversa es la siguiente:
A+=Vw-1U⊤
Esta matriz da una solución al problema del tipo:
x=a+b+v
donde v es una combinación lineal cualquiera de vectores del núcleo definido por la matriz a. Se demuestra que la solución que minimiza la norma del error r≡|ax-b| (mínimos cuadrados) es precisamente la que cumple v=0. Acabamos de deducir la forma de cualquier solución hallada por el método de los mínimos cuadrados.
En el apartado 6.1.1.1 analizaremos las propiedades de la SVD en el cálculo de ajustes polinómicos de series de datos.
3.4 Autovalores
- eig
- Calcula los autovalores y autovectores de una matriz cuadrada.
- 1
- El operador llama las rutinas de las bibliotecas ATLAS y LAPACK, una auténtica maravilla de la tecnología.
- 2
- Las matrices sparse son parte de los cursos avanzados de análisis numérico y no suelen entrar en los planes de estudios de las carreras de ingeniería. Son sin embargo una herramienta muy potente y en algunos casos imprescindible para resolver problemas de ingeniería bastante típicos. No podríamos resolver problemas de grandes estructuras sin la existencia del almacenamiento de matrices sparse y su relación con los métodos iterativos. Merece la pena echarle un vistazo a la bibliografía aunque lo más importante es conocer su existencia y los síntomas que fuerzan su uso.
- 3
- Octave informa del tipo de almacenamiento utilizado, Matlab no dice absolutamente nada.
