Hello,

I am trying to use the GNU Scientific Library (GSL) to compute the leading eigenvector (corresponding to the largest eigenvalue) of a matrix M of size N (N ~ 100). M is square and symmetric. All its elements are real and positive, with values between 0 and 1. As such, the Perron–Frobenius theorem dictates that the leading eigenvector has all positive values.

I have tried this with a number of different sets of values for the elements of M. Some of the time, the values of the leading eigenvector printed out are all positive, as expected. They seem to correspond with what I would expect the leading eigenvector to be. However, the rest of the time, all of the elements in the eigenvector are negative, between 0 and -1. They do not seem to correspond at all to what I would expect, eg. if I made all the elements positive, they still seem to be inconsistent to what I would expect.

As such, this program is working for some sets of data for M, but not for others. I cannot see a pattern in the datasets between the correct solutions and the incorrect solutions.

As an example of one of the incorrect solutions, if I set element (i, j) of M to be equal to (i * j), then the leading eigenvector has entirely negative values; but they should all be positive! Please see the code below for this implementation.

Any ideas on what I am doing wrong?

Thanks.

Code:
// Initialize the size of the matrix
int N = 100;

// Create the matrix M
gsl_matrix* M = gsl_matrix_alloc(N, N);

// Assign the elements of M
for (int i = 0; i < N; i++)
{
	for (int j = 0; j < N; j++)
	{
		gsl_matrix_set(M, i, j, i * j);
	}
}

// Create the vectors to store the eigenvalues and eigenvectors
gsl_vector* evals = gsl_vector_alloc(N);
gsl_matrix* evecs = gsl_matrix_alloc(N, N);

// Allocate an eigenvector / eigenvalue workspace
gsl_eigen_symmv_workspace* workspace = gsl_eigen_symmv_alloc(N);

// Compute the eigenvalues and eigenvectors
gsl_eigen_symmv(M, evals, evecs, workspace);

// Find the index of the leading eigenvalue
int index = gsl_vector_max_index(evals);

/// Print out the leading eigenvector
for (int i = 0; i < N; i++)
{
	cout << gsl_matrix_get(evecs, i, index) << endl;
}

// Free the gsl memory
gsl_eigen_symmv_free(workspace);
gsl_matrix_free(M);
gsl_vector_free(evals);
gsl_matrix_free(evecs);