function update_subspace(obj, x_in)
% Update the low dimensional subspace by adding a new gradient direction x_in to be incorporated into the subspace.;
    if obj.K_current >= obj.M
        return; % no need to update the subspace if it spans the full space;
    end
    if sum(~isfinite(x_in)) > 0
        return; % bad vector! 
    end
    x_in_length = sqrt(sum(x_in.^2));
    if x_in_length < obj.eps
        return; % nothing to do with a short new vector
    end
    xnew = x_in/x_in_length; % make x unit length;
    % Find the component of x pointing out of the existing subspace.;
    % This is done multiple times for numerical stability.;
    % $$ q_{orth}=f_{j}^{'}(x^{t})-P^{t-1}(P^{t-1})^{T} f_{j}^{'}(x^{t}) $$  Equation (7) 
    for i = 1:3
        xnew = xnew - obj.P * (obj.P.'*xnew);
        ss = sqrt(sum(xnew.^2));
        if ss < obj.eps % barely points out of the existing subspace
            return; % no need to add a new direction to the subspace; 
        end 
        xnew = xnew / ss; % make it unit length;
        if ss > 0.1 % largely orthogonal >> good numerical stability
            break;
        end
    end
    % update the subspace with a new column containing the new direction;
    % $$ P^{t}=\left [P^{t-1} \frac{q_{orth}}{\left \|q_{orth}\right\|}\right] $$ Equation (8)
    obj.P(:,obj.K_current+1) = xnew;
    obj.K_current = obj.K_current + 1;

% Restricting the size of the subspace: In case the subspace dimensionality K^{t} exceeds K_{max}, it is collapsed to a new subspace
% that includes only the most recent gradient and position measurements from each subfunction through QR decomposition
    if obj.K_current >= obj.K_max
    % collapse the subspace when it exceeded its maximum allowed size
    % xl is passed explicity to make sure it's used if not in
    % history yet.
        xl = (obj.P.') * (x_in);
        obj.collapse_subspace(xl); % use the function collapse_subspace below
    end
end