Asynchronous parallel simulation algorithm of large-scale cortical spiking neural network based on GPU

Abstract

The application belongs to the technical field of neural network simulation and analog, and particularly relates to a large-scale cortex pulse neural network asynchronous parallel simulation algorithm based on GPU. The application utilizes the advantages of multi-thread and texture memory of a computing graphics card, combines the general form of a biological brain receiving external stimulation and the general connection mode between neurons in the cortex, designs an asynchronous parallel algorithm framework of GPU and CPU, GPU is responsible for parallel evolution of neuron dynamics equations and block parallel calculation of isotropic connection in a local network, CPU is responsible for processing anisotropic long-range connection, and different neuron dynamics equations and plasticity learning rules can be compatible. Compared with the prior art, the application can effectively improve simulation speed, provides a tool for simulating a biological brain cortex in a single computing node, and is suitable for a single node multi-graphics card and a multi-node multi-graphics card distributed operation model.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of neural network simulation and specifically relates to a large-scale pulse neural network parallel simulation algorithm based on a computing graphics card. BACKGROUND

[0002] Existing large-scale neural network simulation of the cerebral cortex mostly deploys NEST, PyNN, Brian and other pulse neuron network simulation software on a CPU cluster. These software has a unified algorithm framework, that is, the dynamics equation of a single neuron model is solved in parallel on each CPU core thread, and then the transmission of pulse signals between different neurons is realized through parallel programming (MPI) based on message passing. For example, the European Human Brain Project (HBP), the Swiss Blue Brain Project (BBP) and the American Brain Project (BRAIN) all use this technology. Recently, the mouse visual cortex model of the Allen Institute used 5 nodes of 24-core CPUs (Xeon E5-2620 v2, 2.10GHz) for 1 hour of computing time (totaling 600G of memory) to simulate the response of 35,000 integrated firing (LIF) and 10,000 multi-compartment pulse neuron networks for 1 second (Arkhipov, et al., 2018). However, the number of neurons in the primary visual cortex of mammals is actually on the order of 10 8 Therefore, it is still a great challenge for the current CPU parallel algorithm in terms of computing power and memory requirements to complete the simulation of the primary visual cortex even if all simplified LIF neuron models are used.

[0003] On the other hand, the branch processing capability of the CPU core is very strong relative to the core (CUDA-Core) of the computing graphics card (GPU), but in the simulation of the pulse neuron network, only the firing event involves branching, and the average firing rate of the cortical neuron is 1-10 Hz, which cannot fully utilize the branch processing advantage of the CPU. At the same time, the connections between cortical neurons have certain regularity in space, and the linear memory under the CPU framework and the MPI-based parallel processing cannot utilize this topological property of the neural network.

[0004] At present, the number of cores of a single GPU (such as the H100 of NVIDIA Company) has reached 16,896 (1.4GHz), and its block-type processing unit architecture and two-dimensional data reading and built-in interpolation characteristics of the texture memory have great advantages over the linear memory of the CPU when processing local connection matrices and spatially continuous input stimuli. Therefore, how to fully utilize the GPU to design a new algorithm is a breakthrough to reduce the energy consumption of large-scale pulse neural network simulation and lower the application threshold of large-scale pulse neural network. SUMMARY

[0005] The purpose of this invention is to provide an asynchronous parallel simulation algorithm for large-scale cortical spiking neural networks that can reduce simulation energy consumption and improve computing power.

[0006] The asynchronous parallel simulation algorithm for large-scale spiking neural networks provided by this invention is based on the hardware architecture of existing computing graphics cards (GPUs) (such as those from NVIDIA). It leverages the advantages of multi-threading and textured memory in GPUs, combined with the general forms of external stimulation received by the biological brain and the common connection patterns between neurons in the cortex, to design and implement an asynchronous parallel algorithm framework between the GPU and the central processing unit (CPU). The GPU is responsible for the parallel evolution of neuronal dynamics equations and the block-based parallel computation of isotropic connections within the local network, while the CPU handles anisotropic long-range connections. It is compatible with different neuronal dynamics equations and plasticity learning rules, improving the storage and computational efficiency of connection matrices and spiking signal transmission in GPU memory. The specific steps are as follows:

[0007] (1) Preprocessing of external input stimuli, specifically including:

