PDL / Vector and matrix algebra / Sort Eigenvalues of symmetric matrix

Question

I'm trying to compute eigenvectors of symmetric matrices with the great eigens_sym function of the great PDL Perl module.

use strict; # Be pessimistic and careful with interpretations in the compilation stage, before the very execution of this program . Doubts may abort execution, probably with error messages to stdout 

use PDL ;   # Load the PDL extension of Perl 



my $matrix ; # declare this variable 

# Set $matrix to my matrix A 

$matrix = new PDL(
                [
                    [  2 ,  1 ,  9 ,  2 , -6 ,  5 , ] ,
                    [  0 ,  5 ,  3 ,  0 , -8 , -3 , ] ,
                    [  1 ,  0 ,  5 ,  1 , -6 ,  5 , ] ,
                    [ -4 , -6 , -7 ,  2 , -1 ,  5 , ] ,
                    [  5 , -2 ,  3 ,  1 , -8 ,  3 , ] ,
                    [  1 , -6 ,  9 , -2 ,  9 ,  5 , ]

                ]
               ) ; 

print "\nA: \n" ; 

print $matrix ; 

# Make and display $matrix2 = The transposition of A 

    my $matrix2 =  $matrix -> transpose ;

    print "\nA* : \n" ; 

    print $matrix2 ; 

my $matrix3 = $matrix2 x $matrix ; # Here $matrix3 is set to the matrix product A* A : 

print "\nA*A : \n" ; 

print $matrix3 ; 

The mathematics says, for any square matrix A made of real numbers, A* A is always symmetric (and also semidefinite positive, but this is not necessary for the PDL's eigens_sym function), and can help to tell us very important informations for understanding A and all it means for us

my ($ev,$e) = eigens_sym $matrix3 ;  # X-rated matrix, id est:  Compute the eigenvalues and the eigenvectors of A* A   !

In other words : diagonalize the matrix product A* A . Mathematics says, it is always possible to diagonalize any symmetric matrix of real numbers, but it can be very difficult as the dimension of the matrix increases .

print "\nEigenvalues with eigenvectors of A*A : \n\n" ; 

print $e ;  # The eigenvalues 

print $ev ; # " The eigenvectors are returned in COLUMNS of the returned PDL "

prints :

A: 

[
 [ 2  1  9  2 -6  5]
 [ 0  5  3  0 -8 -3]
 [ 1  0  5  1 -6  5]
 [-4 -6 -7  2 -1  5]
 [ 5 -2  3  1 -8  3]
 [ 1 -6  9 -2  9  5]
]


A*A : 

[
 [ 47  10  75   0 -45  15]
 [ 10 102   6   0 -78 -76]
 [ 75   6 254  -6 -44  80]
 [  0   0  -6  14 -46  18]
 [-45 -78 -44 -46 282 -20]
 [ 15 -76  80  18 -20 118]
]


Eigenvalues with eigenvectors of A*A : 

[ 23.285687  4.3552833  277.10958 0.79058644  367.06761  144.39125]

[
 [  0.86502676   0.39257997  0.090239841 -0.070634235  -0.25470872  -0.14000255]
 [ -0.36255636   0.53399064  -0.40332417  -0.31499006  -0.16574127  -0.54226155]
 [ -0.19994359  -0.20028032   0.58246948      0.19444  -0.60114242  -0.42598184]
 [ -0.14771441   0.50816552 -0.082459683   0.81070439 -0.092150685   0.21775057]
 [-0.037584019   0.29620226   0.55789823  0.032583693   0.69935299  -0.33082462]
 [ -0.23889757   0.41791719   0.41456683  -0.44685776  -0.22066884   0.58994146]
]

. So far so good, beware that here the "column" eigenvectors appear as rows [ ... , ... , etc etc ] , I' m very happy but have some problems :

1. sometimes I need to have the result with the eigenvalues [ 23.285687 4.3552833 277.10958 0.79058644 367.06761 144.39125] sorted in ascending or descending order, how can I have this ?

2. How can I free the memory used for the objects $matrix2 and $matrix3 after the above job ?

3. For further vector/matrix processing : are there a good tutorial and a manual for Perl PDL ?

BTW / Just a comment on the installation of the PDL module in your Perl system / I installed PDL after Perl: at first run I had a dependency problem due to many missing modules in my perl system (PDL is a big module), as of Nov 2023 I solved this serious problem by running from command line

#     $ sudo perl -MCPAN -e 'shell' 
#     <etc etc>
#     cpan[1]> install PDL 

# This automatic installation took several minutes to run and was apparently successful, although some error messages were issued

Answer 1

To answer your first question how to sort the eigenvalues: You can use the qsort function:

# Sort eigenvalues in ascending order
my $sorted_ev_asc = qsort($e);

# Sorted eigenvalues in descending order
my $sorted_ev_desc = $sorted_ev_asc->slice("-1:0");

For the second question: How to free the memory of $matrix3? Perl automatically manages memory for variables, and you generally don't need to explicitly free memory. However, if you want to release the memory occupied by a PDL variable explicitly, you can use the null function.

$matrix3 = PDL->null;

For the third question: Where to find a tutorial on PDL matrix operations? Have a look at this page for example: https://metacpan.org/pod/PDL::MatrixOps

Update:

To also sort the eigenvectors, you could use the qsorti function:

my ($evec,$eval) = eigens_sym $matrix3 ;  # Make the eigenvector 
# Get the indices of the sorted elements
my $indices = qsorti($eval);
# Sort eigenvectors based on the permutation indices
my $sorted_eigenvalues = $eval->slice($indices);
# Sort eigenvectors based on the permutation indices
my $sorted_eigenvectors = $evec->slice(':', $indices);

Edit2:

To change the number of significant digits in the output, there is a package variable called doubleformat that can be set locally (or globally). For example:

local $PDL::doubleformat = "%10.16g";  # Set to the desired precision
# Print sorted eigenvalues
say "Eigenvalues in ascending order: $sorted_eigenvalues";

Which gives output like this:

Eigenvalues in ascending order: 

  [-104.0242117201316  -50.37924528190175  -4.799244776313021 
     55.46970789893607 117.6355876358277  225.0974062435826
  ]

Answer 2

Thank you HÅKON HÆGLAND you SOLVED all problems , I'm completely happy.

Here is your solution in my arrangement:

 my ($evectors,$evalues) = eigens_sym $matrix3 ;  # Make eigenvector matrix ; here as  $matrix3  any real and (square) symmetric, non necessarily semidefinite positive matrix, is allowed 
  # " The eigenvectors are returned in COLUMNS " , but these columns may appear as rows such as [ 2 -3 5 ]  if you print the object

#  Great. But probably eigen-values and eigen-vectors are not sorted now! 

my $indices ; # create variable

  $indices = qsorti($evalues); # For eigenvalues in ASCENDING order 
# $indices = $indices->slice("-1:0"); # Add this line if you want eigenvalues in DESCENDING order

# Sort eigenvalues and eigenvectors
    $evalues  = $evalues  -> slice($indices); 
    $evectors = $evectors -> slice(':', $indices); 

# Done, now the eigen-values and -vectors are sorted! For example you can print the result

# You can change the number of  sig fig  aka significant figures in the floating-point number print formats 
 local $PDL::doubleformat = "% 10.16g";  # Set to the desired precision
# Available also Perl's new hex-bin float format, example "% 18.13a" , where MAYBE the .13 could give us the perfect and complete representation of the hardware thingy behind Double-precision floating-point format ) 

 print "\nEigenvalues with eigenvectors : \n\n" ; 

 print $evalues ;  # The eigenvalues 

 print $evectors ; # The eigenvectors