Putting the "MAX" back into NARMAX

The library doesn't just raises exceptions; it raises eyebrows and occasionally, the bar for what is assumed to be possible...

Anyways, this GPU-accelerated Python package implements the series of papers that constitute my PhD in mathematics for AI / ML where I present new algorithms for learning symbolic NARMAX machine learning models.

It currently fully implements the machine learning algorithms described my first paper "Arborescent Orthogonal Least Squares Regression (AOrLSR)" and already partially implements elements from the upcoming papers, which include NARMAX-expression compiler (Symbolic Oscillator) and "Dictionary Morphing Orthogonal Least Squares Regression (DMOrLSR)" both based on my "Recursive Forward Orthogonal Least Squares Regression (rFOrLSR)" to fit "Non-Linear Auto-Regressive Moving-Average Exogenous input systems (NARMAX)". So, now that we have covered all the fancy acronyms, we might get into some explanations.
Otherwise jump straight into the library examples/ tutorials.

Note 1 (Github's poor LaTex): Github's LaTex engine is unreliable, so please forgive that certain expressions (especially sums and underlines) are not rendered properly or not at all. All $x$ and $\chi$ in matrix equations are of course vectors and should be underlined (check the readme.md if in doubt).

Note 2 (rFOrLSR): You might ask yourself "how am I even supposed to pronounce rFOrLSR"?
Imagine you're a French pirate trying to pronounce "Airforce". Being French, you'll ignore the last letter in the word, making it "rFOrLS" and being a pirate, you'll say "ARRRRRRRforce" which fully suffices.

NARMAX...who?

• N for Non-linear: Any non-linearity (functions or piecewise definitions) or non-linear combination method (products, divisions, conjunctions, disjunctions, etc.) applied to any of the below types of terms.

• AR for Auto-Regressive: Any type of feeding the system-output $y$ back into the system, via feedback or recursion. Thus, anything containing temporal terms (with $y[k-j]$ terms such as $\tanh(0.5y[k-1] + 0.3y[k-2]y[k-3])$ or $y^a[k-1]x^b[k-3]$) or spatial $y$ terms ($y_1y_2y_3$) or any spatio-temporal combination.

• MA for Moving Average: Any type and distribution of internal noise terms $e[k]$, which accounts for fitting error or any variables not present in the system input. Can also represent internal random states spanning multiple time steps or spatial distributions.

• X for exogenous input: Any type of system-input terms $x$ which can be temporal ($x^n[k-j]$ or $e^{x[k]y[k-2]}$) or spatial ($x_1x_2x_3$, etc) or spatio-temporal (a mixture of both).

A NARMAX system contains thus any arbitrary mixture of the above terms. The provided library examples (being those in the AOrLSR paper) demonstrate the fitting of the below dynamic (=temporal) systems. Those have no other special properties than being stable for an input of amplitude 1 (or potentially more) and being the most ridiculously non-linear systems I came up with at that time to demonstrate the fitting power and flexibility of the algorithms.

1. Linear-in-the-parameter Example

$y[k] = 0.2x[k] + 0.3x^3[k-1] + 0.7|x[k-2]x^2[k-1]| +0.5e^{x[k-3]x[k-2]} - \cos(y[k-1]x[k-2]) -0.4|x[k-1]y^2[k-2]| - y^3[k-3]$

This is essentially a monomial expansion of IIR terms ($y[k-j]$ and $x[k-j]$) also passed through common non-linearities such as abs, cos and exp, yielding a heavily non-linear NARX system.
Tutorial for this example

2. Rational Example

$y[k]=\frac{0.6|x[k]|-0.35x^3[k]-0.3x[k-1]y[k-2]+0.1|y[k-1]|}{1-0.4|x[k]|+0.3|x[k-1]x[k]|-0.2x^3[k-1]+ 0.3y[k-1]x[k-2]}$

This demonstrates that (for NARX systems) rational non-linear models can be fitted by linearizing the terms: $y[k]=\frac{A}{1+B}\iff y[k](1+B)=A⟺y[k]=A-y[k]B$, $A$ and $B$ being linear-in-the-parameter systems such as system 1 in the above example.
Tutorial for this example

3. Expansion-In-An-Expression Example

$y = \text{sgn}(x)(1-\frac{1}{1+|x|A})$ with $A≔\Sigma_{j\in J}\theta _j |x|^j$ and $J\subseteq \mathbb{N}$

This is a memoryless NX (Non-linear exogenous input) system aka a normal scalar function, depending only on $x$. This system shows that NARMAX expansions can be inserted into expressions to impose constraints or system properties (here quick convergence to $\text{sgn(x)}$ and low error around the origin) or obtain complex fitting. This specific expansion is designed to emulate tanh with another continuous rational function. The provided code also demonstrates how to choose the number of terms in such an expansion and how to create a custom validation function.
Tutorial for this example

4. MIMO & MISO Example

$y_1[k] = 0.2 x_1[k] + 0.3 x_2^3[k] + 0.7 |x_3[k]| + 0.5 x_2[k-3] x_1[k-2] - 0.3 y_2[k-1] x_2^2[k-2] - 0.8 |x_3[k-1] y_1[k-2]| - 0.7 x_1[k-1] x_2^2[k-1]$ $y_2[k] = 0.3 x_1[k-1] + 0.5 x_3^3[k] + 0.7 |y_1[k-1]| + 0.6 y_1[k-3] x_1[k-2] - 0.4 y_1[k-1] x_3^2[k-2] - 0.9 |x_3[k-1] y_2[k-2]| - 0.7 x_3[k-1] x_2^2[k-1]$

This is a MIMO (Multiple Input Multiple Output) system / function with 3 input channels / variables and 2 output channels / variables, which is constituted by two MISO (Multiple Input Single Output) systems / functions: one per output. This demonstrates that the rFOrLSR can fit systems / functions with an arbitrary input and output dimensionality: $\mathbb{R}^n \rightarrow \mathbb{R}^m$ (in this example $\mathbb{R}^3 \rightarrow \mathbb{R}^2$).
Tutorial for this example

The NARMAX fitting steps:

1. Expansion Type Selection: As usual in machine learning, one must first choose an appropriate expansion type (FIR, IIR, Monomially-Expanded IIR, RBF, Wavelet, arbitrary non-linearities, etc.). As expected, the model quality strongly depends on how relevant the chosen expansion is. The advantage of this library's rFOrLSR is that any type of expansion and any mix of expansions is supported, as the rFOrLSR is based on vector matching methods.
   This is achieved by creating a fitting "dictionary" $D_C \in \mathbb{R}^{p \times n_C}$ (Pytorch matrix) containing the candidate regressors $\underline{\varphi}_k \in \mathbb{R}^{p}$ stacked column-wise and passing it to the library.

2. Model Structure Detection: The first real challenge is to find the correct regressors from all those present in the user-defined dictionary $D_C$, as most system behaviors can be sufficiently well described with very few terms.
   To illustrate, the first system above contains a single cosine term which needs to be retrieved from the set of cosines with relevant monomial expansions as arguments.

3. Model Parameter Estimation: Finally, once the correct expression terms are selected, their regression (= scaling) coefficients must be chosen. The rFOrLSR's optimization criterion is least squares, however my further papers will switch to more robust (noise and glitch resistant) optimization methods.

NARMAX...why?

This section is dedicated to all the people who asked me something along the lines of "but AI is currently the thing, why not use neural networks like everyone else?".

First of all, I'd like to point out that artificial neural networks (ANNs), including our friends Qwen, ChatGPT, DeepSeek, etc. are some subclass of NARMAXes. LLMs based on Auto-regressive Transformers are certainly non-linear (N) and auto-regressive (AR), have some internal random states affecting the computations (MA) and take exogenous inputs (X) being whatever you ask them.

Thus, this section is really about symbolic fitting vs black-box networks (which includes neural networks and other algorithm classes this library can fit such as RBF networks and to some extend Wavelet networks, etc.)

Interpretability & Provable Correctness: Symbolic models often result in equations having a clear physical or at least mathematical interpretation. This is important in situations where understanding the underlying relationships between inputs and outputs is required, such as physics, biology, or engineering. Black-box networks, however, do not provide easily interpretable models, which in turn can't be proved to be correct.
To illustrate, NARMAX models are used, amongst others, in the field of robotics, as they allow for example to determine a) which part of the system is affected by b) the input of which sensor c) at what time lag and d) how.

