最小2乗法 - その11

たまたま人づてに前の会社で僕が最後にやっていたテーマが今どうなっているか聞いた。評価用サンプルを作る段階まで到達して、うまく行っているらしい。それはいいんだけど、最後に僕がやってた近似計算は結局無視されて、高価な照明系シミュレーションのソフトウェアを導入して計算しているという。

照明系の計算ソフトは光線追跡に基づいていて、一番の問題となる回折効果（SWS的な効果も含む）は別の手段を使わないと盛り込めないので、光の利用効率は悪くないという結果にしかならない、従って無駄であってだから近似的な計算を僕がやる、と言ってその計算を実行して、後のために詳細なレポートにまとめた。それが去年の暮れから年明けにかけての話だった。

僕がいなくなったあとそのソフトウェアのライセンスを買って光学的なモデルを作って計算すると、効率はまったく劣化しないという結果になって、試作することになった、試作品を測定すると計算よりは悪いと言う結果だったが、そのまま行けということになった、という。ふーん、それはよかったですねえ、と笑顔で聞いていたけど、腹の底では、まったくもって僕が言ったことを無視して、しかも僕が苦労してやった計算まで無かったことにしやがって、とムカムカした。

つまりようするに、近似ではあるが本質的な身内の計算より、結果に影響しない複雑なモデリングをして高価なソフトを使った計算のほうを信用する、ということだ。そしてこれが技術を売りにしている会社の、技術の現場で行われていることだ（これまでの僕の言動にも問題がある、という指摘は甘んじて受けよう。しかしそれはそれこれはこれだ）。おそらくそのうちこの会社では、節電のためにブレーカにシールを貼ったり、謹製のマイナスイオン応用製品が出たりするんだろう。

まあ、今さらどうなろうが僕の知ったことではないけど、なんだかおも〜い無力感を感じた。それでは好きなだけ無駄な金と時間と工数を費やして無意味な結果を出せばいい。ほんとに僕の知ったこっちゃない。

．．．．さて、（はぁ。ため息が出てしまった）気を取り直して最小2乗法の実装。

5  Objective-Cでの実装

Householder変換でQR分解して不能の方程式の最小2乗解を求めるObjective-Cのクラスを書いてみる。

Objective-Cとは関係ないんだけど、OS Xでこういう数値計算をするときはvDSPライブラリを使ったほうがいい場合がある（単精度計算では非常に有利。倍精度ではとんとんぐらい）。

Householder変換ではベクトルの内積の計算をいっぱい使うことになるけど、内積計算はvDSPにある。

void vDSP_dotpr (
   const float __vDSP_input1[],
   vDSP_Stride __vDSP_stride1,
   const float __vDSP_input2[],
   vDSP_Stride __vDSP_stride2,
   float *__vDSP_result,
   vDSP_Length __vDSP_size
);

void vDSP_dotprD (
   const double __vDSP_input1[],
   vDSP_Stride __vDSP_stride1,
   const double __vDSP_input2[],
   vDSP_Stride __vDSP_stride2,
   double *__vDSP_result,
   vDSP_Length __vDSP_size
);

vDSP_dotpr()は単精度用、vDSP_dotprD()は倍精度用の内積関数である。

また、Objective-CではCとコンパチで行列すなわち2次元配列は行（横方向）が連続に並んでいる（row major）とみなすのが自然である。しかしHouseholder変換の計算では行列演算よりむしろ列ベクトルを扱うことのほうが多い。データのメモリ上の局所性と言う観点からは、FORTRANのような列方向が連続なほう（column major）が有利である。しかし、とちゅうで転置するなんていうことをやっているとそのオーバーヘッドが馬鹿にならなくなるので、行連続のまま列ベクトルを扱う必要がある。

vDSPではベクトル（1次元配列）の要素は連続に並んでいなくてもいい。その点でもvDSPを使うのが便利である。ほんとはCBLASを使うのが効率上はいちばんいい、ということがわかってるんだけど、CBLASのベクトルは連続でないといけないので、残念ながら使えない（もちろん列連続なデータにすればCBLASが使える。しかしそこまでするんならLAPACKを使えよ、という感じになって独自に実装する意味は「ボケ防止」以外になくなってしまう）。

数値計算の実装でいつも悩ましいのは、ベクトルや行列を単なるCの構造体、あるいはもっとただの配列へのポインタにするか、あるいはObjective-Cのオブジェクトにするか、ということ。これまでいろんなパターンを書いてきたけど、いまだにどっちがいいか決められない。

