Numerical simulations suggest that the expression $$ A=G^\dagger (PGG^\dagger P)^+G, $$ where $^+$ denotes the Moore-Pensore pseudoinverse, $P$ is the projector $$ P=I-\frac{1}{c^\dagger G^\dagger G c} Gcc^\dagger G^\dagger $$ where $c$ is some nonzero vector, and $G$ has orthogonal columns, seems to simplify to $$ A=I-\frac{1}{c^\dagger c}cc^\dagger,$$ but all my attempts failed in proving it. I discuss below some ideas.
First, numerically it seems that $A$ does not depend on the norm of the columns of $G$; if we were able to prove this, then we could restrict to unitary matrices, $G^\dagger G=I$. But actually I didn't succeed in proving that. Second, assume that all elements of $c$ are nonzero; then I observe that the rank of $PGG^\dagger P$ is $n-1$, where $\mathrm{rank\,}G=n$. So the pseudoinverse is “almost an inverse.” Third, if $P$ were $I$, then $A=I$, therefore the term $-(cc^\dagger)/\|c\|^2$ is there because of $P$.
Any suggestion is very appreciated.
Answer: Denote $G=QD^{1/2}$ and rewrite $PGG^\dagger P=QAQ^\dagger$ for some $A$. From the properties of MP pseudoinverse follows that $$ G^\dagger( PGG^\dagger P )^+ G = D^{1/2} A^+ D^{1/2}. $$ The pseudoinverse of $A$ is given by $$ A^+ = D^{-1/2}\left(I-\frac{1}{\|c\|^2}cc^\dagger \right) D^{-1/2}, $$ that satisfies the four properties defining the MP pseudoinverse.