Qingyuan Wang
email: qingyuan.wang@ucdconnect.ie

Antoine Frappé
email: †antoine.frappe@junia.com

Benoit Larras
email: benoit.larras@junia.com

Barry Cardiff
email: ‡barry.cardiff@ucd.ie

Deepu John
email: deepu.john@ucd.ie

Complexity Reduction of CNNs using Multi-Scale Group Convolution for IoT Edge Sensors

In this paper, we propose Multi-Scale Group Convolution (MSGC) an optimization to the conventional convolutional layer, to address the high computational complexity issue in deploying convolutional neural networks (CNN) on the Internet of Things (IoT) enabled edge sensors. The proposed method reduces complexity by grouping input channels of a convolution layer into smaller groups, thereby reducing the number of intermediate connections and complexity of matrix computations in a CNN. This approach results in a minor performance loss, which is compensated by utilizing a characteristic of group convolution to extract multi-scale features. The proposed technique is applied for detecting cardiac arrhythmias from electrocardiogram (ECG) data using CNNs to be deployed in edge sensors. For the binary classification of ECG into Normal or Anomalous beats, the proposed MSGC-based CNN achieved an average 30% reduction in computations while achieving similar or better performance compared to the conventional CNNs. We used the Physionet MIT-BIH Arrhythmia database for performance evaluation, and in the best scenario, our approach increases accuracy by 0.47%, F1 score by 1.87% while only using 64.41% MACs and 83.62% parameters. This optimization strategy can be extended to other CNN models where computational complexity reduction is critical for deployment in edge devices.

I. INTRODUCTION

Cardiovascular diseases (CVDs) are one of the leading healthcare challenges in terms of risks, costs, and mortality [START_REF] Organization | The World Health Report 2002: Reducing Risks, Promoting Healthy Life[END_REF]. The symptoms of CVDs can be sporadic, sudden, and can go unnoticed. Therefore continuous monitoring and a fast diagnosis are needed for CVD patients even outside the hospital premises. Deploying machine learning (ML) techniques on wearable sensors is widely considered a promising approach for continuous cardiac arrhythmia monitoring and detection. On-device edge intelligence using deep neural networks (DNNs) can enable long-term continuous health monitoring and early warnings so that proactive actions can be taken without causing privacy concerns.

Convolutional neural networks (CNN) are commonly used for CVD monitoring and detection [START_REF] Ebrahimi | A review on deep learning methods for ECG arrhythmia classification[END_REF] due to their high accuracies. CNN's inbuilt signal filtering and shift invariance feature make it more robust than perceptron networks. CNNs are also friendly to parallel computation and thus hardware implementations. Many studies employed 1D-CNN approaches with good detection performance [START_REF] Oh | Automated diagnosis of arrhythmia using combination of CNN and LSTM techniques with variable length heart beats[END_REF][START_REF] Yildirim | Arrhythmia detection using deep convolutional neural network with long duration ECG signals[END_REF]. Several 2D-CNN based techniques that process images generated from This publication has emanated from research supported in part by 1) a grant from Science Foundation Ireland under Grant number 18/CRT/6183 and 2) CHIST-ERA grant JEDAI CHIST-ERA-18-ACAI-003.

Conv Groups

with differnt sized kernels Regular Conv Layer MSGC Layer

Conv Operation Channel Shuffling

le n g th channels Fig. 1: A Typical CNN layer vs an MSGC layer transform domain features of ECG signals are also reported [START_REF]ECG arrhythmia classification using a 2-D convolutional neural network[END_REF][START_REF] Salem | ECG Arrhythmia Classification Using Transfer Learning from 2-Dimensional Deep CNN Features[END_REF]. All these works claim overall detection accuracy of above 90%. However, their models are usually too large and require enormous computational resources, preventing them from efficiently deploying in low-cost, low-power embedded IoT sensors. It is possible to reduce the model complexity by directly reducing the model size. Nevertheless, this may result in a significant performance drop and defeat the purpose. Model compression techniques are often used for the complexity reduction of DNNs so that they can be deployed on resource-limited devices. Common techniques include 1) Model pruning, 2) Quantization and 3) Knowledge distillation etc. Model pruning removes redundant and less significant weights from a trained neural network, reducing computations and model parameters. Quantization replaces trained weights with fewer-bit representations to limit the complexity of arithmetic operations. Knowledge distillation helps to train a simplified (smaller) model using the original (larger) model so that the smaller model can replace the larger one [START_REF] Cheng | A Survey of Model Compression and Acceleration for Deep Neural Networks[END_REF]. These methods sacrifice performance for efficiency and can be applied universally to any model. Domain-specific knowledge also could be used for model compression. For instance, MobileNet [START_REF] Howard | MobileNets: Efficient Convolutional Neural Networks for Mobile Vision Applications[END_REF] is an optimized CNN architecture designed for imaging tasks on mobile devices. Neural architecture search is a technique that automatically searches for the best architecture [START_REF]Neural Architecture Search: A Survey[END_REF]. These techniques use prior knowledge or simulation results to optimize the model and do not necessarily reduce performance.

In this paper, we propose an architecture optimization technique for CNNs and apply it for binary classification of ECG for CVD event detection. We present Multi-Scale Group Convolution (MSGC), an optimization for the convolutional layer in CNNs, for reducing the computational complexity for deployment in resource-limited devices.

The difference between a typical Convolutional and MSGC layer is illustrated in Fig. 1. MSGC is a combination of group convolution [START_REF] Krizhevsky | ImageNet Classification with Deep Convolutional Neural Networks[END_REF][START_REF] Xie | Aggregated Residual Transformations for Deep Neural Networks[END_REF], multi-scale feature extraction [START_REF] Dempster | ROCKET: Exceptionally fast and accurate time series classification using random convolutional kernels[END_REF], and channel shuffling [START_REF] Zhang | ShuffleNet: An Extremely Efficient Convolutional Neural Network for Mobile Devices[END_REF]. Group convolution reduces convolutional layer complexity by grouping input channels into smaller groups, thereby reducing the number of intermediate connections and complexity of matrix computations. This results in a minor performance loss, which is then compensated by the latter two techniques. While reducing complexity and model size, the proposed technique maintains a similar detection performance to that of a conventional CNNbased model. This technique is inspired by CNN optimization methods used in computer vision. The effect of this method is evaluated using the Physionet MIT-BIH Arrhythmia database. The details of MSGC are presented in Sec.II. Experimental details, Complexity comparisons with conventional CNNs, and Results are discussed in Sec.III.

II. MULTI-SCALE GROUP CONVOLUTION

The primary motivation of MSGC is to minimize computations while retaining the ability for feature extraction from the input. As Fig. 2 shows, MSGC combines 1) Group convolution, 2) Multi-scale feature extraction, and 3) Channel shuffling into a single layer to achieve this. Group convolution reduces computation by eliminating connections between different groups. This also allows for different kernel (filter) sizes in different groups, encouraging the model to view the input one different scales. After the convolutions, channel shuffling (a computation-free operation) is done to exchange information among multiple groups and create combinations of different kernel sizes with more layers stacking together.

Fig. 2: Components of Multi-Scale Group Convolution

A. Group Convolution

In group convolution, the input channels are divided into smaller groups, and convolution operation is done within each group. Fig. 2a shows an example of separated input channels, each of which has two independent output channels. This method was initially proposed with AlexNet [START_REF] Krizhevsky | ImageNet Classification with Deep Convolutional Neural Networks[END_REF] to distribute the model on two GPUs. While, the main advantage in our application is that it can reduce the convolution layer complexity and may not result in a performance drop [START_REF] Xie | Aggregated Residual Transformations for Deep Neural Networks[END_REF]. This is because group convolution works like an internal micro ensemble of multiple independent submodels, and this can compensate for losses compared to a single large model. However, this advantage requires that every submodel has enough ability to make an independent prediction. For very small models, the submodels from group convolution are less effective in making independent predictions, and therefore it may not improve and instead possibly reduce the performance. The following methods are used to compensate for any performance loss.

B. Multi-Scale Feature Extraction

The characteristics of useful features could vary slightly in every input sample. For ECG, the signal morphology can vary 1) temporally due to user activity, 2) measurement methods, 3) anatomical variations etc. These variations in signal characteristics make feature extraction less effective if we only use fixed kernels. DNNs are able to address this issue if the model size is big enough and with a sufficiently large training dataset. However, it is difficult for a size-limited model to achieve this by training only. Artificially introducing more prior knowledge, also known as applying inductive bias, is beneficial for feature extraction in this case. In our proposed method, we assign different kernel sizes for convolution groups, as this will enable the model to view the input signal on different scales. There are other successful precedents reported for this technique. In Rocket [START_REF] Dempster | ROCKET: Exceptionally fast and accurate time series classification using random convolutional kernels[END_REF], a large number of untrained convolution kernels were used to cover feature extraction from all possible input variations. This could be considered an extreme case of pre-defined multi-scale feature extraction. In our model, we assigned different kernel sizes to different groups of output channels to enable multi-scale feature extraction. Note that these multi-scale groups are not necessarily the same as the groups in the group convolution. This means different kernel sizes can be used even when the group convolution is not applied. Fig. 2b illustrates a case of one input channel being convolved in three different ways.

C. Channel Shuffling

The sparse channel connections due to group convolution operation make every group independent. However, after stacking layers, each channel will never communicate with other channels outside its group as they are not connected. One complexity-friendly method to establish communication between channels is shuffling them after group convolution, which was first employed in ShuffleNet [START_REF] Zhang | ShuffleNet: An Extremely Efficient Convolutional Neural Network for Mobile Devices[END_REF]. The shuffling operation can be implemented by re-indexing the channels and does not involve any computation. Furthermore, channel shuffling will create a mixture of various kernel sizes. In our implementation, channel shuffling is designed to make all combinations among groups. Hence, the shuffling is fixed and the number of channels must be the multiple of g 2 , where g is the number of groups. Fig. 2c shows a shuffling between four channels of two groups.

III. EXPERIMENTS

The MSGC-based model is evaluated by comparing it with a typical CNN model. We built several low-complexity networks based on the architecture in Fig. 3. We denote the size of input signal as [l in , c in ]. The network will change the channel numbers to c h while decreasing the length gradually. At the end of these models, a small fully-connected network with n h neurons will make the binary prediction. The complexity of these models are estimated in the number of multiply-accumulate (MAC) operation and number of parameters while their performances are evaluated in accuracy (ACC) and F1 score [START_REF] Ebrahimi | A review on deep learning methods for ECG arrhythmia classification[END_REF]. We used the MIT-BIH Arrhythmia database for performance evaluation. It contains 48 half-hours, two lead ECG records, sampled at 360 Hz. We used a single lead (MLII) signal as the input and extracted 240 sample points centered at the Rpeak location as a beat sample. The signal was normalized to follow the standard normal distribution to reduce any impact of amplitude variations across multiple records. The beat annotations are converted into two classes, i.e. normal and abnormal, to support binary classification. We divided the dataset into DS1 and DS2 as per the recommendation of the Association for the Advancement of Medical Instrumentation (AAMI) [START_REF] Ec57 | Testing and reporting performance results of cardiac rhythm and ST segment measurement algorithms[END_REF]. The records in DS1 are used for global training. From DS2, the first 5 minutes of each record is used for patient-specific training, and the rest are used for testing as per [START_REF] Ec57 | Testing and reporting performance results of cardiac rhythm and ST segment measurement algorithms[END_REF]. After the pre-processing steps, our experiments involve a total of 93,391 beats, about 11% of them are abnormal beats. ∼47k (50%), ∼7.4k (8%) and ∼3.9k (42%) beats are used for global training, patient-specific training and testing, respectively.

