<jupyter_start><jupyter_text>Tensor Nuclear Norm minimization with Discrete Cosine Transform (TNN-DCT)This notebook shows how to implement a TNN-DCT imputer on some real-world data sets. For an in-depth discussion of TNN-DCT, please see [1].[1] <NAME>, <NAME>, <NAME> (2019). Low-rank Tensor Completion with a New Tensor Nuclear Norm Induced by Invertible Linear Transforms. CVPR 2019 [PDF]<jupyter_code>import numpy as np
from numpy.linalg import inv as inv<jupyter_output><empty_output><jupyter_text>Define kernelsWe start by introducing some necessary functions that relies on `Numpy`.<jupyter_code>def ten2mat(tensor, mode):
    return np.reshape(np.moveaxis(tensor, mode, 0), (tensor.shape[mode], -1), order = 'F')

def mat2ten(mat, dim, mode):
    index = list()
    index.append(mode)
    for i in range(dim.shape[0]):
        if i != mode:
            index.append(i)
    return np.moveaxis(np.reshape(mat, list(dim[index]), order = 'F'), 0, mode)
from scipy.fftpack import dctn, idctn

def tsvt_tnn(tensor, tau):
    dim = tensor.shape
    X = np.zeros(dim)
    tensor = dctn(tensor, axes = (2, ), norm = 'ortho')
    for t in range(dim[2]):
        u, s, v = np.linalg.svd(tensor[:, :, t], full_matrices = False)
        s = s - tau
        s[s < 0] = 0
        X[:, :, t] = u @ np.diag(s) @ v
    return idctn(X, axes = (2, ), norm = 'ortho')<jupyter_output><empty_output><jupyter_text>compute_mape: Compute the value of Mean Absolute Percentage Error (MAPE).compute_rmse: Compute the value of Root Mean Square Error (RMSE).> Note that $$\mathrm{MAPE}=\frac{1}{n} \sum_{i=1}^{n} \frac{\left|y_{i}-\hat{y}_{i}\right|}{y_{i}} \times 100, \quad\mathrm{RMSE}=\sqrt{\frac{1}{n} \sum_{i=1}^{n}\left(y_{i}-\hat{y}_{i}\right)^{2}},$$ where $n$ is the total number of estimated values, and $y_i$ and $\hat{y}_i$ are the actual value and its estimation, respectively.<jupyter_code>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])<jupyter_output><empty_output><jupyter_text>The main idea behind TNN-DCT is to implement tensor singular value thresholding on discrete cosine transformed data.<jupyter_code>def imputer(dense_tensor, sparse_tensor, rho0, epsilon, maxiter):
    """TNN-DCT-imputer."""
    dim = np.array(sparse_tensor.shape)
    pos_missing = np.where(sparse_tensor == 0)
    pos_test = np.where((dense_tensor != 0) & (sparse_tensor == 0))

    T = np.zeros(dim)                         # \boldsymbol{\mathcal{T}}
    Z = sparse_tensor.copy()                     # \boldsymbol{Z}
    Z[pos_missing] = np.mean(sparse_tensor[sparse_tensor != 0])
    it = 0
    last_ten = sparse_tensor.copy()
    snorm = np.linalg.norm(sparse_tensor)
    rho = rho0
    while True:
        rho = min(rho * 1.05, 1e5)
        X = tsvt_tnn(Z - T / rho, 1 / rho)
        temp2 = rho * X + T
        Z[pos_missing] = temp2[pos_missing] / rho
        T = T + rho * (X - Z)
        tol = np.linalg.norm((X - last_ten)) / snorm
        last_ten = X.copy()
        it += 1
        if it % 100 == 0:
            print('Iter: {}'.format(it))
            print('Tolerance: {:.6}'.format(tol))
            print('MAPE: {:.6}'.format(compute_mape(dense_tensor[pos_test], X[pos_test])))
            print('RMSE: {:.6}'.format(compute_rmse(dense_tensor[pos_test], X[pos_test])))
            print()
        if (tol < epsilon) or (it >= maxiter):
            break

    print('Total iteration: {}'.format(it))
    print('Tolerance: {:.6}'.format(tol))
    print('Imputation MAPE: {:.6}'.format(compute_mape(dense_tensor[pos_test], X[pos_test])))
    print('Imputation RMSE: {:.6}'.format(compute_rmse(dense_tensor[pos_test], X[pos_test])))
    print()

    return X<jupyter_output><empty_output><jupyter_text>Guangzhou dataWe generate **random missing (RM)** values on Guangzhou traffic speed data set.<jupyter_code>import scipy.io

