Reference: 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

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}.$$

Temporal Regularized Matrix Factorization(TRMF)

Temporal Regularized Matrix Factorization (TRMF) framework is an approach to incorporate temporal dependencies into matrix factorization models which use well-studied time series models to describe temporal dependencies among ${\boldsymbol{x}_t}$ explicitly.Such models take 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.

Matrix Computation Concepts

Kronecker product

  • Definition:

Given two matrices $A\in\mathbb{R}^{m_1\times n_1}$ and $B\in\mathbb{R}^{m_2\times n_2}$, then, the Kronecker product between these two matrices is defined as

$$A\otimes B=\left[ \begin{array}{cccc} a_{11}B & a_{12}B & \cdots & a_{1m_2}B \\ a_{21}B & a_{22}B & \cdots & a_{2m_2}B \\ \vdots & \vdots & \ddots & \vdots \\ a_{m_11}B & a_{m_12}B & \cdots & a_{m_1m_2}B \\ \end{array} \right]$$

where the symbol $\otimes$ denotes Kronecker product, and the size of resulted $A\otimes B$ is $(m_1m_2)\times (n_1n_2)$ (i.e., $m_1\times m_2$ columns and $n_1\times n_2$ rows).

  • Example:

If $A=\left[ \begin{array}{cc} 1 & 2 \\ 3 & 4 \\ \end{array} \right]$ and $B=\left[ \begin{array}{ccc} 5 & 6 & 7\\ 8 & 9 & 10 \\ \end{array} \right]$, then, we have

$$A\otimes B=\left[ \begin{array}{cc} 1\times \left[ \begin{array}{ccc} 5 & 6 & 7\\ 8 & 9 & 10\\ \end{array} \right] & 2\times \left[ \begin{array}{ccc} 5 & 6 & 7\\ 8 & 9 & 10\\ \end{array} \right] \\ 3\times \left[ \begin{array}{ccc} 5 & 6 & 7\\ 8 & 9 & 10\\ \end{array} \right] & 4\times \left[ \begin{array}{ccc} 5 & 6 & 7\\ 8 & 9 & 10\\ \end{array} \right] \\ \end{array} \right]$$$$=\left[ \begin{array}{cccccc} 5 & 6 & 7 & 10 & 12 & 14 \\ 8 & 9 & 10 & 16 & 18 & 20 \\ 15 & 18 & 21 & 20 & 24 & 28 \\ 24 & 27 & 30 & 32 & 36 & 40 \\ \end{array} \right]\in\mathbb{R}^{4\times 6}.$$

Khatri-Rao product (kr_prod)

  • Definition:

Given two matrices $A=\left( \boldsymbol{a}_1,\boldsymbol{a}_2,...,\boldsymbol{a}_r \right)\in\mathbb{R}^{m\times r}$ and $B=\left( \boldsymbol{b}_1,\boldsymbol{b}_2,...,\boldsymbol{b}_r \right)\in\mathbb{R}^{n\times r}$ with same number of columns, then, the Khatri-Rao product (or column-wise Kronecker product) between $A$ and $B$ is given as follows,

$$A\odot B=\left( \boldsymbol{a}_1\otimes \boldsymbol{b}_1,\boldsymbol{a}_2\otimes \boldsymbol{b}_2,...,\boldsymbol{a}_r\otimes \boldsymbol{b}_r \right)\in\mathbb{R}^{(mn)\times r}$$

where the symbol $\odot$ denotes Khatri-Rao product, and $\otimes$ denotes Kronecker product.

  • Example:

If $A=\left[ \begin{array}{cc} 1 & 2 \\ 3 & 4 \\ \end{array} \right]=\left( \boldsymbol{a}_1,\boldsymbol{a}_2 \right) $ and $B=\left[ \begin{array}{cc} 5 & 6 \\ 7 & 8 \\ 9 & 10 \\ \end{array} \right]=\left( \boldsymbol{b}_1,\boldsymbol{b}_2 \right) $, then, we have

