Search

Graphical models have found widespread applications in many areas of modern statistics and machine learning. Iterative Proportional Fitting (IPF) and its variants have become the default method for undirected graphical model estimation, and are thus ubiquitous in the field. As the IPF is an iterative approach, it is not always readily scalable to modern high-dimensional data regimes. In this paper we propose a novel and fast non-iterative method for positive definite graphical model estimation in high dimensions, one that directly addresses the shortcomings of IPF and its variants. In addition, the proposed method has a number of other attractive properties. First, we show formally that as the dimension p grows, the proportion of graphs for which the proposed method will outperform the state-of-the-art in terms of computational complexity and performance tends to 1, affirming its efficacy in modern settings. Second, the proposed approach can be readily combined with scalable non-iterative thresholding-based methods for high-dimensional sparsity selection. Third, the proposed method has high-dimensional statistical guarantees. Moreover, our numerical experiments also show that the proposed method achieves scalability without compromising on statistical precision. Fourth, unlike the IPF, which depends on the Gaussian likelihood, the proposed method is much more robust.

Knowledge graphs (KGs) are known for their large scale and knowledge inference ability, but are also notorious for the incompleteness associated with them. Due to the long-tail distribution of the relations in KGs, few-shot KG completion has been proposed as a solution to alleviate incompleteness and expand the coverage of KGs. It aims to make predictions for triplets involving novel relations when only a few training triplets are provided as reference. Previous methods have mostly focused on designing local neighbor aggregators to learn entity-level information and/or imposing sequential dependency assumption at the triplet level to learn meta relation information. However, valuable pairwise triplet-level interactions and context-level relational information have been largely overlooked for learning meta representations of few-shot relations. In this paper, we propose a hierarchical relational learning method (HiRe) for few-shot KG completion. By jointly capturing three levels of relational information (entity-level, triplet-level and context-level), HiRe can effectively learn and refine the meta representation of few-shot relations, and consequently generalize very well to new unseen relations. Extensive experiments on two benchmark datasets validate the superiority of HiRe against other state-of-the-art methods., Comment: 10 pages, 5 figures

Many applications benefit from theory relevant to the identification of variables having large correlations or partial correlations in high dimension. Recently there has been progress in the ultra-high dimensional setting when the sample size $n$ is fixed and the dimension $p$ tends to infinity. Despite these advances, the correlation screening framework suffers from practical, methodological and theoretical deficiencies. For instance, previous correlation screening theory requires that the population covariance matrix be sparse and block diagonal. This block sparsity assumption is however restrictive in practical applications. As a second example, correlation and partial correlation screening requires the estimation of dependence measures, which can be computationally prohibitive. In this paper, we propose a unifying approach to correlation and partial correlation mining that is not restricted to block diagonal correlation structure, thus yielding a methodology that is suitable for modern applications. By making connections to random geometric graphs, the number of highly correlated or partial correlated variables are shown to have compound Poisson finite-sample characterizations, which hold for both the finite $p$ case and when $p$ tends to infinity. The unifying framework also demonstrates a duality between correlation and partial correlation screening with theoretical and practical consequences., Comment: 70 pages, 9 figures

A surprising result of FitzGerald and Horn (1977) shows that $A^{\circ \alpha} := (a_{ij}^\alpha)$ is positive semidefinite (p.s.d.) for every entrywise nonnegative $n \times n$ p.s.d. matrix $A = (a_{ij})$ if and only if $\alpha$ is a positive integer or $\alpha \geq n-2$. Given a graph $G$, we consider the refined problem of characterizing the set $\mathcal{H}_G$ of entrywise powers preserving positivity for matrices with a zero pattern encoded by $G$. Using algebraic and combinatorial methods, we study how the geometry of $G$ influences the set $\mathcal{H}_G$. Our treatment provides new and exciting connections between combinatorics and analysis, and leads us to introduce and compute a new graph invariant called the critical exponent., Comment: 12 pages, final version. This is an extended abstract of arXiv:1504.04069 in FPSAC 2017

Simple random walks are a basic staple of the foundation of probability theory and form the building block of many useful and complex stochastic processes. In this paper we study a natural generalization of the random walk to a process in which the allowed step sizes take values in the set $\{\pm1,\pm2,\ldots,\pm k\}$, a process we call a random leap. The need to analyze such models arises naturally in modern-day data science and so-called "big data" applications. We provide closed-form expressions for quantities associated with first passage times and absorption events of random leaps. These expressions are formulated in terms of the roots of the characteristic polynomial of a certain recurrence relation associated with the transition probabilities. Our analysis shows that the expressions for absorption probabilities for the classical simple random walk are a special case of a universal result that is very elegant. We also consider an important variant of a random leap: the reflecting random leap. We demonstrate that the reflecting random leap exhibits more interesting behavior in regard to the existence of a stationary distribution and properties thereof. Questions relating to recurrence/transience are also addressed, as well as an application of the random leap.

Bayesian shrinkage methods have generated a lot of recent interest as tools for high-dimensional regression and model selection. These methods naturally facilitate tractable uncertainty quantification and incorporation of prior information. A common feature of these models, including the Bayesian lasso, global-local shrinkage priors, and spike-and-slab priors is that the corresponding priors on the regression coefficients can be expressed as scale mixture of normals. While the three-step Gibbs sampler used to sample from the often intractable associated posterior density has been shown to be geometrically ergodic for several of these models (Khare and Hobert, 2013; Pal and Khare, 2014), it has been demonstrated recently that convergence of this sampler can still be quite slow in modern high-dimensional settings despite this apparent theoretical safeguard. We propose a new method to draw from the same posterior via a tractable two-step blocked Gibbs sampler. We demonstrate that our proposed two-step blocked sampler exhibits vastly superior convergence behavior compared to the original three- step sampler in high-dimensional regimes on both real and simulated data. We also provide a detailed theoretical underpinning to the new method in the context of the Bayesian lasso. First, we derive explicit upper bounds for the (geometric) rate of convergence. Furthermore, we demonstrate theoretically that while the original Bayesian lasso chain is not Hilbert-Schmidt, the proposed chain is trace class (and hence Hilbert-Schmidt). The trace class property has useful theoretical and practical implications. It implies that the corresponding Markov operator is compact, and its eigenvalues are summable. It also facilitates a rigorous comparison of the two-step blocked chain with "sandwich" algorithms which aim to improve performance of the two-step chain by inserting an inexpensive extra step., Comment: This paper is already available at arXiv with identifier arXiv:1509.03697, and was submitted again due to some confusion

