矩阵因式分解（LU矩阵分解）与GSL实现

来源：互联网发布：java第三方接口开发编辑：程序博客网时间：2024/06/08 06:31

矩阵的因式分解是把一个矩阵A表示为两个或更多个矩阵的乘积，是将复杂的数据进行分解，其中有多种方法，例如：LU分解，秩分解，QR分解，奇异值分解，谱分解等。这里主要介绍对LU分解的认识。

根据参考的书籍，这里的LU分解只限于一系列具有相同系数矩阵的线性方程：

Ax=b₁, Ax=b₂, Ax=b_{p （1）}

当A为可逆矩阵时，可计算A^-1，然后计算A^-1 b₁，A^-1 b₂，等等。但是，真正在社会实践的运用中，又是如何计算并使用的呢？实际而言，（1）中的第一个方程是由行变换解出的，并同时得出矩阵A的LU分解。

设A为m×n阶矩阵，则A_m_×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>

2 0