Why does $A^TAx = A^Tb$ have infinitely many solution algebraically when $A$ has dependent columns? 
*

*This is a problem from least square approximation, where we solve the equation $A^TAx = A^Tb$ when $Ax = b$ is unsolvable.

*The case I am dealing with is when A has dependent columns, i.e. A is an m by n matrix where the rank r is smaller than n. In this case $A^TA$ is a singular n by n matrix with dependent columns, the rank of which is also r (Rank($A^TA$)=Rank($A$)).

*Now in the book Introduction to Linear Algebra of the legendary Gilbert Strang, he says and I quote, when A is singular, $A^TA$ is also singular, and the equation $A^TAx = A^Tb$ had infinitely many solutions, the pseudoinverse gives us a way to choose a "best solution" $x^+=A^+b$.

*I understand why the equation has infinitely many solutions geometrically:
Because what the equation asserts geometrically is to find the projection of b, denoted by p, in the column space of A, then solve the new equation $A\hat x$ = p. Because we can always project b onto the column space of A, whether it's singular or not, we know there must be a solution to the equation $A\hat x$=p and if there is a solution, there are infinitely many because A is singular.

*My question is how do we know that the equation have infinitely many solution algebraically, to make it clearer, I don't understand why the equation has at least one solution. I do understand that once it has at least one solution, it has infinitely many.

*Algebraically, I understands that $A^Tb$ will take us to $C(A^T)$, and it will take away the part of b that lies in $N(A)$. But what does it has to do with   $C(A^TA)$ ? My hypothesis is there is some formula regarding $C(A^TA)$ and $C(A^T)$ that I am not aware of. For example, if $C(A^TA) = C(A^T)$, then my problem is solved.

*Also, I found this How come least square can have many solutions?, I know what  $\hat x^TA\hat x$ in sums will looks like, but I don't know where it came from, but assuming that this is actually correct, I understand the arguments made in that thread. 

*Any instruction will be appreciated. 

 A: Let $A$ be an $n \times m$ matrix. 
Consider 
$$T_{A^T} : \mathbb R^n \to \mathbb R^m , T_{A^T}=A^Tx \\
T_{A} : \mathbb R^m \to \mathbb R^n , T_{A}=Ax \\
T_{A^TA} : \mathbb R^m \to \mathbb R^m , T_{A^TA}=A^TAx $$
Since $T_{A^TA}=T_{A^T}\circ T_A$ we have
$$C(A^TA)=range(T_{A^TA}) \subset range(T_{A^T})= C(A^T)$$
If we can show that these two spaces have the same dimension we are done.
But this comes for free from the rank nullity theorem and the fact that 
$$ker(T_{A^TA})=\ker(T_A)$$
Indeed, the $\supset$ inclusion is obvious. For $\subset$ let $x \in ker(T_{A^TA})$, then 
$$A^TAx=0 \Rightarrow x^TA^TAx=0 \Rightarrow (Ax)^T(Ax)=0 \Rightarrow (Ax)\cdot(Ax)=0 \Rightarrow \| Ax\|^2=0 \Rightarrow Ax=0$$ 
A: Existence
First and foremost, we must have existence before we can talk about uniqueness. To have existence, we require that the data vector $b$ is not in the $\color{red}{null}$ space:
$$
 b\notin\color{red}{\mathcal{N} \left( \mathbf{A} \right)}
$$
Uniqueness
When the matrix $\mathbf{A}$ has dependent columns, the $\color{red}{null}$ space $\color{red}{\mathcal{N} \left( \mathbf{A} \right)}$ is nontrivial. Let $\color{red}{z}$ be any vector in the $\color{red}{null}$ space, and let $x_{LS}$ be a known least squares minimizer.
$$
\mathbf{A} \left( x_{LS} + \color{red}{z} \right) = \mathbf{A} x_{LS} +  \mathbf{A} \color{red}{z} = \mathbf{A} x_{LS} + \mathbf{0} = \mathbf{A} x_{LS}
$$
Under the action of $\mathbf{A}$, the solutions $\left( x_{LS} + \color{red}{z} \right)$ and $x_{LS}$ are equivalent.
Geometrically, the issue is shown below. The general least squares solution is the affine space represented by the red, dashed line. 


