#1.はじめに
　RでもPytorchが使えるようになりました。しかも、pythonバックエンドで動かすkerasのようなパッケージではなく、C++ライブラリで実装された本格的なパッケージです。しかし、Quitaにもあまり情報がないので、Rstudo AI blog から紹介してみたいと思います。
基本的にDeepLで訳しました。ブログはGPUを用いていますが、CPU版で実行しました。ブログとは結果が若干異なっていますが、ご了承ください。

#2.実行環境

R version 4.1.0 (2021-05-18) 
RStudio Version 1.4.1717

#3.Please allow me to introduce myself

　昨年1月のrstudio::confでは、カンファレンスがまだ物理的な場所で行われていた遠い過去のことですが、私の同僚であるDanielは、tensorflowエコシステムの新機能と進行中の開発を紹介する講演を行いました。質疑応答では、意外なことを聞かれました。彼はためらった。彼は躊躇しました。実際にそれは計画されていましたし、彼は以前にtorchテンソルをネイティブに実装して遊んだことがありましたが、「それ」がどれほどうまくいくかは完全にはわからなかったのです。

　「つまり、PyTorchホイールをインストールしてreticulate経由でインポートするのではなく、Python Torchにバインドしない実装です。代わりに、テンソル計算と自動微分については基盤となるC++ライブラリlibtorchに委ね、ニューラルネットワーク機能（レイヤー、アクティベーション、オプティマイザ）はRで直接実装します。1つは、ソフトウェアスタックがスリムになることで、インストール時に起こりうる問題が少なくなり、トラブルシューティングの際に探す場所が少なくなります。第二に、Pythonに依存しないことで、torchはユーザーが適切なPython環境をインストールして維持する必要がありません。これは、OSや環境によって大きな違いがあります。例えば、多くの組織では、従業員が自分のラップトップにインストールされた特権的なソフトウェアを操作することは許されていません。

