ARD Reference

def fit(self,X,y):
    '''
    Fit the ARD regression model to the data.

    Parameters
    ----------
    X : {array-like, sparse matrix}, shape (n_samples, n_features)
        Training data, matrix of explanatory variables.

    y : array-like, shape (n_samples,)
        Target values.

    Returns
    -------
    self : object
        Returns the instance itself.
    '''
    X, y = check_X_y(X, y, dtype=np.float64, y_numeric=True)
    # ── Route to the appropriate algorithm ─────────────────
    if self.algorithm == 'em':
        return self._fit_em(X, y)
    # ── Sequential SBL (Tipping & Faul 2003) ──────────────
    # Save original-scale data for cross-validation (prevents leakage
    # from full-dataset centering into per-fold evaluation).
    X_orig, y_orig = X.copy(), y.copy()
    X, y, X_mean, y_mean, X_std = self._center_data(X, y)
    n_samples, n_features = X.shape

    # Initialize history storage
    self.history_ = {
        'iterations': [],
        'cv_scores': [],
        'n_features': [],
        'states': [] if self.store_history else None
    }
    self.best_iteration_ = None
    self.best_cv_score_ = None

    # For backward compatibility
    cv_list = []
    current_r = 0

    #  precompute X'*Y , X'*X for faster iterations & allocate memory for
    #  sparsity & quality vectors
    XY     = np.dot(X.T,y)
    XX     = np.dot(X.T,X)
    XXd    = np.diag(XX)

    #  initialise precision of noise & and coefficients
    var_y  = np.var(y)

    # check that variance is non zero !!!
    if var_y == 0 :
        beta = 1e-2
    else:
        beta = 1. / np.var(y)

    A      = np.inf * np.ones(n_features)
    active = np.zeros(n_features , dtype = bool)

    # in case of almost perfect multicollinearity between some features
    # start from feature 0
    if np.sum( XXd - X_mean**2 < np.finfo(np.float32).eps ) > 0:
        A[0]       = np.finfo(np.float16).eps
        active[0]  = True
    else:
        # start from a single basis vector with largest projection on targets
        proj  = XY**2 / XXd
        start = np.argmax(proj)
        active[start] = True
        A[start]      = XXd[start]/( proj[start] - var_y)

    warning_flag = 0

    # Store best model state for retrospective selection
    best_state = None
    best_cv_score = -np.inf
    best_iteration = -1

    for i in range(self.n_iter):
        XXa     = XX[active,:][:,active]
        XYa     = XY[active]
        Aa      =  A[active]

        # mean & covariance of posterior distribution
        Mn,Ri,cholesky  = self._posterior_dist(Aa,beta,XXa,XYa)
        if cholesky:
            Sdiag  = np.sum(Ri**2,0)
        else:
            Sdiag  = np.copy(np.diag(Ri))
            warning_flag += 1

        # raise warning in case cholesky failes
        if warning_flag == 1:
            warnings.warn(("Cholesky decomposition failed ! Algorithm uses pinvh, "
                           "which is significantly slower, if you use RVR it "
                           "is advised to change parameters of kernel"))

        # compute quality & sparsity parameters
        s,q,S,Q = self._sparsity_quality(XX,XXd,XY,XYa,Aa,Ri,active,beta,cholesky)

        # update precision parameter for noise distribution
        rss     = np.sum( ( y - np.dot(X[:,active] , Mn) )**2 )
        beta    = n_samples - np.sum(active) + np.sum(Aa * Sdiag )
        beta   /= ( rss + np.finfo(np.float32).eps )

        # update precision parameters of coefficients
        A,converged  = update_precisions(Q,S,q,s,A,active,self.tol,
                                         n_samples,False)

        # --- Cross-validation scoring (if enabled) ---
        cv_score = None
        if self.cv:
            cv_score = self._compute_cv_score(X_orig, y_orig, active, beta, A, XX, XY,
                                             X_mean, y_mean, X_std)

            # Store in history
            self.history_['iterations'].append(i)
            self.history_['cv_scores'].append(cv_score)
            self.history_['n_features'].append(np.sum(active))

            if self.store_history:
                # Store model state
                state = {
                    'active': active.copy(),
                    'coef': np.zeros(n_features),
                    'coef_active': Mn.copy(),
                    'lambda': A.copy(),
                    'alpha': beta,
                    'sigma': Ri.copy() if not cholesky else None,
                    'cholesky': cholesky
                }
                state['coef'][active] = Mn
                self.history_['states'].append(state)

            # Update best model for retrospective selection
            if cv_score is not None and cv_score > best_cv_score:
                best_cv_score = cv_score
                best_iteration = i
                # Store best state
                best_state = {
                    'active': active.copy(),
                    'A': A.copy(),
                    'beta': beta,
                    'XXa': XXa.copy(),
                    'XYa': XYa.copy(),
                    'Aa': Aa.copy(),
                    'X_mean': X_mean.copy(),
                    'y_mean': y_mean,
                    'X_std': X_std.copy()
                }

            # For backward compatibility
            cv_list.append(cv_score)
            if i == 0:
                current_r = cv_score if cv_score is not None else 0

            # Print CV status
            if self.verbose:
                print(f'Iteration: {i:<4}  CV Score: {cv_score:.6f}  '
                      f'Active features: {np.sum(active)}')

        # --- Verbose output for main iteration progress ---
        if self.verbose and not self.cv:
            # Use f-string for consistency and clarity
            print(f'Iteration: {i:<5}, number of features remaining: {np.sum(active)}')

        # --- Check for convergence (ARD criteria only, no CV early stopping) ---
        # Note: CV-based early stopping is disabled when retrospective_selection=True
        if converged or i == self.n_iter - 1:
            if self.verbose:
                print(f'Finished ARD iterations at iteration {i+1}.')
                if converged:
                    print('Algorithm converged (ARD criteria).')
                elif i == self.n_iter - 1:
                    print('Reached maximum number of iterations without full convergence.')

            # If using retrospective selection and we have a best state, restore it
            if self.cv and self.retrospective_selection and best_state is not None:
                if self.verbose:
                    print(f'Restoring best model from iteration {best_iteration} '
                          f'with CV score: {best_cv_score:.6f}')

                # Restore best state
                active = best_state['active']
                A = best_state['A']
                beta = best_state['beta']
                XXa = best_state['XXa']
                XYa = best_state['XYa']
                Aa = best_state['Aa']
                X_mean = best_state['X_mean']
                y_mean = best_state['y_mean']
                X_std = best_state['X_std']

                # Update best iteration attributes
                self.best_iteration_ = best_iteration
                self.best_cv_score_ = best_cv_score

            # Break only if not using retrospective selection, if we're at the last iteration, or if converged
            if not self.retrospective_selection or i == self.n_iter - 1 or converged:
                break


    # after last update of alpha & beta update parameters
    # of posterior distribution
    XXa,XYa,Aa         = XX[active,:][:,active],XY[active],A[active]
    Mn, Sn, cholesky   = self._posterior_dist(Aa,beta,XXa,XYa,True)
    self.coef_         = np.zeros(n_features)
    self.coef_[active] = Mn
    self.sigma_        = Sn
    self.active_       = active
    self.lambda_       = A
    self.alpha_        = beta

    # Post-fit pruning: zero out coefficients whose estimated
    # precision exceeds threshold_lambda (mirrors sklearn ARDRegression).
    if np.isfinite(self.threshold_lambda):
        prune_mask = (self.lambda_ > self.threshold_lambda) & self.active_
        self.coef_[prune_mask] = 0.0
        self.active_[prune_mask] = False
        if self.verbose and np.any(prune_mask):
            print(f"  Pruned {np.sum(prune_mask)} features with "
                  f"lambda > {self.threshold_lambda:.0f}")

    self._set_intercept(X_mean,y_mean,X_std)

    # Store scores for backward compatibility
    self.scores_ = cv_list if self.cv else []

    if self.cv and self.verbose:
        print(('Number of features in the model: {0}').format(np.sum(active)))
    return self

def predict_dist(self,X):
    '''
    Computes predictive distribution for test set.
    Predictive distribution for each data point is one dimensional
    Gaussian and therefore is characterised by mean and variance.

    Parameters
    ----------
    X : {array-like, sparse matrix}, shape (n_samples_test, n_features)
        Test data, matrix of explanatory variables.

    Returns
    -------
    y_hat : array, shape (n_samples_test,)
        Estimated values of targets on the test set (mean of the predictive distribution).

    var_hat : array, shape (n_samples_test,)
        Variance of the predictive distribution.
    '''
    y_hat     = self._decision_function(X)
    var_hat   = 1./self.alpha_
    var_hat  += np.sum( np.dot(X[:,self.active_],self.sigma_) * X[:,self.active_], axis = 1)
    return y_hat, var_hat