1

事前定義されたサイズ N x N の複雑な double 行列 "A" を割り当てるアルゴリズムがあります。要素は最初はゼロです。また、サイズ N x N の行列を割り当てて、逆行列 "A_inv" を格納しました。アルゴリズム中に、「A」の要素が満たされます。各反復 i で、最終的にサイズ ix i の部分行列になります。したがって、N=4 の 2 回目の反復では次のようになります。

| x00 x01 0.0 0.0 |
| x10 x11 0.0 0.0 |
| 0.0 0.0 0.0 0.0 |
| 0.0 0.0 0.0 0.0 |

ここで、x はゼロ以外の値を示します。ここで、行列の非ゼロ部分 (この例では 2x2 行列) を反転したいと考えています。これまでのところ、私は次の方法でこれを行ってきました。

  1. "A" の非ゼロ要素を 2x2 gsl 行列にコピーします
  2. gsl LU 分解を使用して 2x2 gsl 行列を反転します
  3. 2x2 反転行列を A_inv にコピーします

このアプローチの問題は、各反復で行列を 2 回コピーする必要があることです。1 回はより小さい nxn gsl 行列に、もう 1 回は結果の逆 nxn gsl 行列を A_inv にコピーします。

誰かがもっと直接的な方法を知っているかどうか疑問に思っていました。gsl 関数を使用して行列の一部のみを反転し、ゼロ要素を無視する方法はありますか? 次のように言います。

A = NxN matrix
A_inv = invert_submatrix(A,n,n)

ここで、n < N です。ここでinvert_submatrix()は、A の nxn 部分のみを考慮します。さらに、元の行列 "A" は、この反転によって変更されてはなりません。たぶん、最後の要求により、とにかくマトリックスをコピーする必要が生じます。その場合、私が今行っていることよりも効率的ではありません。とは言っても、gsl アルゴリズムは、私が通常思いつくものよりもはるかに効率的である傾向があります。したがって、これに関する考えは大歓迎です。

4

1 に答える 1

0

悲しいことに、GSL では LU 分解が行われています。部分行列を変更しAないようにする必要がある場合、最初から部分行列をコピーせずにそれを実行できるかどうかはわかりません。ただし、行列ビューを使用して、LU 分解から直接逆を構築することができます。逆を構築してからコピーする必要はありません。

gsl_matrix *invert_submatrix( const gsl_matrix *m, size_t sub_size )
{
    gsl_matrix *inv = gsl_matrix_calloc( m->size1, m->size2 );

    // Create views onto the submatrices we care about 
    gsl_matrix_const_view m_sub_view = 
        gsl_matrix_const_submatrix(m, 0, 0, sub_size, sub_size);

    gsl_matrix_view inv_sub_view = 
        gsl_matrix_submatrix(inv, 0, 0, sub_size, sub_size);

    const gsl_matrix *m_sub = &m_sub_view.matrix;
    gsl_matrix *inv_sub = &inv_sub_view.matrix;

    // Create a matrix for the LU decomposition as GSL does this inplace.
    gsl_permutation *perm = gsl_permutation_alloc(sub_size);
    gsl_matrix *LU = gsl_matrix_alloc(sub_size, sub_size);
    gsl_matrix_memcpy(LU, m_sub);

    int s;
    gsl_linalg_LU_decomp(LU, perm, &s);
    gsl_linalg_LU_invert(LU, perm, inv_sub);

    gsl_matrix_free(LU);
    gsl_permutation_free(perm);

    return inv;
}
于 2015-01-15T12:07:44.573 に答える