How does the complexity in Matrix Product states ansatz drop from $D^N$ to $ND$? I have just started to read about DMRG and MPS.
It is said that in case of simple 1D chain with spins states $|\uparrow\rangle$; $|\downarrow\rangle$ and any state in the complete Hilbert space of such a system could be written as :
$|\Psi\rangle=\sum\limits_{i_1,...,i_n} 
C^{i_1,...,i_n}|i_1,...,i_n\rangle$
Where each index runs over the local basis at each site. This presents a complexity of the order $2^N$ and by writing it in the form of matrix product states we would making a simplification a mean field simplification and reduce the complexity to $2N$.
$C^{i_1,...,i_n}=C^{i_1}C^{i_2}...C^{i_n}$
EDIT:
A General case:
Suppose I have Heisenberg chain kind of a 1D Lattice problem with N sites, where each site lives on a Hilbert space of dimensionality D.
Given the fact that the dimensionality of Hilbert space of the entire system scales exponentially as $D^N$ because of the entanglement, it becomes a very complex problem to solve.
Questions:

*

*How does writing the system in terms of MPS reduce the complexity of the problem and bring it down to $ND$ ?

*What is the link between using MPS ansatz and mean field approximation?

 A: $\newcommand{\ket}[1]{\left|#1\right>}$
A general quantum state can be written as
\begin{equation}
    \ket{\psi} = \sum\limits_{\sigma_1,...,\sigma_N}\Psi_{\sigma_1,...,\sigma_N}\ket{\sigma_1}...\ket{\sigma_N} \equiv \Psi_{\sigma_1,...,\sigma_N}\ket{\sigma_1...\sigma_N} 
\end{equation}
where $\Psi_{\sigma_1...\sigma_N}$ are the coefficients that completely describe the state in the standard spin basis $\{\ket{\sigma_i}\}$. In the equations that follow we repeatedly perform SVD on $\Psi_{\sigma_1,...,\sigma_N}$ to convert the multi-dimensional tensor into a product of 3-dimensional tensors. The object that gets SVD is always a matrix and this matrix is the multi-dimensional tensor we want to SVD but reshaped. The indices found in parenthesis below are merged into one index so that a multi-dimensional tensor is reshaped into a two index matrix. The SVD takes a matrix M and gives out $M = USV^\dagger$.
\begin{align}
&\Psi_{\sigma_1,...,\sigma_N} = \Psi_{\sigma_1,(\sigma_2...\sigma_N)} = \sum\limits_{a_1 = 1}^{D_1 = min(dim(\sigma_1),dim((\sigma_2...\sigma_N)))}U_{\sigma_1,a_1}S_{a_1,a_1}V_{a_1,(\sigma_2,...,\sigma_N)}^\dagger =\\ 
& = \sum\limits_{a_1 = 1}^{D_1}U_{\sigma_1,a_1}\tilde{V}_{a_1,(\sigma_2,...,\sigma_N)}^\dagger \label{svd_repeat} = \sum\limits_{a_1 = 1}^{D_1}U_{\sigma_1,a_1}\Psi_{(a_1,\sigma_2),(\sigma_3,...,\sigma_N)} = \\
& = \sum\limits_{a_1=1}^{D_1}\sum\limits_{a_2=1}^{D_2 = min(dim(a_1,\sigma_2),dim(\sigma_3,...,\sigma_N))} U_{\sigma_1,a_1}U_{(a_1,\sigma_2),a_2}S_{a_2,a_2}V_{a_2,(\sigma_3,...,\sigma_N)}^\dagger = \\
&= \sum\limits_{a_1=1}^{D_1}\sum\limits_{a_2=1}^{D_2} U_{\sigma_1,a_1}U_{(a_1,\sigma_2),a_2}\tilde{V}_{a_2,(\sigma_3,...,\sigma_N)}^\dagger = ... = \\
& = \sum\limits_{a_1 = 1}^{D_1}...\sum\limits_{a_{N-1}=1}^{D_{N-1}}U_{\sigma_1,a_1}U_{a_1,\sigma_2,a_2}...U_{a_{N-2}\sigma_{N-1}a_{N-1}}U_{a_{N-1}\sigma_N} = \\
&\equiv \sum\limits_{a_1 = 1}^{D_1}...\sum\limits_{a_{N-1}=1}^{D_{N-1}}A_{1,a_1,\sigma_1}A_{a_1,a_2,\sigma_2}...A_{a_{N-2}a_{N-1}\sigma_{N-1}}A_{a_{N-1},1,\sigma_N} \label{svd_repeat_last}
\end{align}
The last line manifestly shows that this representation of the $\Psi_{\sigma_1,...,\sigma_N}$ coefficients forms an MPS. In this case the MPS is in left canonical form, but we could have done similar manipulations to end up with a right or mixed canonical form. Up to this point everything is exact and the MPS stores all the information contained in the coefficients $\Psi_{\sigma_1,...,\sigma_N}$. However, we can truncate the number of singular values in each $S$ matrix and keep only the $D$ largest ones. This truncation achieves the best approximation for a matrix of rank $R_{before}$ by a matrix of rank $R_{after}$ such that $R_{after} < R_{before}$. Then the last line above becomes
\begin{equation}
\label{approx_rep}
    \Psi_{\sigma_1,...,\sigma_N} \approx \sum\limits_{a_1 = 1}^{D}...\sum\limits_{a_{N-1}=1}^{D}A_{1,a_1,\sigma_1}A_{a_1,a_2,\sigma_2}...A_{a_{N-2}a_{N-1}\sigma_{N-1}}A_{a_{N-1},1,\sigma_N}
\end{equation}
If we count the number of data on each side of the above equation, on the left we have $d^N$ components and on the right $\sim N\cdot dD^2$. We have managed to reduce an exponential number of components into a polynomial number of components, thus avoiding the curse of dimensionality.
The above answers your initial question.