B. Experiment Details

The models are implemented with PyTorch, and their complexity is estimated using ptflops [START_REF] Sovrasov | Flops counter for convolutional networks in pytorch framework[END_REF]. All models are firstly trained with the global training set for 200 epochs. Then, for every record in DS2, the global model will be further trained using the first 5 minutes of data for 50 epochs and is evaluated on the rest of the record. The patient-specific model will be reset to the global trained state after evaluating each DS2 record. We applied some common training boosters such as an advanced optimizer, learning rate scheduler and data augmentation method to improve performance and stability during both training phases. AdamW [START_REF]Decoupled Weight Decay Regularization[END_REF] is used as the optimizer with a learning rate of 0.0003, and a cosine annealing scheduler is used to adjust the learning rate. Mixup [START_REF] Zhang | mixup: Beyond Empirical Risk Minimization[END_REF] is also applied to reduce the over-fitting issue.

Our experiments include six cases of different tiny models. The experiment starts from a model with 3 Conv layers and two fully-connected layers, as illustrated in Fig 3 . It has 9 hidden channels (denoted by c h ) for all Conv layers, and 64 hidden neurons (denoted by n h ) for fully-connected layers, i.e. c h = 9 and n h = 64. The purpose of the experiment is to validate if MSGC can be a replacement for conventional convolution in our application. We also tested some network variations with more channels and layers to evaluate MSGC's generalisability. Specifically, the width of networks can be extended to 18 hidden layers and 128 hidden neurons, i.e. c h = 18 and n h = 128. The number of Conv layers (denoted by M ) can also be increased to 5 or 7.

We also tested the components of MSGC separately to verify our hypothesis that the three methods can work together. Specifically, we tested four types of convolution: 1) group convolution alone, 2) MSGC without channel shuffling, 3) complete MSGC, and 4) conventional convolution. Deep learning algorithms are not deterministic and may sometimes be unstable when the model size and training dataset are small. Therefore, we repeated all experiments with different seeds three times to increase the results' reliability. The metrics in Tab.I are the average of those three repeats. Fig. 4: Results comparison between regular CNN and MSGC Fig. 4 visualizes four test cases comparing between MSGC and conventional convolution. The blue bar shows the relative MACs between displayed cases. A considerable saving of MACs is observed in every case. The green and red bars represent the accuracy and F1 score. The performance is nearly unchanged between the two convolution types. Meanwhile, we can see a trend of more complexity results in better F1 score. Tab.I compares the simulation results of six test cases. The proposed method achieved a significant reduction in MAC operations for all test cases. Increasing the number of convolution layers does not result in superior results but increasing the number of hidden channels improves performance. Regarding different convolution types, using group convolution only has the lowest complexity. It reduces the complexity significantly with a slight drop in performance. The MSGC layer can achieve similar or even better performances to the regular CNN layer with reduced complexities. In the best case (M = 5, c h = 9, n h = 64), 26.11% of MACs and 7.05% of parameters are removed while the accuracy and F1 score are improved by 0.50% and 1.75%. Applying MSGC with channel shuffling improves the ACC and F1 scores in almost all cases, with no increase in complexity. However, we note one exception for the network with M = 7, c h = 18 and n h = 128. In this case, MSGC without channel shuffling outperforms the complete MSGC, achieving an accuracy of 98.41% and an F1 score of 93.01%. We think this is because this configuration coincidentally suits the dataset. The performance improvement of MSGC is from the multi-scale operation that captures features on different scales. In practice, enabling or disabling channel shuffling, varying the number of Conv layers and hidden channels will make the model sensitive to different scales.

C. Results

85

IV. CONCLUSION

This paper proposes a novel optimization to the conventional convolution layer which can be an efficient replacement for convolutional layers in low complexity models. The experimental results show that the MSGC version of the network consistently outperformed the original version of the network. Furthermore, this improvement in complexity is achieved without sacrificing performance, and this model optimization method can co-exist with other approaches for further model compression. It provides a candidate component for building an AI-based system for low-power, resource-limited edge devices.

Fig. 3 :

 3 Fig. 3: Model Architecture used in our experiments

  MSGC No CS refers to Multi-Scale Group Convolution without Channel Shuffling.

	Channels	Conv Layers	Conv Type	MACs	Parameters	ACC	F1
			Regular Conv	111,757	6,411	97.51%	89.21%
		7	MSGC	78,853 (70.56%)	5,781 (90.17%)	97.65% (+0.14%)	89.94% (+0.73%)
			MSGC No CS *	78,853 (70.56%)	5,781 (90.17%)	97.03% (-0.49%)	87.03% (-2.19%)
			Group Conv	57,361 (51.33%)	5,439 (84.84%)	97.33% (-0.19%)	88.32% (-0.90%)
			Regular Conv	82,183	5,871	97.22%	88.31%
	9	5	MSGC	60,727 (73.89%)	5,457 (92.95%)	97.72% (+0.50%)	90.06% (+1.75%)
			MSGC No CS	60,727 (73.89%)	5,457 (92.95%)	97.63% (+0.41%)	89.77% (+1.46%)
			Group Conv	44,959 (54.71%)	5,223 (88.96%)	97.28% (+0.06%)	87.92% (-0.39%)
			Regular Conv	52,609	5,331	97.50%	89.20%
		3	MSGC	45,481 (86.45%)	5,133 (96.29%)	97.45% (-0.05%)	89.00% (-0.20%)
			MSGC No CS	45,481 (86.45%)	5,133 (96.29%)	97.45% (-0.05%)	88.87% (-0.32%)
			Group Conv	35,437 (67.36%)	5,007 (93.92%)	97.15% (-0.35%)	87.63% (-1.56%)
			Regular Conv	377,869	15,609	97.94%	91.14%
		7	MSGC	243,373 (64.41%)	13,053 (83.62%)	98.06% (+0.12%)	91.61% (+0.47%)
			MSGC No CS	243,373 (64.41%)	13,053 (83.62%)	98.41% (+0.47%)	93.01% (+1.87%)
			Group Conv	166,045 (43.94%)	11,721 (75.09%)	97.85% (-0.09%)	90.66% (-0.48%)
			Regular Conv	267,205	13,557	98.01%	91.45%
	18	5	MSGC	178,501 (66.80%)	11,865 (87.52%)	98.01% (-0.01%)	91.43% (-0.02%)
			MSGC No CS	178,501 (66.80%)	11,865 (87.52%)	97.95% (-0.07%)	91.19% (-0.25%)
			Group Conv	124,069 (46.43%)	10,965 (80.88%)	97.93% (-0.08%)	90.99% (-0.46%)
			Regular Conv	156,541	11,505	97.87%	90.99%
		3	MSGC	113,629 (72.59%)	10,677 (92.80%)	97.95% (+0.08%)	91.26% (+0.26%)
			MSGC No CS	113,629 (72.59%)	10,677 (92.80%)	97.80% (-0.07%)	90.59% (-0.41%)
			Group Conv	82,093 (52.44%)	10,209 (88.74%)	98.01% (+0.14%)	91.39% (+0.40%)

*

TABLE I :

 I Simulation Results of model variants using different convolutional layers.


Maryam Saeed
email: maryam.saeed@ucdconnect.ie

Olev Märtens
email: olev.martens@taltech.ee

Benoit Larras
email: ‡benoit.larras@junia.com

Antoine Frappé
email: antoine.frappe@junia.com

Barry Cardiff
email: barry.cardiff@ucd.ie

Event-Driven ECG Classification using Functional Approximation and Chebyshev Polynomials

Keywords: level-crossing ADC, event-driven, electrocardiograms, functional approximation, chebyshev polynomials, artificial neural networks, arrhythmia

Level-crossing ADCs reduce the size of data streams in wearable devices. However, in the context of electrocardiogram (ECG) signals, such an event-driven data source results in a variable length two-dimensional (time-amplitude tuples) data vector for each ECG beat. It is difficult to apply many standard signal processing techniques to this data making classifiers more complex. In this paper we resolve these difficulties by mapping the variable length 2D vectors to a fixed length feature vector comprising the first 80 coefficients of a Chebyshev polynomial expansion of the ECG beat. We show that, by using these 80 coefficients, the average percentage root-mean-square error is only ≈3.08%. Using this feature set we constructed a simple three-layered ANN binary (Normal / Abnormal) ECG classifier and we demonstrate 98.15% average accuracy and 96.07% average sensitivity. We also constructed a 4-class ANN, using the same ANN structure and we achieved 98.80% average accuracy and 91.5% average sensitivity. Both these networks have only 20k parameters and outperform the state-of-the-art classifiers, enabling low-power edge computing.

I. INTRODUCTION

Long-term monitoring of electrocardiogram (ECG) signals effectively identifies individuals with high risks of cardiovascular disease. Morphological changes in ECG beats can indicate abnormalities of the heart and the early detection of these changes is important for preventing life-threatening conditions. A recent study showed the significance of monitoring ECG signals for longer than 14 days, in both cardiovascular disease diagnoses and the study of the efficacy of treatment [START_REF] Mahajan | Efficacy of long-term monitoring on detecting critical cardiac arrhythmias[END_REF]. Low-power arrhythmia classifiers that prolong battery life in wearable devices [START_REF] Leong | A low-power vlsi arrhythmia classifier[END_REF]- [START_REF] Sivapalan | Annet: a lightweight neural network for ecg anomaly detection in iot edge sensors[END_REF] have been proposed for longterm ECG classification tasks. Many such methods use highdimensional ECG data with computationally expensive deep learning algorithms to classify arrhythmias. Some methods transmit ECG data off-chip to classify on a less powerconstrained devices e.g. a mobile phone [START_REF] Gradl | Real-time ecg monitoring and arrhythmia detection using android-based mobile devices[END_REF] or cloud-based systems [START_REF] Huda | A low-cost, low-energy wearable ecg system with cloud-based arrhythmia detection[END_REF].

Recent development in level-crossing analogue-to-digital converters (ADCs) has shown a promising future for lowpower wearable devices. Level-crossing ADCs can implicitly compress ECG data by 3x while maintaining signal quality in comparison with SAR ADCs [START_REF] Saeed | Evaluation of level-crossing adcs for event-driven ecg classification[END_REF], however the complexity of many non-uniform signal processing techniques make the realization of simple classifiers difficult on a low-powered device. For example, each beat in a event-driven ECG system comprises a variable number of (non-uniform) two dimension (2D) time-amplitude samples, making it impossible to feed into conventional deep learning networks. A recent study [START_REF] Melchert | Functional approximation for the classification of smooth time series[END_REF] explored functional approximation coefficients as features on several time-series datasets. They showed that there is only negligible or no change in classification accuracy using approximation coefficients as features. However, they achieved a large reduction in input dimensions. This study [START_REF] Melchert | Functional approximation for the classification of smooth time series[END_REF] is the motivation for our work. In [START_REF] Prakash | Detection of cardiac disease with less number of electrocardiogram sensor samples using chebyshev[END_REF], Chebyshev polynomials were used to classify entire ECG records in a dataset, however the focus of our work is the classification of individual ECG beats and so direct comparison with [START_REF] Prakash | Detection of cardiac disease with less number of electrocardiogram sensor samples using chebyshev[END_REF] is not possible.

