Eigenvalue bug?

[Michael Meyer <spyqqqdia at yahoo dot com>]

Hello,

Let b[i], i=1,...,n be an increasing sequence (here
b[i]=sqrt(i+1)).
The matrix A defined as A_{ij}=b(min{i,j})/b(max{i,j})
is then known to be positive semidefinite.

Let lambda[i] be the sequence of eigenvalues of A and
U (= evec in gsl notation) the matrix whose columns
are the eigenvectors of A.

Let D=diag(sqrt(lambda[i])) denote the
diagonal matrix with the square root of the
eigenvalues of A on its diagonal and set

  L=U*D (matrix product).

One can then easily show that LL'=A (the prime denotes
transpose).

The following code sets out to check this equality,
but the equality works only in dimension n=2 and fails
in higher dimensions.

Is this a bug or can I simply not code this right?
Sorry for the akward placement of comments, I had to
keep 
the lines short or else yahoo mail messes these up.

C++ CODE:


#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <iostream>
#include <stdio.h>
#include <math.h>
#include <unistd.h>

using namespace std;

 int main(int argc, char* argv[])
 {
   int n=10; //all matrices n by n
   
   //sequence b[i]=sqrt(i+1):
   double* b=new double[n];
   for(int i=0;i<n;i++) b[i]=sqrt(i+1);             

   //initialize matrix A=0
   gsl_matrix* A=gsl_matrix_calloc(n,n);
   
   //set A_{ij}=b(min{i,j})/b(max{i,j}):
   for(int i=0;i<n;i++)
   for(int j=0;j<n;j++)
   {
     //x=b(min{i,j})/b(max{i,j})
     double x=(i<j)? b[i]/b[j] : b[j]/b[i];
     gsl_matrix_set(A,i,j,x);                         
 
   }

   //prepare eigenvector computation
   gsl_matrix* evec=gsl_matrix_calloc(n,n);
   gsl_vector* eval=gsl_vector_calloc(n);
   gsl_eigen_symmv_workspace* 
   w=gsl_eigen_symmv_alloc(n);


   //compute eigenvalues and eigenvectors of A 
   gsl_eigen_symmv(A,eval,evec,w);     

   //eigenvalues are now in vector eval,
   //eigenvectors the columns of matrix evec


   //eigenvectors before sorting
   for(int i=0;i<n;i++) 
   cerr<<gsl_vector_get(eval,i)<<endl;
   cerr<<endl;
  
   //sort eigenvalues descending, 
   //eigenvectors correspondingly
   gsl_eigen_symmv_sort(eval,evec,
                        GSL_EIGEN_SORT_VAL_DESC);  
                                                      
     

   //load eigenvalues into lambda[]
   //lambda[i]=eigenvalue_i:
   double* lambda=new double[n];
   for(int i=0;i<n;i++) 
   lambda[i]=gsl_vector_get(eval,i);   

   
   //allocate matrix diag(sqrt(lambda[j]))
   //diagonal matrix with sqrt(eigenvalue_i) 
   //down the diagonal:
   gsl_matrix* diagSqrtLambda=gsl_matrix_calloc(n,n);
   for(int i=0;i<n;i++)
   for(int j=0;j<n;j++)
   {
     double x = (!(i==j))? 0 : sqrt(lambda[j]);
     gsl_matrix_set(diagSqrtLambda,i,j,x);
   }

    //allocate L=evec*diagSqrtLambda
    //L initialized as L=0:
    gsl_matrix* L=gsl_matrix_calloc(n,n);
    //L=L+evec*diagLambda:   
    gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,
                   1.0,evec,diagSqrtLambda,1.0,L);    

    //allocate B=LL'
    //B initialized as B=0:
    gsl_matrix* B=gsl_matrix_calloc(n,n);
    //B=B+LL':                                
    gsl_blas_dgemm(CblasNoTrans,CblasTrans,
                   1.0,L,L,1.0,B);               

    //check if A=B, i.e. A-B=0
    cerr<<endl<<endl<<"Matrix A-LL':"<<endl<<endl;
    //A=A-B:
    gsl_matrix_sub(A,B);
    //print each entry of A on a separate line:       
                                         
    gsl_matrix_fprintf(stderr,A,"%f");
                                       

  return 0;
}
  


Thanks,

Mike


__________________________________________________
Do You Yahoo!?
NEW from Yahoo! GeoCities - quick and easy web site hosting, just $8.95/month.
http://geocities.yahoo.com/ps/info1