見出し画像

excelで試す、重回帰分析と主成分分析。(組み合わせて属性分析を担うことができるのか？)

アールル

今秋、機会あって、比較的大きめなアンケートデータを頂くことができました(30000円：高額注意)。
大きな買い物でしたが、その分本格的に分析ツールを動かして、データを整形するのにかなり適しており、分析ツールを動かすと、いろいろ面白そうな結果や分析方法を考えることができました。
データの提供者様ご本人の分析記事はこちらになります。
そこで、遊びでいろいろ解析した結果を、とりあえずここに残すことにします。SPSSは持っていないため、全部VBAを動かしつつ、エクセルの機能を使って計算しております。正確性についてはご容赦を。
また、正統なアカデミックな分析の手法については一切ここには載っていませんので、転用の際には十分注意して使うようにしてください。
エクセルの計算力に驚かされた何日間かでした。本題に入る前に、前説明として、統計分析の手法の簡単な説明を2種類書きます。ご存じでしたら読み飛ばしてください。
文字数がすごいことになってる。

追記:えっと、、一行かけて言うのはあれですが、Excelを含めて、コンピュータ計算による逆行列の計算は、規格(IEEE 754、だそうです。)の問題からかなり行数の少ない段階から、誤差の生まれる計算方法になっています。Excelの逆行列は、計算しやすいですが、正確性と厳密性には非常に難のある計算方法でありますので、よく注意されてください。すみません。
逆行列の計算部分については、マクロで掃き出し法の自動化を行うことで多少精度の改善はできますが、それでも十分な精度を保てないことを確認しています。掃き出し法で逆行列を作ることを介さず、行列の等式の両辺に掃き出し法に相当する同じ行列を掛けてみたところ、計算結果にある程度の精度を保てている可能性を感じましたので、この方法をとりあえず推奨します。最後のLU分解での固有ベクトルの求め方も、同じ方法で大丈夫だと思いますが、計算がやや複雑になりそうなので、別に記事を作ろうと考えております。

重回帰分析

一つの目的変数に対して、複数の説明変数で目的変数を予測できるようなモデルを考えることが、重回帰分析の主な目的です。
ある程度の数のサンプルを観測し、その結果を複数の要素を表す変数の値で置き換えたデータがあるとする。
例えば、ある一つの変数を以下のような他の変数の一次式の形で表すことで、ある変数の値を予測する関係式が作れないかを考える時に使われるそうです。
データから目的の指標の変化を予測したり、他の指標から予測したいデータの分類を判別する式を推定するのに使われるようです。

まずここのページにある解説を借りて、最初に基本的な回帰方程式と、最小二乗法による偏回帰係数の求め方を簡単に示します。

重回帰方程式

b:偏回帰係数　x:説明変数(nは説明変数の番号)が、実測値

に対する残差、

ができる限り0に近くなるように重回帰方程式を作ります。
後で確認しておかないとめんどくさくなるため、先に注意書きを書いておきますが、上の式で言うxの部分はデータに応じて値を変化させる変数です。
全てのxに掛け算でかかっているbのアルファベットは、データによって変動をすることのない、変数xの係数ということなります。

最小二乗法による係数の推定方法

残差平方和

この値が最小になるように、偏回帰係数を推定します。それぞれの説明変数の偏回帰係数bに対して、δS/δb=0になるようにbを決めるので、

となります。

を代入して、整理する。

となります。xiとxjの平均偏差の積和を、

xiとyの平均偏差の積和を、

とすると、

…

に変形できます。（正規方程式）
bをn行1列のベクトルb(b1,b2,…,bn)、Sijをn×n次元の(i,j)成分をSijに持つ正方行列S、Siyをn行1列のベクトルyと表現すれば、

Sb=y

と表現することができます。Sについての逆行列を両辺に掛けて、

と変形したベクトルbが、最小二乗法での偏回帰係数の解になります。(切片は平均についての式から求められる。)

最小二乗法による重回帰分析を扱うに当たって、今回主眼を置きたいこと

先のページに依れば、重回帰モデルの仮定として、
・重回帰式は線形の一次式であること
・n個の誤差同士には系列相関がない
・n個の誤差の平均値は0で、分散は等しく、その分布は正規分布に従う
・説明変数と誤差は互いに独立である
とされるそうです。この仮定に従う限りは、説明変数と目的変数のそれぞれの積和を使い、最小二乗法にて推定される、偏回帰係数を求めることができると考えられるようです。
本来は母集団の分布の推定を含めて、厳密に推定するもので、回帰関数の当てはまりを良くするために変数の数を絞る必要があるみたいですが、この変数の値を少し強めたら、最終的に調べたい変数がどれだけ上がるか、関数的に捉えたいことに主眼を置きたいと思います。
最小二乗法で求める回帰係数は偏微分(他の変数を固定して、一つの変数だけを仮に動かしたら、元の関数はどう変化するか？)で求まる係数であるため、偏微分の性質になるべくイメージを近づけて理解をしたいという大まかな見立てを立てて、この記事を書いています。

主成分分析

相関のある多数の変数から、相関のない少数で全体のばらつきをよく表す主成分と呼ばれる変数を合成する多変量解析の一手法。とされる(Wikipediaより)。データの次元を削減するために主に行われるとされます。

観測データから共分散行列もしくは相関行列を取り出し、この行列の固有ベクトルと固有値を計算します。

この共分散(or相関)行列から計算した固有ベクトルと固有値は、それぞれ、データの分散を最大化した上で、互いに直交するそれぞれの大きさが1のベクトル(固有ベクトル)と、それぞれのベクトルで説明できる分散(固有値)の大きさを表しているということになります。

つまり、それぞれの固有ベクトルは、全て互いに直交しており、データの相関構造を分散の強さとして最もよく表すことのできる組み合わせを取っているということになります。
たぶん、この理解であってるはず。。
また、それぞれの固有ベクトルに対応する固有値は、全てのデータをそのベクトルの向く方向を軸にした時に、そのベクトル方向でのデータの分散の強さを数値として表しているそうです。

