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

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

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

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

Фон

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

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

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

Пример

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

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

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

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

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

Позиция	0	1	2	3	4	5	6	7
Исходное хранилище	11	12	13	14	21	22	23	24
Транспонированное хранилище	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 ) такую, что:

Nm+n=P(Mn+m)\,

для всех

(n,m)\in [0,N-1]\times [0,M-1]\,.

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

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

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

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

Доказательство
If 0 ≤ a = Mn + m < MN − 1, then Na mod (MN−1) = MN n + Nm mod (MN − 1) = n + Nm. ^{[ProofNote 1]}^{[ProofNote 2]}

Во-вторых, обратная перестановка определяется выражением:

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

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

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

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

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

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

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

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

Доказательство (Бреннер, 1973)
Let s be the smallest element of the cycle, and $d=\gcd(s,MN-1)$ . From the definition of the permutation P above, every other element x of the cycle is obtained by repeatedly multiplying s by N modulo MN−1, and therefore every other element is divisible by d. But, since N and $MN - 1$ are coprime, x cannot be divisible by any factor of $MN - 1$ larger than d, and hence $d=\gcd(x,MN-1)$ .

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

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

Доказательство по определению P выше.
The length k of the cycle containing s is the smallest k > 0 such that $sN^{k}=s{\bmod {(}}MN-1)$ . Clearly, this is the same as the smallest k>0 such that $(-s)N^{k}=-s{\bmod {(}}MN-1)$ , since we are just multiplying both sides by −1, and $MN-1-s=-s{\bmod {(}}MN-1)$ .

Примечание доказательств
^ MN x mod (MN−1) = (MN − 1) x + x mod (MN−1) = x for 0 ≤ x < MN − 1. ^ The first (a = 0) and last (a = MN−1) elements are always left invariant under transposition.

Алгоритмы

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

Аксессор транспонировать

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

Квадратные матрицы

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

for n = 0 to N - 1
    for m = n + 1 to N
        swap A(n,m) with A(m,n)

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

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

Неквадратные матрицы: по циклам

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

for each length>1 cycle C of the permutation
    pick a starting address s in C
    let D = data at s
    let x = predecessor of s in the cycle
    while x ≠ s
        move data from x to successor of x
        let x = predecessor of x
    move data from D to successor of s

Различия между алгоритмами заключаются главным образом в том, как они находят циклы, как они находят начальные адреса в каждом цикле и как они гарантируют, что каждый цикл перемещается ровно один раз. Обычно, как обсуждалось выше, циклы перемещаются парами, поскольку s и MN -1- находятся в циклах одинаковой длины (возможно, в одном и том же цикле). Иногда небольшой временной массив, обычно длиной 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 в квадрате $NM/d^{2}$ блоков (что эффективно, поскольку перемещаемые блоки большие и последовательные, а длина циклов не превышает 2). Это еще больше упрощается, если N кратно M (или наоборот), поскольку требуется только один из двух неуместных проходов.

Другой алгоритм для невзаимно простых измерений, включающий множественные вспомогательные транспозиции, был описан Catanzaro et al. (2014). Для случая, когда $| Н - М |$ мал, Доу (1995) описывает другой алгоритм, требующий $| Н - М | \cdot min(N, M)$ дополнительное хранилище, включающее $квадратное транспонирование min(N, M) \cdot min(N, M)$ , которому предшествует или за которым следует небольшое транспонирование неуместно. Фриго и Джонсон (2005) описывают адаптацию этих алгоритмов для использования методов без учета кэша для процессоров общего назначения, полагающихся на строки кэша для использования пространственной локальности.

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

Ссылки

^ «numpy.swapaxes — Руководство по NumPy v1.15» . docs.scipy.org . Проверено 22 января 2019 г.
^ Харрис, Марк (18 февраля 2013 г.). «Эффективное транспонирование матрицы в CUDA C/C++» . Блог разработчиков NVIDIA .

