User Tools

Site Tools


Calculating Eigenvalues

Calculating all eigenvalues and eigenvectors of a symmetric matrix in one swoop results in a highly non linear and non-convex model and NLP solvers most likely will have problems solving these. However, one may compute the eigenvalues one by one by maximizing a norm and iteratively introducing a constraint that makes the new eigenvector orthogonal to the ones already found. A little trick in the generation of initial values removes components parallel to already found eigenvectors and ensures that the initial point to each solve is feasible. All constraints except one are linear and the objective is quadratic, so the model is easy to solve.

Alternatively, the tool eigenvector that is distributed with GAMS can be called. The model below (contributed by Arne Drud) gives an example of this approach for a Covariance matrix.

$Title Covariance Matrix and its Eigenvalues and -vectors.
    j    'variables'         /j1*j25/
    n    'observations'      /n1*n40/

Alias (j,jj,jjj);

    jn1(j,n)    'observations of the variables'
    meanj(j)    'mean value of the observations'
    covj(j,jj)  'First covariance matrix';

* forming a covariance matrix
jn1(j,n)   = normal(0,1);
meanj(j)   = sum(n, jn1(j,n))/card(n);
covj(j,jj) = sum(n, (jn1(j,n) - meanj(j))*(jn1(jj,n) - meanj(jj)))/card(n);

    evecj(j,j)  'Eigenvectors found so far'
    evalj(j)    'Eigenvalues found so far';

*  Define a model that finds the largest eigen value and eigen vector:
Set       jev(j)   'Subset already having an associated eigenvector-eigenvalue pair';
Scalar    normrv   'Norm of initial random vector';

jev(j) = no;
evecj(j,j) = 0;

    evec(j)    'Eigenvector'
    eval       'Eigenvalue';

    eq_obj     'Definition of Objective'
    eq_orth(j) 'The new Eigenvector must be orthogonal to earlier Eigenvectors'
    eq_norm    'The norm of the eigen vector must be 1';

eq_obj        .. eval =e= sum((j,jj), evec(j)*covj(j,jj)*evec(jj) );

eq_orth(jev)  .. sum( j, evecj(j,jev) * evec(j) ) =e= 0;

eq_norm       .. sum(j, sqr(evec(j)) ) =e= 1;

Model mod_ev / all /;

parameter helper(j);

option limRow=0, limCol=0, solPrint=silent, solveLink=%solveLink.loadLibrary%;
*  We initialize the eigenvector with random numbers. To prevent using
*  a lot of time in phase 1 we make sure the orthogonality constraints
*  are satisfied by subtracting a multiple of the already found
*  eigenvectors. We do not check that the remaining vector still is
*  nonzero.
   evec.l(j)     = uniform(-1,1);
   helper(jev)   = sum(j, evec.l(j)*evecj(j,jev));
   evec.l(j)     = evec.l(j) - sum(jev,helper(jev)*evecj(j,jev));
   normrv        = sqrt( sum(j, sqr(evec.l(j)) ) );
   evec.l(j)     = evec.l(j) / normrv;
   solve mod_ev using nlp maximizing eval;
   abort$(mod_ev.solvestat ne 1) 'solvestat not 1', mod_ev.solvestat;
   abort$(mod_ev.modelstat ne 1 and mod_ev.modelstat ne 2) 'modelstat not 1 or 2';
   jev(jjj)      = Yes;
   evecj(j,jjj)  = evec.l(j);
   evalj(jjj)    = eval.l;

    testj(j,j) 'New covariance matrix'
    difj(j,j)  'Difference between original and new covariance';

testj(j,jj) = sum( jjj, evecj(j,jjj) * evalj(jjj) * evecj(jj,jjj) );
difj(j,jj)  = covj(j,jj) - testj(j,jj);
abort$(smax((j,jj), abs(difj(j,jj)))>1e-6) difj;

* Use eigenvector tool to calculate Eigenvectors and Eigenvalues
execute_unload 'a.gdx', j, covj;
execute 'eigenvector.exe a.gdx j covj b.gdx evalj evecj';
execute_load 'b.gdx', evalj, evecj;
testj(j,jj) = sum( jjj, evecj(j,jjj) * evalj(jjj) * evecj(jj,jjj) );
difj(j,jj)  = covj(j,jj) - testj(j,jj);
abort$(smax((j,jj), abs(difj(j,jj)))>1e-6) difj;
IMPRESSUM / LEGAL NOTICEPRIVACY POLICY gams/calculating_eigenvalues.txt · Last modified: 2020/05/28 07:57 by Michael Bussieck