Generalized Score Functions for Causal Discovery

Huang et al. (2018). Generalized score functions for causal discovery.

A score-based algorithm requires two components: the score function and the search procedure. In addition to Markov assumptions, an ideal score function should meet the basic criteria as following

• characterization of conditional dependency
• local consistency and score equivalence
• robustness to multiple, if not all, data types and distributions, including
• non-linear relations
• non-Gaussian distributions
• (infinite) multi-dimensionality
• heterogeneity (i.e., continuous/discreet, stationary/time-series data)
• tractability and computational efficiency

While the previous algorithms are limited in assuming a certain model class or distribution to begin with, Generalized Score Functions (GSF) is a novel effort to unify score-based approaches into a generalized framework that can work on a wide class of data types and distributions.

Methods (refs)	Non linearity	Non Gaussian	Multi dimensionality	Heterogeneous data
BIC, MDL BGe, BDeu BDe, BCCD	No	No	No	No
Mixed BIC	No	No	No	Stationary only
LiNGAM ICA-LiNGAM Direct-LiNGAM	No	Yes	No	Yes
ANM CAM PNL	Yes	Yes	No	Yes
GSF	Yes	Yes	Yes	Stationary only

Notations

Let $X$ be a random variable with domain $\mathcal{X}$ in any dimension.
${\mathcal{H}}_X$ is defined as a Reproducing Kernel Hilbert Space (RKHS) on $\mathcal{X}$, in which

1. ${\varphi }_X:\mathcal{X}⟼{\mathcal{H}}_X$ is a continuous feature map projecting a point $x\in \mathcal{X}$ onto ${\mathcal{H}}_X$


2. $k_X:\mathcal{X}\times \mathcal{X}⟼\mathbb{R}$ is a positive definite kernel function such that $$k_X(x,x^{\prime })=⟨{\varphi }_X(x),{\varphi }_X(x^{\prime })⟩=⟨k_X(.,x),k_X(.,x^{\prime })⟩$$
3. probabilty law of $X$ is $P_X$ and ${\mathcal{H}}_X\subset L^2(P_X)$ - the space of square integrable functions
$\to$ This is a property of completeness, basically saying that with the probability of existence $P_X$, ${\mathcal{H}}_X$ contain functions that are square integrable. In ${\mathcal{H}}_X$, $L^2$ norm is defined.


4. for any point $x$ in a sample of $n$ observations,
   • the reproducing property yields $$\mathrm{\forall }x\in \mathcal{X},\mathrm{\forall }f\in {\mathcal{H}}_{\mathcal{X}},⟨f(.),k_X(.,x)⟩=f(x)=\sum _{i=1}^n{k}_X(x_i,x)c_i$$ with $c_i$ being the weight
   • its empirical feature map $k_x=[k_X(x_1,x),...,k_X(x_n,x){]}^T$


5. $K_X$ represents the $n\times n$ kernel matrix with each element $K_{ij}=k_X(x_i,x_j)$, and $\widetilde{K_X}=H{K}_{X}H$ is the centralized kernel matrix where $H=I-\frac{1}{n}\mathbb{1}{\mathbb{1}}^T$ with $n\times n$ idenity matrix $I$ and vector of 1's $\mathbb{1}$

Similar notations are applied to other variables $Y$ and $Z$

1. Characterization of Conditional Independence

The score function must act as a proxy for conditional independence test, which, in the general case, is characterized by conditional covariance operators. The kernel-based conditional independence test is described as follow

where ${\mathrm{\Sigma }}_{XX|Z}={\mathrm{\Sigma }}_{XX}-{\mathrm{\Sigma }}_{XZ}{\mathrm{\Sigma }}_{ZZ}^{-1}{\mathrm{\Sigma }}_{ZX}$.

