SILVar: Single Index Latent Variable Models

TL;DR

SILVar introduces a semi-parametric, non-linear regression model that jointly estimates system non-linearities, direct interactions, and unmeasured effects using a convex optimization framework, applicable to complex networked systems.

Contribution

It proposes a novel semi-parametric model with a regularized empirical risk minimization approach, connecting it to existing paradigms like Robust PCA and VAR, and provides an efficient learning algorithm.

Findings

Effective in modeling complex networked systems with latent variables

Demonstrates strong performance on simulated and real data

Bridges gaps between different modeling paradigms

Abstract

A semi-parametric, non-linear regression model in the presence of latent variables is introduced. These latent variables can correspond to unmodeled phenomena or unmeasured agents in a complex networked system. This new formulation allows joint estimation of certain non-linearities in the system, the direct interactions between measured variables, and the effects of unmodeled elements on the observed system. The particular form of the model adopted is justified, and learning is posed as a regularized empirical risk minimization. This leads to classes of structured convex optimization problems with a "sparse plus low-rank" flavor. Relations between the proposed model and several common model paradigms, such as those of Robust Principal Component Analysis (PCA) and Vector Autoregression (VAR), are established. Particularly in the VAR setting, the low-rank contributions can come from broad trends exhibited in the time series. Details of the algorithm for learning the model are presented. Experiments demonstrate the performance of the model and the estimation algorithm on simulated and real data.

Equations

$$D_{\varphi }(\mathbf{y}\Vert \mathbf{x})=\varphi (\mathbf{y})-\varphi (\mathbf{x})-\nabla \varphi (\mathbf{x}{)}^{\top }(\mathbf{y}-\mathbf{x}).$$

$${\varphi }_{\ast }(\mathbf{x})\stackrel{\Delta }{=}\underset{\mathbf{y}}{\sup }{\mathbf{y}}^{\top }\mathbf{x}-\varphi (\mathbf{y}).$$

$${\varphi }_{\ast }(x)=x(\nabla \varphi {)}^{-1}(x)-\varphi ((\nabla \varphi {)}^{-1}(x))$$

$$(\nabla \varphi {)}^{-1}=\nabla {\varphi }_{\ast }({\varphi }_{\ast }{)}_{\ast }=\varphi$$

$$\mathbb{E}[y_i|{\mathbf{x}}_i]=g\left({\mathbf{a}}^{\top }{\mathbf{x}}_i\right),$$

$$P(y_i|{\mathbf{x}}_i)=\exp \left\{-D_{\varphi }\left(y_i\Vert g\left({\mathbf{a}}^{\top }{\mathbf{x}}_i\right)\right)\right\}$$

$$\hat{\mathbf{a}}$$

$$=\underset{\mathbf{a}}{\mathrm{argmin}}\sum _{i=1}^n{D}_{\varphi }\left(y_i\Vert g\left({\mathbf{a}}^{\top }{\mathbf{x}}_i\right)\right)$$

$$=\underset{\mathbf{a}}{\mathrm{argmin}}\sum _{i=1}^n\begin{array}{rl}& [\varphi \left(y_i\right)-\varphi \left(g\left({\mathbf{a}}^{\top }{\mathbf{x}}_i\right)\right)\\ & -\nabla \varphi \left(g\left({\mathbf{a}}^{\top }{\mathbf{x}}_i\right)\right)\left(y_i-g\left({\mathbf{a}}^{\top }{\mathbf{x}}_i\right)\right)]\end{array}$$

$$\stackrel{(a)}{=}\underset{\mathbf{a}}{\mathrm{argmin}}\sum _{i=1}^n\begin{array}{rl}& [G_{\ast }\left(y_i\right)-y_i\left({\mathbf{a}}^{\top }{\mathbf{x}}_i\right)-\varphi \left(g\left({\mathbf{a}}^{\top }{\mathbf{x}}_i\right)\right)\\ & +\left({\mathbf{a}}^{\top }{\mathbf{x}}_i\right)g\left({\mathbf{a}}^{\top }{\mathbf{x}}_i\right)]\end{array}$$

$$\stackrel{(b)}{=}\underset{\mathbf{a}}{\mathrm{argmin}}\sum _{i=1}^n\left[G_{\ast }\left(y_i\right)+G\left({\mathbf{a}}^{\top }{\mathbf{x}}_i\right)-y_i\left({\mathbf{a}}^{\top }{\mathbf{x}}_i\right)\right]$$

$$=\underset{\mathbf{a}}{\mathrm{argmin}}{\hat{F}}_1\left(\mathbf{y},\mathbf{X},g,\mathbf{a}\right)$$

$${\hat{F}}_1(\mathbf{y},\mathbf{X},g,\mathbf{a})\stackrel{\Delta }{=}\frac{1}{n}\sum _{i=1}^n\left[G_{\ast }(y_i)+G\left({\mathbf{a}}^{\top }{\mathbf{x}}_i\right)-y_i\left({\mathbf{a}}^{\top }{\mathbf{x}}_i\right)\right]$$

$$(\hat{g},\hat{\mathbf{a}})$$

$$\text{s.t. }g=\nabla G\in {\mathcal{C}}^1.$$

$$y_i=g(x_i)+w_i,$$

$$\hat{\mathbf{g}}\stackrel{\Delta }{=}\text{LMR}(\mathbf{y},\mathbf{x})$$

$$\text{s.t. }0\le g\left(x_{[j+1]}\right)-g\left(x_{[j]}\right)\le x_{[j+1]}-x_{[j]}$$

$$\text{for }j=1,\dots ,n-1.$$

$${\nabla }_{\mathbf{a}}{F}_1$$

$$\mathbb{E}\left[y_{ji}|{\mathbf{x}}_i\right]$$

$$\Rightarrow \mathbb{E}\left[{\mathbf{y}}_i|{\mathbf{x}}_i\right]$$

$$\left(\overline{g},\overline{\mathbf{A}},\overline{\mathbf{B}}\right)=\begin{array}{rl}& \underset{g,\mathbf{A},\mathbf{B}}{\mathrm{argmin}}{F}_2\left({\mathbf{y}}_i,{\mathbf{x}}_i,{\mathbf{z}}_i,g,\mathbf{A},\mathbf{B}\right)\\ & \text{s.t. }g=\nabla G\in {\mathcal{C}}^1,\end{array}$$

$$F_2({\mathbf{y}}_i,{\mathbf{x}}_i,{\mathbf{z}}_i,g,\mathbf{A},\mathbf{B})\stackrel{\Delta }{=}\mathbb{E}$$

$$-{\mathbf{y}}^{\top }(\mathbf{A}{\mathbf{x}}_i+\mathbf{B}{\mathbf{z}}_i)].$$

$$(\hat{g},\hat{\mathbf{A}})=$$

$$\text{s.t. }g=\nabla G\in {\mathcal{C}}^1.$$

$$\max \left(\mathbb{E}[\Vert \mathbf{x}{\Vert }_2\Vert \mathbf{x}{\Vert }_{\infty }^2],\mathbb{E}[\Vert \mathbf{x}{\Vert }_2\Vert \mathbf{x}{\Vert }_{\infty }\Vert \mathbf{z}{\Vert }_2],\mathbb{E}[\Vert \mathbf{x}{\Vert }_2\Vert \mathbf{z}{\Vert }_2^2]\right)\le s_{Nr},$$

$$\hat{\mathbf{A}}=q(\overline{\mathbf{A}}+\mathbf{L})+\mathbf{E},$$

$$\displaystyle\mathbf{L}=\Big{(}\overline{\mathbf{B}}\mathbb{E}[\mathbf{z}\mathbf{x}^{\top}]+\left(\overline{\mathbf{g}}(\mathbf{0})-\widehat{\mathbf{g}}(\mathbf{0})\right)\boldsymbol{\mu}_{\mathbf{x}}^{\top}\Big{)}\left(\mathbb{E}[\mathbf{x}\mathbf{x}^{\top}]\right)^{\dagger}$$

$$\Rightarrow$$

$$\begin{array}{rl}\frac{1}{MN}\Vert \mathbf{E}{\Vert }_F& \le \frac{2J{\sigma }_{\ell }\sqrt{N}}{{\hat{g}}^{\prime }(0)M}{s}_{Nr}{k}^2,\end{array}$$

