A Smart Monitoring of a Water Quality Detector System

Received Jan 17, 2018 Revised Mar 28, 2018 Accepted Apr 12, 2018 The importance to monitor the water quality level is undeniable due to significant impact to human health and ecosystem. The project aims to develop a wireless water quality monitoring system that aids in continuous measurements of water conditions based on pH and turbidity measurements. These two sensors are connected to microprocessor and transmitted to the database by using a Wi-Fi module as a bridge. The developed system was successfully detect both the pH and turbidity values hence updating in IoT platform. Based on the results obtained, the test water sample can be classified to class IIB which is suitable for water recreational used body contact. Overall, the developed system offers fast and easy monitoring of pH and turbidity levels with IoT application for continuous maintenance of clean water. The work is just concern on the physical water parameters hence further extend to chemical parameter for verifying a better result in measuring the WQI value.

INTRODUCTION
The water quality index (WQI) is always used in order to justify or classify the level of water quality. WQI acts as mean indication of water quality measurement through the determination of physico-chemical parameters of surface water [1]. It is important to measure the level of the water quality before consuming or safe use for other purposes. For water quality classification, several physical, biological and chemical parameters that have significant impact on the water quality have to be identified. The measurements obtained are then compared to water quality standards of Malaysia.
However, water samples are normally collected at regular periods continued by laboratory analysis thus ask for larger time consumption. Besides, there are certain chemical and biological processes such as oxidation-reduction potential that require for on-site measurement to ensure accuracy and other aspects to be considered as discussed in [2]. Therefore, online water quality monitoring systems which potential to deliver continuous data is highly demanded. The objective of this work is to develop a low-cos wireless water quality monitoring system that aids in continuous measurements of water conditions. The Department of Environment (DOE) water quality parameters and index including its classes is presented in Table 1. According to Table 1, there are six parameters to be measured for water quality control which are Ammoniacal Nitrogen, Biochemical Oxygen Demand (BOD), Chemical Oxygen Demand (COD), Dissolved Oxygen (DO), pH, and Total Suspended Solid (TSS) [3]. However, since this is an initial study of water project, two parameters of DOE stated which are pH and total suspended solid (tss) are just emphasized. There are various studies of water quality monitoring system such discussed in [4] and [5]. In [4], a system which enables remote probing and real-time monitoring of the water quality parameters and enables observation of current and historical water quality status are developed. A ZigBee based on wiresless sensor network monitoring is applied. [5] implements water monitoring system which cover the analyses of pH, conductivity, dissolved oxygen and temperature's level. Sounds an alarm will be triggered if there is a water contamination or change in water quality. The parameters are measured with off-the shelf sensors and data is sent to a base station via GPRS. Meanwhile, the aim of [2] is to develop a low-cost wireless water quality monitoring system that aids in continuous measurements of water conditions. A prototype sensor as one component of the Autonomous Live Animal Response Monitor (ALARM) has been developed. The system measure temperature, light intensity, pH, electrical conductivity, total dissolved solids, salinity, dissolved oxygen and oxidation reduction potential as the parameters. On top of that, several previous works has been reviewed in terms of hardware and software selection such as microcontroller, pH and turbidity sensors.
Based on [6], the microcontroller used for the system is a Raspberry Pi (RPi) model B. This model has 256 Mb of RAM memory and 700 MHz processor and its architecture allows the development of several, different embedded system projects through the GPIO pins. This makes less demanding the collaboration with other equipment. The RPi is a very low-cost computer, with the compact size which supports Unix platforms. The difference for the proposed project is the used of WeMos D1 R2. The WeMos D1 R2 have almost same layout with the Arduino UNO except the number of analog I/O pins provided. The Arduino UNO offers six analog I/O pins while WeMos D1 only provided with one analog I/O pin. On the bright side, the WeMos is already embedded with a Wi-Fi module ESP8266 and cheaper compared to Arduino UNO. With Wi-fi module, the project is specifically use the ESP-12F that was embedded on the WeMos D1 R2.
In terms of pH sensor, [2] uses pH sensor from Phidgets. It measures the full pH range from 0 to 14 and operates in the temperature range of 0•−80• C. The Phidgets also supplies adaptor to convert BNC to analog voltage after sensing. Using this adaptor, the pH sensor data is acquired at Arduino analog pin. Meanwhile, [7] measures pH through the use of a conventional glass electrode with a reference electrode setup, the other is using an Ion-Selective-Field-Effect-Transistor (ISFET). The pH sensor will consist of a conventional glass electrode as these electrodes are more reliable and economical for long term monitoring. However, this project will use an analog pH sensor SKU:SEN0161 by DFRobot. It also measure the pH range from 0 to 14 and operates in the temperature range from 0 o C-60 o C. It has an LED which works as the Power Indicator, a BNC connector and PH2.0 sensor interface.
Meanwhile, turbidity sensor detects suspended particles in water by measuring the light transmittance and scattering rate which changes with the amount of total suspended solids (tss) in water. As the tss increases, the liquid turbidity level increases. [2] uses a light sensor to measure water turbidity which operate at 3.3V and 5V. The range of light intensity can be measured from 0−1000 lx and sensor output is of nonratiometric type. However, the project will use a turbidity sensor SKU: SEN0189 by DFRobot. Other research work studies the water monitoring dan control can be further review in [8][9][10]. The paper is organized as follows. The project implementation is explained in Section 2 while the result and analysis obtained is discussed in Section 3. The conclusion of the project drawn is next presented in Section 4.

RESEARCH METHOD
The implementation of this project can be divided into two major sections; hardware and software implementation.

Hardware Implementation
The hardware involves are microcontroller, pH and turbidity sensors and Wi-Fi module. In general, the microcontroller processes the data taken by pH and turbidity sensors hence sending the data to the database by using the Wi-Fi module as a bridge. Figure 1 shows a WeMos D1 R2 Arduino Compatible layout or commonly referred as WeMos D1. As discussed, WeMos D1 microcontroller has the same layout with Arduino UNO except a Wi-Fi module embedded onto the board.  Figure 2 shows an analog pH sensor SKU:SEN0161 by DFRobot. The sensor is connected to the first WeMos D1 analog pin, A0 as shown in Figure 2. Turbidity sensor SKU: SEN0189 by DFRobot is shown in Figure 3. This sensor provides analog and digital signal output modes while the threshold is adjustable when in digital signal mode. This sensor is next connected to the second WeMos D1 analog pin, A0.

Wi-Fi module
The ESP8266 module type ESP-12F is used as a Wi-Fi Module as shown in Figure 4. This Wi-Fi Module is a self-contained SOC with integrated TCP/IP protocol stack that can give any microcontroller access to Wi-Fi network. The ESP8266 is capable of either hosting an application or offloading all Wi-Fi networking functions from another application processor. Each ESP8266 module comes pre-programmed with an AT command set firmware, which can simply plug-and-play to Arduino device and get about as much Wi-Fi-ability as a Wi-Fi Shield offers. The data collected were next archived in online database, Ubidots. Ubidots are specializing in hardware and software development for IoT projects. It is an open source software that facilitates developers to easily capture sensor data and turn it into useful information. The Ubidots platform can be used to send data to the cloud from any internet-enabled device by configuring actions and alerts based on the real-time data and disengage the value of data through visual tools.
Finally, the data was displayed at webpage created by open software, Wix. The wix.com is one of the most popular HTML5 website builders (cloud-based web development). The Wix introduced a mobile editor that the size or the style better in mobile viewing. The concept of building the websites by using their online drag and drop tools is very easy to use. The block diagram of developed complete system is presented in Figure 5. The pH sensor and turbidity sensor were connected to the Arduino. The Arduino will then process the data and send it to database through Wi-Fi module ESP8266. The user will be able to receive the updated data hence allowing for water quality level classification based on DOE measurement.
The pH and turbidity data were recorded in both offline and online measurement. For offline action, the data were recorded for five days and maintained to be collected on 5pm to 6pm. The data were next continuously collected by online measurement in IoT platform.

