Principal Component Analysis (PCA)

September 2026 · 17 min read

PCA finds directions along which the training data varies most. It gives each row new coordinates along those directions. Keeping only a few coordinates compresses the row, while discarding variation in the remaining directions.

The central picture is a cloud of points and a line through its center. Project the points onto the line. PCA chooses the line that retains the most variance; equivalently, it minimizes squared reconstruction distance among such one-dimensional linear projections.

Start with four points and remove their mean

Our training rows are A = (3, 2), B = (4, 1), C = (6, 5), and D = (7, 4). Treat the two features as comparable measurements in fixed units. Their mean is $\mu=(5,3)$. Subtracting it places the cloud around the origin:

One coordinate along a line

Let $u$ be a unit-length direction vector. To get a row's scalar coordinate $z_i$, multiply matching entries of $u$ and the centered row, then add them: this is the dot product. To reconstruct an approximation in the original feature space, move $z_i$ units along $u$ and add the mean back:

The projected coordinate is a score. It is not the same object as the direction $u$. One score replaces two input values here, and reconstruction shows what information that compression loses.

At 45°, $u=(1,1)/\sqrt2$. Row A's score is $(-2-1)/\sqrt2\approx-2.121$. Its reconstructed position is (3.5, 1.5), half a unit away on each feature, so its squared reconstruction distance is $0.5^2+0.5^2=0.5$. Row B lands at the same reconstructed point.

Keeping the horizontal axis alone retains 50% of the variance and has average reconstruction error 2.5. The 45° PCA axis retains 90% and lowers that error to 0.5. Keeping both orthogonal axes reconstructs every training point exactly, but performs no dimensionality reduction.

How does PCA find the axis?

The sample covariance matrix summarizes how the centered features vary. Its diagonal entries are the individual feature variances; an off-diagonal entry measures how two features move together. With $N$ training rows, form the matrix by averaging their centered outer products using denominator $N-1$:

For any unit direction $u$, the variance of its projected scores is $u^\top S u$. PCA maximizes this quantity. The maximizing direction is an eigenvector of $S$: applying $S$ to it only scales it. That scale, the eigenvalue, is the variance along the direction.

Here the first direction $(1,1)/\sqrt2$ has variance 6. The perpendicular direction has variance $2/3$. The first therefore retains $6/(6+2/3)=0.9$, or 90%, of the total variance. Later principal directions maximize the remaining variance while staying orthogonal to earlier ones.

Explained variance is not prediction accuracy

Consider a separate four-row dataset with inputs $(x_1,x_2)$ taking all combinations of $x_1\in\{-10,10\}$ and $x_2\in\{-1,1\}$. Let the label be determined entirely by whether $x_2$ is positive. The features have zero covariance, and the variance of $x_1$ is 100 times larger.

Explained variance measures spread in the representation you supplied. It does not certify useful prediction, causal relevance or removal of noise. If reduction feeds a supervised model, choose and evaluate the number of components using that model's actual task metric.

Scaling changes which variation matters

Scikit-learn's PCA centers features automatically but does not standardize their scales. Converting one measurement to a different unit can therefore change the leading direction. Standardization asks PCA to work with equalized marginal variances; keeping original scales asks a different question.

Choose according to what those units mean. Equal treatment may help when scales are arbitrary, while standardization can remove meaningful amplitude differences. Fit any scaler on training data and save it with PCA. Whitening is a further operation that rescales retained scores; it changes their relative distance scales and does not restore discarded directions.

Fit once, then transform new rows

Save the training mean and retained directions. A new point (8, 4) is centered using the same mean (5, 3), gets score $4/\sqrt2\approx2.828$ on the chosen 45° axis, and reconstructs to (7, 5). Recomputing PCA for that point or batch would change the coordinate system expected by a downstream model.

import numpy as np
from sklearn.decomposition import PCA

