1. Dimensionality, variables, and relationships

A table with many columns provides numerous descriptions of its observations, but this does not guarantee that every column contributes a distinct relationship. Several measurements may vary together; others may introduce specific information or mainly reflect noise. Dimensionality reduction constructs a representation with fewer coordinates while seeking to preserve properties relevant to a specified question.

A data matrix $X\in {ℝ}^{n\times p}$ contains n observations in its rows and p quantitative variables in its columns. The original representation space has dimension p; the rank of the centered data cannot exceed min(n−1, p). Intrinsic dimension expresses, more broadly, how many degrees of freedom are needed to describe the structure of the data at the required level of precision.

Two related operations are involved: selecting original variables and extracting combinations of them. PCA directly performs the second and provides a mediated route to assessing the relevance of the first. This operational distinction does not remove their conceptual connection: the derived coordinates are constituted by the original variables and their relationships of covariation.

The curse of dimensionality encompasses difficulties that arise under certain conditions as the number of dimensions increases: greater data requirements, concentration of distances, or instability of estimates. It is not a law rendering every wide table useless. Sample size, effective structure, distribution, metric, and task all matter.

The problem of fitting lines and planes to systems of points appears in Pearson (1901). Contemporary machine learning motivations extend its applications, but should be distinguished from the history of the procedure. PCA and the later theory of random projections address related problems through different guarantees.

2. What a transformation preserves

Centering subtracts the mean of each variable. Translation changes the origin and preserves distances between observations. Standardization additionally divides each column by its standard deviation, changing the weight of its differences. This second operation generally changes the geometry expressed in the original units.

Neither centering nor standardization, by itself, transforms an arbitrary distribution into a Gaussian distribution. Standardizing a column is also different from dividing a direction vector by its length: the latter produces a unit vector.

Let V be a full orthogonal matrix. For any difference d between two observations, $\Vert dV{\Vert }_2=\Vert d{\Vert }_2$ holds. Splitting the basis into retained and discarded directions, the Pythagorean identity gives:

A reduced PCA representation is therefore not necessarily an isometric embedding. It may even map two distinct observations to the same point. Exact preservation in a lower-dimensional subspace is possible when all differences between the data points belong to the retained subspace; otherwise, distortion must be evaluated.

3. The mathematical construction of PCA

We use a consistent convention: X_c is the centered matrix, S is its covariance matrix with divisor n−1, and v_j is a unit direction vector. When the variables are standardized, their covariance matrix is the correlation matrix of the original data.

The first direction maximizes projected variance among all unit directions. Subsequent directions maximize the remaining variance subject to orthogonality to the previous ones. Without the norm constraint, multiplying the coefficients by an arbitrarily large constant would make the maximum undefined.

Singular value decomposition provides a direct and numerically convenient calculation. If $X_c=U\Sigma V^{\mathsf{T}}$, the columns of V are the principal directions, the scores are T=X_cV=UΣ, and the eigenvalues of S are ${\lambda }_j={\sigma }_j^2/(n-1)$. A singular value of the data matrix is not simply the square root of an eigenvalue of its covariance matrix: the divisor matters. The function prcomp() uses SVD and this sample convention (R Core Team, n.d.-c).

Maximizing retained variance and minimizing this squared error are equivalent formulations of the same problem. The equivalence follows because total variation decomposes into orthogonal projected and residual parts. It does not depend on separating a supposedly “analytical” methodology from a “numerical” one.

Descriptive PCA does not require normally distributed data (Jolliffe & Cadima, 2016). Additional assumptions arise when making inferences, specifying probabilistic models, or setting detection thresholds. With repeated eigenvalues, the individual directions within their subspace are not unique; even without repetition, each direction may reverse sign without changing the solution.

4. Eigenvalues, averages, and variance

In linear algebra, an eigenvalue is the factor by which a transformation multiplies an eigenvector: $Av=\lambda v$. In PCA, the transformation is the covariance or correlation matrix, and that factor acquires a precise statistical interpretation: the variance of the projection onto the corresponding unit direction.

An average therefore underlies the eigenvalue: the expectation of the squared, centered projections of the observations. Dispersion is measured along a direction defined by the eigenvector. This relationship allows the eigenvalue to be understood as a quantitative summary of variation in that direction.

