For a given data matrix $A$ (with variables in columns and data points in rows), it seems like $A^TA$ plays an important role in statistics. For example, it is an important part of the analytical solution of ordinary least squares. Or, for PCA, its eigenvectors are the principal components of the data.

I understand how to calculate $A^TA$, but I was wondering if there's an intuitive interpretation of what this matrix represents, which leads to its important role?

Geometrically, matrix $\bf A'A$ is called matrix of scalar products (= dot products, = inner products). Algebraically, it is called sum-of-squares-and-cross-products matrix (SSCP).

Its $i$-th diagonal element is equal to $\sum a_{(i)}^2$, where $a_{(i)}$ denotes values in the $i$-th column of $\bf A$ and $\sum$ is the sum across rows. The $ij$-th off-diagonal element therein is $\sum a_{(i)}a_{(j)}$.

There is a number of important association coefficients and their square matrices  are called angular similarities or SSCP-type similarities:


Dividing SSCP matrix by $n$, the sample size or number of rows of $\bf A$, you get MSCP (mean-square-and-cross-product) matrix. The pairwise formula of this association measure is hence $\frac{\sum xy}{n}$ (with vectors $x$ and $y$ being a pair of columns from $\bf A$).
If you center columns (variables) of $\bf A$, then $\bf A'A$ is the scatter (or co-scatter, if to be rigorous) matrix and $\mathbf {A'A}/(n-1)$ is the covariance matrix. Pairwise formula of covariance is $\frac{\sum c_xc_y}{n-1}$ with $c_x$ and $c_y$ denoting centerted columns.
If you z-standardize columns of $\bf A$ (subtract the column mean and divide by the standard deviation), then $\mathbf {A'A}/(n-1)$ is the Pearson correlation matrix: correlation is covariance for standardized variables. Pairwise formula of correlation is $\frac{\sum z_xz_y}{n-1}$ with $z_x$ and $z_y$ denoting standardized columns. The correlation is also called  coefficient of linearity.
If you unit-scale columns of $\bf A$ (bring their SS, sum-of-squares, to 1), then $\bf A'A$ is the cosine similarity matrix. The equivalent pairwise formula thus appears to be $\sum u_xu_y = \frac{\sum{xy}}{\sqrt{\sum x^2}\sqrt{\sum y^2}}$ with $u_x$ and $u_y$ denoting L2-normalized columns. Cosine similarity is also called coefficient of proportionality.
If you center and then unit-scale columns of $\bf A$, then $\bf A'A$ is again the Pearson correlation matrix, because correlation is cosine for centered variables$^{1,2}$: $\sum cu_xcu_y = \frac{\sum{c_xc_y}}{\sqrt{\sum c_x^2}\sqrt{\sum c_y^2}}$


Alongside these four principal association measures let us also mention some other, also based on of $\bf A'A$, to top it off. They can be seen as measures alternative to cosine similarity because they adopt different from it normalization, the denominator in the formula:


Coefficient of identity [Zegers & ten Berge, 1985] has its denominator in the form of arithmetic mean rather than geometric mean: $\frac{\sum{xy}}{(\sum x^2+\sum y^2)/2}$. It can be 1 if and only if the being compared columns of $\bf A$ are identical.
Another usable coefficient like it is called similarity ratio: $\frac{\sum{xy}}{\sum x^2 + \sum y^2 -\sum {xy}} = \frac{\sum{xy}}{\sum {xy} + \sum {(x-y)^2}}$.
Finally, if values in $\bf A$ are nonnegative and their sum within the columns is 1 (e.g. they are proportions), then $\bf \sqrt {A}'\sqrt A$ is the matrix of fidelity or Bhattacharyya coefficient.




$^1$ One way also to compute correlation or covariance matrix, used by many statistical packages, bypasses centering the data and departs straight from SSCP matrix $\bf A'A$ this way. Let $\bf s$ be the row vector of column sums of data $\bf A$ while $n$ is the number of rows in the data. Then (1) compute the scatter matrix as $\bf C = A'A-s's/ \it n$ [thence, $\mathbf C/(n-1)$ will be the covariance matrix]; (2) the diagonal of $\bf C$ is the sums of squared deviations, row vector $\bf d$; (3) compute correlation matrix $\bf R=C/\sqrt{d'd}$.

$^2$ An acute but statistically novice reader might find it difficult reconciling the two definitions of correlation - as "covariance" (which includes averaging by sample size, the division by df="n-1") and as "cosine" (which implies no such averaging). But in fact no real averaging in the first formula of correlation takes place. The thing is that st. deviation, by which z-standardization was achieved, had been in turn computed with the division by that same df; and so the denominator "n-1" in the formula of correlation-as-covariance entirely cancels if you unwrap the formula: the formula turns into the formula of cosine. To compute empirical correlation value you really need not to know $n$ (except when computing the mean, to center).
