インプレース行列転置

インプレース行列転置（ in-situ行列転置とも呼ばれる）は、N × M行列をコンピュータメモリ内でインプレース転置する問題であり、理想的にはO ( 1 )（制限付き）の追加ストレージ、または最大でもNMよりはるかに少ない追加ストレージを使用します。通常、行列は行優先または列優先の順序（つまり、連続する行または列がそれぞれ連続して配置されている）で格納されていると想定されます。

インプレース転置 (in-situ 転置) の実行は、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}}.}$

行優先形式では、これはコンピュータのメモリにシーケンス (11、12、13、14、21、22、23、24) として保存されます。つまり、2 つの行が連続して保存されます。これを転置すると、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 までの記憶場所の番号を付ける場合、この順列は 4 つのサイクルで構成されます。

(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行列が 0 ベースのインデックスを使用して行優先順に格納されていると仮定します。つまり、n = 0,..., N −1 およびm = 0,..., M −1 の ( n , m ) 要素は、アドレス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{otherwise}}.\end{cases}}}$

(これは、 N × M転置の逆はM × N転置であるという事実の結果に過ぎませんが、 P ^{−1 を}Pと合成すると恒等行列になることを明示的に示すことも簡単です。)

Cate & Twigg (1977) によって証明されたように、順列の固定点(長さ 1 のサイクル) の数は正確に1 + gcd( N −1, M −1)であり、ここで gcd は最大公約数です。たとえば、N = Mの場合、固定点の数は単にN (行列の対角線) です。一方、 N  − 1とM  − 1が互いに素である場合、固定点は行列の左上隅と右下隅の 2 つだけです。

任意の長さk >1のサイクル数は次のように与えられます（Cate & Twigg、1977）。

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

ここで、 μ はメビウス関数であり、その和はkの約数 dについて求められます。

さらに、 a =1（つまり、行列の最初の行の2番目の要素）を含むサイクルは常に最大長Lのサイクルであり、他のすべてのサイクルの長さkはLの約数でなければなりません （Cate＆Twigg、1977）。

与えられたサイクルCでは、すべての要素は同じ最大公約数を持ちます。 ${\displaystyle x\in C}$${\displaystyle d=\gcd(x,MN-1)}$

この定理は、順列のサイクルの検索に役立ちます。効率的な検索では、MN −1 の約数の倍数のみを調べることができるためです (Brenner、1973)。

Laflin & Brebner (1970) は、サイクルはペアで現れることが多いことを指摘しました。これは、サイクルのペアを一度に並べ替えるいくつかのアルゴリズムによって利用されます。特に、s を長さkのサイクルCの最小要素とします。したがって、MN −1− sも長さkのサイクルの要素です(同じサイクルである可能性もあります)。

以下は、インプレース行列転置を実行するために公開されているアルゴリズムを簡単にまとめたものです。 これらのアルゴリズムの一部を実装したソース コードは、以下の参考文献に記載されています。