$$\hat{\mathbf{y}}=\hat{\mathbf{g}}\left(\left(\hat{\mathbf{A}}+\hat{\mathbf{L}}\right)\mathbf{x}\right),$$

$$(\hat{g},\hat{\mathbf{A}},\hat{\mathbf{L}})=$$

$$\text{s.t. }g=\nabla G\in {\mathcal{C}}^1,$$

$${\hat{F}}_3(\mathbf{Y},\mathbf{X},g,\mathbf{A})=\frac{1}{n}\sum _{i=1}^n\left[\sum _{j=1}^m\left[G_{\ast }\left(y_{ji}\right)+G\left({\mathbf{a}}_j^{\top }{\mathbf{x}}_i\right)\right]-{\mathbf{y}}_i^{\top }\left(\mathbf{A}{\mathbf{x}}_i\right)\right],$$

$$\hat{\mathbf{Y}}=\hat{\mathbf{g}}(\hat{\mathbf{A}}+\hat{\mathbf{L}}).$$

Full text

I Introduction

How real is this relationship? This is a ubiquitous question that presents itself not only in judging interpersonal connections but also in evaluating correlations and causality throughout science and engineering. Two reasons for reaching incorrect conclusions based on observed relationships in collected data are chance and outside influences. For example, we can flip two coins that both show heads, or observe that today’s temperature measurements on the west coast of the continental USA seem to correlate with tomorrow’s on the east coast throughout the year. In the first case, we might not immediately conclude that coins are related, since the number of flips we observe is not very large relative to the possible variance of the process, and the apparent link we observed is up to chance. In the second case, we still may hesitate to use west coast weather to understand and predict east coast weather, since in reality both are closely following a seasonal trend.

Establishing interpretable relationships between entities while mitigating the effects of chance can be achieved via sparse optimization methods, such as regression (Lasso) [1] and inverse covariance estimation [2]. In addition, the extension to time series via vector autoregression [3, 4] yields interpretations related to Granger causality [5]. In each of these settings, estimated nonzero values correspond to actual relations, while zeros correspond to absence of relations.

However, we are often unable to collect data to observe all relevant variables, and this leads to observing relationships that may be caused by common links with those unobserved variables. The hidden variables in this model are fairly general; they can possibly model underlying trends in the data, or the effects of a larger network on an observed subnetwork. For example, one year of daily temperature measurements across a country could be related through a graph based on geographical and meteorological features, but all exhibit the same significant trend due to the changing seasons. We have no single sensor that directly measures this trend. In the literature, a standard pipeline is to de-trend the data as a preprocessing step, and then estimate or use a graph to describe the variations of the data on top of the underlying trends [6, 7, 8].

Alternatively, attempts have been made to capture the effects of hidden variables via sparse plus low-rank optimization [9]. This has been extended to time series [10], and even to a non-linear setting via Generalized Linear Models (GLMs) [11]. What if the form of the non-linearity (link function) is not known? Regression using a GLM with an unknown link function is also known as a Single Index Model (SIM). Recent results have shown good performance when using SIMs for sparse regression [12].

So far, when choosing a model, current methods will impose a fixed form for the (non-)linearity, assume the absence of any underlying trends, perform separate pre-processing or partitioning in an attempt to remove or otherwise explicitly handle such trends, or take some combination of these steps. To address all of these issues, we propose a model with a non-linear function applied to a linear argument that captures the effects of latent variables, which manifest as unmodeled trends in the data. Thus, we introduce the Single Index Latent Variable (SILVar) model, which uses the SIM in a sparse plus low-rank optimization setting to enable general, interpretable multi-task regression in the presence of unknown non-linearities and unobserved variables. That is, we propose the SILVar model not only to use for regression in difficult settings, but also as a tool for uncovering hidden relationships buried in data.

First, we establish notation and review prerequisites in Section II. Next, we introduce the SILVar model and discuss several paradigms in which it can be applied in Section III. Then, we detail the numerical procedure for learning the SILVar model in Section IV. Finally, we demonstrate the performance via experiments on synthetic and real data in Section VI.

II Background

In this section, we introduce the background concepts and notation used throughout the remainder of the paper.

II-A Bregman Divergence and Convex Conjugates

For a given convex function $\phi$, the Bregman Divergence [13] associated with $\phi$ between $\mathbf{y}$ and $\mathbf{x}$ is denoted

[TABLE]

The Bregman Divergence is a non-negative (asymmetric) quasi-metric. Two familiar special cases of the Bregman Divergence are Euclidean Distance when $\phi(\mathbf{x})=\frac{1}{2}\|\mathbf{x}\|_{2}^{2}$ and Kullback-Liebler Divergence when $\phi(\mathbf{x})=\sum\limits_{i}x_{i}\log x_{i}$ in the case that $\mathbf{x}$ is a valid probability distribution (i.e., $\mathbf{x}\geq 0$ and $\sum\limits_{i}x_{i}=1$).

The convex conjugate $\phi_{*}$ of a function $\phi$ is given by

[TABLE]

The convex conjugate arises naturally in optimization problems when deriving a dual form for the original (primal) problem. For closed, convex, differentiable, 1-D function $\phi$ with invertible gradient, the following properties hold

[TABLE]

where $(\cdot)^{-1}$ denotes the inverse function, not the multiplicative inverse. In words, these properties give an alternate form of the conjugate in terms of gradients of the original function, state that the function inverse of the gradient is equal to the gradient of the conjugate, and state that the conjugate of the conjugate is the original function.

II-B Generalized Linear and Single-Index Models

The Generalized Linear Model (GLM) can be described using several parameterizations. We adopt the one based on the Bregman Divergence [14]. For observations $y_{i}\in\mathbb{R}$ and $\mathbf{x}_{i}\in\mathbb{R}^{p}$, let $\mathbf{y}=(y_{1}\ldots y_{n})^{\top}$, $\mathbf{X}=(\mathbf{x}_{1}\ldots\mathbf{x}_{n})$. The model is parameterized by 1) a non-linear link function $g=(\nabla\phi)^{-1}$ where $\phi$ is a closed, convex, differentiable, invertible function; and 2) a vector $\mathbf{a}\in\mathbb{R}^{p}$. We have the model written as

[TABLE]

(note that some references use $g^{-1}$ as the link function where we use $g$) and the corresponding likelihood function written as

[TABLE]

where the likelihood is expressed with respect to an appropriate base measure [15], which can be omitted for notational clarity.

Let $G=\phi_{*}$ and $g=\nabla G=(\nabla\phi)^{-1}$. Then, for data $\{\mathbf{x}_{i},y_{i}\}$ with conditionally independent $y_{i}$ given $\mathbf{x}_{i}$ (note that this is not necessarily assuming that $\mathbf{x}_{i}$ are independent), learning the model $\mathbf{a}$ assuming $g$ is known can be achieved via empirical risk minimization,

[TABLE]

where equality $(a)$ arises from the second property in (3), equality $(b)$ arises from the first property, and we introduce

[TABLE]

for notational compactness.

The Single Index Model (SIM) [16] takes the same form as the GLM. The crucial difference is in the estimation of the models. When learning a GLM, the link function $g$ is known and the linear parameter $\mathbf{a}$ is estimated; however when learning a SIM, the link function $g$ needs to be estimated along with the linear parameter $\mathbf{a}$.

Recently, it has been shown that, when the function $g$ is restricted to be monotonic increasing and Lipschitz, learning SIMs becomes computationally tractable [15] with performance guarantees in high-dimensional settings [12]. Thus, with scalar $u$ defining the set $\mathcal{C}^{u}=\{g:\forall y>x,\;0\leq g(y)-g(x)\leq u(y-x)\}$ of monotonic increasing $u$-Lipschitz functions, this leads to the optimization problem,

[TABLE]

II-C Lipschitz Monotonic Regression and Estimating SIMs

