用manim实现Gram-Schmidt正交化过程

news2024/10/6 18:27:43

在线性代数中，正交基有许多美丽的性质。例如，由正交列向量组成的矩阵(又称正交矩阵)可以通过矩阵的转置很容易地进行反转。此外，例如：在由彼此正交的向量张成的子空间上投影向量也更容易。Gram-Schmidt过程是一个重要的算法，它允许我们将任意基转换为生成同一子空间的正交基。在这篇文章中，我们将使用一个流行的开源库。manim在3D中实现和可视化这个算法

Gram-Schmidt过程是一种用于将一组线性无关的向量转化为一组正交（或正交归一化）的向量的算法。这个过程在数学和工程中广泛应用，特别是在计算机图形学、信号处理和统计分析中。

Gram-Schmidt 正交化的步骤

假设我们有一组线性无关的向量 ${v_{1},v_{2},...,v_{n}}$ ，Gram-Schmidt过程的步骤如下：

初始化：给定向量 $v_{1}$ ，定它为第一个正交向量 $u_{1}$ :

$u_{1}=v_{1}$

2.后续向量的正交化：对于每个 $i=1,2,3,...,n$ ，计算新的正交向量 $u_{1}$ ：

$u_{1}=v_{1}- \sum_{j=1}^{i-1}proj_{u_{i}}(v_{i})$

其中，投影 $proj_{u_{i}}(v_{i})$ 是：

$proj_{u_{i}}(v_{i})=\frac{v_{i}.u_{j}}{v_{j}.u_{j}}.u_{j}$

这里表示向量的点积。

3.归一化（可选）：如果需要将向量归一化，使其单位长度，可以通过以下公式进行：

$q_{i}=\frac{u_{i}}{||u_{i}||}$

其中 $||u_{i}||$ 是向量 $u_{i}$ 的模。

应用

正交基：Gram-Schmidt过程的输出结果是正交基，可以用于简化内积空间中的计算。
数值稳定性：在计算中使用正交向量能提高数值稳定性，特别是在进行矩阵分解（如QR分解）时。
计算简化：在许多应用中（如最小二乘法），正交向量可以简化计算，使得更容易处理高维数据。

用manim实现向量的表示

代码为了实现 Gram-Schmidt 正交化的可视化过程而构建，使用了 Manim 库。

from manim import *  # 导入 Manim 库  
import numpy as np  # 导入 NumPy 库以进行数值计算  
from enum import Enum  # 导入枚举模块以创建动作类型  

# 定义颜色  
basis_i_color = GREEN  # 基向量 i 的颜色  
basis_j_color = RED    # 基向量 j 的颜色  
basis_k_color = GOLD   # 基向量 k 的颜色  
q_color = PURPLE       # 向量 q 的颜色  
q_shifted_color = PINK # 移动的向量 q 的颜色  
projection_color = BLUE # 投影向量的颜色  

# 定义动作枚举，便于管理不同的操作  
class Action(Enum):  
    UPDATE_MATRIX_REMOVE_Q = 1  # 更新矩阵并移除 q 向量  
    ADD_PROJECTION = 2           # 添加一个投影向量  
    REMOVE_PROJECTIONS_SET_Q = 3 # 移除投影向量并设置新的 q 向量  
    NORMALIZE_Q = 4              # 规范化 q 向量  

# 实现 Gram-Schmidt 过程  
def gram_schmidt(A):  
    (n, m) = A.shape  # 获取矩阵 A 的形状  
    
    for i in range(m):  
        q = A[:, i]  # 选取矩阵 A 的第 i 列  
        
        for j in range(i):  
            # 计算并从 q 中减去之前列的投影  
            projection = np.dot(A[:, j], A[:, i]) * A[:, j]  
            q = q - projection  
        
        # 归一化 q 向量  
        q /= np.linalg.norm(q)  
        yield Action.NORMALIZE_Q, q  # 返回规范化的 q 向量  