In [START_REF] Moody | The impact of the mit-bih arrhythmia database[END_REF] the authors show that it is possible to represent continuous time functions with a weighted sum of Chebyshev polynomials. In our we apply this concept to transform the non-uniformly sampled ECG beats to a fixed number of features (namely the 80 Chebyshev weights/coefficients) and use these to classify arrhythmias within the MIT-BIH dataset. In this context, the major contributions of this paper are as follows: a) This work uses functional approximation to estimate a fixed number of features from non-uniform levelcrossing samples, and b) presents a low-power ANN (binary and 4-class) to classify four types of beats in the MIT-BIH dataset with results better than the state-of-the-art. The rest of the article is organized as follows. Section II presents the level-crossing ADC principles, introduces the functional approximation technique, and proposes a three-layered ANN classifier. Section III presents results and a comparison with the state-of-the-art. Finally, in Section IV, conclusions are drawn.

II. METHODOLOGY

A. Level-Crossing Sampling

In contrast to uniform sampling, level-crossing sampling does not sample at every clock cycle of the ADC. A levelcrossing ADC produces an output only when the amplitude of the input signal crosses either an upper or a lower threshold at which point the threshold are updated. If there are no level crossings within a clock roll-over period, the ADC will produce an output at that time. Fig. 1 shows the 7-bit levelcrossing ADC design used in this study with a 2385 Hz clock frequency and 6-bit clock timer (corresponding to a roll-over time of 26.8 ms). These parameters were chosen based on the analysis in [START_REF] Saeed | Evaluation of level-crossing adcs for event-driven ecg classification[END_REF].

At every level-crossing the output is a tuple comprising an ECG amplitude ECG out , and corresponding timing information T I n , where n is an index indicating the n th sample. T I n is an N bit word (in our case 6-bit) and indicates the time between two consecutive samples w.r.t the clock frequency of 2385 Hz. Fig. 2 shows an example of the level-crossing sampling of an ECG beat from the MIT-BIH dataset.

B. Functional Approximation Features

Functional approximation is a method for estimating an underlying unknown function from a set of observations. Functional approximation can also be used where samples exist at irregular intervals in time [START_REF] Melchert | Functional approximation for the classification of smooth time series[END_REF] assuming these samples are from a continuous signal [START_REF] Ramsay | Functional data analysis[END_REF].

We define the variable length vector ⃗ t j,i as the vector of time values (in seconds) where a variable number of level crossing events associated with the j th beat occurred. These time vectors are formed by only selecting those level crossing events that lie with a specific window about the R-peak of the j th beat, denoted here as occurring at time t R j , as follows:

-260 ms < t j,i -t R j < 400 ms

For the purpose of Chebyshev polynomial approximation we need to normalize these time values to the range [-1, +1] so we perform the following mapping:

⃗ x j,i = 2 ⃗ t j,i -t R j + 260 ms 260 ms + 400 ms -1
Let f j (x) denote the continuous function representing the j th beat. Then, the amplitudes of the level-crossings at the times t j,i can be represented as f j (x j,i ).

Given a suitable set of finite basis functions, f j (x) can be represented as a weighted sum of basis functions defined on the interval [-1, 1] as:

f j (x) = n k=0 c j,k g k (x)
where, c j,k are the functional approximation coefficients of the j th beat and they represent the amount of the respective basis function present in the signal f j (x). Motivated by [START_REF] Melchert | Functional approximation for the classification of smooth time series[END_REF], we choose Chebyshev polynomials of the first kind as the set of basis functions in this study. They can be represented by the following recurrence relations:

T 0 (x) = 1 T 1 (x) = x T n+1 (x) = 2xT n (x) -T n-1 (x)
Then, we compute the c j,k by a dot product of the points corresponding to the roots of the (n + 1) th polynomial [START_REF] Melchert | Functional approximation for the classification of smooth time series[END_REF]:

c j,k = 2 n + 1 n l=0 f j (x l )T k (x l ) (1)
where, x l time points are defined by:

where, x l = cos l + 1 2

π n + 1 (2)
However, we do not know exactly the points f j (x l ), therefore we use linear interpolation to approximate the f j (x l ) from the set of points f j (x j,i ) throughout this study.

1) Raised Cosine Windowing and Rotation: Fig. 3a shows the distribution of xl points on an ECG beat using [START_REF] Leong | A low-power vlsi arrhythmia classifier[END_REF]. It can be seen that the distribution of these points is mostly concentrated towards the edges of the ECG beat, whereas, the highest activity of the QRS complex lies around the R peak. We observed that using these x l time points resulted in high error when the ECG beat was reconstructed using the coefficients. To address this issue, we introduce a preprocessing step, which applies raised cosine windowing to the ECG beat followed by a rotation around the R peak.

The raised cosine window is used to taper of the edges of the ECG beat as follows: where, β is the roll-off factor and signifies the percentage of the signal being tapered off. Fig. 4 shows the effect of windowing the ECG beat with a roll-off factor of 0.3. Next, the beat is rotated around the R-peak such that the Rpeak is now concentrated around the edges as shown in Fig. 3b and the new time points x l capture the R peak more uniformly. Fig. 3c shows a reconstructed ECG beat using windowing and rotation method and 50 approximation coefficients. The percentage root-mean-square difference (PRD) between the original ECG beat and the reconstructed beat is 2.034%, which is well under the signal quality limit of 10% [START_REF] Saeed | Event-driven ecg classification using an open-source, lcadc based non-uniformly sampled dataset[END_REF].

w(t) =        1 |t| < (1 -β) 0 |t| > 1 1/2 1 + cos π |t|-(1-β) β elsewhere
2) Choosing the right number of coefficients: To understand the effect of windowing and rotation on the ECG beat, we ran simulations on the entire MIT-BIH dataset with and without the technique described in Section II-B1. We observed that without the windowing and rotation technique, 200 more x l time points and approximation coefficients were required for the PRD to converge. Similarly, we observed that when using windowing and rotation, the PRD converged to ≈3.08% at n = 200 interpolation time points, x l , and for k = 80 approximation coefficients. Therefore, for the rest of this article, we use n = 200 and k = 80 for estimating the approximation coefficients. Due to the page limits of this article, this comparative study is not presented in detail here.

C. Classification

For every beat in the MIT-BIH dataset, we generate a feature vector of 81 values using functional approximation with n = 200 and k = 80. This feature vector is fed into a three-layered ANN as shown in Fig. 5. There are 128, 64 and 32 neurons in the three hidden layers. Finally, each output is classified into one of the following four classes: normal (N), supraventricular (S), ventricular (V), and fusion (F) beats. The number of hidden layers and the respective number of neurons were chosen from several simulations of one to five hidden layers and 8-256 neurons per layer. Here, the dataset is split into 62%, 13%, and 25% of training, validation and test set. Synthetic minority over-sampling [START_REF] Brownlee | Smote for imbalanced classification with python[END_REF] is used to handle the class-imbalance. Adam optimizer, sparse categorical cross entropy loss and a batch size of 64 is used to train the model. We also tested the same classifier for binary classification performance, where each beat was classified into two classes: normal (N) and abnormal(A). Accuracy (ACC), sensitivity (SEN), positive predictivity (+PV or precision), and false positive rate (FPR) are used to assess the classifier performance as defined in [START_REF] Saeed | Evaluation of level-crossing adcs for event-driven ecg classification[END_REF].

III. RESULTS AND DISCUSSION

We also tested the classification performance using SVM with Gaussian and Cubic kernels as well as an ensemble of bagged trees. However, the ANN showed consistently superior performance for lower computational complexity. Table I shows the results for the four-class and two-class ANN. The 4class ANN has an average 98.80% accuracy, 91.5% sensitivity, 83.28% precision, and 1.81% false positive rate. The 2-class ANN has an average 98.15% accuracy, 96.07% sensitivity, 94.48% precision, and 3.88% false positive rate.

Table II shows a comparison with previous arrhythmia classifiers in literature using level-crossing ADCs. This study and [START_REF] Saeed | Evaluation of level-crossing adcs for event-driven ecg classification[END_REF] use a class-oriented evaluation scheme whereas, [START_REF] Zhao | A 13.34 µw event-driven patientspecific ann cardiac arrhythmia classifier for wearable ecg sensors[END_REF] uses a patient-specific scheme. The performance in [START_REF] Saeed | Evaluation of level-crossing adcs for event-driven ecg classification[END_REF] is comparable to this work but requires 6x more computational resources. The classifier in [START_REF] Zhao | A 13.34 µw event-driven patientspecific ann cardiac arrhythmia classifier for wearable ecg sensors[END_REF] outperforms this work for the S class, specifically in precision and requires much less computations. However, the classifier performs poorly for the minority F class with 49.03% sensitivity. Additionally, as it is a patient-specific design it requires labeled data from each patient to train the classifier before analyzing real-time data. This is generally not possible in real-time analysis. In [START_REF] Duforest | Slope-based event-driven feature extraction for cardiac arrhythmia classification[END_REF], a patient-specific 2-class ANN was presented using event-driven ECG using 296 parametes. They reported an average accuracy of 98.4% and an average sensitivity of 98.4% using selected records from the database.

Table III shows a comparison of this work with previous arrhythmia classifiers using uniform-sampling. In [START_REF] Janveja | A dnn-based low power ecg co-processor architecture to classify cardiac arrhythmia for wearable devices[END_REF], the authors used a patient-specific approach to classify S and V classes only. Their ANN model required a smaller amount of computational resources, however, the accuracy and sensitivity performance was much lower than this work. Similarly, in [START_REF] Tang | A real-time arrhythmia heartbeats classification algorithm using parallel delta modulations and rotated linear-kernel support vector machines[END_REF], a patient specific approach was used to classify S and V type beats using a very low power linear SVM that only requires 32 multiplications and 32 additions. However, the performance was much lower compared to this work.

A. Complexity considerations

As a measure of complexity of the pre-processing step, i.e. the conversion from the 2D variable length vectors to fixed length Chebyshev coefficients, we count the number of multiplications (or equivalently divisions) required. According the linear interpolation process to generate the 200 interpolated values, f j (x l ), requires approximately 400 multiplication operations. The subsequent windowing function requires an additional 200 multiplications and the rotation in time can be implemented for free by careful pre-rotation of the stored Chebyshev basis functions 1 . Each of the 80 coefficients is computed by a length 200 dot product each requiring 200 multiplications as per [START_REF] Mahajan | Efficacy of long-term monitoring on detecting critical cardiac arrhythmias[END_REF]. In total the pre-processing steps require 400+200+80*200 ≈ 17k multiplications per ECG beat. Similarly, we can take the number of multiplications in the ANN itself as being approximately equal to the number of trained parameters as enumerated in Tables II andIII.

IV. CONCLUSION