普通の主成分分析では、一度共分散行列を作ったら、その固有ベクトルのうち、元のデータの分散を最もよく説明できる固有ベクトルを第一主成分とします。
また、第一主成分と直交するベクトルのうち、分散の最も大きいもの、つまり、固有ベクトルのうち、固有値の大きい二つのベクトルを選んで、因子分析に繋げることになるようです。
これは、共分散(相関)行列の固有値と固有ベクトルから、固有値の大きい固有ベクトル2つを選んで、それらによりデータの主な相関を説明することを目指すもののはずです。

主成分分析は、データの内部構造を分析するのに主眼置いた上に、それ単体で何かを分析されることは少なく、データ分析の下ごしらえ的な分析手法のようです。少なくとも、心理学統計の分野においては。
しかし、時節柄、Twitterではデータの構造そのものよりも、その構造から何を言えるかについて話題になっていたので(暴力性とモテ議論、が11月ごろしばらくTwitterの一部界隈で盛り上がっていました。)、気ままに話題から発想を得て、データの構造を下ごしらえできる主成分分析の手法を崩して、いろいろと遊んでみました。

データ構造を表す固有ベクトルを全て回帰分析の従属変数に放り込んでみたらどうなる？

SPSSというアプリケーションの機能の強さはわかりませんが、SPSS、エクセル、共に単体で分析できる範囲は広く、組み合わせてもいろんな分析ができたような記憶があります。(ちなみに、随分前にSPSSは使った記憶があります。)

SPSSの性能の限界の話をすると、重回帰分析は重回帰分析として、かなり多めの変数について従属変数にして分析ができたかな？と思います。ただ、回帰分析を扱う上での統計学的な問題として、変数が多めだと重回帰関数としての説明力に問題があるそうなのです。(昔、データの変数の削減がよくわからず苦心したことがあったように思います。)主成分分析は、あまり使ったことはなかったはずですが、主成分2つから3つくらいに変数を絞るために使う技法だと、今のところは捉えています。

それはともかく。。

さて、エクセルでは、t検定などの統計的有意性を検証するためのデータは出してくれないのですが(エクセル統計という追加機能パッケージは、別にあるそうです。)、基本的な数値セットの整形には手軽に使えるツールが多く、手作りで基本分析のツールを作ることはできます。

手作りする際に、敢えて方法論を踏み外して、オリジナルのデータを作ってみましょうと考えて、データをいろいろ整形してみていました。元データは、モテに関するインターネット調査(すももさん作成)で、サンプル数は1000(n＝1000)、質問項目数は100項目程度です。モテを表す、交際経験があるかどうか(1：交際経験あり。0：交際経験なし。)のダミー変数を目的変数とした、全変数の回帰分析を試しにやってみました。

超多変数(200変数以上)による回帰分析です。

①データを整形します。

量的変数(年収項目：1→年収なし。2→0〜200万未満。3→200万以上〜400万未満。4→400万以上)のデータは、そのまま。もしくは、数値を関連が期待できるスケールに変換します。

質的変数は、選択なし、以外の全ての選択肢をダミー変数に置き換えます。
例えば、血液型AorBorOorABなら、1：A、0：Aではないの変数を一つ。1：B、0：Bではないの変数を一つ、という具合に4つの選択肢の全てについてダミー変数を用意します。
選択肢なしについてダミー変数を用意するかどうかについては、後で目指したい解釈によって決めることにします。
一部の質的変数も、変数の質的変化に従った、直線的な相関を示さないことが考えられますので、一部の質的変数もダミー変数に置き換えます。

今回は結果として、交際経験を説明したい変数を、説明変数とし、全部で212個の説明変数を設定することになりました。目的変数である交際経験のダミー変数を1個として、今回は計213個の変数を用いて、重回帰関数を求めます。

②全ての変数につき、平均値を求めます。2種類の変数同士で全てのサンプルについて、平均値との差の積を求め、足し合わせます。全ての変数の組み合わせについて、これを行います。

変数の全てについて、平均値を求め、平均値との差を求めるシートを作成します。変数同士の積和は、SUMPRODUCT関数にて、達成することができます。2種類の変数同士でそれぞれ、全てのサンプルについて、平均値との差の積を求め、足し合わせます。
足し合わせた数値を2つ変数の積和といいます。
行番号と列番号が一致する場所については、同じ行もしくは列の変数から、平均の差の2乗和を全て求めて入力します。
正方形の行列を、正方行列と言います。行番号と列番号をひっくり返しても同じ数字が並ぶようにして、正方行列を完成させます。
マクロを作るか、頑張って数式をコピーしましょう。
説明変数がm個あれば、説明変数の個数の2乗個の数値でできた、m行×m列の正方形の行列ができるはずです。
m行×m列の正方行列です。今回は212個の変数を作ったので、212行×212列です。

また、目的変数である、交際経験の変数と他の変数との積和も全て求めて、別にセルを作って記入しておきましょう。

交際経験の変数とその他の変数との積和がm個。また、説明変数同士の積和がmの2乗個。計mの２乗＋m個の数値が、新たに用意されるはずです。
作った行列は、サンプル数(n＝1000)で割れば、共分散行列になります。

③逆行列を求めます。

正方行列から、逆行列を作ります。(この時、交際経験の変数を含めて逆行列を作ると、最後に回帰係数が求まりません。エラーになります。)

200以上の行と列がある正方行列でも、Excelなら一瞬で逆行列を求めることができます。できるます！(とりあえず強調。)
行列を表示したい範囲を選択した状態で、MINVERSE関数を入力します。そして、選択範囲を動かさないまま、ctrl+shift+enterを押します。
すると、選択した範囲が表示範囲となり逆行列が表示されます。

