Covid-19 Pandemic Data Analysis using Tensor Methods ^††thanks: Submitted to the editors December 29, 2022. \fundingThis work was funded by National Science Foundation under Grant No. DMS-1439786 and Grant No. MCB-2126374.

Tensor decomposition is a powerful tool in data analysis, computer vision, scientific computing, machine learning and many other fields. In fact, tensor models for dimensionality reduction have been highly effective in machine learning applications like classification and regression. Tensor decomposition on its own have been successful in many modern big data analysis [11, 8, 24, 31].

In this work, we focus on analyzing Covid-19 pandemic data [28, 1] by using tensor decomposition models. Our goals are to predict future infection, locate and identify hotspots from data-set coming from multiple sources across different spatial regions over time. One is interested in determining where and when changes occur in the pattern. The strategy is set up an optimization which separates the data in following: let $\mathcal{Y}=\mathcal{L}+\mathcal{S}$ where $\mathcal{Y}$ is the given tensor, $\mathcal{L}$ is a low rank reconstructed tensor of $\mathcal{Y}$ and $\mathcal{S}$ is the sparse tensor. In video processing, the original video is separated into its background and foreground subspace to detect anomalous activities [11, 32]. The tensor $\mathcal{L}$ is the background, and $\mathcal{S}$ is the foreground. The sparse tensor $\mathcal{S}$ can provide anomalous activities. Similarly, $\mathcal{S}$ will contain hotspot occurrence. To achieve this separation, we formulate the following:

where $\mathcal{L}=\sum^{R}_{r=1}\bm{\alpha}_{r}\bf a_{r}\circ\bf b_{r}\circ\bf c_{r}$, $\|\cdot\|_{F}$ is the tensor Frobenius norm, and $\|\cdot\|_{\ell_{1}}$ is the vector one norm. The optimization problem (1) is neither convex nor differentiable [12]. In addition, this formulation is amenable to tensor completion problems where one can find missing data. This is basis for predicting future Covid-19 infection cases in presence of previous pandemic data.

Most of the algorithms suggested for solving (1) are based on ALS (alternating least squares) [12]. ALS is fast, easy to program and effective. However, ALS has some drawbacks. There is a need for more efficient and accurate methods. Thus, we propose a sampling method for ALS to leverage the efficiency of ALS and to allow bigger tensor data.

We denote a vector by a bold lower-case letter $\mathbf{a}$. The bold upper-case letter $\mathbf{A}$ represents a matrix and the symbol of a tensor is a calligraphic letter $\mathcal{A}$. Throughout this paper, we focus on third-order tensors $\mathcal{A}=(a_{ijk})\in\mathbb{R}^{I\times J\times K}$ of three indices $1\leq i\leq I,1\leq j\leq J$ and $1\leq k\leq K$, but these can be easily extended to tensors of arbitrary order greater or equal to three.

The three kinds of matricization for third-order $\mathcal{A}$ are $\mathbf{A}_{(1)},\mathbf{A}_{(2)}$ and $\mathbf{A}_{(3)}$, according to respectively arranging the column, row, and tube fibers to be columns of matrices, which are defined by fixing every index but one and denoted by $\mathbf{a}_{:jk}$, $\mathbf{a}_{i:k}$ and $\mathbf{a}_{ij:}$ respectively. We also consider the vectorization for $\mathcal{A}$ to obtain a row vector $\mathbf{a}$ such the elements of $\mathcal{A}$ are arranged according to $k$ varying faster than $j$ and $j$ varying faster than $i$, i.e., $\mathbf{a}=(a_{111},\cdots,a_{11K},a_{121},\cdots,a_{12K},\cdots,a_{1J1},\cdots,a_{1JK},\cdots)$. Kronecker Product of two matrices $\mathbf{A}\in\mathbb{R}^{m\times n}$ and $\mathbf{B}\in\mathbb{R}^{p\times q}$ is denoted as $\mathbf{A}\otimes\mathbf{B}\in\mathbb{R}^{mp\times nq}$ and obtained as the product of each element of $\mathbf{A}$ and the matrix$\mathbf{B}$. Khatri-Rao Product of two matrices $\mathbf{A}$ and $\mathbf{B}$ with same columns is the column-wise Kronecker product: $\mathbf{A\odot B}=[\mathbf{a}_{1}\otimes\mathbf{b}_{1}\hskip 4.26773pt\ldots\hskip 4.26773pt\mathbf{a}_{R}\otimes\mathbf{b}_{R}].$ The outer product of a rank-one third order tensor is denoted as $\mathbf{a}\circ\mathbf{b}\circ\mathbf{c}\in\mathbb{R}^{I\times J\times K}$ of three nonzero vectors $\mathbf{a},\mathbf{b}$ and $\mathbf{c}$ is a rank-one tensor with elements $a_{i}b_{j}c_{k}$ for all the indices, i.e. the matricizations of $\mathbf{a}\circ\mathbf{b}\circ\mathbf{c}\in\mathbb{R}^{I\times J\times K}$ are rank-one matrices. A canonical polyadic (CP) decomposition of $\mathcal{A}\in\mathbb{R}^{I\times J\times K}$ expresses $\mathcal{A}$ as a sum of rank-one outer products:

