◆ Zgglse()

Sub Zgglse	(	M As	Long,
		N As	Long,
		P As	Long,
		A() As	Complex,
		B() As	Complex,
		C() As	Complex,
		D() As	Complex,
		X() As	Complex,
		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 を最小化する. ただし,
( -0.82+0.83i 0.18-0.94i -0.18-0.12i )

A = ( -0.76-0.24i 0.57-0.16i -0.08-0.27i )

( 1.90+0.26i -0.98+0.54i 0.21+0.28i )

( 0.50-0.30i -0.31+0.37i 0.22+0.19i )

( 0.57-0.91i -0.28-0.45i 0.25+0.91i )

B = ( 0.83+0.46i 0.63-0.19i -0.69+0.09i )

( 0.24-1.33i -0.56-0.67i 0.90+1.25i )

( 1.7126-0.6648i )

c = ( 0.8697+0.7604i )

( -2.1048-1.6171i )

( -0.9297+0.1252i )

( -1.5111+0.3107i )

d = ( -0.0941-1.2737i )

( -1.5579+1.0462i )

とする.
Sub Ex_Zgglse()

Const M = 4, N = 3, P = 3

Dim A(M - 1, N - 1) As Complex, B(P - 1, N - 1) As Complex

Dim C(M - 1) As Complex, D(P - 1) As Complex, X(N - 1) As Complex

Dim S As Double, Info As Long

A(0, 0) = Cmplx(-0.82, 0.83): A(0, 1) = Cmplx(0.18, -0.94): A(0, 2) = Cmplx(-0.18, -0.12)

A(1, 0) = Cmplx(-0.76, -0.24): A(1, 1) = Cmplx(0.57, -0.16): A(1, 2) = Cmplx(-0.08, -0.27)

A(2, 0) = Cmplx(1.9, 0.26): A(2, 1) = Cmplx(-0.98, 0.54): A(2, 2) = Cmplx(0.21, 0.28)

A(3, 0) = Cmplx(0.5, -0.3): A(3, 1) = Cmplx(-0.31, 0.37): A(3, 2) = Cmplx(0.22, 0.19)

B(0, 0) = Cmplx(0.57, -0.91): B(0, 1) = Cmplx(-0.28, -0.45): B(0, 2) = Cmplx(0.25, 0.91)

B(1, 0) = Cmplx(0.83, 0.46): B(1, 1) = Cmplx(0.63, -0.19): B(1, 2) = Cmplx(-0.69, 0.09)

B(2, 0) = Cmplx(0.24, -1.33): B(2, 1) = Cmplx(-0.56, -0.67): B(2, 2) = Cmplx(0.9, 1.25)

C(0) = Cmplx(1.7126, -0.6648): C(1) = Cmplx(0.8697, 0.7604)

C(2) = Cmplx(-2.1048, -1.6171): C(3) = Cmplx(-0.9297, 0.1252)

D(0) = Cmplx(-1.5111, 0.3107): D(1) = Cmplx(-0.0941, -1.2737)

D(2) = Cmplx(-1.5579, 1.0462)

Call Zgglse(M, N, P, A(), B(), C(), D(), X(), Info)

S = Dznrm2(M - N + P, C(N - P))

Debug.Print "X ="

Debug.Print Creal(X(0)), Cimag(X(0)), Creal(X(1)), Cimag(X(1))

Debug.Print Creal(X(2)), Cimag(X(2))

Debug.Print "SumSq =", S, "Info =", Info

End Sub

実行結果: X =

-0.82 -0.940000000000001 0.739999999999999 0.2

0.48 0.209999999999999

SumSq = 5.2372425074627E-15 Info = 0