This average of squared projections is not the average of the eigenvector’s coordinates. The coordinates are variable weights; the scores are values for observations; variance summarizes the dispersion of those scores. Confusing these three objects changes the meaning of the result.

Furthermore, v and −v share an eigenvalue, while their coordinate averages change sign. Multiplying the entire matrix S by ten preserves its eigenvectors and multiplies its eigenvalues by ten. An eigenvalue depends jointly on the transformation and the direction; the direction alone does not determine its magnitude.

The eigenvector determines how to project; the eigenvalue expresses how much variation the data contain in that projection.

5. A mediated assessment of variable relevance

An account of PCA limited to saying that it “creates new coordinates” leaves its interpretive role incomplete. The new coordinates are constructed from the original variables; by ranking the components and studying their composition, we can identify, through these mediations, which variables participate in the retained structure and to what extent.

Three questions should be distinguished: which variables define a direction; what fraction of each variable’s variation is represented; and how much that direction contributes to total variation. Answering them jointly requires examining coefficients, variable–component correlations, contributions, and eigenvalues. A variable may have little relevance to PC1 and considerable relevance to PC2 or PC3.

For standardized variables, the coordinates in the correlation circle are $\sqrt{{\lambda }_j}{v}_{aj}$. Their squares describe how well the variable is represented on each axis. In this standardized PCA, the conventional contribution of a variable to an axis is $100v_{aj}^2$ and sums to 100% across variables.

These expressions formalize the mediation: both the ranking of components, given by their eigenvalues, and the participation of the variable in them are involved. Since $S_{aa}={\sum }_{j=1}^p{\lambda }_j{v}_{aj}^2$, we have $0\le Q_a(k)\le 1$. For standardized variables, S_aa=1 and R_a(k)=Q_a(k).

The term loading follows different conventions across programs: it may designate eigenvector coefficients or coordinates scaled by the square root of the eigenvalue. This article always specifies the object being used. An observation’s squared cosine, in turn, concerns the representation of that observation rather than a variable’s contribution.

6. A fully worked PCA example

Consider four observations on two variables, with mean (3, 3). The example is constructed so that every step can be checked without relying on software output. We work with covariances, without standardizing.

The characteristic polynomial is $(10/3-\lambda {)}^2-(8/3{)}^2=0$. It gives ${\lambda }_1=6$ and ${\lambda }_2=2/3$. We choose the directions $v_1=(1,1{)}^{\mathsf{T}}/\sqrt{2}$ and $v_2=(-1,1{)}^{\mathsf{T}}/\sqrt{2}$. They have unit norm, are orthogonal, and satisfy the eigenvector equation.

The PC1 scores sum to zero. Their squares sum to 18, so their sample variance is 18/3=6. The squared PC2 scores sum to 2, and their variance is 2/3. This shows directly where the average of squared scores defining each eigenvalue enters the calculation.

For both variables, the reconstructed variance is $6(1/\sqrt{2}{)}^2=3$. Since their original variance is 10/3, each is represented at 90%. Each contributes 50% to PC1. These percentages differ because they answer different questions: the composition of an axis and the representation of a variable.

X <- rbind(A=c(1,2), B=c(2,1), C=c(4,5), D=c(5,4))
colnames(X) <- c("X1", "X2")
Xc <- scale(X, center=TRUE, scale=FALSE)
S <- cov(X)
V <- cbind(PC1=c(1,1), PC2=c(-1,1))/sqrt(2)
T <- Xc %*% V
lambda <- diag(t(V) %*% S %*% V)
Xc_hat <- T[,1,drop=FALSE] %*% t(V[,1,drop=FALSE])
X_hat <- sweep(Xc_hat, 2, colMeans(X), "+")
lambda
T
X_hat
sum((Xc-Xc_hat)^2)
stopifnot(max(abs(lambda-c(6,2/3))) < 1e-10,
          abs(sum((Xc-Xc_hat)^2)-2) < 1e-10,
          max(abs(t(V)%*%V-diag(2))) < 1e-10)

7. How many components should be retained?

The proportion of variance explained by component j is obtained by dividing its eigenvalue by the sum of all eigenvalues. The cumulative proportion measures how much a representation with k dimensions retains:

