LAPACKのDGEEV関数メモ

スポンサーリンク

非対称行列の固有値計算関数DGEEVについてのメモ

作成途中なので、随時更新

外部関数

名前説明使用場面
DGEBAL非対称行列をスケーリングおよび置換してバランス化する。固有値計算の前処理
DGEBACKバランス化された行列を元の形に戻す。固有ベクトル計算後の後処理
DGEHRD一般行列をヘッセンベルク行列に変換する。固有値計算の簡約化
DHSEQRシュール形式に基づく固有値計算を行う。固有値とシュールベクトルの計算
DLACPY行列の要素をコピーする。行列のバックアップや複製
DLARTG2つの値を使ってGivens回転行列を生成する。QR分解や回転行列操作
DLASCL行列全体をスケーリングする。行列の条件数改善
DORGHRヘッセンベルク行列から直交行列を生成する。QR分解や直交行列の利用
DROT2つのベクトルにGivens回転を適用する。QR分解や回転操作
DSCALベクトル全体にスカラーを掛ける。ベクトルのスケーリング
DTREVC3右固有ベクトルを計算するための高精度な改良版。固有値計算後の固有ベクトル抽出
XERBLAエラーを報告するための標準サブルーチン。不正な引数やエラー条件の検出

外部サブルーチン

名前型説明
使用場面
LSAMELogical2つの文字が同じかどうかを大小区別なく比較する。サブルーチンの引数検証。
LAPACK関数内でオプション文字（例: 'N', 'T', 'V'）を検証する際に使用。
IDAMAXIntegerベクトル中で最大絶対値を持つ要素のインデックスを返す。ベクトル操作の際に最も大きい値を特定するために使用。BLASレベル1関数の一部。
ILAENVIntegerパフォーマンスを最適化するためのブロックサイズや他のパラメータを決定する。LAPACK内の自動チューニング。
LAPACKのサブルーチンが内部で最適なブロックサイズや並列化パラメータを取得する際に利用。
DLAMCHDouble precision浮動小数点数の特性（マシンイプシロンや最大値など）を返す。数値計算の精度や範囲を考慮。
（例: 非正規化数の扱い、丸め誤差の評価）。
DLANGEDouble precision行列のノルム（最大絶対値、フロベニウスノルムなど）を計算する。行列の性質を評価。
条件数の推定や前処理として行列のスケールを評価する際に使用。
DLAPY2Double precision2つの値\(x\) と\(y\)の平方和の平方根（\(\sqrt{x^2+y^2}\) ）を安全に計算する。オーバーフローやアンダーフローを避ける安全な方法で三平方の計算を行う際に使用。
DNRM2Double precisionベクトルの2ノルム（ユークリッド長さ）を計算する。ベクトルの正規化や長さを評価する場面。例えば、QR分解やランク計算などで使用されることが多い。

スケーリング

340～353行

*
*     Scale A if max element outside range [SMLNUM,BIGNUM]
*
      anrm = dlange( 'M', n, n, a, lda, dum )
      scalea = .false.
      IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
         scalea = .true.
         cscale = smlnum
      ELSE IF( anrm.GT.bignum ) THEN
         scalea = .true.
         cscale = bignum
      END IF
      IF( scalea )
     $   CALL dlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )

この部分は行列 \(A\) のスケーリングを行う処理で、数値安定性を高めるための重要な前処理。

スケーリングの目的

行列のスケーリングは、以下のような目的で行われます：

  1. 数値安定性の向上:
    • 数値計算では、非常に小さい値（アンダーフロー）や非常に大きい値（オーバーフロー）があると計算結果の精度が低下します。
    • 行列の要素が適切な範囲に収まるようにスケーリングすることで、計算の精度を保ちます。
  2. 演算の効率化:
    • 行列の要素が適切にスケーリングされていると、アルゴリズムがより効率的に動作します。

コードの詳細

スケーリング条件のチェック

anrm = dlange( 'M', n, n, a, lda, dum )
  • dlange('M', ...):
    • 行列 \(A\) の最大絶対値（max norm）を計算します。
    • これにより、行列\(A\) のスケーリングが必要かどうかを判定します。
IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
   scalea = .true.
   cscale = smlnum
ELSE IF( anrm.GT.bignum ) THEN
   scalea = .true.
   cscale = bignum