X = np.array([[3,2], [4,1], [6,5], [7,4]], dtype=float)
full = PCA(n_components=2, svd_solver="full").fit(X)
one = PCA(n_components=1, svd_solver="full").fit(X)
scores = one.transform(X)
reconstructed = one.inverse_transform(scores)
error = np.mean(np.sum((X - reconstructed) ** 2, axis=1))
print(full.mean_)                       # [5, 3]
print(full.explained_variance_)         # [6, 2/3]
print(full.explained_variance_ratio_)   # [0.9, 0.1]
print(one.components_)                  # [1/sqrt(2), 1/sqrt(2)], up to sign.
print(scores.ravel())                   # [-2.1213, -2.1213, 2.1213, 2.1213], up to sign.
print(reconstructed)                    # [3.5,1.5] twice, [6.5,4.5] twice.
print(error)                            # 0.5: squared Euclidean error per row.

new = np.array([[8, 4]], dtype=float)
print(one.transform(new))              # 2.8284, with the orientation above.
print(one.inverse_transform(one.transform(new)))  # [[7, 5]]

# A separate counterexample: the low-variance coordinate determines the label.
X_signal = np.array([[-10,-1], [-10,1], [10,-1], [10,1]], dtype=float)
labels = np.array([0,1,0,1])
signal_pca = PCA(n_components=1, svd_solver="full").fit(X_signal)
print(signal_pca.explained_variance_ratio_)  # [100/101], about 99.01%.
print(signal_pca.transform(X_signal).ravel())
# Each retained coordinate value appears with both labels.

Put scaling, PCA and the downstream estimator inside the pipeline passed to validation. Even though PCA does not use labels, fitting it on evaluation rows lets those rows influence the learned representation. Use validation for component-count choices and a separate evaluation set or nested cross-validation for the final comparison.

PCA fits a linear subspace. A curved low-dimensional structure may need a nonlinear representation; kernel PCA is the next step. Outliers can also dominate covariance, so inspect the measurements before interpreting a principal direction.

Optional implementation: fit a two-feature PCA model

The C++ version solves the two-feature case directly from its sample covariance. It saves the mean and leading axis so the same projection can be applied to new rows.

#include <array>
#include <cmath>
#include <stdexcept>
#include <vector>
using Point = std::array<double, 2>;
struct PCA2 { Point mean, axis; double first_variance, total_variance; };

PCA2 fit_pca2(const std::vector<Point>& X) {
    if (X.size() < 2) throw std::invalid_argument("at least two training rows required");
    Point mean{0,0};
    for (const auto& x : X) {
        for (int j = 0; j < 2; ++j) {
            if (!std::isfinite(x[j])) throw std::invalid_argument("finite inputs required");
            mean[j] += x[j] / X.size();
        }
    }
    double a = 0, b = 0, c = 0;
    for (const auto& x : X) {
        double u = x[0] - mean[0], v = x[1] - mean[1];
        a += u*u / (X.size()-1);
        b += u*v / (X.size()-1);
        c += v*v / (X.size()-1);
    }
    double total = a + c;
    if (!std::isfinite(total) || !std::isfinite(b))
        throw std::overflow_error("covariance overflow");
    if (!(total > 0)) throw std::domain_error("no variance to explain");
    // Leading direction of the symmetric covariance [[a,b],[b,c]].
    double angle = .5 * std::atan2(2*b, a-c);
    Point axis{std::cos(angle), std::sin(angle)};
    double first = .5 * (total + std::hypot(a-c, 2*b));
    return {mean, axis, first, total};
}

double project(const PCA2& model, const Point& x) {
    if (!std::isfinite(x[0]) || !std::isfinite(x[1]))
        throw std::invalid_argument("finite prediction input required");
    return (x[0]-model.mean[0])*model.axis[0]
         + (x[1]-model.mean[1])*model.axis[1];
}
Point reconstruct(const PCA2& model, double score) {
    return {model.mean[0]+score*model.axis[0], model.mean[1]+score*model.axis[1]};
}

References and further reading

• Shlens: A Tutorial on Principal Component Analysis — geometric intuition, covariance, change of basis and the connection to SVD.
• Scikit-learn: PCA — centering, fitted outputs, explained variance and solver behavior.
• Scikit-learn: Data leakage — keeping learned transformations inside training boundaries.