完全情報最尤推定法（FIML）で操作変数のある連立方程式モデルを推定

計量経済学では同時方程式を考えることが多く、そのために一般最小二乗法（OLS）の発展である操作変数法（IVM）とそのバリエーションである二段階最小二乗法（2SLS）*1や三段階最小二乗法（3SLS）、モーメント法の発展である一般化モーメント法（GMM）がよく使われているのですが、これらは尤度函数を用いないのでベイズ推定にできません。頻度主義者がベイジアンになるのは難しいです。
しかし、完全情報最尤推定法（FIML）と言うのがあり、実際のところは簡単に最尤法で同時方程式を推定できます。FIMLはRでは欠損値処理アルゴリズムのように扱われている気がするのですが、教科書的な手法なので確認しておきましょう*2。最尤法だけに不均一分散や欠損値処理の対処も工夫できますし、ベイズ推定しない場合でも、便利に使える局面は多そうです。

理論モデルとOLS推定量に入る同時性バイアス

推定の例をつくるわけですが、需要曲線と供給曲線を考えましょう*3。

需要関数：
供給関数：
均衡条件：

は需要、は供給、は価格、は天候などの操作変数、は係数、は誤差項です。ここで均衡価格を計算すると、

均衡価格：

となります。観測データは均衡価格になるわけですが、例えば需要関数の推定でを説明変数に用いると、の第3項にがあるため、となりOLSの前提を満たしません。需要関数の価格の係数の推定量はとバイアスが入ります。

完全情報最尤推定法

需要関数と均衡価格式を同時推定することを考えましょう。推定する係数を、誤差項をと置きなおし、方程式一本ごとに誤差項があると考え、二本の方程式を行列にまとめて書くと、次のようになります。は観測数です。

前節の係数とは、やと言った対応関係があります。
右辺の内生変数と外生変数を分離します。

内生変数を左辺にまとめます。

左辺第1項と第2項をまとめましょう。

この行列をいちいち書いていられないので、書き直します。

FIMLでは、この構造方程式を推定します。こう書くとOLSで入るバイアスが入らない理由が分かりづらいですが、この構造方程式の尤度函数は説明変数に内生変数が入らない誘導形の尤度函数から導かれます*4。また、数理的な特性を使って整理すると以下のようなconcentrated loglikelihood function*5にできることが知られています*6。

は内生変数（もしくは左辺にまとめられた被説明変）の数で、ここでは2です。の4つのパラメーターを上の対数尤度函数で推定すればよいです。

データ

理論モデルに沿った式で、乱数を使ってデータを生成します。

set.seed(1231)
n <- 300
a10 <- 100
a11 <- -1
a20 <- 2
a21 <- 3
a22 <- 4
z <- runif(n, min=0, max=3)
mu <- rnorm(n, mean=0, sd=2)
nu <- rnorm(n, mean=0, sd=1)
p <- (a20 - a10 + a22*z + nu - mu)/(a11 - a21)
d <- a10 + a11*p + mu
s <- a20 + a21*p + a22*z + nu

推定

concentrated loglikelihood functionをoptim関数で推定します。

Y <- cbind(d, p)
W <- cbind(rep(1, n), z)
g <- 2

llf <- function(p){
    B <- matrix(c(p[1], 0, p[2], p[3]), 2, 2)
    G <- matrix(c(1, p[4], 0, 1), 2, 2) # Γ
    E <- Y%*%G - W%*%B
    S <- t(E)%*%E/n # ∑
    -g*n/2*(log(2*pi) + 1) + n*log(det(G)) - n/2*log(det(S))
}

r_optim <- optim(c(1, 0, 0, 0), function(p){ -llf(p) }, method="BFGS")

r_optim$parが推定されたになりますが、内生変数の係数となる4番目のパラメーターの符号が逆になることに注意してください。なお、optim関数の引数にhessian=TRUEをつけるとヘッセ行列がとれるので、標準誤差を計算することができます*7。

OLSとIVとの比較

推定結果の比較をしておきましょう。FIMLはOLSよりはかなりマシ、IVよりは僅かに劣るといった感じでしょうか。

係数	真の値	OLS	IV	FIML
	100	75.95501	100.11595	100.1228135
	-1	0.03949	-1.01101	-1.0113085

なお、最尤法は漸近的に有効だけに、サンプルサイズが小さいとIVとの差が広がる傾向があります。IV推定量は以下のように計算できるので、nを変えて試してみてください。

zm <- matrix(c(rep.int(1, n), c(z)), n, 2)
xm <- matrix(c(rep.int(1, n), c(p)), n, 2)
solve(t(zm) %*% xm) %*% (t(zm) %*% d)

*1:操作変数が1つの場合はIVMと同値なのですが、一応、別なので。

*2:昔読んだGreene (2003)やWooldridge (2002)では最尤法で連立方程式モデルの推定の仕方は書いていなかったのですが、Davidson and MacKinnon (2003)では第12章で詳しく説明されています。