where $\mathbf{a}_{r}\in\mathbb{R}^{I},\mathbf{b}_{r}\in\mathbb{R}^{J},\mathbf{c}_{r}\in\mathbb{R}^{K}$ for $1\leq r\leq R$ and and $\circ$ is the outer product. The outer product $\mathbf{a}_{r}\circ\mathbf{b}_{r}\circ\mathbf{c}_{r}$ is a rank-one component and the integer $R$ is the number of rank-one components in tensor $\mathcal{A}$. The minimal number $R$ such that the decomposition (2) holds is the rank of tensor $\mathcal{A}$, which is denoted by $\mbox{rank}(\mathcal{A})$. For any tensor $\mathcal{A}\in\mathbb{R}^{I\times J\times K}$, $\mbox{rank}(\mathcal{A})$ has an upper bound $\min\{IJ,JK,IK\}$ [13]. In fact, tensor rank is NP-hard over $\mathbb{R}$ and $\mathbb{C}$ [7]. Another standard tensor decomposition is higher-order Tucker (HOT) decomposition. HOT is a generalization of matrix SVD where $m\times n$ matrix $\mathbf{M}$ has a factorizaion $\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{T}$, where $\mathbf{U}$ and $\mathbf{V}$ are orthogonal matrices and $\Sigma$ is a diagonal matrix with $\sigma_{ij}=0$ if $i\neq 0$ and otherwise $\sigma_{ii}>=0$. Its generalization in third-order tensors is

where $\mathcal{T}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}}$ is the given tensor, $\mathcal{G}\in\mathbb{R}^{R_{1}\times R_{2}\times R_{3}}$ is the core tensor and $\mathbf{U^{(i)}}\in\mathbb{R}^{I_{i}\times R_{i}}$ for $i=1,2,3$ is an orthogonal matrix. The Tucker contracted product $\times_{1}$ is defined as $\mathcal{G}\times_{1}\mathbf{U}^{(1)}=\sum_{r_{1}}\mathcal{G}_{r_{1}r_{2}r_{3}}\mathbf{U}^{(1)}_{i_{1}r_{1}}\in\mathbb{R}^{I_{1}\times R_{2}\times R_{3}}$.

Given a tensor $\mathcal{X}\in\mathbb{R}^{I\times J\times K}$ and a fixed positive integer $R$, the ALS algorithm solves three independent least squares problems alternatively. The least squares problems can be solved via various methods such as QR factorization, Cholesky factorization and etc. However, it requires of an $\mathcal{O}(R^{3})$ floating point operations per second. To reduce this we let $S\subset\{1,\ldots R\}$ be the set of sample indices and $A_{s},B_{s}$ and $C_{s}$ represent the sub-matrices obtained by choosing the columns of $A,B$, and $C$ according to the index set $S$ respectively. The partial derivatives of the objective function $f$ with respect to the blocks $A_{s},B_{s}$ and $C_{s}$ are

and

the stationary points of the above equations can be obtained by setting each gradient equal to zero. For instance $\nabla_{A_{S}}f=0$ implies the following normal equation

where $A^{S}$ represents the sub matrix of $A$ obtained by sampling the rows of $A$ corresponding to the sampling set $S$. Similar results can be derived by solving the equations $\nabla_{B_{S}}f=0$ and $\nabla_{C_{S}}f=0$. This reduces the latter computational complexity to $\mathcal{O}(\max\{|S|^{3}\})$. When $I,J\,\text{and}\,K$ are relatively large, the reduction could be significant. The sampling set $S$ is chosen based on the performance of each block variable in each iteration. For instance, if the update for the block $A_{i}$ yields more decrease in the objective function compared to the block $A_{j}$, the index $i$ is replaced by $j$ in the next iteration. The differentiation of ALS can be found here [26]. The derivation of SMALS and further analysis can be found [10].

In [31], an iterative method based on proximal algorithms called low-rank approximation of tensors (LRAT) was proposed to solve the minimum rank optimization:

where $\mathcal{L}=\sum^{R}_{r=1}\alpha_{r}\mathbf{a_{r}}\circ\mathbf{b_{r}}\circ\mathbf{c_{r}}$, $\mathcal{C}$, a given data tensor. The implementation is in Algorithm 4. In a recent work, [8], a more practical choice of the parameter $\sigma$ led to a more efficient and accurate algorithm. The practical regularization method is based on a flexible Golub-Kahan (fGK) process. In Figures 14 and 15, we implemented the LRAT + fGK algorithm with the convergence plot; the sparse model predicts the general phenomena of the Covid-19 cases.

.
In figure $7$, the error plot shows the difference between the original data and estimation. The horizontal line detects an anomaly and this spot coincides with a spike in the number of infections.

In this work, we apply various tensor models and tensor algorithms to analyze Covid-19 data. The standard tensor models, CP and HOOI, with their off-the-shelves algorithms, ALS and HOOI, are tested against a new sampling method for ALS (SMALS). The numerical results are very promising as it cuts the down time while keeping the Frobenious norm errors relatively consistent with ALS. Here tensor sparse model [31, 8] is used as the model for prediction of future Covid-19 infection cases as a tensor completion problem. The numerical results are impressive as the tensor completion algorithm can predict infection a week ahead as well as a quarter ahead. Moreover, the tensor sparse model can locate which counties exhibit hotspots. The sparse tensor model is based on proximal algorithms and flexible hybrid method by Golub-Kahan methods for efficient practical implementation.

In our future work, we would like to have a more mathematical and methodical technique in locating and detecting hotspots. We have used $l_{1}$ minimization in the tensor completion; we will work on the efficacy of $l_{0}$ minimization of low-rank approximation of CP decomposition and tensor completion algorithm in the proximal framework.

This material is based upon work supported by the National Science Foundation under Grant No. DMS-1439786 while the author, C. Navasca, was in residence at the Institute for Computational and Experimental Research in Mathematics in Providence, RI, during the Model and Dimension Reduction in Uncertain and Dynamic Systems Program. C. Navasca is also in part supported by National Science Foundation No. MCB-2126374.