PCA (Principal Component Analysis) in Machine Learning

import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

# Load your feature matrix (n_samples, n_features)
X = np.random.randn(10000, 200)  # Replace with real data

# Step 1: Always standardize before PCA
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Step 2: Fit PCA with all components to analyze variance
pca_full = PCA()  # Keep all components
pca_full.fit(X_scaled)

# Step 3: Scree plot — explained variance per component
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Individual explained variance
axes[0].bar(range(1, 51), pca_full.explained_variance_ratio_[:50])
axes[0].set_xlabel('Principal Component')
axes[0].set_ylabel('Explained Variance Ratio')
axes[0].set_title('Scree Plot')

# Cumulative explained variance
cumulative_var = np.cumsum(pca_full.explained_variance_ratio_)
axes[1].plot(range(1, len(cumulative_var) + 1), cumulative_var)
axes[1].axhline(y=0.95, color='r', linestyle='--', label='95% threshold')
axes[1].set_xlabel('Number of Components')
axes[1].set_ylabel('Cumulative Explained Variance')
axes[1].set_title('Cumulative Variance Plot')
axes[1].legend()
plt.tight_layout()
plt.savefig('scree_plot.png', dpi=150)

# Step 4: Choose k where cumulative variance >= 95%
k = np.argmax(cumulative_var >= 0.95) + 1
print(f'Components needed for 95% variance: {k}')

# Step 5: Fit final PCA with selected k
pca = PCA(n_components=k)
X_reduced = pca.fit_transform(X_scaled)
print(f'Original shape: {X_scaled.shape}')
print(f'Reduced shape:  {X_reduced.shape}')
print(f'Compression ratio: {X_scaled.shape[1] / X_reduced.shape[1]:.1f}x')

import numpy as np
from sklearn.decomposition import IncrementalPCA
from sklearn.preprocessing import StandardScaler
import joblib

# Configuration
N_COMPONENTS = 50
BATCH_SIZE = 5000
N_TOTAL = 1_000_000  # Total samples (too large for memory)
N_FEATURES = 500

# Initialize incremental PCA
ipca = IncrementalPCA(n_components=N_COMPONENTS, batch_size=BATCH_SIZE)

# Simulated data generator (replace with your data loader)
def data_generator(total_samples, n_features, batch_size):
    """Yields batches from disk/database without loading all into memory."""
    for start in range(0, total_samples, batch_size):
        end = min(start + batch_size, total_samples)
        # In production: load from Parquet, HDF5, or database
        batch = np.random.randn(end - start, n_features)
        yield batch

# Pass 1: Fit the scaler incrementally
scaler = StandardScaler()
for batch in data_generator(N_TOTAL, N_FEATURES, BATCH_SIZE):
    scaler.partial_fit(batch)

# Pass 2: Fit IncrementalPCA on scaled data
for batch in data_generator(N_TOTAL, N_FEATURES, BATCH_SIZE):
    batch_scaled = scaler.transform(batch)
    ipca.partial_fit(batch_scaled)

print(f'Explained variance (top {N_COMPONENTS}): '
      f'{sum(ipca.explained_variance_ratio_):.3f}')

# Save artifacts for production serving
joblib.dump({
    'scaler': scaler,
    'ipca': ipca,
    'n_components': N_COMPONENTS,
}, 'pca_artifacts.joblib')

# Inference: transform a single batch
artifacts = joblib.load('pca_artifacts.joblib')
X_new = np.random.randn(100, N_FEATURES)
X_transformed = artifacts['ipca'].transform(
    artifacts['scaler'].transform(X_new)
)
print(f'Inference output shape: {X_transformed.shape}')

import numpy as np
from sklearn.decomposition import KernelPCA
from sklearn.datasets import make_circles
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt

# Generate nonlinear data (concentric circles)
X, y = make_circles(n_samples=2000, factor=0.3, noise=0.05, random_state=42)

# Standard PCA fails on nonlinear structure
from sklearn.decomposition import PCA
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X)

# Kernel PCA with RBF kernel captures nonlinear structure
kpca = KernelPCA(
    n_components=2,
    kernel='rbf',
    gamma=15.0,        # RBF kernel bandwidth
    fit_inverse_transform=True,  # Enable reconstruction
    alpha=0.01,         # Ridge regularization for inverse
    random_state=42
)
X_kpca = kpca.fit_transform(X)

# Visualize comparison
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

axes[0].scatter(X[:, 0], X[:, 1], c=y, cmap='viridis', s=10)
axes[0].set_title('Original Data')

