A 3D rotation around an arbitrary point $(x_0, y_0, z_0)$ is described by
$$\left[ \begin{matrix} x^\prime \\ y^\prime \\ z^\prime \end{matrix} \right] =
\left[ \begin{matrix} X_x & Y_x & Z_x \\ X_y & Y_y & Z_y \\ X_z & Y_z & Z_z \end{matrix} \right] \left[ \begin{matrix} x - x_0 \\ y - y_0 \\ z - z_0 \end{matrix} \right] + \left[ \begin{matrix} x_0 \\ y_0 \\ z_0 \end{matrix} \right]$$
which, as OP noted, first subtracts the center point, rotates around the origin, then adds back the center point; equivalently written as
$$\vec{p}^\prime = \mathbf{R} \left( \vec{p} - \vec{p}_0 \right) + \vec{p}_0$$
We can combine the two translations, saving three subtractions per point – not much, but might help in a computer program. This is because
$$\vec{p}^\prime = \mathbf{R} \vec{p} + \left( \vec{p}_0 - \mathbf{R} \vec{p}_0 \right)$$
In other words, you can use the simple form, either
$$\left[ \begin{matrix} x^\prime \\ y^\prime \\ z^\prime \end{matrix} \right] =
\left[ \begin{matrix} X_x & Y_x & Z_x \\ X_y & Y_y & Z_y \\ X_z & Y_z & Z_z \end{matrix} \right] \left[ \begin{matrix} x \\ y \\ z \end{matrix} \right] + \left[ \begin{matrix} T_x \\ T_y \\ T_z \end{matrix} \right]$$
or, equivalently,
$$\left[ \begin{matrix} x^\prime \\ y^\prime \\ z^\prime \\ 1 \end{matrix} \right] = \left[ \begin{matrix} X_x & Y_x & Z_x & T_x \\ X_y & Y_y & Z_y & T_y \\ X_z & Y_z & Z_z & T_z \\ 0 & 0 & 0 & 1 \end{matrix} \right] \left[ \begin{matrix} x \\ y \\ z \\ 1 \end{matrix} \right]$$
where
$$\left[ \begin{matrix} T_x \\ T_y \\ T_z \end{matrix} \right] =
\left[ \begin{matrix} x_0 \\ y_0 \\ z_0 \end{matrix} \right] -
\left[ \begin{matrix} X_x & Y_x & Z_x \\ X_y & Y_y & Z_y \\ X_z & Y_z & Z_z \end{matrix} \right] \left[ \begin{matrix} x_0 \\ y_0 \\ z_0 \end{matrix} \right]$$
to apply a rotation $\mathbf{R}$ around a centerpoint $(x_0, y_0, z_0)$.
The 4×4 matrix form is particularly useful if you use SIMD, like SSE or AVX, with four-component vectors; that's one reason why many 3D libraries use it. Another is that the same form can be used for projection.
It will turn out that the rotation in polar coordinates is nothing else than a shift
of coordinate functions by an angle $\alpha$. Note that we usually shift a function
$g$ on $\mathbb R$ by $\alpha$ by going to the function $x\mapsto g(x\color{red}{-}\alpha)$.
A very convenient method to formalize vector fields in different coordinates is to treat them as derivations. This is very common in differential geometry. In $2d$ Cartesian coordinates a vector field is then
$$
X=X_1\frac{\partial}{\partial x_1}+X_2\frac{\partial}{\partial x_2}
$$
where $X_1,X_2$ are scalar functions on $\mathbb R^2$.
From the Cartesian-polar transformation
\begin{align}
r&=\sqrt{x_1^2+x_2^2}\,,&
\theta&=\arctan\frac{x_2}{x_1}\,\\
x_1&=r\cos\theta\,,&x_2&=r\sin\theta
\end{align}
we get
\begin{align}
\frac{\partial x_1}{\partial r}&=\cos\theta\,,&\frac{\partial x_1}{\partial \theta}=-r\sin\theta\,,\\
\frac{\partial x_2}{\partial r}&=\sin\theta\,,&\frac{\partial x_2}{\partial \theta}=r\cos\theta\,.
\end{align}
By the chain rule we have for any differentiable function $f$ on $\mathbb R^2$
\begin{align}
\frac{\partial f}{\partial r}&=\frac{\partial f}{\partial x_1}\frac{\partial x_1}{\partial r}+\frac{\partial f}{\partial x_2}\frac{\partial x_2}{\partial r}
=\frac{\partial f}{\partial x_1}\cos\theta+\frac{\partial f}{\partial x_2}\sin\theta\,,\\[3mm]
\frac{\partial f}{\partial \theta}&=\frac{\partial f}{\partial x_1}\frac{\partial x_1}{\partial \theta}+\frac{\partial f}{\partial x_2}\frac{\partial x_2}{\partial \theta}
=-\frac{\partial f}{\partial x_1}r\sin\theta+r\frac{\partial f}{\partial x_2}\cos\theta\,.
\end{align}
In matrix notation and dropping the test function this is the transformation rule for the basis vector fields:
$$
\begin{pmatrix}\frac{\partial}{\partial r}\\\frac{\partial}{\partial \theta} \end{pmatrix}=\begin{pmatrix}\cos\theta &\sin\theta\\-r\sin\theta&r\cos\theta\end{pmatrix}\begin{pmatrix}\frac{\partial}{\partial x_1}\\\frac{\partial}{\partial x_2} \end{pmatrix}.
$$
Clearly,
$$
\begin{pmatrix}\frac{\partial}{\partial x_1}\\\frac{\partial}{\partial x_2} \end{pmatrix}=\begin{pmatrix}\cos\theta &-\frac{1}{r}\sin\theta\\\sin\theta&\frac{1}{r}\cos\theta\end{pmatrix}\begin{pmatrix}\frac{\partial}{\partial r}\\\frac{\partial}{\partial \theta} \end{pmatrix}.
$$
The vector field $X$ in polar coordinates is therefore,
\begin{align}
X'&=X_1'\Big(\cos\theta\frac{\partial}{\partial r}-\frac{1}{r}\sin\theta\frac{\partial}{\partial \theta}\Big)+X_2'\Big(\sin\theta\frac{\partial}{\partial r}+\frac{1}{r}\cos\theta\frac{\partial}{\partial \theta}\Big)\\[3mm]
&=\Big(X_1'\cos\theta+X_2'\sin\theta\Big)\frac{\partial}{\partial r}+
\Big(-X_1'\frac{1}{r}\sin\theta+X_2'\frac{1}{r}\cos\theta\Big)\frac{\partial}{\partial \theta}
\end{align}
where $X'_1(r,\theta)=X_1(x_2,x_2)$ and $X'_2(r,\theta)=X_2(x_2,x_2)\,.$
The whole point is now that nothing is simpler than rotating the vector field $X'$ by an angle $\alpha$ around the origin.
This rotation is a shift which gives
\begin{align}
X''
&=\Big(X_1''\cos(\theta-\alpha)+X_2'\sin(\theta-\alpha)\Big)\frac{\partial}{\partial r}+
\Big(-X_1''\frac{1}{r}\sin(\theta-\alpha)+X_2''\frac{1}{r}\cos(\theta-\alpha)\Big)\frac{\partial}{\partial \theta}
\end{align}
where $X''_1(r,\theta)=X'_1(r,\theta-\alpha)$ and $X''_2(r,\theta)=X_2'(r,\theta-\alpha)\,.$
Best Answer
Active rotation of the vector $x = \vec{OP}$ to a new vector $x'$: $$ x' = A x $$ where $A$ is the matrix of the rotation.
We want to get the same new coordinates $x_i'$ by a change of point of view, the vector stays the same, by a linear transformation of the coordinate system, where the base vectors $e_i$ get transformed into new base vectors $e_i'$ according to $$ e_i' = B e_i $$ for $i \in \{ 1, \ldots, n \}$.
\begin{align} \underbrace{A^{-1} x'}_{(4)} = x &= \underbrace{\sum_{i=1}^n x_i e_i}_{(1)} \\ &=^! \underbrace{\sum_{i=1}^n x_i' e_i'}_{(2)} = \sum_{i=1}^n x_i' B e_i = \sum_{i=1}^n B x_i' e_i = B \sum_{i=1}^n x_i' e_i = \underbrace{B x'}_{(3)} \end{align} which leads to $B = A^{-1}$.
$(1)$ vector $x$ as linear combination of original base vectors, having coordinates $x_i$
$(2)$ the same vector $x$ expressed as linear combination of new base vectors, having new coordinates $x_i'$
$(3)$ the effect of $B$ (interpreted as active transformation within the old coordinate system) is $B x' = x$
$(4)$ this should be consistent with the active transformation $A x = x'$