Progress in Data Transmission and Storage for Long Pulse Data Acquisition System

Since the study in the field of fusion has gradually developed toward the long-pulse experiment mode, long-pulse data has gradually become one of the main data types for pulsed experiments in the field of fusion. For long-pulse data, which is a kind of pulse-type data, it will be more difficult to transmit and store than short-pulse data because of its significant characteristics. In addition, in the design of data acquisition and control system (DACS) in fusion field, Experimental Physics and Industrial Control System (EPICS) has now gradually become the main framework of experimental control system to meet the diversity of devices and complexity of subsystems in large experimental system. However, due to the limitation of EPICS, its effectiveness in handling data transmission and storage under high speed data acquisition is not satisfactory. To solve the data transmission and storage under high-speed sampling, this paper proposes a data transmission and storage solution based on TCP/IP protocol and MDSplus database, which is designed with the concept of segmentation, i.e., data generated from experiments longer than 100 seconds are uploaded and stored in a segmented form. Currently, this system has been tested and applied, and the test result shows that the solution is feasible and the overall test system operates stably and reliably.

Introduction
At present, long pulse experiments have become an important study direction in the field of fusion. Through long pulse experiments, experimenters can more accurately analyse the operating state of each subsystem and the behaviour of large experimental devices during the experiments, which will help to study the effects of different experimental parameters on the experimental results and further upgrade the experimental devices. Compared with short-pulse experiments, the duration of long-pulse experiments usually exceeds 100s, so the experimenter usually needs to analyse the data only after the experiments are finished. Therefore, it is necessary to explore a method to synchronize the transmission, storage and visualization of long pulse data. With the study development of the experimental system, the control system of large fusion hosts is now usually developed using EPICS [1], which is first applied to the KSTAR system in Korea [2], and then the ITER organization and Huazhong University of Science and Technology have adopted EPICS to improve their own systems [3][4]. Compared to some data acquisition frameworks, EPICS provides input output controller (IOC) and channel access (CA) protocol, making it possible to support most current bus architectures on the one hand, and on the other hand, those control systems developed based on EPICS can precisely control subsystems through reliable CA protocol. However, nothing is perfect, the characteristics of EPICS make it ineffective in terms of data storage and data transmission for high speed data acquisition. In fact, EPICS does not provide reliable data retrieval and storage, so for data 2 storage processing, an additional database is usually required. For data transmission, EPICS can achieve data transmission at low acquisition rate through channel archiver [5], which creates a predefined buffer and archives the collected data in real time according to the scanning rate of EPICS. However, when the data acquisition rate is higher than 1kHz, EPICS will have a certain delay in scanning the data. Therefore, based on the above data transmission mechanism, it is easy to conclude that it is difficult to achieve better data transmission through channel archiver, which means that the system designed based on EPICS needs to consider data transmission at high speed data acquisition to avoid data loss.
Inspired by the ITER CODAC system, the current DACS has been upgraded based on the EPICS framework [6], and the architecture of it is shown in Figure 1. DACS is composed of three layers, namely the central control system, the server control system, and the sub-control system. The central control system is the human-machine interface (HMI), which includes both remote and local HMIs; the server control layer is mainly composed of a data server and a network server, the data server is used for long-term storage for experimental data, while the network server is used to monitor the data server and provide remote services for users; the sub-control system runs on an industrial personal computer (IPC), which can be composed of multiple IPCs to form a workgroup. The overall workflow of DACS can be briefly described as follows: the experimenter issues a trigger command from the HMI, the relevant configuration data will be obtained by the sub-control system, and at the same time, the sub-control system will start the data acquisition task, and the collected data will be asynchronously uploaded to the data server for long-term storage and visualization on the HMI. This paper focuses on the design of the DACS for data transmission and data storage in the case of high-speed sampling.

The overall design of data transmission and storage
As mentioned above, EPICS is deficient in data transmission and storage. Therefore, MDSplus [7] is used for data storage. MDSplus database is a software that integrates data storage, data processing, and data visualization which is specially designed for long pulse experiments by MIT. It has been widely used for long pulse data processing and storage due to its easy operation and user-friendly interface. As for data transmission, the data transmission process here refers to the real-time transmission of the experimental data from the IOC to the data server for long-term storage. As mentioned earlier, MDSplus provides channel archiver to archive data from the EPICS IOC to MDSplus. However, limited by the working mechanism of the channel archiver, this data transmission method is not suitable for high speed data acquisition. In view of the fact that data transmission under low speed data 3 acquisition can be completed by channel archiver, this paper only discusses data transmission under high speed data acquisition. The collected data is finally transferred to MDSplus in the data server for long-term storage for later analysis and utilization by the researchers. The process of data transmission and storage under high speed data acquisition is shown in Figure 2. In the design and implementation about the scheme, Data acquisition tasks at high sampling speeds are performed by IPC and then the collected data will be scanned by EPICS IOC and read in the process variable (PV). In order to reduce the network transmission load and save storage space, the raw experimental data will be compressed by a data compression algorithm. After that， the compressed data will be integrated into a defined and standardized data structure, which is then transmitted to the data server via TCP/IP based data transmission structure [8]. The packet will then be parsed and the data will be stored in MDSplus using data segmentation techniques. When the experimenter needs to access the experimental data, it can be obtained based on the specific sampling time and channel number of a particular experiment. Finally, the waveform of the raw data will be visualized on the HMI.
It is worth noting that the storage format needs to be specified in advance. Data obtained by different data acquisition channels in the same experiment number will be stored in different data files to prevent data confusion; while experimental data from different experiment numbers but under the same channel should be stored by experiment number to prevent data overwriting. In order to easily display data waveforms to the user, the data segmentation requirements must be met in a specific data format. The requirements of this system under high-speed sampling can be summarized as follows.
 Implementation of communication between the IPC and the data server for data transmission and reception at high speed data acquisition.  Correct parsing of the data structure to facilitate access to the original data and prevent data loss.  Long pulse data segmentation processing.  Experimental data segmentation storage.  Accessing segmented data from the MDSplus database according to pre-set parameters.  A clear, simple and user-friendly data visualization interface.

High speed data acquisition
Given that the overall framework of EPICS is implemented through C and C++, the data acquisition program is also implemented using the C programming language in order to maximize compatibility between the data acquisition program and EPICS. In addition, since EPICS provides a variety of PVs for implementing different functions, we establish a generic device model based on 'ai (analog input)', 'bo (binary output)', and 'calc (calculation)' in order to enhance the compatibility of the system with different DAQ cards. The 'ai' is used to read the collected data, the 'bo' is used to control the experimental devices, and the 'calc' is used to monitor the status of those devices. This device model unifies the abstract programming standard of experimental devices, especially for DAQ cards, so it provides great convenience for the subsequent maintenance, update and expansion of DACS. Before the data acquisition task starts, the experimenter needs to configure the experimental devices in detail through HMI, such as the number of acquisition channels, the selection of voltage range, the setting of acquisition frequency, and the acquisition time. At the same time, the program will retrieve the specified directory to determine the current experiment number to prevent the new experiment data from overwriting the original experiment data. When the experimental system is on standby, the IPC will be in a waiting state, ready to receive the 'START' command from the central control system, after which the data acquisition task will be started with the obtained configuration. Considering the need for real-time access in long pulse experiments, the steps of data acquisition, data transmission and data segment storage will be performed synchronously. Through the high speed data acquisition program, the data collected by the DAQ devices will be continuously scanned into the EPICS IOC, and the next step is to transfer and store the collected data for further processing.
In this study, PCI-9112 DAQ card is used on the IPC for data acquisition [9]. The card has 12-bit sampling accuracy and 16 channels, the highest sampling rate can reach 110 kHz for each independent channel. In addition, the minimum hardware and software environment required for the overall system operation should meet the following requirements: CPU performance of 2.0 GHz or above; memory performance of 4 GB or above; hard disk storage capacity of at least 500 GB. The software platform is a C/C++ programming environment based on the Linux operating system. Moreover, MDSplus database and Linux shell programming language are required for the successful implementation of the system.

Data transmission at high speed sampling
Once data acquisition is finished, one thing that needs to be considered is data transmission. Data acquisition at high speed sampling temporarily stores the collected data in a file named after the channel number that is used for data acquisition in this experiment. In order to reduce network transmission load and save data storage space, this data file will be compressed at this time using a data compression algorithm. The LZO (Lempel-Ziv-Oberhumer) data compression algorithm [10] is chosen for data compression after considering the support for data segmentation and the rate of compression and decompression and space size. After testing, the size of original data file processed by LZO can be reduced by 30%.
Lossless transmission and reception of raw data is crucial, which is why we adopt the TCP/IP protocol to achieve reliable data transmission under high-speed sampling. The data transmission structure based on TCP/IP protocol is shown in Figure 3. As can be seen from the figure, the data structure contains three parts, which are TCP header, data header and compressed data file. The TCP header includes the IP address of the IPC host, the IP address of the data server, the size of this segment of the data structure and CRC checksum. In this design, the IP of IPC is used to mark the source address of data, and the IP of the data server is used to indicate the address to which this data will be sent, and finally, the CRC is used to verify the correctness of the packet after transmission. The data header mainly includes some configuration information about the acquisition of this data, including the acquisition channel number, signal name, sampling frequency, sampling duration, sampling accuracy and voltage range. This information is used to accurately parse the data structure in the data parsing phase and store the data in the correct node in the MDSplus tree. And the last part is the compressed data file. After comparison, the collected data can be transmitted correctly in this way.  Figure 3. The data transmission structure based on TCP/IP Protocol