　では、なぜダニエルは躊躇し、私の記憶が正しければ、あまり結論の出ない答えを出したのでしょうか。一方で、libtorchに対するコンパイルが、ある種のオペレーティングシステム上で、深刻な問題を引き起こすかどうかは明らかではありませんでした。(1 一方で、PyTorchをRで再実装するための膨大な作業は、すべてではありませんが、かなりの量になり、気が遠くなりそうでした。今日、やるべきことはまだたくさんありますが(このスレッドは最後に取り上げます)、主な障害は克服され、torchがRコミュニティの役に立つことができる十分なコンポーネントが利用可能になりました。それでは早速、ニューラルネットワークを学習してみましょう。

※torch,torchvisionもCRANから簡単にインストールできます。

library(torch)

インストールと、GPUサポートが問題なく動作しているかどうか（CUDA対応のNVidia GPUがあると仮定して）を素早くチェックするために、CUDAデバイス上にテンソルを作成します。
※ここではCPU版です


torch_tensor(1)

output

torch_tensor
1
[CPUFloatType{1}]

torchには、テンソル、ネットワークモジュール、一般的なデータ読み込み機能がありますが、データタイプに特化した機能は、専用のパッケージで提供されています（または提供される予定です）。一般的に、これらの機能は、データセット、データの前処理と読み込みのためのツール、学習済みモデルの3つのタイプで構成されています。

この記事を書いている時点では、PyTorchには、ビジョン、テキスト、オーディオの3つの分野の専用ライブラリがあります。Rの場合も同様に進めていく予定ですが、「予定」というのは、torchtextとtorchaudioがまだ作られていないからです。現時点では、torchvisionがあれば十分です。

library(torchvision)

#3.1 データの読み込みと前処理

PyTorchにバンドルされているビジョンデータセットのリストは長く、torchvisionに継続的に追加されています。

私たちが今必要としているものはすでに利用可能で、それは - MNIST? ...とまではいきませんが。私のお気に入りの「MNIST dropin」、Kuzushiji-MNIST (Clanuwat et al. 2018)です。MNISTを置き換えるために明示的に作成された他のデータセットと同様に、このデータセットは10個のクラスを持っています。この場合、文字は解像度28x28のグレースケール画像として描かれています。

ここでは、最初の32文字を紹介します。

データをインポートします。

train_ds <- kmnist_dataset(
  ".",
  download = TRUE,
  train = TRUE,
  transform = transform_to_tensor
)

test_ds <- kmnist_dataset(
  ".",
  download = TRUE,
  train = FALSE,
  transform = transform_to_tensor
)

　transformの引数に注目してください。 transform_to_tensorは、画像を受け取り、2つの変換を行います。まず，ピクセルを0から1の範囲に正規化します．次に，もう1つの次元を前面に追加しています。

　今までkerasを使っていた人の予想に反して、この追加次元はバッチ次元ではありません。バッチは次に紹介するdataloaderで処理されます。その代わりに、torchではデフォルトで幅と高さの前にチャンネル数が来ます。

　torchを使っていて非常に便利だと感じたのは、オブジェクトの検査が簡単にできることです。Rの配列やtorchのテンソルではなく、データセットやカスタムオブジェクトを扱っているにも関わらず、簡単に中身を覗くことができます。torchのインデックスは1ベースで、Rユーザーの直感に従ったものになっています。
次のコードを実行します。その結果は、

train_ds[1]

output

$x
torch_tensor
(1,.,.) = 
 Columns 1 to 9  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.2000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0431  0.8314
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.5804  1.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.8039  1.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0275  0.9608  1.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0902  1.0000  0.9373
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.1882  1.0000  0.7529
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.2078  1.0000  0.5804
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.3725  1.0000  0.5608
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.5843  1.0000  0.4745
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.7922  1.0000  0.8353
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0627  0.9647  1.0000  1.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.3490  1.0000  1.0000  1.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.7529  1.0000  1.0000  0.9961
  0.0000  0.0000  0.0000  0.0000  0.1765  0.9882  1.0000  1.0000  0.9882
  0.0000  0.0000  0.0000  0.0000  0.6314  1.0000  1.0000  1.0000  0.9843
  0.0000  0.0000  0.0000  0.2196  0.9882  1.0000  1.0000  1.0000  0.9765
  0.0000  0.0000  0.0000  0.7059  1.0000  1.0000  1.0000  1.0000  0.9804
  0.0000  0.0000  0.0588  0.9882  1.0000  1.0000  1.0000  1.0000  0.9294
  0.0000  0.0000  0.0667  0.9961  0.9529  0.8784  1.0000  1.0000  0.9373
  0.0000  0.0000  0.0196  0.5882  0.2745  0.4039  1.0000  1.0000  0.6588
  0.0000  0.0000  0.0000  0.0000  0.0000  0.2549  1.0000  0.9961  0.2431
  0.0000  0.0000  0.0000  0.0000  0.0000  0.1176  0.9922  0.7294  0.0078
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0314  0.5647  0.1020  0.0000

... [the output was truncated (use n=-1 to disable)]
[ CPUFloatType{1,28,28} ]

$y
[1] 9

※$xが入力　$yがターゲット

データセットの最初の要素である，inputとtargetにそれぞれ対応する2つのテンソルのRリストを与えます
それでは、テンソルの次元を見てみましょう。

train_ds[1][[1]]$size()

output

[1]  1 28 28

※チャンネルはGrayの1チャンネル、28x28ピクセル

データが手に入ったら、それをきれいにバッチ処理して深層学習モデルに供給する人が必要です。
torchでは、これがデータローダの仕事です。

#3.2 データローダ

トレーニングセットとテストセットには、それぞれデータローダがあります。

train_dl <- dataloader(train_ds, batch_size = 32, shuffle = TRUE)
test_dl <- dataloader(test_ds, batch_size = 32)

※バッチサイズ32をデータローダにセット

ここでもtorchは、正しいことをしたかどうかを簡単に確認することができます。最初のバッチの内容を見るには、次のようにします。

train_iter <- train_dl$.iter()
train_iter$.next()

output

$x
torch_tensor
(1,1,.,.) = 
 Columns 1 to 9  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.1098
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.1882  0.6824
  0.0000  0.0000  0.0000  0.0000  0.0000  0.1490  0.5725  0.8510  1.0000
  0.0000  0.0000  0.0000  0.0000  0.1294  0.8588  1.0000  1.0000  1.0000
  0.0000  0.0000  0.0000  0.0000  0.1529  0.9176  1.0000  1.0000  0.9765
  0.0000  0.0000  0.0000  0.0000  0.2549  0.9569  0.9725  1.0000  0.9216
  0.0000  0.0000  0.0000  0.0000  0.0471  0.5490  0.9647  1.0000  1.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0275  0.3490  0.9451  1.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0196  0.4588  0.9725  1.0000
  0.0000  0.0000  0.0000  0.0000  0.0902  0.6196  0.9765  0.9922  0.8353
  0.0000  0.0000  0.0000  0.1569  0.7765  1.0000  1.0000  0.9176  0.2392
  0.0000  0.0000  0.1294  0.8392  1.0000  1.0000  1.0000  0.9059  0.2118
  0.0000  0.0510  0.7373  1.0000  1.0000  1.0000  1.0000  0.8902  0.1882
  0.0000  0.1686  0.9412  1.0000  1.0000  1.0000  1.0000  0.8980  0.4196
  0.0000  0.0902  0.8784  1.0000  1.0000  0.9608  0.8980  1.0000  0.9608
  0.0000  0.0118  0.6039  1.0000  0.9961  0.6431  0.2863  0.9333  1.0000
  0.0000  0.0000  0.0863  0.7804  0.6980  0.1255  0.1373  0.8784  1.0000
  0.0000  0.0000  0.0000  0.0745  0.0471  0.0784  0.7569  0.9922  1.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.2824  0.9882  1.0000  1.0000
  0.0000  0.0000  0.0000  0.0000  0.0078  0.6431  1.0000  1.0000  1.0000
  0.0000  0.0000  0.0000  0.0000  0.2118  0.9255  1.0000  1.0000  1.0000
  0.0000  0.0000  0.0000  0.0275  0.6431  1.0000  1.0000  1.0000  1.0000
  0.0000  0.0000  0.0000  0.2980  0.9725  1.0000  1.0000  1.0000  1.0000
  0.0000  0.0000  0.0000  0.4039  1.0000  1.0000  1.0000  1.0000  1.0000
  0.0000  0.0000  0.0000  0.2275  0.9451  1.0000  1.0000  1.0000  0.9608
  0.0000  0.0000  0.0000  0.2118  0.8902  1.0000  1.0000  0.9647  0.4471
  0.0000  0.0000  0.1098  0.8078  1.0000  1.0000  1.0000  0.7647  0.0510
  0.0000  0.0000  0.2863  0.9961  1.0000  1.0000  0.9961  0.5059  0.0039

... [the output was truncated (use n=-1 to disable)]
[ CPUFloatType{32,1,28,28} ]

$y
torch_tensor
  2
  1
  6
  7
  9
  2
  1
  4
  5
  1
  8
  8
  3
  5
  8
  3
  8
  4
  5
  5
 10
  4
  8
  4
  8
  6
  9
  1
  4
  2
... [the output was truncated (use n=-1 to disable)]
[ CPULongType{32} ]

このような機能は、よく知られたデータセットを扱うときには必要ないと思われるかもしれませんが、ドメイン固有の前処理がたくさん必要なときには非常に役に立つことがわかります。
データを読み込む方法がわかったところで、データを視覚化するための前提条件がすべて整いました。ここでは、上の最初の文字を表示するためのコードを紹介します。

par(mfrow = c(4,8), mar = rep(0, 4))
images <- train_dl$.iter()$.next()[[1]][1:32, 1, , ] 
images %>%
  purrr::array_tree(1) %>%
  purrr::map(as.raster) %>%
  purrr::iwalk(~{plot(.x)})

これでネットワークを定義する準備ができました。

#3.3 ネットワーク

kerasのカスタムモデルを使ったことがある人（またはPyTorchを使ったことがある人）は、以下のようにネットワークを定義する方法にあまり違和感を覚えないかもしれません。

net <- nn_module(
  
  "KMNIST-CNN",
  
  initialize = function() {
    # in_channels, out_channels, kernel_size, stride = 1, padding = 0
    self$conv1 <- nn_conv2d(1, 32, 3)
    self$conv2 <- nn_conv2d(32, 64, 3)
    self$dropout1 <- nn_dropout2d(0.25)
    self$dropout2 <- nn_dropout2d(0.5)
    self$fc1 <- nn_linear(9216, 128)
    self$fc2 <- nn_linear(128, 10)
  },
  
  forward = function(x) {
    x %>% 
      self$conv1() %>%
      nnf_relu() %>%
      self$conv2() %>%
      nnf_relu() %>%
      nnf_max_pool2d(2) %>%
      self$dropout1() %>%
      torch_flatten(start_dim = 2) %>%
      self$fc1() %>%
      nnf_relu() %>%
      self$dropout2() %>%
      self$fc2()
  　　　　　　　　　#最後の層にrelu活性化関数がありませんが、後ほど説明があります。　　　　　　　　　　　
　}
)

nn_module()を使って、ネットワークのコンポーネントを格納するR6クラスを定義します。そのレイヤーは initialize() で作成され、forward() はネットワークのフォワードパスで何が起こるかを記述します。用語の説明をします。torchでは、レイヤーはネットワークと同様にモジュールと呼ばれています。これには意味があります。torchでは、レイヤーをネットワークと同様にモジュールと呼びます。

レイヤ（申し訳ありませんが、モジュール）自体は見覚えがあるかもしれません。当然のことながら、nn_conv2d()は2次元の畳み込みを行い、nn_linear()は重み行列の乗算とバイアスのベクトルを加えます。しかし、nn_linear(128, 10)のような数字は何でしょうか？

torchでは、レイヤーのユニット数ではなく、レイヤーを通過する「データ」の入出力次元を指定します。つまり、nn_linear(128, 10)は、128個の入力接続を持ち、10個の値を出力します（各クラスに1個ずつ）。このケースのように、次元を指定するのが簡単な場合もあります。入力辺の数（つまり、前の層の出力辺の数と同じ）はわかっていますし、必要な出力値の数もわかっています。しかし、前のモジュールではどうでしょうか。どうやって9216個の入力コネクションにたどり着くのでしょうか？

ここでは、ちょっとした計算が必要です。形状に影響を与えるものであれば、その変換を記録し、そうでなければ無視します。

そこで、まず、batch_size x 1 x 28 x 28 という形状の入力テンソルを用意します。すると,

1)nn_conv2d(1, 32, 3)、または nn_conv2d(in_channels = 1, out_channels = 32, kernel_size = 3)は、カーネルサイズ3、ストライド1（デフォルト）、パディングなし（デフォルト）の畳み込みを行います。ドキュメントを参照して結果の出力サイズを調べることもできますが，直感的に，カーネルサイズ 3 でパディングなしの場合，画像は各方向に 1 ピクセルずつ縮小され，26 x 26 の空間解像度になると考えることもできます．32のチャンネルごとにです。したがって、実際の出力形状は、 batch_size x 32 x 26 x 26 となります。
※最初の画像サイズ28 -2=26サイズ

