# How to reconstruct the density matrix $\rho$ from the overlap matrix $T_{a,a'}={\rm Tr}(M^{(a)}M^{(a')})$?

Suppose we have $$N$$-qubit POVM $${\bf M} = \{M^{(a_1)} \otimes M^{(a_2)} \otimes \cdots \otimes M^{(a_N)}\}_{a_1, \ldots, a_N}.$$ Given an $$N$$-qubit state $$\rho$$, the measurement outcome $${\bf a} = (a_1, a_2, \ldots, a_N)$$ occures with probability $$P({\bf a})$$ such that $$\sum_{{\bf a}}P({\bf a})=1$$. Equivalently, we can write $$P({\bf a}) = Tr(M^{({\bf a})} \rho).$$

Provided that the measurement is informationally complete, the density matrix can be unambiguously inferred from the probability distribution of measurement outcomes: $$\tag{1} \rho = \sum_{{\bf a, a'}} P({\bf a}) T^{-1}_{{\bf a, a'}} M^{({\bf a'})},$$ where the matrix $$T$$ is called an overlap matrix defined as follows: $$T_{{\bf a, a'}} := Tr(M^{({\bf a})}M^{({\bf a'})}).$$

I would like to know how (1) can be understood from intuitive and mathematical points of view. My main point of confusion is the matrix $$T^{-1}_{{\bf a,a'}}$$.

Let $$\rho$$ be an arbitrary state, and $$\{\mu_b\}_b$$ some POVM, so that the outcome probabilities are $$p_b(\rho)=\operatorname{tr}(\mu_b \rho)$$. You're asking what are the (Hermitian) operators $$M_b$$ such that you get the decomposition $$\rho=\sum_b p_b(\rho) M_b$$ for all $$\rho$$.
This is completely analogous to the following question: given a finite-dimensional vector space $$V$$ and a finite subset $$\{v_k\}\subset V$$ that spans $$V$$, is there some set $$\{w_k\}$$ such that for all $$v\in V$$ we have $$v=\sum_k \langle v_k,v\rangle w_k$$? One way to answer this question is to observe this relation is equivalent to $$v=WV^\dagger v$$, with $$V$$ and $$W$$ matrices whose columns are the vectors $$\{v_k\}$$ and $$\{w_k\}$$, respectively. If we want this for all $$v\in V$$, that means we're asking for vectors $$w_k$$ such that $$WV^\dagger =I$$. This is an inhomogeneous linear system, a solution of which can always be written $$W=(V^\dagger)^+\equiv (V V^\dagger)^{-1} V,$$ assuming $$VV^\dagger$$ is invertible, that is, $$V$$ is surjective (which is true iff $$\{v_k\}$$ span $$V$$). You might now observe that $$VV^\dagger=\sum_k v_kv_k^\dagger$$. Defining $$S\equiv VV^\dagger$$, we can restate this result as saying that a viable set of "dual" vectors $$w_k$$ is given by $$w_k = S^{-1}v_k$$, and therefore $$v = \sum_k \langle v_k,v\rangle S^{-1}(v_k). %= \sum_k \langle S^{-1}(v_k),v\rangle v_k.$$
All of this can be applied to decompose states in terms of POVM elements, just replacing $$v\to \rho$$ and $$v_k\to \mu_b$$, and using the trace inner product: $$\langle\mu_b,\rho\rangle\equiv \operatorname{tr}(\mu_b\rho)$$. You'll then end up with the decomposition $$\rho = \sum_b \langle\mu_b,\rho\rangle S^{-1}(\mu_b), \qquad S \equiv \sum_b |\mu_b\rangle\!\rangle\langle\!\langle\mu_b|,$$ where $$S$$ is the superoperator defined as $$S(X)\equiv \sum_b \mu_b \langle\mu_b ,X\rangle$$ for all $$X$$. I think this is essentially the expression you mentioned, though your $$T$$ should not be the "overlap matrix" but rather the $$S$$ operator above. Unless you're making some additional assumption (the expressions might be identical for minimal IC-POVMs, though I haven't really checked).