*3:Hayashi (2000)にある例を踏襲しています。

*4:Davidson and MacKinnon (2003)もしくはオンラインで配布している2021年版のpp.503–505とpp.544–545を参照。誘導形はSURと同じ多変量正規分布の尤度函数になるので、による変形の影響を整理すれば構造形の尤度函数が導出できます。

*5:定訳不明

*6:一階条件の節の一階条件に直接関係のないExercise 12.24をに注意しての成分を考えて解くことで得ることができる･･･ハズ。

*7:Rと手作業で覚える最尤法 - 餡子付゛録゛

MCMCpackのMCMCmetrop1RでBayes Factorを雑に計算

マルコフ連鎖モンテカルロ法（MCMC）を使いたいベイジアンの統計ユーザーはStanを使っている人が多く、R使いもフロントエンドRStanを通して利用していることがほとんどだと思うのですが、MCMCpackと言うのもあります。このMCMCpackでベイズ因子（Bayes Factor）を使ってモデル選択をしてみましょう。

MCMCpackにはよく使われるらしい特定の共役事前分布*1を置ける線形回帰やポアソン回帰などの関数が用意されており、それらはベイズ因子を簡単に出力できるように設計されています。しかし、今回は正規分布と対数正規分布のどちらが当てはまりがよいかモデル選択をするので使えません。そして、汎用的なベイズ因子の計算関数は無いようです。

詰まった感があるのですが、推定モデルごとに周辺尤度を計算できればベイズ因子は簡単に計算できます。そして、事後分布の推定が終わっていれば、ベイズの定理から周辺尤度を簡単に計算できます。精度を考えなければ。
ベイズの定理の事後確率と周辺尤度を入れ替えてみましょう。

はデータと言うか観測値のベクトル、はパラメーターのベクトル、は確率、は条件付き確率です。は任意のですが、精度の都合でMAP推定量など、なるべく大きな値を選びます。はで計算した尤度、はの事前確率です。

1. データセット

今回は以下の対数正規分布からの乱数を使います。

set.seed(20211204)
n <- 15
y <- rlnorm(n, meanlog=1.2, sdlog=0.8)

真のパラメーターが分かっている方が、推定がうまくいっているか分かりやすいので。

ヒストグラムを描くと、こんな感じです。

2. 対数尤度関数

正規分布と対数正規分布の対数尤度関数は以下になります。

# 正規分布の対数尤度関数
llf_nd <- function(theta){
  if(any(theta < 0.0)){
    return(-Inf)
  }
  sum(dnorm(y, mean=theta[1], sd=theta[2], log=TRUE))
}

# 対数正規分布の対数尤度関数
llf_lnd <- function(theta){
  if(any(theta < 0.0)){
    return(-Inf)
  }
  sum(dlnorm(y, meanlog=theta[1], sdlog=theta[2], log=TRUE))
}

3. MCMCで推定

MCMCpackをインストールしてあったら、簡単にできます。なお、MCMCmetrop1R関数は初期値をoptimでMAP推定してからMCMCを開始するため、最適化アルゴリズムがMAP推定に失敗するとエラーを出します。そしてデフォルト指定のBFGSがエラーを出してきたので、L-BFGS-Bに代えました。

init_theta <- c(mean(y), sd(y))
library(MCMCpack)

# 正規分布のモデルの事前分布は(0 7)と(0 5)の一様分布
prior.nd <- function(theta){
    dunif(theta[1], min=0, max=7, log=TRUE) + dunif(theta[2], min=0, max=5, log=TRUE)
}
post.samp.nd <- MCMCmetrop1R(function(theta){
  llf_nd(theta) + prior.nd(theta)
}, theta.init=log(init_theta), optim.method="L-BFGS-B")

# 対数正規分布のモデルの事前分布は(0 log(7))と(0 log(5))の一様分布
prior.lnd <- function(theta){
    dunif(theta[1], min=0, max=log(7), log=TRUE) + dunif(theta[2], min=0, max=log(5), log=TRUE)
}
post.samp.lnd <- MCMCmetrop1R(function(theta){
  llf_lnd(theta) + prior.lnd(theta)
}, theta.init=log(init_theta))

4. 中央値から周辺尤度を計算

今回は手軽に中央値を使います。ノンパラメトリックな確率密度関数の推定法…とかっこよく言いたいところですが、MCMCを繰り返して精度をあげるChib (1995)のような事はしていないので雑です。

# 中央値を得る
medians.nd <- summary(post.samp.nd)$quantiles[, "50%"]
medians.lnd <- summary(post.samp.lnd)$quantiles[, "50%"]

