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

Highlights 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 • 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 pp; selected GNLSE errors drop up to 2500x and 325x.

arXiv:2610.13663v1 [cs.AI] 9 Oct 2026

• 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 Moloshnikova,∗ , Dmitry Kunitsyna , Alexander Sboeva,b,∗∗ , Maksim Skorokhodova , Roman Rybkaa,c and Yuanchao Liud a National Research Centre “Kurchatov Institute”, 1 Akademika Kurchatova Square, Moscow, 123182, Moscow, Russia b National Research Nuclear University MEPhI, 31 Kashirskoye Shosse, Moscow, 115409, Moscow, Russia c Moscow Institute of Physics and Technology, 9 Institutskiy Lane, Dolgoprudny, 141701, Moscow, Russia d Harbin Institute of Technology, No. 92 West Dazhi Street, Nangang District, Harbin, 150001, Heilongjiang, China

ARTICLE INFO

ABSTRACT

Keywords: Large language models Multi-agent systems Automated machine learning Physics-informed neural networks Spiking neural networks

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 invariantpreservation 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.

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, Chen, Guo et al. (2023); Yao, Zhao, Yu et al. (2023b). Such agents are already used in programming, ML engineering, and individual scientific disciplines Jimenez, Yang, Wettig et al. (2024); Chan, Chowdhury, Jaffe et al. (2024); Lu, Lu, Lange et al. (2024); Schmidgall, Su, Wang 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. ∗ Corresponding author

∗∗ Principal corresponding author

ivan-rus@yandex.ru (I. Moloshnikov); Kunitsyn_DE@nrcki.ru (D. Kunitsyn); sag111@mail.ru (A. Sboev); Skorokhodov_MS@nrcki.ru (M. Skorokhodov); Rybka_RB@nrcki.ru (R. Rybka); ycliuharbin@163.com (Y. Liu) ORCID (s): 0000-0002-5000-7257 (I. Moloshnikov); 0009-0002-9640-6787 (D. Kunitsyn); 0000-0002-6921-4133 (A. Sboev); 0000-0001-9399-9332 (M. Skorokhodov); 0000-0002-5595-6398 (R. Rybka); 0009-0007-0430-7658 (Y. Liu)

Moloshnikov et al.: Preprint submitted to Elsevier

Page 1 of 40

A Multi-Agent LLM System for Automated Tuning

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, Chen, Wang et al. (2024); Huang, Liu, Chen et al. (2024) identify typical components: a planning module (Chain-of-Thought Wei, Wang, Schuurmans et al. (2022), Tree-of-Thoughts Yao, Yu, Zhao et al. (2023a)), dynamic memory and self-reflection Shinn, Cassano, Berman et al. (2023), calls to external tools, and multi-agent designs with role separation Chen, Liu, Han, Zhang and Liu (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, Yoon, Chen, Shin, Arik and Pfister (2026), including 36 % gold medals) and to approximately 35 % across all difficulty levels (R&D-Agent on GPT-5) Yang, Yang, Fang, Zhang, Wang, Xian, Li, Li, Xu, Li et al. (2025). Similar multi-agent designs are used in materials discovery and biochemistry Takahara, Mizoguchi and Liu (2025); Shi, Li, Wang et al. (2026); Ghafarollahi and Buehler (2024); Swanson, Wu, Bulaong et al. (2025). Voyager Wang, Xie, Jiang, Mandlekar, Xiao, Zhu, Fan and Anandkumar (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, Lange, Lu et al. (2025), Agent Laboratory Schmidgall et al. (2025)) and agent-based robot control Skorokhodov, Latalin, Rybka 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, Sboev, Kutukov et al. (2025b); separate heuristics have been proposed for each component, but no unified automatic procedure has been established. LangPINN He, You, Tian 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, Li, Park et al. (2022); Na, Mok, Park et al. (2022); formally, neural architecture search (NAS) also includes the generation of individual functional modules within a predefined space of Moloshnikov et al.: Preprint submitted to Elsevier

Page 2 of 40

A Multi-Agent LLM System for Automated Tuning

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.

Moloshnikov et al.: Preprint submitted to Elsevier

Page 3 of 40

A Multi-Agent LLM System for Automated Tuning

Designation of connections

Human Researcher sets a task once

required when using detailed prompt only

Detailed-prompt extension: 12 method classes + coverage map

required steps regardless of the prompt if an error occurs during code testing (≤ 10 attempts)

One iteration 250 iterations per run (SNN) / 50 (PINNs) Prompt Task Experiment log

Analyst output method class Technical Specification

LLM Coder codes under a fixed software contract

LLM Analyst hypothesis generation

Traceback and error message

Sandboxed task pipeline execution (Docker)

Metric

1. Code testing on a small subset 2. Run on the full training and test datasets

LLM Recorder describers the performed experiment as JSON

Code for editable module(s) of pipeline

An example of a pipeline with 2 editable blocks: Input → [Fixed] → [Editable] → [Fixed] → [Editable] → Output → Metric

JSON record name, description, params, metric, status, method class

Update Experiment log → Next iteration

Experiment log table of JSON records (external memory)

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.

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.

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. Moloshnikov et al.: Preprint submitted to Elsevier

Page 4 of 40

A Multi-Agent LLM System for Automated Tuning

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 spikingnetwork 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, Kunitsyn, Davydov, Vlasov, Serenko, Rybka and Liu (2025); Kunitsyn, Sboev, Davydov 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 𝐫 ∈ ℝ𝐾 ; the agent-modifiable part is the Decoder(𝐫) 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, Bottou, Bengio and Haffner (1998) (28×28 images, 𝑁 = 784 spike generators) and the scikit-learn digits dataset Pedregosa, Varoquaux, Gramfort, Michel, Thirion, Grisel, Blondel, Prettenhofer, Weiss, Dubourg et al. (2011) (8 × 8, 𝑁 = 64). Both contain 10 classes, so the network has 𝐾 = 10 output LIF neurons. Accuracy is used as the target metric. As Moloshnikov et al.: Preprint submitted to Elsevier

Page 5 of 40

A Multi-Agent LLM System for Automated Tuning

Figure 2: SNN pipeline. The input encoding, the SNN trained with STDP, and the output activity vector 𝐫 (left and middle) are fixed. Only the decoder module Decoder(𝐫) (right) is generated and optimized by the framework.

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.

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, Moloshnikov, Sboev et al. (2024); Moloshnikov, Kuvakin and Sboev (2025a), and the second is the generalized Kuramoto–Sivashinsky equation (6) from Kudryashov (2026). ( ) 𝑖𝑞𝑡 + 𝑞𝑥𝑥 + |𝑞|2 𝑞 1 − 𝛼|𝑞|2 + 𝛽|𝑞|4 = 0,

(1)

where the complex field 𝑞(𝑥, 𝑡) = 𝑢(𝑥, 𝑡)+𝑖𝑣(𝑥, 𝑡) is approximated by a fully connected network, while the equation (𝑝𝑑𝑒 ), initial-condition (𝑖𝑐 ), and boundary-condition (𝑏𝑐 ) residuals enter the loss functional () according to the Physics-Informed Neural Networks framework Raissi, Perdikaris and Karniadakis (2019), a detailed description is given in the appendix A. In all configurations, 𝛽 = 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 𝐹𝐼 , which maps the coordinates (𝑥, 𝑡) to the MLP input features, and the output transformation 𝐹𝑂 , which maps the MLP outputs and coordinates to the predicted physical field. For the GNLSE problems, the prediction consists of the two real components (𝑢, 𝑣); for the Kuramoto–Sivashinsky problem, it is the single real field 𝑢. The agents may change the mathematical forms of 𝐹𝐼 and 𝐹𝑂 and may increase the intermediate MLP output dimension when required by 𝐹𝑂 , but cannot modify the remaining architecture, loss functional, or training procedure. In the baseline, 𝐹𝐼 passes (𝑥, 𝑡) to the MLP unchanged and 𝐹𝑂 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).

Moloshnikov et al.: Preprint submitted to Elsevier

