最小2乗について

片山 泰男 (Yasuo Katayama) 2011/11/16



データに合わせてエラーを最小にする、目的をもって何かを設計する最小2乗法といわれる方法がある。まずは、最も単純な例を挙げよう。

x に対する y という、(x, y) の多くの測定があって、これらを比例 y= ax で近似する係数 a を求めたい。誤差を E= Σ (y- ax)^2 とすると、 E は a を変数としてみるとき2次関数になる。E(a) は a のある値 a0 で極値をとり、そこは、∂E/∂a= 0 という、aの変化がEに影響しない点に なっている。多項式で極とは微分係数が0でその前後で符号が反転し、しなければ変曲点である。極値は、最大又は最小であり得るが、aが+∞や-∞ では E も+∞になるから、E(a0)は、最小値である。E を展開し E = Σy^2 + a^2 Σx^2 - 2aΣxy であり、E の a についての微分、∂E/∂a= 2aΣx^2 - 2Σxy を =0 とおき、a= Σxy/Σx^2 がすぐに得られる。偏微分=0 が E の最小値を与える a の式を導くのである。

いくつかの仮定が、これを支えていた。(1)エラーが(根拠なく)2乗誤差であること。エラーは、なぜ2乗誤差か、誤差の定義に正解はないのだろうか。 ここでは、最小2乗法を使うためだけに2乗誤差が使われたのである。例えば、2乗誤差でなく絶対値誤差Σ|y-ax|がよい場合がある。例えば、動き ベクトル検索では画像誤差を、E= Σ|I(x, y)- P(x+mvx, y+mvy)| として、2次元の動きベクトルの関数として E を最小にする動きベクトルを全数 探索することが多い。勿論、全数探索以外にも多くの方法があるが、フルサーチは、探索範囲内で最小のエラーをもつずれの位置を捜し出す。 それ以上に良い結果は、動き場(動きベクトルの空間分布)にLPFしてスムースな動きにしてその符号量を削減することや、動きベクトルの符号量を エラーに導入する方法以外に考えられない。

そのとき、2乗誤差よりも絶対値誤差の方が圧倒的に使われる。その理由は、探索は単純に探索であって、最小2乗法でないから2乗を使う利点がなく、 そして2乗は、2乗が値のビット数を倍にし、絶対値回路と比べ2乗回路はその規模が数倍に大きい。ソフトウエアでは絶対値と2乗は複雑さに差が小 さく、テーブル索引では差はもっと小さい。両者は実用上の動き検出能力に大きな差はないが、2乗がSNRの定義に合致するから SNR的に最適なのは 2乗誤差である。絶対値が 2乗の仕方なく使う簡略した代替という面はあるが、それだけではない。2乗は画素に大きく広がる低域の誤差よりも高域 又はパルス的誤差を優先して検出する欠点がある。2乗和とは誤差の絶対値にそれ自身を重み付けた和であるからである。2乗平均は、平均とは大きく 違う。例えば、4箇所に1がある(1,1,1,1)を1箇所に集中する(0,0,0,4)と平均は1のままで2乗平均は√(16/4)=2になる。2乗誤差最小は、誤差を分散さ せるだけかもしれない。



逆に利点としては、ソフトウエアでは、誤差がマッチング窓内の画素間の誤差を累積して行われ、部分和が暫定の最小値を超えた場合、すぐに累積を 中止して他の候補に向かうという省略手法を利用することが多いが、2乗が値の範囲を広げるために、2乗探索の最小値と平均値との比率は、絶対値に よる探索より小さく、誤差の累積過程で暫定の最小値の超過の検出を早めることもよく知られている。各候補ベクトルでの処理として、暫定の最小値 を超過しない間は誤差の累積を繰り返し、暫定最小値を超えると(最小値の可能性を失うため)中止する。それ以前に終了すると、暫定最小値を更新する。

