Métodos de subespacios de Krylov

y sus aplicaciones a la Dinámica Browniana

Edwin Armando Bedolla Montiel

12 de febrero, 2021, Grupo de Materia Blanda.

Contenido

  1. ¿Qué son los métodos de subespacios de Krylov?
    • ¿En qué problemas se pueden emplearse?
    • ¿Cómo se implementan?
  2. Ejemplos: sistemas lineales y regresión lineal.
  3. Interacciones hidrodinámicas en dinámica Browniana.

Definición

Sea \( A \) una matriz invertible de \( n \times n \), y sea \( b \) un vector de dimensión \( n \), entonces, un subespacio de Krylov de orden \( r \) está definido como \[\begin{equation} \mathcal{K}_r(A, b) = \text{span} \{ b, Ab, A^2 b, A^3 b, \dots, A^{r-1} b\} \end{equation} \]

Aproximaciones

Esto significa que podemos crear un subespacio lineal usando solamente el producto \( Ab. \)
Para visualizar esto, podemos pensar en las series de Taylor \[ f(x-a)=\sum_{n=0}^{\infty}\frac{f^{(n)}(a)}{n!}(x-a)^n \]
Las series de Taylor nos dan una buena aproximación de la función \( f(x). \)
\[ e^x = \sum_{n=0}^{\infty} \frac{x^n}{n!} \]
\[ e^1 = \sum_{n=0}^{N=10} \frac{1^n}{n!}=2.718281801146385 \] \[ e = 2.718281828459045 \]
De la misma forma, el subespacio de Krylov \( \mathcal{K}_r \) nos permite construir un espacio donde, con pocos vectores base, podemos encontrar una buena aproximación.
¿En qué tipo de problemas puedo aplicar estos métodos?

Regresión lineal

Se quiere resolver el problema \[ \mathbf{y}=\mathbf{X \beta} + \alpha + \mathbf{\epsilon} \] donde \( \beta \) es el vector de coeficientes y \( \alpha \) es conocido como el sesgo u ordenada al origen.
Se reorganiza el problema a un sistema lineal \[ \mathbf{X^{\mathsf{T}} y} = \mathbf{X^{\mathsf{T}} X \beta} \] donde ahora se multiplica por la matriz transpuesta para que la matriz \( X^{\mathsf{T}} X \) sea una matriz cuadrada.
Para problemas pequeños, la solución al problema es simple. Uno de lo más conocidos es mediante la descomposición \( LU \).
Toda matriz cuadrada tiene descomposición \[ PA = LU \] donde \( P \) es una matriz de permutaciones entre las filas; \( L \) es una matriz triangular inferior; y \( U \) es una matriz tringular superior.
La solución es, entonces \[ P \mathbf{y} = P \mathbf{X^{\mathsf{T}} X \beta} = LU \beta \] Es decir, primero se resuelve el problema \[ Lz = P\mathbf{y} \] y luego \[ U \mathbf{\beta} = z \]
Pero ¡ojo!, nótese que la operación \[ \mathbf{X^{\mathsf{T}} \beta} \] da como resultado a un vector, y aún más, ¡se parece a nuestra definición de \( \mathcal{K}_r \)!

Método del gradiente conjugado

