Dimensionality Reduction

This code demonstrates the complete PCA workflow on the Iris dataset, reducing it from 4 features down to just 2. The first critical step is data standardization using StandardScaler().fit_transform(X) — this scales all features to have a mean of 0 and standard deviation of 1. This step is absolutely essential for PCA because without it, features with naturally larger values (like measurements in centimeters vs millimeters) would dominate the principal components purely due to their scale, not their actual importance. Next, we create a PCA object with PCA(n_components=2) to specify that we want exactly 2 output dimensions. The fit_transform() method does two things in one step: it first "fits" PCA by learning the principal components (the directions of maximum variance in the data), then "transforms" the data by projecting it onto these new axes. The explained_variance_ratio_ attribute is incredibly useful — it tells us what percentage of the original information (variance) each component captures. For example, if you see [0.73, 0.22], it means the first principal component captures 73% of the variance and the second captures 22%, for a total of 95% of the original information preserved in just 2 dimensions instead of 4!

This code shows two approaches for determining the optimal number of principal components to keep. First, we fit PCA with no n_components argument, which means it calculates ALL possible components (64 for the digits dataset). We then use np.cumsum() to calculate the cumulative sum of explained variance ratios — this transforms something like [0.12, 0.22, 0.10...] into a running total [0.12, 0.34, 0.44...], showing how much total variance is explained as we add more components. Using np.argmax(cumsum >= 0.95), we find the first index where the cumulative variance reaches or exceeds 95%, telling us exactly how many components we need to preserve that much information. But there's an even easier way! You can pass a decimal value directly to PCA using PCA(n_components=0.95), and scikit-learn will automatically determine how many components are needed to retain 95% of the variance. For the 64-feature digits dataset, you typically only need about 29 components to keep 95% of the variance — that's a reduction of more than half the features while losing only 5% of the information! This technique helps you balance dimensionality reduction against information loss.

We start by importing all necessary libraries: PCA for dimensionality reduction, StandardScaler for normalizing our data, load_digits to get the handwritten digits dataset, numpy for numerical operations, and matplotlib.pyplot for creating visualizations. After loading the digits dataset, we immediately scale it using StandardScaler because PCA is sensitive to feature scales — without scaling, features with larger numerical values would unfairly dominate the principal components. The digits dataset has 1797 samples (handwritten digit images) with 64 features each (8×8 pixels).

Here we create a PCA object without specifying n_components, which tells scikit-learn to keep ALL principal components (64 in this case). This is essential for scree plots because we need to see the variance explained by every component to make an informed decision about how many to keep. The explained_variance_ratio_ attribute is an array where each element tells us what fraction of the total variance that component captures. For example, if the first element is 0.12, it means PC1 explains 12% of the total variance in the data. The components are automatically sorted from highest to lowest variance.

This creates the classic scree plot as a bar chart. The name "scree plot" comes from geology — scree is the pile of broken rocks at the base of a cliff, and the shape of this plot resembles that slope. We use plt.subplots(1, 2) to create a figure with two side-by-side plots. The first plot shows the first 20 components as bars, where the height of each bar represents how much variance that single component explains. You'll typically see the first few bars being tall (PC1 might explain 12%, PC2 might explain 10%), then a steep drop-off, followed by a long "tail" of short bars. The "elbow" is where the bars stop dropping quickly and become uniformly short — that's your sweet spot for choosing components.

The second plot shows cumulative explained variance, which is often more practical for decision-making. We use np.cumsum() to calculate a running total — if PC1 explains 12% and PC2 explains 10%, the cumsum would be [0.12, 0.22, ...]. This tells us "with N components, how much total variance do I capture?" We add horizontal reference lines at 90% (orange) and 95% (red) thresholds using axhline() — these are common targets in practice. You can simply look at where your curve crosses these lines to determine how many components you need. For instance, if the curve crosses 0.95 at x=29, you need 29 components to retain 95% of the variance. The plt.tight_layout() ensures the plots don't overlap.