Covariance estimation for high-dimensional datasets is a fundamental problem in modern day statistics with numerous applications. In these high dimensional datasets, the number of variables p is typically larger than the sample size n. A popular way of tackling this challenge is to induce sparsity in the covariance matrix, its inverse or a relevant transformation. In particular, methods inducing sparsity in the Cholesky pa- rameter of the inverse covariance matrix can be useful as they are guaranteed to give a positive definite estimate of the covariance matrix. Also, the estimated sparsity pattern corresponds to a Directed Acyclic Graph (DAG) model for Gaussian data. In recent years, two useful penalized likelihood methods for sparse estimation of this Cholesky parameter (with no restrictions on the sparsity pattern) have been developed. How- ever, these methods either consider a non-convex optimization problem which can lead to convergence issues and singular estimates of the covariance matrix when p > n, or achieve a convex formulation by placing a strict constraint on the conditional variance parameters. In this paper, we propose a new penalized likelihood method for sparse estimation of the inverse covariance Cholesky parameter that aims to overcome some of the shortcomings of current methods, but retains their respective strengths. We ob- tain a jointly convex formulation for our objective function, which leads to convergence guarantees, even when p > n. The approach always leads to a positive definite and symmetric estimator of the covariance matrix. We establish high-dimensional estima- tion and graph selection consistency, and also demonstrate finite sample performance on simulated/real data.

We prove a refinement of the inequality by Hoffmann-Jorgensen that is significant for three reasons. First, our result improves on the state-of-the-art even for real-valued random variables. Second, the result unifies several versions in the Banach space literature, including those by Johnson and Schechtman [Ann. Probab. 17 (1989)], Klass and Nowicki [Ann. Probab. 28 (2000)], and Hitczenko and Montgomery-Smith [Ann. Probab. 29 (2001)]. Finally, we show that the Hoffmann-Jorgensen inequality (including our generalized version) holds not only in Banach spaces but more generally, in a very primitive mathematical framework required to state the inequality: a metric semigroup $\mathscr{G}$. This includes normed linear spaces as well as all compact, discrete, or (connected) abelian Lie groups., Comment: 11 pages, published in the Annals of Probability. The Introduction section shares motivating examples with arXiv:1506.02605

In this paper, we exploit the theory of dense graph limits to provide a new framework to study the stability of graph partitioning methods, which we call structural consistency. Both stability under perturbation as well as asymptotic consistency (i.e., convergence with probability $1$ as the sample size goes to infinity under a fixed probability model) follow from our notion of structural consistency. By formulating structural consistency as a continuity result on the graphon space, we obtain robust results that are completely independent of the data generating mechanism. In particular, our results apply in settings where observations are not independent, thereby significantly generalizing the common probabilistic approach where data are assumed to be i.i.d. In order to make precise the notion of structural consistency of graph partitioning, we begin by extending the theory of graph limits to include vertex colored graphons. We then define continuous node-level statistics and prove that graph partitioning based on such statistics is consistent. Finally, we derive the structural consistency of commonly used clustering algorithms in a general model-free setting. These include clustering based on local graph statistics such as homomorphism densities, as well as the popular spectral clustering using the normalized Laplacian. We posit that proving the continuity of clustering algorithms in the graph limit topology can stand on its own as a more robust form of model-free consistency. We also believe that the mathematical framework developed in this paper goes beyond the study of clustering algorithms, and will guide the development of similar model-free frameworks to analyze other procedures in the broader mathematical sciences., Comment: 33 pages

We introduce PseudoNet, a new pseudolikelihood-based estimator of the inverse covariance matrix, that has a number of useful statistical and computational properties. We show, through detailed experiments with synthetic and also real-world finance as well as wind power data, that PseudoNet outperforms related methods in terms of estimation error and support recovery, making it well-suited for use in a downstream application, where obtaining low estimation error can be important. We also show, under regularity conditions, that PseudoNet is consistent. Our proof assumes the existence of accurate estimates of the diagonal entries of the underlying inverse covariance matrix; we additionally provide a two-step method to obtain these estimates, even in a high-dimensional setting, going beyond the proofs for related methods. Unlike other pseudolikelihood-based methods, we also show that PseudoNet does not saturate, i.e., in high dimensions, there is no hard limit on the number of nonzero entries in the PseudoNet estimate. We present a fast algorithm as well as screening rules that make computing the PseudoNet estimate over a range of tuning parameters tractable.

Bayesian shrinkage methods have generated a lot of recent interest as tools for high-dimensional regression and model selection. These methods naturally facilitate tractable uncertainty quantification and incorporation of prior information. This benefit has led to extensive use of the Bayesian shrinkage methods across diverse applications. A common feature of these models is that the corresponding priors on the regression coefficients can be expressed as scale mixture of normals. While the three-step Gibbs sampler used to sample from the often intractable associated posterior density has been shown to be geometrically ergodic for several of these models, it has been demonstrated recently that convergence of this sampler can still be quite slow in modern high-dimensional settings despite this apparent theoretical safeguard. We propose a new method to draw from the same posterior via a tractable two-step blocked Gibbs sampler. We demonstrate that our proposed two-step blocked sampler exhibits vastly superior convergence behavior compared to the original three-step sampler in high-dimensional regimes on both real and simulated data. We also provide a detailed theoretical underpinning to the new method in the context of the Bayesian lasso. First, we prove that the proposed two-step sampler is geometrically ergodic, and derive explicit upper bounds for the (geometric) rate of convergence. Furthermore, we demonstrate theoretically that while the original Bayesian lasso chain is not Hilbert-Schmidt, the proposed chain is trace class (and hence Hilbert-Schmidt). The trace class property implies that the corresponding Markov operator is compact, and its (countably many) eigenvalues are summable. It also facilitates a rigorous comparison of the two-step blocked chain with "sandwich" algorithms which aim to improve performance of the two-step chain by inserting an inexpensive extra step.

