📐 线性代数基础¶

import numpy as np

# 向量定义
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

# 内积（点积）: 衡量相似度
dot = np.dot(a, b)  # 1*4 + 2*5 + 3*6 = 32
print(f"内积: {dot}")

# 外积: 生成矩阵
outer = np.outer(a, b)  # 3x3矩阵
print(f"外积:\n{outer}")

# 范数
l2_norm = np.linalg.norm(a)       # L2范数 = sqrt(14)
l1_norm = np.linalg.norm(a, 1)    # L1范数 = 6
linf_norm = np.linalg.norm(a, np.inf)  # L∞范数 = 3
print(f"L2范数: {l2_norm:.4f}, L1范数: {l1_norm}, L∞范数: {linf_norm}")

# 余弦相似度: AI中最常用的相似度度量
cos_sim = np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b))
print(f"余弦相似度: {cos_sim:.4f}")

# 线性变换可视化
import matplotlib.pyplot as plt

# 原始单位圆上的点
theta = np.linspace(0, 2*np.pi, 100)
circle = np.array([np.cos(theta), np.sin(theta)])

# 变换矩阵（拉伸+旋转）
A = np.array([[2, 1], [0, 1.5]])
transformed = A @ circle

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].plot(circle[0], circle[1]); axes[0].set_title("原始")
axes[1].plot(transformed[0], transformed[1]); axes[1].set_title("变换后")
plt.show()

A = np.array([[1, 2], [3, 4]])

# 转置
print(f"转置:\n{A.T}")

# 行列式: det(A) = ad - bc
det = np.linalg.det(A)  # 1*4 - 2*3 = -2
print(f"行列式: {det:.4f}")

# 逆矩阵: AA^{-1} = I
A_inv = np.linalg.inv(A)
print(f"逆矩阵:\n{A_inv}")
print(f"验证 A @ A^-1 = I:\n{A @ A_inv}")

# 秩
rank = np.linalg.matrix_rank(A)
print(f"秩: {rank}")

# 迹: tr(A) = 对角元素之和 = 特征值之和
trace = np.trace(A)  # 1+4 = 5
print(f"迹: {trace}")

A = np.array([[4, 2], [1, 3]])

# 特征分解
eigenvalues, eigenvectors = np.linalg.eig(A)
print(f"特征值: {eigenvalues}")      # [5, 2]
print(f"特征向量:\n{eigenvectors}")

# 验证: Av = λv
for i in range(len(eigenvalues)):
    v = eigenvectors[:, i]
    lam = eigenvalues[i]
    print(f"Av = {A @ v}, λv = {lam * v}")  # 应该相等

# 对称矩阵特征分解
S = np.array([[3, 1], [1, 3]])
eigenvalues, Q = np.linalg.eigh(S)  # eigh专用于对称矩阵
Lambda = np.diag(eigenvalues)

# 验证 S = QΛQ^T
S_reconstructed = Q @ Lambda @ Q.T
print(f"重构误差: {np.max(np.abs(S - S_reconstructed)):.2e}")

# 用SVD压缩图像
def svd_compress(image_matrix, k):
    """保留前k个奇异值"""
    U, s, Vt = np.linalg.svd(image_matrix, full_matrices=False)
    # 截断
    U_k = U[:, :k]
    s_k = np.diag(s[:k])
    Vt_k = Vt[:k, :]
    return U_k @ s_k @ Vt_k

# 原始: m*n个存储, 压缩后: k*(m+n+1)个存储
m, n = 512, 512
for k in [5, 20, 50, 100]:
    ratio = k * (m + n + 1) / (m * n)
    print(f"k={k}: 压缩率 {ratio:.2%}")

def pca_from_scratch(X, n_components):
    """手写PCA实现"""
    # 1. 中心化
    X_centered = X - X.mean(axis=0)

    # 2. 计算协方差矩阵
    cov_matrix = np.cov(X_centered.T)

    # 3. 特征分解
    eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)

    # 4. 按特征值降序排列
    idx = np.argsort(eigenvalues)[::-1]
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]

    # 5. 选取前n_components个主成分
    W = eigenvectors[:, :n_components]

    # 6. 投影
    X_projected = X_centered @ W

    # 方差解释比例
    explained_var_ratio = eigenvalues[:n_components] / eigenvalues.sum()

    return X_projected, explained_var_ratio