2)nnf_relu()はReLUの起動を行いますが、形状には一切触れません。

3)nn_conv2d(32, 64, 3)では、パディングゼロ、カーネルサイズ3の畳み込みを行います。このときの出力サイズは batch_size x 64 x 24 x 24 です。2回目のnnf_relu() は出力形状には何もしません
※画像サイズ26-2＝24

4)nnf_max_pool2d(2) (等価: nnf_max_pool2d(kernel_size = 2))では、拡張子2×2の領域に最大プーリングを適用して、出力をbatch_size x 64 x 12 x 12のフォーマットに縮小します。
※画像サイズ24/2=12

5)nn_dropout2d(0.25)は形状的には問題ありませんが、後でリニアレイヤーを適用したい場合は、すべてのチャンネル、高さ、幅の軸を1つの次元に統合する必要があります。

6)torch_flatten(start_dim = 2). 出力形状は batch_size * 9216 となり、64 * 12 * 12 = 9216 となります。このようにして、9216個の入力接続がであがります。

7)nn_linear(9216, 128) に入力されます。9216から128にアウトプットされます。

8)nnf_relu() と nn_dropout2d(0.5) は寸法は変えません

9)nn_linear(128, 10)は、10個のクラスのそれぞれに1つずつ、望ましい出力スコアを与えます。

