1

我有一个算法，它分配一个具有预定义大小 N x N 的复数双矩阵“A”。元素最初为零。我还分配了大小为 N x N 的矩阵来存储逆矩阵“A_inv”。在算法期间，“A”的元素被填充。在每次迭代 i 时，我都会得到一个大小为 ix i 的子矩阵。所以第二次迭代看起来像这样，其中 N=4：

| 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

这种方法的问题是我必须在每次迭代中复制一个矩阵两次。一次到较小的 nxn gsl 矩阵，一次将生成的反向 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 回答