Prior Knowledge & Constraints: Symbolic models allow the incorporation of domain knowledge and constraints into the model structure.
To illustrate, if the system is known to be linear, only linear terms are added to the dictionary or if the system is oscillatory in nature one fills the dictionary with sine and cosine regressors.
My rFOrLSR even allows to impose regressors to further constraint the model and impose user knowledge.

Data Efficiency: Symbolic models require very little data for training, which allows to fit them in scarce data scenarios or even keep them updated in real-time. This is a clear advantage when data acquisition is expensive (like for biological processes) or when computational resources for fitting are limited. To illustrate, many NARMAX papers fit their example models with a 500 floats dataset, whereas neural networks require many GB or TB of data.

Noise and Outliers Handling: Symbolic models can be more robust to noise and outliers in the data. In particular, NARMAX systems including MA regressors can create a model of the noise structure to further stabilize the model and guarantee bias-free parameter estimation. Black-box networks can be sensitive to noisy data and can overfit to outliers, leading to poor generalization.

Computational Efficiency: Symbolic models often comprise very few terms and are thus computationally cheap, which can be of interest for real-time applications or embedded systems.

Extrapolation / generalization: Symbolic models may perform better in extrapolation tasks, where predictions need to be made beyond the training data's value range than black-box networks. Indeed, if the model contains the correct regressors, the correct or almost correct equation can be obtained. The fitting error would be very low and the model would perform well on unseen data and data ranges.

