関連情報

線形最小二乗問題

Fortranによるサンプルソースコード
使用ルーチン名:e04ncf

Keyword: 線形最小二乗問題

概要

本サンプルは線形最小二乗問題を解くFortranによるサンプルプログラムです。 本サンプルは以下に示される二次関数を最小化する解を求めて、出力します。

線形最小二乗問題のデータ 

※本サンプルはNAG Fortranライブラリに含まれるルーチン e04ncf() のExampleコードです。本サンプル及びルーチンの詳細情報は e04ncf のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで

入力データ

(本ルーチンの詳細はe04ncf のマニュアルページを参照)

このデータをダウンロード
E04NCF Example Program Data
  10   9   3                                           :Values of M, N and NCLIN
  1.0   1.0   1.0   1.0   1.0   1.0   1.0   1.0   1.0
  1.0   2.0   1.0   1.0   1.0   1.0   2.0   0.0   0.0
  1.0   1.0   3.0   1.0   1.0   1.0  -1.0  -1.0  -3.0
  1.0   1.0   1.0   4.0   1.0   1.0   1.0   1.0   1.0
  1.0   1.0   1.0   3.0   1.0   1.0   1.0   1.0   1.0
  1.0   1.0   2.0   1.0   1.0   0.0   0.0   0.0  -1.0
  1.0   1.0   1.0   1.0   0.0   1.0   1.0   1.0   1.0
  1.0   1.0   1.0   0.0   1.0   1.0   1.0   1.0   1.0
  1.0   1.0   0.0   1.0   1.0   1.0   2.0   2.0   3.0
  1.0   0.0   1.0   1.0   1.0   1.0   0.0   2.0   2.0        :End of matrix A
  1.0   1.0   1.0   1.0   1.0   1.0   1.0   1.0   1.0   1.0  :End of B
  1.0   1.0   1.0   1.0   1.0   1.0   1.0   1.0   4.0
  1.0   2.0   3.0   4.0  -2.0   1.0   1.0   1.0   1.0
  1.0  -1.0   1.0  -1.0   1.0   1.0   1.0   1.0   1.0        :End of matrix C
  0.0       0.0      -1.0E+25   0.0   0.0   0.0   0.0   0.0   0.0
  2.0      -1.0E+25   1.0                                             :End of BL
  2.0       2.0       2.0       2.0   2.0   2.0   2.0   2.0   2.0
  1.0E+25   2.0       4.0                                             :End of BU
  1.0   0.5   0.3333   0.25   0.2   0.1667   0.1428   0.125   0.1111  :End of X 

  • 1行目はタイトル行で読み飛ばされます。
  • 2行目に係数行列Aの行数(m=10)、変数の数(n=9)、線形制約の数(nclin=3)を指定しています。
  • 3〜12行目には係数行列Aの要素を指定しています。
  • 13行目には係数行列Bを指定しています。
  • 15〜16行目には線形制約行列Cの要素を指定しています。
  • 17〜18行目には変数と制約の下限(bl)を指定しています。
  • 19〜20行目には変数と制約の上限(bu)を指定しています。
  • 21行目にはxの初期値の推定値(x)を指定しています。

出力結果

(本ルーチンの詳細はe04ncf のマニュアルページを参照)

