// // Peter Richter (prichter@princeton.edu) // Princeton University Undergraduate Mathematics Laboratory // May 11, 2001 // // walk.cpp - compute the second-largest eigenvalue of each bipartite graph // of order N and regularity K along a random walk of length M. // #include #include #include #include #define FALSE 0 #define TRUE 1 #define K 3 /* number of edges meeting each vertex */ #define N 500 /* even number of vertices in the graph */ #define M 1000 /* number of steps in the random walk */ #define EPS 1e-8 /* tolerance for power method convergence */ //----------------------------------------------------------------------------- double error_bound(double e, double x[N], int g[N][K]) { /* bound the error in the eigenvalue e using the eigenvector x of e^2 */ // compute the 2-norm of x double x_norm = 0; for (int i = 0; i < N; i++) x_norm += x[i] * x[i]; x_norm = sqrt(x_norm); // normalize x double y[N]; for (int i = 0; i < N; i++) y[i] = x[i] / x_norm; // apply the adjacency matrix to y double t[N]; for (int i = 0; i < N; i++) { t[i] = 0; for (int j = 0; j < K; j++) t[i] += y[g[i][j]]; } // apply the adjacency matrix to t double z[N]; for (int i = 0; i < N; i++) { z[i] = 0; for (int j = 0; j < K; j++) z[i] += t[g[i][j]]; } // compute the error vector zeta double zeta[N]; for (int i = 0; i < N; i++) zeta[i] = z[i] - (e * e) * y[i]; // compute the 2-norm of zeta double zeta_norm = 0; for (int i = 0; i < N; i++) zeta_norm += zeta[i] * zeta[i]; zeta_norm = sqrt(zeta_norm); // return the error bound on e return sqrt(zeta_norm); } //----------------------------------------------------------------------------- double eig_two(int g[N][K], double p) { /* compute the second largest eigenvalue of the adjacency matrix for g */ // choose an initial vector x if this routine has not been called before static int first_call = TRUE; static double x[N]; if (first_call == TRUE) { for (int i = 0; i < N; i++) x[i] = ((double) rand()) / RAND_MAX; first_call = FALSE; } // calculate the second eigenvalue using the power method double new_approx = -K-1, old_approx = -K-2; while (fabs(new_approx - old_approx) >= EPS || error_bound(new_approx, x, g) >= fabs(new_approx - p)) { // compute the orthogonality coefficients double c = 0, d = 0; for (int i = 0; i < N; i++) { c += x[i]; d += (i < N/2) ? x[i] : -x[i]; } c /= N; d /= N; // project out the two principal eigenvectors from x double y[N]; for (int i = 0; i < N; i++) y[i] = (i < N/2) ? x[i]-c-d : x[i]-c+d; // compute the 2-norm of y double y_norm = 0; for (int i = 0; i < N; i++) y_norm += y[i] * y[i]; y_norm = sqrt(y_norm); // normalize y double z[N]; for (int i = 0; i < N; i++) z[i] = y[i] / y_norm; // apply the adjacency matrix to z double t[N]; for (int i = 0; i < N; i++) { t[i] = 0; for (int j = 0; j < K; j++) t[i] += z[g[i][j]]; } // apply the adjacency matrix to t for (int i = 0; i < N; i++) { x[i] = 0; for (int j = 0; j < K; j++) x[i] += t[g[i][j]]; } // compute the square root of (x,z) double lambda = 0; for (int i = 0; i < N; i++) lambda += x[i] * z[i]; lambda = sqrt(lambda); // store this estimate old_approx = new_approx; new_approx = lambda; } return new_approx; } //----------------------------------------------------------------------------- void shuffle(int a[]) { /* permute randomly the elements of the array a */ // do the fisher-yates shuffle! for (int i = (N/2)-1; i >= 0; i--) { int j = rand() % (i+1); int t = a[i]; a[i] = a[j]; a[j] = t; } } //----------------------------------------------------------------------------- void walk_init(int g[N][K]) { /* begin the random walk along the space of bipartite K-regular graphs */ // seed the random number generator long seed = time(NULL); srand(seed); int okay; // flag to indicate the graph is simple and connected int v[K][N/2]; // arrays to hold K permutations of vertices N/2...N-1 do { // reset flag to true okay = TRUE; // fill each array v[i] with the independent set of vertices N/2...N-1 for (int i = 0; i < K; i++) for (int j = 0; j < N/2; j++) v[i][j] = N/2 + j; // shuffle each array v[i] randomly; each (j, v[i][j]) represents an edge for (int i = 0; i < K; i++) shuffle(v[i]); // check that the graph is simple for (int i = 0; i < N/2; i++) for (int j = 0; j < K; j++) for (int r = 0; r < K; r++) if (r != j && v[j][i] == v[r][i]) okay = FALSE; // re-shuffle if the graph is not simple if (okay == FALSE) continue; // store the graph in the N-by-K table g for (int i = 0; i < K; i++) for (int j = 0; j < N/2; j++) g[j][i] = v[i][j]; for (int i = N/2; i < N; i++) for (int j = 0, s = 0; j < N/2; j++) for (int r = 0; r < K; r++) if (g[j][r] == i) g[i][s++] = j; // check that the graph is connected; loop again if it is not if (fabs(eig_two(g, 1/EPS) - K) < sqrt(EPS)) okay = FALSE; } while (okay == FALSE); } //----------------------------------------------------------------------------- void walk_step(int g[N][K]) { /* step to the next graph in the random walk */ int okay; // flag to indicate the next step preserves structure do { // reset flag to true okay = TRUE; // pick two distinct vertices from independent set of vertices 0...(N/2)-1 int v2, v1 = rand() % (N/2); do { v2 = rand() % (N/2); } while (v2 == v1); // pick distinct vertices connected to v1 (but not v2) and v2 (but not v1) int i1 = rand() % K, w1 = g[v1][i1]; int i2 = rand() % K, w2 = g[v2][i2]; for (int i = 0; i < K; i++) if (w1 == g[v2][i] || w2 == g[v1][i]) okay = FALSE; // loop again if this was not accomplished if (okay == FALSE) continue; // make the next step in the walk g[v1][i1] = w2; g[v2][i2] = w1; for (int i = 0; i < K; i++) { if (g[w2][i] == v2) g[w2][i] = v1; if (g[w1][i] == v1) g[w1][i] = v2; } } while (okay == FALSE); } //----------------------------------------------------------------------------- void main() { int g[N][K]; // the graph; each (i, g[i][j]) represents an edge double p = 2*sqrt(K-1); // critical value for characterizing ramanujan graphs double e; // second largest eigenvalue of the adjacency matrix // perform the random walk, recording the eigenvalue of each matrix walk_init(g); for (int i = 0; i < M; i++) { if ((e = eig_two(g, p)) <= p) // <- g is ramanujan cout << e << '\t' << '*' << endl; else // <- g is not ramanujan cout << e << endl; walk_step(g); } } //-----------------------------------------------------------------------------