高次の最小2乗法の計算

高次の最小2乗法の計算に関してですが、
https://imagingsolution.net/wordpress/wp-content …

この数式のように、行列計算により求めることができます。

VB.NETに行列計算のライブラリであるMathNet.Numericsを使ってプログラムを書いたのですが
10次くらいまでは計算できるのですが、
それ以上になると、要素数が多すぎるせいか、
計算できずエラーとなります。

どのようにすれば、10次以上の計算をすることができますか？

行列を小さく分解する方法などあるのでしょうか？

質問者からの補足コメント

• ありがとうございます。
  
  たくさん質問したいことがあるのですが、
  
  > という推定式を求めたいのですか?
  仰る通りです。
  
  
  > ①高次項を作る時は、必ずxを中心化してから累乗すること。
  中心化、というのは、オフセットをさっ引くという意味ですか？
  y＝a＋f
  の式のaで求まるのではないのでしょうか？
  
  
  
  > ②20次にもなるとサンプル数は200くらいないと、線形制約が入り逆行列が解けなくなります。
  何次以上だと何サンプル必要といった式はありますか？
  あと、線形制約とはどのようなものでしょうか？
  
  
  > 高次回帰の他にも色々あります。
  例えば、どういったものがありますでしょうか？

  No.1の回答に寄せられた補足コメントです。 補足日時：2019/06/07 21:50

• ありがとうございます。
  
  中心化に関してなのですが、
  これってつまり、
  y＝a＋f（x,x^2,x^3,x^4,x^5,・・・）
  のうち、f(x)の項だけで先にフィッティングを行い
  さっ引いた後に、それ以外の項をフィッティングすれば
  集束しやすいという意味ですか？
  
  もし、そうだとして、0次項（a）は先に求めなくて良いのでしょうか？
  
  何となく0→1→2→3→・・・
  と順番にフィッティングしていけば
  集束しやすいのではないかと予想するのですが、
  そういうことはないのでしょうか？

  No.2の回答に寄せられた補足コメントです。 補足日時：2019/06/14 21:14

No.2 回答者： kamiyasiro 回答日時： 2019/06/09 11:22

#1です。

なにか事例のデータはありますか。それで説明するのが早いのですが、概念で説明します。

・線形制約とは、
もし、各列のデータをA,B,C,Dとすると、A＋B＋C＝Dのような線形の等式が成立することを線形制約と言い、行列のランクが落ちます。その時は行列式の値（デターミナント値）が0になります。
もとの行列をXとすると、そのランクは、
rank(X)＝rank(XTX)＝rank(XXT)・・・Tは転置
なので、XTX、一般的には相関係数行列のdet値を調べます。ランク落ちするとdet値がほぼ0になります。このときは、逆行列が解けないため重回帰分析が破綻します。

・中心化とは、
1,2,3の2次項をそのまま作成すると、1,4,9となり、１次項と2次項は強い相関を持ちます。しかし、１次項を一旦中心化し（平均を引き）、-1,0,1として2次項を作ると1,0,1となり、１次項と２次項の内積は0になります。すなわち直交化します。このように、高次項が１次項と相関を持つことを避けることによって、高次の効果が解けるようになります。奇数次項の相関は無くなりませんが、強い相関は無くなります。
中心化あり／なしで、高次項まで作成したXの相関係数行列のdet値を見れば、中心化の効果が分かります。

・逆行列が解けないとは、
（逆行列）＝１／det(元の行列)・（余因子行列）
なので、det値が0あるいはニアリィ0だと、割り算の項が無限大となり、逆行列が存在できません。

・予期せぬ線形制約とは、何次だと何個サンプルが必要か、
次数が200次元もあると、相関係数行列の下半分は200Ｃ2（200から2を取る組み合わせ）で、19900もの数値があります。もし、サンプルが1000個しかないと、18900もの数値は全て従属となり、線形制約のないデータであってもdet値はほぼ0になります。例えば、10次元でも、10C2=45ですから50個程度のサンプルがないと解けません。それでもオーバーフィッティングになっています。ですから、200個程度のサンプルが必要です。

・高次回帰（多項式回帰）以外の代表的な手法は、
よく用いられるのが、gam（一般化加法モデル）です。ガウシアンの重畳としてモデルを考えます。そのほかにはBスプライン、罰則付きスプライン、ハンペルフィルタ、カルマンフィルタなどがあります。

No.1 回答者： kamiyasiro 回答日時： 2019/06/06 00:06

企業で統計を推進する立場の者です。

y＝a＋f（x,x^2,x^3,x^4,x^5,・・・）+ε

という推定式を求めたいのですか?
通常、説明変数行列をX、係数ベクトルをβとすると、

β＝(XTX)^-1・XT・ｙ　ここでTは転置

という式になるんですが・・・
まあ、それは良いとして、「～が特異です」のような逆行列の部分のエラーですよね。

①高次項を作る時は、必ずxを中心化してから累乗すること。そうしないと、容易に線形制約が入り逆行列が解けなくなります。
②20次にもなるとサンプル数は200くらいないと、線形制約が入り逆行列が解けなくなります。

逆行列を解く前の、大きなカッコのdet値（行列式の値）を取ってみて、10のマイナス5乗より小さいと、逆行列が解けたとしてもかなり怪しい解になります。

対策としては、
①グネグネ曲線の近似方法は、高次回帰の他にも色々あります。
②また、det値が0に漸近し、X行列がランク落ちに近い時でも、「正則化」（スパースモデリング）を使えば解けます。