Podemos entonces construir un subespacio de Krylov tal que \[ \mathcal{K}_r (\mathbf{X}^{\mathsf{T}} \mathbf{X}, \mathbf{X}^{\mathsf{T}} \beta) \equiv \text{span} \{ \mathbf{X}^{\mathsf{T}} \beta, (\mathbf{X}^{\mathsf{T}} \mathbf{X}) \mathbf{X}^{\mathsf{T}} \beta, \\ (\mathbf{X}^{\mathsf{T}} \mathbf{X})^2 \mathbf{X}^{\mathsf{T}} \beta, \\ (\mathbf{X}^{\mathsf{T}} \mathbf{X})^3 \mathbf{X}^{\mathsf{T}} \beta, \dots, (\mathbf{X}^{\mathsf{T}} \mathbf{X})^{r-1} \mathbf{X}^{\mathsf{T}} \beta \} \]
Dado que es un espacio lineal nos preguntamos, entonces, ¿cuál es la mejor combinación lineal que resuelve el problema?
El método del gradiente conjugado nos dice que la mejor combinación es aquella donde se cumpla que \( \mathbf{r} \) definido como \[ \mathbf{r}=\mathbf{y} - \mathbf{X} \beta \] y conocido como residuo sea ortogonal a \( \mathcal{K}_r \).
Sin embargo, es impráctico crear todo el espacio \( \mathcal{K}_r \) pues no se necesitan todos los vectores base.
La solución es simple, de forma iterativa, vamos ortogonalizando (semejante al proceso de Gram-Schmidt) hasta que nuestra base tenga suficientes vectores tal que el resultado es suficientemente bueno.
Esta base tiene una propiedad, todo vector residuo es un vector conjugado respecto a los vectores del espacio \( \mathcal{K}_r \).
Esto significa que, mientras ortogonalizamos \( \mathcal{K}_r \) iremos ortogonalizando el nuevo espacio \( \mathcal{K}_{residuos} \).

Implementación

Este algoritmo es rápido, y tiene la gran ventaja de que es fácil de implementar.

Definiciones y asignaciones

Sea \( A=\mathbf{X^{\mathsf{T}} X} \) y \( \mathbf{b} = \mathbf{X^{\mathsf{T}} y} \), entonces \[ \begin{align} & x_0 := \mathbf{0} \text{ solución inicial} \\ & \mathbf{r}_0 := \mathbf{b} - \mathbf{A x}_0 \\ & \mathbf{p}_0 := \mathbf{r}_0 \\ & k := 0 \\ \end{align} \] Aquí, \( \mathbf{p} \) es la notación para los vectores base que formarán el espacio lineal de los residuos, \( \mathcal{K}_{residuos} \).
\[ \begin{align} & \text{inicia bucle} \\ & \qquad \alpha_k := \frac{\mathbf{r}_k^\mathsf{T} \mathbf{r}_k}{\mathbf{p}_k^\mathsf{T} \mathbf{A p}_k} \\ & \qquad \mathbf{x}_{k+1} := \mathbf{x}_k + \alpha_k \mathbf{p}_k \\ & \qquad \mathbf{r}_{k+1} := \mathbf{r}_k - \alpha_k \mathbf{A p}_k \\ & \qquad \hbox{si } \mathbf{r}_{k+1} \text{ es suficientemente pequeño, salir del bucle} \\ & \qquad \beta_k := \frac{\mathbf{r}_{k+1}^\mathsf{T} \mathbf{r}_{k+1}}{\mathbf{r}_k^\mathsf{T} \mathbf{r}_k} \\ & \qquad \mathbf{p}_{k+1} := \mathbf{r}_{k+1} + \beta_k \mathbf{p}_k \\ & \qquad k := k + 1 \\ & \text{termina bucle} \\ & \text{regresa } \mathbf{x}_{k+1} \text{ como el resultado} \end{align} \]

\( \alpha_k \) son los coeficientes de la combinación lineal en \( \mathcal{K}_{residuos} \).

\( \beta_k \) son los coeficientes de la combinación lineal en \( \mathcal{K}_r \).

Condiciones

El método del gradiente conjugado solamente funciona si la matriz \( A \) es simétrica y definida positiva.

Si no lo es, otro método debe emplearse.

Alternativas

  1. GMRES (Generalized minimal residual method)
  2. BiCG (Biconjugate gradient method)
  3. BiCGSTAB (Biconjugate gradient stabilized method)

Dinámica Browniana

