向量正交的艺术

import numpy as np

# 正交向量示例
v1 = np.array([1, 0, 0])  # x轴方向
v2 = np.array([0, 1, 0])  # y轴方向
v3 = np.array([0, 0, 1])  # z轴方向

print(f"v1 · v2 = {np.dot(v1, v2)}")  # 0
print(f"v1 · v3 = {np.dot(v1, v3)}")  # 0
print(f"v2 · v3 = {np.dot(v2, v3)}")  # 0

def project_onto_orthonormal_basis(vector, basis):
    """在标准正交基上的投影"""
    projections = []
    for base_vector in basis:
        # 对于标准正交基，投影系数就是内积！
        coefficient = np.dot(vector, base_vector)
        projections.append(coefficient * base_vector)
    return projections

# 示例：将向量投影到标准正交基上
vector = np.array([3, 4, 5])
standard_basis = [
    np.array([1, 0, 0]),
    np.array([0, 1, 0]),
    np.array([0, 0, 1])
]

projections = project_onto_orthonormal_basis(vector, standard_basis)
print("投影结果:", projections)
# 结果：每个投影就是原向量的对应分量！

# 几乎线性相关的向量（病态情况）
bad_basis = np.array([
    [1.0, 1.0],
    [1.0, 1.000001]  # 几乎相同的方向
])

print("条件数:", np.linalg.cond(bad_basis))  # 非常大的数！

# 正交基
good_basis = np.array([
    [1.0, 0.0],
    [0.0, 1.0]
])

print("条件数:", np.linalg.cond(good_basis))  # 接近1！

def demonstrate_numerical_stability():
    """演示正交基的数值稳定性"""
    
    # 创建一个向量
    original = np.array([1.0, 2.0, 3.0])
    
    # 非正交基
    non_orthogonal = np.array([
        [1.0, 0.1, 0.1],
        [0.1, 1.0, 0.1],
        [0.1, 0.1, 1.0]
    ])
    
    # 正交基（单位矩阵）
    orthogonal = np.eye(3)
    
    # 多次变换后再变换回来
    result_non_orth = original.copy()
    result_orth = original.copy()
    
    for _ in range(100):
        # 非正交基的累积误差
        result_non_orth = non_orthogonal @ result_non_orth
        result_non_orth = np.linalg.inv(non_orthogonal) @ result_non_orth
        
        # 正交基的稳定性
        result_orth = orthogonal @ result_orth
        result_orth = orthogonal.T @ result_orth  # 正交矩阵的逆就是转置
    
    print(f"原向量: {original}")
    print(f"非正交基结果: {result_non_orth}")
    print(f"正交基结果: {result_orth}")
    print(f"非正交基误差: {np.linalg.norm(result_non_orth - original)}")
    print(f"正交基误差: {np.linalg.norm(result_orth - original)}")

demonstrate_numerical_stability()

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(12, 5))

# 2D正交向量
ax1 = fig.add_subplot(121)
ax1.arrow(0, 0, 1, 0, head_width=0.1, head_length=0.1, fc='red', ec='red', label='x轴')
ax1.arrow(0, 0, 0, 1, head_width=0.1, head_length=0.1, fc='blue', ec='blue', label='y轴')
ax1.set_xlim(-0.5, 1.5)
ax1.set_ylim(-0.5, 1.5)
ax1.grid(True)
ax1.set_aspect('equal')
ax1.legend()
ax1.set_title('2D正交基：完全独立的方向')

# 3D正交向量
ax2 = fig.add_subplot(122, projection='3d')
ax2.quiver(0, 0, 0, 1, 0, 0, color='red', label='x轴')
ax2.quiver(0, 0, 0, 0, 1, 0, color='blue', label='y轴')
ax2.quiver(0, 0, 0, 0, 0, 1, color='green', label='z轴')
ax2.set_xlim(0, 1)
ax2.set_ylim(0, 1)
ax2.set_zlim(0, 1)
ax2.legend()
ax2.set_title('3D正交基：三个独立方向')

plt.tight_layout()
plt.show()

