Documentation for AOLS¶

class AOLS(BaseMSS):
    r"""Accelerated Orthogonal Least Squares Algorithm.

    Build Polynomial NARMAX model using the Accelerated Orthogonal Least-Squares ([1]_).
    This algorithm is based on the Matlab code available on:
    https://github.com/realabolfazl/AOLS/

    The NARMAX model is described as:

    $$
        y_k= F^\ell[y_{k-1}, \dotsc, y_{k-n_y},x_{k-d}, x_{k-d-1},
        \dotsc, x_{k-d-n_x}, e_{k-1}, \dotsc, e_{k-n_e}] + e_k
    $$

    where $n_y\in \mathbb{N}^*$, $n_x \in \mathbb{N}$, $n_e \in \mathbb{N}$,
    are the maximum lags for the system output and input respectively;
    $x_k \in \mathbb{R}^{n_x}$ is the system input and $y_k \in \mathbb{R}^{n_y}$
    is the system output at discrete time $k \in \mathbb{N}^n$;
    $e_k \in \mathbb{R}^{n_e}$ stands for uncertainties and possible noise
    at discrete time $k$. In this case, $\mathcal{F}^\ell$ is some nonlinear function
    of the input and output regressors with nonlinearity degree $\ell \in \mathbb{N}$
    and $d$ is a time delay typically set to $d=1$.

    Parameters
    ----------
    ylag : int, default=2
        The maximum lag of the output.
    xlag : int, default=2
        The maximum lag of the input.
    k : int, default=1
        The sparsity level.
    L : int, default=1
        Number of selected indices per iteration.
    threshold : float, default=1e-9
        The desired accuracy used to stop the iterative selection.

    Examples
    --------
    >>> import numpy as np
    >>> import matplotlib.pyplot as plt
    >>> from sysidentpy.model_structure_selection import AOLS
    >>> from sysidentpy.basis_function import Polynomial
    >>> from sysidentpy.utils.display_results import results
    >>> from sysidentpy.metrics import root_relative_squared_error
    >>> from sysidentpy.utils.generate_data import get_miso_data, get_siso_data
    >>> x_train, x_valid, y_train, y_valid = get_siso_data(n=1000,
    ...                                                    colored_noise=True,
    ...                                                    sigma=0.2,
    ...                                                    train_percentage=90)
    >>> basis_function = Polynomial(degree=2)
    >>> model = AOLS(basis_function=basis_function,
    ...              ylag=2, xlag=2
    ...              )
    >>> model.fit(x_train, y_train)
    >>> yhat = model.predict(x_valid, y_valid)
    >>> rrse = root_relative_squared_error(y_valid, yhat)
    >>> print(rrse)
    0.001993603325328823
    >>> r = pd.DataFrame(
    ...     results(
    ...         model.final_model, model.theta, model.err,
    ...         model.n_terms, err_precision=8, dtype='sci'
    ...         ),
    ...     columns=['Regressors', 'Parameters', 'ERR'])
    >>> print(r)
        Regressors Parameters         ERR
    0        x1(k-2)     0.9000       0.0
    1         y(k-1)     0.1999       0.0
    2  x1(k-1)y(k-1)     0.1000       0.0

    References
    ----------
    - Manuscript: Accelerated Orthogonal Least-Squares for Large-Scale
       Sparse Reconstruction
       https://www.sciencedirect.com/science/article/abs/pii/S1051200418305311
    - Code:
       https://github.com/realabolfazl/AOLS/

    """

    def __init__(
        self,
        *,
        ylag: Union[int, list] = 2,
        xlag: Union[int, list] = 2,
        k: int = 1,
        L: int = 1,
        threshold: float = 10e-10,
        model_type: str = "NARMAX",
        estimator: Estimators = LeastSquares(),
        basis_function: Union[Polynomial, Fourier] = Polynomial(),
    ):
        self.basis_function = basis_function
        self.model_type = model_type
        self.xlag = xlag
        self.ylag = ylag
        self.max_lag = self._get_max_lag()
        self.k = k
        self.L = L
        self.estimator = estimator
        self.threshold = threshold
        self.res = None
        self.n_inputs = None
        self.theta = None
        self.regressor_code = None
        self.pivv = None
        self.final_model = None
        self.n_terms = None
        self.err = None
        self._validate_params()

    def _validate_params(self):
        """Validate input params."""
        if isinstance(self.ylag, int) and self.ylag < 1:
            raise ValueError(f"ylag must be integer and > zero. Got {self.ylag}")

        if isinstance(self.xlag, int) and self.xlag < 1:
            raise ValueError(f"xlag must be integer and > zero. Got {self.xlag}")

        if not isinstance(self.xlag, (int, list)):
            raise ValueError(f"xlag must be integer and > zero. Got {self.xlag}")

        if not isinstance(self.ylag, (int, list)):
            raise ValueError(f"ylag must be integer and > zero. Got {self.ylag}")

        if not isinstance(self.k, int) or self.k < 1:
            raise ValueError(f"k must be integer and > zero. Got {self.k}")

        if not isinstance(self.L, int) or self.L < 1:
            raise ValueError(f"L must be integer and > zero. Got {self.L}")

        if not isinstance(self.threshold, (int, float)) or self.threshold < 0:
            raise ValueError(
                f"threshold must be integer and > zero. Got {self.threshold}"
            )

    def aols(
        self, psi: np.ndarray, y: np.ndarray
    ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
        """Perform the Accelerated Orthogonal Least-Squares algorithm.

        Parameters
        ----------
        psi : ndarray of floats
            The information matrix of the model.
        y : array-like of shape = n_samples
            The target data used in the identification process.

        Returns
        -------
        theta : array-like of shape = number_of_model_elements
            The estimated coefficients for the selected regressors.
        piv : array-like of shape = number_of_model_elements
            Contains the index to put the regressors in the correct order
            based on err values.
        residual_norm : float
            The final residual norm.

        References
        ----------
        - Manuscript: Accelerated Orthogonal Least-Squares for Large-Scale
           Sparse Reconstruction
           https://www.sciencedirect.com/science/article/abs/pii/S1051200418305311

        """
        xp = get_namespace(psi, y)
        target_device = _device(psi, y)
        n, m = psi.shape
        theta = _zeros(xp, (m, 1), dtype=psi.dtype, target_device=target_device)
        r = _copy(xp, xp.reshape(y[self.max_lag :, :], (-1, 1)))
        it = 0
        max_iter = int(min(self.k, np.floor(n / self.L)))
        # selected_indices is a small bookkeeping array; keep as numpy
        selected_indices = np.full(max_iter * self.L, -1, dtype=np.int64)
        basis_matrix = _zeros(
            xp,
            (n, max_iter * self.L),
            dtype=psi.dtype,
            target_device=target_device,
        )
        transformed_psi = _copy(xp, psi)
        eps = np.finfo(float).eps
        numerator = xp.reshape(r.T @ psi, (-1,))
        denominator = xp.sum(psi * transformed_psi, axis=0)
        mask = xp.abs(denominator) > eps
        q = xp.where(mask, numerator / denominator, xp.zeros_like(numerator))
        while (
            float(_to_numpy(_vector_norm(xp, xp.reshape(r, (-1,))))) > self.threshold
            and it < max_iter
        ):
            it = it + 1
            offset = (it - 1) * self.L
            if it > 1:
                basis_vec = xp.reshape(basis_matrix[:, offset], (-1, 1))
                transformed_psi = transformed_psi - basis_vec @ (basis_vec.T @ psi)

            contribution = xp.sum(transformed_psi * transformed_psi, axis=0) * (q**2)
            previous = selected_indices[:offset]
            sub_ind = list(previous[previous >= 0].astype(int))
            # Zero out already-selected and non-finite contributions
            if _is_numpy_namespace(xp):
                contribution[sub_ind] = 0
                contribution[~np.isfinite(contribution)] = 0
            else:
                for si in sub_ind:
                    contribution = _set_element(xp, contribution, si, 0.0)
                contribution = xp.where(
                    xp.isfinite(contribution), contribution, xp.zeros_like(contribution)
                )
            block_size = min(self.L, contribution.shape[0])
            if block_size == 0:
                break
            # argpartition is numpy-specific; convert for this small operation
            if _is_numpy_namespace(xp):
                top_candidates = np.argpartition(contribution, -block_size)[
                    -block_size:
                ]
                current_indices = top_candidates[
                    contribution[top_candidates].argsort()[::-1]
                ]
            else:
                # Fallback: full argsort for non-numpy backends
                sorted_idx = xp.flip(xp.argsort(contribution))
                current_indices = sorted_idx[:block_size]
            current_indices_np = (
                np.asarray(_to_numpy(current_indices), dtype=np.intp)
                if not _is_numpy_namespace(xp)
                else current_indices
            )
            selected_indices[offset : offset + block_size] = current_indices_np
            for i, idx in enumerate(current_indices_np):
                col = int(idx)
                temp = xp.reshape(transformed_psi[:, col], (-1, 1)) * q[col]
                temp_norm = float(_to_numpy(_vector_norm(xp, xp.reshape(temp, (-1,)))))
                if temp_norm <= eps:
                    continue
                basis_matrix[:, offset + i] = xp.reshape(temp / temp_norm, (-1,))
                r = r - temp

                basis_vec = xp.reshape(basis_matrix[:, offset + i], (-1, 1))
                transformed_psi = transformed_psi - basis_vec @ (basis_vec.T @ psi)

                numerator = xp.reshape(r.T @ psi, (-1,))
                denominator = xp.sum(psi * transformed_psi, axis=0)
                mask = xp.abs(denominator) > eps
                q = xp.where(mask, numerator / denominator, xp.zeros_like(numerator))

        selected_indices = selected_indices[selected_indices >= 0].ravel().astype(int)
        residual_norm = float(_to_numpy(_vector_norm(xp, xp.reshape(r, (-1,)))))
        theta[selected_indices] = self.estimator.optimize(
            psi[:, selected_indices], xp.reshape(y[self.max_lag :, 0], (-1, 1))
        )
        if self.L > 1 and len(selected_indices) > self.k:
            if _is_numpy_namespace(xp):
                sorted_local = np.argsort(np.abs(theta[selected_indices]).ravel())[
                    ::-1
                ][: self.k]
            else:
                abs_vals = xp.abs(xp.reshape(theta[selected_indices], (-1,)))
                sorted_local = np.asarray(
                    _to_numpy(xp.flip(xp.argsort(abs_vals))), dtype=np.intp
                )[: self.k]
            top_indices = selected_indices[sorted_local]
            theta_filtered = xp.zeros_like(theta)
            theta_filtered[top_indices] = self.estimator.optimize(
                psi[:, top_indices],
                xp.reshape(y[self.max_lag :, 0], (-1, 1)),
            )
            theta = theta_filtered
            selected_indices = top_indices
            residual_norm = float(
                _to_numpy(
                    _vector_norm(
                        xp,
                        xp.reshape(
                            xp.reshape(y[self.max_lag :, :], (-1, 1))
                            - psi[:, selected_indices] @ theta[selected_indices],
                            (-1,),
                        ),
                    )
                )
            )

        if _is_numpy_namespace(xp):
            pivv = np.argwhere(theta.ravel() != 0).ravel()
        else:
            flat = xp.reshape(theta, (-1,))
            pivv = np.asarray(_to_numpy(xp.nonzero(flat != 0)[0]), dtype=np.intp)
        theta_vals = theta[theta != 0]
        return xp.reshape(theta_vals, (-1, 1)), pivv, residual_norm

    def fit(self, *, X: Optional[np.ndarray] = None, y: Optional[np.ndarray] = None):
        """Fit polynomial NARMAX model using AOLS algorithm.

        The 'fit' function allows a friendly usage by the user.
        Given two arguments, x and y, fit training data.

        Parameters
        ----------
        X : ndarray of floats
            The input data to be used in the training process.
        y : ndarray of floats
            The output data to be used in the training process.

        Returns
        -------
        model : ndarray of int
            The model code representation.
        piv : array-like of shape = number_of_model_elements
            Contains the index to put the regressors in the correct order
            based on err values.
        theta : array-like of shape = number_of_model_elements
            The estimated parameters of the model.
        err : array-like of shape = number_of_model_elements
            The respective ERR calculated for each regressor.
        info_values : array-like of shape = n_regressor
            Vector with values of akaike's information criterion
            for models with N terms (where N is the
            vector position + 1).

        """
        if y is None:
            raise ValueError("y cannot be None")

        self.max_lag = self._get_max_lag()
        lagged_data = build_lagged_matrix(X, y, self.xlag, self.ylag, self.model_type)
        reg_matrix = self.basis_function.fit(
            lagged_data,
            self.max_lag,
            self.ylag,
            self.xlag,
            self.model_type,
            predefined_regressors=None,
        )

        if X is not None:
            self.n_inputs = num_features(X)
        else:
            self.n_inputs = 1  # just to create the regressor space base

        self.regressor_code = self._regressor_space_for_feature_matrix(
            self.n_inputs, n_features=reg_matrix.shape[1]
        )
        self.theta, self.pivv, self.res = self.aols(reg_matrix, y)
        self.pivv = np.asarray(_to_numpy(self.pivv), dtype=np.intp).reshape(-1)
        self.final_model = self.regressor_code[self.pivv, :].copy()

        self.n_terms = self.theta.shape[
            0
        ]  # the number of terms we selected (necessary in the 'results' methods)
        self.err = self.n_terms * [
            0
        ]  # just to use the `results` method. Will be changed in future updates.
        return self

    def predict(
        self,
        *,
        X: Optional[np.ndarray] = None,
        y: Optional[np.ndarray] = None,
        steps_ahead: Optional[int] = None,
        forecast_horizon: int = 0,
    ) -> np.ndarray:
        """Return the predicted values given an input.

        The predict function allows a friendly usage by the user.
        Given a previously trained model, predict values given
        a new set of data.

        Parameters
        ----------
        X : ndarray of floats
            The input data to be used in the prediction process.
        y : ndarray of floats
            The output data to be used in the prediction process.
        steps_ahead : int, optional
            ``None`` selects free-run simulation, 1 selects one-step-ahead
            prediction, and values greater than 1 select n-step-ahead prediction.
        forecast_horizon : int, default=0
            Number of values predicted beyond the initial conditions for a NAR
            free-run prediction when ``X`` is ``None``.

        Returns
        -------
        yhat : ndarray of floats
            The predicted values of the model.

        """
        return super().predict(
            X=X,
            y=y,
            steps_ahead=steps_ahead,
            forecast_horizon=forecast_horizon,
        )

    def _one_step_ahead_prediction(
        self, x_base: Optional[np.ndarray], y: Optional[np.ndarray] = None
    ) -> np.ndarray:
        """Perform the 1-step-ahead prediction of a model.

        Parameters
        ----------
        y : array-like of shape = max_lag
            Initial conditions values of the model
            to start recursive process.
        x : ndarray of floats of shape = n_samples
            Vector with input values to be used in model simulation.

        Returns
        -------
        yhat : ndarray of floats
               The 1-step-ahead predicted values of the model.

        """
        lagged_data = build_lagged_matrix(
            x_base, y, self.xlag, self.ylag, self.model_type
        )
        x_base = self.basis_function.transform(
            lagged_data,
            self.max_lag,
            self.ylag,
            self.xlag,
            self.model_type,
            predefined_regressors=self.pivv[: len(self.final_model)],
        )

        yhat = super()._one_step_ahead_prediction(x_base)
        return get_namespace(yhat).reshape(yhat, (-1, 1))