The estimation of $g$ with the objective function including terms $G$ and $G_{*}$ at first appears to be an intractably difficult calculus of variations problem. However, there is a marginalization technique that cleverly avoids directly estimating functional gradients with respect to $G$ and $G_{*}$ [15] and admits gradient-based optimization algorithms for learning. The marginalization utilizes Lipschitz monotonic regression (LMR) as a subproblem. Thus, before introducing this marginalization, we first review LMR.

II-C1 LMR

Given ordered pairs $\{x_{i},y_{i}\}$ and independent Gaussian $w_{i}$, consider the model

[TABLE]

which intuitively treats $\{y_{i}\}$ as noisy observations of a function $g$ indexed by $x$, sampled at points $\{x_{i}\}$. Let $\widehat{g}_{i}=g(x_{i})$, an estimate of the function value, with $\widehat{\mathbf{g}}=(\widehat{g}_{1}\ldots\widehat{g}_{n})^{\top}$, and $x_{[j]}$ denote the $j^{th}$ element of the $\{x_{i}\}$ sorted in ascending order. Then LMR is described by the problem,

[TABLE]

While there may be in fact infinitely many possible (continuous) monotonic increasing Lipschitz functions $g$ that pass through the points $\widehat{\mathbf{g}}$, the solution vector $\widehat{\mathbf{g}}$ is unique. We will introduce a simple yet effective algorithm for solving this problem later in Section IV-A.

II-C2 Estimating SIMs

We now return to the objective function (8). Let $\widehat{\mathbf{g}}=\textrm{LMR}(\mathbf{y},\mathbf{X}^{\top}\mathbf{a})$. Then the gradient w.r.t. $\mathbf{a}$ can be expressed in terms of an estimate of $g$ without explicit computation of $G$ or $G_{*}$,

[TABLE]

This allows us to apply gradient or quasi-Newton methods to solve the minimization in $\mathbf{a}$, which is itself a convex problem since the original problem was jointly convex in $g$ and $\mathbf{a}$.

III Single Index Latent Variable Models

In this section, we build the Single Index Latent Variable (SILVar) model from fundamental concepts.

III-A Multitask regression and latent variables

First, we extend the SIM to the multivariate case and then examine how latent variables can affect learning of the linear parameter. Let $\mathbf{y}_{i}=\left(y_{1i}\;\ldots\;y_{mi}\right)^{\top}$, $\mathbf{g}(\mathbf{x})=\left(g_{1}(x_{1})\;\ldots\;g_{m}(x_{m})\right)^{\top}$, and $\mathbf{A}=\left(\mathbf{a}_{1}\;\ldots\;\mathbf{a}_{m}\right)^{\top}$. Consider the vectorization,

[TABLE]

For the remainder of this paper, we make an assumption that all $g_{j}=g$ for notational simplicity, though the same analysis readily extends to the case where $g_{j}$ are distinct.

Now, let us introduce a set of latent variables $\mathbf{z}_{i}\in\mathbb{R}^{r}$ with $r\ll p$ and the corresponding linear parameter $\mathbf{B}=(\mathbf{b}_{1}\ldots\mathbf{b}_{m})^{\top}\in\mathbb{R}^{m\times r}$ (note we can incorporate a linear offset by augmenting the variable $\mathbf{z}\leftarrow(\mathbf{z}^{\top}\;1)^{\top}$ and adding the linear offset as a column of $\mathbf{B}$). This leads to the asymptotic maximum likelihood estimate,

[TABLE]

where

[TABLE]

Now consider the case in which the true distribution remains the same, but we only observe $\mathbf{x}_{i}$ and not $\mathbf{z}_{i}$,

[TABLE]

We now propose a relation between the two models in (13) and (III-A), which will finally lead to the SILVar model. Here we present the abridged theorem, and relegate the full expressions and derivation to Appendix A. To establish notation, let primes ( ′) denote derivatives, hats ( $\widehat{}$ ) denote a parameter estimate with only observed variables, overbars ($\overline{\vphantom{\mathbf{a}}\hphantom{\mathbf{a}}}$) denote an underlying parameter estimate when we have access to both observed and latent variables, and we drop the subscripts from the random variables $\mathbf{x}$ and $\mathbf{z}$ to reduce clutter.

Theorem 1**.**

Assume that $\widehat{g}^{\prime}(0)\neq 0$ and that $|\widehat{g}^{\prime\prime}|\leq J$ and $|\overline{g}^{\prime\prime}|\leq J$ for some $J<\infty$. Furthermore, assume in models (III-A) and (III-A) that $\max_{j}\left(\|\widehat{\mathbf{a}}_{j}\|_{1},\|\overline{\mathbf{a}}_{j}\|_{1}+\|\overline{\mathbf{b}}_{j}\|_{2}\right)\leq k,$

[TABLE]

where subscripts in $s_{Nr}$ indicate that the bounds may grow with the values in the subscript. Then, the parameters $\widehat{\mathbf{A}}$ and $\overline{\mathbf{A}}$ from models (III-A) and (III-A) are related as

[TABLE]

where $q=\frac{\overline{g}^{\prime}(0)}{\widehat{g}^{\prime}(0)}$, $\boldsymbol{\mu}_{\mathbf{x}}=\mathbb{E}[\mathbf{x}_{i}]$,

[TABLE]

and $\mathbf{E}=\widehat{\mathbf{A}}-q(\overline{\mathbf{A}}+\mathbf{L})$ is bounded as

[TABLE]

where $\sigma_{\ell}=\left\|\left(\mathbb{E}\left[\mathbf{x}\mathbf{x}^{\top}\right]\right)^{\dagger}\right\|_{2}$, the largest singular value of the pseudo-inverse of the covariance.

The proof is given in Appendix A. The assumptions require that the 2nd order Taylor series expansion is accurate around the point $\mathbf{0}$, that the model parameters $\overline{\mathbf{A}}$ and $\overline{\mathbf{B}}$ are bounded, and that the distributions generating the data are not too spread out (beyond $\mathbf{0}$ where the Taylor series expansion is performed). These are all intuitively reasonable and unsurprising assumptions. Though the theorem poses a hard constraint on $|\widehat{g}^{\prime\prime}|$, we hypothesize that this is a rather strong condition that can be weakened to be in line with similar models.

The theorem determines certain scaling regimes under which the sparse plus low-rank approximation remains reasonable. For instance, consider the case where $M\sim N$ and the moments scale as $s_{Nr}\sim\sqrt{N}$, which is reasonable given their form (i.e., very loosely, $\|\mathbf{x}\|_{\infty}\sim 1$ and $\|\mathbf{x}\|_{2}\sim\sqrt{N}$ and small latent power $\|\mathbf{z}\|_{2}\sim 1$ relative to $N$). Then, to keep the error of constant order, the power of each row of the matrix would need to stay constant $k\sim 1$. If we see $\mathbf{A}$ as a graph adjacency matrix, then, in rough terms, this can correspond intuitively to the case in which the node in-degrees stay constant so that the local complexity remains the same even as the global complexity increases while the network grows. Again, we hypothesize that this overall bound can be tightened given more assumptions (e.g., on the network topology). Thus, we propose the SILVar model,

[TABLE]

and learn it using the optimization problem,

[TABLE]

where

[TABLE]

the empirical version of $F_{3}$, and $h_{1}$ and $h_{2}$ are regularizers on $\mathbf{A}$ and $\mathbf{L}$ respectively. Two natural choices for $h_{2}$ would be $h_{2}(\mathbf{L})=\lambda_{2}\|\mathbf{L}\|_{*}$ the nuclear norm and $h_{2}(\mathbf{L})=\mathbb{I}\{\|\mathbf{L}\|_{*}\leq\lambda_{2}\}$ the indicator of the nuclear norm ball, both relating to the nuclear norm of $\mathbf{L}$, since $\mathbf{L}$ is approximately low rank due to the influence of a relatively small number of latent variables. We may choose different forms for $h_{1}$ depending on our assumptions about the structure of $\mathbf{A}$. For example, if $\mathbf{A}$ is assumed sparse, we may use $h_{1}(\mathbf{A})=\lambda_{1}\|\textrm{v}(\mathbf{A})\|_{1}$, the $\ell_{1}$ norm applied element-wise to the vectorized $\mathbf{A}$ matrix. These examples are extensions to the “sparse and low-rank” models, which have been shown under certain geometric incoherence conditions to be identifiable [9]. In other words, if the sparse component is not too low-rank, and if the low-rank component is not too sparse, then $\mathbf{A}$ and $\mathbf{L}$ can be recovered uniquely.