Markov chain Monte Carlo (MCMC) lies at the core of modern Bayesian methodology, much of which would be impossible without it. Thus, the convergence properties of MCMCs have received significant attention, and in particular, proving (geometric) ergodicity is of critical interest. Trust in the ability of MCMCs to sample from modern-day high-dimensional posteriors, however, has been limited by a widespread perception that these chains typically experience serious convergence problems. In this paper, we first demonstrate that contemporary methods for obtaining convergence rates have serious limitations when the dimension grows. We then propose a framework for rigorously establishing the convergence behavior of commonly used high-dimensional MCMCs. In particular, we demonstrate theoretically the precise nature and severity of the convergence problems of popular MCMCs when implemented in high dimensions, including phase transitions in the convergence rates in various $n$ and $p$ regimes, and a universality result across an entire spectrum of models. We also show that convergence problems effectively eliminate the apparent safeguard of geometric ergodicity. We then demonstrate theoretical principles by which MCMCs can be constructed and analyzed to yield bounded geometric convergence rates even as the dimension $p$ grows without bound. Additionally, we propose a diagnostic tool for establishing convergence.

We study probability inequalities leading to tail estimates in a general semigroup $\mathscr{G}$ with a translation-invariant metric $d_{\mathscr{G}}$. (An important and central example of this in the functional analysis literature is that of $\mathscr{G}$ a Banach space.) Using our prior work [Ann. Prob. 2017] that extends the Hoffmann-Jorgensen inequality to all metric semigroups, we obtain tail estimates and approximate bounds for sums of independent semigroup-valued random variables, their moments, and decreasing rearrangements. In particular, we obtain the "correct" universal constants in several cases, extending results in the Banach space literature by Johnson-Schechtman-Zinn [Ann. Prob. 1985], Hitczenko [Ann. Prob. 1994], and Hitczenko and Montgomery-Smith [Ann. Prob. 2001]. Our results also hold more generally, in a very primitive mathematical framework required to state them: metric semigroups $\mathscr{G}$. This includes all compact, discrete, or (connected) abelian Lie groups., Comment: 13 pages, final version, published in Advances in Operator Theory

In this paper we develop a rigorous foundation for the study of integration and measures on the space $\mathscr{G}(V)$ of all graphs defined on a countable labelled vertex set $V$. We first study several interrelated $\sigma$-algebras and a large family of probability measures on graph space. We then focus on a "dyadic" Hamming distance function $\left\| \cdot \right\|_{\psi,2}$, which was very useful in the study of differentiation on $\mathscr{G}(V)$. The function $\left\| \cdot \right\|_{\psi,2}$ is shown to be a Haar measure-preserving bijection from the subset of infinite graphs to the circle (with the Haar/Lebesgue measure), thereby naturally identifying the two spaces. As a consequence, we establish a "change of variables" formula that enables the transfer of the Riemann-Lebesgue theory on $\mathbb{R}$ to graph space $\mathscr{G}(V)$. This also complements previous work in which a theory of Newton-Leibnitz differentiation was transferred from the real line to $\mathscr{G}(V)$ for countable $V$. Finally, we identify the Pontryagin dual of $\mathscr{G}(V)$, and characterize the positive definite functions on $\mathscr{G}(V)$., Comment: 15 pages, LaTeX

When can reliable inference be drawn in the "Big Data" context? This paper presents a framework for answering this fundamental question in the context of correlation mining, with implications for general large scale inference. In large scale data applications like genomics, connectomics, and eco-informatics the dataset is often variable-rich but sample-starved: a regime where the number $n$ of acquired samples (statistical replicates) is far fewer than the number $p$ of observed variables (genes, neurons, voxels, or chemical constituents). Much of recent work has focused on understanding the computational complexity of proposed methods for "Big Data." Sample complexity however has received relatively less attention, especially in the setting when the sample size $n$ is fixed, and the dimension $p$ grows without bound. To address this gap, we develop a unified statistical framework that explicitly quantifies the sample complexity of various inferential tasks. Sampling regimes can be divided into several categories: 1) the classical asymptotic regime where the variable dimension is fixed and the sample size goes to infinity; 2) the mixed asymptotic regime where both variable dimension and sample size go to infinity at comparable rates; 3) the purely high dimensional asymptotic regime where the variable dimension goes to infinity and the sample size is fixed. Each regime has its niche but only the latter regime applies to exa-scale data dimension. We illustrate this high dimensional framework for the problem of correlation mining, where it is the matrix of pairwise and partial correlations among the variables that are of interest. We demonstrate various regimes of correlation mining based on the unifying perspective of high dimensional learning rates and sample complexity for different structured covariance models and different inference tasks.

Bayesian inference for graphical models has received much attention in the literature in recent years. It is well known that when the graph G is decomposable, Bayesian inference is significantly more tractable than in the general non-decomposable setting. Penalized likelihood inference on the other hand has made tremendous gains in the past few years in terms of scalability and tractability. Bayesian inference, however, has not had the same level of success, though a scalable Bayesian approach has its respective strengths, especially in terms of quantifying uncertainty. To address this gap, we propose a scalable and flexible novel Bayesian approach for estimation and model selection in Gaussian undirected graphical models. We first develop a class of generalized G-Wishart distributions with multiple shape parameters for an arbitrary underlying graph. This class contains the G-Wishart distribution as a special case. We then introduce the class of Generalized Bartlett (GB) graphs, and derive an efficient Gibbs sampling algorithm to obtain posterior draws from generalized G-Wishart distributions corresponding to a GB graph. The class of Generalized Bartlett graphs contains the class of decomposable graphs as a special case, but is substantially larger than the class of decomposable graphs. We proceed to derive theoretical properties of the proposed Gibbs sampler. We then demonstrate that the proposed Gibbs sampler is scalable to significantly higher dimensional problems as compared to using an accept-reject or a Metropolis-Hasting algorithm. Finally, we show the efficacy of the proposed approach on simulated and real data.