END IF
  • anrm<smlnum: 行列の要素が小さすぎる場合（アンダーフローの可能性）
    • スケーリング係数 cscale を smlnum（最小有効値）に設定。
  • anrm>bignum: 行列の要素が大きすぎる場合（オーバーフローの可能性）。
    • スケーリング係数 cscale を bignum（最大有効値）に設定。

スケーリング処理の実行

IF( scalea )
$   CALL dlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )

dlascl:

  • 行列 \(A\) の全要素をスケーリングします。
  • スケーリングの形式は 'G'（一般行列用）で、行列の全要素を anrm/cscale に基づいて調整します。
  • 引数の意味:
    • 'G': スケーリングの種類（一般行列）。
    • anrm: 元のスケーリング係数（最大絶対値）。
    • cscale: 新しいスケーリング係数（smlnum または bignum）。
    • n, n: 行列のサイズ。
    • a: スケーリングされる行列。

具体例

行列 \(A\) の例

初期状態:
$$
A = \begin{bmatrix}
10^{300} & 2 \\
3 & 10^{290}
\end{bmatrix}
$$

NumPyにおける smlnum と bignum 値：
smlnum（最小正規化値）: smlnum=最小正規化値（マシンの最小値）÷マシンのイプシロン
bignum（最大値）:bignum=1÷smlnum

  1. dlange による行列要素の絶対値の最大値anrmの取得:
    • 最大絶対値 \(\text{anrm} = 10^{300}\)
  2. スケーリング要否の確認
    • 条件1: anrm>bignum
      → 行列の最大値が計算の安全範囲を超えている場合。
    • 条件2: anrm<smlnum
      → 行列の最大値が非常に小さくアンダーフローの可能性がある場合。
    • 例では、bignum=\(9.9792015476736×10^{291}\)だとすると
      anrm>bignumになり条件１となる。
  3. スケーリング係数の決定:
    • anrm>bignumの場合：
      $$
      \text{cscale} = \frac{\text{bignum}}{\text{anrm}}$$
      → 行列 \(A\) の要素がすべて bignum 以下になるようにスケーリング。
    • anrm<smlnumの場合：
      $$
      \text{cscale} = \frac{\text{smlnum}}{\text{anrm}}
      $$
      → 行列 \(A\) の要素がすべて smlnum 以上になるようにスケーリング。
    • 例の場合
      $$
      \text{cscale}=\frac{\text{bignum}}{\text{anrm}}​=9.9792015476736×10^{−9}
      $$
  4. dlascl によるスケーリング:
    • \(A\) の全要素を \(\text{cscale}=9.9792015476736×10^{-9}\) にスケーリング。

スケーリング後の行列

$$
A’ =
\begin{bmatrix} 10^{300} × 9.9792015476736 × 10^{-9} & 2 × 9.9792015476736 × 10^{-9} \\
3 × 9.9792015476736 × 10^{-9} & 10^{290} × 9.9792015476736 × 10^{-9}
\end{bmatrix}
$$

結果：
$$
A’ =
\begin{bmatrix} 9.9792015476736 × 10^{291} & 1.99584030953472 × 10^{-8} \\
2.99376046430208 × 10^{-8} & 9.9792015476736 × 10^{281}
\end{bmatrix}
$$

スケーリングの影響

  • スケーリングにより指数部は調整されるものの、少数部（浮動小数点の有効桁）が操作されるため、以下の問題が生じる可能性があります：
    1. 丸め誤差の増加:
      • スケーリング後、元の値に戻す際に丸め誤差が累積する可能性があります。
      • 特に行列の要素間でスケーリング比が大きい場合、元のスケールを戻す際の精度が悪化します。
    2. 有効桁数の減少:
      • 浮動小数点数の有効桁が一定であるため、桁の移動によって元の有効数字が減る場合があります。

スケーリングのメリットとデメリット

メリットデメリット
数値範囲を調整してオーバーフローやアンダーフローを防ぐ。少数部の桁が増えることで丸め誤差が増える可能性がある。
計算の安定性を向上させ、特に固有値計算での精度が向上する。有効桁数が減少し、元に戻す際の復元精度が低下する可能性がある。

LAPACKがスケーリングを使う理由