For all functions $g\in {\mathcal{H}}_{\mathcal{X}}$ and $f\in {\mathcal{H}}_{\mathcal{Z}}$, the cross-covariance operator ${\mathrm{\Sigma }}_{ZX}:{\mathcal{H}}_{\mathcal{X}}⟼{\mathcal{H}}_{\mathcal{Z}}$ of a random vector ($X,Z$) on $\mathcal{X}\times \mathcal{Z}$ defined in Fukumizu, Bach & Jordan (2004) is $$⟨f,{\mathrm{\Sigma }}_{ZX}g{⟩}_{{\mathcal{H}}_{\mathcal{Z}}}=E_{XZ}[g(X)f(Z)]-E_X[g(X)]E_Z[f(Z)]$$

Additionally, the operator ${\mathrm{\Sigma }}_{XX|Z}$ captures conditional variance of a random variable $X$ through $$⟨g,{\mathrm{\Sigma }}_{XX|Z}g{⟩}_{{\mathcal{H}}_{\mathcal{X}}}=E_Z[Va{r}_{X|Z}[g(X)|Z]]$$

Let $\ddot{Z}:=[Y,Z]$. Similarly, we have $$⟨g,{\mathrm{\Sigma }}_{XX|\ddot{Z}}g{⟩}_{{\mathcal{H}}_{\mathcal{X}}}=E_{\ddot{Z}}[Va{r}_{X|\ddot{Z}}[g(X)|\ddot{Z}]]$$

2. Regression in RKHS

Given random variables $X,Z$ over measureable spaces $\mathcal{X},\mathcal{Z}$, the authors exploit regression framework in RKHS to project $Z$ onto finite-dimensional ${\mathcal{H}}_{\mathcal{X}}$, and encode the dependency by regressing ${\varphi }_X(X)$ on $Z$ $${\varphi }_X(X)=F(Z)+U,F:\mathcal{Z}⟼{\mathcal{H}}_X(1)$$ Consider the two regression models in RKHS $${\varphi }_X(X)=F_1(Z)+U_1,{\varphi }_X(X)=F_2(\ddot{Z})+U_2(2)$$ From Lemma 1, we can derive the equality between the conditional variance operators as

${\varphi }_X(X)$ is a vector in infinite-dimensional space, and can be written as ${\varphi }_X(X)=[{\varphi }_1(X),...,{\varphi }_i(X),...{]}^T$. Since ${\varphi }_X(X)$ is a feature map in ${\mathcal{H}}_X$, there exists a function $g\in {\mathcal{H}}_X$ such that $g={\varphi }_i(X)$. Thus, Equation 4 holds for any ${\varphi }_i(X)$.

Furthermore, each component ${\varphi }_i(X)$ and ${\varphi }_j(X)$ is orthogonal, and $Cov({\varphi }_i(X),{\varphi }_j(X))=0$ for $i\ne j$. Notice that ${\varphi }_X(X)|Z$ is also a random vector with ${\varphi }_X(X)|Z=[{\varphi }_1(X)|Z,...,{\varphi }_i(X)|Z,...{]}^T$, the covariance matrix of $\mathsf{Var}[{\varphi }_X(X)|Z]$ is a diagonal matrix in which each non-zero element corresponds to $\mathsf{Var}[{\varphi }_i(X)]$ Thus, we can derive

Given $Z$, $Y$ gives no information to predict $X$. Graphically, if $Z$ is a parent of $X$ and $X$ and $Y$ are independent conditioning on $Z$, there should not be any arrow directly pointing to $X$ from $Y$. Thus, the first structural equation in (2) is a better fit. Here, causal learning can be viewed as a problem of model selection, and the task becomes how to construct a score function for selection regression model in RKHS.

The most straightforward test for goodness-of-fit is the maximum log-likelihood for the regression problem in (1). Unfortunately, we do not know the distribution of ${\varphi }_X(X)$ since $\mathcal{X}$ is infinite-dimensional. The strategy is to reformulate the problem so that it is expressed only in kernel functions, which is tractable in finite space.

3. Formulation

