Efficient dynamics simulation analysis method based on Fourier neural operator

Abstract

The invention belongs to the technical field of model simulation, and particularly relates to an efficient dynamics simulation analysis method based on a Fourier neural operator. The present invention proposes FNO-Speed, and a series of comprehensive solutions for inefficient operations that the FNO solver does not fully utilize hardware. According to the method, two unique optimization methods are adopted, and comprise a multi-level parallel implicit image-to-column general matrix multiplication optimization strategy and a user-defined size high-frequency signal filtering algorithm. According to the method, efficient general matrix multiplication is achieved through an implicit image-to-column and data division strategy to replace pointwise convolution, and fragmentary calculation of frequency domain local linear transformation is eliminated through the latter. The FNO-Speed makes full use of the memory bandwidth, improves the calculation efficiency, and aims to solve the problems of low utilization rate of calculation resources and delay influence caused by large-scale data access calculation.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of model simulation, and particularly relates to an efficient dynamic simulation analysis method based on a Fourier neural operator. BACKGROUND

[0002] In recent years, with the rise of the concept of minimally invasive treatment, local treatment represented by local ablation has developed rapidly and has become another important treatment method after surgical resection. Radiofrequency ablation has important application value in multi-target ablation treatment of localized malignant liver cancer. Ablation refers to a method of directly killing tumors locally under the guidance of imaging, which introduces physical energy into the tumor to increase or decrease the temperature of the target area, causing irreversible damage to the tumor, including radiofrequency, microwave, freezing, intratumoral anhydrous alcohol injection, etc., among which radiofrequency ablation is the most widely used local ablation method. Radiofrequency refers to radio frequency. When the current frequency of radio frequency is high to a certain value, it causes the movement of positive and negative ions in the tissue, frictional heat, and the formation of a 3-5c diameter spherical necrotic lesion in the tumor tissue. The accuracy of the radiofrequency ablation parameterization scheme determines the clinical diagnosis and treatment effect of radiofrequency ablation surgery. Precise dynamic analysis of liver cancer lesions quickly and accurately is the key to the scientific setting of the radiofrequency ablation parameterization scheme, and the current dynamic simulation analysis of liver cancer lesions still has problems such as "low precision", "slow speed" and "much experience intervention", which leads to difficulties in precise guidance of liver cancer radiofrequency ablation and greatly affects the surgical treatment effect. How to quickly and accurately simulate the dynamic behavior of liver cancer lesions and provide scientific basis for the formulation of clinical parameterization scheme has become one of the core challenges of this technology.

[0003] The dynamic analysis method of liver cancer lesions mainly includes methods relying on traditional blood flow network, electric heat pulse, fluid mechanics, electric heat boundary condition, and traditional partial differential equation solution. The dynamic solution analysis of the traditional partial differential equation solution usually relies on numerical solution method, mainly including finite difference method, finite element method and spectral method, etc. These methods can handle very complex partial differential equations and obtain high-precision solutions. However, these methods face many problems in actual solution. When the partial differential equation is high-dimensional or nonlinear, especially in complex systems and multi-scale and multi-physical properties, the numerical method is required to have high precision, which leads to heavy computational burden. The precision and timeliness of dynamic simulation analysis of liver cancer lesion model are one of the key factors affecting liver cancer radiofrequency ablation technology. How to construct a high-precision and fast dynamic simulation analysis method of liver cancer lesion model and improve the precision and timeliness of radiofrequency ablation simulation calculation is a key scientific problem to be solved in the process of precise implementation and clinical application of liver cancer radiofrequency ablation technology.

[0004] To cope with these challenges, researchers urgently need to find new methods that are more efficient than traditional numerical solvers. Artificial intelligence-enabled science is an emerging field that uses artificial intelligence technologies to accelerate and improve scientific discovery, the fifth paradigm, and is becoming a major trend in scientific research. In recent years, deep learning methods for solving partial differential equations have been widely used to replace numerical solvers for scientific computing problems. Artificial intelligence is helping scientists generate hypotheses, design experiments, collect and interpret large data sets, and gain insights that may not be possible through traditional scientific methods. Deep learning solver models can be used to efficiently analyze massive scientific data and solve complex physical equations, saving time required by traditional numerical solvers. Currently, there are data-driven and physics information-driven related works, which are widely used in the field of physics and practical engineering problems.

[0005] With the advent of neural operators, the method of training solvers has gradually shifted to focusing on approximation operators. Neural operator learning is one of the main methods currently applied to the solution of partial differential equations. Neural operator learning is based on the universal approximation theorem, learning the mapping of the solution function from one space to another, so that it can approximate any given nonlinear continuous operation, and then achieve the corresponding output for any input point. It aims to solve a class of problems and opens up new directions for the application of deep learning in the field of partial differential equations, greatly improving the efficiency of solving complex partial differential equations. The method of neural operator learning effectively loads the computational burden in advance to the training time. In the inference time, the speed of the neural operator solver can be several orders of magnitude faster than the corresponding numerical solver, and it avoids the expensive computational cost of repeated iterative solving caused by changes in problem conditions. Among them, the architecture based on the Fourier neural operator has been proven to be an efficient and high-precision method. It parameterizes the kernel in the Fourier domain and uses fast Fourier transform calculations to achieve a time complexity of O(J log J), which is three orders of magnitude faster than traditional PDE solvers. The PDE solver based on FNO performs well in liver cancer lesion dynamics simulation, and can meet the dual demands of clinical computing speed and solving accuracy, providing efficient and reliable technical support for precise planning of radiofrequency ablation surgery.

[0006] However, as the scale and dimensionality of real-world problems increase, Fourier neural operator solvers struggle with high-dimensional, nonlinear, or complex geometric problems. This is especially true when time-dependent problems, such as the Penners biological heat conduction equation used in liver cancer lesion dynamics simulation, extend beyond small-scale 2D time-varying scenarios—that is, problems involving two or three spatial dimensions plus a time dimension. This requires processing large amounts of high-dimensional data and thousands of training simulations, resulting in significant computational costs. Therefore, it is not feasible for complex time-varying problems in real-world scenarios. Due to the heterogeneity of tumor tissue and the diversity of blood vessel distribution, constructing FNO-based liver cancer lesion dynamics simulations faces challenges in high-dimensional solution domains and adapting to irregular and complex geometric boundary conditions, leading to long computation times and high computational costs. To address the clinical characteristics and real-time requirements of liver cancer radiofrequency ablation surgery, further optimization of the FNO solver's computational performance is necessary.

[0007] Existing deep learning frameworks (such as PyTorch) have deeply optimized the implementation of basic operations like matrix multiplication for routine batch data processing, effectively improving computational utilization. However, during the training and inference of Fourier Neural Operator (FNO) solver models, hardware efficiency is low due to the small-batch data processing mode and suboptimal implementation strategies for key operations. For example, trivial small-scale linear transformations in the Fourier Neural Operator layer cannot maximize the effective utilization of hardware cores and memory bandwidth, while introducing additional kernel startup overhead. As the dimensionality and complexity of the problem domain increase, redundant memory accesses and data movement of operators on large-scale data have a significant impact on time consumption; for example, special-shaped convolutions in the Fourier Neural Operator layer lead to increased training and inference time. These performance bottlenecks affect the computational efficiency of the FNO solver in high-resolution and time-sensitive problem scenarios such as the dynamic simulation analysis of liver cancer lesions, failing to provide effective solutions. Although much excellent work has been done on accelerating operators such as matrix multiplication on GPUs and parallel strategies on heterogeneous platforms, these general acceleration methods cannot maximize optimization effects on specific models. Summary of the Invention

