Matrix product operators for sequence-to-sequence learning

The method of choice to study one-dimensional strongly interacting many-body quantum systems is based on matrix product states and operators. Such a method allows one to explore the most relevant and numerically manageable portion of an exponentially large space. It also allows one to accurately describe correlations between distant parts of a system, an important ingredient to account for the context in machine learning tasks. Here we introduce a machine learning model in which matrix product operators are trained to implement sequence-to-sequence prediction, i.e., given a sequence at a time step, it allows one to predict the next sequence. We then apply our algorithm to cellular automata (for which we show exact analytical solutions in terms of matrix product operators) and to nonlinear coupled maps. We show advantages of the proposed algorithm when compared to conditional random fields and a bidirectional, long, short-term-memory neural network. To highlight the flexibility of the algorithm, we also show that it can readily perform classification tasks.

Figure 1

In the training phase a pair of vectors ${\vec{x}}_i$ and ${\vec{y}}_i$ is converted to a pair of MPSs $X_i^{\vec{\sigma }}$ and $Y_i^{\vec{\tau }}$ and used to train an MPO $W_{\vec{\sigma }}^{\vec{\tau }}$. In the prediction phase an input vector ${\vec{x}}_k$ is converted to an MPS $X_k^{\vec{\sigma }}$ and then multiplied by the trained MPO $W_{\vec{\sigma }}^{\vec{\tau }}$. The resulting MPS ${\overline{Y}}_i^{\vec{\tau }}$ is then converted to a vector ${\vec{\overline{y}}}_k$.

Figure 2

Evolution of the long-range rule 153 cellular automata with different random initial conditions from (c) and (e) the trained MPO and (d) and (f) the respective exact evolutions. Also shown are samples of 40 (a) input and (b) output training data. The system is chosen to have $L=40$ and interaction distance $j=3$. The training parameters are training sample size $M=5000$, a bond dimension of the MPO $D_W=8$, an error rate in the training ${\epsilon }_r=0$, a regularizing coefficient $\alpha =0.001$, a maximum number of sweeps $k_{\text{max}}=20$, and convergence tolerance ${\epsilon }_t={10}^{-5}$. (g) Predicted evolution (which perfectly matches the exact one) for rule 18 with fixed boundary conditions. In this case $L=60$, the MPO bond dimension $D_W=4$, and we have used $M=50000$ training pairs of inputs and outputs. The other parameters are $\alpha =0.001$, $k_{\text{max}}=20$, and ${\epsilon }_t={10}^{-5}$.

Figure 3

Average error $\varepsilon$ between exact and predicted evolution, computed using Eq. (14), versus time for (a) $M=7000$ and MPO bond dimensions $D_W=5$ (blue solid line), $D_W=6$ (red dot-dashed line), $D_W=7$ (yellow dashed line), and $D_W=8$ (purple dotted line) and (b) $D_W=8$ with different errors ${\epsilon }_r$ in the training data and sample size $M$: ${\epsilon }_r=0.2$ and $M=3000$ (blue solid line), ${\epsilon }_r=0.2$ and $M=7000$ (red dashed line), and ${\epsilon }_r=0.4$ and $M=7000$ (yellow dot-dashed line). The evolution is given by the long-range rule 153 with $j=3$ as in Fig. 2. The other parameters are $N=100$, $k_{\text{max}}=20$, ${\epsilon }_t={10}^{-5}$, and $\alpha =0.001$.

Figure 4

Portion of the training data for (a) input and (b) output. Each of the $n$ rows is (a) an input or (b) the corresponding output sequence. The input sequence is chosen from uniformly distributed random numbers between 0 and 1 such that their sum is 1. The output sequence is computed using Eq. (15). (c)–(h) Representation of the evolution in time of an initial sequence: (c), (e), and (g) the predicted evolution from the MPO algorithm and (d), (f), and (h) the exact evolution. Each of the pairs (c) and (d), (e) and (f), and (g) and (h) corresponds to a different initial condition. The other parameters are $D_W=20$, $M=20000$, $k_{\text{max}}=20$, ${\varepsilon }_t={10}^{-5}$, and $\alpha =0.001$, while for the evolution we have used $m_1=3$, $m_2=2$, $g_1=-0.1$, and $g_2=0.5$.