スケーリングによるデメリット（丸め誤差や有効桁数の減少）はあるものの、特に固有値計算や行列分解のアルゴリズムでは、以下の理由からスケーリングが重要視されています：

  1. 計算精度の向上:
    • オーバーフローやアンダーフローが発生する場面では、スケーリングを行わない場合よりも精度が大幅に低下するリスクがあります。
    • スケーリングを行うことでそのリスクを軽減します。
  2. アルゴリズムの数値安定性:
    • 特に、ヘッセンベルク変換やシュール形式への変換では、スケーリングがないとアルゴリズム自体が不安定になることがあります。

バランシング

354-359行

*
*     Balance the matrix
*     (Workspace: need N)
*
      ibal = 1
      CALL dgebal( 'B', n, a, lda, ilo, ihi, work( ibal ), ierr )

バランシング（dgebal）の役割

目的

  • 固有値計算の精度と数値安定性を向上させる。
  • 条件数を改善し、行列をより扱いやすい形に変換する。
  • 行列のスケーリングや並べ替えを行い、大きすぎるまたは小さすぎる要素の影響を軽減する。

具体的な処理内容

  1. 行列のスケーリング:
    • 行列の行や列をスケーリングし、要素の絶対値の範囲を狭めます。
    • 行列の条件数を改善します。
  2. 行列の並べ替え（縮約）:
    • 行列のゼロ行またはゼロ列を識別し、非ゼロ部分だけを操作するように行列を縮約します。
    • 縮約後の有効範囲が ilo から ihi で示されます。

コードの詳細解説

dgebal の呼び出し

CALL dgebal( 'B', n, a, lda, ilo, ihi, work(ibal), ierr )
  • 引数の説明:
    • 'B': 処理オプションを指定。
      • 'B': 行列のスケーリングと並べ替えの両方を実行。
      • 'S': スケーリングのみを実行。
      • 'P': 並べ替え（行・列の縮約）のみを実行。
    • n: 行列の次元（\(n \times n\)）。
    • a: 入力行列。バランシング後、変換された行列が上書きされる。
    • lda: 行列のリーディングディメンション。
    • ilo, ihi: バランシング後の行列の有効範囲（縮約された部分）。
    • work(ibal): 作業領域として使用される配列。
    • ierr: エラーフラグ。

バランシング後の情報

  • ilo と ihi:
    • 行列の有効範囲を示します。
    • \(i = 1 \text{ から } ilo-1\) と \(i = ihi+1 \text{ から } n\)の範囲はゼロ行またはゼロ列として扱われ、固有値計算に影響しない部分です。
  • work(ibal):
    • スケーリング係数が保存されます。
    • 固有ベクトルを求める際、この情報を用いて元のスケールに戻します。

処理の流れ

  1. 行列のスケーリング:
    • 各行や列が適切な係数でスケーリングされます。
    • 例えば、ある行の要素がすべて非常に大きい場合、その行を適切にスケールダウン。
  2. 行列の並べ替え:
    • 行列のゼロ行またはゼロ列が検出され、縮約されます。
    • 有効な非ゼロ部分（iloからihi）のみが後続の計算対象。
  3. 結果の出力:
    • 行列 a がバランス化された形に更新されます。
    • スケーリング係数が work(ibal) に保存され、行列の元のスケールに戻すために使用されます。

バランシングの利点と影響

利点

  1. 計算精度の向上:
    • 行列の要素が非常に大きい/小さい場合でも、適切なスケーリングにより数値安定性が向上します。
  2. 効率の向上:
    • ゼロ行やゼロ列を縮約することで、計算量が減少します。

潜在的な課題

  1. 丸め誤差:
    • 行列を元のスケールに戻す際に丸め誤差が蓄積する可能性があります。
  2. 固有ベクトルの整合性:
    • 固有ベクトル計算時にはスケーリング情報を適切に適用して元のスケールに戻す必要があります。

具体例

入力行列

$$
A = \begin{bmatrix} 10 & 1 & 0 \\
1000 & 10 & 0 \\
0 & 0 & 1
\end{bmatrix}
$$

スケーリングと縮約

1. 各行と列のノルム計算

  • 行ノルム（最大絶対値）: \(\text{Row norms} = [10, 1000, 1]\)
  • 列ノルム（最大絶対値）: \(\text{Column norms} = [1000, 10, 1]\)
  • この例では第1行は 1/10、第2行は 1/1000 のスケーリングとなります。

