pagerank.m

基于matlab的pagerank代码

for j=1:k-1
   beta=norm(f(:,j));
   V(:,j+1)=f(:,j)/beta;
   ejt=[zeros(1,j-1) beta];
   Hhat=[H; ejt];
   w=A(V(:,j+1));
   h=V(:,1:j+1)'*w;
   f(:,j+1)=w-V(:,1:j+1)*h;
   H=[Hhat h];
end

% Extend Arnoldi factorization
beta=norm(f(:,k));
V(:,k+1) = f(:,k)/beta;
ejt=[zeros(1,k-1) beta];
H=[H ;ejt];

% ===================================
% pagerank_approx
% ===================================

function [x flag hist dt] = pagerank_approx(A, options)
% use the power iteration algorithm

if (options.verbose > 0)
   fprintf('approximate computation...\n');
end;
 

tol = options.tol;
v = options.v;
maxiter = options.maxiter;
c = options.c;
bp = options.approx_bp;
subiter = options.approx_subiter;
boundary = options.approx_boundary;


n = size(A,1);

%x = v;
%if (isfield(options, 'x0'))
%    x = options.x0;
%end;

if (length(find(v)) ~= n)
    global_pr = 0;
else
    global_pr = 1;
    error('pagerank:invalidParameter',...
        'approximation computations are not implemented for global pagerank yet');
end;

    
hist = zeros(maxiter,1);
delta = 1;
iter = 0;
dt = 0;

% set the initial set of seed pages
if (global_pr)
    if (isfield(options, 'x0'))
        % the seed pages come from the x0 vector if provided
        p = find(options.x0);
        x = x0(p);
    else
        % the seed pages come from the x0 vector (otherwise, choose random)
        p = unique(ceil(rand(250,1)*size(P,1)));
        x = ones(length(p),1)./length(p);
    end;
else
    % the seed pages come from the x0 vector
    p = find(v);
    x = ones(length(p),1)./length(p);
    v = v(p);
end;

local = [];
active = p;
frontier = p;

