C#による 線形最小二乗問題

C#によるサンプルソースコード : 使用関数名:e04nc

Keyword: 線形最小二乗問題

概要

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

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

※本サンプルはnAG Library for .NETに含まれる関数 e04nc() のExampleコードです。本サンプル及び関数の詳細情報は e04nc のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで

入力データ

(本関数の詳細はe04nc のマニュアルページを参照)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

このデータをダウンロード
e04nc 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)、変数の数(n)、線形制約の数(nclin)を指定しています。
  • 3~12行目には目的関数行列Aの要素を指定しています。
  • 13行目には目的関数の係数ベクトル(b)を指定しています。
  • 14~16行目には線形制約行列Cの要素を指定しています。
  • 17行目には変数の下限(bl(1,..,9))を指定しています。
  • 18行目には制約の下限(bl(10,11,12))を指定しています。
  • 19行目には変数の上限(bu(1,..,9))を指定しています。
  • 20行目には制約の上限(bu(10,11,12))を指定しています。
  • 21行目にはxの初期推定値を指定しています。

出力結果

(本関数の詳細はe04nc のマニュアルページを参照)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93

この出力例をダウンロード
e04nc Example Program Results
      Print level = 10

 *** 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.0E-01    0   5.342505E+00  1.9E+01
   4  7.1E-02    0   4.616975E+00  2.2E+00
   5  1.0E-01    0   4.558492E+00  1.3E+00
   6  1.0E+00    0   4.523485E+00  9.8E-16
   7  3.5E-01    0   1.934106E+00  6.9E+00
   8  2.1E-01    0   1.323283E+00  5.1E+00
   9  1.4E-02    0   1.307479E+00  0.0E+00
  10  1.0E+00    0   9.153991E-01  5.3E-15
  11  6.1E-01    0   2.190278E-01  5.9E-01
  12  1.0E+00    0   1.652065E-01  2.2E-15
  13  1.0E+00    0   9.605160E-02  2.2E-15
  14  3.0E-02    0   9.236999E-02  4.5E-01
  15  1.0E+00    0   8.134082E-02  8.3E-16

 Exit from LS problem after    15 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     -4.4409E-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      4.4409E-16

 Exit E04NCF - Optimal LS solution.

 Final LS objective value =   0.8134082E-01

 ** e04nc returned with ifail =   0

  Varbl    Istate      Value                  Lagr Mult

  V    1    1                     0          0.1572
  V    2    0             0.0415261               0
  V    3    0              0.587176               0
  V    4    1                     0          0.8782
  V    5    0             0.0996432               0
  V    6    1                     0          0.1473
  V    7    0             0.0490578               0
  V    8    1                     0          0.8603
  V    9    0              0.305649               0


  L Con    Istate      Value                  Lagr Mult

  L    1    1                     2          0.3777
  L    2    2                     2        -0.05791
  L    3    1                     1          0.1075


  Final objective value =      0.08134082

  • 9~22行目に以下に示すプログラムのオプション引数が出力されています。
    Problem type 最小化される目的関数の種類。"LS1"は標準の線形最小二乗(LS)問題を意味します。
    Hessian 上三角行列 R の内容を制御します。"NO"は変換され並べ替えられた行列HQのコレスキー因子を含むことを意味しています。
    Linear constraints 線形制約の数。
    Feasibility tolerance 実行可能解での許容可能な制約違反の最大値。
    Variables 変数の数。
    Crash tolerance 初期のワーキングセットに含まれる不等式制約が下限/上限のこの範囲内にあることを示しています。
    Objective matirx rows 目的関数行列の行数。
    Rank tolerance 三角因子の推定値を制御します。
    Infinite bound size 無限の境界値。この値以上の上限は+∞と見なされます。またこの値を負に反転させた値以下の下限は−∞と見なされます。
    COLD start コールドスタート。変数や制約の値に基づいて初期のワーキングセットが選ばれていることを意味します。
    Infinite step size 無限解へのステップと見なされる変数の変化の大きさ。
    EPS (machine precision) マシンの精度。
    Print level 結果出力のレベル。"10"は最終解と各反復の結果を出力することを意味します。
    Feasibility phase itns. 実行可能フェーズの最大反復数。
    Monitoring file モニタリング情報の生成を制御します。"-1"はモニタリング情報が生成されないことを意味します。
    Optimality phase itns. 最適化フェーズの最大反復数 。
    Workspace provided is 用意されているワークスペース。
    To solve problem we need 問題を解くのに必要なワークスペース。
  • 24行目に目的関数行列のランクが出力されています。
  • 27~43行目には以下が出力されています。
    Itn 反復回数。
    Step 探索方向に進んだステップ幅。
    Ninf 違反した制約の数。
    Sinf/Objective 制約違反の大きさの加重和(xが実行可能でない場合)もしくは目的関数の値(xが実行可能な場合)。
    Norm Gz 縮小勾配のユークリッド・ノルム。
  • 45行目にはLS問題を解くのに15回の反復が行われたことが示されています。
  • 48~58行目には以下が出力されています。
    Varbl 変数を示す "V"、インデックス。
    State 変数の状態。"LL" は下限であることを意味し、"FR" は下限も上限もワーキングセットにはないことを意味します。
    Value 最後の反復での変数の値。
    Lower Bound 変数の下限。
    Upper Bound 変数の上限。
    Lagr Mult 関連する境界制約のラグランジュ乗数。
    Slack スラック(変数の値と近いほうの下限/上限との差)。
  • 61~65行目には以下が出力されています。
    L Con 線形制約を示す "L"、インデックス。
    State 制約の状態。"LL" は下限であることを意味し、"UL" は上限であることを意味します。
    Value 最後の反復での制約の値。
    Lower Bound 制約の下限。
    Upper Bound 制約の上限。
    Lagr Mult 関連する境界制約のラグランジュ乗数。
    Slack スラック(制約の値と近いほうの下限/上限との差)。
  • 69行目に線形最小二乗(LS)の最終的な目的関数値が出力されています。
  • 本関数 e04nc がエラーを検知せずに終了したことが示されています。
  • 75~83行目には以下が出力されています。
    Varbl 変数を示す "V"、インデックス。
    Istate 変数の状態。"1" は不等式制約の下限がワーキングセットに含まれることを意味し、"0" はFeasibility toleranceの範囲内にあるがワーキングセットにはないことを意味します。
    Value 最後の反復での変数の値。
    Lagr Mult 関連する境界制約のラグランジュ乗数。
  • 86~90行目には以下が出力されています。
    L Con 線形制約を示す "L"、インデックス。
    Istate 制約の状態。"1" は不等式制約の下限がワーキングセットに含まれることを意味し、"2" は不等式制約の上限がワーキングセットに含まれることことを意味します。
    Value 最後の反復での制約の値。
    Lagr Mult 関連する境界制約のラグランジュ乗数。
  • 93行目に最終的な目的関数値が出力されています。

