How to implement a 2x2 operator in Mathematica?

I am trying to write a Mathematica program to compute the following:

For a given Hermitian matrix $\rho$, the operator $L_\theta$ with respect to a parameter $\theta$ is defined as:

\begin{equation} L_\theta = -\vec{y} \cdot \vec{r} \, \mathbb{I} + \vec{y} \cdot \vec{\sigma}, \quad \vec{y} = (\mathbb{I} - \vec{r}\vec{r}^T)^{-1} \, \partial_\theta \vec{r} \end{equation}

Here, $\vec{r}$ is a vector defined as:

\begin{equation} \vec{r} = (\mathrm{Tr}[\rho \sigma_x], \mathrm{Tr}[\rho \sigma_y], \mathrm{Tr}[\rho \sigma_z]), \end{equation}

$\vec{\sigma} = (\sigma_x, \sigma_y, \sigma_z)$ are the Pauli matrices, and $\partial_\theta \vec{r}$ is the derivative with respect to $\theta$.

I would like to write a Mathematica program that, given matrix $\rho[\theta]$, compute vector $\vec{r}$, its derivative with respect to $\theta$, and then constructs $L_\theta$ automatically.

I tried something like this:

 sigmaX = {{0, 1}, {1, 0}}; sigmaY = {{0, -I}, {I, 0}}; sigmaZ = {{1, 0}, {0, -1}}; Vector[rho_] := {Tr[rho.sigmaX], Tr[rho.sigmaY], Tr[rho.sigmaZ]}; rho[theta_] := {{(1 + Cos[theta])/2, Sin[theta]/2}, {Sin[theta]/2, (1 - Cos[theta])/2}}; Operator[rhoFunc_, theta_] := Module[{rho, r, dr, Id, P, L}, rho = rhoFunc[theta]; r = Vector[rho]; dr = D[r, theta]; Id = IdentityMatrix[3]; P = Inverse[Id - Outer[Times, r, r]].dr; L = -P.r IdentityMatrix[2] + P[[1]] sigmaX + P[[2]] sigmaY + P[[3]] sigmaZ; L ]; 

But I am not sure if this is the best way to implement it, and I get errors sometimes with defining $\rho[\theta]$ or computing derivatives.

I would appreciate a solution that works both symbolically and numerically.