A symbolic regression method for network dynamics based on LLM and Bayesian optimization

Abstract

The present invention is applicable to the field of computer technology and provides a network dynamics symbolic regression method based on LLM and Bayesian optimization. This method does not require a predefined basis function library, but instead uses natural language to integrate physical priors to guide the LLM to generate equations, further improving the quality of the generated equations and accelerating the search process. At the same time, it utilizes Bayesian optimization, a serialized optimization method, to efficiently search for global optimal solutions, improving the accuracy of the generated equations and providing an efficient and intelligent solution for symbolic model search. Experiments in various network dynamics scenarios have demonstrated that, without the need for an additional predefined basis function library, the method effectively incorporates prior knowledge by generating candidate equations through LLM and combining it with the efficient search of Bayesian optimization, significantly improving the accuracy of the discovered equations and the efficiency of the search. The resulting equations are far more accurate than current methods based on GNN combined with GP (GG) and sparse regression with a predefined function library (TPSINDy).

Description

Technical Field

[0001] The present invention belongs to the field of computer technology, and in particular relates to a network dynamics symbolic regression method based on LLM and Bayesian optimization. Background Art

[0002] Network dynamics are common in natural and social phenomena, such as electrical signal transmission in neural networks, epidemic spread, and interactions between ecological populations. The behavior of these systems is driven by nonlinear interactions between nodes, and their dynamics can be modeled as differential equations on the network. For example, in scenarios such as epidemic spread, neuronal activity, and interactions between ecological populations, system dynamics are often determined by interactions between node states and topology. Symbolic regression aims to be data-driven, extracting the underlying equations and mathematical laws from observational data to accurately derive dynamic models of complex systems and help uncover the governing laws. For example, establishing an accurate dynamic model of epidemic spread can predict the spread of an epidemic and facilitate its prevention and control. However, in real-world scenarios, the governing equations of the system are often unknown, and high-dimensional, missing, and noisy observational data further complicates modeling. Traditional methods rely on expert experience or simplifying assumptions, making it difficult to accurately capture the underlying laws of complex dynamics.

[0003] Currently, there are two main approaches for symbolic regression of network dynamics. The first utilizes graph neural networks (GNNs) to model spatial dependencies in network dynamics. GNNs model spatial dependencies between nodes through a message-passing mechanism, effectively extracting network topological features. Genetic programming (GP) is then used to solve symbolic equations within the neural network. The second approach involves sparse regression methods based on predefined function libraries, such as TPSINDy. The core idea is to model the dynamic equations as linear combinations of predefined basis functions (such as polynomials, trigonometric functions, and exponential functions), and then filter key terms through sparse optimization. For example, TPSINDy improves efficiency through a two-stage regression approach (first filtering function terms, then optimizing coefficients).

[0004] The Large Language Model (LLM), with its powerful natural language understanding and text generation capabilities, can generate a diverse set of candidate equations to address various network dynamics scenarios. Furthermore, the LLM can leverage prior knowledge (such as physical laws) to guide the generation process, further improving the quality and interpretability of the generated equations.

[0005] However, existing methods still have many flaws: Although the method based on GNN combined with GP can obtain the topological characteristics of the network through the message passing mechanism, its black box nature limits the interpretability of the model. There is a semantic gap between the nonlinear mapping of the neural network and the symbolic expression, which may cause the analytical results to lose some physical meaning. At the same time, the search efficiency of GP is low, especially when the equation complexity is high, the calculation time may exceed the actual requirements. The sparse regression method based on the predefined function library selects key terms through the predefined basis function library and sparse regression, but its basis function library needs to be manually designed, and its performance is highly dependent on the completeness of the function library. If the library does not cover the key terms of the real dynamics, such as nonlinear coupling terms, it cannot be accurately modeled. The symbolic regression technology based on LLM requires the generation of a large number of candidate equations, which leads to a sharp expansion of the search space. There is a lack of effective methods to guide LLM to generate more accurate candidate equations.

[0006] To this end, the present invention proposes a network dynamics symbolic regression method based on LLM and Bayesian optimization, aiming to solve the above problems and achieve accurate and efficient mining of interpretable dynamic equations from data without relying on strong prior knowledge. Summary of the Invention

