function [beta_ridge] = fit_ridge(X, y, lambda)

    [n, p] = size(X); % number of observations
    
    if (p < n)
        hSigma = X' * X / n; % sample covariance matrix
        beta_ridge = (hSigma + lambda * eye(p)) \ (X' * y) / n; % ridge estimator using primal form
    else
        gram_mat = X * X' / n; % gram matrix
        beta_ridge = X' * ((gram_mat + lambda * eye(n)) \ y) / n; % ridge estimator using dual form
    end

    
end