RESULTS AND ANALYSIS
A complete system of Smart Monitoring Water Quality Detector is presented in Figure 6. As observed, the system is floated dan tested at Tasik UTeM.

Experiment Procedure of pH and Turbidity Sensors
Upon completion, the system has been tested and calibrated to achieve stability and accuracy of the measurements. The result of pH troubleshoots's are discussed:

pH sensor
The pH experiment is divided into three categories which are alkaline, acidic, or neutral. These experiments will be conducted in same order but with different test pH solution. The pH sensor is dip into the pH solutions while the reading is recorded. To ensure the reading accuracy, adjustment the "Offset" value in the codes has been done to ±0.10. Table 2-4 show the results of alkaline, acidic and neutral test pH solution using the modeled developed.
According to Table 2-4, it was proved that the applied pH sensor is accurately read the pH level; with the differences between the actual reading and the experimental value is not more than 0.9% on average of percentage error. It is noticed that the pH electrode must be initially washed with distilled water for each measurement of pH test solution for a reliable reading.

Turbidity sensor
The system was tested with filtered water that has actual reading of 0.1. Again, adjustment of "Offset" value in the codes has been done to ±0.10 to ensure the reading accuracy.

System Results
The system is next setup at "Tasik UTeM". Two assessments have been done where the data were collected in offline and online measurements.

Offline measurements
The data of offline measuremet was recorded for five days in two weeks; and maintained to be collected at every 5pm to 6pm. The results obtained are shown in Figure 7.

Online measurements
The data of online measurement was collected with the function of the IoT application. The system was connected to the Service Set Identifier (SSID) provided and transmitted the data to the online database, Ubidots. The Figure 8 shows a result collected by the online database, Ubidots.

Figure 8. pH and tss results on Ubidots platform
The pH results indicate that the tested water is in the neutral range with slightly acidic. The results also indicate a reading at pH 7 and above. The differences in reading is observed due to weather changing which in rainy or sunny days. It was studied that that weather has obviously effects to water pH level. Acid rain is comprehensively acknowledged to be accountable for acidifying soil and water in particular area, mostly in urban area [12]. The rain that fell on soil was mix with fertilizer and then flow down to the lake. This mixture are acidic, is then mix together with the water hence, may turns the water slightly acidic. Meanwhile, there were differences reading by 7 mg/l using offline and online approaches for tss measurement. It was studied that the rain gave an impact on water turbidity especially during heavy rain. It was found that rain events of the early rainy season through mid-tempestuous season are conclusive in initiating turbidity increment [13]. The rain transported the sediment into the water and turns the water murky hence increasing the tss.
The results collected can be summarized to average pH of 6.85 and tss of 50 mg/l; which reflect to class II of WQI (refer Table 1). According to WQI, the IIA needs a conventional water treatment if it uses as a water supply, while IIB is suitable for water recreational used body contact.

CONCLUSION
The project aims to develop a smart system to monitor the water quality level based on Water Quality Index (WQI) as a standard practice by Department of Environment (DOE) Malaysia. The developed system was successfully detect both pH and turbidity levels hence updating in Wix IoT platform. Overall, the system developed offers fast and easy monitoring of pH and turbidity level in ensuring clean water is continuously maintained. Based on the results obtained, the water level of "Tasik UTeM" can be classified into class IIB which is suitable for water recreational used body contact. However, the presented work was just concern on the physical water parameters hence further extend to chemical parameter for verifying a better result in measuring the WQI value.


extfieldsofstudy:
Computer Science
Mathematics

text:
Identification and Visualization of the Underlying Independent Causes of the Diagnostic of Diabetic Retinopathy made by a Deep Learning Classifier

Interpretability is a key factor in the design of automatic classifiers for medical diagnosis. Deep learning models have been proven to be a very effective classification algorithm when trained in a supervised way with enough data. The main concern is the difficulty of inferring rationale interpretations from them. Different attempts have been done in last years in order to convert deep learning classifiers from high confidence statistical black box machines into self-explanatory models. In this paper we go forward into the generation of explanations by identifying the independent causes that use a deep learning model for classifying an image into a certain class. We use a combination of Independent Component Analysis with a Score Visualization technique. In this paper we study the medical problem of classifying an eye fundus image into 5 levels of Diabetic Retinopathy. We conclude that only 3 independent components are enough for the differentiation and correct classification between the 5 disease standard classes. We propose a method for visualizing them and detecting lesions from the generated visual maps.

Introduction
Diabetic Retinopathy (DR) is a leading disabling chronic disease and one of the main causes of blindness and visual impairment in developed countries for diabetic patients. Studies reported that 90% of the cases can be prevented through early detection and treatment. Eye screening through retinal images is used by physicians to detect the lesions related with this disease. Due to the increasing number of diabetic patients, the amount of images to be manually analyzed is becoming unaffordable. Moreover, training new personnel for this type of image-based diagnosis is long, because it requires to acquire expertise by daily practice [1], [2].
Deep Learning (DL) [3], [4], is a subfield of Machine Learning that allow the automatic model construction of very effective image classifiers using a parametric model. These models are able to identify and extract the statistical regularities present in data that are important for optimizing a defined loss function, with the final objective of mapping a high-multidimensional input into a smaller multi-dimensional output (f: R n → R m , n m). This mapping allows the classification of multi-dimensional objects into a small number of categories. The model is composed by many neurons that are organized in layers in a hierarchical way. Every neuron receives the input from a predefined set of neurons. Every connection has a parameter that corresponds to the weight of the connection. The function of every neuron is to make a transformation of the received inputs into a calculated output value. For every incoming connection, the weight is multiplied by the input value received by the neuron. The aggregation of all inputs passes through an activation function that calculates the neuron output. Parameters are usually optimized using a stochastic gradient descent algorithm that minimizes a predefined loss function. These hierarchical models are able to learn multiple levels of representation that correspond to different levels of abstraction, which enables the representation of complex concepts in a compressed way [5], [6].
Convolutional neural networks (CNN) are a deep learning subfield that proved to be very effective for image classification, detection and segmentation. The first successful CNN was presented in [7], and it was designed for hand-written digit recognition. This early CNN implementation used a combination of convolution, pooling and non-linearity that has been the key feature of DL since now. The DL breakthrough took place with the publication of [8], where for the first time a CNN won the Imagenet [9] classification competition by a large margin. In that paper a set of innovative techniques where introduced including data augmentation, the use of rectified linear units (ReLUs) as activation function, the use of dropout for avoiding overfitting, overlapping maxpooling for avoiding the averaging effects of avg-pooling and the use of GPUs for speeding up the training time. Later on, different CNN improvements where also published. In [10], it was the first time that small 3x3 convolution filters where used and combined as a sequence of convolutions. In [11] the authors introduced residual networks. This networks used a combination of 3x3 convolutional layers with a by-pass of the input every two layers that was summed up to the output. Such bypasses improved the gradient propagation through the network allowing the design of deeper networks and improving the classification capabilities.
DL models have been also successfully applied in many medical classification tasks. In [12] a DL classifier was designed achieving dermatologist-level accuracy for skin cancer detection. In [13] a 3D CNN for automated pulmonary nodule detection and classification was designed. In [14] a CNN Alzheimer's disease classifier with high performance was also described. In [15] the authors presented a diabetic retinopathy classifier with better performance than human experts in the detection of the most severe cases of the disease.
Ordinal regression is a term used for multi-class classification for the cases where a underlying property can be used for prestablishing an ordering of the classes. Several quality measures exist in the literature of machine learning and statistics for ordinal regression [16]. Kappa is a well-known statistic coefficient defined by Cohen [17] to measure inter-rater agreement in such cases. Weighted Kappa [18] is another index used for measuring the goodness of a ordinal regression classification. In this last coefficient, disagreements are penalized proportionally to a power of the distance between classes. The penalization most commonly used is the quadratic. In such cases, the index is commonly referred as Quadratic Weighted Kappa (referred as QWK or κ indistinctly in this paper).
In medical diagnosis tasks is important not only the accuracy of predictions but also the reasons behind decisions. Self-explainable models enable the physicians to contrast the information reported with their own knowledge, increasing the probability of a good diagnostic, which may have a significant influence in patient's treatment.
Different attempts have been done for interpreting the results reported by neural networks. In [19] a network propagation technique is used for input-space feature visualization. After this work, [20] used a pixel-wise decomposition for classification decision. This decomposition could be done in two ways: considering the network as a global function, disregarding its topology (functional approach) or using the natural properties of the inherent function topology for applying a message passing technique, propagating back into the pixel space the last layer output values. After this work, in [21] they used a so named Deep Taylor decomposition technique to replace the inherently intractable standard Taylor decomposition using a multitude of simpler analytically tractable Taylor decompositions.
Independent Component Analysis (ICA) [22] is a statistical method for the separation of a multi-dimensional random signal into a linear combination of components that are statistically as independent from each other as possible. The theoretical foundation of ICA is based on the Central Limit Theorem, which establishes that the distribution of the sum (average or linear combination) of N independent random variables approaches a gaussian as N → ∞. When ICA method is applied, it is assumed that such separation exist, ie. that is possible to express the signal as a linear combination of independent components (IC). Perfect independence between random variables is achieved when mutual information between them is zero. Mutual information can be expressed as the Kullback-Leibler divergence between the joint distribution and the product of the distributions of each variable. Mutual information can be decomposed, under linear transforms, as the sum of two terms: one term expressing the decorrelation of the components and one expressing their non-Gaussianity [23]. ICA uses optimization to calculate its components. Two types of optimization objectives can be used: minimize the mutual information or maximize the nongaussianess of each component.
In this paper we study a technique that allows the identification, separation and visualization in the input and hidden space of the IC responsible of a particular DR classification decision taken by a DL classifier given a certain eye fundus image. This is done by calculating the minimum number of IC that are able to encode the maximum information about the particular classification. Identifying such components, we reduce the redundancy of feature space and identify the components that share the minimum mutual information between each other. In that way, under the supposition that the feature extraction phase of the deep learning model has been able to disentangle completely the feature information, we are able to separate the independent elements causing the disease.
We use the pixel-wise visualization method to identify in pixel space such independent causes. We use the assessment of experts clinicians for interpreting such IC.
The paper is structured as follows: in Section 2 the current work on DL applied to DR is briefly introduced, then, the main works on interpretation of DL are presented. Section 3 we present the methods, Section 4 presents the results showing samples of the kind of visual interpretations given by the model and finally Section 5 present the final conclusions of our work.

