Is a matrix multiplied with its transpose something special?

The main thing is presumably that $AA^T$ is symmetric. Indeed $(AA^T)^T=(A^T)^TA^T=AA^T$. For symmetric matrices one has the Spectral Theorem which says that we have a basis of eigenvectors and every eigenvalue is real.

Moreover if $A$ is invertible, then $AA^T$ is also positive definite, since $$x^TAA^Tx=(A^Tx)^T(A^Tx)> 0$$

Then we have: A matrix is positive definite if and only if it's the Gram matrix of a linear independent set of vectors.

Last but not least if one is interested in how much the linear map represented by $A$ changes the norm of a vector one can compute

$$\sqrt{\left<Ax,Ax\right>}=\sqrt{\left<A^TAx,x\right>}$$

which simplifies for eigenvectors $x$ to the eigenvalue $\lambda$ to

$$\sqrt{\left<Ax,Ax\right>}=\sqrt \lambda\sqrt{\left<x,x\right>},$$

The determinant is just the product of these eigenvalues.

$AA^T$ is positive semi-definite, and in a case in which $A$ is a column matrix, it will be a rank 1 matrix and have only one non-zero eigenvalue which equal to $A^TA$ and its corresponding eigenvector is $A$. The rest of the eigenvectors are the null space of $A$ i.e. $\lambda^TA = 0$.

Indeed, independent of the size of $A$, there is a useful relation in the eigenvectors of $AA^T$ to the eigenvectors of $A^TA$; based on the property that $rank(AA^T)=rank(A^TA)$. That the rank is identical implies that the number of non-zero eigenvectors is identical. Moreover, we can infer the eigenvectors of $A^TA$ from $AA^T$ and vice versa. The eigenvector decomposition of $AA^T$ is given by $AA^Tv_i = \lambda_i v_i$. In case $A$ is not a square matrix and $AA^T$ is too large to efficiently compute the eigenvectors (like it frequently occurs in covariance matrix computation), then it's easier to compute the eigenvectors of $A^TA$ given by $A^TAu_i = \lambda_i u_i$. Pre-multiplying both sides of this equation with $A$ yields

$AA^TAu_i=\lambda_iAu_i$.

Now, the originally searched eigenvectors $v_i$ of $AA^T$ can easily be obtained by $v_i:=Au_i$. Note, that the resulted eigenvectors are not yet normalized.

One could name some properties, like if $B=AA^T$ then

$$B^T=(AA^T)^T=(A^T)^TA^T=AA^T=B,$$

so

$$\langle v,Bw\rangle=\langle Bv,w\rangle=\langle A^Tv,A^Tw\rangle.$$

Special? Yes! The matrix $A^TA$ is abundant with the Least Squares Finite Element Method in Numerical Analysis:

• Scale vector in scaled pivoting (numerical methods)
• Solving for streamlines from numerical velocity field
• What is the difference between Finite Difference Methods,
  Finite Element Methods and Finite Volume Methods for solving PDEs?
• Any employment for the Varignon parallelogram?

1. Multiply $A$ on the left with $A^T$, giving $B = A^T A$ .
2. The inverse can of $B$ can be determined by employing our special matrix inversion routine.
   The reason is that the pivots of $B$ are always at the main diagonal: see the first reference.
3. The inverse matrix is $B^{-1} = (A^TA)^{-1} = A^{-1}A^{-T}$ .
   Therefore multiply on the right with $A^T$ , giving : $A^{-1}A^{-T}A^T = A^{-1}$ .
   The inverse that has been sought for is recovered herewith.

program Demo;
type
matrix = array of array of double;
procedure Inverteren(var b : matrix);
{
Matrix inversion
pivots on diagonal
}
var
pe : double;
NDM,i,j,k : integer;
begin
NDM := Length(b);
for k := 0 to NDM-1 do
begin
pe := b[k,k];
b[k,k] := 1;
for i := 0 to NDM-1 do
begin
b[i,k] := b[i,k]/pe;
if (i = k) then Continue;
for j := 0 to NDM-1 do
begin
if (j = k) then Continue;
b[i,j] := b[i,j] - b[i,k]*b[k,j];
end;
end;
for j := 0 to NDM-1 do
b[k,j] := - b[k,j]/pe;
b[k,k] := 1/pe;
end;
end;
procedure inverse(b : matrix; var q : matrix);
{
Matrix inversion
without pivoting
}
var
NDM,i,j,k : integer;
p : matrix;
h : double;
begin
NDM := Length(b);
SetLength(p,NDM,NDM);
{ Transpose * original }
for i := 0 to NDM-1 do
begin
for j := 0 to NDM-1 do
begin
h := 0;
for k := 0 TO NDM-1 do
h := h + b[k,i]b[k,j];
p[i,j] := h;
end;
end;
Inverteren(p);
SetLength(q,NDM,NDM);
{ (B^TB)^(-1)*B^T = B^(-1) }
for i := 0 to NDM-1 do
begin
for j := 0 to NDM-1 do
begin
h := 0;
for k := 0 TO NDM-1 do
h := h + p[i,k]*b[j,k];
q[i,j] := h;
end;
end;
end;
procedure test(NDM : integer);
var
i,j,k : integer;
b,q : matrix;
h : double;
begin
SetLength(b,NDM,NDM);
Random; Random;
for i := 0 to NDM-1 do
begin
for j := 0 to NDM-1 do
begin
b[i,j] := Random;
end;
end;
inverse(b,q);
for i := 0 to NDM-1 do
begin
for j := 0 to NDM-1 do
begin
h := 0;
for k := 0 TO NDM-1 do
h := h + q[i,k]*b[k,j];
Write(h:5:2,' ');
end;
Writeln;
end;
end;
BEGIN
test(9);
END.

 1.00  0.00  0.00 -0.00  0.00  0.00  0.00  0.00 -0.00 
 0.00  1.00  0.00 -0.00  0.00  0.00  0.00  0.00 -0.00 
-0.00  0.00  1.00  0.00 -0.00 -0.00 -0.00 -0.00  0.00 
-0.00 -0.00 -0.00  1.00 -0.00 -0.00 -0.00 -0.00  0.00 
 0.00  0.00  0.00  0.00  1.00  0.00  0.00  0.00  0.00 
 0.00  0.00  0.00 -0.00 -0.00  1.00  0.00  0.00 -0.00 
-0.00 -0.00 -0.00 -0.00 -0.00 -0.00  1.00 -0.00  0.00 
-0.00 -0.00 -0.00  0.00 -0.00 -0.00 -0.00  1.00  0.00 
-0.00 -0.00 -0.00 -0.00 -0.00 -0.00 -0.00 -0.00  1.00 

LATE EDIT.
Theory of Polar Decomposition is described in Wikipedia. If attention is restricted to real-valued (non-singular square invertible) matrices, then an appropriate question and some answers are found in Polar decomposition of real matrices. Especially the following formula over there leaves no doubt that a matrix multiplied with its transpose IS something special: $$ B = \left[B\left(B^TB\right)^{-1/2}\right]\left[\left(B^TB\right)^{1/2}\right] $$ Least Squares methods (employing a matrix multiplied with its transpose) are also very useful with Automated Balancing of Chemical Equations