◆ Zggglm()

Sub Zggglm	(	N As	Long,
		M As	Long,
		P As	Long,
		A() As	Complex,
		B() As	Complex,
		D() As	Complex,
		X() As	Complex,
		Y() As	Complex,
		Info As	Long
	)

一般Gauss-Markov線形モデル(GLM)問題 (複素行列)

目的: 本ルーチンは一般Gauss-Markov線形モデル(GLM)問題を解く.
d = A*x + B*y の制約条件のもとで || y ||_2 を最小化するxを求める.

ここで, AはN×M行列, BはN×P行列, dは与えられたNベクトルである. ただし, M <= N <= M + P とする. また, 次式が成り立つものとする.
rank(A) = M かつ rank(A B) = N

これらの条件により, 制約付き方程式は常に正しく, 一意の解xおよび最小2ノルム解yを持つ. 解は, 次式で与えられる行列のペア(A, B)の一般化QR分解を用いて得られる.
A = Q*(R), B = Q*T*Z

(0)

特に行列Bが正方で非特異ならば, GLM問題は次の重み付き線形最小二乗問題に等しい.
|| inv(B)*(d - A*x) ||_2 を最小化するxを求める.

引数

[in]	N	行列 A および B の行数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る)
[in]	M	行列 A の列数. (0 <= M <= N)
[in]	P	行列 B の列数. (P >= N - M)
[in,out]	A()	配列 A(LA1 - 1, LA2 - 1) (LA1 >= N, LA2 >= M) [in] N×M行列 A. [out] 配列A()の上三角要部分にM×M上三角行列Rが入る.
[in,out]	B()	配列 B(LB1 - 1, LB2 - 1) (LB1 >= N, LB2 >= P) [in] N×P行列 B. [out] N <= P であれば, 配列の一部 B(0〜N-1, P-N〜P-1)の上三角部分にN×N上三角行列Tが入る. N > P であれば, (N-P)番目の対角およびその上三角要素にN×P上台形行列Tが入る.
[in,out]	D()	配列 D(LD - 1) (LD >= N) [in] D()にGLM方程式の左辺ベクトルを入れる. [out] D()は壊される.
[out]	X()	配列 X(LX - 1) (LX >= M)
[out]	Y()	配列 Y(LY - 1) (LY >= P) X()およびY()はGLM問題の解を返す.
[out]	Info	= 0: 正常終了. = -1: パラメータ N の誤り. (N < 0) = -2: パラメータ M の誤り. (M < 0 または M > N) = -3: パラメータ P の誤り. (P < 0 または P < N - M) = -4: パラメータ A() の誤り. = -5: パラメータ B() の誤り. = -6: パラメータ D() の誤り. = -7: パラメータ X() の誤り. = -8: パラメータ Y() の誤り. = 1: (A, B)ペアの一般化QR分解のAに関連する上三角行列Rが特異で, rank(A) < Mである. 最小二乗解を求めることができなかった. = 2: (A, B)ペアの一般化QR分解のBに関する上台形行列Tの底辺側N-M×N-M部分が特異で, rank(A B) < Nである. 最小二乗解を求めることができなかった.

出典: LAPACK

使用例: 一般ガウス・マルコフ線形モデル(GLM)問題を解く. すなわち, d = A*x + B*y の制約条件のもとで || y ||_2 を最小化するxを求める. ただし,
( -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 )

( 1 0 0 0 )

B = ( 0 1 0 0 )

( 0 0 1 0 )

( 0 0 0 1 )

( 1.7126-0.6648i )

d = ( 0.8697+0.7604i )

( -2.1048-1.6171i )

( -0.9297+0.1252i )

とする.
Sub EX_Zggglm()

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

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

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

Dim X(M - 1) As Complex, Y(P - 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(1): B(0, 1) = Cmplx(0): B(0, 2) = Cmplx(0): B(0, 3) = Cmplx(0)

B(1, 0) = Cmplx(0): B(1, 1) = Cmplx(1): B(1, 2) = Cmplx(0): B(1, 3) = Cmplx(0)

B(2, 0) = Cmplx(0): B(2, 1) = Cmplx(0): B(2, 2) = Cmplx(1): B(2, 3) = Cmplx(0)

B(3, 0) = Cmplx(0): B(3, 1) = Cmplx(0): B(3, 2) = Cmplx(0): B(3, 3) = Cmplx(1)

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

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

Call Zggglm(N, M, P, A(), B(), D(), X(), Y(), Info)

S = Dznrm2(P, Y(0))

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.819999999999999 -0.939999999999999 0.740000000000001 0.200000000000002

0.480000000000002 0.21

SumSq = 1.06144666256319E-16 Info = 0