Schmidt decomposition#

Let us assume a quantum state with two subsystems, labeled by physical indices $i$ and $j$, running over dimensions $d_1$ and $d_2$. In quantum notation

\[|\psi\rangle = \sum_{i=1}^{d_i} \sum_{j=1}^{d_j} \psi_{i,j} |i\rangle|j\rangle\]

This quantum state has a minimal decomposition, called the Schmidt decomposition, which expresses the state as a convex combination of states in two orthogonal basis:

\[|\psi\rangle = \sum_{k=1}^d s_k |\phi_k\rangle|\xi_k\rangle\]

such that $\langle\phi_k|\phi_j\rangle=\delta_{jk}$ and $\langle\xi_k|\xi_j\rangle=\delta_{jk}$, with non-negative Schmidt weights $s_kgeq 0$.

From a tensor’s perspective, what we have achieved is a decomposition of the form

\[\psi_{ij} = \sum_k \phi_{ki} s_k \xi_{kj}\]

If the number of Schmidt vectors d is small, this decomposition is advantageous. This may be the case also when the weights $|s_k|^2$ are small enough that they can be dropped above a certain size.

How are the Schmidt vectors computed? The orthogonal basis for the i and j subsystems are actually the eigenstates of the reduced density matrices. For instance

\[\rho_1 = \sum_k s_k^2 |\phi_k\rangle\langle\phi_k| = \mathrm{tr}_2|\psi\rangle\langle\psi|\]

We could thus compute those matrices, diagonalize them and obtain the tensors $\phi_{ki}$ and $\xi_{kj}$.

However, a more efficient approach is to use the singular-value decomposition or SVD, provided in Python by scipy.linalg.svd(). This algorithm recreates a matrix $A$ as the product of two isometries $U$ and $V$, and a diagonal matrix of weights $\Sigma_{ij}=s_i\delta_{ij}$. thus

\[A = U \Sigma V^T\]

Using this decomposition

\[A_{ij} = U_{ik} s_k V_{jk}\]

Now, if we identify the matrix $A$ with our original tensor $\psi$, we can identify the isometries with the Schmidt basis

\[\psi_{ij} = U_{ik} s_k V_{jk} = \phi_{ki} s_k \xi_{kj}\]