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
