Как эффективно инвертировать только часть выделенной матрицы

Question

У меня есть алгоритм, который выделяет сложную двойную матрицу «A» с заранее определенным размером N x N. Элементы изначально равны нулю. Я также выделил матрицу размера N x N для хранения обратного «A_inv». Во время алгоритма элементы «А» заполняются. На каждой итерации 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 | 

где иксы показывают некоторое ненулевое значение. Теперь я хочу инвертировать ненулевую часть матрицы (матрица 2x2 в этом примере). До сих пор я делал это следующим образом:

1. копию ненулевые элементы «А» в матрице 2х2 GSL
2. использование разложения GSL LU, чтобы инвертировать матрицу 2х2 GSL
3. скопируйте преобразованную матрицу 2x2 в A_inv

Проблема с этим подходом заключается в том, что мне приходится копировать по матрице дважды на каждой итерации. Однажды на меньшую n x n gsl-матрицу и один раз, чтобы скопировать полученную инвертированную матрицу nxn gsl в A_inv.

Мне было интересно, знает ли кто-нибудь о более прямом пути. Есть ли способ использовать некоторую функцию gsl только для инвертирования части матрицы и игнорирования любых нулевых элементов? сказать что-то вроде этого:

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

где п < Н. Здесь invert_submatrix() будет рассматривать только п х п часть А. Кроме того, исходная матрица «А» не должны быть изменены с помощью этой инверсии. Возможно, последнее требование требует скопировать по матрице в любом случае, и в этом случае оно будет неэффективным, чем то, что я делаю сейчас. Тем не менее, алгоритмы 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; 
}