tensor = scipy.io.loadmat('../datasets/Guangzhou-data-set/tensor.mat')
dense_tensor = tensor['tensor']
random_tensor = scipy.io.loadmat('../datasets/Guangzhou-data-set/random_tensor.mat')
random_tensor = random_tensor['random_tensor']

missing_rate = 0.2

### Random missing (RM) scenario:
binary_tensor = np.round(random_tensor + 0.5 - missing_rate)
sparse_tensor = np.multiply(dense_tensor, binary_tensor)

dense_tensor = np.transpose(dense_tensor, [0, 2, 1])
sparse_tensor = np.transpose(sparse_tensor, [0, 2, 1])

del tensor, random_tensor,binary_tensor
import time
start = time.time()
rho = 1e-4
epsilon = 1e-4
maxiter = 100
tensor_hat = imputer(dense_tensor, sparse_tensor, rho, epsilon, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
import scipy.io

tensor = scipy.io.loadmat('../datasets/Guangzhou-data-set/tensor.mat')
dense_tensor = tensor['tensor']
random_tensor = scipy.io.loadmat('../datasets/Guangzhou-data-set/random_tensor.mat')
random_tensor = random_tensor['random_tensor']

missing_rate = 0.4

### Random missing (RM) scenario:
binary_tensor = np.round(random_tensor + 0.5 - missing_rate)
sparse_tensor = np.multiply(dense_tensor, binary_tensor)

dense_tensor = np.transpose(dense_tensor, [0, 2, 1])
sparse_tensor = np.transpose(sparse_tensor, [0, 2, 1])

del tensor, random_tensor,binary_tensor
import time
start = time.time()
rho = 1e-4
epsilon = 1e-4
maxiter = 100
tensor_hat = imputer(dense_tensor, sparse_tensor, rho, epsilon, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))<jupyter_output>Iter: 100
Tolerance: 0.000288978
MAPE: 0.0791689
RMSE: 3.32679

Total iteration: 100
Tolerance: 0.000288978
Imputation MAPE: 0.0791689
Imputation RMSE: 3.32679

Running time: 63 seconds<jupyter_text>We generate **non-random missing (NM)** values on Guangzhou traffic speed data set. Then, we conduct the imputation experiment.<jupyter_code>import scipy.io

tensor = scipy.io.loadmat('../datasets/Guangzhou-data-set/tensor.mat')
dense_tensor = tensor['tensor']
random_matrix = scipy.io.loadmat('../datasets/Guangzhou-data-set/random_matrix.mat')
random_matrix = random_matrix['random_matrix']

missing_rate = 0.2

### Non-random missing (NM) scenario:
binary_tensor = np.zeros(dense_tensor.shape)
for i1 in range(dense_tensor.shape[0]):
    for i2 in range(dense_tensor.shape[1]):
        binary_tensor[i1, i2, :] = np.round(random_matrix[i1, i2] + 0.5 - missing_rate)
sparse_tensor = np.multiply(dense_tensor, binary_tensor)

dense_tensor = np.transpose(dense_tensor, [0, 2, 1])
sparse_tensor = np.transpose(sparse_tensor, [0, 2, 1])

del tensor, random_matrix, binary_tensor
import time
start = time.time()
rho = 1e-4
epsilon = 1e-4
maxiter = 100
tensor_hat = imputer(dense_tensor, sparse_tensor, rho, epsilon, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
import scipy.io

tensor = scipy.io.loadmat('../datasets/Guangzhou-data-set/tensor.mat')
dense_tensor = tensor['tensor']
random_matrix = scipy.io.loadmat('../datasets/Guangzhou-data-set/random_matrix.mat')
random_matrix = random_matrix['random_matrix']

missing_rate = 0.4

### Non-random missing (NM) scenario:
binary_tensor = np.zeros(dense_tensor.shape)
for i1 in range(dense_tensor.shape[0]):
    for i2 in range(dense_tensor.shape[1]):
        binary_tensor[i1, i2, :] = np.round(random_matrix[i1, i2] + 0.5 - missing_rate)
sparse_tensor = np.multiply(dense_tensor, binary_tensor)

dense_tensor = np.transpose(dense_tensor, [0, 2, 1])
sparse_tensor = np.transpose(sparse_tensor, [0, 2, 1])

del tensor, random_matrix, binary_tensor
import time
start = time.time()
rho = 1e-4
epsilon = 1e-4
maxiter = 100
tensor_hat = imputer(dense_tensor, sparse_tensor, rho, epsilon, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))<jupyter_output>Iter: 100
Tolerance: 0.000342577
MAPE: 0.116677
RMSE: 4.78242