A rank-deficient matrix has zero eigenvalues. For centered data, the number of positive-variance components cannot exceed min(n−1, p). When p is large relative to n, a good description of the sample can coexist with unstable estimates of population directions. Consistency requires conditions on signal, noise, and the dimensional regime; comparing two integers alone does not settle it.

Parallel analysis has an antecedent in Horn (1965). Sparse PCA imposes or encourages zero coefficients to obtain combinations with few active variables. “Sparse” here refers to the coefficient structure, not simply to having little data. It can support interpretation or regularization under specified assumptions; it does not by itself guarantee recovery of the true structure.

8. Application: the ten decathlon events

The decathlon2 dataset contains 27 athletes, ten events, and three additional columns: rank, overall points, and competition. Following the conventional example in the documentation, we use the first 23 athletes and first ten columns as active data; the remaining four rows can subsequently be projected as supplementary observations. The columns Rank, Points, and Competition do not enter the fit (Kassambara & Mundt, n.d.-a).

The events use seconds, meters, and other performance measurements. Standardization prevents numerical scale from determining each column’s initial weight. In races, a lower time generally indicates better performance; in jumps and throws, a greater distance generally indicates better performance. The signs must be interpreted with this difference in mind.

# First run only, if these packages are missing:
# install.packages(c("factoextra", "cluster"))
data("decathlon2", package="factoextra")
deca_active <- decathlon2[1:23, 1:10]
stopifnot(all(vapply(deca_active, is.numeric, logical(1))),
          !anyNA(deca_active), all(vapply(deca_active, sd, 0) > 0))
pca_deca <- prcomp(deca_active, center=TRUE, scale.=TRUE)

# Fix an orientation to simplify comparisons of signs.
orient_pca <- function(p, anchor=1, sign=-1) {
  for (j in seq_len(ncol(p$rotation))) {
    a <- if (j==1) anchor else which.max(abs(p$rotation[,j]))
    target_sign <- if (j==1) sign else 1
    if (p$rotation[a,j]*target_sign < 0) {
      p$rotation[,j] <- -p$rotation[,j]
      p$x[,j] <- -p$x[,j]
    }
  }
  p
}
pca_deca <- orient_pca(pca_deca)
lambda_deca <- pca_deca$sdev^2
deca_eigenvalues <- data.frame(
  lambda=lambda_deca,
  percentage=100*lambda_deca/sum(lambda_deca),
  cumulative=100*cumsum(lambda_deca)/sum(lambda_deca))
variable_pc_cor <- sweep(pca_deca$rotation, 2, pca_deca$sdev, "*")
contribution <- 100*pca_deca$rotation^2
representation_3 <- rowSums(variable_pc_cor[,1:3,drop=FALSE]^2)
round(deca_eigenvalues, 4)
round(cbind(variable_pc_cor[,1:3], Q3=representation_3), 4)

# Supplementary projection using the fitted centering and scaling:
supplementary_scores <- predict(pca_deca,
                                newdata=decathlon2[24:27,1:10])
supplementary_scores

The eigenvalues 4.1242, 1.8385, and 1.2391 correspond to PC1, PC2, and PC3. Each component combines the ten events; its name is not replaced by that of one of the first three columns. In the chosen orientation, PC1 contrasts short-race times with several jumps and throws. PC2 combines, among other characteristics, pole vault and 1500 meters; because a longer time and a greater height have different sporting meanings, it should not automatically be converted into a single scale of “better performance.”

The 100-meter event contributes 17.54% to PC1 and is represented at 84.67% when three components are retained. The 1500-meter event contributes almost nothing to PC1, but its correlation with PC2 is approximately 0.784 and its contribution to that axis reaches 33.47%. Pole vault contributes approximately 35.46% to PC2. These results show why a mediated assessment of variable relevance must consider more than one axis.

set.seed(20261113)
null_eigenvalues <- replicate(500, {
  permuted <- apply(deca_active, 2, sample)
  eigen(cor(permuted), symmetric=TRUE, only.values=TRUE)$values
})
threshold95 <- apply(null_eigenvalues, 1, quantile, probs=.95)
data.frame(observed=lambda_deca, threshold95=threshold95,
           exceeds=lambda_deca > threshold95)
