2019-03-31

閉殻Hartree-Fock法によるHeH+のエネルギー計算

量子化学計算を勉強するために新しい量子化学―電子構造の理論入門〈上〉を読んでいます。
サンプルコードがFortranだったので勉強がてらC言語に移植しました。

新しい量子化学―電子構造の理論入門〈上〉

作者:Attila Szabo,Neil S. Ostlund
出版社/メーカー: 東京大学出版会
発売日: 1987/07/01
メディア: 単行本

理論

本で200ページくらいかけて説明されている内容をブログに書くのは大変なので省略します。

ちなみに、以下のツイートの「100ページかけて説明してる別の本」はこの本を指します。
難しい内容を易しく簡潔に説明しようとするのは大変で、難しい内容は難しいまま長い時間をかけて説明する必要があるのだなぁという学びがありました。

ある本で30ページで説明されてて何も分からなかった内容を100ページで説明してる別の本読んだら完全に理解して何が分からなかったのか分からなくなった
— きばん (@kivantium) 2019年3月21日

実装

本のFortran実装を素直にC言語で書き直しました。
本のコードを打ち込むのが面倒だったので、検索して見つけたページにあったコードと見比べながら書きました。

最終的に得られたエネルギーはFortran実装の結果と小数第5位まで一致しました。
この差がどこで生じているのか謎ですが、実装は正しいんじゃないでしょうか……

#include <math.h>
#include <stdio.h>

// プロトタイプ宣言
void HFCALC(int, int, double, double, double, double, double);
void INTGRL(int, int, double, double, double, double, double);
void COLECT(int, int, double, double, double, double, double);
void SCF(int, int, double, double, double, double, double);
double calcS(double, double, double);
double calcT(double, double, double);
double F0(double);
double calcV(double, double, double, double, double);
double TWOE(double, double, double, double, double, double, double);
void FORMG(void);
void DIAG(double [][2], double [][2], double [][2]);
void MULT(double [][2], double [][2], double [][2], int, int);
void MATOUT(double [][2], int, int, int, int, char*);

double S12, T11, T12, T22, V11A, V12A, V22A, V11B, V12B, V22B,
       V1111, V2111, V2121, V2211, V2221, V2222;

// Slater型関数に対する最小2乗近似の短縮係数と軌道指数
double COEF[3][3] = {{1.0, 0.678914, 0.444635}, 
                     {0.0, 0.430129, 0.535328},
                     {0.0, 0.0, 0.154329}};
double EXPON[3][3] = {{0.27095, 0.151623, 0.109818},
                      {0.0, 0.851819, 0.405771},
                      {0.0, 0.0, 2.22766}};
double D1[3], A1[3], D2[3], A2[3];
double PI = 3.1415926535898;

double S[2][2], X[2][2], XT[2][2], H[2][2], F[2][2], G[2][2], C[2][2],
       FPRIME[2][2], CPRIME[2][2], P[2][2], OLDP[2][2], TT[2][2][2][2], E[2][2];