Combinations: Also, both (symbolic and blackbox) models are not exclusive. A small and efficient symbolic model can be put first in the processing chain to allow a drastic size reduction of the following neural network, which can fit the symbolic model's residuals.
To illustrate, a small symbolic model explaining 70% of the data variance, could allow to reduce the number of layers and neurons by X amount.

About the library

As described above, the library makes optimal sparse least squares fitting of any vectors / regressors $\underline{\varphi}_k \in \mathbb{R}^{p}$ passed by the user for the given system response $\underline{y}\in \mathbb{R}^{p}$. All provided examples are with severely non-linear auto-regressive systems, however, the user can pass dictionaries containing any type of functions (multivariate, discontinuous, discrete, acausal, etc. - whatever that outputs vectors with numbers really). This is thus a very general mathematical framework supporting for example FIR/MA, IIR/AR, wavelet-, RBF, polynomial, etc. fitting.

Installation

Clone the repository and install locally with pip:

git clone https://github.com/Stee-T/NARMAX.git
pip install ./NARMAX

Library content

Regressor Constructors

The Regressor-CTors allow to easily generate regressors to be horizontally stacked in the dictionary $D_C$ passed to the arborescence. They cover the most common types of expansions like FIR/MA, IIR/AR, Power-series, oscillations, Radial basis functions, Wavelet, etc. Any mixture and any custom user function/regressor can be put in the dictionary to be used for fitting.

