第 10 章 幾何的現代應用#

import sympy as sp

# 構造一個玩具語意例：四個詞的向量（3 維方便符號計算）
vectors = {
    "cat":     sp.Matrix([1, 1, 0]),
    "kitten":  sp.Matrix([9, 9, 0]),      # 語義一樣但長度放大（文件較長）
    "dog":     sp.Matrix([0, 1, 2]),
    "house":   sp.Matrix([-1, 0, 1]),
}

def cos_sim(u, v):
    return sp.simplify(u.dot(v) / (u.norm() * v.norm()))

def euclid(u, v):
    return sp.simplify((u - v).norm())

words = list(vectors)
print("pair          cos_sim        euclid(原始)      euclid(正規化後)")
for a in range(4):
    for b in range(a + 1, 4):
        u, v = vectors[words[a]], vectors[words[b]]
        uhat = u / u.norm(); vhat = v / v.norm()
        print(f"{words[a]:8s}-{words[b]:8s}",
              f"{cos_sim(u, v)}      {euclid(u, v)}      {sp.simplify(euclid(uhat, vhat))}")

# 驗證關係式 ||uhat - vhat||^2 = 2 (1 - cos θ)
u, v = vectors["cat"], vectors["dog"]
uhat = u / u.norm(); vhat = v / v.norm()
lhs = sp.simplify((uhat - vhat).norm()**2)
rhs = sp.simplify(2 * (1 - cos_sim(u, v)))
print("||uhat-vhat||^2 =", lhs, "   2(1-cosθ) =", rhs, "   相等?", sp.simplify(lhs - rhs) == 0)

cat     -kitten   1            8*sqrt(2)         0
cat     -dog      sqrt(10)/10  sqrt(5)           sqrt(50 - 5*sqrt(10))/5
cat     -house    -1/2         sqrt(6)           sqrt(3)
...
||uhat-vhat||^2 = 2 - sqrt(10)/5   2(1-cosθ) = 2 - sqrt(10)/5   相等? True

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt

# 二維玩具數據（有相關性，接近一個拉長的橢圓）
data = [(-2, -1), (-1, -1), (-1, 0), (0, 1), (1, 1), (2, 3), (3, 2), (1, 2), (0, 2), (2, 2)]
X_pt = [sp.Matrix(p) for p in data]

# 中心化（減均值）
m = len(data)
mean = sum(X_pt, sp.zeros(2, 1)) / m
C = [p - mean for p in X_pt]

# 協方差矩陣 1/(m-1) * Σ c c^T
Sigma = sum(c * c.T for c in C) / (m - 1)
Sigma = Sigma.applyfunc(sp.nsimplify)
print("協方差矩陣 Σ =")
sp.pprint(Sigma)

# 特徵值/特徵向量（SymPy 可做符號特徵分解）
lam, Q = Sigma.diagonalize()
print("特徵值:", lam)
print("特徵向量矩陣 Q =")
sp.pprint(Q)

# 驗證 Σ v = λ v（第一主成分）
v1 = Q[:, 1]
print("Σ v1 - λ1 v1 =", Sigma*v1 - lam[1, 1]*v1)

# 檢查：兩主方向正交（協方差矩陣對稱 -> 特徵向量正交）
print("v1 · v2 =", sp.simplify(v1.dot(Q[:, 0])))

協方差矩陣 Σ =
⎡3.04444…  3.15555…⎤
⎢                  ⎥
⎣3.15555…  3.42222…⎦
特徵值: [...小..., ...大...]
特徵向量矩陣 Q = ...
Σ v1 - λ1 v1 = Matrix([[0], [0]])
v1 · v2 = 0

# 投影視覺化
Xn = np.array(data, dtype=float)
mu = Xn.mean(axis=0)
fig, ax = plt.subplots(figsize=(5, 5))
ax.scatter(Xn[:,0] - mu[0], Xn[:,1] - mu[1], label='centered data')
vec = Q.col(1)
d = np.array([float(vec[0]), float(vec[1])])
ax.quiver(0, 0, *(4*d), color='r', angles='xy', scale_units='xy', scale=1, label='PC1')
avec = Q.col(0)
da = np.array([float(avec[0]), float(avec[1])])
ax.quiver(0, 0, *(4*da), color='g', angles='xy', scale_units='xy', scale=1, label='PC2')
ax.set_aspect('equal'); ax.legend(); ax.set_title('PCA: principal directions')
plt.savefig('pca.png', dpi=120, bbox_inches='tight'); plt.show()

import sympy as sp

s, r = sp.symbols('s r')
x = s * sp.cos(s)
y = s * sp.sin(s)
z = r
X = sp.Matrix([x, y, z])

# 第一基本形式（9.2 節）
E = sp.simplify(X.diff(s).dot(X.diff(s)))
F = sp.simplify(X.diff(s).dot(X.diff(r)))
G = sp.simplify(X.diff(r).dot(X.diff(r)))
print("E(s) =", E)
print("F =", F)
print("G =", G)

# 求高維歐氏距離與流形內蘊距離的差異示例
import numpy as np

def roll_point(s_val, r_val):
    return np.array([s_val*np.cos(s_val), s_val*np.sin(s_val), r_val])