For a particular observation $(x,z)$, we map ${\varphi }_X(x)$ into the empirical feature map $k_X$ that is in $n$dimensional space $$k_x=\left[\begin{array}{c}⟨k_X(x_1,x),{\varphi }_X(x){⟩}_{{\mathcal{H}}_X}\\ ...\\ ⟨k_X(x_n,x),{\varphi }_X(x){⟩}_{{\mathcal{H}}_X}\end{array}\right]$$

According the reproducing property, $k_x$ has the same vector structure as ${\varphi }_X(X)$ wherein each element is equivalent to a function $f_i(x)$ given by the linear combination. This means that regressing $k_X$ on $Z$ does not cause loss of information, that is the log-likelihood score reflects model compatibility while remains tractable.

We now turn our attention to the reformulated problem $$k_x=\tilde{F}(z)+\tilde{U}(5)$$ Assuming Gaussian errors, we fit a kernel ridge regression on finite sample size (regression with $L2$ regularization) to estimate $${\tilde{\hat{F}}}_Z(z)\approx \sum _{i=1}^n{k}_Z(z_i,z)c_i$$ where the coefficient matrix $c\in R^n$ can be found by solving $$min\frac{1}{n}||Y-c{K}_Z|{|}_{R^n}^2+\lambda c^T{K}_{Z}c\lambda :\text{regularization parameter}$$ This gives $$\hat{c}=Y(K_Z+\lambda nI{)}^{-1}$$ The fitted function is $${\tilde{\hat{F}}}_Z(z)=\hat{c}{K}_Z=\widetilde{K_X}(K_Z+\lambda nI{)}^{-1}{K}_Z$$ The estimated covariance matrix is $$\hat{\mathrm{\Sigma }}=\frac{1}{n}({\tilde{K}}_X-{\tilde{\hat{F}}}_Z)({\tilde{K}}_X-{\tilde{\hat{F}}}_Z{)}^T=n{\lambda }^2{\tilde{K}}_X({\tilde{K}}_Z+n\lambda I{)}^{-2}{\tilde{K}}_X$$

The maximal value of log-likelihood is represented as

From line 2 to line 3, $$\mathsf{\text{trace}}\{(\widetilde{K_X}-\hat{\tilde{W}}{\mathrm{\Phi }}_Z{)}^T{\hat{\mathrm{\Sigma }}}^{-1}(\widetilde{K_X}-\hat{\tilde{W}}{\mathrm{\Phi }}_Z)\}=\mathsf{\text{trace}}\{{\hat{\mathrm{\Sigma }}}^{-1}(\widetilde{K_X}-\hat{\tilde{W}}{\mathrm{\Phi }}_Z{)}^T(\widetilde{K_X}-\hat{\tilde{W}}{\mathrm{\Phi }}_Z)\}$$ $$=\mathsf{\text{trace}}\{{\hat{\mathrm{\Sigma }}}^{-1}(\widetilde{K_X}-{\tilde{\hat{F}}}_Z)(\widetilde{K_X}-{\tilde{\hat{F}}}_Z{)}^T\}=\mathsf{\text{trace}}\{\hat{\mathrm{\Sigma }}{\hat{\mathrm{\Sigma }}}^{-1}\}=\mathsf{\text{trace}}\{I\}=n$$

4. Score Functions

Maximizing the likelihood can cause over-fitting. Thus, the proposed score functions are Cross-Validation Log-Likelihood and Marginal Log-Likelihood.

4.1. CV Log-Likelihood

Given a dataset $D$ with $m$ variables $X_1,X_2,...,X_m$, we split $D$ into a training set of size $n_1$ and a test set of size $n_0$, then repeat the procedure $Q$ times. Let $D_{1,i}^{(q)}$ and $D_{0,i}^{(q)}$ be the $q^{th}$ training set and the $q^{th}$ training set for variable $X_i$.