Related Work
Many DL based DR classifiers have been published in the last years [24], [25], [15].
In [15] a binary DL classifier was published for the detection of the most severe cases of DR (grouping classes 0 and 1 of DR on one side, and classes 2, 3 and 4 on another). This model was trained using an extended version of the EyePACS dataset mentioned before with a total of 128,175 images and improving the proper tagging of the images using a set of 3 to 7 experts chosen from a panel of 54 US expert Ophtalmologists. This model surpassed the human expert capabilities, reaching at the final operating point approximately 97% sensitivity and 93.5% specificity in the test sets of about 10,000 images for detecting the worse cases of DR. The strength of this model was its ability to predict the more severe cases with a sensitivity and specificity greater than human experts.
In our previous work [25] a DL classifier was published for the prediction of five different disease grades. This model was trained using the public available EyePACS dataset. The training set had 35,126 images and the test set 53,576. The quadratic weighted kappa (QWK) evaluation metric [18] over the test set using a unique DL model without ensembling was close to the reported by human experts.
The drawback of all these implementations, as of many other DL based models, is its lack of interpretability. In last years different approximations have been proposed for converting the initial DL black box classifiers into interpretable ones. Between the more successful interpretation models existing today we have the sensitivity maps [26], layer-wise relevance propagation [20] and Taylor type decomposition models [21].
All these methods use different strategies to backpropagate the classification scores to the input space and distribute the value of the final classification between input pixels. With this score distribution is possible to identify the most relevant pixels for a particular classification.
Layer-wise relevance propagation models allow not only to report the predicted class but also to evaluate the importance of every input pixel of the image in the final classification decision. In such a way, it is possible to determine which pixels are more important in the final decision and facilitate the human experts the construction of rationale explanations based on the interpretation of such maps.

Methods
DL models are organized in layers, being the inputs of each one a combination of the outputs of previous ones. We design the output layer to be a linear combination of last layer feature space components. In this way we are forcing the model to disentangle the important features that, combined in a linear way, allow the achievement of a maximum possible classification score. These components (or other obtained as a linear combination of them, like with ICA analysis), are easy to analyze due to the linear nature of its relationship with the classification scores.

ICA based interpretation procedure
In this paper we go a step forward in the interpretation of score maps. Instead of generating directly the pixel maps associated with a particular class, we try to identify and separate independent elements associated with the disease. Our new contribution comes from the identification, separation and visualization of the IC that explain a particular classification decision. We hypothesize that in order to achieve high performance classification scores, the network has to encode the information required to make the classification. Human experts base its decisions in the number and types of lesions present in the image. That's why in some way, such information has to be present in a disentangled form in last layer feature space, previous to the output layer. Instead of directly visualizing the more important pixels under a classification decision, we split the information of such last feature layer into independent features using a Independent Component Analysis (ICA). A posteriori we use a pixel-wise relevance propagation derived method to visualize such independent component in input space. In this way, we can, not only generate importance pixel maps, but also differentiate between the underlying independent causes of the disease.
We study the last layer feature space, previous to the output layer linear combination in order to identify its properties and try to isolate the independent elements that are causing a particular classification. For this purpose we use a principal component analysis (PCA) [27] to appraise the redundancy of this space and a ICA [22] using different number of components to identify the minimum of them required to achieve a classification score close enough to the achieved without such a dimensional reduction. ICA allows to find a linear representation of non-Gaussian data so that the components are statistically independent, or as independent as possible. Such a representation seems to capture the essential structure of the data in many applications, including feature extraction [22]. When the data is not Gaussian, there are higher order statistics beyond variance that are not being taken into account by PCA. While PCA captures only uncorrelated components, these uncorrelated components are not necessarily independent for general distributions. ICA minimizes the mutual information (or relative Kullback-Leibler divergence) of non-Gaussian data because two distributions with zero mutual information are statistically independent [28].
For finding the optimal number of IC, we apply ICA to the last layer feature space training set vectors. Using different number of IC and comparing the classification performance of the original model with the obtained using a linear combination of the reduced number of calculated components, it is possible to find the optimal number of components (n) that does not significantly reduce the classification performance of the original model. We modify the original model adding a new layer after the last layer feature space to calculate online the components of every analyzed image. This layer is a linear transformation and acts as a dimensionality reduction layer (see fig. 1). The final classification is achieved linearly combining the low dimensional IC layer.   After identifying the optimal n, we use the receptive field and pixel-wise explanation model [29] to visualize the independent scores in the input space. In this way we are visualizing not only a score map explaining a classification but also differentiating and visualizing the mathematically IC responsible of a particular classification.

Mathematical formalization
Let F train be the set of all feature vectors of training set: being T the number of elements of the training set, and m the dimension of feature space vector.
Let S train the set of IC calculated from F train : being n the number of IC, n < m. Every s (i) can be expressed as a linear combination of f (i) : Where W is calculated using a optimization method, minimizing the mutual information between n IC (min S I(S)) [30].
The ordinal regression problem solved using F train can be expressed as: being: being c (i) = Af (i) the predicted class vector and κ val the evaluation function used in the training process of the neural network calculated for the validation set.
We want to solve the same problem using the reduced ICA components (feature space compressed version): being c (i) = Bs (i) the predicted class vector using n IC. The optimal number of IC to select is the one that minimizes the difference in performance between both models: 4. Results

