Microsoft Excelを用いたケモメトリクス計算(2)−反復計算−

吉村 季織, 茂谷 明宏, 高柳 正夫


Return

1 はじめに

実験的に得られた化学データを数学的・統計的手法により解析する「ケモメトリクス」が、近年盛んに用いられるようになってきた。しかし、ケモメトリクスに関する教科書類[1 - 4]がいくつか出版されている一方で、大学等の授業ではほとんど扱われていないのが実情である。そこで本シリーズは、より多くの学生や研究者がケモメトリクスを用いることができるようにするために、最も普及している計算ソフトMicrosoft Excel(以後Excelと記す)によるケモメトリクス計算について解説することを目的としている。
今回は反復計算をExcelのワークシート上で実行する方法を解説する。反復計算では、所定の計算を指定した条件が成り立つまで繰り返す。科学技術計算では、しばしば反復計算が必要となるが、ケモメトリクスでは主成分分析(principal component analysis : PCA)を実行する上で、固有値問題の求解に反復計算が必要となる。しかし、Excelのワークシート上で反復計算を行うことは難しいと考えられており、Visual Basic for Applications (VBA)などの利用[5 - 8]が一般的である。研究の場では、効率の面からも完成されたプログラムを使うことも適当であるといえる。しかし教育の場では、計算のアルゴリズムや数学的な意味合いを理解することも重要であり、Excelはその強力な手助けとなる。計算内容を理解していれば、VBAプログラムや数値計算用ソフトウェアを利用した計算結果をより深く考察することができるようになると思われる。
Excelの基本機能で反復計算を行うには、滴定曲線のシミュレーション[9]や固有値問題の求解[10]のように反復計算そのものをワークシート上で展開する方法がある。ただし、反復回数が予め固定されてしまうこと、計算量が多いと大量のセルを消費し、PCに対する負荷が大きくなるなどの欠点がある。また、循環参照(例えば、A1に"=B1"、B1に"=A1"と入力したときのように、参照関係が無限ループになってしまう場合)を利用することも考えられるが、設定変更が必要なため一般には知られていないこと、入力や削除といった作業をするたびに計算が行われるのでPCに負荷がかかることなどが問題として挙げられる。
我々は、「コピーを行った範囲が、貼り付けを行った後でも有効になっている」というExcelに特徴的な動作に着目した。本稿では、まず1回分の計算をExcelワークシート上に組み立てる。そして「コピー」、「値貼り付け」および「直前の動作を繰り返す」コマンド、この3つの動作を利用することで、これまで難しいと考えられていたExcelワークシート上での反復計算が可能であることを示す。
本稿の内容は、Microsoft Windows上のExcelを想定している。そのため、Mac OS上のExcelにより計算を行う場合には、キー操作などを適宜読み替えてほしい。また、本稿中の図では、Excelのワークシートは列幅を初期状態より狭い5.0 (45 ピクセル)にしてあるので、その点に注意して計算値(最終桁の値など)の桁数を見てほしい。

2 固有値問題とその数値解法

固有値問題の式は、行列X、固有値h、固有ベクトルv(列ベクトル)を用いて、

と書くことができる。PCA用いる場合、実行列とその転置行列の積によりXが得られるため、Xが実対称行列となること、hが正となることが知られている。Xの階数をrとすれば、r個の固有値と固有ベクトルが見つかることになる。これらの固有値のうち、絶対値が最大のものを見つける方法として累乗法(power method)[11]が知られている。累乗法では、まず適当なベクトルvi(||vi|| = 1, ||.||はノルム)を用いて次式により計算して

vi+1を求める。このvi+1viと考えて、式(2)の計算をする。この作業を反復する。理論上、この反復計算を無限回繰り返すことで、固有ベクトル(v¥)と最大の固有値(h = v¥TXv¥、ただし、今回の例では固有値が正であることが分かっているのでh = ||Xv||)を得ることができる。ただし、現実問題として無限回の反復を行えないこと、コンピュータで扱う数値は無限の精度を持っていないことから、適当な収束条件を設けて有限回の反復で中断する。本稿では後述するように、収束する様子が視覚的に観察できるため、機械的な収束の判断を行わず、目視によって判断することとした。

3 Excelによる累乗法の実行

3. 1 ワークシートの準備

式(2)で示される反復計算をExcelの基本機能のみ行う例を挙げる。Figure 1に示すワークシートを用意した。また、このワークシートの作成手順をTable 1に記した。


Figure 1. Preparation for calculation of power method on a worksheet. The ranges surrounded by dashed lines have to be filled by values following this figure. The values in the ranges surrounded by bold solid lines are calculated by formulas shown in Table 1.

Table 1. Procedures to construct worksheet for power method.
addressformula or valueinput style, etc.referred equationdefined name
1A3:E7fill these ranges with values following Figure 1.matrix XX
2H3:H8initial value of vivi
3I3:I7=MMULT(X,vi)array*2(2)Xvi
4J1=SQRT(SUMSQ(Xvi))normal*1(2)norm
5J3:J7=Xvi/normarray*2(2)
6J8=H8+1normal*1(2)
*1 normal input: when inputting a formula in a cell, hit only [Enter] key.*2 array input: when inputting an array formula in cell(s), keep pressing [Ctrl] and [Shift] key, and then hit [Enter] key.