[0008] This invention proposes an efficient dynamic simulation analysis method based on Fourier neural operators. This method addresses the performance bottleneck of Fourier neural operators and solves the problem of latency caused by low utilization of computing resources.

[0009] The technical solution of this invention is implemented as follows: An efficient dynamic simulation analysis method based on Fourier neural operators includes the following steps: S1. Perform multi-level parallel implicit image-to-column general matrix multiplication optimization: take advantage of data reuse opportunities to implicitly transform convolution calculations into matrix multiplications, and improve device memory bandwidth utilization by using vectorized memory access and matrix multiplication based on vector outer product. S2. Perform operator fusion optimization: fuse convolution calculation and subsequent element-wise addition to reduce the overhead of repeated reading and writing and storage of intermediate calculation results; S3. Perform custom-sized high-frequency signal filtering optimization: Adopt a data relocation scheme to perform customized splicing and transposition of low-frequency data in the frequency domain, merge subsequent linear transformation operations, avoid fragmented step memory access, make full use of memory bandwidth, and simplify user operations.

[0010] Through the above technical solutions, the three operation steps effectively solve the key bottlenecks of Fourier neural operators in dynamic simulation, such as computational density, inefficient memory access, and redundancy of high-frequency filtering operations. They significantly improve computational throughput and memory bandwidth utilization, reduce intermediate result storage and repeated read / write overhead, avoid fragmented memory access and simplify user operations, and achieve efficient and accurate dynamic simulation analysis.

[0011] Optionally, in step S1, in order to transform the convolution calculation into matrix multiplication, the input is converted into a large-scale matrix, where each column contains the input elements required for calculation at each convolution kernel position during the convolution process. The elements of each convolution kernel are then used to generate another matrix according to their corresponding positions during the convolution process. Finally, the result of the convolution is obtained by matrix multiplication of the two matrices.

[0012] Through the above technical solution, this operation transforms the originally spatially irregular memory access pattern into accessing and calculating two consecutive, densely stored large matrices, providing a unified computing carrier for subsequent optimization methods such as vectorized memory access and hierarchical storage allocation, and realizing the regularization and efficiency of convolution calculation.

[0013] Optionally, in step S1, to adapt to the GPU's thread hierarchy and fully utilize the special structure of the convolution kernel, the liver tissue temperature field data (including spatial dimension H×W and channel dimension in_C) is expanded along the spatial dimension, and the input is flattened along the H×W dimension. Thus, each row of the input matrix corresponds to all data within each input channel, and the data in each batch is expanded into in_C×( On the one hand, the kernel is a two-dimensional matrix in the shape of 1×1; on the other hand, since the kernel size is 1×1, each element of the kernel matrix corresponds to a kernel, and each column corresponds to the tiling of the kernel in each input channel.

[0014] Through the above technical solution, this operation achieves a triple match between data layout, GPU thread scheduling, and convolution characteristics, which facilitates efficient allocation of computing resources by the hardware, maximizes the utilization of thread parallelism, and further reduces latency.

[0015] Optionally, in step S1, in order to fully utilize the GPU's hardware resources, the matrix multiplication is performed along in_C and (... The computation task of the result matrix is ​​divided into several threads for parallel execution, with each thread block arranged into a three-dimensional grid. The first and second dimensions represent the division along two dimensions of the result matrix, and the third dimension represents the division along the batch size. Taking advantage of the fact that the two dimensions of the convolution kernel matrix are equal, the input matrix is ​​divided along (…). The dimension is divided into several square matrices with the same dimensions as the convolution kernel matrix. Each thread block is assigned a computational task of matrix multiplication between square matrices of size len, and is responsible for calculating the temperature field within a specific channel range.

[0016] Through the above technical solution, this operation achieves precise mapping between computing tasks and GPU thread levels and maximizes parallel efficiency, ensuring full utilization of GPU shared memory and thread resources, further improving cache utilization, and facilitating kernel integration with subsequent operations.

[0017] Optionally, in step S1, the threads are arranged into two-dimensional thread blocks, with the first and second dimensions of the thread blocks tiled along the two dimensions of the block result matrix, and the size len is further tiled along (…). The computational workload is divided into slice_len segments based on the out_C dimension and the out_C dimension, so that each thread is responsible for calculating a portion of the computational results for a small sub-region in a batch. Each thread in a thread block reads sub-blocks of the input matrix and convolution kernel matrix from global memory into shared memory using a merging and vectorized memory access method. For the input matrix, each thread block performs vectorized read and write operations on the input matrix in each batch to the shared memory along (…). The square matrix sub-blocks are divided into dimensions. For the convolution kernel matrix, the transpose operation is performed simultaneously when writing to shared memory. When performing matrix multiplication based on vector outer product, data is loaded from shared memory into registers through vectorized memory access and the block calculation results are cached in registers, further reducing the overhead of memory access and data transfer.

[0018] Through the above technical solution, this operation optimizes memory access and computation mode through thread-level design, improves bandwidth utilization, increases computational throughput, solves memory bottlenecks and parallel efficiency issues in convolutional computation, and enhances the performance of convolutional computation on GPUs.

[0019] Optionally, in step S2, a computational fusion strategy is adopted to integrate the convolution computation and its subsequent bias terms with the superposition of the frequency domain transformation operation output into a single kernel function. The output of the frequency domain transformation operation and the output of the convolution computation adopt a unified tiling strategy, and the computation tasks are allocated to each thread block based on a consistent partitioning scheme. The data is cached in shared memory and added to the matrix multiplication results stored in the register. For the superposition operation of the bias terms, it is completed synchronously when the final result is written back to global memory.

[0020] Through the above technical solution, this operation reduces the number of kernel startups by using shared memory and registers in a coordinated manner, allowing GPU computing resources to continuously focus on effective computing, thus making fuller use of the GPU's computing power and improving end-to-end computing efficiency.