# 测试
np.random.seed(42)
X = np.random.randn(100, 5)  # 100个样本，5个特征
X_pca, ratios = pca_from_scratch(X, 2)
print(f"前2个主成分解释方差比例: {ratios}")
print(f"投影后形状: {X_pca.shape}")  # (100, 2)

import torch
import torch.nn.functional as F

def scaled_dot_product_attention(Q, K, V, mask=None):
    """
    Q: [batch, heads, seq_len, d_k]
    K: [batch, heads, seq_len, d_k]
    V: [batch, heads, seq_len, d_v]
    """
    d_k = Q.size(-1)

    # QK^T: [batch, heads, seq_len, seq_len]
    scores = torch.matmul(Q, K.transpose(-2, -1)) / (d_k ** 0.5)

    if mask is not None:
        scores = scores.masked_fill(mask == 0, -1e9)

    # softmax归一化
    attention_weights = F.softmax(scores, dim=-1)

    # 加权求和: [batch, heads, seq_len, d_v]
    output = torch.matmul(attention_weights, V)

    return output, attention_weights

# 示例
batch, heads, seq_len, d_k = 2, 8, 10, 64
Q = torch.randn(batch, heads, seq_len, d_k)
K = torch.randn(batch, heads, seq_len, d_k)
V = torch.randn(batch, heads, seq_len, d_k)
output, weights = scaled_dot_product_attention(Q, K, V)
print(f"输出形状: {output.shape}")     # [2, 8, 10, 64]
print(f"注意力权重形状: {weights.shape}") # [2, 8, 10, 10]

# Word2Vec经典类比: king - man + woman ≈ queen
# 这依赖于嵌入空间中的线性关系
embeddings = {
    "king": np.array([0.5, 0.8, 0.1]),
    "queen": np.array([0.5, 0.1, 0.8]),
    "man": np.array([0.2, 0.7, 0.1]),
    "woman": np.array([0.2, 0.1, 0.7]),
}

result = embeddings["king"] - embeddings["man"] + embeddings["woman"]
print(f"king - man + woman = {result}")
print(f"queen = {embeddings['queen']}")
# 两者应该非常接近

import numpy as np
from scipy.sparse.linalg import svds

def svd_recommend(ratings_matrix, k=20, top_n=10):
    """
    基于截断SVD的推荐系统

    Args:
        ratings_matrix: 用户-物品评分矩阵 (m x n)，0表示未评分
        k: 隐因子维度
        top_n: 推荐数量
    """
    # 用户均值归一化（处理评分偏置，只对已评分项取均值）
    mask = ratings_matrix != 0
    user_ratings_mean = np.where(
        mask.sum(axis=1, keepdims=True) > 0,
        (ratings_matrix * mask).sum(axis=1, keepdims=True) / mask.sum(axis=1, keepdims=True),
        0
    )
    R_normalized = (ratings_matrix - user_ratings_mean) * mask  # 未评分项保持为0

    # 截断SVD分解
    U, sigma, Vt = svds(R_normalized, k=k)
    sigma = np.diag(sigma)

    # 预测评分矩阵
    predicted_ratings = U @ sigma @ Vt + user_ratings_mean

    # 为每个用户生成推荐
    recommendations = {}
    for user_idx in range(ratings_matrix.shape[0]):
        # 找到用户未评分的物品
        unrated = np.where(ratings_matrix[user_idx, :] == 0)[0]
        # 按预测评分排序
        pred_scores = predicted_ratings[user_idx, unrated]
        top_items = unrated[np.argsort(pred_scores)[::-1][:top_n]]
        recommendations[user_idx] = top_items

    return predicted_ratings, recommendations

# 示例
np.random.seed(42)
R = np.random.randint(0, 6, (50, 30)).astype(float)  # 50用户30物品
R[R < 2] = 0  # 模拟稀疏矩阵

pred, recs = svd_recommend(R, k=5, top_n=5)
print(f"用户0的Top-5推荐物品: {recs[0]}")
print(f"预测评分: {pred[0, recs[0]].round(2)}")

import torch
import torch.nn as nn