III-B Connection to Related Problems

In this section, we show how the SILVar model can be used in various problem settings commonly considered throughout the literature.

III-B1 Generalized Robust PCA

Though we posed our initial problem as a regression problem, if we take our measurement vectors to be $\mathbf{x}_{i}=\mathbf{e}_{i}$ the canonical basis vectors (i.e., so that the overall data matrix $\mathbf{X}=\mathbf{I}$), then we arrive at

[TABLE]

This is precisely the model behind Generalized Robust PCA [17], but with the twist of estimating the link function as well [18]. What is worth noting is that although we arrived at our model via a regression problem with latent variables, the model itself is also able to describe a problem that arises from very different assumptions on how the data is generated and structured.

We also note that the SILVar model can be modified to share a space with the Generalized Low-Rank (GLR) Models [19]. The GLR framework is able to describe arbitrary types of data with an appropriate choice of convex link function $g$ determined a priori, while the SILVar model is restricted to a certain continuous class of convex link functions but aims to learn this function. The modification is simply a matrix factorization $\mathbf{L}=\mathbf{U}\mathbf{V}$ (and “infinite” penalization on $\mathbf{A}$). The explicit factorization makes the problem non-convex but instead block convex, which still allows for alternating convex steps in $\mathbf{U}$ with fixed $\mathbf{V}$ (and vice versa) to tractably reach local minima under certain conditions. Nonetheless, due to the non-convexity, further discussion of this particular extension will be beyond the scope of this paper.

III-B2 Extension to Autoregressive Models

We can apply the SILVar model to learn from time series as well. Consider a set of $N$ time series each of length $K$, $\mathbf{X}\in\mathbb{R}^{N\times K}$. We assume the noise at each time step is independent (note that, with this assumption, the time series are still dependent across time), and take in our previous formation, $\mathbf{y}_{i}\leftarrow\mathbf{x}_{k}$ and $\mathbf{x}_{i}\leftarrow\mathbf{x}_{k-1:k-M}=(\mathbf{x}_{k-1}^{\top}\ldots\mathbf{x}_{k-M}^{\top})^{\top}$ so that the model of order $M$ takes the form,

[TABLE]

and learn it using the optimization problem,

[TABLE]

where $\mathbf{A}=\left(\mathbf{A}^{(1)}\ldots\mathbf{A}^{(M)}\right)$ and $\mathbf{L}=\left(\mathbf{L}^{(1)}\ldots\mathbf{L}^{(M)}\right)$ and

[TABLE]

where $\mathbf{A}^{(i)}\!=\!\left(\mathbf{a}_{1}^{(i)}\ldots\mathbf{a}_{N}^{(i)}\right)^{\top}$, similarly to before. Note that the analysis in the previous section follows naturally in this setting, so that here $\textrm{rank}(\mathbf{L}_{i})\leq r+1$. Then, the matrix $\mathbf{A}$ may be assumed to be group sparse, relating to generalized notions of Granger Causality [20, 3], and one possible corresponding choice of regularizer taking the form $h_{1}(\mathbf{A})=\lambda_{1}\sum\limits_{i,j}\left\|\left(a^{(1)}_{ij}\ldots a^{(M)}_{ij}\right)\right\|_{2}$.

Another structural assumption could be that of the Causal Graph Process model [21], inspired by Signal Processing on Graphs [6], in which $\mathbf{A}^{(i)}$ are matrix polynomials in one underlying matrix $\widetilde{\mathbf{A}}$. This framework utilizes the nonconvex regularizer $h_{1}(\mathbf{A})=\lambda_{1}\|\textrm{v}(\mathbf{A}^{(1)})\|_{1}+\lambda_{2}\sum\limits_{i\neq j}\|\mathbf{A}^{(i)}\mathbf{A}^{(j)}-\mathbf{A}^{(j)}\mathbf{A}^{(i)}\|_{F}^{2}$ to encourage both sparsity and commutativity, which is satisfied if $\mathbf{A}^{(i)}$ are all matrix polynomials in the same matrix. Since this particular regularization is again block convex, convex steps can still be taken in each $\mathbf{A}^{(i)}$ with all other blocks $\mathbf{A}^{(j)}$ for $j\neq i$ fixed, for a computationally tractable algorithm to reach a local minimum under certain conditions. However, further detailed discussion will remain outside the scope of this paper.

The hidden variables in this time series model can even possibly model underlying trends in the data. For example, one year of daily temperature measurements across a country could be related through a graph based on geographical and meteorological features, but all exhibit the same significant trend due to the changing seasons. In previous work, a standard pipeline is to detrend the data as a preprocessing step, and then estimate or use a graph to describe the variations of the data on top of the underlying trends [6, 7, 8]. Instead, the time series can also be modeled as a modified autoregressive process, depending on a low-rank trend $\mathbf{L}^{\prime}=\left(\boldsymbol{\ell}^{\prime}_{1}\;\ldots\;\boldsymbol{\ell}^{\prime}_{K}\right)\in\mathbb{R}^{N\times K}$ and the variations of the process about that trend,

[TABLE]

Substituting this into Equation (20) yields

[TABLE]

Thus we estimate the trend using ridge regression for numerical stability purposes but without enforcing $\mathbf{L}^{\prime}$ to be low rank. We can accomplish this via the simple optimization,

[TABLE]

with $\lambda$ being the regularization parameter. In this way, the extension of SILVar to autoregressive models can allow joint estimation of the effects of the trend and of these variations supported on the graph, as will be demonstrated via experiments in Section VI.

IV Efficiently Learning SILVar Models

In this section, we describe the formulation and corresponding algorithm for learning the SILVar model. Surprisingly, in the single-task setting, learning a SIM is jointly convex in $g$ and $\mathbf{a}$ as demonstrated in [15]. The pseudo-likelihood functional $\widehat{F}_{3}$ used for learning the SILVar model in (21) is thus also jointly convex in $g$, $\mathbf{A}$, and $\mathbf{L}$ by a simple extension from the single-task regression setting.

Lemma 2** (Corollary of Theorem 2 of [15]).**

The $\widehat{F}_{3}$ in the SILVar model learning problem (16) is jointly convex in $g$, $\mathbf{A}$, and $\mathbf{L}$.

This convexity is enough to ensure that the learning can converge and be computationally efficient. Before describing the full algorithm, one detail remains: the implementation of the LMR introduced in Section II-C.

IV-A Lipschitz Monotonic Regression

To tackle LMR, we first introduce the related simpler problem of monotonic regression, which is solved by a well-known algorithm, the pooled adjacent violators (PAV) [22]. The monotonic regression problem is formulated as

[TABLE]

The PAV algorithm has a complexity of $O(n)$, which is due to a single complete sweep of the vector $\mathbf{y}$. The monotonic regression problem can also be seen as a standard $\ell_{2}$ projection onto the convex set of monotonic functions.

We introduce a simple generalization to the monotonic regression,

[TABLE]

This modification allows weaker (if $t_{i}<0$) or stronger (if $t_{i}>0$) local monotonicity conditions to be placed on the estimated function $g$. Let $t_{[1]}=0$ and $s_{[i]}=\sum\limits_{j=1}^{i}t_{[j]}$, and in vectorized form,

[TABLE]

Then, we can rewrite,

[TABLE]

which leads us to recognize that

[TABLE]

Returning to LMR (10), we pose it as the projection onto the intersection of 2 convex sets, the monotonic functions and the functions with upper bounded differences, for which each projection can be solved efficiently using GPAV (for the second set, we can use GPAV on the negated function and negate the result). Thus to perform this optimization efficiently utilizing PAV as a subroutine, we can use an accelerated Dykstra’s Algorithm [23], which is a general method to find the projection onto the non-empty intersection of any number of convex sets. The accelerated algorithm can has a geometric convergence rate for $O(\log n)$ passes, with each pass utilizing PAV with $O(n)$ cost, for a total cost of $O(n\log n)$. We make a brief sidenote that this is better than the $O(n^{2})$ in [24] and simpler than the $O(n\log n)$ in [25] achieved using a complicated tree-like data structure, and has no learning rate parameter to tune compared to an ADMM-based implementation that would also yield $O(n\log n)$. We provide the steps of LMR in Algorithm 1 as a straightforward application of the the accelerated Dykstra’s algorithm, but leave the details of the derivation of the acceleration to [23]. In addition, the numerical stability parameter is included to handle small denominators, but in practice we did not observe the demoninator falling below $\varepsilon=10^{-9}$ in our simulations before convergence.

