my_spec_cut.m

利用加速随机游走计算初始点确定的分类算法

function [cRe] = my_spec_cut(W,k,para)

%recursive do 2-way partition on W by some spectral method

if nargin<3
    para = [0.2 1];
end
p = para(1);%sparsification prob

flag = para(2);

[id1 id2] = my_spec_cut_once(W,p,flag);
cRe(id1,1) = 1;
cRe(id2,1) = 2;

% [Ad,Aod] = MatSpt(W,cRe);
% disp(num2str(nnz(Aod)/nnz(Ad)))

for i=2:k
    cRe = my_spec_cut_more(W,cRe,p,flag);
    
%     [Ad,Aod] = MatSpt(W,cRe);
%     disp(num2str(nnz(Aod)/nnz(Ad)))
end


cRe = cRe - 1;




function cRe1 = my_spec_cut_more(W,cRe0,p,flag)
%cut each partition in cRe0 into 2
cRe1 = cRe0;
n1 = max(cRe0);%# of partitions
for i=1:n1
    pos = find(cRe0==i);
    W0 = W(pos,pos);
    [id1 id2] = my_spec_cut_once(W0,p,flag);
    cRe1(pos(id2)) = i + n1;
end

function [id1 id2] = my_spec_cut_once(W,p,flag)
%cut W into two balance partitions

%by ncut
%[v,val] = ncut(W,2);
%by spectral clustering
[v,val] = my_spec_cut_eigcomp(W,2,p,flag);
[y,I] = sort(full(v(:,2)),1,'ascend');
n = length(I);
id1 = I(1:round(n/2));
id2 = I(1+round(n/2):end);

function [Eigenvectors,Eigenvalues] = my_spec_cut_eigcomp(W,nbEigenValues,p,flag)

%compute desired eigenvector
%flag 1: 2nd smallest of D-W
%flag 2: 2nd largest of S
%flag 3: 2nd largest of S * d^(-1/2)

if nargin < 4
    flag = 3;
end
if nargin < 3
    p = 0.1;%sparsification rate
end
if nargin < 2
    nbEigenValues = 2;
end
%if nargin < 4
    dataNcut.offset = 5e-1;
    dataNcut.verbose = 0;
    dataNcut.maxiterations = 100;
    dataNcut.eigsErrorTolerance = 1e-6;
    dataNcut.valeurMin=1e-6;
%end