本物の回帰分析では、変数同士の内生性をできるだけ低めるために、VIF(分散拡散係数)を求めることがあるそうですが、この逆行列の対角要素からも、擬似的に求めることが可能なはずです。
VIFとは、回帰係数を求めるのに使う、説明変数同士の相関の強さを表すスコアだそうです。
このスコアが高ければ、従属変数同士で共通するファクターが存在する可能性が伺われます。
同じファクターが重複して計上されている可能性があるので、回帰関数が目的の変数から大きく離れた数値になる可能性が高くなるようです。
(とりあえず僕の理解ではそんな感じです。統計については素人なので、間違ってたらごめんなさい。) 
VIFは相関係数行列の逆行列の対角要素、だそうです。
相関係数行列は、それぞれの変数を標準偏差で割った共分散行列です。
また、逆行列にする前の行列は、共分散行列の逆行列の各数値を1000分の1倍した行列に、数値的には等しいです。
逆行列は、実数に対する逆数のような概念であるため、元の行列の要素が実数倍されれば、逆行列は同じ実数の逆数をそれぞれの要素に掛けた行列となります。
従って、逆行列の対角要素から、1000を掛けて変数の標準偏差の2乗で割ればVIFと同じ数値が取り出せることになります。
(共分散行列の対角要素は、同じ変数の標準偏差の2乗で掛ければ、相関行列の逆行列の対角要素になります。)
VIFは10以上あれば多重共線性ありと判定されるそうなので、そういう見方で見るとVIFについてはなんとなく怪しい数値が並んでいるように見えます。

今回は予測精度の高い回帰関数を作るのが目的ではありませんので、VIFは無視します。

④求めた逆行列と、交際経験とその他の変数との間の積和を、行列として掛け合わせて、回帰係数のリストを計算します。

交際経験のダミー変数と他の変数の間の積和を、説明変数同士の積和の正方行列(m行×m列)を作った時の、説明変数の順番と同じ順番に、m行×1列の行列に縦に並べます。
今回は、説明変数が212個ありましたので、212行×1列の行列になるはずです。

③で作った逆行列を、左から掛けます。つまり、③の逆行列に、交際経験の関わる積和の行列を右から掛けることと同じです。
行列計算は右から掛けるか、左からかけるか、どちらの順番で計算するかによって、答えが変わっていきますので、よく注意して計算してください。

Excelでは、MMULT関数が、行列の積を導く関数に該当します。
Excelのセルを範囲選択をして、MMULT関数を入力します。
MMULT関数を行列として表示するために、縦にm列横に1列選択します。
第一変数欄に逆行列の位置を入力し、第二変数に交際経験の変数の関わる積和の行列について、位置を入力します。
m行×m列の正方行列にm行×1列の行列を右から掛けた場合、計算後の行列は、m行×1列の、縦並びの行列になります。

この結果できた行列は、全ての変数を説明変数に入れた、交際経験を予測する重回帰分析の結果に一致します。
積和行列を作った時の変数の並びと同じ順に、上から係数が当てはまっていきます。
切片は、全ての変数の平均値を重回帰関数に代入することで、求めることができます。

この計算は、最初に説明した重回帰分析の説明の最後にある、b＝S^ｰ1yが今求めた計算式と、同じ式の形をしていることから、説明ができます。

今まで説明した方法論については、以下のサイトが詳しいです。

⑤それぞれの変数の重回帰係数同士を見比べます。

計算した重回帰関数自体は、背景因子が重複して計算されている可能性が考えられるため、あまり重要ではありません。

まあ、ここのツリーで扱っている評価のためのスコアは、目的変数と説明変数のそれぞれの標準偏差で標準化された、標準化偏回帰係数の話ですが、標準化しなくても、ここの変数が1大くなると、交際経験スコアがこれだけ上がりやすくなるよという推定を、他の変数との比較で進めることができるかもしれないと思います。

また、全変数の積和を求めると、あくまで目的変数との関係の強さを比較する文脈の上では、変数同士を自由に選びながら、重回帰分析ができるようになると考えられます。
互いに積和のわかる変数同士であれば、今回重回帰分析を実行した同じ手順で、多数の変数を自由に選んで、変数同士の相関を自由に推定したり、背景因子の有無を考察したりする余地が生まれます。
自由な発想の考察や研究の糧になる…かもしれないです。

全変数回帰分析の結果を一部紹介。

それぞれの変数と、回帰係数の間に、弱い相関が想定されるような回帰係数の内容になっています。
スケールの違いや、標準偏差の違いによって、直接変数同士の関連を比較することはできませんが、大雑把な関係の強い弱いの可能性を見て取ることができるように思います。
女性の好みについて、都会的な人が好みな人は、少し交際経験からは遠い可能性がある。
落ち着いた人が好みな人や、結婚相手に相手の学歴を重視する人は、少し交際経験から近い可能性がある。みたいな感じです。

ほんとは多項ロジスティック回帰分析の方が適しているらしいのですが、今回はExcelで確実なロジスティック回帰分析を行う方法が見つかりませんでした。
なので、単純なる重回帰分析の結果で勘弁を。。
これだけの変数を使って、多項ロジスティック回帰分析って、組めるのだろうか？？

ーーメインテーマ1、とりあえず終わり。。。長かった。。ーーー

主成分分析で用いた、固有ベクトルと固有値群を、捨てずに全て回帰分析の変数の一部として使ったらどうなるか？

主成分分析。基本的に、固有値と固有ベクトルの群が互いに直交するベクトル群で、尚且つ共分散行列の元となったデータの分散を最もよく説明できるパラメータになる方向の組み合わせであることを利用した、元データの構造分析に該当する分析です。
が、これを次元削減のために使わず、そのまま回帰分析に掛けて、軸同士の目的変数の相関を比較するために使ったらどうなるか。これが二つ目のお遊びです。

一つ目のテーマで用いた積和行列をそのまま転用して、固有値固有ベクトルから、新しい変数の積和行列まで、全ての計算をExcel計算で算出することができます。
固有ベクトルの性質として、互いのベクトル同士は直交ベクトルであり、尚且つ、新しいベクトル同士に相関関係はないです。
新しい情報が付加されるような計算というよりは、積和行列をもう少し細かく構造化したい時に便利そうといった印象の計算です。積和行列、共分散行列、もしくは、相関行列を、もう少し見やすく表すことができるのではないかなとは思っています。とりあえず、やってみましょう。