この出力例をダウンロード
 E04NCF Example Program Results

 *** E04NCF

 Parameters
 ----------

 Problem type...........       LS1       Hessian................        NO

 Linear constraints.....         3       Feasibility tolerance..  1.05E-08
 Variables..............         9       Crash tolerance........  1.00E-02
 Objective matrix rows..        10       Rank tolerance.........  1.11E-14

 Infinite bound size....  1.00E+20       COLD start.............
 Infinite step size.....  1.00E+20       EPS (machine precision)  1.11E-16

 Print level............        10       Feasibility phase itns.        60
 Monitoring file........        -1       Optimality  phase itns.        60

 Workspace provided is     IWORK(       9),  WORK(     261).
 To solve problem we need  IWORK(       9),  WORK(     261).

 Rank of the objective function data matrix =     6


 Itn     Step Ninf Sinf/Objective  Norm Gz
   0  0.0E+00    1   2.145500E+00  0.0E+00
   1  2.5E-01    1   1.145500E+00  0.0E+00
   2  3.8E-01    0   6.595685E+00  2.3E+01
   3  1.6E-01    0   4.602812E+00  2.2E+00
   4  2.6E-01    0   4.283027E+00  1.6E+00
   5  4.4E-01    0   4.172700E+00  6.2E-01
   6  1.0E+00    0   4.165246E+00  2.7E-15
   7  5.8E-01    0   8.186698E-01  4.2E+00
   8  7.1E-01    0   2.575168E-01  1.2E+00
   9  1.0E+00    0   2.279362E-01  8.9E-16
  10  6.6E-01    0   9.629363E-02  8.9E-16
  11  0.0E+00    0   9.629363E-02  6.5E-16
  12  1.0E+00    0   8.134082E-02  2.4E-15

 Exit from LS problem after    12 iterations.


 Varbl State     Value       Lower Bound   Upper Bound    Lagr Mult   Slack

 V   1    LL    0.00000           .         2.00000      0.1572         .
 V   2    FR   4.152607E-02       .         2.00000           .      4.1526E-02
 V   3    FR   0.587176          None       2.00000           .       1.413
 V   4    LL    0.00000           .         2.00000      0.8782         .
 V   5    FR   9.964323E-02       .         2.00000           .      9.9643E-02
 V   6    LL    0.00000           .         2.00000      0.1473         .
 V   7    FR   4.905781E-02       .         2.00000           .      4.9058E-02
 V   8    LL    0.00000           .         2.00000      0.8603         .
 V   9    FR   0.305649           .         2.00000           .      0.3056


 L Con State     Value       Lower Bound   Upper Bound    Lagr Mult   Slack

 L   1    LL    2.00000       2.00000          None      0.3777     -2.2204E-16
 L   2    UL    2.00000          None       2.00000     -5.7914E-02  2.2204E-16
 L   3    LL    1.00000       1.00000       4.00000      0.1075      2.2204E-16

 Exit E04NCF - Optimal LS solution.

 Final LS objective value =   0.8134082E-01

  • 6〜7行目に線形制約の数、変数の数、目的関数の行列の行数が出力されています。
  • 8〜18行目に以下に示すプログラムのオプション引数が出力されています。
    Problem type 最小化される目的関数の種類。"LS1"は標準的な最小二乗問題を意味します。
    Hessian プログラムの引数hの内容を制御します。"NO"は引数hが変換され並べ替えられた行列HQのコレスキー因子を含むことを意味しています。
    Linear constraints 線形制約の数。
    Feasibility tolerance 実行可能解での許容可能な制約違反の最大値。
    Variables 変数の数。
    Crash tolerance 初期のワーキングセットに含まれる不等式制約が下限/上限のこの範囲内にあることを示しています。
    Objective matrix rows 目的関数行列の行数。
    Rank tolerance 三角因子の推定値を制御します。
    Infinite bound size 無限の境界値。この値以上の上限は+∞と見なされます。またこの値を負に反転させた値以下の下限は−∞と見なされます。
    Cold start 初期のワーキングセットがどのように選ばれるかを示しています。"Cold"は変数や制約の値に基づいて初期のワーキングセットが選ばれていることを意味します。
    Infinite step size 無限解へのステップと見なされる変数の変化の大きさ。
    EPS (machine precision) マシンの精度。
    Print level 結果出力のレベル。"10"は最終解と各反復の結果を出力することを意味します。
    Feasibility phase itns. 実行可能解生成フェーズの最大反復数。
    Monitoring file モニタリング情報がファイルに送られるかどうかを示しています。"-1"はモニタリング情報が生成されないことを意味します。
    Optimality phase itns. 最適化フェーズの最大反復数。
  • 20〜21行目には、与えられたワークスペースの大きさと問題を解くのに要したワークスペースの大きさが出力されています。
  • 23行目に目的関数行列のランクが出力されています。
  • 26〜39行目には以下が出力されています。
    Itn 反復回数。
    Step 探索方向に進んだステップ幅。
    Ninf 違反した制約の数。
    Sinf/Objective 制約違反の大きさの加重和(xが実行可能でない場合)もしくは目的関数の値(xが実行可能な場合)。
    Norm Gz 縮小勾配のユークリッド・ノルム。
  • 41行目にはLS問題を解くのに12回の反復が行われたことが示されています。
  • 44〜54行目には以下が出力されています。
    Varbl 変数を示す"V"、インデックス。
    State 変数の状態。LLは下限であることを意味し、FRは下限も上限もワーキングセットにはないことを意味します。
    Value 最後の反復での変数の値。
    Lower Bound 下限。
    Upper Bound 上限。
    Lagr Mult 関連する境界制約のラグランジュ乗数。
    Slack スラック(変数の値と近いほうの下限/上限との差)。
  • 57〜61行目には以下が出力されています。
    L Con 線形制約を示す"L"、インデックス。
    State 制約の状態。ULは上限であることを意味し、LLは下限であることを意味します。
    Value 最後の反復での制約の値。
    Lower Bound 下限。
    Upper Bound 上限。
    Lagr Mult 関連する境界制約のラグランジュ乗数。
    Slack スラック(変数の値と近いほうの下限/上限との差)。
  • 63行目にLSの最適解が見つかったことが示されています。
  • 65行目にLSの最終的な目的関数値が出力されています。