sum(lambda_deca > 1)         # Kaiser: 3
sum(lambda_deca > threshold95)  # Permutation parallel analysis: 1

9. Reading the plots and Minitab output

Software output represents different mathematical objects. Before interpreting it, identify what each point, arrow, percentage, and unit represents. This applies equally to R, Minitab, and other programs.

In Minitab (Minitab, n.d.), select the active quantitative columns and choose between covariance and correlation analysis. To replicate the geometry of the decathlon example, use the correlation matrix and the same 23 rows. Request eigenvalues, coefficients, scores, and plots, and compare percentages and subspaces. Signs can reverse across programs without changing the solution; divisors and score scales should also be checked.

In the financial example in the Minitab documentation, the first three eigenvalues are 3.5476, 2.1320, and 1.0447, accounting for approximately 84.1% cumulatively. Associations of PC1 with age, residence, employment, and savings support a proposed interpretation in terms of long-term financial stability; this label is a substantive reading of the coefficients and requires examining the variable definitions. Eigenvalues alone do not assign that meaning.

Some Spanish Minitab output headings use “influencias” for loadings. Here this concerns component composition rather than a regression influence measure. Nor is a loading generally a causal derivative of a variable with respect to a response phenomenon. Ordinary PCA does not define a dependent variable.

Mahalanobis distance weights coordinates by their dispersion and accounts for covariation. Retaining only some components changes this distance and can hide anomalies in discarded directions. If the covariance matrix is singular, the ordinary inverse does not exist; using a pseudoinverse or regularization changes the criterion and should be reported.

PCA and least-squares regression also differ in their objectives. A simple regression of Y on X minimizes residuals in the response direction; the PCA line fit minimizes distances orthogonal to the line in the chosen geometry. Principal component regression adds a supervised stage after PCA and must be validated as a predictive model.

10. Clusters and the k-means objective

Clustering organizes observations using a similarity criterion and a research purpose. A statistical group depends on the variables, their scale, the dissimilarity, and the algorithm. The resulting organization should be justified through these decisions and knowledge of the object being studied.

K-means seeks a partition of the observations into k groups and a centroid for each group. Its ordinary objective minimizes the sum of squared Euclidean distances to the centroids. The mean minimizes this criterion; absolute distances or other dissimilarities change the problem and do not automatically yield the same algorithm.

Common algorithms alternate assignments and updates or use other local improvements. Multiple starts reduce the risk of accepting a poor local solution, but do not certify a global minimum. In R, nstart controls how many initializations are tried when a number of centers is supplied (R Core Team, n.d.-b).

K-means is unsupervised learning: labels are constructed from the data rather than supplied as known training classes. Groups are disjoint and exhaustive when every observation receives a label. Hierarchical methods organize partitions at different levels; their groups are also disjoint when the tree is cut at a particular level.

The connection with PCA arises because both can be formulated through quadratic approximations: PCA reconstructs through a linear subspace, and k-means through centroids and discrete assignments. Ding and He (2004) study this relationship through a continuous relaxation of membership indicators. The connection does not guarantee that PC1 separates the groups or that the methods produce the same solution.

11. Application: USArrests and validation

USArrests contains 50 US states and four variables. Murder, Assault, and Rape record arrests per 100,000 residents in 1973; UrbanPop gives the urban population percentage. These are aggregate records of arrests and urbanization, not a direct, complete measurement of every crime committed or of individual people’s characteristics (R Core Team, n.d.-d).

We use the dataset included in R as distributed so that results are comparable with the historical example. Its documentation identifies a transcription error in Maryland’s urbanization value and provides additional rounding corrections. A substantive study should examine the corrected version and document sensitivity; the numbers below explicitly refer to the uncorrected built-in dataset.

library(cluster)
data("USArrests", package="datasets")
usa_z <- scale(USArrests)
dist_usa <- dist(usa_z, method="euclidean")
pca_usa <- orient_pca(prcomp(usa_z, center=FALSE, scale.=FALSE))

