|
- function [Z,theta,itr] = naml(Kmix,K,lambda)
-
- Kint = zeros(size(Kmix{1}));
- for i=1:1:size(Kmix,2)
- Kint = Kint + Kmix{i};
- end
- Kint = Kint / size(Kmix,2);
- %%Kint = center_kernel(Kint);
- Kint = kcenter(Kint);
-
-
- [Z,error] = kernel_kmeans(Kint,K,50);
-
- tol = 0.01;
-
- olde = -1;
- maxe = 0;
- convergence = 0;
- [N,D]=size(Kint);
-
-
- L = Z*((Z'*Z)^-0.5);
-
- itr =1;
- Lnew = L;
-
- while (~convergence)
- disp '***********************************';
- % obtaining Q given L and G
- Sk1= Kint *L*L'*Kint;
- [L Lnew]
- q = rank(Sk1);
- Sk2= Kint*Kint + lambda*Kint;
-
- pinvSk2 = pinv(Sk2'*Sk2)*Sk2';
- Qsolve = pinvSk2*Sk1;
- [Sq,Dq,Vq] = svd(Qsolve);
- Q = Sq(:,1:q);
-
- % obtaining G given Q and L
- [obj,B,t,w,stq] =qoqc_LSSVM_multiclass_correct(Kmix,L,lambda)
- %[w,B,S,time]=silp5(Kmix,L,lambda);
- %[w,Knew,weights] =qoqc2(Kmix,L,lambda,theta);
-
- Knew = zeros(size(Kmix{1}));
- for loop=1:1:size(Kmix,2)
- Knew =Knew+ Kmix{loop}*w(loop);
- end
- Kint = Knew;
-
- Lcond = Q'*(Knew*Knew+lambda*Knew)*Q;
-
-
- % cond(Knew*Knew+lambda*Knew);
- % if abs(cond(Knew*Knew+lambda*Knew)-1)<tol
- % G = Knew *inv(Knew*Knew+lambda.*Knew)*Knew;
- % else
- % G = Knew *pinv(Knew*Knew+lambda.*Knew)*Knew;
- % end
- % [Sw,Dw,Vw]=svd(G);
- % Sw = Sw(:,1:K);
-
- % Lnew = Sw;
-
- %cond(Knew*Knew+lambda*Knew)
- if abs(cond(Lcond)-1)<0.1
- G = Knew*Q*inv(Lcond)*Q'*Knew;
- else
- G = Knew*Q*pinv(Lcond)*Q'*Knew;
- end
- [Sw,Dw,Vw]=svd(G);
-
- diag(Dw)
- Sw = Sw(:,1:K);
-
-
- Lnew = Sw;
-
- e = trace(Lnew'*G*Lnew);
- %e = trace(L'*G*L);
- maxe = max(e,maxe);
-
- theta = w;
- if (abs(olde -e)<tol && abs(maxe-e)<tol)
- convergence = 1;
-
- [Q,R]= qr(Lnew');
- R11 = R(1:K,1:K);
- R11 = R11 + 0.00001*eye(K,K);
- Rnew = inv(R11) * R;
- [rmax,IDX]= max(Rnew);
- IDX = IDX';
- Znew = zeros(size(Z));
- for i=1:1:size(Znew,2)
- Znew(find(IDX==i),i)=1;
- end
- Zspec = Znew;
- return;
-
- else
- olde = e;
- fprintf('Error = %f\n', e);
- L = Lnew;
- itr = itr+1;
- end
- end
-
-
- function [obj,B,t,mu,stq] =qoqc_LSSVM_multiclass_correct(Kmix,L,lambda)
-
- % read how many kernels
- numP = size(Kmix,2);
-
- % read how many clusters
- numK = size(L,2);
-
- % read how many data
- numD = size(Kmix{1},1);
-
- % creat cluster indicator matrix L from partition matrix P
- %L = P*(P'*P)^-0.5;
-
-
-
- %variables
- % B_1 B_2 ,..., B_K, t
-
-
- % init the mosek problem
- clear prob;
-
-
- % c vector .
- c = -ones(numD*numK+1,1);
- c(end)=0.5;
- prob.c =c;
- clear c;
-
-
- % the quadratic term in objective function
-
- prob.qosubi = [];
- prob.qosubj = [];
- prob.qoval = [];
-
-
- for kernelLoop = 1:1:numK
-
- prob.qosubi = [prob.qosubi numD*(kernelLoop-1)+1:numD*(kernelLoop-1)+numD];
- prob.qosubj = [prob.qosubj numD*(kernelLoop-1)+1:numD*(kernelLoop-1)+numD];
- prob.qoval = [prob.qoval 1/lambda*ones(1,numD)];
- end
- prob.qosubi = [prob.qosubi numK*numD+1]';
- prob.qosubj = [prob.qosubj numK*numD+1]';
- prob.qoval = [prob.qoval 0]';
-
-
- % Next quadratic terms in the constraints .
- % constraint of 1st kernel
- prob.qcsubk =[];
- prob.qcsubi =[];
- prob.qcsubj = [];
- prob.qcval =[];
- prob.a = [];
- prob.buc = [];
- prob.blc=[];
- prob.blx = -inf*ones(numD*numK+1,1); % the dual variables are unconstrained, the dummy variable is also unconstrained (actually >=0)
- prob.bux = +inf*ones(numD*numK,1); % the dual variables are unconstrained, the dummy variable is also unconstrained (actually >=0)
- prob.bux = [prob.bux; +inf]; % the dual variables are unconstrained, the dummy variable is also unconstrained (actually >=0)
-
- yacc =[];
- for clusterLoop = 1:1:numK
- yacc= [yacc; L(:,clusterLoop)];
- end
- yacc =[yacc; 0];
-
-
- for kernelLoop = 1:1:numP
- Qk=[];
-
- for clusterLoop = 1:1:numK
- clusterLoop
- Y = diag(L(:,clusterLoop));
- Qk=[Qk;sparse(numD,(clusterLoop-1)*numD) 2*Y*Kmix{kernelLoop}*Y sparse(numD,(numK-clusterLoop)*numD+1)];
- end
-
- Qk=[Qk; sparse(1,numD*numK+1)];
-
- trQk = tril(Qk);
- trQk = sparse(trQk);
- [i,j,s] = find(trQk);
-
- clear Qk;
-
- prob.qcsubk = [prob.qcsubk; kernelLoop*ones(size(i))];
- prob.qcsubi = [prob.qcsubi i'];
- prob.qcsubj = [prob.qcsubj j'];
- prob.qcval = [prob.qcval s'];
-
- prob.a = [prob.a; sparse(1,numK*numD) -1]; % a'*Y*K*Y*a - 2*gamma <=0
- prob.buc = [prob.buc; 0];
- prob.blc = [prob.blc; -inf];
-
-
- end
-
-
- % linear constraints (Ya)'*\vec 1 = 0
-
- size(yacc')
- size(prob.a)
- prob.a = [prob.a; yacc'];
-
- prob.buc = [prob.buc; 0];
- prob.blc = [prob.blc; 0];
-
- prob.qcsubi = prob.qcsubi';
- prob.qcsubj = prob.qcsubj';
- prob.qcval = prob.qcval';
-
- size(prob.c)
- size(prob.qcsubi)
-
- tic;
- [r,res] = mosekopt ('minimize', prob);
- stq =toc;
-
- % get the optimal dual variable solution
-
- size(res.sol.itr.xx)
-
- beta = res.sol.itr.xx(1:numK*numD);
-
-
- B = [];
-
- for loop=1:1:numK
- B = [B beta((loop-1)*numD+1:loop*numD)];
- end
-
- t = res.sol.itr.xx(end);
-
- mu = res.sol.itr.suc(1:numP);
- mu = mu.*2;
-
- obj = 0.5 * t + 0.5/lambda* (beta'*beta)- sum(beta);
-
- function [Z,Y] = kernel_kmeans(Kernel,K,Max_Its,Zold)
-
- %This is a simple implementation of Kernel K-means clustering - an
- %interesting paper which proposed kernel based Kmeans clustering is [1]
- %Girolami, M, Mercer Kernel-Based Clustering in Feature Space,
- %IEEE Trans Neural Networks, 13(3),780 - 784, 2002.
- %INPUTS:
- % Kernel = the kernel matrix for clustering
- % K = the number of clusters
- % Max_Its = maximal number of iterations
- % Zold = an optional parameter to use the Zold as a starting point
- %
- %OUTPUTS:
- % Z = the NxK binary cluster assignment matrix
- % Y = the matrix of disimilarities
-
-
- [N,D]=size(Kernel);
-
-
- if nargin < 4
- %initialise the indictaor matrix to a random segmentation of the data
- Z = zeros(N,K);
- for n = 1:N
- Z(n,rand_int(K)) = 1;
- end
- else
- Z = Zold;
- end
-
- for k=1:1:K
- if (sum(Z(:,k))==0)
- c = rand_int(N);
- Z(c,:) = zeros(1,K);
- Z(c,k)==1;
- end
- end
-
- %main loop
- for its = 1:Max_Its
- %compute the similarity of each data point to each cluster mean in
- %feature space - note we do not need to compute store or update a mean
- %vector s we are using the kernel-trick - cool eh?
-
- for k=1:K
- Nk = sum(Z(:,k));
- Y(:,k) = diag(Kernel) - 2*Kernel*Z(:,k)./Nk + Z(:,k)'*Kernel*Z(:,k)./(Nk^2);
- end
-
- %Now we find the cluster assignment for each point based on the minimum
- %distance of the point from the mean centres in feature space using the
- %Y matrix of dissimilarities
- %for yLoop = 1:N
- % Y(yLoop,:)=Y(yLoop,:)./max(Y(yLoop,:));
- %end
- %Y
- [i,j]=min(Y,[],2);
-
- %this simply updates the indictor matrix Z refleting the new
- %allocations of data points to clusters
- Z = zeros(N,K);
- for n=1:N
- Z(n,j(n)) = 1;
- end
-
- end
-
-
- %return the clutsers that each data point has been allocated to
- for n=1:N
- z(n) = find(Z(n,:));
- end
-
- function u = rand_int(Max_Int)
- u=ceil(Max_Int*rand);
|