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

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 矩阵）。到目前为止，我一直在通过以下方式进行此操作：

将"A"的非零元素复制到2x2 gsl矩阵
使用 gsl LU 分解来反转 2x2 gsl 矩阵
将2x2倒矩阵复制到A_inv

这种方法的问题是我必须在每次迭代中复制矩阵两次。一次是一个较小的 n x n gsl 矩阵，一次是将生成的倒置 nxn gsl 矩阵复制到 A_inv.

我想知道是否有人知道更直接的方法。有没有办法使用一些 gsl 函数只反转矩阵的一部分并忽略任何零元素？像这样说：

A = NxN matrix
A_inv = invert_submatrix(A,n,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;
}

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

How to efficiently invert only part of an allocated matrix

c

matrix

inversion

gsl

submatrix