Finally, we print a numerical summary showing both individual and cumulative variance for the first 10 components. This gives you exact numbers instead of reading from the plot. For the digits dataset, you might see output like: PC1 explains 0.120 (12.0% cumulative), PC2 explains 0.097 (21.7% cumulative), and so on. By PC10, you might already have ~60% cumulative variance. This table helps you make precise decisions — if you need exactly 95% variance, you can see exactly which component crosses that threshold. The combination of visual plots and numerical output gives you a complete picture for choosing the optimal number of components.

This code reveals what's happening "under the hood" of PCA by examining the pca.components_ attribute, which contains the "recipe" for each principal component. Think of each component as a weighted combination of your original features. The components_ attribute is a matrix where each row represents one principal component, and each column shows how much the corresponding original feature contributes to that component (these weights are called "loadings"). For example, if the loading for "petal length" is 0.5, it means petal length contributes positively to that component — samples with longer petals will have higher values on this component. A negative loading like -0.3 means the feature contributes in the opposite direction. The magnitude (absolute value) of each loading tells you how important that feature is: loadings close to 0 mean the feature barely influences the component, while loadings close to 1 or -1 mean the feature is very influential. This interpretability is one of PCA's biggest advantages over methods like t-SNE and UMAP! You can actually explain your components in human terms, like "PC1 seems to represent overall flower size because all measurement features have high positive loadings."

We start by importing TSNE from scikit-learn's manifold module — t-SNE is classified as a "manifold learning" technique because it tries to unfold the complex curved surface (manifold) on which high-dimensional data lives. We load the digits dataset which contains 1797 images of handwritten digits, each represented as 64 pixel values (8×8 images). Just like with PCA, we apply StandardScaler to normalize the data so all features have similar scales. We also keep the target labels y (digits 0-9) so we can color-code our visualization later to see if t-SNE correctly groups similar digits together.

Here we configure t-SNE with its most important parameters. The n_components=2 means we want a 2D output for plotting — unlike PCA, t-SNE is almost never used for more than 3 dimensions because its purpose is visualization. The perplexity=30 is the key hyperparameter: it roughly means "each point considers about 30 neighbors when deciding where to position itself." Think of perplexity as the "attention span" of each point. Lower values (5-10) make points focus only on their immediate neighbors, creating tight local clusters. Higher values (30-50) make points aware of a broader neighborhood, showing more global structure. We set random_state=42 because t-SNE uses random initialization and stochastic optimization — without this, you'd get different results each run! The n_iter=1000 controls how many optimization steps t-SNE takes to arrange the points; more iterations usually mean better results but take longer.

Now we apply t-SNE using fit_transform(X). This is a crucial difference from PCA: t-SNE only has fit_transform(), NOT separate fit() and transform() methods. This means you cannot train t-SNE on one dataset and then embed new points later — if you get new data, you must re-run t-SNE on the entire dataset including the new points. The output X_tsne has shape (1797, 2), meaning each of our 1797 digit images now has just 2 coordinates that we can plot on a 2D scatter plot. The actual coordinate values are arbitrary (they might range from -50 to +50 or any other range) — what matters is the relative positions of points: similar digits should cluster together. When you plot this, you'll see 10 distinct clusters corresponding to digits 0-9!

We set up an experiment to see how different perplexity values affect t-SNE's output. We import NearestNeighbors which we'll use later to measure how spread out the resulting clusters are. Perplexity is a crucial hyperparameter that essentially controls the "scale" at which t-SNE looks at your data — should it focus on very local neighborhoods (low perplexity) or consider broader regions (high perplexity)? There's no single correct answer, which is why we'll test multiple values and compare the results.

We test four different perplexity values: 5, 15, 30, and 50. With perplexity=5, each point only "cares about" its 5 closest neighbors — this creates many small, tight, well-separated clusters, but you might miss the bigger picture of how clusters relate to each other. With perplexity=15-30, you get a balanced view that works well for most datasets. With perplexity=50, points consider a much wider neighborhood, which can cause smaller clusters to merge together but gives a better sense of overall structure. A good rule of thumb: perplexity should be less than the number of points in your smallest expected cluster, and definitely less than the total number of samples.

