Finding the stabilizer group given a state Consider general pure state
$|\psi\rangle$ 
in some hilbert space $\mathcal{H}$ (which could be a tensor product of other Hilbert spaces)
I would like to know whether there is a way to  systematically find the group of all operators on $\mathcal{H}$ under which $|\psi\rangle$ is invariant.
A little more precisely: the group I'm looking for is 
$$
G=\{A \in H(\mathcal{H}): A|\psi\rangle = |\psi\rangle\}
$$
Where $ H(\mathcal{H})$ is the set of Hermitian operators acting on the Hilbert space.
This group I defined here is a stabilizer group in the usual sense in mathematics, and a slight generalization of the ones popularized by Gottesman in quantum information as we're considering a general state. 
 A: First of all, the problem is technically difficult due to the fact that generally unbounded self-adjoint operators like those used in general QM have domain smaller that the whole Hilbert space. 
For this reason I will consider here only bounded self-adjoint operators whose domain, as is well known, is the full Hilbert space. 
Proposition. The elements of $G$ are all the operators of this form:
$$A= P^\perp B P^\perp + P \tag{1}$$
where:
(i) $B=B^\dagger$ is every bounded self-adjoint operator
(ii) $P$ is an orthogonal projector such that $P \geq |\psi\rangle \langle \psi|$ (I am assuming that $|\psi\rangle$ has norm $1$) and $P^\perp := I -P$.
(In other words $P$ is the orthogonal projector on a subspace including $|\psi\rangle$ and $P^\perp$ is the orthogonal projector on the orthogonal subspace to that space.)
PROOF. If $A= P^\perp B P^\perp + P$ then $A|\psi\rangle = |\psi\rangle$ by construction (notice that $P^\perp |\psi\rangle=0$) and $A= A^\dagger$ since $B,P,P^\perp$ are self-adjoint so that
$$A^\dagger= (P^\perp)^\dagger B^\dagger (P^\perp)^\dagger + P^\dagger
= P^\perp B P^\perp + P =A\:.$$
Conversely, if $A\in G$, as it is self-adjoint, let us consider its spectral decomposition:
$$A = \int_{\sigma(A)} \lambda dP^{(A)}(\lambda)\:.$$
As $1 \in \sigma_p(A)$, because $|\psi\rangle$ is an eigenvector with eigenvalue $1$, the spectral measure satisfies $P^{(A)}(\{ 1\}) \geq |\psi\rangle \langle \psi|$ and the integral can be decomposed as:
$$A = \int_{\sigma(A)\setminus\{1\}} \lambda dP^{(A)}(\lambda)  + 1P^{(A)}(\{ 1\})\:.$$
where 
$$C := \int_{\sigma(A)\setminus\{1\}} \lambda dP^{(A)}(\lambda)$$
admits $(P^{(A)})^{\perp}(\cal H)$ as invariant space and vanishes on $P^{(A)}(\cal H)$, just for the additive property of the spectral measure $P^{(A)}$. Therefore, it holds, for $P:=P^{(A)}$,
$$P^\perp C P^\perp =C\:.$$
Defining $B:=C$ we have again, 
 $A= P^\perp B P^\perp + P$.
QED
Notice that, for every $A\in G$, there are many couples  $(B,P)$ associated to it via (1). However running $B$ throughout the real space of self-adjoint operators (1) and $P$ throughout the set of orthogonal projectors on subspaces including $|\psi\rangle$  reproduces all elements of $G$.
You also see that $G$ is not a group, with the usual definition in math, with respect to the composition of operators, because if $\cal H$ has sufficiently large dimension you can find two non-commuting elements of $G$ so that their product does not belong to $G$ as it is not self-adjoint. It is not a group eferhring to the sum of operators because $G$ does not include the null operator.
The problem with unbounded operators is that $BP^{\perp}$ may be undefined because the range of $P^\perp$ may not belong to the domain of $B$. However it is possible to re-arrange a similar construction if specifying some detail.
