分支定界算法的步骤：
1 （初始化）： B= ∞，L=0（当前水平） ,p=0（当前结点）
2（当前结点展开）：将当前结点的所有直接后继结点放入一个目录表（活动表）中，对它们计算并存储D(x, M_p)
3 （规则1检验）：对活动表中的每个结点，若D(x, M_p)>B+r_p，则将该结点从活动表中去掉；
4（回溯）：如果活动表为空（当前结点的孩子结点被全部剪掉），则回溯到上一级，即L=L-1；如果L=0，那么输出X_NN和D(X_NN, x)=B，否则至5 ；
5 （选择最近结点）：在活动表中选择使D(x, M_p)最小的结点p’作为当前结点，若L为最终水平，则转6，否则，置L＝L＋1，转2；
6 （规则2检验）：对当前结点p’中的每个点x_i，检验如果D(x, M_p)>B+D(x_i, M_p)，那么x_i为非最近邻；否则计算D(x_i, x)，如果D(x_i, x)< B，那么NN=i，B=D(x_i, x)

matlab实现代码如下：

X = [randn(200,2)+ones(200,2);randn(200,2)-1.5*ones(200,2);randn(200,2)+2.5*ones(200,2);];
[row,col]=size(X);
all_idx=0;L=3;l=3;   % 每个水平均划分为3个子集
%计算总节点的数目
for i=1:L
    all_idx=all_idx+l^i;
end
Xp=cell(all_idx,1);
Mp=zeros(all_idx,col);
Rp=zeros(all_idx,1);
p=1;
for i=1:L
   if i==1
       [IDX,C,sumd,D] = kmeans(X,l);
        for j=1:l
            Xp(p)={X((IDX==j),:)};
            Mp(p,:)=C(j,:);
            Rp(p)=max(D((IDX==j),j));
            p=p+1;
        end  
   else
       endk=p-1;begink=endk-l^(i-1)+1;      
       for k=begink:endk          
           [IDX,C,sumd,D] = kmeans(Xp{k,1},l);
           X1=Xp{k,1};
           for j=1:l               
                Xp(p)={X1((IDX==j),:)};
                Mp(p,:)=C(j,:);
                Rp(p)=max(D((IDX==j),j));
                p=p+1;
           end  
       end
   end
end
x=randn(1,2);       %待判样本
B=inf;CurL=1;p=0;TT=1;
while TT==1         %步骤2
    Xcurp=cell(1);
    CurTable=cell(l,1);
    CurPinT=zeros(l,1);
    Dx=zeros(l,1);
    RpCur=zeros(l,1);
    %当前节点的直接后继放入目录表  
    for i=1:l  
        CurTable(i,1)=Xp(i+p*l,1);
        CurPinT(i)=i+p*l;
        Dx(i)=norm(x-Mp(i+p*l,:))^2;
        RpCur(i)=Rp(i+p*l);
    end   
    while 1         %步骤3
        [rowT,colT]=size(CurTable);
        for i=1:rowT                  
            if Dx(i)>B+RpCur(i)+eps%从目录表中去掉当前节点p
                CurTable(i,:)=[];
                CurPinT(i)=[];
                Dx(i)=[];
                RpCur(i)=[];
                break;
            end
        end
        [CurRowT,CurColT]=size(CurTable);
        if CurRowT==0
           CurL=CurL-1;p=floor((p-1)/3);
           if CurL==0
              TT=0; break; 
           else
               %转步骤3
           end
        elseif CurRowT>0
            [Dxx,Dxind]=sort(Dx,'ascend');
            p1=CurPinT(Dxind(1));
            p=p1;

            for j=1:CurRowT      %从当前目录表去掉p1
                if CurPinT(j)==p1
                    Xcurp(1,1)=CurTable(j,1);
                    CurTable(j,:)=[];
                    CurPinT(j)=[];
                    CurD=Dx(j);  %记录D(x,Mp)
                    Dx(j)=[];
                    RpCur(j)=[];                   
                    break;
                end
            end
            if CurL==L
                XcurpMat=cell2mat(Xcurp);
                [CurpRow,CurpCol]=size(XcurpMat);
                CurpMean=Mp(p,:);
                for k=1:CurpRow
                    Dxi=norm((XcurpMat(k,:)-CurpMean))^2;
                    if CurD>Dxi+B+eps

                    else
                        Dxxi=norm((x-XcurpMat(k,:)))^2;
                        if Dxxi<B+eps
                            B=Dxxi;Xnn=XcurpMat(k,:);
                        end
                    end
                end
            else
                CurL=CurL+1;
                break;
            end
        end
    end
end
B,Xnn,NN=find(X(:,1)==Xnn(1))
figure, plot(X(1:200,1),X(1:200,2),'m.');
hold on,plot(X(201:400,1),X(201:400,2),'b.');
hold on,plot(X(401:600,1),X(401:600,2),'g.');  
hold on,plot(Xnn(1,1),Xnn(1,2),'kx ','MarkerSize',10,'LineWidth',2);
hold on,plot(x(1),x(2),'r+','MarkerSize',10,'LineWidth',2);
legend('Cluster 1','Cluster 2','Cluster 3','NN','x','Location','NW');