Zgejsv Zgesdd Zgesvd Zgesvdq Zgesvdx Zggsvd3

Zgejsv() Sub Zgejsv ( Joba As  String, Jobu As  String, Jobv As  String, Jobr As  String, Jobt As  String, Jobp As  String, M As  Long, N As  Long, A() As  Complex, S() As  Double, U() As  Complex, V() As  Complex, Info As  Long, Optional Scale0 As  Double, Optional Scale1 As  Double, Optional Ns As  Long, Optional SCond As  Double, Optional Nrank As  Long  ) 特異値分解 (SVD) (前処理付きヤコビ SVDアルゴリズム) (複素行列) 目的 本ルーチンはm >= nであるm×n複素行列Aの特異値分解(SVD)を求める. AのSVDは次のように表される. A = U * SIGMA * V^H, ここで, SIGMAはn個の対角成分を除き0のn×n(m×n)行列である. Uはm×n(またはm×m)ユニタリ行列, Vはn×nユニタリ行列である. SIGMAの対角要素はAの特異値である. UおよびVの列はそれぞれAの左および右特異ベクトルである. 求められた行列UおよびVはそれぞれ配列U()およびV()に格納される. 求められたSIGMAの対角要素は配列S()に格納される. 引数 [in] Joba 精度のレベルを指定する. = "C": A = B * D と表せる場合, 本オプションが有効(相対的に高精度)である. ここで, Bは条件数の良い行列, Dは任意の対角行列である. 列のスケーリングにより精度が落ちることはない. 結果の精度はBの条件数に依存し, 本ルーチンは理論上の最高精度を目指す. 相対誤差 max{i = 1〜N}|d σi|/σi の上限は f(M, N)*epsilon*cond(B) となり, Dに依存しない.   入力行列は, 列のピボット選択付きのQR分解により前処理される. Rank Revealing QR分解による初期前処理およびプリコンディショニングはJobaのすべての値に共通である. 以下では追加動作を指定する. = "E": "C"と同様に計算するが, Bの条件数の推定値も求める. これは現実的な誤差限界を与える. = "F": A = D1 * C * D2 と表せる場合, 本オプションは"C"よりも高精度である. ここで, D1およびD2は条件数の悪い対角スケーリング, Cは条件数の良い行列である. 入力行列の構造が不明であるが高精度を望む場合に本オプションが推奨される. 入力行列Aはフル(行および列の)ピボット選択付きのQR分解により前処理される. = "G": "F"と同様に計算するが, Bの条件数も推定する. ただし, A = D*B である. Aの行が強く重みづけされている場合, この条件数を使用すると安全側すぎる誤差限界を与える. = "A": 小さな特異値はノイズとされ, 行列は数値的にランク落ちとして扱われる. 特異値の計算誤差限界は f(M, N)*epsilon*||A|| である. 求められたSVD A = U*S*V^H は f(M, N)*epsilon*||A|| まで復元できる. これにより, N*epsilon*||A|| より小さい特異値を捨てる(0にする)ことができるようにする. = "R": "A"と同様であるが, 最初のQR分解のRank Revealingプロパティを利用して σr+1 < ε*σr となるrを(三角行列を使い)見つける. rはランク数とされる. SVDの計算は絶対誤差限界を用いて行われるが, "A"の場合より精度が高い. [in] Jobu Uの列を計算するかどうか指定する. = "U": UのN列を配列U()に返す. = "F": M個の左特異ベクトル全部を配列U()に返す. = "W": 配列U()をM*Nの作業領域として使用する. U()の説明を参照せよ. = "N": Uを計算しない. [in] Jobv 行列Vを計算するかどうか指定する. = "V": VのN列を配列V()に返す. ヤコビ回転は明示的には累算されない. = "J": Jobt = "N"の場合, Vのn列を配列V()に返すが, これらはヤコビ回転の積として求められる. = "W": 配列V()をN*Nの作業領域として使用する. V()の説明を参照せよ. = "N": Vを計算しない. [in] Jobr 特異値の範囲を指定する. 指定された範囲外であれば小さな正の特異値を0に設定することができるようにする. c*Aの最大の特異値がおおよそsqrt(big)になるようにA(<> 0)がスケーリングされていれば, c*Aのノルムがsqrt(sfmin)以下(Jobr = "R"の場合)あるいはsmall = sfmin/epsln以下であるAの列をJobrにより無効にすることができる. ただし, big = Dlamch("O"), sfmin = Dlamch("S"), epsln = Dlamch("E")である. = "N": c*Aの小さな列を無効にしない. このオプションは, BLAS, QR分解および三角行列の解を求めるルーチンがこの範囲で動作するように実装されていることを前提とする. Aの条件数がbigよりも大きければZgesvjを使用せよ. = "R": sigma(c*A)の制限された範囲は [sqrt(sfmin), sqrt(big)] である(大ざっぱには上の説明のとおり). このオプションが推奨される. フル範囲 [sfmin, big] で特異値を計算するためには, Zgesvjを使用せよ. [in] Jobt 入力が正方行列のとき, Aの転置のほうが収束に関して良ければA^Hの使用を選択するかもしれない. 正方行列でない場合, Jobtは無視される. これは, A^H*Aの随伴軌道上の2種のエントロピー値に基づいて決定される. = "T": エントロピーテストによりA^Hを入力として使った方がヤコビプロセスの収束が速い可能性が示されれば転置を行う. = "N": 転置を考えない. オプション"T"は, 特異値だけを求める場合, または, SVD全部(U, SIGMA, V)を求める場合にのみ使用できる. 片方の特異ベクトル(UあるいはV)だけの場合, U()とV()の両方を与えて呼び出さなければならない. これは, Aが転置された場合に一方の配列は作業領域として使われるからである. U()およびV()の説明を参照せよ. [in] Jobp 非正規化数を捨てるために構造化された摂動を加えるかどうか指定する. 非正規化数が不完全に実装されている場合, 計算が遅くなってしまう場合, 特に収束が速い場合, にこれを使用すべきである. 簡単のため, SVD全部または特異値だけが必要なときにのみ摂動が加えられる. 実装者/ユーザーは, 1セットの特異ベクトルを計算する場合には, 容易に摂動を加えることができる. = "P": 摂動を加える. = "N": 摂動を加えない. [in] M 入力行列Aの行数. (M >= 0) (M = 0 の場合, 処理を行わずに戻る) [in] N 入力行列Aの列数. (M >= N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in,out] A() 配列 A(LA1 - 1, LA2 - 1) (LA1 >= M, LA2 >= N) [in] M×N 行列 A. [out] A()の内容は壊される. [out] S() 配列 S(LS - 1) (LS >= N) Scale0/Scale1 = 1: Aの特異値. 計算中は配列A()の中の反復行列の列の2-ノルムがS()に入る. Scale0 <> Scale1: Aの特異値は (Scale0/Scale1)*S(0〜N-1) である. この分解形は, σmax(A)がオーバーフローするか, 入力行列Aをスケーリングすることにより小さな特異値がアンダーフローから救われたときに使われる. Jobr = "R"の場合, いくつかの特異値が0と返されるかもしれない. これは, ランク判定値より小さいか非正規化数であるために「0に設定」されたことによるものである. [out] U() 配列 U(LU1 - 1, LU2 - 1) (LU1 >= M, LU2 >= M (Jobu = "F") または N (その他)) Jobu = "U": U()に左特異ベクトルのM×N行列が入る. Jobu = "F": U()に左特異ベクトルのM×M行列が入る. Aの値域の直交補空間の正規直交基底を含む. Jobu = "W": Jobv = "V"かつJobt = "T"かつM = Nの場合, Aの代わりにA^Hが使われたときにU()は作業領域として使用される. この場合, VはA^Hの左特異ベクトルとしてU()の中で計算された後, 配列V()にコピーされる. この"W"オプションは, U()が長さN*Nの作業領域として予約されていることを単に呼び出し元に示しているだけである. Jobu = "N": Jobt = "T"でない限りU()は参照されない. [out] V() 配列 V(LV1 - 1, LV2 - 1) (LV1 >= N, LV2 >= N) Jobv = "V"または"J": V()に右特異ベクトルのN×N行列が入る. Jobv = "W": Jobu = "U"かつJobt = "T"かつM = Nの場合, Aの代わりにA^Hが使われたときにV()は作業領域として使用される. この場合, UはA^Hの右特異ベクトルとしてV()の中で計算された後, 配列U()にコピーされる. この"W"オプションは, V()が長さN*Nの作業領域として予約されていることを単に呼び出し元に示しているだけである. Jobv = "N": Jobt = "T"でない限りV()は参照されない. [out] Info = 0: 正常終了. = -1: パラメータ Jobu の誤り. (Joba <> "C", "E", "F", "G", "A" および "R") = -2: パラメータ Jobu の誤り. (Jobu <> "U", "F", "W" および "N") = -3: パラメータ Jobv の誤り. (Jobv <> "V", "J", "W" および "N") = -4: パラメータ Jobr の誤り. (Jobr <> "N" および "R") = -5: パラメータ Jobt の誤り. (Jobt <> "T" および "N") = -6: パラメータ Jobp の誤り. (Jobp <> "P" および "N") = -7: パラメータ M の誤り. (M < 0) = -8: パラメータ N の誤り. (N < 0 または N > M) = -9: パラメータ A() の誤り. = -10: パラメータ S() の誤り. = -11: パラメータ U() の誤り. = -12: パラメータ V() の誤り. = 1: スケーリング係数が1でない(Scale0 <> Scale1). (警告) = 2: Aの列ノルムのいくつかが非正規化数であった. このデータでは要求精度が保証されない. (警告) = 29: 最大スイープ回数(30)で収束しなかった. [out] Scale0 (Optional) [out] Scale1 (Optional) scale = Scale1/Scale0 はスケーリング係数で, 求められたAの特異値はscale*S(0〜N-1)である(S()の説明を参照せよ). [out] Ns (省略可) 求められた0でない特異値の数. [out] SCond (省略可) 列の均衡化後のAの条件数の推定値. (Joba = "E"または"G"の場合) Scondaはsqrt(||(R^H*R)^(-1)||_1) の推定値である. これはDpoconを使って計算され, N^(-1/4)*Sconda <= ||R^(-1)||_2 <= N^(1/4)*Sconda となる. ここで, RはAのQR分解で得られる三角行列である. しかし, Rが不完全で, ランク数がNより小さいと判定された場合, Sconda = -1 が返される. これは, 最小の特異値が失われた可能性があることを示す. [out] Nrank (省略可) 最初のピボット選択付きQR分解後に判定されたランク数. JobaおよびJobrの説明を参照せよ. 出典 LAPACK 使用例 行列Aの特異値, 左および右特異ベクトルを求める. ただし, ( 0.20-0.11i -0.93-0.32i 0.81+0.37i ) A = ( -0.80-0.92i -0.29+0.86i 0.64+0.51i ) ( 0.71+0.59i -0.15+0.19i 0.20+0.94i ) とする. Sub Ex_Zgejsv() Const M = 3, N = 3 Dim A(M - 1, N - 1) As Complex, U(M - 1, N - 1) As Complex, V(N - 1, N - 1) As Complex Dim S(N - 1) As Double, Ns As Long, Nrank As Long, Info As Long A(0, 0) = Cmplx(0.2, -0.11): A(0, 1) = Cmplx(-0.93, -0.32): A(0, 2) = Cmplx(0.81, 0.37) A(1, 0) = Cmplx(-0.8, -0.92): A(1, 1) = Cmplx(-0.29, 0.86): A(1, 2) = Cmplx(0.64, 0.51) A(2, 0) = Cmplx(0.71, 0.59): A(2, 1) = Cmplx(-0.15, 0.19): A(2, 2) = Cmplx(0.2, 0.94) Call Zgejsv("C", "U", "V", "N", "N", "N", M, N, A(), S(), U(), V(), Info, , , Ns, , Nrank) Debug.Print "Singular values =", S(0), S(1), S(2) Debug.Print "U =" Debug.Print Creal(U(0, 0)), Cimag(U(0, 0)), Creal(U(0, 1)), Cimag(U(0, 1)) Debug.Print Creal(U(1, 0)), Cimag(U(1, 0)), Creal(U(1, 1)), Cimag(U(1, 1)) Debug.Print Creal(U(2, 0)), Cimag(U(2, 0)), Creal(U(2, 1)), Cimag(U(2, 1)) Debug.Print Creal(U(0, 2)), Cimag(U(0, 2)) Debug.Print Creal(U(1, 2)), Cimag(U(1, 2)) Debug.Print Creal(U(2, 2)), Cimag(U(2, 2)) Debug.Print "V^H =" Debug.Print Creal(V(0, 0)), Cimag(V(0, 0)), Creal(V(1, 0)), Cimag(V(1, 0)) Debug.Print Creal(V(0, 1)), Cimag(V(0, 1)), Creal(V(1, 1)), Cimag(V(1, 1)) Debug.Print Creal(V(0, 2)), Cimag(V(0, 2)), Creal(V(1, 2)), Cimag(V(1, 2)) Debug.Print Creal(V(2, 0)), Cimag(V(2, 0)) Debug.Print Creal(V(2, 1)), Cimag(V(2, 1)) Debug.Print Creal(V(2, 2)), Cimag(V(2, 2)) Debug.Print "Ns =", Ns, "Nrank =", Nrank, "Info =", Info End Sub 実行結果 Singular values = 2.07084030821889 1.23513084760163 0.901483337149816 U = 0.322218195271111 -0.402632553976182 0.323580768676013 0.032766875846889 -0.255509191701108 -0.690011678142581 -0.219080816933185 0.593565694755161 0.438566422474147 -1.79134356330225E-02 0.612649925763086 0.34433304824613 -0.583412100357814 -0.536576742191836 0.241341326648261 6.05497287146372E-03 0.108590300074888 0.549219053729666 V^H = 0.603023103819567 1.86227915744298E-04 -0.366673081969094 0.394409352108013 0.265910254090572 -0.608232432266715 0.191162102396799 -6.37222693955672E-02 0.160663801530462 0.412183562948795 0.81816174251754 -6.93889390390723E-18 -0.160579494721751 -0.566188521062472 0.714850576679605 -8.79216760506474E-02 2.71033959614138E-02 -0.366290352480977 Ns = 3 Nrank = 3 Info = 0