[0007] The purpose of the present invention is to provide a network dynamics symbolic regression method based on LLM and Bayesian optimization, aiming to solve the problems raised in the above background technology.

[0008] The purpose of the present invention is achieved through the following technical solutions:

[0009] A symbolic regression method for network dynamics based on LLM and Bayesian optimization, comprising the following steps:

[0010] Step 1: LLM generates candidate equations;

[0011] Collect network dynamics domain knowledge as LLM prompt words, call LLM to generate candidate equations until the number of generated equations reaches a preset value, form the initial search population and add it to the candidate equation set;

[0012] Step 2: Bayesian optimization searches for the optimal solution, including:

[0013] Sub-step 2.1: Build and train a proxy model to convert the equation text description into a vector using an embedding model and use a neural network to predict the fitness value;

[0014] Sub-step 2.2: Design the acquisition function and select the most promising equation;

[0015] Sub-step 2.3: Iteratively update the proxy model until the optimal solution is found or the preset number of iterations is reached;

[0016] Step 3: Evaluate the candidate equations and calculate the fitness value;

[0017] Preprocess the network dynamics observation data into structured data and extract node states and neighbor relationships;

[0018] For the equation combination to be evaluated, the constant term is optimized by the quasi-Newton method, and the normalized mean square error is used as the fitness value;

[0019] Step 4: Determine whether there is a true equation;

[0020] If it exists, its equation form is output and the process ends; if it does not exist, the most promising equation that has been evaluated is added as knowledge to the prompt word of LLM, and the process returns to step 1 to continue generating candidate equations.

[0021] Furthermore, in step 1, the network dynamics domain knowledge includes background information of network dynamics, operation symbols in equations, and the number of variables.

[0022] Furthermore, in the sub-step 2.1, the proxy model is an integrated regression model based on a feedforward neural network, using a snapshot integration technique, specifically including:

[0023] During the training process, model snapshots are saved and the prediction results of the snapshots are averaged to obtain different predictions of the proxy model for the same input; the formula is as follows:

[0024] ;

[0025] in: is the embedding vector of the equation, Represents the proxy model embedding vector of the equation The predicted fitness value of Indicates the The fitness prediction value of each base model; is the number of base models; j is the index value of the base model;

[0026] The formulas for the predicted mean and variance of the surrogate model are as follows:

[0027] ;

[0028] ;

[0029] in: is the embedding vector of the equation, Represents the equation embedding vector obtained by the proxy model The predicted mean of Represents the equation embedding vector obtained by the proxy model The prediction variance of Represents the equation embedding vector obtained by the proxy model The predicted standard deviation of .

[0030] Furthermore, in the sub-step 2.2, the acquisition function adopts the expected improvement, and its calculation formula is:

[0031] ;

[0032] in: is the embedding vector of the equation, represents the value of the acquisition function; Represents the best fitness value of the known dynamic equation; represents a positive number; represents the cumulative distribution function of the standard normal distribution; and are the equation embedding vectors obtained by the proxy model The predicted mean and predicted standard deviation of ; represents the probability density function of the standard normal distribution; Indicates the complexity of the equation;

[0033] By calculation The maximum value of to obtain the embedding vector of the most potential equation to be evaluated next time ,Right now:

[0034] .

[0035] Furthermore, in the sub-step 2.3, the specific process of iteratively updating the proxy model includes:

[0036] The newly evaluated most promising equation and its fitness value are added to the evaluated dataset to form an updated training sample set; the proxy model is retrained based on the updated dataset, and the fitness value of the equation is predicted using the trained proxy model; the optimal equation in the next iteration cycle is screened by maximizing the acquisition function, and the iteration is repeated until the termination condition is met.

[0037] Furthermore, in step 3, the combination of equations to be evaluated includes the equations describing the node's own dynamics. equations and describes the interactions between nodes The equation is calculated by the following formula:

[0038] ;

[0039] in: for and The fitness value of the node itself is dynamic Equations and interactions between nodes The loss value between the equation and the true equation; is the normalized mean square error; express and By minimizing the constant term in get; It is a finite difference calculation; is the preprocessed network dynamics observation data; is a variable The set of neighbors of Operations to obtain neighbor information in observation data; The neighbor node index in the network topology; It is the status of neighbor nodes in the network topology.