• Lagger: Creates delayed $\underline{x}, \underline{y}$ and optionally $\underline{e}$ terms up to the desired order.
  This is the only CTor needed for linear FIR and IIR systems.
  Example: $\underline{x} → [\ [\underline{x}[k-j]\ ]_{j=0}^{n} = [\ \underline{x}[k], \underline{x}[k-1], \underline{x}[k-2],...\ ]$

• Expander: Monomially expands the input by creating all product combinations up to a certain combined power of the passed input.
  Example: Set of two terms with 3rd order expansion: $[\underline{x}_1,\ \underline{x}_2 ]→[\ \underline{x}_1,\ \underline{x}_2,\ \underline{x}_1^2,\ \underline{x}_2^2,\ \underline{x}_1\cdot\underline{x}_2,\ \underline{x}_1^3,\ \underline{x}_2^3,\ \underline{x}_1^2\cdot\underline{x}_2,\ \underline{x}_1\cdot\underline{x}_2^2\ ]$

• NonLinearizer: Applies non-linearities to the input in addition to allowing terms to be in the denominator for rational functions.
  Note that this is limited to elementwise functions, other fancy things must currently be done by hand.
  Example: $\underline{x}→[\ \underline{x}, f_1(\underline{x}), f_2(\underline{x}), \frac{1}{f_3(\underline{x})},...\ ]$

• Booler: TODO Example: $[\underline{x}_1,\ \underline{x}_2\ ]→[\ \underline{x}_1,\ \underline{x}_2,\ !\underline{x}_1,\ !\underline{x}_2,\ \underline{x}_1 && \underline{x}_2,\ \underline{x}_1 || \ \underline{x}_2,\ \underline{x}_1^ \underline{x}_2 ]$

• SmoothedDeriver: (already implemented but not compatible with the updated API) Constructor for derivatives of arbitrary order with arbitrarily strong smoothing. This can be either used to generate derivative regressors (see differential equations) or for Ultra-Orthogonal Least Squares fitting in dedicated Sobolev spaces. Idea: Concatenating the regressor and their derivatives reduces overfitting as there is more information on the regressors as is classically done in solving of differential equations (see fitting in Sobolev spaces, etc).

• Oscillator: (Coming soon) Constructor for common oscillations like cos, sin with options for different frequencies, phases and amplitude decays.

• RBF: (Coming soon) Constructor for common Radial Basis Functions

• Wavelet: (Coming soon) Constructor for common wavelet functions

• Input Signal Generators (Coming soon)
  Probably something facilitating the construction of input signals $\underline{x}$ like Maximum Length Sequences (MLS), multi-level sequences, etc.

• Whatever anyone feels like contributing like Chirps or orthogonal bases for specific types of analysis :)

Analysis tools

• (linear) IIR Analysis Tools: The library contains functions to transform the rFOrLSR output into a standard $\underline{b}$, $\underline{a}$ coefficent-form for further production use. Additionally, several convenience plotting functions are provided such as magnitude and phase response, and pole-zero plots.

• MaxLagPlotter: Tool designed to estimate the maximum lags and the expansion order CTor for $\underline{x}$ and $\underline{y}$ terms via polynomial NARMAXes.
  Uses:

  • Allows automating the selection of the parameters passed to the Lagger and Expander CTors for polynomial fitting.

  • Helps eliminating unnecessary regressors from the dictionary before starting the arborescence for great speed-ups.
    (Note that this is still under development and in some cases much slower than running the arborescence with a large dictionary.)

  • Allows to analyze the percentage of explained variance per lag and per expansion order.

• NonLinearSpectrum: (Coming soon).
  Essentially a Fourier transform for Nonlinear systems based on polynomial NARMAXes.

rFOrLSR & Arborescence

The library's main object is the arborescence which takes in the user-defined dictionaries (imposed regressors $D_S$ and candidate regressors $D_C$) and the desired system input $\underline{y}$ for the fitting. The arborescence, being a breadth-first tree search in the solution space spanned by the dictionary $D_C$ columns, repeatedly calls the rFOrLSR to iteratively refine the solution. The more "levels" the arborescence has, the more search space is traversed and the higher the probability of finding the optimal solution.
For more on that, see the below "More on the Algorithms" section.

Symbolic Oscillator

It's a python-object, internally compiled by the library to allow the use of user-defined symbolic NARMAX systems for data generation and processing. In Machine Learning / Statistic-terms, it's basically the "executable" object containing the fitted model.

The NARMAX.SymbolicOscillator object, thus, allows the user to create a symbolic NARMAX system from the provided lists of regressors (list of strings), regression coefficients (float array) and regression non-linearities (lists of function pointers). The SymbOsc object has an internal parsing engine that analyzes the provided lists, checks for correctness, and forwards it to the internal compiler, which will generate the necessary code for the symbolic NARMAX system with optimal CPU/GPU dispatching (currently only partially implemented, being part of my next paper).

Dictionary Morphing

→ COMING SOON
Essentially, an extension of the framework where the passed regressors are adapted to the system in real time by the rFOrLSR under very mild conditions. To illustrate, the user can pass any elementwise function $f_⦿$ and the arborescence will morph it by adding scaling coefficients $\xi_j$ and arbitrary linear combinations of dictionary terms as arguments. The morphing is thus of the form: $f_⦿ \left(\underline{\chi}_0 \right) \rightarrow f_⦿ \left(\Sigma_{j=0}^r \xi_j \underline{\chi}_j\right)$

Examples:

• Frequency and phase estimation:
  $\cos(x) → \cos(2.1x + 0.2)$ with $\xi_0 = 2.1$, $\xi_1 = 0.2$, $\underline{\chi}_1 = \underline{1}$

• Wavelet translation and dilation:
  $\psi(x) → \psi(ax + b)$ with dilation coefficient $\xi_0 = a$ and translation coefficient $\xi_1 = b$, $\underline{\chi}_1 = \underline{1}$

• Radial Basis Function arguments:
  $\varphi(\left|x_0 - c_0\right|) → \varphi\left(\lVert \Sigma_{j=0}^r (a_j x_j-b_jc_j) \rVert \right)$ with scaling factors $\xi_j = a_j$ and $\xi_j = b_j$, $\underline{c}_j = \underline{1}$

• Arbitrary elementwise functions
  Arbitrary functions where changing coefficients or adding arguments is of interest. Examples include oscillations where the decay / frequency needs to be varied, chirps whose speed / start or end frequency must be varied etc.

Library & PhD Roadmap

Symbolic NARMAX systems (such as Polynomials, Wavelets, RBFs (like GMMs), FIRs, IIRs, etc.) are the main alternative to artificial neural networks (ANNs) for black-box fitting or data and stochastic process analysis and processing. For system fitting, they are, however, used essentially for small to mid-size systems (or in situation where interpretability and provable correctness is paramount), since SotA methods are either very limited in fitting capabilities (FOrLSR-based methods) or prohibitively slow (genetic algorithms). Since NARMAX models have many advantages over ANNs such as interpretability, ease of modification and analysis, computational efficiency, and marginal data requirements compared to ANNs (500 or fewer floats vs GB/TB/PB of data), it’s of interest to push the limits of the FOrLSR (Forward Orthogonal Least Squares Regression) framework. This PhD’s goal is to make the FOrLSR framework capable of reliably fitting increasingly more complex systems with each paper, such that the more efficient and expressive symbolic NARMAX systems can be used instead of the black-box neural networks in an increased number of applications. The focus is on supporting general symbolic NARMAX expansions to allow users to fit their domain-specific functions and expressions, to customize fitting quality and interpretability.

Paper 1 (AOrLSR: Arborescent Orthogonal Least Squares Regression)

Modernization of the mathematical framework, greatly increasing fitting quality and speed, and allowing the use with “big-data” on GPUs and distributed accelerator systems. Status: Math: 100%. Code: 100%. Paper link here. Equivalence to Neural Networks: Much better model strucutre detection, optimization for GPU

Paper 2 (Symbolic Oscillator)

JIT-compiler for generalized NARMAX systems, allowing to measure the domain-specific user-defined functions during fitting and validation to further increase the complexity of reliably fitted systems, in addition to allowing their integration in ML-pipelines with an optimal dispatch between CPU and GPU. Status: Math: 80%. Code: 70% Equivalence to Neural Networks: Forward Pass (Evaluating the model)

Paper 3 (Dictionary Morphing)

Optimization pipeline allowing the framework to adapt the user-given domain-specific functions to the system (or data or stochastic process) for more precise and smaller models. (Paper 1 only allows to use the user-given functions as-is, while this allows the library to optimize them further) Status: Math: 90%. Code: 70% Equivalence to Neural Networks: (Fine-tuning the model parameters)

Paper 4 (Currenlty Unnamed)

Batching and Robust fitting. Batching is the continuation of the first paper, where the framework is pushed to support big data fitting like modern AI requires. Robust fitting: Least-squares is optimal in the case of gaussian noise but is very sensitive to large outliers (Sudden unplanned events, sensor glitches, etc.), robust fitting makes the fitting stable against such problems (often used in robotics and self-driving cars). Status: Math: 20%. Code: 0% Equivalence to Neural Networks: Batching (Fine tuning) + Data-preprocessing (outlier removal and heavy noise compensation built into the fitting algorithm)

Paper 5 (LASAGNE: Lean Arborescent Search Algorithm for Generalized NARMAX Expressions)

Merge all the previous papers, Inserts the algorithms from all the previous papers into the 1st one to tremendously increase the complexity of reliably fitted systems in addition to providing supplementary sparsity theorems and optimizations to accelerate fitting. Status: Math: 30%. Code: 0% Equivalence to Neural Networks: Putting all previous papers together into one cohesive algorithm.