The study of entrywise powers of matrices was originated by Loewner in the pursuit of the Bieberbach conjecture. Since the work of FitzGerald and Horn (1977), it is known that $A^{\circ \alpha} := (a_{ij}^\alpha)$ is positive semidefinite for every entrywise nonnegative $n \times n$ positive semidefinite matrix $A = (a_{ij})$ if and only if $\alpha$ is a positive integer or $\alpha \geq n-2$. This surprising result naturally extends the Schur product theorem, and demonstrates the existence of a sharp phase transition in preserving positivity. In this paper, we study when entrywise powers preserve positivity for matrices with structure of zeros encoded by graphs. To each graph is associated an invariant called its "critical exponent", beyond which every power preserves positivity. In our main result, we determine the critical exponents of all chordal/decomposable graphs, and relate them to the geometry of the underlying graphs. We then examine the critical exponent of important families of non-chordal graphs such as cycles and bipartite graphs. Surprisingly, large families of dense graphs have small critical exponents that do not depend on the number of vertices of the graphs., Comment: 21 pages; added Proposition 4.15 on correlation matrices. Final version, to appear in Journal of Combinatorial Theory, Series A

This paper proposes a general adaptive procedure for budget-limited predictor design in high dimensions called two-stage Sampling, Prediction and Adaptive Regression via Correlation Screening (SPARCS). SPARCS can be applied to high dimensional prediction problems in experimental science, medicine, finance, and engineering, as illustrated by the following. Suppose one wishes to run a sequence of experiments to learn a sparse multivariate predictor of a dependent variable $Y$ (disease prognosis for instance) based on a $p$ dimensional set of independent variables $\mathbf X=[X_1,\ldots, X_p]^T$ (assayed biomarkers). Assume that the cost of acquiring the full set of variables $\mathbf X$ increases linearly in its dimension. SPARCS breaks the data collection into two stages in order to achieve an optimal tradeoff between sampling cost and predictor performance. In the first stage we collect a few ($n$) expensive samples $\{y_i,\mathbf x_i\}_{i=1}^n$, at the full dimension $p\gg n$ of $\mathbf X$, winnowing the number of variables down to a smaller dimension $l < p$ using a type of cross-correlation or regression coefficient screening. In the second stage we collect a larger number $(t-n)$ of cheaper samples of the $l$ variables that passed the screening of the first stage. At the second stage, a low dimensional predictor is constructed by solving the standard regression problem using all $t$ samples of the selected variables. SPARCS is an adaptive online algorithm that implements false positive control on the selected variables, is well suited to small sample sizes, and is scalable to high dimensions. We establish asymptotic bounds for the Familywise Error Rate (FWER), specify high dimensional convergence rates for support recovery, and establish optimal sample allocation rules to the first and second stages., Comment: To appear in IEEE Transactions on Information Theory. 40 Pages. arXiv admin note: text overlap with arXiv:1303.2378

Concurrent time series commonly arise in various applications, including when monitoring the environment such as in air quality measurement networks, weather stations, oceanographic buoys, or in paleo form such as lake sediments, tree rings, ice cores, or coral isotopes, with each monitoring or sampling site providing one of the time series. The goal in such applications is to extract a common time trend or signal in the observed data. Other examples where the goal is to extract a common time trend for multiple time series are in stock price time series, neurological time series, and quality control time series. For this purpose we develop properties of MAF [Maximum Autocorrelation Factors] that linearly combines time series in order to maximize the resulting SNR [signal-to-noise-ratio] where there are multiple smooth signals present in the data. Equivalence is established in a regression setting between MAF and CCA [Canonical Correlation Analysis] even though MAF does not require specific signal knowledge as opposed to CCA. We proceed to derive the theoretical properties of MAF and quantify the SNR advantages of MAF in comparison with PCA [Principal Components Analysis], a commonly used method for linearly combining time series, and compare their statistical sample properties. MAF and PCA are then applied to real and simulated data sets to illustrate MAFs efficacy., Comment: 38 pages, 12 figures

High dimensional covariance estimation and graphical models is a contemporary topic in statistics and machine learning having widespread applications. An important line of research in this regard is to shrink the extreme spectrum of the covariance matrix estimators. A separate line of research in the literature has considered sparse inverse covariance estimation which in turn gives rise to graphical models. In practice, however, a sparse covariance or inverse covariance matrix which is simultaneously well-conditioned and at the same time computationally tractable is desired. There has been little research at the confluence of these three topics. In this paper we consider imposing a condition number constraint to various types of losses used in covariance and inverse covariance matrix estimation. When the loss function can be decomposed as a sum of an orthogonally invariant function of the estimate and its inner product with a function of the sample covariance matrix, we show that a solution path algorithm can be derived, involving a series of ordinary differential equations. The path algorithm is attractive because it provides the entire family of estimates for all possible values of the condition number bound, at the same computational cost of a single estimate with a fixed upper bound. An important finding is that the proximal operator for the condition number constraint, which turns out to be very useful in regularizing loss functions that are not orthogonally invariant and may yield non-positive-definite estimates, can be efficiently computed by this path algorithm. As a concrete illustration of its practical importance, we develop an operator-splitting algorithm that imposes a guarantee of well-conditioning as well as positive definiteness to recently proposed convex pseudo-likelihood based graphical model selection methods., Comment: 19 pages; 1 figure

