Total Variation Regularized Tensor-on-scalar Regression

TL;DR

This paper introduces TVTR, a new tensor-on-scalar regression method with total variation regularization, designed for efficient, scalable estimation of associations in high-dimensional tensor data, with proven consistency and superior performance in simulations and real data.

Contribution

The paper presents a novel TV regularized tensor-on-scalar regression method with an efficient algorithm, convergence proof, and demonstrated superior performance over existing methods.

Findings

TVTR outperforms existing methods in simulations

The algorithm is efficient and scalable in distributed platforms

The method is applicable to brain imaging and activity tracking data

Abstract

In this paper, we propose Total Variation Regularized Tensor-on-scalar Regression(TVTR), a novel method for estimating the association between a tensor outcome (a one dimensional or multidimensional array) and scalar predictors. While the statistical developments proposed here were motivated by the brain mapping and activity tracking, the methodology is designed and presented in generality and is applicable to many other areas of scientific research. The estimator is the solution of a penalized regression problem where the objective is the sum of square error plus a total variation (TV) regularization on the predicted mean across all subjects. We propose an algorithm for the parameter estimation, which is efficient and scalable in distributed computing platform. Proof of the algorithm convergence is provided and the statistical consistency of the estimator is presented via an oracle inequality. We presented 1D and 2D simulation results, and demonstrate that TVTR outperforms existing methods in most cases. We also demonstrate the general applicability of the method by two real data examples including the analysis of the 1D accelerometry subsample of a large community-based study for mood disorders and the analysis of the 3D MRI data from the attention deficient/hyperactive deficient (ADHD) 200 consortium.

Tables

Table 1. Table 1 : Mean Deviation of Coefficients for 1D settings

	sample size	TVTR	periodic_TVTR	FoSR_bspline	FoSR_fourier
Setting 1	25	0.045(0.008)	-	0.054(0.015)	0.053(0.015)
	50	0.024(0.004)	-	0.024(0.006)	0.024(0.006)
	100	0.015(0.002)	-	0.012(0.003)	0.011(0.003)
Setting 2	25	0.076(0.009)	0.033(0.007)	0.081(0.015)	0.078(0.015)
	50	0.041(0.004)	0.020(0.004)	0.051(0.006)	0.049(0.006)
	100	0.020(0.002)	0.012(0.002)	0.039(0.003)	0.037(0.004)

Table 2. Table 2 : Mean Deviation of Coefficients for 2D settings

	sample size	TVTR	Envelope	FoSR	TV_OLS	OLS_TV
Setting 1	25	0.298 (0.038)	0.349(0.028)	0.357 (0.019)	0.501 (0.044)	0.424 (0.048)
	50	0.217 (0.022)	0.241(0.013)	0.339 (0.010)	0.361 (0.028)	0.313(0.033)
	100	0.200 (0.016)	0.169(0.007)	0.330 (0.004)	0.292 (0.017)	0.230 (0.013)
Setting 2	25	0.265 (0.008)	0.561(0.045)	0.321 (0.003)	0.275 (0.006)	0.919 (0.087)
	50	0.192 (0.007)	0.387(0.019)	0.318 (0.001)	0.266 (0.003)	0.545 (0.033)
	100	0.126 (0.006)	0.267(0.010)	0.318 (0.001)	0.259 (0.003)	0.277 (0.012)

Equations

$$\Vert Y-\theta {\Vert }^2+\lambda \sum _{i=1}^{M-1}|{\theta }_{i+1}-{\theta }_i{\Vert }_{{\ell }_1},$$

$$\Vert Y-\Theta {\Vert }^2+\lambda \Vert D\text{vec}(\Theta ){\Vert }_{{\ell }_1},$$

$$Y_i=$$

$$\text{vec}(Y^T)=\text{vec}((X\Gamma {)}^T)+\mathbf{ϵ},$$

$$I_{nM}-H_v=(I-H)\otimes I_M\in {\mathbb{R}}^{nM\times nM}.$$

$$(I-H_v)\text{vec}((X\Gamma {)}^T)=$$

