如何有效地仅反转分配矩阵的一部分

Question

我有一个算法，它分配一个具有预定义大小 N x N 的复数双矩阵“A”。元素最初为零。我还分配了大小为 N x N 的矩阵来存储逆矩阵“A_inv”。 在算法期间，“A”的元素被填充。在每次迭代 i 中，我最终得到一个大小为 i x 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 invert_submatrix() 只考虑A 的n x n 部分。此外，原始矩阵“A”不能被这种反转改变。 也许最后的需求使得无论如何都必须复制矩阵，在这种情况下，它不会比我现在做的更有效率。也就是说，gsl 算法往往比我通常想出的任何方法都高效得多。因此，对此的想法仍然非常受欢迎。

Answer 1

遗憾的是，正如 GSL 所做的那样，它的 LU 分解已经到位，如果您要求不修改它，我不确定是否可以在不先从 A 复制子矩阵的情况下做到这一点。但是，您可以使用 matrix view 直接从 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;
}