机器学习——无监督学习中的降维算法(PCA，流形学习)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
import numpy as np
from numpy.linalg import svd, eig

class PCA:
    """主成分分析实现"""
    
    def __init__(self, n_components=None, explained_variance_ratio=None):
        """
        Parameters:
        -----------
        n_components : int, optional
            降维后的维数
        explained_variance_ratio : float, optional
            保留的累积方差比例 (0, 1)
            如果指定，则自动计算 n_components
        """
        self.n_components = n_components
        self.explained_variance_ratio_threshold = explained_variance_ratio
        self.mean_ = None
        self.components_ = None
        self.explained_variance_ = None
        self.explained_variance_ratio_ = None
        
    def fit(self, X):
        """
        拟合 PCA 模型
        
        Parameters:
        -----------
        X : array-like, shape (n_samples, n_features)
        """
        n_samples, n_features = X.shape
        
        # 中心化
        self.mean_ = np.mean(X, axis=0)
        X_centered = X - self.mean_
        
        # 方法：SVD
        U, S, Vt = svd(X_centered, full_matrices=False)
        
        # 方差（特征值）
        self.explained_variance_ = (S ** 2) / (n_samples - 1)
        
        # 方差比例
        total_variance = np.sum(self.explained_variance_)
        self.explained_variance_ratio_ = self.explained_variance_ / total_variance
        
        # 确定降维维数
        if self.explained_variance_ratio_threshold is not None:
            cumsum = np.cumsum(self.explained_variance_ratio_)
            self.n_components = np.argmax(cumsum >= self.explained_variance_ratio_threshold) + 1
        
        if self.n_components is None:
            self.n_components = min(n_samples, n_features)
        
        # 主成分（特征向量）
        self.components_ = Vt[:self.n_components, :].T
        
        return self
    
    def transform(self, X):
        """
        降维投影
        
        Parameters:
        -----------
        X : array-like, shape (n_samples, n_features)
        
        Returns:
        --------
        X_transformed : array, shape (n_samples, n_components)
        """
        X_centered = X - self.mean_
        return X_centered @ self.components_
    
    def inverse_transform(self, X_transformed):
        """
        重建原始维数
        
        Parameters:
        -----------
        X_transformed : array, shape (n_samples, n_components)
        
        Returns:
        --------
        X_original : array, shape (n_samples, n_features)
        """
        return X_transformed @ self.components_.T + self.mean_
    
    def fit_transform(self, X):
        """拟合并转换"""
        return self.fit(X).transform(X)
    
    def get_info(self):
        """获取 PCA 信息"""
        cumsum_ratio = np.cumsum(self.explained_variance_ratio_)
        print(f"选择的主成分数: {self.n_components}")
        print(f"累积解释方差: {cumsum_ratio[self.n_components-1]:.4f}")
        print(f"\n前 10 个主成分的解释方差:")
        for i in range(min(10, self.n_components)):
            print(f"  PC{i+1}: {self.explained_variance_ratio_[i]:.4f} "
                  f"(累积: {cumsum_ratio[i]:.4f})")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