axes[1].scatter(X_pca[:, 0], X_pca[:, 1], c=y, cmap='viridis', s=10)
axes[1].set_title('Linear PCA (cannot separate)')

axes[2].scatter(X_kpca[:, 0], X_kpca[:, 1], c=y, cmap='viridis', s=10)
axes[2].set_title('Kernel PCA (RBF) — separated!')

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

# Available kernels
for kernel in ['linear', 'rbf', 'poly', 'sigmoid', 'cosine']:
    kpca_k = KernelPCA(n_components=2, kernel=kernel, random_state=42)
    try:
        X_k = kpca_k.fit_transform(X)
        print(f'{kernel:>10}: shape={X_k.shape}')
    except Exception as e:
        print(f'{kernel:>10}: failed — {e}')

import numpy as np
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
import joblib
import time
import json

class PCAFeatureCompressor:
    """Production-ready PCA feature compression for inference pipelines.

    Designed for real-time serving with sub-ms latency.
    Typical use: compress 768-dim BERT embeddings to 50-dim for
    recommendation serving at Flipkart, Myntra, or Swiggy.
    """

    def __init__(self, n_components=50, variance_threshold=None):
        self.n_components = n_components
        self.variance_threshold = variance_threshold
        self.pipeline = None
        self.metadata = {}

    def fit(self, X_train):
        """Fit the compression pipeline on training data."""
        if self.variance_threshold:
            pca = PCA(n_components=self.variance_threshold, svd_solver='full')
        else:
            pca = PCA(
                n_components=self.n_components,
                svd_solver='randomized',  # Fast for large p
                random_state=42
            )

        self.pipeline = Pipeline([
            ('scaler', StandardScaler()),
            ('pca', pca),
        ])
        self.pipeline.fit(X_train)

        # Store metadata for monitoring
        fitted_pca = self.pipeline.named_steps['pca']
        self.metadata = {
            'input_dim': X_train.shape[1],
            'output_dim': fitted_pca.n_components_,
            'compression_ratio': X_train.shape[1] / fitted_pca.n_components_,
            'explained_variance': float(sum(fitted_pca.explained_variance_ratio_)),
            'n_training_samples': X_train.shape[0],
        }
        return self

    def transform(self, X):
        """Transform features at inference time."""
        return self.pipeline.transform(X)

    def benchmark(self, X_sample, n_runs=1000):
        """Measure inference latency."""
        # Warmup
        for _ in range(10):
            self.transform(X_sample[:1])

        # Benchmark single-sample latency
        latencies = []
        for _ in range(n_runs):
            start = time.perf_counter_ns()
            self.transform(X_sample[:1])
            latencies.append((time.perf_counter_ns() - start) / 1e6)

        return {
            'p50_ms': np.percentile(latencies, 50),
            'p95_ms': np.percentile(latencies, 95),
            'p99_ms': np.percentile(latencies, 99),
        }

    def save(self, path):
        """Serialize for deployment."""
        joblib.dump({
            'pipeline': self.pipeline,
            'metadata': self.metadata,
        }, path)

    @classmethod
    def load(cls, path):
        """Load from serialized artifact."""
        artifacts = joblib.load(path)
        instance = cls()
        instance.pipeline = artifacts['pipeline']
        instance.metadata = artifacts['metadata']
        return instance


# Example usage
X_train = np.random.randn(50000, 768)  # BERT embeddings
X_test = np.random.randn(1000, 768)

compressor = PCAFeatureCompressor(n_components=50)
compressor.fit(X_train)

print('Metadata:', json.dumps(compressor.metadata, indent=2))
print('Latency:', compressor.benchmark(X_test))

# Save and reload
compressor.save('pca_compressor.joblib')
loaded = PCAFeatureCompressor.load('pca_compressor.joblib')
X_compressed = loaded.transform(X_test)
print(f'Output: {X_compressed.shape}')  # (1000, 50)

# PCA configuration for a production recommendation pipeline
pca_config:
  n_components: 50
  svd_solver: randomized
  random_state: 42
  whiten: false          # Set true if downstream model is k-means
  # Preprocessing
  scaler: standard       # Options: standard, minmax, robust, none
  # Component selection
  variance_threshold: null  # Set to 0.95 to auto-select k
  # Incremental PCA (for large datasets)
  incremental: false
  batch_size: 5000
  # Serialization
  artifact_path: /models/pca/v2/
  artifact_format: joblib  # Options: joblib, onnx, pickle
  # Monitoring
  log_explained_variance: true
  alert_if_variance_below: 0.90
  reconstruction_error_threshold: 0.05