The score of DAG ${\mathcal{G}}^{\mathcal{h}}$ is represented by $$S_{CV}({\mathcal{G}}^{\mathcal{h}};D)=\sum _{i=1}^m{S}_{CV}(X_i,P{A}_i^{{\mathcal{G}}^{\mathcal{h}}})$$ $$S_{CV}(X_i,P{A}_i^{{\mathcal{G}}^{\mathcal{h}}})=\frac{1}{Q}\sum _{q=1}^{Q}l(\widehat{\tilde{F_i^{(q)}}}|D_{0,i}^{(q)})$$ $X_i$ is the target variable and $P{A}_i^{{\mathcal{G}}^{\mathcal{h}}}$ is the set of predictors in the regression function in Equation 5. The individual likelihood function is represented with kernel matrices as follow

The formula can be derived using matrix manipulation. Notice that line 2 approximates ${\hat{\mathrm{\Sigma }}}^{(q)}$ by adding $\lambda$ to give rise to kernel functions. The paper also proposes removing the fitted values and using only $\widetilde{K_Z}$ to reduce overfitting. This is equivalent to regressing $k_X$ on the empty set. Please refer to the paper for more details.

4.2. Marginal Likelihood

The marginal likelihood score is estimated as $$S_M({\mathcal{G}}^{\mathcal{h}};D)=\sum _{i=1}^m{S}_M(X_i,P{A}_i^{{\mathcal{G}}^{\mathcal{h}}})$$ with $S_M(X_i,P{A}_i^{{\mathcal{G}}^{\mathcal{h}}})=\mathsf{log}p(X_i|P{A}_i^{{\mathcal{G}}^{\mathcal{h}}},\widehat{{\sigma }_i})$, and the noise variance $\widehat{{\sigma }_i}$ is learned by maximizing marginal likelihood with gradient methods.

Reliability

A score function must also satisfy Local Consistency and Score Equivalence.

Local Consistency

The score function is expected to give lower scores to models assigning with incorrect indendency constraints. It has been shown that under mild conditions, CV likelihood score is locally consistent according to Lemma 2.

The model selected by CV likelihood has the same performance with the model selected by the optimal benchmark that is fitted on the entire sample. That is, CV likelihood score supports learning causal relations that hold in the true data generating distribution. This conclusion at the same time supports $K$-fold cross validation with reasonably large $K$.

Marginal likelihood score is also proven to yield local consistency according to Lemma 3.

The score can then be expressed in terms of Bayesian information criterion (BIC), plus a constant, in which BIC is locally consistent.

Score Equivalence

On the other hand, the previous studies provide evidence that non-linear regression gives different scores to models in an equivalence class. Thus, neither CV or Marginal likelihood meet Score Equivalence.

However, the authors show that we can exploit local consistency and local score change to compare the reliability of any two DAGs. It is observed that the learned DAG with the same structure and orientation as the data generative distribution gives the highest scores. However, since there have not been any theoretical proof for this yet, we resort to the problem of searching the optimal Markov equivalence class.

In general, if conditions in Lemma 2 and Lemma 3 hold, the CV likelihood and Marginal likelihood, under RKHS regression framework, are guaranteed to find the Markov equivalence class that is closely consistent with the data generating distribution by using Greedy Equivalence Search algorithm (GES). Details can be found in Proposition 1 and Proposition 2.

Previous score-based methods

• BIC, MDL : Estimating the Dimension of a Model
• BGe : Learning Gaussian Networks
• BDeu, BDe : Learning Bayesian networks: The combination of knowledge and statistical data.
• BCCD : A Bayesian Approach to Constraint Based Causal Inference
• Mixed BIC : Causal discovery from databases with discrete and continuous variables.
• LiNGAM : LiNGAM: non-Gaussian methods for estimating causal structure
• ICA-LiNGAM : A linear non-Gaussian acyclic model for causal discovery
• Direct-LiNGAM : DirectLiNGAM: A direct method for learning a linear non-Gaussian structural equation model
• ANM : Nonlinear causal discovery with additive noise models
• CAM : CAM: Causal Additive Models, highdimensional order search and penalized regression
• PNL : On the identifiability of the post-nonlinear causal model