function update_hessian(obj, indx)
%% BFGS initialization:           
    gd = find(sum(obj.hist_deltatheta(:,:,indx).^2, 1)>0);% test if there is any history for this subfunction
    num_gd = length(gd);
    %% CASE(I): NO sufficient history to run BFGS
    if num_gd == 0 % if no history, initialize with the median eigenvalue from full Hessian of other subfunctions;
        obj.b(:,:,indx) = 0.;
        H = obj.get_full_H_with_diagonal();
        [U, ~] = obj.eigh_wrapper(H);
        obj.min_eig_sub(indx) = median(U)/sum(obj.active); % median eigenvalue of avg Hessian of other active subfunctions
        obj.max_eig_sub(indx) = obj.min_eig_sub(indx);
        if obj.eval_count(indx) > 2
            if sum(obj.eval_count) < 5 % Subfunction was evaluated multiple times, but has no stored history.
                obj.min_eig_sub(indx) = obj.min_eig_sub(indx) / 10.; % Need to initialize SFO with a smaller hessian_init value >> scale down the hessian_init value
            end
        end
        return;
    end
    %% CASE(II): enough history to run BFGS
    % Working in the subspace defined by this subfunction's history (P_hist);
    [P_hist, ~] = qr([obj.hist_deltatheta(:,gd,indx),obj.hist_deltadf(:,gd,indx)], 0);
    deltatheta_P = P_hist.' * obj.hist_deltatheta(:,gd,indx);
    deltadf_P = P_hist.' * obj.hist_deltadf(:,gd,indx);
    %% approximate the smallest eigenvalue \beta by solving the squared secant equation for full history;
    %% \beta will be used as diagonal initialization for BFGS, (i.e. B_0=\beta I);
    % calculate Hessian using pinv & squared equation:   
    % \Delta f_{s}^{'}=B_s \Delta x_s (secant equation >>(df = H dx))
    % (df^T df = dx^T H^T H dx = dx^T H^2 dx)
    % Q=\left [(\Delta x)^{+^{T}} (\Delta f^{'})^{T} \Delta f^{'} (\Delta x)^{+}\right ]^{1/2} Equation (12)
    pdelthet = pinv(deltatheta_P);
    dd = (deltadf_P) * (pdelthet);
    H2 = (dd.') * (dd);
    [H2w, ~] = obj.eigh_wrapper(H2);
    H2w = sqrt(abs(H2w));
    H2w = sort(H2w, 'descend'); % only the top (~num_gd) eigenvalues are expected to be well defined;
    H2w = H2w(1:num_gd);

% Enforcing positive definiteness for the estimated Hessian: Each estimated $ {H_{i}}^t$ is constrained to be +ve definite by explicit eigendecomposition and seeting any
% too-small eigen values to the median positive eigenvalue.
    if min(H2w) == 0 || sum(~isfinite(H2w)) > 0 % ~zero eigenvalue or close (poor estimate of curvature)
    % a failure using this history due to degenerate deltadf (0 case) or deltatheta (non-finite case).;
        H2w(:) = max(obj.min_eig_sub(obj.active)); % Initialize using other subfunctions;
    end
    obj.min_eig_sub(indx) = min(H2w);
    obj.max_eig_sub(indx) = max(H2w);
    % constraining Hessian initialization
    if obj.min_eig_sub(indx) < obj.max_eig_sub(indx)/obj.hess_max_dev % eigenvalues less than (\gamma \lambda_{max})
        obj.min_eig_sub(indx) = obj.max_eig_sub(indx)/obj.hess_max_dev; % constrain using the allowed ratio;
    end

    %% Recalculate Hessian;
    num_hist = size(deltatheta_P, 2); % number of history terms (columns);
    % the new hessian will be (b_p)*(b_p.') + eye().*obj.min_eig_sub(indx);
    b_p = zeros(size(P_hist,2), num_hist*2);
    % step through the history;
    for hist_i = num_hist:-1:1
        s = deltatheta_P(:,hist_i);
        y = deltadf_P(:,hist_i);
        % for numerical stability;
        rscl = sqrt(sum(s.^2));
        s = s/rscl;
        y = y/rscl;
        % the BFGS step proper
        Hs = s.*obj.min_eig_sub(indx) + b_p*((b_p.')*(s));
        term1 = y/sqrt(sum(y.*s));
        sHs = sum(s.*Hs);
        term2 = sqrt(complex(-1.)).*Hs/sqrt(sHs);
        if sum(~isfinite(term1)) > 0 || sum(~isfinite(term2)) > 0
            obj.min_eig_sub(indx) = max(H2w); % invalid bfgs history term
            continue;
        end
        b_p(:,2*(hist_i-1)+2) = term1;
        b_p(:,2*(hist_i-1)+1) = term2;
    end

    H = real((b_p) * (b_p.')) + eye(size(b_p, 1))*obj.min_eig_sub(indx);
    % constrain the Hessian to be positive definite;
    [U, V] = obj.eigh_wrapper(H);
    if max(U) <= 0. %if there aren't any positive eigenvalues after BFGS
        U(:) = obj.max_eig_sub(indx); % set them all to be the same conservative diagonal value;
    end
    % set any too-small eigenvalues to the median positive eigenvalue;
    U_median = median(U(U>0)); 
    U(U<(max(abs(U))/obj.hess_max_dev)) = U_median;
    % Hessian after it's been forced to be positive definite;
    H_posdef = bsxfun(@times, V, U') * V.';
    % break it apart into matrices b & a diagonal term again;
    B_pos = H_posdef - eye(size(b_p, 1))*obj.min_eig_sub(indx);
    [U, V] = obj.eigh_wrapper(B_pos);
    b_p = bsxfun(@times, V, sqrt(obj.reshape_wrapper(U, [1,-1])));
    obj.b(:,:,indx) = 0.;
    obj.b(:,1:size(b_p,2),indx) = (P_hist)*(b_p);
end