First, I think your $H_{SOC}$ is incorrect and should read:
$$
H_{SOC}=\frac{\alpha}{2}\left(L_+\sigma_- + L_-\sigma_+ + L_z\sigma_z\right)\, .
\tag{1}
$$
The twist with your problem is that it uses the real spherical harmonics, rather than the complex exponential form more common in physics and better adapted to the evaluation of matrix elements. The conversion between the two can be found on this wiki page.
Thus:
\begin{align}
p_y&=\frac{i}{\sqrt{2}}\left(Y_{1}^{-1}+Y_1^1\right)\quad\Rightarrow \quad \vert p_y\rangle = \frac{i}{\sqrt{2}}\left(\vert 1,-1\rangle + \vert 1,1\rangle\right)\, ,\tag{2}\\
p_x&=\frac{1}{\sqrt{2}}\left(Y_{1}^{-1}-Y_1^1\right)
\quad\Rightarrow \quad \vert p_y\rangle = \frac{i}{\sqrt{2}}\left(\vert 1,-1\rangle + \vert 1,1\rangle\right)\, ,\tag{3}\\
p_z&=Y_1^0\qquad\qquad \qquad\quad\Rightarrow \quad \vert p_z\rangle = \vert 1,0\rangle\, , \tag{4}
\end{align}
where the arguments $(\theta,\varphi)$ of the $Y_{\ell}^{m}$ have been omitted for clarity.
The next step is to recall the action of the various operators on $Y_{\ell m}$:
\begin{align}
\hat L_\pm Y_{\ell}^{m}&=\sqrt{(\ell\mp m)(\ell\pm m+1)}Y_{\ell}^{m\pm 1}\, ,\tag{5}\\
\hat L_z Y_{\ell}^{m}&= mY_{\ell}^{m}\tag{6}
\end{align}
from which one obtains for instance (if I did not make obvious mistakes), the matrix representation of $\hat L_+$ in the basis $\{\vert p_x\rangle, \vert p_y\rangle, \vert p_z\rangle\}$
\begin{align}
\hat L_+\vert p_x\rangle &= \vert p_z\rangle \qquad
\hat L_+\vert p_y\rangle = i \vert p_z\rangle \qquad
\hat L_+\vert p_z\rangle =-\vert p_x\rangle - i\vert p_y\rangle\, , \tag{7}
\end{align}
so that
\begin{align}
\hat L_+\to \left(\begin{array}{ccc}
0&0&-1\\
0&0&-i\\
1&i&0\end{array}\right)\, ,\qquad \hat L_-\to
\left(\begin{array}{ccc}
0&0&1\\
0&0&-i\\
-1&i&0\end{array}\right)\, . \tag{8}
\end{align}
since the matrix representation for $\hat L_-$ will be transpose conjugate of the matrix representation for $\hat L_+$. Likewise:
\begin{align}
\hat \sigma_+\vert \uparrow \rangle &=0\, , \quad
\hat \sigma_+\vert \downarrow \rangle =\vert \uparrow\rangle\, ,\quad
\hat \sigma_-\vert \uparrow \rangle =\vert\downarrow\rangle\, , \quad
\hat \sigma_-\vert \downarrow \rangle =0\, . \tag{9}
\end{align}
Hence:
\begin{align}
\hat L_+\hat \sigma_- \vert p_x;\uparrow\rangle &= \Bigl[\hat L_+\vert p_x\rangle \Bigr]
\Bigl[\hat \sigma_-\vert \uparrow\rangle\Bigr]= \vert p_z\rangle \vert\downarrow\rangle= \vert p_z;\downarrow\rangle\, , \tag{10}\\
\hat L_-\hat \sigma_+ \vert p_x;\uparrow\rangle &= 0 \, ,\tag{11} \\
\hat L_z\hat\sigma_z\vert p_x;\uparrow\rangle &=
\Bigl[\hat L_z\vert p_x\rangle \Bigr]
\Bigl[\hat \sigma_z\vert \uparrow\rangle\Bigr]=i\vert p_y\rangle \vert \uparrow\rangle=
i\vert p_y;\uparrow\rangle \tag{12}
\end{align}
This is not quite the first column of your matrix, but as I understand it this matrix is not quite correct either since it ought to be hermitian but is not: for instance the matrix element in position $(2,6)$ should be the complex conjugate of the matrix element in position $(6,2)$.
Anyways, my calculation would give the elements in the first column as
$$
(0,0,i,0,0,1)^T \tag{13}
$$
This might be because the basis elements are ordered as
$\{\vert p_x\uparrow\rangle,\vert p_y\uparrow\rangle,\vert p_z\uparrow\rangle,\vert p_x\downarrow\rangle,\vert p_y\downarrow\rangle,\vert p_z\downarrow\rangle\}$ rather than the ordering you give (or else it's quite possible I made an error somewhere). With this ordering my first column would come out as $(0,i,0,0,0,1)^T$, which is identical to yours. Note that this does not fix the hermiticity problem mentioned earlier.
The other columns are found in the same manner.
Just a formal wisecrack: you know how to Lagrange-shift functions,
$$
e^{a\partial_x + b\partial_y} f(x,y)=f(x+a,y+b),
$$
a mere rewriting of the Taylor expansion.
You may further perform two evaluations in an important order,
$$
f(x+a,y+b) \Big|^{x=0}_{y=0}= f(a,b),
$$
and
$$
f(a,b) \Big|^{a=y}_{b=x}= f(y,x),
$$
so that
$$
P^r_{12} f(x,y)= \left( \left((e^{a\partial_x + b\partial_y} f(x,y) ) \Big|^{x=0}_{y=0} \right ) \Bigg|^{a=y}_{b=x} \right) =f(y,x)
$$
so, hyper-formally,
$$
P^r_{12}=\Big|^{a=y}_{b=x} \circ \Big|^{x=0}_{y=0} \circ e^{a\partial_x + b\partial_y} .
$$
This is an infinite-dimensional representation, so its exchange operator is predictably messier than for the other two finite-dimensional representations you exchanged. Not clear to me what you might do with it...
Best Answer
it's a bit difficult to answer, as "coming up" with an expression is very personal and depends on intuition etc. It is especially so if one knows the right answer already, and then "how to come up" with it is plagued by hindsight. Having said that, you can look at the spin exchange as wanting to do the following things:
i) if one of the spins is up and the other is down, take the down up and the up down. So this is achieved by a term like $\sigma_1^{+}\sigma_2^{-} + \sigma_1^{-}\sigma_2^{+}$. Note that acting with this on a state where both spins are identical will lead to zero, which is good for us -- it won't get us into trouble in the case that they are equal.
ii) if both spins are identical, just leave it be, which we can write as $1$. But we want it to operate only if they are equal. So let's subtract $1$ in case they are not equal, which means that $\sigma_1^z\sigma_2^z=-1$. Then we write this part as $(1+\sigma_1^z\sigma_2^z)/2$ which does the trick.
Now we just have to add everything up, and as $\sigma^{\pm} = (\sigma^x \pm i\sigma^y)/2$, you get the expression you wanted.