Representing the conditional independences present in a multivariate random vector via graphs has found widespread use in applications, and such representations are popularly known as graphical models or Markov random fields. These models have many useful properties, but their fundamental attractive feature is their ability to reflect conditional independences between blocks of variables through graph separation, a consequence of the equivalence of the pairwise, local and global Markov properties demonstrated by Pearl and Paz (1985). Modern day applications often necessitate working with either an infinite collection of variables (such as in a spatial-temporal field) or approximating a large high-dimensional finite stochastic system with an infinite-dimensional system. However, it is unclear whether the conditional independences present in an infinite-dimensional random vector or stochastic process can still be represented by separation criteria in an infinite graph. In light of the advantages of using graphs as tools to represent stochastic relationships, we undertake in this paper a general study of infinite graphical models. First, we demonstrate that naive extensions of the assumptions required for the finite case results do not yield equivalence of the Markov properties in the infinite-dimensional setting, thus calling for a more in-depth analysis. To this end, we proceed to derive general conditions which do allow representing the conditional independence in an infinite-dimensional random system by means of graphs, and our results render the result of Pearl and Paz as a special case of a more general phenomenon. We conclude by demonstrating the applicability of our theory through concrete examples of infinite-dimensional graphical models., Comment: 43 pages, LaTeX

The study of very large graphs is a prominent theme in modern-day mathematics. In this paper we develop a rigorous foundation for studying the space of finite labelled graphs and their limits. These limiting objects are naturally countable graphs, and the completed graph space $\mathscr{G}(V)$ is identified with the 2-adic integers as well as the Cantor set. The goal of this paper is to develop a model for differentiation on graph space in the spirit of the Newton-Leibnitz calculus. To this end, we first study the space of all finite labelled graphs and their limiting objects, and establish analogues of left-convergence, homomorphism densities, a Counting Lemma, and a large family of topologically equivalent metrics on labelled graph space. We then establish results akin to the First and Second Derivative Tests for real-valued functions on countable graphs, and completely classify the permutation automorphisms of graph space that preserve its topological and differential structures., Comment: 21 pages. Final version, to appear in Advances in Operator Theory (the journal version is longer)

Sparse high dimensional graphical model selection is a popular topic in contemporary machine learning. To this end, various useful approaches have been proposed in the context of $\ell_1$-penalized estimation in the Gaussian framework. Though many of these inverse covariance estimation approaches are demonstrably scalable and have leveraged recent advances in convex optimization, they still depend on the Gaussian functional form. To address this gap, a convex pseudo-likelihood based partial correlation graph estimation method (CONCORD) has been recently proposed. This method uses coordinate-wise minimization of a regression based pseudo-likelihood, and has been shown to have robust model selection properties in comparison with the Gaussian approach. In direct contrast to the parallel work in the Gaussian setting however, this new convex pseudo-likelihood framework has not leveraged the extensive array of methods that have been proposed in the machine learning literature for convex optimization. In this paper, we address this crucial gap by proposing two proximal gradient methods (CONCORD-ISTA and CONCORD-FISTA) for performing $\ell_1$-regularized inverse covariance matrix estimation in the pseudo-likelihood framework. We present timing comparisons with coordinate-wise minimization and demonstrate that our approach yields tremendous payoffs for $\ell_1$-penalized partial correlation graph estimation outside the Gaussian setting, thus yielding the fastest and most scalable approach for such problems. We undertake a theoretical analysis of our approach and rigorously demonstrate convergence, and also derive rates thereof., Comment: NIPS accepted version

In recent years, a variety of useful extensions of the Wishart have been proposed in the literature for the purposes of studying Markov random fields/graphical models. In particular, generalizations of the Wishart, referred to as Type I and Type II Wishart distributions, have been introduced by Letac and Massam (\emph{Annals of Statistics} 2006) and play important roles in both frequentist and Bayesian inference for Gaussian graphical models. These distributions have been especially useful in high-dimensional settings due to the flexibility offered by their multiple shape parameters. The domain of In this paper we resolve a long-standing conjecture of Letac and Massam (LM) concerning the domains of the multi-parameters of graphical Wishart type distributions. This conjecture, posed in \emph{Annals of Statistics}, also relates fundamentally to the existence of Bayes estimators corresponding to these high dimensional priors. To achieve our goal, we first develop novel theory in the context of probabilistic analysis of graphical models. Using these tools, and a recently introduced class of Wishart distributions for directed acyclic graph (DAG) models, we proceed to give counterexamples to the LM conjecture, thus completely resolving the problem. Our analysis also proceeds to give useful insights on graphical Wishart distributions with implications for Bayesian inference for such models., Comment: 24 pages, 15 figures

Functions preserving Loewner positivity when applied entrywise to positive semidefinite matrices have been widely studied in the literature. Following the work of Schoenberg [Duke Math. J. 9], Rudin [Duke Math. J. 26], and others, it is well-known that functions preserving positivity for matrices of all dimensions are absolutely monotonic (i.e., analytic with nonnegative Taylor coefficients). In this paper, we study functions preserving positivity when applied entrywise to sparse matrices, with zeros encoded by a graph $G$ or a family of graphs $G_n$. Our results generalize Schoenberg and Rudin's results to a modern setting, where functions are often applied entrywise to sparse matrices in order to improve their properties (e.g. better conditioning). The only such result known in the literature is for the complete graph $K_2$. We provide the first such characterization result for a large family of non-complete graphs. Specifically, we characterize functions preserving Loewner positivity on matrices with zeros according to a tree. These functions are multiplicatively midpoint-convex and super-additive. Leveraging the underlying sparsity in matrices thus admits the use of functions which are not necessarily analytic nor absolutely monotonic. We further show that analytic functions preserving positivity on matrices with zeros according to trees can contain arbitrarily long sequences of negative coefficients, thus obviating the need for absolute monotonicity in a very strong sense. This result leads to the question of exactly when absolute monotonicity is necessary when preserving positivity for an arbitrary class of graphs. We then provide a stronger condition in terms of the numerical range of all symmetric matrices, such that functions satisfying this condition on matrices with zeros according to any family of graphs with unbounded degrees are necessarily absolutely monotonic., Comment: Final version (updated references + minor edits); 21 pages, to appear in Transactions of the AMS

