The usual laplacian in two dimensions takes the form,
$$\Delta f = f_{xx} + f_{yy} = f_{rr} + \frac{f_r}{r} + \frac{f_{\theta \theta}}{r^2}$$
when written in polar coordinates. However if I write the hessian of $f$ in polar coordinates, then
\begin{align} \nabla^2 f = \begin{bmatrix} \frac{\partial^2 f}{\partial r^2} & \frac{\partial^2 f}{\partial r\partial \theta} -\frac{1}{r}\frac{\partial f}{\partial \theta} \\ \frac{\partial^2 f}{\partial r\partial \theta} -\frac{1}{r}\frac{\partial f}{\partial \theta} & \frac{\partial^2 f}{\partial \theta^2} + r \frac{\partial f}{\partial r} \end{bmatrix} \end{align}
Taking the trace of this matrix yields $$ f_{rr} + f_{\theta\theta} + rf_r \neq \Delta f \text{ (in radial coordinates)}. $$
Why is this happening and shouldn't the trace remain invariant under a change of coordinates?
NB: We will for now work with holonomic polar coordinates, characterized by the basis $\boldsymbol e_r=\partial_r$, $\boldsymbol e_{\theta}=\partial_\theta$. This is following the formalism of differential geometry. It will not affect the final result.
The Laplacian is defined as the trace of the Hessian, that is $$\Delta :=\operatorname {tr}\mathbf H=H^i{}_i=g^{ij}H_{ij}$$ The Hessian can be viewed as the tensor product of two covariant derivatives: $$\mathbf H=\nabla \nabla \\ H_{ij}=\nabla_i\nabla_j$$
For a scalar $f$, the covariant derivative is the same as the partial derivative, that is $$\nabla_if=\partial_if$$ The object $\nabla f$ is now a covariant vector ($(0,1)$ tensor). For a covariant vector $\boldsymbol \omega$, its covariant derivative has components $$(\nabla\omega)_{ij}=\partial_j\omega_i-\Gamma^{k}_{ij}\omega_k$$ Where $\Gamma$ are our Christoffel symbols. In polar coordinates they are $$\Gamma^\theta_{\theta r}=\Gamma^{\theta}_{r\theta}=\frac{1}{r} \\ \Gamma^r_{\theta\theta}=-r$$ With all others zero. So $$(\nabla\nabla f)_{ij}=(Hf)_{ij}=\partial_i\partial_j f-\Gamma^{k}_{ij}\partial_k f$$
Hence, $$(Hf)_{rr}=\partial_r^2 f \\ (Hf)_{r\theta}=\partial_{r}\partial_{\theta }f-\frac{2}{r}\partial_{\theta }f =(H f)_{\theta r} \\ (H f)_{\theta\theta}=\partial_{\theta}^2 f+r\partial_r f$$ Finally, the inverse metric for polar coordinates has components $$g^{rr}=1~~~,~~~g^{\theta\theta}=\frac{1}{r^2}$$ So $$\Delta f=g^{ij}(Hf)_{ij}=(Hf)_{rr}+\frac{1}{r^2}(Hf)_{\theta\theta} \\ \boxed{\Delta f=\partial_r^2 f+\frac{1}{r}\partial_r f+\frac{1}{r^2}\partial_{\theta}^2 f}$$
As desired!
To answer your question, you are taking the trace of the Hessian by computing $\sum_k H_{kk}$. This is wrong, you need to instead take the sum as $\sum_k H^k{}_k = \sum_{i,j}g^{ij}H_{ij}$. Einstein index notation is often abused - repeated indices should only imply summation if they are paired up and down. Of course in the standard basis the distinction between up and down is irrelevant, but it will lead to wrong answers in curvilinear coordinates.