遅い - なぜこの単純なRループを「ベクトル化」すると結果が変わるのですか?




r リスト (3)

おそらく非常にばかな質問です。

私は次のループを "ベクトル化"しようとしています:

set.seed(0)
x <- round(runif(10), 2)
# [1] 0.90 0.27 0.37 0.57 0.91 0.20 0.90 0.94 0.66 0.63
sig <- sample.int(10)
# [1]  1  2  9  5  3  4  8  6  7 10
for (i in seq_along(sig)) x[i] <- x[sig[i]]
x
# [1] 0.90 0.27 0.66 0.91 0.66 0.91 0.94 0.91 0.94 0.63

私はそれが単にx[sig]だと思うが、結果は一致しない。

set.seed(0)
x <- round(runif(10), 2)
x[] <- x[sig]
x
# [1] 0.90 0.27 0.66 0.91 0.37 0.57 0.94 0.20 0.90 0.63

どうしましたか?

リマーク

明らかに出力から、 forループとx[sig]が異なることが分かります。 後者の意味は明白です:順列、したがって多くの人々はループがちょうど間違ったことをやっていると信じがちです。 しかし、それほど確実ではない。 それはいくつかの明確に定義された動的プロセスである可能性があります。 このQ&Aの目的は、どちらが正しいかを判断するのではなく、なぜそれらが同等でないのかを説明することです。 うまくいけば、「ベクトル化」を理解するための堅実なケーススタディが提供されることを望みます。


準備し始める

ウォームアップとして、より簡単な2つの例を考えてみましょう。

## example 1
x <- 1:11
for (i in 1:10) x[i] <- x[i + 1]
x
# [1]  2  3  4  5  6  7  8  9 10 11 11

x <- 1:11
x[1:10] <- x[2:11]
x
# [1]  2  3  4  5  6  7  8  9 10 11 11

## example 2
x <- 1:11
for (i in 1:10) x[i + 1] <- x[i]
x
# [1] 1 1 1 1 1 1 1 1 1 1 1

x <- 1:11
x[2:11] <- x[1:10]
x
# [1]  1  1  2  3  4  5  6  7  8  9 10

第1の例では「ベクトル化」は成功したが、第2の例では成功しなかった。 どうして?

ここで慎重な分析です。 「ベクトル化」はループアンローリングによって開始され、次に複数の命令を並列に実行します。 ループが "ベクトル化"できるかどうかは、ループによって運ばれるデータ依存性によって決まります。

例1のループを展開すると、

x[1]  <- x[2]
x[2]  <- x[3]
x[3]  <- x[4]
x[4]  <- x[5]
x[5]  <- x[6]
x[6]  <- x[7]
x[7]  <- x[8]
x[8]  <- x[9]
x[9]  <- x[10]
x[10] <- x[11]

これらの命令を1つずつ実行して同時に実行すると、同じ結果が得られます。 したがって、このループは「ベクトル化」することができます。

例2のループは次のとおりです。

x[2]  <- x[1]
x[3]  <- x[2]
x[4]  <- x[3]
x[5]  <- x[4]
x[6]  <- x[5]
x[7]  <- x[6]
x[8]  <- x[7]
x[9]  <- x[8]
x[10] <- x[9]
x[11] <- x[10]

残念ながら、これらの命令を1つずつ実行して同時に実行すると、同じ結果が得られません。 たとえば、1つずつ実行する場合、 x[2]が第1命令で変更され、この変更された値が第2命令のx[3]に渡されます。 したがって、 x[3]x[1]と同じ値を持ちます。 しかし、並列実行では、 x[3]x[2]等しくなります。 その結果、このループを「ベクトル化」することはできません。

「ベクトル化」理論では、

  • 例1は、データに「write-after-read」依存関係を持ちますx[i]は読み込み後に変更されます。
  • 例2は、データの中に「read-after-write」依存関係がありますx[i]は、変更後に読み込まれます。

"write-after-read"データ依存性を持つループは "ベクトル化"することができますが、 "read-after-write"データ依存性を持つループはできません。

深さ

おそらく多くの人々が今混乱しているでしょう。 「ベクトル化」は「並列処理」ですか?

はい。 ハイスピードコンピューティングのためにどのような並列処理コンピュータを設計するのか、人々が疑問を呈した1960年代、Flynnはデザインアイデアを4つのタイプに分類しました。 カテゴリー「SIMD」(単一命令、複数データ)は、「ベクトル化」され、SIMDキャビタビリティを有するコンピュータは、「ベクトルプロセッサ」または「アレイプロセッサ」と呼ばれる。

1960年代には多くのプログラミング言語がありませんでした。 人々は、CPUレジスタを直接プログラムするためにアセンブリ(コンパイラが発明されたときはFORTRAN)を書きました。 「SIMD」コンピュータは、1つの命令でベクトルレジスタに複数のデータをロードし、同時にそれらのデータに対して同じ演算を実行することができます。 したがって、データ処理は実際に並行して行われます。 もう一度例1を考えてみましょう。 ベクトルレジスタが2つのベクトル要素を保持できるとすると、スカラー処理のように10回の反復ではなくベクトル処理を使用して5回の反復でループを実行できます。