# パラメーターyに対応する事後確率を求める
get_posterior_probability <- function(post.samp, y, nos=100){
    p <- length(y) # 次元
    nr <- length(post.samp[, 1])
    count <- function(a){
        v <- 1
        flag <- rep(TRUE, nr)
        for(column in 1:p){
            s <- post.samp[, column]
            range <- max(s) - min(s)
            # y±range*aに含まれる点を選ぶ
            flag <- flag & s <= y[column] + range*a & s >= y[column] - range*a
            v <- v * (2*range*a)
        }
        c("nflag"=sum(flag), "volume"=v)
    }
    # 区間に入る点がnos個になるように、区間を調整する
    l <- 0
    h <- 1
    m <- 0.5
    nz <- h
    for(i in 1:64){
        m <- (l + h)/2
        n <- count(m)[["nflag"]] - nos
        if(n == 0){
            break;
        } else if(n < 0){
            l <- m
        } else {
            h <- m
            nz <- m
        }
    }
    r <- count(ifelse(0<=n, m, nz))
    # r["nflag"]/length(post.samp)は確率（累積分布）になるので、区間の大きさr["volume"]で割って（平均）確率密度にする
    r[["nflag"]]/length(post.samp)/r[["volume"]]
}

# 中央値の事後確率
posterior_p.nd <- get_posterior_probability(post.samp.nd, medians.nd)
posterior_p.lnd <- get_posterior_probability(post.samp.lnd, medians.lnd)

# 対数周辺尤度を作る
mll.nd <- llf_nd(medians.nd)-log(posterior_p.nd) + prior.nd(medians.nd)
mll.lnd <- llf_lnd(medians.lnd)-log(posterior_p.lnd) + prior.lnd(medians.lnd)

5. ベイズ因子を計算

割り算と言うか、引き算して指数をとって終わりです。

sprintf("Bayes Factor: %.3f", exp(mll.lnd-mll.nd))

8.967でSubstantialと言えるような値になりますが、ヒストグラムの見た目にあった結果ではないでしょうか。なお、シャピロ・ウィルク検定をかけると10%有意で正規分布が棄却されて、対数正規分布は棄却されません。

shapiro.test(x=y)$p.value
shapiro.test(x=log(y))$p.value

6. 注意事項

ほぼ思いつきで書いた胡散臭い方法なので、実用の参考にはしないでください*2。

*1: MCMCregress関数だと、b0, B0, c0, d0, sigma.mu, sigma.varで事前分布になる多変量正規分布、逆ガンマ分布のパラメーターを指定できます。

*2: Chib (1995)の方法で計算しなおしたら7.154になりました。8.967は目安として使えなくもないですが、論文に書くと怒られると思います。

Rの主成分回帰で多重共線性を制御

一般線形回帰（OLS）でコスト分析をしたら係数がマイナスの項目が出てきた、工数をかけるほど費用が減っていく･･･と言う嘆きを見かけたのですが、説明変数の相関が高いと推定される係数の標準誤差が大きくなって推定量がおかしくなる、多重共線性による問題だと考えられます*1。

こういう問題を解決する素朴な方法として、主成分回帰（Principal Component Regression）で多重共線性を制御する技があるので紹介します。試してみると望んだ結果が得られるかも知れないので、学部生で卒論に使う回帰分析の係数が気に入らない人にp-hacking的にお勧めです。なお、べたべたとコードを書いていきますが、plsパッケージで手軽に実行でき、さらに同様の目的に使えるPLS回帰も試せるので、基本的にはplsパッケージを使う方が賢いです。

さて、まず解決すべき多重共線性をつくってみましょう。最後に比較するときに、真のモデルの切片項と係数はそれぞれ1であることに注意してください。

set.seed(20201115)
n <- 150 # 十分なサンプルサイズ（e.g. 1500）にすると多重共線性が問題なくなる
x1 <- runif(n)
x2 <- x1 + rnorm(n, sd=0.1)
y <- 1 + x1 + x2 + rnorm(n)

VIFは8.02675です。
計算結果を格納する行列をつくります。

beta <- matrix(0, 3, 2)
se <- matrix(0, 3, 2)
colnames(beta) <- colnames(se) <- c("OLS", "PCR")
rownames(beta) <- rownames(se) <- c("Const.", "x1", "x2")

比較のためにOLSで回帰します。

lm_r <- lm(y ~ x1 + x2)
beta[, "OLS"] <- coef(summary(lm_r))[, 1]
se[, "OLS"] <- coef(summary(lm_r))[, 2]

主成分分析を行います。scale=TRUEを忘れないようにしましょう。

r_prcomp <- prcomp(~ x1 + x2, scale=TRUE)

寄与率を低い主成分を省略し、回帰を行います。今回は、summary(r_prcomp)してCumulative Proportionが0.9を超えているので、第1主成分だけで行います*2。回帰に用いる主成分の数と説明変数の数が一致すると、OLSと同じ予測力になります。

npc <- 1 # 回帰に用いる主成分の数
lm_prc_r <- lm(y ~ predict(r_prcomp)[,1:npc])

