Inverse of a matrix - 使用 LU 演算法

在 GSL 中已經實作了矩陣的反轉換函式

使用 GSL 定義的矩陣
首先應該知道怎麼在 GSL 中使用其定義的矩陣，這部份應參考 GSL Reference Manual 的 Vectors and Matrices 一節的 Matrices 小節。必須注意的是，GSL採用的是 C 語言的架構，即使用 struct 來實作矩陣 gsl_matrix ，再透過許多函式來操作之。

1. 宣告及配置記憶體給矩陣

/* 宣告及配置記憶體 */
gsl_matrix * m1 = gsl_matrix_alloc(m, n);

/* 宣告及配置記憶體，並將所有元素初始化為0 */
gsl_matrix * m2 = gsl_matrix_calloc(m, n);

　ps. gsl_matrix 宣告的矩陣，其元素的型態為 double ，若要使用其他型態，
　　 可使用 gsl_matrix_float, gsl_matrix_int, .. 等。

2.初始化矩陣

/* 設定矩陣 m1 內所有元素的值為 1.0 */
gsl_matrix_set_all (m1, 1.0);

/* 設定矩陣 m1 內所有元素的值為 0 */
gsl_matrix_set_zero (m1);

/* 將矩陣 m1 內對角線的元素值均設為 1 ，其餘為 0，即為單位矩陣 */
gsl_matrix_set_identity (gsl_matrix * m);

3.讀取及寫入矩陣元素

/* 讀出矩陣m1的第i列(row)第j行(column)的值 */
double v1 = gsl_matrix_get(m1, i, j);

/* 將矩陣m1的第i列(row)第j行(column)的值設為 5 */
gsl_matrix_set(m1, i, j, 5);

/* 直接傳回矩陣m1的第i列(row)第j行(column)元素的指標 */
double * ptr = gsl_matrix_ptr(m1, i, j);

　ps. 預設的情況下，這3個函式會檢查傳入的 i, j 值範圍，在大量運算的時候會拖慢系統執行速度，此時可以加上 preprocessor definition GSL_RANGE_CHECK_OFF來關掉此功能。

4.矩陣的操作

　GSL 在矩陣操作上提供了矩陣相加(gsl_matrix_add)、相減(gsl_matrix_sub)、乘上一個常數(gsl_matrix_scale)、加上一個常數(gsl_matrix_add_constant)等函式，比較令人不解的是，此處提供的 gsl_matrix_mul_elements 函式，並非正常的矩陣相乘運算，因此矩陣相乘的工作必須自行另外撰寫程式碼運作。這個程式碼很簡單，以下是我實作的矩陣相乘︰

#include <gsl/gsl_matrix.h>
#include <assert.h>

void gsl_matrix_mul_matrix(const gsl_matrix * a,
                      const gsl_matrix * b,
                      gsl_matrix * c)
{
assert(a->size2 == b->size1
   && a->size1 == c->size1
   && b->size2 == c->size2);

/* set all elements in c to zero */
gsl_matrix_set_zero(c);

/* calculate mul */
int i, j, k;
for(i = 0; i < a->size1; i++)
{
   for(j = 0; j < b->size2; j++)
   {
       for(k = 0; k < a->size2; k++)
       {
           gsl_matrix_set(c, i, j, gsl_matrix_get(c, i, j)
               + gsl_matrix_get(a, i, k)
               * gsl_matrix_get(b, k, j));
       }
   }
}
}

5.使用完陣列記得釋放記憶體

/* 釋放掉為矩陣 m1 配置的記憶體 */
gsl_matrix_free(m1);

計算矩陣的反轉換

矩陣的反轉換函式定義在 Linear Algebra 一節的 LU Decomposition 小節中。假設要計算反矩陣的是 n x n 的矩陣 A，則︰

1.計算 LU Decomposition

gsl_permutation * p = gsl_permutation_alloc(n);
int sign = 1;
gsl_linalg_LU_decomp(A, p, &sign);

　此時 LU Decomposition 的結果儲存在 A 及 p 中。

2.計算原矩陣 A 的反矩陣

gsl_matrix * inverse = gsl_matrix_alloc(n, n);
gsl_linalg_LU_invert(A, p, inverse);

　此時矩陣 inverse 即為原矩陣 A 的反矩陣。然而必須注意的是，經過 gsl_linalg_LU_decomp的計算，A中儲存的已經不再是原矩陣，故使用 gsl_linalg_LU_decomp 函式前應先使用gsl_matrix_memcpy 函式備份原矩陣。

gsl_matrix * A_c = gsl_matrix_alloc(n, n);
gsl_matrix_memcpy(A_c, A);

Reference
GSL Reference Manual - HTML
Matrix Inverse
Identity Matrix