Reputation: 226

Compute the eigenvalues of v ⊗ v

I need to compute the eigenvalues of a symmetric matrix

M = v ⊗ v,

where v is a known, one dimensional, array (a vector) and ⊗ denotes the outer product.

I currently compute v ⊗ v explicitly via a self defined function and then call the lapack routine dsyev to get the desired eigenvalues.

Is there a more efficient method to calculate those eigenvalues (maybe directly from v?) using lapack, blas or a self written algorithm in fortran 2003?

Upvotes: 1

Answers (1)

Ian Bush

Reputation: 7432

It is shown on the sister maths site that

The eigenvalues for this problem are all zero, except one which has value the length squared of v
The eigenvector corresponding to the unique value is the (normalized) vector itself
As a result of 1, 2 and the symmetry of the matrix the N-1 eigenvectors with eigenvalue zero span a space orthogonal to the original vector

Below is a naive implementation of the method. I use a simple Modified Gram-Schmidt to generate the evecs, there is almost certainly a better way. It checks the results and compares the time with LAPACK. The accuracy is comparable to LAPACK, and on single thread on my laptop using openblas and gfortran 9.4 it is about twice as quick as LAPACK, as shown by the runs at the end. Note the first run is compiled with all run time checks turned on to demonstrate code correctness. Again you can almost certainly do better than this.

ijb@ijb-Latitude-5410:~/work/stack$ cat outer_eigen.f90
Module numbers_module

  Use, Intrinsic :: iso_fortran_env, Only : wp => real64
  
  Implicit None

  Public :: wp

  Private

End Module numbers_module

Module lapack_interfaces_module

  Interface
     
     Subroutine dsyev( jobz, uplo, n, a, lda, w, work, lwork, info )

       Use numbers_module, Only : wp

       Implicit None

       Character  :: jobz, uplo
       Integer    :: info, lda, lwork, n
       Real( wp ) :: a( lda, * ), w( * ), work( * )

     End Subroutine dsyev

  End Interface

  Public :: dsyev

  Private
  
End Module lapack_interfaces_module

Module outer_product_eigen_module

  Use numbers_module, Only : wp

  Implicit None
  
  Public :: outer_product_eigen

  Private

Contains

  Subroutine outer_product_eigen( v, evals, evecs )

    ! Returns the evals and evecs of the symmetirx matrix formed by taking the outer product of v with itself

    Use numbers_module, Only : wp
    
    Implicit None

    Real( wp ), Dimension( :    ),              Intent( In    ) :: v
    Real( wp ), Dimension( :    ), Allocatable, Intent(   Out ) :: evals
    Real( wp ), Dimension( :, : ), Allocatable, Intent(   Out ) :: evecs

    Real( wp ) :: swap
    Real( wp ) :: dot
    
    Integer :: n
    Integer :: i, j, k

    n = Size( v )

    ! The eigenvalues
    Allocate( evals( 1:n ) )
    evals( 1:n - 1 ) = 0.0_wp
    evals(   n     ) = Dot_product( v, v )

    ! The eigenvectors - generate by modified Gram-Schmidt making orthogonal to
    ! the original vector
    Allocate( evecs( 1:n, 1:n ) )

    ! To tax my little brain less initially stick the normalized original vector in the first column. Will swap with
    ! the last column once we have generated all the vectors.
    evecs( :, 1 ) = v / Norm2( v )

    ! Now generate the other vectors by making them orthonormal to the previous vectors
    ! Use a simple modified Gram-Shmidt for this. There must be a better mehtod.
    Do i = 2, n
       ! Inital "guess" at vector - could this fail in a pathological case????
       evecs( :, i ) = 0.0_wp
       evecs( i, i ) = 1.0_wp
       ! Project out the previous vectors
       Do j = 1, i - 1
!!$          evecs( :, i ) = evecs( :, i ) - Dot_product( evecs( :, j ), evecs( :, i ) ) * evecs( :, j )
          dot = 0.0_wp
          Do k = 1, n
             dot = dot + evecs( k, j ) * evecs( k, i )
          End Do
          Do k = 1, n
             evecs( k, i ) = evecs( k, i ) - dot * evecs( k, j )
          End Do
       End Do
       ! And normalise
       evecs( :, i ) = evecs( :, i ) / Norm2( evecs( :, i ) )
    End Do

    ! And swap the columns so the evecs are in the same order as the evals
    Do i = 1, n
       swap          = evecs( i, 1 )
       evecs( i, 1 ) = evecs( i, n )
       evecs( i, n ) = swap
    End Do
    
  End Subroutine outer_product_eigen

