Re: [Help-gsl] Real generalized symmetric-definite eigensystems

[Patrick Alken]

Hello,

  Positive definite does not mean all the matrix entries are positive.
It means all the eigenvalues of the matrix are positive. In your case,
both H and S have negative eigenvalues, so neither matrix is positive
definite. In this case, you need to use the Real Generalized
Nonsymmetric eigensolver (even though your matrices are symmetric).

Alternatively, since your S matrix has det(S) = 1, instead of solving

H x = \lambda S x

You could solve:

S^{-1} H x = \lambda x

in which case you could use the Real Nonsymmetric eigensolver on the
matrix S^{-1} H.

Patrick

On 09/28/2016 11:53 PM, Dimitrova, Maria wrote:
> Hello,
>
> I have yet another question about matrix operations, this time regarding 
> eigensystems. In a sample program I define two matrices and try to solve the 
> eigenvalue problem. Both matrices are symmetric and all of their elements are 
> positive. However, on calling the function gsl_eigen_gensymmv I get an error 
> that they have to be positive definite. And well, they are. The function 
> gsl_matrix_isnonneg returns 1.
>
> Here is the output, followed by the source code. Thank you very much.
>
> $ gcc -Wall -Wextra diagonalization.c -lm -lgsl -lgslcblas && ./a.out
> =======================================================================================
>
> Initialized matrices
> Matrix H
>
>  1.000000   2.000000   4.000000
>  2.000000   3.000000   5.000000
>  4.000000   5.000000   6.000000
>
> Matrix S
>
>  7.000000   8.000000   10.000000
>  8.000000   9.000000   11.000000
>  10.000000   11.000000   12.000000
>
> =======================================================================================
>
> Check for positive definiteness
> "H: "
> stat=1
> "S: "
> stat=1
> =======================================================================================
>
> Generalized eigenvalue problem
> gsl: cholesky.c:157: ERROR: matrix must be positive definite
> Default GSL error handler invoked.
> Aborted (core dumped)
>
>
> /************************************************************************************************
>  *   COMPILED WITH THE LINE:                                                  
>                   *
>  *   gcc -Wall -Wextra diagonalization.c -lm -lgsl -lgslcblas && ./a.out      
>                   *
>  
> ************************************************************************************************/
> #include <stdio.h>
> #include <stdlib.h>
> #include <math.h>
> #include <gsl/gsl_math.h>
> #include <gsl/gsl_blas.h>
> #include <gsl/gsl_vector.h>
> #include <gsl/gsl_matrix.h>
> #include <gsl/gsl_eigen.h>
> #define MatrixOrder 3
> #define MSG(msg) printf( "\n" #msg "\n")
> #define DIV 
> printf("\n=======================================================================================\n\n")
> #define DEBUG(msg, var, fmt) printf( #msg "\n" #var "=%" #fmt "\n", var)
> #define PRINTF(var, fmt) printf("\n**DEBUG: " #var "=%" #fmt "\n", var)
> void printArr(gsl_matrix *array, int nrows, int ncols)
> {
>     int i,j;
>     printf("\n");
>     for (i=0;i<nrows;i++)
>     {
>         printf("\n");
>         for (j=0;j<ncols;j++)
>         {
>             printf("% f  ",gsl_matrix_get(array,i,j));
>         }
>     }
>     printf("\n\n");
> }
> int main()
> {
>     int elem=0; // used in the initialization of the arrays
>     int stat=0;
>     int i,j;
>     gsl_matrix *C, *H, *S;
>     gsl_vector *E;
>     gsl_eigen_gensymmv_workspace *wrkEig; // for generalized eigenvalue 
> problem
>     S = gsl_matrix_calloc(MatrixOrder,MatrixOrder);
>     H = gsl_matrix_calloc(MatrixOrder,MatrixOrder);
>     C = gsl_matrix_calloc(MatrixOrder,MatrixOrder);
>     E = gsl_vector_calloc(MatrixOrder);
>     wrkEig = gsl_eigen_gensymmv_alloc(MatrixOrder); // for NxN matrix wrk is 
> O(4N)
> /*****************************************************************************************/
>     // Initialization of the matrices
>     for (i=0;i<MatrixOrder;i++)
>     {
>         for (j=0;j<=i;j++)
>         {
>             elem++;
>             gsl_matrix_set(H,i,j,elem);
>             gsl_matrix_set(H,j,i,elem);
>         }
>     } // H[i][j] initialized
>     DIV;
>     MSG(Initialized matrices);
>     MSG(Matrix H);
>     printArr(H,MatrixOrder,MatrixOrder);
>     for (i=0;i<MatrixOrder;i++)
>     {
>         for (j=0;j<=i;j++)
>         {
>             elem++;
>             gsl_matrix_set(S,i,j,elem);
>             gsl_matrix_set(S,j,i,elem);
>         }
>     } // S[i][j] initialized
>     MSG(Matrix S);
>     printArr(S,MatrixOrder,MatrixOrder);
> /*****************************************************************************************/
>
>     DIV;
>     MSG(Check for positive definiteness);
>     stat=gsl_matrix_isnonneg(H); //return 1 if all the elements of the matrix 
> m are non-negative, and 0 otherwise
>     DEBUG("H: ",stat,d);
>     DEBUG("S: ",stat,d);
>
>     DIV;
>     MSG(Generalized eigenvalue problem);
>     gsl_eigen_gensymmv(H,S,E,C,wrkEig);
>
> /*****************************************************************************************/
>     gsl_matrix_free(H);
>     gsl_matrix_free(S);
>     gsl_vector_free(E);
>     gsl_matrix_free(C);
>     gsl_eigen_gensymmv_free(wrkEig);
>     return 0;
> }
>
>
>
> Best regards,
>
> Maria Dimitrova
>
> address@hidden