オブジェクトにしたほうが生成破棄が簡単になって便利だけど、実際の計算、例えば行列とベクトルの積をメソッドで書いたからといってわかりやすくはならない。むしろ静的変数を持たない関数として実装したほうが素直な感じがする。それは計算の中身がよくわかっていて、わざわざオブジェクトとしてまとめるのがうるさい、と感じるからであろう。計算のオーバーヘッドはゴミ程度の差しかないはずなので、その「うるささ」は単に思い込みでしかない、ということではあるんだけど。

5.2  オブジェクト設計

実装に先立って、おおよその設計をしておこう。とりあえず今回はベクトルやマトリクスをオブジェクトとして扱ってみよう。理由は特にない。これまで計算効率を気にしてあまりオブジェクトにまとめることをせず、単なるメモリの領域として実装することのほうが多かった。しかし大昔ならいざしらず、GBクラスのマトリクスを扱うのでなければ、計算効率をそれほど気にする必要は無くなってきた。それよりもデバグのしやすさ、見通しの良さを優先するほうがいい。だからといってベクトルやマトリクスをオブジェクトにすればわかりやすくなるか、というと実はそうでもない。

今回はひとつの数値計算の実装の実験として、ベクトルとマトリクスをObjective-Cのオブジェクトとして書いてみることにする。

それだけではつまらないし、オブジェクトに書かずにこれまでやってきたような方法で書いたものと比較するのなら意味はあるけど、そこまでやる元気はないので、今回は要素の値が単精度（float）の場合と倍精度（double）の場合とでなるべく同じメソッドを使うようにしよう。

つまり、オブジェクトを生成するときに単精度か倍精度かを決めると、あとは同じメソッドが使えるようにする。もちろんそれぞれに同じ名前のメソッドを別々に実装すればそうなるんだけど、なるべく共通の作業は一カ所にまとめるようにしよう。

具体的には、例えばベクトルのクラス

@interface HHVector : NSObject {
    NSUInteger  eCount;
    NSInteger   stride;
    void        *data;
}

を抽象クラスとして、

@interface HHFloatVector : HHVector
....
@interface HHDoubleVector : HHVector
....

が実際の作業を行うクラスにする。そして使う側からは抽象クラスとして扱えるようにする。これは今までときどきやってきた「なんちゃってクラスクラスタ」に近いけど、数値を保持するオブジェクトはその型、単精度か倍精度かというのを完全に忘れることはできないので、それほどはメリットが無いかもしれない。でも単精度と倍精度を完全に別々に実装するよりはデバグが簡単になるだろう。

また、単精度と倍精度が混在した計算では自動的に型変換が行われるようにできる。まあ、実際にはマトリクス全体の型変換が発生するような計算はそもそもやるべきではないので、これも実際にはそれほどメリットは無い。まあ、実験としてわりきろう。

抽象クラスと言いながら、共通する作業はなるべくまとめて、型に依存する部分だけを子供クラスにやらせるようにする。そのためにいくつか制限を設けることにする。

• 要素の型はinitで決定されて、変更することはできない
• 要素数やストライド（要素の間隔）も変更することはできない

とする。違う型や要素数の変更は新しいオブジェクトを作って、要素をコピーする、ということにする。

また、それとは別に要素を共有する場合も考慮する。それはつまり、マトリクスの行ベクトルをとりだしたとき、新しいベクトルを作るのではなく、マトリクスの要素を使ったベクトルを作ることにする。そうすることでマトリクス間やマトリクスと行列の演算をベクトルの演算の組み合わせとして実装することができる。

例えば、マトリクスどうしの積ABを

ではなく、行ベクトルと列ベクトルの内積

として計算することができる。ここで^ra_jはマトリクスAのj番目の行ベクトル、^cb_iはBのi番目の列ベクトルである。

これは計算効率の点から言えば不利だけど、ネストしたループを書かずにすむ。とはいうもののマトリクスどうしの積ぐらいではバグを埋め込むほどではないし、vDSPを呼べば直接計算できるので、これもそれほどメリットとは言えないけど、頭の体操と考えよう。

そういう動作を実装するには、データ領域を自分で確保したかどうか、というスイッチを持っている必要がある。自分で確保したらオブジェクトがリリースされるときにその領域も解放するが、マトリクスの列ベクトルだった場合はデータ領域は解放しない、ということになる。

ごちゃごちゃ書いてきたけど、ようするにやってみたい、というだけ。まあ、いいじゃん。