fits <- lapply(2:6, function(k) {
  set.seed(20261113+k)
  kmeans(usa_z, centers=k, nstart=100, iter.max=100)
})
k_comparison <- do.call(rbind, lapply(seq_along(fits), function(j) {
  a <- fits[[j]]
  data.frame(k=j+1, W=a$tot.withinss, B=a$betweenss,
    silhouette=mean(silhouette(a$cluster, dist_usa)[,3]))
}))
km2 <- fits[[1]]
round(k_comparison, 4)
km2$size
round(aggregate(USArrests, list(cluster=km2$cluster), mean), 3)
stopifnot(abs(km2$totss-km2$tot.withinss-km2$betweenss) < 1e-8)
summary(pca_usa)

The first two components explain approximately 62.01% and 24.74%, respectively, totaling 86.75%. PC1 and PC2 are combinations of all four variables. The 62.01% is neither a state’s murder percentage nor an average proportion of crime types: it is the fraction of total standardized variance captured by the first axis.

For k=2, the groups have sizes 20 and 30, a within-cluster sum of 102.8624, and a between-cluster sum of 93.1376. The total sum is 196: each of the four standardized columns contributes n−1=49. The between-cluster proportion is 93.1376/196≈47.52%; this percentage concerns the clustering criterion and should not be confused with the variance explained by PC1.

In this solution, the 20-state group has means of 12.165 for Murder, 255.250 for Assault, 68.400 for UrbanPop, and 29.165 for Rape. The 30-state group has means of 4.870, 114.433, 63.633, and 15.943. This description concerns the observed variables and historical period; it does not authorize causal inferences about urbanization or direct generalization to the present.

11.1. Silhouette: cohesion and separation

For observation i, a(i) is its average dissimilarity to the other observations in its cluster. For each other cluster, calculate the average dissimilarity between i and its members; b(i) is the minimum of these averages. The silhouette compares the two quantities (Rousseeuw, 1987):

A value near one indicates that the observation is comparatively closer to its own group; near zero, assignment is ambiguous under the chosen dissimilarity; a negative value suggests greater average proximity to another group. Their mean is a dimensionless index, not a distance between clusters. The ordinary silhouette is defined for solutions with at least two clusters and fewer clusters than observations.

Among the examined solutions with two to six clusters, the highest mean silhouette is 0.4085 for two clusters. This is a relative comparison among those solutions and within that geometry. A moderate mean does not establish that every state is well assigned; inspect individual silhouettes and differences between groups.

11.2. The gap statistic and its selection rule

The gap statistic, proposed by Tibshirani et al. (2001), compares observed dispersion with dispersion under a reference distribution without cluster structure. Simulating that reference is different from resampling observations to study partition stability. The function clusGap() allows the number of simulations, reference space, and distance power to be specified (Maechler, n.d.).

set.seed(20261113)
gap_usa <- clusGap(usa_z,
  FUNcluster=function(x,k) list(cluster=kmeans(x, centers=k,
                                  nstart=50, iter.max=100)$cluster),
  K.max=6, B=501, d.power=2, spaceH0="scaledPCA", verbose=FALSE)
k_gap <- maxSE(gap_usa$Tab[,"gap"], gap_usa$Tab[,"SE.sim"],
               method="Tibs2001SEmax")
round(gap_usa$Tab, 4)
k_gap  # 2 in the documented run

11.3. What “enhanced analysis” means in eclust()

factoextra::eclust() combines fitting, selection of the number of clusters, and visualization. It returns an object containing results of the chosen method and additional components; this does not make it a second-order tensor. In its documentation, nboot controls gap-statistic simulations when k is estimated. Increasing it can reduce Monte Carlo variability, but does not guarantee consistency, representativeness, or cluster stability (Kassambara & Mundt, n.d.-b).

# Convenient alternative for fitting and visualization:
res_eclust <- factoextra::eclust(usa_z, FUNcluster="kmeans", k=2,
  nstart=100, iter.max=100, graph=FALSE, seed=20261115)
factoextra::fviz_cluster(res_eclust, data=usa_z)
factoextra::fviz_silhouette(res_eclust)
# Here k is fixed: nboot does not enter its selection.
# With k=NULL, specify k.max, nboot, and gap_maxSE, and check
# the implementation settings before comparing results.

