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）