Can the wavefunction be inferred from the expectation values of operators? Preface
This question is motivated by $C^*$ type treatments of quantum mechanics where operators (Basically an operator is an object that has a spectrum) are treated as fundamental and states are functionals on those operators that map operators to values. The Hilbert space then arises from these fundamental objects. This is in contrast to usual textbook treatments in which the Hilbert space is introduced first and observables are introduced as operators on that Hilbert space.
In the $C^*$ approach, there seems to be a duality between operators and states in this treatment (similar to the duality between Heisenberg and Schrodinger pictures) that I'm interested in understanding better.
The Question:
Suppose I have a quantum system defined by either a Hilbert space or a set of operators. For example my system might be a spin-1/2 system that has $\sigma_{x, y, z}$ as a base set of operators. Or it might be a quantum harmonic oscillator with $p$ and $x$ as base operators. Suppose there is some set $\mathcal{O}$ that contains a base set of operators. Now suppose that for a certain quantum state $|\psi \rangle$, I know the expectation values for all products of elements of $\mathcal{O}$. If I know all of that, then is it possible for me to deduce $|\psi\rangle$?
For example, in the spin-1/2 case, suppose I know $\langle \sigma_x\rangle$, $\langle\sigma_y^2\rangle$, $\langle \sigma_x \sigma_y \sigma_z\rangle$ etc. Is it possible for me to deduce $|\psi\rangle$?
Or in the harmonic oscillator case suppose I know $\langle x\rangle$, $\langle p \rangle$, $\langle x^3 p^4 x^2 p^2\rangle$ etc. Is it possible for me to deduce $|\psi\rangle$?

*

*If the answer is yes then how do I do it computationally?

*Also, have I included too many operators that you need to know the expectation value of? Can you determine $|\psi\rangle$ with a smaller set of operators? What is the smaller set?

*What exactly determines the "base set" of operators that we take products of like $\{x, p\}$ or $\{\sigma_x, \sigma_y, \sigma_z\}$?

My Thoughts
I think the answer to this question is related to a quantum generalization of Bochner's Theorem but I'm not sure. On a similar note, it may be related to this question I asked previously: Quantum Probability, what makes quantum characteristic functions quantum?
 A: Intuitively, if you know the expectation values of all possible observables, that should be enough to fix the state of the system. This almost sounds tautological, since the state is just a mathematical device to compute expectation values of observables. If you already have all these expectation values, then you have already known maximal information about the system allowed by the principles of quantum mechanics. In fact from a more algebraic point of view, a state is defined to be a linear functional on the space of operators.
Here is a "proof" for finite-dimensional case.
Consider a $d$-dimensional Hilbert space $\mathcal{H}$ ($d<\infty$), with an orthonormal basis $\{|n\rangle\}, n=1,2,\dots, d$. A complete basis for linear operators over $\mathcal{H}$ is then $\{O_{mn}=|m\rangle\langle n|\}$, for $m,n=1,2,\dots,d$ (it is an orthonormal basis in the Hilbert-Schmidt sense). If we know the expectation value of all these operators $\langle \psi|O_{mn}|\psi\rangle=\langle n|\psi\rangle\langle \psi|m\rangle$, then we essentially have the matrix element of the density operator $\rho=|\psi\rangle\langle \psi|$: $\langle n|\rho|m\rangle=\langle \psi|O_{mn}|\psi\rangle$.
The assumption of a pure state can be easily dropped, since for any state given by a density operator $\rho$, $\text{Tr} (O_{mn}\rho)=\langle n|\rho|m\rangle$, so knowing the expectation values of all operators is the same as knowing the density operator.
A: Just a generalization to Meng Cheng's answer and theoretically how you can construct $O_{mn}$ given known operators like $x$ and $p$.
Given some Hermitian operator $H$ which we have full knowledge about, that is to say, we know its eigenvectors $\{|n\rangle\}$ and corresponding eigenvalues $E_n$'s. Then, consider the operator $O_{mn}=e^{-\alpha(H-E_m)^2}H_1e^{-\alpha(H-E_n)^2}$ where $H_1$ is some operator awaited to be determined and $\alpha>0$. Now if we have some wavefunction $|\psi\rangle=\sum_{j}{\alpha_j|j\rangle}$, when $\alpha$ is large enough,
$$\langle\psi|O_{mn}|\psi\rangle \sim \alpha^{\dagger}_m\alpha_n\langle m|H_1|n \rangle$$
so $H_1$ is some operator which we choose to make $\langle m|H_1|n \rangle \neq 0$.
Take harmonic oscillator as an example. Let $H$ be the hamiltonian of the harmonic oscillator, so it can be expressed with $x$ and $p$ as you request. If $n$ and $m$ have the same parity (they are both even or both odd), then take $H_1=x^2$, and if they have different parity, take $H_1=x$. The Hilbert space is now extented to infinite dimensions, and the approach gives you a systematic way to approximate $|\psi\rangle$.
A: 
Can the wavefunction be inferred from expectation value of operators?

