This is essentially the same calculation I outlined in this other answer (though it might not be immediately obvious).$\newcommand{\ket}[1]{\lvert#1\rangle}\newcommand{\ketbra}[2]{\lvert#1\rangle\!\langle#2\rvert}$
Let us denote with $\Pi_Y$ and $\Pi_N=I-\Pi_Y$ the projectors onto the "yes space" and the "no space".
Given an initial state $\ket\psi$, the goal is getting as close to a state in $\Pi_Y$ as possible, as fast as possible.
Because $\{\Pi_Y,\Pi_N\}$ define a separation of the full space, any state can be decomposed using these operators. In particular, we can write
$$\ket\psi=\cos\theta\ket\alpha+\sin\theta\ket\beta,$$
where $\theta$ is defined via $\cos\theta\equiv\|\Pi_Y\ket\psi\|$, while $\ket\alpha\equiv\Pi_Y\ket\psi/\cos\theta$ and $\ket\beta\equiv\Pi_N\ket\psi/\sin\theta$ (I'm using a slightly different notation with the $\theta$ than the one in the OP, sorry about that).
Intuitively, I'm writing $|\psi\rangle$ as a superposition of its projections onto "yes" and "no" spaces, defining $|\alpha\rangle,|\beta\rangle$ in such a way to absorb potential phase terms.
The Grover iterator is defined as $G=-S_\psi S_Y$, where $S_Y$ and $S_\psi$ are reflections in state space, that is, operators which leave untouched some subspace and change the sign on everything else.
More specifically $S_Y$ flips the "yes space", while $S_\psi$ flips the direction corresponding to the initial state $\ket\psi$ (that is, it leaves the direction of the initial state untouched and flips everything else).
Mathematically, these reflections can be written as
$$S_Y\equiv I - 2\Pi_Y = \Pi_N-\Pi_Y, \qquad
S_\psi\equiv I - 2\ket\psi\!\langle\psi\rvert.$$
It follows that the Grover operator reads
$$G=(I-2\ket\psi\!\langle\psi\rvert)(2\Pi_Y-I).$$
Expanding this product we get
$$G=2\Pi_Y-I-4 \lvert\psi\rangle\!\langle\psi\rvert\Pi_Y + 2\ket\psi\!\langle\psi\rvert.$$
Expanding $\ket\psi$ in terms of $\ket\alpha$ and $\ket\beta$, and remembering the property of $\ket\alpha$ that $\Pi_Y\ket\alpha=\ket\alpha$, you can readily verify that this expression becomes, after a bit of algebra, the following (let me use here the shorthand notation $c\equiv\cos\theta$ and $s\equiv\sin\theta$):
$$G=2\Pi_Y-I+2s^2\ketbra\beta\beta-2c^2 \ketbra\alpha\alpha +2cs (\ketbra\alpha\beta-\ketbra\beta\alpha).$$
Now this represents a rotation in state space of some angle with respect to some axis, but because we are only interested in the action of $G$ on states spanned by the $\ket\alpha$ and $\ket\beta$ states, we need only analyse how $G$ acts on these two states. We then readily get:
\begin{align}
G\ket\alpha&=-\cos(2\theta)\ket\alpha-\sin(2\theta)\ket\beta, \\
G\ket\beta&=\phantom{-}\sin(2\theta)\ket\alpha-\cos(2\theta)\ket\beta. \\
\end{align}
Collecting the corresponding amplitudes in a matrix, we conclude that the action of $G$ in the space spanned by $\ket\alpha$ and $\ket\beta$ can be represented as
$$G\doteq\begin{pmatrix}-\cos(2\theta) &\sin(2\theta)\\-\sin(2\theta) & -\cos(2\theta)\end{pmatrix}.$$
Note that this shows that the result holds in a more general scenario that the one often used when first explaining Grover's algorithm.
You can however easily reduce to the standard situation (which you also use in your post) by having $\ket\psi$ be a balanced superposition of all basis states, and $\Pi_Y$ a trace-1 projector (that is, a projector over a one-dimensional subspace).
Another derivation with different notation
Let us use the shorthand notation $\psi\equiv\ketbra\psi\psi$, $Y\equiv\Pi_Y$ and $N\equiv\Pi_N$. We then have
\begin{align}G = (I-2\psi)(2Y-I)
= &Y(1-2\psi)Y - 2N\psi Y \\
& +2 Y\psi N -N(1-2\psi)N.
\end{align}
Then, specialising to the case $\ket\psi=\cos(\theta)\ket y+\sin(\theta)\ket n$, we get
$$G = Y(I - 2\cos^2(\theta) \ketbra y y)Y + \sin(2\theta)(\ketbra y n - \ketbra n y) \\
- N(I- 2 \sin^2(\theta) \ketbra n n)N,$$
and in particular its action on $\ket y$ and $\ket n$ is
\begin{align}
G\ket y &= -\cos(2\theta)\ket y - \sin(2\theta)\ket n, \\
G\ket n &= \phantom{-}\sin(2\theta) \ket y - \cos(2\theta) \ket n.
\end{align}
