

function collapse_subspace(obj, xl)
% collapse the subspace to its smallest dimensionality. xl is a new direction that may not 
% be in the history yet, hence to be included   
    Pl = zeros(obj.K_max,obj.K_max); % the projection matrix from old to new subspace;
    % yy will hold all the directions to pack into the subspace initialized with random noise 
    % to span K_min dimensions even if not all the subfunctions are active yet;
    yy = randn(obj.K_max,obj.K_min);
    % most recent position & gradient for all active subfunctions and the current position & gradient 
    % that are not saved in the history yet;
    yz = [obj.last_df(:,obj.active), obj.last_theta(:,obj.active), xl, (obj.P.') * (obj.theta)];
    yy(:,1:size(yz, 2)) = yz;
    [Pl(:,1:obj.K_min), ~] = qr(yy, 0);
    obj.P = (obj.P) * (Pl);% update the subspace;
    obj.apply_subspace_transformation(Pl.', Pl);% apply Pl (projection matrix from old to new basis) to all the history terms;
    obj.K_current = obj.K_min; % update the stored subspace size;
    obj.reorthogonalize_subspace();% re-orthogonalize the subspace if small errors are accumulated
end

function reorthogonalize_subspace(obj)
    % check if the subspace has become non-orthogonal
    subspace_eigs = eig(obj.P.'*obj.P);
    if max(subspace_eigs) <= 1 + obj.eps
        return
    end
    % performing QR to the non-orthogonal subspace  
    [Porth,~] = qr(obj.P(:,1:obj.K_current),0); % Orthogonal-triangular decomposition.
    Pl = zeros(obj.K_max, obj.K_max);
    Pl(:,1:obj.K_current) = obj.P.' * Porth;
    obj.P(:,1:obj.K_current) = Porth;% update the subspace;
    obj.apply_subspace_transformation(Pl.',Pl); % apply Pl (projection matrix from old to new basis) to all history terms;
end

function apply_subspace_transformation(obj,T_left,T_right)
% Apply change-of-subspace transformation.  This function is called when the subspace is collapsed to project 
% into the new lower dimensional subspace.;
    % T_left - The covariant subspace projection matrix.;
    % T_right - The contravariant subspace projection matrix.;
    % (Note that for orthogonal subspace T_left = T_right)
    [tt, ss] = size(T_left);
    % project history terms into new subspace;
    obj.last_df = (T_right.') * (obj.last_df);
    obj.last_theta = (T_left) * (obj.last_theta);
    obj.hist_deltadf = obj.reshape_wrapper(T_right.' * obj.reshape_wrapper(obj.hist_deltadf, [ss,-1]), [tt,-1,obj.N]);
    obj.hist_deltatheta = obj.reshape_wrapper(T_left * obj.reshape_wrapper(obj.hist_deltatheta, [ss, -1]), [tt,-1,obj.N]);
    % project stored hessian for each subfunction in to new subspace;
    obj.b = obj.reshape_wrapper(T_right.' * obj.reshape_wrapper(obj.b, [ss,-1]), [tt, 2.*obj.max_history,obj.N]);
    % RECOMPUTE full_H & theta_proj when the subspace is collapsed to avoid slow accumulation of numerical errors 
    obj.theta_proj = (obj.P.') * (obj.theta); % theta projected into current working subspace;
    obj.full_H = real(obj.reshape_wrapper(obj.b,[ss,-1]) * obj.reshape_wrapper(obj.b,[ss,-1]).');% full approximate hessian;
end