While the GPAV allows us to quickly perform Lipschitz monotonic regression with our particular choice of $\mathbf{t}$, it should be noted that using other choices for $\mathbf{t}$ can easily generalize the problem to bi-Hölder regression as well, to find functions in the set $\mathcal{C}_{\ell,\alpha}^{u,\beta}=\{g:\forall y>x,\;\ell|y-x|^{\alpha}\leq g(y)-g(x)\leq u|y-x|^{\beta}\}$. This set may prove interesting for further analysis of the estimator and is left as a topic for future investigation.

IV-B Optimization Algorithms

With LMR in hand, we can outline the algorithm for solving the convex problem (III-A). This procedure can be performed for general $h_{1}$ and $h_{2}$, but proximal mappings are efficient to compute for many common regularizers. Thus, we describe the basic algorithm using gradient-based proximal methods (e.g., accelerated gradient [26], and quasi-Newton [27]), which require the ability to compute a gradient and the proximal mapping.

In addition, we can compute the objective function value for purposes of backtracking [28], evaluating algorithm convergence, or computing error (e.g., for validation or testing purposes). While our optimization outputs an estimate of $g$, the objective function depends on the value of $G$ and its conjugate $G_{*}$. Though we cannot necessarily determine $G$ or $G_{*}$ uniquely since $G+c$ and $G$ both yield $g$ as a gradient for any $c\in\mathbb{R}$, the value of $G(x)+G_{*}(y)$ is unique for a fixed $g$. To see this, consider $\widetilde{G}=G+c$ for some constant $c$. Then,

[TABLE]

This allows us to compute the objective function by performing the cumulative numerical integral of $\widehat{\mathbf{g}}$ on points $\textrm{v}(\boldsymbol{\Theta})$ (e.g., G=cumtrapz(theta,ghat) in Matlab). Then, a discrete convex conjugation (also known as the discrete Legendre transform (DLT) [29]) computes $G_{*}$.

We did not notice significant differences in performance accuracy due to our implementation using quasi-Newton methods and backtracking compared to an accelerated proximal gradient method, possibly because both algorithms were run until convergence. Thus, we show only one set of results. However, we note that the runtime was improved by implementing the quasi-Newton and backtracking method.

Algorithm 2 describes the learning procedure and details the main function and gradient computations while assuming a proximal operator is given. The computation of the gradient and the update vector depends on the particular variation of proximal method utilized; with stochastic gradients, the set $\mathcal{I}\subset\{1,\ldots,n\}$ could be pseudorandomly generated at each iteration, while with standard gradients, $\mathcal{I}=\{1,\ldots,n\}$. The cumtrapz procedure takes as input the coordinates of points describing the function $\widehat{\mathbf{g}}$. The DLT procedure takes as its first two inputs the coordinates of points describing the function $\widehat{\mathbf{G}}$, and as its third input the points at which to evaluate the convex conjugate $\widehat{\mathbf{G}}_{*}$.

Here, an observant reader may notice the subtle point that the function $G_{*}$ may only be defined inside a finite or semi-infinite interval. This can occur if the function $g$ is bounded below and/or above, so that $G$ has a minimum/maximum slope, and $G_{*}$ is infinite for any values below that minimum or above that maximum slope of $G$ (to see this, one may refer back to the definition of conjugation (2)). Fortunately, this does not invalidate our method. It is straightforward to see that $\widehat{g}(x)=x$ is always a feasible solution and that $\widehat{G}_{*}(y)=y^{2}/2$ is defined for all $y$; thus starting from this solution, with appropriately chosen step sizes, gradient-based algorithms will avoid the solutions that make the objective function infinite. Furthermore, even if new data $\mathbf{y}_{j}$ falls outside the valid domain for the learned $\widehat{G}_{*}$ and we incur an “infinite” loss using the model, evaluating $\widehat{\mathbf{g}}\left(\left(\widehat{\mathbf{A}}+\widehat{\mathbf{B}}\right)\mathbf{x}_{j}\right)$ is still well defined. This problem is not unique to SIM’s, as assuming a fixed link function $g$ in a GLM can also incur infinite loss if new data does not conform to modeling assumptions implicit to the link function (e.g., the log-likelihood for a negative $\mathbf{y}$ under a non-negative distribution), and making a prediction using the learned GLM is still well-defined. Practically, the loss for SILVar computed using the DLT may be large, but will not be infinite [29].

V Performance

We now provide conditions under which the optimization recovers the parameters of the model.

First let us establish a few additional notations needed. Let $(\widetilde{g},\widetilde{\mathbf{A}},\widetilde{\mathbf{L}})$ be the best Lipschitz monotonic function, sparse matrix, and low-rank matrix that models the true data generation. Let $\widetilde{\mathbf{A}}$ be sparse on the index set $S$ of size $|S|=s_{\mathbf{A}}$ and $\widetilde{\mathbf{L}}=\mathbf{U}\Lambda\mathbf{V}^{\top}$ be the SVD of $\widetilde{L}$, where $\mathbf{U}\in\mathbb{R}^{M\times r_{\mathbf{L}}}$, $\Lambda\in\mathbb{R}^{r_{\mathbf{L}}\times r_{\mathbf{L}}}$, and $\mathbf{V}\in\mathbb{R}^{N\times r_{\mathbf{L}}}$ where $r_{\mathbf{L}}$ is the rank of $\widetilde{\mathbf{L}}$. Then let $\mathbf{A}_{S}$ be equal to $\mathbf{A}$ on $S$ and [math] on $S^{c}$ (so that $\widetilde{\mathbf{A}}_{S}=\widetilde{\mathbf{A}}$ and $\mathbf{A}_{S}+\mathbf{A}_{S^{c}}=\mathbf{A}$), and $\mathbf{L}_{R}=\mathbf{L}-(\mathbf{I}-\mathbf{U}\mathbf{U}^{\top})\mathbf{L}(\mathbf{I}-\mathbf{V}\mathbf{V}^{\top})$ and $\mathbf{L}_{R^{c}}=\mathbf{L}-\mathbf{L}_{R}$ (so that $\widetilde{\mathbf{L}}_{R}=\widetilde{\mathbf{L}}$ and $\mathbf{L}_{R^{c}}+\mathbf{L}_{R}=\mathbf{L}$). Consider the set of approximately $S$-sparse and $R$-low-rank matrices $\mathcal{B}_{\gamma}(S,R)=\{(\boldsymbol{\Phi},\boldsymbol{\Psi}):\gamma\|\boldsymbol{\Phi}_{S^{c}}\|_{1}+\|\boldsymbol{\Psi}_{R^{c}}\|_{*}\leq 3(\gamma\|\boldsymbol{\Phi}_{S}\|_{1}+\|\boldsymbol{\Psi}_{R}\|_{*})\}$. Intuitively, this is the set of matrices for which the energy of $\boldsymbol{\Phi}$ is on the same sparsity set as $\widetilde{\mathbf{A}}$ and similarly the energy of $\boldsymbol{\Psi}$ is in the same low-rank space as $\widetilde{\mathbf{L}}$. Finally, let the marginalization of the loss functional w.r.t. $g$ be $\widehat{m}(\mathbf{Y},\mathbf{X},\mathbf{A})=\underset{g\in\mathcal{C}^{1}}{\min}\widehat{F}_{3}(\mathbf{Y},\mathbf{X},g,\mathbf{A})$, the value of the marginalized function at the true matrix parameters be $\widehat{\widetilde{g}}=\underset{g\in\mathcal{C}^{1}}{\min}\widehat{F}_{3}(\mathbf{Y},\mathbf{X},g,\widetilde{\mathbf{A}}+\widetilde{\mathbf{L}})$, the Hessian of the marginalized functional w.r.t. the matrix parameter be $\widetilde{\mathcal{I}}=\nabla^{2}_{\textrm{v}(\mathbf{A})}\widehat{m}(\mathbf{Y},\mathbf{X},\widetilde{\mathbf{A}}+\widetilde{\mathbf{L}})$, and denote the quadratic form in shorthand $\|\mathbf{v}\|_{\widetilde{\mathcal{I}}}=\mathbf{v}^{\top}\widetilde{\mathcal{I}}\mathbf{v}$.

