title: Keras RNN with LSTM cells for predicting multiple output time series based on multiple intput time series
body: <p>I would like to model RNN with LSTM cells in order to predict multiple output time series based on multiple input time series. To be specific, I have 4 output time series, y1[t], y2[t], y3[t], y4[t], each has a length 3,000 (t=0,...,2999). I also have 3 input time series, x1[t], x2[t], x3[t], and each has a length 3,000 sec (t=0,...,2999). The goal is to predict y1[t],.. y4[t] using all the input time series up to this current time point i.e.:</p>
  
  <pre><code>  y1[t] = f1(x1[k],x2[k],x3[k], k = 0,...,t)
    y2[t] = f2(x1[k],x2[k],x3[k], k = 0,...,t)
    y3[t] = f3(x1[k],x2[k],x3[k], k = 0,...,t)
    y4[t] = f3(x1[k],x2[k],x3[k], k = 0,...,t)
  </code></pre>
  
  <p>For a model to have a long term memory, I created a stateful RNN model by following. <a href="http://philipperemy.github.io/keras-stateful-lstm/" rel="noreferrer">keras-stateful-lstme</a>. The main difference between my case and <a href="http://philipperemy.github.io/keras-stateful-lstm/" rel="noreferrer">keras-stateful-lstme</a> is that I have:</p>
  
  <ul>
  <li>more than 1 output time series</li>
  <li>more than 1 input time series</li>
  <li>the goal is the prediction of continuous time series</li>
  </ul>
  
  <p>My code is running. However the model's prediction result is bad even with a simple data. So I would like to ask you if I am getting anything wrong.</p>
  
  <p>Here is my code with a toy example.</p>
  
  <p>In toy example, our input time series are simple cosign and sign waves:</p>
  
  <pre><code>import numpy as np
  def random_sample(len_timeseries=3000):
      Nchoice = 600
      x1 = np.cos(np.arange(0,len_timeseries)/float(1.0 + np.random.choice(Nchoice)))
      x2 = np.cos(np.arange(0,len_timeseries)/float(1.0 + np.random.choice(Nchoice)))
      x3 = np.sin(np.arange(0,len_timeseries)/float(1.0 + np.random.choice(Nchoice)))
      x4 = np.sin(np.arange(0,len_timeseries)/float(1.0 + np.random.choice(Nchoice)))
      y1 = np.random.random(len_timeseries)
      y2 = np.random.random(len_timeseries)
      y3 = np.random.random(len_timeseries)
      for t in range(3,len_timeseries):
          ## the output time series depend on input as follows: 
          y1[t] = x1[t-2] 
          y2[t] = x2[t-1]*x3[t-2]
          y3[t] = x4[t-3]
      y = np.array([y1,y2,y3]).T
      X = np.array([x1,x2,x3,x4]).T
      return y, X
  def generate_data(Nsequence = 1000):
      X_train = []
      y_train = []
      for isequence in range(Nsequence):
          y, X = random_sample()
          X_train.append(X)
          y_train.append(y)
      return np.array(X_train),np.array(y_train)
  </code></pre>
  
  <p>Please notice that y1 at time point t is simply the value of x1 at t - 2.
  Please also notice that y3 at time point t is simply the value of x1 in the two previous step.</p>
  
  <p>Using these functions, I generated 100 sets of time series y1,y2,y3,x1,x2,x3,x4. Half of them go to training data and the remaining half go to testing data.</p>
  
  <pre><code>Nsequence = 100
  prop = 0.5
  Ntrain = Nsequence*prop
  X, y = generate_data(Nsequence)
  X_train = X[:Ntrain,:,:]
  X_test  = X[Ntrain:,:,:]
  y_train = y[:Ntrain,:,:]
  y_test  = y[Ntrain:,:,:] 
  </code></pre>
  
  <p>X, y are both 3 dimensional and each contains:</p>
  
  <pre><code>#X.shape = (N sequence, length of time series, N input features)
  #y.shape = (N sequence, length of time series, N targets)
  print X.shape, y.shape
  &gt; (100, 3000, 4) (100, 3000, 3)
  </code></pre>
  
  <p>The example of the time series y1, .. y4 and x1, .., x3 are shown as below:</p>
  
  <p><a href="https://i.stack.imgur.com/u100N.png" rel="noreferrer"><img src="https://i.stack.imgur.com/u100N.png" alt="enter image description here"></a>
  <a href="https://i.stack.imgur.com/PtOvt.png" rel="noreferrer"><img src="https://i.stack.imgur.com/PtOvt.png" alt="enter image description here"></a></p>
  
  <p>I standardize these data as:</p>
  
  <pre><code>def standardize(X_train,stat=None):
      ## X_train is 3 dimentional e.g. (Nsample,len_timeseries, Nfeature)
      ## standardization is done with respect to the 3rd dimention
      if stat is None:
          featmean = np.array([np.nanmean(X_train[:,:,itrain]) for itrain in range(X_train.shape[2])]).reshape(1,1,X_train.shape[2])
          featstd = np.array([np.nanstd(X_train[:,:,itrain]) for itrain in range(X_train.shape[2])]).reshape(1,1,X_train.shape[2])
          stat = {"featmean":featmean,"featstd":featstd}
      else:
          featmean = stat["featmean"]
          featstd = stat["featstd"]
      X_train_s = (X_train - featmean)/featstd
      return X_train_s, stat 
  
  X_train_s, X_stat = standardize(X_train,stat=None)
  X_test_s, _ = standardize(X_test,stat=X_stat)
  y_train_s, y_stat = standardize(y_train,stat=None)
  y_test_s, _ = standardize(y_test,stat=y_stat)
  </code></pre>
  
  <p><strong>Create a stateful RNN model with 10 LSTM hidden neurons</strong></p>
  
  <pre><code>from keras.models import Sequential
  from keras.layers.core import Dense, Activation, Dropout
  from keras.layers.recurrent import LSTM
  def create_stateful_model(hidden_neurons):
      # create and fit the LSTM network
  
      model = Sequential()
      model.add(LSTM(hidden_neurons, 
                     batch_input_shape=(1, 1, X_train.shape[2]), 
                     return_sequences=False, 
                     stateful=True))
      model.add(Dropout(0.5))
      model.add(Dense(y_train.shape[2]))
      model.add(Activation("linear"))
      model.compile(loss='mean_squared_error', optimizer="rmsprop",metrics=['mean_squared_error'])
      return model
   model = create_stateful_model(10)
  </code></pre>
  
  <p>Now following code is used to train and validate the RNN model:</p>
  
  <pre><code>def get_R2(y_pred,y_test):
          ## y_pred_s_batch: (Nsample, len_timeseries, Noutput)
          ## the relative percentage error is computed for each output
          overall_mean = np.nanmean(y_test)
          SSres = np.nanmean( (y_pred - y_test)**2 ,axis=0).mean(axis=0)
          SStot = np.nanmean( (y_test - overall_mean)**2 ,axis=0).mean(axis=0)
          R2 = 1 - SSres / SStot 
          print "&lt;R2 testing&gt; target 1:",R2[0],"target 2:",R2[1],"target 3:",R2[2]
          return R2
  
  
  def reshape_batch_input(X_t,y_t=None):
      X_t = np.array(X_t).reshape(1,1,len(X_t)) ## (1,1,4) dimention
      if y_t is not None:
          y_t = np.array([y_t]) ## (1,3)
      return X_t,y_t
  def fit_stateful(model,X_train,y_train,X_test,y_test,nb_epoch=8):
      '''
      reference: http://philipperemy.github.io/keras-stateful-lstm/
  
      X_train: (N_time_series, len_time_series, N_features) = (10,000, 3,600 (max), 2), 
      y_train: (N_time_series, len_time_series, N_output) =   (10,000, 3,600 (max), 4)
  
      '''
      max_len = X_train.shape[1]
  
      print "X_train.shape(Nsequence =",X_train.shape[0],"len_timeseries =",X_train.shape[1],"Nfeats =",X_train.shape[2],")"
      print "y_train.shape(Nsequence =",y_train.shape[0],"len_timeseries =",y_train.shape[1],"Ntargets =",y_train.shape[2],")"
      print('Train...')
      for epoch in range(nb_epoch):
          print('___________________________________')
          print "epoch", epoch+1, "out of ",nb_epoch
          ## ---------- ##
          ##  training  ##
          ## ---------- ##
          mean_tr_acc = []
          mean_tr_loss = []
          for s in range(X_train.shape[0]):
              for t in range(max_len):
                  X_st = X_train[s][t]
                  y_st = y_train[s][t]
                  if np.any(np.isnan(y_st)):
                      break
                  X_st,y_st = reshape_batch_input(X_st,y_st)
                  tr_loss, tr_acc = model.train_on_batch(X_st,y_st)
                  mean_tr_acc.append(tr_acc)
                  mean_tr_loss.append(tr_loss)
              model.reset_states()
  
          ##print('accuracy training = {}'.format(np.mean(mean_tr_acc)))
          print('&lt;loss (mse) training&gt; {}'.format(np.mean(mean_tr_loss)))
          ## ---------- ##
          ##  testing   ##
          ## ---------- ##
          y_pred = predict_stateful(model,X_test)
          eva =  get_R2(y_pred,y_test)
      return model, eva, y_pred
  
  def predict_stateful(model,X_test):
      y_pred = []
      max_len = X_test.shape[1]
      for s in range(X_test.shape[0]):
          y_s_pred = []
          for t in range(max_len):
              X_st = X_test[s][t]
              if np.any(np.isnan(X_st)):
                  ## the rest of y is NA
                  y_s_pred.extend([np.NaN]*(max_len-len(y_s_pred)))
                  break
              X_st,_ = reshape_batch_input(X_st)
              y_st_pred = model.predict_on_batch(X_st)
              y_s_pred.append(y_st_pred[0].tolist())
  
          y_pred.append(y_s_pred)
          model.reset_states()
  
      y_pred = np.array(y_pred)
      return y_pred
  
  
  
  
    model, train_metric, y_pred = fit_stateful(model,
                                          X_train_s,y_train_s,
                                          X_test_s,y_test_s,nb_epoch=15)
  </code></pre>
  
  <p>The output is the following:</p>
  
  <pre><code>X_train.shape(Nsequence = 15 len_timeseries = 3000 Nfeats = 4 )
  y_train.shape(Nsequence = 15 len_timeseries = 3000 Ntargets = 3 )
  Train...
  ___________________________________
  epoch 1 out of  15
  &lt;loss (mse) training&gt; 0.414115458727
  &lt;R2 testing&gt; target 1: 0.664464304688 target 2: -0.574523052322 target 3: 0.526447813052
  ___________________________________
  epoch 2 out of  15
  &lt;loss (mse) training&gt; 0.394549429417
  &lt;R2 testing&gt; target 1: 0.361516087033 target 2: -0.724583671831 target 3: 0.795566178787
  ___________________________________
  epoch 3 out of  15
  &lt;loss (mse) training&gt; 0.403199136257
  &lt;R2 testing&gt; target 1: 0.09610702779 target 2: -0.468219774909 target 3: 0.69419269042
  ___________________________________
  epoch 4 out of  15
  &lt;loss (mse) training&gt; 0.406423777342
  &lt;R2 testing&gt; target 1: 0.469149270848 target 2: -0.725592048946 target 3: 0.732963522766
  ___________________________________
  epoch 5 out of  15
  &lt;loss (mse) training&gt; 0.408153116703
  &lt;R2 testing&gt; target 1: 0.400821776652 target 2: -0.329415365214 target 3: 0.2578432553
  ___________________________________
  epoch 6 out of  15
  &lt;loss (mse) training&gt; 0.421062678099
  &lt;R2 testing&gt; target 1: -0.100464591586 target 2: -0.232403824523 target 3: 0.570606489959
  ___________________________________
  epoch 7 out of  15
  &lt;loss (mse) training&gt; 0.417774856091
  &lt;R2 testing&gt; target 1: 0.320094445321 target 2: -0.606375769083 target 3: 0.349876223119
  ___________________________________
  epoch 8 out of  15
  &lt;loss (mse) training&gt; 0.427440851927
  &lt;R2 testing&gt; target 1: 0.489543715713 target 2: -0.445328806611 target 3: 0.236463139804
  ___________________________________
  epoch 9 out of  15
  &lt;loss (mse) training&gt; 0.422931671143
  &lt;R2 testing&gt; target 1: -0.31006468223 target 2: -0.322621276474 target 3: 0.122573123871
  ___________________________________
  epoch 10 out of  15
  &lt;loss (mse) training&gt; 0.43609803915
  &lt;R2 testing&gt; target 1: 0.459111316554 target 2: -0.313382405804 target 3: 0.636854743292
  ___________________________________
  epoch 11 out of  15
  &lt;loss (mse) training&gt; 0.433844655752
  &lt;R2 testing&gt; target 1: -0.0161015052703 target 2: -0.237462995323 target 3: 0.271788109459
  ___________________________________
  epoch 12 out of  15
  &lt;loss (mse) training&gt; 0.437297314405
  &lt;R2 testing&gt; target 1: -0.493665758658 target 2: -0.234236263092 target 3: 0.047264439493
  ___________________________________
  epoch 13 out of  15
  &lt;loss (mse) training&gt; 0.470605045557
  &lt;R2 testing&gt; target 1: 0.144443089961 target 2: -0.333210874982 target 3: -0.00432615142135
  ___________________________________
  epoch 14 out of  15
  &lt;loss (mse) training&gt; 0.444566756487
  &lt;R2 testing&gt; target 1: -0.053982119103 target 2: -0.0676577449316 target 3: -0.12678037186
  ___________________________________
  epoch 15 out of  15
  &lt;loss (mse) training&gt; 0.482106208801
  &lt;R2 testing&gt; target 1: 0.208482181828 target 2: -0.402982670798 target 3: 0.366757778713
  </code></pre>
  
  <p>As you can see, the training loss is NOT decreasing!!</p>
  
  <p>As the target time series 1 and 3 have very simple relations with the input time series (y1[t] = x1[t-2] , y3[t] = x4[t-3]), I would expect perfect prediction performance. However, testing R2 at every epoch shows that that is not the case. R2 at the final epoch is just about 0.2 and 0.36. Clearly, the algorithm is not converging. I am very puzzled with this result. Please do let me know what I am missing, and why the algorithm is not converging. </p>