C++求矩阵eigenvalues与eigenvectors：无需库的PCA实现需求问询

Got it, let's break this down for you—since you're building PCA step-by-step and need to implement eigenvalue/eigenvector calculation from scratch (no external libs, no Matlab switch), here's a straightforward, easy-to-follow approach with actionable code and explanations tailored to PCA's needs.

Since PCA works with symmetric covariance matrices (a key detail that simplifies things), the best beginner-friendly algorithms are:

• Power Method: Great for finding the largest eigenvalue/eigenvector first (perfect for PCA's top components), easy to grasp and code.
• Power Method + Deflation: Extends the above to find subsequent eigenvalues/eigenvectors by "removing" the influence of already found components.

These avoid the complexity of full QR decomposition or other advanced methods, while aligning perfectly with what PCA actually needs (you usually care about the top N components, not all of them).

Let's walk through this with pure Python (no numpy/pandas tricks—though you can adapt this to any language):

1. First, Prep Your Data (PCA Prerequisite) #

Before calculating eigenvalues, you need to center your data (subtract the mean of each feature) and compute the covariance matrix. Here's how to do that manually:

def center_data(data):
    # data is a 2D list: rows = samples, columns = features
    num_samples = len(data)
    num_features = len(data[0])
    
    # Calculate mean for each feature
    feature_means = [sum(col) / num_samples for col in zip(*data)]
    
    # Subtract mean from each sample's feature value
    centered = []
    for sample in data:
        centered_sample = [sample[i] - feature_means[i] for i in range(num_features)]
        centered.append(centered_sample)
    return centered

def compute_covariance_matrix(centered_data):
    num_samples = len(centered_data)
    num_features = len(centered_data[0])
    
    # Initialize covariance matrix (square, size = num_features x num_features)
    cov_matrix = [[0.0 for _ in range(num_features)] for _ in range(num_features)]
    
    # Calculate covariance between each pair of features
    for i in range(num_features):
        for j in range(num_features):
            # Covariance = E[(X_i - μ_i)(X_j - μ_j)]
            cov_sum = sum(centered_data[k][i] * centered_data[k][j] for k in range(num_samples))
            cov_matrix[i][j] = cov_sum / (num_samples - 1)  # Unbiased estimate
    return cov_matrix

2. Power Method for Largest Eigenvalue/Eigenvector #

The power method works by iteratively multiplying a random vector by the covariance matrix, normalizing it, until the vector stops changing (converges). Here's the code:

def power_method(matrix, max_iterations=1000, tolerance=1e-6):
    num_rows = len(matrix)
    # Start with a random initial vector (all 1s works too, random avoids edge cases)
    eigenvec = [1.0 for _ in range(num_rows)]
    
    for _ in range(max_iterations):
        # Multiply matrix by eigenvec: new_vec = matrix * eigenvec
        new_vec = [0.0 for _ in range(num_rows)]
        for i in range(num_rows):
            new_vec[i] = sum(matrix[i][j] * eigenvec[j] for j in range(num_rows))
        
        # Normalize the new vector
        norm = sum(x**2 for x in new_vec)**0.5
        new_vec_normalized = [x / norm for x in new_vec]
        
        # Check for convergence: if change in eigenvec is below tolerance
        diff = sum(abs(new_vec_normalized[i] - eigenvec[i]) for i in range(num_rows))
        if diff < tolerance:
            break
        
        eigenvec = new_vec_normalized
    
    # Calculate eigenvalue: λ = (matrix * eigenvec) · eigenvec
    eigenvalue = sum(matrix[i][j] * eigenvec[j] * eigenvec[i] for i in range(num_rows) for j in range(num_rows))
    return eigenvalue, eigenvec

3. Extend to Multiple Eigenvalues/Eigenvectors (Deflation) #

To get the next largest eigenvalue, we "deflate" the covariance matrix by removing the contribution of the already found eigenvector. For symmetric matrices, this is simple using outer products:

def deflate_matrix(matrix, eigenvalue, eigenvec):
    num_rows = len(matrix)
    # Create outer product of eigenvec with itself
    outer_product = [[eigenvec[i] * eigenvec[j] for j in range(num_rows)] for i in range(num_rows)]
    # Deflate: subtract λ * (v * v^T) from the original matrix
    deflated = [[matrix[i][j] - eigenvalue * outer_product[i][j] for j in range(num_rows)] for i in range(num_rows)]
    return deflated

def get_top_n_eigenpairs(matrix, n):
    eigenpairs = []
    current_matrix = [row.copy() for row in matrix]  # Avoid modifying original
    
    for _ in range(n):
        val, vec = power_method(current_matrix)
        eigenpairs.append((val, vec))
        current_matrix = deflate_matrix(current_matrix, val, vec)
    
    # Sort eigenpairs in descending order of eigenvalue (since deflation can sometimes mix order)
    eigenpairs.sort(key=lambda x: x[0], reverse=True)
    return eigenpairs

4. Put It All Together for PCA #

Once you have the top N eigenpairs, you can project your centered data onto these components—this is the core of PCA. Here's a quick snippet to tie it together:

# Example usage
if __name__ == "__main__":
    # Sample data: 5 samples, 3 features
    data = [
        [1.2, 3.1, 2.0],
        [1.5, 3.2, 2.1],
        [1.7, 3.0, 1.9],
        [1.3, 2.9, 2.2],
        [1.6, 3.3, 2.0]
    ]
    
    centered = center_data(data)
    cov_matrix = compute_covariance_matrix(centered)
    
    # Get top 2 eigenpairs (for 2D PCA projection)
    top_eigenpairs = get_top_n_eigenpairs(cov_matrix, 2)
    
    print("Top Eigenvalues and Eigenvectors:")
    for idx, (val, vec) in enumerate(top_eigenpairs):
        print(f"Component {idx+1}:")
        print(f"  Eigenvalue: {val:.4f}")
        print(f"  Eigenvector: {[round(x,4) for x in vec]}")

• Convergence: The tolerance (1e-6) controls how precise you want the eigenvectors to be—adjust based on your needs.
• Symmetric Matrix Advantage: Since covariance matrices are symmetric, the power method converges faster, and deflation works reliably without extra steps.
• Normalization: We normalize eigenvectors at each step to prevent the vector values from blowing up or shrinking to zero.

内容的提问来源于stack exchange，提问作者Daniel