int main(void) {
  int IOP = 2;
  // 基底関数STO-NG
  int N = 3;
  // 基底状態の平衡結合長(a.u.)
  double R = 1.4632;
  // Heに対するSTO-3G軌道指数値
  // He原子に対する最良の軌道指数は1.6875
  // HeH+では電子雲が縮んでいることを反映して1.24倍
  double ZETA1 = 2.0925;
  // Hに対するSTO-3G軌道指数値
  // H原子に対する最良の軌道指数は1.0
  // HeH+では電子雲が縮んでいることを反映して1.24倍
  double ZETA2 = 1.24;
  double ZA = 2.0;
  double ZB = 1.0;
  HFCALC(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
}

/*
 * 1s基底関数に対応するSTO-NG最小基底を用いて2電子系2原子分子に対する
 * Hartree-Fock計算を行う関数
 *
 * IOP=0 何も表示しない
 * IOP=1 収束の結果のみ表示
 * IOP=2 全繰り返しの結果を表示
 * N STO-NG を使う (N=1, 2, 3)
 * R 結合長 (a.u.)
 * ZETA1 原子Aの基底関数のSlater軌道指数値
 * ZETA2 原子Bの基底関数のSlater軌道指数値
 * ZA 原子Aの原子番号
 * ZB 原子Bの原子番号
 */
void HFCALC(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
            double ZB) {
  if (IOP != 0) {
    printf("   STO-%dG for atomic numbers %5.2lf and %5.2lf\n\n\n", N, ZA, ZB);
  }
  INTGRL(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
  COLECT(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
  SCF(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
}

/* SCF計算に必要な積分を計算する関数*/
void INTGRL(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
            double ZB) {
  double R2 = R * R;
  // 軌道指数のスケーリングと短縮係数の規格化
  for (int i = 0; i < N; ++i) {
    A1[i] = EXPON[i][N-1] * pow(ZETA1, 2.0);
    D1[i] = COEF[i][N-1] * pow(2.0 * A1[i] / PI, 0.75);
    A2[i] = EXPON[i][N-1] * pow(ZETA2, 2.0);
    D2[i] = COEF[i][N-1] * pow(2.0 * A2[i] / PI, 0.75);
  }
  S12 = 0.0;
  T11 = 0.0;
  T12 = 0.0;
  T22 = 0.0;
  V11A = 0.0;
  V12A = 0.0;
  V22A = 0.0;
  V11B = 0.0;
  V12B = 0.0;
  V22B = 0.0;
  V1111 = 0.0;
  V2111 = 0.0;
  V2121 = 0.0;
  V2211 = 0.0;
  V2221 = 0.0;
  V2222 = 0.0;

  // 1電子積分の計算
  // V12A = 原子核Aへの各引力の非対角項
  // RAP2 = 中心Aと中心Pの距離の2乗
  for (int i = 0; i < N; ++i) {
    for (int j = 0; j < N; ++j) {
      double RAP = A2[j] * R / (A1[i] + A2[j]);
      double RAP2 = pow(RAP, 2.0);
      double RBP2 = pow(R - RAP, 2.0);
      S12 += calcS(A1[i], A2[j], R2) * D1[i] * D2[j];
      T11 += calcT(A1[i], A1[j], 0.0) * D1[i] * D1[j];
      T12 += calcT(A1[i], A2[j], R2) * D1[i] * D2[j];
      T22 += calcT(A2[i], A2[j], 0.0) * D2[i] * D2[j];
      V11A += calcV(A1[i], A1[j], 0.0, 0.0, ZA) * D1[i] * D1[j];
      V12A += calcV(A1[i], A2[j], R2, RAP2, ZA) * D1[i] * D2[j];
      V22A += calcV(A2[i], A2[j], 0.0, R2, ZA) * D2[i] * D2[j];
      V11B += calcV(A1[i], A1[j], 0.0, R2, ZB) * D1[i] * D1[j];
      V12B += calcV(A1[i], A2[j], R2, RBP2, ZB) * D1[i] * D2[j];
      V22B += calcV(A2[i], A2[j], 0.0, 0.0, ZB) * D2[i] * D2[j];
    }
  }
  // 2電子積分の計算
  for (int i = 0; i < N; ++i) {
    for (int j = 0; j < N; ++j) {
      for (int k = 0; k < N; ++k) {
        for (int l = 0; l < N; ++l) {
          double RAP = A2[i] * R / (A2[i] + A1[j]);
          //double RBP = R - RAP;
          double RAQ = A2[k] * R / (A2[k] + A1[l]);
          double RBQ = R - RAQ;
          double RPQ = RAP - RAQ;
          double RAP2 = RAP * RAP;
          //double RBP2 = RBP * RBP;
          //double RAQ2 = RAQ * RAQ;
          double RBQ2 = RBQ * RBQ;
          double RPQ2 = RPQ * RPQ;
          V1111 += TWOE(A1[i], A1[j], A1[k], A1[l], 0.0, 0.0, 0.0) 
                   * D1[i] * D1[j] * D1[k] * D1[l];
          V2111 += TWOE(A2[i], A1[j], A1[k], A1[l], R2, 0.0, RAP2)
                   * D2[i] * D1[j] * D1[k] * D1[l];
          V2121 += TWOE(A2[i], A1[j], A2[k], A1[l], R2, R2, RPQ2)
                   * D2[i] * D1[j] * D2[k] * D1[l];
          V2211 += TWOE(A2[i], A2[j], A1[k], A1[l], 0.0, 0.0, R2)
                   * D2[i] * D2[j] * D1[k] * D1[l];
          V2221 += TWOE(A2[i], A2[j], A2[k], A1[l], 0.0, R2, RBQ2)
                   * D2[i] * D2[j] * D2[k] * D1[l];
          V2222 += TWOE(A2[i], A2[j], A2[k], A2[l], 0.0, 0.0, 0.0) * D2[i] *
                   D2[j] * D2[k] * D2[l];
        }
      }
    }
  }
  if (IOP != 0) {
    printf("   R          ZETA1      ZETA2      S12       T11\n\n");
    printf("%11.6lf%11.6lf%11.6lf%11.6lf%11.6lf\n\n\n", R, ZETA1, ZETA2, S12, T11);
    printf("   T12        T22        V11A       V12A        V22A\n\n");
    printf("%11.6lf%11.6lf%11.6lf%11.6lf%11.6lf\n\n\n", T12, T22, V11A, V12A, V22A);
    printf("   V11B       V12B       V22B       V1111        V2111\n\n");
    printf("%11.6lf%11.6lf%11.6lf%11.6lf%11.6lf\n\n\n", V11B, V12B, V22B, V1111,
           V2111);
    printf("   V2121      V2211       V2221     V2222\n\n");
    printf("%11.6lf%11.6lf%11.6lf%11.6lf\n", V2121, V2211, V2221, V2222);
  }
}

// 重なり積分
double calcS(double A, double B, double RAB2) {
  return pow(PI / (A + B), 1.5) * exp(-A * B * RAB2 / (A + B));
}

// 運動エネルギー積分
double calcT(double A, double B, double RAB2) {
  return A * B / (A + B) * (3.0 - 2.0 * A * B * RAB2 / (A + B)) *
         pow(PI / (A + B), 1.5) * exp(-A*B*RAB2/(A+B));
}

// F0関数
double F0(double ARG) {
  double ret;
  if (ARG < 1.0e-6) {
    ret = 1.0 - ARG / 3.0;
  } else {
    ret = sqrt(PI / ARG) * erf(sqrt(ARG)) / 2.0;
  }
  return ret;
}

// 核引力積分
double calcV(double A, double B, double RAB2, double RCP2, double ZC) {
  double v = 2.0 * PI / (A + B) * F0((A + B) * RCP2) * exp(-A * B * RAB2 / (A + B));
  return -v * ZC;
}

// 2電子積分
double TWOE(double A, double B, double C, double D, double RAB2, double RCD2,
            double RPQ2) {
  return 2.0 * pow(PI, 2.5) / ((A + B) * (C + D) * sqrt(A + B + C + D))
         * F0((A + B) * (C + D) * RPQ2 / (A + B + C + D))
         * exp(-A * B * RAB2 / (A + B) - C * D * RCD2 / (C + D));
}

// 積分の結果から行列を求める
void COLECT(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
            double ZB) {
  // 核ハミルトニアン cf. 式(3.153)
  H[0][0] = T11 + V11A + V11B;
  H[0][1] = T12 + V12A + V12B;
  H[1][0] = H[0][1];
  H[1][1] = T22 + V22A + V22B;

  // 重なり行列 cf. 式(3.136)
  S[0][0] = 1.0;
  S[0][1] = S12;
  S[1][0] = S[0][1];
  S[1][1] = 1.0;

  // 正準直交化 cf. 式(3.169), p.191
  X[0][0] = 1.0 / sqrt(2.0 * (1.0 + S12));
  X[1][0] = X[0][0];
  X[0][1] = 1.0 / sqrt(2.0 * (1.0 - S12));
  X[1][1] = -X[0][1];

  // 変換行列の転置
  XT[0][0] = X[0][0];
  XT[0][1] = X[1][0];
  XT[1][0] = X[0][1];
  XT[1][1] = X[1][1];

  // 2電子積分の行列
  TT[0][0][0][0] = V1111;
  TT[1][0][0][0] = V2111;
  TT[0][1][0][0] = V2111;
  TT[0][0][1][0] = V2111;
  TT[0][0][0][1] = V2111;
  TT[1][0][1][0] = V2121;
  TT[0][1][1][0] = V2121;
  TT[1][0][0][1] = V2121;
  TT[0][1][0][1] = V2121;
  TT[1][1][0][0] = V2211;
  TT[0][0][1][1] = V2211;
  TT[1][1][1][0] = V2221;
  TT[1][1][0][1] = V2221;
  TT[1][0][1][1] = V2221;
  TT[0][1][1][1] = V2221;
  TT[1][1][1][1] = V2222;

  if (IOP != 0) {
    MATOUT(S, 2, 2, 2, 2, "S   ");
    MATOUT(X, 2, 2, 2, 2, "X   ");
    MATOUT(H, 2, 2, 2, 2, "H   ");
    printf("\n");
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        for (int k = 0; k < 2; ++k) {
          for (int l = 0; l < 2; ++l) {
            printf("   (%2d%2d%2d%2d )%10.6lf\n", i+1, j+1, k+1, l+1, TT[i][j][k][l]);
          }
        }
      }
    }
  }
}

// SCFを実行する
void SCF(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
         double ZB) {
  // 密度行列の収束規準
  double CRIT = 1.0e-4;
  // 繰り返しの最大回数
  int MAXIT = 25;
  // 核ハミルトニアンを最初の推定に用いる（つまりPを零行列で初期化）
  for (int i = 0; i < 2; ++i) {
    for (int j = 0; j < 2; ++j) {
      P[i][j] = 0.0;
    }
  }
  if (IOP == 2) {
    MATOUT(P, 2, 2, 2, 2, "P    ");
  }
  int ITER = 0;
  double EN;
  for (;;) {
    ITER++;
    if (IOP == 2) {
      printf("\n    START OF ITERATION NUMBER = %2d\n", ITER);
    }
    // Fock行列の2電子部分を計算
    FORMG();

    if (IOP == 2) {
      MATOUT(G, 2, 2, 2, 2, "G    ");
    }

    // 核ハミルトニアンを足してFock行列を得る
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        F[i][j] = H[i][j] + G[i][j];
      }
    }

    // 電子エネルギーを計算する
    EN = 0.0;
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        EN += 0.5 * P[i][j] * (H[i][j] + F[i][j]);
      }
    }

    if (IOP == 2) {
      MATOUT(F, 2, 2, 2, 2, "F   ");
      printf("\n\n\n    ELECTRONIC ENERGY = %20.12lf\n", EN);
    }

    // Gを用いてFock行列を変換
    MULT(F, X, G, 2, 2);
    MULT(XT, G, FPRIME, 2, 2);
    // Fock行列の対角化
    DIAG(FPRIME, CPRIME, E);
    // 固有ベクトルから行列Cを得る
    MULT(X, CPRIME, C, 2, 2);
    // 新しい密度行列を計算する
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        OLDP[i][j] = P[i][j];
        P[i][j] = 0.0;
        for (int k = 0; k < 1; ++k) { // ここはk<1で本当にいいのか不明
          P[i][j] += 2.0 * C[i][k] * C[j][k];
        }
      }
    }
    if (IOP == 2) {
      MATOUT(FPRIME, 2, 2, 2, 2, "F'   ");
      MATOUT(CPRIME, 2, 2, 2, 2, "C'   ");
      MATOUT(E, 2, 2, 2, 2, "E   ");
      MATOUT(C, 2, 2, 2, 2, "C   ");
      MATOUT(P, 2, 2, 2, 2, "P   ");
    }
    // deltaを計算する
    double DELTA = 0.0;
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        DELTA += pow(P[i][j] - OLDP[i][j], 2.0);
      }
    }
    DELTA = sqrt(DELTA / 4.0);
    if (IOP != 0) {
      printf("\n    DELTA(CONVERGENCE OF DENSITY MATRIX) = %10.6lf\n", DELTA);
    }

    // 収束判定
    if (DELTA < CRIT) break;
    if (ITER < MAXIT) continue;
    // 収束しなかった場合
    printf("    NO CONVERGENCE IN SCF");
    return;
  }

  double ENT = EN + ZA * ZB / R;
  if (IOP != 0) {
    printf("\n\n    CALCULATION CONVERGED"
           "\n\n    ELECTORONIC ENERGY = %20.12lf"
           "\n\n    TOTAL ENERGY =       %20.12f", EN, ENT);
  }
  if (IOP == 1) {
    MATOUT(G, 2, 2, 2, 2, "G   ");
    MATOUT(F, 2, 2, 2, 2, "F   ");
    MATOUT(E, 2, 2, 2, 2, "E   ");
    MATOUT(C, 2, 2, 2, 2, "C   ");
    MATOUT(P, 2, 2, 2, 2, "P   ");
  }
  MULT(P, S, OLDP, 2, 2);
  if (IOP != 0) {
    MATOUT(OLDP, 2, 2, 2, 2, "PS  ");
  }
}