$$D_{i,j}=\left\{\begin{aligned} -1&\mbox{ if the edge $e_{j}$ leaves vertex $v_{i}$;}\\ 1&\mbox{ if it enters vertex $v_{i}$;}\\ 0&\mbox{ otherwise}\end{aligned}\right.$$

$$\underset{\Gamma }{\min }\frac{1}{2}\Vert Y-X\Gamma {\Vert }_F^2+\lambda \Vert X\Gamma D{\Vert }_{{\ell }_1}.$$

$$\begin{array}{rl}\underset{\theta }{\min }& \frac{1}{2}\Vert \text{vec}(Y^T)-\theta {\Vert }_F^2+\lambda \Vert D_v^T\theta {\Vert }_{{\ell }_1};\\ \text{s.t.}& (I-H_v)\theta =0.\end{array}$$

$$\min$$

$$s.t.$$

$${\kappa }_T:=\underset{\theta \in {\mathbb{R}}^T}{\inf }\frac{\sqrt{|T|}\Vert \theta \Vert }{\Vert (\theta D{)}_T{\Vert }_{{\ell }_1}};\kappa =\underset{T\subset [m]}{\inf }{\kappa }_T$$

$$\Vert {\theta }^{\ast }-\hat{\theta }{\Vert }_F^2\le$$

$$+64{\sigma }^2\log \left(\frac{2enM}{\delta }\right)+8{\rho }^2{\sigma }^2\log \left(\frac{mnM}{\delta }\right){\kappa }_T^{-2}|T|$$

$$\Vert {\theta }^{\ast }-\hat{\theta }{\Vert }_F^2\le$$

$$+64{\sigma }^2\log \left(\frac{2enM}{\delta }\right)+\frac{2{\rho }^2{\sigma }^2\log \left(\frac{mnM}{\delta }\right)|T|}{\min \{d,|T|\}}$$

$$\frac{1}{Mp}\Vert \hat{\Gamma }-{\Gamma }^{\ast }{\Vert }_F^2\le$$

$$+\frac{64{\sigma }^2\log \left(\frac{2enM}{\delta }\right)}{nMp}+\frac{8{\rho }^2{\sigma }^2\log \left(\frac{mnM}{\delta }\right)|T|}{nMp{\kappa }_T^2}$$

$$Y=0.3\sin (\pi t/100)+0.5X_1\cos (\pi t/100)-0.3X_2\sin (\pi t/50)+0.5\cos (\pi t/25)+2\mathcal{N}(0,1).$$

$$Y_{ij}=\mathbb{I}(1\le j\le 20)+\mathbb{I}(101\le j\le 120)+0.5X_{1i}\left(\mathbb{I}(31\le j\le 70)+\mathbb{I}(131\le j\le 170)\right)$$

$$-X_{2i}\left(\mathbb{I}(71\le j\le 80)+\mathbb{I}(171\le j\le 180)\right)+X_{3i}\left(\mathbb{I}(61\le j\le 100)+\mathbb{I}(161\le j\le 200)\right)+2\mathcal{N}(0,1).$$

$$Y=X_1{\Gamma }_1+X_2{\Gamma }_2+X_3{\Gamma }_3+ϵ$$

$$\min$$

$$s.t.$$

$$f(\theta )=\frac{1}{2}\Vert y-\theta {\Vert }^2,g_1(\eta )=\delta ((I-H_v)\eta =0),g_2(\mu )=\Vert D_v^T\mu {\Vert }_1$$

$$A_1=\left[\begin{array}{c}-I\\ -I\end{array}\right],A_2=\left[\begin{array}{c}I\\ 0\end{array}\right],A_3=\left[\begin{array}{c}0\\ I\end{array}\right]\in {\mathbb{R}}^{2nM\times nM}$$

$$\mathcal{L}(\theta ,\eta ,\mu ,U,V,\rho )=f(\theta )+g_1(\eta )+g_2(\mu )+{\left[\begin{array}{c}U\\ V\end{array}\right]}^T(A_1\theta +A_2\eta +A_3\mu )+\frac{\rho }{2}\Vert A_1\theta +A_2\eta +A_3\mu {\Vert }^2.$$

$${\theta }^{(k+1)}=\arg \underset{\theta }{\min }f(\theta )+{\left[\begin{array}{c}U\\ V\end{array}\right]}^T(A_1\theta )+\frac{\rho }{2}\Vert A_1\theta +A_2{\eta }^{(k)}+A_3{\mu }^{(k)}{\Vert }^2$$

$${\eta }^{(k+1)}=\arg \underset{\eta }{\min }{g}_1(\eta )+{\left[\begin{array}{c}U\\ V\end{array}\right]}^T(A_2\eta )+\frac{\rho }{2}\Vert A_1{\theta }^{(k+1)}+A_2\eta +A_3{\mu }^{(k)}{\Vert }^2$$

$${\mu }^{(k+1)}=\arg \underset{\mu }{\min }{g}_2(\mu )+{\left[\begin{array}{c}U\\ V\end{array}\right]}^T(A_3\mu )+\frac{\rho }{2}\Vert A_1{\theta }^{(k+1)}+A_2{\eta }^{(k+1)}+A_3\mu {\Vert }^2$$

$$\left[\begin{array}{c}{U}^{(k+1)}\\ V^{(k+1)}\end{array}\right]=\left[\begin{array}{c}{U}^{(k)}\\ V^{(k)}\end{array}\right]+\rho (A_1{\theta }^{(k+1)}+A_2{\eta }^{(k+1)}+A_3\mu )$$

$$f(\theta )-f({\theta }^{(k+1)})+(\theta -{\theta }^{(k+1)})\left(-A_1^T({\lambda }^k-\rho (A_1\theta +A_2\eta +A_3\mu ))\right)\ge 0,\forall \theta \in {\mathbb{R}}^{nM}$$

$$g_1(\eta )-g_1({\eta }^{(k+1)})+(\eta -{\eta }^{(k+1)})\left(-A_2^T({\lambda }^k-\rho (A_1\theta +A_2\eta +A_3\mu ))\right)\ge 0,\forall \eta \in {\mathbb{R}}^{nM}$$

$$g_2(\mu )-g_2({\mu }^{(k+1)})+(\mu -{\mu }^{(k+1)})\left(-A_3^T({\lambda }^k-\rho (A_1\theta +A_2\eta +A_3\mu ))\right)\ge 0,\forall \mu \in {\mathbb{R}}^{nM}$$

$$f(\theta )-f({\theta }^{(k+1)})+(\theta -{\theta }^{(k+1)})\left(-A_1^T({\lambda }^k-\rho (A_1\theta +A_2\eta +A_3\mu ))\right)\ge 0,\forall \theta \in {\mathbb{R}}^{nM}$$

$$g_1(\eta )-g_1({\eta }^{(k+1)})+(\eta -{\eta }^{(k+1)})\left(-A_2^T({\lambda }^k-\rho (A_1\theta +A_2\eta ))\right)\ge 0,\forall \eta \in {\mathbb{R}}^{nM}$$

$$g_2(\mu )-g_2({\mu }^{(k+1)})+(\mu -{\mu }^{(k+1)})\left(-A_3^T({\lambda }^k-\rho (A_1\theta +A_3\mu ))\right)\ge 0,\forall \mu \in {\mathbb{R}}^{nM}$$

Full text

Keywords: Penalized regression, Tensor-on-scalar regression, Multivariate linear regression, Total variation denoising

1 Introduction

High-dimensional data continue to proliferate in modern biomedical research where many of outcomes are one-dimensional or multidimensional arrays (tensors) with complex spatial and temporal structures. For example, wearable computing sensors such as accelerometers and smart phones, have been increasingly used in epidemiological cohort studies such as National Health and Nutrition Examination Survey (NHANES) and UK Biobank Naslund et al. (2015); Pantelopoulos & Bourbakis (2010); Patel et al. (2015). Instead of collecting subjects’ physical activity levels through self-reported survey questionnaire that are scarcely sampled and subject to reporting biases, these devices have allowed us to observe physical activity intensities continuously in high frequencies over an extended period and have shown to be a more powerful way to track activity markers and circadian rhythm patterns that link with health-related behaviors and clinical outcomes such as mood, obesity and aging. Meanwhile in neuroscience, sustained growth in computing power and the increasing affordability to collect high-throughput medical images has made various modalities of brain imaging such as MRI, fMRI and PET available for researchers to search for promising brain markers that are associated with or predict the patient-specific diagnosis and prognosis. To untangle such associations, appropriate statistical methods that acknowledge the multidimensional structures of the outcome data and correlations induced by time or spatial continuity are needed. In particular, we were motivated by two specific studies that collect 1D curve and 3D array as outcome measures. The first example came from the accelerometry subsample of a large community-based study for mood disorders, the National Institute of Mental Health Family Study of Affective Spectrum Disorder (Merikangas et al., 2014). The study consist minute-by-minute physical activity intensity measures (i.e., activity counts) from a set of about 350 participants using the Philips Respronics 24 hours per day over 2 weeks. Despite that the observed activity data are highly noisy as shown in Figure 3, we hypothesize that the underlying activity status (inactive, sedentary, moderate and vigorous) transit smoothly from one to another. The second study was from attention deficient/hyperactive deficient (ADHD) 200 consortium as part of the Human Connectome Initiative. The study included both structural and BOLD functional magnetic resonance imaging (MRI) scans of 776 subjects, including 491 typically developing controls (TDs) and 285 ADHD children from multiple centers. We were interested in investigating regions where gray matter concentration differs between TD and AD children by taking the voxel-based morphometry images as outcome in the model. In addition to the low signal-to-noise ratio, the images were also formulated in a three-dimensional spatially coherent space where gray matter exists.

As previously noted, the challenges of developing a general and robust analytic framework for tensor response data are three-fold: the high dimensionality of the response, the low signal-to-noise ratio, and the complex dependency structures. While there exist an enormous body of literature tackling regression with high- or ultrahigh-dimensional data, they often deal with scalar outcomes and high-dimensional data as predictors through various regularizations and sparsity assumptions on the associations with the outcome variable. There have been relatively few works on regression with multivariate response, and even fewer works on regression with tensor responses. One line of research treats the multivariate observations as random realizations of functional objects and conduct function-on-scalar regression (FoSR) (Reiss et al., 2010; Goldsmith & Kitago, 2016; Scheipl et al., 2015). FoSR assumes both the outcomes and the regression coefficients are smooth over time or space, and can be represented in the space expanded by a set of pre-specified basis functions. When the outcome is a multi-dimensional array, one could reform the data into a vector and assumes smoothness within local time intervals or neighborhood of voxels. Such assumptions might violate the spatial continuity embedded in the original data that is beyond Euclidean distance. Another line of research is reduced-rank multivariate regression (Chen & Huang, 2012; Chen et al., 2013, 2012; Peng et al., 2010). Reduced-rank regression take advantage of inter-relationship within the response variables to improve prediction accuracy. Our method, unlike the aforementioned two categories, model the response data as a tensor object, where the vector response can be viewed as a special case of a 1D tensor. To the best of our knowledge, few papers have tackled tensor response with dimensional flexibility except the envelope methods (Li & Zhang, 2017; Cook et al., 2010).Envelope method assumes a generalized sparsity in the data, and use BIC for selecting the latent dimension. However, the envelope methods (Li & Zhang, 2017) aims at multi-dimensional tensor outcomes not including 1D cases.

It is worth pointing out that our questions of interest are different from a class of tensor/functional regression models that quantify the relationship between a scalar outcome and a functional/tensor regressor (Goldsmith et al., 2011; Zhou et al., 2013; Goldsmith et al., 2014; Wang et al., 2014). By contrast, our proposal reverses the role by treating the tensor object as response and the vector of covariates as predictors. These two kinds of problems have different interpretations and applications. The former problem where the tensor object are predictors, focuses on assessing the effect on a clinical outcome as the biological signals (e.g. image) changes, and hence could be used for disease classification and prognosis. The latter where the tensor object is the response variable in regression, aims to identify patterns in the observed data and discover novel markers that are indicative of the disease status and demographics. Hence in this work, we focuses on the interpretability of the identified patterns and scientific conclusions are drawn from the insights gained by the estimated multi-dimensional patterns .

Secondly, biomedical images or sensor data often suffer from low signal-to-noise ratio and it is common practice to denoise the observed signals before taking it into further processing. One of the most widely used techniques in signal and image processing is Total Variation (TV) denoising. It recovers blurred signals or images by penalizing the average intensity difference of adjacent records/pixels using $L_{1}$ norm. Steidl et al. (2006) generalized the approach to penalize higher-order gradient difference and pointed out a connection between the dual formulation of TV denoising and spline smoothing of images. TV denoising can be viewed as a support vector regression problem in the discrete counterpart of the Sobolev space. Efficient algorithms have been proposed to implement the approach (see e.g. Beck & Teboulle (2009); Condat (2013); Padilla et al. (2016)). Total variation penalty has also been used in the statistics community to handle regression problems where the parameter has certain smooth structure. For example, fused lasso (Tibshirani et al., 2005) was proposed by adding penalty between adjacent sites on genomes when assessing the association between phenotypes and genotypes.

Inspired by the idea of TV de-noising, we propose a regularization term in the regression problem to encourage smoothness in the outcome space. It is important to differentiate it from the fused lasso problem, which assumes the regression coefficients are smooth. Regularization directly imposed on the outcome image and tensor space simultaneously de-noises the complex structured input data in the model fitting step without either data pre-smoothing. We will see in the experiments that this simultaneous treatment is crucial in preserving weak but consistent signals, that could be blurred out by off-the-shell de-noising algorithms.

Thirdly, the proposed method can also easily handle complex dependency structure and preserve the inherited adjacency structure of the voxels or time points from the observed data. The smoothness constraint can go beyond adjacent voxels and be customized by any adjacency structures. This could include, for example, the periodicity in longitudinal data such as seasonal effects in physical activity or temporal brain activation during a task fMRI scan.

The proposed method is mainly designed for continuous outcomes and does not require explicit assumptions on the noise distributions. For convenience of theoretical analysis, we present an oracle inequality of the estimator, and show that it is statistical consistent when the noise is i.i.d. Gaussian. To solve the aforementioned regularized regression problem, we propose an Alternating Direction Method of Multipliers (ADMM) algorithm, which, our analysis proves, is globally convergent. In both simulation studies and real data analysis, we compare our method with related work, and the proposed method exhibit strong performance especially when the sample size is small. Statistical inference on the estimated regression coefficients could be achieved using Bootstrap resampling that help the interpretability of the model in the real data analysis.

The rest of the paper is organized as follows. Section 2 describes the proposed method and introduces the model with some notations. In Section 3.2, we propose a fast and convergent algorithm to achieve the estimation. We develop the consistency of the estimator via an oracle inequality in Section 4. We demonstrate the performance in 1D and 2D settings by simulation follows in Section 5. We present the data analysis of an activity data of our proposed methods and compared to some other existing methods the result in 6. We conclude the paper with some discussion in Section 8. Detailed proofs are outlined in the web appendices.

2 Preparations

2.1 Notations and Operations

We begin with the introduction for the notations used in the paper. Throughout the paper we use lower-case and upper-case letters to represent vectors and matrices. Multidimensional array $A\in\mathbb{R}^{r_{1}\times\dots\times r_{m}}$ is called an mth-order-tensor. We will use $\|\cdot\|_{F}$ for Frobenius norm of a tensor, $\|A\|^{2}_{F}=\sum_{i_{1},\dots,i_{m}}A_{i_{1},\dots,i_{m}}^{2}$. We denote the $\ell_{1}$ norm of a as the $\ell_{1}$ norm of its vectorization: $\|A\|_{\ell_{1}}=\|\text{vec}(A)\|_{\ell_{1}}=\sum_{i_{1},\dots,i_{m}}|A_{ij}|$. $\otimes$ represents the Kronecker product between two tensors. The $\text{vec}(\cdot)$ operator denote the vectorization operation that stack the entries of a tensor into column vector, so that an entry $a_{i_{1},\dots,i_{m}}$ of $A$ maps to the $j$th entry of $\text{vec}(\cdot)$, where $j=1+\sum_{k=1}^{m}(i_{k}-1)\Pi_{k^{\prime}=1}^{k-1}r_{k}^{\prime}$. and $\text{mat}(\cdot)$ denote the the inverse operator. For an integer $n$, denote $[n]=\{1,2,\dots,n\}$.

2.2 TV de-noising

Next we review the total variation denoising objective where our method inherited the $L_{1}$ penalty for smoothing. The TV-denoising method for a 1D case,

[TABLE]

where $Y\in\mathbb{R}^{n}$ is a response vector, $\theta=(\theta_{1},\dots,\theta_{n})\in\mathbb{R}^{n}=E(Y)$ is the mean of vector of $Y$, where we often assume $\theta_{i}=\theta_{i+1}$.

In a m-dimensional setting, $Y$ is an $m^{\rm th}$ order tensor,

[TABLE]

where $D$ is an incidence matrix that we will introduce in the next session.

3 Model and Inference

3.1 Regression Model

Without loss of generality, for a more intuitive presentation, we now introduce our method for 2D tensor outcomes. Suppose we observe the covariates and 2D tensor outcomes (for example 2D images) of size $m_{1}\times m_{2}$ from $n$ subjects $(X_{i},\bar{Y}_{i})_{i=1}^{n}$, where $X_{i}\in\mathbb{R}^{p},\bar{Y}_{i}\in\mathbb{R}^{m1\times m2}$, $M=m_{1}m_{2}$ is the total number of outcome entries and $Y_{i}=\text{vec}({\bar{Y}_{i}})\in\mathbb{R}^{M}$ Consider the following model

[TABLE]

where $\Gamma\in\mathbb{R}^{p\times M}$ is the coefficient matrix of interest, and each row $\Gamma_{t\bullet}\in\mathbb{R}^{M}$ is an image, representing the coefficient map corresponds to the $t$th feature. While our theory is derived under the case of i.i.d. Gaussian noise, the method can be applied to other noise structure with real-value responses. Let $X\in\mathbb{R}^{n\times p}$ be the design matrix and $Y\in\mathbb{R}^{n\times M}$ be the matrix of outcomes. Now we can stack all the observations and reformulate the problem into the following matrix form:

[TABLE]

where \text{vec}(Y^{T})=\left(\begin{array}[]{c}Y_{1}\\ Y_{2}\\ \vdots\\ Y_{n}\\ \end{array}\right)\in\mathbb{R}^{nM} is a long vector.

Define the hat matrix as $H:=X(X^{T}X)^{-1}X^{T}$, which is consistent with that used in the classical regression analysis. We have the projection matrix to $\text{span}(X)^{\perp}$ as $I-H$. In the vectorization form, we can define the extended projection matrix to project measurements from each voxel to the linear space $\text{span}(X)^{\perp}$.

[TABLE]

Notice that for any $\Gamma\in\mathbb{R}^{p\times M}$,

[TABLE]

Now let us introduce the smoothing graph in this problem. In a 2D or higher dimensional image outcome, it is usually the intrinsic nature that the adjacent voxels are similar, except for a small number of edges. This can be summarized by defining a graph $G$, whose nodes are all the voxels in the image, and edges are the pairs of voxels which should have values that are close. For example, the grid graph is the most commonly used, where each voxel is connected to the four voxels adjacent to it in the image. The graph $G$ can also be chosen to reflect some sophisticated affinity relationship.

Let $D$ be the incidence matrix for the graph $G$ whose edges represent the smoothing affinity. To be specific, let $n$ and $m$ be the number of vertices’s and edges respectively, $D\in\mathbb{R}^{M\times m}$ where

[TABLE]

Note the orientation of the edge does not matter, since it corresponds to a negation of $D$, which does not change the $\ell_{1}$ norm. Define the extended incidence matrix $D_{v}=I_{n}\otimes D\in\mathbb{R}^{nM\times nm}$, then $\|X\Gamma D\|_{\ell_{1}}=\sum_{i=1}^{n}\|X_{i}\Gamma D\|_{\ell_{1}}$. We propose the Total Variation Regularized Tensor-on-scalar Regression (TVTR) by minimizing the following loss.

[TABLE]

where $\lambda$ is a tuning parameter controlling the tradeoffs between the linear correlation and the smoothness of the fitted image. It has been observed in Steidl et al. (2006) that total variation regularization bears strong relationship with functional analysis with splines. To be more specific, it is shown in a strictly discrete setting that thin plate splines (Duchon, 1977) of degree $m-1$ solve also a minimization problem with quadratic data term and $m$-th order TV regularization term.

3.2 Algorithm

Our goal is to estimate the multi-dimensional coefficients $\Gamma$ in the objective function (P). Although it is a convex function, the non-smoothness of the total variation regularization makes traditional algorithms suffer from slow convergence. In many applications, the optimization problem involves a large number of variables, and cannot be efficiently handled by generic optimization solvers.

Instead, we propose to use the Alternating Directional Method of Multipliers (ADMM). ADMM was originally introduced by Glowinski & Marroco (1975); Gabay & Mercier (1976) and has been frequently used in distributed settings (see Boyd et al. (2011) for an overview). The idea of using ADMM to solve total variation regularized problems were first derived by Wahlberg et al. (2012), which introduced an efficient and scalable optimization algorithm for problems satisfying a TV-separability condition. The condition requires the cost function is the sum of two terms, one that is separable in the variable blocks, and a second that is separable in the difference between consecutive variable blocks. In each iteration of our method, the first step involves separately optimizing over each variable block, which can be carried out in parallel. The second step is not separable in the variables, but can be carried out by applying ADMM to $L_{1}$ mean filtering very efficiently. This algorithm has been successfully applied to fused lasso and total variation denoising problems and exhibited competitive performance. It is also very easy to extend to distributed computing setting.

To make our problem satisfy the above TV-separability condition, we consider a reformulation of the original problem. Define $\theta=X\Gamma\in\mathbb{R}^{nM}$, the problem is equivalent to,

[TABLE]

The solution of (P) and that of (CP) have the following relationship.

Lemma 1**.**

Let $\hat{\theta}$ be the optimal solution of (CP), and $\hat{\Gamma}$ be the optimal solution of (P). If rank($X)=p$, then $\hat{\Gamma}=(X^{T}X)^{-1}X^{T}\text{mat}(\hat{\theta})_{n\times M}$.

Now to solve (CP), we introduce two auxiliary variables and write the original problem in the following form.

[TABLE]

where $\delta(\cdot)$ is the characteristic function, which takes 0 if the condition in the parenthesis is satisfied and infinity otherwise. The objective function is separable for $\theta,\eta$ and $\mu$, we solve it with the Alternating Direction Method of Multiplier (ADMM), which is summarized in Algorithm 1.

In Algorithm 1, line 4-6 are the primal variable updates and line 7-8 are the dual variable updates. It is common practice to set $\rho=1$. The update of $\theta$ in line 4 and projection step in line 5 both have analytical form and only consist of matrix multiplication. The computational bottleneck is the subproblem in line 6, which can be recognized as a graph fused lasso (GFL) problem that takes observation $\theta_{k+1}+V_{k}$, regularization parameter $\lambda/\rho$, and graph incidence matrix $D_{v}^{T}$. GFL is a generalization of fused lasso. It penalizes the first differences of the signal across edges and can be efficiently solved with many off-the-shell GFL solvers. Below we highlight several solvers for GFL, and compare their efficiency under different assumptions.

When the graph is a chain graph, the GFL reduces to a 1D fused lasso problem, and can be solved in $O(n)$ time, by, for example, the “taut string” algorithm derived by Davies & Kovac (2001) and the dynamic programming based algorithm by Johnson (2013). These methods are fast in practice. For 2D cases, Kolmogorov et al. (2016) generalizes the dynamic programming idea to solve the fused lasso problem on a tree in $O(n\log n)$ time. Barbero & Sra (2014) uses operator splitting techniques like Douglas-Rachford splitting and extends fast 1D fused lasso optimizers to work over grid graphs. Over general graphs structure, numerous algorithms have been proposed in recent years: Chambolle & Darbon (2009) described a direct algorithm based on a reduction to parametric max flow programming; Chambolle & Pock (2011) described what can be seen as a kind of preconditioned ADMM-style algorithm; Kovac & Smith (2011) described an active set approach; Landrieu & Obozinski (2016) derived a method based on graph cuts. In our experiments, we use the one proposed in Tansey & Scott (2015) which leverages fast 1D fused lasso solvers in an ADMM decomposition over trails of the graph. One can choose the graph fused lasso solver that best suits the target graph structure.

Despite being a popular algorithm for many constrained optimization problems, the convergence of ADMM does not have a simple and unified answer. We present the convergence of Algorithm 1 in the following theorem.

Theorem 1**.**

Algorithm 1 converges linearly to the unique global optimum of problem (CP).

To obtain the result, we first convert problem CP into a two-block problem; then apply recent result in general two-block problems Nishihara et al. (2015) to complete the proof. The detailed proof is deferred to supplementary materials.

3.3 Software

We provide the open source python code for implementation of Algorithm 1 on https://github.com/summeryingliu/imagereg.

4 Theoretical Results

In this section, we analyze the statistical property of the estimator given by (P). We present the result via an oracle inequality on the prediction error of the regression. Before we state the theoretical result, we first introduce two quantities which are widely used in the analysis of sparse recovery problems.

Definition 1** (Compatibility factor).**

Let $D\in\mathbb{R}^{M\times m}$ be an incidence matrix. The compatibility factor of $D$ for a set $\emptyset\subsetneq T\subset[m]$ is defined as

[TABLE]

Compatibility factor gets its name based on the idea that, on the subset of edges indicated by set $T$, we require the $\ell_{1}$-norm and the $\ell_{2}$-norm to be somehow compatible. Compared with other conditions used to derive sparsity oracle inequalities, such as restricted eigenvalue conditions or irrepresentable conditions, the compatibility factor greater than 0 is shown to be weaker by Van De Geer et al. (2009). More discussion about the relationship between different conditions can be found in Van De Geer et al. (2009). For graphs with bounded degree, it is shown that the compatibility condition is always satisfied.

Lemma 2** (Hütter & Rigollet (2016)).**

Let $D$ be the incidence matrix of a graph $G$ with maximal degree $d$ and $\emptyset\neq T\subset E$. Then, $\kappa_{T}\geq\frac{1}{2\min\{\sqrt{d},\sqrt{|T|}\}}$.

Definition 2** (Inverse scaling factor).**

The inverse scaling factor of an incidence matrix $D$ is defined as $\rho(D):=\max_{j\in[m]}\|s_{j}\|$, where $S=(D^{T})^{\dagger}=[s_{1}^{T},\cdots,s_{m}^{T}]^{T}$ is the pseudo inverse of $D^{T}$.

By design of $D_{v}$, it is clear that $\rho(D_{v})=\rho(D)$. Now we present the main result.

Theorem 2** (Oracle Inequality for Projected TV Regression).**

Under model (1) with $\epsilon_{i}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\mathcal{N}(0,\sigma^{2}I_{M})$, define $\theta^{*}=\text{vec}(\Gamma^{*T}X^{T})$, $\hat{\theta}$ is the solution for (CP). For any $\delta>0$, if $\lambda=\rho\sigma\sqrt{\log(mnM/\delta)}$, then with probability at least $1-\delta$,

[TABLE]

The proof of Theorem 2 is inspired by that in Hütter & Rigollet (2016) and is deferred to Section B in supplementary materials.

Remark 1**.**

*The length of $\theta$ scales as $nM$. Hence we care about the mean recovery error $\frac{1}{nM}\|\theta^{*}-\hat{\theta}\|^{2}$, which converges at rate $O\left(\frac{\log(mnM)}{nM}\right)$. The upper bound exhibits the trade-off between two quantities: the number of changing point $|T|$ and the total variation of the “smooth” part $\|(D\theta)_{T^{c}}\|_{\ell_{1}}$. To be more specific, given the total variation of the entire data fixed, when the model is piecewise smooth but have drastic change at the changing point, $|T|$ will dominate $\|(D\theta)_{T^{c}}\|_{\ell_{1}}$; on the other hand, if the data has few changing point but fluctuate a lot in each piece, the total variation of $\theta$ in $T^{c}$ will be large compared to $|T|$. *

If the graph has bounded degree, by Lemma 2 we immediately have the following corollary.

Corollary 1**.**

If the maximal degree of the penalty graph $G$ is $d$, $\lambda=\rho\sigma\sqrt{\log(mnM/\delta)}$, then with probability at least $1-\delta$,

[TABLE]

In particular, if the features has isotropic covariance matrix, then we will have the following bound on the parameter estimation.

Corollary 2**.**

If $\frac{1}{n}X^{T}X=I_{p}$, and $\lambda=\rho\sigma\sqrt{\log(mnM/\delta)}$, then with probability at least $1-\delta$,

[TABLE]

Remark 2**.**

The condition in Corollary 2 can be achieved by normalizing the input design matrix. When the covariance matrix for the features is not isotropic, the error term should be represented in Mahalanobis distance instead, that is, the error along different axis needs to be reweighted by the variance in that direction.

5 Simulation Studies

In this section, we present the simulation results on recovering signal $\Gamma$ with various methods. Since the goal is to accurately estimate $\Gamma$, we use the mean devision (square root of the mean squared error) of the coefficients as the performance metric, which is defined as $\frac{1}{\sqrt{Mp}}\|\hat{\Gamma}-\Gamma^{*}\|_{F}$. The standard deviations of the metric over replicates are reported in the parenthesis.

We conduct experiments on 1D and 2D synthetic data. For 1D data, we compare our method with the state-of-the-art function-on-scalar regression with B-spline and Fourier basis, which is implemented in R package refund; The setting and results are summarized in Section 5.1. For 2D data, the function-on-scalar method is applied on the vectorized image. We also compare with Parsimonious Tensor Response Regression (Li & Zhang, 2017), and some two-step procedures which we discuss in Section 5.2.

5.1 One Dimensional Simulation

We generate 2 simulation scenarios, each with 200 replications. For each scenario, we present results for sample sizes of 25, 50 and 100. In the first setting we generate model that satisfies the assumptions of function-on-scalar regression, i.e., the signals are sparse if expanded in Fourier basis. The second setting is generated not following this model assumption.

For the first setting, the predictors are generated as following: $(X_{1},X_{2})$ are drawn from categorical distribution that takes value $(1,0)$ with probability $1/4$, $(0,1)$ with probability $1/4$, and $(0,0)$ with probability $1/2$. $X_{3}$ is drawn from standard normal distribution. The true signals are from Fourier basis functions, with index $t$ ranges from $1$ to $200$, and

[TABLE]

We compared our proposed method with Functional-on-Scalar regression methods (FoSR) using B-spline and Fourier basis functions with $L_{2}$ penalty, where the tuning parameter is selected through cross validation from a grid of $(2^{-5},2^{-3},2^{-1},2,2^{3},2^{5})$. It is noticeable that scenario 1 is in favor of the FoSR methods since it is generated from the true model assumption with Fourier basis. The Parsimonious Tensor Response Regression (Li & Zhang, 2017) is proposed only for tensor regression problem and not intended for 1D function-on-scalar regression problem.

Table 1 presents the simulation results for scenario 1. When sample size is $25$, the proposed method shows a better performance than all the competitors in terms of smaller mean and standard devision. Our method performs similarly with the FoSR methods when sample size is 50. With large sample size 100, the FoSR method with the correct basis has the best performance. The Fourier basis function performs slightly better then b-spline basis, however the difference is not phenomenal.

The second simulation setting showed the advantage of our methods where the function on scalar regression model are misspecified. It also demonstrates the advantage of our proposed method in the flexibility for defining arbitrary adjacency structure with self defined incidence matrix. We generate piecewise constant signals, and in addition the signals for indexes 1 to 100 are identical with signals in 101:200. This setting is motivated by the scenario in time series observations, where one may have observations of the same pattern in repeated time periods. For example, periodicity that correspond to day of the week or season of the year might be observed in subjects’ levels of physical activity. Similar brain activation time series might also manifest during a task fMRI experiment when the stimulations were given repeatedly. We used this 1D simulation to demonstrate the performance of incorporating some prior knowledge compared with not using this knowledge.

We consider two regimes in defining the edges. For vanilla TVTR, we only consider edges connecting adjacent observations: $E=\{(i,i+1),i=1,\dots,199\}$. For the periodic_TVTR, we add in 100 more edges that connects the $i$th and $100+i$th observations, which encourages the periodicity of the signal. Here the edge set is $E=\{(i,i+1),i=1,\dots,199\}\cup\{(i,i+100),i=1,\dots,100\}$. The true signal is generated by the following model:

[TABLE]

The results are shown in Table 1. We can see TVTR with or without the added edges outperforms the other methods. As sample size increases all the methods improves with smaller mean and SD for the mean deviation. By adding the additional edges, our proposed method shows better performance, and the advantage is more phenomenal with small sample size. The proposed method with added edges to incorporate prior knowledge and encourage similar patterns reduced the mean deviation to less than half for the one estimated without added edges.

5.2 Two Dimensional Simulation

In this section, we present the experimental results for 2D synthetic data. All tuning parameters in the proposed method are chosen by 4-fold cross validation with a grid search in $[0.1,0.25,0.5,1,1.5,2,3]$. Some methods such as the one proposed in Reiss et al. (2010) fail to adapt to the scale we consider in this simulation due to the large memory cost. We demonstrate the performance of the function-on-scalar framework using penalized splines (Goldsmith & Kitago, 2016), with an efficient variational Bayes implementation (bayes_fosr in R package refund). For envelope method proposed by Li & Zhang (2017), we use the matlab package provided by the authors. We also consider two two-stage methods: TV-smooth before and after the regression, where the regression model is the voxel-wise ordinary linear regression.

All methods used in the simulation are summarized as following.

• •

TVTR: Total Variation regularized Tensor-on-scalar Regression by Algorithm 1.

• •

Envelope: Parsimonious Tensor Response Regression by Li & Zhang (2017);

• •

FoSR: A variational Bayes implementation for penalized splines.

• •

TV_OLS : Each image is denoised individually by total variation regularization, and the estimator is achieved by conducting a voxel-wise OLS on the denoised images.

• •

OLS_TV: A TV denoised estimator from voxel-wise OLS regression.

We consider two different settings in the simulation.

In the first setting, we simulate 3 predictors as indicators of 3 disease groups and an integer value predictor as age. We generated 3 patient-level variables: $X_{1}$ and $X_{2}$ are disease group dummy variables, $X_{1}=0,X_{2}=0$ denote the control group, $X_{1}=1,X_{2}=0$ represent the disease group 1, and $X_{1}=0,X_{2}=1$ represent the disease group 2. A integer value variable $X_{3}$ is generated uniformly over integers from $56$ to $75$. The image size is $m_{1}=40$ by $m_{2}=40$, the true coefficient maps ($\Gamma_{1},\Gamma_{2},\Gamma_{3}$) are block-wise constant (as shown in top left panel of Figures 1). And the outcome is generated by

[TABLE]

where $\epsilon\in\mathbb{R}^{m_{1}\times m_{2}}$, and $\epsilon_{ij}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\mathcal{N}(0,2)$.

In the second setting, two binary variables $X_{1}$ and $X_{2}$ are generated same as in setting 1. The coefficient map of $X_{1}$ is generated to have active regions with different sizes: $1$, $4$, $25$ pixels, and the true coefficients are $2$, $1.5$ and $1$, respectively. Please refer to Figure 2 for the actual arrangements of these block-wise true signals.

In Figure 1 and 2, we present one example of the estimation for various methods with training sample size of 100. The FoSR is designed for one-dimensional functionals, so we implemented it with the input as the vectorization of the image, this is the main reason why it only maintains the continuity along horizontal directions and fails to detect other types of structure. There are no off-the-shelf implementation of 2D FoSR as far as we know.

According to Table 2, with increasing sample sizes, all methods perform better, but our method outperforms the others in all scenarios. Our proposed method outperformed the other methods except for setting 1 with sample size 100. Envelope method demonstrates better performance with the large sample size for setting 1 which also has larger ’block size’. Envelope method is outperformed by our proposed method and TV_OLS (regression after TV denoising) for setting 2 with all three sample sizes, when the block has various sizes.

According to Figures 1 and 2 that our proposed method outperforms the other methods in recovering the true structure of the parameters, and achieves a cleaner cut on the boundary. Envelope method is able to capture and display the main pattern, however does not encourage sparsity of the edges, thus the solution is not “smooth” anywhere. TV_OLS (regression after TV denoising) estimated coefficient maps are more ’vague’ at the boundaries, since the the smoothing step is done separately for individual images, the estimated change points/edges will not be consistent across the smoothed images. For setting 2, TV_OLS performs well when sample size is small, but it does not improve much when sample size increases, since the TV denoising for single images can only recover piecewise constant signals for the moderate-size blocks. The signal from small regions will be smoothed out regardless of increasing the sample size because the denoising step only involves data from one image. In comparison, our proposed method gains much performance improvement with increasing sample size, as the signal is strengthened by more observations. By simultaneously conducting Image-on-scalar Regression and Total Variation Smoothing, our proposed method distinguishes signal on small regions (which is shared in all observations) versus random noise (which has no consistent behavior across different images), while competing methods fail to do so. OLS_TV performed well for setting 1 and has a “smoother” estimation for the coefficient map than our method. But it fails to detect the small regions in setting 2 and performed worse than TV_OLS.

To summarize, TVTR is a universal method for recovering piecewise constant signals in small and large areas. It shows evident advantages of recovering signal in small areas compared with the two-stage methods and envelope method.

6 Application to Activity Data

As discussed in the introduction, for the NIMH family study of spectrum disorders, we are interested in removing noises from the minute-level physical activity measures, and meanwhile accounting for heterogeneity in physical activity due to differences in demographic characteristics as well as health-related factors such as disease groups. Hence the outcome measures in our example are one-dimensional time series of activity intensities assessed over two weeks for each participant. Given the skewness of the observed activity counts (Shou et al., 2017), we conducted a log (counts+1) transformation beforehand and pre-smoothed the outcome data using moving average filter over a $60$-minute window, that is, averaged across an hour. As shown in Figure 3, the raw activity counts were fairly noisy, but we assume that the underlying activity status for a particular participant (e.g., being inactive, moderately or vigorously active) would transit smoothly and continuously over time. Therefore such data might benefit from total variation smoothing that would restrict dramatic fluctuations from minute to minute. The covariates that we accounted in the model include diagnosis, age, gender, body mass index (BMI) and day of week. The associations between time-varying physical activity intensities and each of these scalar variables were estimated using our proposed model. To make the results more clinically interpretable, we categorized BMI into underweight ($<18.5$), normal ($18.5\leq$BMI$\leq 25$), overweight ($25<$BMI$\leq 30$) and obese (BMI$>30$), with the normal participants being the reference group. Similarly, age was also stratified into four categories: adolescence (under 18), adulthood (18 to 40), middle age (40 to 60) and elderly ($\geq 60$). There were 5 diagnostic groups, including healthy control, type I bipolar (BPI), type II bipolar (BPII), major depressive disorder (MDD), and other disorders. The sample included 339 subjects with 14 days of data and our final analysis was based on 301 subjects with non-missing covariates. In total, there were 4214 daily activity curves and the observations were balanced for the 7 days of the week (Monday through Sunday).

The estimated time-varying effects for each covariate are shown as the red curve in Figure 4. The point-wise confidence bands in red were obtained as the upper and lower $2.5\%$ quantiles over $100$ bootstrap samples. We observed that compared to the middle-age group, seniors tend to have an overall lower activity intensity throughout the day, especially from late afternoon to night. Both adolescent and adult groups have shown higher night-time activity and lower early morning activity, but the difference is more prominent between adolescent and middle-aged group. Participants who were obese were observed to have lower day-time activity and higher nighttime activity as compared to subjects with normal weights. The overweight subjects showed a similar trend but with a wider confidence band and smaller magnitude, so the trend were not significant. The underweight subjects were observed to have a reversed trend as compared to normal. In terms of diagnostic groups, BPI patients showed lower activity intensity later of the day (after 12pm) as compared with healthy controls, which is consistent with the findings reported in (Shou et al., 2017). Females were also observed to have higher average activity intensities during day time and lower intensities at night. Interestingly, when investigating the activity patterns over 7 days a week with Sunday as the reference, we independently observed similar patterns across weekdays (Monday to Friday), where the morning activity intensities are much higher than Sunday, but no difference were shown for the rest of the day. Note that on Friday, we observed a rising level of activity towards midnight, indicating that subjects were engaged in more late-night activities on Friday. Saturday demonstrated a non-differentiable pattern compared with Sunday. But a similar rising pattern were also observed towards midnight.

To further evaluate our proposed method, we equally split the two weeks of activity data per subject into the training and testing datasets, where each dataset contained about 2107 days of measures respectively from the same set of subjects. We compared our proposed model with several competing methods including the minute-wise OLS regression, the function-on-scalar regression (FoSR) with B-spline basis functions, FoSR with $L_{2}$ penalty and the envelope method. The models were fit based on the training data and evaluated on the testing dataset by calculating the corresponding mean square errors (MSE) between the predicted activity counts and the observed counts. The prediction error (MSE) of our proposed method is 2.76, which is the smaller than both the minute-wise OLS (2.79), the FoSR with B-splines (2.80), as well as FoSR with $L_{2}$ penalty (2.79). The envelope method that had similar performance in the simulation studies demonstrated large MSE and is not suitable for the 1D application. As a reference, the total sample variance in the testing set is 2.90 which indicates that the activity data were very noisy. Although neither of the methods did a perfect job predicting activity counts, our proposed method were able to obtain interpretable estimates of for the association between activity and covariates and meanwhile demonstrate some improvement in prediction than the other methods.

7 Application to 3D Brain Imaging Analysis

We now demonstrate the proposed method in analyzing 3D brain imaging. Our data came from a public data source of attention hyperactivity disorder (ADHD), the ADHD-200 sample, which was an initiative effort to understand neuro basis of the disease by ADHD consortium led by Dr. Michael Milham. The dataset contains structural and BOLD functional magnetic resonance imaging (MRI) scans from 776 subjects, including 491 typically developing controls (TDs) and 285 ADHD children.

We focused on the gray matter voxel-based morphometry (VBM) maps preprocessed via the Burner pipeline using SPM8, downloaded from the Neuro Bureau (https://www.nitrc.org/plugins/mwiki/index.php?title=neurobureau:BurnerPipeline). The intensities of the images quantifies the normalized gray matter concentration changes at specific brain locations. The same dataset was used in the tensor envelope based regression method by (Li & Zhang, 2017), which is one of the state-of-the-art methods for the 3D regression problem. Thus we compare the results of our proposed methods with the tensor envelope based regression and OLS estimates presented in (Li & Zhang, 2017). The images are of dimension 121$\times$ 145$\times$ 121. Following (Li & Zhang, 2017), we downsized the images to 30$\times$ 36$\times$ 30 using TensorReg. The covariates in the regression model include ADHD, age, gender, handedness. We also included the site indicators in the multivariate analysis to adjust for the systematic site effects. The estimated regression coefficients were mapped back to the 3D population-average template of the Burner VBM images. We also registered Harvard-Oxford Atlas (Bohland et al., 2009) to the same template for references of the cortical regions.

The computation cost of 3D data raised a challenge to many algorithms including ours. This specific example consists images from 770 subjects, each consists 32400 nodes and 94140 edges for a grid smooth structure. Algorithm 1 offers an iterative procedure of solving a TV denoising problem. In this real data example and our simulation studies, the algorithm converges in 10-20 iterations, when the parameters are initialized with random noise following normal distribution. The most computational expensive step of Algorithm 1 is line 6, where one is solving a TV denoising problem for a graph with $32400\times 770$ nodes and $94140\times 770$ edges. Solving TV denoising in a single step with this scale will cost memory failure on a personal computer. Fortunately, this task can be distributively solved. Since there are no edges across subjects, one can solve it by divide and conquer. i.e. within each iteration, at line 6, one could divide the samples to several subgroups, solve for the $\mu^{(k+1)}$ for each subgroup, and combine the results together. The combined estimates of $\mu^{(k+1)}$ is equivalent to the solutions from computing the combined TV denoising problem. For this dataset, we solved line 6 by ‘divide and conquer’ with mini-batch size of 8 samples, each with the tuning parameter $\lambda=0.05$. The computation was conducted on a Windows 10 desktop computer with Intel(R) Core(TM) i7-7820U CPU@3.60GHz processor (16 cores), 32.00 GB installed memory(RAM). It takes approximately 10 hours for obtaining the estimates under a chosen tuning parameter. Admittedly, the computational challenge caused problem in conducting cross validation for tuning parameters and for obtaining the bootstrap confidence intervals on a single machine. However, the algorithm is highly distributed and bootstrap would be practical using parallel computing resources.

Here we demonstrate the estimated regression coefficients in 3D representation based on our proposed method. The coefficient maps correspond to the association between voxel-wise gray matter morphometry with ADHD, age, gender and handedness (1 for left-handedness and 0 for right-handedness). We only demonstrated the areas where the coefficients’ magnitudes are in the 10%, 5%, 2.5% and 0.5% tails across the whole brain.

As shown in the Figure 5, we observe that after adjusting for age, gender, handedness and site effects, ADHD showed the largest reduction in the gray matter volume in the frontal lobe such as the anterior and posterior cingulate gyrus and frontal operculum. Such regions have been known to be associated with sensory and motor responses, as well as cognitive and emotional functions. Significant cortical thinning in anterior cingulate gyrus has been reported among children with ADHD (van Rooij et al., 2015). These findings were novel in our analysis methods, while in the results of the envelope method (Li & Zhang, 2017) the effect of frontal lope were not captured. Instead, their methods identified superior temporal gyrus, pyramid and uvula in cerebellum. Additionally, we observed that cortical thickness was increased in frontal operculum with older children, which is linked with better task controls(Higo et al., 2011). Females are shown to have thicker temporal fusiform cortex and frontal medial cortex. Interestingly, although the magnitude of the handedness effects are small, we observed an asymmetric association where the left-handed children tend to have greater gray matter thickness on the left hemisphere of brain in regions such as posterior parahippocampal gyrus, lingual gyrus and temporal occipital fusiform cortex. Such findings are consistent with several previous studies (Cuzzocreo et al., 2009).

8 Conclusion and Discussion

We propose a new method for tensor-on-scalar regression with a TV regularization term to encourage sparsity of adjacent outcome estimates. The proposed method can be universally applied to tensor-on-scaler regression problem with any dimensions. The proposed algorithm is showed to converge, it is both scalable and consistent. Our method also has the appealing new feature of incorporating prior knowledge for the smoothness by user-defined adjacency matrix. A extension of application of the method is in brain connectivity studies, where the smoothness can be defined as the connected area.

Although the method focuses on providing point estimation. In medical translational researches, in addition to signal recovering accuracy, proper inference of variability of the estimation is desirable. Thus we demonstrate bootstrap confidence interval can be computed to help interpreting the result in the real data example of activity data analysis. However, as dimension of the outcome increases from 1D to 3D, the computational burden increases tremendously that bootstrap CI is might be too costly to compute without utilizing powerful parallel computing cluster, which is a limitation of our current work.

Our current method only estimate the mean association not taking into account correlation of tensor-outcomes. It remains for future research how to more efficiently estimate the mean, with multiple outcome tensors per subject by adding random effect in the model.

Acknowlegements

The authors would like to thank Drs. Xin Zhang and Lexin Li for providing codes and guidance of their envelope tensor regression algorithms. We also would like to acknowledge Junchi Li and Michael Daniels for discussion and advices to improve the quality of the work, and Jinshan Zeng for his comments on the convergence properties of the algorithm. HS’s research was partially supported by the Intramural Research Program of the NIMH. The NIMH family study was conducted under the support of by the Intramural Research Program (Merikangas; grant number Z-01-MH002804) and under clinical protocol 03-M-0211 (NCT00071786).

Disclaimer

The views and opinions expressed in this article are those of the authors and should not be construed to represent the views of any of the sponsoring organizations, agencies or US Government.

SUPPLEMENTARY MATERIAL

Appendix A Proofs in Section 3

A.1 Proof of Lemma 1

Proof of Lemma 1.

Recall that $H_{X}$ is the projection matrix onto $\text{span}(X)^{\perp}$, therefore $\theta\in\text{span}(X)$ is equivalent to $H_{X}\theta=0$. By a variable transformation $\theta=X\Gamma$, solving (P) is equivalent to solving (CP), and $\hat{\theta}=X\hat{\Gamma}$. When $X^{T}X$ is invertible, we solve $\hat{\Gamma}=(X^{T}X)^{-1}X^{T}\hat{\theta}$. ∎

A.2 Proof of Theorem 1

Proof of Theorem 1.

We rewrite problem (2) as

[TABLE]

where

[TABLE]

Algorithm 1 is presented in the scaled form of ADMM. As is shown in Boyd et al. (2011), it is equivalent to the unscaled form as below. Denote the Lagrangian as

[TABLE]

The unscaled form is

[TABLE]

We first convert the three-block problem into a two-block problem with the technique in chen2016direct.

By first-order optimality condition in iteration $k$, we have

[TABLE]

Note that $A_{2}^{T}A_{3}=0$, we have

[TABLE]

which is also the first order optimality condition for the regime:

[TABLE]

Clearly, this is a specific application of the two-block ADMM by regarding $(\eta,\mu)$ as one variable, $B:=[A_{2},A_{3}]$ as one matrix, and $g(\eta,\mu):=g_{1}(\eta)+g_{2}(\mu)$ as one function. Existing convergence results for the two-block ADMM thus hold for our case.

Now we note that $f$ is Lipschitz differentiable and strongly convex, both $f(\theta)$ and $g(\eta,\mu)$ are closed (their sublevel sets are closed) and proper (they neither take on the value $-\infty$ nor are they uniformly equal to $\infty$)). Also $A$ and $B$ both have full column rank. Therefore by Theorem 7 in Nishihara et al. (2015), Algorithm 1 converges to the unique global solution of (P) linearly. ∎

Appendix B Proofs in Section 4

B.1 Proof of Theorem 2

Proof of Theorem 2.

Define $H_{X}=I-H_{v}$ is the projection matrix projecting each voxel to $\text{span}(X)$. Let $\hat{\theta}\in\mathbb{R}^{nM}$ be the optimal solution for the following constrained optimization problem.

[TABLE]

When $\text{rank}(X)=p$, by Lemma 1, $\hat{\Gamma}=(X^{T}X)^{-1}X^{T}\text{mat}(\theta)_{n\times M}$. We first prove an oracle inequality for $\hat{\theta}$. By the KKT condition, $\exists z\in\text{sign}(D_{v}^{T}\hat{\theta}),\alpha\in\mathbb{R}^{nM}$ such that

[TABLE]

Multiplying $H_{X}$ on the left of Eq. (3), we have $2H_{X}(\hat{\theta}-y)+\lambda H_{X}D_{v}^{T}z=0$. Equivalently,

[TABLE]

Note that $H_{X}$ is the projection matrix to the column space of $\text{span}(X)$, it suffices to consider $\forall\bar{\theta}=\text{vec}(\bar{\Gamma}^{T}X^{T})$. By definition of the sign operator, the following holds:

[TABLE]

Subtracting the former from the latter, and replacing $\text{vec}(Y^{T})$ with $\theta^{*}+\epsilon$, we get

[TABLE]

Note $(\bar{\theta}-\hat{\theta})^{T}(\theta^{*}-\hat{\theta})=\frac{1}{4}\left(\|\bar{\theta}-\hat{\theta}\|^{2}+\|\theta^{*}-\hat{\theta}\|^{2}-\|\bar{\theta}-\theta^{*}\|^{2}\right)$,

[TABLE]

To bound $(\hat{\theta}-\bar{\theta})^{T}\epsilon$, note $DD^{T}$ is the graph Laplacian, therefore when the graph is connected, we have $\text{ker}(D^{T})=\text{ker}(DD^{T})=\text{span}\{1_{M}\}$. Define $D^{\dagger}$ be the pseudo inverse of $D$, then $I-(D^{\dagger})^{T}D^{T}$ is the projection matrix onto ker($D^{T}$),

[TABLE]

In view of the fact that $\epsilon_{i}\sim\mathcal{N}(0,\sigma^{2}I_{M})$ and $(I-(D^{\dagger})^{T}D^{T})$ being a projection matrix to a one dimensional space, by the tail bound for Gaussian random variables, $\forall i\in[n],\forall\delta>0$,

[TABLE]

For the second part, by the maximal inequality for Gaussian random variables (massart2007concentration Thm 3.12), and the variance of elements of $(D^{\dagger})^{T})\epsilon_{i}$ is upper bounded by $\rho^{2}\sigma^{2}$,

