Gauss Elimination算法分析與實現

高斯消去法分爲兩個過程：第一步是前向消元（forward elimination），也就是將係數矩陣轉化成上三角矩陣的過程；第二步是回代（back substitution）過程，自底向上求解方程組的過程。

選擇主元（pivoting），主元（pivot）必定不能爲0.且主元應當儘量大。因此通常狀況下，咱們採起部分選主元（partial pivoting）。選擇主元爲所在各列中絕對值最大的元素（或者非0亦可）。而後將對應的該行與主元行進行交換，而後消元。

function x = gausselimination(A,b)
ab = [A,b];
[r,c] = size(ab);
for k=1:r-1
    [rm,im]=max(abs(ab(k:r,k)));//im是最大元素對應的索引
    im = im +k-1;
    if(ab(im,k)~=0)
        if(im~=k)
            ab([k im],:)=ab([im k],:);
        end
    end
    for i=k+1:r
        ab(i,k:c)=ab(i,k:c)-ab(i,k)/ab(k,k)*ab(k,k:c);
    end
end

x=zeros(r,1);
x(r)=ab(r,c)/ab(r,r);
for i=r-1:-1:1
    x(i)=(ab(i,c)-ab(i,i+1:r)*x(i+1:r))/ab(i,i);
end

以上是matlab代碼實現（保存爲gausselimination.m）。如下是運行結果：

>> A = [10 -7 0;-3 2 6;5 -1 5]

A =

    10    -7     0

    -3     2     6

     5    -1     5

>> b = [7;4;6]

b =

     7

     4

     6

>> gausselimination(A,b)

ans =

    0.0000

   -1.0000

    1.0000

In linear algebra, Gaussian elimination (also known as row reduction) is an algorithm for solving systems of linear equations. It is usually understood as a sequence of operations performed on the associated matrix of coefficients. This method can also be used to find the rank of a matrix, to calculate the determinant of a matrix, and to calculate the inverse of an invertible square matrix. The method is named after Carl Friedrich Gauss (1777–1855), although it was known to Chinese mathematicians as early as 179 CE (see History section).

The process of row reduction makes use of elementary row operations, and can be divided into two parts. The first part (sometimes called Forward Elimination) reduces a given system to row echelon form, from which one can tell whether there are no solutions, a unique solution, or infinitely many solutions. The second part (sometimes called back substitution) continues to use row operations until the solution is found; in other words, it puts the matrix into reduced row echelon form.

算法解釋以下：

Another point of view, which turns out to be very useful to analyze the algorithm, is that row reduction produces a matrix decomposition of the original matrix. The elementary row operations may be viewed as the multiplication on the left of the original matrix by elementary matrices. Alternatively, a sequence of elementary operations that reduces a single row may be viewed as multiplication by a Frobenius matrix. Then the first part of the algorithm computes an LU decomposition, while the second part writes the original matrix as the product of a uniquely determined invertible matrix and a uniquely determined reduced row echelon matrix.

接下來一塊兒看一個小程序，它的目的不過是爲了求解線性方程組Ax+b=x,化簡以後，能夠獲得(I -A)x=b。請注意在MATLAB中的語法是怎樣的。爲此，首先咱們初始化矩陣A和列向量b

>> A=[0.85 0.04;-0.04 0.85]

A =

    0.8500    0.0400

   -0.0400    0.8500

>> b=[0;1.6]

b =

         0

    1.6000

注意觀察接下來的操做

>> I = eye(2,2)

I =

     1     0
     0     1

>> x=(I-A)\b

x =

    2.6556
    9.9585

這樣就解得了這個方程的解。最關鍵的是這個符號「\」，那麼咱們不由好奇，這個符號的功能是怎樣實現的。

\ --「It’s amazing how much numerical linear algebra is contained in that single character」

爲此介紹一個函數bslashtx

function x = bslashtx(A,b)
% BSLASHTX  Solve linear system (backslash)
% x = bslashtx(A,b) solves A*x = b


[n,n] = size(A);
if isequal(triu(A,1),zeros(n,n))
   % Lower triangular
   x = forward(A,b);
   return
elseif isequal(tril(A,-1),zeros(n,n))
   % Upper triangular
   x = backsubs(A,b);
   return
elseif isequal(A,A')
   [R,fail] = chol(A);
   if ~fail
      % Positive definite
      y = forward(R',b);
      x = backsubs(R,y);
      return
   end
end

% Triangular factorization
[L,U,p] = lutx(A);

% Permutation and forward elimination
y = forward(L,b(p));

% Back substitution
x = backsubs(U,y);


% ------------------------------

function x = forward(L,x)
% FORWARD. Forward elimination.
% For lower triangular L, x = forward(L,b) solves L*x = b.
[n,n] = size(L);
x(1) = x(1)/L(1,1);
for k = 2:n
   j = 1:k-1;
   x(k) = (x(k) - L(k,j)*x(j))/L(k,k);
end


% ------------------------------

function x = backsubs(U,x)
% BACKSUBS.  Back substitution.
% For upper triangular U, x = backsubs(U,b) solves U*x = b.
[n,n] = size(U);
x(n) = x(n)/U(n,n);
for k = n-1:-1:1
   j = k+1:n;
   x(k) = (x(k) - U(k,j)*x(j))/U(k,k);
end