Paper deep dive

A Multi-Agent LLM System for Automated Tuning of Specialized Neural Network Models

Ivan Moloshnikov, Dmitry Kunitsyn, Alexander Sboev, Maksim Skorokhodov, Roman Rybka, Yuanchao Liu

Year: 2026Venue: arXiv preprintArea: cs.AIType: PreprintEmbeddings: 154

Intelligence

Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%

Last extracted: 10/13/2026, 2:53:46 AM

Summary

This paper introduces a multi-agent Large Language Model (LLM) framework for the automated tuning of specific functional blocks within specialized neural network models, rather than optimizing entire architectures. The system employs three distinct agent roles—an Analyst, a Coder, and a Recorder—operating within a sandboxed environment with an experiment log. The framework was validated across five tasks: spiking neural network (SNN) classification on MNIST and Digits datasets, and solving high-order differential equations (GNLSE and KS) using Physics-Informed Neural Networks (PINN). Results demonstrate significant improvements over baselines, including increased classification accuracy and reduced error metrics for physical simulations.

Entities (11)

MNIST · dataset · 99%Multi-agent LLM framework · method · 98%digits · dataset · 95%LLM Analyst · agentrole · 95%LLM Coder · agentrole · 95%LLM Recorder · agentrole · 95%Physics-Informed Neural Network · modeltype · 92%Spiking Neural Network · modeltype · 92%GNLSE · equation · 90%KS Equation · equation · 90%Correlation Classwise Network · modelarchitecture · 88%

Relation Signals (10)

Multi-agent LLM framework → evaluatedon → digits

confidence 95% · spiking classification on MNIST and Digits

Multi-agent LLM framework → evaluatedon → MNIST

confidence 95% · spiking classification on MNIST and Digits

Multi-agent LLM framework → usesagentrole → LLM Recorder

confidence 95% · The multi-agent large language model (LLM) framework has three roles---analyst, coder, and recorder

Multi-agent LLM framework → usesagentrole → LLM Analyst

confidence 95% · The multi-agent large language model (LLM) framework has three roles---analyst, coder, and recorder

Multi-agent LLM framework → usesagentrole → LLM Coder

confidence 95% · The multi-agent large language model (LLM) framework has three roles---analyst, coder, and recorder

Multi-agent LLM framework → appliedto → Physics-Informed Neural Network

confidence 90% · solving high-order differential equations with Physics-Informed Neural Networks (PINN)

Multi-agent LLM framework → appliedto → Spiking Neural Network

confidence 90% · tuning a classification spiking neural network with correlation learning

Spiking Neural Network → usesarchitecture → Correlation Classwise Network

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

Abstract:Many machine learning approaches are two-level: an expert selects a class of methods, within this class the choice of key components is critical and usually manual. This paper addresses the automatic tuning of these components and proposes a framework transferable across models of different nature. The multi-agent large language model (LLM) framework has three roles---analyst, coder, and recorder~---with an experiment log and a sandboxed execution environment. Unlike systems that optimize the entire ML script or search architectures in a discrete space, the framework operates over an open but semantically narrow space: the agent generates code only for the specified components, leaving the overall network structure unchanged. The approach is validated in two directions: tuning a classification spiking neural network with correlation learning, and solving high-order differential equations with Physics-Informed Neural Networks (PINN). Across five tasks---spiking classification on MNIST and Digits, and boundary-value problems for the generalized nonlinear Schrödinger equation (GNLSE) and Kuramoto--Sivashinsky (KS) equations---the agent outperforms the baselines: classification accuracy is 13 percentage points (MNIST) and 6 percentage points (Digits) higher; the error against the analytical single-soliton solution is reduced approximately 2491-fold; the invariant-preservation error drops approximately 325-fold in the selected GNLSE run with a Gaussian perturbation and by more than 100-fold on average over repeated initializations; and the mass drift drops 14-fold for the KS problem with a Gaussian perturbation. The framework thus tunes key components on tasks of different nature without changing its internal structure.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

151,044 characters extracted from source content.

Expand or collapse full text