In this study, a low-power feature extraction technique using functional approximation and Chebyshev polynmials is presented. Using 200 Chebyshev time points and 80 approximation coefficients, event-driven ECG beats can be accurately reconstructed with only an average ≈3.08% PRD. A lowcomputation artificial neural network with three-hidden layers was also presented. The 2-class ANN classifies the dataset with 98.15% average accuracy, 96.07% average sensitivity. The 4-class ANN classifies the dataset with 98.80% average accuracy and 91.5% average sensitivity. Both these networks take only 20k parameters and outperforms the state-of-the-art classifiers. The combination of event driven ADCs and our presented algorithms represent a power efficient solution for ECG recording and classification at the edge.

Fig. 1 .

 1 Fig. 1. The 7-bit level-crossing design with a 2385Mz clock frequency and 6-bit clock timer.

Fig. 2 .

 2 Fig. 2. Level-crossing sampling of ECG from the MIT-BIH dataset.

Fig. 3 .Fig. 4 .

 34 Fig. 3. (a) An example of the x l time points using linear interpolation on an ECG beat for n=50, (b) after linear interpolation, windowing and a rotation for n=50, and (c) Reconstructed ECG beat using 50 approximation coefficients (PRD of 2.034%).

Fig. 5 .

 5 Fig. 5. The three-layered artificial neural network used to classify the ECG beats in the MIT-BIH dataset.

The windowing can be absorbed into the stored Chebyshev basis functions.

This work is supported by 1) JEDAI project under the Chist-Era Program; 2) Schlumberger Foundation's Faculty for the Future Program 3) Irish Research Council under the New Foundations Scheme and 4) Microelectronic Circuit Centre Ireland.


Alexandre Brouste

Christophe Dutang

Lilit Hovsepyan

Tom Rohmer
email: tom.rohmer@inrae.fr

One-step closed-form estimator for generalized linear model with categorical explanatory variables

Keywords: Generalized linear models, explicit estimators, categorical explanatory variables, one-step procedure

The parameters of generalized linear models are generally estimated by the maximum likelihood estimator (MLE), computed using a Newton-Raphson type algorithm that can be time-consuming for a large number of variables or modalities, or a large sample size. Explicit estimators exist for these models but they are not always asymptotically efficient, especially for simple effects models, although they are fast to calculate compared to the MLE. The article proposes a fast and asymptotically efficient estimation of the parameters of generalized linear models with categorical explanatory variables. It is based on a one-step procedure where a single step of the gradient descent is performed on the log-likelihood function initialized from the explicit estimators. This work presents the theoretical results obtained, the simulations carried out and an application to car insurance pricing.

Introduction