Paper 6 (Multi-layer LASAGNE for Deep NARMAX-Networks)

Reinterprets the generalized symbolic NARMAX model as a single-layer network and generalizes it to be multi-layer, with the advantage over traditional artificial neural networks (ANNs) that the network choses its optimal structure (size, connections and activation functions) by itself, while remaining interpretable (human readable-equations). Status: Math: 0%. Code: 0% Equivalence to Neural Networks: Multi-layer Neural Network with build-in architecture selection and Data-cleaning

More about the Algorithms

Recursive Forward Orthogonal Least Squares Regression - rFOrLSR

The rFOrLSR (described in the first section of the AOrLSR paper) is a recursive matrix form of the classical FOrLSR (Forward Orthogonal Least Squares Regression) algorithm. Thus, the double for-loop inside the FOrLSR is replaced by a single matrix operation, which is recursively updated. The matrix-form allows the large-scale BLAS-like or GPU-based optimizations offered by this library, while the recursion flattens the FOrLSR’s complexity from quadratic in model length $n_r$ ( $O(n_r^2)$ ) to linear ( $O(n_r)$ ). Then, the remaining FOrLSR procedure is rearranged and vectorized, in addition to allowing to efficiently impose regressors (which is not a native FOrLSR feature). These optimizations greatly reduce its computational cost, allowing use it as node processing procedure in large search trees like the proposed arborescence.

The (r)FOrLSR is a greedy algorithm selecting iteratively the regressors the most similar to the desired system output $\underline{y}$ to form the NARMAX model. After each regressor is added, to avoid describing the same system output information (aka $\underline{y}$ variance) repeatedly, the information is removed (via orthogonalization) from the remaining regressors in the dictionary of candidates $D_C$.
The (r)FOrLSR is thus a greedy algorithm, since at each iteration, it takes the locally best choice, maximizing the increase in explained output variance without considering the global search space, which could contain shorter and/or lower validation error models. When a suboptimal choice is made, the following regressors are often chosen to compensate for that error, leading to poor models. For this reason, an arborescence spanning a large search space is proposed by my papers/library, see below.

Arborescent Orthogonal Least Squares Regression - AOrLSR

The AOrLSR (described in the second section of the AOrLSR paper) mitigates the (r)FOrLSR's greediness by tactically imposing regressors in the model and aborting or skipping regressions once their results become predictable. My paper proposes several theorems and corollaries dictating which regressions must be computed and which not. Most of the time, only a marginal fraction (<0.5-5%) of the search space must be computed to be traversed. This allows to find expansions with fewer regressors (sparse fitting) and/or lower error by choosing different regressors.

The arborescence also provides a validation mechanism to select the best model. Once the search space is traversed (or if the search is ended manually), the validation procedure iterates over all computed models and runs the user-defined validation function on the models with the smallest number of regressors (sparsity criterion), diminishing overfitting and model computational expense. The arborescence object accepts a function pointer and a dictionary for the validation, allowing arbitrary validation methods on arbitrary validation data, such that any user-defined criterion can be used, such as signal mean-square error, frequency-domain similarities, computational expense estimations, etc.

The AOrLSR as sparsifying solver

The (r)FOrLSR and the AOrLSR are, on an abstract level, linear equation system solvers, maximizing the solution-vector sparsity by choosing the most relevant LHS columns, while minimizing the mean squared error with respect to the RHS $\underline{y}$. The LHS is $M≔D_S⫲D_C$, where $D_S$ allows to impose columns to the selection while the AOrLSR is free to choose any desired columns from $D_C$. $M$’s dimensions being arbitrary, the system can be under-, well- or over-determined.

To illustrate, be the following system with $p \in \mathbb{N}^*$: $M \underline{x} = \underline{y}$

