如何通过正交变换将二次型化为标准型？

Question

在利用正交变换将二次型化为标准型时，一个常见技术问题是：如何确保所求的正交矩阵由实对称矩阵的特征向量构成且保持正交归一性？实践中，当特征值重数大于1时，需验证对应特征向量是否线性无关并进行施密特正交化，否则可能导致变换矩阵不正交，从而无法实现合同对角化。此外，忽略特征向量归一化步骤也会使最终变换失效。

Answer 1

1. 问题引入：正交变换与二次型标准型的数学基础

在IT领域，尤其是机器学习、计算机图形学和优化算法中，二次型的标准化处理是一个核心代数操作。其本质是通过正交变换将实对称矩阵 \( A \) 对角化，即寻找一个正交矩阵 \( P \)，使得：

\[ P^TAP = D \]

其中 \( D \) 为对角矩阵，包含 \( A \) 的特征值。这一过程依赖于谱定理——任意实对称矩阵可被正交对角化。然而，在实际计算中，若不严格保证特征向量的正交性与归一性，会导致变换失效。

2. 常见技术问题分析

• 特征向量非正交：当特征值重数大于1时，对应特征空间的基可能线性无关但不正交。
• 未归一化向量：即使正交，若未单位化，构造的矩阵 \( P \) 不满足 \( P^TP = I \)。
• 数值稳定性差：浮点运算误差累积可能导致正交性退化。
• 施密特正交化遗漏：多维特征子空间中未执行Gram-Schmidt过程。
• 特征向量顺序错乱：影响最终标准型中特征值排列一致性。

3. 数学原理与实现流程（含Mermaid流程图）

1. 输入实对称矩阵 \( A \in \mathbb{R}^{n \times n} \)
2. 求解特征值：解特征方程 \( \det(A - \lambda I) = 0 \)
3. 对每个特征值 \( \lambda_i \)，求解齐次线性方程组 \( (A - \lambda_i I)\mathbf{x} = 0 \)
4. 若代数重数 > 1，则检查几何重数是否匹配，并提取线性无关特征向量
5. 对重特征值对应的特征向量集执行施密特正交化
6. 将所有正交向量单位化（归一化）
7. 按特征值顺序排列特征向量构成矩阵 \( P \)
8. 验证 \( P^TP \approx I \) 和 \( P^TAP \approx D \)

import numpy as np

def orthogonal_diagonalize(A):
    # 验证实对称性
    assert np.allclose(A, A.T), "Matrix must be symmetric"
    
    # 求特征值与特征向量
    eigenvals, eigenvecs = np.linalg.eigh(A)
    
    # eigh 返回已正交归一化的特征向量（推荐使用）
    P = eigenvecs
    D = np.diag(eigenvals)
    
    # 验证正交性
    assert np.allclose(P.T @ P, np.eye(P.shape[1])), "Orthogonality failed"
    assert np.allclose(P.T @ A @ P, D), "Diagonalization failed"
    
    return P, D

graph TD
    A[开始] --> B{输入实对称矩阵A}
    B --> C[计算特征值与特征向量]
    C --> D{是否存在重特征值?}
    D -- 是 --> E[对该特征子空间进行施密特正交化]
    D -- 否 --> F[继续]
    E --> G[所有特征向量归一化]
    F --> G
    G --> H[构造正交矩阵P]
    H --> I[验证P^T P = I 且 P^T A P = D]
    I --> J[输出P和D]

4. 施密特正交化详解与代码实现

步骤	公式	说明
1	\( \mathbf{u}_1 = \mathbf{v}_1 \)	保留第一个向量
2	\( \mathbf{u}_2 = \mathbf{v}_2 - \frac{\langle\mathbf{v}_2,\mathbf{u}_1\rangle}{\|\mathbf{u}_1\|^2}\mathbf{u}_1 \)	减去投影分量
3	\( \mathbf{u}_k = \mathbf{v}_k - \sum_{j=1}^{k-1} \frac{\langle\mathbf{v}_k,\mathbf{u}_j\rangle}{\|\mathbf{u}_j\|^2}\mathbf{u}_j \)	递推正交化
4	\( \mathbf{e}_k = \frac{\mathbf{u}_k}{\|\mathbf{u}_k\|} \)	单位化

def gram_schmidt(vectors):
    U = []
    for v in vectors:
        u = v.copy()
        for ui in U:
            proj = np.dot(v, ui) / np.dot(ui, ui) * ui
            u -= proj
        if np.linalg.norm(u) > 1e-10:  # 排除零向量
            u = u / np.linalg.norm(u)
            U.append(u)
    return np.array(U).T