①固有値と固有ベクトルを全て求める。

Excelには、固有値と固有ベクトルを求める関数は定義されていないのですが、全ての固有値と固有ベクトルを求める計算式を組むことはできます。VBAのマクロを利用しますが、そんなに難しくはありません。
ハウスホルダー行列を用いたQR法による求め方を記事の最後に乗せます。一旦は本筋とは関係ないテクニックになりますので、割愛します。

②固有ベクトルの軸から求められる新しいスコアのそれぞれついて、積和に相当する数値を求める。

例えば、固有ベクトルをz{c1,c2,…,cm}とし、zの軸方向の大きさを表すスコアをrと定めると、rのスコアはr＝c1X1＋c2X2＋…＋cmXmと表すことができます。
平均値との差に限ってはzが元の変数Xとその係数から求められるため、元の全ての変数の積和からzの積和を求めることができます。
ここからの計算で、切片は平均値との引き算をするときに全て相殺されて、最終的な式には残らないように計算できます。従って切片は全て省略して記載します。

とすると、

のようになります。
新しい変数は、全て、元の変数に固有ベクトルの各要素を係数として掛け算した式で表すことができます。
これはすなわち、元の変数の積和で新しい変数の積和を全て表現できることになります。

目的変数である、交際経験の変数との積和も、元の変数との積和のリストが全て用意されていれば、同じように求めることができます。

Excelで新しい変数の積和を求める場合、元の変数の積和の行列の上と左に、積和を求めたい二つの固有ベクトルの数値を横と縦にそれぞれ並べます。
上に並べた数値は下へ下へと同じ数値を掛けて、左に並べた数値は右へ右へと同じ数値を掛けて、新しい行列を作ります。
この行列の要素を全てそのまま足し合わせれば、新しい変数の積和を求めることができます。
200も変数があると、全ての変数の組み合わせについて積和を求めるのはたいへん手間がかかり、まず手作業では終わりません。
そのため、ループ処理をかけて、VBAのマクロで求めました。
全ての計算が終わるまでに2時間程度かかりました。

③新しいベクトルで重回帰係数を求める。

積和が全て求まれば、後は、1つ目のメインテーマで行った方法と、同じ方法を使い、全てのベクトルを使って、回帰係数を求めることができます。
固有ベクトル同士の回帰係数も、自由自在に求めることはできます。
しかし、積和行列や、共分散行列から求めた固有ベクトルは、互いに直交する性質があります。
そのため、回帰係数は0になるか、エラーになると思います。

固有ベクトルを全て使って求めた回帰係数を一部紹介

ほとんどのベクトルが、一つの変数のみ、成分として強いものを持っているように見える図表です。
ただ、一部、固有ベクトルの要素としての絶対値が大きい変数同士が同居しているベクトルがあるように見えます。
同じ固有ベクトルに強い成分を持つ変数同士で、やや強い相関があるのかもしれません。
また、主成分分析で扱うベクトルと、回帰分析の係数は本来別物であるので、固有ベクトルの固有値(=説明できる分散の量)と、それぞれの固有ベクトルの回帰係数との間には、相関はありません。

中には、一つの固有ベクトルとの強い相関をあまり持たない変数も、あるのかもしれません。

応用してみるために考えてみる

固有ベクトルに分解して見てみたところ、重回帰係数は案外大きいものと小さいものの差が大きかったです。
固有値が表す分散と合わせて、ある程度の集団の傾向を説明したり、交際経験との相関をベクトル毎にみて取ることができるように見えます。
固有ベクトル分解で得られる互いに直交するベクトルの群は、見方を変えれば、個々のデータを多次元の直交座標で表すための、ベクトルの方向を与えてくれるとみることもできます。
固有ベクトルによる新しいスコアへの変換を通して、新しい回帰関数に当てはめた場合、普通の方法とは違った視点から、データを眺めることができるかもしれないと思います。

(固有値と固有ベクトルへの分解は、元のデータに最もよくフィットするモデルのダイヤモンドを、探してくれています。)

また、それぞれの固有ベクトルの数値に注意して見てみた場合、案外似たような変数が、同じ固有ベクトルの中で高い数値を取っていることに気づくかもしれません。
積和による行列の固有値分解は、元の集団の構造(分散)を最もよく説明できる直交ベクトルに分解できるとされております。
この数値の傾向は、元の集団構造の相関の違いを疑わせる情報の一つかもしれないと考えています。

固有ベクトルから導く重回帰関数を、属性の強さや、ある属性の成分の強さとして解釈した場合、相関関係の強い変数同士の相関を説明するために有用な情報を与えているようにも見えます。

主成分分析と次元削減

主成分分析の応用機能として、次元削減というテクニックがあるそうです。

固有ベクトル分解を行った後、次元削減という手法を使うと、固有値の大きい第一主成分や第二主成分以外の成分を削ぎ落して、情報を圧縮して表すことができるというものだそうです。
次元削減を軽く数学的に説明すると、固有ベクトルのうち、主な成分2つを使ってｍ行×2列の行列を作り、主成分分析の元として使った共分散行列(もしくは相関行列)に右から掛けると、m個の変数からなるデータを、2つの新しい変数からなる、新しいスコアを作ることができます。
これは、数学的には、射像というものだそうです。

ただ、この射像。主成分分析では、固有値の大きな、主要なベクトル以外のベクトルの情報を全部削ぎ落すそうですが、逆に、主要なベクトルの情報だけを削ぎ落す手法も逆に考えれば作れるような気がします。
主成分分析から、因子分析に持ち込む際に、互いに直交しているような、第一主成分や第二主成分を回転させるだけで本当にいいのか？
主成分を削ぎ落した残りの情報から新たに固有値の大きい固有ベクトルを見出して、そのベクトルについて因子回転に持ち込んだ方がいい結果が出るのか。。
もうちょっと検討を加えてみたいと思っております。
持ち合わせてる知識と技術的には、比較は可能のように思いますので、もう少し調べてみたいと思っています。

