Transposición de matriz en el lugar

La transposición de matriz in situ , también llamada transposición de matriz in situ , es el problema de transponer una matriz N × M in situ en la memoria de la computadora , idealmente con O (1) (limitado) almacenamiento adicional, o como máximo con almacenamiento adicional mucho menor que NM . Por lo general, se supone que la matriz se almacena en orden de fila o columna principal (es decir, filas o columnas contiguas, respectivamente, dispuestas consecutivamente).

Realizar una transposición en el lugar (transposición in situ) es más difícil cuando N ≠ M , es decir, para una matriz no cuadrada (rectangular), donde implica una permutación compleja de los elementos de datos, con muchos ciclos de longitud mayor que 2. En contraste, para una matriz cuadrada ( N = M ), todos los ciclos tienen una longitud de 1 o 2, y la transposición se puede lograr mediante un simple bucle para intercambiar el triángulo superior de la matriz con el triángulo inferior. Surgen más complicaciones si se desea maximizar la localidad de memoria para mejorar la utilización de la línea de caché o para operar fuera del núcleo (donde la matriz no cabe en la memoria principal), ya que las transposiciones implican inherentemente accesos a la memoria no consecutivos.

El problema de la transposición in situ no cuadrada se ha estudiado al menos desde fines de la década de 1950 y se conocen varios algoritmos, incluidos varios que intentan optimizar la localidad para contextos de memoria caché, fuera del núcleo o similares relacionados con la memoria.

En una computadora , a menudo se puede evitar la transposición explícita de una matriz en la memoria simplemente accediendo a los mismos datos en un orden diferente. Por ejemplo, las bibliotecas de software para álgebra lineal , como BLAS , suelen proporcionar opciones para especificar que ciertas matrices se deben interpretar en orden transpuesto para evitar el movimiento de datos.

Sin embargo, sigue habiendo una serie de circunstancias en las que es necesario o deseable reordenar físicamente una matriz en la memoria a su orden transpuesto. Por ejemplo, con una matriz almacenada en orden de fila principal , las filas de la matriz son contiguas en la memoria y las columnas son discontinuas. Si se deben realizar operaciones repetidas en las columnas, por ejemplo en un algoritmo de transformada rápida de Fourier (por ejemplo, Frigo y Johnson, 2005), la transposición de la matriz en la memoria (para hacer que las columnas sean contiguas) puede mejorar el rendimiento al aumentar la localidad de la memoria . Dado que estas situaciones normalmente coinciden con el caso de matrices muy grandes (que exceden el tamaño de la caché), se vuelve deseable realizar la transposición en el lugar con un almacenamiento adicional mínimo.

Además, como problema puramente matemático, la transposición en el lugar implica una serie de interesantes problemas de teoría de números que se han resuelto a lo largo de varias décadas.

Por ejemplo, considere la matriz 2×4:

${\displaystyle {\begin{bmatrix}11&12&13&14\\21&22&23&24\end{bmatrix}}.}$

En formato de fila principal, esto se almacenaría en la memoria de la computadora como la secuencia (11, 12, 13, 14, 21, 22, 23, 24), es decir, las dos filas almacenadas consecutivamente. Si transponemos esto, obtenemos la matriz 4×2:

${\displaystyle {\begin{bmatrix}11&21\\12&22\\13&23\\14&24\end{bmatrix}}}$

que se almacena en la memoria de la computadora como la secuencia (11, 21, 12, 22, 13, 23, 14, 24).

Si numeramos las ubicaciones de almacenamiento del 0 al 7, de izquierda a derecha, entonces esta permutación consta de cuatro ciclos:

(0), (1 2 4), (3 6 5), (7)