Entrywise functions preserving the cone of positive semidefinite matrices have been studied by many authors, most notably by Schoenberg [Duke Math. J. 9, 1942] and Rudin [Duke Math. J. 26, 1959]. Following their work, it is well-known that entrywise functions preserving Loewner positivity in all dimensions are precisely the absolutely monotonic functions. However, there are strong theoretical and practical motivations to study functions preserving positivity in a fixed dimension $n$. Such characterizations for a fixed value of $n$ are difficult to obtain, and in fact are only known in the $2 \times 2$ case. In this paper, using a novel and intuitive approach, we study entrywise functions preserving positivity on distinguished submanifolds inside the cone obtained by imposing rank constraints. These rank constraints are prevalent in applications, and provide a natural way to relax the elusive original problem of preserving positivity in a fixed dimension. In our main result, we characterize entrywise functions mapping $n \times n$ positive semidefinite matrices of rank at most $l$ into positive semidefinite matrices of rank at most $k$ for $1 \leq l \leq n$ and $1 \leq k < n$. We also demonstrate how an important necessary condition for preserving positivity by Horn and Loewner [Trans. Amer. Math. Soc. 136, 1969] can be significantly generalized by adding rank constraints. Finally, our techniques allow us to obtain an elementary proof of the classical characterization of functions preserving positivity in all dimensions obtained by Schoenberg and Rudin., Comment: Published in Trans. Amer. Math. Soc. This later version has a few changes to Theorem B and related results (mostly to the statements). The edits made are listed on the final page

Several methods have been recently proposed for estimating sparse Gaussian graphical models using $\ell_{1}$ regularization on the inverse covariance matrix. Despite recent advances, contemporary applications require methods that are even faster in order to handle ill-conditioned high dimensional modern day datasets. In this paper, we propose a new method, G-AMA, to solve the sparse inverse covariance estimation problem using Alternating Minimization Algorithm (AMA), that effectively works as a proximal gradient algorithm on the dual problem. Our approach has several novel advantages over existing methods. First, we demonstrate that G-AMA is faster than the previous best algorithms by many orders of magnitude and is thus an ideal approach for modern high throughput applications. Second, global linear convergence of G-AMA is demonstrated rigorously, underscoring its good theoretical properties. Third, the dual algorithm operates on the covariance matrix, and thus easily facilitates incorporating additional constraints on pairwise/marginal relationships between feature pairs based on domain specific knowledge. Over and above estimating a sparse inverse covariance matrix, we also illustrate how to (1) incorporate constraints on the (bivariate) correlations and, (2) incorporate equality (equisparsity) or linear constraints between individual inverse covariance elements. Fourth, we also show that G-AMA is better adept at handling extremely ill-conditioned problems, as is often the case with real data. The methodology is demonstrated on both simulated and real datasets to illustrate its superior performance over recently proposed methods., Comment: 21 pages, 3 figures

Entrywise powers of matrices have been well-studied in the literature, and have recently received renewed attention in the regularization of high-dimensional correlation matrices. In this paper, we study powers of positive semidefinite block matrices $(H_{st})_{s,t=1}^n$ with complex entries. We first characterize the powers $\alpha\in\mathbb{R}$ such that the blockwise power map $(H_{st}) \mapsto (H_{st}^\alpha)$ preserves Loewner positivity. The characterization is obtained by exploiting connections with the theory of matrix monotone functions developed by Loewner. Second, we revisit previous work by Choudhury [Proc. AMS 108] who had provided a lower bound on $\alpha$ for preserving positivity when the blocks $H_{st}$ pairwise commute. We completely settle this problem by characterizing the full set of powers preserving positivity in this setting. Our characterizations generalize previous work by FitzGerald-Horn, Bhatia-Elsner, and Hiai from scalars to arbitrary block size, and in particular, generalize the Schur Product Theorem. Finally, a natural and unifying framework for studying the case of diagonalizable blocks consists of replacing real powers by general characters of the complex plane. We thus classify such characters, and generalize our results to this more general setting. In the course of our work, given $\beta\in\mathbb{Z}$, we provide lower and upper bounds for the threshold power $\alpha >0$ above which the complex characters $re^{i\theta}\mapsto r^\alpha e^{i\beta\theta}$ preserve positivity when applied entrywise to positive semidefinite matrices. In particular, we completely resolve the $n=3$ case of a question raised in 2001 by Xingzhi Zhan. As an application, we extend previous work by de Pillis [Duke Math. J. 36] by classifying the characters $K$ of the complex plane for which the map $(H_{st})_{s,t=1}^n \mapsto (K({\rm tr}(H_{st})))_{s,t=1}^n$ preserves positivity., Comment: 17 pages, LaTeX

We consider the general problem of minimizing an objective function which is the sum of a convex function (not strictly convex) and absolute values of a subset of variables (or equivalently the l1-norm of the variables). This problem appears exten- sively in modern statistical applications associated with high-dimensional data or "big data", and corresponds to optimizing l1-regularized likelihoods in the context of model selection. In such applications, cyclic coordinatewise minimization (CCM), where the objective function is sequentially minimized with respect to each individual coordi- nate, is often employed as it offers a computationally cheap and effective optimization method. Consequently, it is crucial to obtain theoretical guarantees of convergence for the sequence of iterates produced by the cyclic coordinatewise minimization in this setting. Moreover, as the objective corresponds to at l1-regularized likelihoods of many variables, it is important to obtain convergence of the iterates themselves, and not just the function values. Previous results in the literature only establish either, (i) that every limit point of the sequence of iterates is a stationary point of the objective function, or (ii) establish convergence under special assumptions, or (iii) establish con- vergence for a different minimization approach (which uses quadratic approximation based gradient descent followed by an inexact line search), (iv) establish convergence of only the function values of the sequence of iterates produced by random coordinatewise minimization (a variant of CCM). In this paper, a rigorous general proof of convergence for the cyclic coordinatewise minimization algorithm is provided. We demonstrate the usefulness of our general results in contemporary applications.