2. スケーリング適用

  • スケーリングの要否判定
    • ノルム比=\(\text{ノルム比} = \frac{\text{max(row, col)}}{\text{min(row, col)}}\)
    • スレッショルド（例: \(10^2\)）より大きい場合にスケーリングを適用。
    • 列ノルム：
      • 第一列: 1000
      • 第二列: 10
      • 第三列: 1
    • ノルム比（例: 第一列と第二列）：
      $$
      \frac{\text{第一列ノルム}}{\text{第二列ノルム}} = \frac{1000}{10} = 100
      $$
      → スレッショルド（例: \(10^2\)）を超えるため、スケーリングが適用。
    • 第二列と第三列：
      $$
      \frac{\text{第二列ノルム}}{\text{第三列ノルム}} = \frac{10}{1} = 10
      $$
      → スレッショルド（例: \(10^2\)）を満たさないため、スケーリング不要。
    • 偏りが起こらないようにノルム比の底は、シフトして行きます。
  • 各行および列を、その最大絶対値で正規化するようにスケーリングを適用します。
    • 第1行: すべての要素を \(\frac{1}{10}\)​ 倍。
    • 第2行: すべての要素を \(\frac{1}{1000}\) 倍。
    • 第3行: 変化なし（ノルムが1）。
    • 第1列: \(\frac{1}{1000}\) 倍。
    • 第2列: 変化なし（ノルムが10）。
    • 第3列: 変化なし（ノルムが1）。

3. 縮約（ilo と ihi）

  • 行列 AAA にゼロ行やゼロ列は存在しないため、ilo = 1 と ihi = 3 になります。
  • つまり、有効な範囲は行列全体です。

4. スケーリング情報の格納

  • work(ibal) に各行および列のスケーリング係数が格納されます。
    $$\text{work(ibal)} = [10, 1000, 1, 1000, 10, 1]$$
    • 前半3つが行のスケーリング係数。
    • 後半3つが列のスケーリング係数。
  • ibalはスケーリング係数の格納開始位置のため、1を指定します。（FORTRANは0からでは無く１から）

結果:
$$
A’ =
\begin{bmatrix} 1 & 0.1 & 0 \\
1 & 0.01 & 0 \\
0 & 0 & 1
\end{bmatrix}
$$

  • ilo = 1, ihi = 3（全範囲有効）。
  • work(ibal) にスケール係数が保存。\(\text{work(ibal)} = [10, 1000, 1, 1000, 10, 1]\)

デスケーリング（work配列の使途）

スケーリングによって行列式や固有ベクトルが変化しますが、work 配列を利用することで元のスケールに戻せます。

work 配列によるデスケーリング

  1. スケーリング係数の格納
    • dgebal は行と列のスケーリング係数を work 配列に格納します。
    • 例えば：
      • 行ごとのスケーリング係数：work(1:n)
      • 列ごとのスケーリング係数：work(n+1:2n)
  2. デスケーリングに使用
    • 行列式(Determinant)や固有ベクトルを元のスケールに戻す際、work に保存されたスケーリング係数が利用されます。
  3. 具体例: 行列式の復元
    • 元の行列式を復元するには、work に保存されたスケーリング係数を逆に適用します。
    • 行列式の復元：
      $$
      \text{det}(A) = \frac{\text{det}(A’)}{\text{det}(D_1) \cdot \text{det}(D_2)}
      $$
    • ここで、\(\text{det}(D_1) \) と \(\text{det}(D_2) \) は work のスケーリング係数を掛け合わせたものです。
  4. 固有ベクトルの復元
    • 固有ベクトルもスケーリングによって影響を受けますが、work に基づいて元のスケールに戻せます。
    • 固有ベクトル xxx を復元する場合：
      $$
      x_{\text{original}} = D_2^{-1} \cdot x
      $$

入力行列

$$
A =
\begin{bmatrix} 10 & 1 & 0 \\
1000 & 10 & 0 \\
0 & 0 & 1
\end{bmatrix}
$$

スケーリング

  1. 行ごとにスケーリング：
    • 第1行：\(\frac{1}{10}\)
    • 第2行：\(\frac{1}{1000}\)
    • 第3行：そのまま。
  2. 列ごとにスケーリング：
    • 第1列：\(\frac{1}{1000}\)
    • 第2列：そのまま。
    • 第3列：そのまま。

結果：
$$
A’ =
\begin{bmatrix}
0.01 & 0.1 & 0 \\
1 & 0.01 & 0 \\
0 & 0 & 1
\end{bmatrix}
$$

スケーリング後の行列式

