矩陣因式分解（LU矩陣分解）與GSL實作

   矩陣的因式分解是把一個矩陣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")

int
main (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>