Background
I am a PhD student in Physics and I am currently developing quite refined computer codes that allow to simulate many-body quantum systems living on a lattice. The difficulty resides in the fact that the corresponding Hilbert space, being the tensor product of many local spaces, has a huge dimension. To those who are interested in details I can say that I use the so-called Density Matrix Renormalization Group (DMRG) algorithm in its more evolved formulation using Matrix Product States (MPSs). Quite some time ago (five months at least) I stumble upon a problem which is very practical from my point of view (I cannot proceed any farther with my code if I don't solve it), but at the same time of a conceptual and strictly mathematical nature. I believe this to be a very interesting issue for its implication on my field of research but at the same time quite an interesting and challenging mathematical problem in (multi-)linear algebra. In a special case it states the existence of the Schmidt decomposition of a vector living on a tensor product space, a well-known mathematical tool in the field of quantum information theory, the latter being the theoretical framework behind MPSs. Despite my considerable effort and the help of another PhD student in numerical analysis and linear algebra the problem is still open.
The Problem
Suppose we are given a set of $n^2$ linearly independent vectors $\lbrace v_{ij}\rbrace_{i,j\in[1,\dots,n]}$ in a complex vector space $V$ endowed with a positive definite hermitian scalar product. All possible scalar products between vectors in the above set can be collected in the Gramian matrix
$[G]_{kl,ij} = \langle v_{kl}\,;v_{ij}\rangle\quad$.
$[A]_{i,j}$ is the entry in the $i$-th row and $j$-th column of matrix $A$. Rows and columns can be indexed by pairs of indices as in the Gramian matrix $G$ above. The goal is to find a $n\times n$ matrix $X$ generating a set of transformed vectors $\lbrace\tilde{v}_{kl}\rbrace_{k,l\in[1,\dots,n]}$ as follows
$\tilde{v}_{kl} = > \sum_{i,j}X_{ki}X_{lj}^{-\mathrm{T}}v_{ij}\quad$,
( $X^{-\mathrm{T}} = (X^{-1})^\mathrm{T}$, the inverse transpose of $X$ ) such that the $\tilde{v}_{kl}$ are orthogonal for $k=l$ and arbitrary otherwise, namely
$\langle \tilde{v}_{aa};\tilde{v}_{bb} \rangle = \sigma_{a}\delta_{ab}\quad$.
Here $\delta_{ab}$ is the Kronecker delta and $\sigma_a$ are positive real numbers for $a = 1,\dots,n$.
Since the Gramian matrix is indexed by pairs of integers it is endowed with a tensor product structure and we can restate the problem as follows: find a $n \times n$ matrix $X$ such that
$\mathcal{P}(X^\mathrm{T}\otimes > X^{-1})^{\dagger}G(X^\mathrm{T}\otimes > X^{-1})\mathcal{P} = \sigma\qquad > \text{with} \qquad [\sigma]_{a,b} = > \sigma_{a}\delta_{ab}\,$.
Here $\mathcal{P}$ is just a projector on rows and columns indexed by pairs of the form $ii$. So, for instance, $\mathcal{P}G\mathcal{P}$ is a $n\times n$ matrix with entries $[ \mathcal{P}G\mathcal{P}]_{i,j}= [G]_{ii,jj}$. Also the dagger $\dagger$ superscripts stands for the transponse conjugate of a matrix.
Comments
A few comments are in order. The definition of the transformed vectors $\tilde{v}_{kl}$ implies that $w = \sum_k \tilde{v}_{kk} = \sum_i v_{ii}$ is left invariant by the transformation and in fact this is the reason behind the choice of such a definition.
The problem admits a solution in a special case. Let $V = V^{(L)} \otimes V^{(R)}$ and $v_{ij} = v_{i}^{(L)}\otimes v_{j}^{(R)}$ with $\mathcal{B}_L = \lbrace v_i^{(L)}\rbrace_{i\in[1,\dots,n]}$ and $\mathcal{B}_R = \lbrace v_j^{(R)}\rbrace_{j\in[1,\dots,n]}$ linearly independent sets of vectors in $V^{(L)}$ and $V^{(R)}$ respectively. The scalar product on $V$ is defined in the standard way from scalar products on $V^{(L)}$ and $V^{(R)}$, so
$\langle v_{kl};v_{ij}\rangle = \langle v_k^{(L)}\otimes v_l^{(R)};v_i^{(L)}\otimes v_j^{(R)}\rangle = \langle v_k^{(L)};v_i^{(L)}\rangle_L\langle v_l^{(R)};v_j^{(R)}\rangle_R\quad$.
Then it is not difficult to find a matrix $X$ such that $\tilde{v}^{(L)}_k = \sum_{i}X_{ki}v_{i}^{(L)}$ and $\tilde{v}^{(R)}_l = \sum_{l}X_{lj}^{-\mathrm{T}}v_{j}^{(R)}$ are orthogonal in their respective spaces (for all the details you can have a look at an extended summary of the problem here ). Therefore
$w = \sum_i v_i^{(L)}\otimes v_i^{(R)} = \sum_k \tilde{v}_{k}^{(L)} \otimes \tilde{v}_{k}^{(R)}\quad$
and on the rightmost side one can recognize the Schmidt decomposition of $w$.
I am interested in knowing if such a matrix $X$ exists and under which conditions. Moreover if it exists, I am interested in a practical algorithm to find it since the matrix $X$ is to be calculated several times in a numerical code.