Data
In this study we use the EyePACS dataset of the Diabetic Retinopathy Detection competition hosted on the internet Kaggle Platform. For every patient right and left eye images are reported. All images are classified by ophthalmologists according to the standard severity scale presented before in [31]. The images are taken in variable conditions: by different cameras, illumination conditions and resolutions. The

Baseline Model
Our baseline model uses a 3x640x640 input image obtained from a minimal preprocessing step where only the external background borders are trimmed and later resized to the required input size. It is a CNN of 391,325 parameters, divided in 17 layers. Layers are divided in two groups: the feature extractor and the classifier. The feature extraction has 7 blocks of 2 layers. Every layer is a stack of a 3x3 convolution with stride 1x1 and padding 1x1 followed by a batch normalization and a ReLU activation function. Between every block a 2x2 max-pooling operation of stride 2x2 is applied. Afterwards, the classification phase takes place using a 2x2 convolution. A 4x4 average-pooling reduces the dimensionality to get a final 64 feature vector that are linearly combined to obtain the output scores of every class. A soft-max function allows the conversion of scores to probabilities to feed the values to the proper cost function during the optimization process. The feature extractor has 16 filters in the first block, 32 in the second and 64 in all the other.

Model Modifications
For all the training set we calculate the last layer feature space, obtaining a 64-dimensional vector as a representation of each image. We observe that this vector is highly redundant. After applying a PCA, with only 10 components it is possible to explain 99% of the variance. Using the 64-dimensional feature-space vector (m = 64) of all the training set, we calculate a set of ICAs using different n values. With each one, we train a linear classifier to calculate the evaluation metric obtained over a validation set ( fig. 1b). We choose the minimal n that allows achieving maximum performance. The optimal n for this problem is 3, achieving a QW K val = 0.790 not far from the achieved by the original model without dimensionality reduction (QW K val = 0.800). Fig. 2 shows the contribution of each component to the score of each class. We can see that ICA values are differentiating between extreme classes 0, 2 and 4. As we are training the network using as a loss function QW K, the optimization takes place as an ordinal regression. We expect the network to find the underlying causes present in the image that produce the predefined sorting of the classes, ie. the types of lesions and the number of them. Class 0 score contributions come from ICA 0 > 0, ICA 1 > 0 and ICA 2 > 0; being the class markers of the presence of disease ICA 0 < 0, ICA 1 < 0 and ICA 2 < 0. Analyzing the pixels with higher negative signals in the three components will give us the points that are contributing the most to the signaling of a possible presence of the disease. Backpropagating the scores of each one of this negative components will give a richer explanation with distinction between three possible independent causes of the final diagnostic given by the model.
We use a two-dimensional t-SNE visualization [32] for qualitative evaluation of the difference between the quality of the separation using the 64 feature vector and of the reduced version with only 3 ICA components. In fig. 3 we show the 2D t-SNE visualization for the original feature space and for the ICA reduced space. We can see how for both spaces class 0, 2 and 4 classes are clearly separated. Class 0 and 1 are not properly separated and in the case of 3 and 4, although the separation is not perfect, it is possible to distinguish a different location of both classes for both spaces. As the quantitative index QWK also show, there is not a qualitative separation capability difference between the t-SNE map of the original features and the map of the ICA 3 dimensional compressed feature map. For such 1.25% difference in the classification performance index (κ orig = 0.800 vs κ ICA = 0.790) we can conclude that the independent component analysis has been able to find a adequate compressed expression of the information encoded in the network.

Score components contribution for a test sample
We use a pixel-wise relevance propagation derived method for visualizing each ICA component independently. In this way, it is possible to visualize the mathematically independent contributions, enhancing the localization of different types of primary elements causing the disease. Figure 4 shows the intermediate score maps generated using a receptive field of 61x61. Figure 5 shows the same intermediate scores for a class 4 image. These hidden layer maps are useful when a general map of the lesion locations is enough.
Input-space pixel scores are the most suitable when pixel detail is required for detecting the individual lesions causing the disease.    Figure 6 show the input space final scores for ICA 2 . We show the points having contributions higher than a prestablished limit, in this case with values higher than two standard deviations. ICA 2 score maps have an almost perfect match with real lesions having a very low rate of false positives and false negatives. ICA 0 and ICA 1 include statistical regularities not directly related with relevant lesions. From such interpretation we conclude that ICA acts as a filter separating lesion information present in images from blink artifacts, noise and other statistical regularities present in images.

Conclusions
In this paper we studied some of the feature space properties of a already trained diabetic retinopathy deep learning classifier of 17 layers. We saw how, even being a small network, the redundancy of the feature space is very high. We hypothesized that if the network is able to achieve human-level classifications, in some way, it has been able to identify the important statistical regularities present in the image that are important for the classification. As last layer is a linear classifier, such properties has to be disentagled, ie. expressed in a linear way in such last layer. We applied a independent component analysis over such last layer with the method explained in the paper in order to find the optimal number of components that maximize the classification capabilities of such compressed version of the features. We experimentally proved that reducing to only 3 components is possible to achieve almost 99% of the evaluation metric. Such value experimentally prove that the ICA compression has been able to extract all the important elements required for the classification. With such elements, we applied a visualization technique for identifying the elements in the image of each one of the components, using a pixel-wise derived visualization technique. We are able to generate three maps for every image each one identifying independent statistical regularities important for the classification.
Our method allows not only the classification of retinographies but also the identification and localization in the image of the mathematically independent signs of the disease. The presented ICA score model is of general applicability and can be easily adapted for the usage in other image classification tasks thus, as future work we plan to test it in other types medical images (eg. cancer detection in mammograms).


text:
A Modern Take on the Bias-Variance Tradeoff in Neural Networks

The bias-variance tradeoff tells us that as model complexity increases, bias falls and variances increases, leading to a U-shaped test error curve. However, recent empirical results with over-parameterized neural networks are marked by a striking absence of the classic U-shaped test error curve: test error keeps decreasing in wider networks. This suggests that there might not be a bias-variance tradeoff in neural networks with respect to network width, unlike was originally claimed by, e.g., Geman et al. (1992). Motivated by the shaky evidence used to support this claim in neural networks, we measure bias and variance in the modern setting. We find that both bias and variance can decrease as the number of parameters grows. To better understand this, we introduce a new decomposition of the variance to disentangle the effects of optimization and data sampling. We also provide theoretical analysis in a simplified setting that is consistent with our empirical findings.