% % make W matrix sparse
% W = sparsifyc(W,dataNcut.valeurMin);
% 
% % check for matrix symmetry
% if max(max(abs(W-W'))) > 1e-10 %voir (-12) 
%     %disp(max(max(abs(W-W'))));
%     error('W not symmetric');
% end


n = size(W,1);
nbEigenValues = min(nbEigenValues,n);
offset = dataNcut.offset;


% degrees and regularization
d = sum(abs(W),2);
dr = 0.5 * (d - sum(W,2));
d = d + offset * 2;
dr = dr + offset;
W = W + spdiags(dr,0,n,n);

Dinvsqrt = 1./sqrt(d+eps);
%P = spmtimesd(W,Dinvsqrt,Dinvsqrt);
if flag==2|flag==3
    P=BLin_W2P(W,3);%normalization by sqrt of row and col
%     tmp = [1:n]';
%     tmp(:,2) = tmp(:,1);
%     tmp(:,3) = 1;
%     P = spconvert(tmp) - P;
    %P = (P+P')/2;
else
    tmp = [1:n]';
    tmp(:,2) = tmp(:,1);
    tmp(:,3) = d;
    P = spconvert(tmp) - W;
    %P = (P+P')/2;
    clear tmp;    
end
clear W;
if p<1
    P = mysparsify(P,p);    
end
P = (P+P')/2;

options.issym = 1;
     
if dataNcut.verbose
    options.disp = 3; 
else
    options.disp = 0; 
end
options.maxit = dataNcut.maxiterations;
options.tol = dataNcut.eigsErrorTolerance;

options.v0 = ones(size(P,1),1);
options.p = max(35,2*nbEigenValues); %voir
options.p = min(options.p,n);

%warning off
if flag==2|flag==3
    [vbar,s,convergence] = eigs2(P,nbEigenValues,'LA',options); 
else
    [vbar,s,convergence] = eigs2(P,nbEigenValues,'SA',options); 
end
%warning on
if flag==1
    %s = real(diag(s));
    [x,y] = sort(-s); 
    %[x,y] = sort(s); 
    Eigenvalues = -x;
    vbar = vbar(:,y);
end
if flag==3|flag==2
    s = real(diag(s));
    [x,y] = sort(-s); 
    %[x,y] = sort(s); 
    Eigenvalues = -x;
    vbar = vbar(:,y);
    if flag==3
        Eigenvectors = spdiags(Dinvsqrt,0,n,n) * vbar;
    else
        Eigenvectors = vbar;
    end
%     for  i=1:size(Eigenvectors,2)
%         Eigenvectors(:,i) = (Eigenvectors(:,i) / norm(Eigenvectors(:,i))  )*norm(ones(n,1));
%         if Eigenvectors(1,i)~=0
%             Eigenvectors(:,i) = - Eigenvectors(:,i) * sign(Eigenvectors(1,i));
%         end
%     end
else
    Eigenvectors = vbar;
end
    
function W0 = mysparsify(W,p)

%sparsify matrix W
[a(:,1) a(:,2) a(:,3)] = find(W);
b = rand(size(a,1),1);
I = find(b>p);
a(I,3) = 0;
I = find(b<=p);
a(I,3) = a(I,3)/p;
W0 = spconvert(a);
[s0,t0] = size(W0);
[s,t] = size(W);
if s0<s|t0<t
    W0(s,t) = 0;
end



function  varargout = eigs2(varargin)
%
% Slightly modified version from matlab eigs, Timothee Cour, 2004
%
%   EIGS  Find a few eigenvalues and eigenvectors of a matrix using ARPACK.
%   D = EIGS(A) returns a vector of A's 6 largest magnitude eigenvalues.
%   A must be square and should be large and sparse.
%
%   [V,D] = EIGS(A) returns a diagonal matrix D of A's 6 largest magnitude
%   eigenvalues and a matrix V whose columns are the corresponding eigenvectors.
%
%   [V,D,FLAG] = EIGS(A) also returns a convergence flag.  If FLAG is 0
%   then all the eigenvalues converged; otherwise not all converged.
%
%   EIGS(A,B) solves the generalized eigenvalue problem A*V == B*V*D.  B must
%   be symmetric (or Hermitian) positive definite and the same size as A.
%   EIGS(A,[],...) indicates the standard eigenvalue problem A*V == V*D.
%
%   EIGS(A,K) and EIGS(A,B,K) return the K largest magnitude eigenvalues.
%
%   EIGS(A,K,SIGMA) and EIGS(A,B,K,SIGMA) return K eigenvalues based on SIGMA:
%      'LM' or 'SM' - Largest or Smallest Magnitude
%   For real symmetric problems, SIGMA may also be:
%      'LA' or 'SA' - Largest or Smallest Algebraic
%      'BE' - Both Ends, one more from high end if K is odd
%   For nonsymmetric and complex problems, SIGMA may also be:
%      'LR' or 'SR' - Largest or Smallest Real part
%      'LI' or 'SI' - Largest or Smallest Imaginary part
%   If SIGMA is a real or complex scalar including 0, EIGS finds the eigenvalues
%   closest to SIGMA.  For scalar SIGMA, and also when SIGMA = 'SM' which uses
%   the same algorithm as SIGMA = 0, B need only be symmetric (or Hermitian)
%   positive semi-definite since it is not Cholesky factored as in the other cases.
%
%   EIGS(A,K,SIGMA,OPTS) and EIGS(A,B,K,SIGMA,OPTS) specify options:
%   OPTS.issym: symmetry of A or A-SIGMA*B represented by AFUN [{0} | 1]
%   OPTS.isreal: complexity of A or A-SIGMA*B represented by AFUN [0 | {1}]
%   OPTS.tol: convergence: Ritz estimate residual <= tol*NORM(A) [scalar | {eps}]
%   OPTS.maxit: maximum number of iterations [integer | {300}]
%   OPTS.p: number of Lanczos vectors: K+1<p<=N [integer | {2K}]
%   OPTS.v0: starting vector [N-by-1 vector | {randomly generated by ARPACK}]
%   OPTS.disp: diagnostic information display level [0 | {1} | 2]
%   OPTS.cholB: B is actually its Cholesky factor CHOL(B) [{0} | 1]
%   OPTS.permB: sparse B is actually CHOL(B(permB,permB)) [permB | {1:N}]
%
%   EIGS(AFUN,N) accepts the function AFUN instead of the matrix A.
%   Y = AFUN(X) should return
%      A*X            if SIGMA is not specified, or is a string other than 'SM'
%      A\X            if SIGMA is 0 or 'SM'
%      (A-SIGMA*I)\X  if SIGMA is a nonzero scalar (standard eigenvalue problem)
%      (A-SIGMA*B)\X  if SIGMA is a nonzero scalar (generalized eigenvalue problem)
%   N is the size of A. The matrix A, A-SIGMA*I or A-SIGMA*B represented by AFUN is
%   assumed to be real and nonsymmetric unless specified otherwise by OPTS.isreal
%   and OPTS.issym. In all these EIGS syntaxes, EIGS(A,...) may be replaced by
%   EIGS(AFUN,N,...).
%
%   EIGS(AFUN,N,K,SIGMA,OPTS,P1,P2,...) and EIGS(AFUN,N,B,K,SIGMA,OPTS,P1,P2,...)
%   provide for additional arguments which are passed to AFUN(X,P1,P2,...).
%
%   Examples:
%      A = delsq(numgrid('C',15));  d1 = eigs(A,5,'SM');
%   Equivalently, if dnRk is the following one-line function:
%      function y = dnRk(x,R,k)
%      y = (delsq(numgrid(R,k))) \ x;
%   then pass dnRk's additional arguments, 'C' and 15, to EIGS:
%      n = size(A,1);  opts.issym = 1;  d2 = eigs(@dnRk,n,5,'SM',opts,'C',15);
%
%   See also EIG, SVDS, ARPACKC.
%
%   Copyright 1984-2002 The MathWorks, Inc.
%   $Revision: 1.45 $  $Date: 2002/05/14 18:50:58 $


% arguments_Afun = varargin{7-Amatrix-Bnotthere:end};
%(pour l'instant : n'accepte que 2 arguments dans le cas de Afun : Afun(W,X))
%permet d'aller plus vite en minimisant les acces a varargin 
%



nombre_A_times_X = 0; 
nombreIterations = 0; 

cputms = zeros(5,1);
t0 = cputime; % start timing pre-processing

if (nargout > 3)
    error('Too many output arguments.')
end

% Process inputs and do error-checking
if isa(varargin{1},'double')
    A = varargin{1};
    Amatrix = 1;
else
    A = fcnchk(varargin{1});
    Amatrix = 0;
end

isrealprob = 1; % isrealprob = isreal(A) & isreal(B) & isreal(sigma)
if Amatrix
    isrealprob = isreal(A);
end

issymA = 0;
if Amatrix
    issymA = isequal(A,A');
end

if Amatrix
    [m,n] = size(A);
    if (m ~= n)
        error('A must be a square matrix or a function which computes A*x.')
    end
else
    n = varargin{2};
    nstr = 'Size of problem, ''n'', must be a positive integer.';
    if ~isequal(size(n),[1,1]) | ~isreal(n)
        error(nstr)
    end
    if (round(n) ~= n)
        warning('MATLAB:eigs:NonIntegerSize',['%s\n         ' ...
            'Rounding input size.'],nstr)
        n = round(n);
    end
    if issparse(n)
        n = full(n);
    end
end

Bnotthere = 0;
Bstr = sprintf(['Generalized matrix B must be the same size as A and' ...
        ' either a symmetric positive (semi-)definite matrix or' ...
        ' its Cholesky factor.']);
if (nargin < (3-Amatrix-Bnotthere))
    B = [];
    Bnotthere = 1;
else
    Bk = varargin{3-Amatrix-Bnotthere};
    if isempty(Bk) % allow eigs(A,[],k,sigma,opts);
        B = Bk;
    else
        if isequal(size(Bk),[1,1]) & (n ~= 1)
            B = [];
            k = Bk;
            Bnotthere = 1;
        else % eigs(9,8,...) assumes A=9, B=8, ... NOT A=9, k=8, ...
            B = Bk;
            if ~isa(B,'double') | ~isequal(size(B),[n,n])
                error(Bstr)
            end
            isrealprob = isrealprob & isreal(B);
        end
    end
end

if Amatrix & ((nargin - ~Bnotthere)>4)
    error('Too many inputs.')
end

if (nargin < (4-Amatrix-Bnotthere))
    k = min(n,6);
else
    k = varargin{4-Amatrix-Bnotthere};
end

kstr = ['Number of eigenvalues requested, k, must be a' ...
        ' positive integer <= n.'];
if ~isa(k,'double') | ~isequal(size(k),[1,1]) | ~isreal(k) | (k>n)
    error(kstr)
end
if issparse(k)
    k = full(k);
end
if (round(k) ~= k)
    warning('MATLAB:eigs:NonIntegerEigQty',['%s\n         ' ...
            'Rounding number of eigenvalues.'],kstr)
    k = round(k);
end

whchstr = sprintf(['Eigenvalue range sigma must be a valid 2-element string.']);
if (nargin < (5-Amatrix-Bnotthere))
    % default: eigs('LM') => ARPACK which='LM', sigma=0
    eigs_sigma = 'LM';
    whch = 'LM';
    sigma = 0;
else
    eigs_sigma = varargin{5-Amatrix-Bnotthere};
    if isstr(eigs_sigma)
        % eigs(string) => ARPACK which=string, sigma=0
        if ~isequal(size(eigs_sigma),[1,2])
            whchstr = [whchstr sprintf(['\nFor real symmetric A, the choices are ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.'], ...
                    'LM','SM','LA','SA','BE')];
            whchstr = [whchstr sprintf(['\nFor non-symmetric or complex A, the choices are ''%s'', ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.\n'], ...
                    'LM','SM','LR','SR','LI','SI')];
            error(whchstr)
        end
        eigs_sigma = upper(eigs_sigma);
        if isequal(eigs_sigma,'SM')
            % eigs('SM') => ARPACK which='LM', sigma=0
            whch = 'LM';
        else
            % eigs(string), where string~='SM' => ARPACK which=string, sigma=0
            whch = eigs_sigma;
        end
        sigma = 0;
    else
        % eigs(scalar) => ARPACK which='LM', sigma=scalar
        if ~isa(eigs_sigma,'double') | ~isequal(size(eigs_sigma),[1,1])
            error('Eigenvalue shift sigma must be a scalar.')
        end
        sigma = eigs_sigma;
        if issparse(sigma)
            sigma = full(sigma);
        end
        isrealprob = isrealprob & isreal(sigma);
        whch = 'LM';
    end
end

tol = eps; % ARPACK's minimum tolerance is eps/2 (DLAMCH's EPS)
maxit = [];
p = [];
info = int32(0); % use a random starting vector
display = 1;
cholB = 0;

if (nargin >= (6-Amatrix-Bnotthere))
    opts = varargin{6-Amatrix-Bnotthere};
    if ~isa(opts,'struct')
        error('Options argument must be a structure.')
    end
    
    if isfield(opts,'issym') & ~Amatrix
        issymA = opts.issym;
        if (issymA ~= 0) & (issymA ~= 1)
            error('opts.issym must be 0 or 1.')
        end