Given a matrix of connection weights $W$ (and bias vectors $\vec{a}$ and $\vec{b}$), a restricted Boltzmann machine generates a probability distribution over a discrete set of inputs $\vec{v} = \{v_1, \dots, v_N\}$ by $$ p(\vec{v}) = \sum_{h_1, \dots h_M} \exp(\sum_{j}^{N}a_j v_j + \sum_{i}^{M}b_i h_i + \sum_{i j}W_{i j}h_i v_j) $$
Here, $v_i$ and $h_i$ can take the values $+1$ or $-1$.
For my application, $p(\vec{v})$ can be complex, so the RBM is defined with complex network weights $W_{i j}$ (idem. bias vectors). Of course, $p(\vec{v})$ can no longer be interpreted as probability distribution.
I would like to know if, given a complex-valued RBM, I can write the distribution over the squared moduli $|p(\vec{v})|^2$ as an RBM with real weights. I know it should be possible, given that an RBM can approximate any probability distribution, but can it be done efficiently?
The new hidden units of this RBM do not necessarily have to take the values +1 or -1, or even be binary-valued.
What I tried
The complex-valued RBM has the alternate form $$ p(\vec{v}) = \exp(\sum_{j}^{N}a_j v_j)\prod_{i}^{M}\sum_{h_i}\exp(h_i(b_i + \sum_j W_{i j} v_j)) $$
The objective is to find a similar form for $p(\vec{v})p(\vec{v})^{*}$, with real weight matrix $V$, and real bias vectors $A$ and $B$. It is quickly seen that $$ A_i = a_j + a_j^{*} $$
If I write out the RBM-form of $|p(\vec{v})|^2$, forgetting about the bias vector $A$, I get $$ |p(\vec{v})|^2 = \prod_{i}^{M}\sum_{h_i, h_i'}\exp(h_i(b_i + \alpha_i) + h_i'(b_i^{*} + \alpha_i^{*})) $$ where $$ \alpha_i = \sum_j^N W_{i j} v_j $$ I am looking for an efficient representation $$ |p(\vec{v})|^2 = \prod_{i}^{M'}\sum_{H_i}\exp(H_i(B_i + \sum_j^N V_{i j} v_j)) $$ with $B_i$, $V_{i j}$ real, where $H_i$ does not have to be binary-valued, and the number of new hidden units $M'$ can be different from $M$.
How can I start to solve this? Help is very much appreciated.