Find the shortest distance of a point to the solution set of a linear equation Given an affine set
$$
\mathbf{S}=\left\{ \mathbf{x|Ax}=\mathbf{y} \right\} \subseteq \mathbb{R}^N
$$
where the matrix $\mathbf{A}\in\mathbb{R}^{M\times N} (M\ll N)$ satisfying the orthogonal form: $\mathbf{AA}^\top=\mathbf{I}_M,$ and vector $\mathbf{y}\in\mathbb{R}^{M}$ are both known.
Question: For any point $\mathbf{z}\in\mathbb{R}^N$, find $\mathbf{r}\in\mathbf{S}$ with a shortest Euclidean distance $l$ from $\mathbf{z}$.

My Efforts:
My main concern is the non-trivial case where the linear equation $\mathbf{Ax}=\mathbf{y}$ has infinite solutions.
Firstly, I turn the above question into the following optimization problem:
$$
\underset{\mathbf{r}}{\min}f\left( \mathbf{r};\mathbf{A},\mathbf{y},\mathbf{z}, \lambda \right),
$$
where
$$
f\left( \mathbf{r};\mathbf{A},\mathbf{y},\mathbf{z}, \lambda \right) =\frac{1}{2}\lVert \mathbf{r}-\mathbf{z} \rVert _{2}^{2}+\frac{\lambda}{2}\lVert \mathbf{Ar}-\mathbf{y} \rVert _{2}^{2}, \lambda\in \mathbb{R}^+.
$$
Let
$$
\frac{\partial f}{\partial \mathbf{r}}=\left( \mathbf{r}-\mathbf{z} \right) +\lambda \mathbf{A}^{\top}\left( \mathbf{Ar}-\mathbf{y} \right) =\mathbf{0}_N,
$$
and I got:
$$
\mathbf{r}=\left( \mathbf{I}_N+\lambda \mathbf{A}^{\top}\mathbf{A} \right) ^{-1}\left( \lambda \mathbf{A}^{\top}\mathbf{y}+\mathbf{z} \right),
$$
then I actually do not know what to do next, and if my idea was correct.

Update:
I notice that:
(1) when $\lambda \rightarrow 0$, we have $\mathbf{r} \rightarrow \mathbf{z}$, and
(2) when $\lambda \rightarrow +\infty$, we have $\mathbf{r} \rightarrow \left( \mathbf{A}^{\top}\mathbf{A} \right) ^{-1} \mathbf{A}^{\top}\mathbf{y}$.
I can understand (1), but am confused by (2).
After a long struggle, I give up solving this problem and guess the solution is:
$$
\mathbf{r}=\mathbf{A}^\top\mathbf{y}+(\mathbf{I}_N-\mathbf{A}^{\top}\mathbf{A})\mathbf{z}.
$$
 A: It seems that I solved this problem by myself.

My Solution:
For any $\mathbf{x}\in \mathbb{R}^N$, there exists a unique range-nullspace decomposition (RND):
$$
\mathbf{x}=\mathbf{A^\top Ax}+(\mathbf{I}_N-\mathbf{A^\top A})\mathbf{x},
$$
so for any $\mathbf{s}\in \mathbf{S}$, we have $\mathbf{s}=\mathbf{A^\top y}+(\mathbf{I}_N-\mathbf{A^\top A})\mathbf{s}$, and $\mathbf{z}=\mathbf{A^\top Az}+(\mathbf{I}_N-\mathbf{A^\top A})\mathbf{z}$.
Let
$$
f(\mathbf{s};\mathbf{z})=\frac{1}{2}\lVert \mathbf{s}-\mathbf{z}\rVert^2_2\\
=\frac{1}{2}\lVert \mathbf{A}^\top(\mathbf{y}-\mathbf{Az})+(\mathbf{I}_N-\mathbf{A^\top A})(\mathbf{s}-\mathbf{z})\rVert^2_2,
$$
and
$$
\frac{\partial f}{\partial \mathbf{s}}=(\mathbf{I}_N-\mathbf{A^\top A})(\mathbf{s}-\mathbf{z}) = \mathbf{0}_N.
$$
The above equation is the necessary condition for an optimal point, so assume $\mathbf{r}$ is a solution, we have
$$
\mathbf{r}=\mathbf{A^\top y}+(\mathbf{I}_N-\mathbf{A^\top A})\mathbf{r}\\
=\mathbf{A^\top y}+(\mathbf{I}_N-\mathbf{A^\top A})(\mathbf{r}-\mathbf{z}+\mathbf{z})\\
=\mathbf{A^\top y}+(\mathbf{I}_N-\mathbf{A^\top A})\mathbf{z},
$$
and the shortest Euclidean distance
$$
l=\lVert \mathbf{r}-\mathbf{z}\rVert_2=\lVert \mathbf{y}-\mathbf{Az}\rVert_2,
$$
i.e., the final unique solution is: $\mathbf{r}=\mathbf{A^\top y}+(\mathbf{I}_N-\mathbf{A^\top A})\mathbf{z}$, and $l=\lVert \mathbf{y}-\mathbf{Az}\rVert_2$.