End Module outer_product_eigen_module

Program testit

  Use, Intrinsic :: iso_fortran_env, Only : stdio => output_unit
  
  Use numbers_module            , Only : wp
  Use outer_product_eigen_module, Only : outer_product_eigen
  Use lapack_interfaces_module  , Only : dsyev
  
  Implicit None

  Real( wp ), Dimension( :, : ), Allocatable :: m
  Real( wp ), Dimension( :, : ), Allocatable :: evecs

  Real( wp ), Dimension( : ), Allocatable :: v
  Real( wp ), Dimension( : ), Allocatable :: evals
  Real( wp ), Dimension( : ), Allocatable :: work

  Real( wp ), Dimension( 1:1 ) :: size_work
  
  Integer :: n
  Integer :: info
  Integer :: start, finish, rate

  ! Size of problem
  Write( stdio, '( "Order of problem?" )' )
  Read ( *, * ) n
  Write( stdio, * )

  ! Insitalise v
  Allocate( v( 1:n ) )
  Call Random_Number( v )
  v = v - 0.5_wp

  ! Solve via method in https://math.stackexchange.com/questions/403614/eigenvalues-of-outer-product-matrix-of-two-n-dimensional-vectors
  Call system_clock( start, rate )
  Call outer_product_eigen( v, evals, evecs )
  Call system_clock( finish, rate )
  Call check_results( 'Outer eigen', v, evals, evecs )
  Write( stdio, '( "Time for outer eigen: ", f0.3, " seconds" )' ) Real( finish - start ) / rate
  Write( stdio, * )
  
  ! Find the results via LAPACK
  Deallocate( evals )
  Call system_clock( start, rate )
  Call generate_matrix( v, m )
  ! LAPACK overwtites the matrix with the evecs
  evecs = m
  Allocate( evals( 1:n ) )
  Call dsyev( 'V', 'U', n, evecs, Size( evecs, Dim = 1 ), evals, size_work, -1, info )
  Allocate( work( 1:Nint( size_work( 1 ) ) ) )
  Call dsyev( 'V', 'U', n, evecs, Size( evecs, Dim = 1 ), evals, work, Size( work ), info )
  Call system_clock( finish, rate )
  Write( stdio, '( "LAPACK info = ", i0 )' ) info

  If( info == 0 ) Then
     Call check_results( 'LAPACK', v, evals, evecs )
     Write( stdio, '( "Time for LAPACK: ", f0.3, " seconds" )' ) Real( finish - start ) / rate
  End If

Contains

  Pure Subroutine generate_matrix( v, m )

    ! Generate the matrix M = V .outer_product. V

    Use numbers_module, Only : wp
    
    Implicit None

    Real( wp ), Dimension( :    )             , Intent( In    ) :: v
    Real( wp ), Dimension( :, : ), Allocatable, Intent(   Out ) :: m

    Integer :: n
    Integer :: i, j

    n = Size( v )
    Allocate( m( 1:n, 1:n ) )

    Do j = 1, n
       Do i = 1, n
          m( i, j ) = v( i ) * v( j )
       End Do
    End Do

  End Subroutine generate_matrix

  Subroutine check_results( method, v, evals, evecs )

    ! Check the results of the given method

    Use, Intrinsic :: iso_fortran_env, Only : stdio => output_unit

    Use numbers_module, Only : wp

    Implicit None

    Character( Len = * )         , Intent( In ) :: method
    Real( wp ), Dimension( :    ), Intent( In ) :: v
    Real( wp ), Dimension( :    ), Intent( In ) :: evals
    Real( wp ), Dimension( :, : ), Intent( In ) :: evecs

    Real( wp ), Dimension( :, : ), Allocatable :: m
    Real( wp ), Dimension( :, : ), Allocatable :: tmp

    Integer :: i
    
    Write( stdio, '( "Checking results from ", a )' ) method

    ! Check orthonormality
    tmp = Matmul( Transpose( evecs ), evecs )
    Do i = 1, Size( tmp, Dim = 1 )
       tmp( i, i ) = tmp( i, i ) - 1.0_wp
    End Do
    Write( stdio, '( "Maximum error in orthonormality: ", g20.12 )' ) Maxval( Abs( tmp ) )

    ! Check is solution to eval problem
    Call generate_matrix( v, m )
    tmp = Matmul( Transpose( evecs ), Matmul( m, evecs ) )
    Do i = 1, Size( tmp, Dim = 1 )
       tmp( i, i ) = tmp( i, i ) - evals( i )
    End Do
    Write( stdio, '( "Maximum error in solution      : ", g20.12 )' ) Maxval( Abs( tmp ) )

  End Subroutine check_results

