◆ Zgelss()

Sub Zgelss	(	M As	Long,
		N As	Long,
		A() As	Complex,
		B() As	Complex,
		S() As	Double,
		RCond As	Double,
		Rank As	Long,
		Info As	Long,
		Optional Nrhs As	Long = `1`
	)

優決定または劣決定系連立一次方程式 Ax = b の解 (特異値分解 (SVD)) (複素行列)

目的: 本ルーチンは複素線形最小二乗問題の最小ノルム解を, Aの特異値分解(SVD)を使って求める. すなわち,
|| A * X - B || を最小化する.

AはM×N行列でランク落ちしていてもよい.
いくつかの右辺ベクトル b および解ベクトル x を1回の呼び出しで扱うことができる. これらのベクトルは, M×Nrhs右辺行列BおよびN×Nrhs解行列Xの列として格納される.

Aの有効ランク数は, 最大の特異値のrcond倍よりも小さな特異値を0として扱うことにより決められる.

引数

[in]	M	行列 A の行数. (M >= 0) (M = 0 の場合, Rank = 0 として戻る)
[in]	N	行列 A の列数. (N >= 0) (N = 0 の場合, Rank = 0 として戻る)
[in,out]	A()	配列 A(LA1 - 1, LA2 - 1) (LA1 >= M, LA2 >= N) [in] M×N行列 A. [out] A()の最初のmin(M, N)行はその右特異ベクトル(行ごとに格納)により上書きされる.
[in,out]	B()	配列 B(LB1 - 1, LB2 - 1) (LB1 >= max(M, N), LB2 >= Nrhs) (2次元配列) または B(LB - 1) (LB >= max(M, N), Nrhs = 1) (1次元配列) [in] M×Nrhs右辺行列 B. [out] B()は N×Nrhs解行列 X により上書きされる. M >= N かつ Rank = N の場合, i列の解の残差二乗和は同じ列のN〜M-1番要素の二乗和で与えられる.
[out]	S()	配列 S(LS - 1) (LS >= min(M, N)) 行列Aの特異値 (降順). Aの2-ノルムによる条件数は S(0)/S(min(M, N)-1) である.
[in]	RCond	RCondはAの有効ランク数を決めるために使われる. S(i) <= RCond*S(0) となる特異値は0として扱われる. RCond < 0の場合, 代わりにマシンイプシロンが使われる.
[out]	Rank	Aの有効ランク数, すなわち, RCond*S(0)より大きい特異値の数.
[out]	Info	= 0: 正常終了. = -1: パラメータ M の誤り. (M < 0) = -2: パラメータ N の誤り. (N < 0) = -3: パラメータ A() の誤り. = -4: パラメータ B() の誤り. = -5: パラメータ S() の誤り. = -9: パラメータ Nrhs の誤り. (Nrhs < 0) = i > 0: SVDの計算アルゴリズムが収束しなかった. 中間結果の二重対角形の副対角要素のうちi個が0に収束しなかった.
[in]	Nrhs	(省略可) 右辺の数, すなわち, 行列BおよびXの列数. (Nrhs >= 0) (Nrhs = 0 の場合, Rank = 0 として戻る) (省略時 = 1)

出典: LAPACK

使用例: 優決定系連立1次方程式 Ax = B の最小二乗解を求める. また, 分散を求める. ただし,
( -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.7126-0.6648i )

B = ( 0.8697+0.7604i )

( -2.1048-1.6171i )

( -0.9297+0.1252i )

とする.
Sub Ex_Zgelss()

Const M = 4, N = 3

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

Dim Sigma(N - 1) As Double, RCond As Double, Rank As Long, Info As Long, I 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) = Cmplx(1.7126, -0.6648): B(1) = Cmplx(0.8697, 0.7604)

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

RCond = 0.0001

Call Zgelss(M, N, A(), B(), Sigma(), RCond, Rank, Info)

If Info <> 0 Then

Debug.Print "Error in Zgelss: Info =", Info

Exit Sub

End If

Debug.Print "X ="

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

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

Call Zgecovs(0, N, A(), Sigma(), Ci(), Info)

Debug.Print "Var ="

Debug.Print Creal(Ci(0)), Creal(Ci(1)), Creal(Ci(2))

Debug.Print "Rank =", Rank, "Info =", Info

End Sub

実行結果: X =

-0.82 -0.94 0.740000000000001 0.200000000000003

0.480000000000002 0.21

Var =

5.46169501938982 15.1880504061464 21.4241290120714

Rank = 3 Info = 0