[0040] Compared with the prior art, the present invention has the following beneficial effects:

[0041] This method, without requiring a predefined basis function library, can integrate physical prior knowledge into LLM prompts through a concise natural language description. This guides the LLM in generating candidate equations, further improving the quality of the generated equations and accelerating the search process. Furthermore, it leverages Bayesian optimization, a serialized optimization method, to efficiently search for global optimal solutions through an iterative mechanism of surrogate models and acquisition functions, thereby improving the accuracy of the generated equations and providing an efficient and intelligent solution for symbolic model search. Experiments conducted in various network dynamics scenarios demonstrate that, without requiring an additional predefined basis function library, this method effectively incorporates prior knowledge by combining LLM generation of candidate equations with the efficient search of Bayesian optimization. This significantly improves the accuracy of the discovered equations while also enhancing search efficiency. The resulting equations achieve an accuracy far superior to existing methods based on GNN combined with GP (GG) and sparse regression with a predefined function library (TPSINDy). BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 It is a framework diagram of the method of the present invention.

[0043] Figure 2 Flowchart of the method of the present invention.

[0044] Figure 3 Generate a flow chart of candidate equations for LLM.

[0045] Figure 4 Flowchart for Bayesian optimization. DETAILED DESCRIPTION

[0046] In order to have a clearer understanding of the technical features, objectives and beneficial effects of the present invention, the technical solution of the present invention is now described in detail below, but it should not be understood as limiting the scope of implementation of the present invention.

[0047] The present invention provides a network dynamics symbolic regression method based on LLM and Bayesian optimization, the framework diagram and flow chart of which are as follows: Figure 1 and Figure 2 As shown, the method includes the following steps:

[0048] Step 1: LLM generates candidate equations;

[0049] The process of LLM generating candidate equations is as follows Figure 3 As shown in the figure, relevant knowledge about network dynamics is first collected, including background information about network dynamics, operational symbols in equations, and the number of variables. This information is used as prompts for the LLM, which is then invoked to generate candidate equations. When the number of generated equations reaches a preset value, n, these equations form the initial search population and are added to the set of candidate equations. If this number is not reached, the LLM is called again until the required number is met.

[0050] Step 2: Bayesian optimization (BO) searches for the optimal solution;

[0051] Bayesian optimization iteratively searches for the optimal equation through the surrogate model and acquisition function. The process is as follows Figure 4 As shown. First, the evaluated candidate equations and their fitness values ​​are mapped into vectors through the embedding model. The mapped vectors are used to train the proxy model. The trained proxy model is used to predict the fitness value of the equation and output the predicted mean and predicted variance. In the first iteration, all candidate equations generated by the LLM are evaluated to initialize the training proxy model. Based on the acquisition function that maximizes the predicted mean and variance, the optimal equation is screened and evaluated. Finally, the selected equation is added to the set of evaluated candidate equations for subsequent training of the proxy model. The specific steps are as follows:

[0052] Sub-step 2.1: Proxy model construction and training;

[0053] The surrogate model in Bayesian optimization is mainly used to approximate the objective function to be optimized. Its design quality directly affects the optimization effect. This paper uses an integrated regression model based on a feedforward neural network as a surrogate model to predict the fitness value of the equation. The specific process is as follows:

[0054] The equation text description is converted into a vector through an embedding model (such as NV-Embed-v1), that is:

[0055] ;

[0056] in: Represents the predicted fitness value; NN is a neural network, output The predicted fitness value of Represents the text description of the equation; It is the embedding vector obtained by the embedding model of the equation text description.

[0057] Using a multi-layer perceptron (MLP) as the base model, ensemble learning combines the predictions of several base models, NN1, NN2, etc., to improve model performance and robustness. To improve the computational efficiency of surrogate model training and inference, snapshot ensembles are employed. Specifically, while training the neural network, M model snapshots are saved, each corresponding to a different state of the surrogate model. The predictions from these snapshots are then averaged. This snapshot ensemble yields different predictions from the surrogate model for the same input, providing an estimate of the uncertainty of the model's predictions. The formula is as follows:

[0058] ;

[0059] in: is the embedding vector of the equation, Represents the proxy model embedding vector of the equation The predicted fitness value of Indicates the The fitness prediction value of each base model; is the number of base models; j is the index value of the base model.

