餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

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]*coef(lm_prc_r)[2:(npc + 1)]/sd(x1) * coef(summary(lm_prc_r))[, 2][2:(npc + 1)])^2))
se["x2", "PCR"] <- sqrt(sum((r_prcomp$rotation[2, ][1:npc]*coef(lm_prc_r)[2:(npc + 1)]/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]*coef(lm_prc_r)[2:(npc + 1)]*mean(x1)/sd(x1)*coef(summary(lm_prc_r))[, 2][2:(npc + 1)])^2)
    + sum((r_prcomp$rotation[2,][1:npc]*coef(lm_prc_r)[2:(npc + 1)]*mean(x2)/sd(x2)*coef(summary(lm_prc_r))[, 2][2:(npc + 1)])^2)
)

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: C_{pcr}を主成分回帰の切片工、 \beta_{pcr}を主成分回帰の回帰式の係数、 r_iを主成分の説明変数 x_iの係数、主成分回帰の誤差項を \etaとすると、主成分回帰は  y = C_{pcr} + \beta_{pcr} \bigg( \sum^2_{i=1} r_i \frac{x_i - \bar{x_i}}{sd(x_i)} \bigg) + \eta と書けるので、式を x_iの項で分けて整理すると、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

f:id:uncorrelated:20201023143052p:plain

これで.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>は要らないので消しました。
(これはチートペーパーに記述の問題はありませんが)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とエラーが出て止まります。

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

コンパイルが終わって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 })

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

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

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

1. Makevarsの編集

gccOpenMPを使うためには、引数に-fopenmpフラグを加える必要があります。R CMDでこの引数をつけるためには、LinuxMacOS Xでは$HOME/.R/Makevarsに、MS-Windowsでは%HOME%/.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と比較して速度的なアドバンテージはほぼ無いようです。

疎な巨大対称行列の一般化固有値問題の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で解消すると高速化できます。

EigenC++用の線形代数ライブラリなのですが、名前が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を使う方法が紹介されていました。

FASTGLを使ってC++で積分

数値解析における数値積分法の代表としてガウス求積と言う技法があり、それに必要なガウス点と重みの計算でJuliaがC++より桁違いに速いと言う話が流れていたのですが、Juliaで使われているFastGaussQuadrature.jlが使っているアルゴリズムC++版を使ったところ、案の定、だいたい同じの速度になりました。

Fast computation of Gauss-Legendre Quadrature Points and Weights(FASTGL)のアルゴリズム開発者のページからfastgl.cppとfastgl.hppをダウンロードしてきて、積分に使うだけです。サブディレクトリfastglにこの2つのファイルを配置したところで、以下のソースコードのファイルを作成してコンパイルしたら簡単に実行できます。-O3 -lmオプションはいると思いますが。

#include <iostream>
#include <chrono>
#include <iomanip>
#include <functional>
#include "fastgl/fastgl.hpp"

/*
  FASTGLを使って積分をする関数
  args: ラムダ関数f(x), xの最小値, xの最大値, 節点の数
*/
double integral(std::function<double(double)> f, double s, double e, int NumberOfKnots){
  double Res = 0.0;
  double xm = 0.5 * (s + e);
  double xr = 0.5 * (e - s);

  for ( int k = 1 ; k <= NumberOfKnots ; ++k ){
    fastgl::QuadPair p = fastgl::GLPair ( NumberOfKnots, k );
    Res += p.weight * f( xm + xr * p.x());
  }

  return Res*xr;
}

int main(){
  std::chrono::system_clock::time_point stime, etime;
  stime = std::chrono::system_clock::now();

  // boost::math::quadratureあたりと異なり、節点が100万個でも大した時間にならない
  double Res = integral([](const double& t) { return 1.0 / 2 / std::sqrt(t); }, 1, 4, 1000000);

  etime = std::chrono::system_clock::now();

  std::cout << std::setprecision(16);
  std::cout << "Results: " << Res << std::endl;
  std::cout << "Elapsted Time: " << std::chrono::duration_cast<std::chrono::milliseconds>(etime - stime).count() << " ms" << std::endl;

  return 0;
}

モンテカルロ法で円周率を計算するベンチマーク