Generalized linear models (GLMs) are regression models where the distribution of the response variable belongs to the exponential family with a natural parameter which is a linear combination of explanatory variables (up to a link function, see e.g. [START_REF] Mccullagh | Generalized linear models[END_REF] for an introduction). The unknown parameters of such models are generally estimated by the maximum likelihood estimator (MLE) for which the asymptotic normality (and efficiency) have been established by [START_REF] Fahrmeir | Consistency and asymptotic normality of the maximum likelihood estimator in generalized linear models[END_REF]. Since the MLE has no closed-form formula in general, it is numerically computed by a Newton-type method (known as the Fisher scoring algorithm which can be rewritten as an iteratively re-weighted least square method (IWLS), see [START_REF] Mccullagh | Generalized linear models[END_REF]. This computation method is time-consuming for large datasets or for numerous explanatory variables.

The setting of sole categorical explanatory variables is singular due to the non-identifiability of the model and linear identifiability conditions are usually imposed. In this setting, for some specific models, closed-form MLE has been proposed by [START_REF] Brouste | Closed-form maximum likelihood estimator for generalized linear models in the case of categorical explanatory variables: application to insurance loss modeling[END_REF] which is asymptotically efficient. For other models, as the single effect models, closed-form alternative estimator which is consistent, asymptotically normal but not asymptotically efficient has also been built (see e.g. [START_REF] Brouste | A closed-form alternative estimator for GLM with categorical explanatory variables[END_REF].

Dealing with categorical explanatory variables is of particular interest in practical applications. A finite number of risk groups relying on categorical variables is used for the pricing of guarantees [START_REF] Denuit | Effective Statistical Learning Methods for Actuaries I: GLMs and Extensions[END_REF] or for modeling disease, fertility and milk production in dairy cattle in [START_REF] Kadarmideen | Linear and threshold model genetic parameters for disease, fertility and milk production in dairy cattle[END_REF]. The parameters of GLMs can be fastly estimated with the closed formula and large datasets and/or large number of modalities can be handled.

But in most situations, explanatory variables are used as single effect ((McCullagh & Nelder 1989, Chapters 4 and 6), [START_REF] Lindsey | Applying Generalized Linear Models[END_REF], (Denuit et al. 2020, Chapter 4), (Wuethrich & Merz 2021, Chapter 5)). Consequently, the loss of asymptotical efficiency for the gain of computation speed could be questioned. We propose in this paper a fast and asymptotically efficient procedure to estimate the parameters of GLM with categorical variables based on the one-step procedure.

The one-step procedure was initially considered for the estimation of parameters in independent and identically distributed (i.i.d.) samples (see Le [START_REF] Cam | On the asymptotic theory of estimation and testing hypothesis[END_REF]). In such procedure, an initial guess estimator is proposed which is fast to be computed but not asymptotically efficient. Then, a single step of the gradient descent method is done on the log-likelihood function in order to correct the initial estimation and reach asymptotic efficiency, see [START_REF] Brouste | Onestep -Le Cam's onestep estimation procedure[END_REF]. With some recent developments, the one-step procedure has been successfully generalized to more sophisticated statistical experiments as diffusion processes in [START_REF] Kamatani | Hybrid multi-step estimators for stochastic differential equations based on sampled data[END_REF], [START_REF] Gloter | Adaptive estimation for degenerate diffusion processes[END_REF], ergodic Markov chains in [START_REF] Kutoyants | On multi-step MLE-process for markov sequences[END_REF], inhomogeneous Poisson counting processes in [START_REF] Dabye | Method of moments estimators and multi-step mle for poisson processes[END_REF], fractional Gaussian and stable noises observed at high frequency in [START_REF] Brouste | Efficient estimation of stable Lévy process with symmetric jumps[END_REF], [START_REF] Brouste | Onestep estimation for the fractional gaussian noise model at high-frequency[END_REF].

The paper falls into 5 parts. The notations for GLMs are given in Section 2. The notion of restricted parameter is described in Section 3, which is new from our previous studies in [START_REF] Brouste | Closed-form maximum likelihood estimator for generalized linear models in the case of categorical explanatory variables: application to insurance loss modeling[END_REF], [START_REF] Brouste | A closed-form alternative estimator for GLM with categorical explanatory variables[END_REF], with the asymptotic properties of the MLE and the closed-form estimator. The fast one-step estimation procedure is also presented in Section 3 with a convergence result. Monte-Carlo simulations on samples of finite size are done in Section 4 and an application to car insurance is given in Section 5. Technical lemmas and the proof of the main are postponed in Appendix.

Preliminaries on GLMs

Notation for GLMs

Bold notations are given for the vectors. The index i ∈ I = {1, . . . , n} is reserved for the observations, while the indexes j, k are used for the explanatory variables.

The observation sample Y = (Y 1 , . . . , Y n ) is composed of independent random variables valued in Y ⊂ R where for i ∈ I, Y i belongs to the one-parameter exponential family of probability measures valued in Λ ⊂ R. In this setting, the log-likelihood log L of the sample is The parameters λ 1 , . . . , λ n in Equation ( 1) depend on a p-dimensional parameter ϑ ∈ Θ ⊂ R p . For i = 1, . . . , n, denote µ i = b (λ i (ϑ)). Theoretical moments of Y i are explicitly given as a function of a and derivatives of b

log L(ϑ, φ | Y ) = n i=1 λ i (ϑ)Y i -b (λ i (ϑ)) a(φ) + n i=1 c(Y i , φ), (1)
E ϑ Y i = b (λ i (ϑ)) = µ i and Var ϑ Y i = b (λ i (ϑ))a(φ) = V (µ i )a(φ).
(2)

The function

V : µ → V (µ) = b • (b ) -1 (µ)
is known as the variance function of the expectation µ. Using a twice continuously differentiable and bijective function g from b (Λ) to R, GLMs are defined by assuming the following relation between the expectation E ϑ Y i and the predictor g(µ i ) = x T i ϑ = η i , for all ϑ ∈ Θ, where η i are the linear predictors. The parameter ϑ ∈ Θ ⊂ R p and the parameter φ > 0 are unknown and are to be estimated. The function g is called the link function in the regression framework. In other words, the bijective function = (b ) -1 • g -1 is settled; then we have λ i (ϑ) = (η i ). We summarize with the following relations

X × Θ -→ D -1 Λ,
where D is the space of linear predictor and X the possible set of value of x i for i ∈ I. Here is chosen and, consecutively Θ, Λ and D must be set. We talk of canonical link function, when is the identity function, that is to say g = (b ) -1 .

Score and Fisher information for GLMs

We introduce T as the part of the log-likelihood depending on ϑ only in order to express the first-order condition verified by the maximum likelihood estimator (MLE)

T (ϑ | Y ) = n i=1 Y i (η i ) -b ( (η i )) .
Hence, the log-likelihood (1) can be rewritten as

log L(ϑ, φ | Y ) = T (ϑ | Y ) a(φ) + n i=1 c(Y i , φ).
The gradient of log L is

∇ (ϑ,φ) log L(ϑ, φ | Y ) =   ∇ ϑ T (ϑ | Y )/a(φ) -a (φ) a(φ) 2 T (ϑ | Y ) + n i=1 ∂ ∂φ c(Y i , φ)   =: U (ϑ)/a(φ) V(ϑ, φ) .
The MLE ( ϑ n , φ n ) for (ϑ, φ) satisfies

U ( ϑ n ) = 0 and V( ϑ n , φ n ) = 0. (3)
Using (2), we define the weight matrix and the explanatory variables matrix by

W (ϑ) =    1 (g (µ 1 )) 2 V (µ 1 ) . . . 1 (g (µn)) 2 V (µn)    , X =    x (1) 1 . . . x (p) 1 . . . x (1) n . . . x (p) n    , (4)
and Z the vector of mean deviations

Z(ϑ) =    (y 1 -µ 1 )g (µ 1 ) . . . (y n -µ n )g (µ n )    . (5)
Note that all diagonal terms W are positive so that W (ϑ) is invertible. U can be rewritten as a matrix multiplication U (ϑ) = X T W (ϑ)Z(ϑ), cf. (25) in Appendix A.2, where for a matrix M , M T denotes the transposition of the matrix M . We have the following expression for the score vector and the Fisher information

S n (ϑ) = X T W (ϑ)Z(ϑ) a(φ) , I n (ϑ) = X T W (ϑ)X a(φ) .
Since [START_REF] Fahrmeir | Consistency and asymptotic normality of the maximum likelihood estimator in generalized linear models[END_REF], under regularity conditions, the MLE ϑ n of ϑ asymptotically exists and as soon as the MLE is unique, that is to say there is no over-parametrization in the model, we have

I T /2 n (ϑ)( ϑ n -ϑ) L -→ n→+∞ N p (0 p , I p ) , with I 1/2 n I
T /2 n = I n where I p is the identity matrix of R p×p .

Estimation of GLMs with categorical variables

Consider the case where all m explanatory variables are categorical, that is for j = 1, . . . , m every observations (x

) i take values in a finite set {v j,1 , . . . , v j,d j } and x

(1) i = 1 is the intercept. Assuming values are unordered, x (j+1) i needs to be encoded using binary dummies as

x (j+1),k i = 1 {x (j+1) i =v j,k } , k ∈ {1, . . . , d j }.
These binary dummies can be used both in single-effect models or with cross-effect models.

In the following, we introduce the notion of restricted parameter which is new from our previous studies [START_REF] Brouste | Closed-form maximum likelihood estimator for generalized linear models in the case of categorical explanatory variables: application to insurance loss modeling[END_REF], [START_REF] Brouste | A closed-form alternative estimator for GLM with categorical explanatory variables[END_REF] and recall the asymptotic properties of the MLE of this parameter in Section 3.1. The properties of the closed-form estimator (CFE) are given in Section 3.2.

It is worth mentioning that, on the one hand, the MLE will be shown to be asymptotically efficient but time-consuming in the general setting. On the other hand, the CFE is fast to be computed but is not asymptotically efficient.

Notion of restricted parameter

For the sake of presentation, we firstly consider the single-effect model and secondly present the general case.

Model with single effects only

First let consider the GLMs with single effect only

g (E ϑ Y i ) =ϑ 1 + m+1 j=2 d j k=1 x (j),k i ϑ (j) k . (6)
For j = 2, . . . , m + 1, denote ϑ (j) = (ϑ

(j)
k ) k=1,...,d j . Hence the parameter vector is written as ϑ = (ϑ 1 , ϑ (2) , . . . , ϑ (m+1) ) T . The linear predictor vector η = (η i ) i=1...,n can be rewritten as η = Xϑ, with

X = 1 n x (2) . . . x (m+1) , x (j) =    x (j),1 1 . . . x (j),d j 1 . . . . . . x (j),1 n . . . x (j),d j n    .
Because the redundancies of the matrix X going from

d j k=1 x (j),k i
= 1 for all i, j, at least m linear conditions are needed to be imposed to get the identifiability of the model. Let contrast vectors R j such that R T j ϑ (j) = 0, j = 1, . . . , m. This is equivalent to Rϑ = 0 with

R =    0 R T 1 . . . 0 . . . . . . 0 0 . . . R T m    .
Note that the matrix R is a particular case of the constrast matrix considered in [START_REF] Brouste | A closed-form alternative estimator for GLM with categorical explanatory variables[END_REF]. In other words, the identifiability conditions assume that one component of ϑ (j) can be rewritten as a linear combination of the others. Therefore, we define a restricted parameter θ(j+1) of size d j -1 and a matrix B j of size d j × (d j -1) such that

ϑ (j+1) = B j θ(j+1) . (7)
In Table 1, two examples of such vectors R j are given as well as the corresponding restricted parameter and the corresponding coding matrix B j . Note that for a given R j , the choice of the restricted parameter is not necessarily unique. We refer to Venables (2023) for further examples of coding matrices B j and associated contrasts matrices in the R statistical software (R Core Team 2023).

Name

R j B j implication θ(j+1) R code zero-sum (1, . . . , 1) I d j -1 -1 T d j -1 ϑ (j+1) d j = - d j -1 k=1 ϑ (j+1) k    ϑ (j+1) 1 . . . ϑ (j+1) d j -1    contr.sum ref. category (1, 0, . . . , 0) 0 T d j -1 I d j -1 ϑ (j+1) 1 = 0    ϑ (j+1) 2 . . . ϑ (j+1) d j    contr.treatment
Table 1: Contrast examples and restricted parameters

Using the restricted parameters and the matrices B 1 , . . . , B m , the linear predictor can be rewritten as η = X θ, ( 8)

with X = 1 n x (2) B 1 . . . x (m+1) B m , θ = (ϑ 1 , θ(2) , . . . , θ(m+1) ).
Using Appendix A.3, the score vector and the Fisher information write in an analogous way as the previous section

S(ϑ) = XT W (ϑ)Z(ϑ) a(φ) , Ĩn (ϑ) = XT W (ϑ) X a(φ) . (9)
Note that the matrix W can be big for large datasets, and the Information matrix can be time-consuming. To reduce the computing time, note that, it can be rewritten as Ĩn (ϑ) = ( X S) T X/a(φ), where is Hadamard's product and S is the vector constituted with the diagonal elements of W .

General case

To take all possible GLM settings into account, we consider a GLM with predictor defined as

g (E ϑ Y i ) =ϑ 1 + m+1 j=2 d j k=1 x (j),k i ϑ (j) k
Intercept and single effect

+ j 2 <j 3 k 2 ,k 3 x (j 2 ),k 2 i x (j 3 ),k 3 i ϑ (j 2 ,j 3 ) k 2 ,k 3 Double effect + j 2 <j 3 <j 4 k 2 ,k 3 ,k 4 x (j 2 ),k 2 i x (j 3 ),k 3 i x (j 4 ),k 4 i ϑ (j 2 ,j 3 ,j 4 ) k 2 ,k 3 ,k 4 Triple effect + . . . . . . + k 2 ,...,k m+1 x (2),k 2 i . . . x (m+1),k m+1 i ϑ (2,...,m+1) k 2 ,...,k m+1 , All crossed effect ( 10
)
where g is the link function and indexes j i are in {2, . . . , m + 1} and k j are in {1, . . . , d j } for j = 2, . . . , m + 1. The linear predictor η = (η i ) i=1...,n can be rewritten as η = Xϑ, with m+1) x (2,3) . . . x (m,m+1) x (2,3,4) . . . x (m-1,m,m+1) . . . x (2,...,m+1) , and with for 2 ≤ j 2 < j 3 < . . . ,

X = 1 n x (2) . . . x (
x (j 2 ,j 3 ,...) =    x (j 2 ,j 3 ,...,),1 1 . . . x (j 2 ,j 3 ,...),d j 2 d j 2 ... 1 . . . . . . x (j 2 ,j 3 ,...),1 n . . . x (j 2 ,j 3 ,...),d j 2 d j 2 ... n    , x (j 2 ,j 3 ,...,),k 1 k 2 ,... i = x (j 2 ),k 1 i x (j 3 ),k 2 i . . .
The unknown parameter vector is ϑ = ϑ 1 , (ϑ (j) ) j , (ϑ (j 2 ,j 3 ) ) j 2 <j 3 , (ϑ (j 2 ,j 3 ,j 4 ) ) j 2 <j 3 <j 4 , . . . , (ϑ (2,...,m+1) ) T , with ϑ (j) = (ϑ

(j) k ) k , ϑ (j 2 ,j 3 ) = (ϑ (j 2 ,j 3 ) k 2 ,k 3 ) k 2 ,k 3 ,. . . , ϑ (2,...,m+1) = (ϑ (2,...,m+1) k 2 ,...,k m+1 ) k 2 ,...,k m+1
). The total number of parameters of the model in ( 10) is

p = 1 + m+1 j=2 d j + j 2 <j 3 d j 2 d j 3 + • • • + m+1 j=2 d j .
Because the model in ( 10) is not identifiable, we need to impose at least p -m+1 j=2 d j linear conditions. Here, we propose to reduce the dimension of the parameters. name parameter size of the restricted param. contrast number Intercept

ϑ 1 1 0 single effect ϑ (j) d j -1 d j -(d j -1) double effect ϑ (j 2 ,j 3 ) (d j 3 -1)(d j 2 -1) d j 3 d j 2 -(d j 3 -1)(d j 2 -1) triple effect ϑ (j 2 ,j 3 ,j 4 ) (d j 4 -1)(d j 3 -1)(d j 2 -1) d j 4 d j 3 d j 2 -(d j 4 -1)(d j 3 -1)(d j 2 -1) . . . . . . . . . all crossed-effect ϑ (2,...,m+1) m+1 j=2 (d j -1) m+1 j=2 d j - m+1 j=2 (d j -1)
Table 2: size of the restricted parameters and contrast numbers

In Table 2, the number of considered linear contrasts and the size of the restricted parameter θ(j) are given for each elements of the parameter vector. Note that this decomposition is not unique. For example in the single-effect part, one could impose a null intercept or a zero-condition on one of the parameter.

Let R j , R j 2 ,j 3 , . . . , R j 2 ,j 3 ,...,j m+1 , 2 ≤ j 2 < j 3 < . . . < j m+1 ≤ m + 1, matrices such that R T j ϑ (j) = 0, R T j 2 ,j 3 ϑ (j 2 ,j 3 ) , . . . , R T 2,...,m+1 ϑ (2,...,m+1) = 0. That is to say, Rϑ -1 = 0, with R the diagonal block matrix R = diag ((R j ) j , (R j 2 ,j 3 ) j 2 ,j 3 , . . . , R 2,...,m+1 ) , and ϑ -1 the parameter vector without the intercept ϑ 1 .

For all m-effect parts of the model, we can define a restricted parameter θ(j) , θ(j 2 ,j 3 ) , . . ., θ(2,...,m+1) of respective size dj = d j -1, dj 2 ,j 3 = (d j 3 -1)(d j 2 -1), . . ., d2,...,m+1 = m+1 j=2 (d j -1) and the corresponding coding matrix B j , B j 2 ,j 3 , . . ., B 2,...,m+1 respectively of size d j × dj , d j 2 d j 3 × dj 2 ,j 3 , . . . , m+1 j=2 d j × d2,...,m+1 such that ϑ (j) = B j θ(j) , ϑ (j 2 ,j 3 ) = B j 2 j 3 θ(j 2 ,j 3 ) , . . . , ϑ (2,...,m+1) = B 2,...,m+1 θ(2,...,m+1) .

The coding matrices can be rewritten in terms of the Kronecker product. Consider B 2 , B 3 and

B 4 of respective dimension d 2 × (d 2 -1), d 3 × (d 3 -1) and d 4 × (d 4 -1). Then the Kronecker product B 2,3 = B 3 ⊗ B 2 has dimension d 3 d 2 × (d 3 -1)(d 2 -1). In the same way, B 2,3,4 = B 4 ⊗ B 2,3 has dimension d 4 d 3 d 2 × (d 4 -1)(d 3 -1)(d 2 -1). More generally, for B j of dimension d j × (d j -1), 2 j=m+1 B j has dimension 2 j=m+1 d j × 2 j=m+1 (d j -1
). Using the restricted parameters and the coding matrices B j , B j 2 j 3 , . . . , B 2,...,m the linear predictor can be rewritten as

η = X θ, ( 12
)
with the new regressor matrix

X = 1 n x (2) B 1 . . . x (m+1) B m x (2,3) B 2,3 . . . x (m,m+1) B m,m+1 . . . x (2,...,m+1) B 2,...,m+1
Note that the linear predictor η i in ( 12) is identical whether we use ϑ or θ, leading to the same expectation µ i , henceforth the same matrices W (.) and Z(.).

Asymptotic properties of the MLE and the closed-form estimator

Since the covariates are supposed to be categorical in this paper, the vector of linear predictors η defined in (12) takes d = j d j distinct values namely h 1 , . . . , h d . We consider the unique d × p matrix Q composed of binary dummies such that η = Qϑ, η = (h j ) j=1,...,d .

For example, in the single-effect model ( 6

), Q = (M () m , M (2) m , . . . , M (m+1) m ) with M () m = 1 d m+1 ...d 2 ,
and for 2 < j < m + 1

M (2) m = 1 d m+1 ...d 3 ⊗ I d 2 , M (m+1) m = I d m+1 ⊗ 1 dm...d 2 , M (j) m = 1 d m+1 ⊗ I d j ⊗ 1 d j-1 ...d 2 .
The general form of Q for the model ( 10) is given in [START_REF] Brouste | A closed-form alternative estimator for GLM with categorical explanatory variables[END_REF].

Using the same restricted parameter as (11), and considering vectors R j of size d j+1 such that R j ϑ (j+1) = 0, the restricted linear predictors η can be rewritten as previously as

η = Q θ.
For example in the single-effect model

Q = M () m M (2) m B 1 . . . M (m+1) m B m , θ = (ϑ 1 , θ(2) , . . . , θ(m+1) ).
Note the the score S and the information matrix Ĩn (ϑ) can be rewritten in term of the Q matrix as following a(φ)

S(ϑ) = n QT Σ n (ϑ)Z (ϑ), a(φ) Ĩn (ϑ) = n QT Σ n (ϑ) Q,
with Σ n is the diagonal matrix whose the diagonal elements are

Σ n,j,j (ϑ) = m j n (g (µ j )) 2 V (µ j ) -1 , µ j = g -1 (h j ), j = 1, . . . , d, (13)
and

Z (ϑ) = ((ȳ (j) n -µ j )g (µ j )) j=1,...,d , (14)
see Appendix A.3. Hence, when p j satisfies m j n → p j as n → ∞, the asymptotic distribution of the MLE converges in distribution

√ n( θn -θ) L -→ n→+∞ N p 0 p , Ĩ-1 (ϑ) , (15)
with p the dimension of the restricted parameter and the information matrix expressed as

Ĩ(ϑ) = QΣ QT a(φ) , Σ(ϑ) = diag(v 1 , . . . , v d ), (16)
with diagonal elements are asymptotic variances v j = p j ((g (µ j )) 2 V (µ j )) -1 . Because the MLE can be time-consuming for large datasets or for a large number of explanatory variables or modalities for each variables, as soon as Q T Q + R T R is definite positive, an alternative closed-form estimator (CFE) which avoids the using of a time-consuming IWLS algorithm has been defined in [START_REF] Brouste | A closed-form alternative estimator for GLM with categorical explanatory variables[END_REF].

The closed-form estimator of the restricted parameter is

θCFE n = ( QT Q) -1 QT g(Y n ), where g(Y n ) =         g(Y 1 n ) . . . g(Y k n ) . . . g(Y d n )         , Y k n = n i=1;η i =h k Y i m k , m k = #{i ∈ {1, . . . , n}; η i = h k }. (17)
From Theorem 1 in [START_REF] Brouste | A closed-form alternative estimator for GLM with categorical explanatory variables[END_REF], we know the asymptotic distribution of θCFE

n : √ n( θCFE n -θ) L -→ n→+∞ N p 0 p , a (φ) ( QT Q) -1 QT Σ -1 ( θ) Q( QT Q) -1 , (18)
with Σ(ϑ) -1 the diagonal matrix whose diagonal elements are 1/v j and p j satisfies

m j n → p j as n → ∞.
As mentioned in [START_REF] Brouste | A closed-form alternative estimator for GLM with categorical explanatory variables[END_REF], as soon as the closed-form estimator coincides with the MLE (in this case the Q is a square matrix), θCFE n = Q-1 g( Ȳn ), and consequently using the delta method we get

√ n( θCFE n -θ) L -→ n→+∞ N p 0 p , a (φ) Q-1 Σ -1 (ϑ)( Q-1 ) T . (19)
That is the asymptotic variance of √ n( θCFE n -θ) coincides with Ĩ-1 (ϑ) defined in (15). In general, for instance in single effect models, the CFE is not the MLE and is not asymptotically efficient. For this reason, it is desirable to add a single step of the gradient descent method in order to reach asymptotic efficiency. Hence we consider a one-step version of θCFE n in the next subsection.

One-Step Closed-form Estimator

The One-Step Closed-form Estimator (OS-CFE) of θ is defined as

θOS-CFE n = θCFE n + Ĩn ( θCFE n ) -1 S( θCFE n ) = ( Q Q) -1 QT g(Y n ) + ( QT Σ n ( θCFE n ) Q) -1 QT Σ n ( θCFE n )Z ( θCFE n ),
where Z is defined in (14). It is worth emphasizing that the OS-CFE of θ does not depend on the dispersion parameter φ by simplification. The main result is that the OS-CFE of the restricted parameter θ is asymptotically equivalent in probability to the MLE. The proof of Theorem 1 is postponed in Appendix B.

Theorem 1. Under some regular conditions, as soon as for all j = 1, . . . , d the frequencies

m j n → p j as n → ∞, √ n( θOS-CFE n -θn ) P -→ n→+∞ 0.
It also means that the OS-CFE is asymptotically normal with an optimal asymptotic variance, i.e.

√ n( θOS-CFE

n -θ) L -→ n→+∞ N p 0 p , Ĩ-1 (ϑ) ( 20
)
where Ĩ(ϑ) is defined in ( 16).

Monte Carlo illustrations

The performances on finite size samples of the aforementioned estimators (MLE, CFE, OS-CFE), in terms of computation times and asymptotic variances, are assessed with numerical examples. Monte Carlo simulations of samples for Poisson and Gamma GLMs are done. All computations are carried out with the R statistical software (R Core Team 2023). More precisely, for the Poisson GLMs, the sample Y = (Y 1 , . . . , Y n ) is composed of independent Poisson-distributed random variables with respective distributions

f (y i , µ i ) = µ y i i y i ! exp(-µ i y i ), µ i > 0, y i ∈ N, i = 1, . . . , n.
In other words, Y i is characterized in Equation ( 1) by

λ i = µ i , a(φ) = 1, b(λ) = exp(λ) and c(y i , φ) = -log(y i !).
For the Gamma GLMs, the sample Y = (Y 1 , . . . , Y n ) is composed of independent Gammadistributed random variables with respective distributions

f (y i , (α, β i )) = y α-1 i Γ(α) β α i exp(-β i y i ), α > 0, β i > 0, y i > 0, i = 1, . . . , n.
In other words, Y i is characterized in Equation ( 1) by

λ i = - β i α , a(φ) = φ = 1 α , b(λ) = -log(-λ) and c(y i , φ) = 1 φ -1 log(y i ) -log Γ 1 φ + 1 φ log 1 φ .
In this numerical example the canonical setting is used ( is the identity function) leading to a log link function for the Poisson distribution (g(x) = log(x)) and the inverse link function for the Gamma distribution (g(x) = 1/x).

Simulations with a fixed sample size and a modality number

In our simulations, we consider two explanatory variables x

(2) i , x

(3) i and d 2 , d 3 modalities. As a starting point d 2 and d 3 have respectively been taken equal to 2 and 3, and the true parameters are chosen arbitrarily for each case as represented in Table 3. For the Gamma distribution, the dispersion parameter φ is also chosen arbitrary: φ = 8. We consider B = 10 4 Monte Carlo simulations of samples with the size of n = 10 4 .

Distribution Link function intercept Variable 2

Variable 3 g

(x) ϑ 1 ϑ (2) 1 ϑ (2) 2 ϑ (3) 1 ϑ (3) 2 ϑ (3) 3 Gamma 1/x 10 1 -1 2 3 -5 Poisson log(x) 0.05 1 -1 0.5 0.5 -1
Table 3: True parameter values Further, for the sake of ease, histograms for Poisson and Gamma GLMs are drawn for only ϑ 1 and ϑ

(2) 2 parameters. The sequence of OS-CFE naturaly overcomes the performance of CFE in terms of asymptotic variance (see Figures 1 and2). According to the other comparison in terms of computation time which is highlighted in Table 4, OS-CFE is almost 50 times faster than MLE to be computed for the dataset with the size of n = 10 4 .

Simulations with an increasing sample size

In this section for both Poisson and Gamma GLMs keeping the number of modalities constant (d 2 = 2 and d 3 = 3), the sample sizes are increasing from n = 10 3 to n = 10 7 . The number of Monte Carlo simulations are set to B = 100. The total computation time of the three estimators are computed. The comparative gain of the OS-CFE over the MLE is increasing with the sample size for both distributions, from 4 to 80 times faster as shown in Tables 5 and6.

Dataset

Simulations with an increasing number of modalities

Table 7 summarizes the results of the Monte Carlo simulations with an increasing number of modalities when the number of simulations are fixed to B = 100 and sample size is n = 10 5 . Simulations have been done for 2 categorical variables with equal number of modalities d 2 = d 3 = d varying from d = 5 to 40. In this particular case, the true parameter value for Gamma distribution are chosen by simply taking the ϑ 1 = 3d + 1, the ϑ In addition, for the Poisson distribution case, ϑ 1 has been chosen to be equal to 0.5, ϑ

1 to ϑ (j) d-1 are 1, 2, . . . , 1 or 2 each divided by the sum of all 1 and 2's, and the ϑ (j) d is the minus sum of all modalities for each of the two variables.

The data presented in Table 7 demonstrates that as the number of modalities increases, the computation time for all estimators also increases for both Poisson and Gamma GLMs. The OS-CFE method demonstrates faster computation times when compared to the MLE method, even at the highest number of modalities. For confidentiality reasons, the modality values are not revealed. The single effect models

g (E ϑ Y i ) = ϑ 1 + m+1 j=2 d j k=1 x (j),k i ϑ (j) k ,
is generally used by the insurers to model the claim amounts with (non-canonical) Gamma GLMs with a log link function (g(x) = log(x)). The "reference category" linear contrast has been used. It is worth recalling that in this setting the MLE has no closed-form and the closed-form estimator is not efficient. In order to compute the log-likelihood, we use Equation (3) to fit the dispersion parameter.

The one-step estimator has been applied to the Covea dataset:

Conclusion

Generalized linear models with single effects and sole categorical explanatory variables are widely used in different applications (e.g., insurance, agriculture, biology). On the one hand, the classical iteratively re-weighted least-square calibration algorithm is asymptotically efficient but can be time consuming for large datasets and/or large number of variables. On the other hand, closed-form estimators proposed in [START_REF] Brouste | A closed-form alternative estimator for GLM with categorical explanatory variables[END_REF] are fast to be computed but are not asymptotically efficient.

In this paper, we proposed a fast and asymptotically efficient method for the calibration of GLMs with categorical explanatory variables based on the one-step procedure that can also be applied to single effect models. It is 30 times faster than the classical method on both simulated and real datasets.

The proposed estimator could be further used for variable and/or model selection. For instance, the model selection with fastly computable Akaike Information Criterion could be treated in a further work. The one-step procedure can be extended in many directions: e.g., for multivariate regression models. R Core Team ( 2023 A Gradient and Hessian of the log-likelihood

Let b = (b ) -1 . We recall that X =    x (1) 1 . . . x (p) 1 . . . x (1) n . . . x (p) n    , Y =    Y 1 . . . Y n    .

A.1 Gradient and Hessian of a single observation

For ease of notation, we consider a single observation in (1). We recall that µ i = g -1 (η i ) and

V (µ i ) = b (λ i (ϑ)). The ith contribution is l i (ϑ) = λ i (ϑ)y i -b (λ i (ϑ)) a(φ) + c(y i , φ). (21)
To derive the first and second order derivatives, we first compute

∂λ i ∂ϑ j = ( b) (g -1 (η i )) × (g -1 ) (η i ) × x (j) i = 1 b ( b(g -1 (η i )) × x (j) i g (g -1 (η i )) = x (j) i g (µ i )V (µ i ) . Hence ∂l i (ϑ) ∂ϑ j = y i -b (λ i (ϑ)) a(φ) ∂λ i ∂ϑ j = x (j) i y i -µ i a(φ)g (µ i )V (µ i ) . ( 22
) Furthermore, ∂ 2 l i (ϑ) ∂ϑ j ∂ϑ l = ∂ ∂ϑ l y i -b (λ i ) a(φ i ) 1 b (λ i ) x (j) i g (µ i ) + ∂ ∂ϑ l 1 b (λ i ) y i -b (λ i ) a(φ) x (j) i g (µ i ) + ∂ ∂ϑ l x (j) i g (µ i ) y i -b (λ i ) a(φ) 1 b (λ i ) .
The first derivative term is

∂ ∂ϑ l y i -b (λ i ) a(φ) = -b (λ i ) a(φ) × ∂λ i ∂ϑ l = -b (λ i ) a(φ) × 1 b (λ i ) x (l) i g (µ i ) = -x (l) i a(φ)g (µ i ) .
The second derivative term is

∂ ∂ϑ l 1 b (λ i ) = - b (λ i ) (b (λ i )) 2 ∂λ i ∂ϑ l = - b (λ i ) (b (λ i )) 3 x (l) i g (µ i ) = - b (λ i ) (V (µ i )) 3 x (l) i g (µ i ) .
The third derivative term is

∂ ∂ϑ l x (j) i g (µ i ) = - x (j) i (g (µ i )) 2 ∂(g (µ i )) ∂ϑ l = - x (j) i g (µ i ) (g (µ i )) 2 ∂µ i ∂ϑ l = - x (j) i x (l) i g (µ i ) (g (µ i )) 3 , since ∂µ i ∂ϑ l = ∂(b (λ i )) ∂ϑ l = b (λ i ) ∂λ i ∂ϑ l = b (λ i ) × 1 b (λ i ) x (l) i g (µ i ) = x (l) i g (µ i ) .
This leads to

∂ 2 l i (ϑ) ∂ϑ j ∂ϑ l = - x (l) i x (j) i a(φ i )V (µ i )g (µ i ) 2 - y i -b (λ i ) a(φ i ) x (j) i x (l) i b (λ i ) g (µ i ) 2 (V (µ i )) 3 - y i -b (λ i ) a(φ i ) x (j) i x (l) i g (µ i ) (g (µ i )) 3 V (µ i ) . (23)

A.2 Score and Fisher information

The score component is obtained by summing ( 22) over observations

S j (ϑ) = ∂l i (ϑ) ∂ϑ j = 1 a(φ) n i=1 x (j) i y i -µ i g (µ i )V (µ i ) . (24)
Using the matrices of weights, covariables (4), and the mean deviation vector (5), the score is obtained by

S(ϑ) = 1 a(φ) X T W (ϑ)Z(ϑ). ( 25
)
Therefore, the information matrix is obtained by summing (23) over observations and taking the expectation with respect to Y i (using E(Y i ) -µ i = 0)

I n,l,j (ϑ) = -E i ∂ 2 l i (ϑ) ∂ϑ j ∂ϑ l = n i=1 x (l) i x (j) i a(φ)V (µ i )g (µ i ) 2 . ( 26
)
So that the Fisher information matrix is

I n (ϑ) = (I n,l,j (ϑ)) l,j = X T W (ϑ)X a(φ) . (27)
A.3 Information matrix and score in term of Q Define q(k) j the element of the jth row and kth column of Q. Note that the ith row of X is equal to the jth row of Q for all i such that η i = h j . Using (26), the information matrix is Ĩn (ϑ) = ( Ĩn,k,l (ϑ)) k,l=1,...,p with Ĩn,k,l (ϑ) = 1 a(φ)

n i=1 x(k) i x(l) i V (µ i )(g (µ i )) 2 = 1 a(φ) d j=1 n i=1;η i =h j x(k) i x(l) i V (µ i )(g (µ i )) 2 = 1 a(φ) d j=1 1 V (µ j )(g (µ j )) 2 n i=1;η i =h j x(k) i x(l) i = 1 a(φ) d j=1
m j V (µ j )(g (µ j )) 2 q(k) j q(l) j .

Hence the information matrix rewrites similarly to the non-restricted case Ĩn (ϑ) = n QT Σ n (ϑ) Q a(φ) . Now, using ( 22) the score component is

Ũk (ϑ) = n i=1 x(k) i (y i -µ i ) V (µ i )g (µ i ) = d j=1 1 V (µ j )g (µ j ) n i=1;η i =h j q(k) j (y i -µ j ) = d j=1 m j q(k) j (ȳ (j)
n -µ j ) V (µ j )g (µ j ) .

This yields to Ũ (ϑ) = n QT Σ n (ϑ)Z (ϑ).

B Proof of Theorem 1

Let us recall that the One-Step Closed-form Estimator (OS-CFE) is defined as θOS-CFE where I p is the p × p identity matrix.

In our setting, since the sequence of initial guess estimators is √ n-consistent (see Equation ( 19)), we get the asymptotic equivalence by showing that the quantity

  where a : R → R, b : Λ → R and c : Y × R → R are fixed real-valued measurable functions and φ is the dispersion parameter, e.g. McCullagh & Nelder (1989, Section 2.2).

Figure 1 :

 1 Figure 1: Histograms for the B = 10 4 simulations of the renormalized statistical errors of MLE, CFE, OS-CFE for the Poisson distribution with 2 categorical variables with d 2 = 2, d 3 = 3 for θ 1 = ϑ 1 and θ 2 = ϑ (2) 2 . Red and blue lines are the theoretical Gaussian asymptotic densities respectively of the MLE (in red) and CFE (in blue).

k

  are equal to k/d for k = 1, . . . , d -1, and the ϑ

  n (ϑ) = log L(ϑ | Y ). The mean-value theorem gives, for the initial sequence of guess estimators

Table 4 :

 4 Total computation time (s) based on B = 10 4 runs for Poisson and Gamma distributions

	Computation time	MLE	CFE OS-CFE
	Poisson	848.07 9.05	17.73
	Gamma	1601.44 10.65	31.61

Table 5 :

 5 Total computation time (s) based on B = 100 runs for Poisson distribution.

	size 10 3 10 4	10 5	10 6	10 7
	MLE	0.33 2.35 24.86 257.15 2678.38
	CFE	0.07 0.06 0.26	2.50	26.77
	OS-CFE	0.08 0.11 0.50	5.14	56.10
	Dataset size 10 3 10 4	10 5	10 6	10 7
	MLE	0.41 3.64 30.63 389.88 5120.97
	CFE	0.06 0.06 0.29	2.96	38.29
	OS-CFE	0.07 0.11 0.52	5.30	64.98

Table 6 :

 6 Total computation time (s) based on B = 100 runs for Gamma distribution.

Table 7 :

 7 Total computation times based on B = 100 runs for Poisson (canonical link) and Gamma (canonical link) GLMs.5 Application to claim amounts in car insuranceThe Covea Affinity dataset under study is composed of 76,446 claim amounts ranging from 4 to 33,531 EUR. Three covariates have been selected from the 124 available for the pricing of the guarantee

	d	Poisson		Gamma
		MLE CFE OS-CFE MLE CFE OS-CFE
	5 36.49 0.31	0.71	50.59 0.41	0.67
	10 71.95 0.58	0.92	109.39 0.60	0.91
	15 135.12 0.97	1.56	175.90 0.79	1.67
	20 150.19 1.15	2.42	181.79 0.92	1.92
	25 173.01 1.42	3.97	346.64 2.23	5.50
	30 264.51 2.61	10.27	429.85 2.51	9.58
	35 254.86 3.23	14.10	618.42 3.97	19.96
	40 343.16 3.85	25.74	675.58 4.60	28.12
	vehicle brand with d 2 = 2 modalities,
	pricing segment with d 3 = 6 modalities,
	age class with d 4 = 8 modalities.

Table 8 :

 8 Table 8 give parameter estimates for MLE, CFE and OS-CFE. The CFE and OS-CFE were almost 30 times faster to obtain than the MLE, with similar estimate and similar fitted log-likelihood. Values of ϑ n , log-likelihood and total computation time (s) for CFE, OS-CFE and MLE

		CFE OS-CFE	MLE
	ϑ 1	6.23	6.04	6.03
	ϑ (2) 2	0.24	0.08	0.03
	ϑ (3) 2	0.18	0.22	0.22
	ϑ (3) 3	-0.48	0.04	-0.01
	ϑ (3) 4	-0.07	0.08	0.09
	ϑ (3) 5	0.06	0.18	0.19
	ϑ (3) 6	0.20	0.21	0.22
	ϑ (4) 2	-0.07	0.00	0.01
	ϑ (4) 3	0.06	0.16	0.16
	ϑ (4) 4	0.17	0.18	0.18
	ϑ (4) 5	0.34	0.41	0.40
	ϑ (4) 6	0.11	0.44	0.42
	ϑ (4) 7	0.16	0.25	0.26
	ϑ (4) 8	-0.01	0.34	0.33
	log L -554,868 -553,708 -553,685
	Time (s)	0.01	0.01	0.30

  ), R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria. URL: https://www.R-project.org/

	Venables, B. (2023), codingMatrices: Alternative Factor Coding Matrices for Linear Model For-
	mulae. R package version 0.4.0.
	URL: https://CRAN.R-project.org/package=codingMatrices
	Wuethrich, M. & Merz, M. (2021), Statistical foundations of actuarial learning and its applications.
	SSRN papers.
	URL: https://ssrn.com/abstract=3822407

Acknowledgments We would like to thank Héloïse Bertrand and Nicolas Villette for fruitful discussions on the topic. This research benefited from the support of the ANR project 'Efficient inference for large and high-frequency data' (ANR-21-CE40-0021), the 'Chair Risques Émergents ou Atypiques en Assurance', under the aegis of Fondation du Risque, a joint initiative by Le Mans Université and MMA company, member of Covea group and the 'Chair Impact de la Transition Climatique en Assurance', under the aegis of Fondation du Risque, a joint initiative by Le Mans Université and Groupama Centre-Manche company, member of Groupama group. This preprint

has not undergone peer review or any post-submission improvements or corrections. The Version of Record of this article is published in Statistics and Computing, and is available online at https://10.1007/s11222-023-10313-4.

converges to zero as n → ∞.

For the quantity (B), we get from ( 15) and ( 18)

is bounded in probability.

For the quantity (A), by Markov's law of large number and Equation ( 23) we have

Using the consistency of the initial guess estimator and the uniform continuity of the Fisher information matrix ( 16), we get the convergence to zero of the quantity (A).


Alexandre Brouste

Christophe Dutang

Lilit Hovsepyan

Tom Rohmer

One-step closed-form estimator for generalized linear model with categorical explanatory variables

Keywords: modèles linéaires généralisés, procédure one-step, assurance generalized linear models, one-step procedure, insurance

L'article propose une procédure d'estimation rapide et asymptotiquement efficace des paramètres des modèles linéaires généralisés à variables catégorielles. Des estimateurs explicites existent [2] pour ces modèles mais ils ne sont pas toujours asympotiquement efficaces, notamment pour les modèles à effets simples. La procédure proposée dans cet article est basée sur une approche one-step où une unique étape de la descente de gradient est effectuée sur la fonction de log-vraisemblance initialisée à partir des estimateurs explicites. Ce travail présente de manière succincte les résultats théoriques obtenus, les simulations effectuées et une application à la tarification en assurance automobile.

Notations for GLMs with categorical variables

The observation sample Y = (Y 1 , . . . , Y n ) is composed of independent random variables valued in Y ⊂ R where for i ∈ I, Y i belongs to the one-parameter exponential family of probability measures valued in Λ ⊂ R. In this setting, the log-likelihood log L of the sample is

log L(ϑ, ϕ | Y ) = n i=1 λ i (ϑ)Y i -b (λ i (ϑ)) a(ϕ) + n i=1 c(Y i , ϕ), (1)
where a : R → R, b : Λ → R and c : Y × R → R are fixed real-valued measurable functions and ϕ is the dispersion parameter. Here, λ 1 , λ 2 , . . . , λ n are related to exogenous explanatory variables x i , i = 1, . . . , n by the relations

g(E ϑ Y i ) = x T i ϑ =: η i , for all ϑ ∈ Θ, i = 1, . . . , n, (2)
or, equivalently, b ′ (λ i ) = g -1 (η i ), i = 1, . . . , n.

The parameter ϑ ∈ Θ ⊂ R p and ϕ in (1) are to be estimated. The function g is called the link function in the regression framework.

Let us consider m categorical explanatory variables (z

(j)
i ) i=1,...,n , j = 2, . . . , m + 1, that take values in a finite set {v j,1 , . . . , v j,d j }. Assuming values are unordered, we encode our explanatory variables using binary dummies as

x (j),k i = 1 {z (j) i =v j,k } , k ∈ {1, .
. . , d j }. These binary dummies can be used both in single-effect models or in cross-effect models. The equation ( 2) can be generally reconstructed in the following way:

g (E ϑ Y i ) =ϑ 1 + m+1 j=2 d j k=1 x (j),k i ϑ (j) k
Intercept and single effect +

j 2 <j 3 k 2 ,k 3 x (j 2 ),k 2 i x (j 3 ),k 3 i ϑ (j 2 ,j 3 ) k 2 ,k 3 Double effect + j 2 <j 3 <j 4 k 2 ,k 3 ,k 4 x (j 2 ),k 2 i x (j 3 ),k 3 i x (j 4 ),k 4 i ϑ (j 2 ,j 3 ,j 4 ) k 2 ,k 3 ,k 4 Triple effect + . . . + k 2 ,...,k m+1 x (2),k 2 i . . . x (m+1),k m+1 i ϑ (2,...,m+1) k 2 ,...,k m+1 , All crossed effect (3)
where indexes j i are in {2, . . . , m + 1} and k j are in {1, . . . , d j } for j = 2, . . . , m + 1.

In the setting of categorical explanatory variables, the model (3) is not identifiable and linear constraints have to be imposed. For this reason, we consider a restricted parameter θ which is new from our previous studies [START_REF] Brouste | Closed-form maximum likelihood estimator for generalized linear models in the case of categorical explanatory variables: application to insurance loss modeling[END_REF][START_REF] Brouste | A closed-form alternative estimator for glm with categorical explanatory variables[END_REF].

The MLE ( θn , ϕ n ) for ( θ, ϕ) satisfies ( θn , ϕ n ) = arg max ( θ,ϕ)∈ Θ×R + * log L( θ, ϕ | Y ). ( 4
)
It is worth mentioning that the MLE is asymptotically efficient generally but not in a closed form and its computation can be time-consuming for large samples or numerous explanatory variables.

2 One-step closed-form estimator

Closed-form Estimator

The closed-form estimator (CFE) of the restricted parameter is defined in [START_REF] Brouste | Closed-form maximum likelihood estimator for generalized linear models in the case of categorical explanatory variables: application to insurance loss modeling[END_REF][START_REF] Brouste | A closed-form alternative estimator for glm with categorical explanatory variables[END_REF] by

θCFE n = ( QT Q) -1 QT g(Y n ).
where

g(Y n ) =         g(Y 1 n ) . . . g(Y k n ) . . . g(Y d n )         , Y k n = n i=1;η i =h k Y i m k , m k = #{i ∈ {1, . . . , n}; η i = h k } (5)
where h k are the different possible values of the linear predictors defined in [START_REF] Brouste | A closed-form alternative estimator for glm with categorical explanatory variables[END_REF]. Here, the matrix Q is related to the identifiability condition and to the structure of the model (3).

It is known that the CFE is asymptotically normal, namely

√ n( θCFE n -θ) L -→ n→+∞ N p ⋆ 0 p ⋆ , a (ϕ) ( QT Q) -1 QT Σ -1 ( θ) Q( QT Q) -1 (6)
where Σ -1 (ϑ) is defined in [START_REF] Brouste | A closed-form alternative estimator for glm with categorical explanatory variables[END_REF].

In general, for instance in single effect models, the CFE is not the MLE and is not asymptotically efficient. Hence, we consider a one-step version of θCFE n in the next subsection.

One-Step Closed-form Estimator

The One-Step Closed-form Estimator (OS-CFE) of θ is defined as

θOS-CFE n = θCFE n + Ĩn ( θCFE n ) -1 S( θCFE n )
where S stands for the score S( θ) = ∇ θL(( θ, ϕ) | Y ) and Ĩn stands for the Fisher Information matrix Ĩn ( θ) = -E ∇ 2 θL( θ, ϕ | Y ) which can be described in terms of the matrix Q.

It is worth emphasizing that the OS-CFE of θ does not depend on the dispersion parameter ϕ by simplification. The main result is that the OS-CFE of the restricted parameter θ is asymptotically equivalent to the MLE. Theorem 1. Under the regularity conditions, as soon as for all j = 1, . . . , d the frequencies

m j n → p j as n → ∞, √ n( θOS-CFE n -θn ) P -→ n→+∞ 0.
It also means that the OS-CFE is asymptotically normal with an optimal asymptotic variance.

Monte Carlo illustrations

The performances on finite size samples of the aforementioned estimators (MLE, CFE, OS-CFE), in terms of computation times and asymptotic variances, are assessed with numerical examples. Monte Carlo simulations of samples for Poisson and Gamma GLMs are done.

In this numerical example the canonical setting is used (ℓ is the identity function) leading to a log link function for the Poisson distribution (g(x) = log(x)) and the inverse link function for the Gamma distribution (g(x) = 1/x).

The sequence of OS-CFE naturally overcomes the performance of CFE in terms of asymptotic variance (see Figures 1 and2). According to the other comparison in terms of computation time which is highlighted in Table

Application to claim amounts in car insurance

The Covea Affinity dataset under study is composed of 76,446 claim amounts ranging from 4 to 33,531 EUR. Three covariates have been selected from the 124 available for the pricing of the guarantee

• vehicle brand with d 2 = 2 modalities,
• pricing segment with d 3 = 6 modalities,

• age class with d 4 = 8 modalities.

For confidentiality reasons, the modality values are not revealed.

The single effect models

g (E ϑ Y i ) = ϑ 1 + m+1 j=2 d j k=1 x (j),k i ϑ (j) k ,
is generally used by the insurers to model the claim amounts with (non-canonical) Gamma GLMs with a log link function (g(x) = log(x)). The "reference category" linear contrast has been used. It is worth recalling that in this setting the MLE has no closed form and the closed-form estimator is not efficient. In order to compute the log-likelihood, we use Equation (4) to fit the dispersion parameter.

Conclusion

Generalized linear models with single effects and single categorical explanatory variables are widely used in different applications (insurance, agriculture, biology). The classical iterative re-weighted least square calibration method is asymptotically efficient, but can be time consuming for large datasets. On the other hand, closed-form estimators proposed in [START_REF] Brouste | A closed-form alternative estimator for glm with categorical explanatory variables[END_REF] are faster to compute, but they are not asymptotically efficient.

In this paper we proposed fast and asymptotically efficient method for the calibration of GLMs with categorical explanatory variables based on the one-step procedure that can also be applied to single effect models. It is 30 times faster on the simulated and the Covea Affinity datasets compared with the classical methods. cient inference for large and high-frequency data' (ANR-21-CE40-0021), the 'Chair Risques Émergents ou Atypiques en Assurance', under the aegis of Fondation du Risque, a joint initiative by Le Mans Université and MMA company, member of Covea group and the 'Chair Impact de la Transition Climatique en Assurance', under the aegis of Fondation du Risque, a joint initiative by Le Mans Université and Groupama Centre-Manche company, member of Groupama group.

Figure 1 : 2 ) 2 .

 122 Figure 1: Histograms for the B = 10 4 simulations of the renormalized statistical errors of MLE, CFE, OS-CFE for the Poisson distribution with 2 categorical variables with d 2 = 2, d 3 = 3 for ϑ 1 and ϑ (2) 2 . Red and blue lines are the theoretical Gaussian asymptotic densities respectively of the MLE (in red) and CFE (in blue).

Figure 2 :

 2 Figure 2: Histograms for the B = 10 4 Monte Carlo simulations of the renormalized statistical errors of MLE, CFE, OS-CFE for the Gamma distribution (canonical link) with 2 categorical variables with d 2 = 2, d 3 = 3 for ϑ 1 and ϑ (2) 2 and fixed ϕ = 8. Red and blue lines are the theoretical Gaussian asymptotic densities respectively of the MLE (in red) and CFE (in blue).

Table 1 :

 1 1, OS-CFE is almost 50 times faster than MLE to be computed for the dataset with the size of n = 10 4 .

	Computation time	MLE	CFE OS-CFE
	Poisson	848.07 9.05	17.73
	Gamma	1601.44 10.65	31.61

Total computation time (s) for Poisson and Gamma distributions

Table 2 :

 2 The one-step estimator has been applied to the Covea dataset. It is almost 30 times faster than the MLE with similar estimation. Values of ϑ n , log-likelihood and total computation time (s) for CFE, OSCFE and MLE

		CFE OSCFE	MLE
	ϑ 1	6.23	6.04	6.03
	ϑ (2) 2	0.24	0.08	0.03
	ϑ (3) 2	0.18	0.22	0.22
	ϑ (3) 3	-0.48	0.04	-0.01
	ϑ (3) 4	-0.07	0.08	0.09
	ϑ (3) 5	0.06	0.18	0.19
	ϑ (3) 6	0.20	0.21	0.22
	ϑ (4) 2	-0.07	0.00	0.01
	ϑ (4) 3	0.06	0.16	0.16
	ϑ (4) 4	0.17	0.18	0.18
	ϑ (4) 5	0.34	0.41	0.40
	ϑ (4) 6	0.11	0.44	0.42
	ϑ (4) 7	0.16	0.25	0.26
	ϑ (4) 8	-0.01	0.34	0.33
	log L -554,868 -553,708 -553,685
	Time (s)	0.01	0.01	0.30

Acknowledgments We would like to thank Héloïse Bertrand and Nicolas Villette for fruitful discussions on the topic. This research benefited from the support of the ANR project 'Effi-