基于半监督支持向量机算法，能得到全局最优解，很好用

function [Y,alpha,b,obj] = tsvm_bb(Y,params) % Unlabeled points have Y=0 global K A cnt % K = entire kernel matrix cnt = 0; % cnt = count the number of recursive calls lab = find(Y~=0); unlab = find(Y==0); % Compute the objective function corresponding to % the "true moon". Only for debugging. % Y2 = Y; Y2(3:1+params.np)=-1; Y2(2+params.np:end)=1; % obj0 = obj_fun(Y2,params); % Set the upper bound to some solution. % Does not really matter: could be skipped % Or can use any other solution Y0 = Y; np = params.np - sum(Y==1); invK = inv([[K+eye(size(K))/params.C ones(length(K),1)]; ... [ones(1,length(K)) 0]]); out = -invK(unlab,unlab) \ (invK(unlab,lab)*Y(lab)); [foo,ind] = sort(-out); Y0(unlab(ind(1:np))) = 1; Y0(unlab(ind(np+1:end))) = -1; ub = obj_fun(Y0,params); set(0,'RecursionLimit',10+length(Y)); % Let's start ! [obj, Y] = bb(Y,params,ub); % Compute the alpha corresponding to the solution [obj, alpha, b] = obj_fun(Y,params); function [obj,Yb] = bb(Y,params,ub) % ub is the current upper bound % returns: obj is smallest objective function in the subtree % and Yb the corresponding solution global cnt cnt = cnt + 1; Yb = Y; % The balancing constraint isn't (or won't) be satisfied -> exit if (sum(Y==1) > params.np) | (sum(Y==-1) > length(Y)-params.np) obj = +Inf; return; end; % Train with current labeled set. This gives a lower bound on % the objective value in this subtree. [obj, alpha, b, next] = obj_fun(Y,params); % fprintf('%d %d %f %f\r',sum(Y~=0),sum(alpha~=0),obj,ub); % If we are larger than the current best solution -> exit if obj > ub return; end; % Leaf -> exit if all(Y~=0) % fprintf(' %d %f\n',cnt,obj); return; end; % Go down in the tree. Start by the point for which we are the % most confident about the label (such that the upper bound get % reduced quickly) Y2 = Y; Y2(next.ind) = next.sign; params.alpha = zeros(length(Y),1); params.alpha(Y~=0) = alpha; params.b = b; [obj, Yb] = bb(Y2,params,ub); % Try the other sign (hopefully this part of the tree will not % be explored too much Y2(next.ind) = -Y2(next.ind); [obj2, Yb2] = bb(Y2,params,min(obj,ub)); % Keep the best solution if obj2<obj Yb = Yb2; obj = obj2; end; function [obj, alpha, b, next] = obj_fun(Y,params) % Do the SVM training % Next contains the unlabeled point which that we should try to label next global K lab = find(Y~=0); unlab = find(Y==0); if isfield(params,'alpha') alpha = params.alpha; b = params.b; out = K(lab,alpha~=0)*alpha(alpha~=0)+b ; svset = find(Y(lab).*out <= 1); else svset = [1:length(lab)]'; end; old_svset = []; iter = 0; itermax = 20; % Primal SVM training (Newton steps) while iter<itermax if length(svset)==length(old_svset), if all(svset==old_svset) break; end; end; old_svset = svset; H = K(lab(svset),lab(svset)) + eye(length(svset))/params.C; H(end+1,:) = 1; % To take the bias into account H(:,end+1) = 1; H(end,end) = 0; % Find the parameters if var(Y(lab)) == 0 alpha = zeros(length(lab),1); b = Y(lab(1)); else par = H\[Y(lab(svset));0]; alpha = par(1:end-1); b = par(end); end; % Compute the new set of support vectors out = K(lab,lab(svset))*alpha+b ; svset = find(Y(lab).*out <= 1); iter = iter + 1; end; alphasv = alpha; alpha = zeros(length(lab),1); alpha(old_svset) = alphasv; if iter == itermax warning('Maximum number of Newton steps reached'); elseif any(alpha.*Y(lab)<0) warning('Bug in learn'); end; % Quick way to compute the objective function % (can be derived from KKT conditions) obj = sum(abs(alpha)); % Try to see which point should be labeled next. % Naive implementation = try each unlabeled point, with both signs and % take the largest objective function. % Approximate solution based on rank one update of the kernel matrix. % Would be exact if the set of SV does not change. % Other possibility would be to take the closet point from the labeled % set (in feature space). if isempty(unlab) next = []; return; end; K2 = [K(lab(svset),unlab); ones(1,length(unlab))]; out = K2'*[alphasv; b]; % Output of the unlabeled points span = 1./(diag(K(unlab,unlab))+1/params.C-sum(K2.*(inv(H)*K2),1)'); % The increase of the objective function is estimated as span*(1-out*y) [foo,ind] = max((1+abs(out)).^2.*span); % NB: max(1-out,1+out) = 1+abs(out) next.ind = unlab(ind); % We are almost sure about the label of this point % because if add it in with the wrong label, the objective function % will increase a lot. next.sign = sign(out(ind)); % We will start by the "correct" label