To objectively compare embeddings, we measure the average distance between neighboring points using NearestNeighbors. We use n_neighbors=6 because the first neighbor is always the point itself (with distance 0), so distances[:, 1:] gives us the 5 actual nearest neighbors. Lower average distance means points are packed more tightly together. You'll typically see that low perplexity values create tighter local clusters (lower average distance), while high perplexity values spread points out more evenly (higher average distance). This metric helps you quantify what you're seeing visually in the plots. The best practice is to try 2-3 different perplexity values, visualize each one, and choose the one that best reveals the structure you're interested in.

We're setting up a speed comparison experiment to demonstrate why PCA before t-SNE is a professional best practice. The digits dataset has 64 features (8×8 pixel images), which is high enough that t-SNE will take noticeable time. We import the time module to measure execution time for each approach. The question we're answering: Is it faster to run t-SNE directly on all 64 dimensions, or to first reduce to fewer dimensions with PCA?

This is the straightforward approach that most beginners use: apply t-SNE directly to the full 64-dimensional data. We record the time before and after using time.time() to measure how long it takes. The reason this is slow is that t-SNE must compute distances between every pair of points in 64-dimensional space — that's a lot of calculations! For N samples with D dimensions, the complexity scales with both N² (all pairs) and D (dimensions). This approach works fine for small datasets, but becomes painfully slow as data size or dimensionality increases.

Here we apply PCA first, reducing from 64 dimensions to 50. This is extremely fast — typically just milliseconds even for large datasets — because PCA is a linear operation that just multiplies matrices. The explained_variance_ratio_ tells us how much information we kept; with 50 components, we typically retain 95%+ of the original variance, meaning we've lost almost nothing important while making the data much easier to work with. This preprocessing step removes noise and redundancy, leaving only the most important structure.

Now we run t-SNE on the PCA-reduced data instead of the original. Even though we're still computing the same embedding, t-SNE now only needs to calculate distances in 50 dimensions instead of 64 — and more importantly, those 50 dimensions contain cleaner, more meaningful information without noise. The result is usually just as good (often better!) than direct t-SNE, but computed faster.

We compare the total time of both approaches. The PCA + t-SNE combination is typically 2-10x faster, and the speedup becomes even more dramatic with higher-dimensional data (like images with thousands of pixels or text embeddings with hundreds of dimensions). The original t-SNE authors themselves recommend this preprocessing step for any dataset with more than about 50 features. Beyond speed, PCA removes noise that might confuse t-SNE, often producing cleaner visualizations. This two-step approach is considered the professional standard and should be your default workflow for any serious dimensionality reduction task!

This code applies UMAP (Uniform Manifold Approximation and Projection) to reduce the 64-dimensional digits dataset to 2D for visualization. Note that UMAP is a separate package — you need to install it with pip install umap-learn (not just "umap"). Like t-SNE, we typically use n_components=2 for visualization purposes. The n_neighbors=15 parameter is similar to t-SNE's perplexity — it controls how many neighbors UMAP considers for each point when building its high-dimensional graph. Lower values focus on very local structure (tight clusters), while higher values capture more global relationships (smoother, more connected embeddings). The min_dist=0.1 parameter controls how tightly UMAP packs points together in the output. With min_dist=0.0, points can be placed right on top of each other, creating very dense, compact clusters. With higher values like 0.5 or 1.0, points are spread out more evenly. Like t-SNE, UMAP is stochastic (involves randomness), so we set random_state=42 for reproducible results. The major advantages of UMAP over t-SNE are speed (10-100x faster on large datasets), better preservation of global structure (distances between clusters are more meaningful), and the ability to transform new data (which we'll see next). For most modern applications, UMAP has become the go-to choice for non-linear dimensionality reduction.