[0008] (1.1) Within the i-th time step (t) i ≤t <t i+1 Input stimuli with one, two, or three-dimensional signals Mapped onto a two-dimensional, three-channel texture memory; since texture memory can have different layers, the third dimension can be used as a discrete dimension, or as data of temporal input stimuli at different times; therefore, texture memory can process basic external input signals, including: vision, hearing (two-dimensional signal), touch (two-dimensional curved surface signal), and limb angle sensation (each joint is a set of two-dimensional signals).

[0009] (1.2) Create a feedforward model of the brain nuclei that process the corresponding senses, with its output pulse signals serving as inputs to cortical neurons:

[0010] (1.2.1) Utilizing the linear interpolation properties of texture memory, spatiotemporal input is continuously sampled within the receptive field corresponding to each nucleus neuron to obtain the sampled input stimulus intensity. Refers to the discretized two-dimensional or three-dimensional coordinate values;

[0011] (1.2.2) Utilizing the Block Reduce parallelism of the GPU (Harris, 2018), the input samples are processed by the convolution kernel of each nucleus neuron. Perform parallel convolution operations:

[0012]

[0013] (1.2.3) The convolution value c is calculated by the nonlinear gain curve function G(c) of the nuclear neuron measured in electrophysiological experiments (arbitrary form, commonly fitted by S-type function), and the change of its discharge rate with time is obtained as the frequency of the time-varying Poisson process;

[0014] (1.2.4) The conversion of the time-varying discharge rate G(c) to the discharge time t sp of the Poisson process is completed; the two-dimensional reading property of the (discrete) texture memory is used to map the corresponding position vector of the nuclear neuron in the external input space to the texture memory, record its two-dimensional serial number, and update the t sp of the corresponding nuclear neuron at the beginning of each time step as the feedforward signal input of the cortical neuron.

[0015] (2) Pre-generation of network connections, specifically including:

[0016] (2.1) Feedforward connection: based on anatomical experimental data, the feedforward connection of nuclear to cortical neurons is directly generated, or the connection is first initialized and a learning rule is added in the feedforward synapse to form a plastic feedforward connection;

[0017] (2.2) Intracortical connection: based on anatomical experimental data, each brain region in the cerebral cortex is generally divided into 6 layers, the intralayer connection of neurons in each layer and the connection mode with other layers in the same plane include both short-range and long-range, specifically:

[0018] (2.2.1) The position of each layer of neurons in the two-dimensional space is determined, and the two-dimensional space is locally blocked, with a total of n local blocks, and the number of neurons in each local block is m, which is a multiple of 32 (usually 1024), to ensure that the parallel scale m x n is equal to the total number of neurons; the size of the GPU in actual operation is determined by the number of hardware units such as cores (CUDA-Core) and SM (streaming multiprocessor);

[0019] (2.2.2) The short-range and long-range connection ranges of neurons in each layer and between layers are determined, the local block center distance and maximum radius are calculated, and the serial numbers of each adjacent local block within the short-range and long-range connection ranges are recorded;

[0020] (2.2.3) Each local block and other adjacent local blocks within the short-range connection range form an adjacency matrix, which is used to store the short-range connection strength in the GPU memory; for the relatively small number of long-range connection strengths, they are stored in the CPU memory in the form of neuron serial number vector and connection strength vector.

[0021] (3) Parallel evolution of single neuron model, pulse signal is stored by block, including:

[0022] (3.1) Confirm the single neuron model used in large-scale neural network, such as Leaky integrate-and-fire, Adaptive exponential integrate-and-fire, multi-compartment model, etc. If multi-compartment model is used (no more than 1024 compartments), the row number of the neuron in the local block can be marked to represent the number of neurons in the local block, and the column number represents the number of compartments of each neuron. At the same time, the numerical method for evolving the model dynamics equation (Runge-Kutta, etc.) needs to be fixed, so that its precision is consistent with the time precision of the pulse firing time algorithm;

[0023] (3.2) Confirm the form of the pulse signal based on the time-varying synaptic conductance or current;

[0024] (3.3) For input layer neurons, first read the pulse firing output of the corresponding kernel front-end neurons according to (2.1), and calculate the input conductance from the kernel part according to (3.2). Otherwise, the dynamics equation in (3.1) can be directly evolved, and it is determined whether the corresponding neuron has a pulse firing in this step, and the firing time is estimated. The estimation time precision needs to be consistent with the numerical method for evolving the dynamics equation;