void FORMG(void) {
  for (int i = 0; i < 2; ++i) {
    for (int j = 0; j < 2; ++j) {
      G[i][j] = 0.0;
      for (int k = 0; k < 2; ++k) {
        for (int l = 0; l < 2; ++l) {
          G[i][j] += P[k][l] * (TT[i][j][k][l] - 0.5 * TT[i][l][k][j]);
        }
      }
    }
  }
}

// Fを対角化し、Cに固有ベクトル、Eに固有値を入れる
void DIAG(double F[][2], double C[][2], double E[][2]) {
  double theta;
  if (fabs(F[0][0] - F[1][1]) > 1.0e-20) {
    theta = 0.5 * atan(2.0 * F[0][1] / (F[0][0] - F[1][1]));
  } else {
    theta = PI / 4.0;
  }
  C[0][0] = cos(theta);
  C[1][0] = sin(theta);
  C[0][1] = sin(theta);
  C[1][1] = -cos(theta);
  E[0][0] = F[0][0] * pow(cos(theta), 2.0) + F[1][1] * pow(sin(theta), 2.0)
            + F[0][1] * sin(2.0 * theta);
  E[1][1] = F[1][1] * pow(cos(theta), 2.0) + F[0][0] * pow(sin(theta), 2.0)
            - F[0][1] * sin(2.0 * theta);
  E[1][0] = 0.0;
  E[0][1] = 0.0;
  // 固有値と固有ベクトルを整列する
  if (E[1][1] > E[0][0]) return;
  double temp = E[1][1];
  E[1][1] = E[0][0];
  E[0][0] = temp;
  temp = C[0][1];
  C[0][1] = C[0][0];
  C[0][0] = temp;
  temp = C[1][1];
  C[1][1] = C[1][0];
  C[1][0] = temp;
}

