Gram-Schmidt正交化方法及其程序實現（Matlab）

Gram-Schmidt正交化方法

參考文獻：http://blog.sciencenet.cn/home.php?mod=space&uid=315774&do=blog&id=383334

問題：設有向量組α1, α2, . . . , αk，求一組與之等價的規范正交向量組v1, v2, . . . , vk（稱作原向量組的規范正交化向量組）

將上述問題轉換為下面三個問題來解釋：

（1）設v1為單位向量，alpha為v1線性尤關的任意向量.求與 v1正交的單位向量v2。(要求用v1和alpha線性表示v2）
（2）設alpha1和alpha2是兩個線性尤關的向量.求兩個正交的單位向量v1與v2:(要求用alpha1和alpha2線性表示）
（3）設alpha1,alpha2,  alpha3是三個線性尤關的向量.求三個正交的單位向量v1，v2，v3（要求用alpha1,alpha2,  alpha3線性表示）
（4）采用遞推，最后回到原問題。
具體實施步驟：

1）將alpha_j+1看作總量，並計算J已在每個v_i上的投影向量:(alpha_j+1  v_i)*v_i;

2)從從總量中減去投影向量的和，得到補向量:

3)對補向量單位化，得到

偽碼表示如下：

（貌似算法表達有一點點問題）

用Matlab描述如下：

function [v] = GS(A)
v(:,1)=A(:,1)/norm(A(:,1))
Num=0;%迭代步數
[Ahang,Alie]=size(A)  %矩陣的行和列
for j=2:Alie
    for i=1:j-1
        sum(:,1)=zeros(1,Ahang); %臨時變量 存儲后面的求和
        for k=1:j-1
            %disp('---A----')
            %disp(A(:,j));
            %disp(v(:,k));
            h(k,j)=A(:,j)'*v(:,k);
            %disp(sum(:,1))
            %disp(k)
            %disp(j)
            %disp(h(k,j))
            %disp(v(:,k))
            sum(:,1) =sum(:,1)+h(k,j)*v(:,k);
            Num=Num+1;
        end
        v(:,j)=A(:,j)-sum(:,1);
        v(:,j)=v(:,j)/norm(v(:,j));
    end
end
disp('Num=');disp(Num)
%上述某些冗余行，主要方便調試

改進版

function [v] = GS(A)
v(:,1)=A(:,1)/norm(A(:,1))
Num=0;%迭代步數
[Ahang,Alie]=size(A)  %矩陣的行和列
for j=2:Alie
    sum(:,1)=zeros(1,Ahang); %臨時變量 存儲后面的求和
    for i=1:j-1
        h(i,j)=A(:,j)'*v(:,i);
         sum(:,1) =sum(:,1)+h(i,j)*v(:,i);
            Num=Num+1;
    end
    v(:,j)=A(:,j)-sum(:,1);
    v(:,j)=v(:,j)/norm(v(:,j));
end
disp('Num=');disp(Num)

 

該程序可以實現方陣，或者是行列不等的情況。

在matlab主窗口，輸入：

>> A=[1 2 3;2 1 3]

其顯示結果為：（部分）

Num=
     4


v =

    0.4472    0.8944    0.4472
    0.8944   -0.4472    0.8944

輸入：

>> A=[1 2 3;2 1 3;1 1 2]

其顯示結果為：

Num=
     4


v =

    0.4082    0.8616    0.4082
    0.8165   -0.4924    0.8165
    0.4082    0.1231    0.4082

輸入：

>> A=[1 2 3;2 1 3;1 1 2;2 3 1]

其結果顯示如下：

Num=
     4


v =

    0.3162    0.5285    0.6454
    0.6325   -0.7047    0.6001
    0.3162   -0.0587    0.4152
    0.6325    0.4698    0.2259

可以通過norm(v(:,k))  (其中k為1~Ahang)查看是否為單位向量 或者是v(:,i)'*v(:,1)是否為1

通過v(:,i)'*v(:,j)測試是否為正交（i!=j）