この単純なRループを「ベクトル化」すると異なる結果が得られるのはなぜですか？

準備し始める

ウォームアップとして、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

「ベクトル化」は最初の例では成功しますが、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番目の命令で変更され、この変更された値がx[3]2番目の命令で渡されます。したがってx[3]、と同じ値になりx[1]ます。ただし、並列実行でx[3]は、はになりx[2]ます。その結果、このループを「ベクトル化」することはできません。

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

• 例1には、データに「読み取り後の書き込み」依存関係x[i]があります。読み取り後に変更されます。
• 例2には、データに「書き込み後の読み取り」依存関係がありますx[i]。変更後に読み取られます。

「読み取り後の書き込み」データ依存性を持つループは「ベクトル化」できますが、「書き込み後読み取り」データ依存性を持つループは「ベクトル化」できません。

深く

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

はい。1960年代、高性能コンピューティング用にどのような並列処理コンピュータを設計するのか疑問に思ったとき、フリンは設計のアイデアを4つのタイプに分類しました。カテゴリ「SIMD」（単一命令、複数データ）は「ベクトル化」と呼ばれ、「SIMD」機能を備えたコンピュータは「ベクトルプロセッサ」または「アレイプロセッサ」と呼ばれます。

1960年代には、プログラミング言語はあまりありませんでした。人々はCPUレジスターを直接プログラムするためにアセンブリー（そしてコンパイラーが発明されたときはFORTRAN）を書きました。「SIMD」コンピュータは、単一の命令で複数のデータをベクトルレジスタにロードし、それらのデータに対して同時に同じ演算を実行できます。したがって、データ処理は確かに並列です。例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]

ここでa、bそしてcメモリにエイリアスされていない、すなわち、メモリブロックは、ベクトルを格納しa、bかつ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セッションでは、アドレス<0x28f7340>が指すメモリブロックが割り当てられておりx、コードを実行すると異なる値が表示される場合があります。ただし、の出力はtracememループ後に変更されませんx。つまり、のコピーは作成されません。したがって、ループは実際に余分なメモリを使用せずに「インプレース」変更を行っています。

ただし、ループは「インプレース」順列を実行していません。「インプレース」順列は、より複雑な操作です。の要素をxループに沿って交換する必要があるだけでなく、の要素sigも交換する必要があります（最終的にsigはそうなります1:10）。