ソースコード

(本関数の詳細はe04nc のマニュアルページを参照)

※本サンプルソースコードは .NET環境用の科学技術・統計計算ライブラリである「nAG Library for .NET」の関数を呼び出します。
サンプルのコンパイル及び実行方法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162

このソースコードをダウンロード
//      e04nc Example Program Text
//      C# version, nAG Copyright 2008
using System;
using NagLibrary;
using System.IO;
namespace NagDotNetExamples
{
  public class E04NCE
  {
    static bool defaultdata = true;
    static string datafile = "";
    // as a command line argument. It defaults to the named file specified below otherwise.
    //
    static void Main(String[] args)
    {
      if (args.Length == 1)
      {
        defaultdata = false;
        datafile = args[0];
      }
      StartExample();
    }
    public static void StartExample()
    {
      try
      {
        DataReader sr = null;
        if(defaultdata)
        {
          sr = new DataReader("exampledata/e04nce.d");
        }
        else
        {
          sr = new DataReader(datafile);
        }
        double obj = 0.0;
        int i,  ifail,  iter,  j,  m,  n,  nclin;
        Console.WriteLine("e04nc Example Program Results");
        //      Skip heading in data file
        sr.Reset();
        sr.Reset();
        m = int.Parse(sr.Next());
        n = int.Parse(sr.Next());
        nclin = int.Parse(sr.Next());
        //
        //         Read a, b, c, bl, bu and x from data file
        //
        double[,] a = new double[Math.Max(1,m), n];
        double[] b = new double[m];
        double[] bl = new double[n + nclin];
        double[] bu = new double[n + nclin];
        double[,] c = new double[nclin, n];
        double[] clamda = new double[n + nclin];
        double[] cvec = new double[n];
        double[] work = new double[nclin];
        double[] x = new double[n];
        int[] istate = new int[n + nclin];
        int[] kx = new int[n];
        sr.Reset();
        for (i = 1; i <= m; i++)
        {
          for (j = 1; j <= n; j++)
          {
            a[i - 1, j - 1] = double.Parse(sr.Next());
          }
        }
        sr.Reset();
        for (i = 1; i <= m; i++)
        {
          b[i - 1] = double.Parse(sr.Next());
        }
        sr.Reset();
        for (i = 1; i <= nclin; i++)
        {
          for (j = 1; j <= n; j++)
          {
            c[i - 1, j - 1] = double.Parse(sr.Next());
          }
        }
        sr.Reset();
        for (i = 1; i <= n + nclin; i++)
        {
          bl[i - 1] = double.Parse(sr.Next());
        }
        sr.Reset();
        for (i = 1; i <= n + nclin; i++)
        {
          bu[i - 1] = double.Parse(sr.Next());
        }
        sr.Reset();
        for (i = 1; i <= n; i++)
        {
          x[i - 1] = double.Parse(sr.Next());
        }
        //
        //         Initialise E04.e04nc and check for error exits
        //
        E04.e04ncOptions options = new E04.e04ncOptions();
        //
        //            Solve the problem
        //
        options.Set("Defaults");
        options.Set("List");
        options.Set("Print level = 10");
        //
        E04.e04nc(m, n, nclin, c, bl, bu, cvec, istate, kx, x, a, b, out iter, out obj, clamda,
        options, out ifail);
        //
        //            A valid licence was found: continue.
        //            Check for error exits
        //
        Console.WriteLine("");
        if (ifail == 6)
        {
          Console.WriteLine("   ** An input parameter is invalid");
        }
        else if (ifail >= 0)
        {
          Console.WriteLine(" ** e04nc returned with ifail = {0, 3}", ifail);
          Console.WriteLine("");
          Console.WriteLine("  Varbl    Istate      Value                  Lagr Mult");
          Console.WriteLine("");
          for (i = 1; i <= n; i++)
          {
            Console.WriteLine("  V  {0,3}  {1,3}        {2,14:g6}    {3,12:g4}", i, istate[i - 1], x[i - 1], clamda[i - 1]);
          }
          if (nclin > 0)
          {
            //
            //  This performs the matrix vector multiplication C*X
            // (linear constraint values) and puts the result in
            // the first nclin locations of work.
            //
            F06.f06pa("N", nclin, n, 1.00e0, c, x, 1, 0.00e0, work, 1, out ifail);
            Console.WriteLine("");
            Console.WriteLine("");
            Console.WriteLine("  L Con    Istate      Value                  Lagr Mult");
            Console.WriteLine("");
            for (i = n + 1; i <= n + nclin; i++)
            {
              j = i - n;
              Console.WriteLine("  L  {0,3}  {1,3}        {2,14:g6}    {3,12:g4}", j, istate[i - 1], work[j - 1], clamda[i - 1]);
            }
          }
          Console.WriteLine("");
          Console.WriteLine("");
          Console.WriteLine("  Final objective value = {0,15:g7}", obj);
        }
        else
        {
          Console.WriteLine(" ** e04nc returned with ifail = {0, 3}", ifail);
        }
        //
      }
      catch (Exception e)
      {
        Console.WriteLine(e.Message);
        Console.WriteLine( "Exception Raised");
      }
    }
  }
}


関連情報
© 日本ニューメリカルアルゴリズムズグループ株式会社 2025
Privacy Policy  /  Trademarks