Eligibility vector for softmax policy with policy gradients
Calculation of gradient
\begin{align}
\nabla_{\theta} \log(\pi_{\theta}(a|s)) &amp;= \phi(s,a) - \mathbb E[\phi (s, \cdot)]\\
&amp;= \phi(s,a) - \sum_{a'} \pi(a'|s) \phi(s,a')
\end{align}
is only valid for linear function approximation with action preferences of form
\begin{equation}
h(s, a, \theta) = \theta^T \phi(s,a)
\end{equation}
and softmax policy
\begin{equation}
\pi(a|s) = \frac{e^{h(s,a,\theta)}}{\sum_{a'} e^{h(s,a',\theta)}}
\end{equation}
The gradient would be calculated as it is written. For example, if your current state is $s = (1, 1)$ and in that state you have actions $a_0 = 0$ and $a_1 = 1$ and probabilities for those actions are $\pi(a_0|s) = 0.7$, $\pi(a_1|s) = 0.3$ then gradient for action $a_0$ would be
\begin{equation}
\nabla_{\theta} \log(\pi_{\theta}(a_0|s)) = (1, 1, 0)^T - (0.7 \cdot (1, 1, 0)^T + 0.3\cdot (1,1,1)^T) = (0, 0, -0.3)^T
\end{equation}
Feature vector $\phi$ can be basically anything you want. For example you could stack state feature and action (like I did in small example), you could use polynomials, radial basis functions, tile coding, etc.  

If you're using multilayered neural network you would have to propagate gradients through all layers, usually done with backpropagation algorithm. Easiest way is to use automatic differentiation software (e.g. Tensorflow) which can do that for you so you don't have to write your implementation. All you have to do is define your objective function that you want to optimize
\begin{equation}
J_\theta = \sum_t \log(\pi(a_t|s_t, \theta)) G_t
\end{equation}
and software will automatically calculate gradient $\nabla J_{\theta}$ and update weights.