Now, consider the following assumptions:

There exists some $\alpha>0$ such that

[TABLE]

This is essentially a Restricted Strong Convexity condition standard for structured recovery [30]. 2. 2.

Let $\widehat{\widetilde{\boldsymbol{\Gamma}}}\in\mathbb{R}^{M\times K}$ be the matrix with columns given by

[TABLE]

Then for some $\lambda>0$ and $\gamma>0$,

[TABLE]

where $\|\mathbf{A}\|_{\infty}$ is the largest element of $\mathbf{A}$ in magnitude, and $\|\mathbf{A}\|_{2}$ is the largest singular value of $\mathbf{A}$.

This states that the error is not too powerful and is not too correlated with the regression variables, and it is also similar to standard conditions for structured recovery. 3. 3.

Let $\tau=\gamma\sqrt{\frac{s_{\mathbf{A}}}{r_{\mathbf{L}}}}$ and $\mu=\frac{1}{32r_{\mathbf{L}}(1+\tau^{2})}$, then

[TABLE]

This is a condition on the incoherence between the sparsity set $S$ and low rank subspace $R$ (see [31, 9, 32, 17] for other similar variations on rank-sparsity incoherence considered previously). Basically, this states that $S$ should be such that $S$-sparse matrices are not approximately low-rank, while $R$ should be such that $R$-low-rank matrices are not approximately sparse.

Then we have the following result,

Theorem 3**.**

Under Assumptions 1-3, the solution to the optimization problem (21) satisfies

[TABLE]

We give the proof for this theorem in the Appendix. The theorem relates the performance of the optimization to conditions on the problem parameters. Of course, in practice these conditions are difficult to check, since they require knowledge of ground truth. Even given ground truth, verifying Assumption 3 requires solving another separate non-convex optimization problem. Thus, future directions for work could include investigation into which kinds of sparse matrix (or network) models for $\mathbf{A}$ and low-rank models $\mathbf{L}$ produce favorable scaling regimes in the problem parameters with high probability.

VI Experiments

We study the performance of the algorithm via simulations on synthetically generated data as well as real data. In these experiments, we show the different regression settings, as discussed in Section III-B, under which the SILVar model can be applied.

VI-A Synthetic Data

The synthetic data corresponds to a multi-task regression setting. The data was generated by first creating random sparse matrix $\mathbf{A}_{f}=(\mathbf{A}\;\mathbf{A}_{h})$ where $\mathbf{B}\in\mathbb{R}^{p\times H}$ and $H=\lfloor hp\rfloor$ the number of hidden variables, and $h$ is the proportion of hidden variables. The elements of $\mathbf{A}_{f}$ were first generated as i.i.d. Normal variables, and then a sparse mask was applied to $\mathbf{A}$ choose $\lfloor m\log_{10}(p)\rfloor$ non-zeros, and $\mathbf{B}$ was scaled by a factor of $\frac{1}{3\sqrt{H}}$. Next, the data $\mathbf{X}_{f}=(\mathbf{X}^{\top}\;\mathbf{Z}^{\top})^{\top}$ were generated with each vector drawn i.i.d. from a multivariate Gaussian with mean $\mathbf{0}$ and covariance matrix $\boldsymbol{\Sigma}_{f}\in\mathbb{R}^{(p+H)\times(p+H)}$. This was achieved by creating $\boldsymbol{\Sigma}_{f}^{1/2}$ by thresholding a matrix with diagonal entries of $1$ and off-diagonal entries drawn from an i.i.d. uniform distribution in the interval $[-0.5,0.5]$ to be greater than 0.35 in magnitude. Then $\boldsymbol{\Sigma}_{f}^{1/2}$ was pre-multiplied with a matrix of i.i.d. Normal variables. Finally, we generated $\mathbf{Y}=g(\mathbf{A}_{f}\mathbf{X}_{f})+\mathbf{W}$ for 2 different link functions $g_{1}(x)=\log(1+e^{c_{1}x})$ and $g_{2}(x)=\frac{2}{1+e^{-c_{2}x}}-1$, and $\mathbf{W}$ is added i.i.d. Gaussian noise. Then, the task is to learn $(\widehat{g},\widehat{\mathbf{A}},\widehat{\mathbf{L}})$ the SILVar model (16) from $(\mathbf{Y},\mathbf{X})$ without access to $\mathbf{Z}$. As comparison, we use an Oracle GLM model, in which the true $g$ is given but the parameters $(\widehat{\mathbf{A}}_{o},\widehat{\mathbf{L}}_{o})$ still need to be learned. We note that while there is a slight mismatch between the true and assumed noise distributions, the task is nonetheless difficult, yet our estimation using the SILVar model can still exhibit good performance with respect to the Oracle.

The experiments were carried out across a range of problem sizes. The dimension of $\mathbf{y}_{i}$ was fixed at $m=25$. The dimension of the observed $\mathbf{x}_{i}$ was set to $p\in\{25,50\}$, and the proportion of the dimension of hidden $\mathbf{z}_{i}$ was set to $h\in\{0.1,0.2\}$. The number of data samples was varied in $k\in\{25,50,100,150,200\}$. Validation was performed using a set of samples generated independently but using the same $\mathbf{A}_{f}$ and $g$, and using a grid of (\lambda_{S},\lambda_{L})\in\left\{10^{i/4}\big{|}i\in\{-8,-7,\ldots,7,8\}\right\}^{2}. The whole process of data generation and model learning is repeated 20 times for each experimental condition. The average $\ell_{1}$ errors between the true $\mathbf{A}$ and the best estimated $\widehat{\mathbf{A}}$ (determined via validation) are shown in Figure 1. The first row is for the function $g_{1}$ and the second row is for the function $g_{2}$ for the various experimental conditions. The empirical results show that the SILVar model and learning algorithm lose out slightly on performance for estimating the link function $g$ but overall still match the Oracle fairly well in most cases. We also see that SILVar performs better than just the sparse SIM [12], but that the purely sparse Oracle GLM can still provide an improvement over SILVar for not knowing the link function. We also see several other intuitive behaviors for both the SILVar and Oracle models: as we increase the amount of data $n$, the performance increases. Note that due to the scaling of weights $\mathbf{B}$ by the factor of $1/\sqrt{H}$ to keep the average power $\|\mathbf{B}\|_{F}$ roughly constant, increasing the proportion of hidden variables (from $h=0.1$ to $h=0.2$) did not directly lead to significant performance decreases. We also note that $g_{2}$ seems to be somehow more difficult to estimate than $g_{1}$, and additionally that the performance of the SILVar model w.r.t. the Oracle model is worse with $g_{2}$ than with $g_{1}$. This could have something to do with the 2 saturations in $g_{2}$ as opposed to the 1 saturation in $g_{1}$, corresponding in an intuitive sense to a more non-linear setting that contains more information needing to be captured while estimating $\widehat{g}$.

VI-B Temperature Data

In this setting, we wish to learn the graph capturing relations between the weather patterns at different cities. The data is a real world multivariate time series consisting of daily temperature measurements (in ∘F) for $365$ consecutive days during the year $2011$ taken at each of $150$ different cities across the continental USA.

Previously, the analysis on this dataset has been performed by first fitting with a $4$th order polynomial and then estimating a sparse graph from an autoregressive model using a known link function $g(x)=x$ assuming Gaussian noise [8].