Figure 5

(a) Average position $\langle l\rangle$, (b) standard deviation ${\mathrm{\Delta }}_l$, and (c) total probability $P_T$ versus time for $D_W=5$ and ${\varepsilon }_t={10}^{-5}$ (green dashed line), $D_W=20$ and ${\varepsilon }_t={10}^{-5}$ (red dot-dashed line), and $D_W=20$ and ${\varepsilon }_t={10}^{-7}$ (blue solid line). The exact solutions, given by Eq. (15), are represented by the black dotted lines. The other parameters are $N=100$, $M=20000$, $k_{\text{max}}=20$, and $\alpha =0.001$, while for the evolution we have used $m_1=3$, $m_2=2$, $g_1=-0.1$, and $g_2=0.5$.

Figure 6

Average error $\varepsilon$ between exact and predicted evolution, computed using Eq. (14), versus time for (a) $M=20000$, $\alpha =0.001$, and different bond dimensions $D_W=5$ (blue solid line), $D_W=10$ (red dashed line), and $D_W=20$ (yellow dot-dashed line); (b) $D_W=20$, $\alpha =0.001$, and different numbers of training data $M=5000$ (blue solid line), $M=10000$ (red dashed line), and $M=20000$ (yellow dot-dashed line; and (c) $D_W=20$, $M=20000$, and different regularization constants $\alpha =0.1$ (blue solid line), $\alpha =0.001$ (red dashed line), $\alpha ={10}^{-5}$ (yellow dot-dashed line), and $\alpha ={10}^{-7}$ (purple dotted line). The other parameters are $N=100$, $k_{\text{max}}=20$, and ${\epsilon }_t={10}^{-7}$. The exact evolution is given by Eq. (15) with the parameters $m_1=3$, $m_2=2$, $g_1=-0.1$, and $g_2=0.5$. (d) Comparison with BiLSTM model. The red dashed curve is the error for the BiLSTM model and the blue solid curve is for the MPO algorithm with $D_W=20$, $\alpha =0.001$, $N=100$, $k_{\text{max}}=20$, and ${\epsilon }_t={10}^{-7}$. Both algorithms have been trained with $M=20000$ samples and we have used $N=100$ initial conditions. For the BiLSTM model we have used a batch size $d_B=1$, a hidden size of LSTM cell $d_{\text{LSTM}}=100$, and learning rate $l_R=0.001$.

Figure 7

(a) Graphical representation of the linear equation (A1) in which every tensor is a structure with a different number of legs depending on the number of tensor indices (for lighter notation we have removed the indices $l$, ${\sigma }_l$, and ${\tau }_l$). The vertical lines correspond to physical indices (e.g., ${\sigma }_l$ or ${\tau }_l$) and the horizontal lines to auxiliary indices (e.g., $a_l$, $b_l$, or $c_l$). Joined legs are summed over. (b) Elementary tensors $W_{b_{l-1},b_l}^{{\sigma }_l,{\tau }_l}$, $X_{i,a_{l-1},a_l}^{{\sigma }_l}$, and $Y_{i,c_{l-1},c_l}^{{\tau }_l}$. (c) Graphical representation of intermediate tensors which are computed to evaluate more efficiently Eq. (A1).

Figure 8

Bidirectional LSTM neural network architecture. The model is composed of two rows of LSTM cells with information propagating from the left in the bottom row and from the right in the top row. The information from the LSTM cells is concatenated in $h_{i,l}^{fb}$ and then processed in the multilayer perceptrons (MLP) to return the predicted output ${\overline{y}}_{i,l}$.