# 创建一个新的场景，用于展示 Gram-Schmidt  
class GramSchmidt(Scene):  
    def construct(self):  
        # 创建基础向量的颜色，并初始化向量和矩阵 M  
        M = np.random.rand(3, 3)  # 随机初始化 3x3 矩阵 M  
        projection_vectors = []     # 初始化用于存储投影向量的列表  
        q = None                    # 初始化 q 向量为 None  

        # 初始化基向量箭头  
        i_vec = Vector(M[:, 0], color=basis_i_color)  # i 向量  
        j_vec = Vector(M[:, 1], color=basis_j_color)  # j 向量  
        k_vec = Vector(M[:, 2], color=basis_k_color)  # k 向量  

        # 播放箭头和矩阵的动画  
        self.play(GrowArrow(i_vec), GrowArrow(j_vec), GrowArrow(k_vec))  
        self.wait()  

        # 遍历 Gram-Schmidt 的步骤  
        for (action, payload) in gram_schmidt(M):  
            if action == Action.UPDATE_MATRIX_REMOVE_Q:  
                # 此步骤用于更新矩阵，并移除旧的 q 向量  
                assert q is not None  
                M_rounded = np.round(M.copy(), 2)  # 对矩阵 M 进行四舍五入  
                matrix = self.create_matrix(M_rounded)  # 创建新的矩阵对象  
                self.remove(matrix)  # 移除旧的矩阵  

                # 更新基向量  
                i_vec_new = Vector(M[:, 0], color=basis_i_color)  
                j_vec_new = Vector(M[:, 1], color=basis_j_color)  
                k_vec_new = Vector(M[:, 2], color=basis_k_color)  

                animation_time = 2.0  # 动画时间  
                
                # 播放更新基向量的动画  
                self.play(  
                    FadeOut(q, run_time=animation_time * 0.75),  
                    ReplacementTransform(i_vec, i_vec_new, run_time=animation_time),  
                    ReplacementTransform(j_vec, j_vec_new, run_time=animation_time),  
                    ReplacementTransform(k_vec, k_vec_new, run_time=animation_time)  
                )  
                
                self.wait()  # 等待一段时间  
                # 更新现有向量引用  
                i_vec, j_vec, k_vec = i_vec_new, j_vec_new, k_vec_new  

            elif action == Action.ADD_PROJECTION:  
                # 添加新的投影向量  
                p = Vector(payload, color=projection_color)  
                projection_vectors.append(p)  # 将投影向量添加到列表中  
                self.play(GrowArrow(p))  # 播放投影箭头的动画  

                self.wait()  # 等待  

                if len(projection_vectors) == 2:  
                    # 当有两个投影向量时，更新其显示效果  
                    first_projection_end = projection_vectors[0].get_end()  
                    p_shifted = Arrow(first_projection_end, first_projection_end + payload, buff=0, color=projection_color)  
                    projection_vectors[1] = p_shifted  

                    self.play(ReplacementTransform(p, p_shifted))  # 替换动画  
                    self.wait()  # 等待  

            elif action == Action.REMOVE_PROJECTIONS_SET_Q:  
                # 移除投影并设置新的 q 向量  
                if not projection_vectors:  
                    q = Vector(payload, color=q_color)  
                    self.play(GrowArrow(q))  # 播放 q 向量的动画  
                    self.wait()  
                else:  
                    last_projection_end = projection_vectors[-1].get_end()  
                    q_shifted = Arrow(last_projection_end, last_projection_end + payload, buff=0, color=q_shifted_color)  
                    self.play(GrowArrow(q_shifted))  
                    self.wait()  

                    q = Vector(payload, color=q_color)  # 设置新的 q 向量  
                    self.play(  
                        ReplacementTransform(q_shifted, q),  
                        *[FadeOut(p) for p in projection_vectors]  # 移除所有投影  
                    )  
                    self.wait()  
                    projection_vectors = []  # 清空投影向量列表  

            elif action == Action.NORMALIZE_Q:  
                # 规范化 q 向量  
                q_normalized = Vector(payload, color=q_color)  
                self.play(ReplacementTransform(q, q_normalized))  # 播放规范化动画  
                self.wait()  
                q = q_normalized  # 更新 q 向量  

            else:  
                assert False  # 确保没有未识别的动作  

            self.wait(1)  # 在每轮操作间等待  

        # 验证结果  
        assert np.allclose(M.T @ M, np.identity(3))  # 确保 M 的转置乘以 M 为单位矩阵  
        self.wait(15)  # 等待更长时间以查看最终结果