第1主成分を、例えば平均全体工数のような抽象的な概念に対応するとして分析を進めてもよいのですが、そのような概念をあまり説明したくない場合は、主成分の係数を、元の説明変数の係数に分解しましょう。第1主成分の回帰係数と第1主成分の固有ベクトル（説明変数を標準化した値の係数）から、元の説明変数の回帰係数を求めます*3。

beta["x1", "PCR"] <- sum(r_prcomp$rotation[1,][1:npc]*coef(lm_prc_r)[2:(npc + 1)]/sd(x1))
beta["x2", "PCR"] <- sum(r_prcomp$rotation[2,][1:npc]*coef(lm_prc_r)[2:(npc + 1)]/sd(x2))
beta["Const.", "PCR"] <- coef(lm_prc_r)[1] - beta["x1", "PCR"]*mean(x1) - beta["x2", "PCR"]*mean(x2)

続けて分散の加法公式を使って、上の計算手順をもとに、元の説明変数の擬似的な標準誤差を求めます。

se["x1", "PCR"] <- sqrt(sum((r_prcomp$rotation[1, ][1:npc]/sd(x1) * coef(summary(lm_prc_r))[, 2][2:(npc + 1)])^2))
se["x2", "PCR"] <- sqrt(sum((r_prcomp$rotation[2, ][1:npc]/sd(x2) * coef(summary(lm_prc_r))[, 2][2:(npc + 1)])^2))
se["Const.", "PCR"] <- sqrt(
    coef(summary(lm_prc_r))[, 2][1]^2
    + sum(((
        r_prcomp$rotation[1,][1:npc]*mean(x1)/sd(x1)
         + r_prcomp$rotation[2,][1:npc]*mean(x2)/sd(x2)
        ) * coef(summary(lm_prc_r))[, 2][2:(npc + 1)])^2)
)

t検定に使う場合の自由度は、標本サイズから、回帰に使った主成分の数と、切片項の分の1を引いたものになるのに注意してください。
OLSとPCRの推定量と残差平方和を比較しておきましょう。

# 推定量を比較する
cat("\nestimated coeffcients:\n")
print(beta)

cat("\nstandard errors:\n")
print(se)

# 残差平方和を比較する
cat(
    sprintf("\nOLSの残差平方和: %.3f\n主成分回帰の残差平方和: %.3f\n", 
        sum((y - predict(lm_r))^2), 
        sum((y - model.matrix(~ x1 + x2) %*% beta[, 2])^2)
    )
)

係数	真のモデル	OLS	PCR
Const.	1	1.10562969	1.0806319
x1	1	-0.01407916	0.9303908
x2	1	1.80575174	0.8960958

標準誤差	OLS	PCR
Const.	0.1775728	0.09022349
x1	0.8366178	0.05424751
x2	0.8057794	0.05224790

OLSの残差平方和は143.823、PCRの残差平方和は145.11と予測力は僅かに悪化しますが、係数を見るとPCRの推定量の方が真のモデルに近い事がわかります。

*1:さらにomitted variable biasで、見えない成分を代理してしまっている可能性もあります。例えば、常連客との取引では低コストになる仕組みがあり、常連客によく提供するサービス項目の工数であったりすると、常連客のコスト低減効果を代理することもありえます。

*2:省略する成分は、偶発的に生じた見かけ上の、もしくは想定外だが従属変数への影響が大きくない構造（e.g. 上述の常連客のコスト低減効果）になります。

*3:を主成分回帰の切片工、を主成分回帰の回帰式の係数、を主成分の説明変数の係数、主成分回帰の誤差項をとすると、主成分回帰は と書けるので、式をの項で分けて整理すると、OLSと同様の説明変数とその係数の式に直せます。

VSCodeで最小限のLaTeX環境を整備する

Visual Studio Code（以下、VSCode）1.5でLaTeX環境を整備したのですが、参考にしたページの設定が心理的な参入障壁になりそうなほど頑張っていたので、備忘録をかねてWindows 10でplatex/pbibtexを使うだらっとした最小限の設定を紹介します。

1. LaTeXのインストールと設定

4GBあるTeX Live 2020のISOイメージをダウンロードしてきてインストールします。Windows 10であれば、ダウンロードしたファイルを右クリックしてマウントすると、ドライブが割り当てられるので、ドライブを開いてinstall-tl-windows.batを実行します。

イマドキはlatexmkコマンドを使ってTeXのコンパイルをするので、その設定になる%USERPROFILE%\.latexmkrcを作ります。notepadに、

$latex = 'platex -synctex=1 %O %S';
$bibtex = 'pbibtex %O %B';
$dvipdf = 'dvipdfmx %O -o %D %S';
$makeindex = 'mendex %O -o %D %S';
$max_repeat = 10;

と書いて、ファイル種類を*.*、ファイル名を%USERPROFILE%\.latexmkrcで、名前をつけて保存してください。これでlatexmk -pdfdvi example.texで、example.pdfが作れるようになります。
分野ごとの参考文献の表示スタイル（e.g. jecon.bst）などもこの段階で追加しておきましょう*1。