[TABLE]

Applying union bound, with probability at least $1-\delta$,

[TABLE]

By the triangle inequality,

[TABLE]

Hence

[TABLE]

By Definition 1 , $\|D^{T}(\hat{\theta}_{i}-\bar{\theta}_{i})_{T}\|_{\ell_{1}}\leq\kappa_{T}^{-1}\sqrt{|T|}\|\hat{\theta}_{i}-\bar{\theta}_{i}\|.$ We now plug above and (5) back to (4), and take $\lambda=\rho\sigma\sqrt{\log(mnM/\delta)}$, then with probability at least $1-c_{7}n^{-1}$,

[TABLE]

Use Young’s inequality,

[TABLE]

Canceling out $\|\hat{\theta}-\bar{\theta}\|^{2}$ on both sides of (6),

[TABLE]

Taking infimum on the right and plugging in $\lambda$ we have

[TABLE]

∎

B.2 Proof of Corollary 2

Corollary 2 can be proved by combining Lemma 1 and Theorem 2.

Proof of Corollary 2.

When $\frac{1}{n}X^{T}X=I_{p}$, we will be able to control the error in $\Gamma$, notice that $\|\text{vec}(A)\|_{2}^{2}=\|A\|_{F}^{2}$, we have from (B.1) that

[TABLE]

∎