A causal discovery method based on conditional independence testing and continuous optimization

Abstract

This invention belongs to the field of machine learning technology, specifically a causal discovery method based on conditional independence testing and continuous optimization. The method includes: establishing a masking matrix using conditional independence testing; designing constraint terms (hard / soft constraints) in a continuous optimization loss function based on the masking matrix; establishing a variational graph autoencoder model using a combination of variational autoencoders and structural causal models, and deriving the corresponding evidence lower bound objective function; and optimizing the equality constraint optimization problem using the quadratic penalty function method to obtain the weight matrix, thereby discovering causal relationships. The method has been validated on virtual datasets and real biological datasets. This invention overcomes the limitations of existing causal discovery methods based on continuous optimization in discovering causal relationships in scenarios with limited samples and heterogeneous noise.

Description

Technical Field

[0001] This invention belongs to the field of machine learning technology, specifically relating to a causal discovery method based on continuous optimization. Background Technology

[0002] In AI-driven legislative public opinion analysis, legislation, from a causal perspective, is an intervention. Therefore, essentially, all our public opinion analysis is based on the distribution of interventions, which falls within the scope of causal science. Furthermore, when we need to explain the results of public opinion predictions, analyzing from a causal rather than a correlational perspective can greatly improve the interpretability of the model. Making reliable, stable, and highly interpretable predictions requires us to characterize the causal relationships underlying the data; in other words, we need to identify causal relationships from these observational data, i.e., achieve causal discovery. In the field of artificial intelligence, causal discovery is becoming a very important problem and is beginning to demonstrate its value in multiple interdisciplinary fields (such as bioinformatics). Existing causal discovery methods can generally be divided into the following categories: constraint-based methods, scoring-based methods, function-based causal model methods, and hybrid methods. Here, we focus on scoring-based methods. The idea is to select a specific scoring function, score all candidate causal graphs, and select the one with the highest score as the output causal graph. The commonly chosen scoring function is the Bayesian information criterion score. However, as the number of candidate causal graphs increases exponentially with the number of vertices, theory confirms that finding the optimally scored causal graph is an NP-hard problem. Therefore, directly searching for the optimally scored graph is impractical when the number of vertices is large. Early solutions employed a greedy approach, seeking local optima by adding, deleting, or reversing edges. These traditional score-based causal discovery methods typically model the problem as a combinatorial optimization problem. In recent years, researchers have proposed a differentiable equivalent expression for discrete acyclic constraints, allowing the causal discovery problem to be solved in the continuous domain using numerical optimization algorithms, rather than being limited to the discrete domain.

[0003] F(W)

[0004] stG(W)∈DAGs

[0005] F(W)

[0006] sth(W) = 0

[0007] Where h(W)=tr(e) W°WThe formula ∫t-d (tr(·) represents the trace of the matrix, ° represents the Hadamard product of the matrices, and d is the size of the causal graph) is used to measure whether a graph is a directed acyclic graph (DAG). This transformation technique reduces the time cost of causal discovery algorithms and has broad integration potential with models that utilize gradient information for optimization, such as neural networks. Such methods are often referred to as causal discovery methods based on continuous optimization. The typical implementation involves designing a specific loss function F(W) and optimization method, and using a numerical optimization algorithm to solve an equality optimization problem, i.e., finding the solution W that minimizes F(W) while satisfying h(W) = 0. ★ The corresponding causal graph G(W) ★ (Generally, this is achieved by setting a specific threshold τ.) An edge from i to j is considered to exist if the absolute value of the edge is greater than a threshold τ (the threshold τ depends on the specific scenario, but is generally set to 0.3). However, this type of method still has many problems. First, its performance drops significantly when the sample size is small. Second, it performs poorly when dealing with heterogeneous noise. Furthermore, continuous optimization methods have been shown to be susceptible to the influence of numerical range; scaling the original data can alter the final optimization result, which greatly affects its usability in real-world scenarios. Summary of the Invention

