GNU Scientific Library eigenvectors

[karnavor]

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.

// 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);

[superbonzo]

... dictates that the leading eigenvector has all positive values.

this cannot be the case, just consider that if v is an eigenvector then -v is an eigenvector as well. Quoting wikipedia, the theorem says that "There exists an eigenvector v [...] such that all components of v are positive". Therefore, your program is correct as long as the found leading eigenvector has non-zero components with the same sign.

[karnavor]

Hi, thanks for your reply. The Wikipedia article states that:

"the Perron–Frobenius theorem ... asserts that a real square matrix with positive entries has a unique largest real eigenvalue and that the corresponding eigenvector has strictly positive components."

My program finds the largest eigenvalue (I have checked this value to verify it is the largest), and then finds the associated eigenvector, which it prints out. So according to the theory, these should all be positive values. Or are you saying that a vector with all negative values is identical to a vector with the same values, but all positive?

I have been playing around with the code, and I've noticed two strange behaviours:
1. If N <= 10, then the leading eigenvector has all positive values, as expected. But if N > 10, then all the values are negative.
2. If I add a very small value to each element in M, then the leading eigenvector has all positive values, as expected. ie. The matrix assignment becomes: gsl_matrix_set(M, i, j, 0.0001 + i * j);

Thanks for any help.

[Lindley]

Or are you saying that a vector with all negative values is identical to a vector with the same values, but all positive?

Not identical, but equivalent in this sense. Consider a 2D Gaussian. This can be represented by an ellipse. The first eigenvector tells you the direction of the ellipse's major axis. It doesn't matter if you reverse that vector; it's still describing the direction of the ellipse's major axis. (The first eigenvalue would in this case be proportional to the length of the major axis. The second eigenvector and eigenvalue describe the minor axis similarly.)

[superbonzo]

Hi, thanks for your reply. The Wikipedia article states that:

"the Perron–Frobenius theorem ... asserts that a real square matrix with positive entries has a unique largest real eigenvalue and that the corresponding eigenvector has strictly positive components."

My program finds the largest eigenvalue (I have checked this value to verify it is the largest), and then finds the associated eigenvector, which it prints out. So according to the theory, these should all be positive values. Or are you saying that a vector with all negative values is identical to a vector with the same values, but all positive?

I have been playing around with the code, and I've noticed two strange behaviours:
1. If N <= 10, then the leading eigenvector has all positive values, as expected. But if N > 10, then all the values are negative.
2. If I add a very small value to each element in M, then the leading eigenvector has all positive values, as expected. ie. The matrix assignment becomes: gsl_matrix_set(M, i, j, 0.0001 + i * j);

Thanks for any help.

they are not "identical", they are both eigenvectors. That is, if your program outputs a vector "v" with negative components then "-v" is the eigenvector with positive component you're looking for.

whether your algorithm outputs a positive or negative vector depends on the specific way that algorithm diagonalizes your input matrix; that is, unless specified otherwise by the algorithm documentation, it will "randomly" give one or the other case ( although apparently "strange" patterns can occur, of course ).

BTW, note that if you're working with floating point numbers you could as well observe inconsistencies in the sign of the output vector, when dealing with big matrices and/or small matrix components ...