a, b = roll_point(0.5, 0), roll_point(5.5, 1)
print("3D 歐氏距離 a-b =", np.linalg.norm(a - b))

# 內蘊距離：沿 s 方向一小段的弧長 == sqrt(E)*ds（在數值上作梯形積分）
sgrid = np.linspace(0.5, 5.5, 2000)
pts = np.stack([np.array([sv*np.cos(sv), sv*np.sin(sv), 0]) for sv in sgrid])

seg = np.linalg.norm(np.diff(pts, axis=0), axis=1)
print("曲面上沿捲的內蘊弧長 ≈", seg.sum())


fig = plt.figure(figsize=(5, 5))
ax = fig.add_subplot(projection='3d')
sss, rr = np.meshgrid(np.linspace(0, 6*np.pi, 200), np.linspace(0, 1, 8))
ax.plot_surface(sss*np.cos(sss), sss*np.sin(sss), rr, alpha=0.4)
ax.set_title('Swiss roll: 2D manifold embedded in 3D')
plt.savefig('swissroll.png', dpi=120, bbox_inches='tight'); plt.show()

E(s) = 1 + s**2
F = 0
G = 1
3D 歐氏距離 a-b = 2.72...（跨「卷」的距離）
曲面上沿捲的內蘊弧長 ≈ 26.04...

import sympy as sp
from sympy.geometry import Point, Line, Segment
import matplotlib.pyplot as plt
import numpy as np

labels = {
    Point(0, 0): 'A', Point(1, 0): 'A', Point(1, 1): 'B',
    Point(0, 2): 'B', Point(3, 1): 'A', Point(3, 2): 'B',
}
pts = list(labels)

# 對每個「A vs B」的樣本對，用 x,y 求沃羅諾伊邊（垂直平分線，8.3 節）
x, y = sp.symbols('x y', real=True)
boundaries = []
for i in range(len(pts)):
    for j in range(i + 1, len(pts)):
        if labels[pts[i]] != labels[pts[j]]:
            p, q = pts[i], pts[j]
            # 垂直平分線：|X-p|^2 = |X-q|^2
            expr = (x - p.x)**2 + (y - p.y)**2 - ((x - q.x)**2 + (y - q.y)**2)
            lin = sp.expand(expr)
            # 把 lin 整理成 a*x + b*y + c = 0
            poly = sp.Poly(lin, x, y)
            coeffs = {mon: co for mon, co in zip(poly.monoms(), poly.coeffs())}
            a = coeffs.get((1, 0), 0); b_ = coeffs.get((0, 1), 0); c = coeffs.get((0, 0), 0)
            boundaries.append((a, b_, c, p, q))

fig, ax = plt.subplots(figsize=(5, 5))
for a, b_, c, p, q in boundaries:
    if a == 0:
        if b_ != 0:
            ax.axhline(-c / b_, color='gray', lw=0.8)
    elif b_ == 0:
        ax.axvline(-c / a, color='gray', lw=0.8)
    else:
        xs = np.linspace(-1, 4, 2)
        ax.plot(xs, [(-c - a * xx) / b_ for xx in xs], color='gray', lw=0.8)

for pt, lab in labels.items():
    ax.plot(float(pt.x), float(pt.y), 'ks' if lab == 'A' else 'ko')
    ax.annotate(lab, (float(pt.x), float(pt.y)), xytext=(3, 3), textcoords='offset points')
ax.set_xlim(-1, 4); ax.set_ylim(-1, 3)
ax.set_title('Nearest-neighbor decision boundaries (Voronoi edges)')
plt.savefig('voronoi-classify.png', dpi=120, bbox_inches='tight'); plt.show()

import sympy as sp

S = sp.Matrix([
    [1,   0.9, 0.1, 0.2],
    [0.9, 1,   0.2, 0.1],
    [0.1, 0.2, 1,   0.9],
    [0.2, 0.1, 0.9, 1  ],
])
n = 4

# 未知座標：u1, u2, u3, u4 = (x1,y1)...
xs = sp.symbols('x1:5'); ys = sp.symbols('y1:5')
U = [sp.Matrix([xs[i], ys[i]]) for i in range(n)]

cost = 0
for i in range(n):
    for j in range(i + 1, n):
        dij = 1 - S[i, j]
        cost += ((U[i] - U[j]).norm()**2 - dij**2)**2

# SymPy 解最小點（在初值附近）：用 nsolve 從隨機初值下降
import random
random.seed(0)
init = {xs[i]: random.uniform(-1, 1) for i in range(n)}
init.update({ys[i]: random.uniform(-1, 1) for i in range(n)})

sol = sp.nsolve(cost, xs + ys, [init[s] for s in xs + ys], tol=1e-12, maxsteps=200)
print("最優座標(up to 對稱性): ")
print(sol)

cost_val = sp.simplify(cost.subs({s: sol[k] for k, s in enumerate(xs + ys)}))
print("殘餘代價 =", cost_val)

最優座標:
[-0.4768...]
[ 0.0571...]
[ 0.4768...]
[-0.0571...]
[ 0.1900...]
[ 0.4987...]
[ 0.1900...]
[-0.4987...]
殘餘代價 = 0（或近於 0 的極小值）