[0008] The purpose of this invention is to provide a causal discovery method based on conditional independence testing and continuous optimization, which can achieve reliable causal discovery at a faster speed, especially avoiding the performance degradation of existing methods in scenarios with fewer samples and heterogeneous noise.

[0009] The causal discovery method based on conditional independence testing and continuous optimization provided by this invention employs several novel techniques, including: designing a masking matrix based on conditional independence (CI) testing and deriving a continuous optimization regularization term from it; designing a variational graph autoencoder model based on a combination of variational autoencoders and structural causal models to model the causal generation process; and reducing the entire problem to an equality-constrained optimization problem based on continuous optimization and optimizing it using a quadratic penalty function method. The specific steps are as follows:

[0010] (I) Modeling the causal generation process; Modeling the causal generation process using a variational graph autoencoder model; In this process, variational inference techniques are introduced, and variational parameters and true causal graph parameters are learned by minimizing the KL divergence between the variational distribution and the posterior distribution of the true latent variables; In order to ensure that the optimized weight matrix can derive an acyclic causal graph, a differentiable constraint is introduced to restrict the search space to a directed acyclic graph space.

[0011] (ii) Construction of conditional independence constraints: The image is preprocessed to obtain smaller image patches, and the image patch features are extracted using a pre-trained image feature extraction model, and finally a set of image patch representations are output.

[0012] (III) Equality constraint optimization: The objective function of the variational graph autoencoder can be combined with the conditional independence constraint, and can be expressed as an equality constraint optimization problem. The parameters are optimized by using the quadratic penalty function method to obtain the final causal graph.

[0013] The causal generation process modeling described in step (I) is as follows:

[0014] Suppose that the structural causal model behind the observation data satisfies the following form:

[0015] X=f2((IW T ) -1 f1(Z)), (1)

[0016] Where X is the observed variable, and W is the causal weight matrix (W0). i,j ≠0 corresponds to an edge i→j in the causal graph G, i.e. Z represents the noise latent variable, and I represents the identity matrix. By introducing a deep neural network as a general approximator to fit the true causal mechanism, Z can be expressed as a function of the observed data X through inverse function transformation, as shown in the framework below. Figure 1 As shown:

[0017] Z=f4((IW T (2)f3(X))

[0018] In this method, f1(·) and f2(·), f3(·) and f4(·) are set to be the same multilayer perceptron (MLP) to fit two square integrable functions in a high-dimensional space. This can achieve a high accuracy and greatly improve the algorithm's running speed.

[0019] For the following log-likelihood function:

[0020]

[0021] Here, n represents the sample size. Direct optimization requires calculating a high-dimensional integral over the noisy latent variable Z, which generally cannot yield an analytical solution unless very strict distributional assumptions are made. Therefore, this invention introduces a variational posterior distribution q(X|Z) to approximate the true posterior distribution p(X|Z), and uses the Kullback-Leibler divergence to constrain the distance between them, thus deriving the log-likelihood lower bound, i.e., the evidence lower bound (ELBO):

[0022]

[0023]

[0024] D KL This represents the KL divergence. Given a sample X k At that time, the distribution of the latent variable Z generated by the encoder follows q(Z|X). k The decoder attempts to base its algorithm on probability density p(X). k |Z) Reconstruct X through the latent variable Z k .

[0025] Essentially, the parameters to be determined are the weight matrix W representing the causal graph (i.e., W0). i,j ≠0 corresponds to the existence of an edge from i to j) and variational posterior distribution q θ The variational parameter θ of (Z|X) (represented as an unknown parameter in the variational distribution family) is used to construct a variational graph autoencoder. Since the weight matrix W corresponds to causal relationships, its corresponding graph needs to be restricted to a directed acyclic graph. Therefore, this invention introduces a differentiability constraint:

[0026]

[0027] Similar implementations used to control the acyclic property of a graph include:

[0028]

[0029] Here, tr(·) represents the trace of the matrix, ° represents the Hadamard product of the matrices, d is the size of the causal graph, I represents the identity matrix, and μ is a constant; the overall framework is as follows: Under the constraints, W and θ are obtained by maximizing the lower bound of evidence ELBO through an optimization algorithm. This belongs to the standard causal discovery method based on continuous optimization.