2. VSCodeのインストールと設定

VSCodeの公式サイトからダウンロードしてきて、インストールします。続けて、LaTeX WorkshopとJapanese Language Pack for Visual Studio Codeの2つの拡張機能を、ショートカットCTRL + SHIFT + Xを押すか、以下のアイコンをクリックしてVS Code Marketplaceのペーンを呼び出し、検索してインストールします *2。

これで.texファイルと.bibファイルの編集でシンタックスハイライトと言うか色分け表示してくれるようになり、コード補間が効くようになります。

次に、VSCodeで.texファイルを編集中にウィンドウ内右上の三角形のRunボタンを押したときの挙動を設定します。notepadかVSCodeで%APPDATA%\Code\User\settings.jsonを開いて、

    // for LaTeX
    "latex-workshop.latex.recipes": [
        {
            "name": "latexmk",
            "tools": [
                "latexmk"
            ]
        },
    ],
    "latex-workshop.latex.tools": [
        {
            "name": "latexmk",
            "command": "latexmk",
            "args": [
                "-pdfdvi",
                "%DOC%",
            ],
        },
    ],

をファイル内先頭の{と最後の}の間に挟みます。なお、このファイルはフォントサイズやタブ幅などVSCodeのエディターの諸々の設定を書いておける場所で作りこむことが可能です*3が、上のようなコマンド変数以外は、VSCodeでCTRL + ,を押してダイアログ的な画面を使った方がよいかも知れません。

*1: jecon.bstの場合は、jecon.bstを落としてきて、C:\texlive\2020\texmf-dist\pbibtex\bstに放り込むだけになると思います。

*2: Managing Extensions in Visual Studio Code

*3: vscodeがいい感じになる設定を継ぎ足していく（俺流vscode秘伝のタレ） - Qiita

OpenBLASをリンクしたWindows版R 4.1 PR

ループでスカラ演算を繰り返すよりもベクトル演算をした方が速いRですが、標準構成ではNetlib BLASと言う線形代数ライブラリBLASのリファレンス実装を使っていて、行列演算、ベクトル演算も速いとは言えないものとなっています。チューニングされたOpenBLASに差し替えることで、高速化をしてみましょう。

Ubuntu Linuxでは、

sudo apt install libopenblas0-pthread

で、OpenBLASがインストールされ、aptでインストールしたRはOpenBLASを使うようになります。簡単ですね。しかし、Windowsだと時間のかかる作業が必要です。
基本的にはチートペーパーを見ながら、Rの開発版をコンパイルできる環境を揃えて、ビルドに使う設定ファイルを書き換えれば終わりなのですが、微妙に書いてある通りではコンパイルできませんでした。

まず、Step 4: Create blas patchのところですが、現行のMakefile.winがホワイトスペースに空白ではなくタブを使っているせいか、示されているblas.diffではパッチがあたりません。blas.diffをつくりなおす必要がありました。作り直したblas.diffをr-base-master.zipを展開してできるフォルダーC:\r-base-masterの中に入れました。

次に、Step 5: Adjust existing filesのところですが、示されているPKGBUILDの末尾の</jeroen@berkeley.edu>は要らないので消しました。見比べるとmakedependsに"${MINGW_PACKAGE_PREFIX}-openblas"、sourceにblas.diff、sha256sumsに'SKIP'、# Add your patches hereにpatch -Np1 -i "${srcdir}/blas.diff"をそれぞれ追加しているだけのようなので、チートペーパーのPKGBUILDをコピペする必要は無いかもです。（これはチートペーパーに記述の問題はありませんが）MkRules.local.inのQPDFへのパスは、

QPDF = C:/qpdf-10.0.1

と、インストール先のQPDFのホームにしました。
export PATH="$PATH:/c/progra~1/MiKTeX 2.9/miktex/bin/x64"の行は、

export PATH="$PATH:/c/progra~1/MikTeX/miktex/bin/x64"

と、（自明ですが）MikTeXのインストール先にあわせました。

設定は以上で終わりですが、別のチートペーパーを参考にビルドする前にMikTeXでフォントを作成しておかないとFont ts1-zi4r at 540 not foundとエラーが出て止まります*1。

Step 6: Build Rはマイナーなところですが、rtools40のmsys.exeを使えと書いてありますが、手元の環境ではC:\rtools40\msys2.exeでした。
なお、コンパイル途中でQPDFがインストール確認をしてくるので、少なくとも1回は了解ボタンを押すまでは、放置しておいてもコンパイルは終了しません。また、マイナーな計算結果差のようですが、微妙に1つのテストでエラーが出ました。これは4.1 PRのせいかも知れません*2。

コンパイルが終わってR-devel-win.exeができたらすぐにインストールしてみたのですが、Windows版ではsessionInfo()をしても使っているBLASの種類を教えてくれないです。以下のように行列演算をさせ、4コアCPUで35倍以上の経過時間差を出して、差し替えた感を実感しました。

