Dggglm Dgglse

Dgglse() Sub Dgglse ( M As  Long, N As  Long, P As  Long, A() As  Double, B() As  Double, C() As  Double, D() As  Double, X() As  Double, Info As  Long  ) 線形等式制約最小二乗(LSE)問題 目的 本ルーチンは次の線形等式制約最小二乗(LSE)問題を解く. B*x = d の制約条件のもとで || c - Ax ||_2 を最小化する. ここで, Aはm×n行列, Bはp×n行列, cは与えられたmベクトル, dは与えられたpベクトルである. ただし, p <= n <= m + p とする. また, 次式が成り立つものとする. rank(B) = p かつ rank( (A) ) = n ( (B) ) これらの条件により, LSE問題が一意の解を持つことが保証される. 解は次式で与えられる行列のペア(B, A)の一般化RQ分解を用いて得られる. B = (0 R)*Q, A = Z*T*Q 引数 [in] M 行列 A の行数. (M >= 0) [in] N 行列 A および B の列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in] P 行列 B の行数. (0 <= P <= N <= M + P) [in,out] A() 配列 A(LA1 - 1, LA2 - 1) (LA1 >= M, LA2 >= N) [in] M×N行列 A. [out] 配列の対角および上三角要素にmin(M, N)×N上台形行列Tが入る. [in,out] B() 配列 B(LB1 - 1, LB2 - 1) (LB1 >= P, LB2 >= N) [in] P×N行列 B. [out] 配列の一部 B(0〜P-1, N-P〜N-1)の上三角部分にP×P上三角行列Rが入る. [in,out] C() 配列 C(LC - 1) (LC >= M) [in] LSE問題の最小二乗部分の右辺ベクトル. [out] 解の残差二乗和が C(N-P)〜C(M-1) の二乗和で与えられる. [in,out] D() 配列 D(LD - 1) (LD >= P) [in] 制約方程式の右辺ベクトル. [out] D()は上書きされる. [out] X() 配列 X(LX - 1) (LX >= N) LSE問題の解. [out] Info = 0: 正常終了. = -1: パラメータ M の誤り. (M < 0) = -2: パラメータ N の誤り. (N < 0) = -3: パラメータ P の誤り. (P < 0 または P > N または P < N - M) = -4: パラメータ A() の誤り. = -5: パラメータ B() の誤り. = -6: パラメータ C() の誤り. = -7: パラメータ D() の誤り. = -8: パラメータ X() の誤り. = 1: (B, A)ペアの一般化RQ分解のBに関連する上三角行列Rが特異で, rank(B) < P である. 最小二乗解を求めることができなかった. = 2: (B, A)ペアの一般化RQ分解のAに関する上台形行列TのN-P×N-P部分が特異で, rank((A^T B^T)^T) < N である. 最小二乗解を求めることができなかった. 出典 LAPACK 使用例 線形等式制約最小二乗(LSE)問題を解く. すなわち, B*x = d の制約条件のもとで || c - Ax ||_2 を最小化する. ただし, ( -1.06 0.48 -0.04 ) A = ( -1.19 0.73 -0.24 ) ( 1.97 -0.89 0.56 ) ( 0.68 -0.53 0.08 ) ( -0.58 -0.79 0.82 ) B = ( 0.77 0.71 -0.55 ) ( -1.36 -1.22 1.66 ) ( 0.3884 ) c = ( 0.1120 ) ( -0.3644 ) ( -0.0002 ) ( 1.8250 ) d = ( -1.7058 ) ( 3.4904 ) とする. Sub Ex_Dgglse() Const M = 4, N = 3, P = 3 Dim A(M - 1, N - 1) As Double, B(P - 1, N - 1) As Double Dim C(M - 1) As Double, D(P - 1) As Double, X(N - 1) As Double Dim S As Double, Info As Long A(0, 0) = -1.06: A(0, 1) = 0.48: A(0, 2) = -0.04 A(1, 0) = -1.19: A(1, 1) = 0.73: A(1, 2) = -0.24 A(2, 0) = 1.97: A(2, 1) = -0.89: A(2, 2) = 0.56 A(3, 0) = 0.68: A(3, 1) = -0.53: A(3, 2) = 0.08 B(0, 0) = -0.58: B(0, 1) = -0.79: B(0, 2) = 0.82 B(1, 0) = 0.77: B(1, 1) = 0.71: B(1, 2) = -0.55 B(2, 0) = -1.36: B(2, 1) = -1.22: B(2, 2) = 1.66 C(0) = 0.3884: C(1) = 0.112: C(2) = -0.3644: C(3) = -0.0002 D(0) = 1.825: D(1) = -1.7058: D(2) = 3.4904 Call Dgglse(M, N, P, A(), B(), C(), D(), X(), Info) S = Dnrm2(M - N + P, C(N - P)) Debug.Print "X =", X(0), X(1), X(2) Debug.Print "SumSq =", S, "Info =", Info End Sub 実行結果 X = -0.820000000000002 -0.939999999999998 0.74 SumSq = 7.1082457242971E-15 Info = 0