П.Ф. Уиндли, «Транспонирование матриц в цифровом компьютере», Computer Journal 2 , стр. 47-48 (1959).
Г. Палл и Э. Зейден, «Проблема абелевых групп с применением к транспонированию матрицы на электронном компьютере», Math. Комп. 14 , с. 189–192 (1960).
Дж. Бутройд, «Алгоритм 302: Транспонирование векторного хранимого массива», Транзакции ACM в математическом программном обеспечении 10 (5), стр. 292-293 (1967). дои : 10.1145/363282.363304
Сьюзен Лафлин и М.А. Бребнер, «Алгоритм 380: транспонирование прямоугольной матрицы на месте», Транзакции ACM в математическом программном обеспечении 13 (5), стр. 324-326 (1970). doi : 10.1145/362349.362368 Исходный код .
Норман Бреннер, «Алгоритм 467: транспонирование матрицы на месте», Транзакции ACM в математическом программном обеспечении 16 (11), стр. 692-694 (1973). doi : 10.1145/355611.362542 Исходный код .
WO Alltop, «Компьютерный алгоритм транспонирования неквадратных матриц», IEEE Trans. Вычислить. 24 (10), с. 1038-1040 (1975).
Эско Г. Кейт и Дэвид В. Твигг, «Алгоритм 513: Анализ транспозиции на месте», Транзакции ACM в математическом программном обеспечении 3 (1), стр. 104-110 (1977). doi : 10.1145/355719.355729 Исходный код .
Брайан Катандзаро, Александр Келлер и Майкл Гарланд, «Декомпозиция для транспонирования матрицы на месте», Труды 19-го симпозиума ACM SIGPLAN по принципам и практике параллельного программирования (PPoPP '14), стр. 193–206 (2014). дои : 10.1145/2555243.2555253
Мюррей Доу, «Транспонирование матрицы на векторном компьютере», Parallel Computing 21 (12), стр. 1997-2005 (1995).
Дональд Э. Кнут, Искусство компьютерного программирования, том 1: Фундаментальные алгоритмы , третье издание, раздел 1.3.3, упражнение 12 (Аддисон-Уэсли: Нью-Йорк, 1997).
М. Фриго, К. Э. Лейзерсон, Х. Прокоп и С. Рамачандран, «Алгоритмы, не обращающие внимания на кэш», в материалах 40-го симпозиума IEEE по основам компьютерных наук (FOCS 99), стр. 285-297 (1999). два : 10.1109/SFCS.1999.814600
Дж. Су и В. К. Прасанна, «Эффективный алгоритм транспонирования матриц вне ядра», IEEE Trans. Компьютеры 51 (4), с. 420-438 (2002). дои : 10.1109/12.995452
С. Кришнамурти, Г. Баумгартнер, Д. Кочорва, К.-К. Лам и П. Садаяппан, « Эффективное параллельное транспонирование матрицы вне ядра », Международный журнал высокопроизводительных вычислений и сетей 2 (2-4), стр. 110-119 (2004).
М. Фриго и С.Г. Джонсон, « Проектирование и реализация FFTW3 », Труды IEEE 93 (2), 216–231 (2005). Исходный код библиотеки FFTW , которая включает в себя оптимизированные последовательные и параллельные квадратичные и неквадратные транспонирования в дополнение к БПФ .
Фейт Э. Фич , Дж. Ян Манро и Патрисио В. Поблете, «Перестановка на месте», SIAM Journal on Computing 24 (2), стр. 266-278 (1995).
Фред Г. Густавсон и Тадеуш Свирщ, «Транспонирование прямоугольных матриц на месте», Конспекты лекций по информатике 4699 , стр. 560-569 (2007), из материалов семинара 2006 года по новейшим достижениям в области научных и параллельных вычислений (PARA 2006) (Умео, Швеция, июнь 2006 г.).
Слоан, Нью-Джерси (ред.). «Последовательность A093055 (Количество неодноэлементных циклов при транспонировании на месте прямоугольной матрицы j X k)» . Электронная энциклопедия целочисленных последовательностей . Фонд ОЭИС.
Слоан, Нью-Джерси (ред.). «Последовательность A093056 (Длина самого длинного цикла в транспозиции in-situ прямоугольной матрицы j X k)» . Электронная энциклопедия целочисленных последовательностей . Фонд ОЭИС.
Слоан, Нью-Джерси (ред.). «Последовательность A093057 (Количество матричных элементов, остающихся в фиксированном положении при транспонировании на месте прямоугольной матрицы j X k)» . Электронная энциклопедия целочисленных последовательностей . Фонд ОЭИС.

Внешние ссылки

Исходный код

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

[1] MN x mod (MN−1) = (MN − 1) x + x mod (MN−1) = x for 0 ≤ x < MN − 1.

[2] The first (a = 0) and last (a = MN−1) elements are always left invariant under transposition.

[3] «numpy.swapaxes — Руководство по NumPy v1.15» . docs.scipy.org . Проверено 22 января 2019 г.

[4] Харрис, Марк (18 февраля 2013 г.). «Эффективное транспонирование матрицы в CUDA C/C++» . Блог разработчиков NVIDIA .

[ProofNote 1]

[ProofNote 2]

[1]

[2]