$$A\odot B=\left( \boldsymbol{a}_1\otimes \boldsymbol{b}_1,\boldsymbol{a}_2\otimes \boldsymbol{b}_2 \right) $$$$=\left[ \begin{array}{cc} \left[ \begin{array}{c} 1 \\ 3 \\ \end{array} \right]\otimes \left[ \begin{array}{c} 5 \\ 7 \\ 9 \\ \end{array} \right] & \left[ \begin{array}{c} 2 \\ 4 \\ \end{array} \right]\otimes \left[ \begin{array}{c} 6 \\ 8 \\ 10 \\ \end{array} \right] \\ \end{array} \right]$$$$=\left[ \begin{array}{cc} 5 & 12 \\ 7 & 16 \\ 9 & 20 \\ 15 & 24 \\ 21 & 32 \\ 27 & 40 \\ \end{array} \right]\in\mathbb{R}^{6\times 2}.$$
In [1]:
import numpy as np
from numpy.linalg import inv as inv
In [2]:
def kr_prod(a, b):
    return np.einsum('ir, jr -> ijr', a, b).reshape(a.shape[0] * b.shape[0], -1)
In [3]:
import numpy as np
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8], [9, 10]])
print(kr_prod(A, B))
[[ 5 12]
 [ 7 16]
 [ 9 20]
 [15 24]
 [21 32]
 [27 40]]
In [4]:
def TRMF(dense_mat, sparse_mat, W, X, theta, time_lags, lambda_w, lambda_x, lambda_theta, eta, maxiter):
    dim1, dim2 = sparse_mat.shape
    binary_mat = np.zeros((dim1, dim2))
    position = np.where((sparse_mat != 0))
    binary_mat[position] = 1
    pos = np.where((dense_mat != 0) & (sparse_mat == 0))
    d, r = theta.shape

    for iters in range(maxiter):
        var1 = X.T
        var2 = kr_prod(var1, var1)
        var3 = np.matmul(var2, binary_mat.T).reshape([r, r, dim1]) + np.dstack([lambda_w * np.eye(r)] * dim1)
        var4 = np.matmul(var1, sparse_mat.T)
        for i in range(dim1):
            W[i, :] = np.matmul(inv(var3[:, :, i]), var4[:, i])

        var1 = W.T
        var2 = kr_prod(var1, var1)
        var3 = np.matmul(var2, binary_mat).reshape([r, r, dim2]) + np.dstack([lambda_x * eta * np.eye(r)] * dim2)
        var4 = np.matmul(var1, sparse_mat)
        for t in range(dim2):
            Mt = np.zeros((r, r))
            Nt = np.zeros(r)
            if t < np.max(time_lags):
                Pt = np.zeros((r, r))
                Qt = np.zeros(r)
            else:
                Pt = np.eye(r)
                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
                    var5 = (X[t + time_lags[k], :] 
                            - np.einsum('ij, ij -> j', theta0, X[t + time_lags[k] - time_lags, :]))
                    Nt += np.multiply(Ak, var5)
            var_mu = var4[:, t] + lambda_x * Nt + lambda_x * Qt
            var_Lambda = var3[:, :, t] + lambda_x * Mt + lambda_x * Pt
            X[t, :] = np.matmul(inv(var_Lambda), var_mu)
        
        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(r)/lambda_x), var2)

        mat_hat = np.matmul(W, X.T)
        mape = np.sum(np.abs(dense_mat[pos] - mat_hat[pos])/dense_mat[pos])/dense_mat[pos].shape[0]
        rmse = np.sqrt(np.sum((dense_mat[pos] - mat_hat[pos]) ** 2)/dense_mat[pos].shape[0])
        
        if (iters + 1) % 200 == 0:
            print('Iter: {}'.format(iters + 1))
            print('Imputation MAPE: {:.6}'.format(mape))
            print('Imputation RMSE: {:.6}'.format(rmse))
            print()
In [5]:
import scipy.io