Page 6 of 40

A Multi-Agent LLM System for Automated Tuning

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: ) ( ⎧𝑖 𝑞𝑡 + 𝑞𝑥𝑥 + |𝑞|2 𝑞 1 − 𝛼|𝑞|2 + 𝛽|𝑞|4 = 0, ⎪ ̃ 0), ⎪𝑞(𝑥, 0) = 𝑞(𝑥, ⎨ ⎪𝑞(−25, 𝑡) = 0, 𝑞(25, 𝑡) = 0, ⎪(𝑥, 𝑡) ∈ [−25, 25] × [0, 3], ⎩

(2)

where 𝛼 = 0.1, 𝛽 = 0, and the initial profile 𝑞̃ is specified by the analytical soliton (12) with parameters 𝑘 = 1, 𝑤 = 0.88, 𝑥0 = 0, 𝜃0 = 0, and 𝜇 = 4(𝑘2 − 𝑤). Because an analytical solution is known for√this problem, quality is evaluated directly using the relative amplitude error Relℎ = ‖ℎtrue − ℎpred ‖2 ∕‖ℎtrue ‖2 , ℎ = 𝑢2 + 𝑣2 . This problem provides a transparent test: the agent can rely on the explicit structure of the solution. Soliton problem with a Gaussian perturbation: ( ) ⎧𝑖 𝑞𝑡 + 𝑞𝑥𝑥 + |𝑞|2 𝑞 1 − 0.3|𝑞|2 = 0, ⎪ 2 ⎪𝑞(𝑥, 0) = 𝑞2 (𝑥, 0) + 0.2 𝑒−𝑥 ∕50 , ⎨ ⎪𝑞(−50, 𝑡) = 0, 𝑞(50, 𝑡) = 0, ⎪(𝑥, 𝑡) ∈ [−50, 50] × [0, 25], ⎩ Moloshnikov et al.: Preprint submitted to Elsevier

(3)

Page 7 of 40

A Multi-Agent LLM System for Automated Tuning

where the soliton component 𝑞2 has the same form as (12) with 𝛼 = 0.3, 𝑘 = 1, 𝑤 = 0.88, 𝑥0 = −30, and 𝜃0 = 0. The Gaussian packet is 𝐴𝑔 exp[−(𝑥 − 𝐵𝑔 )2 ∕(2𝐶𝑔2 )] with 𝐴𝑔 = 0.2, 𝐵𝑔 = 0, and 𝐶𝑔 = 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 ( ) (4) Resbvp = 31 ⟨|GNLSE |⟩ + ⟨|𝑞pred (𝑥, 0) − 𝑞(𝑥, 0)|⟩ + ⟨|𝑞pred |𝜕 ⟩ , where GNLSE = 𝑖𝑞𝑡 +𝑞𝑥𝑥 +|𝑞|2 𝑞(1−0.3|𝑞|2 ), the averages are taken over a uniform evaluation grid, and 𝜕 denotes the Dirichlet boundaries 𝑥 = ±50. The second primary metric is the relative conservation error of the integrals of motion 𝐼1 = ∫ |𝑞|2 𝑑𝑥 and 𝐼2 = ∫ (𝑞 ∗ 𝑞𝑥 − 𝑞𝑞𝑥∗ ) 𝑑𝑥. For each invariant, we compute its maximum relative drift from the initial value over the considered time interval, | | max𝑡 |𝐼𝑗 (𝑡) − 𝐼𝑗 (0)| | |, Rel𝐼𝑗 = | | |𝐼𝑗 (0)| | |

𝑗 ∈ {1, 2}.

These individual errors are combined into the overall conservation metric Rel𝐼 = Rel𝐼1 + Rel𝐼2 . A segmented PINN solution from Zavertyaev et al. (2024) is retained as an auxiliary amplitude reference Relℎ = ‖ℎref −ℎpred ‖2 ∕‖ℎref ‖2 . On the agent domain that field itself satisfies Resbvp = 9.11⋅10−5 and Rel𝑖 = 9.87⋅10−4 . Relℎ has lowest priority: if a hypothesis improves both Resbvp and Rel𝑖 relative to this reference, it is regarded as a more accurate solution of the boundary-value problem than the stored PINN field, and Relℎ 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, Rosenbluth and Tang (1976); Sudan and Rosenbluth (1976); Chen, Berger, Lominadze, Rosenbluth and Rutherford (1977). We write the source equation for the scalar field 𝑣(𝑥, 𝑡) in the normalization 𝑣𝑡 + 𝑣𝑥𝑥 + 𝜎𝑣𝑥𝑥𝑥 + 𝑣𝑥𝑥𝑥𝑥 + 𝛾𝑣𝑣𝑥 + 𝜒(𝑣2 )𝑥𝑥 = 0.

(5)

Here 𝛾 = 2𝛽 is introduced to avoid confusion with the coefficient 𝛽 in (1). For 𝛽 = 3𝜎 and 𝜎 = −1, the coefficients used in this study are 𝛾 = 2𝛽 = 6𝜎 = −6 and 𝜒 = 3. The corresponding exact traveling wave has the nonzero background −1∕3. We therefore remove this background by introducing 1 𝑣(𝑥, 𝑡) = 𝑢(𝑥, 𝑡) − , 3 so that the neural network approximates the zero-background scalar field 𝑢(𝑥, 𝑡). 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 𝑞(𝑥, 𝑡), 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 𝐹𝐼 and 𝐹𝑂 as in the GNLSE problems. Traveling-wave problem with a Gaussian perturbation: ⎧ + 3(𝑢2 )𝑥𝑥 = 0, ⎪𝑢𝑡 − 𝑢𝑥𝑥 − 𝑢𝑥𝑥𝑥 + 𝑢𝑥𝑥𝑥𝑥 − 6𝑢𝑢𝑥 + 2𝑢 [𝑥 ] ⎪ (𝑥 + 6)2 1 , 𝑢(𝑥, 0) = + 0.02 exp − ⎪ 8 2 cosh2 (𝑥∕2) ⎨ ⎪𝑢(−30, 𝑡) = 𝑢 (−30, 𝑡) = 0, 𝑢(40, 𝑡) = 𝑢 (40, 𝑡) = 0, 𝑥𝑥 𝑥𝑥 ⎪ ⎪(𝑥, 𝑡) ∈ [−30, 40] × [0, 10]. ⎩ Moloshnikov et al.: Preprint submitted to Elsevier

(6)

Page 8 of 40

A Multi-Agent LLM System for Automated Tuning

Figure 4: PINN architecture. The coordinates (𝑥, 𝑡) pass through the transformation layer 𝐹𝐼 , the fully connected network, and the transformation layer 𝐹𝑂 , producing predictions 𝑢, ̂ 𝑣. ̂ Automatic differentiation computes the derivatives for the loss function. The loss function () consists of the residual in the differential equation (𝑝𝑑𝑒 ) and the errors in the initial (𝑖𝑐 ) and boundary conditions (𝑏𝑐 ). The optimizer updates the network weights.

The unperturbed pulse has the form 𝑢ex (𝑥, 𝑡) =

1 (speed 1, amplitude 1∕2, starting at 𝑥0 = 0); the Gaussian 2 cosh2 ((𝑥−𝑡)∕2)

packet is located to the left (𝑥𝑔 = −6, width 𝑤 = 2) and catches up with the pulse at approximately 𝑡 ≈ 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 = 41 ⟨|KS |⟩ + ⟨|𝑢pred (𝑥, 0) − 𝑢(𝑥, 0)|⟩ + ⟨|𝑢pred |𝜕 ⟩ + ⟨|𝜕𝑥𝑥 𝑢pred |𝜕 ⟩ , (7)

where KS is the residual (22), the averages are taken over a uniform 1000 × 10 grid in the domain, and 𝜕 denotes the two zero boundaries. The physical diagnostic is the relative mass drift of the predicted solution: | | |𝑀pred (𝑡𝑘 ) − 𝑀pred (0)| |, | | |𝑀pred (0)| | |

| Δ(𝑢) 𝑚 = max 𝑘

𝑀pred (𝑡) =

∫

𝑢pred (𝑥, 𝑡) 𝑑𝑥,

(8)

with the maximum taken over ten uniformly spaced time slices 𝑡𝑘 ∈ [0, 10].

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.

Moloshnikov et al.: Preprint submitted to Elsevier

Page 9 of 40

A Multi-Agent LLM System for Automated Tuning

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 %.

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, Bartley, Urbanowicz and Moore (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) = 250 evaluated candidates. • TPOT sequential uses a population of one individual and 249 consecutive generations, 1 × (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 Moloshnikov et al.: Preprint submitted to Elsevier

Page 10 of 40

A Multi-Agent LLM System for Automated Tuning

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 Relℎ relative to the baseline for the single-soliton GNLSE problem, a decrease in Resbvp and Rel𝑖 for the GNLSE problem with a Gaussian perturbation, and a decrease in Resbvp 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 𝐹𝐼 , 𝐹𝑂 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 ⟨language 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, Relℎ is the final characteristic of a hypothesis; the best Relℎ (mean over three internal seeds) is reported. In the GNLSE problem with a Gaussian perturbation, hypotheses are ranked by the conservation error Rel𝑖 , while the boundary-value-problem residual Resbvp (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 Resbvp (7), and the mass drift Δ(𝑢) 𝑚 (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 Relℎ = 4.69 ⋅ 10−2 ,

5.51 ⋅ 10−2 , and 4.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 Rel𝑖 = 1.28 ⋅ 10−1 , 2.15 ⋅ 10−1 , and 1.85 ⋅ 10−1 ; over 𝑛 = 10 additional baseline trainings, mean is Rel𝑖 = (1.65 ± 0.18) ⋅ 10−1 (Table 11). In the Kuramoto–Sivashinsky −2 problem, the baseline gives Resbvp = 1.41 ⋅ 10−3 and Δ(𝑢) 𝑚 = 9.31 ⋅ 10 . 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 Moloshnikov et al.: Preprint submitted to Elsevier

Page 11 of 40

A Multi-Agent LLM System for Automated Tuning Table 1 Baseline and best agent results across five tasks. Task (metric) SNN, MNIST (Acc) SNN, scikit-learn digits (Acc) PINN, GNLSE, single soliton (Relℎ ) PINN, GNLSE, soliton + Gaussian (Rel𝑖 ) ) PINN, KS, traveling wave + Gaussian (Δ(𝑢) 𝑚

Baseline (single run)

Best agent result (best in series)

Improvement over baseline (pp / factor)

80 % 91 %

93 % 97 %

+13 pp +6 pp

5.51 ⋅ 10−2 2.15 ⋅ 10−1 9.31 ⋅ 10−2

2.21 ⋅ 10−5 6.61 ⋅ 10−4 6.56 ⋅ 10−3

2491× 325× 14×

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 𝑊 and 𝑝 values. The paired-analysis methodology and rationale for selecting the test are provided in Supplementary Section E. The unadjusted and Holm-adjusted 𝑝 values, together with the 𝑊 statistic for all six tests, are reported in Supplementary Section E; the difference is statistically significant in all six tests at 𝛼 = 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 Resbvp and Rel𝑖 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 Rel𝑖 and its two components but not Resbvp ; 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 − Accbase is reported in percentage points, whereas for PINN tasks the errorreduction factor Err 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 Rel𝑖 . For the Kuramoto–Sivashinsky problem, it reports the mass drift Δ(𝑢) 𝑚 , an independent physical diagnostic recorded alongside the optimized boundary-value-problem residual Resbvp ; 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 Rel𝑖 = (1.65 ± 0.18) ⋅ 10−1 for the baseline, (1.42 ± 0.59) ⋅ 10−3 for the configuration selected by Resbvp , and (1.60 ± 0.81) ⋅ 10−3 for the configuration selected by Rel𝑖 . The corresponding reductions relative to the baseline mean are 116× and 103×, respectively (Table 11). Overall, the best agent configuration outperforms the baseline on the reported metric in all five tasks.

Moloshnikov et al.: Preprint submitted to Elsevier

Page 12 of 40

A Multi-Agent LLM System for Automated Tuning 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.

Mean maximum accuracy ± std, %

Best accuracy, %

Runs with accuracy ≥ 90 %

simple detailed detailed without coverage map

91.13 ± 0.44 92.15 ± 0.62

91.65 92.83

10 out of 10 10 out of 10

91.02 ± 0.25

91.45

10 out of 10

simple detailed detailed without coverage map

79.88 ± 0.00 91.76 ± 0.56

79.88 92.70

0 out of 10 10 out of 10

82.34 ± 0.70

83.18

0 out of 10

simple detailed detailed without coverage map

79.88 ± 0.00 85.88 ± 2.42

79.88 89.16

0 out of 10 0 out of 10

80.18 ± 0.45

81.15

0 out of 10

Model

Prompt

Claude Opus 4.6 Claude Opus 4.6 Claude Opus 4.6 Perplexity Sonar Perplexity Sonar Perplexity Sonar Qwen3-30B-A3B-Instruct-2507 Qwen3-30B-A3B-Instruct-2507 Qwen3-30B-A3B-Instruct-2507

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 % to 92.15 ± 0.62 %, and the best result from 91.65 % to 92.83 %. It gives the higher run maximum in eight of the ten pairs, and the two-sided exact Wilcoxon signed-rank test yields 𝑊 = 5, 𝑝 = 0.0195, which remains significant after the Holm correction over the six tests of Section 4.6. The paired gain averages 1.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 % in all ten runs (std = 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 % to 92.70 % with a mean of 91.76 ± 0.56 %, so the 90 % reference point was exceeded in all ten runs, whereas for Qwen3-30B-A3BInstruct-2507 the mean maximum was 85.88 ± 2.42 % and the best run reached 89.16 % without crossing it. For both models the prompt effect is significant by the same paired test (𝑊 = 0, 𝑝 = 0.00195, the minimum attainable value at 𝑛 = 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 % to 96.12 ± 0.66 %, with all ten runs reaching 93 %. For Perplexity Sonar, the mean changed from 93.44 ± 0.57 % to 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-30BA3B-Instruct-2507, the mean increased from 91.63 ± 0.74 % to 93.66 ± 0.79 %, and the number of runs reaching 93 % increased from 0/10 to 7/10.

5.3. Cost of the agent search The positions 𝑘min and 𝑘max 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 𝐶run , and the cost 𝐶 ∗ accumulated up to the earliest position 𝑘min . 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 Moloshnikov et al.: Preprint submitted to Elsevier

Page 13 of 40

A Multi-Agent LLM System for Automated Tuning 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

TPOT TPOT Random Choice

classic sequential LLM sampler

89.60 89.38 89.48

90.15 89.93 89.93

201 243 67 210 -

detailed

82.28

89.16

15

209

detailed detailed

90.93 91.22

92.70 92.83

8 2

147 43

Qwen3-30BA3B-Instruct-2507 Perplexity Sonar Claude Opus 4.6

Acc𝑚𝑖𝑛 , % Acc𝑚𝑎𝑥 , % k𝑚𝑖𝑛 k𝑚𝑎𝑥

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: • Acc𝑚𝑖𝑛 , Acc𝑚𝑎𝑥 — the smallest and the largest best accuracies. • k𝑚𝑖𝑛 , k𝑚𝑎𝑥 — 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 𝑘min and 𝑘max are not reported.

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.

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 𝑟, out out in 𝐶𝑟 = 𝑝in 𝑟 𝑇𝑟 + 𝑝𝑟 𝑇𝑟 ,

𝐶run = 𝐶an + 𝐶pr + 𝐶rec ,

(9)

where 𝑝𝑟 are the per-token rates and 𝑇𝑟 are the token counts taken from the run logs. The Analyst input at iteration 𝑖 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), 𝐹 = 4085 tokens, and the experiment log of 𝑖 records of 𝑑log 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 𝑑log in the prompt exceeds the length 𝑑 of the description text itself, which is what the Recorder returns. Over the 𝑁 = 250 accepted iterations of a run the Analyst also makes 𝑁an −𝑁 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 𝑑log 𝑁∕2 — the average log length over the run — the total Analyst input is in 𝑇an = 𝑁an 𝐹 + 𝑑log

(

𝑁(𝑁 + 1) 𝑁 + (𝑁an − 𝑁) 2 2

) .

(10)

The first term is the fixed prompt part paid at each of the 𝑁an Analyst calls; the second term collects the log tokens read during the 𝑁 accepted iterations (exact sum) and during the 𝑁an − 𝑁 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 𝑑log = 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 𝑁an = 265, 252, and 255 Moloshnikov et al.: Preprint submitted to Elsevier

Page 14 of 40

A Multi-Agent LLM System for Automated Tuning 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 Performance Accmin , % Accmax , % 𝑘min 𝑘max

Claude det.

Sonar det.

Qwen det.

Random Choice

91.22 92.83 2 43

90.93 92.70 8 147

82.28 89.16 15 209

89.48 89.93 — —

Resources per run, tokens in millions (M) Evaluated decoders 250 250 250 250 Analyst calls, 𝑁an 265 252 255 347 Programmer calls / iter 3.31 1.70 2.05 — Analyst in / out, M 11.10 / 1.17 10.89 / 0.58 11.20 / 0.58 0.28 / 0.042 Programmer in / out, M 12.98 / 4.66 2.40 / 1.73 2.79 / 1.32 — Recorder in / out, M 2.61 / 0.058 1.71 / 0.061 1.38 / 0.063 — Cost Analyst rate, $/M in / out Programmer, Recorder rate, $/M 𝐶an , $ 𝐶pr , $ 𝐶rec , $ 𝐶run , $ 𝐶 ∗ up to 𝑘min , $

5 / 25 1/1 0.15 / 0.60 0.25 / 2.00 (gpt-5-mini) 84.66 11.47 2.03 12.56 4.06 3.35 0.77 0.55 0.47 97.99 16.08 5.85 0.37 0.21 0.26

0.15 / 0.60 — 0.07 — — 0.07 0.022

in = 11.10, 10.89, and 11.20 M input tokens per run, of which the (Table 4). Equation (10) with these values gives 𝑇an fixed-prompt contribution 𝑁an 𝐹 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 out = 𝑁 ⋅ 𝑇 out = 1.17, 0.58, and 0.58 M per run. Qwen3-30B-A3B-Instruct-2507, so that 𝑇an an an,call 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-A3BInstruct-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 in = 2.61, 1.71, and 1.38 M and 𝑇 out = 𝑁𝑑 = 0.058, 0.061, and 0.063 M returns the description of length 𝑑, so that 𝑇rec rec 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 𝑘min 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 % Moloshnikov et al.: Preprint submitted to Elsevier

Page 15 of 40

A Multi-Agent LLM System for Automated Tuning Table 5 PINN results for the single soliton (Relℎ 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 Relℎ

Best Relℎ

Change relative to baseline

sonar-pro anthropic-claude-opus-4.6 qwen-qwen3-30b-a3b-instruct-2507

4.69 ⋅ 10−2 5.51 ⋅ 10−2 4.85 ⋅ 10−2

1.97 ⋅ 10−4 2.21 ⋅ 10−5 8.86 ⋅ 10−3

238× better 2491× better 5.5× better

threshold. When Qwen3-30B-A3B-Instruct-2507 is used free of charge, 𝐶run = 0 for both configurations, and they differ only in token usage and in 𝐶 ∗ . 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) 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 𝐹𝐼 , 𝐹𝑂 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 Relℎ = 2.21 ⋅ 10−5 , a 2491-fold reduction from 5.51 ⋅ 10−2 . The best sonar-pro hypothesis reduces Relℎ from 4.69 ⋅ 10−2 to 1.97 ⋅ 10−4 (238×), while Qwen3-30B-A3B-Instruct-2507 reduces it from 4.85 ⋅ 10−2 to 8.86 ⋅ 10−3 (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. The best sonar-pro configuration is tw_nonlinear_pinn_core. Its 𝐹𝐼 supplies a normalized traveling-wave coordinate, normalized time, and a local nonlinearity-density feature; 𝐹𝑂 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 𝐹𝑂 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; 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 Resbvp (4) and Rel𝑖 ; the reruns reported here do not include Relℎ . The human-designed PINN from prior work has Resbvp = 9.11 ⋅ 10−5 and Rel𝑖 = 9.87 ⋅ 10−4 . The singlerun results are summarized in Table 6. The shared identity baseline combines the lowest Rel𝑖 and the lowest Resbvp observed across the analyst runs, and each reduction factor is computed relative to its Rel𝑖 . The 325× improvement is obtained relative to the identity baseline from the same Claude Opus 4.6 run, whereas the 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 Rel𝑖 from 1.28 ⋅ 10−1 to 6.61 ⋅ 10−4 (194×) while giving Resbvp = 3.95 ⋅ 10−5 . A different Claude configuration, SplitAlpha_ExpPhase_HighP, gives the lowest residual, Moloshnikov et al.: Preprint submitted to Elsevier

Page 16 of 40

A Multi-Agent LLM System for Automated Tuning 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 Rel𝑖 are computed relative to this shared baseline. Model

Result

Human-designed PINN from prior work Zavertyaev et al. (2024)

Rel𝑖

Reported result

Identity baseline (best across runs)

Change in Rel𝑖 relative to baseline

Resbvp −4

−5

130× better

9.87 ⋅ 10

9.11 ⋅ 10

1.28 ⋅ 10−1

1.04 ⋅ 10−3

—

sonar-pro

Best Rel𝑖

6.35 ⋅ 10−3

2.37 ⋅ 10−2

20× better

anthropic-claude-opus-4.6 anthropic-claude-opus-4.6

Best Rel𝑖 Best Resbvp

6.61 ⋅ 10−4 1.40 ⋅ 10−3

3.95 ⋅ 10−5 2.08 ⋅ 10−5

194× better 91× better

qwen-qwen3-30b-a3b-instruct-2507

Best Rel𝑖

9.43 ⋅ 10−2

1.38 ⋅ 10−3

1.36× better

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). Resbvp is the metric the analyst optimized; Δ(𝑢) 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 Baseline (𝐹𝐼 , 𝐹𝑂 identity) Best Resbvp (9× better) Best Δ(𝑢) (14× better) 𝑚

Resbvp

Δ(𝑢) 𝑚

1.41 ⋅ 10−3 1.51 ⋅ 10−4 1.70 ⋅ 10−4

9.31 ⋅ 10−2 3.47 ⋅ 10−2 6.56 ⋅ 10−3

Hypothesis — MultiplicativeIC_DualAdditive_DualZoom BoundedMultAdd_SechZoom_DualPhase

Resbvp = 2.08 ⋅ 10−5 , with Rel𝑖 = 1.40 ⋅ 10−3 . The best sonar-pro hypothesis by Rel𝑖 improves this metric 20×, but its Resbvp = 2.37 ⋅ 10−2 is worse than the shared baseline; the best Qwen hypothesis improves Rel𝑖 only 1.36× and likewise does not improve the baseline residual. Thus, only the Claude hypotheses improve both reported criteria substantially. Ten additional trainings from random initializations were performed for the identity baseline and for both Claude configurations. The configuration selected by Resbvp gives Rel𝑖 = (1.42 ± 0.59) ⋅ 10−3 , a 116× reduction from the repeated baseline mean, while the configuration selected by Rel𝑖 gives (1.60 ± 0.81) ⋅ 10−3 , a 103× reduction. The repeated-run export does not contain Resbvp , 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 𝐹𝐼 , 𝐹𝑂 transformations while minimizing the boundary-value-problem residual Resbvp . The mass drift Δ(𝑢) 𝑚 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 (Resbvp ), the best architecture is MultiplicativeIC_DualAdditive_DualZoom, which reduces Resbvp from the baseline value of 1.41 ⋅ 10−3 to 1.511 ⋅ 10−4 , a factor of 9.3. On the physical diagnostic, the −2 to 6.56 ⋅ 10−3 , an best architecture is BoundedMultAdd_SechZoom_DualPhase, which reduces Δ(𝑢) 𝑚 from 9.31 ⋅ 10 approximately 14-fold reduction, while also improving Resbvp by a factor of 8.3. ( ) Both winning hypotheses share a multiplicative–additive decomposition 𝑢̂ = 𝑢(𝑥, 0) 𝑒−𝛼𝑡 1 + 𝑡 𝛾𝑚 𝑁1 + ( ) 𝑡 𝑊 (𝑥) 𝛾1 𝑁2 + 𝛾2 𝑁3 , in which 𝑁1,2,3 are the three MLP outputs, the initial profile is embedded exactly at 𝑡 = 0, the window 𝑊 (𝑥) = 𝑠𝐿 𝑠𝑅 built from the two normalized distances to the boundaries vanishes at 𝑥 = −30 and 𝑥 = 40, and the input block supplies features in the co-moving coordinate 𝜉 = 𝑥 − 𝑡 — a rescaled tanh(𝜉∕𝜆), the profile sech2 (𝜉∕2) of the unperturbed pulse, and a phase marker sin(𝜋𝑡∕6) of the collision. The hypothesis that wins on mass drift differs in one detail: the multiplicative channel is bounded, 𝑁1 ↦ tanh(𝑁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 𝑀 = ∫ 𝑢 𝑑𝑥 an order of magnitude closer to its initial value, which illustrates that the two criteria are not interchangeable.

Moloshnikov et al.: Preprint submitted to Elsevier

Page 17 of 40

A Multi-Agent LLM System for Automated Tuning

Figure 5: Predicted solution of the Kuramoto–Sivashinsky problem at 𝑡 = 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 𝑥, the vertical axis is 𝑢. The baseline pulse loses roughly half of its amplitude by 𝑡 = 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∕2.

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 % to 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 signedrank test (𝑊 = 5, 𝑝 = 0.0195), with the detailed prompt yielding higher accuracy in eight of ten pairs; the mean paired difference is 1.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-Instruct2507. 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 %, so the 90 % threshold was exceeded in all ten runs. For Qwen3-30B-A3B-Instruct-2507, the mean maximum was 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.

Moloshnikov et al.: Preprint submitted to Elsevier

Page 18 of 40

A Multi-Agent LLM System for Automated Tuning

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 % in all ten runs. The detailed prompt provided gains ranging from 2.40 to 9.28 percentage points in individual runs. The mean maximum accuracy was 85.88 ± 2.42 %, and the best result was 89.16 %. The high variability (std = 2.42 percentage points at 𝑛 = 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 %, which is 1.14 percentage points below the detailed prompt (92.15 %) and nearly matches the level of the simple prompt (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 % (Δ = 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 % (Δ = 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-A3BInstruct-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 %) 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 𝑘min = 2 to 𝑘max = 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 rederives relations between experiments from the log at every step, which limits the depth of multi-step reasoning that can be sustained over long sessions. Moloshnikov et al.: Preprint submitted to Elsevier

Page 19 of 40

A Multi-Agent LLM System for Automated Tuning

• Bounded number of iterations. All reported runs are capped at 250 iterations for SNN tasks and 50 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 LLMdriven 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-A3BInstruct-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.

Moloshnikov et al.: Preprint submitted to Elsevier

Page 20 of 40

A Multi-Agent LLM System for Automated Tuning

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.

CRediT authorship contribution statement Ivan Moloshnikov: Conceptualization, Methodology, Software, Investigation, Data curation, Writing – original draft. Dmitry Kunitsyn: Methodology, Investigation, Software, Writing – original draft. Alexander Sboev: Conceptualization, Methodology, Investigation, Data curation, Supervision, Validation, Writing – original draft. Maksim Skorokhodov: Software, Validation, Investigation, Visualization, Writing – original draft. Roman Rybka: Supervision, Project administration. Yuanchao Liu: Formal analysis, Methodology.

A. Baseline PINN model The model problem is a boundary-value problem for the generalized nonlinear Schrödinger equation (GNLSE) with cubic–quintic nonlinearity: ( ) ⎧𝑖 𝑞𝑡 + 𝑞𝑥𝑥 + |𝑞|2 𝑞 1 − 𝛼|𝑞|2 + 𝛽|𝑞|4 = 0, ⎪ ⎪𝑞(𝑥, 0) = 𝑞(𝑥, ̃ 0), ⎨ ⎪𝑞(𝑥min , 𝑡) = 0, 𝑞(𝑥max , 𝑡) = 0, ⎪ ⎩(𝑥, 𝑡) ∈ [𝑥min , 𝑥max ] × [𝑡min , 𝑡max ].

(11)

In all configurations, 𝛽 = 0; the quintic term is retained in the equation for generality. Two configurations are considered. In the single-soliton problem, 𝛼 = 0.1 and the domain is (𝑥, 𝑡) ∈ [−25, 25] × [0, 3]; the initial profile is specified by the analytical soliton solution: √ √ 𝜇𝑧 ⋅ 𝑒𝑖(𝑘𝑥+𝜃0 ) , 𝑧 = 𝑒 𝜇(𝑥−𝑥0 ) , (12) 𝑞(𝑥, ̃ 0) = ( ) 2 𝛼𝜇 2 1 𝑧 + 1 − 𝑧 2 3 Moloshnikov et al.: Preprint submitted to Elsevier

Page 21 of 40

A Multi-Agent LLM System for Automated Tuning 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 Neurons per layer Activation (hidden layers) Activation (output layer) Optimizer Learning rate Number of training iterations Collocation points (𝑁𝑓 ) Initial-condition points (𝑁0 ) Boundary-condition points (𝑁𝑏 ) Collocation-point sampling Weights 𝑤pde , 𝑤ic , 𝑤bc

2 32 sin linear Adam 10−3 10 000 1 000 100 50 fixed random 1, 1, 1

4 64 sin linear Adam 10−3 30 000 30 000 3 000 300 adaptive (v3) 1, 1, 1

where 𝑘 = 1, 𝑤 = 0.88, 𝑥0 = 0, 𝜃0 = 0 and 𝜇 = 4(𝑘2 −𝑤). The complex field is decomposed into two real components: 𝑞(𝑥, 𝑡) = 𝑢(𝑥, 𝑡) + 𝑖 𝑣(𝑥, 𝑡). The analytical soliton is used as the reference when calculating the error. In the second configuration (a soliton with a Gaussian perturbation), 𝛼 = 0.3, the domain is extended to (𝑥, 𝑡) ∈ [−50, 50] × [0, 25], and the initial condition is the sum of the soliton profile (12) (with 𝑥0 = −30) and a real-valued Gaussian packet: 2

(13)

𝑞(𝑥, 0) = 𝑞2 (𝑥, 0) + 0.2 𝑒−𝑥 ∕50 ,

where 𝑞2 coincides with (12) for the specified parameters, equivalently 𝐴𝑔 exp[−(𝑥 − 𝐵𝑔 )2 ∕(2𝐶𝑔2 )] with 𝐴𝑔 = 0.2, 𝐵𝑔 = 0, and 𝐶𝑔 = 5 (the printed 𝐶𝑔2 in Zavertyaev et al. (2024) is a typographical error, must be 2𝐶𝑔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] that field has Resbvp = 9.11 ⋅ 10−5 and Rel𝑖 = 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 in 𝑃𝑟 𝑙 in Eq. 6 of Moloshnikov et al. (2025b)). The primary metrics of the search are Resbvp (4) and the relative conservation error of the integrals of motion Rel𝑖 = Rel𝑖1 + Rel𝑖2 ; Relℎ has lowest priority. The baseline is a fully connected network MLP ∶ ℝ2 → ℝ2 , (𝑥, 𝑡) ↦ (𝑢, 𝑣), with two hidden layers of 32 neurons each and sin activation. In the modified architecture, two tunable blocks, 𝐹𝐼 (𝑥, 𝑡) and 𝐹𝑂 (𝑢, 𝑣, 𝑥, 𝑡), 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: (14)

 = pde + ic + bc ,

where pde is the mean squared residual of Eq. (11) at the collocation points, ic is the initial-condition error, and bc is the boundary-condition error. The derivatives 𝑞𝑡 , 𝑞𝑥𝑥 are obtained by TensorFlow automatic differentiation. For the single-soliton problem the metric is: Relℎ =

‖ℎtrue − ℎpred ‖2 ‖ℎtrue ‖2

,

ℎ=

√ 𝑢2 + 𝑣2 .

(15)

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 𝑅0 = 10 before squaring, the learning rate is lowered Moloshnikov et al.: Preprint submitted to Elsevier

Page 22 of 40

A Multi-Agent LLM System for Automated Tuning Table 9 Hyperparameters of the baseline PINN for the Kuramoto–Sivashinsky traveling-wave problem with a Gaussian perturbation. Parameter

Value

Number of hidden layers Neurons per layer Network output Activation (hidden layers) Activation (output layer) Optimizer Learning rate Residual scale 𝑅0 Number of training iterations Collocation points (𝑁𝑓 ) Initial-condition points (𝑁0 ) Boundary-condition points (𝑁𝑏 ) Collocation-point sampling Adaptive-sampling parameters Evaluation grid (𝑛𝑥 × 𝑛𝑡 ) Weights 𝑤pde , 𝑤ic , 𝑤bc

4 32 𝑢 (one real scalar) sin linear Adam, clipnorm = 1 5 ⋅ 10−5 10 50 000 20 000 1 000 200 adaptive (v3), every 1 000 iterations 𝜆2 = 0.01, oversampling ×5 1000 × 10 1, 1, 1

to 5 ⋅ 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 𝑢 and 𝑢𝑥𝑥 at the two boundaries, so the loss reads ⟩ ⟩ ⟨ ⟩ ⟨ ⟨  = (KS ∕𝑅0 )2 + (𝑢pred (𝑥, 0) − 𝑢(𝑥, 0))2 + 𝑢2 + 𝑢2𝑥𝑥 𝜕 .

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 (16)

𝑣𝑡 + 𝑣𝑥𝑥 + 𝜎𝑣𝑥𝑥𝑥 + 𝑣𝑥𝑥𝑥𝑥 + 𝛾𝑣𝑣𝑥 + 𝜒(𝑣2 )𝑥𝑥 = 0,

where 𝛾 = 2𝛽, so that the nonlinear convection term integrates to 𝛽𝑣2 under a traveling-wave reduction. For 𝛽 = 3𝜎, substitution of 𝜎 = −1, 𝛾 = 2𝛽 = 6𝜎 = −6, and 𝜒 = 3 gives (17)

𝑣𝑡 + 𝑣𝑥𝑥 − 𝑣𝑥𝑥𝑥 + 𝑣𝑥𝑥𝑥𝑥 − 6𝑣𝑣𝑥 + 3(𝑣2 )𝑥𝑥 = 0.

The exact wave from Kudryashov (2026) is expressed on the background −1∕3. To make the learned field decay to zero away from the pulse, we set 1 𝑣=𝑢− , 3

1 𝑢=𝑣+ . 3

(18)

Because the shift is constant, all differentiated linear terms are unchanged: 𝑣𝑡 = 𝑢𝑡 ,

𝑣𝑥 = 𝑢𝑥 ,

𝑣𝑥𝑥 = 𝑢𝑥𝑥 ,

𝑣𝑥𝑥𝑥 = 𝑢𝑥𝑥𝑥 ,

The nonlinear convection term transforms as ( ) 1 −6𝑣𝑣𝑥 = −6 𝑢 − 𝑢 = −6𝑢𝑢𝑥 + 2𝑢𝑥 , 3 𝑥 whereas the nonlinear diffusion term becomes ( ) 1 2 2 1 𝑣2 = 𝑢 − = 𝑢2 − 𝑢 + , 3 3 9 Moloshnikov et al.: Preprint submitted to Elsevier

𝑣𝑥𝑥𝑥𝑥 = 𝑢𝑥𝑥𝑥𝑥 .

(19)

(20) Page 23 of 40

A Multi-Agent LLM System for Automated Tuning

(21)

3(𝑣2 )𝑥𝑥 = 3(𝑢2 )𝑥𝑥 − 2𝑢𝑥𝑥 . Inserting these expressions into Eq. (17) yields 𝑢𝑡 + 𝑢𝑥𝑥 − 𝑢𝑥𝑥𝑥 + 𝑢𝑥𝑥𝑥𝑥 − 6𝑢𝑢𝑥 + 2𝑢𝑥 + 3(𝑢2 )𝑥𝑥 − 2𝑢𝑥𝑥 = 0. Combining the two 𝑢𝑥𝑥 terms gives the residual used in the implementation:

(22)

KS = 𝑢𝑡 − 𝑢𝑥𝑥 − 𝑢𝑥𝑥𝑥 + 𝑢𝑥𝑥𝑥𝑥 − 6𝑢𝑢𝑥 + 2𝑢𝑥 + 3(𝑢2 )𝑥𝑥 .

Thus, the negative coefficient of 𝑢𝑥𝑥 and the additional linear-advection term 2𝑢𝑥 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 1 1 , 𝑣ex (𝑥, 𝑡) = − + 3 2 cosh2 ((𝑥 − 𝑡 − 𝑥0 )∕2)

𝑢ex (𝑥, 𝑡) =

1 . 2 cosh ((𝑥 − 𝑡 − 𝑥0 )∕2) 2

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, Davydov, Sboev, Vlasov and Serenko (2024) are not used because preliminary experiments on MNIST showed no appreciable gain. The approach is applied to two datasets: MNIST (28 × 28, 𝑁 = 784 spike generators) and the scikit-learn digits Pedregosa et al. (2011) (8 × 8, 𝑁 = 64). Both have ten classes, pixel values are normalized by the 𝐿1 norm of the input vector, and the number of output LIF neurons is 𝐾 = 10. The network has two layers: the input layer contains 𝑁 spike generators, and the output layer contains 𝐾 Leaky-Integrate-and-Fire (LIF) neurons; the connections are fully connected and trainable, with no lateral connections. The input vector 𝑉 = (𝑣1 , … , 𝑣𝑁 ) is converted into an array of correlated spike sequences {𝑆𝑖 }: a reference sequence 𝑆0 of length 𝑁𝑡 is constructed with constant spacing 𝑁ref ; component 𝑣𝑖 is interpreted as the target Pearson correlation coefficient between 𝑆𝑖 and 𝑆0 , and 𝑆𝑖 is obtained by randomly removing some spikes from 𝑆0 . Supervision is implemented using two currents applied to the output LIF neurons: for a sample of class 𝑘, a positive current 𝐼exc = 107 pA, synchronized with the 𝑆0 spikes with a shift 𝜏𝑠 , is applied to the corresponding neuron; the other neurons receive an inhibitory current 𝐼inh = −103 pA. The weights are adjusted using a local multiplicative STDP rule with a bounding symmetric nearest-neighbor scheme: { 𝐴+ ⋅ 𝑤𝜇+ ⋅ 𝑒−Δ𝑡∕𝜏+ , Δ𝑡 > 0, Δ𝑤(Δ𝑡) = 𝐴− ⋅ 𝑤𝜇− ⋅ 𝑒Δ𝑡∕𝜏− , Δ𝑡 < 0, where Δ𝑡 = 𝑡post − 𝑡pre , 𝑤 is the current weight, 𝐴± 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.

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. Moloshnikov et al.: Preprint submitted to Elsevier

Page 24 of 40

A Multi-Agent LLM System for Automated Tuning Table 10 Hyperparameters of the baseline SNN and its training, common to all configurations. Parameter

Value

𝑉𝑟𝑒𝑠𝑡 (resting potential) 𝑉𝑟𝑒𝑠𝑒𝑡 (post-spike potential) 𝑉𝑡ℎ (threshold potential) 𝜏𝑚 (membrane time constant) 𝜏𝑠 (teacher-spike shift) 𝑁𝑟𝑒𝑓 (spacing between spikes in 𝑆0 ) 𝑁𝑡 (number of time steps in 𝑆0 ) 𝑇𝑠𝑖𝑚 (simulation time, 𝑑𝑡 = 0.1 ms) 𝐴+ (potentiation amplitude) 𝐴− (depression amplitude) 𝜏+ (potentiation window) 𝜏− (depression window) 𝜇+ (weight dependence, potentiation) 𝜇− (weight dependence, depression)

−70 mV −70 mV −56 mV 103 ms 0.1 ms 3 3000 300 ms 1.2 ⋅ 10−3 10−3 200 ms 160 ms 0.96 0.26

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(𝑑) =

5 ) ( 1∑ Acc 𝑑 (−𝑘) , (𝑘) , 5 𝑘=1

(23)

where 𝑑 (−𝑘) denotes the decoder trained on four folds, and (𝑘) is the fifth fold excluded from training. The analyst used CVAcc(𝑑) when formulating the next hypothesis. After the specified search budget was exhausted, the decoder with the highest CVAcc(𝑑) 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.

C. Analysis of the best PINN hypotheses The single-soliton comparison includes Perplexity Sonar (sonar-pro), Claude Opus 4.6, and Qwen3-30B-A3BInstruct-2507 (Table 5). The agent controls the input transformation 𝐹𝐼 and output transformation 𝐹𝑂 , 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 Relℎ = 2.21 ⋅ 10−5 , a 2491-fold reduction from its identity baseline of 5.51 ⋅ 10−2 . It uses seven analytical features, including normalized and logarithmic envelopes and two amplitude–time cross-features, while 𝐹𝑂 applies time-gated exponential amplitude and additive phase corrections with separate envelope weights. The best sonar-pro configuration, tw_nonlinear_pinn_core, gives Relℎ = 1.97 ⋅ 10−4 , a 238-fold reduction from its identity baseline of 4.69 ⋅ 10−2 . Its 𝐹𝐼 supplies a normalized traveling-wave coordinate, normalized time, and a local nonlinearity-density feature; 𝐹𝑂 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 Relℎ = 8.86 ⋅ 10−3 , a 5.48-fold reduction from its identity baseline of 4.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.

Moloshnikov et al.: Preprint submitted to Elsevier

Page 25 of 40

A Multi-Agent LLM System for Automated Tuning 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 Rel𝑖 . Configuration Baseline (𝑛 = 10) Residual-selected configuration (𝑛 = 10) Invariant-selected configuration (𝑛 = 10)

Mean ± std

Median

Min–max

CV

Reduction factor in Rel𝑖 relative to baseline mean

(1.65 ± 0.18) ⋅ 10−1 (1.42 ± 0.59) ⋅ 10−3 (1.60 ± 0.81) ⋅ 10−3

1.66 ⋅ 10−1 1.50 ⋅ 10−3 1.59 ⋅ 10−3

1.29 ⋅ 10−1 –1.94 ⋅ 10−1 5.40 ⋅ 10−4 –2.41 ⋅ 10−3 7.57 ⋅ 10−4 –3.31 ⋅ 10−3

11.1 % 41.8 % 50.4 %

1× 116.1× 103.3×

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.

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 agentdesigned 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. 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−4 , lower than in any of its ten subsequent retrainings. By contrast, the original error of the residualselected configuration, 1.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.

E. Additional classification results Method-class coverage map. The aggregate coverage map of decoder method classes from the complete 250-

iteration log for the ⟨Perplexity 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.

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. When the simple and detailed prompts were compared, the mean maximum accuracy for Qwen3-30B-A3BInstruct-2507 increased from 79.88 ± 0.00 % to 85.88 ± 2.42 %. A positive difference was also observed in all ten pairs and averaged 6.00 percentage points. For Claude Opus 4.6, the mean maximum increased from 91.13 ± 0.44 % to 92.15 ± 0.62 %. The mean paired difference was 1.02 percentage points: a positive effect was observed in 8 of 10 pairs, whereas accuracy decreased by 0.07 and 0.11 percentage points in two pairs. For Perplexity Sonar, switching to the Moloshnikov et al.: Preprint submitted to Elsevier

Page 26 of 40

A Multi-Agent LLM System for Automated Tuning Table 12 Coverage map of decoder method classes for one 250-iteration run (Perplexity Sonar, detailed prompt, MNIST). Method class

Methods

Best acc, %

154 32 13 11 11 9 8 5 4 3 0 0

92.70 87.90 82.32 88.56 83.73 86.05 81.68 80.43 85.75 73.93 — —

Ensembles LDA / QDA Softmax / calibration 𝑘-nearest neighbors Naive Bayes Mahalanobis Rank-based Prototypes / centroids KDE / density Threshold methods / rules Correlation / cosine Margin-based

Status explored explored explored explored explored limited data limited data limited data limited data limited data not explored not explored

Table 13 Maximum accuracy in ten matched pairs of complete agent-system runs with the simple and detailed prompts for Qwen330B-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, %

Δ, pp

1 2 3 4 5 6 7 8 9 10

79.88 79.88 79.88 79.88 79.88 79.88 79.88 79.88 79.88 79.88

82.28 83.35 88.57 86.48 89.16 84.10 87.20 85.50 83.90 88.30

2.40 3.47 8.69 6.60 9.28 4.22 7.32 5.62 4.02 8.42

Mean Std. dev.

79.88 0.00

85.88 2.42

6.00 2.42

detailed prompt increased the mean maximum accuracy from 79.88 ± 0.00 % to 91.76 ± 0.56 %. A positive difference was observed in all ten pairs, with a mean of 11.89 percentage points. The coverage-map ablation showed an additional effect of this component of the detailed prompt. For Qwen3-30BA3B-Instruct-2507, mean accuracy increased from 80.18 ± 0.45 % to 85.88 ± 2.42 %, i.e., by 5.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 % without the map and 91.76 ± 0.56 % with the map, with a mean increase of 9.42 percentage points. For Claude Opus 4.6, the presence of the map increased the mean maximum accuracy from 91.02 ± 0.25 % to 92.15 ± 0.62 %; the mean paired increase was 1.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 𝑊 statistic was defined as the smaller of the rank sums for positive and negative nonzero differences. The significance level was set to 𝛼 = 0.05. For the comparison of the simple and detailed prompts for Claude Opus 4.6, 𝑊 = 5 and 𝑝 = 0.0195 were obtained. For Perplexity Sonar and Qwen3-30B-A3B-Instruct-2507, all ten paired differences were positive. Therefore, in both

Moloshnikov et al.: Preprint submitted to Elsevier

Page 27 of 40

A Multi-Agent LLM System for Automated Tuning 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, %

Δ, pp

1 2 3 4 5 6 7 8 9 10

82.28 83.35 88.57 86.48 89.16 84.10 87.20 85.50 83.90 88.30

79.88 80.42 79.88 81.15 79.88 80.73 79.88 79.88 80.25 79.88

2.40 2.93 8.69 5.33 9.28 3.37 7.32 5.62 3.65 8.42

Mean Std. dev.

85.88 2.42

80.18 0.45

5.70 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, %

Δ, pp

1 2 3 4 5 6 7 8 9 10

79.88 79.88 79.88 79.88 79.88 79.88 79.88 79.88 79.88 79.88

92.70 92.57 91.95 91.82 91.77 91.74 91.64 91.45 91.08 90.93

12.82 12.69 12.07 11.94 11.89 11.86 11.76 11.57 11.20 11.05

Mean Std. dev.

79.88 0.00

91.76 0.56

11.89 0.56

cases the two-sided exact Wilcoxon signed-rank test gave the same result: 𝑊 = 0, 𝑝 = 0.001953. The null hypothesis was rejected for all three models at the 𝛼 = 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 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, i.e., approximately 26.49 %. To ensure control of the family-wise probability of at least one type I error at 0.05, the Holm–Bonferroni correction was applied. This procedure controls the family-wise error even when the statistical tests are dependent. The original

Moloshnikov et al.: Preprint submitted to Elsevier

Page 28 of 40

A Multi-Agent LLM System for Automated Tuning 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, %

Δ, pp

1 2 3 4 5 6 7 8 9 10

92.70 92.57 91.95 91.82 91.77 91.74 91.64 91.45 91.08 90.93

83.18 82.92 82.77 82.66 82.65 82.47 82.15 82.04 81.87 80.71

9.52 9.65 9.18 9.16 9.12 9.27 9.49 9.41 9.21 10.22

Mean Std. dev.

91.76 0.56

82.34 0.70

9.42 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, %

Δ, pp

1 2 3 4 5 6 7 8 9 10

90.20 91.15 91.65 91.60 91.33 90.71 91.42 91.17 90.88 91.22

92.83 91.88 91.58 92.49 91.22 92.35 92.69 92.57 92.66 91.24

2.63 0.73 -0.07 0.89 -0.11 1.64 1.27 1.40 1.78 0.02

Mean Std. dev.

91.13 0.44

92.15 0.62

1.02 0.90

𝑝 values were ordered in ascending order and sequentially compared with the thresholds 𝛼 𝛼 , , … , 𝛼. 6 5 For the five comparisons with an original 𝑝 = 0.001953, the adjusted value was 𝑝Holm = 0.0117. For the comparison of the simple and detailed prompts for Claude Opus 4.6, with an original 𝑝 = 0.0195, 𝑝Holm = 0.0195 was obtained. All six adjusted 𝑝 values were below 0.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, 𝑊 = 0, 𝑝 = 0.001953. After applying the Holm–Bonferroni correction across all six comparisons, the adjusted value was 𝑝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.02 percentage points for Claude Opus 4.6, 11.89 percentage points for Perplexity Sonar, and 6.00 percentage points for Qwen3-30B-A3B-Instruct2507. For Claude Opus 4.6, the two-sided exact Wilcoxon signed-rank test gave 𝑊 = 5, an original 𝑝 = 0.0195, and Moloshnikov et al.: Preprint submitted to Elsevier

Page 29 of 40

A Multi-Agent LLM System for Automated Tuning 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, %

Δ, pp

1 2 3 4 5 6 7 8 9 10

91.22 91.24 91.58 91.88 92.35 92.49 92.57 92.66 92.69 92.83

90.66 90.73 90.89 90.91 90.95 91.01 91.04 91.20 91.32 91.45

0.56 0.51 0.69 0.97 1.40 1.48 1.53 1.46 1.37 1.38

Mean Std. dev.

92.15 0.62

91.02 0.25

1.14 0.41

Table 19 Comparison of analyst models and prompt variants on the digits 8×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

𝑛

Claude Opus 4.6 Claude Opus 4.6 Perplexity Sonar Perplexity Sonar Qwen3-30B-A3B-Instruct-2507 Qwen3-30B-A3B-Instruct-2507

simple detailed simple detailed simple detailed

10 10 10 10 10 10

Mean maximum accuracy ± std, % 95.07 ± 0.52 96.12 ± 0.66 93.44 ± 0.57 93.31 ± 0.93 91.63 ± 0.74 93.66 ± 0.79

Best accuracy, % 95.56 96.67 93.89 94.17 92.78 94.17

Runs with accuracy ≥ 93 % 10 out of 10 10 out of 10 9 out of 10 8 out of 10 0 out of 10 7 out of 10

𝑝Holm = 0.0195. For Perplexity Sonar and Qwen3-30B-A3B-Instruct-2507, positive differences were obtained in all ten pairs; for each model, 𝑊 = 0, the original 𝑝 = 0.001953, and 𝑝Holm = 0.0117.

Results on the digits dataset. The same framework was applied to digits 8 × 8; the run design was the same as for MNIST. The baseline argmax decoder gives an accuracy of 90.83 %, whereas the best ⟨Claude 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 𝑐̂ = arg max𝑐 [log 𝑝(𝑓 ∣ 𝑐) + log 𝑝(𝑐)]; hyperparameters were selected using 5-fold cross-validation.

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. 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.

Moloshnikov et al.: Preprint submitted to Elsevier

Page 30 of 40

A Multi-Agent LLM System for Automated Tuning

(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.

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.9654 versus the digits baseline 0.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-30BA3B-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.

Moloshnikov et al.: Preprint submitted to Elsevier

Page 31 of 40

A Multi-Agent LLM System for Automated Tuning

(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.

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),

Moloshnikov et al.: Preprint submitted to Elsevier

Page 32 of 40

A Multi-Agent LLM System for Automated Tuning 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

𝑘-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

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

Method combinations

StackingEstimator, VotingClassifier

Covariance methods, LDA, and QDA

},

Linear decision rule without changing the SNN parameters Multistage pipelines and prediction aggregation

"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",

Moloshnikov et al.: Preprint submitted to Elsevier

Page 33 of 40

A Multi-Agent LLM System for Automated Tuning

},

}

],

"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.

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": {}}`,

Moloshnikov et al.: Preprint submitted to Elsevier

Page 34 of 40

A Multi-Agent LLM System for Automated Tuning 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_neurons}0`: 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?

Moloshnikov et al.: Preprint submitted to Elsevier

Page 35 of 40

A Multi-Agent LLM System for Automated Tuning ### 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 \mu_c and \sigma_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: \sigma^2 = max(\sigma^2, \epsilon). ### 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.

Moloshnikov et al.: Preprint submitted to Elsevier

Page 36 of 40

A Multi-Agent LLM System for Automated Tuning ### 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.

Moloshnikov et al.: Preprint submitted to Elsevier

Page 37 of 40

A Multi-Agent LLM System for Automated Tuning 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": {} } ```

Moloshnikov et al.: Preprint submitted to Elsevier

Page 38 of 40

A Multi-Agent LLM System for Automated Tuning

References 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, 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, 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, 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/dd/ d4dd00013g, doi:10.1039/D4DD00013G. 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, 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, 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, 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, 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, N.A., 2026. Traveling wave solutions of the Kuramoto–Sivashinsky equation with nonlinear convection. Physics Letters A 585, 131644. 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, 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, 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, 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, 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, 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, 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, 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, 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. pp. 485–492. URL: https://doi.org/10.1145/2908812.2908918, doi:10.1145/2908812.2908918. 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 API documentation. https://docs.perplexity.ai. Accessed: 2026-05-08. Qwen Team, 2025. Qwen3-30B-A3B-Instruct-2507. URL: https://huggingface.co/Qwen/Qwen3-30B-A3B-Instruct-2507. 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, 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. pp. 741–749. doi:10.1007/ 978-3-031-50381-8_80. 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. pp. 227–241. doi:10.1007/ 978-981-96-6579-2_16. 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, 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.

Moloshnikov et al.: Preprint submitted to Elsevier

Page 39 of 40

A Multi-Agent LLM System for Automated Tuning 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, 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, 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, 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://www.nature.com/articles/s41586-025-09442-9, doi:10.1038/s41586-025-09442-9. 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, 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, 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, pp. 24824–24837. URL: https://arxiv.org/abs/2201.11903, doi:10.48550/arXiv.2201.11903. 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, 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, 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, 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, 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, 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.

Moloshnikov et al.: Preprint submitted to Elsevier

Page 40 of 40

Related documents

• A Multi-Agent LLM System for Automated Tuning of Specialized Neural Network Models#1298159
• Learning Reusable Agent Subroutines from Experience#1344605
• Autonomous Assessment of Generalizability of AI Agent Capabilities#1110150
• Skill-MAS: Evolving Meta-Skill for Automatic Multi-Agent Systems#1131578
• AgentTune: A Unified Training Lifecycle for Tool-Using LLM Agents#1297687