$$
\text{det}(A’) = \text{det}(A) \cdot (1/10) \cdot (1/1000) \cdot (1/1000)
$$

行列式の復元

  • work に格納されたスケーリング係数： \(\text{work} = [1/10, 1/1000, 1, 1/1000, 1, 1]\)
  • 行列式の復元： \(\text{det}(A) = \text{det}(A’) \cdot 10 \cdot 1000 \cdot 1000\)

スケーリング関数とバランシング関数内のスケーリングの役割の違い

dlasclが粗い均しで、dgebalがより滑らかな均しというイメージです。

dlascl (スケーリング関数)と dgebal（バランシング関数） の役割の違い

1. dlascl の役割

  • 目的:
    • 行列全体をスケーリングし、要素を適切な範囲（[smlnum,bignum]）に収める。
    • オーバーフローやアンダーフローを防ぐための前処理として使用される。
  • 適用範囲:
    • 行列全体に一様にスケーリング係数を適用します。
    • 行や列ごとのスケーリングは行わず、行列の最大値や最小値に基づいて調整します。

2. dgebal の役割

  • 目的:
    • 行や列ごとにスケーリングを適用し、行列の条件数を改善。
    • シフトによる縮約によってゼロ行やゼロ列を効率的に扱う。
  • 適用範囲:
    • 各行や列に対して異なるスケーリングを適用します。
    • また、ゼロ行やゼロ列を検出して右下隅に追いやる縮約も同時に行います。

両者の違いと補完関係

特徴dlascldgebal
スケーリングの範囲行列全体を均一にスケーリング行や列ごとに個別のスケーリングを適用
スケーリングの目的演算範囲（数値範囲）を安全に調整条件数の改善、計算の安定性向上
ゼロ行・列の処理行わないゼロ行やゼロ列を検出し縮約（シフト）
適用タイミング固有値計算の全体的な安定性を確保する前処理固有値計算に特化した安定性向上の前処理

なぜ dlascl が必要なのか

  1. 数値範囲の調整
    • dlascl は、行列全体をスケーリングして「安全な数値範囲」に収める役割を担います。
    • 特に、行列の要素が非常に大きい（オーバーフロー）または非常に小さい（アンダーフロー）場合に、計算全体の安定性を確保します。
  2. dgebal では補えない部分
    • dgebal は行や列ごとのスケーリングを行いますが、行列全体の数値範囲を調整することを目的としていません。
    • 例えば、行列全体が非常に大きな値（例えば \(10^{300}\)）を持つ場合、dgebal だけではその値を範囲内に収めることはしません。
  3. 後続処理のための準備
    • dlascl は、固有値計算や他のアルゴリズムが過剰なスケールの行列で不安定にならないようにするための基盤を提供します。

なぜ dlascl と dgebal の両方が使われるのか

  1. dlascl は全体の数値範囲を調整
    • 例: 行列全体のスケールが \(10^{-300}\) ～ \(10^{300}\) のように極端な場合、これを \([10^{-292}, 10^{292}]\)に収めます。
    • スケーリングは全体的に適用されるため、行列の個々の行や列を均一に縮小または拡大します。
  2. dgebal は行列の内部構造を改善
    • スケーリング後、行や列ごとのスケールのばらつきをさらに均一化します。
    • ゼロ行やゼロ列を検出して縮約することで、計算量を削減します。

例を用いた違いの説明

行列

$$
A = \begin{bmatrix}
10^{300} & 2 \\
3 & 10^{-300}
\end{bmatrix}
$$

1. dlascl によるスケーリング

  • 最大値 \(\text{anrm} = 10^{300}\) を基準にスケーリング。
  • 結果：
    $$
    A’ = \begin{bmatrix}
    10^{8} & 2 \times 10^{-292} \\
    3 \times 10^{-292} & 10^{-292}
    \end{bmatrix}
    $$

2. dgebal によるバランシング

  • 行と列のノルムに基づいてさらにスケーリング。
  • 例えば、第1列は \(10^{8}\) と \(3 \times 10^{-292}\) のばらつきがあるためスケール調整。
  • ゼロ行列部分（もしあれば）を右下に移動。
  • 結果：
    $$
    A” = \begin{bmatrix}
    1 & 0.2 \times 10^{-292} \\
    1 \times 10^{-292} & 1
    \end{bmatrix}
    $$
スポンサーリンク

コメント

タイトルとURLをコピーしました