function optimization_step(obj)
% Perform a single optimization step. 
    time_pass_start = tic();
    indx = obj.get_target_index(); %% choose an index to update;
    [f, df_proj] = obj.f_df_wrapper(obj.theta, indx); % evaluate subfunction value & gradient at new position;
    [step_failure, f, df_proj] = obj.handle_step_failure(f, df_proj, indx); % check for a failed update step, & adjust f, df, & obj.theta;
    % as appropriate if one occurs.;
    obj.update_history(indx, obj.theta_proj, f, df_proj); % add the change in theta & the change in gradient to the history for this subfunction;
    obj.total_distance = obj.total_distance + sqrt(sum((obj.theta - obj.theta_prior_step).^2)); % increment the total distance traveled using the last update;
    H_pre_update = real(obj.b(:,:,indx) * obj.b(:,:,indx).'); % the current contribution from this subfunction to the total Hessian approximation;
    obj.update_hessian(indx);%% update this subfunction's Hessian estimate;
    H_new = real(obj.b(:,:,indx) * obj.b(:,:,indx).');% the new contribution from this subfunction to the total approximate hessian;
    obj.full_H = obj.full_H + H_new - H_pre_update;% update total Hessian using this subfunction's updated contribution;

    % calculate the total gradient, total Hessian, & total function value at the current location;
    full_df = 0.;
    for i = 1:obj.N
        dtheta = obj.theta_proj - obj.last_theta(:,i);
        bdtheta = obj.b(:,:,i).' * dtheta;
        Hdtheta = real(obj.b(:,:,i) * bdtheta);
        Hdtheta = Hdtheta + dtheta.*obj.min_eig_sub(i); % the diagonal contribution;
        full_df = full_df + Hdtheta + obj.last_df(:,i);
    end
    full_H_combined = obj.get_full_H_with_diagonal();
    full_H_inv = inv(full_H_combined);
    % calculate an update step;
    dtheta_proj = -(full_H_inv) * (full_df) .* obj.step_scale;% calculate the Newton step described in Equation (5)
    dtheta_proj_length = sqrt(sum(dtheta_proj(:).^2));
    if dtheta_proj_length < obj.minimum_step_length
        dtheta_proj = dtheta_proj*obj.minimum_step_length/dtheta_proj_length;% forcing the update to have minimum step length to avoid numerical errors
        dtheta_proj_length = obj.minimum_step_length;% assign the minimum step length to the update step length
    end
    if sum(obj.eval_count) > obj.N && dtheta_proj_length > obj.eps
    % only allow a step to be up to a factor of max_step_length_ratio longer than the average step length
        avg_length = obj.total_distance / sum(obj.eval_count);
        length_ratio = dtheta_proj_length / avg_length;
        ratio_scale = obj.max_step_length_ratio;
        if length_ratio > ratio_scale
            % truncating step length from (dtheta_proj_length) to (ratio_scale*avg_length) 
            dtheta_proj_length = dtheta_proj_length/(length_ratio/ratio_scale);
            dtheta_proj = dtheta_proj/(length_ratio/ratio_scale);
        end
    end

    dtheta = (obj.P) * (dtheta_proj);% the update to theta, in the full dimensional space;
    obj.theta_prior_step = obj.theta;% backup the prior position, in case this is a failed step;
    obj.theta = obj.theta + dtheta;% update theta to the new location;
    obj.theta_proj = obj.theta_proj + dtheta_proj;
    % the predicted improvement from this update step;
    obj.f_predicted_total_improvement = 0.5 .* dtheta_proj.' * (full_H_combined * dtheta_proj);
    obj.expand_active_subfunctions(full_H_inv, step_failure); %% expand the set of active subfunctions as appropriate;
    % record how much time was taken by this learning step;
    time_diff = toc(time_pass_start);
    obj.time_pass = obj.time_pass + time_diff;
end