class KernelPCA:
    """核主成分分析"""
    
    def __init__(self, n_components=2, kernel='rbf', gamma=None, c=None, degree=3):
        """
        Parameters:
        -----------
        kernel : str
            'linear', 'rbf', 'poly', 'sigmoid'
        gamma : float
            RBF 核的参数
        """
        self.n_components = n_components
        self.kernel = kernel
        self.gamma = gamma
        self.c = c
        self.degree = degree
        self.X_fit = None
        self.K_fit = None
        self.alphas_ = None
        self.lambdas_ = None
        
    def _kernel_matrix(self, X, Y=None):
        """计算核矩阵"""
        if Y is None:
            Y = X
        
        if self.kernel == 'linear':
            K = X @ Y.T
        
        elif self.kernel == 'rbf':
            # ||x-y||² = ||x||² + ||y||² - 2<x,y>
            X_norm_sq = np.sum(X**2, axis=1, keepdims=True)
            Y_norm_sq = np.sum(Y**2, axis=1, keepdims=True).T
            K = X_norm_sq + Y_norm_sq - 2 * (X @ Y.T)
            gamma = self.gamma if self.gamma is not None else 1.0 / X.shape[1]
            K = np.exp(-gamma * K)
        
        elif self.kernel == 'poly':
            K = (X @ Y.T + 1) ** self.degree
        
        elif self.kernel == 'sigmoid':
            K = np.tanh(self.c * (X @ Y.T) + 1)
        
        return K
    
    def fit(self, X):
        """拟合 KPCA"""
        n_samples = X.shape[0]
        self.X_fit = X
        
        # 计算核矩阵
        K = self._kernel_matrix(X)
        
        # 中心化核矩阵
        one_n = np.ones((n_samples, n_samples)) / n_samples
        K_centered = K - one_n @ K - K @ one_n + one_n @ K @ one_n
        
        # 特征值分解
        eigenvalues, eigenvectors = eig(K_centered)
        
        # 排序（降序）
        idx = np.argsort(eigenvalues)[::-1]
        self.lambdas_ = eigenvalues[idx]
        self.alphas_ = eigenvectors[:, idx[:self.n_components]]
        
        # 正交化
        for i in range(self.n_components):
            self.alphas_[:, i] /= np.sqrt(self.lambdas_[i])
        
        self.K_fit = K_centered
        
        return self
    
    def transform(self, X):
        """降维投影"""
        K = self._kernel_matrix(X, self.X_fit)
        
        # 中心化
        n_samples_fit = self.X_fit.shape[0]
        one_m_n = np.ones((X.shape[0], n_samples_fit)) / n_samples_fit
        K_centered = K - one_m_n @ self.K_fit - K @ np.ones((n_samples_fit, n_samples_fit)) / n_samples_fit
        
        return K_centered @ self.alphas_
    
    def fit_transform(self, X):
        """拟合并转换"""
        self.fit(X)
        return self.K_fit @ self.alphas_

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
import cv2
import numpy as np
from sklearn.metrics.pairwise import euclidean_distances
from pathlib import Path

class FaceRecognitionSystem:
    """人脸识别系统"""
    
    def __init__(self, pca_components=1024, use_kpca=False):
        self.pca = PCA(n_components=pca_components)
        self.kpca = KernelPCA(n_components=512, kernel='rbf') if use_kpca else None
        self.face_features = []
        self.face_names = []
        self.face_ids = []
        self.verification_threshold = 0.6
        
    def load_features(self, features_file):
        """
        加载预提取的特征
        
        Features 应该来自 ResNet-50 倒数第二层
        """
        features = np.load(features_file)  # shape: (n_faces, 4096)
        names = np.load(features_file.replace('features', 'names'))
        
        # 拟合 PCA
        self.pca.fit(features)
        reduced_features = self.pca.transform(features)
        
        # 可选：拟合 KPCA
        if self.kpca is not None:
            reduced_features = self.kpca.fit_transform(reduced_features)
        
        self.face_features = reduced_features
        self.face_names = names
        self.face_ids = np.arange(len(names))
        
        print(f"已加载 {len(names)} 个人脸")
        print(f"特征维数: {reduced_features.shape[1]}")
        
    def verify_face(self, query_features):
        """
        人脸验证 (1:1 匹配)
        
        Returns:
        --------
        (is_match, confidence, distance)
        """
        # 特征降维
        query_reduced = self.pca.transform(query_features[np.newaxis, :])[0]
        
        if self.kpca is not None:
            query_reduced = self.kpca.transform(query_reduced[np.newaxis, :])[0]
        
        # 计算距离
        distances = euclidean_distances(
            query_reduced.reshape(1, -1),
            self.face_features.reshape(len(self.face_features), -1)
        )[0]
        
        min_distance = np.min(distances)
        is_match = min_distance < self.verification_threshold
        confidence = 1.0 - (min_distance / 2.0)  # 转换为置信度
        
        return is_match, confidence, min_distance
    
    def identify_face(self, query_features, top_k=5):
        """
        人脸识别 (1:N 搜索)
        
        Returns:
        --------
        [(name, distance, confidence), ...]
        """
        # 特征降维
        query_reduced = self.pca.transform(query_features[np.newaxis, :])[0]
        
        if self.kpca is not None:
            query_reduced = self.kpca.transform(query_reduced[np.newaxis, :])[0]
        
        # 计算距离
        distances = euclidean_distances(
            query_reduced.reshape(1, -1),
            self.face_features.reshape(len(self.face_features), -1)
        )[0]
        
        # Top-K
        top_indices = np.argsort(distances)[:top_k]
        results = [
            (self.face_names[i], distances[i], 1.0 - distances[i] / 2.0)
            for i in top_indices
        ]
        
        return results
    
    def set_verification_threshold(self, threshold):
        """设置验证阈值"""
        self.verification_threshold = threshold

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
# 示例：合成数据
np.random.seed(42)

