Both the Heisenberg picture and the Schrödinger picture
are rigorously equivalent when they are defined rigorously.
First, we define time evolution as a projective representation ##\mathbb{R}\ni\tau\longmapsto\gamma_{\tau}##, where the maps ##\rho\longmapsto\gamma_{\tau}(\rho)## are Kadison symmetries (##\rho## is a positive trace class operator of trace 1; Kadison symmetries preserve the convexity of the space of states; when reduced to pure states, it's simply a Wigner symmetry). By Kadison's Theorem and using the fact that the group here is ##\mathbb{R}##, we get that there's a unitary one parameter group ##\mathbb{R}\ni\tau\longmapsto U_{\tau}## such that ##\rho_{_{\tau}}=\gamma_{\tau}(\rho)=U_{\tau}\rho U_{\tau}^{-1}##. This gives the Schrödinger evolution for states (all states, not only pure).
Now, taking the symmetries ##\gamma_{\tau}##, characterized by the unitary one parameter group ##U_{\tau}##, we can define the
dual action, ##\gamma_{\tau}^{*}##, of the symmetry, which acts on the lattice of orthogonal projectors and preserves its structure: ##\gamma_{\tau}^{*}(P)\doteq U_{\tau}^{-1}PU_{\tau}## (all operators here are bounded, so there are no domain issues). An observable ##A## is represented by a projector valued measure ##P^{(A)}## over the real line. So, we define the dual action on observables by: ##P^{(A_{H}(\tau))}(E)\doteq\gamma_{\tau}^{*}(P^{(A)}(E))##, for any ##E## in the Borel algebra of the real line. We call ##A_{H}(\tau)## the observable ##A## in the Heisenberg picture, characterizaed by the PVM ##P^{(A_{H}(\tau))}##.
It's easy to show that: ##tr(\rho\gamma^{*}(P))=tr(\gamma(\rho)P)##. But, since ##tr(\rho P^{(A)}(E))## represents the probability that the measuring of ##A## in state ##\rho## falls in ##E##, the mentioned identity shows the equivalence of both pictures in a rigorous way.
And that's it.
By the spectral theorem, the dual action is obviously extended to self-adjoint operators and we get: ##A_{H}(\tau)=\gamma_{\tau}^{*}(A)= U_{\tau}^{-1}AU_{\tau}##. Also: ##\sigma(A_{H}(\tau))=\sigma(A)## and, defining the expectation value as ##<A>_{\rho}=\int_{\mathbb{R}}\lambda\mathrm{d}\mu_{\rho}^{(A)}(\lambda)## (where ##\mu_{\rho}^{(A)}## is the probability measure ##\mu_{\rho}^{(A)}(E)=tr(\rho P^{(A)}(E))##), we also get ##<A>_{\rho_{_{t}}}=<A_{H}(t)>_{\rho}##, which are the usual standard results. (notice that we only need the trace identity mentioned above for proving this result, there's no need to invoke the operators ##A##, that's the advantage of introducing this probability measure).
We could start to introduce the (most likely, unbounded) self-adjoint operators associated with the PVM ##P^{(A)}## and also to write the time evolution in infinitesimal or differential form, but in doing that you introduce a lot of spurious technical problems related to the domain of the operators involved.
The above definitions and results also hold when there's an explicit time dependence. In that case, the definition is: ##P^{(A_{Ht}(\tau))}(E)\doteq\gamma_{\tau}^{*}(P^{(A_{t})}(E))##. This can also be extended to non-homogeneous time evolution (i.e., for "time dependent hamiltonians"), though I guess some of the details may need some work, since the time evolution is not given by the one used here.