MDSplus-based data segmentation storage
Now, the collected data has been transmitted to the data server through the above data structure. As mentioned before, the collected data needs to be compressed by LZO, and the original data file will be changed to a lzo file after being compressed. Since MDSplus does not support LZO files for storage, in order to ensure that the data files can be stored successfully, it is necessary to perform decompression operation after the data has been transmitted to the data server. As for decompression, LZO's library provides decompression function that can be called directly at the time of use. After decompression, the compressed file will be converted to a txt file and the txt file is the one that will be stored in MDSplus. The next step is data segmentation and storing the segmented data into the MDSplus. MDSplus stores all experimental data with a tree structure [11]. On the one hand, it unifies the data storage format; on the other hand, it can clearly express the logical relationship between tree nodes. There are two types of trees in the MDSplus: the model tree and the pulse tree. The model tree mainly stores the structure information of the tree, such as the node name and other specific tree structures. The pulse tree will be built after the experimental data was stored in the model tree. After this, if the structure of the model tree is changed again afterwards, the logical structure of the created pulse tree will not change again [12]. In other words, the change of the logical structure of the model tree does not affect the structure of the pulse tree. The relationship between them is shown in Figure 4.  the structure of the model tree, which can save a lot of development time. In this paper, in order to better present the data storage structure of the MDSplus after data transmission, a four-level structure of the MDSplus tree was created in the test server, consisting of the experiment number, device number, channel number and data. The final MDSplus tree is shown in Figure 5.
In this MDSplus tree, 'EXPNUM X' represents the experiment number, 'DEVICE X' represents the device number, 'CH X' represents the channel number and the 'DATA X' represents the segmented data file. As mentioned before, the DAQ card has 16 data acquisition channels, so, there are 16 one-to-one corresponding subtree nodes named after the channel number of the DAQ card in this tree structure. In this design, the model tree of this MDSplus consists of four branches, including parameter, timing, raw and process. The parameter branch stores some constant data, such as amplification, data type, etc. The timing branch is used to store the trigger time, data acquisition time and the frequency of this experiment. The raw branch stores the raw experimental data. The processing branch stores the corresponding signal information of the transmitted packets after processing. As mentioned before, the pulse tree will be built after store the data into the model tree. Therefore, if the experiment starts, the data server needs to stay connected to each IPC all the time when the data server receives a data packet from the IPC, the server program will first parse the packet and write the data to the corresponding MDSplus tree node according to the parsed result. In fact, the parsed data is saved in three files according to the file type in the pulse tree, including tree file, data file and characteristic file. The tree file stores information about the structure of the MDSplus tree; the data file stores the real data. In theory, all data from a long pulse experiment will be written to the node of the pulse tree, while the actual data will only be stored in the data file of the tree node, and some information about this experiment, such as the offset address, will be stored in the access log of the node. MDSplus only maps the '*.tree' file to memory [13], and as for the actual data, it is stored in the hard disk, in other words, the data and the tree information are stored separately, which has the advantage of greatly increasing the data transfer speed.

Data segmentation storage and accessing.
In the specific MDSplus data storage format, segmented storage is used to store data. The core design concept of segmented storage is the "data segment", by which MDSplus can achieve continuous data reads and writes. The key technology of data segmentation is the segmented reading and segmented writing of the collected data. Segmented write means that each segment data is appended to the data node, instead of overwriting the original data, as storage space allows; segmented read means that each segment can be read from MDSplus in a separate form, without being read only with the overall data. By reading and writing in segments, the collected data can be stored in MDSplus in the form of segments without waiting for the end of the experiment. In addition, the experimenter can read the stored experimental data in time during the experiment. The principle of this technique is that when retrieving data from segmented records, MDSplus will select the corresponding segment based on the timestamp, then read all segments that match the retrieval time and join these segments to finally construct the final result that satisfies the retrieval criteria and return it to the user. Through data segmentation storage, the experimenter can read data from MDSplus that is in the experimental period but has been collected and stored. The associated program can then visualize the target data on HMI. The process of accessing data from MDSplus is shown in Figure 6.

Testing
In order to verify the integrity of the experimental data and the reliability of the segmented storage after processing by this system at high speed sampling, after the segmented experimental data is transmitted and stored in MDSplus, NBWave [14] is used for accessing the same experimental data in segment and in whole respectively, and the segmented data is compared with the original data to check whether the data is correctly stored in MDSplus. The test experiment is a long pulse experiment of up to 100 seconds with a square wave and 100 kHz data sampling frequency. In this test, 10 seconds is chosen as the segmentation basis. Due to the high speed data acquisition, for better display of the comparison results, a segment of data is taken here, whose time length is 10s, and then it is divided into 10 segments on the basis of 1 second. The raw data waveform for 10 seconds under this configuration is shown in Figure 7. Then the 2nd to 3rd seconds data stored in MDSplus is randomly selected and then visualized in the other data visualization window and compared with the original data waveform. The result of comparison between them is shown in Figure 8. In this figure, above the middle red line is the original data waveform, and below is the data of the 2nd to the 3th seconds. It can be seen that the upper part and the lower part are exactly the same in terms of the y-axis value that

Conclusion
This paper presents the progress of data transmission and storage for long pulse data acquisition system. It proposes a scheme to cope with high speed data acquisition based on EPICS and designs a reasonable data transmission structure based on TCP/IP protocol, and also proposes a segmented storage scheme based on MDSplus for the data storage in long-pulse experimental mode. It integrates many tools, including EPICS system, MDSplus database, LZO data compression algorithm, TCP/IP network transmission protocol and data segmentation processing technology. The application of the program makes it possible to access real-time data for high speed data acquisition in long-pulse experimental mode. After testing, the results show that the method is effective and fully meets the requirements of practical applications.


Remote Recording of Emotional and Activity Data: A Methodological Study

The impact of physical exercise on emotional well-being is one of the most important factors that drive sustained physical activity engagement. It is also one of the least studied topics on account of the rather elaborated setups required to quantify emotional expression during exercise. The wide adoption of at-home, physical exercise solutions has compounded this problem due to the secluded nature of these activities. We propose here a new methodology that would allow for mass emotional expression and physical activity data gathering using publicly available sources such as at-home exercise videos. We have shown that the methodology is robust enough to extract high resolution, reliable data from home videos posted on popular video share sites. The source-code and instructions for practical use are published online such that researchers can access data pertinent to emotional response during exertion in real world situations with applicability for prospective and retrospective studies.

Introduction
During the last decade, a vast number of physical activity quantification, data recording, and feedback devices and applications were developed and commercialized by virtually every major technology company (Shin et al., 2019). The explicit and implicit goals of most of these devices/methodologies is to increase engagement in physical activity and induce healthy behaviors in an otherwise sedentary population (Mercer et al., 2016). Despite the ubiquity of physical activity data availability, the actual number of people engaged in regular physical exercise has not changed significantly. This leads to the inevitable conclusion that biofeedback in the form of raw or processed physiological and/or biomechanical data does not seem to have the desired impact on activity related behaviors. A number of reasons were proposed for this rather unexpected lack of effect. The most preeminent one is the lack of emotional impact of the biofeedback data (Coen et al., 2020) and, more broadly, our limited understanding of physical exercise on the acute and longterm emotional response (Wenzler et al., 2016).
The ability to quantify the impact of physical activity under various conditions on the emotional state of the athlete would open new possibilities to manipulate the environmental, exercise, and feedback conditions in order to increase physical activity engagement. However, emotional response to exercise has been notoriously difficult to quantify and, as a consequence, very little information is available in the scientific literature (Khanal et al., 2018). The primary reason for this lack of information is the fact that the methodologies used to quantify emotional states are rarely ever compatible with intense physical activity. Furthermore, in the rare occasions where these types of measurements were attempted, the elaborate setups required to collect the data have severely limited the number of subjects that can be tested. Recently, however, a number 3 of data rich sources have become publicly available in the form of recorded (streamed) physical exercise sessions (McGloin et al., 2020). These recordings contain information about physical activity as measured by various instruments used by the athletes as well as a number of features that are related to the emotional state of the user. Thus, should one be able to extract the information pertinent to physical activity and emotional expression, our understanding of the impact of physical activity on emotions (and vice-versa) would be greatly improved. Given that this data is publicly available and does not require specialized, expensive, and time-consuming setups, this approach opens the possibility for scientists from fields unrelated to physical exercise to test various hypotheses related to the relation between physical activity, emotions, behaviors, and overall engagement in physical activity.
It is therefore the aim of this study to create a set of tools that would allow one to extract all the information available in exercise video streams and make it available for further processing.
In particular, we focused here on videos involving a stationary bike but the methodologies can be expanded/extrapolated to other forms of exercise that have similar setups.

