Friday, September 22, 2006

Inverse of a matrix - 使用 LU 演算法

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

使用 GSL 定義的矩陣
首先應該知道怎麼在 GSL 中使用其定義的矩陣,這部份應參考 GSL Reference ManualVectors 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

5 comments:

Anonymous said...

您好
我想請教一下
關於complex matrix的inverse
該如何求解呢?

Chinson said...
This comment has been removed by the author.
Chinson said...

您好,complex的求解應該可以參考 GSL 的 Reference :
http://www.gnu.org/software/gsl/manual/html_node/LU-Decomposition.html

GSL有提供各種相對於 complex 的函式可以使用

Anonymous said...

您好
我有看到那些Reference了
我會再試試看
謝謝您的指教

by the way~
方便跟您留個mail嗎?
這是我的信箱 s96325511@ncnu.edu.tw

真的謝謝您的回覆~^^

Chinson said...

不客氣^^
GSL我目前為止用過的函式還很少,
能幫上的忙其實很有限,呵呵

我的email: chinsonyeh@gmail.com