The question is
Find matrix $R$ so that $R\vec{a} = \vec{a} - \frac{\vec{a} \cdot \vec{b}}{\|\vec{b}\|^2} \vec{b} = \vec{a} - \vec{a}_{\parallel\vec{b}} = \vec{a}_{\bot \vec{b}}$.
I want to check if my solution is correct.
For start, $\vec{a}_{\bot \vec{b}} = \frac{\vec{b} \times \vec{a} \times \vec{b}}{\|b\|^2}$ is something I proved previously using the vector triple product and its property according to which $\vec{b} \times (\vec{a} \times \vec{b}) = (\vec{b} \times \vec{a}) \times \vec{b}$; another property is $(\vec{b} \times \vec{a}) \times \vec{c} = -\vec{c} \times (\vec{b} \times \vec{a})$, which in our case gives $-\vec{b} \times (\vec{b} \times \vec{a})$. Knowing that $\vec{u} \times \vec{v}$ can be written under the form of a matrix as $[u]_{\times}\vec{v}$ with $[u]_{\times} = \begin{bmatrix} 0 & -u_z & u_y \\ u_z & 0 & -u_x \\ -u_y & u_x & 0 \\ \end{bmatrix}$ we get $-\vec{b} \times (\vec{b} \times \vec{a}) = [-b]_{\times}([b]_{\times}\vec{a}) = [-b]_{\times}[b]_{\times}\vec{a}$. Considering the fact that $\vec{a}_{\bot \vec{b}} = \frac{\vec{b} \times \vec{a} \times \vec{b}}{\|b\|^2}$, we get $R = \frac{[-b]_{\times}[b]_{\times}}{\|b\|^2}$. Is my reasoning correct here or am I missing something?
First, the quoted formula is incorrect. It should have $\|\vec b\|^2$ in the denominator—you have to normalize $\vec b$ to compute the projection of $\vec a$ onto $\vec b$. Your identity $$\vec a_{\perp\vec b} = {\vec b\times\vec a\times\vec b \over \|\vec b\|}$$ suffers from the same problem: it should also have $\|\vec b\|^2$ in the denominator, and for the same reason: you need to use the normalized $\vec b/\|\vec b\|$ in both places.
Anyway, once you’ve corrected for the missing normalization, your derivation works (and in fact turns up in the derivation of Rodrigues’ rotation formula). However, it’s only valid in spaces that have a cross product, which limits its applicability. A more general solution can be written directly from the given (corrected) formula, after a bit of algebraic manipulation: $$R\vec a = \vec a - {\vec a \cdot \vec b \over \|\vec b\|^2} \vec b = \vec a - {1\over \|\vec b\|^2}\left(\vec b^T\vec a\right)\vec b = \vec a - {\vec b \vec b^T \over \|\vec b\|^2} \vec a = \left (I-{\vec b \vec b^T \over \|\vec b\|^2}\right) \vec a.$$ You’ll often see the denominator written as $\vec b^T\vec b$ instead. This relies on the fact that $\vec a\cdot\vec b = \vec a^T\vec b = \vec b^T\vec a$ and that, since $\vec a\cdot\vec b$ is a scalar, you can move it from right to left in the term $\left(\vec b^T\vec a\right)\vec b$. Note, too that the dot product and it’s expression as a matrix product requires working in an orthononormal basis, but you we’re doing that anyway when you converted cross products into matrix multiplications.