tensor = scipy.io.loadmat('../datasets/Guangzhou-data-set/tensor.mat')
tensor = tensor['tensor']
random_matrix = scipy.io.loadmat('../datasets/Guangzhou-data-set/random_matrix.mat')
random_matrix = random_matrix['random_matrix']
random_tensor = scipy.io.loadmat('../datasets/Guangzhou-data-set/random_tensor.mat')
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
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
time_lags = np.array([1, 2, 144])
d = time_lags.shape[0]
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
maxiter = 1000
TRMF(dense_mat, sparse_mat, W, X, theta, time_lags, lambda_w, lambda_x, lambda_theta, eta, 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 = scipy.io.loadmat('../datasets/Guangzhou-data-set/tensor.mat')
tensor = tensor['tensor']
random_matrix = scipy.io.loadmat('../datasets/Guangzhou-data-set/random_matrix.mat')
random_matrix = random_matrix['random_matrix']
random_tensor = scipy.io.loadmat('../datasets/Guangzhou-data-set/random_tensor.mat')
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
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
time_lags = np.array([1, 2, 144])
d = time_lags.shape[0]
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
maxiter = 1000
TRMF(dense_mat, sparse_mat, W, X, theta, time_lags, lambda_w, lambda_x, lambda_theta, eta, 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 = scipy.io.loadmat('../datasets/Guangzhou-data-set/tensor.mat')
tensor = tensor['tensor']
random_matrix = scipy.io.loadmat('../datasets/Guangzhou-data-set/random_matrix.mat')
random_matrix = random_matrix['random_matrix']
random_tensor = scipy.io.loadmat('../datasets/Guangzhou-data-set/random_tensor.mat')
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
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
time_lags = np.array([1, 2, 144])
d = time_lags.shape[0]
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
maxiter = 1000
TRMF(dense_mat, sparse_mat, W, X, theta, time_lags, lambda_w, lambda_x, lambda_theta, eta, 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 = scipy.io.loadmat('../datasets/Guangzhou-data-set/tensor.mat')
tensor = tensor['tensor']
random_matrix = scipy.io.loadmat('../datasets/Guangzhou-data-set/random_matrix.mat')
random_matrix = random_matrix['random_matrix']
random_tensor = scipy.io.loadmat('../datasets/Guangzhou-data-set/random_tensor.mat')
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
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
time_lags = np.array([1, 2, 144])
d = time_lags.shape[0]
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
maxiter = 1000
TRMF(dense_mat, sparse_mat, W, X, theta, time_lags, lambda_w, lambda_x, lambda_theta, eta, 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
In [13]:
import scipy.io

tensor = scipy.io.loadmat('../datasets/Birmingham-data-set/tensor.mat')
tensor = tensor['tensor']
random_matrix = scipy.io.loadmat('../datasets/Birmingham-data-set/random_matrix.mat')
random_matrix = random_matrix['random_matrix']
random_tensor = scipy.io.loadmat('../datasets/Birmingham-data-set/random_tensor.mat')
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
lambda_w = 100
lambda_x = 100
lambda_theta = 100
eta = 0.01
time_lags = np.array([1, 2, 18])
d = time_lags.shape[0]
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
maxiter = 1000
TRMF(dense_mat, sparse_mat, W, X, theta, time_lags, lambda_w, lambda_x, lambda_theta, eta, 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 = scipy.io.loadmat('../datasets/Birmingham-data-set/tensor.mat')
tensor = tensor['tensor']
random_matrix = scipy.io.loadmat('../datasets/Birmingham-data-set/random_matrix.mat')
random_matrix = random_matrix['random_matrix']
random_tensor = scipy.io.loadmat('../datasets/Birmingham-data-set/random_tensor.mat')
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
lambda_w = 100
lambda_x = 100
lambda_theta = 100
eta = 0.01
time_lags = np.array([1, 2, 18])
d = time_lags.shape[0]
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
maxiter = 1000
TRMF(dense_mat, sparse_mat, W, X, theta, time_lags, lambda_w, lambda_x, lambda_theta, eta, 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 = scipy.io.loadmat('../datasets/Birmingham-data-set/tensor.mat')
tensor = tensor['tensor']
random_matrix = scipy.io.loadmat('../datasets/Birmingham-data-set/random_matrix.mat')
random_matrix = random_matrix['random_matrix']
random_tensor = scipy.io.loadmat('../datasets/Birmingham-data-set/random_tensor.mat')
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
lambda_w = 100
lambda_x = 100
lambda_theta = 100
eta = 0.01
time_lags = np.array([1, 2, 18])
d = time_lags.shape[0]
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
maxiter = 1000
TRMF(dense_mat, sparse_mat, W, X, theta, time_lags, lambda_w, lambda_x, lambda_theta, eta, 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 = scipy.io.loadmat('../datasets/Birmingham-data-set/tensor.mat')
tensor = tensor['tensor']
random_matrix = scipy.io.loadmat('../datasets/Birmingham-data-set/random_matrix.mat')
random_matrix = random_matrix['random_matrix']
random_tensor = scipy.io.loadmat('../datasets/Birmingham-data-set/random_tensor.mat')
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
lambda_w = 100
lambda_x = 100
lambda_theta = 100
eta = 0.01
time_lags = np.array([1, 2, 18])
d = time_lags.shape[0]
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
maxiter = 1000
TRMF(dense_mat, sparse_mat, W, X, theta, time_lags, lambda_w, lambda_x, lambda_theta, eta, 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
In [21]:
import scipy.io

tensor = scipy.io.loadmat('../datasets/Hangzhou-data-set/tensor.mat')
tensor = tensor['tensor']
random_matrix = scipy.io.loadmat('../datasets/Hangzhou-data-set/random_matrix.mat')
random_matrix = random_matrix['random_matrix']
random_tensor = scipy.io.loadmat('../datasets/Hangzhou-data-set/random_tensor.mat')
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
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
time_lags = np.array([1, 2, 108])
d = time_lags.shape[0]
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
maxiter = 1000
TRMF(dense_mat, sparse_mat, W, X, theta, time_lags, lambda_w, lambda_x, lambda_theta, eta, 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 = scipy.io.loadmat('../datasets/Hangzhou-data-set/tensor.mat')
tensor = tensor['tensor']
random_matrix = scipy.io.loadmat('../datasets/Hangzhou-data-set/random_matrix.mat')
random_matrix = random_matrix['random_matrix']
random_tensor = scipy.io.loadmat('../datasets/Hangzhou-data-set/random_tensor.mat')
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
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
time_lags = np.array([1, 2, 108])
d = time_lags.shape[0]
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
maxiter = 1000
TRMF(dense_mat, sparse_mat, W, X, theta, time_lags, lambda_w, lambda_x, lambda_theta, eta, 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 = scipy.io.loadmat('../datasets/Hangzhou-data-set/tensor.mat')
tensor = tensor['tensor']
random_matrix = scipy.io.loadmat('../datasets/Hangzhou-data-set/random_matrix.mat')
random_matrix = random_matrix['random_matrix']
random_tensor = scipy.io.loadmat('../datasets/Hangzhou-data-set/random_tensor.mat')
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
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
time_lags = np.array([1, 2, 108])
d = time_lags.shape[0]
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
maxiter = 1000
TRMF(dense_mat, sparse_mat, W, X, theta, time_lags, lambda_w, lambda_x, lambda_theta, eta, 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 = scipy.io.loadmat('../datasets/Hangzhou-data-set/tensor.mat')
tensor = tensor['tensor']
random_matrix = scipy.io.loadmat('../datasets/Hangzhou-data-set/random_matrix.mat')
random_matrix = random_matrix['random_matrix']
random_tensor = scipy.io.loadmat('../datasets/Hangzhou-data-set/random_tensor.mat')
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
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
time_lags = np.array([1, 2, 108])
d = time_lags.shape[0]
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
maxiter = 1000
TRMF(dense_mat, sparse_mat, W, X, theta, time_lags, lambda_w, lambda_x, lambda_theta, eta, 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
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
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
time_lags = np.array([1, 2, 288])
d = time_lags.shape[0]
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
maxiter = 1000
TRMF(dense_mat, sparse_mat, W, X, theta, time_lags, lambda_w, lambda_x, lambda_theta, eta, 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
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
time_lags = np.array([1, 2, 288])
d = time_lags.shape[0]
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
maxiter = 1000
TRMF(dense_mat, sparse_mat, W, X, theta, time_lags, lambda_w, lambda_x, lambda_theta, eta, 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
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
time_lags = np.array([1, 2, 288])
d = time_lags.shape[0]
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
maxiter = 1000
TRMF(dense_mat, sparse_mat, W, X, theta, time_lags, lambda_w, lambda_x, lambda_theta, eta, 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
lambda_w = 500
lambda_x = 500
lambda_theta = 500
eta = 0.03
time_lags = np.array([1, 2, 288])
d = time_lags.shape[0]
W = 0.1 * np.random.rand(dim1, rank)
X = 0.1 * np.random.rand(dim2, rank)
theta = 0.1 * np.random.rand(d, rank)
maxiter = 1000
TRMF(dense_mat, sparse_mat, W, X, theta, time_lags, lambda_w, lambda_x, lambda_theta, eta, 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