From 0db8e68a8431633f6eec08a41f87878eb7fe1c46 Mon Sep 17 00:00:00 2001
From: Anna Yesypenko <ay256@cornell.edu>
Date: Fri, 9 Dec 2016 18:31:10 -0500
Subject: [PATCH] find sparse eigenvector given eigenvalue and region of nodes
 that support eigenvector

---
 src/eigenbasis.m        | 26 ++++++++++++++++++
 src/freq_cheb.m         | 37 +++++++++++++++++++++++++
 src/recover_eigvec.m    | 60 +++++++++++++++++++++++++++++++++++++++++
 test/test_recover_eig.m | 34 +++++++++++++++++++++++
 4 files changed, 157 insertions(+)
 create mode 100644 src/eigenbasis.m
 create mode 100644 src/freq_cheb.m
 create mode 100644 src/recover_eigvec.m
 create mode 100644 test/test_recover_eig.m

diff --git a/src/eigenbasis.m b/src/eigenbasis.m
new file mode 100644
index 0000000..c11b362
--- /dev/null
+++ b/src/eigenbasis.m
@@ -0,0 +1,26 @@
+% Given a matrix {A}, an eigenvalue {lambda}, and
+% the support {s} of the corresponding eigenbasis V,
+% return V.
+function [L,V] = eigenbasis(A, lambda, s, rtol)
+
+  % Form block matrix reordering rows
+  I = find(s);
+  
+  [Lambda, Q] = eig(full(A(I, I)));
+  
+  M = A(:,I);
+  V = zeros(length(A), 1);
+  L = zeros(1, 1);
+  
+  for j=1:length(Lambda)
+      l = Lambda(j,j);
+      if (abs(lambda - l) < rtol)      
+          v = zeros(length(A),1);
+          v(I) = Q(:,j);
+          resid = norm (M * v(I) - l * v);
+          if (resid < rtol)
+              L = [L; l];
+              V = [V; v];
+          end
+      end
+  end
\ No newline at end of file
diff --git a/src/freq_cheb.m b/src/freq_cheb.m
new file mode 100644
index 0000000..64cb926
--- /dev/null
+++ b/src/freq_cheb.m
@@ -0,0 +1,37 @@
+% [calc_eig] = freq_cheb (m_seq)
+%
+% Given Chebyshev moment sequence {m_seq} for LDoS on node k,
+% return {calc_eig}, the eigenvalues with correponding
+% eigenvectors supported on node k.
+
+function [calc_eig] = freq_cheb (m_seq)
+
+    m = length(m_seq); % number of moments
+    d = m / 2; % degree
+
+    H = hankel(m_seq(1:d+1), m_seq(d+1: m));
+    [~,~,V] = svd(H);
+    p = V(:, end);
+    
+    % Find the coefficients of the polynomial p given by
+    % Chebyshev moments c_{km} and roots of p
+    % p(x) = a_0 + a_1 x + \cdots + a_{d-1} x^{d-1} + a_d
+    r = roots (flipud(p));
+    % Fix roots so they are on the unit circle
+    for j=1:length(r)
+        n = norm([real(r(j)), imag(r(j))]); % should be approx 1
+        r(j) = r(j) / n;
+    end
+    
+    calc_eig = zeros(length(r), 1);
+    
+    for j=1:length(r)
+        % r_j = e^{i w_j}
+        % w_j = arccos (eig_j)
+        w_j = log(r(j)) / sqrt(-1);
+        eig_j = cos(w_j);
+        assert (abs(imag(eig_j)) < 10^(-10));
+        calc_eig(j,1) = real(eig_j);
+    end
+    calc_eig = abs(calc_eig);
+end
\ No newline at end of file
diff --git a/src/recover_eigvec.m b/src/recover_eigvec.m
new file mode 100644
index 0000000..0946e13
--- /dev/null
+++ b/src/recover_eigvec.m
@@ -0,0 +1,60 @@
+% [V] = recover_eigenvec (A, lambda, center, max_hops)
+%
+% Determine the (sparse) eigenvector corresponding to lambda 
+% that it is supported in neighborhood within h hops of specified center.
+% 
+% Inputs:
+% - A:       given symmetric matrix in R^{n x n}
+% - lambda:  eigenvalue of A
+% - center:  a node known to support eigenvector, for center in [1, n]
+% - h:       explore at most h hops from specified center
+
+function [V, s] = recover_eigvec(A, lambda, center, h, rtol)
+    n = length(A);
+    
+    % TODO how do we determine d?
+    num_probe = 10;
+    d = 10; % degree of function
+    m = 2 * d; % number of moments
+    
+    [c,~] = moments_cheb_ldos (A, num_probe, m, 1);
+    
+    % Initialize visited array
+    distance = zeros (1, n);
+    distance(:) = intmax; distance(center) = 0;
+    % Initialize queue for BFS
+    queue = zeros (1, n);
+    head = 1; last = 2;
+    queue(head) = center;
+    % Initialize support for eigenvector correponding to lambda
+    s = zeros (n, 1);
+    
+    % Iterate through nodes in h-hop neighborhood
+    while (head < last)
+        node = queue(head);
+        head = head + 1;
+        d = distance(node);
+        if (d == h)
+            continue;
+        end
+        % Get the one-hop neighborhood
+        neigh = A (node,:);
+        for i=1:n
+            if (neigh(i) ~= 0 && distance(i) == intmax || i == node)
+                distance(i) = d + 1;
+                m_seq = c(:,i);
+                supp_eigs = freq_cheb(m_seq);
+                if (any(abs(supp_eigs - lambda) < rtol))
+                    queue(last) = i;
+                    last = last + 1;
+                    s(i) = 1;
+                end
+            end
+        end
+    end
+    
+    % Power iteration will converge to eigenvector v s.t. Av = lambda * v
+    % and v and s have same support
+    [~,V] = eigenbasis(A, lambda, s, rtol);
+
+end
\ No newline at end of file
diff --git a/test/test_recover_eig.m b/test/test_recover_eig.m
new file mode 100644
index 0000000..c563c09
--- /dev/null
+++ b/test/test_recover_eig.m
@@ -0,0 +1,34 @@
+% K_5 with whiskers
+A = [0 1 1 1 1 1 0 0;
+     1 0 1 1 1 0 1 0;
+     1 1 0 1 1 0 0 1;
+     1 1 1 0 1 0 0 0;
+     1 1 1 1 0 0 0 0;
+     1 0 0 0 0 0 0 0;
+     0 1 0 0 0 0 0 0;
+     0 0 1 0 0 0 0 0];
+ 
+A = matrix_normalize(A, 's');
+[U,D,~] = svd(full(A));
+d = diag(D);
+hops = 4;
+rtol = 1e-10;
+
+sparse_vs = 0;
+solved_vs = 0;
+
+for j=1:length(A)
+    lambda = d(j);
+    v = U(:,j);
+    nonzeros = find(abs(v) > rtol);
+    if (length(nonzeros) > 3/4 * length(A))
+        continue
+    end
+    center = nonzeros(1);
+    sparse_vs = sparse_vs + 1;    
+    [V,s] = recover_eigvec(A, lambda, center, hops, rtol);
+    
+    if (norm(A * V - lambda * V) < rtol && ~all(V))
+        solved_vs = solved_vs + 1;
+    end
+end