\[ \mathbf{r}(t+\Delta t)=\mathbf{r}(t)+\frac{\Delta t}{k_B T} \mathbf{DF}+(\nabla \cdot \mathbf{D}) \Delta t + \mathbf{g} \] donde \( \langle \mathbf{g} \rangle = 0 \) y \( \langle \mathbf{gg^T} \rangle = 2 \mathbf{D}\Delta t .\)
\( \mathbf{D} \) es conocido como el tensor hidrodinámico, y normalmente se escoge uno que cumpla con \(\nabla \cdot \mathbf{D}=0 \). El más común es el tensor de Rotne-Prager y su modificación, el tensor de Rotne-Prager-Yamakawa.
\( \mathbf{g} \) son los desplazamientos estocásticos, que siguen una distribución normal multivariada con vector de medias zero y matriz de covarianza \( 2 \mathbf{D}\Delta t .\)
Normalmente, esto se calcula así \[ \mathbf{g}=\sqrt{2\Delta t}\mathbf{y}=\sqrt{2\Delta t}\mathbf{Bz} \] con \[ \mathbf{D}=\mathbf{BB^{\mathsf{T}}} \] y \( \mathbf{z} \) sigue una distribución normal estándar univariada.
En el algoritmo original de Ermak-McCammon, esto se hace con una descomposición de Cholesky, pero considerando que el tensor \( \mathbf{D} \) es de tamaño \( 3n \times 3n \) con \( n \) el número de partículas, esto se vuelve impráctico rápidamente.

En el 2012, Ando et al. propusieron emplear métodos de subespacios de Krylov.

Se dieron cuenta que el problema es calcular \( \mathbf{D^{1/2}} ,\) esto porque \(\mathbf{D}=\mathbf{D^{1/2}} \mathbf{(D^{1/2})}^{\mathsf{T}} \)

Además, se dieron cuenta que no necesitan saber \( \mathbf{D^{1/2}} \) sino \( \mathbf{D^{1/2} z} .\)

Después de todo, solamente son desplazamientos estocásticos, no deben ser precisos, sino aproximados.

Propuesta

Construir un subespacio de Krylov para realizar esta aproximación, tal que \[ \mathcal{K}_m(\mathbf{D, z})=\text{span } \{ \mathbf{z}, \mathbf{Dz}, \mathbf{D^{1/2}z}, \dots, \mathbf{D^{m-1}z} \} \]

Asignaciones

\[ \mathbf{z} \sim \mathcal{N}(0, \mathbf{I}) \\ j := 0 \text{ contador de bucle} \\ \qquad v_0 := \frac{\mathbf{z}}{\lVert \mathbf{z} \rVert} \\ m := 30 \text{ el tamaño del subespacio de Krylov} \\ \]
\[ \begin{align} & \text{inicia bucle, desde j hasta m} \\ & \qquad \mathbf{w} :=\mathbf{D} \mathbf{v}_j \\ & \qquad \text{si } j > 0 \\ & \qquad \quad \mathbf{w}=\mathbf{w}-h_{j-1,j} \mathbf{v}_{j-1} \\ & \qquad h_{j,j}=\mathbf{w^T} \mathbf{v}_j \\ & \text{si } j < m \\ & \qquad \quad \mathbf{w}=\mathbf{w}-h_{j,j} \mathbf{v}_{j} \\ & \qquad \quad h_{j+1,j}=h_{j,j+1}=\lVert \mathbf{w} \rVert \\ & \qquad \quad \mathbf{v}_{j+1}=\frac{\mathbf{w}}{h_{j+1,j}} \\ \end{align} \]
El resultado aproximado es \[ \widetilde{\mathbf{y}} = \lVert \mathbf{z} \rVert \mathbf{V}_m \mathbf{H}_m^{1/2} \mathbf{e}_1 \]

Se hizo una aproximación de la matriz \( D^{1/2} \); esta aproximación es la matriz \( H^{1/2} .\)

Se va ortognalizando el espacio según se necesiten vectores.