This code demonstrates UMAP's killer feature that sets it apart from t-SNE: the ability to transform new, unseen data points without refitting the entire model. We first split our data using train_test_split into 80% training and 20% test sets to simulate a real-world scenario where you have training data and later receive new data. We then call reducer.fit(X_train) which learns the embedding structure from only the training data. After fitting, we can use reducer.transform(X_train) to get the 2D coordinates for our training points, but here's the magic: we can also call reducer.transform(X_test) to embed completely new test data that UMAP never saw during training! The new points are projected into the same 2D space in a way that's consistent with the original embedding. This is impossible with t-SNE — if you want to add even a single new point with t-SNE, you must re-run the entire algorithm on all your data from scratch. UMAP's transform capability is crucial for production machine learning pipelines: you can train UMAP once on your training data, save the fitted model, and then quickly embed new data points as they arrive without expensive recomputation. This makes UMAP suitable for real-time applications, streaming data, and any scenario where you need to process new samples on-the-fly.

This code demonstrates the simplest form of feature selection: removing features that barely vary across samples. The logic is intuitive — if a feature has nearly the same value for every sample, it cannot possibly help distinguish between different samples or classes! We create a VarianceThreshold(threshold=0.5) selector which will remove any feature whose variance is below 0.5. When we call fit_transform(X), scikit-learn calculates the variance of each feature across all samples and keeps only those features with variance above our threshold. The get_support() method returns a boolean array showing which original features survived (True) or were removed (False). For example, imagine a "country" feature where 99% of your samples are from the same country — this feature has almost zero variance and provides virtually no discriminative information, so it should be removed. Variance threshold is an excellent first step in any feature selection pipeline because it's extremely fast, requires no labels (unsupervised), and removes obviously useless features before you apply more sophisticated methods. A practical tip: start with threshold=0 to remove only constant features, then gradually increase the threshold to see how aggressively you can trim features without hurting model performance.

This code identifies and removes redundant features by detecting highly correlated pairs. When two features are strongly correlated, they essentially carry the same information — keeping both just adds noise and computational overhead. We compute the correlation matrix using np.corrcoef(X.T), which calculates the Pearson correlation coefficient between every pair of features. The result is a matrix where the value at position [i,j] represents the correlation between feature i and feature j. Correlation ranges from -1 to +1: a value of +1.0 means the features move perfectly together (when one goes up, the other goes up by a proportional amount), -1.0 means they move perfectly opposite (when one goes up, the other goes down), and 0 means no linear relationship. We set threshold = 0.95, meaning any pair with absolute correlation above 0.95 is considered redundant. For each such pair, we add one feature to our removal set using to_remove.add(j) — we keep the first feature and remove the second. A classic example: if your dataset has both "height_inches" and "height_cm", they're 100% correlated since one is just a unit conversion of the other. There's no reason to keep both! We use 0.95 as a conservative threshold, but you could use 0.8 or 0.9 to be more aggressive in removing features, especially if you suspect multicollinearity issues.

This code uses statistical tests to score each feature based on how well it helps distinguish between classes, then keeps only the top K highest-scoring features. We use f_classif as our scoring function, which performs an ANOVA F-test for each feature. This test essentially asks: "Do the values of this feature differ significantly between classes?" If a feature has very different average values for malignant vs benign tumors, it gets a high F-score because it's useful for classification. We create SelectKBest(score_func=f_classif, k=10) to keep only the 10 features with the highest scores. Notice that unlike PCA, we call fit_transform(X, y) with both features AND labels — this is because SelectKBest is a supervised method that needs to know the target variable to evaluate how useful each feature is for prediction. After fitting, selector.scores_ contains the F-score for each original feature, letting you see exactly how useful each feature is (higher score = more discriminative). Different scoring functions are available for different situations: mutual_info_classif captures non-linear relationships between features and target, chi2 works with non-negative features (like word counts in text), and f_regression is used when your target is continuous instead of categorical. SelectKBest is perfect when you want interpretable results and need to explain which features your model relies on.

