index code

src/index_estrada.m

@@ -0,0 +1,42 @@
+% EE = index_estrada(c, ab )
+% EE = index_estrada(Afuns, n, nZ, N, ab)
+% EE = index_estrada(As, nZ, N, ab)
+%
+% Compute the Estrada index tr(exp(A)) by the Chebyshev expansion of 
+% exponential function and moments of A; the spectrum of As should 
+% already lie in [-1,1], and the rescaling parameters should be given
+% as the last argument.
+%
+% Inputs:
+%    c: A column vector of N Chebyshev moment estimates of A
+%    As: Rescaled matrix or function to apply matrix (to multiple RHS)
+%    n: Dimension of the space (if A is a function)
+%    nZ: Number of probe vectors with which we want to compute moments
+%    N: Number of moments to compute
+%    ab: Rescaling parameters. Original matrix should be ab(1)*A+ab(2)
+%
+% Output:
+%    EE: Estrada index of the form tr(exp(A))
+%
+function EE = index_estrada(varargin)
+    % If rescaling parameters are given
+    if length(varargin{end})==2
+        ab = varargin{end};
+        varargin = varargin(1:end-1);
+    else
+        ab = [1,0];
+    end
+
+    % Use the moments if they are given; otherwise compute them
+    if min(size(varargin{1}))==1 && length(varargin{1})>1
+        c = varargin{1};
+    else
+        [c,~] = moments_cheb_dos(varargin{:});
+    end
+
+    % Get the coefficients of Chebyshev expansion of exponential
+    N = length(c);
+    w = moments_exponential(N,ab(1));
+    EE=exp(ab(2))*w'*c;
+
+end

src/index_inform.m

@@ -0,0 +1,52 @@
+% IC = index_inform(A, nZ, N)
+%
+% Compute the information centrality by the Chebyshev expansion of 
+% 1/(ab(1)*x_ab(2)) and moments of (L+ee^T/n).
+%
+%    A: Adjacency matrix
+%    nZ: Number of probe vectors with which we want to compute moments
+%    N: Number of moments to compute
+%
+% Output:
+%    IC: Information centrality
+%
+function IC = index_inform(A, nZ, N)
+    if nargin<3
+        N = 150;
+    end
+
+    if nargin<2
+        nZ=100;
+    end
+
+    n = size(A,1);
+
+    % Compute the eigenrange of Laplacian
+    L = matrix_laplacian(A);
+    opts = [];
+    opts.isreal = 1;
+    opts.issym = 1;
+    range = sort(eigs(L, 2, 'sm', opts));
+    range(3) = eigs(L, 1, 'lm', opts);
+    range(1) = 1;
+    range = sort(range);
+    range = range([1,end]);
+
+    % Linear rescaling of (L+ee^T/n)
+    Lfun = @(X) L*X;
+    Bfun = @(X) (bsxfun(@plus,mean(X,1),Lfun(X)));
+    [Bfuns,ab] = rescale_mfunc(Bfun,n,range);
+
+    % Compute the ldos moments
+    [c,~] = moments_cheb_ldos(Bfuns,n,nZ,N);
+
+    % Compute the Chebyshev expansion of 1/(ab(1)x+ab(2))
+    z = -ab(2)/ab(1);
+    w = 1./(sqrt(z^2-1)*(z-sqrt(z^2-1)).^(0:N-1))/ab(1);
+    w(2:N)=2*w(2:N);
+
+    % Compute information centrality
+    d = c'*w'-(n-1)/n^2;
+    IC = 1./(d+mean(d)-2/n^2);
+
+end

src/index_sub_exp.m

