Temporal Regularized Matrix Factorization (TRMF) is an effective tool for imputing missing data within a given multivariate time series and forecasting time series with missing values. This approach is from the following literature:

Hsiang-Fu Yu, Nikhil Rao, Inderjit S. Dhillon, 2016. Temporal regularized matrix factorization for high-dimensional time series prediction. 30th Conference on Neural Information Processing Systems (NIPS 2016), Barcelona, Spain.

## Quick Run¶

This notebook is publicly available for any usage at our data imputation project. Please click transdim.

## Data Organization: Matrix Structure¶

In this post, we consider a dataset of $m$ discrete time series $\boldsymbol{y}_{i}\in\mathbb{R}^{f},i\in\left\{1,2,...,m\right\}$. The time series may have missing elements. We express spatio-temporal dataset as a matrix $Y\in\mathbb{R}^{m\times f}$ with $m$ rows (e.g., locations) and $f$ columns (e.g., discrete time intervals),

$$Y=\left[ \begin{array}{cccc} y_{11} & y_{12} & \cdots & y_{1f} \\ y_{21} & y_{22} & \cdots & y_{2f} \\ \vdots & \vdots & \ddots & \vdots \\ y_{m1} & y_{m2} & \cdots & y_{mf} \\ \end{array} \right]\in\mathbb{R}^{m\times f}.$$

## TRMF model¶

Temporal Regularized Matrix Factorization (TRMF) is an approach to incorporate temporal dependencies into commonly-used matrix factorization model. The temporal dependencies are described among ${\boldsymbol{x}_t}$ explicitly. Such approach takes the form:

$$\boldsymbol{x}_{t}\approx\sum_{l\in\mathcal{L}}\boldsymbol{\theta}_{l}\circledast\boldsymbol{x}_{t-l},$$

where this autoregressive (AR) is specialized by a lag set $\mathcal{L}=\left\{l_1,l_2,...,l_d\right\}$ (e.g., $\mathcal{L}=\left\{1,2,144\right\}$) and weights $\boldsymbol{\theta}_{l}\in\mathbb{R}^{r},\forall l$, and we further define

$$\mathcal{R}_{AR}\left(X\mid \mathcal{L},\Theta,\eta\right)=\frac{1}{2}\sum_{t=l_d+1}^{f}\left(\boldsymbol{x}_{t}-\sum_{l\in\mathcal{L}}\boldsymbol{\theta}_{l}\circledast\boldsymbol{x}_{t-l}\right)^\top\left(\boldsymbol{x}_{t}-\sum_{l\in\mathcal{L}}\boldsymbol{\theta}_{l}\circledast\boldsymbol{x}_{t-l}\right)+\frac{\eta}{2}\sum_{t=1}^{f}\boldsymbol{x}_{t}^\top\boldsymbol{x}_{t}.$$

Thus, TRMF-AR is given by solving

$$\min_{W,X,\Theta}\frac{1}{2}\underbrace{\sum_{(i,t)\in\Omega}\left(y_{it}-\boldsymbol{w}_{i}^T\boldsymbol{x}_{t}\right)^2}_{\text{sum of squared residual errors}}+\lambda_{w}\underbrace{\mathcal{R}_{w}\left(W\right)}_{W-\text{regularizer}}+\lambda_{x}\underbrace{\mathcal{R}_{AR}\left(X\mid \mathcal{L},\Theta,\eta\right)}_{\text{AR-regularizer}}+\lambda_{\theta}\underbrace{\mathcal{R}_{\theta}\left(\Theta\right)}_{\Theta-\text{regularizer}}$$

where $\mathcal{R}_{w}\left(W\right)=\frac{1}{2}\sum_{i=1}^{m}\boldsymbol{w}_{i}^\top\boldsymbol{w}_{i}$ and $\mathcal{R}_{\theta}\left(\Theta\right)=\frac{1}{2}\sum_{l\in\mathcal{L}}\boldsymbol{\theta}_{l}^\top\boldsymbol{\theta}_{l}$ are regularization terms.

### Define TRMF model with Numpy¶

Observing the optimization problem of TRMF model as mentioned above, we categorize the parameters within this model as parameters (i.e., init_para in the TRMF function) and hyperparameters (i.e., init_hyper).

• Parameters include spatial matrix $W$, temporal matrix $X$, and AR coefficients $\Theta$.
• Hyperparameters include weight parameters on some regularizers, i.e., $\lambda_w$, $\lambda_x$, $\lambda_\theta$, and $\eta$.
In [1]:
import numpy as np
from numpy.linalg import inv as inv

def TRMF(dense_mat, sparse_mat, init_para, init_hyper, time_lags, maxiter):
"""Temporal Regularized Matrix Factorization, TRMF."""

## Initialize parameters
W = init_para["W"]
X = init_para["X"]
theta = init_para["theta"]

## Set hyperparameters
lambda_w = init_hyper["lambda_w"]
lambda_x = init_hyper["lambda_x"]
lambda_theta = init_hyper["lambda_theta"]
eta = init_hyper["eta"]

