矩陣因式分解（LU矩陣分解）與GSL實現，lugsl

矩陣因式分解（LU矩陣分解）與GSL實現，lugsl

   矩陣的因式分解是把一個矩陣A表示為兩個或更多個矩陣的乘積，是將複雜的資料進行分解，其中有多種方法，例如：LU分解，秩分解，QR分解，奇異值分解，譜分解等。這裡主要介紹對LU分解的認識。

根據參考的書籍，這裡的LU分解只限於一系列具有相同係數矩陣的線性方程：

Ax=b1,   Ax=b2, Ax=bp                  （1）

當A為可逆矩陣時，可計算A-1，然後計算A-1 b1，A-1 b2，等等。但是，真正在社會實踐的運用中，又是如何計算並使用的呢？ 實際而言，（1）中的第一個方程是由行變換解出的，並同時得出矩陣A的LU分解。

設A為m×n階矩陣，則Am×n可進行化簡為階梯形，此時不必行對換，那麼A可寫成形式A=LU，L是m×m下三角矩陣，主對角線元素全是1，U是A的一個等價的m×n階梯形矩陣。如下：

這樣的一個分解稱為LU分解，矩陣L是可逆的，我們稱L為單位下三角矩陣。

由上，我們可知，當A=LU時，方程Ax=b可寫成L（Ux）=b，把Ux寫成y，可以有解下面一對方程來求解x：

                  Ly=b

                  Ux=y

首先解Ly=b然後解Ux=y求得x，如下，每個方程都比較容易解，因和都是三角矩陣。下面，舉出一道例題;

例：求下列矩陣的LU分解：

因為A有4行，故L為4×4矩陣，L的計算方式為第一列是A的第一列除以它的第一行主元元素，L如下：

比較A與L的第一列。把A的第一列的後3個元素變換為零同時也為L的後三列變換，下面是A變為階梯形U:

將上述A到U的行變化結果放入L中：

故得到所求出的L和U滿足LU=A，利用LU分解，我們可以進行線性方程組的計算，簡化這種計算。後我又參考了網路上的最新資訊，得到即使矩陣無法復原，LU仍然可能存在。實際上，如果一個秩為k的矩陣的前k個順序主子式不為零，那麼它就可以進行LU分解，但反之則不然。目前，在任意域上一個方塊矩陣可進行LU分解的充要條件已經被發現，這些充要條件可以用某些特定子矩陣的秩表示。

/*********************************************************************************************************************************************************************************************************************************************************************************************************************************************************/

GNU Scientific Library(GSL)是一個為C和C++程式員提供的科學數值運算庫。該科學計算庫異常強大，提供了如下方面的支援： 
Complex Numbers               Roots of Polynomials          Special Functions Vectors and Matrices           Permutations                     Sorting BLAS Support                      Linear Algebra                    Eigensystems Fast Fourier Transforms      Quadrature                         Random Numbers Quasi-Random Sequences  Random Distributions          Statistics 
Histograms                          N-Tuples                             Monte Carlo Integration Simulated Annealing            Differential Equations         Interpolation Numerical Differentiation     Chebyshev Approximation  Series Acceleration Discrete Hankel Transforms Root-Finding                       Minimization Least-Squares Fitting          Physical Constants             IEEE Floating-Point Discrete Wavelet Transforms                                          Basis splines

對於經常處理複雜數學計算的開發人員來說，無疑是莫大的解脫

~ 

該項目的首頁是：

http://www.gnu.org/software/gsl/gsl.html

 

對於經常處理複雜數學計算的開發人員來說，無疑是莫大的解脫

~ 

該項目的首頁是：

http://www.gnu.org/software/gsl/gsl.html

 

對於經常處理複雜數學計算的開發人員來說，無疑是莫大的解脫

~ 

該項目的首頁是：

http://www.gnu.org/software/gsl/gsl.html

 

對於經常處理複雜數學計算的開發人員來說，無疑是莫大的解脫，其首頁為:

GSL---GNU Scientific Library

GSL實現：

<span style="font-family:FangSong_GB2312;font-size:18px;"><span style="font-family:KaiTi_GB2312;font-size:14px;">#include <stdio.h>#include <gsl/gsl_linalg.h>#pragma comment(lib, "libgsl_d.lib")#pragma comment(lib, "libgslcblas_d.lib")intmain (void){  double a_data[] = { 0.18, 0.60, 0.57, 0.96,                      0.41, 0.24, 0.99, 0.58,                      0.14, 0.30, 0.97, 0.66,                      0.51, 0.13, 0.19, 0.85 };  double b_data[] = { 1.0, 2.0, 3.0, 4.0 };  gsl_matrix_view m     = gsl_matrix_view_array (a_data, 4, 4);  gsl_vector_view b    = gsl_vector_view_array (b_data, 4);  gsl_vector *x = gsl_vector_alloc (4);    int s;  gsl_permutation * p = gsl_permutation_alloc (4);  gsl_linalg_LU_decomp (&m.matrix, p, &s);  gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x);  printf ("x = \n");  gsl_vector_fprintf (stdout, x, "%g");  gsl_permutation_free (p);  return 0;}</span></span>