　もっと複雑なネットワークだったら......と思うでしょう。計算が非常に面倒になるかもしれません。幸い、torchの柔軟性を利用すれば、別の方法があります。各レイヤーは独立して呼び出し可能なので、サンプルデータを作って何が起こるか見てみましょう。

ここに「画像」のサンプルがあります。正確には、「画像」を含む1つのアイテムのバッチです。

x <- torch_randn(c(1, 1, 28, 28))

※torch_randn :　平均0，分散1の正規乱数のテンソルを作成する

その上で、最初のconv2dモジュールを呼び出したらどうでしょうか。

conv1 <- nn_conv2d(1, 32, 3)
conv1(x)$size()

output

[1]  1 32 26 26

といった具合です。これは、torchの柔軟性がいかにニューラルネットの開発を容易にするかを示す一例に過ぎません。
メインスレッドに戻ります。モデルをインスタンス化します。

model <- net()

入力データと出力データについて、これはトレーニングループの中で行われます。

#3.4 トレーニング

torchでは、オプティマイザを作成する際に、モデルのパラメータという操作対象を指定します。

損失関数はどうするか？2つ以上のクラスを持つ分類には、クロスエントロピーを使います。

kerasのcategorical cross entropyとは異なり、torchのnnf_cross_entropy()は、ソフトマックス活性化を適用して得られるような確率を含む予測を想定していますが、生の出力(logits)で動作します。これが、ネットワークの最後の線形層に活性化が行われなかった理由です。