This code demonstrates Recursive Feature Elimination (RFE), which uses an actual machine learning model to iteratively identify and remove the least important features. We use a RandomForestClassifier as our estimator because Random Forest naturally provides feature importance scores based on how much each feature helps split the data and improve predictions. We create RFE(estimator=model, n_features_to_select=10, step=1), which will recursively remove one feature at a time (step=1) until only 10 features remain. Here's how RFE works step by step: first, it trains the model on all features and looks at the feature importances. It removes the single least important feature, then trains the model again on the remaining features. This process repeats — train, remove the worst, train, remove the worst — until only 10 features are left. The rfe.support_ attribute is a boolean array showing which features survived the elimination process (True means kept, False means removed). The rfe.ranking_ attribute is even more informative: features with ranking 1 were kept, ranking 2 means "removed in the last round," ranking 3 means "removed second-to-last," and so on. RFE's major advantage over simple statistical methods like SelectKBest is that it considers feature interactions — since a real model evaluates importance, the removal of one feature might change the importance of others. The downside is speed: RFE trains the model many times (once per elimination round), so it can be slow for large datasets or complex models.

We set up a fair comparison by importing all three dimensionality reduction methods: PCA for linear reduction, TSNE for non-linear visualization, and umap (which needs to be installed separately with pip install umap-learn). We also import silhouette_score from sklearn's metrics module — this will measure how well-separated the resulting 2D clusters are. The score ranges from -1 to +1, where higher values mean clusters are more distinct and well-separated. We use the digits dataset (1797 samples, 64 features) because it has known class labels (digits 0-9), allowing us to objectively measure which method creates the best cluster separation.

PCA is the fastest method by far — often completing in just a few milliseconds even on larger datasets. This is because PCA only involves matrix operations (computing eigenvectors of the covariance matrix), which are highly optimized. However, you'll notice its silhouette score is typically lower than t-SNE or UMAP. This is because PCA is a linear method — it can only find straight directions in the data and cannot capture curved or complex cluster shapes. The upside is that PCA preserves global structure (relative distances are maintained) and can transform new data points using the same learned projection. The explained_variance_ratio_ tells us how much information we retained; with just 2 components from 64, we typically keep only 20-30% of the variance, which explains the lower cluster quality.

t-SNE typically produces the highest silhouette score — meaning the tightest, best-separated clusters — because it's specifically designed to preserve local neighborhood structure. Points that were neighbors in 64D will be neighbors in 2D. However, this quality comes at a cost: t-SNE is significantly slower (seconds instead of milliseconds) because it uses iterative optimization to arrange points. The major limitation is that t-SNE cannot transform new data — if you get a new sample, you must re-run t-SNE on the entire dataset including the new point. This makes t-SNE unsuitable for production ML pipelines but perfect for exploratory data analysis and creating beautiful visualizations for reports or publications. The random_state=42 ensures reproducibility since t-SNE uses random initialization.

UMAP often achieves the best of both worlds: it produces cluster quality comparable to t-SNE (sometimes even better) while being much faster — often only slightly slower than PCA. UMAP's speed advantage becomes even more dramatic on larger datasets where it can be 10-100x faster than t-SNE. Critically, UMAP CAN transform new data using reducer.transform(new_data), making it suitable for production ML systems where you need to embed new samples without re-running the entire algorithm. UMAP also tends to preserve global structure better than t-SNE, meaning the relative positions of clusters in the visualization more accurately reflect their relationships in the original high-dimensional space. For most modern applications, UMAP has become the default choice for non-linear dimensionality reduction.

This summary table gives you a clear decision framework. For ML preprocessing (feeding reduced data into a classifier), use PCA because it's deterministic, fast, interpretable, and can transform new data — plus its linear nature works well as a preprocessing step. For visualization in most cases, use UMAP because it balances speed and quality, and can handle large datasets that would take forever with t-SNE. For publication-quality figures or when you need the absolute best cluster visualization on smaller datasets, use t-SNE because it often produces the most visually striking separation. The "Transform New Data?" column is crucial for production systems — if you're building an application that needs to embed user uploads or streaming data, t-SNE is simply not an option, leaving you with PCA or UMAP depending on whether you need linear or non-linear reduction.