Matlab实现Gram-Schmidt正交化方法

Gram-Schmidt正交化方法**

**
  参考了一些网络资料，利用Matlab工具实现Gram-Schmidt正交化。

1 Gram-Schmidt 理论介绍

  Gram-Schmidt 理论是线性代数中的一项基本理论。它在欧几里得空间的基础上，给出了一种将线性无关向量组正交化的方法。该方法的基本思想是：通过对向量组中的每一个向量进行处理，将其投影到前面所有向量张成的子空间上，并将其减去，使得每个向量都与前面的向量正交。

  具体来说，设 { v 1 , v 2 , ⋯   , v n } \{\mathbf{v}_1,\mathbf{v}_2,\cdots,\mathbf{v}_n\} {v1​,v2​,⋯,vn​} 是欧几里得空间 $V$ 中的一组线性无关向量，则 Gram-Schmidt 过程如下：

1. 令 ${\mathbf{u}}_1={\mathbf{v}}_1$。

2. 对于 $i=2,3,\cdots ,n$，令
   u i = v i − proj ⁡ span ⁡ { u 1 , u 2 , ⋯   , u i − 1 } v i = v i − ⟨ v i , u i − 1 ⟩ ∥ u i − 1 ∥ 2 u i − 1 − ⟨ v i , u i − 2 ⟩ ∥ u i − 2 ∥ 2 u i − 2 − ⋯ − ⟨ v i , u 1 ⟩ ∥ u 1 ∥ 2 u 1 \begin{aligned} \mathbf{u}_i&=\mathbf{v}_i-\operatorname{proj}_{\operatorname{span}\{\mathbf{u}_1,\mathbf{u}_2,\cdots,\mathbf{u}_{i-1}\}}\mathbf{v}_i \\ &=\mathbf{v}_i-\frac{\langle\mathbf{v}_i,\mathbf{u}_{i-1}\rangle}{\|\mathbf{u}_{i-1}\|^2}\mathbf{u}_{i-1}-\frac{\langle\mathbf{v}_i,\mathbf{u}_{i-2}\rangle}{\|\mathbf{u}_{i-2}\|^2}\mathbf{u}_{i-2}-\cdots-\frac{\langle\mathbf{v}_i,\mathbf{u}_1\rangle}{\|\mathbf{u}_1\|^2}\mathbf{u}_1 \end{aligned} ui​​=vi​−projspan{u1​,u2​,⋯,ui−1​}​vi​=vi​−∥ui−1​∥2⟨vi​,ui−1​⟩​ui−1​−∥ui−2​∥2⟨vi​,ui−2​⟩​ui−2​−⋯−∥u1​∥2⟨vi​,u1​⟩​u1​​

这样， { u 1 , u 2 , ⋯   , u n } \{\mathbf{u}_1,\mathbf{u}_2,\cdots,\mathbf{u}_n\} {u1​,u2​,⋯,un​} 就是一组正交向量组，且 span ⁡ { v 1 , v 2 , ⋯   , v n } = span ⁡ { u 1 , u 2 , ⋯   , u n } \operatorname{span}\{\mathbf{v}_1,\mathbf{v}_2,\cdots,\mathbf{v}_n\}=\operatorname{span}\{\mathbf{u}_1,\mathbf{u}_2,\cdots,\mathbf{u}_n\} span{v1​,v2​,⋯,vn​}=span{u1​,u2​,⋯,un​}。

  Gram-Schmidt 理论在计算机科学、物理学、数学等学科中都有广泛的应用，尤其在线性代数、矩阵论、数值分析等领域中常常被使用。

1 Matlab实现及验证

%% Gram-Schmidt 正交化matlab函数
function [ou] = Schmidt(A)
    [N,~] = size(A);
    b = zeros(N,3);
    for i = 1:N
        if i == 1
            b(i,:) = A(i,:);
            continue;
        end
        at = A(i,:);
        at2 = at;
        ind = 0;
        while 1
            ind = ind + 1;
            if ind >= i
                break
            end
            bt = b(ind,:);
            at2 = at2 - ( dot(at,bt)/dot(bt,bt)).*bt;
        end
        b(i,:) = at2;
    end
    ou = zeros(N,3);
    for i = 1:N
        bt2 = b(i,:);
        ou(i,:) = bt2/norm(bt2);
    end
end

%%测试验证
%% 测试主程序
clear;
aa = zeros(3,3);
aa(1,:) = [1,1,1];
aa(2,:) = [1,2,1];
aa(3,:) = [0,-1,1];


bb = Schmidt(aa);
bb(1,:) %out 0.5774    0.5774    0.5774
bb(2,:) %out -0.4082    0.8165   -0.4082
bb(3,:) %out -0.7071         0    0.7071