// AとBの積をCに入れる
void MULT(double A[][2], double B[][2], double C[][2], int IM, int M) {
  for (int i = 0; i < M; ++i) {
    for (int j = 0; j < M; ++j) {
      C[i][j] = 0.0;
      for (int k = 0; k < M; ++k) {
        C[i][j] += A[i][k] * B[k][j];
      }
    }
  }
}

void MATOUT(double A[][2], int IM, int IN, int M, int N, char *LABEL) {
  int IHIGH = 0;
  for (;;) {
    int LOW = IHIGH + 1;
    IHIGH = IHIGH + 5;
    IHIGH = (IHIGH < N) ? IHIGH : N;
    printf("\n\n\n    THE %s ARRAY\n", LABEL);
    printf("               ");
    for (int i = LOW; i <= IHIGH; ++i) {
      printf("          %3d      ", i);
    }
    printf("\n");
    for (int i = 1; i <= M; ++i) {
      printf("%10d     ", i);
      for (int j = LOW; j <= IHIGH; ++j) {
        printf(" %18.10lf", A[i - 1][j - 1]);
      }
      printf("\n");
    }
    if (N <= IHIGH) break;
  }
}

最終的な出力は

    CALCULATION CONVERGED

    ELECTORONIC ENERGY =      -4.227529304014

    TOTAL ENERGY =            -2.860662163500


    THE PS   ARRAY
                           1                  2      
         1            1.5296370381       1.1199277981
         2            0.6424383868       0.4703629619

