Gram-Schmidt Orthogonalization

chimez

2020-12-19 16:11

Gram-Schmidt 正交化方法¶

把一组向量正交化的常用方法就是 Gram-Schmidt 正交化过程。这个过程基于投影 $$ \mathrm{proj}_u(v) = \frac{\langle u,v \rangle}{\langle u , u \rangle} u $$ 每一次迭代得到一个与之前所有矢量都正交的新矢量 $$ \begin{array}{lll} b_1 &=& a_1 \\ b_2 &=& a_2 - \mathrm{proj}_{b_1}(a_2) \\ b_3 &=& a_3 - \mathrm{proj}_{b_1}(a_3) - \mathrm{proj}_{b_2}(a_3) \\ b_k &=& a_k - \sum_{i=1}^{k-1}\mathrm{proj}_{b_i}(a_k) \end{array} $$

In [1]:

import numpy as np

In [2]:

def proj(u, v, do_norm=False):
    if do_norm:
        result = np.vdot(u,v) / np.linalg.norm(u,u) * u
    else:
        result = np.vdot(u,v) * u
    return result.flatten()

In [3]:

def gram_schmidt(A):
    b = A.copy()
    b[:,0] = b[:,0] / np.linalg.norm(b[:,0])
    for i in range(1, A.shape[1]):
        b[:,i] = (b[:,i] 
                  - np.sum(np.apply_along_axis(
                      lambda u: proj(u, b[:,i]), 
                      0, 
                      b[:,:i]),
                           axis=1).flatten())
        b[:,i] = b[:,i] / np.linalg.norm(b[:,i])
    return b

In [4]:

a = np.random.rand(4,4)
b = gram_schmidt(a)

In [5]:

from itertools import combinations
for i,j in combinations(range(4),2):
    print(np.isclose(b[:,i] @ b[:,j], 0))

True
True
True
True
True
True

In [ ]:

Gram-Schmidt 正交化方法¶

评论