Ffmpeg Pre-Processing + Facial Expression Analysis
This study utilized a cohort of 10 publicly available Twitch streams relating to Zwift indoor biking competitions. Videos were selected based on whether the athlete's face and heart rate and power readings were available. Utilizing UnTwitch's online platform, videos were downloaded at 1080p (with the exception of [FR] Workout Time whose highest video quality setting was at 900p). Some videos were split into two separate parts, as UnTwitch can only download 60 minute portions of a video at a time and many videos were above one hour in length.
Each of the videos were subject to pre-processing using version 4.2.2 of FFmpeg. A python script was constructed to easily implement some of FFmpeg's tools as it relates to the study. For any given video, each of its parts were appropriately trimmed to not only capture just the race but also ensure that each frame was entirely unique (not duplicated between parts). Both parts of the video were then set to have the same fps value (4x4, FTP Builder, and GLS were set to 60 fps, others 30 fps). The videos were concatenated together to produce the entire full-length video. For facial expression analysis, the resulting video was cropped to feature only the streaming inset (the athlete).
Emotional response data was extracted from the videos using the iMotions (Version 8.0) biometric software [iMotions]. iMotions uses the Affectiva developed algorithms to derive seven measures of emotion: joy, anger, sadness, contempt, fear, surprise, and disgust.