音声など時系列や画像において、誤差の定義の評価は確定せず、例えば、画像の動きベクトル探索では、対応する画素を独立した次元としてみる空間 の距離の誤差は、画素間の連接を全く無視する点が不満である(*)。

(2)エラー関数の滑らかさを前提にする。最小2乗に必要な誤差関数は、2階微分可能、又は微分可能で微分が連続である。そのとき微分0の場所に極が あると考えられる。しかし、個別の問題においてその条件は、満たされない。例えば任意の画像のずれの誤差である。

(3) 多くの場合、求めたいパラメータcは、1次元でなく多次元であり、例えば動きベクトルでは2次元である。2次元に張るスカラー関数として誤差が ある。それは空間に分布するポテンシャルであり、最適点とは、場所を少々変更しても変わらないポテンシャルをもつ点である、このことを使って最 適点を求めることができる場合がある。物体は、最小作用の原理で静止し、運動は、そのような運動だけが実現する。そのような最適化問題としての 物理がある。絶対値誤差では、誤差最小点は、そこでの微分が2方向からの値が一致しない点になって、式を導くのに不便である。また、最適点で傾き が0であることは、誤差の大きさに対する場所の誤差が大きいという性質を想起させ、いやなものである(**)。誤差の違いが大差ないときに、最小点を求 め続ける意味はあまりない。大差ない誤差には、場所の精度の必要性もないからである。画像自体が平坦な領域ではそのようなことが起き、動きベク トルは不必要にばらつき広がる。その動きベクトルの符号化に必要な符号量は無駄になる。

(*) 現在、2乗誤差(MSE)、絶対値誤差(SAD)ではなく、画像の画素連接を使う誤差、アダマール変換(4x4)後の絶対値誤差(TSAD)が使われる。直交変換は ユークリッド距離を変えないため直交変換後のMSEは変換の意味がない。MSEとは変換の意味を失うような、まずい誤差(歪)の定義であった。アダマール 変換はDCTの代用加減算であり、直交変換前の予測の選択にMSEやSADを使っていた予測の選択、例えば、MPEG-2のフレーム/フィールド予測の選択、B画像 の2予測の4選択等は、TSAD判断すれば符号化効率が向上するだろう。符号量を使うほうが良いが、そこまでしないときの符号量と歪の最適化 RDOである。

(**) 定義域の側の違いは符号量に影響する。RDOの見地から、誤差を少しだけ減らすのに大きな符号量を使うのは無駄である。



もしも、xを誤差とすると、 E= Σ (x- cy)^2 であり、∂E/∂c= 2cΣy^2 - 2Σxy を=0とおくと、c= Σxy/Σy^2。これは最初の例のaの逆数でない。 つまり、同じ結果を与えない。これはサンプル(x, y)から回帰直線への水平の足の長さを最小にする話である。最初の話は、垂直の足の長さを最小に する話であった。傾きを求めるのに、yの誤差Δy^2を使うか、xの誤差Δx^2を使うかによって結果は違い、さらに、x,yに対称であるサンプル(x,y) の座標と回帰直線との距離、d= ΔxΔy/√(Δx^2+Δy^2)の2乗は、また違う解を与える任意性がある。何を最小にしようとしているかが違うなら解も 異なるのである。

最初の簡単な例に定数項 b を付けると、もう複雑になる。 1次式 y= ax + b で (オフセットbも含めて) 近似するには、誤差を E= Σ(y - ax - b)^2 として、
E を展開すれば、E= Σy^2 + a^2Σx^2 + Σb^2 - 2aΣxy - 2bΣy + 2abΣxであり、
∂E/∂a= 2aΣx^2 -2Σxy + 2bΣx、
∂E/∂b= 2Nb - 2Σy + 2aΣx、 (N:測定点数)
これらふたつの偏微分をそれぞれ0とする条件は、丁度、パラメータ次元数の連立方程式になって、