INTRODUCTION
There is a dominant dogma in machine learning: "The price to pay for achieving low bias is high variance" (Geman et al., 1992).
The quantities of interest here are the bias and variance of a learned model's prediction on a new input, where the randomness comes from the sampling of the training data. This idea that bias decreases while variance increases with model capacity, leading to a U-shaped test error curve is commonly known as the bias-variance tradeoff ( Fig. 1 (left)).
There exist experimental evidence and theory that support the idea of a tradeoff. In their landmark paper, Geman et al. (1992) measure bias and variance in various models. They show convincing experimental evidence for the bias-variance tradeoff in nonparametric methods such as kNN (k-nearest neighbor) and kernel regression. They also show experiments on neural networks and claim that bias decreases and variance increases with network width. Statistical learning theory (Vapnik, 1998) successfully predicts these U-shaped test error curves implied by a tradeoff for a number of classic machine learning models. A key element is identifying a notion of model capacity, understood as the main parameter controlling this tradeoff.
A number of different research directions have spawned in response to these findings. Neyshabur et al. (2015) hypothesize the existence of an implicit regularization mechanism. Some study the role that optimization plays (Soudry et al., 2018;Gunasekar et al., 2018). Others suggest new measures of capacity (Liang et al., 2017;Neyshabur et al., 2019). All approaches focus on test error, rather than studying bias and variance directly (Neyshabur et al., 2019;Geiger et al., 2019;Liang et al., 2017;Belkin et al., 2019a).
Test error analysis does not give a definitive answer on the lack of a bias-variance tradeoff. Consider boosting: it is known that its test error often decreases with arXiv:1810.08591v4 [cs.LG] 18 Dec 2019  (Fortmann-Roe, 2012). We find that both bias and variance decrease when we increase network width on MNIST (right) and other datasets (Section 4). These results seem to contradict the traditional intuition of a strict tradeoff.
the number of rounds (Schapire & Singer, 1999, Figures 8-10). In spite of this monotonicity in test error, Bühlmann & Yu (2003) show that variance grows at an exponentially decaying rate, calling this an "exponential bias-variance tradeoff." To study the bias-variance tradeoff, one has to isolate and measure bias and variance individually. To the best of our knowledge, there has not been published work reporting such measurements on neural networks since Geman et al. (1992).
We go back to basics and study bias and variance. We start by taking a closer look at Geman et al. (1992, Figure 16 and Figure 8 (top))'s experiments with neural networks. We notice that their experiments do not support their claim that "bias falls and variance increases with the number of hidden units." The authors attribute this inconsistency to convergence issues and maintain their claim that the bias-variance tradeoff is universal. Motivated by this inconsistency, we perform a set of bias-variance experiments with modern neural networks.
We measure prediction bias and variance of fully connected neural networks. These measurements allow us to reason directly about whether there exists a tradeoff with respect to network width. We find evidence that both bias and variance can decrease at the same time as network width increases in common classification and regression settings ( Fig. 1 and Section 4).
We observe the qualitative lack of a bias-variance tradeoff in network width with a number of gradientbased optimizers. In order to take a closer look at the roles of optimization and data sampling, we propose a simple decomposition of total prediction variance (Section 3.3). We use the law of total variance to get a term that corresponds to average (over data samplings) vari-ance due to optimization and a term that corresponds to variance due to training set sampling of an ensemble of differently initialized networks. Variance due to optimization is significant in the under-parameterized regime and monotonically decreases with width in the over-parameterized regime. There, total variance is much lower and dominated by variance due to sampling (Fig. 2).
We provide theoretical analysis, consistent with our empirical findings, in simplified analysis settings: i) prediction variance does not grow arbitrarily with number of parameters in fixed-design linear models; ii) variance due to optimization diminishes with number of parameters in neural networks under strong assumptions.
Organization The rest of this paper is organized as follows. We discuss relevant related work in Section 2. Section 3 establishes necessary preliminaries, including our variance decomposition. In Section 4, we empirically study the impact of network width on variance.
In Section 5, we present theoretical analysis in support of our findings. Neyshabur et al. (2015); Neyshabur (2017) point out that because increasing network width does not lead to a U-shaped test error curve, there must be some form of implicit regularization controlling capacity. Our work is consistent with this finding, but by approaching the problem from the bias-variance perspective, we gain additional insights: 1) We specifically address the hypothesis that decreased bias must come at the ex- pense of increased variance (see Geman et al. (1992) and Appendix E) by measuring both quantities. 2) Our more fine-grain approach reveals that variance due to optimization vanishes with width, while variances due to sampling increases and levels off. This insight about variance due to sampling is consistent with existing variance results for boosting (Bühlmann & Yu, 2003). To ensure that we are studying networks of increasing capacity, one of the experimental controls we use throughout the paper is to verify that bias is decreasing.

RELATED WORK
In independent concurrent work, Spigler et al. (2018) is the only one to theoretically analyze variance. Their work differs from ours in that we run experiments with neural networks on complex, real data, while they carry out a theoretical analysis of linear models in a simplified teacher (data generating distribution) setting.

Set-up
We consider the typical supervised learning task of predicting an output y ∈ Y from an input x ∈ X , where the pairs (x, y) are drawn from some unknown joint distribution, D. The learning problem consists of learning a function h S : X → Y from a finite training dataset S of m i.i.d. samples from D. The quality of a predictor h can quantified by the expected error, for some loss function : Y × Y → R.
In this paper, predictors h θ are parameterized by the weights θ ∈ R N of neural networks. We consider the average performance over possible training sets (denoted by the random variable S) of size m. This is the same quantity Geman et al. (1992) consider. While S is the only random quantity studied in the traditional bias-variance decomposition, we also study randomness coming from optimization. We denote the on arXiv five months after our paper appeared on arXiv.
random variable for optimization randomness (e.g. initialization) by O.
Formally, given a fixed training set S and fixed optimization randomness O, the learning algorithm A produces θ = A(S, O). Randomness in optimization translates to randomness in A(S, ·). Given a fixed training set, we encode the randomness due to O in a conditional distribution p(θ|S). Marginalizing over the training set S of size m gives a marginal distribution p(θ) = E S p(θ|S) on the weights learned by A from m samples. In this context, the average performance of the learning algorithm using training sets of size m can be expressed in the following ways:

Bias-variance decomposition
We briefly recall the standard bias-variance decomposition in the case of squared-loss. We work in the context of classification, where each class k ∈ {1 · · · K} is represented by a one-hot vector in R K . The predictor outputs a score or probability vector in R K . In this context, the risk in Eq.
(2) decomposes into three sources of error (Geman et al., 1992): The first term is an intrinsic error term independent of the predictor; the second is a bias term: whereȳ(x) denotes the expectation E[y|x] of y given x. The third term is the expected variance of the output predictions: where the expectation over θ can be done as in Eq.
(2). Interpreting this bias-variance decomposition as a bias-variance tradeoff is quite pervasive (see, e.g., Hastie et al. (2001, Chapter 2.9), Goodfellow et al. (2016, Bishop (2006, Chapter 3.2)). It is generally invoked to emphasize that the model selected should be of the complexity that achieves the optimal balance between bias and variance.
Note that risks computed with classification losses (e.g cross-entropy or 0-1 loss) do not have such a clean bias-variance decomposition (Domingos, 2000;James, 2003). However, it is natural to expect that bias and variance are useful indicators of the performance of models that are not assessed with squared error. In fact, we show the classification risk can be bounded as 4 times the regression risk in Appendix D.4. To empirically examine this connection, in all of our graphs that have "test error" or "training error" on some classification task, we plot the 0-1 classification error (see, e.g., Fig. 3b).

Further decomposing variance into its sources
In the set-up of Section 3.1 the prediction is a random variable that depends on two sources of randomness: the randomly drawn training set, S, and any optimization randomness, O, encoded into the conditional p(·|S). In certain regimes, one gets significantly different predictions when using a different initialization. Similarly, the output of a learned predictor changes when we use a different training set. How do we start disentangling variance caused by sampling from variance caused by optimization? There are few different ways; here we describe one of them.
Our goal is to measure prediction variance due to sampling, while controlling for the effect of optimization randomness.
Definition 1 ((Ensemble) Variance due to sampling). We consider the variance of an ensemble of infinitely many predictors with different optimization randomness (e.g. random initializations): A common practice to estimate variance due to optimization effects is to run multiple seeds on a fixed training set. Definition 2 ((Mean) Variance due to optimization). We consider the average (over training sets) variance over optimization randomness for a fixed training set: The law of total variance naturally decomposes variance into these very terms: We use this decomposition to get a finer understanding of our observations (Fig. 2).

EXPERIMENTS
In this section, we study how variance of fully connected single hidden layer networks varies with width. We provide evidence against Geman et al.  Our main finding is that, for all tasks that we study, bias and variance both decrease as we scale network width. We also provide a meaningful decomposition of prediction variance into a variance due to sampling term and a variance due to optimization term.

Common experimental details
We run experiments on different datasets: MNIST, SVHN, CIFAR10, small MNIST, and a sinusoid regression task. Averages over data samples are performed by taking the training set S and creating 50 bootstrap replicate training sets S by sampling with replacement from S. We train 50 different neural networks for each hidden layer size using these different training sets. Then, we estimate E bias 2 and E variance as in Section 3.2, where the population expectation E x is estimated with an average over the test set. To estimate the two terms from the law of total variance (Equation 4), we use 10 random seeds for the outer expectation and 10 for the inner expectation, resulting in a total of 100 neural networks for each hidden layer size. Furthermore, we compute 99% confidence intervals for our bias and variance estimates using the bootstrap (Efron, 1979).
The networks are initialized using PyTorch's default initialization, which scales the variance of the weight initialization distribution inversely proportional to the width (LeCun et al., 1998;Glorot & Bengio, 2010). The networks are trained using SGD with momentum and generally run for long after 100% training set accuracy is reached (e.g. 500 epochs for full data MNIST and 10000 epochs for small data MNIST). The overall trends we find are robust to how long the networks are trained after the training error converges. The step size hyperparameter is specified in each of the sections, and the momentum hyperparameter is always set to 0.9. To make our study as general as possible, we consider networks without regularization bells and whistles such as weight decay, dropout, or data augmentation, which Zhang et al. (2017) found to not be necessary for good generalization.

Decreasing variance in full data setting
We find a clear decreasing trend in variance with width of the network in the full data MNIST setting ( Fig. 1). We also see the same trend with CIFAR10 (Appendix B.1) and SVHN (Appendix B.2). In these experiments, the same step size is used for all networks for a given dataset (0.1 for MNIST and 0.005 for CI-FAR10 and SVHN). The trend is the same with or without early stopping, so early stopping is not necessary to see decreasing variance, similar to how it was not necessary to see better test set performance with width in Neyshabur et al. (2015). Wider ResNets are known to achieve lower test error (Zagoruyko & Komodakis, 2016); this likely translates to decreasing variance with width in convolutional networks as well. Much of the over-parameterization literature focuses on over-parameterization in width; interestingly, the variance trend is not the same when varying depth (Appendix C).

Testing the limits: decreasing variance in the small data setting
Decreasing the size of the dataset can only increase variance. To study the robustness of the above obser- vation, we decrease the size of the training set to just 100 examples. In this small data setting, somewhat surprisingly, we still see that both bias and variance decrease with width ( Fig. 3a). The test error behaves similarly ( Fig. 3b). Because performance is more sensitive to step size in the small data setting, the step size for each network size is tuned using a validation set (see Appendix B.4 for step sizes). The training for tuning is stopped after 1000 epochs, whereas the training for the final models is stopped after 10000 epochs. Note that because we see decreasing bias with width, effective capacity is, indeed, increasing while variance is decreasing.
One control that motivates the experimental design choice of optimal step size is that it leads to the conventional decreasing bias trend ( Fig. 3a) that indicates increasing effective capacity. In fact, in the corresponding experiment where step size is the same 0.01 for all network sizes, we do not see monotonically decreasing bias (Appendix B.5).
This sensitivity to step size in the small data setting is evidence that we are testing the limits of our hypothesis. By looking at the small data setting, we are able to test our hypothesis when the ratio of size of network to dataset size is quite large, and we still find this decreasing trend in variance (Fig. 3a).
To see how dependent this phenomenon is on SGD, we also run these experiments using batch gradient descent and PyTorch's version of LBFGS. Interestingly, we find a decreasing variance trend with those optimizers as well. These experiments are included in Appendix B.6.

Decoupling variance due to sampling from variance due to optimization
In order to better understand this variance phenomenon in neural networks, we separate the variance due to sampling from the variance due to optimiza-tion, according to the law of total variance (Equation 4). Contrary to what traditional bias-variance tradeoff intuition would suggest, we find variance due to sampling increases slowly and levels off, once sufficiently over-parameterized (Fig. 2). Furthermore, we find that variance due to optimization decreases with width, causing the total variance to decrease with width ( Fig. 2).
A body of recent work has provided evidence that over-parameterization (in width) helps gradient descent optimize to global minima in neural networks (Du et al., 2019;Du & Lee, 2018;Soltanolkotabi et al., 2017;Livni et al., 2014;Zhang et al., 2018). Always reaching a global minimum implies low variance due to optimization on the training set. Our observation of decreasing variance on the test set shows that the over-parameterization (in width) effect on optimization seems to extend to generalization, on the data sets we consider.

Visualization with regression on sinusoid
We trained different width neural networks on a noisy sinusoidal distribution with 80 independent training examples. This sinusoid regression setting also exhibits the familiar bias-variance trends (Fig. 3c) and trends of the two components of the variance and the test error ( Fig. 19 of Appendix B.7).
Because this setting is low-dimensional, we can visualize the learned functions. The classic caricature of high capacity models is that they fit the training data in a very erratic way (example in Fig. 15 of Appendix B.7). We find that wider networks learn sinusoidal functions that are much more similar than the functions learned by their narrower counterparts (Fig. 4). We have analogous plots for all of the other widths and ones that visualize the variance similar to how it is commonly visualized for Gaussian processes in Appendix B.7. THEORETICAL INSIGHTS Our empirical results demonstrate that in the practical setting, variance due to optimization decreases with network width while variance due to sampling increases slowly and levels off once sufficiently overparameterized. In Section 5.1, we discuss the simple case of linear models and point out that non-increasing variance can already be seen in the over-parameterized setting. In Section 5.2 we take inspiration from linear models to provide arguments for the behavior of variance in increasingly wide neural networks, and we discuss the assumptions we make.

Insights from linear models
In this section, we review the classic result that the variance of a linear model grows with the number of parameters (Hastie et al., 2009, Section 7.3) and point out that variance behaves differently in the overparameterized setting.
We consider least-squares linear regression in a standard setting which assumes a noisy linear mapping y = θ T x + between input feature vectors x ∈ R N and real outputs, where denotes the noise random variable with E[ ] = 0 and Var( ) = σ 2 . In this context, the over-parameterized setting is when the dimension N of the input space is larger than the number m of examples.
Let X denote the m × N design matrix whose i th row is the training point x T i , let Y denote the corresponding labels, and let Σ = X T X denote the empirical covariance matrix. We consider the fixed-design setting where X is fixed, so all of the randomness due to data sampling comes solely from . A learns weightsθ from (X, Y ), either by a closed-form solution or by gradient descent, using a standard initialization The predictor makes a prediction on x ∼ D: h(x) =θ T x. Then, the quantity we care about is E x Var(h(x)).

Under-parameterized setting
The case where N ≤ m is standard: if X has maximal rank, Σ is invertible; the solution is independent of the initialization and given byθ = Σ −1 X T Y . All of the variance is a result of randomness in the noise . For a fixed x, This grows with the number of parameters N . For example, taking the expected value over the empirical distribution,p, of the sample, we recover that the variance grows with N : We provide a reproduction of the proofs in Appendix D.1.

Over-parameterized setting
The over-parameterized case where N > m is more interesting: even if X has maximal rank, Σ is not invertible. This leads to a subspace of solutions, but gradient descent yields a unique solution from updates that belong to the span of the training points x i (row space of X) (LeCun et al., 1991), which is of dimension r = rank(X) = rank(Σ). Correspondingly, no learning occurs in the null space of X, which is of dimension N − r. Therefore, gradient descent yields the solution that is closest to initialization: where P ⊥ projects onto the null space of X and + denotes the Moore-Penrose inverse.
The variance has two contributions: one due to initialization and one due to sampling (here, the noise ), as in Eq. (4). These are made explicit in Proposition 1.
Proposition 1 (Variance in over-parameterized linear models). Consider the over-parameterized setting where N > m. For a fixed x, the variance decomposition of Eq. (4) yields This does not grow with the number of parameters N . In fact, because Σ −1 is replaced with Σ + , the variance scales as the dimension of the data (i.e the rank of X), as opposed to the number of parameters. For example, taking the expected value over the empirical distribution,p, of the sample, we obtain where r = rank(X). We provide the proofs for overparameterized linear models in Appendix D.2.

A more general result
We will illustrate our arguments in the following simplified setting, where M, M ⊥ , and d(N ) are the more general analogs of rowspace(X), nullspace(X), and r (respectively): Setting. Let N be the dimension of the parameter space. The prediction for a fixed example x, given by a trained network parameterized by θ depends on: (i) a subspace of the parameter space, M ∈ R N with relatively small dimension, d(N ), which depends only on the learning task.
(ii) parameter components corresponding to directions orthogonal to M. The orthogonal M ⊥ of M has dimension, N − d(N ), and is essentially irrelevant to the learning task.
We can write the parameter vector as a sum of these two components θ = θ M + θ M ⊥ . We will further make the following assumptions.
Assumption 1 The optimization of the loss function is invariant with respect to θ M⊥ .
Assumption 2 Regardless of initialization, the optimization method consistently yields a solution with the same θ M component (i.e. the same vector when projected onto M).

Variance due to initialization
Given the above assumptions, the following result shows that the variance from initialization 3 vanishes as we increase N . The full proof, which builds on concentration results for Gaussians (based on Levy's lemma (Ledoux, 2001)), is given in Appendix D.3.
Theorem 1 (Decay of variance due to initialization).
Consider the setting of Section 5.2 Let θ denote the parameters at the end of the learning process. Then, for a fixed data set and parameters initialized as θ 0 ∼ N (0, 1 N I), the variance of the prediction satisfies the inequality, where L is the Lipschitz constant of the prediction with respect to θ, and for some universal constant C > O.
This result guarantees that the variance decreases to zero as N increases, provided the Lipschitz constant L grows more slowly than the square root of dimension,

Variance due to sampling
Under the above assumptions, the parameters at the end of learning take the form θ = θ * M + θ 0M ⊥ . For Discussion on assumptions We made strong assumptions, but there is some support for them in the literature. The existence of a subspace M ⊥ in which no learning occurs was also conjectured by Advani & Saxe (2017) and shown to hold in linear neural networks under a simplifying assumption that decouples the dynamics of the weights in different layers. Li et al. (2018) empirically showed the existence of a critical number d(N ) = d of relevant parameters for a given learning task, independent of the size of the model. Sagun et al. (2017) showed that the spectrum of the Hessian for over-parameterized networks splits into (i) a bulk centered near zero and (ii) a small number of large eigenvalues; and Gur-Ari et al. (2018) recently gave evidence that the small subspace spanned by the Hessian's top eigenvectors is preserved over long periods of training. These results suggest that learning occurs mainly in a small number of directions.

CONCLUSION AND FUTURE WORK
We provide evidence against Geman et al. (1992)'s claim that "the price to pay for achieving low bias is high variance," finding that both bias and variance decrease with network width. Geman et al. (1992)'s claim is found throughout machine learning and is meant to generally apply to all of machine learning (Appendix E), and it is correct in many cases (e.g. kNN, kernel regression, splines). Is this lack of a tradeoff specific to neural networks or is it present in other models as well such as decision trees?
We propose a new decomposition of the variance, finding variance due to sampling (analog of regular variance in simple settings) does not appear to be dependent on width, once sufficiently over-parameterized, and that variance due to optimization decreases with width. By taking inspiration from linear models, we perform a theoretical analysis of the variance that is consistent with our empirical observations.
We view future work that uses the bias-variance lens as promising. For example, a probabilistic notion of effective capacity of a model is natural when studying generalization through this lens (Appendix A). We did not study how bias and variance change over the course of training; that would make an interesting direction for future work. We also see further theoretical treatment of variance as a fruitful direction for better understanding complexity and generalization abilities of neural networks.
where S is a training set drawn from D m , h(A, S) is a random variable drawn from the distribution over learned functions induced by D and the randomness in A; p is the corresponding density. Thinking about a model'shypothesis class can lead to drastically different intuitions for the complexity of a model and its variance (Fig. 5). This is at the core of the intuition for why the traditional view of bias-variance as a tradeoff does not hold in all cases.

Bias and Variance
Bias Variance (a) Variance decreases with width, even in the small data setting (SGD). This figure is in the main paper, but we include it here to compare with the corresponding step sizes used. Note that the U curve shown in Fig. 12 when we do not tune the step size is explained by the fact that the constant step chosen is a "good" step size for some networks and "bad" for others. Results from Keskar et al.  2017) show that a step size that corresponds well to the noise structure in SGD is important for achieving good test set accuracy. Because our networks are different sizes, their stochastic optimization process will have a different landscape and noise structure. By tuning the step size, we are making the experimental design choice to keep optimality of step size constant across networks, rather than keeping step size constant across networks. To us, choosing this control makes much more sense than choosing to control for step size. Note that Park et al. (2019) show that, as long as the network size is not too big, controlling for "optimality of step size" and "keeping step size constant" with increasing network width actually correspond to the same thing, as long as the networks are not too big, which fits well with the fact that we used the same step size across network widths for almost all of our experiments.
B.6 Other optimizers for width experiment on small data MNIST

Average Error
Test Error Train Error Figure 13: Variance decreases with width in the small data setting, even when using batch gradient descent.  Figure 15: Caricature examples of high variance learners on sinusoid task. Below, we find that this does not happen with increasingly wide neural networks ( Fig. 17 and Fig. 18).   Average Test Loss Test Loss Figure 19: We observe the same trends of decomposed variance (left) and test error (right) in the sinusoid regression setting.  Figure 20: We can see that, when using dynamical isometry, bias decreases with depth and variance slowly increases with depth (left). This increase in variance is so small that it only translates to a an increase in test error of about 0.1% for depth 25 to depth 200 (right).

C.2 Discussion on need for careful experimental design
Depth is an important component of deep learning. We study its effect on bias and variance by fixing width and varying depth. However, there are pathological problems associated with training very deep networks such as vanishing/exploding gradient (Hochreiter, 1991;Bengio et al., 1994;Glorot & Bengio, 2010), signal not being able to propagate through the network (Schoenholz et al., 2017), and gradients resembling white noise (Balduzzi et al., 2017). He et al. (2016) pointed out that very deep networks experience high test set error and argued it was due to high training set loss. However, while skip connections (He et al., 2016), better initialization (Glorot & Bengio, 2010), and batch normalization (Ioffe & Szegedy, 2015) have largely served to facilitate low training loss in very deep networks, the problem of high test set error still remains.
The current best practices for achieving low test error in very deep networks arose out of trying to solve the above problems in training. An initial step was to ensure the mean squared singular value of the input-output Jacobian, at initialization, is close to 1 (Glorot & Bengio, 2010). More recently, there has been work on a stronger condition known as dynamical isometry, where all singular values remain close to 1 (Saxe et al., 2014;Pennington et al., 2017). Pennington et al. (2017) also empirically found that dynamical isometry helped achieve low test set error. Furthermore, Xiao et al. (2018, Figure 1) found evidence that test set performance did not degrade with depth when they lifted dynamical isometry to CNNs. This why we settled on dynamical isometry as the best known practice to control for as many confounding factors as possible.
We first ran experiments with vanilla full connected networks (Fig. 21). These have clear training issues where networks of depth more than 20 take very long to train to the target training loss of 5e-5. The bias curve is not even monotonically decreasing. Clearly, there are important confounding factors not controlled for in this simple setting. Still, note that variance increases roughly linearly with depth.
We then study fully connected networks with skip connections between every 2 layers (Fig. 22). While this allows us to train deeper networks than without skip connections, many of the same issues persist (e.g. bias still not monotonically decreasing). The bias, variance, and test error curves are all checkmark-shaped.  Figure 22: While the addition of skip connections (between every other layer) might push the bottom of the U curve in test error out to 10 skip connections (21 layers), which is further than where the bottom was observed without skip connections (3 layers), test error still degrades noticeably in greater depths. Additionally, bias still does not even monotonically decrease with depth. While skip connections appear to have helped control for the factors we want to control, they were not completely satisfying. All networks have training error 0 and are trained to the same training loss of 5e-5.

C.5 Dynamical isometry depth experiments
The figures in this section are included in the main paper, but they are included here for comparison to the above and for completeness. Proof. For a fixed x, we have h(x) = x Tθ . Takingθ = Σ −1 X T Y to be the gradient descent solution, and using Y = Xθ + , we obtain: x T θ, and the variance is, Taking the expected value over the empirical distribution,p, of the sample, we find an explicit increasing dependence on N : Proof. By the law of total variance, Here have h(x) = x Tθ , whereθ the gradient descent solutionθ = P ⊥ (θ 0 ) + Σ + X T Y , and θ 0 ∼ N (0, 1 N I). Then, Since E θ0 (h(x)) = x T Σ + X T Y , the calculation of Var (E θ0 )h(x)) is similar as in D.1, where Σ −1 is replaced by Taking the expected value over the empirical distribution,p, of the sample, we find an explicit dependence on r = rank(X), not N : where I + r denotes the diagonal matrix with 1 for the first r diagonal elements and 0 for the remaining N − r elements.