class LoRALinear(nn.Module):  # 继承nn.Module定义神经网络层
    """LoRA线性层 - SVD低秩思想的直接应用"""
    def __init__(self, in_features, out_features, rank=4, alpha=1.0):  # __init__构造方法，创建对象时自动调用
        super().__init__()  # super()调用父类方法
        self.linear = nn.Linear(in_features, out_features, bias=False)
        self.linear.weight.requires_grad = False  # 冻结原始权重

        # 低秩分解 ΔW = B @ A
        self.lora_A = nn.Parameter(torch.randn(rank, in_features) * 0.01)
        self.lora_B = nn.Parameter(torch.zeros(out_features, rank))
        self.scale = alpha / rank

    def forward(self, x):
        # W @ x + scale * (B @ A) @ x
        return self.linear(x) + self.scale * (x @ self.lora_A.T @ self.lora_B.T)

# LoRA参数量对比
d, k, r = 4096, 4096, 8
print(f"全量参数: {d * k:,} = {d*k/1e6:.1f}M")
print(f"LoRA参数: {d * r + r * k:,} = {(d*r + r*k)/1e6:.3f}M")
print(f"压缩比: {d * k / (d * r + r * k):.0f}x")

import torch
import torch.nn.functional as F

def attention_with_gradient_check(Q, K, V):
    """验证注意力机制的梯度"""
    Q.requires_grad_(True)
    K.requires_grad_(True)
    V.requires_grad_(True)

    d_k = Q.size(-1)
    scores = torch.matmul(Q, K.transpose(-2, -1)) / (d_k ** 0.5)
    attn_weights = F.softmax(scores, dim=-1)
    output = torch.matmul(attn_weights, V)

    # 反向传播
    loss = output.sum()
    loss.backward()  # 反向传播计算梯度

    print(f"Q梯度范围: [{Q.grad.min():.4f}, {Q.grad.max():.4f}]")
    print(f"K梯度范围: [{K.grad.min():.4f}, {K.grad.max():.4f}]")
    print(f"V梯度范围: [{V.grad.min():.4f}, {V.grad.max():.4f}]")
    print(f"注意力权重熵: {-(attn_weights * attn_weights.log()).sum(-1).mean():.4f}")

    return output, attn_weights

# 测试
Q = torch.randn(2, 4, 64)  # (batch, seq_len, d_k)
K = torch.randn(2, 6, 64)
V = torch.randn(2, 6, 64)
out, weights = attention_with_gradient_check(Q, K, V)
print(f"输出形状: {out.shape}, 权重形状: {weights.shape}")

import numpy as np

# 陷阱1: 大数吃小数
a = 1e16
b = 1.0
print(f"(a + b) - a = {(a + b) - a}")  # 期望1.0，实际0.0

# 陷阱2: 矩阵条件数
A = np.array([[1, 1], [1, 1.0001]])
print(f"条件数: {np.linalg.cond(A):.0f}")  # 极大 → 求逆不稳定
b = np.array([2, 2.0001])
x = np.linalg.solve(A, b)
print(f"解: {x}")  # 解对b的微小扰动极其敏感

# 安全的log-sum-exp
def stable_logsumexp(x):
    """数值稳定的log(sum(exp(x)))"""
    c = np.max(x)
    return c + np.log(np.sum(np.exp(x - c)))

x = np.array([1000, 1001, 1002])
print(f"不稳定: {np.log(np.sum(np.exp(x)))}")  # inf
print(f"稳定版: {stable_logsumexp(x):.4f}")  # 正确结果

def cholesky_solve(A, b):
    """用Cholesky分解求解Ax=b（A正定时通常更高效）"""
    L = np.linalg.cholesky(A)
    # Ly = b (前代)
    y = np.linalg.solve(L, b)
    # L^T x = y (回代)
    x = np.linalg.solve(L.T, y)
    return x

# 高斯过程中的应用：从多维高斯采样
def sample_multivariate_gaussian(mu, cov, n_samples=5):
    """利用Cholesky从N(mu, cov)中采样"""
    L = np.linalg.cholesky(cov + 1e-6 * np.eye(len(mu)))  # 加jitter保正定
    z = np.random.randn(len(mu), n_samples)
    return mu[:, None] + L @ z

mu = np.zeros(3)
cov = np.array([[1, 0.8, 0.3], [0.8, 1, 0.5], [0.3, 0.5, 1]])
samples = sample_multivariate_gaussian(mu, cov, 1000)
print(f"采样均值: {samples.mean(axis=1).round(3)}")
print(f"采样协方差:\n{np.cov(samples).round(3)}")