def demonstrate_coordinate_system():
    """演示正交基作为坐标系的自然性"""
    
    # 任意一个点
    point = np.array([3, 4, 5])
    
    # 在标准正交基下的表示
    x_component = point[0]  # x方向的分量
    y_component = point[1]  # y方向的分量
    z_component = point[2]  # z方向的分量
    
    print(f"点 {point} 在正交基下的分解：")
    print(f"x方向: {x_component}")
    print(f"y方向: {y_component}")  
    print(f"z方向: {z_component}")
    
    # 验证：分量的平方和等于长度的平方（勾股定理！）
    length_squared = x_component**2 + y_component**2 + z_component**2
    actual_length_squared = np.dot(point, point)
    
    print(f"长度平方（勾股定理）: {length_squared}")
    print(f"实际长度平方: {actual_length_squared}")

demonstrate_coordinate_system()

def orthogonal_matrix_benefits():
    """演示正交矩阵的计算优势"""
    
    # 创建一个正交矩阵（旋转矩阵）
    theta = np.pi / 4  # 45度
    Q = np.array([
        [np.cos(theta), -np.sin(theta)],
        [np.sin(theta),  np.cos(theta)]
    ])
    
    print("正交矩阵 Q:")
    print(Q)
    
    # 验证正交性质
    print(f"\nQ @ Q.T = \n{Q @ Q.T}")  # 应该是单位矩阵
    print(f"det(Q) = {np.linalg.det(Q)}")  # 应该是1
    
    # 计算逆矩阵：对于正交矩阵，逆 = 转置
    inverse_fast = Q.T  # O(1)操作！
    inverse_slow = np.linalg.inv(Q)  # O(n^3)操作
    
    print(f"\n快速逆矩阵(转置): \n{inverse_fast}")
    print(f"传统逆矩阵计算: \n{inverse_slow}")
    print(f"差异: {np.max(np.abs(inverse_fast - inverse_slow))}")

orthogonal_matrix_benefits()

def solve_with_orthogonal_matrix():
    """使用正交矩阵快速求解线性方程组"""
    
    # 方程组 Qx = b，其中Q是正交矩阵
    Q = np.array([
        [1/np.sqrt(2), -1/np.sqrt(2)],
        [1/np.sqrt(2),  1/np.sqrt(2)]
    ])
    
    b = np.array([1, 2])
    
    # 传统方法：x = Q^(-1) * b
    x_traditional = np.linalg.solve(Q, b)
    
    # 正交矩阵快速方法：x = Q^T * b
    x_fast = Q.T @ b
    
    print(f"传统方法结果: {x_traditional}")
    print(f"快速方法结果: {x_fast}")
    print(f"差异: {np.max(np.abs(x_traditional - x_fast))}")

solve_with_orthogonal_matrix()

def rotation_matrix_parameterization():
    """演示旋转矩阵的紧凑表示"""
    
    # 2D旋转矩阵只需要1个参数（角度）
    theta = np.pi / 3
    R_2d = np.array([
        [np.cos(theta), -np.sin(theta)],
        [np.sin(theta),  np.cos(theta)]
    ])
    
    print("2D旋转矩阵（4个数字，但只有1个自由度）:")
    print(R_2d)
    
    # 3D旋转矩阵只需要3个参数（欧拉角）
    # 而不是9个矩阵元素
    print(f"\n存储效率: 用1个角度参数表示4个矩阵元素")
    print(f"压缩比: {4/1} : 1")

rotation_matrix_parameterization()

def fourier_basis_example():
    """演示傅里叶基函数的正交性"""
    
    # 创建时间序列
    t = np.linspace(0, 2*np.pi, 100)
    
    # 正交的正弦和余弦函数
    cos_1 = np.cos(t)
    sin_1 = np.sin(t)
    cos_2 = np.cos(2*t)
    sin_2 = np.sin(2*t)
    
    # 验证正交性
    print("傅里叶基函数的正交性:")
    print(f"cos(t) · sin(t) = {np.trapz(cos_1 * sin_1, t):.6f}")
    print(f"cos(t) · cos(2t) = {np.trapz(cos_1 * cos_2, t):.6f}")
    print(f"sin(t) · sin(2t) = {np.trapz(sin_1 * sin_2, t):.6f}")
    
    # 信号分解示例
    signal = 3 * cos_1 + 2 * sin_2  # 复合信号
    
    # 提取分量（利用正交性）
    cos_1_coeff = np.trapz(signal * cos_1, t) / np.trapz(cos_1 * cos_1, t)
    sin_2_coeff = np.trapz(signal * sin_2, t) / np.trapz(sin_2 * sin_2, t)
    
    print(f"\n信号分解:")
    print(f"cos(t)的系数: {cos_1_coeff:.2f} (真实值: 3)")
    print(f"sin(2t)的系数: {sin_2_coeff:.2f} (真实值: 2)")