dim1, dim2 = sparse_mat.shape
pos_train = np.where(sparse_mat != 0)
pos_test = np.where((dense_mat != 0) & (sparse_mat == 0))
binary_mat = sparse_mat.copy()
binary_mat[pos_train] = 1
d, rank = theta.shape

for it in range(maxiter):
## Update spatial matrix W
for i in range(dim1):
pos0 = np.where(sparse_mat[i, :] != 0)
Xt = X[pos0[0], :]
vec0 = np.matmul(Xt.T, sparse_mat[i, pos0[0]])
mat0 = inv(np.matmul(Xt.T, Xt) + lambda_w * np.eye(rank))
W[i, :] = np.matmul(mat0, vec0)
## Update temporal matrix X
for t in range(dim2):
pos0 = np.where(sparse_mat[:, t] != 0)
Wt = W[pos0[0], :]
Mt = np.zeros((rank, rank))
Nt = np.zeros(rank)
if t < np.max(time_lags):
Pt = np.zeros((rank, rank))
Qt = np.zeros(rank)
else:
Pt = np.eye(rank)
Qt = np.einsum('ij, ij -> j', theta, X[t - time_lags, :])
if t < dim2 - np.min(time_lags):
if t >= np.max(time_lags) and t < dim2 - np.max(time_lags):
index = list(range(0, d))
else:
index = list(np.where((t + time_lags >= np.max(time_lags)) & (t + time_lags < dim2)))[0]
for k in index:
Ak = theta[k, :]
Mt += np.diag(Ak ** 2)
theta0 = theta.copy()
theta0[k, :] = 0
Nt += np.multiply(Ak, X[t + time_lags[k], :]
- np.einsum('ij, ij -> j', theta0, X[t + time_lags[k] - time_lags, :]))
vec0 = np.matmul(Wt.T, sparse_mat[pos0[0], t]) + lambda_x * Nt + lambda_x * Qt
mat0 = inv(np.matmul(Wt.T, Wt) + lambda_x * Mt + lambda_x * Pt)
X[t, :] = np.matmul(mat0, vec0)
## Update AR coefficients theta
for k in range(d):
theta0 = theta.copy()
theta0[k, :] = 0
mat0 = np.zeros((dim2 - np.max(time_lags), rank))
for L in range(d):
mat0 += np.matmul(X[np.max(time_lags) - time_lags[L] : dim2 - time_lags[L] , :],
np.diag(theta0[L, :]))
VarPi = X[np.max(time_lags) : dim2, :] - mat0
var1 = np.zeros((rank, rank))
var2 = np.zeros(rank)
for t in range(np.max(time_lags), dim2):
B = X[t - time_lags[k], :]
var1 += np.diag(np.multiply(B, B))
var2 += np.matmul(np.diag(B), VarPi[t - np.max(time_lags), :])
theta[k, :] = np.matmul(inv(var1 + lambda_theta * np.eye(rank) / lambda_x), var2)

mat_hat = np.matmul(W, X.T)
mape = np.sum(np.abs(dense_mat[pos_test] - mat_hat[pos_test])
/ dense_mat[pos_test]) / dense_mat[pos_test].shape[0]
rmse = np.sqrt(np.sum((dense_mat[pos_test] - mat_hat[pos_test]) ** 2)/dense_mat[pos_test].shape[0])

if (it + 1) % 200 == 0:
print('Iter: {}'.format(it + 1))
print('Imputation MAPE: {:.6}'.format(mape))
print('Imputation RMSE: {:.6}'.format(rmse))
print()

## Missing Data Imputation¶

In the following, we apply the above defined TRMF function to the task of missing data imputation task on the following spatiotemporal multivariate time series datasets/matrices:

The original data sets have been adapted into our experiments, and it is now available at the fold of datasets.

### Experiments on Guangzhou Data Set¶

In [5]:
import scipy.io

tensor = tensor['tensor']
random_matrix = random_matrix['random_matrix']
random_tensor = random_tensor['random_tensor']

dense_mat = tensor.reshape([tensor.shape[0], tensor.shape[1] * tensor.shape[2]])
missing_rate = 0.2

# =============================================================================
### Random missing (RM) scenario
### Set the RM scenario by:
binary_mat = (np.round(random_tensor + 0.5 - missing_rate)
.reshape([random_tensor.shape[0], random_tensor.shape[1] * random_tensor.shape[2]]))
# =============================================================================

