Los programas pueden ser implementados en cualquiera de
los lenguajes más o menos standard: Matlab, Maple, C++,
Fortran, etc. Importante: todo programa tiene que estar
acompañado de una "documentación", soporte
escrito donde se explique claramente qué se hizo, como
funciona cada una de las componentes, gráficos o tablas de los
experimentos realizados, conclusiones obtenidas, dificultades
encontradas, etc..
El trabajo puede ser hecho tanto
individualmente o como de a dos. Hay flexibilidad en los temas:
ustedes también pueden proponer o discutir conmigo otros temas
de programación que les interese hacer. Ante cualquier duda o
inquietud, no duden en contactarme por e-mail, o mejor aún,
personalmente hasta principios de Feb/2006.
1. Algoritmo de
Strassen para la inversión de matrices densas
Este
algoritmo fue explicado en el curso, muestra que la complejidad
teórica de la inversión es igual a la complejidad de la
multiplicación de matrices, y por lo tanto la inversión
puede hacerse en mejor complejidad que el N^3 del algoritmo de
eliminación. El objetivo es estudiar este algoritmo desde un
punto de vista práctico.
Implementar un algoritmo de multiplicación de matrices. En una primera etapa, se puede usar la multiplicación standard, pero luego hay que implementar alguna multiplicación rapida, por ejemplo el algoritmo clásico de Strassen o bien el del articulo de Kaporin, que tiene complejidad N^2.8 en lugar del N^3 del algoritmo standard.
Implementar el algoritmo de tipo "divide-and-conquer" de Strassen, explicado en el libro de Pan y en las notas del curso.
Determinar el rango de aplicación práctica, es decir la máxima dimensión N que puede resolver en un tiempo razonable (e.g. 10').
Estudiar su estabilidad numérica (= pérdida de precisión por errores de redondeo).
Comparar los resultados (velocidad y estabilidad) con los del algoritmo de eliminación (implementarlo basándose en el pseudo-código en las notas del curso (capítulo 4) o bien usar una versión ya implementada).
Para hacer más: la versión
presentada no es estable ya que no tiene pivotage. Testear formas de
estabilizarlo, por ejemplo usando precondicionadores estructurados,
ver [Pan, S 5.6].
Referencias:
Aho, Hopcroft, Ullman, The design and analysis of computer algorithms. Addison-Wesley 1975.
Pan, Structured matrices and polynomials: unified superfast algorithms (Chapter 4). Birkhauser, 2001.
Kaporin, A practical algorithm for matrix multiplication. Num. Linear Alg. Appl. 1999.
Notas del curso (Capítulo 5)
2. Algoritmo rápido unificado de
descomposición PLU de matrices estructuradas
El
algoritmo rápido presentado en el curso, suponía que
los operadores de desplazamiento S,T eran respectivamente diagonal y
triangular superior. Esto restringe la aplicabilidad del algoritmo,
por ejemplo la estructura tipo Toeplitz no puede ser tratada
directamente, sino que hay que reducirla al tipo Vandermonde usando
un precondicionador circulante.
En el libro de Pan (y
parcialmente en las notas del curso) se explica que esta restricción
no es necesaria, y cómo el algoritmo se adapta a operadores
quasi-diagonales, lo cual cubre todos las estructuras clásicas
en forma unificada. Las ideas son muy similares a las presentadas en
el curso, pero cambian los detalles. El objetivo es implementar la
versión rápida (complejidad O(\alpha*N^2)),
estudiar su velocidad y estabilidad numérica.
Implementar el algoritmo. Usar el algoritmo rápido presentado en el curso, pero usando la actualización de generadores para el caso quasi-diagonal siguiendo [Pan, Ch. 4].
Incorporar pivotage parcial.
Con esta base, programar su aplicación a las cuatro estructuras clásicas (Toeplitz, Hankel, Vandermonde, Cauchy). En cada caso hay que implementar fórmulas de reconstrucción de la primera línea y columna de una matriz a partir de sus generadores, ver [Pan, S 4.4].
Determinar el rango de aplicación práctica, es decir la máxima dimensión N que puede resolver en un tiempo razonable (e.g. 10').
Estudiar su estabilidad numérica (= pérdida de precisión por errores de redondeo).
Comparar los resultados (velocidad y estabilidad) con el algoritmo de eliminación (implementarlo basándose en el pseudo-código en las notas del curso (capítulo 4) o bien usar una versión ya implementada).
Referencias:
Pan, Structured matrices and polynomials: unified superfast algorithms (Chapter 4). Birkhauser, 2001.
Notas del curso (Capítulo 5)
3. Algoritmo super-rápido de decodificación
de códigos de Reed-Solomon
El algoritmo
super-rápido para matrices estructuradas explicado en el curso
puede aplicarse a la decodificación de códigos de
Reed-Solomon de longitud N en O(N*log^3(N)) ops de F_q.
Esta implementación requiere un cierto trabajo a nivel de
detalles (que hasta donde yo sé, no están escritos en
ningún lado) pero es un trabajo accesible y sumamente
interesante.
Implementar el algoritmo super-rapido adaptándolo específicamente para la estructura tipo Vandermonde. Lo que se pide es la parte de álgebra lineal, es decir el cálculo del núcleo de la matriz de Berlekamp-Massey. Se puede obviar la división final de polinomios, para la cual sería necesaria implementar aparte (o conseguir) la división rápida basada en la TFR.
Implementar la aritmética de F_q para q=2^r, o mejor aún, usar alguna implementación ya disponible.
Determinar el rango de aplicación práctica, es decir la máxima dimensión N que puede resolver en un tiempo razonable (e.g. 10').
Déterminar experimentalemente la probabilidad de fracaso del algoritmo (es decir que alguno de los pivots encontrados sea 0).
Referencias:
Olshevsky, Shokrollahi, A displacement approach to decoding algebraic codes. Contemporary Math. 1999.
Pan, Structured matrices and polynomials: unified superfast algorithms (Chapter 5). Birkhauser, 2001.
Notas del curso (Capítulo 5).
4. Multigrilla para la ecuación de Hemholtz
Se trata de resolver numéricamente la ecuación
de Hemholtz en el cuadrado:
-u_xx-u_yy+a*u_x=f
para (x,y) en el cuadrado
[0,1]^2;
u=0 para (x,y)
en el borde;
via el multigrilla. Se tratará de
determinar el comportamiento del problema con respecto a los
distintos parámetros. Se recomienda programar previamente la
ecuación de Hemholtz en dimensión 1
Discretizar la ecuación con diferencias finitas en una malla uniforme.
Implementar las distintas componentes del ciclo multigrilla y del multigrilla completo.
Experimente cambiando el parámetro de regularización, la cantidad de pasos de regularización por ciclo, ciclo V y W. Considere en particular la ecuación -e(u_xx+u_yy)+a*u_x=x(1-x)sin(lx) que tiene solución exacta u=(Ax^2+Bx+C)sin(lx) con A=-1/s, B=(1-2aA)/s, (2eA-aB)/s, y s=1/(el^2). Compare los algorithmos para e=0.01, 0.1, 1; a=0.1, 1, 10.
Referencias:
Briggs, Henson, McCormick: A multigrid tutotial. SIAM 2000.
5. Distribución del condicionamiento de
matrices aleatórias con estructuras particulares
Generar
matrices aleatórias con estructura dada: simétrica,
banda, sparse con soporte fijado de antemano, Toeplitz,
Toeplitz-block-Toeplitz, Hankel, Vandermonde,
Cauchy. Se busca obtener experimentalmente la
distribución del condicionamiento, y compararla con los
resultados de Edelman para matrices densas y coeficientes
siguiendo una distribución gaussiana.
Aprender a generar coefficientes siguiendo una distribución dada, aprender a calcular el condicionamiento de una matriz (se puede y se recomienda usar las rutinas ya implementadas en Matlab o en librerías específicas).
Manejar el programa eficientemente para calcular numerosos ejemplos en la mayor dimensión posible. Sería ideal disponer de máquinas para hacer grandes cálculos, no sé (pero puedo averiguarlo) si hay en la UBA ni son disponibles para los estudiantes.
Comparar las distribuciones encontradas con las obtenidas teóricamente por Edelman. En primer lugar verificar experimentalmente los resultados teóricos: el caso denso (real o complejo) con el condicionamiento de Demmel k_D; para luego pasar a las distintas estructuras.
Referencias:
Edelman, Ph.D. Thesis. MIT 1989.
Schatzman, Les bonnes matrices sont rares, aparecerá en Matapli.