Group method of data handling (GMDH) is a family of inductive algorithms for computer-based mathematical modeling of multi-parametric datasets that features fully automatic structural and parametric optimization of models.

GMDH is used in such fields as data mining, knowledge discovery, prediction, complex systems modeling, optimization and pattern recognition.^[1] GMDH algorithms are characterized by inductive procedure that performs sorting-out of gradually complicated polynomial models and selecting the best solution by means of the external criterion. The last section of ^[2] contains a summary of the applications of GMDH in the 1970s.

Other names include "polynomial feedforward neural network",^[3] or "self-organization of models". It was one of the first deep learning methods, used to train an eight-layer neural net in 1971.^[4][5]

This section is based on.^[2]

This is the general problem of statistical modelling of data: Consider a dataset ${\displaystyle \{(x_{1},...,x_{k};y)_{s}\}_{s=1:n}}$, with ${\displaystyle n}$ points. Each point contains ${\displaystyle x_{1},...,x_{k}}$ observations, and one target ${\displaystyle y}$ to predict. How to best predict the target based on the observations?

First, we split the full dataset into two parts: a training set and a validation set. The training set would be used to fit more and more model parameters, and the validation set would be used to decide which parameters to include, and when to stop fitting completely.

The GMDH starts by considering degree-2 polynomial in 2 variables. Suppose we want to predict the target using just the ${\displaystyle i,j}$ parts of the observation, and using only degree-2 polynomials, then the most we can do is this:$${\displaystyle y\approx f_{a,b,c,d,e,h}(x_{i},x_{j}):=a+bx_{i}+cx_{j}+dx_{i}^{2}+ex_{j}^{2}+fx_{i}x_{j}}$$where the parameters ${\displaystyle a,b,c,d,e,f}$ are computed by linear regression. Now, the parameters ${\displaystyle a,b,c,d,e,f}$ depend on which ${\displaystyle i,j}$ we have chosen, and we do not know which ${\displaystyle i,j}$ we should choose, so we choose all of them. That is, we perform all ${\displaystyle {\frac {1}{2}}k(k-1)}$ such polynomial regressions:$${\displaystyle y\approx f_{(i,j);a,b,c,d,e,h}(x_{i},x_{j}):=a_{i,j}+b_{i,j}x_{i}+c_{i,j}x_{j}+d_{i,j}x_{i}^{2}+e_{i,j}x_{j}^{2}+f_{i,j}x_{i}x_{j}\quad \forall 1\leq i<j\leq k}$$obtaining ${\displaystyle {\frac {1}{2}}k(k-1)}$ polynomial models of the dataset.

We do not want to accept all the polynomial models, since it would contain too many models. To only select the best subset of these models, we run each model ${\displaystyle f_{(i,j);a,b,c,d,e,h}}$ on the validation dataset, and select the models whose mean-square-error is below a threshold. We also write down the smallest mean-square-error achieved as ${\displaystyle minMSE_{1}}$.

Suppose that after this process, we have obtained a set of ${\displaystyle k_{1}}$ models. We now run the models on the training dataset, to obtain a sequence of transformed observations: ${\displaystyle z_{1},z_{2},...,z_{k_{1}}}$. The same algorithm can now be run again.

The algorithm continues, giving us ${\displaystyle minMSE_{1},minMSE_{2},...}$. As long as each ${\displaystyle minMSE}$ is smaller than the previous one, the process continues, giving us increasingly deep models. As soon as some ${\displaystyle minMSE_{L+1}>minMSE_{L}}$, the algorithm terminates. The last layer fitted (layer ${\displaystyle L+1}$) is discarded, as it has overfit the training set. The previous layers are outputted.

More sophisticated methods for deciding when to terminate are possible. For example, one might keep running the algorithm for several more steps, in the hope of passing a temporary rise in ${\displaystyle minMSE}$.

Instead of a degree-2 polynomial in 2 variables, each unit may use higher-degree polynomials in more variables:^[1]

${\displaystyle Y(x_{1},\dots ,x_{n})=a_{0}+\sum \limits _{i=1}^{n}{a_{i}}x_{i}+\sum \limits _{i=1}^{n}{\sum \limits _{j=i}^{n}{a_{ij}}}x_{i}x_{j}+\sum \limits _{i=1}^{n}{\sum \limits _{j=i}^{n}{\sum \limits _{k=j}^{n}{a_{ijk}}}}x_{i}x_{j}x_{k}+\cdots }$

And more generally, a GMDH model with multiple inputs and one output is a subset of components of the base function (1):

${\displaystyle Y(x_{1},\dots ,x_{n})=a_{0}+\sum \limits _{i=1}^{m}a_{i}f_{i}}$

where f_i are elementary functions dependent on different sets of inputs, a_i are coefficients and m is the number of the base function components.

External criteria are optimization objectives for the model, such as minimizing mean-squared error on the validation set, as given above. The most common criteria are:

• Criterion of Regularity (CR) – least mean squares on a validation set.
• Least squares on a cross-validation set.
• Criterion of Minimum bias or Consistency – squared difference between the estimated outputs (or coefficients vectors) of two models fit on the A and B set, divided by squared predictions on the B set.^[1]

Like linear regression, which fits a linear equation over data, GMDH fits arbitrarily high orders of polynomial equations over data.^[6][7]

To choose between models, two or more subsets of a data sample are used, similar to the train-validation-test split.

GMDH combined ideas from:^[8] black box modeling, successive genetic selection of pairwise features,^[9] the Gabor's principle of "freedom of decisions choice",^[10] and the Beer's principle of external additions.^[11]