Pasamos de un espacio de \( 3n \times 3n \) a un espacio \( m \times m \) donde \( m \ll 3n \), en este caso, \( m \equiv 30 .\)


                                    znorm = dnrm2( s,xr,1 )
                                    vm(:,1) = xr / znorm
                                    do i = 1,m
                                        call dgemv( 'n',s,s,1.0_dp,sigma,s,vm(:,i),1,0.0_dp,w,1 )
                                        if ( i > 1 ) then
                                            w = w - ( h(i-1,i) * vm(:,i-1) )
                                        end if
                                        k = size(w, 1)
                                        h(i,i) = ddot( k,w,1,vm(:,i),1 )
                                        if ( i < m ) then
                                            w = w - ( h(i,i) * vm(:,i) )
                                            k = size(w, 1)
                                            h(i+1,i) = dnrm2( k,w,1 )
                                            h(i,i+1) = h(i+1,i)
                                            vm(:,i+1) = w / h(i+1,i)
                                        end if
                                    end do
                                    call sqrt_matrix( h,temp )
                                    call dgemv( 'N',m,m,1.0_dp,temp,m,eid,1,0.0_dp,v,1 )
                                    call dgemv( 'N',s,m,1.0_dp,vm,s,v,1,0.0_dp,res,1 )
                                    r = znorm * res

¿Dónde encontrar estas implementaciones?

Fortran

  • PETSc/TAO, sistemas lineales, sistemas de ecuaciones diferenciales ordinarias y parciales. (Argonne National Laboratory)
  • ARPACK para aplicar métodos de subespacios de Krylov para encontrar valores y vectores propios de sistemas muy grandes. (Rice University y Sandia National Laboratory)

Python

  • Scipy , contiene una colección de estos métodos, algunos son envolturas (e.g. ARPACK), pero ya no es necesario implementarlos.

Julia

  • IterativeSolvers , es un paquete con una colección de rutinas para sistemas lineales.
  • Krylov , acepta no solamente matrices, sino cualquier tipo de operador lineal.

Conclusiones

  1. Los métodos de subespacios de Krylov son métodos iterativos para encontrar resultados aproximados.

  2. Dependiendo de las propiedades de las matrices (u operadores lineales), se deben emplear variaciones de estos métodos.

  3. Se estudió un caso particular aplicado a la dinámica Browniana.

Referencias

  • Strang, G. (1993). Introduction to linear algebra (Vol. 3). Wellesley, MA: Wellesley-Cambridge Press.
  • Demmel, J. W. (1997). Applied numerical linear algebra. Society for Industrial and Applied Mathematics.
  • Van der Vorst, H. A. (2003). Iterative Krylov methods for large linear systems (No. 13). Cambridge University Press.
  • Hestenes, M. R., & Stiefel, E. (1952). Methods of conjugate gradients for solving linear systems (Vol. 49, No. 1). Washington, DC: NBS.

Referencias

  • Cipra, B. A. (2000). The best of the 20th century: Editors name top 10 algorithms. SIAM news, 33(4), 1-2.
  • Gutknecht, M. H. (2007). A brief introduction to Krylov space methods for solving linear systems. In Frontiers of Computational Science (pp. 53-62). Springer, Berlin, Heidelberg.
  • Ermak, D. L., & McCammon, J. A. (1978). Brownian dynamics with hydrodynamic interactions. The Journal of chemical physics, 69(4), 1352-1360.

Referencias

  • Ando, T., Chow, E., Saad, Y., & Skolnick, J. (2012). Krylov subspace methods for computing hydrodynamic interactions in Brownian dynamics simulations. The Journal of chemical physics, 137(6), 064106.
  • Shewchuk, J. (1994). An Introduction to the Conjugate Gradient Method Without the Agonizing Pain.
  • Yamakawa, H. (1970). Transport Properties of Polymer Chains in Dilute Solution: Hydrodynamic Interaction. The Journal of Chemical Physics, 53(1), 436–443. https://doi.org/10.1063/1.1673799

Proyectos

  • dbfort Dinámica Browniana en Fortran 2008+.
  • LeastSquaresSVM Machine Learning con Máquinas de Soporte Vectorial.