基于PCA算法的人脸识别（包含人脸库）

%% This function is modified from Matlab Package: L1-Homotopy % BPDN_homotopy_function.m % % Solves the following basis pursuit denoising (BPDN) problem % min_x \lambda ||x||_1 + 1/2*||y-Ax||_2^2 % % Inputs: % A - m x n measurement matrix % y - measurement vector % lambda - final value of regularization parameter % maxiter - maximum number of homotopy iterations % % Outputs: % x_out - output for BPDN % total_iter - number of homotopy iterations taken by the solver % % Written by: Salman Asif, Georgia Tech % Email: sasif@ece.gatech.edu % %-------------------------------------------+ % Copyright (c) 2007. Muhammad Salman Asif %-------------------------------------------+ function [x_out, e_out, total_iter] = SolveHomotopy_CBM_std(A, y, varargin) global N n gamma_x z_x xk_temp del_x_vec pk_temp dk epsilon isNonnegative t0 = tic ; lambda = 1e-6; maxiter = 100; isNonnegative = false; verbose = false; xk_1 = []; tolerance = 1e-4 ; STOPPING_TIME = -2; STOPPING_GROUND_TRUTH = -1; STOPPING_DUALITY_GAP = 1; STOPPING_SPARSE_SUPPORT = 2; STOPPING_OBJECTIVE_VALUE = 3; STOPPING_SUBGRADIENT = 4; STOPPING_DEFAULT = STOPPING_OBJECTIVE_VALUE; stoppingCriterion = STOPPING_DEFAULT; % Parse the optional inputs. if (mod(length(varargin), 2) ~= 0 ), error(['Extra Parameters passed to the function ''' mfilename ''' must be passed in pairs.']); end parameterCount = length(varargin)/2; for parameterIndex = 1:parameterCount, parameterName = varargin{parameterIndex*2 - 1}; parameterValue = varargin{parameterIndex*2}; switch lower(parameterName) case 'stoppingcriterion' stoppingCriterion = parameterValue; case 'initialization' xk_1 = parameterValue; if ~all(size(xk_1)==[n,1]) error('The dimension of the initial x0 does not match.'); end case 'groundtruth' xG = parameterValue; case 'lambda' lambda = parameterValue; case 'maxiteration' maxiter = parameterValue; case 'isnonnegative' isNonnegative = parameterValue; case 'tolerance' tolerance = parameterValue; case 'verbose' verbose = parameterValue; case 'maxtime' maxTime = parameterValue; otherwise error(['The parameter ''' parameterName ''' is not recognized by the function ''' mfilename '''.']); end end clear varargin [K, n] = size(A); B = [A eye(K)]; At = A'; Bt = B'; N = K + n; timeSteps = nan(1,maxiter) ; % Initialization of primal and dual sign and support z_x = zeros(N,1); gamma_x = []; % Primal support % Initial step Primal_constrk = -[At*y; y]; if isNonnegative [c i] = min(Primal_constrk); c = max(-c, 0); else [c i] = max(abs(Primal_constrk)); end epsilon = c; nz_x = zeros(N,1); if isempty(xk_1) xk_1 = zeros(N,1); gamma_xk = i; else gamma_xk = find(abs(xk_1)>eps*10); nz_x(gamma_xk) = 1; end f = epsilon*norm(xk_1,1) + 1/2*norm(y - A*xk_1(1:n) - xk_1(n+1:end))^2; z_x(gamma_xk) = -sign(Primal_constrk(gamma_xk)); %Primal_constrk(gamma_xk) = sign(Primal_constrk(gamma_xk))*epsilon; z_xk = z_x; % loop parameters iter = 0; out_x = []; old_delta = 0; count_delta_stop = 0; gamma_xkx = gamma_xk(gamma_xk<=n); gamma_xke = gamma_xk(gamma_xk>n)-n; AtgxAgx = [At(gamma_xkx,:)*A(:,gamma_xkx) At(gamma_xkx,gamma_xke); A(gamma_xke,gamma_xkx) eye(length(gamma_xke))]; iAtgxAgx = inv(AtgxAgx); while iter < maxiter iter = iter+1; gamma_x = gamma_xk; gamma_xx = gamma_x(gamma_x<=n); gamma_xe = gamma_x(gamma_x>n)-n; z_x = z_xk; x_k = xk_1; %%%%%%%%%%%%%%%%%%%%% %%%% update on x %%%% %%%%%%%%%%%%%%%%%%%%% % Update direction del_x = iAtgxAgx*z_x(gamma_x); del_x_vec = zeros(N,1); del_x_vec(gamma_x) = del_x; Asupported = B(:,gamma_x); Agdelx = Asupported*del_x; dk = [At*Agdelx; Agdelx]; %%% CONTROL THE MACHINE PRECISION ERROR AT EVERY OPERATION: LIKE BELOW. pk_temp = Primal_constrk; gammaL_temp = find(abs(abs(Primal_constrk)-epsilon)<min(epsilon,2*eps)); pk_temp(gammaL_temp) = sign(Primal_constrk(gammaL_temp))*epsilon; xk_temp = x_k; xk_temp(abs(x_k)<2*eps) = 0; %%%--- % Compute the step size [i_delta, delta, out_x] = update_primal(out_x); if old_delta < 4*eps && delta < 4*eps count_delta_stop = count_delta_stop + 1; if count_delta_stop >= 500 if verbose disp('stuck in some corner'); end break; end else count_delta_stop = 0; end old_delta = delta; xk_1 = x_k+delta*del_x_vec; Primal_constrk = Primal_constrk+delta*dk; epsilon_old = epsilon; epsilon = epsilon-delta; if epsilon <= lambda; % xk_1 = x_k + (epsilon_old-lambda)*del_x_vec; % disp('Reach prescribed lambda in SolveHomotopy_CBM.'); break; end timeSteps(iter) = toc(t0) ; % if mod(iter, 100) == 0 % disp([ 'epsilon = ' num2str(epsilon) ', time =' num2str(timeSteps(iter))]); % end % compute stopping criteria and test for termination keep_going = true; if delta~=0 switch stoppingCriterion case STOPPING_GROUND_TRUTH keep_going = norm(xk_1(1:n)-xG)>tolerance; case STOPPING_SPARSE_SUPPORT nz_x_prev = nz_x; nz_x = (abs(xk_1)>eps*10); num_nz_x = sum(nz_x(:)); num_changes_active = (sum(nz_x(:)~=nz_x_prev(:))); if num_nz_x >= 1 criterionActiveSet = num_changes_active / num_nz_x; keep_going = (criterionActiveSet > tolerance); end case STOPPING_DUALITY_GAP error('Duality gap is not a valid stopping criterion for Homotopy.'); case STOPPING_OBJECTIVE_VALUE % continue if not yeat reached target value tolA prev_f = f; f = lambda*norm(xk_1,1) + 1/2*norm(y-Asupported*xk_1(gamma_x))^2; keep_going = (abs((prev_f-f)/prev_f) > tolerance); case STOPPING_SUBGRADIENT keep_going = norm(delta*del_x_vec)>tolerance; case STOPPING_TIME keep_going = timeSteps(iter) < maxTime ; otherwise, error('Undefined stopping criterion'); end % end of the stopping criteria switch end if ~keep_going %disp('Maximum time reached!'); break; end if ~isempty(out_x) % If an element is removed from gamma_x len_gamma = length(gamma_x); outx_index = find(gamma_x==out_x(1)); gamma_x(outx_index) = gamma_x(end); gamma_x = gamma_x(1:end-1); gamma_xk = gamma_x; rowi = outx_index; % ith row of A is swapped with last row (out_x) colj = outx_index; % jth column of A is swapped with last column (out_lambda) AtgxAgx_ij = AtgxAgx; temp_row = AtgxAgx_ij(rowi,:); AtgxAgx_ij(rowi,:) = AtgxAgx_ij(len_gamma,:); AtgxAgx_ij(len_gamma,:) = temp_row; temp_col = AtgxAgx_ij(:,colj); AtgxAgx_ij(:,colj) = AtgxAgx_ij(:,len_gamma); AtgxAgx_ij(:,len_gamma) = temp_col; iAtgxAgx_ij = iAtgxAgx; temp_row = iAtgxAgx_ij(colj,:); iAtgxAgx_ij(colj,:) = iAtgxAgx_ij(len_gamma,:); iAtgxAgx_ij(len_gamma,:) = temp_row; temp_col = iAtgxAgx_ij(:,rowi); iAtgxAgx_ij(:,rowi) = iAtgxAgx_ij(:,len_gamma); iAtgxAgx_ij(:,len_gamma) = temp_col; AtgxAgx = AtgxAgx_ij(1:l