set.seed(1013)
n <- 5000
M1 <- matrix(runif(n^2, min=1, max=10), n, n)
M2 <- matrix(runif(n^2, min=1, max=10), n, n)
system.time({ M3 <- M1 %*% M2 })

ところで先日書いてみたマクロ経済学の数値解析は、ほとんど速くなりませんでした。速くなりませんでした。

*1:チートペーパーのコマンドをうってinitexmf: security risk: running with elevated privilegesとエラーが出る人は、付属ユーティリティMiKTeX Consoleで操作した方がよいかもです。起動してadministration modeを選択し、メニューのTasksのRefresh file name database、Refresh font map filesを順番に実行します。

*2:export rsource_url='https://cran.r-project.org/src/base-prerelease/R-latest.tar.gz'; ./full-build.shとコンパイルを実行すると、開発版ではなくリリース候補版のソースコードを指定できます。

Rの拡張でOpenMPを使ってみる

プロセッサのマルチコア化が進んだ現代なので、時間のかかる計算では並列処理をするコードを書くのが望ましいです。スレッドやセマフォの制御を直接プログラミングすると骨が折れるのですが、シングルコア向け逐次処理コードを並列処理に変換してくれるOpenMPという仕組みがあります。Cで書いたRの拡張で使ってみましょう。

1. Makevarsの編集

gccでOpenMPを使うためには、引数に-fopenmpフラグを加える必要があります。R CMDでこの引数をつけるためには、LinuxやMacOS Xでは$HOME/.R/Makevarsに、MS-Windowsでは%HOME%*1/.R/Makevars.winと言うテキストファイルに、

CFLAGS= -fopenmp
CXXFLAGS= -fopenmp
# PKG_CXXFLAGSか、PKG_CXXFLAGSに-std=c++11と書いてCXX11FLAGSに-fopenmpをつける（C++14/17でも同様）ほうがよいときもあります。

と言うようにかいておきます。なお、-O3オプションなども同様につきます。

2. Cのソースコードの編集

ループされていて一定以上の処理時間がかかるコードの例として、2つの一行を除けば同じ（無意味な）関数をOMPexample.cとして保存します。

#include<R.h>
#include<Rinternals.h>

SEXP example(SEXP m){
  int i, j, n;
  double s1 = 0, s2 = 0;
  int nor, noc;
  SEXP rv, dim;

  /* 第一引数が行列か確認 */
  if(!isMatrix(m))
    error("A matrix is required for the first argument.");

  /* 行列の行数、列数を得る */
  dim = getAttrib(m, R_DimSymbol);
  nor = INTEGER(dim)[0]; /* 行数 */
  noc = INTEGER(dim)[1]; /* 列数 */

  for(j=0;j<noc;j++){
    for(i=0;i<nor;i++){
      s1 += REAL(m)[i+j*nor]*(i+1)/(j+1); /* 行列と言っても一次元配列 */
      s2 -= REAL(m)[i+j*nor]*(i+1)/(j+1);
    }
  }

  rv = PROTECT(allocVector(REALSXP, 2));
  REAL(rv)[0] = s1;
  REAL(rv)[1] = s2;
  UNPROTECT(1);

  return rv;
}

SEXP OMPexample(SEXP m){
  int i, j, n;
  double s1 = 0, s2 = 0;
  int nor, noc;
  SEXP rv, dim;

  /* 第一引数が行列か確認 */
  if(!isMatrix(m))
    error("A matrix is required for the first argument.");

  /* 行列の行数、列数を得る */
  dim = getAttrib(m, R_DimSymbol);
  nor = INTEGER(dim)[0]; /* 行数 */
  noc = INTEGER(dim)[1]; /* 列数 */

  #pragma omp parallel for private(i, j) reduction(+:s1) reduction(-:s2)
  for(j=0;j<noc;j++){
    for(i=0;i<nor;i++){
      s1 += REAL(m)[i+j*nor]*(i+1)/(j+1); /* 行列と言っても一次元配列 */
      s2 -= REAL(m)[i+j*nor]*(i+1)/(j+1);
    }
  }

  rv = PROTECT(allocVector(REALSXP, 2));
  REAL(rv)[0] = s1;
  REAL(rv)[1] = s2;
  UNPROTECT(1);

  return rv;
}

大学の情報基盤センターなるところが書いたネット上にあるOpenMPの解説を読めばすぐ分かりますが、プリプロセッサ命令

  #pragma omp parallel for private(i, j) reduction(+:s1) reduction(-:s2)

がOpenMPの記述になります。private(i, j)は並列処理で変数iとjは、それぞれ異なる値をシェアしない変数として扱うことを意味し、reduction(+:s1)は、それぞれの並列処理で加算されていくそれぞれの変数s1は、並列処理の終了時に一つに合算されることを意味します。reduction(-:s2)は、加算ではなく減算処理をあわせます。
なお、OpenMPが有効でない場合は、このプリプロセッサは無視されて逐次処理されるコードが出力されます。