ーーーメインテーマ2、終わり。。。ーーー

おまけ：QR法による固有値分解(ハウスホルダー変換を用いた分解の方法、Excel)

Excelで固有値と固有ベクトルを求める方法。最初ソルバーを使っていたのですが、固有値の数が多すぎて求めきれず、何か方法はないかと探していたら、QR法による固有値問題の解決方法が書かれていることを見つけました。

細かい理論の説明はWikipediaに譲りますが、方法をよくよく検討したところ、Excelと、VBAのマクロを用いて実現が可能でしたので、方法を紹介します。

用意するもの

※Excelワークシート6枚、m行×ｍ列の正方行列(固有値と固有ベクトルを求めたいベクトル)

①シート1：xベクトル、ｙベクトル、uベクトル、u×u行列計算用、xベクトルとuベクトルのユークリッド距離計算用
(縦にxベクトル、yベクトル、uベクトルの数値を計算します。
xベクトルとuベクトルについてその隣にMMULT関数で、uベクトルと、TRANSPOSE関数で、縦横を入れ替えたuベクトルとの掛け算の行列を記入します。
また、u×u行列はuベクトルのユークリッド距離の2乗で割っておいてください。xベクトル、yベクトルの初期値については、すべて0とします。)

②シート2：E-2×u×u(=Q：ハウスホルダー行列)行列計算用

③シート3：A'(ダッシュ)行列計算用(シート2のE-2×u×u行列と、シート4のA’行列との行列の掛け算を記入します。)

④シート4：上三角行列完成用、A行列記入用
(初期値として、固有値と固有ベクトルを求めたい行列の値をコピペしておきます。
計算の結果、元の値は全て消失しますので、よく注意してください。
また、行列の計算式を記入した場合、マクロで値の変更ができないので、エラーが出ます。
マクロによって、シート3の行列の値をコピペします。)

⑤シート5：ハウスホルダー行列×A’行列掛け算用(シート2の行列とシート6の行列を掛け算します。)→シート6の行列に、シート2の行列の転置行列を、右から掛けます。

⑥シート6：Q^-1行列完成用(シート5の行列の値をマクロでコピペします。
初期値はif(row()=column,1,0)などの関数を使い、単位行列Eを記入しておいてください。
完成後の行列の逆行列を計算すると、固有値に対応した固有ベクトルになります。)

手順：

上記の通りに、ワークシートを6枚用意してください。
シート1のxベクトルとyベクトルの欄は無記入。もしくはすべてに0を記入します。
uベクトルの欄はxベクトル-yベクトルとして、数式を記入してください。
xベクトルとuベクトルの各値の平方和の平方根の数式を、SUMSQ関数とSQRT関数で用意しておいてください。
uベクトルの平方和については、0にすると途中の計算でエラーが出てしまいますので、if関数で、0の時は0.000000001が出るようにしておいてください。

マクロの記録を立ち上げ、適当にシートの値をコピペして、その動作をマクロに記録してください。
後のことを考えると、のちにマクロでコピペする領域をなぞってコピペした方がいいです。
なるべく貼り付けのオプションで、数値のみの貼り付けを選んで貼り付けてください。

記録したマクロを立ち上げ、編集してください。まず初めに、マクロの最初に、

Application.Calculation = xlCalculationManual

マクロの最後に、

Application.Calculation = xlCalculationAutomatic

を記入してください。マクロ計算中のシートの自動計算を全て止めるためです。

マクロは、シート4の1列目から順に、m列まで1列ずつxベクトルの数値欄に数値をコピペすることで始まり、ｍ回までの反復計算で作業を行います。

マクロの流れは、今が第k回目の反復計算として、次の通りです。

xベクトル、yベクトルの値を0で初期化します。
→シート4のk列目の、k行目から下の数値全部を、xベクトルのk行目から下にコピペします。
→xベクトルの各値の平方和の平方根の値のみを再計算します。(Range("B1").Calculate)→
xベクトルの各値の平方和の平方根の値をｙベクトルのk行目の値にコピペします。
→Excel全体で、1度だけ再計算を実行します。(Application.Calculate)
→シート3の行列の値を、シート4に数値のみコピーします。(シート4のk列目のk行目より一つ下から下全てのセルは、割り切れないために、０に近い小さな数値が入力されますが、お好みで消去して0にして大丈夫です。)
→シート5の行列の値を、シート6に数値のみコピーします。
→最初に戻る

この通りにやってみて、シート4に上三角行列が全て完成したら成功です。
シート4の対角成分(行数=列数)は固有値を表し、シート6の行列の逆行列を求めたら、縦に1列ずつ取り出せば、全て固有値に対応した固有ベクトルの数値に対応します。

行列の固有値と固有ベクトルを全て求めることができました。

このマクロは、一応検証してみたところ、正方行列ではない行列(行数>列数)の特異値を求めるのにも使えるようで、正方行列でない行列を、右下に端を合わせて記入しても、きれいな上三角行列が完成します。
これは、次元削減をした後の行列の特異値分解にも転用できる可能性があります。
もう少し、計算方法を検討してみようと思っております。

追記、

QR分解だけでは、QR法の全ての説明を尽くした訳ではないことが判明しましたので、追記してQR法を完成させます。上のQR分解が全ての固有値と固有ベクトルを尽くしていると思ったけど、冷静に考えたら対角化できてなかったから、同じ値という訳にはいかないやん、、ということです。ちなみに、前回の記事のQR分解についても、間違いを一つ発見しましたので、先に追記します。

シート5：

QR分解後の相似変換。