[0060] The formulas for the predicted mean and variance of the surrogate model are as follows:

[0061] ;

[0062] ;

[0063] in: is the embedding vector of the equation, Represents the equation embedding vector obtained by the proxy model The predicted mean of Represents the equation embedding vector obtained by the proxy model The prediction variance of Represents the equation embedding vector obtained by the proxy model The predicted standard deviation of .

[0064] When training the proxy model, the model parameters are optimized by gradient descent. ,Right now:

[0065] ;

[0066] in: represents the optimal parameters of the model; i Indicates the i group dynamics equations; nrepresents the total number of evaluated equations; Indicates the model parameters Next i The model predicts fitness values ​​for the group dynamics equations; Indicates the The true fitness value of the group dynamics equation.

[0067] Sub-step 2.2: Acquisition function design and optimal equation selection;

[0068] The acquisition function guides the search process and selects the next most promising candidate solution. The goal of the acquisition function is to balance exploration and exploitation, that is, to find a balance between areas known to perform well and areas unknown but potentially even better. Common acquisition functions in BO include probability boost, expectation boost, and confidence bound. This paper uses expectation hints to construct the acquisition function, and its basic formula is:

[0069] ;

[0070] in: represents the value of the acquisition function; It is the expectation calculation of the forecast distribution; Represents the best fitness value of the known dynamic equation; Represents the fitness value of the current equation. Based on the predicted mean and variance of the surrogate model, it is approximately calculated as:

[0071] ;

[0072] in: is the embedding vector of the equation, represents the value of the acquisition function; Represents the best fitness value of the known dynamic equation; Represents a small positive number used to fine-tune the acquisition function results; represents the cumulative distribution function of the standard normal distribution; and are the predicted mean and standard deviation obtained by the surrogate model, respectively; represents the probability density function of the standard normal distribution; Indicates the complexity of the equation;

[0073] By calculation The maximum value of to obtain the embedding vector of the most potential equation to be evaluated next time ,Right now:

[0074] .

[0075] Sub-step 2.3: Iteratively update the proxy model;

[0076] Evaluate the most promising equations screened by the acquisition function and then enter the iterative optimization phase. The specific process is as follows:

[0077] The newly evaluated most promising equation and its true fitness value are added to the evaluated dataset to form an updated training sample set; the proxy model is retrained based on the updated dataset, and the fitness of the equation is predicted using the trained proxy model. The most promising equation in the next iteration cycle is screened by maximizing the acquisition function, and the above "evaluation → update → training → prediction → screening" process is repeated until the optimal solution is found (the fitness value is less than the set target threshold) or the preset number of iterations is reached, and the optimization process is terminated.

[0078] Step 3: Evaluate the candidate equations and calculate the fitness value;

[0079] Data preprocessing step: transform network dynamics observation data Processed into structured data, including node status information and node neighbor relationships. Operation extracts each variable Neighbor set And the corresponding neighbor status , providing an effective data basis for subsequent equation evaluation.

[0080] Constant optimization and fitness calculation: Each combination of equations to be evaluated is composed of a constant describing the node's own dynamics. Equations (including operators) that describe the interactions between nodes The constant term in the equation is optimized by the quasi-Newton method. , the goal is to minimize the normalized mean square error (NMSE) between the equation's predicted value and the true value, and directly use NMSE as the fitness value (the smaller the value, the better the fitting effect). The specific calculation formula is as follows:

[0081] ;

[0082] in: for and The fitness value of the node itself is dynamic Equations and interactions between nodes The loss value between the equation and the true equation; express and The constant term in can be minimized by an optimization algorithm such as the quasi-Newton method get; It is a finite difference calculation; is the preprocessed network dynamics observation data; is a variable The set of neighbors of Operations to obtain neighbor information in observation data; The neighbor node index in the network topology; is the status of neighbor nodes in the network topology; is the normalized mean square error, and the formula is as follows:

[0083] ;

[0084] in: Indicates the number of data to be evaluated, Indicates the The value of the real data, Represents the predicted value of the current equation.

[0085] Step 4: Determine whether there is a true equation;

[0086] Check whether there is a true equation (i.e., an equation that truly describes the network dynamics) among the current candidate equations. If so, output its equation form and end the process. If not, add the most promising equation evaluated as relevant knowledge to the LLM prompt word and return to step 1 to continue generating candidate equations.

