Транспонирование матриц на месте

Транспонирование матрицы на месте , также называемое транспонированием матрицы на месте , представляет собой задачу транспонирования матрицы N × M на месте в памяти компьютера , в идеале с O (1) (ограниченным) дополнительным хранилищем или, самое большее, с дополнительным хранилищем, намного меньшим, чем NM . Обычно предполагается, что матрица хранится в строковом или столбцовом порядке (т. е. смежные строки или столбцы, соответственно, расположены последовательно).

Выполнение транспонирования на месте (транспонирование на месте) наиболее сложно, когда N ≠ M , т. е. для неквадратной (прямоугольной) матрицы, где оно включает в себя сложную перестановку элементов данных со многими циклами длиной больше 2. Напротив, для квадратной матрицы ( N = M ) все циклы имеют длину 1 или 2, и транспонирование может быть достигнуто простым циклом для обмена верхнего треугольника матрицы с нижним треугольником. Дополнительные осложнения возникают, если кто-то хочет максимизировать локальность памяти , чтобы улучшить использование строки кэша или работать вне ядра (когда матрица не помещается в основную память), поскольку транспонирование по своей сути подразумевает непоследовательные доступы к памяти.

Проблема неквадратной транспозиции на месте изучается по крайней мере с конца 1950-х годов, и известно несколько алгоритмов, включая несколько, которые пытаются оптимизировать локальность для кэша, внеядерной памяти или аналогичных контекстов, связанных с памятью.

На компьютере часто можно избежать явного транспонирования матрицы в памяти , просто обратившись к тем же данным в другом порядке. Например, библиотеки программного обеспечения для линейной алгебры , такие как BLAS , обычно предоставляют опции для указания того, что определенные матрицы должны интерпретироваться в транспонированном порядке, чтобы избежать перемещения данных.

Однако остается ряд обстоятельств, при которых необходимо или желательно физически переупорядочить матрицу в памяти в транспонированный порядок. Например, при матрице, хранящейся в порядке по строкам , строки матрицы являются смежными в памяти, а столбцы — несмежными. Если необходимо выполнить повторные операции над столбцами, например, в алгоритме быстрого преобразования Фурье (например, Frigo & Johnson, 2005), транспонирование матрицы в памяти (чтобы сделать столбцы смежными) может улучшить производительность за счет увеличения локальности памяти . Поскольку эти ситуации обычно совпадают со случаем очень больших матриц (которые превышают размер кэша), выполнение транспонирования на месте с минимальным дополнительным хранилищем становится желательным.

Кроме того, как чисто математическая задача, транспонирование на месте включает в себя ряд интересных головоломок теории чисел , которые решались в течение нескольких десятилетий.

Например, рассмотрим матрицу 2×4:

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

В формате row-major это будет храниться в памяти компьютера как последовательность (11, 12, 13, 14, 21, 22, 23, 24), т. е. две строки, хранящиеся последовательно. Если мы транспонируем это, то получим матрицу 4×2:

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

которая хранится в памяти компьютера в виде последовательности (11, 21, 12, 22, 13, 23, 14, 24).

Если пронумеровать ячейки хранения от 0 до 7 слева направо, то эта перестановка будет состоять из четырех циклов:

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

То есть значение в позиции 0 переходит в позицию 0 (цикл длиной 1, без перемещения данных). Далее значение в позиции 1 (в исходном хранилище: 11, 12 , 13, 14, 21, 22, 23, 24) переходит в позицию 2 (в транспонированном хранилище: 11, 21, 12, 22, 13 , 23 , 14, 24), в то время как значение в позиции 2 (11, 12, 13 , 14, 21, 22, 23, 24) переходит в позицию 4 (11, 21, 12, 22, 13 , 23, 14, 24), а позиция 4 (11, 12, 13, 14, 21 , 22, 23, 24) возвращается в позицию 1 (11, 21 , 12, 22, 13, 23, 14, 24). Аналогично для значений в позиции 7 и позициях (3 6 5).

Далее мы предполагаем, что матрица N × M хранится в строчном порядке с индексами, начинающимися с нуля. Это означает, что элемент ( n , m ) для n = 0,..., N −1 и m = 0,..., M −1 хранится по адресу a = Mn + m (плюс некоторое смещение в памяти, которое мы игнорируем). В транспонированной матрице M × N соответствующий элемент ( m , n ) хранится по адресу a' = Nm + n , снова в строчном порядке. Мы определяем перестановку транспозиции как функцию a' = P ( a ) такую, что:

${\displaystyle Nm+n=P(Mn+m)\,}$для всех${\displaystyle (n,m)\in [0,N-1]\times [0,M-1]\,.}$

Это определяет перестановку чисел .${\displaystyle a=0,\ldots ,MN-1}$