[orcid=0000-0002-5000-7257] [orcid=0009-0002-9640-6787] [orcid=0000-0002-6921-4133] [orcid=0000-0001-9399-9332] [orcid=0000-0002-5595-6398] [orcid=0009-0007-0430-7658] †highlights: A multi-agent large language model framework tunes functional blocks of specialized neural networks. The same framework transfers across five spiking classification and physics-informed neural network tasks. MNIST accuracy rises +13 p; selected GNLSE errors drop up to 2500x and 325x. Detailed prompts substantially improve performance of less capable large language models. A Multi-Agent LLM System for Automated Tuning of Specialized Neural Network Models Ivan Moloshnikov ivan-rus@yandex.ru Dmitry Kunitsyn Kunitsyn_DE@nrcki.ru Alexander Sboev sag111@mail.ru Maksim Skorokhodov Skorokhodov_MS@nrcki.ru Roman Rybka Rybka_RB@nrcki.ru Yuanchao Liu ycliuharbin@163.com organization=National Research Centre “Kurchatov Institute”, addressline=1 Akademika Kurchatova Square, city=Moscow, postcode=123182, state=Moscow, country=Russia organization=National Research Nuclear University MEPhI, addressline=31 Kashirskoye Shosse, city=Moscow, postcode=115409, state=Moscow, country=Russia organization=Moscow Institute of Physics and Technology, addressline=9 Institutskiy Lane, city=Dolgoprudny, postcode=141701, state=Moscow, country=Russia organization=Harbin Institute of Technology, addressline=No. 92 West Dazhi Street, Nangang District, city=Harbin, postcode=150001, state=Heilongjiang, country=China Abstract Many machine learning approaches are two-level: an expert selects a class of methods, within this class the choice of key components is critical and usually manual. This paper addresses the automatic tuning of these components and proposes a framework transferable across models of different nature. The multi-agent large language model (LLM) framework has three roles—analyst, coder, and recorder —with an experiment log and a sandboxed execution environment. Unlike systems that optimize the entire ML script or search architectures in a discrete space, the framework operates over an open but semantically narrow space: the agent generates code only for the specified components, leaving the overall network structure unchanged. The approach is validated in two directions: tuning a classification spiking neural network with correlation learning, and solving high-order differential equations with Physics-Informed Neural Networks (PINN). Across five tasks—spiking classification on MNIST and Digits, and boundary-value problems for the generalized nonlinear Schrödinger equation (GNLSE) and Kuramoto–Sivashinsky (KS) equations—the agent outperforms the baselines: classification accuracy is 13 percentage points (MNIST) and 6 percentage points (Digits) higher; the error against the analytical single-soliton solution is reduced approximately 2491-fold; the invariant-preservation error drops approximately 325-fold in the selected GNLSE run with a Gaussian perturbation and by more than 100-fold on average over repeated initializations; and the mass drift drops 14-fold for the KS problem with a Gaussian perturbation. The framework thus tunes key components on tasks of different nature without changing its internal structure. keywords Large language models ,Multi-agent systems ,Automated machine learning ,Physics-informed neural networks ,Spiking neural networks †credit: Conceptualization, Methodology, Software, Investigation, Data curation, Writing – original draft†credit: Methodology, Investigation, Software, Writing – original draft†credit: Conceptualization, Methodology, Investigation, Data curation, Supervision, Validation, Writing – original draft†credit: Software, Validation, Investigation, Visualization, Writing – original draft†credit: Supervision, Project administration†credit: Formal analysis, Methodology†corresponding: Corresponding author†corresponding: Principal corresponding author 1 Introduction Large language models (LLMs) have become the foundation of a new class of tools —agent systems that can not only respond to queries but also plan actions, use tools, execute code, and accumulate experience from previous steps Xi et al. (2023); Yao et al. (2023b). Such agents are already used in programming, ML engineering, and individual scientific disciplines Jimenez et al. (2024); Chan et al. (2024); Lu et al. (2024); Schmidgall et al. (2025); however, their use as a meta-optimizer for specialized neural network models remains understudied. A detailed review of related work is given in Section 2. The quality of specialized neural network models depends strongly on the task: the representation of inputs and outputs, the method for decoding internal activity, the balance of loss-function terms, and other components must be selected anew for each problem, and no universal selection rule exists. In practice, this selection is performed manually, which makes the process labor-intensive and poorly reproducible. At the same time, models of different types can be decomposed into similar modules: input encoding and preprocessing, a trainable feature transformation, a decoding or output transformation, and a loss function. The tuning task is therefore reduced not to redesigning the entire network for every new problem but to modifying its individual modules, which enables a single tuning procedure applicable to models of different types. Studies in which LLM agents are used as a general tool for automatic tuning of individual functional blocks of a model with a fixed architecture and training procedure are almost absent from the open literature (see Section 2). Existing systems focus either on programming tasks and the general ML cycle or on constructing a network from a textual task description. In contrast, the proposed framework does not build a network but modifies its predefined blocks using a single agent interaction pattern that does not depend on the particular application task. The main contributions are as follows. 1. We propose a multi-agent LLM framework for the automatic tuning of individual functional blocks of a neural network model (Section 3). 2. We show that the same framework transfers across five tasks with different types of information representation without changes to its internal structure (Section 5, Table 1). 3. We show that a detailed prompt improves the performance of less capable LLMs, bringing it closer to that of more advanced models (Section 5.2, Table 2). 4. The multi-agent framework improves the baseline models: on the spiking neural network (SNN) tasks (Tables 2 and 19), accuracy increases by 13% on MNIST and by 6% on the scikit-learn digits; on the physics-informed neural network (PINN) tasks (Section 5.4–5.6), the relative error decreases approximately 2500-fold for the generalized nonlinear Schrödinger equation (GNLSE) problem with a known analytical solution, and for problems without analytical solutions, the conservation error decreases approximately 325-fold in the selected GNLSE soliton–Gaussian run (more than 100-fold on average over repeated initializations) and the mass drift decreases 14-fold for the Kuramoto–Sivashinsky (KS) wave with a Gaussian perturbation. 2 Related Work LLM agents can autonomously plan actions, use tools, and accumulate knowledge while solving multi-step tasks Xi et al. (2023). The ReAct (Reasoning + Acting) pattern has become a baseline architecture Yao et al. (2023b). Recent surveys Guo et al. (2024); Huang et al. (2024) identify typical components: a planning module (Chain-of-Thought Wei et al. (2022), Tree-of-Thoughts Yao et al. (2023a)), dynamic memory and self-reflection Shinn et al. (2023), calls to external tools, and multi-agent designs with role separation Chen et al. (2024). When agents act as meta-optimizers of neural network models, three of these components are essential: role separation distinguishes hypothesis formulation from its software implementation; logged memory prevents rejected variants from being proposed again; and a sandboxed execution environment is required to run generated code safely. LLM agents have been applied to a diverse set of tasks. SWE-bench Jimenez et al. (2024) and MLE-Bench Chan et al. (2024) have become standard benchmarks for programming and automated ML: the medal rate of the first published o1-preview configuration with AIDE was 16.9 % Chan et al. (2024), while later agents increased this value to 64 % of medals on MLE-Bench Lite (MLE-STAR Nam et al. (2026), including 36 % gold medals) and to approximately 35 % across all difficulty levels (R&D-Agent on GPT-5) Yang et al. (2025). Similar multi-agent designs are used in materials discovery and biochemistry Takahara et al. (2025); Shi et al. (2026); Ghafarollahi and Buehler (2024); Swanson et al. (2025). Voyager Wang et al. (2023) demonstrates lifelong accumulation of a skill library in the form of executable code, which is close to the idea of logged hypotheses used in this work. Systems for automating the scientific cycle (The AI Scientist Lu et al. (2024); Yamada et al. (2025), Agent Laboratory Schmidgall et al. (2025)) and agent-based robot control Skorokhodov et al. (2025) have been developed in parallel; they show that autonomous research is possible in principle but do not incorporate domain-specific knowledge of PINNs or SNNs. PINN quality depends on the choice of architecture, activations, the balance of loss-function terms, the distribution of collocation points, and the representation of inputs and outputs Moloshnikov et al. (2025b); separate heuristics have been proposed for each component, but no unified automatic procedure has been established. Lang-PINN He et al. (2025), the study closest to the agent-based approach, constructs a PINN from a textual task description but focuses on network construction rather than on selecting its functional blocks for a physical problem. For SNNs, automatic approaches to selecting the architecture and its individual blocks rely on evolutionary algorithms and reinforcement learning Kim et al. (2022); Na et al. (2022); formally, neural architecture search (NAS) also includes the generation of individual functional modules within a predefined space of architectural variants. Unlike NAS, the proposed approach uses an LLM, exploits domain knowledge and knowledge of model properties, locally modifies a selected module of an existing network without searching a fixed space of variants, and stores the history of hypotheses in natural language, so that each subsequent iteration can use a substantive analysis of previous iterations rather than only their numerical results. To the authors’ knowledge, the open literature contains no work in which LLM agents are used as a general tool for automatically tuning individual functional blocks of neural network models without redesigning the network itself. MLE-STAR Nam et al. (2026) is the closest agent system; the similarity is limited to two components: logged experiment memory and a sandboxed execution environment. MLE-STAR applies iterative code modifications with an update only when the metric improves on the outer loop, processes the entire ML script, selects the target block with a separate agent, assigns web search to a separate retriever agent, and involves more than ten agents. 3 Method: Multi-agent LLM framework The framework optimizes a selected functional block while keeping the rest of the task pipeline fixed. Before a run, the Human Researcher specifies the task, the block that the agents may modify, the evaluation metric, the software contract for the generated code, and a baseline. These settings allow the same framework to be applied to different problems without changing its architecture or execution loop. The framework is built on top of an OpenAI-compatible API and separates hypothesis generation, software implementation, and experiment documentation into three LLM roles. A sandboxed execution environment and an external experiment log close the iterative loop (Fig. 1). The framework consists of five main components: • LLM Analyst receives the task prompt and the current experiment log. It proposes one hypothesis for the agent-modifiable block and expresses the idea as a technical specification with a fixed structure. If the Human Researcher provides a taxonomy of method classes, the Analyst also assigns the hypothesis to one of these classes. • LLM Coder translates the technical specification into Python code. In addition to implementing the proposed method, the code must satisfy a task-specific software contract that defines the allowed dependencies, class names, call signatures, input and output shapes, and trainable parameters. • Sandboxed execution environment (Docker) inserts the generated module into the fixed task pipeline, first performs a short functional test, and then evaluates a valid candidate. A runtime error is returned to the LLM Coder for correction, whereas the metric from a successful run is passed to the logging stage. • LLM Recorder summarizes the specification, validated code, and evaluation result in a structured JSON record, including the method name, short description, and final hyperparameters. • Experiment log is a table of structured records that serves as the external memory of the framework. Each record is assembled from the outputs of the pipeline stages and contains: – a run identifier assigned by the framework; – the experiment name and a short description generated by the LLM Recorder from the technical specification and validated code; – the execution status assigned by the control loop from the functional-test and evaluation outcomes; – the resulting target-metric value computed by the task-specific evaluation procedure for a successfully evaluated candidate; – the method class selected by the LLM Analyst when a task-specific taxonomy is used; – the final hyperparameter values extracted by the LLM Recorder from the validated generated code. Records with the With error status remain stored in the complete experiment log. Before the log is included in the Analyst prompt, the framework filters out these records, so the Analyst receives the baseline and successfully evaluated records only. The run budget counts only successfully evaluated records with the OK status. The technical specification describes the method proposed by the Analyst, whereas the software contract defines the interface and constraints for the generated code. The agents may modify only the selected block; the rest of the task pipeline remains fixed. This separation allows generated code to be inserted and tested automatically. Applying the framework to a new task requires a new prompt, software contract, target metric, and baseline, but the agent roles and execution loop remain unchanged. Figure 1: Architecture of the multi-agent framework. The Human Researcher configures the task, agent-modifiable block(s), target metric, and software contract once. At each iteration, the LLM Analyst reads the experiment log and produces a technical specification, which the LLM Coder converts into contract-constrained code. The generated module is first functionally tested and then evaluated in a Docker sandbox; runtime errors are returned to the Coder, while a successful evaluation is passed to the LLM Recorder. The Recorder appends a structured JSON record to the experiment log, which serves as the external memory for the next iteration. The detailed-prompt mode additionally provides the Analyst with a method-class taxonomy and a dynamically updated coverage map. Only task-specific prompts, the software contract, the editable model block, and the target metric change between tasks; the agent roles and loop structure remain unchanged. Baseline record. Each framework run starts with a task-specific baseline record in the experiment log. This record evaluates the fixed pipeline with the simple baseline implementation of the modifiable block and provides the initial reference value of the target metric. The baseline implementations used in this study are defined in Section 4. Standard and detailed Analyst prompts. The framework supports two modes of hypothesis generation. The standard prompt contains the Analyst role, problem statement, target metric, implementation constraints, required format of the technical specification, and complete experiment log. The detailed prompt retains these elements and adds a structured Chain-of-Thought procedure and a list of method classes specified by the Human Researcher for systematic exploration. Each hypothesis is assigned to one principal class, which is also stored in its experiment-log record, and the log is converted into a dynamically updated coverage map. For every class, the map reports the number of evaluated hypotheses, the best metric value, and a coverage status such as “not explored,” “limited data,” or “explored.” The map is included in the next Analyst prompt, which prioritizes suitable classes with little or no coverage and revisits an explored class only when a substantially different hypothesis is proposed. This mechanism structures an open natural-language search without reducing it to a finite list of implementations. The task-specific prompt configurations are described in Section 4; the detailed-mode prompt materials are provided in Supplementary Materials. 3.1 Run protocol One framework run is a series with a task-specific number of successful iterations. At each iteration, one hypothesis is tested in six steps: 1. The Analyst receives the task prompt and the current experiment log. At the first iteration, the log contains only the baseline record. In detailed mode, the current coverage map is generated from the log and included in the prompt. 2. The Analyst creates a technical specification according to the prompt template. In detailed mode, the specification also identifies the method class. 3. The Coder implements the specification subject to the corresponding software contract. 4. The code undergoes an automatic functional test on a small sample in the sandboxed environment. If the test fails, the error message, generated code, and original specification are returned to the Coder. The Coder analyzes the reported error and rewrites the code while preserving the proposed method and software contract; the revised version is then tested again. If valid code cannot be obtained within 10 attempts, the experiment is assigned the With error status; it does not count toward the run budget, and the framework proceeds to a new hypothesis. 5. Code that passes the functional test is evaluated according to the task-specific evaluation protocol, and the resulting metric is assigned the OK status. The budgets, metrics, data-use rules, and repetition schemes used in this study are specified in Section 4.5 and the accompanying appendices. 6. The Recorder combines the specification, validated code, execution result, and fitted hyperparameters into a JSON record and appends it to the log. The updated log becomes the external memory for the next iteration. 4 Experimental setup The same framework is applied to five tasks with different types of information representation: two spiking-network classification tasks, two boundary-value problems for the generalized nonlinear Schrödinger equation, and one boundary-value problem for the generalized Kuramoto–Sivashinsky equation. 4.1 SNN task: Correlation Classwise Network decoder The classifier is based on a Correlation Classwise Network (CCN) spiking neural network Sboev et al. (2025); Kunitsyn et al. (2025), whose weights are trained using a local spike-timing-dependent plasticity (STDP) rule. Image pixels are normalized and encoded as correlated spike trains that serve as the network input. The output spike activity of the leaky integrate-and-fire (LIF) neurons is then translated into a class label by the decoding module. The fixed part, which the framework cannot modify, comprises the input encoding, the CCN, and the recorded output activity vector ∈ℝKr ^K; the agent-modifiable part is the Decoder⁡()Decoder(r) module (Fig. 2). Its software contract permits fitting parameters from training activity and class labels but prohibits modifying the trained SNN, spike encoding, or recorded activity and prohibits error backpropagation. The decoder must return one valid class label per sample. We evaluate the method on two datasets: MNIST LeCun et al. (1998) (28×2828× 28 images, N=784N=784 spike generators) and the scikit-learn digits dataset Pedregosa et al. (2011) (8×88× 8, N=64N=64). Both contain 10 classes, so the network has K=10K=10 output LIF neurons. Accuracy is used as the target metric. As a baseline, we use an arg⁡max decoder: the predicted class corresponds to the output neuron with the highest spike rate. Candidate selection uses only development data; the held-out test sets are unavailable to the agents, as detailed in Section B.1. Full details of the SNN model, the encoding scheme, the STDP rule, and the hyperparameters are given in Appendix B. Figure 2: SNN pipeline. The input encoding, the SNN trained with STDP, and the output activity vector r (left and middle) are fixed. Only the decoder module Decoder⁡()Decoder(r) (right) is generated and optimized by the framework. 4.2 PINN tasks: GNLSE and Kuramoto–Sivashinsky Two equations are investigated in this study. The first is the generalized nonlinear Schrödinger equation with the nonlinear correction (1) from Zavertyaev et al. (2024); Moloshnikov et al. (2025a), and the second is the generalized Kuramoto–Sivashinsky equation (6) from Kudryashov (2026). iqt+qxx+|q|2q(1−α|q|2+β|q|4)=0,iq_t+q_x+|q|^2q (1-α|q|^2+β|q|^4 )=0, (1) where the complex field q⁡(x,t)=u⁡(x,t)+iv(x,t)q(x,t)=u(x,t)+iv(x,t) is approximated by a fully connected network, while the equation (ℒpdeL_pde), initial-condition (ℒicL_ic), and boundary-condition (ℒbcL_bc) residuals enter the loss functional (ℒL) according to the Physics-Informed Neural Networks framework Raissi et al. (2019), a detailed description is given in the appendix A. In all configurations, β=0β=0. The fixed part comprises the multilayer perceptron (MLP) architecture, loss functional, and training procedure. The agent-modifiable part comprises two transformations (Fig. 4): the input transformation FIF_I, which maps the coordinates (x,t)(x,t) to the MLP input features, and the output transformation FOF_O, which maps the MLP outputs and coordinates to the predicted physical field. For the GNLSE problems, the prediction consists of the two real components (u,v)(u,v); for the Kuramoto–Sivashinsky problem, it is the single real field u. The agents may change the mathematical forms of FIF_I and FOF_O and may increase the intermediate MLP output dimension when required by FOF_O, but cannot modify the remaining architecture, loss functional, or training procedure. In the baseline, FIF_I passes (x,t)(x,t) to the MLP unchanged and FOF_O returns the raw MLP outputs. The framework is applied to two boundary-value problems for (1), which differ in the initial condition, domain, and target metric, and to one boundary-value problem for (6). Figure 3: Reference solutions of the two GNLSE boundary-value problems, shown as a surface (top) and a space–time heat map (bottom). Left: the analytical single soliton. Right: the segmented PINN solution of the soliton with a Gaussian perturbation from Zavertyaev et al. (2024) (Figure 5). The Kuramoto–Sivashinsky problem is not shown because it is solved without a reference solution. Single-soliton problem: iqt+qxx+|q|2q(1−α|q|2+β|q|4)=0,q(x, 0)=q~(x, 0),q(−25,t)=0,q(25,t)=0,(x,t)∈[−25,25]×[0,3], casesi\,q_t+q_x+|q|^2q (1-α|q|^2+β|q|^4 )=0,\\[2.0pt] q(x,\,0)= q(x,\,0),\\[2.0pt] q(-25,\,t)=0, q(25,\,t)=0,\\[2.0pt] (x,t)∈[-25,25]×[0,3], cases (2) where α=0.1α=0.1, β=0β=0, and the initial profile q~ q is specified by the analytical soliton (12) with parameters k=1k=1, w=0.88w=0.88, x0=0x_0=0, θ0=0 _0=0, and μ=4(k2−w)μ=4(k^2-w). Because an analytical solution is known for this problem, quality is evaluated directly using the relative amplitude error Relh=‖htrue−hpred‖2/‖htrue‖2Rel_h=\|h_true-h_pred\|_2/\|h_true\|_2, h=u2+v2h= u^2+v^2. This problem provides a transparent test: the agent can rely on the explicit structure of the solution. Soliton problem with a Gaussian perturbation: iqt+qxx+|q|2q(1−0.3|q|2)=0,q(x, 0)=q2(x, 0)+0.2e−x2/50,q(−50,t)=0,q(50,t)=0,(x,t)∈[−50,50]×[0,25], casesi\,q_t+q_x+|q|^2q (1-0.3|q|^2 )=0,\\[2.0pt] q(x,\,0)=q_2(x,\,0)+0.2\,e^-x^2/50,\\[2.0pt] q(-50,\,t)=0, q(50,\,t)=0,\\[2.0pt] (x,t)∈[-50,50]×[0,25], cases (3) where the soliton component q2q_2 has the same form as (12) with α=0.3α=0.3, k=1k=1, w=0.88w=0.88, x0=−30x_0=-30, and θ0=0 _0=0. The Gaussian packet is Agexp[−(x−Bg)2/(2Cg2)]A_g [-(x-B_g)^2/(2C_g^2)] with Ag=0.2A_g=0.2, Bg=0B_g=0, and Cg=5C_g=5. No analytical solution is known, so the search is driven by two physical metrics. The mean absolute residual of the boundary-value problem is Resbvp=13(⟨|ℛGNLSE|⟩+⟨|qpred(x,0)−q⁡(x,0)|⟩+⟨|qpred|∂⟩),Res_bvp= 13 ( |R_GNLSE| + |q_pred(x,0)-q(x,0)| + |q_pred|_∂ ), (4) where ℛGNLSE=iqt+qxx+|q|2q(1−0.3|q|2)R_GNLSE=iq_t+q_x+|q|^2q(1-0.3|q|^2), the averages are taken over a uniform evaluation grid, and ∂ denotes the Dirichlet boundaries x=±50x=± 50. The second primary metric is the relative conservation error of the integrals of motion I1=∫|q|2xI_1= |q|^2\,dx and I2=∫(q∗qx−qqx∗)xI_2= (q^*q_x-q_x^*)\,dx. For each invariant, we compute its maximum relative drift from the initial value over the considered time interval, RelIj=maxt⁡|Ij(t)−Ij(0)||Ij(0)|,j∈1,2.Rel_I_j= _t |I_j(t)-I_j(0) | |I_j(0) |, j∈\1,2\. These individual errors are combined into the overall conservation metric RelI=RelI1+RelI2.Rel_I=Rel_I_1+Rel_I_2. A segmented PINN solution from Zavertyaev et al. (2024) is retained as an auxiliary amplitude reference Relh=‖href−hpred‖2/‖href‖2Rel_h=\|h_ref-h_pred\|_2/\|h_ref\|_2. On the agent domain that field itself satisfies Resbvp=9.11⋅10−5Res_bvp=9.11· 10^-5 and Reli=9.87⋅10−4Rel_i=9.87· 10^-4. RelhRel_h has lowest priority: if a hypothesis improves both ResbvpRes_bvp and ReliRel_i relative to this reference, it is regarded as a more accurate solution of the boundary-value problem than the stored PINN field, and RelhRel_h is disregarded. This problem is the principal substantive test of the method because the agent cannot simply embed a ready-made analytical solution. The complete formulations of both problems, the analytical reference soliton, the collocation-point sampling schemes, and the hyperparameters are provided in Appendix A. The second equation considered is a generalization of the Kuramoto–Sivashinsky equation with nonlinear convection Kudryashov (2026). This generalization is important in high-temperature plasma research, where the nonlinear dependence of thermal conductivity must be taken into account Liu et al. (1976); Sudan and Rosenbluth (1976); Chen et al. (1977). We write the source equation for the scalar field v⁡(x,t)v(x,t) in the normalization vt+vxx+σvxxx+vxxxx+γvvx+χ(v2)xx=0.v_t+v_x+σ v_x+v_x+γ v_x+χ(v^2)_x=0. (5) Here γ=2βγ=2β is introduced to avoid confusion with the coefficient β in (1). For β=3σβ=3σ and σ=−1σ=-1, the coefficients used in this study are γ=2β=6σ=−6γ=2β=6σ=-6 and χ=3χ=3. The corresponding exact traveling wave has the nonzero background −1/3-1/3. We therefore remove this background by introducing v⁡(x,t)=u⁡(x,t)−13,v(x,t)=u(x,t)- 13, so that the neural network approximates the zero-background scalar field u⁡(x,t)u(x,t). Under this change of variables, Eq. (5) becomes the differential equation in (6) (the complete substitution is given in Appendix A.2). Unlike (1), where the neural network approximates the complex field q⁡(x,t)q(x,t), this problem consequently has a single real-valued output, and the boundary conditions are zero for both the field and its second derivative. The general form of the loss functional follows the same three-term PINN scheme as for the GNLSE, but the MLP size, optimizer settings, and training budget differ; they are listed separately in Table 9. The agent-modifiable part comprises the same transformations FIF_I and FOF_O as in the GNLSE problems. Traveling-wave problem with a Gaussian perturbation: ut−uxx−uxxx+uxxxx−6uux+2ux+3(u2)xx=0,u⁡(x,0)=12cosh2⁡(x/2)+0.02exp⁡[−(x+6)28],u(−30,t)=uxx(−30,t)=0,u(40,t)=uxx(40,t)=0,(x,t)∈[−30,40]×[0,10]. casesu_t-u_x-u_x+u_x-6u_x+2u_x+3(u^2)_x=0,\\[2.0pt] u(x,0)= 12 ^2(x/2)+0.02\, \! [- (x+6)^28 ],\\[2.0pt] u(-30,t)=u_x(-30,t)=0, u(40,t)=u_x(40,t)=0,\\[2.0pt] (x,t)∈[-30,40]×[0,10]. cases (6) The unperturbed pulse has the form uex(x,t)=12cosh2⁡((x−t)/2)u_ex(x,t)= 12 ^2((x-t)/2) (speed 11, amplitude 1/21/2, starting at x0=0x_0=0); the Gaussian packet is located to the left (xg=−6x_g=-6, width w=2w=2) and catches up with the pulse at approximately t≈6t≈ 6, after which the perturbation intensifies. The complete analytical solution for this problem is unknown, so quality is assessed only by how well the network satisfies the boundary-value problem itself and by a physical diagnostic. The primary metric of the search is the mean absolute residual of the boundary-value problem, Resbvp=14(⟨|ℛKS|⟩+⟨|upred(x,0)−u⁡(x,0)|⟩+⟨|upred|∂⟩+⟨|∂xxupred|∂⟩),Res_bvp= 14 ( |R_KS| + |u_pred(x,0)-u(x,0)| + |u_pred|_∂ + | _xu_pred|_∂ ), (7) where ℛKSR_KS is the residual (22), the averages are taken over a uniform 1000×101000× 10 grid in the domain, and ∂ denotes the two zero boundaries. The physical diagnostic is the relative mass drift of the predicted solution: Δm(u)=maxk⁡|Mpred(tk)−Mpred(0)||Mpred(0)|,Mpred(t)=∫upred(x,t)x, _m^(u)= _k |M_pred(t_k)-M_pred(0) | |M_pred(0) |, M_pred(t)= u_pred(x,t)\,dx, (8) with the maximum taken over ten uniformly spaced time slices tk∈[0,10]t_k∈[0,10]. Figure 4: PINN architecture. The coordinates (x,t)(x,t) pass through the transformation layer FIF_I, the fully connected network, and the transformation layer FOF_O, producing predictions u^,v u, v. Automatic differentiation computes the derivatives for the loss function. The loss function (ℒL) consists of the residual in the differential equation (ℒpdeL_pde) and the errors in the initial (ℒicL_ic) and boundary conditions (ℒbcL_bc). The optimizer updates the network weights. 4.3 Task-specific prompt configurations The two universal Analyst modes defined in Section 3 are instantiated differently across the experiments. The PINN tasks use only the standard mode, with prompts that specify the differential equation, the modifiable transformations, the software contract, and the physical target metric; they do not use a method-family catalog, Chain-of-Thought procedure, or coverage map. The SNN classification experiments compare the standard mode, denoted the simple prompt in the result tables, with the detailed mode. The simple prompt specifies the Analyst role, classification problem, target accuracy, data-use and implementation constraints, examples of admissible decoder types, required technical-specification format, and complete experiment log. The detailed prompt adds a researcher-defined taxonomy of 12 decoder-method classes, a structured Chain-of-Thought procedure for coverage and failure analysis, an explicit class field in the technical specification, and the dynamically updated coverage map. The 12 classes define broad areas of exploration rather than 12 fixed algorithms; each class may contain multiple feature transformations, decision rules, and hyperparameter settings, and hypotheses may combine methods from different classes. The complete class taxonomy and prompt texts are provided in Supplementary Section F. 4.4 Search baselines and budget: SNN decoder on MNIST The quality reached by the framework depends both on the search performed by the LLM Analyst and on the constraints that the Human Researcher writes into the prompt. To separate these two contributions, we compare the framework with baselines that search the same space of decoders without an Analyst. All comparisons in this subsection are performed on MNIST, where the argmax baseline gives 79.88%79.88\,\%. Shared search space. The decoder catalog is a fixed discrete space of estimators with grids of their hyperparameters, together with the operators that build ensembles and multistage pipelines. It is specified once as a TPOT config_dict that is listed in Supplementary Section E.2 and its classes follow the decoder taxonomy of the detailed Analyst prompt (Section 4.3, Supplementary Materials). The Analyst prompt permits exactly the decoders contained in the catalog. The restriction reaches the systems in different ways: the agent receives it as text and can in principle depart from it, whereas the baselines receive it as a constraint on the sampler. We inspected the Experiment logs of all agent runs and found no record outside the catalog. The space of admissible decoders is therefore common to all systems, and they differ in the search strategy and in the selection protocol described below. Compared configurations. Three baseline configurations are compared with the framework, each with a budget of 250 evaluated decoders per run. • Random Choice is an LLM-based sampler that draws decoder configurations from the same catalog through independent calls to Qwen3-30B-A3B-Instruct-2507 deployed on a rented GPU, without a Chain-of-Thought stage and without a debugging cycle. Each call receives a compressed description of the catalog and returns one decoder configuration in JSON, sampled independently of the previous choices. The configuration is implemented and evaluated by the same sandboxed pipeline as an agent hypothesis. A call that yields a configuration outside the catalog, a configuration that fails the functional test, or a repetition of a configuration already recorded is discarded and does not enter the Experiment log. Calls are issued until the Experiment log contains 250 evaluated decoders. • TPOT Olson et al. (2016) classic is the standard genetic algorithm with a population of 25 individuals and 9 generations, which together with the initial population gives 25×(9+1)=25025×(9+1)=250 evaluated candidates. • TPOT sequential uses a population of one individual and 249 consecutive generations, 1×(249+1)=2501×(249+1)=250 candidates, so that the budget is preserved while population-level exploration is replaced by a single mutating lineage. In both TPOT modes the search may combine models and build multistage processing sequences, because the corresponding operators are part of the catalog. TPOT samples and mutates pipelines according to its internal grammar, whereas Random Choice draws every configuration independently under an instruction to sample uniformly over the catalog. Budget. For the framework, the budget is the number of successfully evaluated records with the OK status. Records with the With error status remain stored in the complete Experiment log but do not count toward the budget and are excluded from the filtered log view passed to the Analyst. For Random Choice, the budget is the number of successfully evaluated decoder records; invalid calls discarded as described above do not count toward it. For TPOT, the budget is the number of fitted pipelines reported by the optimizer. In every configuration the budget is 250 per run, and every configuration was run 10 times. Within one run we take the best accuracy over the 250 evaluated decoders and the position at which it was first reached; Table 3 reports the smallest and the largest of these ten per-run values. Selection and evaluation. The baselines were run on the same trained SNN instances and the same output spike matrices as the prompt comparisons of Table 2, so every baseline run is matched with the corresponding agent run and the comparison is paired. The internal scoring rules of the compared systems are not identical. In the framework and in Random Choice, a candidate is scored on a validation set formed from the training set, and the value entered in the Experiment log (for the Analyst) is the mean over three random seeds. The held-out test set is not available to the Analyst (Section 4.5). TPOT selects its best pipeline by mean accuracy in five-fold cross-validation on the training set, after which the selected pipeline is retrained on the entire training set and evaluated once on the test set. The comparison in Table 3 is therefore made on the test accuracy of the best decoder of a run, which is defined in the same way for all systems. 4.5 Evaluation protocol Run budget and thresholds. The specified run budget (50 for the PINN tasks and 250 for the SNN tasks) counts only successful iterations; see Section 3.1 for the iteration steps. The target metric threshold (accuracy above 90 % for MNIST and 93 % for Digits; a decrease in RelhRel_h relative to the baseline for the single-soliton GNLSE problem, a decrease in ResbvpRes_bvp and ReliRel_i for the GNLSE problem with a Gaussian perturbation, and a decrease in ResbvpRes_bvp for the Kuramoto–Sivashinsky problem) provides a reference point for the analyst but is not a stopping criterion. Each hypothesis is evaluated with three random seeds in the single-soliton problem and with a single seed in the remaining two PINN problems, where one training run takes tens of minutes. Compared configurations. In both GNLSE problems, the baseline with identity FI,FOF_I,F_O transformations is compared with three agent configurations: Perplexity Sonar (sonar-pro) Perplexity AI (2026), Qwen3-30B-A3B-Instruct-2507 (qwen-qwen3-30b-a3b-instruct-2507) Qwen Team (2025), and Claude Opus 4.6 (anthropic-claude-opus-4.6); the Analyst prompt is the same for different LLMs. Because a complete agent run for the Kuramoto–Sivashinsky problem is computationally expensive, we evaluated only Claude Opus 4.6, which produced the best results on the other tasks. The comparison for this problem is therefore limited to the Claude configuration and its baseline; different LLM models are not compared with one another. The Qwen3-30B-A3B-Instruct-2507 weights Qwen Team (2025) were downloaded from the model’s public repository. In the classification tasks, pairs ⟨ model, prompt⟩ are compared: Claude Opus 4.6, Perplexity Sonar, and Qwen3-30B-A3B-Instruct-2507, each with the simple and detailed prompts defined in Section 4.3; the baseline is an argmax decoder based on the output spike rates of the LIF neurons. Ablation runs with the detailed prompt and the coverage-map section removed were also performed for all three models, while retaining the method catalog and Chain-of-Thought. Accuracy comparisons are provided in Tables 13, 14, 15, 16, 17, and 18. Metrics. In the single-soliton problem, RelhRel_h is the final characteristic of a hypothesis; the best RelhRel_h (mean over three internal seeds) is reported. In the GNLSE problem with a Gaussian perturbation, hypotheses are ranked by the conservation error ReliRel_i, while the boundary-value-problem residual ResbvpRes_bvp (4) is reported alongside it. The robustness of the best configurations was additionally tested using repeated runs (Appendix D). In the Kuramoto–Sivashinsky problem, hypotheses are ranked by the boundary-value-problem residual ResbvpRes_bvp (7), and the mass drift Δm(u) _m^(u) (8) is reported as an independent physical diagnostic. In the classification tasks, the target metric is test accuracy, but during a run accuracy is evaluated on a validation set formed from the training set so that test-set accuracy is not entered in the experiment log. Baseline model accuracies. In the single-soliton problem, the identity-map baselines give Relh=4.69⋅10−2Rel_h=4.69· 10^-2, 5.51⋅10−25.51· 10^-2, and 4.85⋅10−24.85· 10^-2 in the Sonar, Claude, and Qwen agent runs, respectively. In the GNLSE problem with a Gaussian perturbation, the corresponding single-run baselines give Reli=1.28⋅10−1Rel_i=1.28· 10^-1, 2.15⋅10−12.15· 10^-1, and 1.85⋅10−11.85· 10^-1; over n=10n=10 additional baseline trainings, mean is Reli=(1.65±0.18)⋅10−1Rel_i=(1.65± 0.18)· 10^-1 (Table 11). In the Kuramoto–Sivashinsky problem, the baseline gives Resbvp=1.41⋅10−3Res_bvp=1.41· 10^-3 and Δm(u)=9.31⋅10−2 _m^(u)=9.31· 10^-2. The SNN baseline (argmax over output spike rates, without Gaussian receptive fields (GRF) or bagging) gives an accuracy of 79.88 % on MNIST and 90.83 % on digits, consistent with previously published results Sboev et al. (2025). The gap between the datasets motivated target accuracies of 90 % on MNIST and 93 % on digits. 4.6 Statistical protocol Repeated MNIST runs and statistical analysis. To assess the effect of prompt structure, 10 complete runs, each comprising 250 iterations of agent search, were performed for every configuration. Two paired comparisons were performed for each of the three models: the simple prompt was compared with the full detailed prompt, and the detailed prompt with a coverage map was compared with its ablated version without the map. Each comparison comprised ten matched pairs of complete agent runs. Within each pair, the same trained SNN instance and identical output spike data were used, so differences between configurations were evaluated from paired differences in the maximum accuracy achieved. For Perplexity Sonar, the difference between prompt variants was consistent across all ten runs: with the simple prompt, the maximum accuracy remained at the baseline level, whereas with the detailed prompt it exceeded 90 % in every case. For Qwen3-30B-A3B-Instruct-2507, the detailed prompt also increased the maxima relative to the baseline in all ten runs, although the magnitude of the improvement varied considerably. A paired statistical analysis over ten matched runs was performed for all three analyst models. The prompt effect was tested using a two-sided exact Wilcoxon signed-rank test; two tests were performed for each of the three analyst models—“detailed versus simple” and “detailed versus detailed without the coverage map”—giving six tests to which the Holm correction for multiple comparisons was applied. Pairs were formed from runs with the same seed, which determines the framework initialization and the order of data presentation; in the Supplementary tables, the rows are ordered by decreasing accuracy of the detailed prompt rather than by seed number. For configurations in which the accuracy of the simple prompt exactly matched the baseline in every run, the paired test degenerates into a sign test, as reflected in the reported W and p values. The paired-analysis methodology and rationale for selecting the test are provided in Supplementary Section E. The unadjusted and Holm-adjusted p values, together with the W statistic for all six tests, are reported in Supplementary Section E; the difference is statistically significant in all six tests at α=0.05α=0.05 after the Holm correction. Repeated GNLSE trainings. For the soliton with a Gaussian perturbation, the identity baseline and the two Claude configurations selected by ResbvpRes_bvp and ReliRel_i were each retrained ten times from independent random initializations. The random seeds were not recorded, so the three series are unpaired. The exported repeated-run data contain ReliRel_i and its two components but not ResbvpRes_bvp; consequently, this experiment characterizes robustness of invariant preservation only and does not establish a repeated-run success rate for the joint criterion. 5 Results 5.1 Summary across five tasks The final results obtained by applying the same agent framework to five tasks are summarized in Table 1. Accuracy (higher is better) is used for the classification tasks, whereas a relative error or physical diagnostic (lower is better) is used for the PINN tasks. The improvement over the baseline is expressed in task-specific units: for classification tasks, the absolute accuracy increase Acc−AccbaseAcc-Acc_base is reported in percentage points, whereas for PINN tasks the error-reduction factor Errbase/ErrErr_base/Err is reported. The accuracy values and improvements shown in the table are rounded to integers. The corresponding unrounded accuracies are 79.88 and 92.83 % for MNIST (Table 2) and 90.83 and 96.67 % for digits (Table 19). For the GNLSE soliton with a Gaussian perturbation, Table 1 reports the conservation error ReliRel_i. For the Kuramoto–Sivashinsky problem, it reports the mass drift Δm(u) _m^(u), an independent physical diagnostic recorded alongside the optimized boundary-value-problem residual ResbvpRes_bvp; the latter decreases by a factor of 9.3 in the same run (Table 7). Each row compares a single best agent result with the corresponding identity baseline from the same analyst run. For the GNLSE soliton with a Gaussian perturbation, robustness was additionally evaluated over ten independent trainings. These runs give Reli=(1.65±0.18)⋅10−1Rel_i=(1.65± 0.18)· 10^-1 for the baseline, (1.42±0.59)⋅10−3(1.42± 0.59)· 10^-3 for the configuration selected by ResbvpRes_bvp, and (1.60±0.81)⋅10−3(1.60± 0.81)· 10^-3 for the configuration selected by ReliRel_i. The corresponding reductions relative to the baseline mean are 116×116× and 103×103×, respectively (Table 11). Overall, the best agent configuration outperforms the baseline on the reported metric in all five tasks. Table 1: Baseline and best agent results across five tasks. Task (metric) Baseline (single run) Best agent result (best in series) Improvement over baseline (p / factor) SNN, MNIST (Acc) 80 % 93 % +13+13 p SNN, scikit-learn digits (Acc) 91 % 97 % +6+6 p PINN, GNLSE, single soliton (RelhRel_h) 5.51⋅10−25.51· 10^-2 2.21⋅10−52.21· 10^-5 2491×2491× PINN, GNLSE, soliton + Gaussian (ReliRel_i) 2.15⋅10−12.15· 10^-1 6.61⋅10−46.61· 10^-4 325×325× PINN, KS, traveling wave + Gaussian (Δm(u) _m^(u)) 9.31⋅10−29.31· 10^-2 6.56⋅10−36.56· 10^-3 14×14× 5.2 SNN classification results The agent configurations on MNIST are compared in Table 2. For Claude Opus 4.6, the detailed prompt raises the mean maximum accuracy over ten runs paired by seed from 91.13±0.44%91.13± 0.44\,\% to 92.15±0.62%92.15± 0.62\,\%, and the best result from 91.65%91.65\,\% to 92.83%92.83\,\%. It gives the higher run maximum in eight of the ten pairs, and the two-sided exact Wilcoxon signed-rank test yields W=5W=5, p=0.0195p=0.0195, which remains significant after the Holm correction over the six tests of Section 4.6. The paired gain averages 1.02±0.901.02± 0.90 percentage points (Table 17), so for this model the effect is comparable to the spread between runs. For Perplexity Sonar and Qwen3-30B-A3B-Instruct-2507, none of the 250 decoders proposed per run under the simple prompt outperformed the baseline, so the run maximum for both models matched the baseline of 79.88%79.88\,\% in all ten runs (std =0.00=0.00). The detailed prompt lifted both models above the baseline, though to different extents: for Perplexity Sonar, the individual-run maxima ranged from 90.93%90.93\,\% to 92.70%92.70\,\% with a mean of 91.76±0.56%91.76± 0.56\,\%, so the 90%90\,\% reference point was exceeded in all ten runs, whereas for Qwen3-30B-A3B-Instruct-2507 the mean maximum was 85.88±2.42%85.88± 2.42\,\% and the best run reached 89.16%89.16\,\% without crossing it. For both models the prompt effect is significant by the same paired test (W=0W=0, p=0.00195p=0.00195, the minimum attainable value at n=10n=10; Tables 15 and 13). On the scikit-learn digits dataset, the effect of the detailed prompt depended on the analyst model. For Claude Opus 4.6, the mean run maximum increased from 95.07±0.52%95.07± 0.52\,\% to 96.12±0.66%96.12± 0.66\,\%, with all ten runs reaching 93 %. For Perplexity Sonar, the mean changed from 93.44±0.57%93.44± 0.57\,\% to 93.31±0.93%93.31± 0.93\,\%, while the number of runs reaching 93 % decreased from 9/10 to 8/10, although the best observed result increased from 93.89 % to 94.17 %. For Qwen3-30B-A3B-Instruct-2507, the mean increased from 91.63±0.74%91.63± 0.74\,\% to 93.66±0.79%93.66± 0.79\,\%, and the number of runs reaching 93 % increased from 0/10 to 7/10. Table 2: Comparison of analyst models and prompt variants on MNIST. For each configuration, the table reports the number of complete agent runs, the mean and standard deviation of the maximum accuracy achieved within an individual run, the highest accuracy across all runs, and the number of runs in which the 90 % threshold was reached. Ten matched runs were performed for every configuration. Model Prompt Mean maximum accuracy ± std, % Best accuracy, % Runs with accuracy ≥90%≥ 90\,\% Claude Opus 4.6 simple 91.13±0.4491.13± 0.44 91.65 10 out of 10 Claude Opus 4.6 detailed 92.15±0.6292.15± 0.62 92.83 10 out of 10 Claude Opus 4.6 detailed without coverage map 91.02±0.2591.02± 0.25 91.45 10 out of 10 Perplexity Sonar simple 79.88±0.0079.88± 0.00 79.88 0 out of 10 Perplexity Sonar detailed 91.76±0.5691.76± 0.56 92.70 10 out of 10 Perplexity Sonar detailed without coverage map 82.34±0.7082.34± 0.70 83.18 0 out of 10 Qwen3-30B-A3B-Instruct-2507 simple 79.88±0.0079.88± 0.00 79.88 0 out of 10 Qwen3-30B-A3B-Instruct-2507 detailed 85.88±2.4285.88± 2.42 89.16 0 out of 10 Qwen3-30B-A3B-Instruct-2507 detailed without coverage map 80.18±0.4580.18± 0.45 81.15 0 out of 10 Table 3: Comparison of the agent system, TPOT evolutionary search (classic and sequential modes), and the Random Choice baseline for the SNN decoder on MNIST with a budget of 250 evaluated models. Method Mode AccminAcc_min, % AccmaxAcc_max, % kmink_min kmaxk_max TPOT classic 89.60 90.15 201 243 TPOT sequential 89.38 89.93 67 210 Random Choice LLM sampler 89.48 89.93 - - Qwen3-30B- A3B-Instruct-2507 detailed 82.28 89.16 15 209 Perplexity Sonar detailed 90.93 92.70 8 147 Claude Opus 4.6 detailed 91.22 92.83 2 43 Every Method was run 10 times. Within one run we take the best accuracy over the 250 evaluated models. The minimum and maximum reported values were derived from these 10 values: • AccminAcc_min, AccmaxAcc_max — the smallest and the largest best accuracies. • kmink_min, kmaxk_max — the smallest and the largest positions of evaluated models at which the best accuracy was first reached within a single run. Random Choice issues independent calls, so the position at which the run maximum appears does not characterise search progress and kmink_ and kmaxk_ are not reported. 5.3 Cost of the agent search The positions kmink_ and kmaxk_ in Table 3 characterize the length of the search, not its computational cost. The framework and the Random Choice baseline are therefore compared by three quantities: the token usage per run, the cost of a complete run CrunC_run, and the cost C∗C^* accumulated up to the earliest position kmink_ . The token usage depends on the structure of the Analyst input, on the number of Coder attempts per iteration (Section 3.1), and on the size of the Recorder input. Rates. The Coder and the Recorder use gpt-5-mini in all configurations, at the provider’s public rate $0.25/$2.00 per M input/output tokens. For the Analyst, the rates in effect in June 2026 were $5/$25 per M input/output tokens for Claude Opus 4.6 and $1/$1 per M for Perplexity Sonar. Qwen3-30B-A3B-Instruct-2507 is deployed on a rented GPU and has no provider price, so its cost and the cost of Random Choice are given at the managed hosting rate for the same model class, $0.15/$0.60 per M tokens. Random Choice uses the same Qwen model in a single-agent configuration with 250 evaluated decoders without a Chain-of-Thought stage and without a debugging cycle. Its token usage was measured from the run logs in the same way as for the agent configurations. Table 4: Performance, resources, and cost of the agent framework and the Random Choice baseline on MNIST (detailed prompt, budget of 250 evaluated decoders per run). The performance block is taken from Table 3. Metric Claude det. Sonar det. Qwen det. Random Choice Performance AccminAcc_ , % 91.2291.22 90.9390.93 82.2882.28 89.4889.48 AccmaxAcc_ , % 92.8392.83 92.7092.70 89.1689.16 89.9389.93 kmink_ 22 88 1515 — kmaxk_ 4343 147147 209209 — Resources per run, tokens in millions (M) Evaluated decoders 250250 250250 250250 250250 Analyst calls, NanN_an 265265 252252 255255 347347 Programmer calls / iter 3.313.31 1.701.70 2.052.05 — Analyst in / out, M 11.1011.10 / 1.171.17 10.8910.89 / 0.580.58 11.2011.20 / 0.580.58 0.280.28 / 0.0420.042 Programmer in / out, M 12.9812.98 / 4.664.66 2.402.40 / 1.731.73 2.792.79 / 1.321.32 — Recorder in / out, M 2.612.61 / 0.0580.058 1.711.71 / 0.0610.061 1.381.38 / 0.0630.063 — Cost Analyst rate, $/M in / out 55 / 2525 11 / 11 0.150.15 / 0.600.60 0.150.15 / 0.600.60 Programmer, Recorder rate, $/M 0.250.25 / 2.002.00 (gpt-5-mini) — CanC_an, $ 84.6684.66 11.4711.47 2.032.03 0.070.07 CprC_pr, $ 12.5612.56 4.064.06 3.353.35 — CrecC_rec, $ 0.770.77 0.550.55 0.470.47 — CrunC_run, $ 97.9997.99 16.0816.08 5.855.85 0.070.07 C∗C^* up to kmink_ , $ 0.370.37 0.210.21 0.260.26 0.0220.022 Token usage. The cost of a run is the sum of the token costs of the three roles of the framework, the Analyst, the Programmer, and the Recorder. For a role r, Cr=prinTrin+proutTrout,Crun=Can+Cpr+Crec,C_r=p_r^in\,T_r^in+p_r^out\,T_r^out, C_run=C_an+C_pr+C_rec, (9) where prp_r are the per-token rates and TrT_r are the token counts taken from the run logs. The Analyst input at iteration i consists of the fixed part of the detailed prompt (the role, problem statement, constraints, catalog of 12 method classes, coverage map, Chain-of-Thought instruction, and decoder examples), F=4085F=4085 tokens, and the experiment log of i records of dlogd_ tokens each. A log record is passed to the Analyst as a JSON object with the fields run_id, name, description, accuracy, and method_class, so its footprint dlogd_ in the prompt exceeds the length d of the description text itself, which is what the Recorder returns. Over the N=250N=250 accepted iterations of a run the Analyst also makes Nan−N_an-N additional calls that failed to produce an accepted decoder; at each such call the Analyst reads the log at its current length. Approximating the mean log length at the moment of a failed call by dlogN/2d_ N/2 — the average log length over the run — the total Analyst input is Tanin=NanF+dlog(N⁡(N+1)2+(Nan−N)N2).T_an^in=N_anF+d_ ( N(N+1)2+(N_an-N)\, N2 ). (10) The first term is the fixed prompt part paid at each of the NanN_an Analyst calls; the second term collects the log tokens read during the N accepted iterations (exact sum) and during the Nan−N_an-N failed calls (approximate). The mean length of a log record in the Analyst prompt, measured with the cl100k_base tokenizer as the increment of the prompt length caused by one accepted record, is dlog=301d_ =301 tokens for Claude Opus 4.6, 312 for Perplexity Sonar, and 317 for Qwen3-30B-A3B-Instruct-2507; the description text alone accounts for 231, 244, and 252 of them. The numbers of Analyst calls per accepted iteration, measured on the same runs, gives Nan=265N_an=265, 252, and 255 (Table 4). Equation (10) with these values gives Tanin=11.10T_an^in=11.10, 10.89, and 11.20 M input tokens per run, of which the fixed-prompt contribution NanFN_anF is 1.08, 1.03, and 1.04 M and the log contribution is 10.02, 9.86, and 10.16 M. The Analyst output contains the Chain-of-Thought reasoning and the technical specification of the method. Its mean length per call, measured on the same runs, is 4399 tokens for Claude Opus 4.6, 2315 for Perplexity Sonar, and 2291 for Qwen3-30B-A3B-Instruct-2507, so that Tanout=Nan⋅Tan,callout=1.17T_an^out=N_an· T_an,call^out=1.17, 0.58, and 0.58 M per run. The Coder receives its fixed prompt of 612 tokens together with the current specification and returns the Python implementation of the Decoding class; the mean length of accepted code is 4982 tokens for Claude Opus 4.6, 3749 for Perplexity Sonar, and 2513 for Qwen3-30B-A3B-Instruct-2507. If the code does not pass the functional test, the failing code and its trace stay in the dialogue and are re-read at every subsequent attempt of the same iteration, so the Coder input grows quadratically with the length of the debugging chain. The debug rate, measured from the runs, is 3.31 Coder calls per accepted iteration for Claude Opus 4.6, 1.70 for Perplexity Sonar, and 2.05 for Qwen3-30B-A3B-Instruct-2507, with mean failed-attempt lengths of 5991, 4641, and 2724 tokens, respectively. Reconstructing every Coder dialogue of the runs token by token gives Coder token counts per run of 12.98, 2.40, and 2.79 M for the input and 4.66, 1.73, and 1.32 M for the output. The Recorder receives its own prompt of 99 tokens together with the Coder prompt and the accepted code, and returns the description of length d, so that Trecin=2.61T_rec^in=2.61, 1.71, and 1.38 M and Trecout=Nd=0.058T_rec^out=Nd=0.058, 0.061, and 0.063 M per run. Summation of (9) over the three roles gives complete-run costs of $97.99, $16.08, and $5.85 for Claude Opus 4.6, Perplexity Sonar, and Qwen3-30B-A3B-Instruct-2507, with the Analyst contributing $84.66, $11.47, and $2.03, and costs up to kmink_ of $0.37, $0.21, and $0.26. Cost of Random Choice. Random Choice gives the lowest token usage among the compared systems. Each mini-call contains the catalog of method classes as input and one JSON decoder configuration as output; at about 800 input and 120 output tokens per call and 347 calls per run, the run requires 0.28 M input and 0.042 M output tokens. These values are 40 and 14 times smaller than the Analyst token counts of the Qwen agent with the detailed prompt (11.20 M and 0.58 M). In accuracy the two configurations are close: the per-run maxima of Random Choice lie between 89.48 and 89.93 %, those of the Qwen agent between 82.28 and 89.16 % (Table 3), and neither configuration reaches the 90 % threshold. When Qwen3-30B-A3B-Instruct-2507 is used free of charge, Crun=0C_run=0 for both configurations, and they differ only in token usage and in C∗C^*. At the managed rate for this model class, the Analyst cost differs by a factor of 29 ($0.07 for Random Choice against $2.03 for the agent framework with Qwen as the LLM Analyst), and with the Coder and the Recorder included by a factor of 84 ($0.07 against $5.85, Table 4). Sampling of the same catalog using Random Choice method is therefore two orders of magnitude cheaper than the Qwen agent framework at comparable accuracy. Both proprietary Analysts with the detailed prompt exceed it: their minimum run maxima are 90.93 % for Perplexity Sonar and 91.22 % for Claude Opus 4.6, against 89.93 % for Random Choice. The framework is thus preferable to sampling of the same catalog when the Analyst model is sufficiently strong; for a weak Analyst, the additional cost and the additional components of the framework are not compensated by the accuracy. 5.4 PINN results: single soliton The first PINN problem (GNLSE, α=0.1α=0.1, β=0β=0) has an analytical solution and therefore provides a transparent test of the framework: the agent can embed the explicit form of the one-soliton solution in the FI,FOF_I,F_O blocks, and the resulting metrics reflect the correctness of such a “hinted” representation rather than an independent search for a nontrivial technique. The results are summarized in Table 5. Because each agent run includes its own three-seed identity baseline, each reduction factor is calculated against the baseline from the same run. The best Claude Opus 4.6 hypothesis gives Relh=2.21⋅10−5Rel_h=2.21· 10^-5, a 24912491-fold reduction from 5.51⋅10−25.51· 10^-2. The best sonar-pro hypothesis reduces RelhRel_h from 4.69⋅10−24.69· 10^-2 to 1.97⋅10−41.97· 10^-4 (238×238×), while Qwen3-30B-A3B-Instruct-2507 reduces it from 4.85⋅10−24.85· 10^-2 to 8.86⋅10−38.86· 10^-3 (5.48×5.48×). Because the problem permits substitution of the analytical formula, these figures should not be interpreted as a measure of a complex search; the substantive test is deferred to the second problem. Table 5: PINN results for the single soliton (RelhRel_h metric; problem with an analytical solution). Each value is the mean over three internal seeds. Because baseline variation is non-negligible, the change for each analyst is calculated relative to the identity baseline from the same agent run. Model Baseline RelhRel_h Best RelhRel_h Change relative to baseline sonar-pro 4.69⋅10−24.69· 10^-2 1.97⋅10−41.97· 10^-4 238×238× better anthropic-claude-opus-4.6 5.51⋅10−25.51· 10^-2 2.21⋅10−52.21· 10^-5 2491×2491× better qwen-qwen3-30b-a3b-instruct-2507 4.85⋅10−24.85· 10^-2 8.86⋅10−38.86· 10^-3 5.5×5.5× better The best sonar-pro configuration is tw_nonlinear_pinn_core. Its FIF_I supplies a normalized traveling-wave coordinate, normalized time, and a local nonlinearity-density feature; FOF_O uses the analytical traveling-wave amplitude and phase and interprets the MLP outputs as bounded amplitude and phase corrections. The best Claude Opus 4.6 configuration, PCDS-LCC-EXP-NLOG-MIXCROSS, uses seven analytical features, including normalized and logarithmic envelopes and two amplitude–time cross-features, while FOF_O applies time-gated exponential amplitude and additive phase corrections with separate envelope weights. The best Qwen configuration, LogPhaseHybridNorm, normalizes the analytical soliton coordinate at the input and reconstructs the field from the analytical amplitude and carrier phase multiplied by learned complex corrections. Thus, all three best configurations exploit the known soliton structure in both the representation and the reconstruction. A detailed analysis of the hypotheses is provided in Appendix C. 5.5 PINN results: soliton with gaussian perturbation The second PINN problem (GNLSE, α=0.3α=0.3, β=0β=0; the initial condition is a soliton with a superimposed Gaussian perturbation) has no analytical solution, so the agent cannot embed a ready-made formula and must search for a substantive representation. The search is driven by ResbvpRes_bvp (4) and ReliRel_i; the reruns reported here do not include RelhRel_h. The human-designed PINN from prior work has Resbvp=9.11⋅10−5Res_bvp=9.11· 10^-5 and Reli=9.87⋅10−4Rel_i=9.87· 10^-4. The single-run results are summarized in Table 6. The shared identity baseline combines the lowest ReliRel_i and the lowest ResbvpRes_bvp observed across the analyst runs, and each reduction factor is computed relative to its ReliRel_i. The 325×325× improvement is obtained relative to the identity baseline from the same Claude Opus 4.6 run, whereas the 194×194× improvement is obtained relative to the shared cross-run reference baseline. The best Claude Opus 4.6 configuration by conservation error, SPM_XPM_CorrWidth_SplitScale_HighP, reduces ReliRel_i from 1.28⋅10−11.28· 10^-1 to 6.61⋅10−46.61· 10^-4 (194×194×) while giving Resbvp=3.95⋅10−5Res_bvp=3.95· 10^-5. A different Claude configuration, SplitAlpha_ExpPhase_HighP, gives the lowest residual, Resbvp=2.08⋅10−5Res_bvp=2.08· 10^-5, with Reli=1.40⋅10−3Rel_i=1.40· 10^-3. The best sonar-pro hypothesis by ReliRel_i improves this metric 20×20×, but its Resbvp=2.37⋅10−2Res_bvp=2.37· 10^-2 is worse than the shared baseline; the best Qwen hypothesis improves ReliRel_i only 1.36×1.36× and likewise does not improve the baseline residual. Thus, only the Claude hypotheses improve both reported criteria substantially. Table 6: Single-run PINN results for the soliton with a Gaussian perturbation (problem without an analytical solution). The shared identity-baseline row contains the best value of each metric across the analyst runs; these values may therefore come from different runs. Improvements in ReliRel_i are computed relative to this shared baseline. Model Result ReliRel_i ResbvpRes_bvp Change in ReliRel_i relative to baseline Human-designed PINN from prior work Zavertyaev et al. (2024) Reported result 9.87⋅10−49.87· 10^-4 9.11⋅10−59.11· 10^-5 130×130× better Identity baseline (best across runs) 1.28⋅10−11.28· 10^-1 1.04⋅10−31.04· 10^-3 — sonar-pro Best ReliRel_i 6.35⋅10−36.35· 10^-3 2.37⋅10−22.37· 10^-2 20×20× better anthropic-claude-opus-4.6 Best ReliRel_i 6.61⋅10−46.61· 10^-4 3.95⋅10−53.95· 10^-5 194×194× better anthropic-claude-opus-4.6 Best ResbvpRes_bvp 1.40⋅10−31.40· 10^-3 2.08⋅10−52.08· 10^-5 91×91× better qwen-qwen3-30b-a3b-instruct-2507 Best ReliRel_i 9.43⋅10−29.43· 10^-2 1.38⋅10−31.38· 10^-3 1.36×1.36× better Ten additional trainings from random initializations were performed for the identity baseline and for both Claude configurations. The configuration selected by ResbvpRes_bvp gives Reli=(1.42±0.59)⋅10−3Rel_i=(1.42± 0.59)· 10^-3, a 116×116× reduction from the repeated baseline mean, while the configuration selected by ReliRel_i gives (1.60±0.81)⋅10−3(1.60± 0.81)· 10^-3, a 103×103× reduction. The repeated-run export does not contain ResbvpRes_bvp, so these series confirm robustness of the invariant error but not of the joint criterion. Details are provided in Appendix D. 5.6 PINN results: Kuramoto–Sivashinsky In the traveling-wave problem with a Gaussian perturbation, the Claude Opus 4.6 agent generated the FI,FOF_I,F_O transformations while minimizing the boundary-value-problem residual ResbvpRes_bvp. The mass drift Δm(u) _m^(u) was recorded alongside it but did not drive the search. The best values are summarized in Table 7, and the two metrics are computed by different hypotheses. On the primary metric (ResbvpRes_bvp), the best architecture is MultiplicativeIC_DualAdditive_DualZoom, which reduces ResbvpRes_bvp from the baseline value of 1.41⋅10−31.41· 10^-3 to 1.511⋅10−41.511· 10^-4, a factor of 9.3. On the physical diagnostic, the best architecture is BoundedMultAdd_SechZoom_DualPhase, which reduces Δm(u) _m^(u) from 9.31⋅10−29.31· 10^-2 to 6.56⋅10−36.56· 10^-3, an approximately 1414-fold reduction, while also improving ResbvpRes_bvp by a factor of 8.3. Table 7: Best single-run PINN results for the Kuramoto–Sivashinsky traveling wave with a Gaussian perturbation (Claude Opus 4.6, one seed per hypothesis). ResbvpRes_bvp is the metric the analyst optimized; Δm(u) _m^(u) is an independent physical diagnostic. The two columns report the best value of each metric over the run; they are attained by different hypotheses, which are named in the last column. Configuration ResbvpRes_bvp Δm(u) _m^(u) Hypothesis Baseline (FI,FOF_I,F_O identity) 1.41⋅10−31.41· 10^-3 9.31⋅10−29.31· 10^-2 — Best ResbvpRes_bvp (9×9× better) 1.51⋅10−41.51· 10^-4 3.47⋅10−23.47· 10^-2 MultiplicativeIC_DualAdditive_DualZoom Best Δm(u) _m^(u) (14×14× better) 1.70⋅10−41.70· 10^-4 6.56⋅10−36.56· 10^-3 BoundedMultAdd_SechZoom_DualPhase Figure 5: Predicted solution of the Kuramoto–Sivashinsky problem at t=0,6,10t=0,6,10 for the baseline (top) and for the hypothesis with the lowest mass drift, BoundedMultAdd_SechZoom_DualPhase (bottom). No numerical reference is available for this problem, so only the PINN prediction is plotted; the initial condition is shown in the left panels for comparison. The horizontal axis is the coordinate x, the vertical axis is u. The baseline pulse loses roughly half of its amplitude by t=10t=10 and develops spurious low-amplitude oscillations across the domain, whereas the agent hypothesis keeps a localized traveling pulse whose amplitude stays close to the initial value of 1/21/2. Both winning hypotheses share a multiplicative–additive decomposition u^=u⁡(x,0)e−αt(1+tγmN1)+tW(x)(γ1N2+γ2N3) u=u(x,0)\,e^-α t (1+t\, _mN_1 )+t\,W(x) ( _1N_2+ _2N_3 ), in which N1,2,3N_1,2,3 are the three MLP outputs, the initial profile is embedded exactly at t=0t=0, the window W⁡(x)=sLsRW(x)=s_Ls_R built from the two normalized distances to the boundaries vanishes at x=−30x=-30 and x=40x=40, and the input block supplies features in the co-moving coordinate ξ=x−tξ=x-t — a rescaled tanh⁡(ξ/λ) (ξ/λ), the profile sech2(ξ/2)sech^2(ξ/2) of the unperturbed pulse, and a phase marker sin⁡(πt/6) (π t/6) of the collision. The hypothesis that wins on mass drift differs in one detail: the multiplicative channel is bounded, N1↦tanh⁡(N1)N_1 (N_1), so the correction cannot rescale the embedded initial profile without limit. This bound costs little in terms of the residual but keeps the integral quantity M=∫uxM= u\,dx an order of magnitude closer to its initial value, which illustrates that the two criteria are not interchangeable. 6 Ablation studies 6.1 Simple and detailed prompts This ablation compares the simple and detailed SNN prompt configurations defined in Section 4.3. The detailed configuration combines the method-class taxonomy, structured Chain-of-Thought, and coverage map described by the general framework in Section 3. The full Analyst and Coder prompts are provided in Supplementary Section F. The MNIST configurations are compared in Table 2. For Claude Opus 4.6, the detailed prompt increased the mean maximum accuracy over ten matched complete runs from 91.13±0.44%91.13± 0.44\,\% to 92.15±0.62%92.15± 0.62\,\%, and the best result from 91.65 % to 92.83 %. The difference between the prompt variants is confirmed by a two-sided exact Wilcoxon signed-rank test (W=5W=5, p=0.0195p=0.0195), with the detailed prompt yielding higher accuracy in eight of ten pairs; the mean paired difference is 1.02±0.901.02± 0.90 percentage points (Table 17), meaning that the effect is comparable in magnitude to the between-run variation. Ten complete runs were performed with each prompt variant for Perplexity Sonar and Qwen3-30B-A3B-Instruct-2507. With the simple prompt, none of the 250 decoders proposed per run outperformed the baseline, so the run maximum for both models exactly matched the baseline of 79.88 % in all ten cases. The detailed prompt changed the results of both models, but to different extents. The 90 % threshold was fixed before the run series as a reference point for MNIST; for the digits dataset, where the baseline is already 90.83 %, an analogous reference point was set at 93 % Section 4.5 and is not used in the tables in this section. For Perplexity Sonar, the individual-run maxima ranged from 90.93 % to 92.70 % (Table 15), with a mean of 91.76±0.56%91.76± 0.56\,\%, so the 90 % threshold was exceeded in all ten runs. For Qwen3-30B-A3B-Instruct-2507, the mean maximum was 85.88±2.42%85.88± 2.42\,\%, and the best run reached 89.16 % without exceeding the specified 90 % threshold. The difference between the prompt variants for both models is confirmed by a two-sided exact paired Wilcoxon signed-rank test. A run-by-run comparison between the simple and detailed prompts for Qwen3-30B-A3B-Instruct-2507 is provided in Table 13. With the simple prompt, accuracy remained at the baseline level of 79.88%79.88\,\% in all ten runs. The detailed prompt provided gains ranging from 2.402.40 to 9.289.28 percentage points in individual runs. The mean maximum accuracy was 85.88±2.42%85.88± 2.42\,\%, and the best result was 89.16%89.16\,\%. The high variability (std=2.42std=2.42 percentage points at n=10n=10) indicates instability of the search when this model is used. 6.2 Coverage-map ablation The ablation experiment excluding the coverage map was performed for all three analyst models (Tables 14, 16, and 18). For Claude Opus 4.6, removing the map reduces the mean maximum to 91.02±0.25%91.02± 0.25\,\%, which is 1.141.14 percentage points below the detailed prompt (92.15%92.15\,\%) and nearly matches the level of the simple prompt (91.13%91.13\,\%). All ten runs nevertheless retain accuracy above 90 %, indicating that the coverage map has little effect on the result for the strong model. For Perplexity Sonar, removing the map reduces the mean maximum to 82.34±0.70%82.34± 0.70\,\% (Δ=9.42 =9.42 percentage points), and no run reaches the threshold. For Qwen3-30B-A3B-Instruct-2507, the mean maximum falls to 80.18±0.45%80.18± 0.45\,\% (Δ=5.70 =5.70 percentage points), effectively returning to the baseline level of 79.88 %: the coverage map proves to be the key component that enabled this model to improve over the baseline argmax decoder (the class label is assigned by the neuron with the highest spike rate). 7 Discussion The experiments show that prompt detail can partially compensate for the limitations of the analyst model, although the degree of compensation depends on the LLM used. Claude Opus 4.6 achieves accuracy above 90 % with both prompt variants. With the simple prompt, Perplexity Sonar remains at the baseline level in all ten complete runs, but with the detailed prompt it consistently exceeds 90 % and reaches the level of Claude Opus 4.6. For Qwen3-30B-A3B-Instruct-2507, the detailed prompt also produces a substantial gain: the best result increases from 79.88 % to 89.16 %. At the same time, the higher variability of the maxima (85.88±2.42%85.88± 2.42\,\%) and the absence of any run with a maximum accuracy above 90 % show that the instruction structure does not fully compensate for the limitations of this model. When search cost is taken into account (Table 4), the advantage of Claude Opus 4.6 also persists in economic terms: despite the highest cost per run ($98 versus $16 for Sonar), it reaches its best decoder after only kmin=2k_ =2 to kmax=43k_ =43 evaluated candidates (Table 3), which makes the cost of reaching the best result comparable to that of Sonar ($0.37 versus $0.21), whereas Qwen3-30B-A3B-Instruct-2507, deployed on a rented GPU, costs less than both proprietary models ($5.90 per run at the managed rate for the same model class, with the Analyst alone contributing $2.00) but still fails to reach the target threshold over the full run. In every task, the agent finds a technique tied to the physics of the problem: class-conditional spike-rate statistics for the CCN, the soliton structure of the solution for the Schrödinger equation, and a phase-separated representation of the pulse–perturbation interaction for the Kuramoto–Sivashinsky equation. The second PINN problem (a soliton with a Gaussian perturbation) is particularly informative: its analytical solution is unknown, so the gain cannot be reduced to substituting a ready-made formula, and the best configuration (Claude Opus 4.6) improves the conservation of the integrals of motion by more than an order of magnitude. In the KS problem, the improvement appears in a different physical diagnostic, namely mass drift. All five experiments are nevertheless handled by the same framework without changes to its internal structure. 8 Limitations Despite the results reported above, the present study has several limitations that should be taken into account when interpreting the findings and when comparing our system with alternative approaches. • Restricted search space. The agent modifies a fixed catalogue of editable blocks (encoders, neuron models, connectivity schemes for SNNs; loss weights, activation functions and sampling strategies for PINNs). Structural transformations that fall outside this catalogue – for example, adding auxiliary networks or replacing the entire training paradigm – cannot be discovered by the current agent and would require an extended action space. • Flat experiment log without a hypothesis graph. The Recorder stores runs as a chronological JSON log; it does not maintain a directed acyclic graph of hypotheses and their dependencies. As a consequence, the Analyst re-derives relations between experiments from the log at every step, which limits the depth of multi-step reasoning that can be sustained over long sessions. • Bounded number of iterations. All reported runs are capped at 250250 iterations for SNN tasks and 5050 iterations for PINN tasks. These bounds were chosen to keep the wall-clock cost reproducible within the budgets reported in Table 4, but they also mean that the accuracy and relative-error values in Sections 5.2, 5.4, 5.5 and 5.6 are lower bounds on what the agent can eventually reach with a larger budget. • Single-seed evaluation in part of the regimes. For the most expensive configurations (Claude Opus 4.6 on the generalised nonlinear Schrödinger equation and on the Kuramoto–Sivashinsky benchmark) we report a single seed per configuration; the ten-run averages for best models are given only for the soliton with a Gaussian perturbation case. A full multi-seed sweep for every LLM/task pair is left for future work. • Dependence on proprietary LLM APIs. Two of the three evaluated models (Claude Opus 4.6 and Perplexity Sonar) are accessed through commercial APIs; their availability, pricing and behaviour may change between the time of the experiments and the time of reproduction. The open-weights Qwen3-30B-A3B run is provided as a partial mitigation, but its cost/quality trade-off is not directly comparable to the proprietary systems. • No head-to-head comparison with LLM-based NAS baselines. The reported baselines are a random-choice controller and the state-of-the-art results reported in the cited references. A direct comparison with recent LLM-driven neural-architecture-search systems on the same tasks and under the same iteration and cost budgets is desirable but was outside the scope of the present paper; we consider it the natural next step of this line of research. 9 Conclusion To date, automatic tuning of neural networks has been carried out mainly within the framework of an allocated single metric of the given competition: agent-based systems optimize the entire ML script, NAS approaches search over architectures according to accuracy, and PINN frameworks train networks to minimize the relative error. For tasks whose quality is defined by a multidimensional signal—by the conservation of a set of physical invariants or by additional structural requirements—frameworks of a different class are required. The results reported in this paper demonstrate the robustness and generalization ability of the proposed multi-agent framework, which consists of an analyst, a programmer, and a recorder, together with a sandboxed execution environment and an experiment log serving as external memory, on tasks of different nature without changes to its internal structure. The same unchanged framework was applied to five tasks with different types of information representation: accuracy increased from 79.88 % to 92.83 % on MNIST and from 90.83 % to 96.67 % on digits; in the single-soliton PINN problem, the relative amplitude error decreased by a factor of approximately 2491; in the more complex problem with a Gaussian perturbation and no known analytical solution, the conservation error of the integrals of motion decreased by a factor of approximately 325 in the selected run and by factors of 116 and 103, respectively, based on the mean errors over ten additional trainings of the residual-selected and invariant-selected configurations (Table 11); and in the Kuramoto–Sivashinsky problem, the boundary-value-problem residual decreased by a factor of 9.3 and mass drift by a factor of 14 (Tables 1 and 7). A separate result concerns the interaction between the analyst model and prompt structure. On MNIST, the detailed prompt enabled Perplexity Sonar to exceed 90 % accuracy in all ten complete runs and to reach a maximum of 92.70 %, close to the result of Claude Opus 4.6. For Qwen3-30B-A3B-Instruct-2507, maximum accuracy increased from the baseline level of 79.88 % to 89.16 %, but the target threshold was not reached. Prompt design can therefore substantially improve the performance of a less capable LLM, but only partially compensates for differences between models. Future work includes increasing the complexity of the PINN problem (the interaction of two solitons and other regimes without analytical solutions), jointly tuning several blocks (the decoder and plasticity rules for the SNN), and constructing models that account for the constraints of neuromorphic systems. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Funding This work was supported by ongoing institutional funding. No additional grants to carry out or direct this particular research were obtained. Acknowledgments The work was carried out according to the state assignment of the National Research Centre ”Kurchatov Institute” using the computing resources of the federal center for collective use ”Complex for modeling and data processing for megascience class installations” of the National Research Center ”Kurchatov Institute”, http://ckp.nrcki.ru/. Data availability The results of the SNN and PINN framework will be made publicly available in the https://github.com/sag111/llm-auto-tuning-results. Additionally, the full texts of the hypotheses proposed by the LLM Analyst and the code written by the LLM Coder for these hypotheses are presented there. The full texts of the prompts are available upon request to the authors of the article. Declaration of generative AI and AI-assisted technologies in the writing process During the preparation of this work, the authors used Perplexity AI to assist with language editing and improving the clarity of the manuscript. After using this tool, the authors reviewed and edited the content as needed and take full responsibility for the content of the manuscript. The large language models evaluated in this study were also used as research tools within the experimental methodology described in Sections 3 and 3.1. This methodological use is distinct from their use in manuscript preparation. Appendix A Baseline PINN model The model problem is a boundary-value problem for the generalized nonlinear Schrödinger equation (GNLSE) with cubic–quintic nonlinearity: iqt+qxx+|q|2q(1−α|q|2+β|q|4)=0,q(x, 0)=q~(x, 0),q(xmin,t)=0,q(xmax,t)=0,(x,t)∈[xmin,xmax]×[tmin,tmax]. casesi\,q_t+q_x+|q|^2q\, (1-α|q|^2+β|q|^4 )=0,\\[4.0pt] q(x,\,0)= q(x,\,0),\\[2.0pt] q(x_ ,\,t)=0, q(x_ ,\,t)=0,\\[2.0pt] (x,t)∈[x_ ,x_ ]×[t_ ,t_ ]. cases (11) In all configurations, β=0β=0; the quintic term is retained in the equation for generality. Two configurations are considered. In the single-soliton problem, α=0.1α=0.1 and the domain is (x,t)∈[−25,25]×[0,3](x,t)∈[-25,25]×[0,3]; the initial profile is specified by the analytical soliton solution: q~(x,0)=μz(12z+1)2−αμ3z2⋅ei⁡(kx+θ0),z=eμ(x−x0), q(x,0)= μ z ( 12z+1 )^2- αμ3z^2· e^i(kx+ _0), z=e μ(x-x_0), (12) where k=1k=1, w=0.88w=0.88, x0=0x_0=0, θ0=0 _0=0 and μ=4(k2−w)μ=4(k^2-w). The complex field is decomposed into two real components: q⁡(x,t)=u⁡(x,t)+iv(x,t)q(x,t)=u(x,t)+i\,v(x,t). The analytical soliton is used as the reference when calculating the error. In the second configuration (a soliton with a Gaussian perturbation), α=0.3α=0.3, the domain is extended to (x,t)∈[−50,50]×[0,25](x,t)∈[-50,50]×[0,25], and the initial condition is the sum of the soliton profile (12) (with x0=−30x_0=-30) and a real-valued Gaussian packet: q(x,0)=q2(x,0)+0.2e−x2/50,q(x,0)=q_2(x,0)+0.2\,e^-x^2/50, (13) where q2q_2 coincides with (12) for the specified parameters, equivalently Agexp[−(x−Bg)2/(2Cg2)]A_g [-(x-B_g)^2/(2C_g^2)] with Ag=0.2A_g=0.2, Bg=0B_g=0, and Cg=5C_g=5 (the printed Cg2C_g^2 in Zavertyaev et al. (2024) is a typographical error, must be 2Cg22C_g^2). The analytical solution is unknown. A PINN field from Zavertyaev et al. (2024) is retained as an auxiliary amplitude reference. On [−50,50]×[0,25][-50,50]×[0,25] that field has Resbvp=9.11⋅10−5Res_bvp=9.11· 10^-5 and Reli=9.87⋅10−4Rel_i=9.87· 10^-4. The baseline network in this configuration is deeper, with 4 hidden layers of 64 neurons each, and collocation-point sampling is adaptive (version v3 from Moloshnikov et al. (2025b), regenerated every 1000 iterations with coefficient λ2=0.01 _2=0.01 in PrlP_rl in Eq. 6 of Moloshnikov et al. (2025b)). The primary metrics of the search are ResbvpRes_bvp (4) and the relative conservation error of the integrals of motion Reli=Reli1+Reli2Rel_i=Rel_i_1+Rel_i_2; RelhRel_h has lowest priority. The baseline is a fully connected network MLP:ℝ2→ℝ2MLP ^2 ^2, (x,t)↦(u,v)(x,t) (u,v), with two hidden layers of 32 neurons each and sin activation. In the modified architecture, two tunable blocks, FI(x,t)F_I(x,t) and FO(u,v,x,t)F_O(u,v,x,t), are attached to the MLP; both blocks are identity mappings in the baseline. The complete set of hyperparameters for both configurations is summarized in Table 8. The loss functional consists of three terms with equal weights: ℒ=ℒpde+ℒic+ℒbc,L=L_pde+L_ic+L_bc, (14) where ℒpdeL_pde is the mean squared residual of Eq. (11) at the collocation points, ℒicL_ic is the initial-condition error, and ℒbcL_bc is the boundary-condition error. The derivatives qt,qxxq_t,q_x are obtained by TensorFlow automatic differentiation. For the single-soliton problem the metric is: Relh=‖htrue−hpred‖2‖htrue‖2,h=u2+v2.Rel_h= \|h_true-h_pred\|_2\|h_true\|_2, h= u^2+v^2. (15) Table 8: Hyperparameters of the baseline PINN for the two GNLSE configurations. The Kuramoto–Sivashinsky configuration is given separately in Table 9. Parameter Single soliton Soliton + Gaussian Number of hidden layers 2 4 Neurons per layer 32 64 Activation (hidden layers) sin sin Activation (output layer) linear linear Optimizer Adam Adam Learning rate 10−310^-3 10−310^-3 Number of training iterations 10 000 30 000 Collocation points (NfN_f) 1 000 30 000 Initial-condition points (N0N_0) 100 3 000 Boundary-condition points (NbN_b) 50 300 Collocation-point sampling fixed random adaptive (v3) Weights wpde,wic,wbcw_pde,\,w_ic,\,w_bc 1,1,11,1,1 1,1,11,1,1 A.1 Kuramoto–Sivashinsky baseline hyperparameters The Kuramoto–Sivashinsky problem uses its own baseline configuration, summarized in Table 9. It differs from the GNLSE configurations in three respects. First, the network output is a single real scalar. Second, the fourth derivative in the residual makes the raw loss stiff, so the residual is divided by R0=10R_0=10 before squaring, the learning rate is lowered to 5⋅10−55· 10^-5, and the gradient norm is clipped at 1. Without these three measures a substantial fraction of initializations diverges to non-finite loss within the first few hundred iterations. Third, the boundary term penalizes both u and uxxu_x at the two boundaries, so the loss reads ℒ=⟨(ℛKS/R0)2⟩+⟨(upred(x,0)−u⁡(x,0))2⟩+⟨u2+uxx2⟩∂.L= (R_KS/R_0)^2 + (u_pred(x,0)-u(x,0))^2 + u^2+u_x^2 _∂. Table 9: Hyperparameters of the baseline PINN for the Kuramoto–Sivashinsky traveling-wave problem with a Gaussian perturbation. Parameter Value Number of hidden layers 4 Neurons per layer 32 Network output u (one real scalar) Activation (hidden layers) sin Activation (output layer) linear Optimizer Adam, clipnorm=1clipnorm=1 Learning rate 5⋅10−55· 10^-5 Residual scale R0R_0 10 Number of training iterations 50 000 Collocation points (NfN_f) 20 000 Initial-condition points (N0N_0) 1 000 Boundary-condition points (NbN_b) 200 Collocation-point sampling adaptive (v3), every 1 000 iterations Adaptive-sampling parameters λ2=0.01 _2=0.01, oversampling ×5× 5 Evaluation grid (nx×ntn_x× n_t) 1000×101000× 10 Weights wpde,wic,wbcw_pde,\,w_ic,\,w_bc 1,1,11,1,1 A.2 Kuramoto–Sivashinsky field shift and PINN residual For completeness, we derive the zero-background form of the generalized Kuramoto–Sivashinsky equation used by the PINN. We use the normalization vt+vxx+σvxxx+vxxxx+γvvx+χ(v2)xx=0,v_t+v_x+σ v_x+v_x+γ v_x+χ(v^2)_x=0, (16) where γ=2βγ=2β, so that the nonlinear convection term integrates to βv2β v^2 under a traveling-wave reduction. For β=3σβ=3σ, substitution of σ=−1σ=-1, γ=2β=6σ=−6γ=2β=6σ=-6, and χ=3χ=3 gives vt+vxx−vxxx+vxxxx−6vvx+3(v2)xx=0.v_t+v_x-v_x+v_x-6v_x+3(v^2)_x=0. (17) The exact wave from Kudryashov (2026) is expressed on the background −1/3-1/3. To make the learned field decay to zero away from the pulse, we set v=u−13,u=v+13.v=u- 13, u=v+ 13. (18) Because the shift is constant, all differentiated linear terms are unchanged: vt=ut,vx=ux,vxx=uxx,vxxx=uxxx,vxxxx=uxxxx.v_t=u_t, v_x=u_x, v_x=u_x, v_x=u_x, v_x=u_x. The nonlinear convection term transforms as −6vvx=−6(u−13)ux=−6uux+2ux,-6v_x=-6 (u- 13 )u_x=-6u_x+2u_x, (19) whereas the nonlinear diffusion term becomes v2 v^2 =(u−13)2=u2−23u+19, = (u- 13 )^2=u^2- 23u+ 19, (20) 3(v2)xx 3(v^2)_x =3(u2)xx−2uxx. =3(u^2)_x-2u_x. (21) Inserting these expressions into Eq. (17) yields ut+uxx−uxxx+uxxxx−6uux+2ux+3(u2)xx−2uxx=0.u_t+u_x-u_x+u_x-6u_x+2u_x+3(u^2)_x-2u_x=0. Combining the two uxxu_x terms gives the residual used in the implementation: ℛKS=ut−uxx−uxxx+uxxxx−6uux+2ux+3(u2)xx.R_KS=u_t-u_x-u_x+u_x-6u_x+2u_x+3(u^2)_x. (22) Thus, the negative coefficient of uxxu_x and the additional linear-advection term 2ux2u_x are both consequences of removing the constant background; they are not independent changes to the source equation. For the unperturbed case, the source-field solution and its zero-background counterpart are vex(x,t)=−13+12cosh2⁡((x−t−x0)/2),uex(x,t)=12cosh2⁡((x−t−x0)/2).v_ex(x,t)=- 13+ 12 ^2((x-t-x_0)/2), u_ex(x,t)= 12 ^2((x-t-x_0)/2). Appendix B Baseline SNN model The Correlation Classwise Network (CCN) Sboev et al. (2025); Kunitsyn et al. (2025) was selected as the baseline model. Its architecture, input encoding, supervisory-signal delivery, and plasticity rule are adopted without modification; the agent’s modifications are restricted to the output-spike-activity decoder. Input expansion using Gaussian receptive fields (GRF) and the bagging ensemble Rybka et al. (2024) are not used because preliminary experiments on MNIST showed no appreciable gain. The approach is applied to two datasets: MNIST (28×2828× 28, N=784N=784 spike generators) and the scikit-learn digits Pedregosa et al. (2011) (8×88× 8, N=64N=64). Both have ten classes, pixel values are normalized by the L1L_1 norm of the input vector, and the number of output LIF neurons is K=10K=10. The network has two layers: the input layer contains N spike generators, and the output layer contains K Leaky-Integrate-and-Fire (LIF) neurons; the connections are fully connected and trainable, with no lateral connections. The input vector V=(v1,…,vN)V=(v_1,…,v_N) is converted into an array of correlated spike sequences Si\S_i\: a reference sequence S0S_0 of length NtN_t is constructed with constant spacing NrefN_ref; component viv_i is interpreted as the target Pearson correlation coefficient between SiS_i and S0S_0, and SiS_i is obtained by randomly removing some spikes from S0S_0. Supervision is implemented using two currents applied to the output LIF neurons: for a sample of class k, a positive current Iexc=107I_exc=10^7 pA, synchronized with the S0S_0 spikes with a shift τs _s, is applied to the corresponding neuron; the other neurons receive an inhibitory current Iinh=−103I_inh=-10^3 pA. The weights are adjusted using a local multiplicative STDP rule with a bounding symmetric nearest-neighbor scheme: Δw(Δt)=A+⋅wμ+⋅e−Δt/τ+,Δt>0,A−⋅wμ−⋅eΔt/τ−,Δt<0, w( t)= casesA^+· w _+· e^- t/ _+,& t>0,\\ A^-· w _-· e t/ _-,& t<0, cases where Δt=tpost−tpre t=t_post-t_pre, w is the current weight, A±A^± are the potentiation and depression amplitudes, τ± _± are the window time constants, and μ± _± are the exponents of the weight dependence. The values of all network and STDP-rule parameters are provided in Table 10. Table 10: Hyperparameters of the baseline SNN and its training, common to all configurations. Parameter Value VrestV_rest (resting potential) −70-70 mV VresetV_reset (post-spike potential) −70-70 mV VthV_th (threshold potential) −56-56 mV τm _m (membrane time constant) 10310^3 ms τs _s (teacher-spike shift) 0.10.1 ms NrefN_ref (spacing between spikes in S0S_0) 3 NtN_t (number of time steps in S0S_0) 3000 TsimT_sim (simulation time, dt=0.1dt=0.1 ms) 300 ms A+A^+ (potentiation amplitude) 1.2⋅10−31.2· 10^-3 A−A^- (depression amplitude) 10−310^-3 τ+ _+ (potentiation window) 200200 ms τ− _- (depression window) 160160 ms μ+ _+ (weight dependence, potentiation) 0.960.96 μ− _- (weight dependence, depression) 0.260.26 B.1 SNN experimental protocol MNIST data split. The standard MNIST test set, containing 10 000 samples, was not used during the agent search. The original training set of 60 000 samples was stratified into training and validation subsets at a 80:20 ratio. The parameters of each proposed decoder were estimated on the training subset, whereas the classification accuracy passed to the agent-search log was calculated only on the validation subset. After the specified search budget was exhausted, the decoder with the highest validation accuracy was selected. The selected decoder was retrained on the complete original MNIST training set and evaluated once on the standard test set. Test samples, their spike representations, and their labels were not provided to the agents and were not used to select the decoder structure or tune its hyperparameters. The same stratified split of the training set was used in all compared configurations. digits data split. Before the agent search began, the digits dataset, containing 1797 samples, was stratified into a development set of 1437 samples and a held-out test set of 360 samples. The test set was not used to calculate metrics passed to the agents, select the decoder structure, or tune its hyperparameters. Stratified five-fold cross-validation was performed on the development set. The five-fold split was fixed once before the experiments and remained unchanged for all decoders, runs, analyst models, prompt variants, and baseline methods. At each iteration, the proposed decoder was trained five times on four folds and evaluated on the remaining fold. The mean classification accuracy was recorded in the experiment log: CVAcc⁡(d)=15∑k=15Acc⁡(d(−k),(k)),CVAcc(d)= 15 _k=1^5Acc (d^(-k),D^(k) ), (23) where d(−k)d^(-k) denotes the decoder trained on four folds, and (k)D^(k) is the fifth fold excluded from training. The analyst used CVAcc⁡(d)CVAcc(d) when formulating the next hypothesis. After the specified search budget was exhausted, the decoder with the highest CVAcc⁡(d)CVAcc(d) was selected. It was then retrained on all 1437 samples in the development set and evaluated once on the 360 samples in the held-out test set. The resulting test accuracy was not recorded in the agent-search log. Appendix C Analysis of the best PINN hypotheses The single-soliton comparison includes Perplexity Sonar (sonar-pro), Claude Opus 4.6, and Qwen3-30B-A3B-Instruct-2507 (Table 5). The agent controls the input transformation FIF_I and output transformation FOF_O, while the MLP, activation function, collocation points, loss functional, and training schedule remain fixed. The best overall result was obtained by the Claude Opus 4.6 configuration PCDS-LCC-EXP-NLOG-MIXCROSS, which gives Relh=2.21⋅10−5Rel_h=2.21· 10^-5, a 24912491-fold reduction from its identity baseline of 5.51⋅10−25.51· 10^-2. It uses seven analytical features, including normalized and logarithmic envelopes and two amplitude–time cross-features, while FOF_O applies time-gated exponential amplitude and additive phase corrections with separate envelope weights. The best sonar-pro configuration, tw_nonlinear_pinn_core, gives Relh=1.97⋅10−4Rel_h=1.97· 10^-4, a 238238-fold reduction from its identity baseline of 4.69⋅10−24.69· 10^-2. Its FIF_I supplies a normalized traveling-wave coordinate, normalized time, and a local nonlinearity-density feature; FOF_O uses the analytical traveling-wave amplitude and phase and interprets the MLP outputs as bounded amplitude and phase corrections. The best Qwen configuration, LogPhaseHybridNorm, gives Relh=8.86⋅10−3Rel_h=8.86· 10^-3, a 5.485.48-fold reduction from its identity baseline of 4.85⋅10−24.85· 10^-2. It normalizes the analytical soliton coordinate at the input and reconstructs the field from the analytical amplitude and carrier phase multiplied by learned complex corrections. All three best configurations exploit the available analytical solution. Their performance should therefore be interpreted as evidence that the framework can discover a useful physics-informed representation, rather than as evidence of performance on problems without a known solution. Appendix D Repeated PINN runs and statistical analysis To assess whether the observed improvements are robust to random initialization, we repeated the training procedure for the GNLSE problem with a Gaussian perturbation. We compared the baseline model with two agent-designed configurations: one selected for its low equation residual and the other for its accurate preservation of the physical invariant. Each configuration was successfully trained ten times. Because the random initializations were not matched across configurations, the resulting samples were treated as independent rather than paired. Table 11 summarizes the results. Table 11: Repeated runs for the GNLSE soliton problem with a Gaussian perturbation. The table reports the mean and sample standard deviation, median, range, and coefficient of variation (CV) of ReliRel_i. Configuration Mean ± std Median Min–max CV Reduction factor in ReliRel_i relative to baseline mean Baseline (n=10n=10) (1.65±0.18)⋅10−1(1.65± 0.18)· 10^-1 1.66⋅10−11.66· 10^-1 1.29⋅10−11.29· 10^-1–1.94⋅10−11.94· 10^-1 11.1%11.1\,\% 1×1× Residual-selected configuration (n=10n=10) (1.42±0.59)⋅10−3(1.42± 0.59)· 10^-3 1.50⋅10−31.50· 10^-3 5.40⋅10−45.40· 10^-4–2.41⋅10−32.41· 10^-3 41.8%41.8\,\% 116.1×116.1× Invariant-selected configuration (n=10n=10) (1.60±0.81)⋅10−3(1.60± 0.81)· 10^-3 1.59⋅10−31.59· 10^-3 7.57⋅10−47.57· 10^-4–3.31⋅10−33.31· 10^-3 50.4%50.4\,\% 103.3×103.3× Both agent-designed configurations reduce the invariant error by more than two orders of magnitude on average relative to the baseline. The residual-selected configuration has a slightly lower mean error and coefficient of variation, although the two distributions overlap substantially. In the original search, the invariant-selected configuration achieved an error of 6.61⋅10−46.61· 10^-4, lower than in any of its ten subsequent retrainings. By contrast, the original error of the residual-selected configuration, 1.40⋅10−31.40· 10^-3, is close to its repeated-run mean. This contrast indicates a selection effect: the best result found during a search can be more favorable than the performance typically obtained when the same configuration is retrained. Appendix E Additional classification results Method-class coverage map. The aggregate coverage map of decoder method classes from the complete 250-iteration log for the ⟨ Sonar, detailed prompt⟩ configuration on MNIST is presented in Table 12. Each valid decoder is assigned to one principal class according to the substance of its decision rule, while decoders that explicitly aggregate several predictions are grouped in the “Ensembles” row. The map shows that the analyst spent most of the budget on hybrid ensemble schemes and discriminant analysis, whereas the correlation, margin-based, and argmax classes remained uncovered in this run. Table 12: Coverage map of decoder method classes for one 250-iteration run (Perplexity Sonar, detailed prompt, MNIST). Method class Methods Best acc, % Status Ensembles 154 92.70 explored LDA / QDA 32 87.90 explored Softmax / calibration 13 82.32 explored k-nearest neighbors 11 88.56 explored Naive Bayes 11 83.73 explored Mahalanobis 9 86.05 limited data Rank-based 8 81.68 limited data Prototypes / centroids 5 80.43 limited data KDE / density 4 85.75 limited data Threshold methods / rules 3 73.93 limited data Correlation / cosine 0 — not explored Margin-based 0 — not explored Best decoders and coverage-map ablation. In the best configuration (Claude Opus 4.6, detailed prompt), the 92.83 % result was obtained by a hybrid method: class-averaged representations serve as centroids, and the decision is based on Mahalanobis distance with class-specific covariance. The pattern is the same for Perplexity Sonar with the detailed prompt: the best decoders belong to LDA/QDA and prototype method classes. These classes use the entire output spike-rate vector, whereas the baseline argmax uses only one of its components. The ablation confirms the role of the coverage map: when it is removed (while the catalog of method classes and Chain-of-Thought are retained), Perplexity Sonar does not reach 90 % during a run, and the best accuracy falls from 92.70 % to 83.18 %; the analyst converges on a narrow set of similar ensemble variations. The results for Qwen3-30B-A3B-Instruct-2507 are provided in Tables 13 and 14, those for Perplexity Sonar in Tables 15 and 16, and those for Claude Opus 4.6 in Tables 17 and 18. Table 13: Maximum accuracy in ten matched pairs of complete agent-system runs with the simple and detailed prompts for Qwen3-30B-A3B-Instruct-2507 on MNIST. Within each pair, the same trained SNN instance and identical output spike data were used. The difference Δ is calculated as the accuracy for the detailed prompt minus the accuracy for the simple prompt. Pair no. Simple, % Detailed, % Δ , p 1 79.88 82.28 2.40 2 79.88 83.35 3.47 3 79.88 88.57 8.69 4 79.88 86.48 6.60 5 79.88 89.16 9.28 6 79.88 84.10 4.22 7 79.88 87.20 7.32 8 79.88 85.50 5.62 9 79.88 83.90 4.02 10 79.88 88.30 8.42 Mean 79.88 85.88 6.00 Std. dev. 0.00 2.42 2.42 Table 14: Comparison of maximum accuracy in ten matched Qwen3-30B-A3B-Instruct-2507 runs on MNIST for the detailed prompt with and without the coverage map. Within each pair, the same SNN instance and identical output spike data were used. The difference Δ is calculated as the accuracy with the coverage map minus the accuracy without the coverage map. Pair no. With map, % Without map, % Δ , p 1 82.28 79.88 2.40 2 83.35 80.42 2.93 3 88.57 79.88 8.69 4 86.48 81.15 5.33 5 89.16 79.88 9.28 6 84.10 80.73 3.37 7 87.20 79.88 7.32 8 85.50 79.88 5.62 9 83.90 80.25 3.65 10 88.30 79.88 8.42 Mean 85.88 80.18 5.70 Std. dev. 2.42 0.45 2.58 Table 15: Maximum accuracy in ten matched pairs of complete agent-system runs with the simple and detailed prompts for Perplexity Sonar on MNIST. Within each pair, the same trained SNN instance and identical output spike data were used. The difference Δ is calculated as the accuracy for the detailed prompt minus the accuracy for the simple prompt. Pair no. Simple, % Detailed, % Δ , p 1 79.88 92.70 12.82 2 79.88 92.57 12.69 3 79.88 91.95 12.07 4 79.88 91.82 11.94 5 79.88 91.77 11.89 6 79.88 91.74 11.86 7 79.88 91.64 11.76 8 79.88 91.45 11.57 9 79.88 91.08 11.20 10 79.88 90.93 11.05 Mean 79.88 91.76 11.89 Std. dev. 0.00 0.56 0.56 Table 16: Comparison of maximum accuracy in ten matched Perplexity Sonar runs on MNIST for the detailed prompt with and without the coverage map. Within each pair, the same SNN instance and identical output spike data were used. The difference Δ is calculated as the accuracy with the coverage map minus the accuracy without the coverage map. Pair no. With map, % Without map, % Δ , p 1 92.70 83.18 9.52 2 92.57 82.92 9.65 3 91.95 82.77 9.18 4 91.82 82.66 9.16 5 91.77 82.65 9.12 6 91.74 82.47 9.27 7 91.64 82.15 9.49 8 91.45 82.04 9.41 9 91.08 81.87 9.21 10 90.93 80.71 10.22 Mean 91.76 82.34 9.42 Std. dev. 0.56 0.70 0.33 Table 17: Maximum accuracy in ten matched pairs of complete agent-system runs with the simple and detailed prompts for Claude Opus 4.6 on MNIST. Within each pair, the same trained SNN instance and identical output spike data were used. The difference Δ is calculated as the accuracy for the detailed prompt minus the accuracy for the simple prompt. Pair no. Simple, % Detailed, % Δ , p 1 90.20 92.83 2.63 2 91.15 91.88 0.73 3 91.65 91.58 -0.07 4 91.60 92.49 0.89 5 91.33 91.22 -0.11 6 90.71 92.35 1.64 7 91.42 92.69 1.27 8 91.17 92.57 1.40 9 90.88 92.66 1.78 10 91.22 91.24 0.02 Mean 91.13 92.15 1.02 Std. dev. 0.44 0.62 0.90 Table 18: Comparison of maximum accuracy in ten matched Claude Opus 4.6 runs on MNIST for the detailed prompt with and without the coverage map. Within each pair, the same SNN instance and identical output spike data were used. The difference Δ is calculated as the accuracy with the coverage map minus the accuracy without the coverage map. Pair no. With map, % Without map, % Δ , p 1 91.22 90.66 0.56 2 91.24 90.73 0.51 3 91.58 90.89 0.69 4 91.88 90.91 0.97 5 92.35 90.95 1.40 6 92.49 91.01 1.48 7 92.57 91.04 1.53 8 92.66 91.20 1.46 9 92.69 91.32 1.37 10 92.83 91.45 1.38 Mean 92.15 91.02 1.14 Std. dev. 0.62 0.25 0.41 When the simple and detailed prompts were compared, the mean maximum accuracy for Qwen3-30B-A3B-Instruct-2507 increased from 79.88±0.00%79.88± 0.00\,\% to 85.88±2.42%85.88± 2.42\,\%. A positive difference was also observed in all ten pairs and averaged 6.006.00 percentage points. For Claude Opus 4.6, the mean maximum increased from 91.13±0.44%91.13± 0.44\,\% to 92.15±0.62%92.15± 0.62\,\%. The mean paired difference was 1.021.02 percentage points: a positive effect was observed in 8 of 10 pairs, whereas accuracy decreased by 0.070.07 and 0.110.11 percentage points in two pairs. For Perplexity Sonar, switching to the detailed prompt increased the mean maximum accuracy from 79.88±0.00%79.88± 0.00\,\% to 91.76±0.56%91.76± 0.56\,\%. A positive difference was observed in all ten pairs, with a mean of 11.8911.89 percentage points. The coverage-map ablation showed an additional effect of this component of the detailed prompt. For Qwen3-30B-A3B-Instruct-2507, mean accuracy increased from 80.18±0.45%80.18± 0.45\,\% to 85.88±2.42%85.88± 2.42\,\%, i.e., by 5.705.70 percentage points. In each of the three ablation comparisons, a positive difference was observed in all ten pairs. For Perplexity Sonar, the corresponding values were 82.34±0.70%82.34± 0.70\,\% without the map and 91.76±0.56%91.76± 0.56\,\% with the map, with a mean increase of 9.429.42 percentage points. For Claude Opus 4.6, the presence of the map increased the mean maximum accuracy from 91.02±0.25%91.02± 0.25\,\% to 92.15±0.62%92.15± 0.62\,\%; the mean paired increase was 1.141.14 percentage points. For each model, the statistical significance of differences between the two configurations was assessed using a two-sided exact Wilcoxon signed-rank test for paired samples. The null hypothesis stated that neither compared configuration had a systematic advantage: the paired accuracy differences were centered at zero. The alternative hypothesis assumed a systematic difference between the configurations. The test was applied under the assumption of a symmetric distribution of paired differences. For the two-sided test, the W statistic was defined as the smaller of the rank sums for positive and negative nonzero differences. The significance level was set to α=0.05α=0.05. For the comparison of the simple and detailed prompts for Claude Opus 4.6, W=5W=5 and p=0.0195p=0.0195 were obtained. For Perplexity Sonar and Qwen3-30B-A3B-Instruct-2507, all ten paired differences were positive. Therefore, in both cases the two-sided exact Wilcoxon signed-rank test gave the same result: W=0W=0, p=0.001953p=0.001953. The null hypothesis was rejected for all three models at the α=0.05α=0.05 significance level. Two research questions were tested for the three models: the effect of switching from the simple to the detailed prompt and the effect of including the coverage map in the detailed prompt. Thus, six statistical comparisons were performed, and their results were used to formulate the overall conclusion about the effectiveness of prompt structure. Testing each hypothesis separately at a significance level of α=0.05α=0.05 does not preserve this level for the full set of comparisons. Specifically, if the six tests were independent and all six null hypotheses were true, the probability of at least one false-positive result would be 1−(1−α)6=1−0.956=0.2649,1-(1-α)^6=1-0.95^6=0.2649, i.e., approximately 26.49%26.49\,\%. To ensure control of the family-wise probability of at least one type I error at 0.050.05, the Holm–Bonferroni correction was applied. This procedure controls the family-wise error even when the statistical tests are dependent. The original p values were ordered in ascending order and sequentially compared with the thresholds α6,α5,…,α. α6,\; α5,\;…,\;α. For the five comparisons with an original p=0.001953p=0.001953, the adjusted value was pHolm=0.0117p_Holm=0.0117. For the comparison of the simple and detailed prompts for Claude Opus 4.6, with an original p=0.0195p=0.0195, pHolm=0.0195p_Holm=0.0195 was obtained. All six adjusted p values were below 0.050.05, so the corresponding null hypotheses were also rejected after correction for multiple comparisons. When the detailed prompt with and without the coverage map was compared, all ten paired differences were positive for each of the three models. The same values, W=0W=0, p=0.001953p=0.001953. After applying the Holm–Bonferroni correction across all six comparisons, the adjusted value was pHolm=0.0117p_Holm=0.0117. Consequently, removing the coverage map led to a statistically significant decrease in the maximum accuracy achieved for each model. After correction for multiple comparisons, the results of all six paired tests remained statistically significant. When switching from the simple to the detailed prompt, the mean accuracy increase was 1.021.02 percentage points for Claude Opus 4.6, 11.8911.89 percentage points for Perplexity Sonar, and 6.006.00 percentage points for Qwen3-30B-A3B-Instruct-2507. For Claude Opus 4.6, the two-sided exact Wilcoxon signed-rank test gave W=5W=5, an original p=0.0195p=0.0195, and pHolm=0.0195p_Holm=0.0195. For Perplexity Sonar and Qwen3-30B-A3B-Instruct-2507, positive differences were obtained in all ten pairs; for each model, W=0W=0, the original p=0.001953p=0.001953, and pHolm=0.0117p_Holm=0.0117. Results on the digits dataset. The same framework was applied to digits 8×88× 8; the run design was the same as for MNIST. The baseline argmax decoder gives an accuracy of 90.83 %, whereas the best ⟨ Opus 4.6, detailed prompt⟩ configuration increases it to 96.67 % (Table 19). The best decoder uses a nonparametric generative classifier based on KDE with a Gaussian kernel in a PCA subspace, and the label is determined using the Bayes rule c^=arg⁡maxc[log⁡p⁡(f∣c)+log⁡p⁡(c)] c= _c[ p(f c)+ p(c)]; hyperparameters were selected using 5-fold cross-validation. Table 19: Comparison of analyst models and prompt variants on the digits 8×88× 8 dataset. For each configuration, the table reports the number of complete agent runs, the mean and standard deviation of the maximum accuracy achieved within an individual run, the highest accuracy across all runs, and the number of runs in which the 93 % threshold was reached. Ten matched runs of 250 iterations were performed for every configuration; the baseline (argmax decoder) gives an accuracy of 90.83 %. Model Prompt n Mean maximum accuracy ± std, % Best accuracy, % Runs with accuracy ≥93%≥ 93\,\% Claude Opus 4.6 simple 10 95.07±0.5295.07± 0.52 95.56 10 out of 10 Claude Opus 4.6 detailed 10 96.12±0.6696.12± 0.66 96.67 10 out of 10 Perplexity Sonar simple 10 93.44±0.5793.44± 0.57 93.89 9 out of 10 Perplexity Sonar detailed 10 93.31±0.9393.31± 0.93 94.17 8 out of 10 Qwen3-30B-A3B-Instruct-2507 simple 10 91.63±0.7491.63± 0.74 92.78 0 out of 10 Qwen3-30B-A3B-Instruct-2507 detailed 10 93.66±0.7993.66± 0.79 94.17 7 out of 10 E.1 Decoder search dynamics In addition to the final accuracy values reported in Tables 2 and 19, the search trajectory itself is of interest: how quickly the agent system finds a decoder that outperforms the baseline and reaches a plateau. Figure 6 shows the maximum accuracy among all decoders proposed up to the current iteration (the cumulative maximum) for each run of the agent system. The MNIST trajectories are collected in Fig. 6, and those for the digits dataset in Fig. 7. (a) Simple prompt, Qwen3-30B (b) Simple prompt, Perplexity Sonar (c) Simple prompt, Claude Opus 4.6 (d) Detailed prompt, Qwen3-30B (e) Detailed prompt, Perplexity Sonar (f) Detailed prompt, Claude Opus 4.6 Figure 6: Search dynamics for the decoding module on the MNIST dataset. The blue line shows the accuracy of the solution obtained at each iteration, the red dashed line shows the best result achieved, and the gray horizontal line shows the baseline accuracy of 79.88 %. The columns correspond to analyst models, and the upper and lower rows correspond to the simple and detailed prompts, respectively. (a) Simple prompt, Qwen3-30B (b) Simple prompt, Perplexity Sonar (c) Simple prompt, Claude Opus 4.6 (d) Detailed prompt, Qwen3-30B (e) Detailed prompt, Perplexity Sonar (f) Detailed prompt, Claude Opus 4.6 Figure 7: Search dynamics for the decoding module on the digits dataset. The blue line shows the accuracy of the solution obtained at each iteration, the red dashed line shows the best result achieved, and the gray horizontal line shows the baseline accuracy of 90.83 %. The columns correspond to analyst models, and the upper and lower rows correspond to the simple and detailed prompts, respectively. The curves demonstrate a pronounced dependence of convergence speed on the analyst model: with Claude Opus 4.6, the system reaches a plateau by iterations 5–10 (Fig. 6, Claude Opus 4.6 panels), whereas Perplexity Sonar in detailed-prompt mode exhibits stepwise dynamics with accuracy gains at iterations 50, 120, and 170; the system periodically returns to revising hypotheses after long plateaus. For the simple prompt in the Perplexity Sonar and Qwen3-30B-A3B-Instruct-2507 configurations (MNIST), the cumulative maximum does not rise above the baseline level over the entire range of 250 iterations: in these runs, the analyst proposed no decoder that outperformed argmax over aggregated spikes. This contrasts with the detailed-prompt version of the same models, where accuracy begins to increase in the first iterations, indicating the sensitivity of the agent system to the level of detail in the instruction provided to the analyst. The largest final gain over the baseline was obtained in the detailed-prompt configuration with Claude Opus 4.6 on digits (best-run accuracy 0.96540.9654 versus the digits baseline 0.90850.9085; Table 19), with the increase occurring almost entirely during the first 10 iterations, while the system spends the remaining more than 200 iterations confirming the maximum found without further improvement. E.2 Decoder catalog and TPOT configuration The same configuration was used in the classic and sequential TPOT modes, which differ in the search strategy but not in the space of admissible solutions, and in the Random Choice sampler described in Section 4.4. The classic mode used a population of 25 and 9 generations, the sequential mode a population of 1 and 249 consecutive generations. The budget of each mode was 250 evaluated candidates per run. Random Choice does not use TPOT as a search strategy: it draws configurations from the same catalog through independent calls to Qwen3-30B-A3B-Instruct-2507 and issues as many calls as are needed to record 250 evaluated decoders. Per run it made 347 calls on average. TPOT selected the best pipeline by mean accuracy in five-fold cross-validation on the training set. After the search was completed, the selected pipeline was retrained on the entire training set and evaluated once on the test set. Test data were not involved in selecting the pipeline structure or its hyperparameters. Mapping between method classes and operators. The catalog is based on the decoder taxonomy from the analyst’s detailed prompt. Some substantive classes map to the same scikit-learn estimator with different metrics or hyperparameters. Therefore, the number of TPOT operators does not match the number of classes in the taxonomy. For convenience, the mapping between decoder method classes and TPOT operators is provided in Table 20. Table 20: Mapping of decoder method classes to TPOT catalog operators. Method class Catalog operator Implementation Prototypes and centroids NearestCentroid Classification by distance to the class centroid k-nearest neighbors KNeighborsClassifier The number of neighbors, weighting rule, and metric are selected from the grid Correlation, cosine, and local-density rules KNeighborsClassifier Implemented using different local-proximity metrics Covariance methods, LDA, and QDA LinearDiscriminant Analysis, QuadraticDiscriminant Analysis Linear and quadratic discriminant rules Naive Bayes classifier GaussianNB Gaussian class-conditional feature model Threshold and tree-based rules DecisionTreeClassifier, ExtraTreeClassifier Threshold partitions of the spike-rate space Tree ensembles RandomForestClassifier, ExtraTreesClassifier Voting by an ensemble of threshold decoders Linear class scores LogisticRegression Linear decision rule without changing the SNN parameters Method combinations StackingEstimator, VotingClassifier Multistage pipelines and prediction aggregation Hyperparameter grids. The complete search-space configuration is provided below. Parameter values were selected only from the listed discrete sets. CONFIG_DETAILED_WITH_ENSEMBLE = "sklearn.neighbors.NearestCentroid": "metric": ["euclidean", "manhattan"], , "sklearn.neighbors.KNeighborsClassifier": "n_neighbors": range(1, 31), "weights": ["uniform", "distance"], "metric": [ "euclidean", "manhattan", "cosine", "chebyshev", ], , "sklearn.discriminant_analysis." "LinearDiscriminantAnalysis": "solver": ["svd", "lsqr", "eigen"], "shrinkage": [None, "auto"], , "sklearn.discriminant_analysis." "QuadraticDiscriminantAnalysis": "reg_param": [0.0, 1e-3, 1e-2, 1e-1, 0.5], , "sklearn.naive_bayes.GaussianNB": , "sklearn.tree.DecisionTreeClassifier": "criterion": ["gini", "entropy"], "max_depth": range(1, 11), "min_samples_split": range(2, 21), "min_samples_leaf": range(1, 21), , "sklearn.tree.ExtraTreeClassifier": "criterion": ["gini", "entropy"], "max_depth": range(1, 11), "min_samples_split": range(2, 21), "min_samples_leaf": range(1, 21), , "sklearn.ensemble.RandomForestClassifier": "n_estimators": [100, 200], "criterion": ["gini", "entropy"], "max_features": [0.5, 0.75, 1.0, "sqrt"], "min_samples_split": range(2, 11), "min_samples_leaf": range(1, 11), "bootstrap": [True, False], , "sklearn.ensemble.ExtraTreesClassifier": "n_estimators": [100, 200], "criterion": ["gini", "entropy"], "max_features": [0.5, 0.75, 1.0, "sqrt"], "min_samples_split": range(2, 11), "min_samples_leaf": range(1, 11), "bootstrap": [True, False], , "sklearn.linear_model.LogisticRegression": "penalty": ["l2"], "C": [0.01, 0.1, 1.0, 10.0, 100.0], "solver": ["lbfgs", "newton-cg"], "max_iter": [1000], , "tpot.builtins.StackingEstimator": "estimator": [ "sklearn.neighbors.NearestCentroid", "sklearn.neighbors.KNeighborsClassifier", "sklearn.discriminant_analysis." "LinearDiscriminantAnalysis", "sklearn.discriminant_analysis." "QuadraticDiscriminantAnalysis", "sklearn.naive_bayes.GaussianNB", "sklearn.tree.DecisionTreeClassifier", "sklearn.tree.ExtraTreeClassifier", "sklearn.ensemble.RandomForestClassifier", "sklearn.ensemble.ExtraTreesClassifier", "sklearn.linear_model.LogisticRegression", ], , "sklearn.ensemble.VotingClassifier": "voting": ["soft"], , The StackingEstimator and VotingClassifier operators were included in the catalog in all modes. Appendix F Multi-Agent Framework Prompts The original language of all prompts is Russian. F.1 Analyst prompt (simple) # Context You are a leading expert in spiking neural networks (SNNs) and in decoding outputs firing rates into class labels for the MNIST task, with over 20 years of practical experience. You are provided with: - Output spike sequences from 10 output neurons for the training and test datasets; - A JSON table containing experiments conducted using various decoding methods. The goal of this research is to develop a new decoding method that achieves an accuracy of 0.9 or higher (90%+) using a vector of firing rates from 10 output neurons. Important feature of the encoding: Input sequences are encoded by randomly selecting spikes without replacement from a given reference spike sequence. Therefore, the temporal structure of individual spike sequences is unstable across runs and should not be used directly. The informative object for decoding is the vector of firing rates of the 10 output neurons. # Your task Develop a method for decoding the firing-rate vector of the output neurons (dimension: 10) that is not included in the provided JSON experiment table, with the goal of achieving an accuracy of 0.9 or higher. **Available data and abilities:** 1. Output spike sequences vector (zeros and ones) for the training dataset - train_spikes. 2. Output spike sequences vector (zeros and ones) for the test dataset - test_spikes. 3. Correct labels for the training dataset - y_train. 4. You may construct new features and new decoding rules based on firing_rates. When developing a method, clearly distinguish between: - the feature-extraction step from the firing rate vector; - the decision rule step based on these features. Examples of acceptable method classes (which may be combined or extended, but are not limited to): - Normalized argmax over firing rates, optionally with neuron-specific weights and biases; - Rules based on the difference between the two highest firing rates (margin-based); - Class calibration using firing rate statistics for each class (e.g., mean firing rates for each neuron and class) and calss-specific score calculations; - Simple threshold-based and piecewise rules (e.g., small if/then/else/ decision trees based on firing rate thresholds and their combinations). **Limitations:** 1. The output spike trains have already been recorded. You cannot modify the network architecture or encoding method to record different sequences. 2. You cannot use backpropagation to train the proposed method. 3. You may fit method parameters only on the training data (train_spikes and y_train). The test data must be used exclusively for evaluation. 4. Your answer must be presented as a developer specification, including: - a clear description of how features are constructed from firing rates (optional); - a description of the decision rule in analytical form: - formulas for intermediate calculations and for the final transformation from firing rates to a class label; - recommended hyperparameters in the following JSON: `"decoding_params": `, where `decoding_params` is a dictionary of parameters, such as thresholds, weights, and search ranges. 5. Do not repeat methods that have already been tested and included in the JSON experiment table. 6. If a parameter is best tuned using the training data (e.g., thresholds, weights, or ranges), mark it with the value `"to_fit"` in `decoding_params` and provide a grid of candidate values for iteration. # Experiment Table (JSON): F.2 Analyst prompt (detailed) You are a leading neuroinformatics expert (20+ years of experience in spiking neural networks (SNNs) and in decoding outputs firing rates into class labels. Your specialty is developing decoding methods that operate WITHOUT backpropagation, relying exclusively on the statistical properties of output-neuron firing rates. # Dataset: dataset_name - Task: task_description - Output neurons: (n_neurons) - Classes: n_classes - Class labels: class_labels - **Target accuracy: target_accuracy+** ## Input data - `train_spikes` - a NumPy array `(N_train, n_neurons)`: the firing rates of the output neurons for the training set; - `test_spikes` - a NumPy array `(N_test, n_neurons0`: the firing rates of the output neurons for the test set; - `y_train` - a NumPy array `(N_train,)`: true class labels. ## Encoding Feature Input sequences are generated by randomly sampling spikes from a refrence sequence. Their temporal structure is **unstable** across runs. **Only** the firing rate vector of the n_neurons output neurons is informative. # Current experiments on the dataset_name Dataset A total of n_experiments experiments were conducted. Below is the table (JSON): ```json current_results_json ``` # Method Class Coverage Below is an automated analysis of which **method classes** have already been tested and which remain unexplored. Your task is to **prioritize methods from unexplored or underexplored classes** whenever they are theoretically applicable to the given problem. method_class_coverage_section **Method class selection rule:** If the table above contains method classes with the status "NOT STUDIED" or "LIMITED DATA", propose a method from one of those classes. Consider an already well-studied class only if you can justify why the proposed variation is fundamntally different from all previously tested methods. # Task Develop a **new** method for decoding firing rate vectors of dimension n_neurons that is **not included** in the experiment table, with the goal of achieving an accuracy of **target_accuracy+**. ## Mandatory Analytical Process Complete the Chain-of-Thought reasoning: ### Step 1 - Class Coverage Analysis Review the method class coverage table above. - Which classes have not been studied or remain underexplored? - Do any of those classes include methods that are theoretically suitable for this task? - What is the highest accuracy achieved by the methods already tested? - Which properties such as normalization, distance metrics, or regulatization may have contributed to thier performance? ### Step 2 - Failure Analysis Which methods degraded to near-chance accuracy (`~random_accuracy:.2f`)? Formulate hypotheses about the possible causes: - overfitting; - matrix degeneracy; - scale instability; - bugs involving array dimensions. **Important:** the failure of a specific implementation does not imply that the entire method class has failed. If kNN or Mahalanobis methods were implemented incorrectly - for example, using a singular matrix or improper normalization - this is a reason to propose **corrected** version rather than abandon the class. ### Step 3 - Selecting a class and synthesizing a method 1. **Priority:** Choose a method from the class with the smallest coverage. 2. If you chose a method from a well-studied class, explain how it fundamentally differs from all previous methofs in that class. 3. Use best practices from successful experiments, such as normalization and regularization, but apply them to a **new type of decision rule**. When developing the method, distinguish between: - **Feature extraction:** constructing features from firing rates. - **Decision rule:** assinging a class label based on these features. ### Step 4 - Self-Check (Required) Before submitting your final answer, verify the following: - [ ] Is the proposed method different from all methods already tested? - [ ] Does the method belong to a **different class** from the one used in most previous experiment? If not, explain why this particular class was selected. - [ ] When working with covariance matrices: was regulatization added (ridge: `Cov + eps * I`, or `np.linalg.pinv`, or diagonal approximation)? - [ ] Does the method scale to n_classes classes? - [ ] Is the output an array of integers from 0 to n_classes_minus_1? ## Method class Catalog Below is the complete catalog of permitted method classes. Methods from different classes may be combined. **All** methods are implementable using pure NumPy. ### 1. k-Nearest Neighbors (kNN) For each test example, find the k nearest neighbors in train_spikes and assign the class by voting. Metrics: Euclidean, Manhattan, Chebyshev, and inverse-distance-weighted voted (1/d). **Tip:** Apply Z-normalization before classification and select k using leave-one-out validation on the training data. ### 2. Mahalanobis Distance Calculate class centroids and covariance matrices, then assign each example to the class with the minimum Mahalanobis distance. Implementable in Numpy via `np.linalg.inv` (or `np.linalg.pinv` for robustness). **Tip:** Apply ridge regularization `Cov_c + eps * I` with eps ~ 1e-4...1e-1. You may also use a pooled covariance matrix instead of the class-specific one. ### 3. Naive Bayes (Gaussian NB) Estimate _c and _c for each neuron and class. Calculate the log-likelihood of each class as the sum of log-Gaussian probabilites, then select the class with the highest log-posterior using argmax. **Tip:** Apply variance smoothing: σ^2 = max(σ^2, ε). ### 4. Discriminant Analysis (LDA / QDA) LDA: general covariance matrix + linear discriminant functions. QDA: class-specific covariances. Fisher LDA: projection onto the axis of maximum interclass variance. All variants can be implemented using NumPy linear algebra. **Tip:** For LDA, solve the generalized eigenvalue problem S_w^-1 @ S_b for example using `np.linalg.eigh`. ### 5. Prototypes / Centroids Calculate a centroid for each class, then assign each example to the class with the minimum Euclidean or cosine distance. Options: L1/L2 normalization, weighted centroids, and median centroids for robustness. ### 6. Correlation Decoding Calculate a class-specific score: Score(x, c) = Pearson(x, template_c) or Score(x, c) = cosine(x, template_c). The class template is the mean or median vector. Assign the class label using Argmax over the class scores. ### 7. KDE / Density Estimation For each class, estimate the density p(x|c) from the training data using the kernel method (Gaussian kernel). Classify each example based on argmax p(x|c) * p(c). **Tip:** Select bandwidth using Silverman’s rule or by maximazing accuracy on the training data. ### 8. Rank-based methods Convert firing rates into ranks, work with the order of neurons. Options: rank-based Z-score, Spearman correlation, margin between ranks. ### 9. Margin-based rules Differences between top-K firing rates. Score = fr[top1] - mean(fr[rest]). Options: normalization, Z-score, temperature scaling. ### 10. Softmax Calibration Softmax with temperature fitted to the training data. Plus: biases by class. Z-normalization can be added before softmax. ### 11. Threshold/Piecewise-Specific Rules Cascade threshold rules (mini-trees), if neuron_k > threshold -> class c. Thresholds are selected on the training data. Options: binary codes, one-versus-rest thresholds. ### 12. Ensembles of heterogeneous methods Combine 2-3 **different** methods from the list above with voting or weighted averaging of scores. Select weights on the training data. # Response format - Technical Assignment (TA) Please formulate your answer **strictly** in the following structure: ## 1. Method name A short, unique name (e.g. "ARMD-v2"). ## 2. Method class Indicate which class from the catalog the method belongs to (for example: "kNN,", "Mahalanobis," "Naive Bayes"). ## 3. Method description One paragraph: what the method does, what mathematical idea it is based on. ## 4. Feature extraction Formulas for constructing features from firing rates: - What intermediate values are calculated? - What normalizations/transformations are applied? ## 5. Decision rule Analytical decision-making formula: - How is the score for each class derived from features? - How is the final label derived from scores? ## 6. Hyperparameters JSON: ```json "decoding_params":... ``` - Fixed parameters -- with specific values. - Parameters for training -- with the "to_fit" value and a grid specification. # Strict constraints 1. The spike sequences have already been recorded; you cannot modify the SNN. 2. **Prohibited:** backpropagation, gradient descent. 3. Parameter fiting is permitted **only** on train_spikes and y_train. Test data may be used exclusively for evaluation. 4. Do not repeat the methods from the experiment table. 5. Output: one-dimensional numpy array of integers (class_labels), of length N_test. F.3 Coder prompt # Task The technical specifications are provided below. Please write Python code for them. Write a class that decodes output spike sequences into class labels for the dataset_name clasification problem. Decoding must be based on firing rate vectors - i.e., spike counts aggregated over time - from n_neurons output neurons. **The class constructor takes exactly 2 parameters:** 1. train_spikes -- a two-dimensional NumPy array of the output firing rates of neurons in the training set of shape (N_train, 10). 2. y_train -- a one-dimensional NumPy array of the class labels of the training set of length (N_train,). **What does the __call__() method accept:** - test_spikes -- a two-dimensional NumPy array of the output firing rates of neurons on the test sample of shape (N_test, 10). **Output of the __call__() method:** y_pred is a one-dimensional Numpy array of length N_test **Restrictions:** - Only standard Python and NumPy are allowed; no other third-party libraries are allowed. - The class name is strictly: Decoding. - The class constructor takes exactly 2 parameters, described above. - If necessary, you can implement a fit method that doesn’t take any arguments and uses only self.train_spikes and self.y_train. - If there are parameters in decoding_params with "to_fit" values, you must implement their tuning on the training data (for example, by selecting thresholds or weights using a grid or another simple method). - After completing the tuning, all hyperparameters (including those marked as "to_fit") should be written to the self.decoding_params attribute as a dictionary with specific values. - The call occurs through the __call__(self, test_spikes) method. - Slice indices must be strictly integer; float/fractional indices are not allowed in slices. - The output is a strictly one-dimensional numpy array of integers from 0 to 9 inclusive, length test_spikes.shape[0]. - By the time the fit method (if implemented) completes, the class must have a self.decoding_params attribute, containing the final values of all the decoding method’s hyperparameters. **Strict withdrawal requirements:** - Return only the decoding code without any comments, text, or additional lines. Blank for this class: ```python class Decoding: def __init__(self, train_spikes, y_train): self.train_spikes = train_spikes self.y_train = y_train def __call__(self, test_spikes): pass # return y_pred ``` F.4 Recorder prompt Write a short description (approximately one paragraph) of the decoding method and provide the hyperparameters for the decoding method. Compose a short name for the experiment. Provide your answer as JSON with the following structure: ```json "name": "Short name", "description": "Description of what the Decoding layer does", "decoding_params": ``` References Chan et al. (2024) Chan, J.S., Chowdhury, N., Jaffe, O., et al., 2024. MLE-bench: Evaluating machine learning agents on machine learning engineering. arXiv preprint arXiv:2410.07095. URL: https://arxiv.org/abs/2410.07095, doi:10.48550/arXiv.2410.07095. Chen et al. (1977) Chen, L., Berger, R.L., Lominadze, J.G., Rosenbluth, M.N., Rutherford, P.H., 1977. Nonlinear saturation of the dissipative trapped-electron instability. Physical Review Letters 39, 754–757. URL: https://doi.org/10.1103/PhysRevLett.39.754, doi:10.1103/PhysRevLett.39.754. Chen et al. (2024) Chen, S., Liu, Y., Han, W., Zhang, W., Liu, T., 2024. A survey on LLM-based multi-agent system: Recent advances and new frontiers in application. arXiv preprint arXiv:2412.17481. URL: https://arxiv.org/abs/2412.17481, doi:10.48550/arXiv.2412.17481. Ghafarollahi and Buehler (2024) Ghafarollahi, A., Buehler, M.J., 2024. ProtAgents: Protein discovery via large language model multi-agent collaborations combining physics and machine learning. Digital Discovery 3, 1389–1409. URL: https://pubs.rsc.org/en/content/articlelanding/2024/d/d4d00013g, doi:10.1039/D4D00013G. Guo et al. (2024) Guo, T., Chen, X., Wang, Y., et al., 2024. Large language model based multi-agents: A survey of progress and challenges. arXiv preprint arXiv:2402.01680. URL: https://arxiv.org/abs/2402.01680, doi:10.48550/arXiv.2402.01680. He et al. (2025) He, X., You, L., Tian, H., et al., 2025. Lang-PINN: From language to physics-informed neural networks via a multi-agent framework. arXiv preprint arXiv:2510.05158. URL: https://arxiv.org/abs/2510.05158, doi:10.48550/arXiv.2510.05158. Huang et al. (2024) Huang, X., Liu, W., Chen, X., et al., 2024. Understanding the planning of LLM agents: A survey. arXiv preprint arXiv:2402.02716. URL: https://arxiv.org/abs/2402.02716, doi:10.48550/arXiv.2402.02716. Jimenez et al. (2024) Jimenez, C.E., Yang, J., Wettig, A., et al., 2024. SWE-bench: Can language models resolve real-world GitHub issues?, in: The Twelfth International Conference on Learning Representations. URL: https://openreview.net/forum?id=VTF8yNQM66. Kim et al. (2022) Kim, Y., Li, Y., Park, H., et al., 2022. Neural architecture search for spiking neural networks, in: European Conference on Computer Vision. URL: https://arxiv.org/abs/2201.10355, doi:10.48550/arXiv.2201.10355. Kudryashov (2026) Kudryashov, N.A., 2026. Traveling wave solutions of the Kuramoto–Sivashinsky equation with nonlinear convection. Physics Letters A 585, 131644. Kunitsyn et al. (2025) Kunitsyn, D., Sboev, A., Davydov, Y., et al., 2025. Probabilistic spiking neural network with correlation-based memristive synaptic update. Moscow University Physics Bulletin 80, S1127–S1133. URL: https://doi.org/10.3103/S002713492570300X, doi:10.3103/S002713492570300X. LeCun et al. (1998) LeCun, Y., Bottou, L., Bengio, Y., Haffner, P., 1998. Gradient-based learning applied to document recognition. Proceedings of the IEEE 86, 2278–2324. doi:10.1109/5.726791. Liu et al. (1976) Liu, C.S., Rosenbluth, M.N., Tang, W.M., 1976. Dissipative universal instability due to trapped electrons in toroidal systems and anomalous diffusion. The Physics of Fluids 19, 1040–1044. URL: https://doi.org/10.1063/1.861574, doi:10.1063/1.861574. Lu et al. (2024) Lu, C., Lu, C., Lange, R.T., et al., 2024. The AI Scientist: Towards fully automated open-ended scientific discovery. arXiv preprint arXiv:2408.06292. URL: https://arxiv.org/abs/2408.06292, doi:10.48550/arXiv.2408.06292. Moloshnikov et al. (2025a) Moloshnikov, I.A., Kuvakin, M.S., Sboev, A.G., 2025a. Neural network modeling of optical solitons described by the generalized nonlinear Schrödinger equation of the sixth order with high nonlinearity. Moscow University Physics Bulletin 80, S1068–S1076. URL: https://doi.org/10.3103/S0027134925703011, doi:10.3103/S0027134925703011. Moloshnikov et al. (2025b) Moloshnikov, I.A., Sboev, A.G., Kutukov, A.A., et al., 2025b. Analysis of neural network methods for obtaining soliton solutions of the nonlinear Schrödinger equation. Chaos, Solitons & Fractals 192, 115943. URL: https://doi.org/10.1016/j.chaos.2024.115943, doi:10.1016/j.chaos.2024.115943. Na et al. (2022) Na, B., Mok, J., Park, S., et al., 2022. AutoSNN: Towards energy-efficient spiking neural networks, in: International Conference on Machine Learning. URL: https://arxiv.org/abs/2201.12738, doi:10.48550/arXiv.2201.12738. Nam et al. (2026) Nam, J., Yoon, J., Chen, J., Shin, J., Arik, S., Pfister, T., 2026. MLE-STAR: Machine learning engineering agent via search and targeted refinement. Advances in Neural Information Processing Systems 38, 116692–116712. doi:10.52202/085713-3889. Olson et al. (2016) Olson, R.S., Bartley, N., Urbanowicz, R.J., Moore, J.H., 2016. Evaluation of a tree-based pipeline optimization tool for automating data science, in: Proceedings of the Genetic and Evolutionary Computation Conference 2016, ACM, New York, NY, USA. p. 485–492. URL: https://doi.org/10.1145/2908812.2908918, doi:10.1145/2908812.2908918. Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., et al., 2011. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research 12, 2825–2830. Perplexity AI (2026) Perplexity AI, 2026. Perplexity API documentation. https://docs.perplexity.ai. Accessed: 2026-05-08. Qwen Team (2025) Qwen Team, 2025. Qwen3-30B-A3B-Instruct-2507. URL: https://huggingface.co/Qwen/Qwen3-30B-A3B-Instruct-2507. Raissi et al. (2019) Raissi, M., Perdikaris, P., Karniadakis, G.E., 2019. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics 378, 686–707. doi:10.1016/j.jcp.2018.10.045. Rybka et al. (2024) Rybka, R., Davydov, Y., Sboev, A., Vlasov, D., Serenko, A., 2024. Image and audio data classification using bagging ensembles of spiking neural networks with memristive plasticity, in: Biologically Inspired Cognitive Architectures 2023, Springer. p. 741–749. doi:10.1007/978-3-031-50381-8_80. Sboev et al. (2025) Sboev, A., Kunitsyn, D., Davydov, Y., Vlasov, D., Serenko, A., Rybka, R., Liu, Y., 2025. Direct correlational spike-timing-dependent plasticity learning applied to classification tasks, in: Neural Information Processing. ICONIP 2024, Springer, Singapore. p. 227–241. doi:10.1007/978-981-96-6579-2_16. Schmidgall et al. (2025) Schmidgall, S., Su, Y., Wang, Z., et al., 2025. Agent Laboratory: Using LLM agents as research assistants. arXiv preprint arXiv:2501.04227. URL: https://arxiv.org/abs/2501.04227, doi:10.48550/arXiv.2501.04227. Shi et al. (2026) Shi, T., Li, Y., Wang, Z., et al., 2026. Knowledge-driven autonomous materials research via collaborative multi-agent and robotic system. Matter 9, 102577. doi:10.1016/j.matt.2025.102577. Shinn et al. (2023) Shinn, N., Cassano, F., Berman, E., et al., 2023. Reflexion: Language agents with verbal reinforcement learning, in: Advances in Neural Information Processing Systems. URL: https://arxiv.org/abs/2303.11366, doi:10.48550/arXiv.2303.11366. Skorokhodov et al. (2025) Skorokhodov, M.S., Latalin, V.E., Rybka, R.B., et al., 2025. Agent-based control of robot behavior using reasoning language models. Automation and Remote Control 86, 265–277. URL: https://doi.org/10.1134/S0005117925601101, doi:10.1134/S0005117925601101. Sudan and Rosenbluth (1976) Sudan, R.N., Rosenbluth, M.N., 1976. Stability of field-reversed ion rings in a background plasma. Physical Review Letters 36, 972–975. URL: https://doi.org/10.1103/PhysRevLett.36.972, doi:10.1103/PhysRevLett.36.972. Swanson et al. (2025) Swanson, K., Wu, W., Bulaong, N.L., et al., 2025. The Virtual Lab of AI agents designs new SARS-CoV-2 nanobodies. Nature 646, 716–723. URL: https://w.nature.com/articles/s41586-025-09442-9, doi:10.1038/s41586-025-09442-9. Takahara et al. (2025) Takahara, I., Mizoguchi, T., Liu, B., 2025. Accelerated inorganic materials design with generative AI agents. arXiv preprint arXiv:2504.00741. URL: https://arxiv.org/abs/2504.00741, doi:10.48550/arXiv.2504.00741. Wang et al. (2023) Wang, G., Xie, Y., Jiang, Y., Mandlekar, A., Xiao, C., Zhu, Y., Fan, L., Anandkumar, A., 2023. Voyager: An open-ended embodied agent with large language models. arXiv preprint arXiv:2305.16291. Wei et al. (2022) Wei, J., Wang, X., Schuurmans, D., et al., 2022. Chain-of-thought prompting elicits reasoning in large language models, in: Advances in Neural Information Processing Systems, p. 24824–24837. URL: https://arxiv.org/abs/2201.11903, doi:10.48550/arXiv.2201.11903. Xi et al. (2023) Xi, Z., Chen, W., Guo, X., et al., 2023. The rise and potential of large language model based agents: A survey. arXiv preprint arXiv:2309.07864. URL: https://arxiv.org/abs/2309.07864, doi:10.48550/arXiv.2309.07864. Yamada et al. (2025) Yamada, Y., Lange, R.T., Lu, C., et al., 2025. The AI Scientist-v2: Workshop-level automated scientific discovery via agentic tree search. arXiv preprint arXiv:2504.08066. URL: https://arxiv.org/abs/2504.08066, doi:10.48550/arXiv.2504.08066. Yang et al. (2025) Yang, X., Yang, X., Fang, S., Zhang, Y., Wang, J., Xian, B., Li, Q., Li, J., Xu, M., Li, Y., et al., 2025. R&D-Agent: An LLM-agent framework towards autonomous data science. arXiv preprint arXiv:2505.14738. Yao et al. (2023a) Yao, S., Yu, D., Zhao, J., et al., 2023a. Tree of thoughts: Deliberate problem solving with large language models, in: Advances in Neural Information Processing Systems. URL: https://arxiv.org/abs/2305.10601, doi:10.48550/arXiv.2305.10601. Yao et al. (2023b) Yao, S., Zhao, J., Yu, D., et al., 2023b. ReAct: Synergizing reasoning and acting in language models, in: International Conference on Learning Representations. URL: https://arxiv.org/abs/2210.03629, doi:10.48550/arXiv.2210.03629. Zavertyaev et al. (2024) Zavertyaev, S.V., Moloshnikov, I.A., Sboev, A.G., et al., 2024. Neural network modeling of optical solitons described by generalized nonlinear Schrödinger equations. Moscow University Physics Bulletin 79, S666–S675. URL: https://doi.org/10.3103/S0027134924702114, doi:10.3103/S0027134924702114.

Instant research discovery

A Multi-Agent LLM System for Automated Tuning of Specialized Neural Network Models

Intelligence

Summary

Entities (11)

Relation Signals (10)

Abstract

Tags

Links

PDF

Full Text

Paper	Published	Area	Tags	Intel	Citations