Total iteration: 100
Tolerance: 0.000342577
Imputation MAPE: 0.116677
Imputation RMSE: 4.78242

Running time: 64 seconds<jupyter_text>PeMS data<jupyter_code>dense_mat = np.load('../datasets/PeMS-data-set/pems.npy')
random_tensor = np.load('../datasets/PeMS-data-set/random_tensor.npy')

missing_rate = 0.2

### Random missing (RM) scenario:
binary_tensor = np.round(random_tensor + 0.5 - missing_rate)
sparse_mat = np.multiply(dense_mat, ten2mat(binary_tensor, 0))

sparse_tensor = mat2ten(sparse_mat, np.array(binary_tensor.shape), 0)
dense_tensor = mat2ten(dense_mat, np.array(binary_tensor.shape), 0)

del dense_mat, random_tensor, binary_tensor
import time
start = time.time()
rho = 1e-4
epsilon = 1e-4
maxiter = 100
tensor_hat = imputer(dense_tensor, sparse_tensor, rho, epsilon, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
dense_mat = np.load('../datasets/PeMS-data-set/pems.npy')
random_tensor = np.load('../datasets/PeMS-data-set/random_tensor.npy')

missing_rate = 0.4

### Random missing (RM) scenario:
binary_tensor = np.round(random_tensor + 0.5 - missing_rate)
sparse_mat = np.multiply(dense_mat, ten2mat(binary_tensor, 0))

sparse_tensor = mat2ten(sparse_mat, np.array(binary_tensor.shape), 0)
dense_tensor = mat2ten(dense_mat, np.array(binary_tensor.shape), 0)

del dense_mat, random_tensor, binary_tensor
import time
start = time.time()
rho = 1e-4
epsilon = 1e-4
maxiter = 100
tensor_hat = imputer(dense_tensor, sparse_tensor, rho, epsilon, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
dense_mat = np.load('../datasets/PeMS-data-set/pems.npy')
random_matrix = np.load('../datasets/PeMS-data-set/random_matrix.npy')

missing_rate = 0.2

### Nonrandom missing (NM) scenario:
binary_tensor = np.zeros((dense_mat.shape[0], 288, 44))
for i1 in range(dense_mat.shape[0]):
    for i2 in range(44):
        binary_tensor[i1,:,i2] = np.round(random_matrix[i1,i2] + 0.5 - missing_rate)
binary_mat = ten2mat(binary_tensor, 0)
sparse_mat = np.multiply(dense_mat, binary_mat)

sparse_tensor = mat2ten(sparse_mat, np.array(binary_tensor.shape), 0)
dense_tensor = mat2ten(dense_mat, np.array(binary_tensor.shape), 0)

del dense_mat, random_matrix, binary_tensor
import time
start = time.time()
rho = 1e-4
epsilon = 1e-4
maxiter = 100
tensor_hat = imputer(dense_tensor, sparse_tensor, rho, epsilon, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
dense_mat = np.load('../datasets/PeMS-data-set/pems.npy')
random_matrix = np.load('../datasets/PeMS-data-set/random_matrix.npy')

missing_rate = 0.4

### Nonrandom missing (NM) scenario:
binary_tensor = np.zeros((dense_mat.shape[0], 288, 44))
for i1 in range(dense_mat.shape[0]):
    for i2 in range(44):
        binary_tensor[i1,:,i2] = np.round(random_matrix[i1,i2] + 0.5 - missing_rate)
binary_mat = ten2mat(binary_tensor, 0)
sparse_mat = np.multiply(dense_mat, binary_mat)

sparse_tensor = mat2ten(sparse_mat, np.array(binary_tensor.shape), 0)
dense_tensor = mat2ten(dense_mat, np.array(binary_tensor.shape), 0)

del dense_mat, random_matrix, binary_tensor
import time
start = time.time()
rho = 1e-4
epsilon = 1e-4
maxiter = 100
tensor_hat = imputer(dense_tensor, sparse_tensor, rho, epsilon, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))<jupyter_output>Iter: 100
Tolerance: 0.000254299
MAPE: 0.105213
RMSE: 6.76111

Total iteration: 100
Tolerance: 0.000254299
Imputation MAPE: 0.105213
Imputation RMSE: 6.76111

Running time: 78 seconds<jupyter_text>Electricity data<jupyter_code>dense_mat = np.load('../datasets/Electricity-data-set/electricity35.npy')
random_tensor = np.load('../datasets/Electricity-data-set/random_tensor.npy')

missing_rate = 0.2

### Random missing (RM) scenario:
binary_tensor = np.round(random_tensor + 0.5 - missing_rate)
sparse_mat = np.multiply(dense_mat, ten2mat(binary_tensor, 0))

sparse_tensor = mat2ten(sparse_mat, np.array(binary_tensor.shape), 0)
dense_tensor = mat2ten(dense_mat, np.array(binary_tensor.shape), 0)

del dense_mat, random_tensor, binary_tensor
import time
start = time.time()
rho = 1e-6
epsilon = 1e-4
maxiter = 100
tensor_hat = imputer(dense_tensor, sparse_tensor, rho, epsilon, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
dense_mat = np.load('../datasets/Electricity-data-set/electricity35.npy')
random_tensor = np.load('../datasets/Electricity-data-set/random_tensor.npy')

missing_rate = 0.4

### Random missing (RM) scenario:
binary_tensor = np.round(random_tensor + 0.5 - missing_rate)
sparse_mat = np.multiply(dense_mat, ten2mat(binary_tensor, 0))

sparse_tensor = mat2ten(sparse_mat, np.array(binary_tensor.shape), 0)
dense_tensor = mat2ten(dense_mat, np.array(binary_tensor.shape), 0)

del dense_mat, random_tensor, binary_tensor
import time
start = time.time()
rho = 1e-6
epsilon = 1e-4
maxiter = 100
tensor_hat = imputer(dense_tensor, sparse_tensor, rho, epsilon, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
dense_mat = np.load('../datasets/Electricity-data-set/electricity35.npy')
random_matrix = np.load('../datasets/Electricity-data-set/random_matrix.npy')

missing_rate = 0.2

### Nonrandom missing (NM) scenario:
binary_tensor = np.zeros((dense_mat.shape[0], 24, 35))
for i1 in range(dense_mat.shape[0]):
    for i2 in range(35):
        binary_tensor[i1,:,i2] = np.round(random_matrix[i1,i2] + 0.5 - missing_rate)
binary_mat = ten2mat(binary_tensor, 0)
sparse_mat = np.multiply(dense_mat, binary_mat)

sparse_tensor = mat2ten(sparse_mat, np.array(binary_tensor.shape), 0)
dense_tensor = mat2ten(dense_mat, np.array(binary_tensor.shape), 0)

del dense_mat, random_matrix, binary_tensor
import time
start = time.time()
rho = 1e-6
epsilon = 1e-4
maxiter = 100
tensor_hat = imputer(dense_tensor, sparse_tensor, rho, epsilon, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
dense_mat = np.load('../datasets/Electricity-data-set/electricity35.npy')
random_matrix = np.load('../datasets/Electricity-data-set/random_matrix.npy')

missing_rate = 0.4

### Nonrandom missing (NM) scenario:
binary_tensor = np.zeros((dense_mat.shape[0], 24, 35))
for i1 in range(dense_mat.shape[0]):
    for i2 in range(35):
        binary_tensor[i1,:,i2] = np.round(random_matrix[i1,i2] + 0.5 - missing_rate)
binary_mat = ten2mat(binary_tensor, 0)
sparse_mat = np.multiply(dense_mat, binary_mat)

sparse_tensor = mat2ten(sparse_mat, np.array(binary_tensor.shape), 0)
dense_tensor = mat2ten(dense_mat, np.array(binary_tensor.shape), 0)

del dense_mat, random_matrix, binary_tensor
import time
start = time.time()
rho = 1e-6
epsilon = 1e-4
maxiter = 100
tensor_hat = imputer(dense_tensor, sparse_tensor, rho, epsilon, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))<jupyter_output>Iter: 100
Tolerance: 0.000179406
MAPE: 0.188558
RMSE: 1024.41

Total iteration: 100
Tolerance: 0.000179406
Imputation MAPE: 0.188558
Imputation RMSE: 1024.41

Running time: 2 seconds
{'fraction_non_alphanumeric': 0.14303063014437997, 'fraction_numerical': 0.0378491431655647, 'mean_word_length': 4.851954204500593, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 34, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '522246', 'license': 'MIT', 'n_tokens_mistral': 6004, 'n_tokens_neox': 5470, 'n_words': 1189, 'path': 'baselines/TNN-DCT-imputer.ipynb', 'repo_name': 'vishalbelsare/tensor-learning', 'size': '14830'}