3. Cのソースコードをコンパイルする

OMPexample.cを用意したら、以下のようにすると、OMPexample.so（もしくはOMPexample.dll）と言うRの拡張が出来上がります。

R CMD SHLIB OMPexample.c

Rtoolsが入っていないと動かないので注意してください。

4. Rから拡張を呼び出して使う

OpenMPを使ったOMPexample関数と、使っていないexample関数の速度比較をRでしてみましょう。

dll <- dyn.load(paste("OMPexample", .Platform$dynlib.ext, sep = ""))
example <- function(m) .Call("example", m)
OMPexample <- function(m) .Call("OMPexample", m)

set.seed(1002)
n <- 5000 # 一瞬で処理が終わる人は増やすと時間をかけられます
m <- matrix(rnorm(n*n), n, n)

print(system.time({
  cat("No OpenMP: ", example(m), "\n")
}))

print(system.time({
  cat("OpenMP: ", OMPexample(m), "\n")
}))

搭載されているプロセッサ数に応じて速くなります。並列化と言っても、手軽に試せますね。

A. Rcpp/Eigen

Rcppの行列はスレッドセーフではないのでOpenMPの並行処理の中で使えませんが、Eigenは最初にEigen::initParallel();をして行列をスレッドセーフにしておくだけで、OpenMPの並列処理の中で行列を使うことができます。つまり、Rcppの行列からEigenの行列を作成して処理を行うことで、RcppからOpenMPを使うことができます。RcppPararellを使うべきと思うかも知れませんが、RcppPararellはWorkerクラスを定義しないといけないので記述が煩雑な一方、Eigen/OpenMPと比較して速度的なアドバンテージはほぼ無いようです。

*1:OSで用意される環境変数ではないので、Rで Sys.getenv("HOME") をして確認してください。

疎な巨大対称行列の一般化固有値問題のC++/Eigenでの効率的な解き方

「水素原子に対するSchrödinger方程式を有限要素法で数値的に解いてみる（C++のソースコード付き）」で紹介されているC++/Eigenのコードを、（ページ著者が）Juliaに移植してみたら12倍くらい速くなったと言うツイートが流れてきて、気になってC++の方のソースコードを見てみたのですが、

• 巨大な疎行列であること
• 必要な固有値が1つであること*1

を活かしておらず、この2点を利用すればもっと時間を短縮できるものでした*2。実際、そう大きくない変更で、手元の環境では25分の1の時間で計算が終わるようになりました。

追記(2021/02/15 01:42):ブログ主が、有限要素法による水素原子のエネルギー固有値が、小さい方から見てどこまで正確かを検証する部分を加筆しました（図1, 図2）。これによって行列の固有値がすべて必要になるので、以下の方法は適切な計算方法にならないです。

なお、すべての固有値を得る場合にJulia比較でC++/Eigenが遅いのは、どうもC++/Eigenが自動的に並列化してマルチコアを活かしてこないところにあるので、その辺を何らかの方法*3で解消すると高速化できます。

EigenはC++用の線形代数ライブラリなのですが、名前がEigenなのに疎行列（Sparse Matrix）用の固有値問題のソルバーが無かったりします。上述のコードでは、GeneralizedSelfAdjointEigenSolverと言う計算する行列式に、疎行列、対称、正定値と言った条件を使わない上に、行列の固有値をすべて用いようとするソルバーが用いられています。汎用ソルバーとしても遅いというような指摘もあるようですが、入力データの特性を活かすのは高速化のイロハです。

この手の数値演算アルゴリズムを実装しようとすると数学の理解やデバッグの手間隙に大変な労力を費やすことになりますが、幸い、手軽に使えるテンプレート・ライブラリでSpectraと言うのがあります。具体的には、一般化固有値問題の計算対象の2つの行列をEigen::SparseMatrixで持ち、SymGEigsSolverと言うソルバーと解きます。なお、使ってみると癖があるなと思う部分があって、特にコンストラクタの第4引数を小さく取りすぎると固有値が収束しなかったとエラーを出して来るところには注意が必要です。

差分は以下になります。コンパイル時にダウンロードして展開したSpectraへインクルード・パスを通さないとコンパイル・エラーになるので気をつけてください。

--- a/src/hydrogen_fem/hydrogen_fem.cpp
+++ b/src/hydrogen_fem/hydrogen_fem.cpp
@@ -12,19 +12,24 @@
#include <memory> // for std::unique_ptr
#include <boost/assert.hpp> // for BOOST_ASSERT
#include <Eigen/Eigenvalues> // for Eigen::GeneralizedSelfAdjointEigenSolver

