|
|
◆ 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は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 "Var ="
Debug.Print "Rank =", Rank, "Info =", Info
End Sub
Function Cmplx(R As Double, Optional I As Double=0) As Complex 複素数の作成 Function Cimag(A As Complex) As Double 複素数の虚数部 Function Creal(A As Complex) As Double 複素数の実数部 Function Ci(X As Double, Optional Info As Long) As Double 余弦積分 Ci(x) 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)) (複素行列) Sub Zgecovs(Job As Long, N As Long, A() As Complex, S() As Double, Ci() As Complex, Info As Long) 線形最小二乗問題の分散・共分散行列 (zgelss用)
- 実行結果
X =
-0.82 -0.94 0.740000000000001 0.200000000000003
0.480000000000002 0.21
Var =
5.46169501938982 15.1880504061464 21.4241290120714
Rank = 3 Info = 0
|
1.9.7