@@ -0,0 +1,54 @@
+% ESC = index_sub_exp(c, ab)
+% ESC = index_sub_exp(beta, c, ab)
+% ESC = index_sub_exp(beta, Afuns, n, nZ, N, ab)
+% ESC = index_sub_exp(beta, As, nZ, N, ab)
+%
+% Compute the exponential subgraph centrality diag(exp(beta*A)) by the
+% Chebyshev expansion of exponential function and moments of A; the 
+% spectrum of As should already lie in [-1,1], and the rescaling 
+% parameters should be given as the last argument.
+%
+% Inputs:
+%    beta: Exponential coefficient
+%    c: An N-by-n matrix of ldos moment estimates of A
+%    As: Rescaled matrix or function to apply matrix (to multiple RHS)
+%    n: Dimension of the space (if A is a function)
+%    nZ: Number of probe vectors with which we want to compute moments
+%    N: Number of moments to compute
+%    ab: Rescaling parameters. Original matrix should be ab(1)*A+ab(2)
+%
+% Output:
+%    ESC: Exponential subgraph centraility of the form diag(exp(beta*A))
+%
+function ESC = index_sub_exp(varargin)
+    % If rescaling parameters are given
+    if length(varargin{end})==2
+        ab = varargin{end};
+        varargin = varargin(1:end-1);
+    else
+        ab = [1,0];
+    end
+
+    % Check if an exponential coefficient is given (default:1)
+    if isnumeric(varargin{1}) && isscalar(varargin{1})
+        beta = varargin{1};
+        varargin = varargin(2:end);
+    else
+        beta = 1;
+    end
+
+    ab = beta*ab;
+
+    % Check if ldos moments are given; otherwise calculate an estimate
+    if size(varargin{1},1)~=size(varargin{1},2)
+        c = varargin{1};
+    else
+        [c,~] = moments_cheb_ldos(varargin{:});
+    end
+
+    % Get the Chebyshev expansi0n of exp(beta*x)
+    N = size(c,1);
+    w = moments_exponential(N,ab(1));
+    ESC = exp(ab(2))*(c'*w);
+
+end

src/index_sub_res.m

@@ -0,0 +1,59 @@
+% ESC = index_sub_res(c, ab)
+% ESC = index_sub_res(alpha, c, ab)
+% ESC = index_sub_res(alpha, Afuns, n, nZ, N, ab)
+% ESC = index_sub_res(alpha, As, nZ, N, ab)
+%
+% Compute the resolvent subgraph centrality diag((I-alpha*A)^-1) by 
+% the Chebyshev expansion of exponential function and moments of A; 
+% the spectrum of As should already lie in [-1,1], and the rescaling 
+% parameters should be given as the last argument. 
+%
+% Inputs:
+%    alpha: Resolvent coefficient
+%    c: An N-by-n matrix of ldos moment estimates of A
+%    As: Rescaled matrix or function to apply matrix (to multiple RHS)
+%    n: Dimension of the space (if A is a function)
+%    nZ: Number of probe vectors with which we want to compute moments
+%    N: Number of moments to compute
+%    ab: Rescaling parameters. Original matrix should be ab(1)*A+ab(2)
+%
+% Output:
+%    RSC: Resolvent subgraph centraility of the form diag((I-alpha*A)^-1)
+%
+function RSC = index_sub_res(varargin)
+    % If rescaling parameters are given
+    if length(varargin{end})==2
+        ab = varargin{end};
+        varargin = varargin(1:end-1);
+    else
+        ab = [1,0];
+    end
+
+    % Check if an exponential coefficient is given (default:0.9)
+    if isnumeric(varargin{1}) && isscalar(varargin{1})
+        alpha = varargin{1};
+        varargin = varargin(2:end);
+    else
+        alpha = 0.9;
+    end
+
+    % Throw a warning if choice of alpha does not guarantee convergence
+    if alpha>=1/sum(ab)
+        warning(['Chebyshev expansion will not converge because'...
+                ' function has pole in [-1,1]. alpha must be less'...
+                ' than 1/sum(ab)']);
+    end
+
+    % Check if ldos moments are given; otherwise calculate an estimate
+    if size(varargin{1},1)~=size(varargin{1},2)
+        c = varargin{1};
+    else
+        [c,~] = moments_cheb_ldos(varargin{:});
+    end
+
+    % Get the Chebyshev expansion of 1/(1-alpha*x)
+    N = size(c,1);
+    w = moments_resolvent(N,alpha,ab);
+    RSC = c'*w;
+
+end

src/moments_exponential.m

@@ -0,0 +1,22 @@
+% [c] = moments_exponential(N)
+% [c] = moments_exponential(N, beta)
+%
+% Compute Chebyshev moments 0 through N-1 of exponential function
+% f(x)=exp(beta*x) on [-1,1]
+%
+% Input:
+%    N: Number of moments
+%    beta: Exponential coefficient
+%
+% Output:
+%    c: A column vector of N moments
+%
+function c = moments_exponential(N,beta)
+    if nargin<2
+        beta = 1;
+    end
+
+    c = besseli(0:N-1,beta);
+    c(2:N) = c(2:N)*2;
+    c=c';
+end

src/moments_resolvent.m

@@ -0,0 +1,30 @@
+% [c] = moments_resolvent(N)
+% [c] = moments_resolvent(N, alpha)
+%
+% Compute Chebyshev moments 0 through N-1 of the function
+% f(x)=1/(1-alpha*x) on [-1,1]
+%
+% Input:
+%    N: Number of moments
+%    alpha: Resolvent coefficient in (0,1)
+%    ab: Rescaling parameters. Original matrix should be ab(1)*A+ab(2)
+%
+% Output:
+%    c: A column vector of N moments
+%
+function c = moments_resolvent(N,alpha,ab)
+    if nargin<3
+        ab = [1,0];
+    end
+    if nargin<2
+        alpha = 0.9;
+    end
+
+    ab = alpha*ab;
+
+    z = (1-ab(2))/ab(1);
+    c = 1./(sqrt(z^2-1)*(z+sqrt(z^2-1)).^(0:N-1))/ab(1);
+    c(2:N) = c(2:N)*2;
+    c=c';
+
+end