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).
The matrix $A^TA$ contains all the inner products of all columns in $A$. The diagonal thus contains the squared norms of columns. If you think about geometry and orthogonal projections onto the column space spanned by the columns in $A$ you may recall that norms and inner products of the vectors spanning this space play a central role in the computation of the projection. Least squares regression as well as principal components can be understood in terms of orthogonal projections.
Also note that if the columns of $A$ are orthonormal, thus forming an orthonormal basis for the column space, then $A^TA = I$ $-$ the identity matrix.
@NRH gave a good technical answer.
If you want something really basic, you can think of $A^TA$ as the matrix equivalent of $A^2$ for a scalar.
Although it has been already discussed that $\textbf{A}^T\textbf{A}$ has the meaning of taking dot products, I would only add a graphical representation of this multiplication.
Indeed, while rows of the matrix $\textbf{A}^T$ (and columns of the matrix $\textbf{A}$) represent variables, we treat each variable measurements as a multidimensional vector. Multiplying the row $row_p$ of $\textbf{A}^T$ with the column $col_p$ of $\textbf{A}$ is equivalent to taking the dot product of two vectors: $dot(row_p, col_p)$ - the result being the entry at position $(p,p)$ inside the matrix $\textbf{A}^T \textbf{A}$.
Similarly, multiplying the row $p$ of $\textbf{A}^T$ with the column $k$ of $\textbf{A}$ is equivalent to the dot product: $dot(row_p, col_k)$, with the result at position $(p,k)$.
The entry $(p, k)$ of the resulting matrix $\textbf{A}^T\textbf{A}$ has the meaning of how much the vector $row_p$ is in the direction of the vector $col_k$. If the dot product of two vectors $row_i$ and $col_j$ is other than zero, some information about a vector $row_i$ is carried by a vector $col_j$, and vice versa.
This idea plays an important role in Principal Component Analysis, where we want to find a new representation of our initial data matrix $\textbf{A}$ such that, there is no more information carried about any column $i$ in any other column $j \neq i$. Studying PCA deeper, you will see that a "new version" of the covariance matrix is computed and it becomes a diagonal matrix which I leave to you to realize that... indeed it means what I expressed in the previous sentence.
An important view of the geometry of $A'A$ is this (the viewpoint strongly stressed in Strang's book on "Linear Algebra and Its Applications"): Suppose A is an $m \times n$-matrix of rank k, representing a linear map $A: R^n \rightarrow R^m$. Let Col(A) and Row(A) be the column and row spaces of $A$. Then
(a) As a real symmetric matrix, $(A'A): R^n \rightarrow R^n$ has a basis $\{e_1,..., e_n\}$ of eigenvectors with non-zero eigenvalues $d_1,\ldots,d_k$. Thus:
$(A'A)(x_1e_1 + \ldots + x_ne_n) = d_1x_1e_1 + ... + d_kx_ke_k$.
(b) Range(A) = Col(A), by definition of Col(A). So A|Row(A) maps Row(A) into Col(A).
(c) Kernel(A) is the orthogonal complement of Row(A). This is because matrix multiplication is defined in terms of the dot products (row i)*(col j). (So $Av'= 0 \iff \text{v is in Kernel(A)} \iff v \text{is in orthogonal complement of Row(A)}$
(d) $A(R^n)=A(\text{Row}(A))$ and $A|\text{Row(A)}:\text{Row(A)} \rightarrow Col(A)$ is an isomorphism.
Reason: If v = r+k (r \in Row(A), k \in Kernel(A),from (c)) then
A(v) = A(r) + 0 = A(r) where A(r) = 0 <==> r = 0$.
[Incidentally gives a proof that Row rank = Column rank!]
(e) Applying (d), $A'|:Col(A)=\text{Row(A)} \rightarrow \text{Col(A')}=\text{Row(A)}$ is an isomorphism
(f)By (d) and (e): $A'A(R^n) = \text{Row(A)}$ and A'A maps Row(A) isomorphically onto Row(A).
I think it is more intuitive to think of $A^\top A$ as a bilinear form than as a linear map.
Say we have a linear transformation $\mathbb{R}^n \to \mathbb{R}^m$ with matrix $A$ (with respect to the standard basis $\vec{e}_1, \vec{e}_2, \dots \vec{e}_n$).
For an example we can visualize, picture $\mathbb{R}^2 \to \mathbb{R}^3$ given by
$$ A = \begin{bmatrix} 1 & 1\\ 1 & 2\\ 1 & 0 \end{bmatrix} $$
Call the columns
$$ \vec{v}_1 = \begin{bmatrix}1 \\ 1 \\ 1 \end{bmatrix} \,\,\textrm{ and }\,\, \vec{v}_2 = \begin{bmatrix}1 \\ 2 \\ 0 \end{bmatrix} $$
We visualize this linear transformation as a way to "stretch" $\mathbb{R}^2$ onto the plane $H = \textrm{span}(\vec{v}_1, \vec{v}_2)$ by sending $\vec{e}_1 \mapsto \vec{v}_1$ and $\vec{e}_2 \mapsto \vec{v}_2$.
$A^\top A$ is a gadget which lets us compute the dot product of the images of vectors from $\mathbb{R}^2$ in the plane $H$.
The computation is simple. Let $\vec{v}, \vec{w} \in \mathbb{R}^2$ be two vectors in the domain of the transformation.
$$ \begin{align*} (A \vec{v}) \cdot (A \vec{w}) & = (A\vec{w})^\top (A \vec{v})\\ &= \vec{w}^\top A^\top A \vec{v} \end{align*} $$
This allow us to interpret many of the uses of $A^\top A$ in linear algebra and statistics in a nice geometric way.
For instance, let's try to understand how to project a vector $\vec{y}$ into the plane $H$, and how this relates to OLS linear regression.
Let the projected vector be called $\textrm{proj}_H(\vec{y})$. We certainly want
$$ \vec{y} = \textrm{proj}_H(\vec{y}) + \vec{y}^\perp $$
for some vector $\vec{y}$ which is perpendicular to $\vec{H}$.
We also want $\textrm{proj}_H(\vec{y})$ to be in $H$, which is to say, it should be in the image of $A$. Say we have $\textrm{proj}_H(\vec{y}) = A \vec{\beta}$ (here I am thinking of $\beta$ as the vector of coefficients of $\vec{v}_1$ and $\vec{v}_2$).
So we have
$$ \vec{y} = A\vec{\beta} + \vec{y}^\perp $$
Taking $A^\top$ of both sides we have
$$ A^\top \vec{y} = A^\top A \vec{\beta} + \vec{0} $$
since $A^\top$ dots all of the columns of $A$ with $\vec{y}^\perp$, and these should all be zero.
We can directly solve this to get
$$ \beta = (A^\top A)^{-1} (A^\top \vec{y}) $$
which is the "normal equation".
So the actual projection is then $\textrm{proj}_H(\vec{y}) = A \vec{\beta} = A (A^\top A)^{-1}(A^\top \vec{y})$
However, I think it is also instructive to back up a step and think about what this means from a geometric perspective.
Starting with $A^\top \vec{y} = A^\top A \vec{\beta}$ we can dot both sides with $\alpha$ where $\alpha$ is some other potential weighting of the columns. Then we have
$$ \begin{align*} \vec{\alpha}^\top A^\top A \vec{\beta} &= \vec{\alpha}^\top A^\top \vec{y}\\ (A \vec{\alpha})\cdot (A \vec{\beta}) & = (A \vec{\alpha}) \cdot \vec{y} \end{align*} $$
This is saying that, of all the vectors in $H$, $A \vec{\beta}$ is the vector whose dot product with any other vector $\vec{w} = A \vec{\alpha} \in H$ is exactly the same as the dot product of $\vec{y}$ with $\vec{w}$. This is a nice rephrasing of what it means to find the "best approximation" of $\vec{y}$ by a vector in $H$, but which avoids talking about any vectors "outside" of $H$ other than $\vec{y}$.
How does this relate to OLS linear regression?
Say we have explanatory variables $x_1, x_2, ..., x_k$. We have $N$ observations which pair $(x_{1i}, x_{2i}, ..., x_{ki})$ with $y_i$. We want to find $\beta_i$ of the form $y = \beta_0 + \sum_{i = 1}^N \beta_i x_i = $ which fits this data while minimilizing MSE.
We can reframe this as attempting to minimize
$$ |\vec{y} - X\beta|^2 $$
where $\vec{y} \in \mathbb{R}^N$ is the vector with coordinates $y_i$ and $X$ is the matrix whose first column is all $1$, and whose subsequent columns have coordinates $x_{ij}$.
In other words, we are attempting to project $\vec{y}$ into the hyperplane which is the span of $X$! This is exactly the problem we solved above.
Another intuition coming from this interpretation:
Given a matrix $A: \mathbb{R}^n \to \mathbb{R}^m$, consider the Rayleigh quotient map $R: S^{n-1} \to \mathbb{R}$ given by
$$ R(\vec{u}) = \vec{u}^\top A^\top A \vec{u} = |(A \vec{u})|^2 $$
The last expression gives this map a direct geometric interpretation: we take a unit vector, apply $A$ and square its length.
When we apply $A$ to the unit sphere we obtain an ellipsoid in the span of the columns of $A$. So $R(\vec{u})$ is keeping track of the squared norm of this point on an ellipsoid.
The image of critical points of this map represent the vertexes of this ellipsoid.
Imagine $\vec{v} \in \mathbb{S}^{n-1}$. Because of the geometry of the sphere, every tangent vector $\vec{w}$ is perpendicular to $\vec{v}$. Notice that $A \vec{w}$ will still be a tangent vector to the ellipse at $A \vec{v}$. It is geometrically obvious to me (think about an ellipse!) that $A(\vec{v}) \perp A(\vec{w})$ is zero for all $\vec{w} \in \vec{v}^\perp$ if and only $A\vec{v}$ is a vertex of the ellipsoid.
In other words, the tangent plane to an ellipsoid is only perpendicular to a radial vector at the vertexes.
So the critical points $\vec{v}$ of $R$ satisfy
$$ A(\vec{v}) \cdot A(\vec{w}) = 0 \textrm{ for all } \vec{w} \in \vec{v}^\perp $$
This allows us to easily see that the critical points of $R$ are eigenvectors of $A^\perp A$. Let $\vec{v}$ be a critical point of $R$. Write $A^\perp A (\vec{v})$ in the form $\lambda \vec{v} + \vec{w}$ for $\vec{w} \in \vec{v}^\perp$. Dotting both sides with $\vec{w}$ gives us:
$$ \begin{align*} \vec{w}^\top A^\top A \vec{v} &= \vec{w}^\top (\lambda \vec{v} + \vec{w})\\ (A \vec{w}) \cdot (A\vec{v}) &= |\vec{w}|^2 \\ |\vec{w}|^2 &= 0 \end{align*} $$
So $\vec{v}$ is an eigenvector of $A^\top A$. This is equivalent to the usual derivation with Lagrange Multipliers. The only difference is that I replaced the differential calculus with the geometric fact that the radial vector of an ellipsoid is perpendicular to the tangent plane only at the vertexes.
So we can interpret the eigenvectors of $A^\top A$ as being the critical points of the map $R$ which has a very natural geometric interpretation. For instance, the interpretation of the eigenvector $\vec{v}_1 $with the greatest eigenvalue $\lambda_1$ is "$A(\vec{v}_1)$ is the vector on the ellipsoid $\textrm{Image}(S^{m-1})$ which is furthest from the origin". This is the first principle component in PCA. The other principle components are the other axes of this ellipsoid!
How does this relate to statistics?
Using the same notation as above let $\vec{X}$ as our design matrix and $\vec{y}$ be our observed outcome vector. Let the orthonormal eigenvectors of $X^\top X$ be $\vec{v}_i$ with eigenvalues $\lambda_i$ with this sequence of eigenvalues decreasing in value. Then we know that we want
$$ X^\top \vec{y} = X^\top X \vec{\beta} $$
Writing $\vec{\beta}$ in the basis of eigenvectors we have
$$ X^\top \vec{y} = \sum_{i=1}^{k+1} \alpha_i \vec{v}_i $$
If we eliminate all of the terms from this sum with "small" eigenvalues $|\alpha_i| \leq \epsilon$ then the new $\vec{\beta}'$ will live in the lower dimensional subspace generated by only the "large" eigenvectors, but as we can see $X^\top X \vec{\beta}'$ will still be "close" to $X^\top \vec{y}$. In this sense, we could consider projecting $\vec{y}$ onto first just the largest eigenvector, then the first two, then the first three, etc. As we do so the projection will approach the true regression, but we have diminishing returns with smaller eigenvalues. It is interesting that the highest eigenvalue terms are carrying "most of the information" about the regression.
Note: To train your geometric intuition a bit more, consider the following three ellipsoids: The image of $S^2$ under:
In (3) we have eigenvalues $3$, $2$, $1$ in the directions $e_1$, $e_2$, $e_3$ and these are our only choices.
In (2) we have eigenvalues $3$ and $1$. Our eigenspace for $\lambda = 1$ is 2 dimensional and is spanned by $e_2$ and $e_3$.
In (1) we have eigenvalue $1$ and the eigenspace is the entire domain!
For a "generic map" $A$ we should expect that the ellipsoid $\textrm{Im}(A)$ will be an ellipsoid whose axis lengths are all distinct, so $A^\top A$ will usually have a diagonalization. However it is possible (though extremely rare for real data) that two of the axes could have the same length leading to higher dimensional eigenspaces. This would give us some choice in the selection of an orthogonal basis for this eigenspace.
There are levels of intuition. For those familiar with matrix notation instatistics the intuition is to think of it as a square of the random variable: $x\to E[x^2]$ vs $A\to A^TA$
In matrix notation a sample of the random variable $x$ observations $x_i$ or a population are represented by a column vector: $$a=\begin{bmatrix} x_1 \\ x_2 \\ \dots \\ x_n \end{bmatrix}$$
So, if you want to get a sample mean of the square of the variable $x$, you simply get a dot product $$\bar{x^2}=\frac{a\cdot a} n$$, which is the same in matrix notation as $A^TA$.
Notice, that if the sample mean of the variable is ZERO, then the variance is equal to the mean of the square: $\sigma^2=E[x^2]$ which is analogous to $A^TA$. This is the reason why in PCA you need the zero mean, and why $A^TA$ shows up, after all PCA is to decompose the variance matrix of the data set.
