Contents
% 2.3.1. Expanding the subspace with a new observation: % The function below is written to update the subspace at each optimization step t to include the most recent % gradient direction as an additional column. This is done by finding the gradient vector component that lies % outside the exsiting subspace and appending it to the current subspace. 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; % 2.3.2. 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 % P^{'}=orth \left(\left[f_{1}^{'}(x^{\tau_{1}^{t}}) \ldots f_{N}^{'}(x^{\tau_{N}^{t}}) % \quad x^{\tau_{1}^{t}} \ldots x^{\tau_{N}^{t}}\right]\right ) Equation (9) if obj.K_current >= obj.K_max % collapse the subspace when it exceeded its maximum allowed size % xl may not be in the history yet, so we pass it in explicitly to make; % sure it's used; xl = (obj.P.') * (x_in); obj.collapse_subspace(xl); % use the function collapse_subspace below end end 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); % TODO(jascha) this may be a stricter cutoff than we need 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 to subspace projection matrix.; % T_right - The contravariant subspace projection matrix.; % (Currently T_left = T_right since the subspace is orthogonal. Will change if eg the code is adapted to also; % incorporate a 'natural gradient' based parameter space transformation.); [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
Not enough input arguments.
Error in subspace_construction (line 7)
if obj.K_current >= obj.M