【线性代数】施密特正交化方法——Python实现

【线性代数】施密特正交化方法——Python实现

文章目录

思想Python代码

思想

施密特正交化方法：
将n维子空间中的任意一组基向量变换成标准正交向量。

假设有两个向量 a ⃗ vec{a} a 和 b ⃗ vec{b} b ，若要使两向量正交，则 a ⃗ vec{a} a 不变， b ⃗ vec{b} b 可分解为 b ⃗ vec{b} b 在 a ⃗ vec{a} a 上的投影 b ′ ⃗ vec{b'} b′ 和误差向量 e ⃗ vec{e} e 。因为 a ⃗ vec{a} a 和 e ⃗ vec{e} e 正交，所以将 a ⃗ vec{a} a 和 e ⃗ vec{e} e 标准化后，即为一组标准正交向量。

设初始矩阵为 A = [ a 1 , a 2 , a 3 , . . . , a n ] A=[a_1,a_2,a_3,...,a_n] A=[a1​,a2​,a3​,...,an​]， a i a_i ai​是 m × 1 small mtimes 1 m×1的列向量 ( i ∈ ( 1 , n ) ) (i in (1,n)) (i∈(1,n))，则 A A A是 m × n small mtimes n m×n的矩阵。
设要求的一组正交向量为 P = [ p 1 , p 2 , p 3 , . . . , p n ] P=[p_1,p_2,p_3,...,p_n] P=[p1​,p2​,p3​,...,pn​]。
首先可知 p 1 = a 1 p_1=a_1 p1​=a1​。
其次求 p 2 p_2 p2​，即求 a 2 a_2 a2​在 p 1 p_1 p1​投影的误差向量。
a 2 a_2 a2​在 p 1 p_1 p1​投影向量的投影系数 x 21 = ( p 1 T p 1 ) − 1 p 1 T a 2 x_{21}=(p_1^{T}p_1)^{-1}p_1^{T}a_2 x21​=(p1T​p1​)−1p1T​a2​。
投影向量 a 2 ′ = p 1 ( p 1 T p 1 ) − 1 p 1 T p 2 a'_2=p_1(p_1^{T}p_1)^{-1}p_1^{T}p_2 a2′​=p1​(p1T​p1​)−1p1T​p2​。
p 2 = a 2 − a 2 ′ = a 2 − p 1 ( p 1 T p 1 ) − 1 p 1 T a 2 p_2=a_2-a'_2=a_2-p_1(p_1^{T}p_1)^{-1}p_1^{T}a_2 p2​=a2​−a2′​=a2​−p1​(p1T​p1​)−1p1T​a2​。
再次求 p 2 p_2 p2​，即求 a 3 a_3 a3​在 p 1 p_1 p1​和 p 2 p_2 p2​构成的平面上投影的误差向量。
a 3 a_3 a3​在 p 1 p_1 p1​和 p 2 p_2 p2​构成的平面上投影向量可以写成 p 1 p_1 p1​和 p 2 p_2 p2​的和。
投影系数为 x 31 = ( p 1 T p 1 ) − 1 p 1 T a 3 x_{31}=(p_1^{T}p_1)^{-1}p_1^{T}a_3 x31​=(p1T​p1​)−1p1T​a3​和 x 32 = ( p 2 T p 2 ) − 1 p 2 T a 3 x_{32}=(p_2^{T}p_2)^{-1}p_2^{T}a_3 x32​=(p2T​p2​)−1p2T​a3​
投影向量 a 3 ′ = p 1 ( p 1 T p 1 ) − 1 p 1 T a 3 + p 2 ( p 2 T p 2 ) − 1 p 2 T a 3 a'_3=p_1(p_1^{T}p_1)^{-1}p_1^{T}a_3+p_2(p_2^{T}p_2)^{-1}p_2^{T}a_3 a3′​=p1​(p1T​p1​)−1p1T​a3​+p2​(p2T​p2​)−1p2T​a3​。
p 3 = a 3 − a 3 ′ = a 3 − [ p 1 ( p 1 T p 1 ) − 1 p 1 T a 3 + p 2 ( p 2 T p 2 ) − 1 p 2 T a 3 ] p_3=a_3-a'_3=a_3-[p_1(p_1^{T}p_1)^{-1}p_1^{T}a_3+p_2(p_2^{T}p_2)^{-1}p_2^{T}a_3] p3​=a3​−a3′​=a3​−[p1​(p1T​p1​)−1p1T​a3​+p2​(p2T​p2​)−1p2T​a3​]。
……
用此方法迭代即可求出这一组正交向量 P P P。
然后将其标准化。
设标准正交矩阵为 Q = [ q 1 , q 2 , q 3 , . . . , q n ] Q=[q_1,q_2,q_3,...,q_n] Q=[q1​,q2​,q3​,...,qn​]。
标准化方法为 q i = p i ∣ p i ∣ ( i ∈ ( 1 , n ) ) q_i=frac{p_i}{|p_i|}(i in (1,n)) qi​=∣pi​∣pi​​(i∈(1,n))

Python代码

import numpy as np
from scipy import linalg

def matmul_mulelms(*matrixs):
    '''
    连乘函数。将输入的矩阵按照输入顺序进行连乘。

    Parameters
    ----------
    *matrixs : 矩阵
        按计算顺序输入参数.

    Raises
    ------
    ValueError
        当参数个数小于2时，不满足乘法的要求.

    Returns
    -------
    res : 矩阵
        返回连乘的结果.

    '''
    if len(matrixs)<2:
        raise ValueError('Please input more than one parameters.')
    res = matrixs[0]
    for i in range(1,len(matrixs)):
        res = np.matmul(res, matrixs[i])
    return res

# 3.3.4 施密特正交化
def One_Col_Matrix(array):
    '''
    确保为列矩阵

    Parameters
    ----------
    array : 矩阵，向量或数组

    Raises
    ------
    ValueError
        获得的参数不是1xn或mx1时，报错.

    Returns
    -------
    TYPE
        返回列矩阵.

    '''
    mat = np.mat(array)
    if mat.shape[0] == 1:
        return mat.T
    elif mat.shape[1] == 1:
        return mat
    else:
        raise ValueError('Please input 1 row array or 1 column array')

def Transfor_Unit_Vector(matrix):
    '''
    将每列都转换为标准列向量，即模等于1

    Parameters
    ----------
    matrix : 矩阵

    Returns
    -------
    unit_mat : 矩阵
        每列模都为1的矩阵.

    '''
    col_num = matrix.shape[1]
    # 初始化为零矩阵
    unit_mat = np.zeros((matrix.shape))
    for col in range(col_num):
        vector = matrix[:,col]
        unit_vector = vector / np.linalg.norm(vector)
        unit_mat[:,col] = unit_vector.T
    return unit_mat

def Gram_Schmidt_Orthogonality(matrix):
    '''
    施密特正交化方法

    Parameters
    ----------
    matrix : 矩阵

    Returns
    -------
    标准正交化矩阵。

    '''
    col_num = matrix.shape[1]
    # 第一列无需变换
    gram_schmidt_mat = One_Col_Matrix(matrix[:,0])
    for col in range(1,col_num):
        raw_vector = One_Col_Matrix(matrix[:,col])
        orthogonal_vector = One_Col_Matrix(matrix[:,col])
        if len(gram_schmidt_mat.shape)==1:
            # 当矩阵为列向量是，shape的返回值为“(row,)”，没有col的值
            gram_schmidt_mat_col_num = 1
        else:
            gram_schmidt_mat_col_num = gram_schmidt_mat.shape[1]
        for base_vector_col in range(gram_schmidt_mat_col_num):
            base_vector = gram_schmidt_mat[:,base_vector_col]
            prejective_vector = matmul_mulelms(base_vector, linalg.inv(np.matmul(base_vector.T,base_vector)), base_vector.T, raw_vector)
            orthogonal_vector = orthogonal_vector - prejective_vector
        gram_schmidt_mat = np.hstack((gram_schmidt_mat,orthogonal_vector))
    #print(gram_schmidt_mat)
    return Transfor_Unit_Vector(gram_schmidt_mat)
### 测试用例
A = np.array([[1,1,1],
              [-1,0,-1],
              [0,-1,1]])
print(Gram_Schmidt_Orthogonality(A))