_dgejsv _dgesdd _dgesvd _dgesvdq _dgesvdx _dggsvd3

_dgejsv() void _dgejsv ( char  joba, char  jobu, char  jobv, char  jobr, char  jobt, char  jobp, int  m, int  n, int  lda, double  a[], double  sva[], int  ldu, double  u[], int  ldv, double  v[], double  work[], int  lwork, int  iwork[], int *  info  ) 特異値分解 (SVD) (前処理付きヤコビ SVDアルゴリズム) 目的 本ルーチンはm >= nであるm×n実行列Aの特異値分解(SVD)を求める. AのSVDは次のように表される. A = U * SIGMA * V^T, ここで, 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の対角要素は配列sva[]に格納される. dgejsvは, 非常に小さな特異値およびその特異ベクトルを他のSVDルーチンよりかなり精度良く求めることができる場合がある. 引数 [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^Tは 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': Vのn列を配列v[][]に返すが, これらはヤコビ回転の積として求められる. 本オプションはjobu != 'N'の場合, すなわち, フルSVDを計算する場合にのみ許される. = '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よりも大きければdgesvjを使用せよ. = 'R': sigma(c*A)の制限された範囲は [sqrt(sfmin), sqrt(big)] である(大ざっぱには上の説明のとおり). このオプションが推奨される. フル範囲 [sfmin, big] で特異値を計算するためには, dgesvjを使用せよ. [in] jobt 入力が正方行列のとき, Aの転置のほうが収束に関して良ければA^Tの使用を選択するかもしれない. 正方行列でない場合, jobtは無視される. これは, A^T*Aの随伴軌道上の2種のエントロピー値に基づいて決定される. work[5]およびwork[6]の説明を参照せよ. = 'T': エントロピーテストによりA^Tを入力として使った方がヤコビプロセスの収束が速い可能性が示されれば転置を行う. = '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] lda 二次元配列a[][]の整合寸法. (lda >= max(1, m)) [in,out] a[][] 配列 a[la][lda] (la >= n) [in] m×n 行列 A. [out] a[][]の内容は壊される. [out] sva[] 配列 sva[lsva] (lsva >= n) work[0]/work[1] = 1: Aの特異値. 計算中は配列a[][]の中の反復行列の列の2-ノルムがsva[]に入る. work[0] != work[1]: Aの特異値は (work[0]/work[1])*sva[0〜n-1] である. この分解形式は, σmax(A)がオーバーフローするか, 入力行列Aをスケーリングすることにより小さな特異値がアンダーフローから救われたときに使われる. jobr = 'R'の場合, いくつかの特異値が0と返されるかもしれない. これは, ランク判定値より小さいか非正規化数であるために「0に設定」されたことによるものである. [in] ldu 二次元配列u[][]の整合寸法. (ldu >= 1 (jobu = 'N'の場合), ldu >= m (jobu = 'U', 'F'または'W'の場合)) [out] u[][] 配列 u[lu][ldu] (lu >= m (jobu = 'F'の場合), lu >= n (その他の場合)) jobu = 'U': u[][]に左特異ベクトルのm×n行列が入る. jobu = 'F': u[][]に左特異ベクトルのm×m行列が入る. Aの値域の直交補空間の正規直交基底を含む. jobu = 'W': jobv = 'V'かつjobt = 'T'かつm = nの場合, Aの代わりにA^Tが使われたときにu[][]は作業領域として使用される. この場合, VはA^Tの左特異ベクトルとしてu[][]の中で計算された後, 配列v[][]にコピーされる. この'W'オプションは, u[][]が長さn*nの作業領域として予約されていることを単に呼び出し元に示しているだけである. jobu = 'N': jobt = 'T'でない限りu[][]は参照されない. [in] ldv 二次元配列v[][]の整合寸法. (ldv >= 1 (jobv = 'N'の場合), ldv >= n (jobv = 'V', 'J'または'W'の場合)) [out] v[][] 配列 v[lv][ldv] (lv >= n) jobv = 'V'または'J': v[][]に右特異ベクトルのn×n行列が入る. jobv = 'W': jobu = 'U'かつjobt = 'T'かつm = nの場合, Aの代わりにA^Tが使われたときにv[][]は作業領域として使用される. この場合, UはA^Tの右特異ベクトルとしてv[][]の中で計算された後, 配列u[][]にコピーされる. この'W'オプションは, v[][]が長さn*nの作業領域として予約されていることを単に呼び出し元に示しているだけである. jobv = 'N': jobt = 'T'でない限りv[][]は参照されない. [out] work[] 配列 work[lwork] 作業領域. work[0], work[1]: scale = work[1]/work[0] はスケーリング係数で, 求められたAの特異値はscale*sva[0〜n-1]である(sva[]の説明を参照せよ). work[2]: = sconda. 列の均衡化後のAの条件数の推定値. (joba = 'E'または'G'の場合) scondaは sqrt(||(R^T*R)^(-1)||_1) の推定値である. これはdpoconを使って計算され, n^(-1/4)*sconda <= ||R^(-1)||_2 <= n^(1/4)*sconda となる. ここで, RはAのQR分解で得られる三角行列である. しかし, Rが不完全で, ランク数がnより小さいと判定された場合, sconda = -1 が返される. これは, 最小の特異値が失われた可能性があることを示す. SVD全部が必要な場合, 次の2つの条件数がアルゴリズムの解析のために有用である. これらは, 解法の詳細に精通した開発者／実装者のために提供される. work[3]: 1回目のQR分解の三角行列のスケーリングされた条件数の推定値. work[4]: 2回目のQR分解の三角行列のスケーリングされた条件数の推定値. 次の2つのパラメータは, jobt = 'T'の場合に計算される. これらは, 解法の詳細に精通した開発者／実装者のために提供される. work[5]: A^T*Aのエントロピー. これは, Probability simplexの点としてとられたdiag(A^T*A)/Trace(A^T*A)のShannonエントロピーである. work[6]: A*A^Tのエントロピー. [in] lwork 配列work[]のサイズ lworkはjobに依存する: 特異値だけが必要なとき(jobu = 'N'かつjobv = 'N') かつ   スケーリングされた条件数の推定値が不要 (!(joba = 'E'または'G')):     最小要求サイズは lwork >= max(2*m+n, 4*n+1, 7).     最適パフォーマンス(ブロックアルゴリズムのコード)のための最適値は lwork >= max(2*m+n, 3*n+(n+1)*nb, 7), ただし, nbは最適ブロックサイズである.   Aのスケーリングされた条件数の推定値が必要(joba = 'E'または'G'):     最小要求サイズは lwork >= max(2*m+n, n*n+4*n, 7).     最適パフォーマンス(ブロックアルゴリズムのコード)のための最適値は lwork >= max(2*m+n, 3*n+(n+1)*nb, n*n+4*n, 7), ただし, nbは最適ブロックサイズである. 特異値および右特異ベクトルが必要なとき ((jobv = 'V'または'J') かつ !(jobu = 'U' または 'F')):   最小要求サイズは lwork >= max(2*m+n, 4*n+1, 7).   最適パフォーマンスのためには lwork >= max(2*m+n, 3*n+(n+1)*nb, 7), ただし, nbは最適ブロックサイズである. 特異値および左特異ベクトルが必要なとき(!(jobv = 'V'または'J') かつ (jobu = 'U' または 'F')):   最小要求サイズは lwork >= max(2*m+n, 4*n+1, 7).   最適パフォーマンスのためには:     jobu = 'U'の場合 lwork >= max(2*m+n, 3*n+(n+1)*nb, 7),     jobu = 'F'の場合 lwork >= max(2*m+n, 3*n+(n+1)*nb, n+m*nb, 7),     ただし, nbは最適ブロックサイズである. 特異値および左右両方の特異ベクトルが必要なとき(jobu = 'U'または'F') かつ   jobv = 'V'の場合, 最小要求サイズは lwork >= max(2*m+n, 6*n+2*n*n, 7),   jobv = 'J'の場合, 最小要求サイズは lwork >= max(2*m+n, 4*n+n*n, 2*n+n*n+6, 7), 最適パフォーマンスのためには, lworkはさらにn+m*nbより大きくなくてはならない. ただし, nbは最適ブロックサイズである. lwork = -1の場合, 作業領域サイズの問い合わせとみなす. その呼び出しで使われたjobパラメータに応じたwork[]の最適サイズをwork[0]に, 最小サイズをwork[1]に返す. [out] iwork[] 配列 iwork[liwork] (liwork >= max(m + 3*n, 3)) 作業領域. iwork[0]: 最初のピボット選択付きQR分解後に判定されたランク数. jobaおよびjobrの説明を参照せよ. iwork[1]: 求められた0でない特異値の数. iwork[2]: 0でなければ警告メッセージである. iwork[2] = 1であれば, Aの列ノルムのいくつかが非正規化数である. このデータでは要求精度が保証されない. [out] info = 0: 正常終了 = -1: 入力パラメータ job1 の誤り (joba != 'C', 'E', 'F', 'G', 'A'および'R') = -2: 入力パラメータ jobu の誤り (jobu != 'U', 'F', 'W'および'N') = -3: 入力パラメータ jobv の誤り (jobv != 'V', 'J', 'W'および'N') = -4: 入力パラメータ jobr の誤り (jobr != 'R'および'N') = -5: 入力パラメータ jobt の誤り (jobt != 'T'および'N') = -6: 入力パラメータ jobp の誤り (jobp != 'P'および'N') = -7: 入力パラメータ m の誤り (m < 0) = -8: 入力パラメータ n の誤り (n < 0 または n > m) = -9: 入力パラメータ lda の誤り (lda < m) = -12: 入力パラメータ ldu の誤り (lduが小さすぎる) = -14: 入力パラメータ ldv の誤り (ldvが小さすぎる) = -17: 入力パラメータ lwork の誤り (lworkが小さすぎる) > 0: dgejsvが最大スイープ回数で収束しなかった. 求められた値は正しくないかもしれない. 出典 LAPACK