*

*If the $\color{red}{null}$ space is trivial, $\color{red}{\mathcal{N} \left( \mathbf{A} \right)} = \mathbf{0}$, the least squares solution is the point $\color{blue}{\mathbf{A}^{+}b}$ in the $\color{blue}{range}$ space $\color{blue}{\mathcal{R} \left( \mathbf{A}^{*} \right)}$.

*If the $\color{red}{null}$ space  is not trivial, $\color{red}{\mathcal{N} \left( \mathbf{A} \right)} \ne \mathbf{0}$, the least squares solution is the affine space going through the point $\color{blue}{\mathbf{A}^{+}b}$ in $\color{blue}{range}$ space $\color{blue}{\mathcal{R} \left( \mathbf{A}^{*} \right)}$ and extending through the $\color{red}{null}$ space.


Least squares solution
Given a matrix $\mathbf{A}\in\mathbb{C}^{m\times n}_{\rho}$, and a data vector $b\in\mathbb{C}^{m}$ such that
$$
 b\notin\color{red}{\mathcal{N} \left( \mathbf{A} \right)}
$$
The least squares solution is defined as 
$$
 x_{LS} = \left\{
 x\in\mathbb{C}^{n} \colon
\lVert
 \mathbf{A} x - b
\rVert_{2}^{2}
\text{ is minimized}
\right\}
$$
The least squares solution is computed using
$$
 x_{LS} = 
\color{blue}{\mathbf{A}^{+} b} +
\color{red}{ 
\left(
\mathbf{I}_{n} - \mathbf{A}^{+} \mathbf{A}
\right) y}, \quad y \in \mathbb{C}^{n}
$$
In this form, it's clear that having a unique solution demands 
$$
 \mathbf{A}^{+} \mathbf{A} = \mathbf{I}_{n}
$$
which happens when $\color{red}{\mathcal{N} \left( \mathbf{A} \right)} = \mathbf{0}$. More details are in the subsequent links.

Explore Stack Exchange:

Existence and uniqueness of least squares solutions: Query about the Moore Penrose pseudoinverse method
Derivation of the SVD solution and conditions on existence and uniqueness: Singular value decomposition proof
Subspace decomposition for least squares: Singular Value Decomposition
How the two null spaces affect the least squares solution: Pseudo-inverse of a matrix that is neither fat nor tall?, What forms does the Moore-Penrose inverse take under systems with full rank, full column rank, and full row rank?
How full column rank changes the inverse: How to find the singular value decomposition of $A^{T}A$ & $\left( A^{T}A \right)^{-1}$
How null spaces affect the pseudoinverse: generalized inverse of a matrix and convergence for singular matrix, When pseudo inverse and general inverse of a invertible square matrix will be equal or not equal?
A: Forgetting how we got the equation (by transposing and multiplying), and just looking at it in its own right, we get this:
$c = A^Tb$ is a column vector and $B = A^TA$ is a square matrix with linearly dependent columns. Lastly, we know that $Bx = c$ has at least one solution.
Algebraically, $Bx$ is a linear combination of the columns of $B$. There are infinitely many of these linear combinations that give $0$, so for each linear combination of columns there are infinitely many others that give the same result (just add one of the infinitely many linear combinations that give $0$). That includes an arbitrary linear combination that gives $c$, i.e. a solution to $Bx = c$.
A: Imagine you have two data points (2,4) and (2,6). Now, if we put the slope formula into matrix form, we will have
\begin{pmatrix}2 & 1 \\2 & 1 \\\end{pmatrix}
as our X matrix. First column is for the input X and second for the slope "b". We set "b" as one because we want to give it freedom to move (otherwise, 0 would kill the "b" out of the formula, mx = y). Then, we have the vector (M B) for our unknown coefficients. These two matrices multiplied will give the output vector (4 6).
Now, try to think of how you would graph through these two data points. The problem is that there are infinite solutions depending on what you set the bias "b" to be equal to. Why do we have infinite least square solutions? Well, The column space is not full rank and we have a dependency.
Now, b is a free variable then. The algebraic intuition would come from setting b with different values to obtain the equivalent solution.
