Dgecov Dgecovs Dgecovy Dgels Dgelsd Dgelss Dgelsy Dgetsls

Dgelsd() Sub Dgelsd ( M As  Long, N As  Long, A() As  Double, B() As  Double, S() As  Double, RCond As  Double, Rank As  Long, Info As  Long, Optional Nrhs As  Long = 1  ) 優決定または劣決定系連立一次方程式 Ax = b の解 (特異値分解 (SVD)) (分割統治法) 目的 本ルーチンは線形最小二乗問題の最小ノルム解を A の特異値分解(SVD)を使って求める. 2-norm(| b - A*x |) を最小化する. A は M x N 行列で, ランク落ちしていてもよい. いくつかの右辺ベクトル b および解ベクトル x を 1 回の呼び出しで扱うことができる. これらのベクトルは, M x Nrhs 右辺行列 B および N x Nrhs 解行列 X の列として格納される. 問題を以下の 3 ステップにより解く. (1) ハウスホルダー変換により係数行列 A を二重対角形にする. これにより, 元の問題は「二重対角最小二乗問題」(BLS) に変換される. (2) 分割統治法により BLS を解く. (3) 元の最小二乗問題を解くためにハウスホルダー変換を逆に適用する. 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 x N 行列 A. [out] A() は壊される. [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 x Nrhs 右辺行列 B. [out] B() は N x 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 の最小二乗解を求める. また, 分散を求める. ただし, ( -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.3884 ) B = ( 0.1120 ) ( -0.3644 ) ( -0.0002 ) とする. Sub Ex_Dgelsd() Const M = 4, N = 3 Dim A(M - 1, N - 1) As Double, B(M - 1) As Double, Ci(N - 1) As Double Dim Sigma(N - 1) As Double, RCond As Double, Rank As Long, Info As Long Dim I 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.3884: B(1) = 0.112: B(2) = -0.3644: B(3) = -0.0002 RCond = 0.0001 Call Dgelsd(M, N, A(), B(), Sigma(), RCond, Rank, Info) If Info <> 0 Then Debug.Print "Error in Dgelss: Info =", Info Exit Sub End If Debug.Print "X =", B(0), B(1), B(2) Debug.Print "Rank =", Rank, "Info =", Info End Sub 実行結果 X = -0.82 -0.94 0.740000000000001 Rank = 3 Info = 0