D.3 Proof of Theorem 1
First we state some known concentration results (Ledoux, 2001) that we will use in the proof. Lemma 1 (Levy). Let h : S n R → R be a function on the n-dimensional Euclidean sphere of radius R, with Lipschitz constant L; and θ ∈ S n R chosen uniformly at random for the normalized measure. Then for some universal constant C > 0.
Uniform measures on high dimensional spheres approximate Gaussian distributions (Ledoux, 2001). Using this, Levy's lemma yields an analogous concentration inequality for functions of Gaussian variables: Lemma 2 (Gaussian concentration). Let h : R n → R be a function on the Euclidean space R n , with Lipschitz constant L; and θ ∼ N (0, σI n ) sampled from an isotropic n-dimensional Gaussian. Then: for some universal constant C > 0.
Note that in the Gaussian case, the bound is dimension free.
In turn, concentration inequalities give variance bounds for functions of random variables. Corollary 1. Let h be a function satisfying the conditions of Theorem 2, and . Then Var(h) = Var(g) and Now swapping expectation and integral (by Fubini theorem), and by using the identity E1 |g|>t = P(|g| > t), we obtain We are now ready to prove Theorem 1. We first recall our assumptions: Assumption 1. The optimization of the loss function is invariant with respect to θ M⊥ .
Assumption 2. Along M, optimization yields solutions independently of the initialization θ 0 .
We add the following assumptions.
Assumption 3. The prediction h θ (x) is L-Lipschitz with respect to θ M⊥ .
Assumption 4. The network parameters are initialized as We first prove that the Gaussian concentration theorem translates into concentration of predictions in the setting of Section 5.2.1.
Theorem 2 (Concentration of predictions). Consider the setting of Section 5.2 and Assumptions 1 and 4. Let θ denote the parameters at the end of the learning process. Then, for a fixed data set, S we get concentration of the prediction, under initialization randomness, for some universal constant C > 0.
Proof. In our setting, the parameters at the end of learning can be expressed as where θ * M is independent of the initialization θ 0 . To simplify notation, we will assume that, at least locally around θ * M , M is spanned by the first d(N ) standard basis vectors, and M ⊥ by the remaining N − d(N ). This will allow us, from now on, to use the same variable names for θ M and θ M ⊥ to denote their lower-dimensional representations of dimension d(N ) and N − d(N ) respectively. More generally, we can assume that there is a mapping from θ M and θ M ⊥ to those lower-dimensional representations.
From Assumptions 1 and 4 we get Let g(θ M ⊥ ) h θ * M +θ M ⊥ (x). By Assumption 3, g(·) is L-Lipschitz. Then, by the Gaussian concentration theorem we get, The result of Theorem 1 immediately follows from Theorem 2 and Corollary 1, with σ 2 = 1/N : Provided the Lipschitz constant L of the prediction grows more slowly than the square of dimension, L = o( √ N ), we conclude that the variance vanishes to zero as N grows.