行列を物理的に転置するのは計算コストが高いため、メモリ内の値を移動する代わりに、アクセスパスを転置することができます。イテレータのアクセスパスを交換するだけでよいため、CPUアクセスでこの操作を実行するのは簡単ですが、^[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 が2 の累乗の場合 (連想性が制限されたCPU キャッシュでのキャッシュラインの競合のため)、パフォーマンスが低下する可能性があります。その理由は、内側のループでmが増分されると、 A ( n , m ) またはA ( m , n ) に対応するメモリ アドレスがメモリ内でNずつ不連続にジャンプするためです (配列が列優先形式か行優先形式かによって異なります)。つまり、このアルゴリズムは参照の局所性を活用しません。

キャッシュの利用率を向上させる 1 つの解決策は、キャッシュ ライン サイズで指定されたブロックで、一度に複数の数値を操作するようにアルゴリズムを「ブロック」することです。残念ながら、これはアルゴリズムがキャッシュ ラインのサイズに依存する (「キャッシュ対応」) ことを意味し、複数レベルのキャッシュを備えた最新のコンピューターでは、マシンに依存する複数レベルのブロック化が必要になります。代わりに、再帰アルゴリズムによってパフォーマンスが向上することが示唆されています (Frigo他、1999)。再帰アルゴリズムとは、マトリックスをほぼ同じサイズの 4 つのサブマトリックスに分割し、2 つのサブマトリックスを対角線に沿って再帰的に転置し、2 つのサブマトリックスを対角線の上と下に転置して交換するものです。( Nが十分に小さい場合、上記の単純なアルゴリズムが基本ケースとして使用されます。単純にN =1 まで再帰すると、関数呼び出しのオーバーヘッドが過度に大きくなるためです。) これは、キャッシュ ライン サイズを明示的なパラメーターとして指定せずにキャッシュ ラインを活用できるという意味で、キャッシュを無視するアルゴリズムです。

非正方行列の場合、アルゴリズムはより複雑になります。1980 年以前のアルゴリズムの多くは、「サイクルに従う」アルゴリズムとして説明できます。つまり、サイクルをループし、サイクル内のある場所から次の場所へデータを移動します。擬似コード形式では、次のようになります。

順列のサイクルCの1より大きい長さごとにCの
    開始アドレスs
    を選択するD = sのデータとするサイクル内の
    x = sの前の
        データとする x ≠ s の場合、x から x の後続のデータへ移動する
    x = xの前のデータとするDから
        sの後続の
    データへ移動する 

アルゴリズム間の違いは、主にサイクルの配置方法、各サイクルの開始アドレスの見つけ方、各サイクルが正確に 1 回移動されるようにする方法にあります。通常、上で説明したように、sとMN −1− s は同じ長さのサイクル (同じサイクルの可能性もあります) にあるため、サイクルはペアで移動されます。場合によっては、アルゴリズムを高速化するために、通常は長さM + Nの小さなスクラッチ配列(例: Brenner、1973 年、Cate と Twigg、1977 年) を使用して、配列内の訪問済みの場所のサブセットを追跡します。

特定のサイクルがすでに移動されているかどうかを判断するための最も単純な方式は、O ( MN ) の補助ストレージ（要素ごとに 1ビット）を使用して、特定の要素が移動されているかどうかを示すことです。O ( M + N ) またはO ( log  MN )の補助ストレージのみを使用するには、より複雑なアルゴリズムが必要であり、既知のアルゴリズムは、最悪の場合の線形計算コストがせいぜいO ( MN  log  MN )であり、これはKnuthによって最初に証明されました(Fichら、1995 年、Gustavson と Swirszcz、2007 年)。

このようなアルゴリズムは、各データ要素を正確に 1 回移動するように設計されています。ただし、サイクルを計算するためにかなりの量の演算が必要であり、前述のように、サイクルの隣接する要素がNの乗数だけ異なるため、非常に非連続的なメモリ アクセスが必要になります。

いくつかのアルゴリズムは、より大きなデータ移動と若干のより大きなストレージ要件を犠牲にして、より大きなメモリ局所性を実現するように設計されてきました。つまり、各データ要素を複数回移動するかもしれませんが、より多くの連続したメモリ アクセス (より大きな空間局所性) を伴うため、キャッシュに依存する最新の CPU や、連続したデータ ブロックの処理に最適化されたSIMDアーキテクチャでパフォーマンスが向上します。転置の空間局所性が研究されたと思われる最も古いコンテキストは、マトリックスが大きすぎてメイン メモリ (「コア」)に収まらないアウトオブコア操作 (Alltop、1975 年) です。

たとえば、d = gcd ( N , M ) が小さくない場合、少量 ( NM / d ) の追加ストレージを使用して、配列の最大 3 回のパスで転置を実行できます (Alltop、1975 年、Dow、1995 年)。パスのうち 2 つは、個別の小さな転置のシーケンス (小さなバッファーを使用して効率的にアウトオブプレースで実行できます) を伴い、もう 1 つはブロックのインプレースd × d平方転置を伴います(移動するブロックが大きく連続しており、サイクルの長さが最大 ​​2 であるため効率的です)。N が M の倍数である場合 (またはその逆)、2 つのアウトオブプレース パスのうち 1 つだけが必要になるため、これはさらに簡略化されます。 ${\displaystyle NM/d^{2}}$

複数の補助転置を伴う、互いに素でない次元のための別のアルゴリズムは、Catanzaro ら (2014) によって説明されました。 | N  −  M |が小さい場合、Dow (1995) は、 | N  −  M | ⋅ min( N、M )の追加ストレージを必要とする別のアルゴリズムを説明しています。このアルゴリズムでは、min( N、  M ) ⋅ min( N、  M )の平方転置の前または後に小さなアウトオブプレース転置が伴います。Frigo と Johnson (2005) は、空間的局所性を活用するためにキャッシュ ラインに依存する汎用 CPU でキャッシュ無視手法を使用するためにこれらのアルゴリズムを適応させたことを説明しています。

アウトオブコア行列転置に関する研究は、行列がメインメモリに収まらず、主にハードディスクに保存する必要がある場合、いくつかの例外を除いて、主にN = M正方行列の場合に焦点が当てられてきました(例: Alltop、1975)。特に並列コンピューティングに適用されるアウトオブコアアルゴリズムのレビューは、Suh & Prasanna (2002) や Krishnamoorth 他 (2004) に記載されています。

• OFFT - Fortran での正方行列の再帰ブロックインプレース転置
• Jason Stratos Papadopoulos、 Cにおける正方行列のブロック化されたインプレース転置、sci.math.num-analysisニュースグループ (1998 年 4 月 7 日)。
• 正方行列と非正方行列の両方のインプレース転置を実行するための追加コードについては、上記の参考セクションの「ソース コード」リンクを参照してください。
• libmarshal は、GPU 用の長方形行列のブロックされたインプレース転置を実行します。