[0087] The specific implementation of the present invention is described in detail below with reference to specific embodiments.

[0088] Example 1: This example studies five representative uniform network dynamics, covering fields such as biology, ecology, neuroscience, and epidemiology, to verify the effect of the present invention.

[0089] 1. Biochemical kinetics (Bio): Biochemical processes within living cells are mediated by protein-protein interactions, where proteins bind to form protein complexes. The kinetics of these biochemical reactions can be expressed as:

[0090] ;

[0091] in: Indicates time, Indicates time Protein The concentration over time The derivative of and represent the average influx rate and degradation rate of proteins, respectively; Indicates time Protein concentration; Indicates time Protein The concentration value of Indicates protein and The effective rate constant of the interaction between N Indicates the number of proteins in the network topology. , .

[0092] 2. Gene Regulation Dynamics (Gene): The dynamics of the gene regulatory network can be described by the Michaelis-Menten equation, which is in the form:

[0093] ;

[0094] in: Indicates time, Indicates gene The expression level of The derivative of ; parameter Controlling the decay rate; Indicates gene expression levels; f represents a constant, when When , the first term on the right describes the degradation process; when The first term on the right describes dimerization, i.e., two or more identical molecules combine to form a dimer. The second term on the right represents gene activation, where is the Hill coefficient, which quantifies the saturation rate of the influence of neighboring nodes; Represents different genes in the network topology and The effective rate constant of the interaction between N represents the number of genes in the network topology; Indicates gene This embodiment sets , , .

[0095] 3. Lotka-Volterra Model (LV): The Lotka-Volterra Model (LV) describes the population dynamics of competing species:

[0096] ;

[0097] in: Indicates time, Indicates species The time derivative of the population size, Indicates species The population size, Indicates species Growth parameters, In the experiment and Sample from a uniform distribution in the range [0.5, 1.5]. In this example, and Both are set to 1. Represents different populations in the network topology and The interaction constant between them; N Represents the number of populations in the network topology; Represents population scale.

[0098] 4. Neural dynamics (Neur): The firing rate of a neuron can be described by the Wilson-Cowan model as:

[0099] ;

[0100] In this model, Indicates time, Represents neurons The derivative of the activity level over time; Represents neurons Activity level; parameters and They determine the slope and threshold of the neural activation function respectively; Represents neurons and The effective rate constant of the interaction between N Represents the number of neurons in the network topology; Represents neurons This embodiment sets , .

[0101] 5. Epidemic Dynamics (Epi): Epidemic dynamics can be used to describe outbreaks of infectious diseases:

[0102] ;

[0103] In this model, Indicates time, Represents the derivative of the individual infection probability over time. Each node can represent an individual, and the node state Corresponding node The infection probability of represents the rate at which individuals recover from infection, and its value is set to 1.0; Indicates a node and The effective rate constant of the interaction between NRepresents the number of individuals in the network topology; Representation node probability of infection.

[0104] To demonstrate the advantages of symbolic regression in various network dynamics scenarios, we conducted comparative experiments on mean squared error (MSE) and mean absolute percentage error (MAPE) for the five dynamic scenarios described above, comparing the proposed method with a GNN-based sparse regression method combined with GP (GG) and a predefined library (TPSINDy). All experiments were conducted on an Intel(R) Core(TM) i9-14900HX CPU and two NVIDIA GeForce RTX 4090 server GPUs, running Ubuntu 18.04.6. The LLM used in this paper was a locally deployed LLM with qwen2.5-7B parameters.

[0105] The specific results of MSE comparison are shown in Table 1:

[0106] Table 1 MSE results of network dynamics in different scenarios

[0107]

[0108] The specific results of the MAPE comparison are shown in Table 2:

[0109] Table 2 MAPE results of network dynamics under different scenarios

[0110]

[0111] Experimental data show that the mean square error and mean absolute percentage error of the present invention are much lower than those of the method based on GNN combined with GP (GG) and the sparse regression method with a predefined function library (TPSINDy), verifying its advantages in efficiency and accuracy in symbolic regression of network dynamics, and providing a better solution for dynamic modeling of complex systems.