でした。参照したFortranコードの出力は

    CALCULATION CONVERGED

    ELECTRONIC ENERGY =  -0.422752913203D+01

    TOTAL ENERGY =       -0.286066199152D+01



    THE PS   ARRAY
                           1                  2
         1        0.1529637121D+01   0.1119927796D+01
         2        0.6424383099D+00   0.4703628793D+00

でした。

FortranからCに移植するときのtips

詰まったことを挙げておきます。

配列の添字

Fortranでは添字が1から始まるので配列にアクセスするときに添字をうまく調整する必要があります。

多次元配列の並び順

このページに詳しく書いてありますが、C言語の多次元配列はメモリ上でa[0][0], a[0][1], a[1][0], a[1][1], a[2][0], a[2][1]のような順番で並びますが、Fortranではa(1,1), a(2,1), a(3,1), a(1,2), a(2,2), a(3,2)のような順番に並びます。多次元配列の初期化を書き直すときに気をつける必要があります。

整合配列

Fortranではサブルーチンに配列を渡した配列のサイズをサブルーチン側で定義することができます。（整合配列という）

      dimension a(6)
      a(1) = 1
      a(2) = 2
      a(3) = 3
      a(4) = 4
      a(5) = 5
      a(6) = 6
      call sub(a,3,2)
      end

      subroutine sub(a,n,m)
      dimension a(n,m) <---- a(n,*) でも同じ結果を得る
      do i=1,n
         write(6,*) (a(i,j),j=1,m)
      end do
      return
      end

コードはFortran プログラミングの基礎知識から引用しました。

