Accumulated symmetric matrix eigendecomposition

      #include <gandalf/linalg/symmat_eigen.h>

There is also a specific module in Gandalf to compute the eigenvalues and eigenvectors of a positive semi-definite matrix accumulated as a sum of outer products of vectors. This is useful for instance when solving homogeneous linear equations; for an example see the computation of the fundamental & essential matrix in Section 5.2. The symmetric matrix $S$ is constructed as

$$S = \sum_{i=1}^n \sigma_i {\bf x}_i {\bf x}_i^{\top}$$ (3.2)

given the $n$ vectors ${\bf x}_i$ and weighting factors $s_i$, $i=1,\ldots,n$.

There is a structure to hold the accumulated matrix and the resulting eigendecomposition matrices:

      /* structure to hold accumulated symmetric matrix and resulting
       * eigendecomposition of a sum of vector outer products
       */
      typedef struct
      {
         Gan_SquMatrix SxxT; /* accumulated sum of vector outer products */
         Gan_SquMatrix W;    /* diagonal matrix of eigenvalues */
         Gan_Matrix    Z;    /* matrix of eigenvectors */
         Gan_Vector    work; /* workspace vector for computing eigendecomposition */

         /* whether this structure was dynamically allocated */
         Gan_Bool alloc;
      } Gan_SymMatEigenStruct;

To start the computation, create an instance of the structure using

      Gan_SymMatEigenStruct SymEigen;

      /* create structure for computing eigenvalues and eigenvectors,
         initialising accumulated matrix S (here 5x5) to zero */
      gan_symeigen_form ( &SymEigen, 5 );

This routine creates the structure and its constituent matrices, and initialises the outer product sum matrix, in this case a $5\times 5$ matrix, to zero. You can then add a term to the sum in Equation 3.2 using the routine

      /* increment matrix S by 
      gan_symeigen_increment ( &SymEigen, 2.0, /* weighting factor s */
                               3.4, 1.0, 9.7, 3.4, 2.1 ); /* elements of vector x */

Here the first value $2.0$ is the weighting factor $\sigma$ for this vector, and the elements of the vector ${\bf x}$ are passed into a variable argument list (and hence have to be explicitly double type). You should call gan_symeigen_increment() once for each of the $n$ vectors. Then solve for the eigenvalues & eigenvectors using the routine

      /* solve for eigenvalues and eigenvectors */
      gan_symeigen_solve ( &SymEigen );

after which the eigenvalues can be read back from SymEigen.W and the eigenvectors from SymEigen.Z. If you want to reuse the structure on a new eigendecomposition computation, call the routine

      /* reset accumulated symmetric to zero, optionally changing size of matrix */
      gan_symeigen_reset ( &SymEigen, 5 );

where the last argument allows you to change the size of the matrix if desired. Finally to free the structure use

      /* free eigensystem structure and constituent matrices */
      gan_symeigen_free ( &SymEigen );

Error detection: sym_symeigen_form() returns a pointer to the eigensystem structure (&SymEigen), and so returns NULL on error. sym_symeigen_increment(), sym_symeigen_solve() and sym_symeigen_reset() return boolean values, so GAN_FALSE indicates an error. In all cases the Gandalf error handler is invoked.