Оказывается, можно определить простые формулы для P и его обратных величин (Cate & Twigg, 1977). Первое:

${\displaystyle P(a)={\begin{cases}MN-1&{\text{if}}a=MN-1,\\Na{\bmod {(}}MN-1)&{\text{otherwise}},\end{cases}}}$

где «mod» — операция по модулю .

Во-вторых, обратная перестановка задается формулой:

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

(Это всего лишь следствие того факта, что инверсия транспонирования N × M является транспонированием M × N , хотя также легко явно показать, что P ⁻¹ в сочетании с P дает тождество.)

Как доказали Кейт и Твигг (1977), число неподвижных точек (циклов длины 1) перестановки равно точно 1 + gcd( N −1, M −1) , где gcd — наибольший общий делитель . Например, при N = M число неподвижных точек равно просто N (диагонали матрицы). С другой стороны, если N  − 1 и M  − 1 взаимно просты , то единственными двумя неподвижными точками являются верхний левый и нижний правый углы матрицы.

Число циклов любой длины k >1 определяется по формуле (Кейт и Твигг, 1977):

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

где μ — функция Мёбиуса , а сумма берется по делителям d числа k .

Более того, цикл, содержащий a = 1 (т. е. второй элемент первой строки матрицы), всегда является циклом максимальной длины L , а длины k всех остальных циклов должны быть делителями L (Кейт и Твигг, 1977).

Для данного цикла C каждый элемент имеет один и тот же наибольший общий делитель .${\displaystyle x\in C}$${\displaystyle d=\НОД(x,MN-1)}$

Эта теорема полезна при поиске циклов перестановки, поскольку эффективный поиск может рассматривать только кратные делителям MN −1 (Бреннер, 1973).

Лафлин и Бребнер (1970) указали, что циклы часто встречаются парами, что используется несколькими алгоритмами, которые переставляют пары циклов за раз. В частности, пусть s будет наименьшим элементом некоторого цикла C длины k . Из этого следует, что MN −1− s также является элементом цикла длины k (возможно, того же самого цикла).

Ниже кратко изложены опубликованные алгоритмы для выполнения транспонирования матриц на месте. Исходный код, реализующий некоторые из этих алгоритмов, можно найти в ссылках ниже.

Поскольку физическое транспонирование матрицы является вычислительно затратным, вместо перемещения значений в памяти можно транспонировать путь доступа. Выполнить эту операцию для доступа к ЦП тривиально, поскольку пути доступа итераторов должны быть просто заменены, ^[1] однако аппаратное ускорение может потребовать, чтобы это все еще было физически перестроено. ^[2]

Для квадратной матрицы N × N A _{n , m} = A ( n , m ), транспонирование на месте выполняется легко, поскольку все циклы имеют длину 1 (диагонали A _{n , n} ) или длину 2 (верхний треугольник меняется местами с нижним треугольником). Псевдокод для достижения этого (предполагая, что индексы массива начинаются с нуля ) выглядит следующим образом:

для n = 0 до N - 1 для m = n + 1 до N поменять местами A(n,m) и A(m,n)

Этот тип реализации, хотя и прост, может демонстрировать низкую производительность из-за плохого использования кэш-линий, особенно когда N является степенью двойки (из-за конфликтов кэш-линий в кэше ЦП с ограниченной ассоциативностью). Причина этого в том, что по мере увеличения m во внутреннем цикле адрес памяти, соответствующий A ( n , m ) или A ( m , n ), скачком переходит на N в памяти (в зависимости от того, находится ли массив в формате column-major или row-major соответственно). То есть алгоритм не использует локальность ссылки .

Одним из решений для улучшения использования кэша является «блокировка» алгоритма для работы с несколькими числами одновременно, в блоках, заданных размером кэш-строки; к сожалению, это означает, что алгоритм зависит от размера кэш-строки (он «осведомлен о кэше»), и на современном компьютере с несколькими уровнями кэша он требует нескольких уровней машинно-зависимой блокировки. Вместо этого было предложено (Frigo et al. , 1999), что более высокую производительность можно получить с помощью рекурсивного алгоритма: разделить матрицу на четыре подматрицы примерно одинакового размера, транспонировать две подматрицы по диагонали рекурсивно и транспонировать и менять местами две подматрицы выше и ниже диагонали. (Когда N достаточно мало, простой алгоритм выше используется в качестве базового случая, так как наивное повторение вплоть до N = 1 привело бы к чрезмерным накладным расходам на вызов функций.) Это алгоритм, игнорирующий кэш , в том смысле, что он может использовать кэш-строку без явного параметра размера кэш-строки.