D.4 Bound on classification error in terms of regression error
In this section we give a bound on classification risk R classif in terms of the regression risk R reg .
Notation. Our classifier defines a map h : X → R k , which outputs probability vectors h(x) ∈ R k , with k y=1 h(x) y = 1. The classification loss is defined by where I(a) = 1 if predicate a is true and 0 otherwise. Given trained predictors h S indexed by training dataset S, the classification and regression risks are given by, where Y denotes the one-hot vector representation of the class y.
Proposition 2. The classification risk is bounded by four times the regression risk, R classif ≤ 4R reg .
Proof. First note that, if h(x) ∈ R k is a probability vector, then h(x) y < max y h(x) y =⇒ h(x) y < 1 2 By taking the expectation over x, y, we obtain the inequality L(h) ≤ L(h) where We then have, Appendix E Common intuitions from impactful works E.1 Quotes from influential papers "Neural Networks and the Bias/Variance Dilemma" from (Geman et al., 1992): "How big a network should we employ? A small network, with say one hidden unit, is likely to be biased, since the repertoire of available functions spanned by f (x; w) over allowable weights will in this case be quite limited. If the true regression is poorly approximated within this class, there will necessarily be a substantial bias. On the other hand, if we overparameterize, via a large number of hidden units and associated weights, then the bias will be reduced (indeed, with enough weights and hidden units, the network will interpolate the data), but there is then the danger of a significant variance contribution to the mean-squared error. (This may actually be mitigated by incomplete convergence of the minimization algorithm, as we shall see in Section 3.5.5.)" "An Overview of Statistical Learning Theory" from (Vapnik, 1999): "To avoid over fitting (to get a small confidence interval) one has to construct networks with small VC-dimension." "Stability and Generalization" from Bousquet & Elisseeff (2002): "It has long been known that when trying to estimate an unknown function from data, one needs to find a tradeoff between bias and variance. Indeed, on one hand, it is natural to use the largest model in order to be able to approximate any function, while on the other hand, if the model is too large, then the estimation of the best function in the model will be harder given a restricted amount of data." Footnote: "We deliberately do not provide a precise definition of bias and variance and resort to common intuition about these notions." "Understanding the Bias-Variance Tradeoff" from Fortmann-Roe (2012): "At its root, dealing with bias and variance is really about dealing with over-and under-fitting. Bias is reduced and variance is increased in relation to model complexity. As more and more parameters are added to a model, the complexity of the model rises and variance becomes our primary concern while bias steadily falls. For example, as more polynomial terms are added to a linear regression, the greater the resulting model's complexity will be." Here are two excerpts: • "As a general rule, as we use more flexible methods, the variance will increase and the bias will decrease. The relative rate of change of these two quantities determines whether the test MSE increases or decreases. As we increase the flexibility of a class of methods, the bias tends to initially decrease faster than the variance increases. Consequently, the expected test MSE declines. However, at some point increasing flexibility has little impact on the bias but starts to significantly increase the variance. When this happens the test MSE increases" (James et al., 2014, Chapter 2.2.2).
• "Our goal is to minimize the expected loss, which we have decomposed into the sum of a (squared) bias, a variance, and a constant noise term. As we shall see, there is a trade-off between bias and variance, with very flexible models having low bias and high variance, and relatively rigid models having high bias and low variance" (Bishop, 2006, Chapter 3.2).
• "As the model complexity of our procedure is increased, the variance tends to increase and the squared bias tends to decrease" (Hastie et al., 2001, Chapters 2.9).