Deep Learning with dependent sequences

I have a 18x1 cell array (dataA) containing 4xCols matrices, where Cols vary in length. Each cell is an annual observation of 4 variables. The sum of the first 3 should should be equal to the 4th. I treat them as 4 channels of 18 annual observations. The point is that the 4th variable should depend on the values of the 3 others and I wanted to build a deep network where this dependence would take place. That is when I make future predictions I have 1 or 2 or 3 of the first timeseries from deterministic function and I would like predict the 4th when the others will not be enough. The code is below and the question is this: Are my above intentions fulfilled? Because the results are not that promissing..

numChannels=size(dataA{1},1);

%tilerows=floor(size(dataA,1)/2);
h=figure;
tiledlayout(2,3)
for i=1:6
    nexttile
    stackedplot(dataA{i,1}',DisplayLabels=["AT","SS","Res","TG"])
    xlabel("Time in Hours")
end



% partition train and test data
numObservations = numel(dataA);
idxTrain = 1:floor(0.9*numObservations);
idxTest = floor(0.9*numObservations)+1:numObservations;
dataTrain = dataA(idxTrain);
dataTest = dataA(idxTest);

for n = 1:numel(dataTrain)
    X = dataTrain{n};
    XTrain{n} = X(:,1:end-1);
    TTrain{n} = X(:,2:end);
end

% standarize data
muX = mean(cat(2,XTrain{:}),2);
sigmaX = std(cat(2,XTrain{:}),0,2);

muT = mean(cat(2,TTrain{:}),2);
sigmaT = std(cat(2,TTrain{:}),0,2);

for n = 1:numel(XTrain)
    XTrain{n} = (XTrain{n} - muX) ./ sigmaX;
    TTrain{n} = (TTrain{n} - muT) ./ sigmaT;
end


%% Define LSTM architecture
% layers
layers = [
    sequenceInputLayer(numChannels)
    lstmLayer(28*24)
    fullyConnectedLayer(numChannels)
    regressionLayer];

% options
options = trainingOptions("adam", ...
    MaxEpochs=200, ...
    SequencePaddingDirection="left", ...
    Shuffle="every-epoch", ...
    Plots="training-progress", ...
    Verbose=0);

% train
net = trainNetwork(XTrain,TTrain,layers,options);

disp("After training.Press any key to continue.")
pause()
%%
% Test RNN
for n = 1:size(dataTest,1)
    X = dataTest{n};
    XTest{n} = (X(:,1:end-1) - muX) ./ sigmaX;
    TTest{n} = (X(:,2:end) - muT) ./ sigmaT;
end

YTest = predict(net,XTest,SequencePaddingDirection="left");

for i = 1:size(YTest,1)
    rmse(i) = sqrt(mean((YTest{i} - TTest{i}).^2,"all"));
end

figure
%histogram(rmse)
j=0;
for i=idxTest(1,1):idxTest(1,size(idxTest,2))
    j=j+1;
    TtA(j,1)=year(tA{i}(1,1));
end
bar(TtA',rmse)
xlabel("Test Years")
ylabel("RMSE")

disp('mean RMSE')
mean(rmse)

% forecast future time steps
% Given an input time series or sequence, to forecast the values of 
% multiple future time steps, use the predictAndUpdateState function 
% to predict time steps one at a time and update the RNN state at each prediction. 
% For each prediction, use the previous prediction as the input to the function.
idx = size(XTest,2);
X = XTest{idx};
T = TTest{idx};

figure
stackedplot(tA{idxTest(1,idx)}(1,1:end-1)',X',DisplayLabels=["AT","SS","Res","TG"])
xlabel("Time Step")
title("Test Observation " + year(tA{idxTest(1,idx)}(1,1)))

%% open loop forecasting

% Initialize the RNN by resetting the state 
% make an initial prediction using the first few time steps of the input data (3 days) 
% Update the RNN state using these first time steps of the input data

lfdpred=floor((size(XTest{1,idx},2)/24)/2); % lower than half of the days for initial prediction
firstDays=input("give first days for initial prediction (equal or lower than " + lfdpred + ")");

net = resetState(net);
offset = firstDays*24; 
[net,~] = predictAndUpdateState(net,X(:,1:offset));

% To forecast further predictions, loop over time steps and update the RNN state 
% using the predictAndUpdateState function. Forecast values for the remaining 
% time steps of the test observation by looping over the time steps of the input data 
% and using them as input to the RNN. The first prediction is the value 
% corresponding to the time step offset + 1.
numTimeSteps = size(X,2);
numPredictionTimeSteps = numTimeSteps - offset;
Y = zeros(numChannels,numPredictionTimeSteps);

for t = 1:numPredictionTimeSteps
    Xt = X(:,offset+t);
    [net,Y(:,t)] = predictAndUpdateState(net,Xt);
end

% compare the prediction with the target values
figure
t = tiledlayout(numChannels,1);
title(t,"Open Loop Forecasting")
lab=["AT","SS","Res","TG"];
for i = 1:numChannels
    nexttile
    plot(tA{idxTest(1,idx)}(1,2:end),T(i,:))
    hold on
    plot(tA{idxTest(1,idx)}(1,offset+1:numTimeSteps+1),[T(i,offset) Y(i,:)],'--')
    ylabel(lab(i))
end

xlabel("Time Step")
nexttile(1)
legend(["Input" "Forecasted"])

disp("press any key for closed loop forecasting")
pause()

%% closed loop forecasting

% Initialize the RNN state by first resetting the state using the resetState function, 
% then make an initial prediction Z using the first few time steps of the input data. 
% Update the RNN state using all time steps of the input data.
net = resetState(net);
offset = size(X,2);
[net,Z] = predictAndUpdateState(net,X);

% To forecast further predictions, loop over time steps and update the RNN state 
% using the predictAndUpdateState function. Forecast the next xx (e.g.200) time steps 
% by iteratively passing the previous predicted value to the RNN. 
% Because the RNN does not require the input data to make any further predictions, 
% you can specify any number of time steps to forecast.
numPredictionTimeSteps = 7*24; % 7 days
Xt = Z(:,end);
Y = zeros(numChannels,numPredictionTimeSteps);

for t = 1:numPredictionTimeSteps
    [net,Y(:,t)] = predictAndUpdateState(net,Xt);
    Xt = Y(:,t);
end

numTimeSteps = offset + numPredictionTimeSteps;

figure
t = tiledlayout(numChannels,1);
title(t,"Closed Loop Forecasting")
dt=duration(1,0,0); % increment by an hour
for i = 1:numChannels
    nexttile
    plot(tA{idxTest(1,idx)}(1,2:end),T(i,1:offset))
    hold on
    plot(tA{idxTest(1,idx)}(1,offset):dt:tA{idxTest(1,idx)}(1,offset)+(numPredictionTimeSteps)*dt,[T(i,offset) Y(i,:)],'--')
    ylabel(lab(i))
end

xlabel("Time Step")
nexttile(1)
legend(["Input" "Forecasted"])


%%
pred=[T(1:4,1:offset).*sigmaT+muT Y(1:4,:).*sigmaT+muT];
out_time(1,:)=tA{idxTest(1,idx)}(1,2):dt:tA{idxTest(1,idx)}(1,1)+numTimeSteps*dt;
out_data(1:4,:)=pred;

for i=1:size(out_data,2)
    for j=1:size(alld{21,1},2)
        if out_time(1,i)==allt{21,1}(1,j)
            in_data(1:4,i)=alld{21,1}(1:4,j);
        end
    end
end


figure
t1=tiledlayout(numChannels,1);
title(t1,"Closed loop prediction")
for i=1:numChannels
    nexttile
    plot(out_time(1,offset+1:end),in_data(i,offset+1:end))
    hold on
    plot(out_time(1,offset+1:end),out_data(i,offset+1:end),'--')
    ylabel(lab(i))
end
xlabel("time")
nexttile(1)
legend(["Raw Data" "Prediction"])

[Ben]

My understanding is you have time-series like where i denotes different observations, and . If this is always the case I would only predict with the LSTM, and use a custom layer or function layer to produce .

Here's a toy example:

% create dummy data with x(4) = x(1)+x(2)+x(3);

seqLen=10;

independentData = randn(3,seqLen);

data = cat(1, independentData, sum(independentData,1));

x = data(:,1:end-1);

y = data(:,2:end);

% create network which forces y(4) = y(1) + y(2) + y(3);

net = [sequenceInputLayer(4); lstmLayer(3); functionLayer(@(x) cat(1,x,sum(x,1)); regressionLayer];

opts = trainingOptions("adam",MaxEpochs=1);

net = trainNetwork(x,y,net,opts);

[Nestoras Papadopoulos]

Thank you very much! I'll try this.

Should I cut the sequences to smaller ones (let's say 200-400 timesteps) to produce more training sequences?

[Nestoras Papadopoulos]

In your toy example you put 3 hidden units. Is there a reason for that or is it just an arbitrary number?

I am using much more than that!

Also, in functionLayer you concatenat matrix x with the sum of it. My question is that if this is correct, because x has already the summation of the elements.

Thanks in advance.

[Ben]

The hidden size was chose to be 3 because I was assuming is dependent on 3 variables . You can use any hidden sizes before the functionLayer but the idea is that your model really only predicts and then defines . That way you can enforce this equation exactly.

If you want to model independently, but have a soft constraint that should be equal to then you can add a regularization term in a custom loss . To do that you'd use a network with a 4 dimensional output and add this custom loss term, e.g. l2loss(z, w+x+y).