モンテカルロ法で円周率を計算するベンチマークで、Juliaと比較してgccが悪い言われているのですが、本当なのか検証してみましょう。

f:id:uncorrelated:20200927174314p:plain

CとJuliaでこの分岐が単純なマイクロベンチマークを実行すると、乱数生成の速度が結果を大きく左右しえます。後方互換性の問題で、Cの標準の乱数生成関数は線形合同法になっています。Juliaの方は標準がdouble precision SIMD oriented Fast Mersenne Twister pseudorandom number generator (dSFMT)で、明示的にMersenneTwisterを使った場合もdSFMTを用いているようです。Cの方もdSFMTを用いてベンチマークを取ることにします。

Cのソースコード

dSFMTがよくチューニングされているので、私には工夫の余地は無く、オーソドックスに書いています。なお、生成した実行ファイルをreadelfで調べたところ、dSFMTがヘッダーファイルでinline static宣言している関数はインライン展開されていました。10回計測して、一番短い経過時間を採用しています。

#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <limits.h>
#include "dSFMT.h"
#include <float.h>
#include <math.h>

/*
dSFMT-src-2.2.3.tgz を展開したディレクトリで、ソースコードのファイル名をpi_mc_dsfmt.cとして、

clang -o pi_mc pi_mc_dsfmt.c dSFMT-src-2.2.3/dSFMT.c -IdSFMT-src-2.2.3/ -DDSFMT_MEXP=521 -msse2 -DHAVE_SSE2 -O3 -lm

コンパイルしました。gccとclangはオプションが同じなので、gccとの比較は容易に行えます。
*/

double pi_mc(long n_iter) {
  double  x, y;
  dsfmt_t rng;

  dsfmt_init_gen_rand(&rng, time(NULL));

  int n_in = 0;
  for (long i = 1; i <= n_iter; i++) {

    x = dsfmt_genrand_open_open(&rng);
    y = dsfmt_genrand_open_open(&rng);

    if (x * x + y * y < 1.0)
      n_in++;
  }

  return (double)4.0*n_in/n_iter;
}

int main(){
  struct timespec t1, t2;
  double  dt, sum_dt, min_dt;
  double  r = 0;
  size_t  n;
  int  j, nos = 10;

  printf("log10(n)\telapsed time (seconds)\n");

  for(n=1; n<=1000000*1000; n*=10){
    sum_dt = 0;
    min_dt = DBL_MAX;
    for(j=0; j<nos; j++){
      clock_gettime(CLOCK_THREAD_CPUTIME_ID, &t1); // clock()を使って計測しても、ほぼ同じ数字になります

      r = pi_mc(n);

      clock_gettime(CLOCK_THREAD_CPUTIME_ID, &t2);

//      fprintf(stderr, "estimated pi: %f\n", r);

      dt = (double)(t2.tv_sec - t1.tv_sec) + (double)(t2.tv_nsec - t1.tv_nsec)/1000000000;
      
      sum_dt += dt;
      if(min_dt > dt)
        min_dt = dt;
    }

    printf("%.0f\t%.09f\n", log10(n), min_dt);
  }

  return 0;  
}

Juliaのスクリプト

Cのとだいたい同じですが、グローバル変数へのアクセスが無いように一番外側ループ部分も関数内にしまいました。

using Random
# using Random: default_rng

function pi_mc(n_iter)
n_in = 0
# rng = default_rng()
rng = MersenneTwister(time_ns())
for i in 1:n_iter
x, y = rand(rng), rand(rng)
n_in += ifelse(x^2 + y^2 ≤ 1, 1, 0)
end
4n_in/n_iter
end

function loop(s, e)
  n = s
  while(n <= e)
    sum_dt = 0;
    min_dt = 3600;
    for i in 1:10
      t1 = time_ns()
      r = pi_mc(n)
      t2 = time_ns()
      dt = (t2 - t1)/1000000000
      sum_dt = sum_dt + dt
      if(min_dt > dt)
        min_dt = dt
      end
#      println("estimated pi: ", r)
    end
    println(log10(n), "\t", min_dt)
    n = n*10
  end
end

println("log10(n)", "\t", "elapsed time (seconds)")

loop(1, 10^9)