Es decir, el valor en la posición 0 pasa a la posición 0 (un ciclo de longitud 1, sin movimiento de datos). A continuación, el valor en la posición 1 (en el almacenamiento original: 11, 12 , 13, 14, 21, 22, 23, 24) va a la posición 2 (en el almacenamiento transpuesto 11, 21, 12 , 22, 13, 23, 14, 24), mientras que el valor en la posición 2 (11, 12, 13 , 14, 21, 22, 23, 24) va a la posición 4 (11, 21, 12, 22, 13 , 23, 14, 24), y la posición 4 (11, 12, 13, 14, 21 , 22, 23, 24) vuelve a la posición 1 (11, 21 , 12, 22, 13, 23, 14, 24). Lo mismo ocurre con los valores en la posición 7 y las posiciones (3 6 5).

En lo que sigue, suponemos que la matriz N × M se almacena en orden de fila principal con índices basados ​​en cero. Esto significa que el elemento ( n , m ), para n = 0,..., N −1 y m = 0,..., M −1, se almacena en una dirección a = Mn + m (más algún desplazamiento en la memoria, que ignoramos). En la matriz M × N transpuesta , el elemento ( m , n ) correspondiente se almacena en la dirección a' = Nm + n , nuevamente en orden de fila principal. Definimos la permutación de transposición como la función a' = P ( a ) tal que:

${\displaystyle Nm+n=P(Mn+m)\,}$a pesar de${\displaystyle (n,m)\en [0,N-1]\times [0,M-1]\,.}$

Esto define una permutación de los números . ${\displaystyle a=0,\ldots ,MN-1}$

Resulta que se pueden definir fórmulas simples para P y su inversa (Cate y Twigg, 1977). Primero:

${\displaystyle P(a)={\begin{cases}MN-1&{\text{si }}a=MN-1,\\Na{\bmod {(}}MN-1)&{\text{en caso contrario}},\end{cases}}}$

donde "mod" es la operación módulo .

En segundo lugar, la permutación inversa viene dada por:

${\displaystyle P^{-1}(a')={\begin{cases}MN-1&{\text{si }}a'=MN-1,\\Ma'{\bmod {(}}MN-1)&{\text{en caso contrario}}.\end{cases}}}$

(Esto es sólo una consecuencia del hecho de que la inversa de una transpuesta N × M es una transpuesta M × N , aunque también es fácil demostrar explícitamente que P ⁻¹ compuesto con P da la identidad).

Como demostraron Cate y Twigg (1977), el número de puntos fijos (ciclos de longitud 1) de la permutación es precisamente 1 + mcd( N −1, M −1) , donde mcd es el máximo común divisor . Por ejemplo, con N = M el número de puntos fijos es simplemente N (la diagonal de la matriz). Si N  − 1 y M  − 1 son coprimos , por otro lado, los únicos dos puntos fijos son las esquinas superior izquierda e inferior derecha de la matriz.

El número de ciclos de cualquier longitud k > 1 viene dado por (Cate y Twigg, 1977):

${\displaystyle {\frac {1}{k}}\suma _{d|k}\mu (k/d)\mcd(N^{d}-1,MN-1),}$

donde μ es la función de Möbius y la suma es sobre los divisores d de k .

Además, el ciclo que contiene a = 1 (es decir, el segundo elemento de la primera fila de la matriz) es siempre un ciclo de longitud máxima L , y las longitudes k de todos los demás ciclos deben ser divisores de L (Cate y Twigg, 1977).

Para un ciclo dado C , cada elemento tiene el mismo máximo común divisor . ${\displaystyle x\en C}$${\displaystyle d=\mcd(x,MN-1)}$

Este teorema es útil en la búsqueda de ciclos de la permutación, ya que una búsqueda eficiente sólo puede considerar múltiplos de divisores de MN −1 (Brenner, 1973).

Laflin y Brebner (1970) señalaron que los ciclos suelen presentarse en pares, lo que se aprovecha mediante varios algoritmos que permutan pares de ciclos a la vez. En particular, supongamos que s es el elemento más pequeño de algún ciclo C de longitud k . De ello se deduce que MN −1− s también es un elemento de un ciclo de longitud k (posiblemente el mismo ciclo).

A continuación se resumen brevemente los algoritmos publicados para realizar la transposición de matrices en el lugar. El código fuente que implementa algunos de estos algoritmos se puede encontrar en las referencias que aparecen a continuación.

Debido a que la transposición física de una matriz es costosa en términos computacionales, en lugar de mover valores en la memoria, se puede transponer la ruta de acceso. Es trivial realizar esta operación para el acceso a la CPU, ya que las rutas de acceso de los iteradores simplemente deben intercambiarse, ^[1] sin embargo, la aceleración de hardware puede requerir que aún se realineen físicamente. ^[2]

Para una matriz cuadrada N × N A _{n , m} = A ( n , m ), la transposición en el lugar es fácil porque todos los ciclos tienen una longitud de 1 (las diagonales A _{n , n} ) o de 2 (el triángulo superior se intercambia con el triángulo inferior). El pseudocódigo para lograr esto (asumiendo índices de matriz basados ​​en cero ) es:

para n = 0 a N - 1
     para m = n + 1 a N
        intercambia A(n,m) con A(m,n)

Este tipo de implementación, aunque simple, puede exhibir un desempeño deficiente debido a una mala utilización de la línea de caché, especialmente cuando N es una potencia de dos (debido a conflictos de línea de caché en una caché de CPU con asociatividad limitada). La razón para esto es que, a medida que m se incrementa en el bucle interno, la dirección de memoria correspondiente a A ( n , m ) o A ( m , n ) salta de manera no contigua por N en la memoria (dependiendo de si la matriz está en formato de columna principal o de fila principal, respectivamente). Es decir, el algoritmo no explota la localidad de referencia .

Una solución para mejorar la utilización de la caché es "bloquear" el algoritmo para que opere sobre varios números a la vez, en bloques dados por el tamaño de la línea de caché; desafortunadamente, esto significa que el algoritmo depende del tamaño de la línea de caché (es "consciente de la caché"), y en una computadora moderna con múltiples niveles de caché requiere múltiples niveles de bloqueo dependientes de la máquina. En cambio, se ha sugerido (Frigo et al. , 1999) que se puede obtener un mejor rendimiento mediante un algoritmo recursivo : dividir la matriz en cuatro submatrices de tamaño aproximadamente igual, transponiendo las dos submatrices a lo largo de la diagonal de forma recursiva y transponiendo e intercambiando las dos submatrices por encima y por debajo de la diagonal. (Cuando N es suficientemente pequeño, el algoritmo simple anterior se utiliza como caso base, ya que recurrir ingenuamente hasta N = 1 tendría una sobrecarga excesiva de llamadas de función). Este es un algoritmo ajeno a la caché , en el sentido de que puede explotar la línea de caché sin que el tamaño de la línea de caché sea un parámetro explícito.

En el caso de matrices no cuadradas, los algoritmos son más complejos. Muchos de los algoritmos anteriores a 1980 podrían describirse como algoritmos de "seguimiento de ciclos", es decir, que repiten los ciclos y trasladan los datos de una ubicación a la siguiente en el ciclo. En forma de pseudocódigo:

para cada longitud>1 ciclo C de la permutación
    elija una dirección de inicio s en C 
    sea D = datos en s 
    sea x = predecesor de s en el ciclo
     mientras  x ≠ s 
        mueva datos de x al sucesor de x 
        sea x = predecesor de x 
    mueva datos de D al sucesor de s

Las diferencias entre los algoritmos radican principalmente en cómo localizan los ciclos, cómo encuentran las direcciones iniciales en cada ciclo y cómo garantizan que cada ciclo se mueva exactamente una vez. Normalmente, como se ha comentado anteriormente, los ciclos se mueven en pares, ya que s y MN −1− s están en ciclos de la misma longitud (posiblemente el mismo ciclo). A veces, se utiliza una pequeña matriz de referencia, normalmente de longitud M + N (p. ej., Brenner, 1973; Cate y Twigg, 1977) para realizar un seguimiento de un subconjunto de ubicaciones en la matriz que se han visitado, para acelerar el algoritmo.

Para determinar si un ciclo dado ya se ha movido, el esquema más simple sería usar un almacenamiento auxiliar O ( MN ), un bit por elemento, para indicar si un elemento dado se ha movido. Para usar solo un almacenamiento auxiliar O ( M + N ) o incluso O (log  MN ) , se requieren algoritmos más complejos, y los algoritmos conocidos tienen un costo computacional linealítico en el peor de los casos de O ( MN  log  MN ) en el mejor de los casos, como lo demostró por primera vez Knuth (Fich et al. , 1995; Gustavson y Swirszcz, 2007).

Estos algoritmos están diseñados para mover cada elemento de datos exactamente una vez. Sin embargo, también implican una cantidad considerable de aritmética para calcular los ciclos y requieren accesos a la memoria en gran medida no consecutivos, ya que los elementos adyacentes de los ciclos difieren en factores multiplicativos de N , como se explicó anteriormente.

Se han diseñado varios algoritmos para lograr una mayor localidad de memoria a costa de un mayor movimiento de datos, así como de unos requisitos de almacenamiento ligeramente mayores. Es decir, pueden mover cada elemento de datos más de una vez, pero implican un mayor acceso consecutivo a memoria (mayor localidad espacial), lo que puede mejorar el rendimiento en las CPU modernas que dependen de cachés, así como en arquitecturas SIMD optimizadas para procesar bloques de datos consecutivos. El contexto más antiguo en el que parece haberse estudiado la localidad espacial de transposición es para operaciones fuera del núcleo (por Alltop, 1975), donde la matriz es demasiado grande para caber en la memoria principal (" core ").

Por ejemplo, si d = gcd ( N , M ) no es pequeño, se puede realizar la transposición usando una pequeña cantidad ( NM / d ) de almacenamiento adicional, con como máximo tres pasadas sobre la matriz (Alltop, 1975; Dow, 1995). Dos de las pasadas involucran una secuencia de transposiciones pequeñas separadas (que se pueden realizar de manera eficiente fuera de lugar usando un búfer pequeño) y una involucra una transposición cuadrada de bloques en el lugar d × d (que es eficiente ya que los bloques que se mueven son grandes y consecutivos, y los ciclos tienen una longitud de como máximo 2). Esto se simplifica aún más si N es un múltiplo de M (o viceversa), ya que solo se requiere una de las dos pasadas fuera de lugar. ${\displaystyle NM/d^{2}}$

Catanzaro et al. (2014) describieron otro algoritmo para dimensiones no coprimas , que implica múltiples transposiciones subsidiarias. Para el caso en que | N  −  M | es pequeño, Dow (1995) describe otro algoritmo que requiere | N  −  M | ⋅ min( N , M ) de almacenamiento adicional, que implica una transposición cuadrada min( N ,  M ) ⋅ min( N ,  M ) precedida o seguida por una pequeña transposición fuera de lugar. Frigo y Johnson (2005) describen la adaptación de estos algoritmos para utilizar técnicas ajenas a la caché para CPU de propósito general que dependen de líneas de caché para explotar la localidad espacial.

El trabajo sobre transposición de matrices fuera del núcleo, donde la matriz no cabe en la memoria principal y debe almacenarse en gran parte en un disco duro , se ha centrado principalmente en el caso de matriz cuadrada N = M , con algunas excepciones (por ejemplo, Alltop, 1975). Se pueden encontrar revisiones de algoritmos fuera del núcleo, especialmente aplicados a la computación paralela , en, por ejemplo, Suh y Prasanna (2002) y Krishnamoorth et al. (2004).

• OFFT: transposición recursiva en bloque de matrices cuadradas, en Fortran
• Jason Stratos Papadopoulos, transposición bloqueada en el lugar de matrices cuadradas, en C , grupo de noticias sci.math.num-analysis (7 de abril de 1998).
• Consulte los enlaces de "Código fuente" en la sección de referencias más arriba, para obtener código adicional para realizar transposiciones en el lugar de matrices cuadradas y no cuadradas.
• libmarshal Transposición bloqueada en el lugar de matrices rectangulares para las GPU.