The AOrLSR with a $1 %$ error tolerance could chose the regressor (here columns of $M$) index-set $\mathcal{L}={1,3}$ and $\underline{\hat{\theta}} \in \mathbb{R}^{|\mathcal{L}|=2}$, such that $\underline{x}[\mathcal{L}]:=\hat{\underline{\theta}}\iff \underline{x}=[\hat{\theta}_1, 0, \hat{\theta}_2, 0,0]^T$, while a $0.1 %$ tolerance could yield $\mathcal{L}={1,3,4}$ and $\underline{\hat{\theta}} \in \mathbb{R}^{|\mathcal{L}|=3}$, thus $\underline{x}=\left[\hat{\theta}_1, 0, \hat{\theta}_2, \hat{\theta}_3, 0\right]^T$.

The linear equation system is thus reduced to the most relevant subset of $M$'s columns and $\underline{x}$ entries to describe the RHS $\underline{y}$: Thus, $[ \underline{m} _j]_{j \in \mathcal{L}} \underline{\hat{\theta}}=\underline{y}$ is solved via a QR-like decomposition.

Dictionary Morphing Orthogonal Least Squares Regression - DMOrLSR

COMING SOON

This part of the library is described in my "Dictionary Morphing Orthogonal Least Squares Regression (DMOrLSR)" paper. Many regressors of interest are (elementwise) non-linearities taking one or more arguments $f_⦿(\xi)$. It is, however, very unlikely that the passed candidate dictionary $D_C$ contains the optimal non-linearities of the form $f_⦿\left(\Sigma_{j=1}^r \xi_j \underline{\chi}_j\right)$, as this would require an quasi-infinite number of infinitesimal increments for each $\xi_j$ and a combinatorial order of $\underline{\chi}_j$ combinations for each $r&gt;0$. This would overflow RAM (see example 1 in the paper) and take insane computing times due to the enormous candidate dictionary $D_C$. It is thus more efficient to scan $D_C$ for the most fitting regressor and add arguments and coefficients to adapt the regressor to the system.
(See examples in the morphing section above)

For comparison, classic expansions like Fourier and Laplace fit the entire given dictionary (= imposed terms and expansion order), whereas the (r)FOrLSR selects only the most relevant terms (= system dependent term selection and expansion order). The next logical step is adapting the selected regressors to the system, making the expansion fully system dependent.

The morphing procedure comprises four steps. After the rFOrLSR regressor selection, the regressor is parsed to deduce the non-linearity $f$ and its argument ${\underline{\chi}}_0$. If the function is morphable, the procedure continues with a genetic vector-space generating algorithm (GenVSGen) determining how many and which supplementary arguments ${\underline{\chi}}_j$ improve the fitting. Finally, an infinitesimal optimizer determines the optimal coefficients $\underline{\xi}$. Once the rFOrLSR has selected and morphed sufficiently many regressors to meet its precision threshold, the final infinitesimal optimizer adapts the regression coefficients $\underline{\hat{\theta}}$ and the respective function arguments $\underline{\xi}_j$.

The morphing works with closed-form solutions of the gradient and Hessian, such that no AutoGrad is needed, allowing many orders of magnitude faster fitting procedures on arbitrary functions.

The morphing allows:

• Fitting much more complex expressions (= more expressive models)
• Sparser fitting (= smaller resulting expressions for a desired error tolerance)
• Lower error (= smaller mean squared error for same length expressions)

Contributions & Licensing

Contributions

All contributions to the library are welcome:

• Regressor Constructors for whatever you happen to be fitting!
• Unit Tests
• Documentation and supplementary examples
• Non-linear analysis tools

Please open a pull request with a clear description of the change and relevant tests. We may request revisions for style, correctness, or maintainability.

Licensing

The NARMAX library is licensed under the 3-Clause BSD License with Copyright (c) 2023-2026 to Stéphane J.P.S. Thunus. By submitting a contribution, you agree that your submission may be incorporated into this project under the same license, and that you have the right to submit it.

This library is thus open-source and free in all senses of the word: Free as in usable without any payment and free in the sense of "this is a free country and I do what I want", **eagle screeching in the background as I wipe a tear of joy from my eyes shimmering from hope while the summer wind gently caresses my hair with soft whispers of boundless possibilities**.

Comment on AI Submissions

Any form of AI-assisted or AI-Agent submission are welcome: bug reports, contributing Unit Tests, Code, docs, etc.. However, as for human contributors, the account submitting the PR / bug report is responsible and will be banned if low quality trash is repeatedly submitted.