forum.alglib.net

Hmmm, maybe I've found error. You have 200x200=40000 variables. How you can construct 200x200 covariation matrix from them? You should have 40000x40000 matrix!

Wow... It's correct! Covariance and correlation matrix must be bigger! But... How can i calculate them? :?

Do you really have only two values for each variables? I.e. only two points, each in the space with 200x200 dimension? Or it is just example, and actual problem will have many more points? If only two points, then your correlation matrix will include only +1's, zeros and -1's. Correlation on two-element sample is always +1, 0 or -1.

As for the calculation, you can't use ALGLIB for realtime calculation of eigenvalue of 40000x40000 matrix. Too much data, too much time.
Maybe some iterative solver like ARPACK will help you.
For the case with two points in the sample (two very "heavy" points, I should say), you can use analytic form of the correlation matrix to implicitly calculate C*x products.
For the many-point samples you should do some clever optimizations to calculate such products efficiently. Some analytic form for product of covariance matrix and vector.

Hi, i'm back! The problem still continues...Now i use the covariance matrix but eigenvalues remain negative! The matrix is correct (the number of values for each variable is dynamic, however i follow your advice and i put at least 3 values), but now rmatrixevd and smatrixevd returns different results!
P.S. i believe that the covariance matrix is semidefinite-positive...I'm right?
Thanks
P.P.S. with smatrixevd all the result eigenvalues are negative! (now covariance matrix dimension is dynamic, so i'm trying with small numbers like 100 by 100)

Yes, you are right - it must be semidefinite. Can you post here your code, which solves simplified, reduced problems (now full code, with all initializations) and for which rmatrixevd and smatrixevd return different results?

wait, does it mean that different entries of this matrix were calculated using different number of (different) values? Generally, such matrix is not guaranteed to be semi-definite.

const static int x=3; //numero di immagini
const static int dim=20;//dimensione immagini
double cov[dim*dim][dim*dim];
//calcolo la media delle immagini
double media[dim][dim];
//inizializzo matrice
for(int i=0;i<dim;i++){
for(int j=0;j<dim;j++){
media[i][j]=0;
}
}
//calcolo matrice somma di tutte le matrici immagine
for(int n=0;n<x;n++){
for(int i=0;i<dim;i++){
for(int j=0;j<dim;j++){
media[i][j]+= (double) pixel[i][j][n];
}//fine for j
}//fine for i
}//fine for n
//media aritmetica
for(int i=0;i<dim;i++){
for(int j=0;j<dim;j++){
media[i][j]=media[i][j]/x;
}
}
//inizializzo matrice di covarianza
for(int i=0;i<dim*dim;i++){
for(int j=0;j<dim*dim;j++){
cov[i][j]=0;
}//fine for j
}//fine for i
//calcolo matrice di covarianza
double scarto[dim][dim][x];
//inizializzazione matrice deviazione standard
for(int n=0;n<x;n++){
for(int i=0;i<dim;i++){
for(int j=0;j<dim;j++){
scarto[i][j][n]=0;
}//fine for j
}//fine for i
}//fine for n
//calcolo lo scarto membro a membro
for(int n=0;n<x;n++){
for(int i=0;i<dim;i++){
for(int j=0;j<dim;j++){
scarto[i][j][n]+=pixel[i][j][n]-media[i][j];
}//fine for j
}//fine for i
}//fine for n
//calcolo della matrice di covarianza
double tmp=0;
int l,m,y=0;
for(int n=0;n<x;n++){
l=0;
m=0;
y=0;
for(int k=0;k<dim*dim;k++){
for(int i=0;i<dim;i++){
for(int j=0;j<dim;j++){
tmp=scarto[l][m][n]*scarto[i][j][n];
cov[k][y]+=tmp;
y++;
}//fine for j
}//fine for i
y=0;
if(m==dim-1){
l++;
m=0;
}
else m++;
}//k
}//fine for n
for(int i=0;i<dim*dim;i++){
for(int j=0;j<dim*dim;j++){
cov[i][j]=cov[i][j]/x;
}//fine for j
}//fine for i
//calcolo autovalori e autovettori attraverso libreria ALGLIB
stringstream matrice;
for(int i=0;i<dim*dim;i++){
if(i==0) matrice<<"[[";
else matrice<<"[";
for(int j=0;j<dim*dim;j++){
if(j<dim*dim-1) matrice<<cov[i][j]<<",";
else matrice<<cov[i][j];
}
if(i<dim*dim-1) matrice<<"],";
else matrice<<"]]";
}
string s=matrice.str();
const char *c= s.c_str();
alglib::real_2d_array mat(c);
alglib::real_1d_array autovalRe;
alglib::real_1d_array autovalIm;
alglib::real_2d_array autovettDx;
alglib::real_2d_array autovettSx;
//alglib::real_2d_array matrix=mat;
//procedura per il calcolo degli autovettori e degli autovalori di una matrice simmetrica
bool ok=smatrixevd(mat,dim*dim,1,false,autovalRe,autovettDx);
//procedura per il calcolo degli autovettori e degli autovalori di una matrice reale
//bool ok=rmatrixevd(matrix,dim*dim,1,autovalRe,autovalIm,autovettSx,autovettDx);
if(ok){
double autovalMod[dim*dim];
//trovo il modulo di tutti gli autovalori
for(int i=0;i<dim*dim;i++){
if(autovalRe[i]<0) autovalMod[i]=-autovalRe[i];
else autovalMod[i]=autovalRe[i];
}
//selezione dinamica degli autovalori principali (data una soglia->0.9)
//prender? i primi M che soddisfano la relazione |Sommatoria dei primi M/(sommatoria dei restanti dim-M)|>=0.9
bool ext=false;
int M=0;
while(!ext){
double principali=0;
for(int i=0;i<M+1;i++){
principali+=autovalMod[i];
}
double rimanenti=0;
for(int i=0;i<dim;i++){
rimanenti+=autovalMod[i];
}
if((abs(principali/rimanenti))>=(0.9)) ext=true;
else{
if(M==(dim-1)) break;
M++;
}
}
//normalizzo gli autovettori
for(int i=0;i<M+1;i++){
long double vett=0;
for(int j=0;j<dim*dim;j++){
vett+=autovettDx[j][i]*autovettDx[j][i];
}
vett=sqrt(vett);
for(int j=0;j<dim*dim;j++){
autovettDx[j][i]=autovettDx[j][i]/vett;
}
}
double adjustedData[dim][dim];//matrice temporanea contenente l'immagine normalizzata
for(int i=0;i<dim;i++){
for(int j=0;j<dim;j++){
adjustedData[i][j]= 0;
}
}
//procedura per il calcolo della proiezione di ogni immagine nello spazio degli autovettori
for(int n=0;n<x;n++){
for(int i=0;i<dim;i++){
for(int j=0;j<dim;j++){
adjustedData[i][j]= (double) pixel[i][j][n] - media[i][j];
}
}
//allocazione mat dinamica (matrice delle features, cio? degli autovettori selezionati)
double** featureV;
featureV = new double*[M+1];
for(int i=0; i<M+1; i++){
featureV[i] = new double[dim*dim];
}
for(int i=0;i<dim*dim;i++){
for(int j=0;j<M+1;j++){
featureV[j][i]=autovettDx[i][j];//creazione matrice features gi? trasposta
}
}
//matrice dinamica che conterr? la proiezione delle immagini nell'autospazio
double** immComp;
immComp = new double*[M];
for(int i=0; i<M+1; i++){
immComp[i] = new double[dim*dim];
}
for(int i=0;i<dim*dim;i++){
for(int k=0;k<M+1;k++){
immComp[k][i]=0;
}
}
//prodotto matrice features (Trasp) per matrice immagine (trasposta)
for(int i=0;i<dim*dim;i++){
for(int k=0;k<M+1;k++){
double acc=0;
for(int j=0;j<dim;j++){
acc+=featureV[k][j]*adjustedData[j][i];
}
immComp[k][i]=acc;
}
}
stringstream percorso;
percorso<<""<<"c:\\ImmCompressa"<<n<<".txt";
stringstream output;
for(int i=0;i<M+1;i++){
for(int j=0;j<dim*dim;j++){
output<<immComp[i][j]<<" ";
}
output<<"\n";
}
ofstream testo(percorso.str());
testo<<output.str();
testo.close();
}//fine for n

}//fine if calcolo autovettori ed autovalori
else{
cout<<"Si e' verificato un errore. Impossibile calcolare gli autovalori.";
}

Data is stored in int pixel[dim][dim][x]...With another function i catch data from grayscale bitmaps...As you can see this is a pca algorithm or at least should be so!