使用gsl库实现matlab中矩阵的除法(\)

说明

• 在matlab中 \ 与 .\ 是不一样的计算法则，前者要求逆后再相乘，后者只是对应相除。

---

百度了很多，但是gsl库中好像没有直接实现a \ b 这种矩阵除法的函数。该等式在matlab中等价与 inv(a) * b ，这是在a为方阵的情况下，若a不为方阵，那么等价与pinv(a)*b。但是gsl库中也没发现inv()，pinv() 函数，需要利用gsl库自己现实。

eg:当矩阵是个方阵时，现实复数矩阵的inv()，矩阵的求逆。a\b 等价于inv(a) *b。

void gsl_matrix_complex_inv(gsl_matrix_complex* a) {
    size_t n = a->size1;
    gsl_matrix_complex *temp1 = gsl_matrix_complex_calloc(n, n);
    gsl_matrix_complex_memcpy(temp1, a);
    gsl_permutation *p = gsl_permutation_calloc(n);
    int sign = 0;
    gsl_linalg_complex_LU_decomp(temp1, p, &sign);
    gsl_matrix_complex *inverse = gsl_matrix_complex_calloc(n, n);
    gsl_linalg_complex_LU_invert(temp1, p, inverse);
    gsl_matrix_complex_memcpy(a, inverse);
    gsl_permutation_free(p);
    gsl_matrix_complex_free(temp1);
    gsl_matrix_complex_free(inverse);
}

当矩阵不是方阵时，实现pinv() 伪逆函数，
我使用直接求导的方式，即 pinv(a) = inv(a^ *a)*a^, 其中a^ 为a的转置矩阵，即行列交换。注意这里是叉乘，不是点乘，在matlab中俩种的计算方式也是不同的。a\b 等价于inv(a^ *a)*a^ *b。

• 所以这里我们要先求转置矩阵，转置矩阵叉乘原矩阵，再调用已经实现的gsl_matrix_complex_inv函数，再叉乘转置矩阵。最后叉乘另一个矩阵即可。
• gsl库中求转置的函数为gsl_matrix_complex_transpose_memcpy

void gsl_matrix_complex_pinv(gsl_matrix_complex* p) {
    gsl_matrix_complex *gma_s = gsl_matrix_complex_alloc(p->size2, p->size1);
    gsl_matrix_complex_transpose_memcpy(gma_s, p);  //转置矩阵
    gsl_matrix_complex *gb = gma_s;
    gsl_matrix_complex_sub(gb, p) //转置叉乘原矩阵
    gsl_matrix_complex_inv(gb); //求逆
    gsl_matrix_complex_sub(gb, p) //逆矩阵叉乘转置矩阵
}

gsl_matrix_complex_sub是点乘还是叉乘来着，我是自己实现的叉乘，若是需要也可自己实现。