[0021] Optionally, in step S3, the thread blocks are arranged in a two-dimensional format, with the frequency domain data along (B*in_C)×( The two-dimensional flattening of the thread block along the dimension of () Tasks are divided along ( ) dimensions, and each thread follows ( ) The task of calculating the coordinate offset of elements and extracting data is performed in parallel along the ( ) dimension, where the storage address of the data is along ( ) The dimensional tiling is jumpy. Based on the truncation size in each dimension, the memory address offset of the retained element is calculated. The non-contiguous truncation data is read from memory according to the calculated offset. Then, based on the truncation size, the address offset after rearrangement is calculated. The truncation data is then merged, concatenated, and transposed according to a predetermined pattern. After calculation, the results are inversely transformed along (…). Tile the same length along the dimension, transpose and put it back into the original position.

[0022] Through the above technical solution, this operation achieves coordinated optimization of memory access efficiency and computational continuity in the high-frequency signal filtering process by structured rearrangement and parallel address mapping of frequency domain data. It solves the problems of low efficiency of step-by-step memory access and mismatch between data layout and computation, and improves the execution efficiency of the high-frequency signal filtering module of the Fourier neural operator.

[0023] Optionally, in step S3, after merging and transposing, the low-frequency signals scattered in the frequency domain are spliced ​​together into a shape of ( A tensor of size () × B × in_C can be viewed as a batch size of ( A two-dimensional matrix of shape B×in_C, after data rearrangement optimization, the serial operation of low-frequency signal block linear transformation is performed by ( )×B×in_C and ( Replace the efficient matrix multiplication of ) × in_C × out_C, where ( The whole is considered as a batch dimension. In each batch, matrix multiplication with shape B×in_C and in_C×out_C is performed. The originally scattered low-frequency signal data is integrated into a matrix form with continuous memory addresses.

[0024] Through the above technical solution, this operation realizes the structured integration and computational conversion of low-frequency signals in the frequency domain, achieves dual optimization of storage regularization and computational parallelization, ensures the accuracy of low-frequency signal extraction, and breaks through the performance limitations of the frequency domain module, effectively improving the end-to-end computational efficiency of the Fourier neural operator.

[0025] After adopting the above technical solution, the beneficial effects of the present invention are: This invention proposes an efficient dynamic simulation analysis method based on Fourier neural operators (FNO-Speed), along with a series of comprehensive solutions to address the inefficient operations of the FNO solver that underutilize hardware. The invention employs two unique optimization methods: a multi-level parallel implicit image-to-column general matrix multiplication (im2col GEMM) optimization strategy and a custom-sized high-frequency signal filtering algorithm. The former replaces pointwise convolution with efficient general matrix multiplication (GEMM) through implicit image-to-column (im2col) and data partitioning strategies, while the latter eliminates the fragmented computation of local linear transformations in the frequency domain. FNO-Speed ​​fully utilizes memory bandwidth to improve computational efficiency, aiming to solve the problems of low computational resource utilization and latency caused by large-scale data memory access computation.

[0026] The purpose of this invention is to address the performance bottleneck of Fourier neural operators and solve the latency caused by low utilization of computing resources. It provides a series of comprehensive solutions to address the inefficient operation of Fourier neural operator solvers that underutilizes hardware, enabling efficient training and inference. This allows for the solution of complex partial differential equations in a short time with high-precision results, providing an efficient and accurate dynamic simulation analysis method for liver cancer lesion dynamics analysis and offering a scientific basis for the formulation of clinical parameterization schemes for liver cancer radiofrequency ablation surgery. Attached Figure Description

[0027] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0028] Figure 1 This is a detailed structural diagram of the FNO architecture (taking Geo-FNO as an example) and the FNO layer; Figure 2This is a schematic diagram illustrating the application of the FNO-Speed ​​optimization strategy to the FNO layer of an FNO model; Figure 3 This is a schematic diagram of a multi-level parallel implicit image-to-column general matrix multiplication process; where (a) represents mapping the input data and convolution kernel to a computation matrix through implicit image-to-column transformation; (b) represents dividing the input matrix and convolution matrix into square matrices of size len and distributing them to several thread blocks; (c) represents dividing the len length into slice_len segments along the (HxW) dimension and the out_C dimension, distributing them to multiple threads, and implementing matrix multiplication based on vector outer products; Figure 4 This is a schematic diagram illustrating data storage reconfiguration optimization, which uses two-dimensional thread blocks to move dispersed low-frequency signals and perform data splicing, splitting, and transposition. Figure 5 This is a diagram illustrating the execution time and speedup ratio of the training epochs for the baseline model and the optimized model under different batch sizes; Figure 6 This is a schematic diagram showing the acceleration effects of different acceleration strategies when the batch size is 16 in a training scenario for solving a two-dimensional plastic deformation problem. Figure 7 This is a schematic diagram illustrating the GPU core execution time and speedup ratio of PyTorch implementation with the multi-level parallel implicit im2col GEMM strategy; Figure 8 This is a schematic diagram illustrating the GPU core execution time and speedup ratio for a PyTorch implementation of a custom-sized high-frequency signal filtering strategy with forward propagation. Figure 9 This is a schematic diagram illustrating the GPU core execution time and speedup ratio for backpropagation of a PyTorch implementation and a custom-sized high-frequency signal filtering strategy. Detailed Implementation

[0029] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0030] This application discloses an efficient dynamic simulation analysis method based on Fourier neural operators.

[0031] Example according to Figures 1 to 9 As shown, an efficient dynamic simulation analysis method based on Fourier neural operators is presented, including the following content.

[0032] I. Neural Operators Neural operator learning is one of the mainstream methods for solving partial differential equations and is the only known model that can guarantee both discrete invariance and global approximation. The neural operator architecture is multi-layered, mapping directly between function spaces on bounded domains through an end-to-end combination of linear integral operators and nonlinear activation functions. Furthermore, it can be discretized using various methods and resolution levels without retraining, thus satisfying discrete invariance.

[0033] Several operator designs have been proposed in existing research, including graph-based neural operators, low-rank neural operators, multi-pole graph-based neural operators, and Fourier neural operators. Among them, the Fourier neural operator network starts from the frequency domain formula of the neural operator and applies an efficient algorithm in the fast transform method, which makes the Fourier neural operator network have good performance and high accuracy in solving partial differential equations, achieving excellent results in solving various problems.

[0034] II. Fourier Neural Operator (FNO) Instead of directly applying the kernel function to the solution domain for integration, FNO considers representing and parameterizing the kernel function in Fourier space. This allows FNO to utilize the Fast Fourier Transform method to reduce the computational complexity of complex kernel integral operators to linear levels.

[0035] Taking the Geo-FNO model as an example, the backbone of the Geo-FNO model (Geometrically Aware Fourier Neural Operator) is the Fourier neural operator. The structure of the Fourier neural operator is as follows: Figure 1 As shown, the input is divided into two channels, upper and lower. At the top, a Fourier transform F is applied, followed by a linear transform R to filter out higher-frequency signals, and then an inverse Fourier transform F⁻¹ is applied. At the bottom, a linear transform W is applied. The upper and lower channels are merged and passed through an activation function to output the layer.

[0036] It can be seen that the most important components of the Fourier neural operator include the frequency domain local transformation in the upper channels and the linear transformation in the lower channels. These two parts also account for a significant amount of computational overhead in the Fourier neural operator, limiting its computational speed and becoming performance bottlenecks in the training and inference processes of the Fourier neural operator network. This invention analyzes the performance of the FNO network in solving different problems, and the device-side operation is the main overhead during the training process.

[0037] III. Dynamic Simulation Model of Liver Cancer Lesions Based on FNO This invention addresses the characteristics of radiofrequency ablation of liver cancer lesions by constructing an FNO-based solver based on dynamic models with different initial and boundary conditions, and employing accelerated optimization methods to achieve accurate and rapid simulation of the heat conduction dynamics of liver cancer lesions at full resolution.

[0038] Assuming that during radiofrequency ablation surgery for liver cancer, when the tissue blood perfusion rate is lower than the baseline value... This indicates that the tissue has been damaged. Let the lesion area be... The goal of the surgery is to minimize the sum of the volumes of the following two parts. Firstly, the volume of cancer cells that failed to be eliminated. Secondly, the volume of normal cells affected by side effects. Therefore, under the parameterized constraints of the ablation needle, it is necessary to solve the Penners biological heat transfer equation that satisfies the corresponding boundary and initial conditions, so that... To reach the minimum.

[0039] Based on the characteristics of radiofrequency ablation surgery for liver cancer, the partial differential equation for heat conduction is expressed by the Penners equation and the biological heat transfer equation as follows: ; in, It is based on the product of liver tissue density and heat-melting ratio given by the apparent heat-melting method. It is the thermal conductivity of the tissue as a function of temperature. , , These are blood density, specific heat capacity, and baseline blood perfusion rate of liver tissue. For blood temperature, The temperature of liver tissue. This is the heat source pulse function generated by a single-needle ablation needle. This represents the temperature gradient of liver tissue.

[0040] Based on the actual heat conduction characteristics of liver tissue, the heat source boundary condition at the point where the ablation needle generates the heat source pulse is the Robin boundary condition. Tissues far from the heat source do not generate heat conduction, and initially, the liver muscle temperature is the same as the blood temperature. The boundary conditions and initial conditions are as follows: ; in, It is a spatial position vector. and These are represented as the heat source boundary and the adiabatic boundary, respectively. and These represent the convective heat transfer coefficient between blood and tissues, and the temperature of circulating fluids in the body, respectively. This represents the directional derivative of temperature along the boundary normal. This indicates the tissue temperature at spatial location r at the output time (time t=0).

[0041] Spatial sampling is performed based on the defined dynamic model (a biological heat conduction model constructed based on the Penners equation). The obtained discrete data points are used to construct training and test sets, respectively, as input data. An FNO-based solver is constructed to realize the dynamic simulation of the ablation process of liver cancer lesions. The optimization method explored in this invention is used to accelerate the simulation analysis of heat conduction dynamics of liver cancer lesions.

[0042] IV. FNO-speed optimization strategy In this section, the invention applies the FNO-speed optimization strategy to the FNO model, and its overall architecture is as follows: Figure 2 As shown, this strategy replaces the frequency domain local transformation and linear transformation in the original Fourier neural operator layer with a method and a customized high-frequency signal filtering scheme, making the computation and memory access operations of the Fourier neural operator network in training and inference more adaptable to the hierarchical structure of the GPU, alleviating the inefficient use of hardware, and achieving acceleration.

[0043] V. Multi-level parallel implicit image-to-column general matrix multiplication (multi-level parallel implicit im2col GEMM) In this section, the invention will discuss the convolution operator in the FNO layer, analyze the performance bottleneck caused by it in the simulation of liver cancer lesion dynamics, and introduce an optimization strategy for implicit image-to-column universal matrix multiplication to address this bottleneck.

[0044] 1. Challenges of memory access redundancy Convolution operations involve a convolution kernel performing a sliding window operation on the input feature map in both horizontal and vertical directions with a preset stride. Local features within each sliding window are calculated using a weighted summation method, and the results are mapped to the corresponding positions in the output feature map. Finally, the complete output feature matrix is ​​generated by point-by-point stacking. Therefore, convolution operations are characterized by sparse connectivity, meaning the size of the convolution window is much smaller than the size of the output matrix. Each output unit is only related to the convolution kernel and a portion of the input feature map elements, and is a 2-D cross-correlation result obtained from the convolution kernel and the convolution window of the input feature map. Due to the characteristics of convolution operations, there is a large amount of repeated access to convolution kernel elements, leading to access redundancy and additional access latency. To improve the memory access efficiency during computation, the im2col algorithm and the hierarchical storage structure of the GPU can be used to enable efficient convolution computation with the help of Generalized Matrix Multiplication (GEMM). The im2col algorithm is one of the main methods for convolution optimization. It transforms the input matrix and convolution kernel into matrices adapted to the convolution calculation order based on the convolution window, and then obtains the final output matrix through matrix calculation, making the convolution calculation applicable to the GEMM optimization method. However, explicit matrix transformations require additional memory allocation and data movement overhead, which is particularly significant in the dynamic simulation analysis of liver cancer lesion models. For example, when simulating the temperature distribution of liver tissue, explicit transformations can lead to a surge in memory usage, even exceeding the GPU's memory capacity, severely impacting training speed. Furthermore, the GPU's hierarchical storage structure offers opportunities for memory access optimization of matrix transformations, utilizing implicit caching to avoid the overhead of explicit transformations. However, simply placing data in the GPU's hierarchical storage structure without organizing data access in a fine-grained manner makes it difficult to fully utilize the bandwidth of each storage structure, resulting in significant bandwidth waste.

[0045] To avoid the overhead of explicit im2col and to fully utilize GPU memory access bandwidth, this invention designs a scheme to implicitly execute im2col. Based on the spatial sampling characteristics of liver tissue and the multi-level cache structure of the GPU, the input data is reasonably divided. By merging memory transactions and vectorized memory access, the memory access bandwidth on the GPU is improved. Operator fusion is used to further optimize the computation and memory access overhead.

[0046] 2. Implicit Image-to-Column (im2col) Algorithm and Partitioning Strategy Convolution computation is transformed into matrix multiplication. Specifically, the input is converted into a large-scale matrix, where each column contains the input elements required for computation at each convolution kernel position during the convolution process. The elements of each convolution kernel are then used to generate another matrix according to their corresponding positions during the convolution process. The final convolution result is obtained by matrix multiplication of the two matrices. The convolution operators in Fourier neural operator layers have a special kernel shape, with a size of 1×1, a stride of 1, and zero padding. They do not need to store repeatedly accessed elements or extra elements, making them very easy to convert into the required matrix without requiring complex coordinate transformations.

[0047] To adapt to the GPU's thread hierarchy and fully utilize the special structure of convolution kernels, such as Figure 3 As shown in (a), this invention expands the liver tissue temperature field data (including spatial dimension H×W and channel dimension in_C) along the spatial dimension and flattens the input along the H×W dimension. Thus, each row of the input matrix corresponds to all data within each input channel, and the data within each batch is expanded into a two-dimensional matrix of shape in_C×(H*W). On the other hand, as... Figure 3 As shown in (a), since the convolution kernel size is 1×1, each element of the convolution kernel matrix corresponds to one convolution kernel, and each column corresponds to the tiling of the convolution kernel within each input channel. Through the above transformation, the convolution calculation is converted into matrix multiplication, with each row of the input matrix and each column of the convolution kernel matrix corresponding to a partial calculation of an output channel. In the dimension of the input matrix, (H*W) is usually much larger than in_C, resulting in a short and wide matrix with a special shape. Since the number of input channels and output channels in the Fourier neural operator layer is equal, the complex convolution calculation of the high-resolution temperature field data of liver tissue is finally transformed into matrix multiplication between a square matrix and a short and wide matrix.

[0048] To fully utilize the GPU's hardware resources, the matrix multiplication task described above is divided along the dimensions in_C and (H*W) and assigned to several threads for parallel execution. Correspondingly, as... Figure 3 As shown in (b), we arrange the thread blocks into a three-dimensional grid, where the first and second dimensions are the division along the two dimensions of the result matrix, and the third dimension represents the division along the batch size. Utilizing the property that the two dimensions of the convolution kernel matrix are equal, the input matrix is ​​divided along the (H*W) dimension into several square matrices identical to the convolution kernel matrix. Each thread block is assigned the computational task of matrix multiplication between square matrices of length len, responsible for calculating the temperature field within a specific channel range. When performing matrix multiplication of shapes M×K and K×N, a larger K can cause threads to occupy the GPU for extended periods, resulting in uneven load distribution. Matrix multiplication between square matrices allows the computational task to be distributed to multiple threads for parallel execution, improving GPU core utilization and avoiding uneven load distribution.

[0049] This invention arranges threads into two-dimensional thread blocks. The first and second dimensions of each thread block are tiled along the two dimensions of the block result matrix. Furthermore, the size `len` is divided into computational loads of slice length `slice_len` along the (H*W) and `out_C` dimensions, respectively. This allows each thread to be responsible for calculating a portion of the computational results for a small sub-region within a batch. This fine-grained partitioning also facilitates the calculation of detailed temperature fields at tumor boundaries.

[0050] 3. Vectorized memory access and matrix multiplication based on vector outer product On GPUs, merging transactions and vectorized memory access can significantly reduce the number of memory access requests and instructions, thereby significantly improving memory bandwidth utilization. In real-time simulation of radiofrequency ablation for liver cancer, efficient memory access is one of the key bottlenecks ensuring the performance of the FNO solver. For example... Figure 3 As shown in (b), according to the tiling and partitioning strategy described above, each thread in the thread block reads sub-blocks of the input matrix and convolution kernel matrix from global memory into shared memory using a merging and vectorized memory access method. For the input matrix, each thread block performs vectorized read-write operations to shared memory on the square matrix sub-blocks of the input matrix divided along the (H*W) dimension of each batch. For the convolution kernel matrix, a transpose operation is performed simultaneously when writing to shared memory, thus ensuring that the merging and vectorized memory access method is maintained when each thread calculates the block matrix multiplication of small sub-regions using the vector outer product method. Figure 3 As shown in (c), when performing matrix multiplication based on vector outer product, data is loaded from shared memory into registers through vectorized memory access, and the results of block calculations are cached in registers, further reducing the overhead of memory access and data transfer. The vector outer product-based scheme can take advantage of the locality of loop traversal, maximize data reuse and minimize row and column load, cache some repeated memory accesses in registers, and avoid meaningless memory accesses.

[0051] 4. Further integration and optimization Operator fusion is typically used in deep learning scenarios where preceding and following operators have data dependencies, and the result of the first operator is quickly consumed by the second. Since there are no cross-thread block data dependencies during computation, each thread block can perform block matrix multiplication operations in parallel. However, subsequent computations mainly involve element-wise operations and depend on the results of convolution calculations, and the output of convolution calculations, as intermediate results, is quickly consumed. This invention employs a computation fusion strategy, integrating the convolution calculation and its subsequent bias terms with the output of the frequency domain transform operation into a single kernel function. The outputs of the frequency domain transform operation and the convolution calculation use a unified tiling strategy, and computational tasks are allocated to each thread block based on a consistent partitioning scheme. Data is cached in shared memory and added to the matrix multiplication results stored in registers. The bias term superposition operation is completed synchronously when the final result is written back to global memory. Operator fusion reduces data migration overhead during computation and avoids additional memory consumption of intermediate computational products, significantly reducing GPU memory requirements, especially in simulating large-size liver cancer lesion models.

[0052] VI. Custom-sized high-frequency signal filtering This invention will discuss feature extraction of frequency domain data in the FNO layer, focusing on the inefficient data migration problem caused by high-frequency signal filtering, and introduce an optimization scheme based on data relocation to solve this bottleneck. Furthermore, a high-frequency signal filtering method that supports custom sizes is designed.

[0053] 1. Challenges of accessing memory with distributed low-frequency signals The FNO layer contains two channels, with the top transform channel selectively retaining low-frequency components and filtering out high-frequency components. Specifically, a Fourier transform is applied to the input data, followed by frequency band filtering in the frequency domain, retaining low-frequency components and performing a linear transform while suppressing high-frequency signals. Finally, the processed frequency domain signal is reconstructed to the time domain through an inverse Fourier transform. In the time domain, filtering out high-frequency signals through Fast Fourier Transform can be represented as "a subtle change in a value at a certain point in time," which often corresponds to detailed information about the signal. Based on this characteristic, filtering out high-frequency signals can be seen as an optimization strategy that balances computational efficiency and accuracy: by selectively discarding high-frequency components, computational speed can be improved within an acceptable range of accuracy loss. In radiofrequency ablation surgery simulations, high-frequency signals correspond to subtle local fluctuations in the tissue temperature field, while low-frequency signals reflect the overall thermal diffusion trend. By selectively suppressing high-frequency components, computational efficiency can be significantly improved while ensuring the accuracy of macroscopic temperature field characteristics, enabling real-time optimization of surgical parameters.

[0054] However, when sampling low-frequency signals, the data is dispersed in spatial coordinates. The sampling process requires discontinuous data extraction across multiple dimensions. For example, in a heat conduction model, multiple discontinuous frequency bands need to be extracted along spatial dimensions. This discontinuous memory access pattern inevitably leads to step-by-step memory access problems, severely limiting the efficiency of memory bandwidth utilization. Furthermore, due to the non-centralized distribution of low-frequency signals in the spatial domain, their linear transformation process requires block extraction from the frequency domain data and serial computation on different frequency components. This approach significantly increases the startup frequency of the computational kernel, resulting in additional overhead.

[0055] The simulation of heat conduction dynamics in liver cancer lesions is time-dependent, and its solution domain contains information in three dimensions: two spatial dimensions (X, Y) and a time dimension T. These factors collectively constrain the overall performance of low-frequency signal processing, necessitating optimization of memory access patterns and computational strategies to improve processing efficiency. To address this performance bottleneck, this invention designs a scheme to adjust the data storage layout and merge small-scale computational tasks, integrating low-frequency signals from multiple corner regions of the liver tissue temperature field data into a single computational task.

[0056] 2. Data Restructuring Considering the dynamic simulation analysis scenario of a liver cancer lesion model, a Fast Fourier Transform (FFT) operation is performed along the spatial and temporal dimensions of the input data. The resulting frequency domain signal is a tensor of shape (B×in_C×X×Y×T), which can be viewed as a batch of size B*in_C. This invention arranges thread blocks in a two-dimensional format. The frequency domain data is unfolded in a two-dimensional flattened manner along the (B*in_C)×(X*Y*T) dimension. The thread blocks are divided into tasks along the (X*Y*T) dimension, and each thread performs the task of calculating the coordinate offset of the elements and extracting data in parallel along the (X*Y*T) dimension. The data storage address jumps along the (X*Y*T) dimension, such as... Figure 4 As shown, the memory address offset of the retained elements is calculated based on the truncation size in each dimension. The non-contiguous truncation data is read from memory based on the calculated offset, and the address offset after rearrangement is calculated again based on the truncation size. The truncation data is then merged, concatenated, and transposed according to a predetermined pattern. This data rearrangement strategy integrates the originally scattered small-scale linear transformations in subsequent operations into large-scale matrix multiplication operations, avoiding cross-step memory accesses in subsequent calculations. After calculation, the results are inversely transformed to flatten the data along the (X*Y*T) dimension to the same length, transpose it, and place it back into its original position.

[0057] After merging and transposing, the low-frequency signals scattered in the frequency domain are spliced ​​together into a tensor of shape (X*Y*T)×B×in_C, which can be regarded as a two-dimensional matrix of batch size (X*Y*T) and shape B×in_C.

[0058] Custom size signal filtering After the aforementioned data rearrangement optimization, the serial operation of low-frequency signal block linear transformation is replaced by efficient matrix multiplication of (X*Y*T)×B×in_C and (X*Y*T)×in_C×out_C. Here, (X*Y*T) is treated as a batch dimension, and matrix multiplication of shape B×in_C and in_C×out_C is performed within each batch. The originally scattered low-frequency signal data is integrated into a matrix form with contiguous memory addresses. This optimization improves the locality of data access, enabling subsequent matrix multiplication operations to achieve a continuous memory access pattern. Furthermore, since the size of (X*Y*T) is often much larger than B×in_C, the transpose operation replaces the computational task with multiple parallel subtasks, each corresponding to a smaller-scale matrix multiplication operation, resulting in higher computational efficiency and accelerating the execution speed of low-frequency signal linear transformation.

[0059] Based on the above optimization strategy, the high-frequency signal filtering process is significantly simplified: users only need to specify the dimension information to be extracted and the corresponding size parameters, without having to manually perform tedious frequency domain data block extraction and calculation operations. This method automatically integrates low-frequency signals according to the user-defined extraction parameters and obtains the transformation result through GGME calculation. The method automatically allocates a zero-initialized memory space of the same size as the complete frequency domain data and maps the transformation result to the new memory region according to its original spatial location, automatically filtering irrelevant high-frequency noise and focusing on the thermodynamic changes in the core lesion region. This custom-sized high-frequency signal filtering method has two advantages: first, it improves the computational efficiency of low-frequency signal mapping by optimizing the data access mode and calculation process; second, it simplifies user operations and lowers the barrier to entry for medical researchers, enabling them to quickly complete simulation analysis for specific lesion parameters without manually performing complex frequency domain block extraction.

[0060] VII. Verification Experiment 1. Experimental Preparation To verify the performance improvement of the proposed FNO-Speed ​​optimization strategy, performance evaluation was conducted on a platform configured with an NVIDIA GeForce RTX 409 and an Intel(R) Core(TM) i7-12700F. The software environment configured for a series of experiments was CUDA 11.5h and PyTorch 1.11.0. This invention uses the publicly available PyTorch implementation model Geo-FNO as a benchmark and encapsulates each operator module in the model to independently test the performance of each operator under different settings.

[0061] This invention uses the solutions to transonic flow around an airfoil, incompressible flow in a pipe, two-dimensional plastic deformation, and climate simulation scenarios as examples. The solution domains for the transonic flow around an airfoil, incompressible flow in a pipe, and climate simulation are two-dimensional, using Geo-FNO 2d and FNO 2d models as the model backbone. Plastic deformation is a time-dependent problem, including a time dimension in addition to the two spatial dimensions; therefore, the Geo-FNO 3d model is used as the model backbone. Specific model configurations are shown in Table 1.

[0062] Table 1. Model configurations for solving transonic flow around an airfoil, incompressible flow in a pipe, plastic deformation, and climate simulation. The transonic flow around an airfoil, the incompressible flow in a pipe, and the two-dimensional plastic deformation problem use the datasets provided in the paper (Fourier neural operators for learnable deformations of partial differential equations over a general geometric domain), while the dataset for climate simulation is provided in the paper (operator learning based on a coupled attention mechanism).

[0063] The transonic flow around an airfoil was solved using a second-order implicit finite volume solver, generating 1000 training data points and 200 test data points. A C-type mesh was used, with a mesh size of approximately (200 × 50) quadrilateral elements. The mesh point locations and the Mach numbers at these mesh points were used as input and output data.

[0064] The incompressible flow problem within a pipe was solved using a Taylor-Hood Q2-Q1 hybrid element implicit finite element solver, generating 1000 training data points and 200 test data points. The mesh size was (129×129), and the horizontal velocities at these mesh points were used as input and output data.

[0065] The plastic deformation problem was solved using the commercial finite element software ABAQUS, which generated 900 training data points and 80 test data points. The data were presented at 20 time steps on a (101×31) structured grid.

[0066] The climate simulation dataset is derived from Physical Science Laboratory meteorological data from 2000 to 2010. It consists of sample pairs of daily surface air temperature and pressure data, with data in a (72×72) grid as input, and uses 1825 training instances and 1825 test instances.

[0067] This invention uses the aforementioned dataset to train the model and compares the benchmark model with the proposed optimization scheme. To ensure data reliability, the average of five experimental results is used for comparison.

[0068] 2. End-to-end performance acceleration comparison In this experiment, the performance of the baseline model and the optimized model is verified by comparing the end-to-end training time for each training epoch. For example... Figure 5 As shown, the baseline represents the benchmark model, FNO-Speed ​​indicates that the optimization strategy proposed in this paper is applied to the model, speedup ratio represents the speedup ratio achieved after adopting the optimization strategy, and (f) and (e) are used to distinguish the running time of the first training round and the average running time of each subsequent training round.

[0069] from Figure 5 It can be seen that the optimized model exhibits a certain speedup effect in both the first training epoch and each subsequent training epoch. The first training epoch takes longer for both the baseline and optimized models because initializing model parameters, initializing the CUDA driver, loading data, and kernels require additional time when the GPU is first invoked for computation. In the 2D scenario, execution time decreases as the batch size increases. This is because the GPU spends a long time waiting for data during preprocessing, and the smaller data size and faster computation speed of 2D problems mean that the gaps between data loading and GPU core execution mask the computation time. Therefore, increasing the batch size increases computation time but reduces gap time, thus shortening the execution time of each training epoch. The optimization strategy merges and integrates multiple computational operations, reducing the waiting gaps in GPU execution and showing significant effects in training with small batches of data. The speedup is more pronounced when the batch size is less than 16, achieving speedups of up to 1.55x, 1.45x, and 1.27x for the first training epoch, and 1.76x, 1.81x, 1.45x, and [missing value] for subsequent training epochs, respectively. As batch size increases, the speedup effect decreases, achieving a speedup of approximately 1.2x for the first training round and approximately 1.3x for subsequent training rounds.

[0070] The modeling problem involves large datasets, making training significantly longer than for two-dimensional problems. While computation time is relatively long, data transfer latency is not negligible, so batch size variations have a relatively small impact on runtime. The optimization strategy's speedup effect directly benefits the processing of larger datasets. As batch size increases, the speedup increases and then stabilizes, achieving approximately 1.3x speedup for the first training epoch and approximately 1.4x speedup for subsequent epochs.

[0071] This invention conducted ablation experiments on the proposed optimization strategy and analyzed the acceleration effect achieved by different optimization strategies, such as... Figure 6 The figure shows the speedup achieved by different optimization strategies when the batch size is 16 during training to solve a two-dimensional plastic deformation problem. Figure 6It can be seen that the optimization strategies proposed in this invention all achieve certain speedup effects. + Implicit Image-to-Column General Matrix Multiplication: Based on the baseline implementation, a multi-level parallel implicit image-to-column general matrix multiplication optimization strategy is added, such as the content of multi-level parallel implicit image-to-column general matrix multiplication. This optimization achieves a 1.16x speedup for the entire model. + Custom-Size High-Frequency Signal Filtering: Based on the above, a custom-size high-frequency signal filtering method is added, such as the content of custom-size high-frequency signal filtering. This optimization achieves a 1.43x speedup for the model.

[0072] This invention compares the training and testing errors of the benchmark model and the proposed optimization method for solving different problems. As shown in Table 2, both the benchmark model and the optimized model, based on Geo-FNO, converged after 500 training iterations. The accuracy of the airfoil and shaping problems decreased slightly, but within an acceptable range. The pipe flow problem showed a smaller convergence error. In climate simulation, the FNO model converged after 1000 training iterations. Using the optimized scheme, the same number of training iterations were completed in less time, while achieving a comparable error magnitude. The optimization scheme proposed in this invention has a negligible impact on model accuracy and requires less training time.

[0073] Table 2 Comparison of accuracy and training time between the baseline and optimization strategies for solving different problem scenarios. 3. Performance Acceleration Comparison of Multi-Level Parallel Implicit Image to Column-Oriented Matrix Multiplication In this experiment, the performance of the multi-level parallel implicit image-to-column general matrix multiplication optimization strategy is verified, and the GPU core runtime is compared under different problem scenarios. The convolution in the benchmark model is implemented by PyTorch using the CuDNN library, and the element-wise addition operation is implemented by PyTorch's naive tensor addition operation.

[0074] Figure 7 This paper compares the CPU core execution time of PyTorch's implementation of convolution and element-wise addition with a multi-level parallel implicit image-to-column general matrix multiplication optimization strategy as the data batch size increases. It can be seen that the execution time of both increases with the data batch size, but the PyTorch implementation consistently outpaces the optimization strategy, indicating that the proposed optimization strategy improves execution speed for different batch sizes, achieving an average speedup of approximately 2x. Implicit image-to-column conversion reduces the additional time overhead from redundant memory accesses, and operator fusion enables parallel execution of convolution and element-wise addition for reading, writing, and computation, further reducing computational overhead. The proposed optimization strategy shortens the execution time of a series of operations implemented in PyTorch to a level similar to the convolution implementation in the CuDNN library.

[0075] Among them, the acceleration effect on the shaping problem is the most significant, with an average speedup of about 2.3. This indicates that the multi-level parallel implicit image-to-column general matrix multiplication optimization method has a better performance improvement effect on processing large-scale data, demonstrating a clear performance advantage in high-dimensional problem solving scenarios. The acceleration effect increases with the batch size and gradually stabilizes, with a speedup of about 2.4 when the batch size is 64. This is because when the batch size is small, the PyTorch implementation does not fully utilize the GPU. The tiling and partitioning strategy proposed in this invention allows threads to fully utilize GPU performance and execute computational operations in parallel and efficiently. When the batch size increases, the computing resources are fully utilized, and performance cannot be further improved, thus gradually stabilizing.

[0076] 4. Performance Acceleration Comparison of Custom-Sized High-Frequency Signal Filtering In this experiment, the performance of the custom-sized high-frequency signal filtering optimization strategy was verified, and the GPU core runtime was compared under different problem scenarios. Figure 7 and Figure 8 As shown, the baseline is the execution time of data slicing and Einstein implementations for low-frequency signal linear transformation and high-frequency signal filtering, implemented in PyTorch.

[0077] Figure 8 and Figure 9 The paper presents a comparison of the execution time of the PyTorch implementation and the optimization strategy proposed in this invention as the data batch size increases during both the forward and backward propagation processes. Figure 8 This describes the speedup effect of the forward propagation process. For 2D problems with small batch sizes, PyTorch implementation takes a long time, especially noticeable in airfoil and pipe flow problems. In small-scale 2D problem scenarios, GPU resources are not fully utilized during execution, resulting in significant idle time and prolonged waiting times, leading to wasted computing resources and low execution efficiency. The optimization strategy reduces waiting time in these gaps by aggregating data through data rearrangement and merging computations. The highest speedup ratios are achieved at a batch size of 16 for airfoil and pipe flow problems, at 9.6x and 5.4x respectively. Climate simulation achieves the highest speedup ratio of 1.7x at a batch size of 32. As the batch size increases, the computation time increases, and the full utilization of computing resources reduces the waiting gaps in GPU execution, resulting in a decrease in the speedup ratio. This explains why the speedup effect is not as significant in shaping problems compared to 2D problems.

[0078] Figure 9Regarding the acceleration effect of the backpropagation process, it can be observed that the optimization strategy proposed in this invention brings significant performance improvements, achieving average speedups of approximately 7x, 12x, 11x, and 20x. The acceleration for the shaping problem is most prominent; while the execution time increases with the batch size, the speedup ratio shows a linear increasing trend. The optimization scheme proposed in this invention simplifies the gradient calculation process during backpropagation in model training. Furthermore, due to data aggregation, the computational overhead of scattered weight gradients is reduced, exhibiting excellent speedup effects, achieving an average speedup of 20x. The speedup effect reaches its maximum at a batch size of 64, with a speedup ratio of approximately 21.94x.

[0079] This invention analyzes the performance bottleneck of the Geo-FNO (Geometrically Aware Fourier Neural Operator). Taking the Geo-FNO model (Geometrically Aware Fourier Neural Operator) as an example, this invention analyzes the optimization potential of pointwise convolution and local linear transformation in the two channels based on the data characteristics and computational properties of the FNO layer, and identifies new performance opportunities for optimization on GPUs. This invention adopts a multi-level parallel implicit image-to-column general matrix multiplication optimization strategy, utilizing data reuse opportunities to implicitly transform convolution calculations into matrix multiplications, and employing a data partitioning strategy to store the matrix in a hierarchical storage structure on the GPU. By utilizing vectorized memory access and matrix multiplication based on vector outer products, the device's memory bandwidth utilization is improved, and operator fusion is performed on subsequent element-wise additions to reduce the overhead of repeated reads and writes. This invention employs a custom-sized high-frequency signal filtering algorithm. This invention proposes a data rearrangement scheme, which performs customized splicing and transposition of low-frequency data in the frequency domain, merges subsequent linear transformation calculations, avoids fragmented step-by-step memory access, fully utilizes memory bandwidth, and reduces the overhead of frequent kernel startup. This invention designs experimental schemes for various cases based on the proposed optimization strategy. This includes performance analysis of FNO-based models during training for different problems, under varying problem dimensions and data sizes, focusing on each module and the end-to-end training process, as well as performance optimization comparisons of the proposed strategies. Through these optimization strategies, FNO-based models achieve stable acceleration across training with different data batch sizes, saving training time and reducing training costs while maintaining accuracy. Average acceleration is achieved in training for both two-dimensional and three-dimensional problems, reducing computation time for liver cancer lesion dynamics simulation analysis while maintaining accuracy, providing faster and more accurate decision support for clinical liver cancer radiofrequency ablation surgery planning.

[0080] This invention analyzes the performance bottleneck of the Geo-FNO model and proposes FNO-Speed, a series of acceleration optimization schemes for the FNO layer to optimize the end-to-end performance of FNO-type models on GPUs during training. Specifically, to address the redundant access caused by pointwise convolution and element-wise addition in linear transformations, this invention proposes a multi-level parallel implicit im2col GEMM optimization strategy. This strategy utilizes data reuse opportunities and vectorized memory access patterns to improve memory access efficiency and bandwidth utilization, achieving a performance improvement of approximately 2.4 times. To address the fragmented step-by-step memory access and frequent kernel startup overhead caused by low-frequency signal linear transformation and high-frequency signal filtering in frequency domain local transformations, this invention designs a custom-sized high-frequency signal filtering method. Through low-frequency data aggregation and block-based computation optimization, it fully utilizes memory bandwidth, reduces memory access latency and kernel startup overhead, achieving an acceleration of approximately 1.1 times in the forward propagation process and approximately 9.5 times and 20 times in the backpropagation process for two-dimensional and three-dimensional problem scenarios, respectively. By applying the aforementioned optimization strategy to the FNO layer of FNO-type models, the end-to-end training time was accelerated by approximately 1.4 times in experiments solving two-dimensional and three-dimensional problems. Experimental results fully demonstrate the effectiveness of the FNO-Speed ​​optimization strategy in improving the performance of the FNO layer, achieving significant acceleration for solving small-scale problems, high-dimensional complex problems, and real-world problem scenarios. The proposed FNO-Speed ​​optimization strategy is universal and can be widely applied to various deep learning models built on FNO, providing a new solution for improving the training efficiency of FNO-type models.

[0081] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the technical solution of the present invention should be included within the protection scope of the present invention.

Claims

1. An efficient dynamic simulation analysis method based on Fourier neural operators, characterized in that, Includes the following steps: S1. Perform multi-level parallel implicit image-to-column general matrix multiplication optimization: take advantage of data reuse opportunities to implicitly transform convolution calculations into matrix multiplications, and improve device memory bandwidth utilization by using vectorized memory access and matrix multiplication based on vector outer product. S2. Perform operator fusion optimization: fuse convolution calculation and subsequent element-wise addition to reduce the overhead of repeated reading and writing and storage of intermediate calculation results; S3. Perform custom-sized high-frequency signal filtering optimization: Adopt a data relocation scheme to perform customized splicing and transposition of low-frequency data in the frequency domain, merge subsequent linear transformation operations, avoid fragmented step memory access, make full use of memory bandwidth, and simplify user operations.

2. The efficient dynamic simulation analysis method based on Fourier neural operators according to claim 1, characterized in that, In step S1, in order to transform the convolution calculation into matrix multiplication, the input is converted into a large-scale matrix, where each column contains the input elements required for calculation at each convolution kernel position during the convolution process. The elements of each convolution kernel are used to generate another matrix according to their corresponding positions during the convolution process. Finally, the result of the convolution is obtained by matrix multiplication of the two matrices.

3. The efficient dynamic simulation analysis method based on Fourier neural operators according to claim 2, characterized in that, In step S1, to adapt to the GPU's thread hierarchy and fully utilize the special structure of the convolution kernel, the liver tissue temperature field data (including spatial dimension H×W and channel dimension in_C) is expanded along the spatial dimension, and the input is flattened along the H×W dimension. Thus, each row of the input matrix corresponds to all data within each input channel, and the data in each batch is expanded into in_C×( On the one hand, the kernel is a two-dimensional matrix in the shape of 1×1; on the other hand, since the kernel size is 1×1, each element of the kernel matrix corresponds to a kernel, and each column corresponds to the tiling of the kernel in each input channel.

4. The efficient dynamic simulation analysis method based on Fourier neural operators according to claim 3, characterized in that, In step S1, in order to make full use of the GPU's hardware resources, the matrix multiplication is performed along in_C and ( The computation task of the result matrix is ​​divided into several threads for parallel execution, with each thread block arranged into a three-dimensional grid. The first and second dimensions represent the division along two dimensions of the result matrix, and the third dimension represents the division along the batch size. Taking advantage of the fact that the two dimensions of the convolution kernel matrix are equal, the input matrix is ​​divided along (…). The dimension is divided into several square matrices with the same dimensions as the convolution kernel matrix. Each thread block is assigned a computational task of matrix multiplication between square matrices of size len, and is responsible for calculating the temperature field within a specific channel range.

5. The efficient dynamic simulation analysis method based on Fourier neural operators according to claim 4, characterized in that, In step S1, the threads are arranged into two-dimensional thread blocks. The first and second dimensions of the thread blocks are tiled along the two dimensions of the block result matrix. Furthermore, the length len is tiled along (…). The computational workload is divided into slice lengths of 0 and out_C dimensions, so that each thread is responsible for calculating a portion of the computation results of a small sub-region in a batch. Each thread in a thread block reads sub-blocks of the input matrix and convolution kernel matrix from global memory into shared memory using a merging and vectorized memory access method. For the input matrix, each thread block performs vectorized read and write operations on the input matrix in each batch to the shared memory along (…). The square matrix sub-blocks are divided into dimensions. For the convolution kernel matrix, the transpose operation is performed simultaneously when writing to shared memory. When performing matrix multiplication based on vector outer product, data is loaded from shared memory into registers through vectorized memory access and the block calculation results are cached in registers, further reducing the overhead of memory access and data transfer.

6. The efficient dynamic simulation analysis method based on Fourier neural operators according to claim 1, characterized in that, In step S2, a computational fusion strategy is adopted to integrate the convolution calculation and its subsequent bias terms with the superposition of the frequency domain transformation operation output into a single kernel function. The output of the frequency domain transformation operation and the output of the convolution calculation adopt a unified tiling strategy, and the computation tasks are allocated to each thread block based on a consistent partitioning scheme. The data is cached in shared memory and added to the matrix multiplication results stored in the register. For the superposition operation of the bias terms, it is completed synchronously when the final result is written back to global memory.

7. The efficient dynamic simulation analysis method based on Fourier neural operators according to claim 1, characterized in that, In step S3, the thread blocks are arranged in a two-dimensional format, with the frequency domain data along (B*in_C)×( The two-dimensional flattening of the thread block along the dimension of () Tasks are divided along ( ) dimensions, and each thread follows ( ) The task of calculating the coordinate offset of elements and extracting data is performed in parallel along the ( ) dimension, where the storage address of the data is along ( ) The dimensional tiling is jumpy. Based on the truncation size in each dimension, the memory address offset of the retained element is calculated. The non-contiguous truncation data is read from memory according to the calculated offset. Then, based on the truncation size, the address offset after rearrangement is calculated. The truncation data is then merged, concatenated, and transposed according to a predetermined pattern. After calculation, the results are inversely transformed along (…). Tile the same length along the dimension, transpose and put it back into the original position.

8. The efficient dynamic simulation analysis method based on Fourier neural operators according to claim 7, characterized in that, In step S3, after merging and transposing, the low-frequency signals scattered in the frequency domain are spliced ​​together into a shape of ( A tensor of size () × B × in_C can be viewed as a batch size of ( A two-dimensional matrix of shape B×in_C, after data rearrangement optimization, the serial operation of low-frequency signal block linear transformation is performed by ( )×B×in_C and ( Replace the efficient matrix multiplication of ) × in_C × out_C, where ( The whole is considered as a batch dimension. In each batch, matrix multiplication with shape B×in_C and in_C×out_C is performed. The originally scattered low-frequency signal data is integrated into a matrix form with continuous memory addresses.

Images