Figure 1の中で点線に囲まれている範囲は、直接数値を入力しなくてはならない。特にA3:E7は行列Xであるため、Figure 1の通りに入力しなくてはならない。ベクトルviであるH3:H7は、反復計算の初期値(v0に当たる)を入力しておく。v0は一般には、Xの中の最も分散の大きな列が用いられるが、できるだけワークシートを簡素にするために、Figure 1に示すような値に設定した。H8は反復回数を示すiであるため、直接計算には用いられないので省略してもかまわないが、反復回数を明確にするために表示するようにした。
実線で囲まれている範囲は、Table 1に示されている数式を入力する必要がある。XviのノルムであるJ1と、i+1であるJ8は数式を打ち込んだ後、通常の入力方法([Enter]キーを打つ)でかまわないが、Xvivi+1はそれぞれI3:I7とJ3:J7に配列数式を用いて、配列として算出するため、配列入力([Shift]キーと[Ctrl]キーを押したままの状態で、[Enter]キーを打つ方法。詳しくは既報[12]を参照のこと)しなくてはならない。

3. 2 反復計算の実行

Excelの基本的な機能、"コピー"、"値貼り付け"、"直前の操作の繰り返し"のみを用いて反復計算を実行する。この手順を示したのがFigure 2である。行列Xは反復ごとに計算に用いられるが、X自体は変化しないので、Figure 2ではXを省略して示した。以下、手順に関しての解説である。
(1) vi+1i+1であるJ3:J8を選択する(a)。
(2)コピーコマンドを実行する(操作例:[Alt]→[E]→[C]の順にキーを打つ)。するとJ3:J8が点滅点線で囲まれる(b)。これは、この範囲がコピーされていることに加え、cut-copy modeになっていることを示している。
(3) viの先頭行であるH3を選択する(c)。
(4) 値貼り付けを行う(操作例:[Alt]→[E]→[S]→[V]の順にキーを打つ)。通常の貼り付けでは、コピー元のセルに数式が入力されている場合、コピー先のセルには数式が貼り付けられる。値貼り付けを行うと、コピー先には"数式"ではなく、"値"が貼り付けられる。値貼り付けを行うと、直前までJ3:J8に表示されていた値が、H3:H8に値として入力される。貼り付けられた値に従い再計算が行われる(d)。また、H8の初期値は"0"であるため、"=H8+1"が入力されているJ8は、初期状態において"1"と表示される。このJ8をコピーしてH8に値貼り付けをしたために、H8には"1"が入力され、数式に従いJ8は"2"と表示されることとなる。すなわち、i = 2の状態が出来上がる。ここで特に注目すべき点は、J3:J8を囲んでいた点滅点線がいまだに表示されていることで、これはこの範囲がいまだコピーされていることに加え、cut-copy modeになっていることを示している。
(5) ここで、続けて値貼り付けを行えば、i=3の状態を求められることになるが、何回もの反復のたびに前述した手順で値貼り付けを行うのは効率が悪い。そこで、[F4]キーを押すことで"直前の操作の繰り返し"コマンドを実行する(もし使用しているキーボートに[F4]キーがない場合は、[Ctrl]キーを押したままの状態にした状態で、さらに[Y]キーを押すことで同様のコマンドを実行できる)。この時点のJ3:J8の状態がH3:H8に値貼り付けされることとなる(e)。また、ワークシード上部の数式バーは、アクティブセルであるH3の数値が表示されている。
(6) "直前の動作の繰り返し"を繰り返す。9回の反復で、H3:H8とJ3:J8の表示値が等しくなった(f)。しかし、これは列幅の範囲で表示できる桁数で等しいだけで、数式バーを見ながら"直前の動作の繰り返し"を行うことで、セルに表示される桁数を超える範囲で値が変化していることが分かる。
(7) さらに"直前の動作の繰り返し"を繰り返す。多数回反復させる場合、[F4]キーを押したままにするとよい(もしくは[Ctrl]キーを押したままで、かつ[Y]キーも押したままにする)。数式バーの数値により、反復が進むにつれて、数値の変動がより低い桁へと移動していく様子を観察することができる。十分な回数反復を行うと、見かけ上数値が変化しなくなるか、非常に小さい範囲で数値が振動するようになるので、収束に達したと判定できるようになる。今回の例では、32回の反復で収束に達したと判断できる(g)。この方法では、ユーザーがキーを押して反復を進めるので、収束の判定もユーザーが納得できるまで押し続けてかまわない。
収束に至った状態において、J1が行列Xの最大の固有値hであり、これに対応する固有ベクトルvがJ3:J7となる。ここまでが、Excelの基本機能だけで反復計算を行う手順である。


Figure 2. An example of iterative calculation on Excel worksheet for solving eigenvalue problem shown in Figure 1. Following above procedure, an iterative calculation without any VBA macro programs can proceed.*1The blinking dashed lines surrounding J3:J8 show this range is copied and cut-copy mode is valid. *2The value in H3 is shown in the formula bar. One can also confirm visually convergence by observing the formula bar. *3In Excel, "redo" is a command to repeat an operation done just before and is run by hit [F4] key. *4When many times iteration is needed, redo can be repeated by keep pressing [F4] key.