方法論はこちらが詳しいです。こちらのページに従って、QRのRとQを入れ替えて相似変換にて固有ベクトルを求めつつ、逆反復法をなぞって固有ベクトルと固有値を計算していこうと思います。求めていった過程での体感では、QRを入れ替えてもほとんど数値が変わらない対角行列を作ることを目指しつつ、固有値が求まったら、元の行列と固有値を元に固有ベクトルを復元していく感じです。

行列RQを使い、固有値をひたすら正しい値に近似させていく。

用意するもの

※Excelワークシート枚、QR分解で求めた行列Qと行列R(それぞれm行×m列の正方行列)、上で作ったQR分解用のワークシートとVBAマクロ。

①QR分解用のワークシート、マクロ。この過程ではQR分解を何度も繰り返します。早く正確に計算を実行できるように、点検しておくといいかもしれません。このワークシートのうち、シート4で計算したR行列と、シート6で計算したQ行列を掛けてみて(MMULT(シート6,シート4))、元行列とQRが一致することを必ず確認しておいてください。
完全一致させるためには、上三角行列を作る工程での対角成分から下の値は、消去しない方が良さそうです。
シート5の計算は、シート6の行列に、シート2の転置行列(MMULT(シート6,TRANSPOSE(シート2)))を掛け算している必要があります。
また、シート4の対角成分をリストにして、固有値に近似しているかどうかの確認に使いますので、対角成分を1行にまとめて表示するためのマスを用意しておいてください。シート4の下の行にINDEX()関数とROW()関数COLUMN()関数を使って、マスの番地を元に参照していく方法がやりやすいと思います。

②シート7：RQ行列計算用のワークシート
シート4とシート6を掛け算した結果RQを書き出しておきます。(MMULT(シート4,シート6))
また、この時、相似変換をする際に、RQのままで計算を続けてしまうと、行列の収束が悪くなりますので、これを防ぐために、μE(Eは単位行列)を相似変換する前に引き算して、相似変換を終わらせてから足し直す計算を途中で行いながら進めます。

ここで指摘されている通り、そのまま相似変換を続けていくと収束が遅くなるようです。計算の高速化のために、上記の処理をしますが、μの値を保存するためのマスをRQ行列のマスとは別に、2つ用意しておいてください。

③シート8：固有値の変化幅チェック用ワークシート
シート4の対角成分を、値を指定してコピペするようにしてください。また、反復試行を行いますので、前の試行での対角成分と、次の試行での対角成分との差分を計算するセルを作っておいてください。(SUMSQ(対角成分(前の試行)-対角成分(次の試行)))となるように計算するといいと思います。このまま入力すると、引き算の計算ができずエラーになりますので、Control+Shift+Enterで確定して、行列形式の計算をさせるように指定すると、上手く計算ができます。

手順：

上記の通り、QR分解で用いたワークシート6枚と、あと追加で2枚新しいワークシートを用意してください。
シート7には、空マス2つを用意して、それぞれをμ1、μ2と名付けた上で、RQ–μ1E+μ2Eの計算をする行列を用意しておきます。
{MMULT(シート4,シート6)-IF(ROW()=COLUMN(),μ1–μ2,0)}あたりの計算式で、行列形式の計算を指定すると、上手くいくと思います。

求める行列Aに対して、まず1度目のQR分解を行います。
シート8では、最初にQR分解した後の行列のうち、シート4の行列の対角成分を1行にまとめて、書き出しておきます。値を指定してコピーする必要があります。

(※)まず、シート7のRQ行列のうち、1番右下の2×2行列を取り出して、別に固有値を求めておいてください。簡単に計算できると思いますが、どこかのブラウザの行列の固有値を計算するサイトを利用すれば、手間はかかりません。計算して出た2つの固有値のうち、どちらか1つをμ1のマスに描き込んでください。

シート7の行列を、今の段階でこのまま、QR分解のシート4にコピーして書き込みます。値を指定して書き込む必要があります。シート6に単位行列を書き込んでセットした後シート7のμ1の値をμ2のマスにコピーしておいてください。その後、μ1の値は0にしておく必要があります。

セットした値で、QR分解を実行します。

QR分解が終わったら、シート8に、シート4の行列の対角成分を、新たに全てコピーして、記入しておきます。これも、値を指定してコピーする必要があります。前回の対角成分との差分も取っておきます。SUMSQ関数で上手くいくと思います。

シート7では、新しいRQ(+μ2E)行列ができていますので、(※)まで戻って、対角成分の差分が十分に小さくなるまで、上の作業を反復してください。10〜20回の反復で、ある程度小さくなるので、手作業でも十分に収束が期待できると思います。

こうして、ある程度計算を反復すると、R行列とQ行列にあまり変化がなくなるところに収束していくと思います。この時のR行列の対角成分が、元行列の固有値の、限りなく近い近似値となります。(厳密に同じ値ではないようですので、注意してください。)

固有ベクトル

固有値を求めたら、次は固有ベクトルを求めることになると思いますが、上の方の印用ページで、固有値から固有ベクトルを求める方法について解説されているので、これに倣って、次は逆反復法を使って固有ベクトルを求めてみます。

逆反復法については、説明は上で紹介したリンク先に書いてあります。
QR法で求めた、全ての固有値の近似値(λ1,λ2,λ3,…,λn)について、元の行列Aと固有値の近似値λkで、AｰλkEの逆行列を求めて、任意の0ベクトルでないベクトルbに左から掛け続けます。全てのベクトルbは、Aの固有ベクトルに有理数の係数を掛けたベクトルの和で表すことができます。A-λkEの逆行列をbに左から掛け続けると、bをAの固有ベクトルの和に分解した時の固有ベクトルの係数のうち、λkに遠い固有値に対応する固有ベクトルの係数が、掛ける毎にどんどん小さくなります。最終的に、λkから遠い固有ベクトルは、掛け算を繰り返す毎に係数がほぼ0に近くなって消えていき、λkに1番近い固有ベクトルの成分のみ残ることになります。
逆反復法は、この性質を利用して、それぞれの固有値に対応した固有ベクトルを計算していく方法です。

