Julia语言中Gram-Schmidt正交化算法的实现方法

Gram-Schmidt正交化的作用是将矩阵的列向量组转换为两两正交的标准正交向量组，对应QR分解中Q矩阵的构造过程，以下是可适配任意维度输入矩阵的实现，分为经典版和数值稳定性更好的修正版两个版本。

核心逻辑 #

• 逐列遍历输入矩阵的原始向量
• 对当前待处理向量，减去它在所有已完成正交化的向量上的投影分量，得到和之前所有向量正交的中间向量
• 对中间向量做归一化得到标准正交基，同时记录投影系数构成上三角矩阵R
• 加入线性相关列判断，避免除零错误

经典Gram-Schmidt实现 #

function gram_schmidt(A::AbstractMatrix)
    m, n = size(A)
    # 初始化存储，自动适配浮点类型避免整数运算误差
    Q = similar(A, float(eltype(A)))
    R = zeros(float(eltype(A)), n, n)

    for j in 1:n
        v = A[:, j]
        # 减去在已正交向量上的投影
        for i in 1:j-1
            R[i, j] = dot(Q[:, i], A[:, j])
            v -= R[i, j] * Q[:, i]
        end
        R[j, j] = norm(v)
        # 处理线性相关列
        if R[j, j] < eps(real(eltype(A))) * 100
            error("输入矩阵包含线性相关列，无法完成满秩正交化")
        end
        Q[:, j] = v / R[j, j]
    end
    return Q, R
end

修正版Gram-Schmidt（数值稳定推荐） #

经典版本在向量组近似线性相关时，浮点误差会快速累积导致正交性变差，修正版调整了投影计算的顺序，稳定性大幅提升，是实际使用的首选：

function modified_gram_schmidt(A::AbstractMatrix)
    m, n = size(A)
    Q = similar(A, float(eltype(A)))
    R = zeros(float(eltype(A)), n, n)
    V = float.(A)  # 临时存储待处理向量

    for i in 1:n
        R[i, i] = norm(V[:, i])
        if R[i, i] < eps(real(eltype(A))) * 100
            error("输入矩阵包含线性相关列，无法完成满秩正交化")
        end
        Q[:, i] = V[:, i] / R[i, i]
        # 对剩余所有待处理向量，提前减去在当前正交基上的投影
        for j in i+1:n
            R[i, j] = dot(Q[:, i], V[:, j])
            V[:, j] -= R[i, j] * Q[:, i]
        end
    end
    return Q, R
end

调用示例 #

# 构造测试矩阵
A = [1.0 2.0 1.0;
     0.0 1.0 0.0;
     1.0 0.0 2.0]

Q, R = modified_gram_schmidt(A)

# 验证正交性：Q'Q 应为单位矩阵
@assert isapprox(Q' * Q, I, atol=1e-10)
# 验证分解正确性：Q*R 应等于原矩阵
@assert isapprox(Q * R, A, atol=1e-10)

# 如果不需要标准正交基，只需要正交向量组，去掉归一化步骤即可

提示：如果需要处理秩亏矩阵，可以把错误抛出逻辑改为将对应列的Q向量置零，同时记录秩的大小，适配列空间基提取的需求。

内容的提问来源于stack exchange，提问作者James Rider