Published: December 27, 2020
Author: Xinyu Chen [GitHub homepage]
Download: This Jupyter notebook is at our GitHub repository. If you want to evaluate the code, please download the notebook from the transdim repository.
This notebook shows how to implement the Bayesian Temporal Regularized Tensor Factorization (BTRTF), a fully Bayesian matrix factorization model, on some real-world data sets. To overcome the missing data problem in multivariate time series, BTRTF takes into account both low-rank matrix structure and time series autoregression. For an in-depth discussion of BTRTF, please see [1].
import numpy as np
from numpy.linalg import inv as inv
from numpy.random import normal as normrnd
from numpy.random import multivariate_normal as mvnrnd
from scipy.linalg import khatri_rao as kr_prod
from scipy.stats import wishart
from scipy.stats import invwishart
from numpy.linalg import solve as solve
from numpy.linalg import cholesky as cholesky_lower
from scipy.linalg import cholesky as cholesky_upper
from scipy.linalg import solve_triangular as solve_ut
import matplotlib.pyplot as plt
%matplotlib inline
def mvnrnd_pre(mu, Lambda):
src = normrnd(size = (mu.shape[0],))
return solve_ut(cholesky_upper(Lambda, overwrite_a = True, check_finite = False),
src, lower = False, check_finite = False, overwrite_b = True) + mu
def cov_mat(mat, mat_bar):
mat = mat - mat_bar
return mat.T @ mat
def ten2mat(tensor, mode):
return np.reshape(np.moveaxis(tensor, mode, 0), (tensor.shape[mode], -1), order = 'F')
def sample_factor_u(tau_sparse_tensor, tau_ind, U, V, X, beta0 = 1):
"""Sampling M-by-R factor matrix U and its hyperparameters (mu_u, Lambda_u)."""
dim1, rank = U.shape
U_bar = np.mean(U, axis = 0)
temp = dim1 / (dim1 + beta0)
var_mu_hyper = temp * U_bar
var_U_hyper = inv(np.eye(rank) + cov_mat(U, U_bar) + temp * beta0 * np.outer(U_bar, U_bar))
var_Lambda_hyper = wishart.rvs(df = dim1 + rank, scale = var_U_hyper)
var_mu_hyper = mvnrnd_pre(var_mu_hyper, (dim1 + beta0) * var_Lambda_hyper)
var1 = kr_prod(X, V).T
var2 = kr_prod(var1, var1)
var3 = (var2 @ ten2mat(tau_ind, 0).T).reshape([rank, rank, dim1]) + var_Lambda_hyper[:, :, None]
var4 = var1 @ ten2mat(tau_sparse_tensor, 0).T + (var_Lambda_hyper @ var_mu_hyper)[:, None]
for i in range(dim1):
U[i, :] = mvnrnd_pre(solve(var3[:, :, i], var4[:, i]), var3[:, :, i])
return U
def sample_factor_v(tau_sparse_tensor, tau_ind, U, V, X, beta0 = 1):
"""Sampling N-by-R factor matrix V and its hyperparameters (mu_v, Lambda_v)."""
dim2, rank = V.shape
V_bar = np.mean(V, axis = 0)
temp = dim2 / (dim2 + beta0)
var_mu_hyper = temp * V_bar
var_V_hyper = inv(np.eye(rank) + cov_mat(V, V_bar) + temp * beta0 * np.outer(V_bar, V_bar))
var_Lambda_hyper = wishart.rvs(df = dim2 + rank, scale = var_V_hyper)
var_mu_hyper = mvnrnd_pre(var_mu_hyper, (dim2 + beta0) * var_Lambda_hyper)
var1 = kr_prod(X, U).T
var2 = kr_prod(var1, var1)
var3 = (var2 @ ten2mat(tau_ind, 1).T).reshape([rank, rank, dim2]) + var_Lambda_hyper[:, :, None]
var4 = var1 @ ten2mat(tau_sparse_tensor, 1).T + (var_Lambda_hyper @ var_mu_hyper)[:, None]
for j in range(dim2):
V[j, :] = mvnrnd_pre(solve(var3[:, :, j], var4[:, j]), var3[:, :, j])
return V
def sample_theta(X, theta, Lambda_x, time_lags, beta0 = 1):
dim, rank = X.shape
d = time_lags.shape[0]
tmax = np.max(time_lags)
theta_bar = np.mean(theta, axis = 0)
temp = d / (d + beta0)
var_theta_hyper = inv(np.eye(rank) + cov_mat(theta, theta_bar)
+ temp * beta0 * np.outer(theta_bar, theta_bar))
var_Lambda_hyper = wishart.rvs(df = d + rank, scale = var_theta_hyper)
var_mu_hyper = mvnrnd_pre(temp * theta_bar, (d + beta0) * var_Lambda_hyper)
for k in range(d):
theta0 = theta.copy()
theta0[k, :] = 0
mat0 = np.zeros((dim - tmax, rank))
for L in range(d):
mat0 += X[tmax - time_lags[L] : dim - time_lags[L], :] @ np.diag(theta0[L, :])
varPi = X[tmax : dim, :] - mat0
var0 = X[tmax - time_lags[k] : dim - time_lags[k], :]
var = np.einsum('ij, jk, ik -> j', var0, Lambda_x, varPi)
var_Lambda = np.einsum('ti, tj, ij -> ij', var0, var0, Lambda_x) + var_Lambda_hyper
theta[k, :] = mvnrnd_pre(solve(var_Lambda, var + var_Lambda_hyper @ var_mu_hyper), var_Lambda)
return theta
def sample_Lambda_x(X, theta, time_lags):
dim, rank = X.shape
d = time_lags.shape[0]
tmax = np.max(time_lags)
mat = X[: tmax, :].T @ X[: tmax, :]
temp = np.zeros((dim - tmax, rank, d))
for k in range(d):
temp[:, :, k] = X[tmax - time_lags[k] : dim - time_lags[k], :]
new_mat = X[tmax : dim, :] - np.einsum('kr, irk -> ir', theta, temp)
Lambda_x = wishart.rvs(df = dim + rank, scale = inv(np.eye(rank) + mat + new_mat.T @ new_mat))
return Lambda_x
def sample_factor_x(tau_sparse_tensor, tau_ind, time_lags, U, V, X, theta, Lambda_x):
"""Sampling T-by-R factor matrix X."""
dim3, rank = X.shape
tmax = np.max(time_lags)
tmin = np.min(time_lags)
d = time_lags.shape[0]
A = np.zeros((d * rank, rank))
for k in range(d):
A[k * rank : (k + 1) * rank, :] = np.diag(theta[k, :])
A0 = np.dstack([A] * d)
for k in range(d):
A0[k * rank : (k + 1) * rank, :, k] = 0
mat0 = Lambda_x @ A.T
mat1 = np.einsum('kij, jt -> kit', A.reshape([d, rank, rank]), Lambda_x)
mat2 = np.einsum('kit, kjt -> ij', mat1, A.reshape([d, rank, rank]))
var1 = kr_prod(V, U).T
var2 = kr_prod(var1, var1)
var3 = (var2 @ ten2mat(tau_ind, 2).T).reshape([rank, rank, dim3]) + Lambda_x[:, :, None]
var4 = var1 @ ten2mat(tau_sparse_tensor, 2).T
for t in range(dim3):
Mt = np.zeros((rank, rank))
Nt = np.zeros(rank)
Qt = mat0 @ X[t - time_lags, :].reshape(rank * d)
index = list(range(0, d))
if t >= dim3 - tmax and t < dim3 - tmin:
index = list(np.where(t + time_lags < dim3))[0]
elif t < tmax:
Qt = np.zeros(rank)
index = list(np.where(t + time_lags >= tmax))[0]
if t < dim3 - tmin:
Mt = mat2.copy()
temp = np.zeros((rank * d, len(index)))
n = 0
for k in index:
temp[:, n] = X[t + time_lags[k] - time_lags, :].reshape(rank * d)
n += 1
temp0 = X[t + time_lags[index], :].T - np.einsum('ijk, ik -> jk', A0[:, :, index], temp)
Nt = np.einsum('kij, jk -> i', mat1[index, :, :], temp0)
var3[:, :, t] = var3[:, :, t] + Mt
if t < tmax:
var3[:, :, t] = var3[:, :, t] - Lambda_x + np.eye(rank)
X[t, :] = mvnrnd_pre(solve(var3[:, :, t], var4[:, t] + Nt + Qt), var3[:, :, t])
return X
def compute_mape(var, var_hat):
return np.sum(np.abs(var - var_hat) / var) / var.shape[0]
def compute_rmse(var, var_hat):
return np.sqrt(np.sum((var - var_hat) ** 2) / var.shape[0])
def ar4cast(theta, X, Lambda_x, time_lags, multi_step):
dim, rank = X.shape
d = time_lags.shape[0]
X_new = np.append(X, np.zeros((multi_step, rank)), axis = 0)
for t in range(multi_step):
X_new[dim + t, :] = mvnrnd_pre(np.einsum('kr, kr -> r', theta, X_new[dim + t - time_lags, :]), Lambda_x)
return X_new
def BTRTF(dense_tensor, sparse_tensor, init, rank, time_lags, burn_iter, gibbs_iter, multi_step = 1):
"""Bayesian Temporal Regularized Tensor Factorization, BTRTF."""
dim1, dim2, dim3 = sparse_tensor.shape
d = time_lags.shape[0]
U = init["U"]
V = init["V"]
X = init["X"]
theta = 0.01 * np.random.randn(d, rank)
if np.isnan(sparse_tensor).any() == False:
ind = sparse_tensor != 0
pos_obs = np.where(ind)
pos_test = np.where((dense_tensor != 0) & (sparse_tensor == 0))
elif np.isnan(sparse_tensor).any() == True:
pos_test = np.where((dense_tensor != 0) & (np.isnan(sparse_tensor)))
ind = ~np.isnan(sparse_tensor)
pos_obs = np.where(ind)
sparse_tensor[np.isnan(sparse_tensor)] = 0
dense_test = dense_tensor[pos_test]
del dense_tensor
U_plus = np.zeros((dim1, rank, gibbs_iter))
V_plus = np.zeros((dim2, rank, gibbs_iter))
X_plus = np.zeros((dim3 + multi_step, rank, gibbs_iter))
theta_plus = np.zeros((d, rank, gibbs_iter))
tau_plus = np.zeros(gibbs_iter)
Lambda_plus = np.zeros((rank, rank, gibbs_iter))
temp_hat = np.zeros(len(pos_test[0]))
show_iter = 500
tau = 1
tensor_hat_plus = np.zeros(sparse_tensor.shape)
tensor_new_plus = np.zeros((dim1, dim2, multi_step))
for it in range(burn_iter + gibbs_iter):
tau_ind = tau * ind
tau_sparse_tensor = tau * sparse_tensor
U = sample_factor_u(tau_sparse_tensor, tau_ind, U, V, X)
V = sample_factor_v(tau_sparse_tensor, tau_ind, U, V, X)
Lambda_x = sample_Lambda_x(X, theta, time_lags)
theta = sample_theta(X, theta, Lambda_x, time_lags)
X = sample_factor_x(tau_sparse_tensor, tau_ind, time_lags, U, V, X, theta, Lambda_x)
tensor_hat = np.einsum('is, js, ts -> ijt', U, V, X)
tau = np.random.gamma(1e-6 + 0.5 * np.sum(ind),
1 / (1e-6 + 0.5 * np.sum(((sparse_tensor - tensor_hat) ** 2) * ind)))
temp_hat += tensor_hat[pos_test]
if (it + 1) % show_iter == 0 and it < burn_iter:
temp_hat = temp_hat / show_iter
print('Iter: {}'.format(it + 1))
print('MAPE: {:.6}'.format(compute_mape(dense_test, temp_hat)))
print('RMSE: {:.6}'.format(compute_rmse(dense_test, temp_hat)))
temp_hat = np.zeros(len(pos_test[0]))
print()
if it + 1 > burn_iter:
U_plus[:, :, it - burn_iter] = U
V_plus[:, :, it - burn_iter] = V
theta_plus[:, :, it - burn_iter] = theta
Lambda_plus[:, :, it - burn_iter] = Lambda_x
tau_plus[it - burn_iter] = tau
tensor_hat_plus += tensor_hat
X0 = ar4cast(theta, X, Lambda_x, time_lags, multi_step)
X_plus[:, :, it - burn_iter] = X0
tensor_new_plus += np.einsum('is, js, ts -> ijt', U, V, X0[- multi_step :, :])
tensor_hat = tensor_hat_plus / gibbs_iter
print('Imputation MAPE: {:.6}'.format(compute_mape(dense_test, tensor_hat[:, :, : dim3][pos_test])))
print('Imputation RMSE: {:.6}'.format(compute_rmse(dense_test, tensor_hat[:, :, : dim3][pos_test])))
print()
tensor_hat = np.append(tensor_hat, tensor_new_plus / gibbs_iter, axis = 2)
tensor_hat[tensor_hat < 0] = 0
return tensor_hat, U_plus, V_plus, X_plus, theta_plus, Lambda_plus, tau_plus
def sample_factor_x_partial(tau_sparse_tensor, tau_ind, time_lags, U, V, X, theta, Lambda_x, back_step):
"""Sampling T-by-R factor matrix X."""
dim3, rank = X.shape
tmax = np.max(time_lags)
tmin = np.min(time_lags)
d = time_lags.shape[0]
A = np.zeros((d * rank, rank))
for k in range(d):
A[k * rank : (k + 1) * rank, :] = np.diag(theta[k, :])
A0 = np.dstack([A] * d)
for k in range(d):
A0[k * rank : (k + 1) * rank, :, k] = 0
mat0 = Lambda_x @ A.T
mat1 = np.einsum('kij, jt -> kit', A.reshape([d, rank, rank]), Lambda_x)
mat2 = np.einsum('kit, kjt -> ij', mat1, A.reshape([d, rank, rank]))
var1 = kr_prod(V, U).T
var2 = kr_prod(var1, var1)
var3 = (var2 @ ten2mat(tau_ind[:, :, - back_step :], 2).T).reshape([rank, rank, back_step]) + Lambda_x[:, :, None]
var4 = var1 @ ten2mat(tau_sparse_tensor[:, :, - back_step :], 2).T
for t in range(dim3 - back_step, dim3):
Mt = np.zeros((rank, rank))
Nt = np.zeros(rank)
Qt = mat0 @ X[t - time_lags, :].reshape(rank * d)
index = list(range(0, d))
if t >= dim3 - tmax and t < dim3 - tmin:
index = list(np.where(t + time_lags < dim3))[0]
if t < dim3 - tmin:
Mt = mat2.copy()
temp = np.zeros((rank * d, len(index)))
n = 0
for k in index:
temp[:, n] = X[t + time_lags[k] - time_lags, :].reshape(rank * d)
n += 1
temp0 = X[t + time_lags[index], :].T - np.einsum('ijk, ik -> jk', A0[:, :, index], temp)
Nt = np.einsum('kij, jk -> i', mat1[index, :, :], temp0)
var3[:, :, t + back_step - dim3] = var3[:, :, t + back_step - dim3] + Mt
X[t, :] = mvnrnd_pre(solve(var3[:, :, t + back_step - dim3],
var4[:, t + back_step - dim3] + Nt + Qt), var3[:, :, t + back_step - dim3])
return X
def BTRTF_partial(dense_tensor, sparse_tensor, init, rank, time_lags, burn_iter, gibbs_iter, multi_step = 1, gamma = 10):
"""Bayesian Temporal Regularized Tensor Factorization, BTRTF."""
dim1, dim2, dim3 = sparse_tensor.shape
U_plus = init["U_plus"]
V_plus = init["V_plus"]
X_plus = init["X_plus"]
theta_plus = init["theta_plus"]
Lambda_plus = init["Lambda_plus"]
tau_plus = init["tau_plus"]
if np.isnan(sparse_tensor).any() == False:
ind = sparse_tensor != 0
pos_obs = np.where(ind)
elif np.isnan(sparse_tensor).any() == True:
ind = ~np.isnan(sparse_tensor)
pos_obs = np.where(ind)
sparse_tensor[np.isnan(sparse_tensor)] = 0
X_new_plus = np.zeros((dim3 + multi_step, rank, gibbs_iter))
tensor_new_plus = np.zeros((dim1, dim2, multi_step))
back_step = gamma * multi_step
for it in range(gibbs_iter):
tau_ind = tau_plus[it] * ind
tau_sparse_tensor = tau_plus[it] * sparse_tensor
X = sample_factor_x_partial(tau_sparse_tensor, tau_ind, time_lags, U_plus[:, :, it], V_plus[:, :, it],
X_plus[:, :, it], theta_plus[:, :, it], Lambda_plus[:, :, it], back_step)
X0 = ar4cast(theta_plus[:, :, it], X, Lambda_plus[:, :, it], time_lags, multi_step)
X_new_plus[:, :, it] = X0
tensor_new_plus += np.einsum('is, js, ts -> ijt', U_plus[:, :, it], V_plus[:, :, it], X0[- multi_step :, :])
tensor_hat = tensor_new_plus / gibbs_iter
tensor_hat[tensor_hat < 0] = 0
return tensor_hat, U_plus, V_plus, X_new_plus, theta_plus, Lambda_plus, tau_plus
from ipywidgets import IntProgress
from IPython.display import display
def BTRTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter, gamma = 10):
dim1, dim2, T = dense_tensor.shape
start_time = T - pred_step
max_count = int(np.ceil(pred_step / multi_step))
tensor_hat = np.zeros((dim1, dim2, max_count * multi_step))
f = IntProgress(min = 0, max = max_count) # instantiate the bar
display(f) # display the bar
for t in range(max_count):
if t == 0:
init = {"U": 0.1 * np.random.randn(dim1, rank),
"V": 0.1 * np.random.randn(dim2, rank),
"X": 0.1 * np.random.randn(start_time, rank)}
tensor, U, V, X_new, theta, Lambda_x, tau = BTRTF(dense_tensor[:, :, : start_time],
sparse_tensor[:, :, : start_time], init, rank, time_lags, burn_iter, gibbs_iter, multi_step)
else:
init = {"U_plus": U, "V_plus": V, "X_plus": X_new, "theta_plus": theta, "Lambda_plus": Lambda_x, "tau_plus": tau}
tensor, U, V, X_new, theta, Lambda_x, tau = BTRTF_partial(dense_tensor[:, :, : start_time + t * multi_step],
sparse_tensor[:, :, : start_time + t * multi_step], init,
rank, time_lags, burn_iter, gibbs_iter, multi_step, gamma)
tensor_hat[:, :, t * multi_step : (t + 1) * multi_step] = tensor[:, :, - multi_step :]
f.value = t
small_dense_tensor = dense_tensor[:, :, start_time : T]
pos = np.where(small_dense_tensor != 0)
print('Prediction MAPE: {:.6}'.format(compute_mape(small_dense_tensor[pos], tensor_hat[pos])))
print('Prediction RMSE: {:.6}'.format(compute_rmse(small_dense_tensor[pos], tensor_hat[pos])))
print()
return tensor_hat
Scenario setting:
import scipy.io
import warnings
warnings.simplefilter('ignore')
dense_tensor = scipy.io.loadmat('../datasets/NYC-data-set/tensor.mat')['tensor'].astype(np.float32)
sparse_tensor = dense_tensor.copy()
Model setting:
import time
rank = 30
pred_step = 7 * 24
time_lags = np.array([1, 2, 3, 24, 25, 26, 7 * 24, 7 * 24 + 1, 7 * 24 + 2])
burn_iter = 1000
gibbs_iter = 200
for multi_step in [2, 4, 6]:
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTRTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step, rank, time_lags, burn_iter, gibbs_iter)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 2.
IntProgress(value=0, max=84)
Iter: 500 MAPE: nan RMSE: nan Iter: 1000 MAPE: nan RMSE: nan Imputation MAPE: nan Imputation RMSE: nan Prediction MAPE: 0.507437 Prediction RMSE: 5.46837 Running time: 1665 seconds Prediction time horizon (delta) = 4.
IntProgress(value=0, max=42)
Iter: 500 MAPE: nan RMSE: nan Iter: 1000 MAPE: nan RMSE: nan Imputation MAPE: nan Imputation RMSE: nan Prediction MAPE: 0.508875 Prediction RMSE: 5.71892 Running time: 1594 seconds Prediction time horizon (delta) = 6.
IntProgress(value=0, max=28)
Iter: 500 MAPE: nan RMSE: nan Iter: 1000 MAPE: nan RMSE: nan Imputation MAPE: nan Imputation RMSE: nan Prediction MAPE: 0.503395 Prediction RMSE: 5.9437 Running time: 1638 seconds
Scenario setting:
import scipy.io
import numpy as np
np.random.seed(1000)
dense_tensor = scipy.io.loadmat('../datasets/NYC-data-set/tensor.mat')['tensor'].astype(np.float32)
dim = dense_tensor.shape
missing_rate = 0.4 # Random missing (RM)
sparse_tensor = dense_tensor * np.round(np.random.rand(dim[0], dim[1], dim[2]) + 0.5 - missing_rate)
Model setting:
import time
rank = 30
pred_step = 7 * 24
time_lags = np.array([1, 2, 3, 24, 25, 26, 7 * 24, 7 * 24 + 1, 7 * 24 + 2])
burn_iter = 1000
gibbs_iter = 200
for multi_step in [2, 4, 6]:
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTRTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step, rank, time_lags, burn_iter, gibbs_iter)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 2.
IntProgress(value=0, max=84)
Iter: 500 MAPE: 0.478202 RMSE: 4.84207 Iter: 1000 MAPE: 0.484866 RMSE: 4.93907 Imputation MAPE: 0.484998 Imputation RMSE: 4.9549 Prediction MAPE: 0.485314 Prediction RMSE: 6.46555 Running time: 1516 seconds Prediction time horizon (delta) = 4.
IntProgress(value=0, max=42)
Iter: 500 MAPE: 0.511286 RMSE: 4.78762 Iter: 1000 MAPE: 0.52187 RMSE: 4.85227 Imputation MAPE: 0.523542 Imputation RMSE: 4.86575 Prediction MAPE: 0.506951 Prediction RMSE: 5.75299 Running time: 1396 seconds Prediction time horizon (delta) = 6.
IntProgress(value=0, max=28)
Iter: 500 MAPE: 0.50354 RMSE: 4.76326 Iter: 1000 MAPE: 0.513024 RMSE: 4.81782 Imputation MAPE: 0.514024 Imputation RMSE: 4.82617 Prediction MAPE: 0.508809 Prediction RMSE: 5.91545 Running time: 1822 seconds
Scenario setting:
import scipy.io
import numpy as np
np.random.seed(1000)
dense_tensor = scipy.io.loadmat('../datasets/NYC-data-set/tensor.mat')['tensor'].astype(np.float32)
dim = dense_tensor.shape
missing_rate = 0.6 # Random missing (RM)
sparse_tensor = dense_tensor * np.round(np.random.rand(dim[0], dim[1], dim[2]) + 0.5 - missing_rate)
Model setting:
import time
rank = 30
pred_step = 7 * 24
time_lags = np.array([1, 2, 3, 24, 25, 26, 7 * 24, 7 * 24 + 1, 7 * 24 + 2])
burn_iter = 1000
gibbs_iter = 200
for multi_step in [2, 4, 6]:
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTRTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step, rank, time_lags, burn_iter, gibbs_iter)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 2.
IntProgress(value=0, max=84)
Iter: 500 MAPE: 0.486001 RMSE: 5.26879 Iter: 1000 MAPE: 0.491243 RMSE: 5.26529 Imputation MAPE: 0.491764 Imputation RMSE: 5.30838 Prediction MAPE: 0.493182 Prediction RMSE: 7.27568 Running time: 1628 seconds Prediction time horizon (delta) = 4.
IntProgress(value=0, max=42)
Iter: 500 MAPE: 0.51287 RMSE: 5.07456 Iter: 1000 MAPE: 0.52265 RMSE: 5.17151 Imputation MAPE: 0.525277 Imputation RMSE: 5.18165 Prediction MAPE: 0.505567 Prediction RMSE: 5.94686 Running time: 1808 seconds Prediction time horizon (delta) = 6.
IntProgress(value=0, max=28)
Iter: 500 MAPE: 0.516195 RMSE: 5.02277 Iter: 1000 MAPE: 0.522737 RMSE: 5.08686 Imputation MAPE: 0.520651 Imputation RMSE: 5.14705 Prediction MAPE: 0.50992 Prediction RMSE: 6.15616 Running time: 1347 seconds
Scenario setting:
import scipy.io
import numpy as np
np.random.seed(1000)
dense_tensor = scipy.io.loadmat('../datasets/NYC-data-set/tensor.mat')['tensor']
dim = dense_tensor.shape
nm_tensor = np.random.rand(dim[0], dim[1], dim[2])
missing_rate = 0.4 # Non-random missing (NM)
binary_tensor = np.zeros(dense_tensor.shape)
for i1 in range(dim[0]):
for i2 in range(dim[1]):
for i3 in range(61):
binary_tensor[i1, i2, i3 * 24 : (i3 + 1) * 24] = np.round(nm_tensor[i1, i2, i3] + 0.5 - missing_rate)
sparse_tensor = dense_tensor * binary_tensor
Model setting:
import time
rank = 30
pred_step = 7 * 24
time_lags = np.array([1, 2, 3, 24, 25, 26, 7 * 24, 7 * 24 + 1, 7 * 24 + 2])
burn_iter = 1000
gibbs_iter = 200
for multi_step in [2, 4, 6]:
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTRTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step, rank, time_lags, burn_iter, gibbs_iter)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 2.
IntProgress(value=0, max=84)
Iter: 500 MAPE: 0.47176 RMSE: 4.99266 Iter: 1000 MAPE: 0.478566 RMSE: 5.16785 Imputation MAPE: 0.478796 Imputation RMSE: 5.19043 Prediction MAPE: 0.485589 Prediction RMSE: 6.30386 Running time: 1599 seconds Prediction time horizon (delta) = 4.
IntProgress(value=0, max=42)
Iter: 500 MAPE: 0.509326 RMSE: 4.92351 Iter: 1000 MAPE: 0.52759 RMSE: 5.02248 Imputation MAPE: 0.527748 Imputation RMSE: 5.04027 Prediction MAPE: 0.503853 Prediction RMSE: 5.87668 Running time: 1518 seconds Prediction time horizon (delta) = 6.
IntProgress(value=0, max=28)
Iter: 500 MAPE: 0.51339 RMSE: 4.87887 Iter: 1000 MAPE: 0.525974 RMSE: 5.02457 Imputation MAPE: 0.528036 Imputation RMSE: 5.03888 Prediction MAPE: 0.500864 Prediction RMSE: 6.11974 Running time: 1512 seconds
Scenario setting:
import numpy as np
import warnings
warnings.simplefilter('ignore')
dense_tensor = np.load('../datasets/Temperature-data-set/tensor.npy').astype(np.float32)
pos = np.where(dense_tensor[:, 0, :] > 50)
dense_tensor[pos[0], :, pos[1]] = 0
sparse_tensor = dense_tensor.copy()
import time
rank = 30
pred_step = 10 * 12
time_lags = np.array([1, 2, 3, 12, 13, 14, 2 * 12, 2 * 12 + 1, 2 * 12 + 2])
burn_iter = 1000
gibbs_iter = 200
multi_step = 2
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTRTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter, gamma = 30)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 2.
IntProgress(value=0, max=60)
Iter: 500 MAPE: nan RMSE: nan Iter: 1000 MAPE: nan RMSE: nan Imputation MAPE: nan Imputation RMSE: nan Prediction MAPE: 0.163308 Prediction RMSE: 5.58539 Running time: 849 seconds
import time
rank = 30
pred_step = 10 * 12
time_lags = np.array([1, 2, 3, 12, 13, 14, 2 * 12, 2 * 12 + 1, 2 * 12 + 2])
burn_iter = 1000
gibbs_iter = 200
for multi_step in [4, 6]:
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTRTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter, gamma = 30)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 4.
IntProgress(value=0, max=30)
Iter: 500 MAPE: nan RMSE: nan Iter: 1000 MAPE: nan RMSE: nan Imputation MAPE: nan Imputation RMSE: nan Prediction MAPE: 0.219357 Prediction RMSE: 7.52495 Running time: 676 seconds Prediction time horizon (delta) = 6.
IntProgress(value=0, max=20)
Iter: 500 MAPE: nan RMSE: nan Iter: 1000 MAPE: nan RMSE: nan Imputation MAPE: nan Imputation RMSE: nan Prediction MAPE: 0.172542 Prediction RMSE: 5.59173 Running time: 642 seconds
Scenario setting:
import numpy as np
np.random.seed(1000)
dense_tensor = np.load('../datasets/Temperature-data-set/tensor.npy').astype(np.float32)
pos = np.where(dense_tensor[:, 0, :] > 50)
dense_tensor[pos[0], :, pos[1]] = 0
random_tensor = np.random.rand(dense_tensor.shape[0], dense_tensor.shape[1], dense_tensor.shape[2])
missing_rate = 0.4
## Random missing (RM)
binary_tensor = np.round(random_tensor + 0.5 - missing_rate)
sparse_tensor = dense_tensor.copy()
sparse_tensor[binary_tensor == 0] = np.nan
sparse_tensor[sparse_tensor == 0] = np.nan
import time
rank = 30
pred_step = 10 * 12
time_lags = np.array([1, 2, 3, 12, 13, 14, 2 * 12, 2 * 12 + 1, 2 * 12 + 2])
burn_iter = 1000
gibbs_iter = 200
multi_step = 2
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTRTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter, gamma = 30)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 2.
IntProgress(value=0, max=60)
Iter: 500 MAPE: 0.0146505 RMSE: 0.488159 Iter: 1000 MAPE: 0.0146261 RMSE: 0.484394 Imputation MAPE: 0.0145402 Imputation RMSE: 0.480634 Prediction MAPE: 0.100703 Prediction RMSE: 3.29517 Running time: 759 seconds
import time
rank = 30
pred_step = 10 * 12
time_lags = np.array([1, 2, 3, 12, 13, 14, 2 * 12, 2 * 12 + 1, 2 * 12 + 2])
burn_iter = 1000
gibbs_iter = 200
for multi_step in [4, 6]:
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTRTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter, gamma = 30)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 4.
IntProgress(value=0, max=30)
Iter: 500 MAPE: 0.0147999 RMSE: 0.495504 Iter: 1000 MAPE: 0.0145005 RMSE: 0.482804 Imputation MAPE: 0.0144735 Imputation RMSE: 0.481308 Prediction MAPE: 0.146639 Prediction RMSE: 4.71453 Running time: 691 seconds Prediction time horizon (delta) = 6.
IntProgress(value=0, max=20)
Iter: 500 MAPE: 0.0148263 RMSE: 0.492939 Iter: 1000 MAPE: 0.0144807 RMSE: 0.479886 Imputation MAPE: 0.0144055 Imputation RMSE: 0.476987 Prediction MAPE: 0.145679 Prediction RMSE: 4.71419 Running time: 664 seconds
Scenario setting:
import numpy as np
np.random.seed(1000)
dense_tensor = np.load('../datasets/Temperature-data-set/tensor.npy').astype(np.float32)
pos = np.where(dense_tensor[:, 0, :] > 50)
dense_tensor[pos[0], :, pos[1]] = 0
random_tensor = np.random.rand(dense_tensor.shape[0], dense_tensor.shape[1], dense_tensor.shape[2])
missing_rate = 0.6
## Random missing (RM)
binary_tensor = np.round(random_tensor + 0.5 - missing_rate)
sparse_tensor = dense_tensor.copy()
sparse_tensor[binary_tensor == 0] = np.nan
sparse_tensor[sparse_tensor == 0] = np.nan
import time
rank = 30
pred_step = 10 * 12
time_lags = np.array([1, 2, 3, 12, 13, 14, 2 * 12, 2 * 12 + 1, 2 * 12 + 2])
burn_iter = 1000
gibbs_iter = 200
multi_step = 2
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTRTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter, gamma = 30)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 2.
IntProgress(value=0, max=60)
Iter: 500 MAPE: 0.0150097 RMSE: 0.502107 Iter: 1000 MAPE: 0.0148378 RMSE: 0.49437 Imputation MAPE: 0.0147177 Imputation RMSE: 0.489617 Prediction MAPE: 0.175718 Prediction RMSE: 5.57326 Running time: 790 seconds
import time
rank = 30
pred_step = 10 * 12
time_lags = np.array([1, 2, 3, 12, 13, 14, 2 * 12, 2 * 12 + 1, 2 * 12 + 2])
burn_iter = 1000
gibbs_iter = 200
for multi_step in [4, 6]:
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTRTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter, gamma = 30)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 4.
IntProgress(value=0, max=30)
Iter: 500 MAPE: 0.0149364 RMSE: 0.496529 Iter: 1000 MAPE: 0.0146882 RMSE: 0.486221 Imputation MAPE: 0.0146582 Imputation RMSE: 0.485313 Prediction MAPE: 0.133308 Prediction RMSE: 4.16352 Running time: 704 seconds Prediction time horizon (delta) = 6.
IntProgress(value=0, max=20)
Iter: 500 MAPE: 0.0148235 RMSE: 0.493898 Iter: 1000 MAPE: 0.0144327 RMSE: 0.480168 Imputation MAPE: 0.0143945 Imputation RMSE: 0.478618 Prediction MAPE: 0.0907117 Prediction RMSE: 3.1397 Running time: 674 seconds
Scenario setting:
import numpy as np
np.random.seed(1000)
dense_tensor = np.load('../datasets/Temperature-data-set/tensor.npy').astype(np.float32)
pos = np.where(dense_tensor[:, 0, :] > 50)
dense_tensor[pos[0], :, pos[1]] = 0
random_tensor = np.random.rand(dense_tensor.shape[0], dense_tensor.shape[1], int(dense_tensor.shape[2] / 3))
missing_rate = 0.4
## Non-random missing (NM)
binary_tensor = np.zeros(dense_tensor.shape)
for i1 in range(dense_tensor.shape[0]):
for i2 in range(dense_tensor.shape[1]):
for i3 in range(int(dense_tensor.shape[2] / 3)):
binary_tensor[i1, i2, i3 * 3 : (i3 + 1) * 3] = np.round(random_tensor[i1, i2, i3] + 0.5 - missing_rate)
sparse_tensor = dense_tensor.copy()
sparse_tensor[binary_tensor == 0] = np.nan
sparse_tensor[sparse_tensor == 0] = np.nan
import time
rank = 30
pred_step = 10 * 12
time_lags = np.array([1, 2, 3, 12, 13, 14, 2 * 12, 2 * 12 + 1, 2 * 12 + 2])
burn_iter = 1000
gibbs_iter = 200
multi_step = 2
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTRTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter, gamma = 30)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 2.
IntProgress(value=0, max=60)
Iter: 500 MAPE: 0.0144751 RMSE: 0.481766 Iter: 1000 MAPE: 0.0144173 RMSE: 0.477133 Imputation MAPE: 0.0143819 Imputation RMSE: 0.476134 Prediction MAPE: 0.285081 Prediction RMSE: 9.44928 Running time: 788 seconds
import time
rank = 30
pred_step = 10 * 12
time_lags = np.array([1, 2, 3, 12, 13, 14, 2 * 12, 2 * 12 + 1, 2 * 12 + 2])
burn_iter = 1000
gibbs_iter = 200
for multi_step in [4, 6]:
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTRTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter, gamma = 30)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 4.
IntProgress(value=0, max=30)
Iter: 500 MAPE: 0.0147757 RMSE: 0.491789 Iter: 1000 MAPE: 0.0146142 RMSE: 0.484462 Imputation MAPE: 0.0145513 Imputation RMSE: 0.481513 Prediction MAPE: 0.103239 Prediction RMSE: 3.34341 Running time: 708 seconds Prediction time horizon (delta) = 6.
IntProgress(value=0, max=20)
Iter: 500 MAPE: 0.0149642 RMSE: 0.497966 Iter: 1000 MAPE: 0.0147473 RMSE: 0.487566 Imputation MAPE: 0.0146458 Imputation RMSE: 0.484409 Prediction MAPE: 0.296197 Prediction RMSE: 9.49325 Running time: 757 seconds