|Σx^2, Σx| |a| = |Σxy|
|  Σx,   N| |b|   |Σy |
面倒な連立方程式は、ガウスの掃き出し法など機械的に解く方法があり、それに頼ればよい。結果は、 a= (NΣxy - ΣyΣx) / (NΣx^2 - (Σx)^2)、そして b= (Σx^2Σy - ΣxyΣx) / (NΣx^2 - (Σx)^2) である。



Wiener の最適フィルタについて、

劣化信号xの対象画素x(k)の隣接画素x(k-1),x(k+1)からの N(N=3) タップフィルタによって、原信号yに近付けるフィルタ y~(k)= a x(k-1) + b x(k) + c x(k+1) を考え、フィルタ係数 (a, b, c) を求める。

E= Σ{y(k) - a x(k-1) - b x(k) - c x(k+1) }^2
= {Σy(k)^2 + a^2Σx(k-1)^2 + b^2Σx(k)^2 + c^2Σx(k+1)^2 }
-2{aΣy(k)x(k-1) + bΣy(k)x(k) + cΣy(k)x(k+1) }
+2{abΣx(k-1)x(k) + acΣx(k-1)x(k+1) + bcΣx(k)x(k+1)}

∂E/∂a=0 から、aΣx(k-1)^2 + bΣx(k-1)x(k) + cΣx(k-1)x(k+1) = Σy(k)x(k-1)
∂E/∂b=0 から、aΣx(k-1)x(k) + bΣx(k)^2 + cΣx(k)x(k+1) = Σy(k)x(k)
∂E/∂c=0 から、aΣx(k-1)x(k+1) + bΣx(k)x(k+1) + cΣx(k+1)^2 = Σy(k)x(k+1)

3式は、自己相関行列×ベクトル(a,b,c)=相互相関ベクトルという、次の連立方程式になる。

|Σx(k-1)^2,	Σx(k-1)x(k),	Σx(k-1)x(k+1)| |a|   |Σy(k)x(k-1)|
|Σx(k-1)x(k),	Σx(k)^2,	Σx(k)x(k+1)  | |b| = |Σy(k)x(k)  |
|Σx(k-1)x(k+1),Σx(k)x(k+1),	Σx(k+1)^2    | |c|   |Σy(k)x(k+1)|
左辺の大きな行列は、自己相関Σx(k+i)x(k+j) をi,j成分にする NxN 対称行列であり、自己相関行列(テプリッツ行列)という。 これは劣化信号xだけに関係する。右辺のベクトルは、目標フィルタ出力y(k)とベクトル(x(k-1),x(k),x(k+1))との相互相関Σy(k)x(k+j) を成分にするベクトルで相互相関ベクトルという。これは入力又は原信号yと劣化信号xの相互相関である。

この方法のフィルタの設計は、最小2乗法だけを使って、xが受けた劣化過程を全く考慮することなく、求める原信号yとの誤差を最小にする フィルタ設計を可能にする。この方法は、明らかに目的論的な方法論をもち、原信号を知る系は劣化を回復できるという考えは、汎用性がある。 例えば動画像符号化のエンコーダ内部のループ内フィルタには、信号によるフィルタ設計が可能であり、フィルタ係数は、デコーダに伝送できる。

最初の例の a= Σxy/Σx^2 は、DCT係数の量子化マトリックスQM_ij(i,j=0〜7)の最適化に使える。DCT係数をY_ij, 量子化係数をX_ijとし、 それらの統計が採られ、Σは例えば1画像とし、最適な量子化マトリックスは、 QM_ij= ΣX_ijY_ij/Σ(X_ij)^2。各係数位置(i,j)毎の Y_ij とX_ijとの相互相関とX_ijの自己相関との比である。人間の視覚特性(高域は粗い量子化に鈍感)を理由にせずにQMを正当化できる。MPEG-1,-2 ではマクロブロック毎の量子化変更ができ、Y_ijは mquant での量子化後の係数、X_ijはさらにQMでの量子化係数になる。