+#include <Eigen/Sparse>
+#include <Spectra/SymGEigsSolver.h>
+#include <Spectra/MatOp/SparseGenMatProd.h>
+#include <Spectra/MatOp/SparseCholesky.h>
+#include <iostream>
namespace hydrogen_fem {
// #region コンストラクタ

Hydrogen_FEM::Hydrogen_FEM()
- : hg_(Eigen::MatrixXd::Zero(NODE_TOTAL, NODE_TOTAL)),
+ : hg_(Eigen::SparseMatrix<double>(NODE_TOTAL, NODE_TOTAL)),
length_(ELE_TOTAL),
mat_A_ele_(boost::extents[ELE_TOTAL][2][2]),
mat_B_ele_(boost::extents[ELE_TOTAL][2][2]),
nod_num_seg_(boost::extents[ELE_TOTAL][2]),
node_r_ele_(boost::extents[ELE_TOTAL][2]),
node_r_glo_(NODE_TOTAL),
- ug_(Eigen::MatrixXd::Zero(NODE_TOTAL, NODE_TOTAL))
+ ug_(Eigen::SparseMatrix<double>(NODE_TOTAL, NODE_TOTAL))
{
}

@@ -42,20 +47,33 @@ namespace hydrogen_fem {

// 全体行列を生成
make_global_matrix();

// 一般化固有値問題を解く
- Eigen::GeneralizedSelfAdjointEigenSolver<Eigen::MatrixXd> es(hg_, ug_);

+  Spectra::SparseSymMatProd<double> Aop(hg_);
+  Spectra::SparseCholesky<double> Bop(ug_);
+
+  // 第3引数の1は要求する固有値の数，第4引数は収束速度（とメモリー利用量と行列操作の複雑さ）を定める値
+  Spectra::SymGEigsSolver<double, Spectra::SMALLEST_ALGE, Spectra::SparseSymMatProd<double>, Spectra::SparseCholesky<double>, Spectra::GEIGS_CHOLESKY> geigs(&Aop, &Bop, 1, NODE_TOTAL<10 ? NODE_TOTAL : 10 + (NODE_TOTAL - 10)/100);
+
+  geigs.init();
+  int nconv = geigs.compute();
+  if(Spectra::SUCCESSFUL != geigs.info()){
+    std::cerr << "geigs.info(): " << geigs.info() << std::endl;
+    exit(-1);
+  }

// エネルギー固有値Eを取得
- auto const e = es.eigenvalues()[0];
+ Eigen::VectorXcd evalues = geigs.eigenvalues();
+ auto const e = evalues[0].real(); // 計算結果は複素数なことに注意．なお、複数の固有値を求めた場合、なぜか小さい順に固有値が並んでいるわけではない.

// 固有ベクトル（波動関数）を取得
- phi_ = es.eigenvectors().col(0);
+ phi_ = geigs.eigenvectors(1).col(0).real(); // なぜか添字が1からなのと、計算結果は複素数なのに注意

// 固有ベクトル（波動関数）を規格化
normalize();

return e;
}

void Hydrogen_FEM::save_result() const
@@ -190,8 +208,8 @@ namespace hydrogen_fem {
for (auto e = 0; e < ELE_TOTAL; e++) {
for (auto i = 0; i < 2; i++) {
for (auto j = 0; j < 2; j++) {
- hg_(nod_num_seg_[e][i], nod_num_seg_[e][j]) += mat_A_ele_[e][i][j];
- ug_(nod_num_seg_[e][i], nod_num_seg_[e][j]) += mat_B_ele_[e][i][j];
+ hg_.coeffRef(nod_num_seg_[e][i], nod_num_seg_[e][j]) += mat_A_ele_[e][i][j];
+ ug_.coeffRef(nod_num_seg_[e][i], nod_num_seg_[e][j]) += mat_B_ele_[e][i][j];
}
}
}

--- a/src/hydrogen_fem/hydrogen_fem.h
+++ b/src/hydrogen_fem/hydrogen_fem.h
@@ -15,6 +15,9 @@
#include <Eigen/Core> // for Eigen::MatrixXd
#include <boost/multi_array.hpp> // for boost::multi_array

+#include <Eigen/Core>
+#include <Eigen/SparseCore>
+
namespace hydrogen_fem {
//! A class.
/*!
@@ -144,7 +147,7 @@ namespace hydrogen_fem {
/*!
左辺の全体行列
*/
- Eigen::MatrixXd hg_;
+ Eigen::SparseMatrix<double> hg_;

//! A private member variable.
/*!
@@ -192,7 +195,7 @@ namespace hydrogen_fem {
/*!
右辺の全体行列
*/
- Eigen::MatrixXd ug_;
+ Eigen::SparseMatrix<double> ug_;

// #region 禁止されたコンストラクタ・メンバ関数

*1:すべての固有値を取りたい場合もあるわけですが、大学のレジュメなどを拝見する限りは、最小値か最大値のどちらか一つを計算すれば間に合うケースが多いようです。

*2:最初からすべて分かったように書いていますが、ページ著者と試行錯誤して、この結論にたどり着いています。

*3:OpenBlASを使う方法が紹介されていました。