Tesseract Pipeline
In order to extract the information pertinent to the biomechanics and physiology of physical exercise from the available video streams, an open source Optical Character Recognition (OCR) software package was applied on individual frames.
Specifically, each video was segmented into its constituent frames at a rate of one frame per second utilizing FFmpeg. The resulting images were scanned for power and heart rate data using Tesseract OCR (version 4.1.1). Before running Tesseract, each image was subject to preprocessing in order to increase the accuracy of the open-source software. A rudimentary pre-processing pipeline was constructed in Python, involving the following image modifications (Fig.   1): • A. Crop + Resize: Each image was cropped to only feature the corresponding metric of interest (heart rate or power). Every resulting cropped image was resized to the dimensions 90 x 44 at 300 dpi.
• B. Sharpen: Using PIL (Python's Imaging Library), images were sharpened by a factor of four, which removed some of the blur present in the original image.
• C. Binarize: The images were then subject to binarization followed by a slight Gaussian blur, which ensured that the final image was in the best possible format before applying Tesseract OCR.

Post-Collection Data Analysis
In order to remove outliers (or incorrectly predicted samples) from the resulting Tesseract dataset, a z-score analysis was implemented followed by simple thresholding. Z-score values that were greater than 3 were removed from the dataset and excluded from analysis. To account for incorrect heart rate and power readings that were still present after thresholding, an exponential moving average function was computed across the dataset and plotted as a black line (Fig. 2C).

6
The relationship between power, heart rate, and emotional response was measured every minute utilizing the Pearson correlation coefficient ( Fig. 2A, 2D). Furthermore, an athlete's emotional tradeoff was plotted as a function of time (Fig. 2B).

Results
The data mining pipeline was constructed in Python, and it allows one to successfully extract information from an otherwise largely untapped dataset (publicly available workout videos).
Through Optical Character Recognition (OCR), information concerning an athlete's heart rate and energy exerted during the workout can be used to gauge physical performance (Fig. 2).

Discussion
The delimiting lines between physical activity, well-being, and mental health were blurred during the last decade(s) to the point that tools, methodologies, and scientific terminologies are now being used interchangeably between these topics. Furthermore, a dramatic transition towards individual physical exercise performed in secluded environments (e.g. home exercise) has been observed during the last few years. The global pandemic and associated social interventions have recently accelerated this phenomenon (Tison et al., 2020), (Woods et al., 2020). As a result, the emotional landscape of sports and exercise has been significantly altered (Alradhawi et al., 2020).
However, these effects are virtually impossible to study and quantify with the field's current tools and methodologies, given the seclusive nature of these activities.
In this context, the methods proposed here offer a unique opportunity to study the impact of sport and exercise on the physical and emotional well-being of individuals exercising in a natural environment rather than an artificial (lab) setup. Furthermore, it allows one to access thousands of hours of emotional, physiological, and biomechanical data with minimal effort, a feat that even the best funded of labs would find difficult to match. Not only does this method allow one to extract information from pre-recorded (public) sources, but it also offers a new tool to researchers who want to create remote recording studies that look simultaneously at emotional expression and physical exercise.
In conclusion, the methodology proposed here allows researchers to quantify an elusive yet very important metric: emotional response to physical exercise in a natural exercise environment.
Coupled together, physical and emotional performance offer a wider array of metrics that can be 8 used to explore the implications and effects of different workout regiments. They also provide a convenient way to measure an athlete's "performance success factors" and "emotional tradeoff" over time while providing information concerning physical activity stickiness and effort dosage.
Given that the emotional response to physical exercise seems to be one of the main drivers of sustained physical activity engagement, a number of behavioral change techniques have made use of this subtle yet powerful motivator to increase the duration and frequency of physical exercise.
This new methodology should allow non-expert users to also collect significant amounts of data pertinent to the emotional response to exercise and use it for their respective fields of study (e.g. behavioral sciences).


extfieldsofstudy:
Computer Science
Physics
Mathematics

text:
Power-law escape rate of SGD

Stochastic gradient descent (SGD) undergoes complicated multiplicative noise for the mean-square loss. We use this property of SGD noise to derive a stochastic differential equation (SDE) with simpler additive noise by performing a random time change. Using this formalism, we show that the log loss barrier $\Delta\log L=\log[L(\theta^s)/L(\theta^*)]$ between a local minimum $\theta^*$ and a saddle $\theta^s$ determines the escape rate of SGD from the local minimum, contrary to the previous results borrowing from physics that the linear loss barrier $\Delta L=L(\theta^s)-L(\theta^*)$ decides the escape rate. Our escape-rate formula strongly depends on the typical magnitude $h^*$ and the number $n$ of the outlier eigenvalues of the Hessian. This result explains an empirical fact that SGD prefers flat minima with low effective dimensions, giving an insight into implicit biases of SGD.

Introduction
Deep learning has achieved breakthroughs in various applications in artificial intelligence such as image classification [1,2], speech recognition [3], natural language processing [4], and natural sciences [5][6][7]. Such unparalleled success of deep learning hinges crucially on stochastic gradient descent (SGD) and its variants as an efficient training algorithm.
Although the loss landscape is highly nonconvex, the SGD often succeeds in finding a global minimum. It has been argued that the SGD noise plays a key role in escaping from local minima [8][9][10][11][12][13][14]. It has also been suggested that SGD has an implicit bias that is beneficial for generalization. For example, SGD may help the network to find flat minima, which are considered to imply good generalization [9,15,16], although it is also pointed out that flatness alone does not ensure good generalization [17]. How and why does SGD help the network escape from bad local minima and find flat minima? These questions have been addressed in several works, and it is now recognized that the SGD noise strength and structure significantly affect the efficiency of escape from local minima. Our work follows this line of research and adds new theoretical perspectives.
In physics and chemistry, escape from a local minimum of the energy landscape due to thermal noise at temperature T is a thoroughly studied fundamental problem [18,19]. When the energy barrier is given by ∆E, the escape rate is proportional to e −∆E/T , which is known as the Arrhenius law. By analogy, in machine learning, escape from a local minimum of the loss function is considered to be determined by the loss barrier height ∆L = L(θ s ) − L(θ * ), where L(θ) denotes the loss function at the network parameters θ, θ * stands for a local minimum of L(θ), and θ s denotes a saddle point that separates θ * from other minima. Indeed, if we assume that the SGD noise is uniform and isotropic, which is often assumed in machine-learning literature [8], the escape rate is proportional to e −∆L/D , where D denotes the SGD noise strength.
In this paper, we show that inhomogeneity of the SGD noise strength in the parameter space brings about drastic modification for the mean-square loss. We show that the escape rate is determined by the logarithmic loss barrier height ∆ log L = log L(θ s ) − log L(θ * ) = log[L(θ s )/L(θ * )]. In other words, the escape rate is determined not by the difference but by the ratio of L(θ s ) and L(θ * ). This result means that even if the loss barrier height ∆L is the same, minima with smaller values of L(θ * ) are more stable.
Moreover, given the fact that the eigenvalue spectrum of the Hessian at a minimum consists of a bulk of almost zero eigenvalues and outliers [20,21], our escape-rate formula implies that SGD prefers flat minima with a low effective dimension, where the effective dimension is defined as the number of outliers [22] and flatness is measured by a typical magnitude of outlier eigenvalues [15]. The previous theories [8,9,[11][12][13][14] have also successfully explained the fact that SGD prefers flat minima, but not the preference of small effective dimensions. The logarithmic loss naturally explains the latter, and sheds light on implicit biases of SGD.
Main contributions: We obtain the following main results: • We derive an equation for approximating the SGD noise in Eq. (6). Remarkably, the SGD noise strength in the mean-square loss is shown to be proportional to the loss function, which is experimentally verified in section 5.2. A key ingredient in deriving Eq. (6) is the decoupling approximation given in Eq. (7). This is a novel approximate method introduced in our analysis, and hence we experimentally verify it in section 5.1.
• We derive a novel stochastic differential equation (SDE) in Eq. (14) via a random time change introduced in Eq. (13). Although the original SDE (4) has a multiplicative noise, the transformed SDE (14) has a simple additive noise with the gradient of the logarithmic loss. This shows the relevance of the logarithmic loss landscape for understanding SGD.
• We derive a novel form of SGD escape rate from a local minimum in Eq. (15). Remarkably, the escape rate takes the power-law form with respect to the ratio between L(θ * ) and L(θ s ). In section 5.3, we experimentally test the validity of this result for a linear regression and a neural network.
• We show that the escape rate of SGD crucially depends on the flatness and the effective dimension, which implies that SGD has implicit biases towards flat minima with low effective dimension. We also show in Eq. (16) that a local minimum with an effective dimension n greater than a certain critical value n c becomes unstable.
Related works: The role of the SGD noise structure has been discussed in some previous works [11][12][13][14]23]. It was pointed out that the anisotropic nature of the SGD noise is important: the SGD noise covariance matrix is aligned with the Hessian of the loss function, which is beneficial for escape from sharp minima [11,13,14]. These previous works, however, do not take the parameter dependence of the SGD noise strength into account, and consequently, escape rates derived there depend exponentially on the loss barrier height, which differs from our formula.
Compared with the anisotropy of the SGD noise, the inhomogeneity of the SGD noise strength has been less explored. In [12,23], the SGD dynamics under a state-dependent noise is discussed. However, in these previous works, the connection between the noise strength and the loss function was not theoretically established, and the logarithmic loss landscape was not discussed. The instability due to large effective dimensions was also not shown. Another recent work [24] observed that the noise is proportional to the loss for specific simple models. In our paper, such a result is derived for more generic models. Gürbüzbalaban et al. [25] showed that SGD will converge to a heavy-tailed stationary distribution due to a multiplicative nature of the SGD noise in a simple linear regression problem. Our paper strengthens this result: we argue that such a heavy-tailed distribution generically appears for the mean-square loss.

Setup
We consider supervised learning. Let D = {(x (µ) , y (µ) ) : µ = 1, 2, . . . , N } be the training dataset, where x (µ) ∈ R d denotes a data vector and y (µ) ∈ R be its label. The network output for a given input x is denoted by f (θ, x) ∈ R, where θ ∈ R P stands for a set of trainable parameters with P being the number of trainable parameters. In this work, we focus on the mean-square loss The training proceeds through optimization of L(θ). In most machine-learning applications, the optimization is done via SGD or its variants. In SGD, the parameter θ k+1 at the time step k + 1 is determined by where η > 0 is the learning rate, B k ⊂ {1, 2, . . . , N } with |B k | = B is a mini-batch used at the kth time step, and L B k denotes the mini-batch loss. Since the training dataset D is randomly divided into mini-batches, the dynamics defined by Eq. (2) is stochastic. When B = N , the full training data samples are used for every iteration. In this case, the dynamics is deterministic and called gradient descent (GD). SGD is interpreted as GD with stochastic noise. By introducing the SGD noise Obviously, ξ k = 0, where the brackets denote the average over possible choices of minibatches. The noise covariance matrix is defined as Σ(θ k ) := ξ k ξ T k . The covariance structure of the SGD noise is important in analyzing the SGD dynamics, which will be discussed in section 3.1.

Stochastic differential equation for SGD
When the parameter update for each iteration is small, which is typically the case when the learning rate η is sufficiently small, the continuous-time approximation can be used [26,27]. By introducing a continuous time variable t ∈ R and regarding η as an infinitesimal time step dt, we have a SDE where dW t ∼ N (0, I P dt) with I n being the n-by-n identity matrix, and the multiplicative noise ηΣ(θ t ) · dW t is interpreted as Itô since the noise ξ k in Eq. (3) should not depend on θ k+1 . Throughout this work, we consider the continuous-time approximation (4) with Gaussian noise. In machine learning, the gradient Langevin dynamics (GLD) is also considered, in which the isotropic and uniform Gaussian noise is injected into the GD as where D > 0 corresponds to the noise strength (it is also called the diffusion coefficient) [11,28,29]. The GLD yields an escape rate proportional to e −∆L/D [18,19]. We will see in section 4 that the SGD noise structure, which is characterized by Σ(θ), drastically alters the escape rate from a local minimum.

Dynamics of SGD
The fluctuation of SGD and the local landscape are the two most important factors that influence the escaping behavior of learning. In sections 3.1 and 3.2, we establish the form of the noise of SGD and show that under this type of noise, the relevant landscape for consideration should be log L instead of L. We then use these two results to derive the escape rate of SGD in section 4.

Structure of the SGD noise covariance
The SGD noise covariance matrix Σ(θ) significantly affects the dynamics [8,11,27,30]. In this section, under some approximations, we derive the following expression of Σ(θ) for the mean-square loss near a local minimum θ * : where H(θ) = ∇ 2 L(θ) is the Hessian. It should be noted that ∇L(θ * ) = 0 and H(θ * ) is positive semidefinite at any local minima θ * . We give a derivation below and the list of the approximations and their justifications in Appendix A. In particular, the following decoupling approximation is a key assumption in the derivation.
This approximation seems justified for large networks in which ∇f behaves as a random vector. In section 5.1, we experimentally verify the decoupling approximation for the entire training dynamics.
Our formula (6) possesses two important properties. First, the noise is aligned with the Hessian, which is well known and has been pointed out in the literature [8,11,13,14]. If the loss landscape has flat directions, which correspond to the directions of the Hessian eigenvectors belonging to vanishingly small eigenvalues, the SGD noise does not arise along these directions. Consequently, the SGD dynamics is frozen along the flat directions, which effectively reduces the dimension of the parameter space explored by the SGD dynamics. This reduction plays an important role in the escape efficiency. Indeed, we will see that the escape rate crucially depends on the effective dimension of a given local minimum.
Second, the noise is proportional to the loss function, which is indeed experimentally confirmed in section 5.2. This property has not been pointed out and not been taken into account in previous studies [8,11,13,14] and therefore gives new insights into the SGD dynamics. Indeed, this property allows us to formulate the Langevin equation on the logarithmic loss landscape with simple additive noise as discussed in section 3.2. This new formalism yields the power-law escape rate summarized as theorem 4.7, and the importance of the effective dimension of local minima for their stability.
It should be noted that although Eq. (6) is derived for the mean-square loss, the proportionality between the noise strength and the loss function seems to hold for more general loss functions such as the cross-entropy loss (see appendix D for the detail).
Derivation of Eq. (6): We start from an analytic expression of Σ(θ), which reads where B N is assumed in the second equality. The derivation of Eq. (8) is found in Jastrzȩbski et al. [8], Smith and Le [27]. Usually, the gradient noise variance dominates the square of the gradient noise mean, and hence the term ∇L∇L T in Eq. (8) is negligible.
For the mean-square loss, we have where the decoupling approximation (7) is used in the last equality. Equation (9) is directly related to the Hessian of the loss function near a (local or global) minimum.
The Hessian H(θ) = ∇ 2 L(θ) is written as It is shown by Papyan [31] that the last term of Eq. (10) does not contribute to outliers (i.e. large eigenvalues) of the Hessian. The dynamics near a local minimum is governed by outliers, and hence we can ignore this term. At θ = θ * , we therefore obtain Let us assume C f (θ * ) for θ within the valley of a local minimum θ * . We then obtain the desired expression (6) by substituting it into Eq. (9).

Logarithmized loss landscape
Let us consider the Itô SDE (4) with the SGD noise covariance (6) near a local minimum θ * , which is written as Let us consider a stochastic time t(τ ) for τ ≥ 0 as and perform a random time change from t to τ [32]. Correspondingly, we introduce the Wiener process dW τ ∼ N (0, In terms of the notationθ τ = θ t , Eq. (12) is expressed as We should note that at a global minimum with L(θ) = 0, which is realized in an overparameterized regime [33], the random time change through Eq. (13) is ill-defined since τ is frozen at a finite value once the model reaches a global minimum. We can overcome this difficulty by adding an infinitesimal constant > 0 to the loss, which makes the loss function positive without changing the finite-time dynamics like the escape from a local minimum θ * with L(θ * ) > 0.
In this way, the Langevin equation on the loss landscape L(θ) with multiplicative noise is transformed to that on the logarithmic loss landscape U (θ) = log L(θ) with simpler additive noise. This formulation indicates the importance of considering the logarithmic loss landscape U (θ) = log L(θ). In the following, we use Eq. (14) to discuss the escape efficiency from local minima.

Escape rate from local minima
Now we present our main result on the escape from a local minimum θ * . First, we formally define some key concepts regarding the escape problem.
Definition 4.1 (basin of attraction). For a given local minimum θ * , its basin of attraction A θ * (or a valley of θ * ) is defined as the set of all the starting points θ 0 such that θ t → θ * as t → ∞ when there is no noise. Definition 4.2 (escape rate). We denote by P * the total probability within A θ * and by J the total flux of the probability current across the boundary of A θ * . 1 The escape rate κ is then defined as J /P * [18].
We now calculate the asymptotic behavior of the escape rate in the weak-noise limit η/B → +0, which has been investigated in many fields, including physics [18,34], chemistry [19], stochastic processes [35,36], and machine learning [8,13]. In these previous studies, the following two basic assumptions are commonly used, which are expected to hold in the weak-noise limit. Indeed, the Kramers formula [18], which is obtained by using Assumptions 4.3 and 4.4 below, is rigorously justified [35]. Assumption 4.3 (quasi-stationary approximation). The parameters θ obey the quasistationary distribution near θ * , which is identified as the stationary distribution restricted to A θ * [37].
Assumption 4.4 (most probable escape path). The escape from θ * is dominated by the most probable escape path (MPEP) [38,39]. It is known that the direction of the MPEP at one point must be the direction of one eigenvector of the Hessian at that point, and the MPEP passes a saddle point θ s at the boundary of A θ * . Now we introduce further key assumptions in our analysis for P > 1.
Assumption 4.5 (Hessian outliers). The eigenvalue spectrum of the Hessian H(θ * ) at θ * has n nonzero eigenvalues, which are called "outliers". The remaining P − n eigenvalues are vanishingly small. Assumption 4.6 (low-dimensional subspace of SGD). Outlier eigenvectors of the Hessian v 1 , v 2 , . . . , v n ∈ R P do not change within A θ * , and the SGD dynamics is restricted to the n-dimensional subspace spanned by those outlier eigenvectors (n is called the effective dimension of the local minimum). Assumption 4.5 is empirically confirmed [20,21]. Since the SGD dynamics is frozen along flat directions as we pointed out in section 3.1, the restriction to the n-dimensional outlier subspace in Assumption 4.6 is justified. Indeed, it is empirically observed that the SGD dynamics is actually restricted to a low-dimensional subspace spanned by top eigenvectors of the Hessian [40]. Now we are ready to state our main result in the form of the following theorem.
Theorem 4.7. Let the model parameter θ t evolve according to the SDE in Eq. (12). Under Assumptions 4.3 to 4.6, the escape rate κ asymptotically behaves as as η/B → +0 2 , where θ s is the saddle on the MPEP, h * e is the eigenvalue of the Hessian at θ * along the MPEP, and h s e is the negative eigenvalue of the Hessian at θ s along the MPEP.
Proof. The proof is given in Appendix B. increases with h * e and n, which indicates that sharp minima (i.e. minima with large h * e ) or minima with large n are unstable. This fact explains why SGD finds flat minima with a low effective dimension n. Equation (15) also implies that the effective dimension of any stable minima must satisfy The instability due to a large effective dimension is a new insight naturally explained by the picture of the logarithmic loss landscape. Let us summarize novel aspects of the escape rate derived here, which are not found in previous studies [8,9,[11][12][13][14]. First, the escape rate in Eq. (15) takes a power-law form with respect to the ratio between L(θ * ) and L(θ s ), which differs from a conventional exponential form with respect to the difference ∆L = L(θ s )−L(θ * ). Intuitively, escaping behavior must be a result of noise in the gradient, and if there is no noise, gradient descent cannot escape any basin of attraction. Our result agrees with this intuition: at the global minimum where L(θ * ) = 0, there is no noise, and the SGD must be stuck in where it is with no escaping behavior, which is directly reflected by our formula. However, the standard exponential escape rate formula predicts a non-zero escape rate even when L(θ * ) = 0, which cannot be correct in principle. Second, Eqs. (15) and (16) imply that the SGD is biased to small effective dimensions n. As we mentioned, it is empirically confirmed that the SGD dynamics is restricted to a low-dimensional subspace spanned by top eigenvectors of the Hessian [40]. Since the restriction to a low-dimensional subspace greatly reduces the capacity of the network [41], Eq. (16) supports the presence of an implicit regularization via the SGD.
The formula in Eq. (15) is tested in section 5.3. Since it is difficult in practice to identify when the barrier crossing occurs, we instead consider the mean first passage time, which imitates the escape time for a non-convex loss landscape. Let us fix a threshold value of the loss function. The first passage time t p is defined as the shortest time at which the loss exceeds the threshold value. We identify the threshold value as L(θ s ), i.e., the loss at the saddle in the escape problem. It is expected that t p is similar to the escape time and proportional to κ −1 . Indeed, at least when noise is isotropic and uniform as in the GLD (5), the inverse of the mean first passage time 1/ t p is logarithmically equivalent to κ, i.e., log κ ∼ − log t p , in the weak-noise asymptotic limit [36,38].

Experiments
Our key theoretical observation is that the SGD noise strength is proportional to the loss function, which is obtained as a result of the decoupling approximation. This property leads us to the Langevin equation (14) with the logarithmic loss gradient and an additive noise through a random time change (13). Equation (14) implies the escape rate (15). In section 5.1, we show that the decoupling approximation is valid during the entire training dynamics. In section 5.2, we measure the SGD noise strength and confirm that it is indeed proportional to the loss function near a minimum. In section 5.3, we experimentally test the validity of Eq. (15) for the escape rate. We will see that numerical results for a linear regression and for a non-linear neural network agree with our theoretical results.

Experimental verification of the decoupling approximation
Let us compare the eigenvalue distribution of the exact matrix ( We consider a binary classification problem using the first 10 4 samples of the MNIST dataset such that we classify each image into even (its label is y = +1) or odd number (its label is y = −1). The network has two hidden layers, each of which has 100 units and the ReLU activation, followed by the output layer of a single unit with no activation. Starting from the Glorot initialization, the training is performed via SGD with the mean-square loss, where we fix η = 0.01 and B = 100. Figure 1 shows histograms of their eigenvalues at different stages of the training: (a) at initialization, (b) after 50 epochs, and (c) after 500 epochs. We see that the exact matrix and the approximate one have statistically similar eigenvalue distributions except for exponentially small eigenvalues during the training dynamics. This shows that the decoupling approximation holds during the entire training dynamics in this experiment.

Measurements of the SGD noise strength
Now we test whether the SGD noise strength is actually proportional to the loss function, which is predicted by Eq. (6), an essential theoretical result of ours. As a measure of the SGD noise strength, let us consider the norm of the noise vector ξ given by Here we present experimental results for two architectures and datasets. First, we consider training of the Fashion-MNIST dataset by using a fully connected neural network with three hidden layers, each of which has 2 × 10 3 units and the ReLU activation, followed by the output layer of 10 units with no activation (classification labels are given in the one-hot representation). Second, we consider training of the CIFAR-10 dataset by using a convolutional neural network. Following Keskar et al. Starting from the Glorot initialization, the network is trained by SGD of the minibatch size B = 100 and η = 0.1 for the mean-square loss. During the training, we measure the training loss and the noise strength N for every epoch. Numerical results are given in Fig. 2. We see that roughly N ∝ L at a later stage of the training, which agrees with our theoretical prediction.
Although N is not proportional to L at an early stage of training, it does not mean that Eq. (9) is invalid there. Since the decoupling approximation is valid for the entire training dynamics, Eq. (9) always holds. The reason why the SGD noise strength does not decrease with the loss function in the early-stage dynamics is that N ≈ 2L(θ) × (1/N ) N µ=1 ∇f (θ, x (µ) ) T ∇f (θ, x (µ) ), but the quantity (1/N ) N µ=1 ∇f (θ, x (µ) ) T ∇f (θ, x (µ) ) also changes during training.
Although Eq. (9) is derived for the mean-square loss, the relation N ∝ L(θ) holds in more general loss functions; see Appendix D for general argument and experiments on the cross entropy loss.

Experimental test of the escape rate formula
We now experimentally verify our escape-rate formula (15). Below, we first present numerical results for a simple linear regression problem, and then for a nonlinear model, i.e., a neural network.

Linear regression
Let us start with the following linear regression problem: each entry of x (µ) ∈ R d and its label y (µ) ∈ R are i.i.d. Gaussian random variables of zero mean and unit variance. We focus on the case of d N . The output for an input x is given by f (θ, x) = θ T x, where θ ∈ R d is the trainable network parameter. We optimize θ via SGD. The meansquare loss L(θ) = (1/2N ) N µ=1 θx (µ) − y (µ) 2 is quadratic and has a unique minimum at θ ≈ 0. the Hessian H = (1/N ) N µ=1 x (µ) x (µ)T has d nonzero eigenvalues, all of which are close to unity. We can therefore identify h * e = h * = 1 and n = d. Before investigating the escape rate, we remark that the stationary distribution of θ in this case is theoretically obtained because of the isotropy of SGD noise. Indeed, H(θ * ) in Eq. (14) can be replaced by h * I d , where I d is the d-dimensional identity matrix, and then the stationary distribution ofθ τ is simply given by the Gibbs distribution under the potential U (θ) = log L(θ):P s (θ) ∝ e −U (θ)/T = L(θ) −1/T , where the "temperature" T is given by T = ηh * /B. In Appendix C, it is shown that the stationary distribution P s (θ) of θ t under Eq. (12) is related toP s (θ) via P s (θ) ∝ L(θ) −1P s (θ). We therefore have Remarkably, it depends on L(θ) polynomially rather than exponentially as in the standard GLD [8,28,29]. This property is closely related to the power-law form of the escape rate. First, we test Eq. (17), i.e. the stationary distribution, for d = 1 and N = 10 5 . We sampled the value of θ k at every 100 iterations (k = j × 100, j = 1, 2, . . . , 10 4 ) and made a histogram. We then fit the histogram to the form P s (θ) ∝ L(θ) −φ and determine the exponent φ. Our theory predicts φ = 1 + B/(ηh * ). Numerical results for the exponent φ are presented in Fig. 3 (a) against B for three fixed learning rates η. In the same figure, theoretical values of φ are plotted in dashed lines. The agreement between theory and experiment is fairly well. For a large learning rate η = 1, the exponent slightly deviates from its theoretical value. This is due to the effect of a finite learning rate (recall that η is assumed to be small in deriving the continuous-time stochastic differential equation).
Next, we test our formula on the escape rate, Eq. (15). Although the mean-square loss is quadratic and no barrier crossing occurs, we can measure the first passage time, which imitates the escape time for a non-convex loss landscape as is mentioned in section 4.

Neural network
The escape-rate formula (15) is also verified in a non-linear model. As in section 5.1, we consider the binary classification problem using the first 10 4 samples of MNINST dataset such that we classify each image into even or odd. The network has one hidden layer with 10 units activated by ReLU, followed by the output layer of a single unit with no activation. We always use the mean-square loss. Starting from the Glorot initialization, the network is pre-trained via SGD with η = 0.01 and B = 100 for 10 5 iterations. We find that after pre-training, the loss becomes almost stationary around at L(θ) ≈ 0.035. We assume that the pre-trained network is near a local minimum. We then further train the pre-trained network via SGD with a new choice of η and B (here we fix B = 10), and measure the first passage time t p . The mean first-passage time over 100 independent runs is measured for varying threshold values which are specified by c = L(θ s )/L(θ * ) > 1. Experimental results are presented in Fig. 3 (c). We see the power-law behavior, which is consistent with our theory.
To further verify our theoretical formula (15), we also compare the power-law exponent for the mean first-passage time with our theoretical prediction B/(ηh * e )+1−n/2. Here, h * e and n are estimated by the Hessian eigenvalues. In Appendix E, we present a numerical result for eigenvalues of the approximate Hessian given by the right-hand side of Eq. (11). By identifying the largest eigenvalue of the Hessian as h * e , we have h * e ≈ 95.6. On the other hand, it is difficult to precisely determine the effective dimension n, but it seems reasonable to estimate n ∼ 10. It turns out that theory and experiment agree with each other by choosing n = 9. Dashed lines in Fig. 3 (c) correspond to our theoretical prediction κ −1 ∝ c B/(ηh * e )+1−n/2 with h * e = 95.6 and n = 9. This excellent agreement shows that our theoretical formula (15) is also valid in non-linear models.

Conclusion
In this work, we have investigated the SGD dynamics via a Langevin approach. With several approximations listed in Appendix A, we have derived Eq. (6), which shows that the SGD noise strength is proportional to the loss function. This SGD noise covariance structure yields the stochastic differential equation (14) with additive noise near a minimum via a random time change (13). The original multiplicative noise is reduced to simpler additive noise, but instead the gradient of the loss function is replaced by that of the logarithmic loss function U (θ) = log L(θ). This new formalism yields the power-law escape rate formula (15) whose exponent depends on η, B, h * e , and n. The escape-rate formula in Eq. (15) explains an empirical fact that SGD favors flat minima with low effective dimensions. The effective dimension of a minimum must satisfy Eq. (16) for its stability. This result as well as the formulation of the SGD dynamics using the logarithmic loss landscape should help understand more deeply the SGD dynamics and its implicit biases in machine learning problems.
Although the present work focuses on the Gaussian noise, the non-Gaussianity can also play an important role. For example, Şimşekli et al. [42] approximated SGD as a Lévydriven SDE, which explains why SGD finds wide minima. It would be an interesting future problem to take the non-Gaussian effect into account.

A. List of approximations and their justifications
In section 3.1, we made several approximations to derive the result (6). For clarity, we list the approximations made and their justifications below.
• In Eq. (8), we make the approximation of B N , which simplifies the expression but is not essential. The main conclusion is not affected by this approximation.
• We ignore ∇L∇L T in Eq. (8), which is justified near a local minimum.
• In Eq. (9), we make the decoupling approximation (7), which is one of the key heuristic approximations in our work. This approximation is verified experimentally in section 5.1.
• We ignore the last term of Eq. (10), which is a common approximation [20]. This approximation is justified if we are only interested in the outliers of the Hessian eigenvalues. Indeed, outliers of the Hessian eigenvalues, which play dominant roles in escape from local minima, are known to be attributable to the first term of the right-hand side of Eq. (10) [31].
• In Eq. (11), we assume that the matrix (1/N ) N µ=1 ∇f (θ, x (µ) )∇f (θ, x (µ) ) T does not change so much in a valley with a given local minimum. This approximation is indirectly verified in Fig. 2 (the proportionality between the loss and the noise strength implies this matrix is actually constant).

B. Proof of theorem 4.7
We now prove theorem 4.7. Under Assumption 4.6, dynamics of θ is restricted to the n-dimensional subspace spanned by the outlier eigenvectors v 1 , v 2 , . . . , v n ∈ R P of the Hessian. Consequently, θ is parametrized by n variables z 1 , z 2 , . . . , z n ∈ R as The Hessian is a diagonal matrix in the basis of (v 1 , v 2 , . . . , v n ) within A θ * because of Assumption 4.6. Therefore, ∂ 2 U/∂z i ∂z j = 0 for any i = j, where U = log L is the logarithmic loss introduced in Eq. (14). It implies that ∂U/∂z i is a function of z i , and hence U is written in the form It should be noted that where the last equality is obtained by putting U = log L and ∂ 2 L/∂z 2 i z=0 = h * i . The stochastic differential equation for the n-dimensional vector z = (z 1 , z 2 , . . . , z n ) T in the rescaled time variable τ is given by dz = −∇ z U dτ + 2ηĤ(θ * )/BdW τ [see Eq. (14), which is equivalent to the following Fokker-Planck equation for the distribution function P (z, τ ) of z at τ : where the probability current density J(z) = (J 1 (z), J 2 (z), . . . , J n (z)) T is explicitly given by Under Assumption 4.4, the escape rate is evaluated by considering the MPEP. It is known that the MPEP aligns with one of the eigenvectors of the Hessian [13]. Without loss of generality, let us assume that the direction of eth eigenvector v e corresponds to the direction of the MPEP. We denote by z ⊥ ∈ R n−1 the displacement perpendicular to the escape direction, and write z = (z e , z ⊥ ). At the saddle z s , h s e < 0 and h s i > 0 for all i = e, where {h s i } is the set of eigenvalues of the Hessian at z s , and We now derive the escape rate formula for asymptotically weak noise η/B → +0 following Kramers [18]. The steady current J ∈ R n is aligned to the escape direction, and hence J e = 0 and J ⊥ = 0 along the MPEP. Let us denote by P * the total probability within the valley A θ * and by J τ the total current per unit change of τ flowing to the outside of the valley through the saddle θ s = θ * +z s . From Assumption 4.3, z is assumed to be in the quasi-stationary distribution. Since the quasi-stationary distribution is concentrated around θ * in the weak-noise limit, L(θ t ) is approximately equal to L(θ * ) before the escape. Consequently, the rescaled time variable τ in Eq. (13) is approximately equal to L(θ * )t, and hence the total probability current J t per unit change of t is given by L(θ * )J τ . The escape rate κ per unit change of t is therefore given by We now evaluate P * and J τ . First, to obtain P * , we must know about the quasistationary distribution P s (z) near θ * . It is obtained by putting ∂P (z, τ )/∂τ = 0 in Eq. (21), which implies J(z) = 0. We therefore have For convenience, we put P s (z) = 0 for z / ∈ A θ * . Then, P * , which is the total probability within A θ * , is given by In the weak-noise limit η/B → +0, we can evaluate the above integral by using the saddle-point method, which is also called the Laplace method. By using Eq. (20), the saddle-point method yields as η/B → +0, and therefore Next, let us evaluate J τ . Along the MPEP (z ⊥ = 0), the current vector J ⊥ perpendicular to the escape direction is zero, and hence the probability distribution near the saddle is evaluated in a similar way as Eq. (26): By substituting it into we obtain where the current along the MPEP (z ⊥ = 0) is denoted by J path := J(z s e , z ⊥ = 0). The total current through the saddle is then evaluated by using the saddle-point method: (n−1)/2 (33) as η/B → +0. When the distribution function is almost stationary, Eq. (21) yields ∂J e (z e , z ⊥ = 0)/∂z e = 0, and hence the current along the MPEP is constant J e (z e , z ⊥ = 0) = J path . By putting z ⊥ = 0 in Eq. (32), we have ηh * e U P (z e , z ⊥ = 0) .
By multiplying e B ηh * e U in both sides and integrating over z e from 0 to z c e , where z c e defined as U e (θ * + z c e v e ) = U e (θ * ) (see Fig. 4), we obtain as η/B → +0, where we used the fact that the probability at z c e is negligible in the weak-noise limit. By using the saddle-point method, the integral in the left-hand side of Eq. (35) By substituting this result in Eq. (35), we obtain where ∆U = U (θ s ) − U (θ * ). The total current J τ in Eq. (33) is then expressed as ηh * e ∆U P (θ * ).
By using Eqs. (29) and (38), the escape rate κ in Eq. (24) is evaluated as Since ∆U = log[L(θ s )/L(θ * )], we finally obtain which is exactly identical to the escape rate formula given in theorem 4.7.

C. Stationary distribution
Suppose thatθ τ obeys a simple Langevin equation The stationary distribution ofθ τ is then given by the Gibbs distributionP s (θ) ∝ e −U (θ)/T . On the other hand, what we want is the stationary distribution P s (θ) of θ t , where θ t =θ τ with τ = t 0 dt L(θ t ). In this section, we show the relation between the two distributions: P s (θ) ∝ L(θ) −1P s (θ). We express the stationary distributions in terms of the long-time average of the delta function: By using the relation τ = t 0 dt L(θ t ), we have dτ = L(θ t )dt. For a sufficiently large t, we also obtain τ ∼ tL, whereL := lim s→∞ (1/s) s 0 dt L(θ t ) denotes the long-time average of L(θ t ). By using them, P s (θ) is rewritten as We thus obtain the desired relation, P s (θ) ∝ L(θ) −1P s (θ).

D. Other loss functions
In our paper, we mainly focus on the mean-square loss, for which we can analytically derive the relation between the loss L(θ) and the SGD noise covariance Σ(θ). An important observation is that the SGD noise strength N is proportional to the loss, i.e., N ∝ L(θ) (see section 5.2 for the definition of N ).  Here, we argue that the relation N ∝ L(θ) also holds in more general situations. During the training, the value of µ will fluctuate from sample to sample. At a certain time step of SGD, let us suppose that N − M samples in the training dataset are already fit correctly and hence µ ≈ 0, whereas the other M samples are not and hence µ = O(1). The loss function is then given by L(θ) = (1/N ) N µ=1 µ ∝ M/N . When µ is small, ∇ µ will also be small. Therefore, for N − M samples with µ ≈ 0, ∇ µ ≈ 0 also holds. The other M samples will have non-small gradients: ∇ µ 2 = O(1). We thus estimate N as N ≈ (1/N ) N µ=1 ∇ T µ ∇ µ ∝ M/N ∝ L(θ). In this way, N ∝ L(θ) will hold, irrespective of the loss function.
However, we emphasize that this is a crude argument. In particular, the above argument will not hold near a global minimum because all the samples are correctly fit there, which implies that µ is small for all µ, in contrast to the above argument relying on the existence of M samples with µ and ∇ µ 2 of O(1). We now experimentally test the relation N ∝ L(θ) for the cross-entropy loss. We

E. Hessian eigenvalues for a neural network in section 5.3
We present numerical results on Hessian eigenvalues in a pre-trained neural network studied in section 5.3. Instead of the exact Hessian, we consider an approximate Hessian given on the right-hand side of Eq. (11), i.e., ∇f (θ, x (µ) )∇f (θ, x (µ) ) T .
Eigenvalues {h i } are arranged in descending order as h 1 ≥ h 2 ≥ · · · ≥ h P (in our model P = 7861). A histogram of the Hessian eigenvalues is presented in Fig. 6 (a). We see that most eigenvalues are close to zero, which corresponds to the bulk, but there are some large eigenvalues, which correspond to the outliers. The largest eigenvalue is λ 1 = 95.6, which is identified as h * in our theoretical formula (15).
Another important quantity is the effective dimension n corresnding to the number of outliers. Since the outliers and the bulk are not sharply separated, it is difficult to precisely determine n. In Fig. 6 (b), we plot h i up to i = 30. From this figure, it seems reasonable to estimate n ≈ 10.
As a heuristic method of determining n, we can consider the following identification: first we define the weight p i for ith eigenmode as p i = h i / P j=1 h j . We then determine n as which gives n = 12.7 in our case. In section 5.3, our formula (15) with n = 9 explains numerical results on the firstpassage time.


extfieldsofstudy:
Computer Science
Medicine

text:
Adapting Off-the-Shelf Source Segmenter for Target Medical Image Segmentation

Unsupervised domain adaptation (UDA) aims to transfer knowledge learned from a labeled source domain to an unlabeled and unseen target domain, which is usually trained on data from both domains. Access to the source domain data at the adaptation stage, however, is often limited, due to data storage or privacy issues. To alleviate this, in this work, we target source free UDA for segmentation, and propose to adapt an ``off-the-shelf"segmentation model pre-trained in the source domain to the target domain, with an adaptive batch-wise normalization statistics adaptation framework. Specifically, the domain-specific low-order batch statistics, i.e., mean and variance, are gradually adapted with an exponential momentum decay scheme, while the consistency of domain shareable high-order batch statistics, i.e., scaling and shifting parameters, is explicitly enforced by our optimization objective. The transferability of each channel is adaptively measured first from which to balance the contribution of each channel. Moreover, the proposed source free UDA framework is orthogonal to unsupervised learning methods, e.g., self-entropy minimization, which can thus be simply added on top of our framework. Extensive experiments on the BraTS 2018 database show that our source free UDA framework outperformed existing source-relaxed UDA methods for the cross-subtype UDA segmentation task and yielded comparable results for the cross-modality UDA segmentation task, compared with a supervised UDA methods with the source data.

Introduction
Accurate tumor segmentation is a critical step for early tumor detection and intervention, and has been significantly improved with advanced deep neural networks (DNN) [25,18,10,9,17]. A segmentation model trained in a source domain, however, usually cannot generalize well in a target domain, e.g., data acquired from a new scanner or different clinical center, in implementation. Besides, annotating data in the new target domain is costly and even infeasible [11]. To address this, unsupervised domain adaptation (UDA) was proposed to transfer knowledge from a labeled source domain to unlabeled target domains [13].
The typical UDA solutions can be classified into three categories: statistic moment matching, feature/pixel-level adversarial learning [15,14,12], and selftraining [33,16]. These UDA methods assume that the source domain data are available and usually trained together with target data. The source data, however, are often inaccessible, due to data storage or privacy issues, for cross-clinical center implementation [1]. Therefore, it is of great importance to apply an "offthe-shelf" source domain model, without access to the source data. For sourcefree classification UDA, Liang et al. [8] proposed to enforce the diverse predictions, while the diversity of neighboring pixels is not suited for the segmentation purpose. In addition, the class prototype [13] and variational inference methods [11] are not scalable for pixel-wise classification based segmentation. More importantly, without distribution alignment, these methods relied on unreliable noisy pseudo labeling.
Recently, the source relaxed UDA [1] was presented to pre-train an additional class ratio predictor in the source domain, by assuming that the class ratio, i.e., pixel proportion in segmentation, is invariant between source and target domains. At the adaptation stage, the class ratio was used as the only transferable knowledge. However, that work [1] has two limitations. First, the class ratio can be different between the two domains, due to label shift [11,13]. For example, a disease incident rate could vary between different countries, and tumor size could vary between different subtypes and populations. Second, the pre-trained class ratio predictor used in [1] is not typical for medical image segmentation, thereby requiring an additional training step using the data in the source domain.
In this work, to address the aforementioned limitations, we propose a practical UDA framework aimed at the source-free UDA for segmentation, without an additional network trained in the source domain or the unrealistic assumption of class ratio consistency between source and target domains. More specifically, our framework hinges on the batch-wise normalization statistics, which are easy to access and compute. Batch Normalization (BN) [6] has been a default setting in the most of modern DNNs, e.g., ResNet [5] and U-Net [30], for faster and more stable training. Notably, the BN statistics of the source domain are stored in the model itself. The low-order batch statistics, e.g., mean and variance, are domain-specific, due to the discrepancy of input data. To gradually adapt the low-order batch statistics from the source domain to the target domain, we develop a momentum-based progression scheme, where the momentum follows an exponential decay w.r.t. the adaptation iteration. For the domain shareable highorder batch statistics, e.g., scaling and shifting parameters, a high-order batch statistics consistent loss is applied to explicitly enforce the discrepancy minimization. The transferability of each channel is adaptively measured first, from which to balance the contribution of each channel. Moreover, the proposed unsupervised self-entropy minimization can be simply added on top of our framework to boost the performance further.
Our contributions are summarized as follows: • To our knowledge, this is the first source relaxed or source free UDA framework for segmentation. We do not need an additional source domain network, or the Fig. 1: Comparison of (a) conventional UDA [28] and (b) our source-relaxed OS-UDA segmentation framework based on the pre-trained "off-the-shelf" model with BN. We minimize the domain discrepancy based on the adaptively computed batch-wise statistics in each channel. The model consists of a feature encoder (Enc) and a segmentor (Seg) akin to [3,32].
unrealistic assumption of the class ratio consistency [1]. Our method only relies on an "off-the-shelf" pre-trained segmentation model with BN in the source domain.
• The domain-specific and shareable batch-wise statistics are explored via the low-order statistics progression with an exponential momentum decay scheme and transferability adaptive high-order statistics consistency loss, respectively.
• Comprehensive evaluations on both cross-subtype (i.e., HGG to LGG) and cross-modality (i.e., T2 to T1/T1ce/FLAIR) UDA tasks using the BraTS 2018 database demonstrate the validity of our proposed framework and its superiority to conventional source-relaxed/source-based UDA methods.

Methodology
We assume that a segmentation model with BN is pre-trained with source domain data, and the batch statistics are inherently stored in the model itself. At the adaptation stage, we fine-tune the model based on the batch-wise statistics and the self-entropy (SE) of target data prediction. The overview of the different setups of conventional UDA and our "off-the-shelf (OS)" UDA is shown in Fig. 1. Below, we briefly revisit the BN in Subsec. 2.1 first and then introduce our OSUDA in Subsec. 2.2. The added unsupervised SE minimization and the overall training protocol are detailed in Subsec. 2.3.

Preliminaries on Batch Normalization
As a default setting in the most of modern DNNs, e.g., ResNet [5] and U-Net [30], Batch Normalization (BN) [6] normalizes the input feature in the lth layer f l ∈ R B×H l ×W l ×C l within a batch in a channel-wise manner to have zero mean and unit variance. B denotes the number of images in a batch, and H l , W l , and C l are the height, width, and channels of layer l. We have samples in a batch, with index b ∈ {1, · · · , B}, spatial index n ∈ {1, · · · , H l × W l }, and channel index c ∈ {1, · · · , C l }. BN calculates the mean of each channel Then, the input feature is normalized as where ∈ R + is a small scalar for numerical stability. γ l,c and β l,c are learnable scaling and shifting parameters, respectively. In testing, the input is usually a single sample rather than a batch with B samples. Therefore, BN stores the exponentially weighted average of the batch statistics at the training stage and used it in testing. Specifically, the mean and variance over the training are tracked progressively, i.e., where η ∈ [0, 1] is a momentum parameter. After K training iterations, µ K l,c , {σ 2 } K l,c , γ K l,c , and β K l,c are stored and used for testing normalization [6].

Adaptive source-relaxed batch-wise statistics adaptation
Early attempts of BN for UDA simply added BN in the target domain, without the interaction with the source domain [7]. Recent studies [2,20,26,19] indicated that the low-order batch statistics, i.e., mean µ l,c and variance {σ 2 } l,c , are domain-specific, because of the divergence of cross-domain representation distributions. Therefore, brute-forcing the same mean and variance across domains can lead to a loss of expressiveness [29]. In contrast, after the low-order batch statistics discrepancy is partially reduced, with domain-specific mean and variance normalization, the high-order batch statistics, i.e., scaling and shifting parameters γ l,c and β l,c , are shareable across domains [20,26]. However, all of the aforementioned methods [2,20,29,26,19] require the source data at the adaptation stage. To address this, in this work, we propose to mitigate the domain shift via the adaptive low-order batch statistics progression with momentum, and explicitly enforce the consistency of the high-order statistics in a source-relaxed manner. Low-order statistics progression with an exponential momentum decay scheme. In order to gradually learn the target domain-specific mean and variance, we propose an exponential low-order batch statistics decay scheme. We initialize the mean and variance in the target domain with the tracked µ K l,c and {σ 2 } K l,c in the source domain, which is similar to applying a model with BN in testing [6]. Then, we progressively update the mean and variance in the t-th adaptation iteration in the target domain as where η t = η 0 exp(−t) is a target adaptation momentum parameter with an exponential decay w.r.t. the iteration t. µ t l,c and {σ 2 } t l,c are the mean and variance of the current target batch. Therefore, the weight of µ K l,c and {σ 2 } K l,c are smoothly decreased along with the target domain adaptation, while µ t l,c and {σ 2 } t l,c gradually represent the batch-wise low-order statistics of the target data. Transferability adaptive high-order statistics consistency. For the highorder batch statistics, i.e., the learned scaling and shifting parameters, we explicitly encourage its consistency between the two domains with the following high-order batch statistics (HBS) loss: where γ K l,c and β K l,c are the learned scaling and shifting parameters in the last iteration of pre-training in the source domain. γ t l,c and β t l,c are the learned scaling and shifting parameters in the t-th adaptation iteration. α l,c is an adaptive parameter to balance between the channels.
We note that the domain divergence can be different among different layers and channels, and the channels with smaller divergence can be more transferable [22]. Accordingly, we would expect that the channels with higher transferability contribute more to the adaptation. In order to quantify the domain discrepancy in each channel, a possible solution is to measure the difference between batch statistics. In the source-relaxed UDA setting, we define the channel-wise sourcetarget distance in the t-th adaptation iteration as Then, the transferability of each channel can be measured by α l,c = L×C×(1+d l,c ) −1 l c (1+d l,c ) −1 . Therefore, the more transferable channels will be assigned with higher importance, i.e., with larger weight (1 + α l,c ) in L l,c .

Self-entropy minimization and overall training protocol
The training in the unlabeled target domain can also be guided by an unsupervised learning framework. The SE minimization is a widely used objective in modern DNNs to encourage the confident prediction, i.e., the maximum softmax value can be high [4,8,24,1]. SE for pixel segmentation is calculated by the averaged entropy of the classifier's softmax prediction given by where H 0 and W 0 are the height and width of the input, and δ b,n is the histogram distribution of the softmax output of the n-th pixel of the b-th image in a batch. Minimizing L SE leads to the output close to a one-hot distribution. At the source-domain pre-training stage, we follow the standard segmentation network training protocol. At the target domain adaptation stage, the overall training objective can be formulated as L = L HBS + λL SE , where λ is used to balance between the BN statistics matching and SE minimization. We note that a trivial solution of SE minimization is that all unlabeled target data could have the same one-hot encoding [4]. Thus, to stabilize the training, we linearly change the hyper-parameter λ from 10 to 0 in training.
For evaluation, we adopted the widely used Dice similarity coefficient and Hausdorff distance metrics as in [32]. The Dice similarity coefficient (the higher, the better) measures the overlapping part between our prediction results and the ground truth. The Hausdorff distance (the lower, the better) is defined between two sets of points in the metric space.

Cross-subtype HGG to LGG UDA
HGG and LGG have different size and position distributions for tumor regions [23]. Following the standard protocol, we used the HGG training set (source Fig. 2: The comparison with the other UDA methods, and an ablation study of adaptive channel-wise weighting and SE minimization for HGG to LGG UDA. domain) to pre-train the segmentation model and adapted it with the LGG training set (target domain) [23]. The evaluation was implemented in the LGG testing set. We adopted the same 2D U-Net backbone in [23], sliced 3D volumes into 2D axial slices with the size of 128×128, and concatenated all four MRI modalities to get a 4-channel input.
The quantitative evaluation results are shown in Table 1. Since the pixel proportion of each class is different between HGG and LGG domains, the class ratio-based CRUDA [1] only achieved marginal improvements with its unsupervised learning objective. We note that the Dice score of the core tumor was worse than the pre-trained source-only model, which can be the case of negative transfer [27]. Our proposed OSUDA achieved the state-of-the-art performance for source-relaxed UDA segmentation, approaching the performance of SEAT [23] with the source data, which can be seen as an "upper-bound. " We used OSUDA-AC and OSUDA-SE to indicate the OSUDA without the adaptive channel-wise weighting and self-entropy minimization, respectively. The better performance of OSUDA over OSUDA-AC and OSUDA-SE demonstrates the effectiveness of adaptive channel-wise weighting and self-entropy minimization. The illustration of the segmentation results is given in Fig. 2. We can see that the predictions of our proposed OSUDA are better than the no adaptation model. In addition, CRUDA [1] had a tendency to predict a larger area for the tumor; and the tumor core is often predicted for the slices without the core.
3.2 Cross-modality T2 to T1/T1ce/FLAIR UDA Because of large appearance discrepancies between different MRI modalities, we further applied our framework to the cross-modality UDA task. Since clinical annotation of the whole tumor is typically performed on T2-weighted MRI, the typical cross-modality UDA setting is to use T2-weighted MRI as the labeled source domain, and T1/T1ce/FLAIR MRI as the unlabeled target domains [32]. We followed the UDA training (80% subjects) and testing (20% subjects) split as in [32], and adopted the same single-channel input backbone. We note that the data were used in an unpaired manner [32]. The quantitative evaluation results are provided in Table 2. Our proposed OSUDA outperformed CRUDA [1] consistently. In addition, in CRUDA, the additional class ratio prediction model was required to be trained with the source data, which is prohibitive in many real-world cases. Furthermore, our OSUDA outperformed several UDA methods trained with the source data, e.g., Cycle-GAN [31] and SIFA [3], for the two metrics. The visual segmentation results of three target modalities are shown in Fig. 3, showing the superior performance of our framework, compared with the comparison methods.

Discussion and Conclusion
This work presented a practical UDA framework for the tumor segmentation task in the absence of the source domain data, only relying on the "off-the-shelf" pre-trained segmentation model with BN in the source domain. We proposed a low-order statistics progression with an exponential momentum decay scheme to gradually learn the target domain-specific mean and variance. The domain shareable high-order statistics consistency is enforced with our HBS loss, which is adaptively weighted based on the channel-wise transferability. The performance was further boosted with the unsupervised learning objective via self-entropy minimization. Our experimental results on the cross-subtype and cross-modality UDA tasks demonstrated that the proposed framework outperformed the comparison methods, and was robust to the class ratio shift.