The answer is yes.
If we consider "the wavefunction" to be the ground state wavefunction then the first Hohenberg-Kohn theorem tells us that the wavefunction is a unique functional of the density.
The density is an expectation value of the density operator.
So, yes is the answer. But, the explicit form of the functional is unknown in all but the simplest cases. Approximating this functional is the preview of "density functional theory." (Although these days it seems like most of the work on "density functional theory" is just a bunch of chemists applying it to systems for which they have no justification to do so.) ;)
A: It actually suffices to know the expectation values of all projection operators of the form $P_\psi:=|\psi\rangle\langle \psi|$ for $\psi \in H$ (which are of course observables). Indeed, suppose we know
$$E_\psi:=\mathrm{Tr}\,\rho\, P_\psi = \langle \psi,\rho\,\psi\rangle $$
for all $\psi \in H$.$^\ddagger$ Then for $\psi_1,\psi_2 \in H$ we find
$$ 2\,\mathrm{Re} \,\langle \psi_1,\rho\,\psi_2\rangle = E_{\psi_1+\psi_2}- E_{\psi_1} - E_{\psi_2}  $$
and
$$ -2\,\mathrm{Im} \,\langle \psi_1,\rho\,\psi_2\rangle = E_{\psi_1+i\psi_2}- E_{\psi_1} - E_{\psi_2} \quad . $$
This eventually allows us to determine $\langle \psi_1,\rho\,\psi_2\rangle$ for all $\psi_1,\psi_2 \in H$ and thus $\rho$ itself. Moreover, it is easy to see that the obtained density operator $\rho$ is unique: For a density operator $\sigma$ with the same expectation values as $\rho$, the above considerations show that for all $\psi_1,\psi_2\in H$ it holds that $\langle \psi_1,\rho\,\psi_2\rangle = \langle \psi_1,\sigma\,\psi_2\rangle$, which in turn implies $||(\rho-\sigma)\psi||=0$ for all $\psi \in H$, i.e. $\rho=\sigma$.
As a last point, note that if $\rho$ happens to be pure, i.e. there exist a unit vector $\phi \in H$ such that $\rho= P_\phi$, then we can determine $\phi$ up to a phase; indeed, it is the only (normalized) eigenvector of $\rho$ with a corresponding non-zero eigenvalue.

$^\ddagger$ In fact/ Alternatively, we only need the expectation values $E_{\psi_k}$, $E_{\psi_k + \psi_j}$ and $E_{\psi_k+i\psi_j}$ for some complete orthonormal basis $\{\psi_k\}_{k\in I}$, where $I$ is an index set with $|I|=\dim H$ and $k,j \in I$. Even more, due to the self-adjointness of $\rho$ we only have to determine $\langle \psi_k,\rho\,\psi_j\rangle$ for $k\leq j$, so we need $E_{\psi_k+i\psi_j}$ only for $k<j$ (and $E_{\psi_k + \psi_j} = E_{\psi_j + \psi_k}$ anyway). Additionally, the trace condition of $\rho$ yields $1=\sum_{k\in I} E_{\psi_k}$ such that in the case of a finite-dimensional Hilbert space $H$ we only need $\dim H-1$ of the $E_{\psi_k}$ and thus in total $d^2-1$ expectation values.
