models/TRMF-CG.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Temporal Regularized Matrix Factorization\n",
    "\n",
    "The following Python codes reproduce the paper:\n",
    "\n",
    "> H.-F. Yu, N. Rao, and I. S. Dhillon. Temporal regularized matrixfactorization for high-dimensional time series prediction. Advances in neural information processing systems, vol. 29, pp. 847–855, 2016."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def compute_mape(var, var_hat):\n",
    "    return np.sum(np.abs(var - var_hat) / var) / var.shape[0]\n",
    "\n",
    "def compute_rmse(var, var_hat):\n",
    "    return np.sqrt(np.sum((var - var_hat) ** 2) / var.shape[0])\n",
    "\n",
    "def generate_Psi(T, d):\n",
    "    Psi = []\n",
    "    for k in range(0, d + 1):\n",
    "        if k == 0:\n",
    "            Psi.append(np.append(np.zeros((T - d, d)), np.eye(T - d), axis = 1))\n",
    "        else:\n",
    "            Psi.append(np.append(np.append(np.zeros((T - d, d - k)), np.eye(T - d), axis = 1), \n",
    "                                 np.zeros((T - d, k)), axis = 1))\n",
    "    return Psi\n",
    "\n",
    "def update_cg(w, r, q, Aq, rold):\n",
    "    alpha = rold / np.inner(q, Aq)\n",
    "    w = w + alpha * q\n",
    "    r = r - alpha * Aq\n",
    "    rnew = np.inner(r, r)\n",
    "    q = r + (rnew / rold) * q\n",
    "    return w, r, q, rnew\n",
    "\n",
    "def ell_w(ind, W, X, rho):\n",
    "    return X @ ((W.T @ X) * ind).T + rho * W\n",
    "\n",
    "def conj_grad_w(sparse_mat, ind, W, X, rho, maxiter = 5):\n",
    "    rank, dim1 = W.shape\n",
    "    w = np.reshape(W, -1, order = 'F')\n",
    "    r = np.reshape(X @ sparse_mat.T - ell_w(ind, W, X, rho), -1, order = 'F')\n",
    "    q = r.copy()\n",
    "    rold = np.inner(r, r)\n",
    "    for it in range(maxiter):\n",
    "        Q = np.reshape(q, (rank, dim1), order = 'F')\n",
    "        Aq = np.reshape(ell_w(ind, Q, X, rho), -1, order = 'F')\n",
    "        w, r, q, rold = update_cg(w, r, q, Aq, rold)\n",
    "    return np.reshape(w, (rank, dim1), order = 'F')\n",
    "\n",
    "def ell_x(sparse_mat, ind, W, X, A, Psi, d, rho, lmbda):\n",
    "    dim1, dim2 = sparse_mat.shape\n",
    "    rank = W.shape[0]\n",
    "    temp = np.zeros((d * rank, dim2 - d))\n",
    "    for k in range(1, d + 1):\n",
    "        temp[(k - 1) * rank : k * rank, :] = X[:, d - k : dim2 - k]\n",
    "    temp0 = X @ Psi[0].T - A @ temp\n",
    "    temp1 = np.zeros((rank, dim2))\n",
    "    for k in range(d):\n",
    "        temp1 += A[:, k * rank : (k + 1) * rank].T @ temp0 @ Psi[k + 1]\n",
    "    return W @ ((W.T @ X) * ind) + rho * X + lmbda * (temp0 @ Psi[0] - temp1)\n",
    "\n",
    "def conj_grad_x(sparse_mat, ind, W, X, A, Psi, d, rho, lmbda, maxiter = 5):\n",
    "    rank, dim1 = W.shape\n",
    "    dim2 = X.shape[1]\n",
    "    x = np.reshape(X, -1, order = 'F')\n",
    "    r = np.reshape(W @ sparse_mat - ell_x(sparse_mat, ind, W, X, A, Psi, d, rho, lmbda), -1, order = 'F')\n",
    "    q = r.copy()\n",
    "    rold = np.inner(r, r)\n",
    "    for it in range(maxiter):\n",
    "        Q = np.reshape(q, (rank, dim2), order = 'F')\n",
    "        Aq = np.reshape(ell_x(sparse_mat, ind, W, Q, A, Psi, d, rho, lmbda), -1, order = 'F')\n",
    "        x, r, q, rold = update_cg(x, r, q, Aq, rold)\n",
    "    return np.reshape(x, (rank, dim2), order = 'F')\n",
    "\n",
    "def ar(vec, d):\n",
    "    T = vec.shape[0]\n",
    "    Q = np.zeros((T - d, d))\n",
    "    for k in range(d):\n",
    "        Q[:, k] = vec[d - k - 1 : T - k - 1]\n",
    "    return np.linalg.pinv(Q) @ vec[d : T]\n",
    "\n",
    "def notmf(dense_mat, sparse_mat, rank, d, lmbda, rho, maxiter):\n",
    "    dim1, dim2 = sparse_mat.shape\n",
    "    W = 0.01 * np.random.randn(rank, dim1)\n",
    "    X = 0.01 * np.random.randn(rank, dim2)\n",
    "    A = 0.01 * np.random.randn(rank, d)\n",
    "    if np.isnan(sparse_mat).any() == False:\n",
    "        ind = sparse_mat != 0\n",
    "        pos_test = np.where((dense_mat != 0) & (sparse_mat == 0))\n",
    "    elif np.isnan(sparse_mat).any() == True:\n",
    "        pos_test = np.where((dense_mat != 0) & (np.isnan(sparse_mat)))\n",
    "        ind = ~np.isnan(sparse_mat)\n",
    "        sparse_mat[np.isnan(sparse_mat)] = 0\n",
    "    dense_test = dense_mat[pos_test]\n",
    "    del dense_mat\n",
    "    Psi = generate_Psi(dim2, d)\n",
    "    show_iter = 100\n",
    "    for it in range(maxiter):\n",
    "        W = conj_grad_w(sparse_mat, ind, W, X, rho)\n",
    "        A_mat = np.zeros((rank, d * rank))\n",
    "        for k in range(d):\n",
    "            A_mat[:, k * rank : (k + 1) * rank] = np.diag(A[:, k])\n",
    "        X = conj_grad_x(sparse_mat, ind, W, X, A_mat, Psi, d, rho, lmbda)\n",
    "        for r in range(rank):\n",
    "            A[r, :] = ar(X[r, :], d)\n",
    "        mat_hat = W.T @ X\n",
    "        if (it + 1) % show_iter == 0:\n",
    "            temp_hat = mat_hat[pos_test]\n",
    "            print('Iter: {}'.format(it + 1))\n",
    "            print('Loss function: {}'.format(np.linalg.norm((sparse_mat - mat_hat) * ind, 'fro')))\n",
    "            print('MAPE: {:.6}'.format(compute_mape(dense_test, temp_hat)))\n",
    "            print('RMSE: {:.6}'.format(compute_rmse(dense_test, temp_hat)))\n",
    "            print()\n",
    "    return mat_hat, W, X, A\n",
    "\n",
    "def notmf_dic(obs, W, X, A, d, lmbda, rho):\n",
    "    dim1, dim2 = obs.shape\n",
    "    rank = X.shape[0]\n",
    "    if np.isnan(obs).any() == False:\n",
    "        ind = obs != 0\n",
    "    elif np.isnan(obs).any() == True:\n",
    "        ind = ~np.isnan(obs)\n",
    "        obs[np.isnan(obs)] = 0\n",
    "    Psi = generate_Psi(dim2, d)\n",
    "    A_mat = np.zeros((rank, d * rank))\n",
    "    for k in range(d):\n",
    "        A_mat[:, k * rank : (k + 1) * rank] = np.diag(A[:, k])\n",
    "    X = conj_grad_x(obs, ind, W, X, A_mat, Psi, d, rho, lmbda)\n",
    "    for r in range(rank):\n",
    "        A[r, :] = ar(X[r, :], d)\n",
    "    return X, A\n",
    "\n",
    "def var4cast(X, A, d, delta):\n",
    "    dim1, dim2 = X.shape\n",
    "    X_hat = np.append(X, np.zeros((dim1, delta)), axis = 1)\n",
    "    A_mat = np.zeros((rank, d * rank))\n",
    "    for k in range(d):\n",
    "        A_mat[:, k * rank : (k + 1) * rank] = np.diag(A[:, k])\n",
    "    for t in range(delta):\n",
    "        X_hat[:, dim2 + t] = A_mat @ X_hat[:, dim2 + t - np.arange(1, d + 1)].T.reshape(dim1 * d)\n",
    "    return X_hat\n",
    "\n",
    "from ipywidgets import IntProgress\n",
    "from IPython.display import display\n",
    "\n",
    "def rolling4cast(dense_mat, sparse_mat, pred_step, delta, rank, d, lmbda, rho, maxiter):\n",
    "    dim1, T = sparse_mat.shape\n",
    "    start_time = T - pred_step\n",
    "    max_count = int(np.ceil(pred_step / delta))\n",
    "    mat_hat = np.zeros((dim1, max_count * delta))\n",
    "    f = IntProgress(min = 0, max = max_count) # instantiate the bar\n",
    "    display(f) # display the bar\n",
    "    for t in range(max_count):\n",
    "        if t == 0:\n",
    "            _, W, X, A = notmf(dense_mat[:, : start_time], sparse_mat[:, : start_time], \n",
    "                                 rank, d, lmbda, rho, maxiter)\n",
    "        else:\n",
    "            X, A = notmf_dic(sparse_mat[:, : start_time + t * delta], W, X_new, A, d, lmbda, rho)\n",
    "        X_new = var4cast(X, A, d, delta)\n",
    "        mat_hat[:, t * delta : (t + 1) * delta] = W.T @ X_new[:, - delta :]\n",
    "        f.value = t\n",
    "    small_dense_mat = dense_mat[:, start_time : T]\n",
    "    pos = np.where((small_dense_mat != 0) & (np.invert(np.isnan(small_dense_mat))))\n",
    "    mape = compute_mape(small_dense_mat[pos], mat_hat[pos])\n",
    "    rmse = compute_rmse(small_dense_mat[pos], mat_hat[pos])\n",
    "    print('Prediction MAPE: {:.6}'.format(mape))\n",
    "    print('Prediction RMSE: {:.6}'.format(rmse))\n",
    "    print()\n",
    "    return mat_hat, W, X, A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "dense_mat = np.load('../datasets/NYC-movement-data-set/hourly_speed_mat_2019_1.npz')['arr_0']\n",
    "for month in range(2, 4):\n",
    "    dense_mat = np.append(dense_mat, np.load('../datasets/NYC-movement-data-set/hourly_speed_mat_2019_{}.npz'.format(month))['arr_0'], axis = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f2757cd4f7604ad8930239df8894c8e5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "IntProgress(value=0, max=28)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Prediction MAPE: 0.21419\n",
      "Prediction RMSE: 4.86009\n",
      "\n",
      "Running time: 1310 seconds\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "start = time.time()\n",
    "dim1, dim2 = dense_mat.shape\n",
    "rank = 10\n",
    "pred_step = 7 * 24\n",
    "delta = 6\n",
    "d = 12\n",
    "lmbda = 1\n",
    "rho =  5\n",
    "maxiter = 50\n",
    "mat_hat, W, X, A = rolling4cast(dense_mat[:, : 24 * 7 * 10], dense_mat[:, : 24 * 7 * 10], \n",
    "                                pred_step, delta, rank, d, lmbda, rho, maxiter)\n",
    "end = time.time()\n",
    "print('Running time: %d seconds'%(end - start))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### License\n",
    "\n",
    "<div class=\"alert alert-block alert-danger\">\n",
    "<b>This work is released under the MIT license.</b>\n",
    "</div>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}