12. Hierarchical clustering and dissimilarities

An agglomerative procedure begins with one observation per cluster and merges clusters according to a criterion. The dendrogram records this sequence. Sets are nested across levels, and cutting the tree yields a partition of disjoint groups. Such nesting is different from simultaneous membership in overlapping groups at the same level.

In R, hclust() offers several linkages, and ward.D2 implements Ward’s criterion through a squared treatment of dissimilarities (R Core Team, n.d.-a). Here we supply ordinary Euclidean distances when using ward.D2; we do not square them beforehand. The call eclust(..., "hclust") has its own documented default, hc_method="ward.D2". BIRCH is another algorithm and does not describe this call.

hc_ward <- hclust(dist_usa, method="ward.D2")
ward_groups <- cutree(hc_ward, k=2)
plot(hc_ward, cex=.55, hang=-1)
rect.hclust(hc_ward, k=2, border="saddlebrown")
mean(silhouette(ward_groups, dist_usa)[,3])

# Compare profile similarity using average linkage:
profile_corr <- cor(t(usa_z), method="pearson")
d_cor <- as.dist(1-profile_corr)
hc_cor <- hclust(d_cor, method="average")
corr_groups <- cutree(hc_cor, k=2)
mean(silhouette(corr_groups, d_cor)[,3])

Correlations between rows compare the relative shape of four values already standardized by column. This is not the correlation between the table’s variables, nor does it preserve each profile’s overall magnitude. With only four coordinates per state, this similarity warrants careful interpretation. The result may group similar profiles even when their absolute levels differ.

Completeness means that every Cauchy sequence converges within the space; it is not part of the general definition of a metric. The Minkowski norm on an ℓ_p space and the indefinite tensor of pseudo-Riemannian geometry are distinct constructions. An indefinite metric tensor does not automatically generate a nonnegative distance usable by the algorithms above.

13. PCA and clustering: a controlled comparison

To study what changes under dimensionality reduction, we hold the standardized data, two clusters, 100 initializations, and the algorithm fixed. We fit k-means first with all four variables and then with one, two, three, and four components. We compare labels using the adjusted Rand index (ARI) and calculate silhouettes both in the fitting space and in the original variables.

In this run, the five partitions are identical. However, the silhouette changes from 0.4085 in four variables to 0.6468 when measured only on PC1. The index improves because the geometry changes: assignments did not improve, because they did not change. With all components retained, distances and silhouettes are recovered, as required by a full orthogonal change of coordinates.

ari <- function(a,b) {
  m <- table(a,b)
  c2 <- function(z) z*(z-1)/2
  n <- sum(m); A <- sum(c2(m))
  B <- sum(c2(rowSums(m))); C <- sum(c2(colSums(m)))
  expected <- B*C/c2(n)
  denominator <- (B+C)/2-expected
  if (abs(denominator) < 1e-12) return(1)
  (A-expected)/denominator
}
pca_comparison <- do.call(rbind, lapply(1:4, function(m) {
  Tm <- pca_usa$x[,seq_len(m),drop=FALSE]
  set.seed(20261115)
  fit <- kmeans(Tm, centers=2, nstart=100, iter.max=100)
  data.frame(components=m,
    variance=100*sum(pca_usa$sdev[seq_len(m)]^2)/sum(pca_usa$sdev^2),
    silhouette_reduced=mean(silhouette(fit$cluster, dist(Tm))[,3]),
    silhouette_original=mean(silhouette(fit$cluster, dist_usa)[,3]),
    ARI=ari(km2$cluster, fit$cluster))
}))
round(pca_comparison, 4)

13.1. Retaining much variance can conceal a contrast

Now construct two categories defined by Y=−1 and Y=1, both with X values −10, −6, −2, 2, 6, and 10. Cross-covariance is zero. Variation in X is much greater than variation in Y, so covariance-based PCA selects PC1=X and retains 560/572≈97.90% of the variation. The two categories overlap point for point.

bands <- data.frame(
  X=rep(c(-10,-6,-2,2,6,10), 2),
  Y=rep(c(-1,1), each=6))