fourier_basis_example()

from sklearn.decomposition import PCA
from sklearn.datasets import make_blobs

def pca_orthogonality_demo():
    """演示PCA中正交主成分的重要性"""
    
    # 生成示例数据
    X, _ = make_blobs(n_samples=100, centers=1, n_features=2, 
                      cluster_std=2.0, random_state=42)
    
    # 应用PCA
    pca = PCA(n_components=2)
    X_pca = pca.fit_transform(X)
    
    # 获取主成分（正交向量）
    components = pca.components_
    
    print("主成分（正交向量）:")
    print(f"第一主成分: {components[0]}")
    print(f"第二主成分: {components[1]}")
    
    # 验证正交性
    orthogonality = np.dot(components[0], components[1])
    print(f"正交性验证 (应该接近0): {orthogonality:.10f}")
    
    # 解释的方差比例
    print(f"解释的方差比例: {pca.explained_variance_ratio_}")

pca_orthogonality_demo()

def orthogonal_weight_initialization():
    """演示正交权重初始化的好处"""
    
    # 随机权重初始化
    np.random.seed(42)
    random_weights = np.random.randn(100, 100)
    
    # 正交权重初始化
    orthogonal_weights, _ = np.linalg.qr(np.random.randn(100, 100))
    
    # 计算条件数（衡量数值稳定性）
    cond_random = np.linalg.cond(random_weights)
    cond_orthogonal = np.linalg.cond(orthogonal_weights)
    
    print("权重矩阵的条件数比较:")
    print(f"随机初始化: {cond_random:.2f}")
    print(f"正交初始化: {cond_orthogonal:.2f}")
    print(f"改善比例: {cond_random/cond_orthogonal:.2f}x")

orthogonal_weight_initialization()

def error_detection_with_orthogonality():
    """利用正交性进行错误检测"""
    
    def is_orthogonal_matrix(matrix, tolerance=1e-10):
        """检查矩阵是否正交"""
        should_be_identity = matrix @ matrix.T
        identity = np.eye(matrix.shape[0])
        error = np.max(np.abs(should_be_identity - identity))
        return error < tolerance, error
    
    # 测试正确的正交矩阵
    correct_matrix = np.eye(3)
    is_orth_1, error_1 = is_orthogonal_matrix(correct_matrix)
    
    # 测试错误的矩阵
    wrong_matrix = np.array([[1, 0.1], [0, 1]])
    is_orth_2, error_2 = is_orthogonal_matrix(wrong_matrix)
    
    print("正交性检测:")
    print(f"单位矩阵: 正交={is_orth_1}, 误差={error_1:.2e}")
    print(f"错误矩阵: 正交={is_orth_2}, 误差={error_2:.2e}")

error_detection_with_orthogonality()

import concurrent.futures

def parallel_projection_demo():
    """演示正交基上并行投影计算"""
    
    # 大向量和正交基
    vector = np.random.randn(1000)
    orthogonal_basis = [
        np.random.randn(1000) for _ in range(10)
    ]
    
    # 施密特正交化确保基向量正交
    from scipy.linalg import orth
    orthogonal_basis = orth(np.column_stack(orthogonal_basis)).T
    
    def compute_projection(base_vector):
        """计算在单个基向量上的投影"""
        coefficient = np.dot(vector, base_vector)
        return coefficient * base_vector
    
    # 串行计算
    import time
    start_time = time.time()
    serial_projections = [compute_projection(base) for base in orthogonal_basis]
    serial_time = time.time() - start_time
    
    # 并行计算
    start_time = time.time()
    with concurrent.futures.ThreadPoolExecutor() as executor:
        parallel_projections = list(executor.map(compute_projection, orthogonal_basis))
    parallel_time = time.time() - start_time
    
    print(f"串行计算时间: {serial_time:.4f}秒")
    print(f"并行计算时间: {parallel_time:.4f}秒")
    print(f"加速比: {serial_time/parallel_time:.2f}x")

# parallel_projection_demo()  # 取消注释以运行