Для неквадратных матриц алгоритмы более сложные. Многие алгоритмы до 1980 года можно было бы описать как алгоритмы "следования циклам". То есть они циклически обходят циклы, перемещая данные из одного места в следующее в цикле. В форме псевдокода:

для каждой длины>1 цикла C перестановки выберите начальный адрес s в C пусть D = данные в s пусть x = предшественник s в цикле пока  x ≠ s переместите данные из x в последующий элемент x пусть x = предшественник x переместите данные из D в последующий элемент s

Различия между алгоритмами в основном заключаются в том, как они находят циклы, как они находят начальные адреса в каждом цикле и как они обеспечивают, чтобы каждый цикл перемещался ровно один раз. Обычно, как обсуждалось выше, циклы перемещаются парами, поскольку s и MN −1− s находятся в циклах одинаковой длины (возможно, в одном и том же цикле). Иногда небольшой массив скретчей, обычно длиной M + N (например, Brenner, 1973; Cate & Twigg, 1977), используется для отслеживания подмножества посещенных мест в массиве, чтобы ускорить алгоритм.

Чтобы определить, был ли данный цикл уже перемещен, простейшей схемой было бы использование O ( MN ) вспомогательной памяти, по одному биту на элемент, чтобы указать, был ли перемещен данный элемент. Чтобы использовать только O ( M + N ) или даже O (log  MN ) вспомогательной памяти, требуются более сложные алгоритмы, а известные алгоритмы имеют худшую линейную вычислительную стоимость O ( MN  log  MN ) в лучшем случае, как впервые доказал Кнут (Fich et al. , 1995; Gustavson & Swirszcz, 2007).

Такие алгоритмы предназначены для перемещения каждого элемента данных ровно один раз. Однако они также требуют значительного количества арифметических операций для вычисления циклов и требуют большого количества непоследовательных обращений к памяти, поскольку соседние элементы циклов отличаются на мультипликативные множители N , как обсуждалось выше.

Было разработано несколько алгоритмов для достижения большей локальности памяти за счет большего перемещения данных, а также немного больших требований к хранению. То есть, они могут перемещать каждый элемент данных более одного раза, но они подразумевают более последовательный доступ к памяти (большую пространственную локальность), что может улучшить производительность на современных процессорах, которые полагаются на кэши, а также на архитектурах SIMD , оптимизированных для обработки последовательных блоков данных. Самый старый контекст, в котором, по-видимому, изучалась пространственная локальность транспозиции, — это работа вне ядра (Alltop, 1975), когда матрица слишком велика, чтобы поместиться в основную память (« ядро »).

Например, если d = gcd ( N , M ) не мало, можно выполнить транспозицию, используя небольшой объем ( NM / d ) дополнительной памяти, максимум с тремя проходами по массиву (Alltop, 1975; Dow, 1995). Два прохода включают последовательность отдельных небольших транспозиций (которые могут быть эффективно выполнены вне места с использованием небольшого буфера), а один включает квадратную транспозицию блоков d × d на месте (что эффективно, поскольку перемещаемые блоки большие и последовательны, а циклы имеют длину не более 2). Это еще больше упрощается, если N кратно M (или наоборот), поскольку требуется только один из двух проходов вне места.${\displaystyle NM/d^{2}}$

Другой алгоритм для невзаимно простых измерений, включающий множественные вспомогательные транспозиции, был описан Катанзаро и др. (2014). Для случая, когда | N  −  M | мало, Доу (1995) описывает другой алгоритм, требующий | N  −  M | ⋅ min( N , M ) дополнительной памяти, включающий квадратную транспозицию min( N ,  M ) ⋅ min( N ,  M ), которой предшествует или за которой следует небольшая транспозиция не на своем месте. Фриго и Джонсон (2005) описывают адаптацию этих алгоритмов для использования методов, забывающих о кэше, для универсальных ЦП, полагающихся на строки кэша для использования пространственной локальности.

Работа по транспонированию матриц вне ядра, когда матрица не помещается в основную память и должна храниться в основном на жестком диске , была сосредоточена в основном на случае квадратных матриц N = M , за некоторыми исключениями (например, Alltop, 1975). Обзоры алгоритмов вне ядра, особенно применяемых к параллельным вычислениям , можно найти, например, в Suh & Prasanna (2002) и Krishnamoorth et al. (2004).

• OFFT - рекурсивное блочное транспонирование квадратных матриц на Фортране
• Джейсон Стратос Пападопулос, заблокированное транспонирование квадратных матриц на месте, на языке C , группа новостей sci.math.num-analysis (7 апреля 1998 г.).
• Дополнительный код для выполнения транспонирования на месте как квадратных, так и неквадратных матриц см. в ссылках «Исходный код» в разделе ссылок выше.
• libmarshal Заблокировано транспонирование прямоугольных матриц на месте для графических процессоров.