[0030] The specific process for constructing the conditional independence constraints described in step (II) is as follows:

[0031] Based on the aforementioned continuous optimization framework, this invention adds a constraint term based on CI testing. This requires establishing a connection between CI and the weight matrix W in continuous optimization, specifically: given a directed acyclic graph G, vertex x i With x j On G, all connections x are d-separated by set Z, i.e., all connections x are d-separated. i With x j All paths are blocked by the set Z. It can be proven that for any two points x in the vertex set... i ,x j Adjacent if and only if there is no set Z satisfying x i ⊥ G x j |Z. Among them, ⊥ G This indicates d-separation. Generally, d-separation ⊥G There is a one-to-one correspondence between CI and W. Since d-separation can be directly obtained from W, there is also a correspondence between CI and W.

[0032] Then, a causal skeleton can be obtained through CI testing. During subsequent continuous optimization, only the parameters in W covered by this skeleton need to be searched. That is, CI can be regarded as prior knowledge before continuous optimization. When optimizing the weight matrix using continuous optimization methods, the prior CI can be used to infer that the values ​​of the uncovered parts are 0, thereby compressing the parameter search space. This method can be formalized as follows: Define the masking matrix M:

[0033]

[0034] like Figure 2 As shown, a set of (conditional) independence relationships are obtained from the data matrix X using the CI test, which yields the values ​​of each element of the masking matrix M.

[0035] After finding M, a CI hard constraint is defined based on the continuous optimization framework:

[0036]

[0037] in, It is a directed acyclic constraint. The sparsity loss is used to control the complexity of the model, and α and β are hyperparameters that adjust the relative magnitude of the loss. τ is a threshold used to determine the existence of edges; if |W i,j |>τ, meaning there exists an edge x i →x j .

[0038] However, since the sample size cannot actually approach infinity, CI tests are not always accurate. Once a Type II error occurs (please explain this concept), i.e., the null hypothesis is incorrect but the null hypothesis is not rejected, it means that the hard constraint will inevitably incorrectly "cut off" the corresponding edge. This clearly contradicts the original intention of introducing CI tests in this invention. Therefore, this invention introduces a second approach, using the aforementioned hard constraint as a regularization term in the objective function, thereby transforming it into a soft constraint:

[0039]

[0040] in, Represents the Hadamard product of matrices. Let γ represent the Frobenius norm of the matrix, and γ be a hyperparameter used to balance the effects of regularization, which can be seen as the belief degree in the CI test. By solving the above objective function through a numerical optimization algorithm, the optimal weight matrix W can be obtained, and the final causal network graph can be derived.

[0041] The two methods described above that integrate CI testing into continuous optimization are hard constraint terms and soft regularization terms, respectively. The former is suitable for scenarios where CI testing rarely fails, such as sparse graph scenarios; while the latter is more robust and less sensitive to CI testing errors, but is relatively more time-consuming. Both application strategies, from an optimization perspective, are better than directly removing M. i,j Edges (i,j) with a value of 1 are strictly more robust.

[0042] Step (iii) of the equality constraint optimization specifically refers to:

[0043] Equations (9) and (10) both belong to equality optimization problems. The usual solution approach is the quadratic penalty method (QPM) or the augmented Lagrangian method (ALM). Since the properties of the augmented Lagrangian method cannot be guaranteed under the given problem conditions, this invention chooses the QPM algorithm for optimization.

[0044] Specifically, let the quadratic penalty function be denoted as... in This refers to the equality constraints in equations (9) and (10), i.e., the hard constraints. and soft constraints ρ is a penalty factor that increases continuously during the optimization process. Starting from the initial ρ0, the Adam algorithm is used to solve the sub-optimization problem under the current ρ, obtain the current minimum value, and determine whether the optimization convergence condition is met (usually by whether the norm of the gradient is greater than a very small threshold). If it is met, it can be directly used as the final optimization result; otherwise, ρ = β·ρ (β>1 is the scaling factor) is updated, and the current minimum value is used as the starting point for the next optimization sub-optimization problem. This process is repeated until convergence or the maximum number of iterations is reached.