C99で導入された可変長配列の機能を使うと同じことができるそうです。cf. 3.3 整合配列、可変長配列
今回は配列のサイズが決め打ちできたのでvoid MULT(double A[][2], double B[][2], double C[][2], int IM, int M) { のように関数を定義しました。

C言語のabs

特にFortranとは関係がないのですが、C言語のabsはint型に対して使う関数で、double型に対してはfabsという別の関数が用意されています。（これはC言語で関数オーバーロードができないことが原因です。C11で導入されたGenericという機能を使うとそれっぽいことが実現できるらしいです。cf. Generic - C言語入門）

2017-07-16

TwitterにMP4動画をアップロードするにはyuv420pを使う必要がある（らしい）

以前の記事でMP4をアニメーションGIFに変換する方法を紹介しました。
kivantium.hateblo.jp

その後Twitterの機能変更でMP4の動画がそのままアップロードできるようになったのですが、アップロードに失敗することが非常に多かったです。

Twitterの公式ドキュメント（How to share and watch videos on Twitter）には

モバイルアプリではMP4とMOVの動画形式をサポートしています。
ブラウザではMP4（H264形式、AACオーディオ）をサポートしています。アップロードできる動画のサイズは最大512MBです。ただし長さは2分20秒間以下にしてください。

最小解像度: 32 x 32
最大解像度: 1920 x 1200（および1200 x 1900）
縦横比: 1:2.39～2.39:1の範囲（両方の値を含む）
最大フレームレート: 40fps
最大ビットレート: 25Mbps

としか書いていないのですが、実際にはフォーマットをyuv420pにしないとアップロードできないらしいです。[要出典]
（そのように書いてあるTwitter公式のドキュメントは見当たらないのですが、多くのサイトでそう指摘されていました。）

今まで失敗していたのはffmpegのデフォルトだとyuv444になってしまうからのようです。

H.264で作ったのにアップロードに失敗した動画をyuv420pにするのは簡単で、

ffmpeg -i input.mp4 -pix_fmt yuv420p output.mp4

とするだけです。

コーデックの変換も含めて行うには

ffmpeg -i input.webm -vcodec libx264 -pix_fmt yuv420p -strict -2 -acodec aac output.mp4

のようにすればいいようです。

参考

2016-03-25

ボトムアップ型自動微分の実験

自動微分

Deep Learningの関係で自動微分が脚光を浴びつつあるような気がしますが、自動微分を解説したページは少なくまだまだマイナーな分野だと思います。先日ようやく「アルゴリズムの自動微分と応用」を一通り眺めたのでいろいろ実験して遊んでいます。今日は自動微分のうち、ボトムアップ型自動微分をオペレータオーバーロードを用いて実現する方法について書きます。

自動微分とは

自動微分は、アルゴリズムによって定義された関数からその関数の偏導関数値を計算するアルゴリズムを導出するための技術です。一般的には $y=x^2$ のような関数が先にあって、その関数を計算するアルゴリズムやプログラムがあるというように考えますが、自動微分ではどちらかというとアルゴリズムが先にあってアルゴリズムが表現する関数が生まれるというような考え方をします。

プログラムで微分を扱う上でよく知られている技術には数値微分と数式微分があります。

数値微分は、小さな値 $h$ を用いて $\frac{f(a+h)-f(a)}{h}$ や $\frac{f(a+h)-f(a-h)}{2h}$ を計算することで微分の近似値を得る技術です。数値微分はあくまで微分の近似値を計算する技術ですが、自動微分ではアルゴリズムで定義された関数から解析的に得た偏導関数の値を計算します。また、数値微分では摂動幅 $h$ を設定する必要がありますが、 $h$ が大きすぎると近似精度が悪くなり、小さすぎると丸め誤差の影響で精度が悪くなるという問題があります。

この問題を理解するために簡単なプログラムを書きました。 $x=\frac{\pi}{4}$ における $\sin(x)$ の微分を $\frac{f(a+h)-f(a)}{h}$ を計算することで求めます。解析的には $\frac{\sqrt{2}}{2}$ が解なので、この2倍を出力するようにしてよく知っている $\sqrt{2}$ の値にどれくらい近づくかを見てみます。

#include <iostream>
#include <iomanip>
#include <cmath>

using namespace std;

int main(void){
    // 小数点以下15桁表示する
    cout << fixed << setprecision(15);
    
    // 刻み幅を小さくして数値微分
    double h = 0.1;
    for(int i=0; i<15; ++i) {
        cout << 2 * (sin(M_PI/4 + h) - sin(M_PI/4))/h << " (h=" << h << ")" << endl;
        h /= 10;
    }
}

結果は

1.341205945807979 (h=0.100000000000000)
1.407118983378419 (h=0.010000000000000)
1.413506219948735 (h=0.001000000000000)
1.414142849338607 (h=0.000100000000000)
1.414206491290315 (h=0.000010000000000)
1.414212855488372 (h=0.000001000000000)
1.414213490757987 (h=0.000000100000000)
1.414213568473599 (h=0.000000010000000)
1.414213635086980 (h=0.000000001000000)
1.414215411443819 (h=0.000000000100000)
1.414224293228016 (h=0.000000000010000)
1.414202088767524 (h=0.000000000001000)
1.414424133372449 (h=0.000000000000100)
1.421085471520200 (h=0.000000000000010)
1.554312234475219 (h=0.000000000000001)

のようになりました。大きすぎても小さすぎても良くなくて、h=0.000000010000000のときの結果が最も良いようです。

一般に、 $\frac{f(a+h)-f(a)}{h}$ の計算誤差をなるべく小さくするためには関数の計算値に含まれる計算誤差を $\delta$ として $h \simeq 2\sqrt{\frac{\delta}{|f''(a)|}}$ とするのがよいことが知られています。しかし、 $f'$ を求めるのに $f''$ の値が必要になってしまうのでこの式を使うのは現実的ではありません。

また、入力変数が多い関数に対しては入力変数の数だけ関数を計算する必要があり計算量が増えてしまいますが、自動微分を使うことで計算量を抑えることができます。

数式微分

MathematicaやMaximaのようなソフトでは数式を処理して偏導関数を求める数式微分を行いますが、自動微分では繰り返しや分岐を含むようなアルゴリズムに対しても偏導関数を求めるアルゴリズム（偏導関数そのものではない）を求めます。

自動微分の種類

自動微分にはボトムアップ型（フォーワードモードとも）とトップダウン型（リバースモードとも）の2種類があります。どちらも微分の連鎖律を用いて関数を単純な演算に分解することで偏導関数値を求めますが、分解の方法が異なります。

ボトムアップ型自動微分

ボトムアップ型自動微分では一つの入力変数に対して、全ての中間変数の偏導関数値を計算していく手法です。

例としてWikipediaの自動微分のページにある $f(x_1, x_2)=x_1 x_2+\sin{x_1}$ の例を考えます。この関数の計算をステップごとに中間変数に保存するとすれば、

$w_1 \leftarrow x_1$
$w_2 \leftarrow x_2$
$w_3 \leftarrow w_1 \cdot w_2$
$w_4 \leftarrow \sin{w_1}$
$w_5 \leftarrow w_3 + w_4$

のように5つの中間変数を用いて表すことができます。

ここで入力変数 $x_1$ について各中間変数の微分を計算すると、

$\frac{\partial w_1}{\partial x_1} = 1$

$\frac{\partial w_2}{\partial x_1} = 0$

$\frac{\partial w_3}{\partial x_1}= \frac{\partial (w_1 \cdot w_2 )}{\partial x_1} = \frac{\partial w_1}{\partial x_1}\cdot w_2 + w_1 \cdot \frac{\partial w_2}{\partial x_1}$

$\frac{\partial w_4}{\partial x_1} = \frac{\partial \sin{w_1}}{\partial x_1} = \frac{\partial w_1}{\partial x_1} \cdot \cos{w_1}$

$\frac{\partial w_5}{\partial x_1} = \frac{\partial (w_3+w_4)}{\partial x_1} = \frac{\partial w_3}{\partial x_1} + \frac{\partial w_4}{\partial x_1}$

のように分解することができます。これは関数値を求める際に必要な各中間変数の計算と同時に行うことができます。

計算グラフを使って表すと以下のようになります。（Wikipediaから引用）
f:id:kivantium:20160325000216p:plain:w600

トップダウン型自動微分

トップダウン型自動微分では一つの出力変数について、全ての中間変数に対する偏導関数値を計算していく手法です。

先ほどと同じ例を考えます。出力変数 $w_5$ について各中間変数に対する微分を計算すると、

$\frac{\partial w_5}{\partial w_5} = 1$

$\frac{\partial w_5}{\partial w_4} = \frac{\partial w_5}{\partial w_5}\cdot \frac{\partial w_5}{\partial w_4} = \frac{\partial w_5}{\partial w_5}\cdot\frac{\partial (w_3+w_4)}{\partial w_4}=\frac{\partial w_5}{\partial w_5}$

$\frac{\partial w_5}{\partial w_3} = \frac{\partial w_5}{\partial w_5}\cdot \frac{\partial (w_3+w_4)}{\partial w_3} = \frac{\partial w_5}{\partial w_5}$

$\frac{\partial w_5}{\partial w_2} = \frac{\partial w_5}{\partial w_3} \cdot \frac{\partial w_3}{\partial w_2} = \frac{\partial w_5}{\partial w_3} \cdot \frac{\partial (w_1 \cdot w_2)}{\partial w_2} = \frac{\partial w_5}{\partial w_3}\cdot w_1$

$\frac{\partial w_5}{\partial w_1} = \frac{\partial w_5}{\partial w_3} \cdot \frac{\partial w_3}{\partial w_1} + \frac{\partial w_5}{\partial w_4} \cdot \frac{\partial w_4}{\partial w_1}= \frac{\partial w_5}{\partial w_3} \cdot \frac{\partial (w_1 \cdot w_2)}{\partial w_1}+ \frac{\partial w_5}{\partial w_4} \cdot \frac{\partial \sin{w_1}}{\partial w_1} = \frac{\partial w_5}{\partial w_3} \cdot w_2 + \frac{\partial w_5}{\partial w_4} \cdot \cos{w_1}$

のように分解することができます。

$w_1, \ w_2$ についての分解が分かりにくいですが、 $w_5$ の変化は中間変数 $w_3, \ w_4$ の変化を経由してしか起こらないことに注目すると納得できるかもしれません。計算グラフは以下の通りです。
f:id:kivantium:20160325002547p:plain:w600

トップダウン型自動微分では一度計算した後で各中間変数がどのように計算されたのかの履歴を逆向きに辿って偏導関数値を計算する必要があるため、計算グラフを保存するためのメモリが必要になります。

勾配を求める際には、入力変数の方が多いときはトップダウン型、出力変数の方が多い時はボトムアップ型を使うのが計算量の点で有利になります。ニューラルネットワークのバックプロパゲーションはトップダウン型の自動微分の一種です。

自動微分の実現方法

自動微分の実現方法には、ソースコードを解析して偏導関数値を求めるプログラムを作る方法とオペレータオーバーロードによる方法があります。ソースコードを解析する方法はコンパイラを作るのと同じような労力が必要になるので難しいですが、オペレータオーバーロードによる方法は比較的容易なので今回はこちらを採用します。

オペレータオーバーロードによるボトムアップ型自動微分

ボトムアップ型自動微分では関数値を計算するのと同時に偏導関数値を計算することができます。そのためdoubleを自動微分用に拡張したBUdoubleというクラスをつくって計算を行うと同時に偏導関数値を保存するようにします。それを実現したのが以下のコードです。

#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;
 
// ボトムアップ型微分積分用doubleクラス
class BU_double {
    // 変数値
    double val;
    // 偏導関数値
    double d_val;
public:
    // コンストラクタ
    BU_double(double v=0, double dv=0) {
        val = v;
        d_val = dv;
    }
    // 微分する入力変数として選択する関数
    void select(void) {
        d_val = 1.0;
    }
    // 変数値を返す
    double get_value(void) {
        return val;
    }
    // 偏導関数値を返す
    double get_d_value(void) {
        return d_val;
    }

    // 各種演算子の定義
    friend BU_double operator + (BU_double x, BU_double y) {
        return BU_double(x.val+y.val, x.d_val+y.d_val);
    }
    friend BU_double operator - (BU_double x, BU_double y) {
        return BU_double(x.val-y.val, x.d_val-y.d_val);
    }
    friend BU_double operator * (BU_double x, BU_double y) {
        return BU_double(x.val*y.val, x.d_val*y.val+x.val*y.d_val);
    }
    friend BU_double operator / (BU_double x, BU_double y) {
        double w = x.val/y.val;
        return BU_double(w, (x.d_val-w*y.d_val)/y.val);
    }
    friend BU_double operator + (BU_double x) {
        return BU_double(x.val, x.d_val);
    }
    friend BU_double operator - (BU_double x) {
        return BU_double(-x.val, -x.d_val);
    }
    friend bool operator < (BU_double x, BU_double y) {
        return x.val < y.val;
    }
    friend bool operator <= (BU_double x, BU_double y) {
        return x.val <= y.val;
    }
    friend bool operator > (BU_double x, BU_double y) {
        return x.val > y.val;
    }
    friend bool operator >= (BU_double x, BU_double y) {
        return x.val >= y.val;
    }
    
    // 基本関数の定義
    friend BU_double sqrt(BU_double x) {
        return BU_double(sqrt(x.val), 0.5*x.d_val/sqrt(x.val));
    }
    friend BU_double exp(BU_double x) {
        return BU_double(exp(x.val), x.d_val*exp(x.val));
    }
    friend BU_double log(BU_double x) {
        return BU_double(log(x.val), x.d_val/x.val);
    }
    friend BU_double sin(BU_double x) {
        return BU_double(sin(x.val), cos(x.val));
    }
    friend BU_double cos(BU_double x) {
        return BU_double(cos(x.val), -sin(x.val));
    }

    // coutに出力するフォーマットの定義
    friend ostream& operator<<(ostream &s, BU_double x) {
        return s << "BU_double("<< x.val << ", " << x.d_val << ") ";
    }

};

int main(void){
    cout << setprecision(15);
    BU_double x, y;

    // sin(x)の微分 
    x = M_PI / 4;
    x.select();
    y = sin(x);
    cout << 2 * y.get_d_value() << endl; // 2016-07-13に修正

    
    // y = 10*x^2 のfor文を使ったアホな定義
    x = 5;
    y = 0;
    x.select();
    for(int i=0; i<10; ++i) y = y + x*x;
    // x = 5.0のときのyとdy/dxを出力
    cout << y << endl;

    // ニュートン法の実行
    x = 2.0; // 適当な初期値
    for(int i=0; i<10000; ++i) {
        x.select();
        // y = (x-sqrt(2))*(x^3+1)
        y = (x-sqrt(2))*(x*x+x+1);
        // yの値がdoubleの丸め誤差と同程度に小さければ終了
        if (y.get_value() < 1e-15) break;
        // ニュートン法の更新式で次のxを求める
        x = x - y.get_value()/y.get_d_value();
    }
    // x(=sqrt(2))を出力
    cout << x.get_value() << endl;
}

出力は以下の通りです。

1.41421356237309
BU_double(250, 100) 
1.4142135623731

出力の1行目は $\sin{x}$ の微分の $x=\frac{\pi}{4}$ のときの値の2倍（すなわち $\sqrt{2}$ ）を表しています。最初に示した数値微分より高い精度で計算できていることが分かります。
2行目は $y=10x^2$ という関数に $x=5.0$ を入れて計算したときの関数値と導関数値です。for文を使った繰り返し文として定義しているにも関わらずきちんと導関数値が求まっていることが分かります
3行目はニュートン法を用いて求めた $(x-\sqrt{2})(x^2+x+1)$ の解を表しています。ニュートン法は $f(x)=0$ の解を $x_{n+1} = x_n - \frac{f(x)}{f'(x)}$ という更新式を用いて反復することで求めるというアルゴリズムです。導関数の計算が自動的に行われるため簡潔に記述することができます。