[0112] The above are only preferred embodiments of the present invention. It should be pointed out that for those skilled in the art, several variations and improvements can be made without departing from the concept of the present invention. These should also be regarded as the scope of protection of the present invention. These will not affect the effect of the implementation of the present invention and the practicality of the patent.

Claims

1. A network dynamics symbolic regression method based on LLM and Bayesian optimization, characterized by: The following steps are involved: Step 1: LLM generates candidate equations; Collect network dynamics domain knowledge as LLM prompt words, call LLM to generate candidate equations until the number of generated equations reaches a preset value, form the initial search population and add it to the candidate equation set; Step 2: Bayesian optimization searches for the optimal solution, including: Sub-step 2.1: Build and train a proxy model to convert the equation text description into a vector using an embedding model and use a neural network to predict the fitness value; Sub-step 2.2: Design the acquisition function and select the most promising equation; Sub-step 2.3: Iteratively update the proxy model until the optimal solution is found or the preset number of iterations is reached; In the sub-step 2.1, the proxy model is an integrated regression model based on a feedforward neural network, using a snapshot integration technique, specifically including: During the training process, model snapshots are saved and the prediction results of the snapshots are averaged to obtain different predictions of the proxy model for the same input; the formula is as follows: ； in: is the embedding vector of the equation, Represents the proxy model embedding vector of the equation The fitness prediction value of Indicates the The fitness prediction value of each base model; is the number of base models; is the index value of the base model; The formulas for the predicted mean and variance of the surrogate model are as follows: ； ； in: is the embedding vector of the equation, Represents the equation embedding vector obtained by the proxy model The predicted mean of Represents the equation embedding vector obtained by the proxy model The prediction variance of Represents the equation embedding vector obtained by the proxy model The predicted standard deviation of Step 3: Evaluate the candidate equations and calculate the fitness value; Preprocess the network dynamics observation data into structured data and extract node states and neighbor relationships; For the equation combination to be evaluated, the constant term is optimized by the quasi-Newton method, and the normalized mean square error is used as the fitness value; Step 4: Determine whether there is a true equation; If it exists, its equation form is output and the process ends; if it does not exist, the most promising equation that has been evaluated is added as knowledge to the prompt word of LLM, and the process returns to step 1 to continue generating candidate equations.

2. The network dynamics symbolic regression method based on LLM and Bayesian optimization according to claim 1 is characterized in that: In step 1, the network dynamics domain knowledge includes background information of network dynamics, operation symbols in equations, and the number of variables.

3. The network dynamics symbolic regression method based on LLM and Bayesian optimization according to claim 1, characterized in that: In the sub-step 2.2, the acquisition function adopts the expected improvement, and its calculation formula is: ； in: is the embedding vector of the equation, represents the value of the acquisition function; Represents the best fitness value of the known dynamic equation; represents a positive number; represents the cumulative distribution function of the standard normal distribution; and are the equation embedding vectors obtained by the proxy model The predicted mean and predicted standard deviation of represents the probability density function of the standard normal distribution; Indicates the complexity of the equation; By calculation The maximum value of to obtain the embedding vector of the most potential equation to be evaluated next time ,Right now: 。 4. The network dynamics symbolic regression method based on LLM and Bayesian optimization according to claim 3 is characterized in that: In the sub-step 2.3, the specific process of iteratively updating the proxy model includes: The newly evaluated most promising equation and its fitness value are added to the evaluated dataset to form an updated training sample set; the proxy model is retrained based on the updated dataset, and the fitness value of the equation is predicted using the trained proxy model; the most promising equation in the next iteration cycle is screened by maximizing the acquisition function, and the iteration is repeated until the termination condition is met.

5. The network dynamics symbolic regression method based on LLM and Bayesian optimization according to claim 1, characterized in that: In step 3, the equation combination to be evaluated includes the equations describing the node's own dynamics. equations and describes the interactions between nodes The equation is calculated by the following formula: ； in: for and The fitness value of the node itself is dynamic Equations and interactions between nodes The loss value between the equation and the true equation; is the normalized mean square error; express and By minimizing the constant term in get; It is a finite difference calculation; is the preprocessed network dynamics observation data; is a variable The set of neighbors of Operations to obtain neighbor information in observation data; The neighbor node index in the network topology; It is the status of neighbor nodes in the network topology.

Images