[0045] The present invention has at least the following beneficial effects:

[0046] (1) The variational graph autoencoder models the causal generation process, introduces the strong fitting ability of the neural network, and ensures the feasibility of parameter learning through variational inference technology.

[0047] (2) Conditional independence constraints, which include constructing edge constraints using low-order conditional independence tests. Generally, faster conditional independence test methods based on partial correlation are used, which can ensure test quality with a small sample size while satisfying certain assumptions;

[0048] (3) Multiple combinations of hard and soft constraints. When conditional independence constraints are incorporated into continuous optimization methods, they are subject to the Type II error of conditional independence testing, which restricts edges that should exist from not existing. This invention designs hard constraints (requiring only that the constraints are valid at the end of optimization) and soft constraints (regularization terms) to mitigate the impact of Type II error in conditional independence testing;

[0049] (4) Optimize the entire framework using the quadratic penalty function method to ensure the convergence of the entire model under hard constraints.

[0050] Other advantages, objectives and features of the present invention will be apparent in part from the following description, and in part from the understanding of those skilled in the art through study and practice of the invention. Attached Figure Description

[0051] Figure 1 This is a framework diagram of the present invention.

[0052] Figure 2 A schematic diagram of the masking matrix is ​​shown.

[0053] Figure 3 The results of causal discovery on the virtual dataset are shown in comparison. Detailed Implementation

[0054] The present invention will now be described in further detail with reference to the accompanying drawings, so that those skilled in the art can implement it based on the description.

[0055] It should be understood that terms such as “having,” “comprising,” and “including” as used herein do not exclude the presence or addition of one or more other elements or combinations thereof.

[0056] I. Modeling of Causal Generation Processes

[0057] like Figure 1 As shown, this embodiment of the invention provides a method for modeling causal generation processes using a variational graph autoencoder. It is assumed that the structural causal model behind the observed data satisfies: X = f2((IW T ) -1 f1(Z)). Z can be expressed as a function of the observed data through inverse function transformation: Z = f4((IW) T f3(X)). Here, f1(·) and f2(·), f3(·) and f4(·) are parameter-sharing neural networks. The implementation is a single-layer MLP with ReLU activation.

[0058] Directly optimizing the following log-likelihood function requires calculating a high-dimensional integral over the noise latent variable Z, and analytical solutions are generally unavailable unless very strict assumptions are made.

[0059]

[0060] Therefore, this study introduces the variational posterior distribution q(Z|X) to approximate the true posterior distribution p(Z|X), and minimizing the KL ratio between them is equivalent to maximizing the log-likelihood lower bound, i.e., the evidence lower bound (ELBO).

[0061]

[0062]

[0063] Among them, D KL This represents the KL divergence. Given a sample X k At that time, the distribution of the latent variable Z generated by the encoder follows q(Z|X). k The decoder attempts to base its algorithm on probability density p(X). k |Z) Reconstruct X through the latent variable Z k The prior distribution is chosen as the matrix standard normal distribution. (I denotes the identity matrix). The variational posterior distribution q(Z|X) is selected from the factored Gaussian distribution family. Similar to the training of traditional variational autoencoders, this study directly models the mean M of q(Z|X). z S and covariance matrix z That is, q(Z|X)=N(M) z ,S z ).

