This can be proven using index notation (and the Einstein summation convention) pretty easily. Note that in this answer, I put all indices downstairs, and whenever any two are repeated, we sum over them.
\begin{align}
\left(\int_{\partial S} \vec{r} \times d \vec{s} \right)_i &= \int_{\partial S} \epsilon_{ijk} r_j\, ds_k
\end{align}
Where $\epsilon_{ijk}$ is the Levi-Civita symbol (defined such that $\epsilon_{123} = 1$ and such that it is totally alternating in the indices)
Note that the RHS of this expression is can be written as $\vec{F} \cdot d \vec{s}$ if you define $F_k = \epsilon_{ijk} r_j$. So,
\begin{align}
\left(\int_{\partial S} \vec{r} \times d \vec{s} \right)_i &= \int_{\partial S} \vec{F} \cdot d \vec{s} \\
&= \int_S (\nabla \times F) \cdot d\vec{a} \tag{Stokes Theorem}
\end{align}
Now, you can continue carrying out this computation by diligently writing out what the curl of $\vec{F}$ is, taking the dot product with $d\vec{a} = \hat{n} \, da$ ($\hat{n}$ being the unit outward normal vector field to the oriented surface $S$). Or, you can continue with the index manipulations as I show here:
\begin{align}
\int_S (\nabla \times F) \cdot d \vec{a} &= \int_S (\nabla \times F)_{\mu} \, da_{\mu} \\
&= \int_S \epsilon_{\mu \nu k} (\partial_{\nu} F_k) \, da_{\mu} \\
&= \int_S \epsilon_{\mu \nu k} \left(\partial_{\nu} (\epsilon_{ijk} r_j ) \right) \, da_{\mu} \\
&= \int_S \epsilon_{\mu \nu k} \epsilon_{ijk} (\partial_{\nu} r_j) \, da_{\mu} \\
&= \int_S \epsilon_{\mu \nu k} \epsilon_{ijk} \delta_{\nu j} \, da_{\mu} \\
&= \int_S \epsilon_{\mu j k} \epsilon_{ijk} \, da_{\mu} \\
&= \int_S (\delta_{\mu i} \delta_{jj} - \delta_{\mu j} \delta_{j i}) \, da_{\mu} \tag{$*$}
\end{align}
Now, unfortunately $(*)$ is just one of those identities which is completely non obvious, and you just have to perform a massive direct computation to verify it (there are several proofs of this fact online).
Note that $\delta_{jj} = \sum_{j=1}^3 \delta_{jj} = \sum_{j=1}^3 1 = 3$ is the dimension of the space we are working in (i.e the trace of the $3 \times 3$ identity matrix)
So, continuing, we have
\begin{align}
\int_S (\delta_{\mu i} \delta_{jj} - \delta_{\mu j} \delta_{j i}) \, da_{\mu} &= \int_S \left( \delta_{\mu i} \cdot 3 - \delta_{\mu i}\right) \, da_{\mu} \\
&= \int_{S} 2 \delta_{\mu i} \, da_{\mu} \\
&= 2 \int_S da_{i} \\
&= (2 \vec{A})_i
\end{align}
In other words we showed that for all $i \in \{1,2,3\}$,
\begin{align}
\left(\int_{\partial S} \vec{r} \times d \vec{s} \right)_i &= \left(2 \int _S d \vec{a}\right)_i \equiv (2 \vec{A})_i
\end{align}
Since all the components of the two vectors (with respect to the standard basis of $\Bbb{R}^3$) are equal, the vectors themselves agree. This proves your formula.
By the way, there might be some variant of Stokes theorem written in vector notation which may be of direct use in this case, but honestly, I can never remember them so I rather just quickly re-derive things myself. I suggest you do not spend too much time trying to remember all the variants either, because they all require you to think of a "clever" function to apply the standard Stokes formula to.
As much as I hate working with indices and the Einstein summation convention, I think this really is a pretty quick computation once you're familiar with the notation.
