Your proof that $Q$ is an orthogonal projection matrix is correct. One geometric interpretation of the vectors $z=Pb$ and $v=Qb$ is exactly as user84413 stated: orthogonal projections onto the range of P and onto the orthogonal complement of the range of P.
Another geometric interpretation I use for students just starting to learn this topic is as follows.
1. P and Q divide our space into two orthogonal spaces
Any vector in $P$ is orthogonal to any vector in $Q$.
For example, if our space was the Cartesian X-Y plane, perhaps P is analogous to the X axis and Q is the Y axis. The X and Y axes are orthogonal, so any vector on the X axis (e.g., $(3,0)$) is orthogonal to any vector on the Y axis (e.g., $(0,7)$).
2. Any vector in our space can be written as a sum of a vector in P and a vector in Q
For any vector $b$, there are two 'component' vectors in $P$ and $Q$ whose sum is vector $b$. Let $z$ be the vector in space $P$ and $v$ be the vector in space $Q$, then
$$b = z + v$$
which is
$$b = Pb + Qb$$
For example, any vector in the X-Y plane can be represented as a sum of:
- a vector that lies solely on the X axis, and
- a vector that lies solely on the Y axis
In the X-Y plane, the vector $b=(3,7)$ can be written as the sum of $z=(3,0)$ (which is on the X-axis) and $v=(0,7$) (which is on the Y-axis). The vector $b=(3,7)$ is like the hypotenuse of a right triangle that has the two legs $z=(3,0)$ and $v=(0,7)$.
The vectors $z$ and $v$ are unique. In other words, once we are given the orthogonal subspaces P and Q there is only one way to get the two components of the vector $b$ that lie in $P$ and $Q$.
3. The 2-norm of any vector in our space can be calculated from the norms of the component vectors in P and Q
$$ \left \lVert b \right \rVert^2 = \left \lVert z \right \rVert^2 + \left \lVert v \right \rVert^2 $$
which is
$$ \left \lVert b \right \rVert^2 = \left \lVert Pb \right \rVert^2 + \left \lVert Qb \right \rVert^2 $$
For example, in the X-Y plane the vector $b=(3,7)$ has a 'length' that can be found by the Pythogorean theorem:
$$ \left \lVert (3,7) \right \rVert^2 = \left \lVert (3,0) \right \rVert^2 + \left \lVert (0,7) \right \rVert^2 $$
$$ \left \lVert (3,7) \right \rVert^2 = 3^2 + 7^2 $$
Say the columns of $A$ are $A_1,\dots,A_n$.
If $B = (b_{ij}) = A^TA$, then:
$b_{ij} = \langle (A_i)^T,A_j\rangle$, and since these (the $A_k$) are orthonormal relative to our inner product, we have:
$b_{ij} = \delta_{ij}$ (the kronecker delta), which is $1$ when $i = j$, and $0$ otherwise, that is to say $B$ is the identity matrix.
So $A^TA = I$, from which we conclude $A^T = A^{-1}$ (If you insist on showing $A$ has a two-sided inverse, see below).
On the other hand, if $A^TA = I$, then (running our argument in reverse), we see the columns of $A$ (and thus the rows of $A^T$) form an orthonormal basis (they form a basis since $A$ is invertible).
By considering $AA^T = I$ in the same way, we see the columns of $A^T$, and thus the rows of $(A^T)^T = A$ also form an orthonormal basis.
Best Answer
Following up on @Augustin suggestion, write
$$(Px)\cdot (Py)=(Px)^T\cdot (Py)=x^T(P^TP)y$$
for arbitrary $x$ and $y$. So if $P^TP=I$ one has $\forall x,y \,(Px)\cdot (Py)=x\cdot y$
Assume now that the last equality holds. Consider a basis $(e_i)$ and compute
$$e_i^TP^TPe_j=(P^TP)_{ij}=(Pe_i)\cdot (Pe_j)=e_i\cdot e_j=\delta_{ij}$$
And $P^TP=I$