[0064] [M z ,logS z ]=f4((IW T f3(X))

[0065] Based on this, a specific closed-form calculation formula can be derived from the KL divergence term in ELBO:

[0066]

[0067] The reconstruction terms can be obtained through Monte Carlo estimation:

[0068]

[0069] Here, c is a constant term, and Z is sampled from the variational posterior distribution q(Z|X). Furthermore, since Z is considered a latent variable, its dimension may differ from d. Sometimes, if the observed data has a small intrinsic dimension, a number smaller than d can be used as the dimension of Z.

[0070] II. Construction of Conditional Independence Constraints

[0071] Building upon the aforementioned continuous optimization framework, this study adds a constraint term based on the CI test. This requires establishing a connection between CI and the weight matrix W in continuous optimization, specifically: given a directed acyclic graph G, vertex x i With x j On G, all connections x are d-separated by set Z, i.e., all connections x are d-separated. i With x j All paths are blocked by the set Z. It can be proven that for any two points x in the vertex set... i ,x j Adjacent if and only if there is no set Z satisfying x i ⊥ G x j |Z. Among them, ⊥ G This indicates d-separation. Generally, d-separation ⊥ G There is a one-to-one correspondence between CI and W. Since d-separation can be directly obtained from W, there is also a correspondence between CI and W.

[0072] Then, a causal framework can be obtained through CI testing. Theoretically, when the CI test is completely correct, only the parameters in W covered by this framework need to be searched in subsequent continuous optimization; the values ​​of the uncovered parts can be inferred to be 0. Therefore, CI can be regarded as a kind of prior knowledge before continuous optimization. When optimizing the weight matrix using continuous optimization methods, the prior CI is used to determine in advance whether some elements in W are 0, thereby compressing the parameter search space.

[0073] like Figure 2 As shown, this idea can be formalized as follows: Define a masking matrix M

[0074]

[0075] The CI test is used to obtain a set of (conditional) independence relationships from the data matrix X, i.e., the values ​​of each element of the masking matrix M. This study uses a partial correlation test (CI). Its main idea is to calculate the partial correlation ρ using the following recursive calculation method. XY·Z :

[0076]

[0077] Where ρ represents the correlation coefficient. Condition set It degenerates into correlation. Apply Fisher's z-transform. Afterwards, if If the null hypothesis x⊥y|Z is true, then reject the null hypothesis; otherwise, accept the null hypothesis x⊥y|Z. Here, |Z| represents the size of set Z, and Φ... -1This represents the inverse of the cumulative distribution function of the standard normal distribution, and α is the significance level (typically taken as 0.05 or 0.01). Furthermore, we set the condition set size to a maximum of 1 to ensure its reliability with low sample sizes.

[0078] After obtaining M, a CI hard constraint is defined based on the continuous optimization framework:

[0079]

[0080]

[0081] in, It is a directed acyclic constraint. The sparsity loss is used to control the complexity of the model. α and β are hyperparameters that adjust the relative magnitude of the loss. τ is a threshold used to determine the existence of edges; if |W i,j |>τ, meaning there exists an edge x i →x j In this study, the threshold τ is generally set to 0.3. The value of α is selected from [0.25, 0.5, 1.0].

[0082] However, since the sample size cannot actually approach infinity, CI testing may produce errors. As a hypothesis test, a Type II error would mean that hard constraints would inevitably "cut off" the corresponding edges incorrectly. This clearly contradicts the original intention of introducing CI testing in this project. Therefore, this project introduces a second approach, treating the aforementioned hard constraints as regularization terms in the objective function, thus transforming them into soft constraints:

[0083]

[0084] Where ° represents the Hadamard product of matrices, Let represent the Frobenius norm of the matrix, and α and β be hyperparameters that adjust the relative magnitude of the loss. γ is a hyperparameter used to balance the effects of regularization and can be seen as the belief degree in the CI test. By solving the above objective function using a numerical optimization algorithm, the optimal weight matrix W can be obtained, leading to the final causal network graph.

[0085] The two methods described above for integrating CI testing into continuous optimization methods are hard constraint terms and soft regularization terms, respectively. The former is suitable for scenarios where CI testing rarely fails, such as sparse graph scenarios; while the latter is more robust and less sensitive to CI testing errors, but is relatively more time-consuming. Both application strategies, from an optimization perspective, are better than directly removing M. i,j Edges (i,j) with a value of 1 are strictly more robust.

[0086] III. Optimization based on equality constraints

[0087] Equations (9) and (10) both belong to equality optimization problems. The usual solution approach is the quadratic penalty method (QPM) or the augmented Lagrangian method (ALM). Since the properties of the augmented Lagrangian method cannot be guaranteed under the given problem conditions, this invention chooses the QPM method for optimization.

[0088] Specifically, let the quadratic penalty function be denoted as... in This refers to the equality constraints in equations (9) and (10), i.e., the hard constraints. and soft constraints ρ is a penalty factor that increases continuously during the optimization process. Starting from an initial ρ0, the Adam algorithm is used to solve the sub-optimization problem under the current ρ, obtaining the current minimum value. The convergence condition is then checked (usually by checking if the gradient norm is greater than a small threshold). If satisfied, it is directly used as the final optimization result; otherwise, ρ is updated to β·ρ (β>1 is a scaling factor), and the current minimum value is used as the starting point for the next sub-optimization problem. This process is repeated until convergence or the maximum number of iterations is reached. In this study, β is generally taken as 1.5.

[0089] This application also provides verification experiments to further demonstrate the technical effectiveness of this application.

[0090] To verify the performance of our method in causal discovery, we designed a virtual dataset and also selected a publicly available biological dataset, Sachs, for testing.

[0091] Virtual datasets: The causal backbone is composed of the backbones of three public datasets: Alarm, Hailfinder, and Win95pts. The Alarm network has 37 vertices, with an average degree of 2.49 and a maximum degree of 6. The Hailfinder network has 56 vertices, with an average degree of 2.36 and a maximum degree of 17. The Win95pts network has 76 vertices, with an average degree of 1.84 and a maximum degree of 9. The data is generated as follows: x = 2sin(W T (x+0.5·1))+W T (x+0.5·1)+z. Where the elements in W are uniformly sampled from (-2.0, -0.5)∪(0.5, 2.0). If x is the root node, the noise z is uniformly distributed. Mid-sampling, otherwise the noise z from a uniform distribution Sampling.

[0092] Sachs is a protein signaling network dataset that validates a real causal graph through biological experiments. The real causal graph has 11 nodes and 17 edges, and the dataset contains 853 independent and identically distributed samples.

[0093] To verify the superiority of this method, this example is compared with the following existing causal discovery methods on publicly available biological datasets: ICA-LiNGAM (excerpted from "Shimizu S, Hoyer PO, A, et al. A linear non-Gaussian acyclic model for causal discovery. in Journal of Machine Learning Research, 2006, 7(10).”), MMHC (extracted from “Ioannis Tsamardinos, Laura E. Brown, Constantin F. Aliferis: The max-min hill-climbing Bayesian network structure learning algorithm. in Machine Learning, 2006, 31–78.”), CAM (extracted from “Peter Bühlmann, Jonas Peters, Jan Ernest: CAM: Causal additive models, high-dimensional order search and penalized regression. in The Annals of Statistics, 2014, pp. 2526–2556.”), NOTEARS (extracted from “Xun Zheng, Bryon Aragam, Pradeep Ravikumar, Eric P. Xing: DAGs with NOTEARS: Continuous Optimization for Structure Learning. in Advances in Neural Information Processing Systems. 2019, pp. 8024–8035.”), GOLEM (extracted from “Ignavier Ng, Amir Emad Ghassami, Kun Zhang: On the Role of Sparsity and DAG Constraints for Learning Linear DAGs. in Advances in Neural Information Processing Systems, 2020, pp. 17943-17954.”), RL-BIC (extracted from “Zhu Shengyu, Ignavier Ng, Zhitang Chen:Causal discovery with reinforcement learning. in arXiv preprint arXiv:1906.04477(2019)."), DARING (from "Yue He, Peng Cui, Zheyan Shen, Renzhe Xu, Furui Liu, Yong Jiang: Daring: Differentiable causal discovery with residual independence." in Proceedings of the 27th ACM SIGKDD Conference on KnowledgeDiscovery&Data Mining.2021, pp.596-605."), Gran-DAG (from "Se′bastien Lachapelle, Philippe Brouillard, Tristan Deleu, Simon Lacoste-Julien: Gradient-based neuraldag learning.in arXiv preprint arXiv:1906.02226(2019)."), DAG-GNN (from "Yue Yu, Jie Chen, Tian Gao, Mo Yu: DAG-GNN: DAG Structure Learning with Graph NeuralNetworks.in International Conference on Machine Learning.PMLR, 2019:7154-7163. This embodiment uses structural Hamming distance (SHD), precision, recall, and F1-score as evaluation metrics to measure the performance of each algorithm. SHD compares a standard distance between two graphs using adjacency matrices. It involves calculating the difference between two (binary) adjacency matrices: each missing or absent edge in the target graph is counted as an error.

[0094] Experimental results are as follows Figure 3As shown in Table 1, the method of this invention uses DAG-GNN+CIR (soft constraint strategy) and DAG-GNN+CIC (hard constraint strategy). Parentheses following these indicate the dimension of the CI test condition set. On the virtual dataset, it can be seen that the method of this invention has a smaller SHD in all cases. In particular, significant performance improvements are achieved when the sample size is 200. Compared with DAG-GNN, the SHD is reduced by 62.34% / 60.17%, 51.80% / 60.02%, and 72.66% / 71.15% on Alarm, Hailfinder, and Win95pts, respectively. This fully demonstrates the effectiveness of the hard and soft constraint strategies in this invention, especially when the sample size is small. The improvement is most significant for Win95pts because Win95pts has the sparsest structure (lowest average degree). DAG-GNN trains deep neural networks to capture the complex nonlinear structure of the data, but it is prone to overfitting with small sample sizes. Furthermore, it is noted that on the Hailfinder skeleton, the soft-constraint policy has a smaller SHD than the hard-constraint policy. On the other two skeletons, their performance is similar. The main reason why the soft-constraint policy outperforms the hard-constraint policy on the Hailfinder is that our chosen CI testing method has a higher Type II error rate on the Hailfinder, and the soft-constraint policy is more robust to Type II errors in CI than the hard-constraint policy. On the real dataset, both constraint policies achieve the minimum SHD of 11, a result also achieved by RL-BIC, DARING, and NOTEARS-MLP.

[0095] Although embodiments of the present invention have been disclosed above, they are not limited to the applications listed in the specification and embodiments. They can be applied to various fields suitable for the present invention. For those skilled in the art, other modifications can be easily made. Therefore, without departing from the general concept defined by the claims and their equivalents, the present invention is not limited to the specific details and illustrations shown and described herein.

[0096] Table 1. Comparison of causal discovery results on the Sachs biological dataset.

[0097]

Claims

1. A causal discovery method based on conditional independence testing and continuous optimization, characterized in that, The specific steps are as follows: (I) Modeling the causal generation process; Modeling the causal generation process using a variational graph autoencoder model; In this process, variational inference techniques are introduced, and variational parameters and true causal graph parameters are learned by minimizing the KL divergence between the variational distribution and the posterior distribution of the true latent variables; In order to ensure that the optimized weight matrix can derive an acyclic causal graph, a differentiable constraint is introduced to restrict the search space to a directed acyclic graph space. (ii) Construction of conditional independence constraints; including preprocessing the image to obtain smaller image patches, and using a pre-trained image feature extraction model to extract image patch features, and finally outputting a set of image patch representations. (III) Equality-constrained optimization: The objective function of the variational graph autoencoder is combined with the conditional independence constraint and expressed as an equality-constrained optimization problem. The parameters are optimized using the quadratic penalty function method to obtain the final causal graph. The causal generation process modeling described in step (I) is as follows: Suppose that the structural causal model behind the observation data satisfies the following form: X=f2((I-W T ) -1 f1(Z)), (1) Where X is the observed variable, W is the causal weight matrix, and W i,j ≠0 corresponds to an edge i→j in the causal graph G, i.e. Z represents the noise latent variable, and I represents the identity matrix. By introducing a deep neural network as a general approximator to fit the true causal mechanism, Z is expressed as a function of the observed data X through inverse function transformation: Z=f4((I-W T )f3(X)) , (2) In this case, f1(·) and f2(·), f3(·) and f4(·) are set to be completely identical multilayer perceptrons (MLPs) to fit two square integrable functions in a high-dimensional space. For the following log-likelihood function: Where n represents the sample size; the variational posterior distribution q(X|Z) is introduced to approximate the true posterior distribution p(X|Z), and the distance between the two is constrained by the KL divergence, thus deriving the lower bound of log-likelihood, i.e., the lower bound of evidence ELBO: D KL Denotes the KL divergence; given a sample X k At that time, the distribution of the latent variable Z generated by the encoder follows q(Z|X). k The decoder is based on the probability density p(X). k |Z) Reconstruct X through the latent variable Z k ; The parameters to be determined are the weight matrix W representing the causal graph and the variational posterior distribution q. θ The variational parameters θ of (Z|X) are combined to form a variational graph autoencoder. Since the weight matrix W corresponds to causal relationships, its corresponding graph needs to be restricted to a directed acyclic graph. Therefore, a differentiability constraint is introduced: Similar implementations used to control the acyclic property of a graph include: Here, tr(·) represents the trace of the matrix. This represents the Hadamard product of matrices, where d is the size of the causal graph, I represents the identity matrix, and μ is a constant; the overall framework is... Under the constraint of maximizing the lower bound of evidence ELBO through an optimization algorithm, W and θ are obtained. The specific process for constructing the conditional independence constraints described in step (II) is as follows: Based on the aforementioned continuous optimization framework, an additional constraint term based on the CI test is added; to this end, CI is linked to the weight matrix W in continuous optimization: given a directed acyclic graph G, vertex x i With x j On G, all connections x are d-separated by set Z, i.e., all connections x are d-separated. i With x j All paths are blocked by set Z; here, any two points x in the vertex set are blocked by set Z; i x j Adjacent if and only if there is no set Z satisfying x i ⊥ G x j |Z; where ⊥ G This represents d-separation; d-separation can be directly obtained from W, therefore there is a correspondence between CI and W; Then, a causal skeleton is obtained through CI testing. During subsequent continuous optimization, only the parameters in W covered by this skeleton need to be searched; that is, CI is regarded as prior knowledge before continuous optimization. When optimizing the weight matrix using the continuous optimization method, the value of the uncovered part is inferred to be 0 with the help of the prior CI, so as to compress the parameter search space; its formalization is as follows: Define the masking matrix M: Using CI testing, a set of conditional independence relationships are obtained from the data matrix X, which is equivalent to obtaining the values ​​of each element of the masking matrix M. After finding M, a CI hard constraint is defined based on the continuous optimization framework: in, It is a directed acyclic constraint. The sparsity loss is used to control the complexity of the model. α and β are hyperparameters that adjust the relative magnitude of the loss. τ is a threshold used to determine the existence of edges. If |W i，j |>τ, meaning there exists an edge x i →x j ; Furthermore, the aforementioned hard constraint is used as a regularization term in the objective function, thereby transforming it into a soft constraint: in, Represents the Hadamard product of matrices. The Frobenius norm of the matrix is ​​represented by γ, which is a hyperparameter used to balance the effects of regularization and can be regarded as the belief degree of CI test. By solving the above objective function through numerical optimization algorithm, the optimal weight matrix W can be obtained, and the final causal network graph can be derived. The equality constraint optimization described in step (iii) is to optimize equations (9) and (10) of the equality constraint problem, specifically using the quadratic penalty function method for optimization; Denote the quadratic penalty function in This refers to the equality constraints in equations (9) and (10), i.e., the hard constraints. and soft constraints ρ is a penalty factor that increases continuously during the optimization process. Starting from the initial ρ0, the Adam algorithm is used to solve the sub-optimization problem under the current ρ, obtain the current minimum value, and determine whether the optimization convergence condition is met. If it is met, it is directly used as the final optimization result; otherwise, ρ is updated to β·ρ, where β>1 is the scaling factor. The current minimum value is used as the starting point of the next optimization sub-optimization problem. This process is repeated until convergence or the maximum number of iterations is reached.

Images