reg <- x[2:3]  ## load vector register
x[1:2] <- reg  ## store vector register
-------------
reg <- x[4:5]  ## load vector register
x[3:4] <- reg  ## store vector register
-------------
reg <- x[6:7]  ## load vector register
x[5:6] <- reg  ## store vector register
-------------
reg <- x[8:9]  ## load vector register
x[7:8] <- reg  ## store vector register
-------------
reg <- x[10:11] ## load vector register
x[9:10] <- reg  ## store vector register

今日、 Rのような多くのプログラミング言語があります。 「ベクトル化」はもはや「SIMD」を明白に指していません。 RはCPUレジスタをプログラムできる言語ではありません。 Rの「ベクトル化」は、単に「SIMD」に類似しています。 以前のQ&Aでは、 「ベクトル化」という用語は、さまざまな状況で異なることを意味していますか? 私はこれを説明しようとしました。 次の図は、この類推がどのように行われているかを示しています。

single (assembly) instruction    -> single R instruction
CPU vector registers             -> temporary vectors
parallel processing in registers -> C/C++/FORTRAN loops with temporary vectors

したがって、例1のループのRの「ベクトル化」は次のようなものです

## the C-level loop is implemented by function "["
tmp <- x[2:11]  ## load data into a temporary vector
x[1:10] <- tmp  ## fill temporary vector into x

私たちがやっているほとんどの時間

x[1:10] <- x[2:10]

一時的なベクトルを変数に明示的に代入する必要はありません。 作成された一時的なメモリブロックは、R変数によって指されないため、ガベージコレクションの対象となります。

完全な画像

上記において、「ベクトル化」は最も簡単な例では導入されていません。 非常に頻繁に、「ベクトル化」が

a[1] <- b[1] + c[1]
a[2] <- b[2] + c[2]
a[3] <- b[3] + c[3]
a[4] <- b[4] + c[4]

abおよびcはメモリ内でエイリアシングされない。すなわち、ベクトルabおよびc格納するメモリブロックは重ならない。 これは理想的なケースです。メモリエイリアシングがデータ依存性を持たないためです。

"データ依存性"とは別に、 "ベクトル化"の "if ... else ..."を扱う "制御依存"もあります。 しかし、時間と空間の理由から私はこの問題について詳述しません。

質問の例に戻る

今、問題のループを調査する時です。

set.seed(0)
x <- round(runif(10), 2)
sig <- sample.int(10)
# [1]  1  2  9  5  3  4  8  6  7 10
for (i in seq_along(sig)) x[i] <- x[sig[i]]

ループをアンロールすると、

x[1]  <- x[1]
x[2]  <- x[2]
x[3]  <- x[9]   ## 3rd instruction
x[4]  <- x[5]
x[5]  <- x[3]   ## 5th instruction
x[6]  <- x[4]
x[7]  <- x[8]
x[8]  <- x[6]
x[9]  <- x[7]
x[10] <- x[10]

第3命令と第5命令の間には "リード・アフター・ライト"データ依存関係があるため、ループを「ベクトル化」することはできません(参考1を参照)。

さて、 x[] <- x[sig]は何をするのですか? 最初に、一時的なベクトルを明示的に書き出しましょう:

tmp <- x[sig]
x[] <- tmp

"["は2回呼び出されるので、実際にはこの "ベクトル化された"コードの後ろに2つのCレベルのループがあります:

tmp[1]  <- x[1]
tmp[2]  <- x[2]
tmp[3]  <- x[9]
tmp[4]  <- x[5]
tmp[5]  <- x[3]
tmp[6]  <- x[4]
tmp[7]  <- x[8]
tmp[8]  <- x[6]
tmp[9]  <- x[7]
tmp[10] <- x[10]

x[1]  <- tmp[1]
x[2]  <- tmp[2]
x[3]  <- tmp[3]
x[4]  <- tmp[4]
x[5]  <- tmp[5]
x[6]  <- tmp[6]
x[7]  <- tmp[7]
x[8]  <- tmp[8]
x[9]  <- tmp[9]
x[10] <- tmp[10]

だからx[] <- x[sig]

for (i in 1:10) tmp[i] <- x[sig[i]]
for (i in 1:10) x[i] <- tmp[i]
rm(tmp); gc()

この質問に与えられた元のループではありません。

備考1

Rcppでループを実装することを「ベクトル化」と見なした場合、それをそのままにしておきます。 しかし、 "SIMD"を使ってC / C ++ループをさらに「ベクトル化」する機会はありません。

備考2

このQ&Aは、 このQ&Aによって動機付けられています。 OPはもともとループを提示

for (i in 1:num) {
  for (j in 1:num) {
    mat[i, j] <- mat[i, mat[j, "rm"]]
  }
}

それを「ベクトル化する」ことは魅力的です

mat[1:num, 1:num] <- mat[1:num, mat[1:num, "rm"]]

潜在的に間違っています。 後でOPがループを

for (i in 1:num) {
  for (j in 1:num) {
    mat[i, j] <- mat[i, 1 + num + mat[j, "rm"]]
  }
}