sparse_mat = np.multiply(dense_mat, binary_mat)
In [6]:
import time
start = time.time()
dim1, dim2 = sparse_mat.shape
rank = 80
time_lags = np.array([1, 2, 144])
d = time_lags.shape[0]
## Initialize parameters
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
init_para = {"W": W, "X": X, "theta": theta}
## Set hyparameters
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
init_hyper = {"lambda_w": lambda_w, "lambda_x": lambda_x, "lambda_theta": lambda_theta, "eta": eta}
maxiter = 1000
TRMF(dense_mat, sparse_mat, init_para, init_hyper, time_lags, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
Iter: 200
Imputation MAPE: 0.0748007
Imputation RMSE: 3.14413

Iter: 400
Imputation MAPE: 0.0747433
Imputation RMSE: 3.14283

Iter: 600
Imputation MAPE: 0.0747292
Imputation RMSE: 3.14258

Iter: 800
Imputation MAPE: 0.0747239
Imputation RMSE: 3.1425

Iter: 1000
Imputation MAPE: 0.0747206
Imputation RMSE: 3.14241

Running time: 5557 seconds
In [7]:
import scipy.io

tensor = tensor['tensor']
random_matrix = random_matrix['random_matrix']
random_tensor = random_tensor['random_tensor']

dense_mat = tensor.reshape([tensor.shape[0], tensor.shape[1] * tensor.shape[2]])
missing_rate = 0.4

# =============================================================================
### Random missing (RM) scenario
### Set the RM scenario by:
binary_mat = (np.round(random_tensor + 0.5 - missing_rate)
.reshape([random_tensor.shape[0], random_tensor.shape[1] * random_tensor.shape[2]]))
# =============================================================================

sparse_mat = np.multiply(dense_mat, binary_mat)
In [8]:
import time
start = time.time()
dim1, dim2 = sparse_mat.shape
rank = 80
time_lags = np.array([1, 2, 144])
d = time_lags.shape[0]
## Initialize parameters
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
init_para = {"W": W, "X": X, "theta": theta}
## Set hyparameters
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
init_hyper = {"lambda_w": lambda_w, "lambda_x": lambda_x, "lambda_theta": lambda_theta, "eta": eta}
maxiter = 1000
TRMF(dense_mat, sparse_mat, init_para, init_hyper, time_lags, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
Iter: 200
Imputation MAPE: 0.0776843
Imputation RMSE: 3.25726

Iter: 400
Imputation MAPE: 0.0775849
Imputation RMSE: 3.25442

Iter: 600
Imputation MAPE: 0.0775707
Imputation RMSE: 3.25406

Iter: 800
Imputation MAPE: 0.0775633
Imputation RMSE: 3.25383

Iter: 1000
Imputation MAPE: 0.0775572
Imputation RMSE: 3.25359

Running time: 5553 seconds
In [9]:
import scipy.io

tensor = tensor['tensor']
random_matrix = random_matrix['random_matrix']
random_tensor = random_tensor['random_tensor']

dense_mat = tensor.reshape([tensor.shape[0], tensor.shape[1] * tensor.shape[2]])
missing_rate = 0.2

# =============================================================================
### Non-random missing (NM) scenario
### Set the NM scenario by:
binary_tensor = np.zeros(tensor.shape)
for i1 in range(tensor.shape[0]):
for i2 in range(tensor.shape[1]):
binary_tensor[i1, i2, :] = np.round(random_matrix[i1, i2] + 0.5 - missing_rate)
binary_mat = binary_tensor.reshape([binary_tensor.shape[0], binary_tensor.shape[1] * binary_tensor.shape[2]])
# =============================================================================

sparse_mat = np.multiply(dense_mat, binary_mat)
In [10]:
import time
start = time.time()
dim1, dim2 = sparse_mat.shape
rank = 10
time_lags = np.array([1, 2, 144])
d = time_lags.shape[0]
## Initialize parameters
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
init_para = {"W": W, "X": X, "theta": theta}
## Set hyparameters
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
init_hyper = {"lambda_w": lambda_w, "lambda_x": lambda_x, "lambda_theta": lambda_theta, "eta": eta}
maxiter = 1000
TRMF(dense_mat, sparse_mat, init_para, init_hyper, time_lags, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
Iter: 200
Imputation MAPE: 0.102348
Imputation RMSE: 4.26877

Iter: 400
Imputation MAPE: 0.102366
Imputation RMSE: 4.27026

Iter: 600
Imputation MAPE: 0.102373
Imputation RMSE: 4.2707

Iter: 800
Imputation MAPE: 0.102376
Imputation RMSE: 4.27089

Iter: 1000
Imputation MAPE: 0.102378
Imputation RMSE: 4.27098

Running time: 1388 seconds
In [29]:
import scipy.io

tensor = tensor['tensor']
random_matrix = random_matrix['random_matrix']
random_tensor = random_tensor['random_tensor']

dense_mat = tensor.reshape([tensor.shape[0], tensor.shape[1] * tensor.shape[2]])
missing_rate = 0.4

# =============================================================================
### Non-random missing (NM) scenario
### Set the NM scenario by:
binary_tensor = np.zeros(tensor.shape)
for i1 in range(tensor.shape[0]):
for i2 in range(tensor.shape[1]):
binary_tensor[i1, i2, :] = np.round(random_matrix[i1, i2] + 0.5 - missing_rate)
binary_mat = binary_tensor.reshape([binary_tensor.shape[0], binary_tensor.shape[1] * binary_tensor.shape[2]])
# =============================================================================

sparse_mat = np.multiply(dense_mat, binary_mat)
In [30]:
import time
start = time.time()
dim1, dim2 = sparse_mat.shape
rank = 10
time_lags = np.array([1, 2, 144])
d = time_lags.shape[0]
## Initialize parameters
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
init_para = {"W": W, "X": X, "theta": theta}
## Set hyparameters
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
init_hyper = {"lambda_w": lambda_w, "lambda_x": lambda_x, "lambda_theta": lambda_theta, "eta": eta}
maxiter = 1000
TRMF(dense_mat, sparse_mat, init_para, init_hyper, time_lags, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
Iter: 200
Imputation MAPE: 0.103593
Imputation RMSE: 4.36593

Iter: 400
Imputation MAPE: 0.103651
Imputation RMSE: 4.36993

Iter: 600
Imputation MAPE: 0.103664
Imputation RMSE: 4.37084

Iter: 800
Imputation MAPE: 0.103669
Imputation RMSE: 4.37113

Iter: 1000
Imputation MAPE: 0.10367
Imputation RMSE: 4.37126

Running time: 1432 seconds

Experiment results of missing data imputation using TRMF:

scenario rank Lambda_w Lambda_x Lambda_theta eta maxiter mape rmse
20%, RM 80 500 500 500 0.03 1000 0.0747 3.1424
40%, RM 80 500 500 500 0.03 1000 0.0776 3.2536
20%, NM 10 500 500 500 0.03 1000 0.1024 4.2710
40%, NM 10 500 500 500 0.03 1000 0.1037 4.3713

### Experiments on Birmingham Data Set¶

In [13]:
import scipy.io

tensor = tensor['tensor']
random_matrix = random_matrix['random_matrix']
random_tensor = random_tensor['random_tensor']

dense_mat = tensor.reshape([tensor.shape[0], tensor.shape[1] * tensor.shape[2]])
missing_rate = 0.1

# =============================================================================
### Random missing (RM) scenario
### Set the RM scenario by:
binary_mat = (np.round(random_tensor + 0.5 - missing_rate)
.reshape([random_tensor.shape[0], random_tensor.shape[1] * random_tensor.shape[2]]))
# =============================================================================

sparse_mat = np.multiply(dense_mat, binary_mat)
In [14]:
import time
start = time.time()
dim1, dim2 = sparse_mat.shape
rank = 30
time_lags = np.array([1, 2, 18])
d = time_lags.shape[0]
## Initialize parameters
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
init_para = {"W": W, "X": X, "theta": theta}
## Set hyparameters
lambda_w = 100
lambda_x = 100
lambda_theta = 100
eta = 0.01
init_hyper = {"lambda_w": lambda_w, "lambda_x": lambda_x, "lambda_theta": lambda_theta, "eta": eta}
maxiter = 1000
TRMF(dense_mat, sparse_mat, init_para, init_hyper, time_lags, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
Iter: 200
Imputation MAPE: 0.0290578
Imputation RMSE: 11.8703

Iter: 400
Imputation MAPE: 0.0281009
Imputation RMSE: 11.433

Iter: 600
Imputation MAPE: 0.0284944
Imputation RMSE: 10.6691

Iter: 800
Imputation MAPE: 0.0277727
Imputation RMSE: 10.4855

Iter: 1000
Imputation MAPE: 0.0277411
Imputation RMSE: 10.5701

Running time: 287 seconds
In [15]:
import scipy.io

tensor = tensor['tensor']
random_matrix = random_matrix['random_matrix']
random_tensor = random_tensor['random_tensor']

dense_mat = tensor.reshape([tensor.shape[0], tensor.shape[1] * tensor.shape[2]])
missing_rate = 0.3

# =============================================================================
### Random missing (RM) scenario
### Set the RM scenario by:
binary_mat = (np.round(random_tensor + 0.5 - missing_rate)
.reshape([random_tensor.shape[0], random_tensor.shape[1] * random_tensor.shape[2]]))
# =============================================================================

sparse_mat = np.multiply(dense_mat, binary_mat)
In [16]:
import time
start = time.time()
dim1, dim2 = sparse_mat.shape
rank = 30
time_lags = np.array([1, 2, 18])
d = time_lags.shape[0]
## Initialize parameters
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
init_para = {"W": W, "X": X, "theta": theta}
## Set hyparameters
lambda_w = 100
lambda_x = 100
lambda_theta = 100
eta = 0.01
init_hyper = {"lambda_w": lambda_w, "lambda_x": lambda_x, "lambda_theta": lambda_theta, "eta": eta}
maxiter = 1000
TRMF(dense_mat, sparse_mat, init_para, init_hyper, time_lags, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
Iter: 200
Imputation MAPE: 0.0384997
Imputation RMSE: 24.408

Iter: 400
Imputation MAPE: 0.0365896
Imputation RMSE: 21.9966

Iter: 600
Imputation MAPE: 0.036827
Imputation RMSE: 21.8252

Iter: 800
Imputation MAPE: 0.0369684
Imputation RMSE: 21.8012

Iter: 1000
Imputation MAPE: 0.0369032
Imputation RMSE: 21.8022

Running time: 293 seconds
In [17]:
import scipy.io

tensor = tensor['tensor']
random_matrix = random_matrix['random_matrix']
random_tensor = random_tensor['random_tensor']

dense_mat = tensor.reshape([tensor.shape[0], tensor.shape[1] * tensor.shape[2]])
missing_rate = 0.1

# =============================================================================
### Non-random missing (NM) scenario
### Set the NM scenario by:
binary_tensor = np.zeros(tensor.shape)
for i1 in range(tensor.shape[0]):
for i2 in range(tensor.shape[1]):
binary_tensor[i1, i2, :] = np.round(random_matrix[i1, i2] + 0.5 - missing_rate)
binary_mat = binary_tensor.reshape([binary_tensor.shape[0], binary_tensor.shape[1] * binary_tensor.shape[2]])
# =============================================================================

sparse_mat = np.multiply(dense_mat, binary_mat)
In [18]:
import time
start = time.time()
dim1, dim2 = sparse_mat.shape
rank = 10
time_lags = np.array([1, 2, 18])
d = time_lags.shape[0]
## Initialize parameters
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
init_para = {"W": W, "X": X, "theta": theta}
## Set hyparameters
lambda_w = 100
lambda_x = 100
lambda_theta = 100
eta = 0.01
init_hyper = {"lambda_w": lambda_w, "lambda_x": lambda_x, "lambda_theta": lambda_theta, "eta": eta}
maxiter = 1000
TRMF(dense_mat, sparse_mat, init_para, init_hyper, time_lags, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
Iter: 200
Imputation MAPE: 0.128079
Imputation RMSE: 29.9688

Iter: 400
Imputation MAPE: 0.12757
Imputation RMSE: 29.6351

Iter: 600
Imputation MAPE: 0.127128
Imputation RMSE: 29.4746

Iter: 800
Imputation MAPE: 0.127255
Imputation RMSE: 29.3843

Iter: 1000
Imputation MAPE: 0.127435
Imputation RMSE: 29.4629

Running time: 205 seconds
In [19]:
import scipy.io

tensor = tensor['tensor']
random_matrix = random_matrix['random_matrix']
random_tensor = random_tensor['random_tensor']

dense_mat = tensor.reshape([tensor.shape[0], tensor.shape[1] * tensor.shape[2]])
missing_rate = 0.3

# =============================================================================
### Non-random missing (NM) scenario
### Set the NM scenario by:
binary_tensor = np.zeros(tensor.shape)
for i1 in range(tensor.shape[0]):
for i2 in range(tensor.shape[1]):
binary_tensor[i1, i2, :] = np.round(random_matrix[i1, i2] + 0.5 - missing_rate)
binary_mat = binary_tensor.reshape([binary_tensor.shape[0], binary_tensor.shape[1] * binary_tensor.shape[2]])
# =============================================================================

sparse_mat = np.multiply(dense_mat, binary_mat)
In [20]:
import time
start = time.time()
dim1, dim2 = sparse_mat.shape
rank = 10
time_lags = np.array([1, 2, 18])
d = time_lags.shape[0]
## Initialize parameters
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
init_para = {"W": W, "X": X, "theta": theta}
## Set hyparameters
lambda_w = 100
lambda_x = 100
lambda_theta = 100
eta = 0.01
init_hyper = {"lambda_w": lambda_w, "lambda_x": lambda_x, "lambda_theta": lambda_theta, "eta": eta}
maxiter = 1000
TRMF(dense_mat, sparse_mat, init_para, init_hyper, time_lags, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
Iter: 200
Imputation MAPE: 0.191454
Imputation RMSE: 118.747

Iter: 400
Imputation MAPE: 0.162524
Imputation RMSE: 81.9807

Iter: 600
Imputation MAPE: 0.160793
Imputation RMSE: 81.9135

Iter: 800
Imputation MAPE: 0.157669
Imputation RMSE: 82.2151

Iter: 1000
Imputation MAPE: 0.163484
Imputation RMSE: 85.9752

Running time: 207 seconds

Experiment results of missing data imputation using TRMF:

scenario rank Lambda_w Lambda_x Lambda_theta eta maxiter mape rmse
10%, RM 30 100 100 100 0.01 1000 0.0277 10.5701
30%, RM 30 100 100 100 0.01 1000 0.0369 21.8022
10%, NM 10 100 100 100 0.01 1000 0.1274 29.4629
30%, NM 10 100 100 100 0.01 1000 0.1635 85.9752

### Experiments on Hangzhou Data Set¶

In [21]:
import scipy.io

tensor = tensor['tensor']
random_matrix = random_matrix['random_matrix']
random_tensor = random_tensor['random_tensor']

dense_mat = tensor.reshape([tensor.shape[0], tensor.shape[1] * tensor.shape[2]])
missing_rate = 0.2

# =============================================================================
### Random missing (RM) scenario
### Set the RM scenario by:
binary_mat = (np.round(random_tensor + 0.5 - missing_rate)
.reshape([random_tensor.shape[0], random_tensor.shape[1] * random_tensor.shape[2]]))
# =============================================================================

sparse_mat = np.multiply(dense_mat, binary_mat)
In [22]:
import time
start = time.time()
dim1, dim2 = sparse_mat.shape
rank = 50
time_lags = np.array([1, 2, 108])
d = time_lags.shape[0]
## Initialize parameters
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
init_para = {"W": W, "X": X, "theta": theta}
## Set hyparameters
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
init_hyper = {"lambda_w": lambda_w, "lambda_x": lambda_x, "lambda_theta": lambda_theta, "eta": eta}
maxiter = 1000
TRMF(dense_mat, sparse_mat, init_para, init_hyper, time_lags, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
Iter: 200
Imputation MAPE: 0.214013
Imputation RMSE: 35.2672

Iter: 400
Imputation MAPE: 0.214398
Imputation RMSE: 36.3985

Iter: 600
Imputation MAPE: 0.214302
Imputation RMSE: 36.8369

Iter: 800
Imputation MAPE: 0.213989
Imputation RMSE: 36.9683

Iter: 1000
Imputation MAPE: 0.213147
Imputation RMSE: 37.0673

Running time: 937 seconds
In [23]:
import scipy.io

tensor = tensor['tensor']
random_matrix = random_matrix['random_matrix']
random_tensor = random_tensor['random_tensor']

dense_mat = tensor.reshape([tensor.shape[0], tensor.shape[1] * tensor.shape[2]])
missing_rate = 0.4

# =============================================================================
### Random missing (RM) scenario
### Set the RM scenario by:
binary_mat = (np.round(random_tensor + 0.5 - missing_rate)
.reshape([random_tensor.shape[0], random_tensor.shape[1] * random_tensor.shape[2]]))
# =============================================================================

sparse_mat = np.multiply(dense_mat, binary_mat)
In [24]:
import time
start = time.time()
dim1, dim2 = sparse_mat.shape
rank = 50
time_lags = np.array([1, 2, 108])
d = time_lags.shape[0]
## Initialize parameters
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
init_para = {"W": W, "X": X, "theta": theta}
## Set hyparameters
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
init_hyper = {"lambda_w": lambda_w, "lambda_x": lambda_x, "lambda_theta": lambda_theta, "eta": eta}
maxiter = 1000
TRMF(dense_mat, sparse_mat, init_para, init_hyper, time_lags, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
Iter: 200
Imputation MAPE: 0.230955
Imputation RMSE: 36.7006

Iter: 400
Imputation MAPE: 0.230059
Imputation RMSE: 37.6722

Iter: 600
Imputation MAPE: 0.22898
Imputation RMSE: 38.0138

Iter: 800
Imputation MAPE: 0.228664
Imputation RMSE: 38.1187

Iter: 1000
Imputation MAPE: 0.228853
Imputation RMSE: 38.15

Running time: 940 seconds
In [25]:
import scipy.io

tensor = tensor['tensor']
random_matrix = random_matrix['random_matrix']
random_tensor = random_tensor['random_tensor']

dense_mat = tensor.reshape([tensor.shape[0], tensor.shape[1] * tensor.shape[2]])
missing_rate = 0.2

# =============================================================================
### Non-random missing (NM) scenario
### Set the NM scenario by:
binary_tensor = np.zeros(tensor.shape)
for i1 in range(tensor.shape[0]):
for i2 in range(tensor.shape[1]):
binary_tensor[i1, i2, :] = np.round(random_matrix[i1, i2] + 0.5 - missing_rate)
binary_mat = binary_tensor.reshape([binary_tensor.shape[0], binary_tensor.shape[1] * binary_tensor.shape[2]])
# =============================================================================

sparse_mat = np.multiply(dense_mat, binary_mat)
In [26]:
import time
start = time.time()
dim1, dim2 = sparse_mat.shape
rank = 10
time_lags = np.array([1, 2, 108])
d = time_lags.shape[0]
## Initialize parameters
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
init_para = {"W": W, "X": X, "theta": theta}
## Set hyparameters
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
init_hyper = {"lambda_w": lambda_w, "lambda_x": lambda_x, "lambda_theta": lambda_theta, "eta": eta}
maxiter = 1000
TRMF(dense_mat, sparse_mat, init_para, init_hyper, time_lags, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
Iter: 200
Imputation MAPE: 0.255072
Imputation RMSE: 40.6638

Iter: 400
Imputation MAPE: 0.258543
Imputation RMSE: 38.1323

Iter: 600
Imputation MAPE: 0.260439
Imputation RMSE: 39.8118

Iter: 800
Imputation MAPE: 0.260712
Imputation RMSE: 40.0304

Iter: 1000
Imputation MAPE: 0.260727
Imputation RMSE: 40.0598

Running time: 405 seconds
In [27]:
import scipy.io

tensor = tensor['tensor']
random_matrix = random_matrix['random_matrix']
random_tensor = random_tensor['random_tensor']

dense_mat = tensor.reshape([tensor.shape[0], tensor.shape[1] * tensor.shape[2]])
missing_rate = 0.4

# =============================================================================
### Non-random missing (NM) scenario
### Set the NM scenario by:
binary_tensor = np.zeros(tensor.shape)
for i1 in range(tensor.shape[0]):
for i2 in range(tensor.shape[1]):
binary_tensor[i1, i2, :] = np.round(random_matrix[i1, i2] + 0.5 - missing_rate)
binary_mat = binary_tensor.reshape([binary_tensor.shape[0], binary_tensor.shape[1] * binary_tensor.shape[2]])
# =============================================================================

sparse_mat = np.multiply(dense_mat, binary_mat)
In [28]:
import time
start = time.time()
dim1, dim2 = sparse_mat.shape
rank = 10
time_lags = np.array([1, 2, 108])
d = time_lags.shape[0]
## Initialize parameters
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
init_para = {"W": W, "X": X, "theta": theta}
## Set hyparameters
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
init_hyper = {"lambda_w": lambda_w, "lambda_x": lambda_x, "lambda_theta": lambda_theta, "eta": eta}
maxiter = 1000
TRMF(dense_mat, sparse_mat, init_para, init_hyper, time_lags, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
Iter: 200
Imputation MAPE: 0.275611
Imputation RMSE: 40.7617

Iter: 400
Imputation MAPE: 0.276333
Imputation RMSE: 42.7749

Iter: 600
Imputation MAPE: 0.273609
Imputation RMSE: 40.2961

Iter: 800
Imputation MAPE: 0.273325
Imputation RMSE: 39.7922

Iter: 1000
Imputation MAPE: 0.273192
Imputation RMSE: 39.7538

Running time: 406 seconds

Experiment results of missing data imputation using TRMF:

scenario rank Lambda_w Lambda_x Lambda_theta eta maxiter mape rmse
20%, RM 50 1000 1000 1000 0.03 1000 0.2131 37.0673
40%, RM 50 1000 1000 1000 0.03 1000 0.2289 38.15
20%, NM 10 1000 500 500 0.03 1000 0.2607 40.0598
40%, NM 10 1000 500 500 0.03 1000 0.2732 39.7538

### Experiments on Seattle Data Set¶

In [5]:
import pandas as pd

dense_mat = pd.read_csv('../datasets/Seattle-data-set/mat.csv', index_col = 0)
RM_mat = pd.read_csv('../datasets/Seattle-data-set/RM_mat.csv', index_col = 0)
dense_mat = dense_mat.values
RM_mat = RM_mat.values

missing_rate = 0.2

# =============================================================================
### Random missing (RM) scenario
### Set the RM scenario by:
binary_mat = np.round(RM_mat + 0.5 - missing_rate)
# =============================================================================

sparse_mat = np.multiply(dense_mat, binary_mat)
In [6]:
import time
start = time.time()
dim1, dim2 = sparse_mat.shape
rank = 50
time_lags = np.array([1, 2, 288])
d = time_lags.shape[0]
## Initialize parameters
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
init_para = {"W": W, "X": X, "theta": theta}
## Set hyparameters
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
init_hyper = {"lambda_w": lambda_w, "lambda_x": lambda_x, "lambda_theta": lambda_theta, "eta": eta}
maxiter = 1000
TRMF(dense_mat, sparse_mat, init_para, init_hyper, time_lags, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
Iter: 200
Imputation MAPE: 0.0596004
Imputation RMSE: 3.71518

Iter: 400
Imputation MAPE: 0.0595838
Imputation RMSE: 3.7148

Iter: 600
Imputation MAPE: 0.0595782
Imputation RMSE: 3.71476

Iter: 800
Imputation MAPE: 0.059576
Imputation RMSE: 3.71476

Iter: 1000
Imputation MAPE: 0.059575
Imputation RMSE: 3.71477

Running time: 3157 seconds
In [7]:
import pandas as pd

dense_mat = pd.read_csv('../datasets/Seattle-data-set/mat.csv', index_col = 0)
RM_mat = pd.read_csv('../datasets/Seattle-data-set/RM_mat.csv', index_col = 0)
dense_mat = dense_mat.values
RM_mat = RM_mat.values

missing_rate = 0.4

# =============================================================================
### Random missing (RM) scenario
### Set the RM scenario by:
binary_mat = np.round(RM_mat + 0.5 - missing_rate)
# =============================================================================

sparse_mat = np.multiply(dense_mat, binary_mat)
In [8]:
import time
start = time.time()
dim1, dim2 = sparse_mat.shape
rank = 50
time_lags = np.array([1, 2, 288])
d = time_lags.shape[0]
## Initialize parameters
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
init_para = {"W": W, "X": X, "theta": theta}
## Set hyparameters
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
init_hyper = {"lambda_w": lambda_w, "lambda_x": lambda_x, "lambda_theta": lambda_theta, "eta": eta}
maxiter = 1000
TRMF(dense_mat, sparse_mat, init_para, init_hyper, time_lags, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
Iter: 200
Imputation MAPE: 0.0616689
Imputation RMSE: 3.79457

Iter: 400
Imputation MAPE: 0.0616347
Imputation RMSE: 3.79348

Iter: 600
Imputation MAPE: 0.0616226
Imputation RMSE: 3.79312

Iter: 800
Imputation MAPE: 0.0616164
Imputation RMSE: 3.79291

Iter: 1000
Imputation MAPE: 0.0616129
Imputation RMSE: 3.79277

Running time: 3098 seconds
In [7]:
import pandas as pd

dense_mat = pd.read_csv('../datasets/Seattle-data-set/mat.csv', index_col = 0)
NM_mat = pd.read_csv('../datasets/Seattle-data-set/NM_mat.csv', index_col = 0)
dense_mat = dense_mat.values
NM_mat = NM_mat.values

missing_rate = 0.2

# =============================================================================
### Non-random missing (NM) scenario
### Set the NM scenario by:
binary_tensor = np.zeros((dense_mat.shape[0], 28, 288))
for i1 in range(binary_tensor.shape[0]):
for i2 in range(binary_tensor.shape[1]):
binary_tensor[i1, i2, :] = np.round(NM_mat[i1, i2] + 0.5 - missing_rate)
# =============================================================================

sparse_mat = np.multiply(dense_mat, binary_tensor.reshape([dense_mat.shape[0], dense_mat.shape[1]]))
In [8]:
import time
start = time.time()
dim1, dim2 = sparse_mat.shape
rank = 10
time_lags = np.array([1, 2, 288])
d = time_lags.shape[0]
## Initialize parameters
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
init_para = {"W": W, "X": X, "theta": theta}
## Set hyparameters
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
init_hyper = {"lambda_w": lambda_w, "lambda_x": lambda_x, "lambda_theta": lambda_theta, "eta": eta}
maxiter = 1000
TRMF(dense_mat, sparse_mat, init_para, init_hyper, time_lags, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
Iter: 200
Imputation MAPE: 0.0912403
Imputation RMSE: 5.26161

Iter: 400
Imputation MAPE: 0.0912302
Imputation RMSE: 5.26205

Iter: 600
Imputation MAPE: 0.0912316
Imputation RMSE: 5.26238

Iter: 800
Imputation MAPE: 0.0912326
Imputation RMSE: 5.26255

Iter: 1000
Imputation MAPE: 0.0912332
Imputation RMSE: 5.26264

Running time: 1370 seconds
In [9]:
import pandas as pd

dense_mat = pd.read_csv('../datasets/Seattle-data-set/mat.csv', index_col = 0)
NM_mat = pd.read_csv('../datasets/Seattle-data-set/NM_mat.csv', index_col = 0)
dense_mat = dense_mat.values
NM_mat = NM_mat.values

missing_rate = 0.4

# =============================================================================
### Non-random missing (NM) scenario
### Set the NM scenario by:
binary_tensor = np.zeros((dense_mat.shape[0], 28, 288))
for i1 in range(binary_tensor.shape[0]):
for i2 in range(binary_tensor.shape[1]):
binary_tensor[i1, i2, :] = np.round(NM_mat[i1, i2] + 0.5 - missing_rate)
# =============================================================================

sparse_mat = np.multiply(dense_mat, binary_tensor.reshape([dense_mat.shape[0], dense_mat.shape[1]]))
In [10]:
import time
start = time.time()
dim1, dim2 = sparse_mat.shape
rank = 10
time_lags = np.array([1, 2, 288])
d = time_lags.shape[0]
## Initialize parameters
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
init_para = {"W": W, "X": X, "theta": theta}
## Set hyparameters
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
init_hyper = {"lambda_w": lambda_w, "lambda_x": lambda_x, "lambda_theta": lambda_theta, "eta": eta}
maxiter = 1000
TRMF(dense_mat, sparse_mat, init_para, init_hyper, time_lags, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
Iter: 200
Imputation MAPE: 0.0919406
Imputation RMSE: 5.29846

Iter: 400
Imputation MAPE: 0.0919397
Imputation RMSE: 5.2991

Iter: 600
Imputation MAPE: 0.0919385
Imputation RMSE: 5.29936

Iter: 800
Imputation MAPE: 0.0919378
Imputation RMSE: 5.29948

Iter: 1000
Imputation MAPE: 0.0919373
Imputation RMSE: 5.29954

Running time: 1414 seconds

Experiment results of missing data imputation using TRMF:

scenario rank Lambda_w Lambda_x Lambda_theta eta maxiter mape rmse
20%, RM 50 1000 1000 1000 0.03 1000 0.0596 3.7148
40%, RM 50 1000 1000 1000 0.03 1000 0.0616 3.7928
20%, NM 10 1000 500 500 0.03 1000 0.0912 5.2626
40%, NM 10 1000 500 500 0.03 1000 0.0919 5.2995