Recently, the theory of dense graph limits has received attention from multiple disciplines including graph theory, computer science, statistical physics, probability, statistics, and group theory. In this paper we initiate the study of the general structure of differentiable graphon parameters $F$. We derive consistency conditions among the higher G\^ateaux derivatives of $F$ when restricted to the subspace of edge weighted graphs $\mathcal{W}_{\bf p}$. Surprisingly, these constraints are rigid enough to imply that the multilinear functionals $\Lambda: \mathcal{W}_{\bf p}^n \to \mathbb{R}$ satisfying the constraints are determined by a finite set of constants indexed by isomorphism classes of multigraphs with $n$ edges and no isolated vertices. Using this structure theory, we explain the central role that homomorphism densities play in the analysis of graphons, by way of a new combinatorial interpretation of their derivatives. In particular, homomorphism densities serve as the monomials in a polynomial algebra that can be used to approximate differential graphon parameters as Taylor polynomials. These ideas are summarized by our main theorem, which asserts that homomorphism densities $t(H,-)$ where $H$ has at most $N$ edges form a basis for the space of smooth graphon parameters whose $(N+1)$st derivatives vanish. As a consequence of this theory, we also extend and derive new proofs of linear independence of multigraph homomorphism densities, and characterize homomorphism densities. In addition, we develop a theory of series expansions, including Taylor's theorem for graph parameters and a uniqueness principle for series. We use this theory to analyze questions raised by Lov\'asz, including studying infinite quantum algebras and the connection between right- and left-homomorphism densities., Comment: Final version (36 pages), accepted for publication in Journal of Combinatorial Theory, Series A

The application of the lasso is espoused in high-dimensional settings where only a small number of the regression coefficients are believed to be nonzero. Moreover, statistical properties of high-dimensional lasso estimators are often proved under the assumption that the correlation between the predictors is bounded. In this vein, coordinatewise methods, the most common means of computing the lasso solution, work well in the presence of low to moderate multicollinearity. The computational speed of coordinatewise algorithms degrades however as sparsity decreases and multicollinearity increases. Motivated by these limitations, we propose the novel "Deterministic Bayesian Lasso" algorithm for computing the lasso solution. This algorithm is developed by considering a limiting version of the Bayesian lasso. The performance of the Deterministic Bayesian Lasso improves as sparsity decreases and multicollinearity increases, and can offer substantial increases in computational speed. A rigorous theoretical analysis demonstrates that (1) the Deterministic Bayesian Lasso algorithm converges to the lasso solution, and (2) it leads to a representation of the lasso estimator which shows how it achieves both $\ell_1$ and $\ell_2$ types of shrinkage simultaneously. Connections to other algorithms are also provided. The benefits of the Deterministic Bayesian Lasso algorithm are then illustrated on simulated and real data.

Entrywise functions preserving Loewner positivity have been studied by many authors, most notably Schoenberg and Rudin. Following their work, it is known that functions preserving positivity when applied entrywise to positive semidefinite matrices of all dimensions are necessarily analytic with nonnegative Taylor coefficients. When the dimension is fixed, it has been shown by Vasudeva and Horn that such functions are automatically continuous and sufficiently differentiable. A natural refinement of the aforementioned problem consists of characterizing functions preserving positivity under rank constraints. In this paper, we begin this study by characterizing entrywise functions which preserve the cone of positive semidefinite real matrices of rank $1$ with entries in a general interval. Classifying such functions is intimately connected to the classical problem of solving Cauchy's functional equations, which have non-measurable solutions. We demonstrate that under mild local measurability assumptions, such functions are automatically smooth and can be completely characterized. We then extend our results by classifying functions preserving positivity on rank $1$ Hermitian complex matrices., Comment: 12 pages; LaTeX, significantly revised from previous version

Entrywise powers of symmetric matrices preserving positivity, monotonicity or convexity with respect to the Loewner ordering arise in various applications, and have received much attention recently in the literature. Following FitzGerald and Horn [J. Math. Anal. Appl., 1977], it is well-known that there exists a critical exponent beyond which all entrywise powers preserve positive definiteness. Similar phase transition phenomena have also recently been shown by Hiai (2009) to occur for monotonicity and convexity. In this paper, we complete the characterization of all the entrywise powers below and above the critical exponents that are positive, monotone, or convex on the cone of positive semidefinite matrices. We then extend the original problem by fully classifying the positive, monotone, or convex powers in a more general setting where additional rank constraints are imposed on the matrices. We also classify the entrywise powers that are super/sub-additive with respect to the Loewner ordering. Finally, we extend all the previous characterizations to matrices with negative entries. Our analysis consequently allows us to answer a question raised by Bhatia and Elsner (2007) regarding the smallest dimension for which even extensions of the power functions do not preserve Loewner positivity., Comment: 19 pages; LaTeX; added complete classification of powers preserving sub-additivity

Graphical models have proven to be powerful tools for representing high-dimensional systems of random variables. One example of such a model is the undirected graph, in which lack of an edge represents conditional independence between two random variables given the rest. Another example is the bidirected graph, in which absence of edges encodes pairwise marginal independence. Both of these classes of graphical models have been extensively studied, and while they are considered to be dual to one another, except in a few instances this duality has not been thoroughly investigated. In this paper, we demonstrate how duality between undirected and bidirected models can be used to transport results for one class of graphical models to the dual model in a transparent manner. We proceed to apply this technique to extend previously existing results as well as to prove new ones, in three important domains. First, we discuss the pairwise and global Markov properties for undirected and bidirected models, using the pseudographoid and reverse-pseudographoid rules which are weaker conditions than the typically used intersection and composition rules. Second, we investigate these pseudographoid and reverse pseudographoid rules in the context of probability distributions, using the concept of duality in the process. Duality allows us to quickly relate them to the more familiar intersection and composition properties. Third and finally, we apply the dualization method to understand the implications of faithfulness, which in turn leads to a more general form of an existing result.