Here, we fit the time series using a $2$nd order AR SILVar model (20) with regularizers for group sparsity $h_{1}(\mathbf{A})=\lambda_{1}\sum\limits_{i,j}\left\|\left(a^{(1)}_{ij}\ldots a^{(M)}_{ij}\right)\right\|_{2}$ where $a_{ij}^{(m)}$ is the $ij$ entry of matrix $\mathbf{A}^{(m)}$, and nuclear norm $h_{2}(\mathbf{L})=\lambda_{2}\sum\limits_{i=1}^{M}\left\|\mathbf{L}^{(i)}\right\|_{*}$. We also estimate the underlying trend using (24).

In Figure 2, we plot the original time series, estimated trends, the estimated network, and residuals. The plots are all from time steps $3$ to $363$, since those are the indices for the trends that we can reliably estimate using the simple procedure (24). In Figure 2(a), we show the original time series, which clearly exhibit a seasonal trend winter–summer–winter. Figure 2(b) shows the estimated link function $\widehat{g}$, which turns out to be linear, though a priori we might not intuitively expect a process as complicated as weather to behave this way when sampled so sparsely. Figures 2(c) and 2(d) show the full prediction $\widehat{\mathbf{x}}_{k}=\sum\limits_{i=1}^{M}(\widehat{\mathbf{A}}^{(i)}+\widehat{\mathbf{L}}^{(i)})\mathbf{x}_{k-i}$ and the residuals $\mathbf{x}_{k}-\widehat{\mathbf{x}}_{k}$, respectively. Note the much lower vertical scale in Figure 2(d). Figures 2(e) and 2(f) show the trend $\widehat{\mathbf{L}}^{\prime}$ learned using (24) and the time series with the estimated trend removed, respectively. The trend estimation procedure captures the basic shape of the seasonal effects on the temperature. Several of the faster fluctuations in the beginning of the year are captured as well, suggesting that they were caused by some larger scale phenomena. Indeed there were several notable storm systems that affected the entire USA in the beginning of the year in short succession [33, 34].

Figure 3 compares two networks $\widehat{\mathbf{A}}^{\prime}$ estimated using SILVar and using just sparse SIM without accounting for the low-rank trends, both with the same sparsity level of $12\%$ non-zeros for display purposes, and where $\widehat{a}^{\prime}_{ij}=\left\|\left(\widehat{a}^{(1)}_{ij}\ldots\widehat{a}^{(M)}_{ij}\right)\right\|_{2}$. Figure 3(a) shows the network $\widehat{\mathbf{A}}^{\prime}$ that is estimated using SILVar. The connections imply predictive dependencies between the temperatures in cities connected by the graph. It is intuitively pleasing that the patterns discovered match well previously established results based on first de-trending the data and then separately estimating a network [8]. That is, we see the effect of the Rocky Mountain chain around $-110^{\circ}$ to $-105^{\circ}$ longitude and the overall west-to-east direction of the weather patterns, matching the prevailing winds. In contrast to that of SILVar, the graph estimated by the sparse SIM shown in Figure 3(b) on the other hand has many additional connections with no basis in actual weather patterns. Two particularly unsatisfying cities are: sunny Los Angeles, California at $(-118,34)$, with its multiple connections to snowy northern cities including Fargo, North Dakota at $(-97,47)$; and Caribou, Maine at $(-68,47)$, with its multiple connections going far westward against prevailing winds including to Helena, Montana at $(-112,47)$. These do not show in the graph estimated by SILVar and shown in Figure 3(a).

VI-C Bike Traffic Data

The bike traffic data was obtained from HealthyRide Pittsburgh [35]. The dataset contains the timestamps and station locations of departure and arrival (among other information) for each of 127,559 trips taken between 50 stations within the city from May 31, 2015 to September 30, 2016, a total of 489 days.

We consider the task of using the total number of rides departing from and arriving in each location at 6:00AM-11:00AM to predict the number of rides departing from each location during the peak period of 11:00AM-2:00PM for each day. This corresponds to $\mathbf{Y}\in\mathbb{N}_{0}^{50\times 489}$ and $\mathbf{X}\in\mathbb{N}_{0}^{100\times 489}$, where $\mathbb{N}_{0}$ is the set of non-negative integers, and $\mathbf{A},\mathbf{L}\in\mathbb{R}^{50\times 100}$. We estimate the SILVar model (16) and compare its performance against a sparse plus low-rank GLM model with an underlying Poisson distribution and fixed link function $g_{\scriptsize{\textrm{GLM}}}(x)=\log(1+e^{x})$. We use $n\in\{60,120,240,360\}$ training samples and compute errors on validation and test sets of size $48$ each, and learn the model on a grid of (\lambda_{S},\lambda_{L})\in\left\{10^{i/4}\big{|}i\in\{-8,-7,\ldots,11,12\}\right\}^{2}. We repeat this 10 times for each setting, using an independent set of training samples each time. We compute testing errors in these cases for the optimal $(\lambda_{S},\lambda_{L})$ with lowest validation errors for both SILVar and GLM models.

We also demonstrate that the low-rank component of the estimated SILVar model captures something intrinsic to the data. Naturally, we expect people’s behavior and thus traffic to be different on business days and on non-business days. A standard pre-processing step would be to segment the data along this line and learn two different models. However, as we use the full dataset to learn one single model, we hypothesize that the learned low-rank $\widehat{\mathbf{L}}$ captures some aspects of this underlying behavior.

Figure 4(a) shows the test Root Mean Squared Errors (RMSEs) for both SILVar and GLM models for varying training sample sizes, averaged across the 10 trials. We see that the SILVar model outperforms the GLM model by learning the link function in addition to the sparse and low-rank regression matrices. Figure 4(b) shows an example of the link function learned by the SILVar model with $n=360$ training samples. Note that the learned SILVar link function performs non-negative clipping of the output, which is consistent with the count-valued nature of the data.

We also test the hypothesis that the learned $\widehat{\mathbf{L}}$ contains information about whether a day is business or non-business through its relatively low-dimensional effects on the data. We perform the singular value decomposition (SVD) on the optimally learned $\widehat{\mathbf{L}}=\widehat{\mathbf{U}}\widehat{\boldsymbol{\Sigma}}\widehat{\mathbf{V}}^{\top}$ for $n=360$ and project the data onto the $r$ top singular components $\widetilde{\mathbf{X}}_{r}=\widehat{\boldsymbol{\Sigma}}_{r}\widehat{\mathbf{V}}_{r}^{\top}\mathbf{X}$. We then use $\widetilde{\mathbf{X}}_{r}$ to train a linear support vector machine (SVM) to classify each day as either a business day or a non-business day, and compare the performance of this lower dimensional feature to that of using the full vector $\mathbf{X}$ to train a linear SVM. If our hypothesis is true then the performance of the classifier trained on $\widetilde{\mathbf{X}}_{r}$ should be competitive with that of the classifier trained on $\mathbf{X}$. We use 50 training samples of $\widetilde{\mathbf{X}}_{r}$ and of $\mathbf{X}$ and test on the remainder of the data. We repeat this 50 times by drawing a new batch of $50$ samples each time. We then vary the proportion of business to non-business days in the training sample to trace out a receiver operating characteristic (ROC).

In Figure 5, we see the results of training linear SVM on $\widetilde{\mathbf{X}}_{r}$ for $r\in\{1,...,6\}$ and on the full data for classifying business and non-business days. We see that using only the first two singular vectors, the performance is fairly poor. However, by simply taking 3 or 4 singular vectors, the classification performance almost matches that of the full data. Surprisingly, using the top 5 or 6 singular vectors achieves performance greater than that of the full data. This suggests that the projection may even play the role of a de-noising filter in some sense. Since we understand that behavior should correlate well with the day being business or non-business, this competitive performance of the classification using the lower dimensional features strongly suggests that the low-rank $\widehat{\mathbf{L}}$ indeed captures the effects of latent behavioral factors on the data.