tic;
while (iter <= maxiter && delta > tol)
    % expand all pages
    if (boundary == 1)
        
        % if we are running the boundary algorithm...
        
        [ignore sp] = sort(-x);
        cs = cumsum(x(sp));
        spactive = active(sp);
        allexpand_ind = cs < (1-bp);
        % actually, we need to add the first 0 after the last 1 in
        % allexpand_ind because we need cumsum to be larger than 1-bp
        allexpand_ind(min(find(allexpand_ind == 0))) = ~0;
        allexpand = spactive(allexpand_ind);
        toexpand = setdiff(allexpand,local);
    else
        %
        % otherwise, just expand all pages with a sufficient tolerance
        %
        allexpand = active(x > bp);
        toexpand = setdiff(allexpand,local);
    end;
    
    if (length(toexpand) > 0)
        xp = zeros(n,1);
        xp([local frontier]) = x;

        local = [local toexpand];
        frontier = setdiff(find(sum(A(local,:),1)), local);
        active = [local frontier];

        x = xp(local);
    else
        xp = zeros(n,1);
        xp([local frontier]) = x;
        x = xp(local);
    end;
    
    Lp = A(local,active);
    outdegree = full(sum(Lp,2));
    outdegree = [outdegree; zeros(length(frontier),1)];

    siter = 0;
    L = [Lp; sparse(length(frontier),length(active))];
    x2 = [x; xp(frontier)];
    while (siter < subiter)
        y = full(c*L'*(invzero(outdegree).*x2));
        omega = 1 - norm(y,1);
        
        % the ordering of local is preseved, so these are always the 
        % correct vertices
        y(1:length(p)) = y(1:length(p)) + omega*v;
        
        x2 = y;
        
        siter = siter+1;
    end;
    
    
    x2 = [x; xp(frontier)];
    
    delta = norm(y-x2,1);
    hist(iter+1) = delta;
    
    if (options.verbose > 0)
        fprintf('iter=%02i; delta=%0.03e expand=%i\n', iter, delta, length(toexpand));
    end
    
    x = y;
    iter = iter + 1;
end;
dt = toc;
% resize hist
hist = hist(1:iter);

xpartial = x;
x = zeros(n,1);
x([local frontier]) = xpartial;

flag = 0;

if (delta > tol && iter == maxiter)
    warning('pagerank:didNotConverge', ...
        'The PageRank algorithm did not converge after %i iterations', ...
        maxiter);
    flag = 1;
end;

no% ===================================
% pagerank_eval
% ===================================
function [x flag hist dt] = pagerank_eval(P,options)

algs = {'power', 'gs', 'arnoldi', 'linsys', 'linsys'};
extra_opts = {struct(''), struct(''), struct(''), struct(''), ...
    struct('linsys_solver',@(f,v,tol,its) gmres(f,v,8,tol, its))};
names = {'power', 'gs', 'arnoldi8', 'bicgstab', 'gmres8'};

v = options.v;
c = options.c;

x = cell(5,1);
flag = cell(5,1);
hist = cell(5,1);
dt = cell(5,1);

web('text://<html><body>Generating PageRank report...</body></html>','-noaddressbox');

htmlend = '</body></html>';
s = {};
s{1} = '<html><head><title>PageRank runtime report</title></head><body><h1>PageRank Report</h1>';

stemp = s;
stemp{end+1} = '<p>Generating graph statistics...</p>';
stemp{end+1} = htmlend;

A = spones(P);
d = dangling(P);

npages = size(P,1);
nedges = nnz(P);
ndangling = sum(d);
maxindeg = full(max(sum(A,1)));
maxoutdeg = full(max(sum(A,2)));
ncomp = components(A);

s{end+1} = '<h2>Graph statistics</h2>';
s{end+1} = '<table border="0" cellspacing="4">';
s{end+1} = sprintf('<tr><td style="font-weight: bold">%s:</td><td>%i</td></tr>', ...
    'Number of pages', npages);
s{end+1} = sprintf('<tr><td style="font-weight: bold">%s:</td><td>%i</td></tr>', ...
    'Number of edges', nedges);
s{end+1} = sprintf('<tr><td style="font-weight: bold">%s:</td><td>%i</td></tr>', ...
    'Number of dangling nodes', ndangling);
s{end+1} = sprintf('<tr><td style="font-weight: bold">%s:</td><td>%i</td></tr>', ...
    'Max in-degree', maxindeg);
s{end+1} = sprintf('<tr><td style="font-weight: bold">%s:</td><td>%i</td></tr>', ...
    'Max out-degree', maxoutdeg);
s{end+1} = sprintf('<tr><td style="font-weight: bold">%s:</td><td>%i</td></tr>', ...
    'Number of strong components:', ncomp);
s{end+1} = '</table>';

sOut = [stemp{:}];
web(['text://' sOut],'-noaddressbox');

s{end+1} = '<h2>Algorithm performance</h2>';
s{end+1} = '<table border="0">';
s{end+1} = sprintf('<tr><td style="text-align: right">%s</td><td>%0.3f</td></tr>', ...
    'c = ', c);
s{end+1} = sprintf('<tr><td style="text-align: right">%s</td><td>%2.2e</td></tr>', ...
    'tol = ', options.tol);
s{end+1} = sprintf('<tr><td style="text-align: right">%s</td><td>%i</td></tr>', ...
    'maxiter = ', options.maxiter);
s{end+1} = '</table>';

s{end+1} = '<table border="0">';
s{end+1} = ['<tr style="text-align: left">' ...
    '<th style="border-bottom:solid 1px">Algorithm</th>' ...
    '<th style="border-bottom:solid 1px">Time</th>' ...
    '<th style="border-bottom:solid 1px">Iterations</th>' ...
    '<th style="border-bottom:solid 1px">Error</th></tr>'];

for (ii=1:length(algs))
    alg = algs{ii};
    extra_opt = extra_opts{ii};
    name = names{ii};
    
    stemp = s;
    stemp{end+1} = '</table>';
    stemp{end+1} = sprintf('<p>Solving for PageRank with %s...</p>', char(name));
    stemp{end+1} = htmlend;

    sOut = [stemp{:}];
    web(['text://' sOut],'-noaddressbox');
    
    extra_opt = merge_structs(struct('alg',char(alg)),extra_opt);
    
    [pi flagi histi dti] = pagerank(P, merge_structs(extra_opt,options));
    
    p{ii} = pi;
    flag{ii} = flagi;
    hist{ii} = histi;
    dt{ii} = dti;
    
    err = norm(pi - c*(pi'*P)' - c*(d'*pi)*v - (1-c)*v,1);
    
    if (mod(ii,2) == 0)
        s{end+1} = sprintf('<tr style="background-color: #cccccc"><td>%s</td><td>%.2f</td><td>%i</td><td>%2.2e</td></tr>',...
            char(name), dti, length(histi), err);
    else
        s{end+1} = sprintf('<tr><td>%s</td><td>%.2f</td><td>%i</td><td>%2.2e</td></tr>',...
            char(name), dti, length(histi), err);
    end;
end

s{end+1} = '</table>';


s{end+1} = htmlend;
sOut = [s{:}];
web(['text://' sOut],'-noaddressbox');

s{end+1} = sprintf('<tr><td>%s</td><td></td><td></td><td></td></tr>',char(name));

%
% plot the time histogram
%
figure(1);
close(1);
figure(1);

dts = cell2mat(dt);
flags = cell2mat(flag);

h2 = bar(dts.*(flags==0));

    set(h2,'FaceColor',[1 1 1]);
    set(h2,'LineWidth',2.0);
    
    set(gca,'XTick', 1:length(algs));
    set(gca,'XTickLabel',names);
    ylabel('time (sec)');
    

%
% plot the history results
%
figure(2);
close(2);
figure(2);
lso = get(0,'DefaultAxesLineStyleOrder');
lsc = get(0,'DefaultAxesColorOrder');

lso = {'o-', 'x:', '+-.', 's--', 'd-'};

    nlso = length(lso);
    curlso = 0;
    nlsc = length(lsc);
    curlsc = 0;
    
    for ii=1:length(algs)
        histi = hist{ii};
        %legendname = fn{ii};
        %line(1:length(mrval.hist), mrval.hist);
        semilogy(1:length(histi),histi,...
            lso{mod(curlso,nlso)+1}, ...
            'Color',lsc(mod(curlsc,nlsc)+1,:),...
            'MarkerSize',3);
        hold on;
        curlso = curlso+1;
        curlsc = curlsc+1;
    end;

    title('PageRank algorithm convergence (WARNING: DIFFERENT Y-SCALES)');
    xlabel('iteration')';
    ylabel('convergence measure');
    
    legend(names{:});
    

function S = merge_structs(A, B)
% MERGE_STRUCTS Merge two structures.
%
% S = merge_structs(A, B) makes the structure S have all the fields from A
% and B.  Conflicts are resolved by using the value in A.
%

%
% merge_structs.m
% David Gleich
%
% Revision 1.00
% 19 Octoboer 2005
%

S = A;

fn = fieldnames(B);

for ii = 1:length(fn)
    if (~isfield(A, fn{ii}))
        S.(fn{ii}) = B.(fn{ii});
    end;
end;

function P = normout(A)
% NORMOUT Normalize the outdegrees of the matrix A.  
%
% P = normout(A)
%
%   P has the same non-zero structure as A, but is normalized such that the
%   sum of each row is 1, assuming that A has non-negative entries. 
%

%
% normout.m
% David Gleich
%
% Revision 1.00
% 19 Octoboer 2005
%

% compute the row-sums/degrees
d = full(sum(A,2));

% invert the non-zeros in the data
id = invzero(d);

% scale the rows of the matrix
P = diag(sparse(id))*A;

function v = invzero(v)
% INVZERO Compute the inverse elements of a vector with zero entries.
%
% iv = invzero(v) 
%
%   iv is 1./v except where v = 0, in which case it is 0.
%

%
% invzero.m
% David Gleich
%
% Revision 1.00
% 19 Octoboer 2005
%

% sparse input are easy to handle
if (issparse(v))
    [m n] = size(v);
    [i j val] = find(v);
    val = 1./val;
    v = sparse(i,j,val,m,n);
    return;
end;

% so are dense input
 
% compute the 0 mask
zm = abs(v) > eps(1);

% invert all non-zeros
v(zm) = 1./v(zm);

function dmask = dangling(A)
% DANGLING  Compute the indicator vector for dangling links for webgraph A
% 
%  d = dangling(A)
%

d = full(sum(A,2));
dmask = d == 0;

function [k,sizes]=components(A)

% based on components.m from (MESHPART Toolkit)
% which had 
% John Gilbert, Xerox PARC, 8 June 1992.
% Copyright (c) 1990-1996 by Xerox Corporation.  All rights reserved.
% HELP COPYRIGHT for complete copyright and licensing notice

[p,p,r,r] = dmperm(A|speye(size(A)));
sizes = diff(r);
k = length(sizes);