Inspired by an analogy between constructing a model out of noisy data, and sending messages through a noisy channel,^[12] they proposed "noise-immune modelling":^[6] the higher the noise, the less parameters must the optimal model have, since the noisy channel does not allow more bits to be sent through.

The model is structured as a feedforward neural network, but without restrictions on the depth, they had a procedure for automatic models structure generation, which imitates the process of biological selection with pairwise genetic features.

The method was originated in 1968 by Prof. Alexey G. Ivakhnenko in the Institute of Cybernetics in Kyiv.

Period 1968–1971 is characterized by application of only regularity criterion for solving of the problems of identification, pattern recognition and short-term forecasting. As reference functions polynomials, logical nets, fuzzy Zadeh sets and Bayes probability formulas were used. Authors were stimulated by very high accuracy of forecasting with the new approach. Noise immunity was not investigated.

Period 1972–1975. The problem of modeling of noised data and incomplete information basis was solved. Multicriteria selection and utilization of additional priory information for noiseimmunity increasing were proposed. Best experiments showed that with extended definition of the optimal model by additional criterion noise level can be ten times more than signal. Then it was improved using Shannon's Theorem of General Communication theory.

Period 1976–1979. The convergence of multilayered GMDH algorithms was investigated. It was shown that some multilayered algorithms have "multilayerness error" – analogous to static error of control systems. In 1977 a solution of objective systems analysis problems by multilayered GMDH algorithms was proposed. It turned out that sorting-out by criteria ensemble finds the only optimal system of equations and therefore to show complex object elements, their main input and output variables.

Period 1980–1988. Many important theoretical results were received. It became clear that full physical models cannot be used for long-term forecasting. It was proved, that non-physical models of GMDH are more accurate for approximation and forecast than physical models of regression analysis. Two-level algorithms which use two different time scales for modeling were developed.

Since 1989 the new algorithms (AC, OCC, PF) for non-parametric modeling of fuzzy objects and SLP for expert systems were developed and investigated.^[13] Present stage of GMDH development can be described as blossom out of deep learning neuronets and parallel inductive algorithms for multiprocessor computers. Such procedure is currently used in deep learning networks.^[14]

There are many different ways to choose an order for partial models consideration. The very first consideration order used in GMDH and originally called multilayered inductive procedure is the most popular one. It is a sorting-out of gradually complicated models generated from base function. The best model is indicated by the minimum of the external criterion characteristic. Multilayered procedure is equivalent to the Artificial Neural Network with polynomial activation function of neurons. Therefore, the algorithm with such an approach usually referred as GMDH-type Neural Network or Polynomial Neural Network. Li showed that GMDH-type neural network performed better than the classical forecasting algorithms such as Single Exponential Smooth, Double Exponential Smooth, ARIMA and back-propagation neural network.^[15]

Another important approach to partial models consideration that becomes more and more popular is a combinatorial search that is either limited or full. This approach has some advantages against Polynomial Neural Networks, but requires considerable computational power and thus is not effective for objects with a large number of inputs. An important achievement of Combinatorial GMDH is that it fully outperforms linear regression approach if noise level in the input data is greater than zero. It guarantees that the most optimal model will be founded during exhaustive sorting.

Basic Combinatorial algorithm makes the following steps:

• Divides data sample at least into two samples A and B.
• Generates subsamples from A according to partial models with steadily increasing complexity.
• Estimates coefficients of partial models at each layer of models complexity.
• Calculates value of external criterion for models on sample B.
• Chooses the best model (set of models) indicated by minimal value of the criterion.
• For the selected model of optimal complexity recalculate coefficients on a whole data sample.

In contrast to GMDH-type neural networks Combinatorial algorithm usually does not stop at the certain level of complexity because a point of increase of criterion value can be simply a local minimum, see Fig.1.

• Combinatorial (COMBI)
• Multilayered Iterative (MIA)
• GN
• Objective System Analysis (OSA)
• Harmonical
• Two-level (ARIMAD)
• Multiplicative–Additive (MAA)
• Objective Computer Clusterization (OCC);
• Pointing Finger (PF) clusterization algorithm;
• Analogues Complexing (AC)
• Harmonical Rediscretization
• Algorithm on the base of Multilayered Theory of Statistical Decisions (MTSD)
• Group of Adaptive Models Evolution (GAME)

• FAKE GAME Project — Open source. Cross-platform.
• GEvom — Free upon request for academic use. Windows-only.
• GMDH Shell — GMDH-based, predictive analytics and time series forecasting software. Free Academic Licensing and Free Trial version available. Windows-only.
• KnowledgeMiner — Commercial product. Mac OS X-only. Free Demo version available.
• PNN Discovery client — Commercial product.
• Sciengy RPF! — Freeware, Open source.
• wGMDH — Weka plugin, Open source.
• R Package – Open source.
• R Package for regression tasks – Open source.
• Python library of MIA algorithm - Open source.
• Python library of basic GMDH algorithms (COMBI, MULTI, MIA, RIA) - Open source.

• A.G. Ivakhnenko. Heuristic Self-Organization in Problems of Engineering Cybernetics, Automatica, vol.6, 1970 — p. 207-219.
• S.J. Farlow. Self-Organizing Methods in Modelling: GMDH Type Algorithms. New-York, Bazel: Marcel Decker Inc., 1984, 350 p.
• H.R. Madala, A.G. Ivakhnenko. Inductive Learning Algorithms for Complex Systems Modeling. CRC Press, Boca Raton, 1994.

• Library of GMDH books and articles
• Group Method of Data Handling