3 (基本演算)線形最小二乗問題

3.1 概要

CPPLapackでは，線形最小二乗問題を解くための関数が各行列クラスのメンバ関数として用意されています．

      dgematrix.dgels
      dgematrix.dgelss
      zgematrix.zgels
      zgematrix.zgelss

3.2 使い方

ここではdouble型の実一般行列を例に紹介します．
Ax=y形の線形最小二乗問題を解く時，以下の様に記述します． AはM行N列行列であり，解xはyに上書きされます．

      A.dgels(y,r);
      A.dgelss(y,S,RANK,RCOND);

dgelsは $m\geq n$ の場合には最小二乗解を， $m\leq n$ の場合は最小ノルム解を求めます． rには計算残差が入ります．rは引数に取らなくてもでも構いません．
dgelssは最小ノルム解を求めます． Sは特異値を出力するdcovector型， RANKはAのランクを出力するdouble型です．RCONDはAの実行精度を決めるdouble型で，デフォルト引数で-1.0(機種精度)に設定されているので，ここは空欄でもかまいません．
詳しくはマニュアルを御覧下さい．

3.3 注意点

計算後，行列AはdgelsではQR,LQ分解の詳細に， dgelssでは右特異ベクトルへと上書きされ，引数yは解xとして上書きされます．元の値を保持したい場合は，関数を呼び出す前にコピーを作っておいて下さい．
dgelsではフルランクを仮定しています．行列Aがランク落ちの可能性がある場合はdgelssを用いてください．

3.4 例題

連立一次方程式

$\displaystyle \begin{bmatrix} 4 & 1 & 3 \ 5 & 1 & -1\ 2 & 2 & 0 \ -3 ... ...ght) = \left( \begin{array}{c} 8\ 2\ 0\ 5\ \end{array} \right)$

を解く．

コードを以下に示します．

//====================================================================[include]
#include ``cpplapack.h''

//=======================================================================[main]
/*! main */
int main(int argc, char** argv)
{
  //// make dgematrix A  ////
  CPPL::dgematrix A(4,3);
  A(0,0)=4;  A(0,1)=1;  A(0,2)=3;
  A(1,0)=5;  A(1,1)=1;  A(1,2)=-1;
  A(2,0)=2;  A(2,1)=2;  A(2,2)=0;
  A(3,0)=-3; A(3,1)=1;  A(3,2)=5;

  //// make dcovector y ////
  CPPL::dcovector y(4);
  y(0)=8;
  y(1)=2;
  y(2)=0;
  y(3)=5;

  //// solve Ax=y ////
  A.dgels(y);

  //// print ////
  std::cout << ``x=\n'' << y << std::endl;

  return 0;
}

実行結果は以下の様になります．

Masashi UESHIMA