pca_bands <- prcomp(bands, center=TRUE, scale.=FALSE)
100*pca_bands$sdev[1]^2/sum(pca_bands$sdev^2)  # 97.9021
cov(bands)
# Standardization would give both columns variance 1 and change
# the geometry: the choice of scale is part of the question.

13.2. Partition stability

Another question is whether assignment withstands perturbations of the data. In 200 subsamples of 40 states, drawn without replacement while keeping the original scaling, k-means was refitted with two groups and compared with the reference labels restricted to those states. Mean ARI was 0.9995, and the 5th, 50th, and 95th percentiles were 1 in the documented run. This is evidence of stability under that particular perturbation, not a population guarantee or causal validation.

set.seed(20261113)
stability <- replicate(200, {
  ix <- sort(sample(nrow(usa_z), 40, replace=FALSE))
  fit <- kmeans(usa_z[ix,], centers=2, nstart=50, iter.max=100)
  ari(km2$cluster[ix], fit$cluster)
})
mean(stability)
quantile(stability, c(.05,.5,.95))
range(stability)

14. Random projections and Johnson–Lindenstrauss

PCA adapts directions to variation in the data. A normalized random projection asks another question: can all distances between a finite number of points be approximately preserved using far fewer coordinates? The Johnson–Lindenstrauss lemma gives an affirmative answer, with dimension depending on the number of points and the tolerated distortion.

Consider N points in a Euclidean space, , and a failure tolerance . Random linear constructions exist for which a dimension of order $m=O({\epsilon }^{-2}\log (N/\delta ))$ suffices to preserve distances simultaneously with probability at least 1−δ. A common formulation controls squared distances:

The historical formulation with factors 1±ε on ordinary distances has the same dimensional order after adjusting constants. The Gaussian construction uses $f(x)=Rx/\sqrt{m}$, with independent standard-normal entries of R. Randomness belongs to the projection matrix; the original observations need not follow a normal distribution (Johnson & Lindenstrauss, 1984; Matoušek, 2013).

Three distinctions are decisive. Dependence is proportional to ε⁻², not merely to the logarithm of ε⁻². The guarantee holds with high probability for the specified construction and settings; it is not an exact isometry with probability one. Finally, preserving distances in a finite set does not mean preserving every direction of the ambient space when dimension is reduced.

set.seed(20261113)
A <- matrix(rnorm(80*800), nrow=80)
R <- matrix(rnorm(120*800), nrow=120)/sqrt(120)
projected <- A %*% t(R)
d_original <- as.vector(dist(A))
d_projected <- as.vector(dist(projected))
ratio <- d_projected/d_original
range(ratio)
max(abs(ratio-1))
plot(d_original, d_projected, pch=16, cex=.35,
     col="saddlebrown", xlab="Original distance",
     ylab="Projected distance")
abline(0,1,col="gray40")

The experiment illustrates the phenomenon; the proof explains how dimension, error, and probability are related. A sufficient theoretical bound may be conservative. The absence of large distortion in one realization does not replace a guarantee for future random realizations.

15. The scope of scientific interpretation

PCA articulates observed relationships among variables and constructs a representation that ranks their variation. In this sense, it can participate in identifying essential relationships of the object: it locates regularities, contrasts, and combinations whose interpretation develops through further conceptual and empirical mediations. Reducing this role to a simple rotation of axes would leave out its scientific use.

At the same time, variance is a particular criterion of relevance. It may capture structure, differences in measurement scale, measurement errors, or noise. A low-variance relationship can be decisive for a specific question, as the band example showed. Identifying essential determinations requires explaining why the variation examined expresses relationships relevant to the phenomenon and how it connects with its history, mechanisms, and other evidence.

Clusters are likewise constructions mediated by a geometry and a criterion. Their stability shows whether certain decisions withstand perturbations; interpretation requires explaining what observations share and why that similarity matters. A group of states defined by arrests and urbanization does not, by itself, constitute a social class, a causal mechanism, or a complete historical typology.

Dimensionality reduction makes a way of preserving relationships explicit. Scientific interpretation requires explaining why those relationships matter for the object being studied.

16. Proofs and extensions

The following extensions preserve the probabilistic and computational development of the original study with consistent notation. They can be read after the applied sections. Algebraic identities are derived here; general theorems and their antecedents are identified in the references.