# 生成高维数据（瑞士卷）
def swiss_roll(n_samples=1000, noise=0.1):
    t = 3 * np.pi * np.random.rand(n_samples)
    height = 21 * np.random.rand(n_samples)
    x = t * np.cos(t)
    y = height
    z = t * np.sin(t)
    
    if noise:
        x += noise * np.random.randn(n_samples)
        y += noise * np.random.randn(n_samples)
        z += noise * np.random.randn(n_samples)
    
    return np.column_stack([x, y, z])

# 生成数据
X = swiss_roll(n_samples=1000, noise=0.1)

# PCA 降维
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X)

# KPCA 降维
kpca = KernelPCA(n_components=2, kernel='rbf', gamma=0.1)
X_kpca = kpca.fit_transform(X)

# 可视化
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

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

# 原始数据
ax1 = fig.add_subplot(131, projection='3d')
ax1.scatter(X[:, 0], X[:, 1], X[:, 2], c=X[:, 0], cmap='viridis', alpha=0.6)
ax1.set_title('原始数据 (瑞士卷)')
ax1.set_xlabel('X')
ax1.set_ylabel('Y')
ax1.set_zlabel('Z')

# PCA 结果
ax2 = fig.add_subplot(132)
ax2.scatter(X_pca[:, 0], X_pca[:, 1], c=X[:, 0], cmap='viridis', alpha=0.6)
ax2.set_title('PCA 降维结果')
ax2.set_xlabel('PC1')
ax2.set_ylabel('PC2')

# KPCA 结果
ax3 = fig.add_subplot(133)
ax3.scatter(X_kpca[:, 0], X_kpca[:, 1], c=X[:, 0], cmap='viridis', alpha=0.6)
ax3.set_title('KPCA 降维结果 (RBF)')
ax3.set_xlabel('KPCA1')
ax3.set_ylabel('KPCA2')

plt.tight_layout()
plt.savefig('dimensionality_reduction_comparison.png', dpi=150)
plt.show()

# PCA 信息
pca_full = PCA(n_components=3)
pca_full.fit(X)
pca_full.get_info()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
import time

def benchmark_algorithms(X, dimensions=[2, 5, 10, 20]):
    """基准测试不同降维算法"""
    
    results = {}
    
    for d in dimensions:
        print(f"\n降维到 {d} 维:")
        
        # PCA
        start = time.time()
        pca = PCA(n_components=d)
        X_pca = pca.fit_transform(X)
        pca_time = time.time() - start
        print(f"  PCA: {pca_time:.4f}s")
        
        # KPCA (RBF)
        start = time.time()
        kpca = KernelPCA(n_components=d, kernel='rbf')
        X_kpca = kpca.fit_transform(X)
        kpca_time = time.time() - start
        print(f"  KPCA: {kpca_time:.4f}s")
        
        # LLE
        start = time.time()
        lle = LLE(n_components=d)
        X_lle = lle.fit_transform(X)
        lle_time = time.time() - start
        print(f"  LLE: {lle_time:.4f}s")
        
        results[d] = {
            'pca': pca_time,
            'kpca': kpca_time,
            'lle': lle_time
        }
    
    return results