Subsections

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

3.1 概要

CPPLapackでは,線形最小二乗問題を解くための関数が各行列クラスのメンバ関数として用意されています.
      dgematrix.dgels
      dgematrix.dgelss
      zgematrix.zgels
      zgematrix.zgelss

3.2 使い方

ここではdouble型の実一般行列を例に紹介します.
Ax=y形の線形最小二乗問題を解く時,以下の様に記述します. AはM行N列行列であり,解xyに上書きされます.
      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型です.RCONDAの実行精度を決めるdouble型で, デフォルト引数で-1.0(機種精度)に設定されているので,ここは空欄でもかまいません.
詳しくはマニュアルを御覧下さい.

3.3 注意点

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;
}



実行結果は以下の様になります.
x=
 0.928934
 -0.92422
 1.74003



Masashi UESHIMA