3. 3 残りの固有値、固有ベクトルに関して

累乗法を用いて最大固有値と対応する固有ベクトルを求めた。行列Xにはこれ以外にも固有値と固有ベクトルが存在する。Xが今回の例のように実対称行列の場合、固有値を|h(1)| > |h(2)| > |h(3)| > ...、これらのh(j)に対応する固有ベクトルをv(j)として、Xを表すと、

と書ける(上付きのTは転置を表す)。累乗法を使って求まったh(1), v(1)を使って、

を求める。X(1)は式(3)からわかるように、Xからj=1の項を消去した残差行列となっており、h(2)が最大の固有値となる。そこで、X(1)に対して累乗法を適用することで、Xの2番目に大きい固有値h(2)とそれに対応する固有ベクトルv(2)が求まる。今回の例で、式(4)を用いてA10:E14に残差行列を求めるとするならば、この範囲に配列数式
    =X-J1*MMULT(vi,TRANSPOSE(vi))  (数式1)
          =X-J1*vi*TRANSPOSE(vi)  (数式1')
を配列入力する。行ベクトルと列ベクトルの積は、(数式1)が基本的な記述であるが、(数式1')のように配列の掛け算としても計算できる。このように算出されたA10:E14をコピーし、A3:E7に値貼り付けする。貼り付けられた行列の固有値、固有ベクトルはFigures 1, 2に書かれている手順で求めることができる。A10:E14に(数式1)が既に入力されているならば、この範囲にはさらなる残差行列が算出されている(今回の例ではX(2)のすべての要素は0とみなせるようになる)。これらの手順を逐次進めていくことで、すべての固有値と固有ベクトルを求めることができる。

4 終りに

ケモメトリクスで頻繁に用いられるPCAを行う際に必要となる、固有値問題の求解のための反復計算を、コピーと値貼り付けにおけるExcelに特徴的な動作を利用し、基本的な機能のみで解決した。
この方法は、VBAプログラミングにも応用が可能である。通常のプログラム言語で固有値問題を解くためには、行列を格納するための配列を準備し、行列の計算手順などを記述しなくてはならない。この部分をワークシート上にて展開することで、VBAでの記述は、値貼り付けの繰り返し、収束の判断、結果の保存、残差行列の算出といった部分を記述するだけでよく、プログラムコードを簡素化することが可能となる。
また様々な場面への適用が可考えられる。例えば、Brownian運動の簡単なシミュレーションとしてランダムウォークを考えるとした場合、従来の方法でExcelの基本機能だけを用いるなら、粒子の運動の軌跡をグラフに図示するのが一般的である。本稿で紹介した方法ならば、粒子の動きを動的にシミュレーションできるという特徴がある。ケモメトリクスに限らず様々な場面で活用してほしい。

今回例に示した行列Xは、下に示す行列Mより生成される。

Mの階数は2であるため、Xの階数も2となり、固有値と固有ベクトルは、それぞれ2つずつ得ることができる。固有値h(1), h(2)はそれぞれ729, 255であり、これらに対応する固有ベクトルv(1), v(2)

である。また、Xの固有値を求めることが、Mの特異値分解と同義であることに着目して、Y @ MTMとし、Yの固有ベクトルをu(i)(1行3列の行ベクトル。||u(i)|| = 1)とするならば、

という関係があることを覚えておくと、固有値を求める際の行列が小さくでき、計算量を減らすことができるので有用である。

参考文献

[ 1] 長谷川健, スペクトル定量分析, 講談社サイエンティフィク (2005).
[ 2] 三井利幸, ケモメトリックスの基礎と応用−分析化学と多変量解析法−, アイピーシー (2003).
[ 3] 尾崎幸洋, 宇田明史, 赤井俊雄, 化学者のための多変量解析−ケモメトリクス入門−, 講談社サイエンティフィク (2002).
[ 4] 宮下芳勝, 佐々木慎一, ケモメトリックス 化学パターン認識と多変量解析, 共立出版 (1995).
[ 5] 内田治, すぐわかるEXCELによる多変量解析, 東京図書 (2000).
[ 6] 吉村忠与志, 青山義弘, トランジスタ技術Special No.78 技術者のためのExcel活用研究, CQ出版社 (2002).
[ 7] 小椋將弘, Excelで簡単 多変量解析, 講談社サイエンティフィク (2006).
[ 8] 菅民郎, Excelで学ぶ多変量解析入門, オーム社 (2007).
[ 9] 吉村季織, J. Comput. Chem. Jpn., 5, 29 (2006).
[10] 室淳子, 石村 貞夫, Excelでやさしく学ぶ多変量解析, 東京図書 (2004).
[11] 高倉葉子, 数値計算の基礎―解法と誤差―, コロナ社 (2007).
[12] 吉村季織, 高柳正夫, J. Comput. Chem. Jpn., 7, 71 (2008).


Return