計算理論だけは、これだけ説明すれば、後は地道に計算するだけなのですが、この場合、逆行列を求める方法が、いきなりexcel標準関数のINVERSE()に飛びつくと、誤差だらけで元の固有ベクトルの推定もままならないほどに実際の値から乖離した計算になりますので、最後に逆行列を含む計算の方法の紹介をします。

逆行列計算の誤差

ちなみに、メジャーな逆行列の計算方法としては、掃き出し法、余因子行列を使った方法、LU分解(エルユー分解)というものを使った方法と3パターンあります。LU分解が、高速で、誤差も少ないとされています。

この方法なのですが、まだ調べてはいないのですが、計算量がかなり抑えられて、誤差も少ないとはされるのですが、余因子行列を使う方法と同じく、式の途中で割り算と、分数の和をかなり多用するようなのです。コンピューター計算をする際には、それぞれの数値はメモリの制約の関係上、決められた桁数より細かい数値は、丸められて、桁から下の数値情報が落とされることになっています。

逆行列の計算の過程では、どの方法を取っても小さい数で大きい数を割る割り算が非常に多く、こうした数値情報が落とされたことによる誤差が大きく出ることになります。誤差の大きさは、元の行列と計算した逆行列を掛けたら、右から掛けても左から掛けても単位行列になることで検証ができます。LU分解も余因子行列よりは計算量が少ないため、こうした桁落ちによる誤差が少なくなるのですが、それでも、誤差が0になることは望めません。

そこで、もう少し誤差の発生を抑えられることを期待して、全ての計算を掛け算と足し算のみで完遂できるように、掃き出し法をもう少し工夫して計算できるようにしてみました。それが、今から紹介する行列計算による掃き出し法です。

掃き出し法は、計算量が多く、LU分解や余因子行列を用いる方法よりも誤差が大きいとされるため、逆行列計算には不向きだとされた歴史的経緯があります。しかし、コンピューターの発展とともに計算スピードが向上し、今では、計算量の比較的多い回り道を要する計算も、一瞬で計算できるようになっています。

何より、Excelで何百行何百列の行列計算を、一瞬で可能になることが分かりましたので、これを使えば、かなり見やすい計算式にすることができると思います。まだ、LU分解による逆行列の計算の方法と比較して、逆行列の精度を検証してはいないのですが、比較的精度の高い計算法として、有用な可能性があると思いましたので、ここで、Excelによる行列計算掃き出し法(オリジナル名、、)を提案してみようと思います。

理論的裏付けがあまりなく、本当なのか？と疑われる側面もあるとは思いますが、今しばらくお付き合いくださると幸いです。

行列計算掃き出し法による逆行列の計算

これから紹介する方法については、元行列が200行×200列だと、かなり時間がかかることになると思います。

逆行列を計算してから、算出した逆行列を掛けたいベクトルや行列に逆行列を掛ける方法よりも、元の行列と逆行列を掛けたいベクトル(行列)をそれぞれ並べて、掃き出し法の要領で、同じ行列を掛け算して、行列変形していきます。
こう計算すると、逆行列の割り算で整形した値を、メモリに収まる定数値に置き換える(桁落ち)することなく、直接逆行列計算後のベクトル(行列)を計算で算出することができると思います。
今回は、逆行列を左から掛ける計算に限って使うので、並べる行列は左右に並べて、左から列毎に掃き出して、計算することになると思いますが、逆行列を右から掛けてベクトル(行列)を求める場合は、上下に並べて、上から行毎に計算するといいと思います。

そのため、今回は掃き出し法を行列の掛け算にて再現した方法を取るのですが、これが200列もあると、1列の計算にかなり時間がかかるようで、手元のコンピュータを使用しても、5〜6分かかりました。

おそらく全ての固有ベクトルを求めるならば、200回反復作業を行う必要があると思いますが、これについてはマクロを組んで、200回全てを自動化できるように組んだ方がいいかと思います。ただし、時間は1日以上かかるのを覚悟です。。100列くらいまでならおそらく時間は少なくとも1/4には縮まり、ある程度有効な計算になるかと思います。。

用意するもの

※Excelワークシート3枚、元行列A、求めたい固有ベクトルに対応する固有値λk

①シート9：計算結果ペースト用シート
元行列Aと、固有値λkから計算した、AｰλkEの行列を最初に記入します。行列の左端の列の隣の列に、固有ベクトルを計算するために、先に説明した任意のベクトルbに対応するベクトルの初期値を記入します。このマスの初期値は、0ベクトルでなければどんな数値でも構いません。ただ、数値が大きくなりすぎると見にくいので、ユークリッド距離(要素の2乗和の平方根)が1くらいのベクトルにしておくといいかもしれません。最後にユークリッド距離による割り算を行うので、必須ではありません。

②シート10：掛け算用ベクトル書き出し用シート
このシートは、空白にしておきます。マクロで、掃き出し法を再現するための行列を記入することになります。このシートで書いた行列は、シート10の行列に左から掛けて、計算を進めます。反復する毎にシート内の数値を全て空にしてから、使う必要があります。

③シート11：計算結果書き出し用シート
このシートも空白にしておきます。マクロで、シート10とシート9の行列の掛け算の結果を書き出し、全て書き出した後に結果をシート9にコピーすることになります。反復する毎にシート内の数値を全て空にしてから、使う必要があります。Excelの行列計算をマクロで行う際には、計算結果のセル数が全て合わせて5461セル以上になるとエラーになるので、注意して使う必要があります。この行列計算の段階でエラーが出るほど大きな行列から、逆行列計算を行う時には、シート9について1行ずつ縦ベクトルを取り出して、順に計算することになります。Excelのセルに記入するための関数とは少しルールが違うので注意します。行列の積の関数は、WorksheetFunction.MMULT()に対応します。

手順：

上記の通り、ワークシートを3枚用意してください。元行列Aと、求めたい固有値λkから、AｰλkEを計算しておき、シート9に記入しておいてください。