Finally, in Figure 6, we plot the $(i,i)$ entries of the optimal $\widehat{\mathbf{A}}$ at $n=360$. This corresponds to locations for which incoming bike rides at 6:00AM-11:00AM are good predictors of outgoing bike rides at 11:00AM-2:00PM, beyond the effect of latent factors such as day of the week. We may intuitively expect this to correlate with locations that have restaurants open for lunch service, so that people would be likely to ride in for lunch or ride out after lunch. This is confirmed by observing that these stations are in Downtown (-80,40.44), the Strip District (-79.975, 40.45), Lawrenceville (-79.96, 40.47), and Oakland (-79.96, 40.44), known locations of many restaurants in Pittsburgh. It is especially interesting to note that Oakland, sandwiched between the University of Pittsburgh and Carnegie Mellon University, is included. Even though the target demographic is largely within walking distance, there is a high density of restaurants open for lunch, which may explain its non-zero coefficient. The remainder of the locations with non-zero coefficients $a_{ii}$ are also near high densities of lunch spots, while the other locations with coefficients $a_{ii}$ of zero are largely either near more residential areas or near restaurants more known for dinner or nightlife rather than lunch, such as Shadyside ($x\geq-79.95$) and Southside ($y\leq 40.43$)).

VII Conclusion

Data exhibit complex dependencies, and it is often a challenge to deal with non-linearities and unmodeled effects when attempting to uncover meaningful relationships among various interacting entities that generate the data. We introduce the SILVar model for performing semi-parametric sparse regression and estimating sparse graphs from data under the presence of non-linearities and latent factors or trends. The SILVar model estimates a non-linear link function $g$ as well as structured regression matrices $\mathbf{A}$ and $\mathbf{L}$ in a sparse and low-rank fashion. We justify the form of the model and relate it to existing methods for general PCA, multi-task regression, and vector autoregression. We provide computationally tractable algorithms for learning the SILVar model and demonstrate its performance against existing regression models and methods on both simulated and real data sets, namely 2011 US weather sensor network data and 2015-2016 Pittsburgh bike traffic data. We see from the simulated data that the SILVar model matches the performance of an Oracle GLM that knows the true link function and only needs to estimate $\mathbf{A}$ and $\mathbf{L}$; we show empirically on the temperature data that the learned $\mathbf{L}$ can capture the effects of underlying trends in time series while $\mathbf{A}$ represents a graph consistent with US weather patterns; and we see that, in the bike data, SILVar outperforms a GLM with a fixed link function, the learned $\mathbf{L}$ encodes latent behavioral aspects of the data, and $\mathbf{A}$ discovers notable locations consistent with the restaurant landscape of Pittsburgh.

Appendix A Proof of Theorem 1

We first restate the theorem for convenience.

Theorem 1**.**

Assume that $\widehat{g}^{\prime}(0)\neq 0$ and that $|\widehat{g}^{\prime\prime}|\leq J$ and $|\overline{g}^{\prime\prime}|\leq J$ for some $J<\infty$. Furthermore, assume in models (III-A) and (III-A) that $\max_{j}\left(\|\widehat{\mathbf{a}}_{j}\|_{1},\|\overline{\mathbf{a}}_{j}\|_{1}+\|\overline{\mathbf{b}}_{j}\|_{2}\right)\leq k,$

[TABLE]

where subscripts in the parameters on the RHS indicate that the bounds may grow with the values in the subscript. Then, the parameters $\widehat{\mathbf{A}}$ and $\overline{\mathbf{A}}$ from models (III-A) and (III-A) are related as

[TABLE]

where $q=\frac{\overline{g}^{\prime}(0)}{\widehat{g}^{\prime}(0)}$, $\boldsymbol{\mu}_{\mathbf{x}}=\mathbb{E}[\mathbf{x}_{i}]$,

[TABLE]

and $\mathbf{E}=\widehat{\mathbf{A}}-q(\overline{\mathbf{A}}+\mathbf{L})$ is bounded as

[TABLE]

where $\sigma_{\ell}=\left\|\left(\mathbb{E}\left[\mathbf{x}\mathbf{x}^{\top}\right]\right)^{\dagger}\right\|_{2}$, the largest singular value of the pseudo-inverse of the covariance.

Proof. We begin with the Taylor series expansion, using the Lagrange form of the remainder,

[TABLE]

where

[TABLE]

$\overline{\mathbf{B}}\!=\!\left(\overline{\mathbf{b}}_{1}\ldots\overline{\mathbf{b}}_{m}\right)^{\top}$ similarly to before, $\mathbf{x}\odot\mathbf{y}$ denotes the element-wise (Hadamard) product, and $\mathbf{x}.^{\wedge}2=\mathbf{x}\odot\mathbf{x}$ denotes element-wise squaring. First let us consider the quantity

[TABLE]

Then,

[TABLE]

where $(a)$ follows from the fact that $\|\mathbf{x}\|_{4}\leq\|\mathbf{x}\|_{2}\;\forall\mathbf{x}\in\mathbb{R}^{N}$.

Similarly, we have

[TABLE]

To finish off, substituting the definition of $\mathbf{E}$ and $q$ into (33) we have

[TABLE]

∎

Appendix B Proof of Theorem 3

First, we present and prove some intermediate results. For convenience of notation, for the remainder of this section let $\boldsymbol{\Phi}=\widehat{\mathbf{A}}-\widetilde{\mathbf{A}}$ and $\boldsymbol{\Psi}=\widehat{\mathbf{L}}-\widetilde{\mathbf{L}}$.

B-A Propositions

Proposition 4.1**.**

Under Assumption 2, the solution $(\widehat{\mathbf{A}},\widehat{\mathbf{L}})$ to the optimization problem (21) satisfies

[TABLE]

Proof. We start with the convexity of the marginalized objective functional at $(\widetilde{\mathbf{A}},\widetilde{\mathbf{L}})$,

[TABLE]

Then consider the optimality of the full objective functional at $(\widehat{g},\widehat{\mathbf{A}},\widehat{\mathbf{L}})$,

[TABLE]

where the last inequality utilizes convexity of the marginalized objective. Then using Assumption 2,

[TABLE]

where the penultimate inequality arises from decomposability of the norm. Specifically,

[TABLE]

Thus, we have the first inequality

[TABLE]

Finally, for the second inequality,

[TABLE]

∎

Proposition 4.2**.**

Under Assumptions 2 and 3, the solution of $(\widehat{\mathbf{A}},\widehat{\mathbf{L}})$ to the optimization problem (21) satisfies

[TABLE]

Proof. This follows directly from Proposition 4.1, Assumption 3, and applying Cauchy-Schwarz twice

[TABLE]

∎

B-B Theorem

Now we restate the theorem again for convenience.

Theorem 3**.**

Under Assumptions 1-3, the solution to the optimization problem (21) satisfies

[TABLE]

Proof. Starting with Proposition 4.1, since our solution is in this restricted set, we can use the stronger convexity condition implied by Assumption 1,

[TABLE]

Revisiting the objective functional at optimality and skipping repetitive algebra (see proof for Proposition 4.1),

[TABLE]

Now from Proposition 4.2, we have

[TABLE]

Combining the previous two inequalities, we have

[TABLE]

where

[TABLE]

which is a conic in standard form, $\lambda^{\prime}=\lambda\sqrt{r_{\mathbf{L}}}$, and the entries of $\mathbf{Q}$ follow from taking $\mu$ given in Assumption 3.

Checking its discriminant,

[TABLE]

so we have that the conic equation describes an ellipse, and there are bounds on the values of $\|\boldsymbol{\Phi}\|_{F}$ and $\|\boldsymbol{\Psi}\|_{F}$.

For these individual bounds, we consider the points at which the gradients of $\boldsymbol{\zeta}^{\top}\mathbf{Q}\boldsymbol{\zeta}$ vanish w.r.t. each of $\|\boldsymbol{\Phi}\|_{F}$ and $\|\boldsymbol{\Psi}\|_{F}$. For $\|\boldsymbol{\Phi}\|_{F}$,

[TABLE]

Plugging this into the equation defined by $\boldsymbol{\zeta}^{\top}\mathbf{Q}\boldsymbol{\zeta}=0$ yields

[TABLE]

Since we are seeking an upper bound for $\|\boldsymbol{\Phi}\|_{F}$, it can be seen that we take the positive root,

[TABLE]

Similarly, for $\|\boldsymbol{\Psi}\|_{F}$,

[TABLE]

Finally, plugging this into $\boldsymbol{\zeta}^{\top}\mathbf{Q}\boldsymbol{\zeta}=0$ and solving for the upper bound yields

[TABLE]

∎