置き換えられる列が最初のnum列であり、ルックアップされる列が最初のnum列の後にあるため、メモリのエイリアシングの問題が解消されます。

備考3

問題のループがx "インプレース"修正を行っているかどうかについていくつかのコメントがありました。 はい、そうです。 tracememを使うことができます:

set.seed(0)
x <- round(runif(10), 2)
sig <- sample.int(10)
tracemem(x)
#[1] "<0x28f7340>"
for (i in seq_along(sig)) x[i] <- x[sig[i]]
tracemem(x)
#[1] "<0x28f7340>"

私のRセッションは、 xアドレス<0x28f7340>で示されるメモリブロックを割り当てており、コードを実行するときに別の値が表示されることがあります。 しかし、 tracememの出力はループの後も変化しません。つまり、 xのコピーは作成されません。 したがって、ループは実際に余分なメモリを使用せずに "インプレース"変更を行っています。

しかし、このループでは、インプレース置換が行われていません。 「インプレース」置換はより複雑な操作です。 x要素をループに沿ってスワップするだけでなく、 sig要素も入れ替える必要があります(最後にsig1:10なります)。


これはメモリブロックのエイリアシング(私が以前に遭遇したことのない用語)とは関係ありません。 特定の順列の例をとり、C言語またはアセンブリ言語(またはそれに関係なく)の言語レベルでの実装に関係なく発生する割り当てを実行します。 シーケンシャルforループがどのように振る舞い、どのような "真"の順列( x[sig]で何が得られるか)がどのように起こるかは、本質的なものです。

sample(10)
 [1]  3  7  1  5  6  9 10  8  4  2

value at 1 goes to 3, and now there are two of those values
value at 2 goes to 7, and now there are two of those values
value at 3 (which was at 1) now goes back to 1 but the values remain unchanged

...これは続けることができますが、これは通常どのように "真の"順列ではなく、非常に珍しくも値の完全な再分配をもたらす方法を示しています。 私は、完全に順序付けられた順列(そのうちの1つ、すなわち10:1だと思う)だけがユニークであった新しいxのセットをもたらすことができると推測しています。

replicate( 100, {x <- round(runif(10), 2); 
                  sig <- sample.int(10); 
                  for (i in seq_along(sig)){ x[i] <- x[sig[i]]}; 
                  sum(duplicated(x)) } )
 #[1] 4 4 4 5 5 5 4 5 6 5 5 5 4 5 5 6 3 4 2 5 4 4 4 4 3 5 3 5 4 5 5 5 5 5 5 5 4 5 5 5 5 4
 #[43] 5 3 4 6 6 6 3 4 5 3 5 4 6 4 5 5 6 4 4 4 5 3 4 3 4 4 3 6 4 7 6 5 6 6 5 4 7 5 6 3 6 4
 #[85] 8 4 5 5 4 5 5 5 4 5 5 4 4 5 4 5

私は重複カウントの分布が大規模なシリーズになる可能性があるのか​​疑問に思った。 かなり対称に見える:

table( replicate( 1000000, {x <- round(runif(10), 5); 
                            sig <- sample.int(10); 
               for (i in seq_along(sig)){ x[i] <- x[sig[i]]}; 
                            sum(duplicated(x)) } ) )

     0      1      2      3      4      5      6      7      8 
     1    269  13113 126104 360416 360827 125707  13269    294 

興味深いことに、Rのベクトル化はSIMDとは異なりますが(OPがうまく説明されているように)、ループが「ベクトル化可能」であるかどうかを判断する際にも同じロジックが適用されます。 ここでは、OPの自己回答の例を使ってデモを行います(少し修正します)。

"write-after-read"依存関係を持つ例1は、 "ベクトル化可能"です。

// "ex1.c"
#include <stdlib.h>
void ex1 (size_t n, size_t *x) {
  for (size_t i = 1; i < n; i++) x[i - 1] = x[i] + 1;
}

gcc -O2 -c -ftree-vectorize -fopt-info-vec ex1.c
#ex1.c:3:3: note: loop vectorized

例2の "read-after-write"依存関係は "ベクトル化可能"ではありません

// "ex2.c"
#include <stdlib.h>
void ex2 (size_t n, size_t *x) {
  for (size_t i = 1; i < n; i++) x[i] = x[i - 1] + 1;
}

gcc -O2 -c -ftree-vectorize -fopt-info-vec-missed ex2.c
#ex2.c:3:3: note: not vectorized, possible dependence between data-refs
#ex2.c:3:3: note: bad data dependence

C99 restrictキーワードを使用して、3つの配列間のメモリブロックエイリアシングのないコンパイラにヒントを与えます。

// "ex3.c"
#include <stdlib.h>
void ex3 (size_t n, size_t * restrict a, size_t * restrict b, size_t * restrict c) {
  for (size_t i = 0; i < n; i++) a[i] = b[i] + c[i];
}

gcc -O2 -c -ftree-vectorize -fopt-info-vec ex3.c
#ex3.c:3:3: note: loop vectorized




vectorization