ここからはしばらく、VBAによるマクロを使った計算の手順となります。正方行列の行数に対応した回数分だけ、反復して処理を行います。
まずシート10に掃き出し法の計算に対応した行列を作るために、値を書き出します。

現在までの反復回数をmと置いて、現在m回目の反復だとします。元の行列Aがn行×n列だとすれば、シート10のn行×n列分のマスを使って、行列を書き出すことになります。

m−1行目までの対角成分のマスには全て1を、それ以外のマスには、全て0を書き出します。m行目には、シート9のm行目のn個の数値(縦ベクトル)を参照して、それぞれの値から計算した数値を、各列に書き出します。

シート9の対角成分をa、シート9のl(エル)行目m列目の数値をbとします。そして、シート10のm行目m列目の対角成分のマスの数値は1/aを記入します。シート10のl行目m列目の対角成分以外のマスの数値に、-b/aを記入します。

この時、対角成分に0が入る時用に、対角成分が0になったら、シート9のm行目以降の別の行とm行目を入れ替える、ピボット処理を入れてもいいかもしれません。

また、m以外の数iについて、i行目m列目のシート9のマスの数値をc、i以外の行のl行目m列目のシート9の数値をdとした時に、m行目m列目のシート10のマスに1を、i行目m列目のシート10のマスを−1に、l行目m列目のシート10のマスに-d/cを、m列目i列目のシート10のマスに1/cを記入しても、同じ処理が達成されます。

m行目以降のシート10の数値については、全ての対角成分に1を、ピボット処理のために数値を入れたマス以外の、全ての対角成分以外のマスに0を記入します。これで、シート10に書き込むべき行列が全て完成します。

その後、シート11に、マクロでシート10の行列にシート9の行列を掛けた結果を書き出します。マクロのワークシート関数の制限で、結果行列にエラーが出る場合(結果行列が全てのセルを合わせて5461セル以上の場合)は、シート10の行列はそのままにして、シート9を1列ずつ、縦ベクトルを取り出して、行列の掛け算を行うようにしてください。

WorksheetFunction.MMULT(シート10,シート9(のうち1列))

このようにすることで、理想的には、5461行×5461+α列の行列について、ちゃんと行列積を求めることができます。ただし、5461行×5461列もの行列を扱った場合、1日経ってもマクロ計算が終わらないくらいに、長い計算になるかと思います。

シート9について、1列ずつ縦ベクトルを取り出してシート10の行列を掛け算すれば、全てのシート9のマスから構成される行列に、左からシート10の行列を掛け算した結果が出力されるようになります。

これを、シート9について1行ずつ横ベクトルを取り出して、シート10の行列を右から掛け算すれば、シート9の行列に、シート10の行列を右から掛け算した結果が出力されることになります。これは、逆行列を左から掛けるか、右から掛けるかに合わせて、切り替えることができるので、覚えておくといいかもしれません。

最後に、シート11に全ての行列計算の結果が出揃ったら、シート11の値を全てシート9にコピーペーストしてください。そして、シート10とシート11の数値を全てクリアしたら、m回目の反復が終了です。求めるつもりであったベクトルbの行列計算を忘れないようにしてください。

以上までの反復を1からnまで全て行えば、行列計算が終了です。シート9には、AｰλkEが最初に記入してあったマスには、n行×n列の単位行列が、ベクトルbが最初に記入してあったマスには、bにA-λkEの逆行列を左から掛けた計算結果が、縦ベクトルとして残っているはずです。値を取り出して、以降の計算に使用してください。

ただし、この反復を、1セット全て行った後には、シート9の単位行列が記入してあるマスについては、10の-20乗くらいの小さな値がそれぞれ対角成分以外のマスに残っていると思います。これは、おそらく全て、行列積の計算上の誤差のはずですので、もう一回同じ反復をしてみてください。

対角成分以外のほぼ全ての成分がゼロになり、かなり高い精度でベクトルbが求められていると思います。

以上が、掃き出し法を再現するための行列計算の手順となります。行列計算掃き出し法と仮に名付けることにします。計算手順の中にメモリの制限のあるマスを多数使うので、誤差が生じるのは相変わらず免れませんが、シート10に値を書き込む時には割り算を使うものの、肝心のシート9とシート10の行列計算では、全て同じ値を用いて、足し算と掛け算のみを使い、行列計算を行なっています。

大きな検証は未だ行なっていませんが、計算後のベクトルb‘をシート9に書き込んだ元行列に右から掛け算した時に、元ベクトルのbが再現できるところまでは確認しております。

ちなみに、1個注意して欲しいのですが、これは、ベクトルbのマスに、n行×n列の単位行列を用意して計算すれば、逆行列をそのまま求めることができます。しかし、こうして求めた逆行列を元の行列と掛け算した時、誤差が出てることも確認しております。

逆行列の計算を回避して計算するのには実際有用だとは思うのですが、計算した逆行列に対して誤差がないとするのはちょっと信頼性に疑問がありますので、逆行列の計算用には用いないでください。ここはもう少し検証してみます。

逆行列計算から、固有ベクトルの算出まで

ここまで完成したら、後一息です。ただし、固有値全てに対応する固有ベクトルを求める場合は、逆行列計算を少なくとも、行数分×2回は行う必要が出るため、計算時間が非常にかかります。

固有値に逆行列を掛ける計算は、2回3回と反復するほどに固有ベクトルの精度が良くなるので、2回以上反復して、同じ逆行列を、左から掛け算してください。

この計算が終わった後のベクトルb‘’について、自身のユークリッド距離(各要素の2乗和の平方根)を求めて、その値で各要素を割り算すれば、求めたベクトルはユークリッド距離1の単位ベクトルとなり、固有値に対応する固有ベクトル(に限りなく近い近似値)になります。

値ぴったりではないはずですが、これがQR法をベースにした固有値固有ベクトル計算の全てとなります。QR法は、固有値計算としてはかなり計算結果が数値的に安定した計算方法となるそうなので、よければ試してみてください。

この記事が気に入ったらサポートをしてみませんか？