ソースコード

(本ルーチンの詳細はe04ncf のマニュアルページを参照)

※本サンプルソースコードは科学技術・統計計算ライブラリである「NAG Fortranライブラリ」のルーチンを呼び出します。
サンプルのコンパイル及び実行方法


このソースコードをダウンロード
    PROGRAM e04ncfe

!      E04NCF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : e04ncf, nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER              :: nin = 5, nout = 6
!      .. Local Scalars ..
       REAL (KIND=nag_wp)              :: obj
       INTEGER                         :: i, ifail, iter, lda, ldc, liwork,    &
                                          lwork, m, n, nclin, sdc
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE :: a(:,:), b(:), bl(:), bu(:), c(:,:),  &
                                          clamda(:), cvec(:), work(:), x(:)
       INTEGER, ALLOCATABLE            :: istate(:), iwork(:), kx(:)
!      .. Intrinsic Functions ..
       INTRINSIC                          max
!      .. Executable Statements ..
       WRITE (nout,*) 'E04NCF Example Program Results'
       FLUSH (nout)

!      Skip heading in data file
       READ (nin,*)

       READ (nin,*) m, n, nclin
       liwork = n
       ldc = max(1,nclin)
       lda = max(1,m)

       IF (nclin>0) THEN
          sdc = n
       ELSE
          sdc = 1
       END IF

!      This particular example problem is of type LS1, so we allocate
!      A(LDA,N), CVEC(1), B(M) and define LWORK as below

       IF (nclin>0) THEN
          lwork = 2*n**2 + 9*n + 6*nclin
       ELSE
          lwork = 9*n
       END IF

       ALLOCATE (istate(n+nclin),kx(n),iwork(liwork),c(ldc,sdc),bl(n+nclin), &
          bu(n+nclin),cvec(1),x(n),a(lda,n),b(m),clamda(n+nclin),work(lwork))

       READ (nin,*) (a(i,1:n),i=1,m)
       READ (nin,*) b(1:m)
       READ (nin,*) (c(i,1:sdc),i=1,nclin)
       READ (nin,*) bl(1:(n+nclin))
       READ (nin,*) bu(1:(n+nclin))
       READ (nin,*) x(1:n)

!      Solve the problem

       ifail = 0
       CALL e04ncf(m,n,nclin,ldc,lda,c,bl,bu,cvec,istate,kx,x,a,b,iter,obj, &
          clamda,iwork,liwork,work,lwork,ifail)

    END PROGRAM e04ncfe


Results matter. Trust NAG.

Privacy Policy | Trademarks