End Program testit
ijb@ijb-Latitude-5410:~/work/stack$ gfortran  -std=f2018 -Wall -Wextra -fcheck=all -O -g -Wimplicit-procedure -Wuse-without-only outer_eigen.f90 -lopenblas
ijb@ijb-Latitude-5410:~/work/stack$ ./a.out
Order of problem?
250

Checking results from Outer eigen
Maximum error in orthonormality:   0.598794597221E-14
Maximum error in solution      :   0.881524230182E-13
Time for outer eigen: .051 seconds

LAPACK info = 0
Checking results from LAPACK
Maximum error in orthonormality:   0.777156117238E-14
Maximum error in solution      :   0.112331207400E-13
Time for LAPACK: .020 seconds
ijb@ijb-Latitude-5410:~/work/stack$ gfortran  -std=f2018 -Wall -Wextra  -Ofast -g -Wimplicit-procedure -Wuse-without-only outer_eigen.f90 -lopenblas
ijb@ijb-Latitude-5410:~/work/stack$ ./a.out
Order of problem?
250

Checking results from Outer eigen
Maximum error in orthonormality:   0.119649769260E-13
Maximum error in solution      :   0.215376845007E-12
Time for outer eigen: .011 seconds

LAPACK info = 0
Checking results from LAPACK
Maximum error in orthonormality:   0.976996261670E-14
Maximum error in solution      :   0.982076900499E-14
Time for LAPACK: .021 seconds
ijb@ijb-Latitude-5410:~/work/stack$ ./a.out
Order of problem?
500

Checking results from Outer eigen
Maximum error in orthonormality:   0.115890803898E-12
Maximum error in solution      :   0.469814678927E-11
Time for outer eigen: .088 seconds

LAPACK info = 0
Checking results from LAPACK
Maximum error in orthonormality:   0.888178419700E-14
Maximum error in solution      :   0.410210860282E-13
Time for LAPACK: .108 seconds
ijb@ijb-Latitude-5410:~/work/stack$ ./a.out
Order of problem?
1000

Checking results from Outer eigen
Maximum error in orthonormality:   0.357838758624E-13
Maximum error in solution      :   0.291765283970E-11
Time for outer eigen: .406 seconds

LAPACK info = 0
Checking results from LAPACK
Maximum error in orthonormality:   0.150990331349E-13
Maximum error in solution      :   0.124051162678E-12
Time for LAPACK: .925 seconds
ijb@ijb-Latitude-5410:~/work/stack$ gfortran  -std=f2018 -Wall -Wextra  -Ofast -g -Wimplicit-procedure -Wuse-without-only outer_eigen.f90 -lopenblas
ijb@ijb-Latitude-5410:~/work/stack$ ./a.out
Order of problem?
2000

Checking results from Outer eigen
Maximum error in orthonormality:   0.866971425206E-14
Maximum error in solution      :   0.147852563526E-11
Time for outer eigen: 4.076 seconds

LAPACK info = 0
Checking results from LAPACK
Maximum error in orthonormality:   0.264233079861E-13
Maximum error in solution      :   0.373884284743E-12
Time for LAPACK: 8.918 seconds
ijb@ijb-Latitude-5410:~/work/stack$ gfortran --version
GNU Fortran (Ubuntu 9.4.0-1ubuntu1~20.04.1) 9.4.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

ijb@ijb-Latitude-5410:~/work/stack$

Upvotes: 4

Compute the eigenvalues of v ⊗ v

Answers (1)

Related Questions