[0025] (3.4) Record the sequence number and firing time of the firing neurons in the array according to (3.2). Since a time step is usually less than (or equal to) 1 millisecond, the number of pulses in the local block is small, and the total number of pulses is recorded directly using the built-in atomic addition function of CUDA / C++ (NVIDIA GPU programming language), and the sequence number of each neuron in the array is assigned, and the firing time is written in parallel.

[0026] (4) Transmission of pulse signals, specifically:

[0027] (4.1) Short-range connection: calculated by GPU, using the same neuron number as (3.1), first traverse the connection matrix ij (see (2.2.2-3) for the generation of adjacency matrix) between the ith local block and the local blocks within the short-range connection range (including itself), then traverse the read (3.3) in the jth local block in the step of pulse firing neurons k generated by the unit response, and multiply the connection strength of each k column in the adjacency matrix ij in parallel to count the conductance / current of itself (the row label of the adjacency matrix corresponds to the neuron number in the ith local block);

[0028] (4.2) Long-range connection: traverse the neurons with pulse discharge at current step (if calculate the actual propagation speed, add the corresponding delay) and read the corresponding connection strength in the long-range connection vector recorded in (2.2.3) to the conductance / current of the post-synaptic neuron multiplied by the unit response;

[0029] (4.3) Make (4.1) and (4.2) asynchronous, after both are completed, merge the long-range input response of each neuron (CPU memory) into the short-range input response (GPU memory) using the same parallel parameters as (3.1) and update the conductance / current of each neuron at the next time by GPU.

[0030] The overall flow of the asynchronous parallel simulation algorithm of the present application is:

[0031] (5.1) Determine whether the sum of the adjacency matrix memory requirements of all local blocks and their adjacent local blocks can be completely loaded onto the GPU memory at one time;

[0032] (5.2) Start the simulation in time steps;

[0033] (5.2.1) Read or load the external input in real time according to step (1.1);

[0034] (5.2.2) Perform preprocessing of the prefrontal feedforward input in steps (1.2.1)-(1.2.4) and evolve the single neuron model in step (3);

[0035] (5.2.3) If the judgment in step (5.1) is true, perform step (4.1), simultaneously download the neuron discharge time from GPU memory to CPU, perform step (4.2), and finally perform step (4.3); otherwise:

[0036] (5.3.1) According to the number of asynchronous engines of the GPU and the size of the GPU memory, combine the adjacency matrices of every k local blocks into an adjacency matrix local group (there are l local groups, so that kxl=n (number of local blocks), asynchronously upload the first local group from CPU memory to GPU memory, and simultaneously download the pulse discharge time of the neurons in the corresponding first local group;

[0037] Note: The memory size of the adjacency matrix of all local groups is nxp x m 2 (where p is the maximum number of neighbors of the local block), which is much smaller than the size of the traditional adjacency matrix: n 2 xm 2 , solving the problem of insufficient memory of large-scale networks;

[0038] (5.3.2) Upload the 2nd local group to the GPU memory, while the GPU performs step (4.1) in parallel, calculating the short-range connection pulse signal transmission in the 1st local group, and the CPU performs step (4.2) asynchronously to calculate the long-range connection pulse signal transmission in the 1st local group, and uploads the corresponding long-range conductance / current of the neurons;

[0039] (5.3.3) Sequentially upload all local groups and perform the calculations of (4.1) and (4.2) asynchronously, and finally perform step (4.3) to complete a simulation cycle.

[0040] Note: If the GPU model used is older and has smaller memory, it may also encounter the situation that it cannot be completed at a time when performing (5.2.2). The method of transferring memory information to the GPU memory information and GPU parallel calculation asynchronously in (5.2.3) can be referred to.

[0041] Compared with the prior art, the present application can effectively improve the simulation speed, provide a tool for a single computing node to simulate a large-scale pulse neural network of millions of neurons in the biological cerebral cortex, and is suitable for a single node multi-GPU and a multi-node multi-GPU distributed operation model. BRIEF DESCRIPTION OF DRAWINGS

[0042] Figure 1 The algorithm framework diagram of the present application. DETAILED DESCRIPTION

[0043] Embodiment:

[0044] Taking the visual cortex pulse neural network as an example on a single CPU (AMD EPYC 7543 32-Core Processor) single GPU (NVIDIA A40 48G) node:

[0045] (1) Medium-scale pulse neuron network: 32x32 pixel visual input of 120 frames per second, 512 lateral geniculate nucleus (LGN) neurons, 5x1024 primary visual cortex neurons to complete 8000 steps of iteration, step length 0.125 ms, only 6 seconds are needed;

[0046] (2) Large-scale pulse neuron network: 512x512 pixel visual input of 120 frames per second, about 60,000 LGN neurons, 56x1024 primary visual cortex neurons (larger than the scale of Allen Institute 4,5000 neurons), the same 8000 steps of iteration, step length 0.125 ms iteration only takes 2 minutes;

[0047] (3) For a small-scale learnable spiking neural network, 128x128 pixels, 30 frames per second, 256 Lateral Geniculate Neurons, 1x1024 Primary Visual Cortex Neurons, it takes only 54 minutes to complete 1,068,000 steps of iteration (1 ms per step, simulating 17.8 minutes).

[0048] The firing rates of the excited cortical neurons in these examples are all within a reasonable experimental range (5-100 Hz), the image stimuli used are drifting gratings with a temporal frequency of 8 Hz and a spatial frequency of 40 cpd (angle of visual field), and the single neuron models are all Adaptive Exponential Integrate-and-fire models.

[0049] REFERENCES

[0050] Arkhipov, A., Gouwens, N. W., Billeh, Y. N., Gratiy, S., Iyer, R., Wei, Z.,... Koch, C. (2018). Visual physiology of the layer 4 cortical circuit in silico. PLOS Computational Biology, e1006535.

[0051] Harris, M. (2018). Optimizing Parallel Reduction in CUDA. Webinar: NVIDIA corporation.

Claims

1. A GPU-based large-scale cortical spiking neural network asynchronous parallel simulation algorithm, characterized in that, The advantages of multi-thread and texture memory of a computing graphics card are used, combined with the general form of external stimulation accepted by the brain and the common connection mode between neurons in the cortex, to design and implement an algorithm framework of GPU and CPU asynchronous parallelism; wherein, the GPU is responsible for the parallel evolution of the neuron dynamics equation and the block parallel calculation of the isotropic connection in the local network, and the CPU is responsible for processing the anisotropic long-range connection; and it can be compatible with different neuron dynamics equations and plasticity learning rules, to improve the storage and calculation efficiency of the connection matrix and pulse signal transmission in the GPU memory; the specific steps are as follows: (1) preprocessing of external input stimulus, including: (1.1) the first i time step t i ≤ t < t i+1 input stimulus intensity with one, two or three dimensional signals mapped onto a two-dimensional three-channel texture memory; since the texture memory has different layers, the third dimension can be a discrete dimension or data at different time points for time-sequential input stimuli; thus, the texture memory is used to process basic external input signals, including: vision, hearing, touch, limb angle feeling; (1.2) creating a feedforward model of the brain nucleus that processes the corresponding sense, and the output pulse signal of the model is used as the input of the cortical neuron: (1.2.1) using the linear interpolation characteristics of the texture memory, continuously sampling the spatiotemporal input within the receptive field range corresponding to each kernel group neuron to obtain the sampled input stimulus intensity , denotes the discretized two-dimensional or three-dimensional coordinate value; (1.2.2) Block Reduce parallel addition algorithm using GPU, convolution kernel of input sample and each kernel group neuron parallel convolution operation: ， (1.2.3) the convolution value c The non-linear gain curve function of the nucleus neuron obtained by electrophysiological experiment G （c） The change of the discharge rate with time is calculated as the frequency of the time-varying Poisson process. (1.2.4) The completion of the time-varying firing rate by Poisson process G （c） to the discharge moment t sp The corresponding position vector of the kernel neuron in the external input space is mapped to the texture memory by using the two-dimensional reading property of the texture memory, and its two-dimensional serial number is recorded, and the corresponding kernel neuron is updated at the beginning of each time step t sp If there is no discharge, it is 0;​ (2) pre-generation of network connection, specifically including: (2.1) feedforward connection: based on anatomical experimental data, the feedforward connection from the nucleus to the cortical neuron is directly generated, or the connection is initialized first and a learning rule is added in the feedforward synapse to form a plastic feedforward connection; (2.2) intracortical connection: based on anatomical experimental data, each brain region in the cerebral cortex is generally divided into 6 layers, the intralayer connection and the connection mode with other layers in the same plane of each layer of neurons include short-range and long-range, specifically: (2.2.1) determine the position of each neuron in the two-dimensional space in each layer, and locally block the two-dimensional space, assume that it is divided into n local blocks, and the number of neurons in each local block is m, which is a multiple of 32 to ensure that the parallel scale m*n is equal to the total number of neurons; (2.2.2) determine the short-range and long-range connection range of neurons in each layer and between layers, calculate the local block center distance and the maximum radius, and record the serial number of each adjacent local block within the short-range and long-range connection range; (2.2.3) each local block and other adjacent local blocks within the short-range connection range form an adjacency matrix, which is used to store the short-range connection strength in the GPU memory; for the relatively small number of long-range connection strengths, they are stored in the CPU memory in the form of neuron serial number vector and connection strength vector; (3) parallel evolution of single neuron model, pulse signal is stored in blocks, including: (3.1) confirm the single neuron model used in the large-scale neural network; if a multi-compartment model is used, the original neuron row serial number in the local block represents the number of neurons in the local block, and the column serial number represents the number of compartments of each neuron; at the same time, fix the numerical method of evolving the model dynamics equation, so that its precision is consistent with the time precision of the algorithm for determining the pulse firing time; (3.2) confirm the form of pulse signal based on the time-varying form of synaptic conductance or current; (3.3) for input layer neurons, first read the pulse firing output of the corresponding nucleus feedforward neuron according to step (2.1), and calculate the input conductance from the nucleus part according to step (3.2); otherwise, directly evolve the dynamics equation in step (3.1), judge whether the corresponding neuron has pulse firing in this step, and estimate the firing time, the estimation time precision and the numerical method of evolving the dynamics equation are consistent. (3.4) Record the sequence number and discharge time of the spiking neuron into the array: calculate the unit response of each pulse signal at the next time according to step (3.2); since a time step is less than or equal to 1 millisecond, the number of neuron pulses in a single local block is small, so the total pulse number is directly recorded by using the built-in atomic addition function of CUDA, and the sequence number of each neuron discharging at this time step in the array is assigned, and the discharge time is written in parallel; (4) Transmission of pulse signals, specifically including: (4.1) Short-range connections: computed by GPU, using the same neuron numbering as (3.1), first traversing the adjacency matrix of the i local block in question and the local blocks within the short-range connection range ij , then traversing the neurons in the j local block that fired in this step, as recorded in (3.3) k , to multiply their unit responses with the connection strengths in the respective ij column of the adjacency matrix k and add the resulting conductances / currents to their own conductances / currents, with the adjacency matrix row index corresponding to the neuron numbering in the i local block in question; (4.2) Long-range connection: traverse the neurons with pulse discharge at the current step and with long-range connection, read the corresponding connection strength in the long-range connection vector recorded in step (2.2.3), multiply the unit response, and add it to the conductance / current of the post-synaptic neuron; (4.3) Steps (4.1) and (4.2) are performed asynchronously, and after completion, the long-range input response of each neuron is merged into the short-range input response using the same parallel parameters as step (3.1), and the conductance / current of each neuron at the next time is updated in parallel.

2. The GPU-based large-scale cortical spiking neural network asynchronous parallel simulation algorithm according to claim 1, wherein, The overall process is as follows: (5.1) Determine whether the sum of the memory requirements of all local blocks and their adjacent local blocks in the adjacency matrix can be completely loaded into the GPU memory at one time; (5.2) Start the simulation by time step; (5.2.1) Read or load the external input in real time according to step (1.1); (5.2.2) Perform preprocessing of the prefrontal input in steps (1.2.1)-(1.2.4) and evolve the single neuron model in step (3); (5.2.3) If the judgment in step (5.1) is true, perform step (4.1), simultaneously download the neuron discharge time from the GPU memory to the CPU, perform step (4.2), and finally perform step (4.3); otherwise: (5.3.1) According to the number of asynchronous engines of the GPU and the size of the GPU memory, combine the adjacency matrices of every k local blocks into an adjacency matrix local group, asynchronously upload the first local group from the CPU memory to the GPU memory, and simultaneously download the pulse discharge time of the neurons in the corresponding first local group; (5.3.2) Upload the second local group to the GPU memory, and simultaneously perform step (4.1) in parallel to calculate the pulse signal transmission of the short-range connection in the first local group, and perform step (4.2) asynchronously on the CPU to calculate the pulse signal transmission of the long-range connection in the first local group, and upload the corresponding long-range conductance / current of the neurons; (5.3.1) Upload all local groups in sequence and perform asynchronous calculation of (4.1) and (4.2), and finally perform step (4.3) to complete one simulation cycle.

Images