ブロック三角行列の不変部分空間の主角度の計算

目的

ブロック三角行列の 2 つの不変部分空間の相対的な位置に関する情報を得る。

解決方法

部分空間はいくつかのベクトルのスパンとして表現され、部分空間の相対的な位置は部分空間の間の主角度のセットを計算して得ることができると仮定します (2 つの部分空間の間の主角度の計算を参照)。さらに、ブロック三角行列の不変部分空間はシルベスター行列方程式を使用して解くものとします。使用するソルバーは、行列の特性に依存します。

三角行列の対角ブロックがどちらも上三角の場合、LAPACK ?trsyl ルーチンを使用します。
三角行列の対角ブロックがどちらも大きくなく、上三角でない場合、LAPACK 線形ソルバーを使用します。
三角行列の対角ブロックがどちらも大きく、上三角で、疎の場合、oneMKL PARDISO ソルバーを使用します。

ソースコード: サンプル (https://www.intel.com/content/dam/develop/external/us/en/documents/mkl-cookbook-samples-120115.zip (英語)) のファイルを参照してください。

ANGLES/uep_subspace1/main.f
ANGLES/uep_subspace2/main.f
ANGLES/uep_subspace3/main.f

LAPACK ?trsyl を使用してシルベスター行列方程式を解く

   CALL DTRSYL('N', 'N', -1, K, N-K, AA, N, AA(K+1,K+1), N, 
&                AA(1,K+1), N, ALPHA, INFO) 
    IF(INFO.EQ.0) THEN 
        PRINT *,"DTRSYL completed, SCALE=",ALPHA 
    ELSE IF(INFO.EQ.1) THEN 
        PRINT *,"DTRSYL solved perturbed equations」 
    ELSE 
        PRINT *,"DTRSYL failed. INFO=",INFO 
        STOP 
    END IF

LAPACK 線形ソルバーを使用してシルベスター行列方程式を解く

   REAL*8 AA(N,N), FF(K*(N-K)), AAA(K*(N-K),K*(N-K)) 
… 
! シルベスター方程式の密係数行列を形成 
    CALL SYLMAT(K, AA, N, N-K, AA(K+1,K+1), N, -1D0, 1D0, AAA, NK, & INFO) 
! SYLMAT により返された INFO を処理 
… 
! シルベスター方程式に相当する連立線形方程式の 
! 右辺を形成
    DO I = 1, K 
        DO J = 1, N-K 
            FF((J-1)*K+I) = AA(I,J+K) 
        END DO 
    END DO 
! 連立線形方程式を解く 
    CALL DGESV(NK, 1, AAA, NK, IPIV, FF, NK, INFO ) 
! DGESV により返された INFO を処理

oneMKL PARDISO を使用してシルベスター行列方程式を解く

   REAL*8 AA(N,N), FF(K*(N-K)), VAL(K*(N-K)*(N-1)) 
    INTEGER ROWINDEX(K*(N-K)+1), COLS(K*(N-K)*(N-1)) 
… 
! シルベスター方程式の疎係数行列を形成 
    CALL FSYLVOP(K, AA, N, N-K, AA(K+1,K+1), N, -1D0, 1D0, COLS, 
&                 ROWINDEX, VAL, INFO) 
! FSYLVOP により返された INFO を処理 
… 
! シルベスター方程式の右辺を形成 
    DO I=1,K 
        DO J=1,N-K 
            FF((J-1)*K+I) = AA(I,J+K) 
        END DO 
    END DO 

    CALL PARDISOINIT (PT, 1, IPARM) 
    CALL PARDISO (PT, 1, 1, 11, 13, NK, VAL, ROWINDEX, 
&                  COLS, PERM, 1, IPARM, 1, FF, X, IERR) 
! PARDISO により返された IERR を処理 
…

説明

使用するルーチン
タスク	ルーチン	説明
上三角対角ブロックを含む行列についてシルベスター行列方程式を解く	dtrsyl	実数擬似三角または複素数三角行列についてシルベスター方程式を解きます。
小さく上三角でない行列についてシルベスター方程式を解く	dgesv	正方行列 `A` および複数の右辺を含む連立線形方程式の解を計算します。
小さくなく、上三角でない行列についてシルベスター方程式を解く	pardiso	単一または複数の右辺を含む 1 セットの疎線形方程式の解を計算します。

行列の不変部分空間の間の主角度を決定するには、最初に N x N 行列をブロック三角形式で表します。

ここで、対角ブロック A および B はそれぞれ、次数 k および N-k の正方行列です。I_k が k の単位行列を示す場合、次の等式

は、標準基底の最初の k ベクトルのスパンが行列 A の変換に対して不変であることを意味します。

別の不変部分空間は、次の合成行列の列のスパンとして得ることができます

ここで、X は長方形の k x (N - k) 行列です。積を計算します。

X がシルベスター方程式 XB - AX = F の解の場合、最後の方程式の結果は次のようになります。

これは、の列でスパンされた部分空間の不変性を表しています。

2 つ目の不変部分空間の基底を直交させるには、QR 因数分解を使用します。

ここで、C は k x (N-k) 行列で、S は (N-k) x (N-k) 行列です。C および S は方程式 C^TC + S^TS = I_N-k を満たします。ここで、R は次数 N-k の上三角正方行列です。C の SVD を使用してこれらの 2 つの不変部分空間の間の主角度を計算します。

Σ の対角要素は主角度の余弦です。

シルベスター方程式の行列

シルベスター方程式 αAX + βXB = F を考えます。

ここで、正方行列 A および B の次数はそれぞれ M と N で、α と β はスカラーです。F は指定される M x N 行列です。

X は求める M x N 行列です。

この行列方程式は、ベクトル x と右辺ベクトル f の MN 成分が不明な、系の MN 線型方程式と見なすことができます。

x = (x₁₁, x₂₁, ..., x_M1, x₁₂, x₂₂, ..., x_M2, ..., x_1N, x_2N, ..., x_MN)^T

f = (f₁₁, f₂₁, ..., f_M1, f₁₂, f₂₂, ..., f_M2, ..., f_1N, f_2N, ..., f_MN)^T

次数 MN の行列は、2 つの行列の合計として表すことができます。1 つの行列は、行列 X (左から) と行列 A の乗算に相当し、サイズ M x M のブロック形式で表現できます。この行列は、対角線上で N ブロックのブロック対角行列を形成します。

合計の別の行列は、行列 X (右から) と行列 B の乗算に相当します。

ここで、I_M は次数 M の単位行列を表します。つまり、係数行列は次のようになります。

この行列は、MN 要素の行に M + N - 1 の非ゼロ要素を含む疎行列です。このため、oneMKL PARDISO スパースソルバーを効率的に使用できます (oneMKL PARDISO 向けに CSR 形式で係数行列を形成するコード ANGLES/source/fsylvop.f を参照)。しかし、M および N が小さい場合は、oneMKL LAPACK 線形ソルバーがより効率的です (dgesv を使用する密行列として係数行列を形成するコード ANGLES/source/sylmat.f を参照)。