実際、学習ループは二重構造になっています。エポックとバッチをループしています。バッチごとに、入力に対してモデルを呼び出し、損失を計算して、オプティマイザーに重みを更新させます。

optimizer <- optim_adam(model$parameters)

for (epoch in 1:5) {

  l <- c()

  coro::loop(for (b in train_dl) {
    # make sure each batch's gradient updates are calculated from a fresh start
    optimizer$zero_grad()
    # get model predictions
    output <- model(b[[1]])
    # calculate loss
    loss <- nnf_cross_entropy(output, b[[2]])
    # calculate gradient
    loss$backward()
    # apply weight updates
    optimizer$step()
    # track losses
    l <- c(l, loss$item())
  })

  cat(sprintf("Loss at epoch %d: %3f\n", epoch, mean(l)))
}

output

Loss at epoch 1: 0.385767
Loss at epoch 2: 0.187037
Loss at epoch 3: 0.142506
Loss at epoch 4: 0.114186
Loss at epoch 5: 0.097570

メトリクスの計算や検証セットでの性能評価など、できることはたくさんありますが、上記は典型的な（単純ではありますが）torchトレーニングループのテンプレートです。
特にオプティマイザー関連のイディオムは何度も目にすることになるでしょう。
最後に、テストセットでのモデルのパフォーマンスを評価してみましょう。

#3.5 評価

モデルをevalモードにすると、torchはこの後の処理でグラデーションの計算やバックプロップを行わないようにします。
テストセットを反復して、損失とバッチで得られた精度を記録します。

test_losses <- c()
total <- 0
correct <- 0

coro::loop(for (b in test_dl) {
  output <- model(b[[1]])
  labels <- b[[2]]
  loss <- nnf_cross_entropy(output, labels)
  test_losses <- c(test_losses, loss$item())
  # torch_max returns a list, with position 1 containing the values 
  # and position 2 containing the respective indices
  predicted <- torch_max(output$data(), dim = 2)[[2]]
  total <- total + labels$size(1)
  # add number of correct classifications in this batch to the aggregate
  correct <- correct + (predicted == labels)$sum()$item()
})

mean(test_losses)

output

[1] 0.3709392

ここでは，正しい分類の割合として計算される平均精度を示しています。

test_accuracy <-  correct/total
test_accuracy

output

[1] 0.9053

第1回目のtorchの例はここまでです。

#3.6 学ぶ

より詳しく知りたい方は、torchサイトのヴィネットをご覧ください。まずは、以下の項目をご覧ください。

"Getting started "シリーズ。シンプルなニューラルネットワークをゼロから構築します。低レベルのテンソル操作から始めて、自動微分やネットワークモジュールなどの高レベルの機能を徐々に追加していきます。
テンソルの詳細。テンソルの作成と索引付け
torchでのバックプロパゲーション：autograd
ご質問や問題がありましたら、GitHubやRStudioコミュニティフォーラムでお気軽にお尋ねください。

#4 参考
GPU版のコードはブログを参考にしてください。
please allow me to introduce myself: Torch for R
#5 Enjoy!