Understanding centennial scale climate variability requires data sets that are accurate, long, continuous and of broad spatial coverage. Since instrumental measurements are generally only available after 1850, temperature fields must be reconstructed using paleoclimate archives, known as proxies. Various climate field reconstructions (CFR) methods have been proposed to relate past temperature to such proxy networks. In this work, we propose a new CFR method, called GraphEM, based on Gaussian Markov random fields embedded within an EM algorithm. Gaussian Markov random fields provide a natural and flexible framework for modeling high-dimensional spatial fields. At the same time, they provide the parameter reduction necessary for obtaining precise and well-conditioned estimates of the covariance structure, even in the sample-starved setting common in paleoclimate applications. In this paper, we propose and compare the performance of different methods to estimate the graphical structure of climate fields, and demonstrate how the GraphEM algorithm can be used to reconstruct past climate variations. The performance of GraphEM is compared to the widely used CFR method RegEM with regularization via truncated total least squares, using synthetic data. Our results show that GraphEM can yield significant improvements, with uniform gains over space, and far better risk properties. We demonstrate that the spatial structure of temperature fields can be well estimated by graphs where each neighbor is only connected to a few geographically close neighbors, and that the increase in performance is directly related to recovering the underlying sparsity in the covariance of the spatial field. Our work demonstrates how significant improvements can be made in climate reconstruction methods by better modeling the covariance structure of the climate field., Comment: Published at http://dx.doi.org/10.1214/14-AOAS794 in the Annals of Applied Statistics (http://www.imstat.org/aoas/) by the Institute of Mathematical Statistics (http://www.imstat.org)

Great strides have been made in the field of reconstructing past temperatures based on models relating temperature to temperature-sensitive paleoclimate proxies. One of the goals of such reconstructions is to assess if current climate is anomalous in a millennial context. These regression based approaches model the conditional mean of the temperature distribution as a function of paleoclimate proxies (or vice versa). Some of the recent focus in the area has considered methods which help reduce the uncertainty inherent in such statistical paleoclimate reconstructions, with the ultimate goal of improving the confidence that can be attached to such endeavors. A second important scientific focus in the subject area is the area of forward models for proxies, the goal of which is to understand the way paleoclimate proxies are driven by temperature and other environmental variables. In this paper we introduce novel statistical methodology for (1) quantile regression with autoregressive residual structure, (2) estimation of corresponding model parameters, (3) development of a rigorous framework for specifying uncertainty estimates of quantities of interest, yielding (4) statistical byproducts that address the two scientific foci discussed above. Our statistical methodology demonstrably produces a more robust reconstruction than is possible by using conditional-mean-fitting methods. Our reconstruction shares some of the common features of past reconstructions, but also gains useful insights. More importantly, we are able to demonstrate a significantly smaller uncertainty than that from previous regression methods. In addition, the quantile regression component allows us to model, in a more complete and flexible way than least squares, the conditional distribution of temperature given proxies. This relationship can be used to inform forward models relating how proxies are driven by temperature.

Sparse high dimensional graphical model selection is a topic of much interest in modern day statistics. A popular approach is to apply l1-penalties to either (1) parametric likelihoods, or, (2) regularized regression/pseudo-likelihoods, with the latter having the distinct advantage that they do not explicitly assume Gaussianity. As none of the popular methods proposed for solving pseudo-likelihood based objective functions have provable convergence guarantees, it is not clear if corresponding estimators exist or are even computable, or if they actually yield correct partial correlation graphs. This paper proposes a new pseudo-likelihood based graphical model selection method that aims to overcome some of the shortcomings of current methods, but at the same time retain all their respective strengths. In particular, we introduce a novel framework that leads to a convex formulation of the partial covariance regression graph problem, resulting in an objective function comprised of quadratic forms. The objective is then optimized via a coordinate-wise approach. The specific functional form of the objective function facilitates rigorous convergence analysis leading to convergence guarantees; an important property that cannot be established using standard results, when the dimension is larger than the sample size, as is often the case in high dimensional applications. These convergence guarantees ensure that estimators are well-defined under very general conditions, and are always computable. In addition, the approach yields estimators that have good large sample properties and also respect symmetry. Furthermore, application to simulated/real data, timing comparisons and numerical convergence is demonstrated. We also present a novel unifying framework that places all graphical pseudo-likelihood methods as special cases of a more general formulation, leading to important insights.

Doubly non-negative matrices arise naturally in many setting including Markov random fields (positively banded graphical models) and in the convergence analysis of Markov chains. In this short note, we settle a recent conjecture by C.R. Johnson et al. [Linear Algebra Appl. 435 (2011)] by proving that the critical exponent beyond which all continuous conventional powers of $n$-by-$n$ doubly nonnegative matrices are doubly nonnegative is exactly $n-2$. We show that the conjecture follows immediately by applying a general characterization from the literature. We prove a stronger form of the conjecture by classifying all powers preserving doubly nonnegative matrices, and proceed to generalize the conjecture for broad classes of functions. We also provide different approaches for settling the original conjecture., Comment: 6 pages; LaTeX; Accepted for publication in Linear Algebra and its Applications

We introduce a new approach to variable selection, called Predictive Correlation Screening, for predictor design. Predictive Correlation Screening (PCS) implements false positive control on the selected variables, is well suited to small sample sizes, and is scalable to high dimensions. We establish asymptotic bounds for Familywise Error Rate (FWER), and resultant mean square error of a linear predictor on the selected variables. We apply Predictive Correlation Screening to the following two-stage predictor design problem. An experimenter wants to learn a multivariate predictor of gene expressions based on successive biological samples assayed on mRNA arrays. She assays the whole genome on a few samples and from these assays she selects a small number of variables using Predictive Correlation Screening. To reduce assay cost, she subsequently assays only the selected variables on the remaining samples, to learn the predictor coefficients. We show superiority of Predictive Correlation Screening relative to LASSO and correlation learning (sometimes popularly referred to in the literature as marginal regression or simple thresholding) in terms of performance and computational complexity., Comment: 31 pages, 9 figures, Appearing in Proceedings of the 16th International Conference on Artificial Intelligence and Statistics (AISTATS)