Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.forth > #60572 > unrolled thread

fast matrix determinant

Back to article view | Back to comp.lang.forth

Contents

  fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-10 10:29 -0700
    Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-12 05:41 -0700
      Re: fast matrix determinant mhx@iae.nl - 2017-06-12 06:03 -0700
        Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-12 07:14 -0700
          Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-12 07:28 -0700
          Re: fast matrix determinant mhx@iae.nl - 2017-06-12 11:21 -0700
            Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-12 16:39 -0700
          Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-14 18:05 +0300
            Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-14 11:05 -0700
              Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-14 21:37 +0300
            Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-15 09:59 +0300
              Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-15 03:30 -0700
                Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-15 17:57 +0300
                  Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-15 19:06 -0700
                    Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-16 10:21 +0300
                      Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-16 03:38 -0700
                        Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-16 15:58 +0300
                          Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-16 06:04 -0700
                            Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-16 16:19 +0300
                      Re: fast matrix determinant mhx@iae.nl - 2017-06-16 03:58 -0700
                        Re: fast matrix determinant albert@cherry.spenarnc.xs4all.nl (Albert van der Horst) - 2017-06-16 13:39 +0200
                          Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-16 05:58 -0700
                        Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-16 05:12 -0700
                          Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-16 05:40 -0700
                            Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-16 17:50 +0300
                              Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-16 13:39 -0700
                                Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-17 10:01 +0300
                            Re: fast matrix determinant mhx@iae.nl - 2017-06-17 01:31 -0700
                              Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-17 15:24 -0700
              Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-15 17:59 +0300
            Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-16 10:11 +0300
              Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-17 09:20 +0300
        Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-13 16:43 -0700
          Re: fast matrix determinant mhx@iae.nl - 2017-06-13 22:04 -0700
          Re: fast matrix determinant mhx@iae.nl - 2017-06-14 12:23 -0700
            Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-14 18:31 -0700
      Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-15 19:13 -0700

Page 1 of 2  [1] 2  Next page →

#60572 — fast matrix determinant

Our recent discussion about fast matrix multiplication suggests to me the need for fast code to compute matrix determinants. The FSL provides two modules, dets.x and lufact.x, which may be used for computing the determinant of a square NxN matrix: the LU factorization of the matrix is first performed and then the determinant is computed. This approach requires only O(N^3) floating point operations, compared to O(N!) for the standard way we learn to compute determinants in school; see

http://www2.hawaii.edu/~norbert/CompPhys/compphysicsch10.html

If we simply want the determinant of an NxN matrix, then the LU factorization may be doing more than the minimum needed operations (I'm not sure about this, yet). It's also not clear to me that putting a square matrix in upper-triangular form has the capability to be transformed into a parallel computation, easily. Besides the FSL algorithm for determinants, I have one other Forth program, determ.x, translated from a Fortran subroutine. These may be found at:

ftp://ccreweb.org/software/fsl/lufact.fs
ftp://ccreweb.org/software/fsl/dets.fs
ftp://ccreweb.org/software/fsl/extras/determ.fs

Are there some other Forth programs for fast computation of determinants, and benchmarks to go with them?

Krishna

[toc] | [next] | [standalone]

#60653

On Saturday, June 10, 2017 at 12:29:13 PM UTC-5, krishna...org wrote:
> ... Besides the FSL algorithm for determinants, I have one other Forth program, determ.x, translated from a Fortran subroutine. These may be found at:
> 
> ftp://ccreweb.org/software/fsl/lufact.fs
> ftp://ccreweb.org/software/fsl/dets.fs
> ftp://ccreweb.org/software/fsl/extras/determ.fs
> 

A properly factored version of kForth's fsl/extras/determ.4th will allow for optimization, both in the Forth source and by selective replacement with CODE word(s). A revised version, called determ2.4th, is given below. It is considerably faster, even though it doesn't use any CODE words, yet.

Krishna
-------------------

\ Test original, unfactored determ (from determ.4th)
\ on 100x100 random matrix
$ kforth-fast determ-test

\ ... 
Initializing the matrix ... 
The matrix A is 100 x 100 

Computing determinant. -4.72105e+24 
Elapsed [ms]: 781  ok

\ Test new, factored determ (from determ2.4th) on same matrix
$ kforth-fast determ-test
\ ...
Initializing the matrix ... 
The matrix A is 100 x 100 

Computing determinant. -4.72105e+24 
Elapsed [ms]: 465  ok


------------------- Test code
\ determ-test.4th
\
\ Krishna Myneni, 2017-06-12
\

include ans-words
include modules
include fsl/fsl-util
include fsl/horner
include fsl/extras/noise
\ include fsl/extras/determ
include fsl/extras/determ2

100 constant NCOLS
100 constant NROWS

NROWS NCOLS DFLOAT MATRIX a{{ 

: init-matrix ( -- )
    A{{ 0 0 }} NCOLS NROWS * 0 DO  ran0 2 pick F!  dfloat+ LOOP drop
;

cr .( Initializing the matrix ... )
init-matrix
cr .( The matrix A is ) NROWS . .( x ) NCOLS . 
cr cr .( Computing determinant. )
: test ( -- ms)  ms@ a{{ Nrows determ ms@ >r f. r> swap - ;
test cr .( Elapsed [ms]: ) .

------------------- determ2.4th
\
\ determ2.4th
\
\ Return the determinant of a real square matrix.
\ The input matrix is destroyed in the process.

\ Krishna Myneni, Creative Consulting for Research and Education
\ http://ccreweb.org
\
\ Revisions: 
\ 
\   2017-06-11 km; rewrote determ.4th, factoring code for
\                  clarity of algorithm and improving execution
\                  efficiency.
\
\ determ is based on a similar routine from P.R. Bevington,
\   "Data Reduction and Error Analysis for the Physical Sciences",
\   1969, McGraw-Hill.
\

Begin-Module

0 value L
0 value Norder		\ order of matrix 
0 ptr arr{{		\ address of matrix ( ptr is same as VALUE in ANS-Forth)
variable roffs

Public:

: set-determ-params ( 'a N -- )
    to Norder  to arr{{
    Norder dfloats roffs !
;

\ find column with next non-zero element, starting at
\ row = col; column = col+1; return the column number, col',
\ with non-zero element, which will be equal to Norder
\ if no non-zero elements found.
: df_next_nonzero ( col -- col' )
    arr{{ swap 
    dup to L
    BEGIN
      2dup L }} F@ F0= 
      L Norder 1- < and
    WHILE
      L 1+ to L
    REPEAT
    2drop L
;

\ exchange the double floats at a1 and a2
: df_xchg ( a1 a2 -- )
    dup pad dfloat move
    2dup dfloat move
    drop pad swap dfloat move    
;

\ swap columns n and m of matrix
: df_swap_cols ( n m -- )
    2>r
    arr{{ 0 r> }} arr{{ 0 r> }}
    Norder 0 DO
      2dup df_xchg
      roffs @ + swap roffs @ +
    LOOP
    2drop
;

\ subtract scaled row L from rows below to put matrix
\ in diagonal form (only does needed parts of rows).
: df_sub_rows ( l -- )
    to L
    arr{{ L L }} f@
    Norder L 1+ DO
      arr{{ I L }} f@ fover f/ 
      Norder L 1+ DO
	arr{{ L I }} f@ fover f* fnegate
	arr{{ J I }} dup >r f@ f+ r> f!
      LOOP
      fdrop
    LOOP
    fdrop
;

: determ ( 'a n -- fdet | a is the matrix, n is its order )
    set-determ-params
    1e   
    Norder 0 DO
      arr{{ I I }} F@  F0= IF  \ zero on diagonal
	\ Find next element in row which is non-zero
	I df_next_nonzero
	dup Norder = IF drop fdrop 0e LEAVE  THEN
	to L
	arr{{ I L }} F@  F0= IF  fdrop 0e LEAVE  THEN
        L I df_swap_cols
	fnegate
      THEN
	
      \ Cumulative product of diagonal elements
      arr{{ I I }} F@ F*

      \ Subtract scaled row I from lower rows
      I Norder 1- < IF  I df_sub_rows  THEN
    LOOP
;

End-Module

TEST-CODE? [IF]   \ test code ==============================================
[undefined] T{      [IF]  s" ttester.4th" included  [THEN]
BASE @
DECIMAL

1e-15 rel-near F!
1e-15 abs-near F!
set-near

2 2 FLOAT MATRIX m2{{
    3e  1e
    5e  2e
2 2 m2{{ }}fput

3 3 FLOAT MATRIX m3{{
    1e  2e  3e
    2e  1e  1e
    3e  1e  2e
3 3 m3{{ }}fput

4 4 FLOAT MATRIX m4{{
    2e  1e  4e  2e
    0e  2e  1e  2e
    0e  2e  1e  1e
    2e  0e  1e  0e
4 4 m4{{ }}fput

TESTING DETERM
t{ m2{{ 2 determ  ->   1e  r}t
t{ m3{{ 3 determ  ->  -4e  r}t
t{ m4{{ 4 determ  -> -10e  r}t  

BASE !
[THEN]

[toc] | [prev] | [next] | [standalone]

#60654

On Monday, June 12, 2017 at 2:41:29 PM UTC+2, krishna...@ccreweb.org wrote:
> On Saturday, June 10, 2017 at 12:29:13 PM UTC-5, krishna...org wrote:
> > ... Besides the FSL algorithm for determinants, I have one other Forth program, determ.x, translated from a Fortran subroutine. These may be found at:
> > 
> > ftp://ccreweb.org/software/fsl/lufact.fs
> > ftp://ccreweb.org/software/fsl/dets.fs
> > ftp://ccreweb.org/software/fsl/extras/determ.fs
> > 
> 
> A properly factored version of kForth's fsl/extras/
> determ.4th will allow for optimization, both in the
> Forth source and by selective replacement with CODE
> word(s). A revised version, called determ2.4th, is 
> given below. It is considerably faster, even though
> it doesn't use any CODE words, yet.

Hi Krishna,

Just curious: what do you need determinants for? I've
occasional use for them for displaying information,
but never in a tight loop?

In general, apart from ddot and dgemm, eigenvalue and
eigenvector decomposers are the most useful routines
to have. In fact, I have been looking for a PLACE
(Matlab) replacement for years. SLICOT has it in
(Fortran) source, but I didn't get to converting 
(or linking) it yet.

-marcel

[toc] | [prev] | [next] | [standalone]

#60657

Hi Marcel,

I need a fast determinant to write some pedagogical code
for illustrating a molecular physics calculation. It is
to be used for solving the eigenvalue problem,

    D( H - lambda×S ) = 0

where H and S are N×N matrices, and lambda(s) are the
eigenvalue(s) to be found. D() is, of course, the
determinant operator. This form of the eigenvalue
problem is typical in molecular physics calculations,
when the basis set of eigenfunctions are non-orthogonal.

The matrix S, consists of the inner products, of the N
eigenfunctions, and it would simply be the identity
matrix if the N eigenfunctions were orthonormal. 
Typically, the hard part in setting up the problem is
evaluating the elements of H and S, each of the elements
requiring the evaluation of a multi-dimensional integral.
Nevertheless, when N becomes large, efficient evaluation
of the determinants is needed, particularly in searching
for the eigenvalues.

A simple way of solving the equation, once the matrix
elements of H and S have been computed in the basis set,
is to simply use the determinant equation within a root
finder, to search for the values of lambda which satisfy
the equation. This is conceptually simple and straight-
forward to implement, but it's probably not the most
computationally efficient method -- there are probably
much better techniques available in a linear algebra
library. But, for pedagogical use, it's nice. Typically,
we are interested in only the lowest eigenvalue.

Krishna

[toc] | [prev] | [next] | [standalone]

#60659

In my message above, please replace "eigenfunction" with
"basis function". The eigenfunctions are special linear
combinations of the basis functions, and they are not
part of the problem being solved here.

Krishna

[toc] | [prev] | [next] | [standalone]

#60666

On Monday, June 12, 2017 at 4:14:07 PM UTC+2, krishna...@ccreweb.org wrote:
> Hi Marcel,
> 
> I need a fast determinant to write some pedagogical code
> for illustrating a molecular physics calculation. It is
> to be used for solving the eigenvalue problem,
> 
>     D( H - lambda×S ) = 0
> 
> where H and S are N×N matrices, and lambda(s) are the
> eigenvalue(s) to be found. D() is, of course, the
> determinant operator. This form of the eigenvalue
> problem is typical in molecular physics calculations,
> when the basis set of eigenfunctions are non-orthogonal.
> 
> The matrix S, consists of the inner products, of the N
> eigenfunctions, and it would simply be the identity
> matrix if the N eigenfunctions were orthonormal. 
> Typically, the hard part in setting up the problem is
> evaluating the elements of H and S, each of the elements
> requiring the evaluation of a multi-dimensional integral.
> Nevertheless, when N becomes large, efficient evaluation
> of the determinants is needed, particularly in searching
> for the eigenvalues.
> 
> A simple way of solving the equation, once the matrix
> elements of H and S have been computed in the basis set,
> is to simply use the determinant equation within a root
> finder, to search for the values of lambda which satisfy
> the equation. This is conceptually simple and straight-
> forward to implement, but it's probably not the most
> computationally efficient method -- there are probably
> much better techniques available in a linear algebra
> library. But, for pedagogical use, it's nice. Typically,
> we are interested in only the lowest eigenvalue.

Looks like you are interested in the eigenvalues/vectors,
not in the determinants as such?

-marcel

[toc] | [prev] | [next] | [standalone]

#60671

On Monday, June 12, 2017 at 1:21:04 PM UTC-5, m...@iae.nl wrote:
> On Monday, June 12, 2017 at 4:14:07 PM UTC+2, krishna...org wrote:
...
> > for illustrating a molecular physics calculation. It is
> > to be used for solving the eigenvalue problem,
> > 
> >     D( H - lambda×S ) = 0
> > 
> > where H and S are N×N matrices, and lambda(s) are the
> > eigenvalue(s) to be found. D() is, of course, the
> > determinant operator. ...
> 
> Looks like you are interested in the eigenvalues/vectors,
> not in the determinants as such?
> 

Well, that's essentially true, but we should remember that eigenvalues of a matrix are arrived at by a determinant equation, which provides the characteristic equation for which the eigenvalues are roots. Determinants are of fundamental importance in the theory of eigenvalues; it's just that the connection is not always obvious when using a canned library package.

A typical way of obtaining eigenvalues is through transformations of a square matrix into tridiagonal form, for which the determinant is easier to write, and hence the characteristic equation for finding the eigenvalues. However, such a method assumes we are solving the typical eigenvalue equation,

   D( H - lambda*I ) = 0

where I is the identity matrix. The more general form replaces I by S, for a non-orthonormal basis representation. One can use Gram-Schmidt orthogonalization to transform the more general form into the form above, but then a similarity transformation has to be performed on H before using the standard methods. Also, if we change from N to N+1 in the order of the matrix, the entire H matrix has to be recomputed, and not just an additional row and column, as for solving D( H - lambda*S ) = 0. 

For QM calculations in which the representation of a wavefunction (eigenfunction) is being successively approximated by the addition of basis functions, use of the method of computing determinants within a root finder shows the gradual refinement of the eigenvalues, without going through the opaque steps of orthogonalization and similarity transformations. However, I'm not sure about the numerical stability of directly solving the determinant equation for eigenvalues.

Another place that determinants come up in QM is for obtaining a multi-particle wavefunction (for identical particles) with the proper symmetry under exchange of coordinates.

At present we lack useful Forth source libraries for finding eigenvalues and eigenvectors of matrices. While it's possible to call external libraries, it's just more convenient to simply load a source file. Optimizations can be performed as needed for the scale of the problem.

Krishna

[toc] | [prev] | [next] | [standalone]

#60728

On 12.06.2017 17:14, krishna.myneni@ccreweb.org wrote:
> Hi Marcel,
> 
> I need a fast determinant to write some pedagogical code
> for illustrating a molecular physics calculation. It is
> to be used for solving the eigenvalue problem,
> 
>      D( H - lambda×S ) = 0
> 
> where H and S are N×N matrices, and lambda(s) are the
> eigenvalue(s) to be found. D() is, of course, the
> determinant operator. This form of the eigenvalue
> problem is typical in molecular physics calculations,
> when the basis set of eigenfunctions are non-orthogonal.
> 
> The matrix S, consists of the inner products, of the N
> eigenfunctions, and it would simply be the identity
> matrix if the N eigenfunctions were orthonormal.
> Typically, the hard part in setting up the problem is
> evaluating the elements of H and S, each of the elements
> requiring the evaluation of a multi-dimensional integral.
> Nevertheless, when N becomes large, efficient evaluation
> of the determinants is needed, particularly in searching
> for the eigenvalues.
> 
> A simple way of solving the equation, once the matrix
> elements of H and S have been computed in the basis set,
> is to simply use the determinant equation within a root
> finder, to search for the values of lambda which satisfy
> the equation. This is conceptually simple and straight-
> forward to implement, but it's probably not the most
> computationally efficient method -- there are probably
> much better techniques available in a linear algebra
> library. But, for pedagogical use, it's nice. Typically,
> we are interested in only the lowest eigenvalue.
> 
> Krishna
> 

Recently, after your post, I have found something that
could be useful:


A symbol between '[' ']' represents a matrix,
between '{' '}' represents vector,
^T   means transpose
^-1  means inverse of matrix.

Theorem (proposition):

Be a real symmetric matrix [A] NxN  with N distinct nonzero eigenvalues
and [C] its transformation matrix into orthonormal basis of [A]
eigenvectors.

Then, for very small perturbation [dA] (Norm([dA] << Norm([A]))
will result delta eigenvalues dL, for every index k, 1 <= k <= N

dLk = Bkk  ,where [B] = [C]^T * [dA] * [C]



<Proof/

[A] is a symmetric matrix NxN, with N distinct eigenvalues.
That means:

      1) eigenvectors {x}i, 1 <= i <= N are linear independents

      2) for every i j,   1 <= i != j <=N,   {x}i^T * {x}j = 0.

---
3)   [A]*{x}=L*(x).

For a very small perturbation [dA] will follows:

[A + dA] * {x + dx} = (L + dl) * {x + dx} <=>

[A]*{x} + [A]*{dx} + [dA]*{x} + [dA]*{dx}
	=
L*{x} + L*{dx} + dL*{x} + dL*{dx}

using 3)   =>

[A]*{dx} + [dA]*{x} + [dA]*{dx}
	 =
L*{dx} + dL*{x} + dL*{dx}

For Norm([dA]) << Norm([A]) , second order terms are negligible, so

[A]*{dx} + [dA]*{x} = L*{dx} + dL*{x} <=>

([A] - L*[I])*{dx} = ([I]*dL - [dA])*{x} <=>

For every 1 <= k <= N:

4)   ([A] -Lk*[I]) * {dx}k = ([I]*dLk - [dA]) * {x}k

Using 1), vector {dx}k can be expressed in terms of eigenvectors {x}i, so

{dx}k  =  SUM(Fik * {x}i) =>

([A] - Lk*[I]) * {dx}k = ([A] - Lk*[I]) * SUM(Fik * {x}i)
      =
SUM(Fik * ([A] - Lk*[I])*{x}i ) =
SUM(Fik * ([A]*{x}i - Lk*[I]*{x}i) ) =
SUM(Fik * (Li*{x}i  - Lk*{x}i) ) =
SUM(Fik * (Li - Lk)*{x}i )

So 4)   becomes

SUM(Fik * (Li - Lk)*{x}i ) = ([I]*dLk - [dA]) * {x}k

multiplying at left side of equality with Transpose of {x}k results

{x}k^T * SUM(Fik * (Li - Lk)*{x}i  = {x}k^T * ([I]*dLk - [dA]) * {x}k   <=>

SUM(Fik * (Li -Lk) * {x}k^T * {x}i = {x}k^T * ([I]*dLk - [dA]) * {x}k

 From 2) results that for i != k  {x}k^T * {x}i = 0,
also for i=k,  Li - Lk = 0, so left side of equality becames 0:

0 = {x}k^T * ([I]*dLk - [dA]) * {x}k

if we note unity eigenvector {v}k = {x}k/Norm({x}k):

0 = {v}k^T *  ([I]*dLk - [dA]) * {v}k


Transforming bilinear form into orthonormal basis of eigenvectors,
equality becames:

5)   0 = {e}k^T * ([C]^T * ([I]*dLk - [dA]) * [C]) * {e}k

    ,with {e}k versor of orthonormal eigenvalues basis.

Let's note:

[D] = [C]^T * ([I]*dLk - [dA]) * [C]  <=>

[D] = [C]^T*[I]*[C] * dLk  -  [C]^T*[dA]*[C]  <=>

[D] = [I]*dLk - [C]^T*[dA]*[C]


Let's note  [B] = [C]^T*[dA]*[C]  ,so


[D] = [I]*dLk - [B]  ,so 5)   becames:


0 = {e}k^T * ([I]*dLk - [B]) * {e}k  <=>

0 = {e}k^T * [I] * {e}k * dLk  -  {e}k^T * [B] * {e}k  <=>

0= dLk - Bkk <=>

*************
* dLk = Bkk *
*************
/Proof>

Quickly  tested for matrix

0 1
1 0

with small perturbation [dA]

a11 a12
a21 a22

,aij << 1

and it holds :)


Have a nice day,
hum9

[toc] | [prev] | [next] | [standalone]

#60749

On Wednesday, June 14, 2017 at 10:05:08 AM UTC-5, hum9 wrote:
> On 12.06.2017 17:14, krishna...org wrote:
...
> Recently, after your post, I have found something that
> could be useful:
> 
> 
> A symbol between '[' ']' represents a matrix,
> between '{' '}' represents vector,
> ^T   means transpose
> ^-1  means inverse of matrix.
> 
> Theorem (proposition):
> 
> Be a real symmetric matrix [A] NxN  with N distinct nonzero eigenvalues
> and [C] its transformation matrix into orthonormal basis of [A]
> eigenvectors.
> 
> Then, for very small perturbation [dA] (Norm([dA] << Norm([A]))
> will result delta eigenvalues dL, for every index k, 1 <= k <= N
> 
> dLk = Bkk  ,where [B] = [C]^T * [dA] * [C]
> 
> ...

That's an interesting result, but I'm not sure how it can be applied to the problem I described. I don't understand what your transformation matrix, [C], is supposed to do? Does it transform basis vectors into a new basis set? If so, the transformation,

   [B] = [C]^T * [dA] * [C]

is like a similarity transformation, for the special case of

   [C]^-1 = [C]^T

I haven't examined your proof in detail, but if one can obtain a first order approximation to the new eigenvalues, based on perturbations (small changes) to the elements of [A], it seems useful.

In my problem description, the matrix [H] ([A] in your description) is not being perturbed, i.e. small changes are not being made to the matrix elements, but it's rank is being increased by 1, in successive approximations.

Krishna

[toc] | [prev] | [next] | [standalone]

#60755

On 14.06.2017 21:05, krishna.myneni@ccreweb.org wrote:
> On Wednesday, June 14, 2017 at 10:05:08 AM UTC-5, hum9 wrote:
>> On 12.06.2017 17:14, krishna...org wrote:
> ...
>> Recently, after your post, I have found something that
>> could be useful:
>>
>>
>> A symbol between '[' ']' represents a matrix,
>> between '{' '}' represents vector,
>> ^T   means transpose
>> ^-1  means inverse of matrix.
>>
>> Theorem (proposition):
>>
>> Be a real symmetric matrix [A] NxN  with N distinct nonzero eigenvalues
>> and [C] its transformation matrix into orthonormal basis of [A]
>> eigenvectors.
>>
>> Then, for very small perturbation [dA] (Norm([dA] << Norm([A]))
>> will result delta eigenvalues dL, for every index k, 1 <= k <= N
>>
>> dLk = Bkk  ,where [B] = [C]^T * [dA] * [C]
>>
>> ...
> 
> That's an interesting result, but I'm not sure how it can be applied to the problem I described. I don't understand what your transformation matrix, [C], is supposed to do? Does it transform basis vectors into a new basis set? If so, the transformation,
> 
>     [B] = [C]^T * [dA] * [C]
> 
> is like a similarity transformation, for the special case of
> 
>     [C]^-1 = [C]^T
> 
> I haven't examined your proof in detail, but if one can obtain a first order approximation to the new eigenvalues, based on perturbations (small changes) to the elements of [A], it seems useful.
> 
> In my problem description, the matrix [H] ([A] in your description) is not being perturbed, i.e. small changes are not being made to the matrix elements, but it's rank is being increased by 1, in successive approximations.
> 
> Krishna
> 

[C]^T has eigenversors {v} of [A] as rows , [C] has them as columns,

     as {v}i^T * {v}j = 1 for i=j, 0 for i !=j ,

so product [C]^T * [C] = [I]. Yes, [C]^T = [C]^-1 in this case.

I think it does not have direct relation with your problem, but could
be useful for approximating solutions.

Have a nice day,
hum9

[toc] | [prev] | [next] | [standalone]

#60773

On 14.06.2017 18:05, hum9 wrote:
> On 12.06.2017 17:14, krishna.myneni@ccreweb.org wrote:
>> Hi Marcel,
>>
>> I need a fast determinant to write some pedagogical code
>> for illustrating a molecular physics calculation. It is
>> to be used for solving the eigenvalue problem,
>>
>>      D( H - lambda×S ) = 0
>>
>> where H and S are N×N matrices, and lambda(s) are the
>> eigenvalue(s) to be found. D() is, of course, the
>> determinant operator. This form of the eigenvalue
>> problem is typical in molecular physics calculations,
>> when the basis set of eigenfunctions are non-orthogonal.
>>
>> The matrix S, consists of the inner products, of the N
>> eigenfunctions, and it would simply be the identity
>> matrix if the N eigenfunctions were orthonormal.
>> Typically, the hard part in setting up the problem is
>> evaluating the elements of H and S, each of the elements
>> requiring the evaluation of a multi-dimensional integral.
>> Nevertheless, when N becomes large, efficient evaluation
>> of the determinants is needed, particularly in searching
>> for the eigenvalues.
>>
>> A simple way of solving the equation, once the matrix
>> elements of H and S have been computed in the basis set,
>> is to simply use the determinant equation within a root
>> finder, to search for the values of lambda which satisfy
>> the equation. This is conceptually simple and straight-
>> forward to implement, but it's probably not the most
>> computationally efficient method -- there are probably
>> much better techniques available in a linear algebra
>> library. But, for pedagogical use, it's nice. Typically,
>> we are interested in only the lowest eigenvalue.
>>
>> Krishna
>>
> 
> Recently, after your post, I have found something that
> could be useful:
> 
> 
> A symbol between '[' ']' represents a matrix,
> between '{' '}' represents vector,
> ^T   means transpose
> ^-1  means inverse of matrix.
> 
> Theorem (proposition):
> 
> Be a real symmetric matrix [A] NxN  with N distinct nonzero eigenvalues
> and [C] its transformation matrix into orthonormal basis of [A]
> eigenvectors.
> 
> Then, for very small perturbation [dA] (Norm([dA] << Norm([A]))
> will result delta eigenvalues dL, for every index k, 1 <= k <= N
> 
> dLk = Bkk  ,where [B] = [C]^T * [dA] * [C]
> 
> 
> 
> <Proof/
> 
> [A] is a symmetric matrix NxN, with N distinct eigenvalues.
> That means:
> 
>       1) eigenvectors {x}i, 1 <= i <= N are linear independents
> 
>       2) for every i j,   1 <= i != j <=N,   {x}i^T * {x}j = 0.
> 
> ---
> 3)   [A]*{x}=L*(x).
> 
> For a very small perturbation [dA] will follows:
> 
> [A + dA] * {x + dx} = (L + dl) * {x + dx} <=>
> 
> [A]*{x} + [A]*{dx} + [dA]*{x} + [dA]*{dx}
>      =
> L*{x} + L*{dx} + dL*{x} + dL*{dx}
> 
> using 3)   =>
> 
> [A]*{dx} + [dA]*{x} + [dA]*{dx}
>       =
> L*{dx} + dL*{x} + dL*{dx}
> 
> For Norm([dA]) << Norm([A]) , second order terms are negligible, so
> 
> [A]*{dx} + [dA]*{x} = L*{dx} + dL*{x} <=>
> 
> ([A] - L*[I])*{dx} = ([I]*dL - [dA])*{x} <=>
> 
> For every 1 <= k <= N:
> 
> 4)   ([A] -Lk*[I]) * {dx}k = ([I]*dLk - [dA]) * {x}k
> 
> Using 1), vector {dx}k can be expressed in terms of eigenvectors {x}i, so
> 
> {dx}k  =  SUM(Fik * {x}i) =>
> 
> ([A] - Lk*[I]) * {dx}k = ([A] - Lk*[I]) * SUM(Fik * {x}i)
>       =
> SUM(Fik * ([A] - Lk*[I])*{x}i ) =
> SUM(Fik * ([A]*{x}i - Lk*[I]*{x}i) ) =
> SUM(Fik * (Li*{x}i  - Lk*{x}i) ) =
> SUM(Fik * (Li - Lk)*{x}i )
> 
> So 4)   becomes
> 
> SUM(Fik * (Li - Lk)*{x}i ) = ([I]*dLk - [dA]) * {x}k
> 
> multiplying at left side of equality with Transpose of {x}k results
> 
> {x}k^T * SUM(Fik * (Li - Lk)*{x}i  = {x}k^T * ([I]*dLk - [dA]) * {x}k   <=>
> 
> SUM(Fik * (Li -Lk) * {x}k^T * {x}i = {x}k^T * ([I]*dLk - [dA]) * {x}k
> 
>  From 2) results that for i != k  {x}k^T * {x}i = 0,
> also for i=k,  Li - Lk = 0, so left side of equality becames 0:
> 
> 0 = {x}k^T * ([I]*dLk - [dA]) * {x}k
> 
> if we note unity eigenvector {v}k = {x}k/Norm({x}k):
> 
> 0 = {v}k^T *  ([I]*dLk - [dA]) * {v}k
> 
> 
> Transforming bilinear form into orthonormal basis of eigenvectors,
> equality becames:
> 
> 5)   0 = {e}k^T * ([C]^T * ([I]*dLk - [dA]) * [C]) * {e}k
> 
>     ,with {e}k versor of orthonormal eigenvalues basis.
> 
> Let's note:
> 
> [D] = [C]^T * ([I]*dLk - [dA]) * [C]  <=>
> 
> [D] = [C]^T*[I]*[C] * dLk  -  [C]^T*[dA]*[C]  <=>
> 
> [D] = [I]*dLk - [C]^T*[dA]*[C]
> 
> 
> Let's note  [B] = [C]^T*[dA]*[C]  ,so
> 
> 
> [D] = [I]*dLk - [B]  ,so 5)   becames:
> 
> 
> 0 = {e}k^T * ([I]*dLk - [B]) * {e}k  <=>
> 
> 0 = {e}k^T * [I] * {e}k * dLk  -  {e}k^T * [B] * {e}k  <=>
> 
> 0= dLk - Bkk <=>
> 
> *************
> * dLk = Bkk *
> *************
> /Proof>
> 
> Quickly  tested for matrix
> 
> 0 1
> 1 0
> 
> with small perturbation [dA]
> 
> a11 a12
> a21 a22
> 
> ,aij << 1
> 
> and it holds :)
> 
> 
> Have a nice day,
> hum9

I reformulated this proposed theorem with less constraints:

If real NxN matrix [A] has N orthonormal eigenvectors,
for a small perturbation [dA] corresponds small eigenvalues changes:
	dLi = Bii , where [B] = [C]^T * [dA] * [C],
with matrix [C] of orthonormal eigenvectors as columns.

hum9

[toc] | [prev] | [next] | [standalone]

#60774

On Thursday, June 15, 2017 at 1:59:44 AM UTC-5, hum9 wrote:
> On 14.06.2017 18:05, hum9 wrote:
> > On 12.06.2017 17:14, krishna...org wrote:
> >> Hi Marcel,
> >>
> >> I need a fast determinant to write some pedagogical code
> >> for illustrating a molecular physics calculation. It is
> >> to be used for solving the eigenvalue problem,
> >>
> >>      D( H - lambda×S ) = 0
> >>
> >> where H and S are N×N matrices, and lambda(s) are the
> >> eigenvalue(s) to be found. D() is, of course, the
> >> determinant operator. This form of the eigenvalue
> >> problem is typical in molecular physics calculations,
> >> when the basis set of eigenfunctions are non-orthogonal.
> >>
> >> The matrix S, consists of the inner products, of the N
> >> eigenfunctions, and it would simply be the identity
> >> matrix if the N eigenfunctions were orthonormal.
> >> Typically, the hard part in setting up the problem is
> >> evaluating the elements of H and S, each of the elements
> >> requiring the evaluation of a multi-dimensional integral.
> >> Nevertheless, when N becomes large, efficient evaluation
> >> of the determinants is needed, particularly in searching
> >> for the eigenvalues.
> >>
> >> A simple way of solving the equation, once the matrix
> >> elements of H and S have been computed in the basis set,
> >> is to simply use the determinant equation within a root
> >> finder, to search for the values of lambda which satisfy
> >> the equation. This is conceptually simple and straight-
> >> forward to implement, but it's probably not the most
> >> computationally efficient method -- there are probably
> >> much better techniques available in a linear algebra
> >> library. But, for pedagogical use, it's nice. Typically,
> >> we are interested in only the lowest eigenvalue.
> >>
> >> Krishna
> >>
> > 
> > Recently, after your post, I have found something that
> > could be useful:
> > 
> > 
> > A symbol between '[' ']' represents a matrix,
> > between '{' '}' represents vector,
> > ^T   means transpose
> > ^-1  means inverse of matrix.
> > 
> > Theorem (proposition):
> > 
> > Be a real symmetric matrix [A] NxN  with N distinct nonzero eigenvalues
> > and [C] its transformation matrix into orthonormal basis of [A]
> > eigenvectors.
> > 
> > Then, for very small perturbation [dA] (Norm([dA] << Norm([A]))
> > will result delta eigenvalues dL, for every index k, 1 <= k <= N
> > 
> > dLk = Bkk  ,where [B] = [C]^T * [dA] * [C]
> > 
> > 
> > 
> > <Proof/
> > 
> > [A] is a symmetric matrix NxN, with N distinct eigenvalues.
> > That means:
> > 
> >       1) eigenvectors {x}i, 1 <= i <= N are linear independents
> > 
> >       2) for every i j,   1 <= i != j <=N,   {x}i^T * {x}j = 0.
> > 
> > ---
> > 3)   [A]*{x}=L*(x).
> > 
> > For a very small perturbation [dA] will follows:
> > 
> > [A + dA] * {x + dx} = (L + dl) * {x + dx} <=>
> > 
> > [A]*{x} + [A]*{dx} + [dA]*{x} + [dA]*{dx}
> >      =
> > L*{x} + L*{dx} + dL*{x} + dL*{dx}
> > 
> > using 3)   =>
> > 
> > [A]*{dx} + [dA]*{x} + [dA]*{dx}
> >       =
> > L*{dx} + dL*{x} + dL*{dx}
> > 
> > For Norm([dA]) << Norm([A]) , second order terms are negligible, so
> > 
> > [A]*{dx} + [dA]*{x} = L*{dx} + dL*{x} <=>
> > 
> > ([A] - L*[I])*{dx} = ([I]*dL - [dA])*{x} <=>
> > 
> > For every 1 <= k <= N:
> > 
> > 4)   ([A] -Lk*[I]) * {dx}k = ([I]*dLk - [dA]) * {x}k
> > 
> > Using 1), vector {dx}k can be expressed in terms of eigenvectors {x}i, so
> > 
> > {dx}k  =  SUM(Fik * {x}i) =>
> > 
> > ([A] - Lk*[I]) * {dx}k = ([A] - Lk*[I]) * SUM(Fik * {x}i)
> >       =
> > SUM(Fik * ([A] - Lk*[I])*{x}i ) =
> > SUM(Fik * ([A]*{x}i - Lk*[I]*{x}i) ) =
> > SUM(Fik * (Li*{x}i  - Lk*{x}i) ) =
> > SUM(Fik * (Li - Lk)*{x}i )
> > 
> > So 4)   becomes
> > 
> > SUM(Fik * (Li - Lk)*{x}i ) = ([I]*dLk - [dA]) * {x}k
> > 
> > multiplying at left side of equality with Transpose of {x}k results
> > 
> > {x}k^T * SUM(Fik * (Li - Lk)*{x}i  = {x}k^T * ([I]*dLk - [dA]) * {x}k   <=>
> > 
> > SUM(Fik * (Li -Lk) * {x}k^T * {x}i = {x}k^T * ([I]*dLk - [dA]) * {x}k
> > 
> >  From 2) results that for i != k  {x}k^T * {x}i = 0,
> > also for i=k,  Li - Lk = 0, so left side of equality becames 0:
> > 
> > 0 = {x}k^T * ([I]*dLk - [dA]) * {x}k
> > 
> > if we note unity eigenvector {v}k = {x}k/Norm({x}k):
> > 
> > 0 = {v}k^T *  ([I]*dLk - [dA]) * {v}k
> > 
> > 
> > Transforming bilinear form into orthonormal basis of eigenvectors,
> > equality becames:
> > 
> > 5)   0 = {e}k^T * ([C]^T * ([I]*dLk - [dA]) * [C]) * {e}k
> > 
> >     ,with {e}k versor of orthonormal eigenvalues basis.
> > 
> > Let's note:
> > 
> > [D] = [C]^T * ([I]*dLk - [dA]) * [C]  <=>
> > 
> > [D] = [C]^T*[I]*[C] * dLk  -  [C]^T*[dA]*[C]  <=>
> > 
> > [D] = [I]*dLk - [C]^T*[dA]*[C]
> > 
> > 
> > Let's note  [B] = [C]^T*[dA]*[C]  ,so
> > 
> > 
> > [D] = [I]*dLk - [B]  ,so 5)   becames:
> > 
> > 
> > 0 = {e}k^T * ([I]*dLk - [B]) * {e}k  <=>
> > 
> > 0 = {e}k^T * [I] * {e}k * dLk  -  {e}k^T * [B] * {e}k  <=>
> > 
> > 0= dLk - Bkk <=>
> > 
> > *************
> > * dLk = Bkk *
> > *************
> > /Proof>
> > 
> > Quickly  tested for matrix
> > 
> > 0 1
> > 1 0
> > 
> > with small perturbation [dA]
> > 
> > a11 a12
> > a21 a22
> > 
> > ,aij << 1
> > 
> > and it holds :)
> > 
> > 
> > Have a nice day,
> > hum9
> 
> I reformulated this proposed theorem with less constraints:
> 
> If real NxN matrix [A] has N orthonormal eigenvectors,
> for a small perturbation [dA] corresponds small eigenvalues changes:
> 	dLi = Bii , where [B] = [C]^T * [dA] * [C],
> with matrix [C] of orthonormal eigenvectors as columns.
> 

That's a clearer statement, but a definition of "small" is needed. Let's try a numerical experiment, using R, below. The results seem to show that your theorem doesn't work very well, when the perturbations are at the 1e-3 level.

Krishna
---------------------

R version 3.3.3 (2017-03-06) -- "Another Canoe"
Copyright (C) 2017 The R Foundation for Statistical Computing
Platform: x86_64-redhat-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

## Make a square matrix, [A]
> A <- matrix(c(12, -51, 4, 6, 167, -68, -4, 24, -41),
  nrow = 3, byrow = TRUE)
> A
     [,1] [,2] [,3]
[1,]   12  -51    4
[2,]    6  167  -68
[3,]   -4   24  -41


## Obtain the eigenvectors of A as the matrix C.
> C <- eigen(A)$vectors
> C
           [,1]      [,2]       [,3]
[1,]  0.3281474 0.2547577 -0.9905263
[2,] -0.9368814 0.3027934  0.0871754
[3,] -0.1207170 0.9183761  0.1061044
> 


## Set of uniform random numbers between 0 and 1
> p <- runif(length(A))
> p
[1] 0.30973242 0.73219453 0.28703609 0.98814978 0.71599764 0.86366133 0.94912773
[8] 0.50863337 0.03185287

## Make a small perturbation matrix, dA, from
## random numbers.
> dA <- matrix(p, nrow=3, byrow=TRUE)
> dA
          [,1]      [,2]       [,3]
[1,] 0.3097324 0.7321945 0.28703609
[2,] 0.9881498 0.7159976 0.86366133
[3,] 0.9491277 0.5086334 0.03185287

## Scale dA so perturbations are small (relative
## to elements of A.
> dA <- dA / 1000
> dA
             [,1]         [,2]         [,3]
[1,] 0.0003097324 0.0007321945 2.870361e-04
[2,] 0.0009881498 0.0007159976 8.636613e-04
[3,] 0.0009491277 0.0005086334 3.185287e-05

# The transposed matrix, C^T
> CT <- t(C)
> CT
           [,1]       [,2]       [,3]
[1,]  0.3281474 -0.9368814 -0.1207170
[2,]  0.2547577  0.3027934  0.9183761
[3,] -0.9905263  0.0871754  0.1061044

## The perturbed matrix, A' = A + dA

> A_prime <- A + dA
> A_prime
          [,1]      [,2]       [,3]
[1,] 12.000310 -50.99927   4.000287
[2,]  6.000988 167.00072 -67.999136
[3,] -3.999051  24.00051 -40.999968

## Compute eigenvalues for A and A_prime

> L <- eigen(A)$values
> L_prime <- eigen(A_prime)$values

## Compute change in eigenvalues using hum9's theorem:
> B <- CT %*% dA %*% C
> B
              [,1]          [,2]          [,3]
[1,]  0.0002396215 -0.0010482326  8.106709e-04
[2,] -0.0004493261  0.0009161406 -1.123436e-03
[3,]  0.0005559052 -0.0004037818  4.391405e-05

## Compare diagonal elements of above matrix
## with actual change in eigenvalues
> L_prime - L
[1] 0.0004830448 0.0003955184 0.0001790197
> 

[toc] | [prev] | [next] | [standalone]

#60780

On 15.06.2017 13:30, krishna.myneni@ccreweb.org wrote:
> On Thursday, June 15, 2017 at 1:59:44 AM UTC-5, hum9 wrote:
>> On 14.06.2017 18:05, hum9 wrote:
>>> On 12.06.2017 17:14, krishna...org wrote:
>>>> Hi Marcel,
>>>>
>>>> I need a fast determinant to write some pedagogical code
>>>> for illustrating a molecular physics calculation. It is
>>>> to be used for solving the eigenvalue problem,
>>>>
>>>>       D( H - lambda×S ) = 0
>>>>
>>>> where H and S are N×N matrices, and lambda(s) are the
>>>> eigenvalue(s) to be found. D() is, of course, the
>>>> determinant operator. This form of the eigenvalue
>>>> problem is typical in molecular physics calculations,
>>>> when the basis set of eigenfunctions are non-orthogonal.
>>>>
>>>> The matrix S, consists of the inner products, of the N
>>>> eigenfunctions, and it would simply be the identity
>>>> matrix if the N eigenfunctions were orthonormal.
>>>> Typically, the hard part in setting up the problem is
>>>> evaluating the elements of H and S, each of the elements
>>>> requiring the evaluation of a multi-dimensional integral.
>>>> Nevertheless, when N becomes large, efficient evaluation
>>>> of the determinants is needed, particularly in searching
>>>> for the eigenvalues.
>>>>
>>>> A simple way of solving the equation, once the matrix
>>>> elements of H and S have been computed in the basis set,
>>>> is to simply use the determinant equation within a root
>>>> finder, to search for the values of lambda which satisfy
>>>> the equation. This is conceptually simple and straight-
>>>> forward to implement, but it's probably not the most
>>>> computationally efficient method -- there are probably
>>>> much better techniques available in a linear algebra
>>>> library. But, for pedagogical use, it's nice. Typically,
>>>> we are interested in only the lowest eigenvalue.
>>>>
>>>> Krishna
>>>>
>>>
>>> Recently, after your post, I have found something that
>>> could be useful:
>>>
>>>
>>> A symbol between '[' ']' represents a matrix,
>>> between '{' '}' represents vector,
>>> ^T   means transpose
>>> ^-1  means inverse of matrix.
>>>
>>> Theorem (proposition):
>>>
>>> Be a real symmetric matrix [A] NxN  with N distinct nonzero eigenvalues
>>> and [C] its transformation matrix into orthonormal basis of [A]
>>> eigenvectors.
>>>
>>> Then, for very small perturbation [dA] (Norm([dA] << Norm([A]))
>>> will result delta eigenvalues dL, for every index k, 1 <= k <= N
>>>
>>> dLk = Bkk  ,where [B] = [C]^T * [dA] * [C]
>>>
>>>
>>>
>>> <Proof/
>>>
>>> [A] is a symmetric matrix NxN, with N distinct eigenvalues.
>>> That means:
>>>
>>>        1) eigenvectors {x}i, 1 <= i <= N are linear independents
>>>
>>>        2) for every i j,   1 <= i != j <=N,   {x}i^T * {x}j = 0.
>>>
>>> ---
>>> 3)   [A]*{x}=L*(x).
>>>
>>> For a very small perturbation [dA] will follows:
>>>
>>> [A + dA] * {x + dx} = (L + dl) * {x + dx} <=>
>>>
>>> [A]*{x} + [A]*{dx} + [dA]*{x} + [dA]*{dx}
>>>       =
>>> L*{x} + L*{dx} + dL*{x} + dL*{dx}
>>>
>>> using 3)   =>
>>>
>>> [A]*{dx} + [dA]*{x} + [dA]*{dx}
>>>        =
>>> L*{dx} + dL*{x} + dL*{dx}
>>>
>>> For Norm([dA]) << Norm([A]) , second order terms are negligible, so
>>>
>>> [A]*{dx} + [dA]*{x} = L*{dx} + dL*{x} <=>
>>>
>>> ([A] - L*[I])*{dx} = ([I]*dL - [dA])*{x} <=>
>>>
>>> For every 1 <= k <= N:
>>>
>>> 4)   ([A] -Lk*[I]) * {dx}k = ([I]*dLk - [dA]) * {x}k
>>>
>>> Using 1), vector {dx}k can be expressed in terms of eigenvectors {x}i, so
>>>
>>> {dx}k  =  SUM(Fik * {x}i) =>
>>>
>>> ([A] - Lk*[I]) * {dx}k = ([A] - Lk*[I]) * SUM(Fik * {x}i)
>>>        =
>>> SUM(Fik * ([A] - Lk*[I])*{x}i ) =
>>> SUM(Fik * ([A]*{x}i - Lk*[I]*{x}i) ) =
>>> SUM(Fik * (Li*{x}i  - Lk*{x}i) ) =
>>> SUM(Fik * (Li - Lk)*{x}i )
>>>
>>> So 4)   becomes
>>>
>>> SUM(Fik * (Li - Lk)*{x}i ) = ([I]*dLk - [dA]) * {x}k
>>>
>>> multiplying at left side of equality with Transpose of {x}k results
>>>
>>> {x}k^T * SUM(Fik * (Li - Lk)*{x}i  = {x}k^T * ([I]*dLk - [dA]) * {x}k   <=>
>>>
>>> SUM(Fik * (Li -Lk) * {x}k^T * {x}i = {x}k^T * ([I]*dLk - [dA]) * {x}k
>>>
>>>   From 2) results that for i != k  {x}k^T * {x}i = 0,
>>> also for i=k,  Li - Lk = 0, so left side of equality becames 0:
>>>
>>> 0 = {x}k^T * ([I]*dLk - [dA]) * {x}k
>>>
>>> if we note unity eigenvector {v}k = {x}k/Norm({x}k):
>>>
>>> 0 = {v}k^T *  ([I]*dLk - [dA]) * {v}k
>>>
>>>
>>> Transforming bilinear form into orthonormal basis of eigenvectors,
>>> equality becames:
>>>
>>> 5)   0 = {e}k^T * ([C]^T * ([I]*dLk - [dA]) * [C]) * {e}k
>>>
>>>      ,with {e}k versor of orthonormal eigenvalues basis.
>>>
>>> Let's note:
>>>
>>> [D] = [C]^T * ([I]*dLk - [dA]) * [C]  <=>
>>>
>>> [D] = [C]^T*[I]*[C] * dLk  -  [C]^T*[dA]*[C]  <=>
>>>
>>> [D] = [I]*dLk - [C]^T*[dA]*[C]
>>>
>>>
>>> Let's note  [B] = [C]^T*[dA]*[C]  ,so
>>>
>>>
>>> [D] = [I]*dLk - [B]  ,so 5)   becames:
>>>
>>>
>>> 0 = {e}k^T * ([I]*dLk - [B]) * {e}k  <=>
>>>
>>> 0 = {e}k^T * [I] * {e}k * dLk  -  {e}k^T * [B] * {e}k  <=>
>>>
>>> 0= dLk - Bkk <=>
>>>
>>> *************
>>> * dLk = Bkk *
>>> *************
>>> /Proof>
>>>
>>> Quickly  tested for matrix
>>>
>>> 0 1
>>> 1 0
>>>
>>> with small perturbation [dA]
>>>
>>> a11 a12
>>> a21 a22
>>>
>>> ,aij << 1
>>>
>>> and it holds :)
>>>
>>>
>>> Have a nice day,
>>> hum9
>>
>> I reformulated this proposed theorem with less constraints:
>>
>> If real NxN matrix [A] has N orthonormal eigenvectors,
>> for a small perturbation [dA] corresponds small eigenvalues changes:
>> 	dLi = Bii , where [B] = [C]^T * [dA] * [C],
>> with matrix [C] of orthonormal eigenvectors as columns.
>>
> 
> That's a clearer statement, but a definition of "small" is needed. Let's try a numerical experiment, using R, below. The results seem to show that your theorem doesn't work very well, when the perturbations are at the 1e-3 level.
> 
> Krishna
> ---------------------
> 
> R version 3.3.3 (2017-03-06) -- "Another Canoe"
> Copyright (C) 2017 The R Foundation for Statistical Computing
> Platform: x86_64-redhat-linux-gnu (64-bit)
> 
> R is free software and comes with ABSOLUTELY NO WARRANTY.
> You are welcome to redistribute it under certain conditions.
> Type 'license()' or 'licence()' for distribution details.
> 
>    Natural language support but running in an English locale
> 
> R is a collaborative project with many contributors.
> Type 'contributors()' for more information and
> 'citation()' on how to cite R or R packages in publications.
> 
> Type 'demo()' for some demos, 'help()' for on-line help, or
> 'help.start()' for an HTML browser interface to help.
> Type 'q()' to quit R.
> 
> ## Make a square matrix, [A]
>> A <- matrix(c(12, -51, 4, 6, 167, -68, -4, 24, -41),
>    nrow = 3, byrow = TRUE)
>> A
>       [,1] [,2] [,3]
> [1,]   12  -51    4
> [2,]    6  167  -68
> [3,]   -4   24  -41
> 
> 
> ## Obtain the eigenvectors of A as the matrix C.
>> C <- eigen(A)$vectors
>> C
>             [,1]      [,2]       [,3]
> [1,]  0.3281474 0.2547577 -0.9905263
> [2,] -0.9368814 0.3027934  0.0871754
> [3,] -0.1207170 0.9183761  0.1061044
>>
> 
> 
> ## Set of uniform random numbers between 0 and 1
>> p <- runif(length(A))
>> p
> [1] 0.30973242 0.73219453 0.28703609 0.98814978 0.71599764 0.86366133 0.94912773
> [8] 0.50863337 0.03185287
> 
> ## Make a small perturbation matrix, dA, from
> ## random numbers.
>> dA <- matrix(p, nrow=3, byrow=TRUE)
>> dA
>            [,1]      [,2]       [,3]
> [1,] 0.3097324 0.7321945 0.28703609
> [2,] 0.9881498 0.7159976 0.86366133
> [3,] 0.9491277 0.5086334 0.03185287
> 
> ## Scale dA so perturbations are small (relative
> ## to elements of A.
>> dA <- dA / 1000
>> dA
>               [,1]         [,2]         [,3]
> [1,] 0.0003097324 0.0007321945 2.870361e-04
> [2,] 0.0009881498 0.0007159976 8.636613e-04
> [3,] 0.0009491277 0.0005086334 3.185287e-05
> 
> # The transposed matrix, C^T
>> CT <- t(C)
>> CT
>             [,1]       [,2]       [,3]
> [1,]  0.3281474 -0.9368814 -0.1207170
> [2,]  0.2547577  0.3027934  0.9183761
> [3,] -0.9905263  0.0871754  0.1061044
> 
> ## The perturbed matrix, A' = A + dA
> 
>> A_prime <- A + dA
>> A_prime
>            [,1]      [,2]       [,3]
> [1,] 12.000310 -50.99927   4.000287
> [2,]  6.000988 167.00072 -67.999136
> [3,] -3.999051  24.00051 -40.999968
> 
> ## Compute eigenvalues for A and A_prime
> 
>> L <- eigen(A)$values
>> L_prime <- eigen(A_prime)$values
> 
> ## Compute change in eigenvalues using hum9's theorem:
>> B <- CT %*% dA %*% C
>> B
>                [,1]          [,2]          [,3]
> [1,]  0.0002396215 -0.0010482326  8.106709e-04
> [2,] -0.0004493261  0.0009161406 -1.123436e-03
> [3,]  0.0005559052 -0.0004037818  4.391405e-05
> 
> ## Compare diagonal elements of above matrix
> ## with actual change in eigenvalues
>> L_prime - L
> [1] 0.0004830448 0.0003955184 0.0001790197
>>

[C]^T * [C] should be [I], but is:

1.00000006788172 -0.31094703528148 -0.41952024562898
-0.31094703528148 0.99999998984406 -0.12850432115031
-0.41952024562898 -0.12850432115031 1.00000004505621

Looks that eigenvectors are not orthonormal.

Have a nice day,
hum9

[toc] | [prev] | [next] | [standalone]

#60799

On Thursday, June 15, 2017 at 9:57:56 AM UTC-5, hum9 wrote:
> On 15.06.2017 13:30, krishna...org wrote:
...
> > 
> > That's a clearer statement, but a definition of "small" is needed. Let's try a numerical experiment, using R, below. The results seem to show that your theorem doesn't work very well, when the perturbations are at the 1e-3 level.
> > 
...
> > ## Make a square matrix, [A]
> >> A <- matrix(c(12, -51, 4, 6, 167, -68, -4, 24, -41),
> >    nrow = 3, byrow = TRUE)
> >> A
> >       [,1] [,2] [,3]
> > [1,]   12  -51    4
> > [2,]    6  167  -68
> > [3,]   -4   24  -41
> > 
> > 
> > ## Obtain the eigenvectors of A as the matrix C.
> >> C <- eigen(A)$vectors
> >> C
> >             [,1]      [,2]       [,3]
> > [1,]  0.3281474 0.2547577 -0.9905263
> > [2,] -0.9368814 0.3027934  0.0871754
> > [3,] -0.1207170 0.9183761  0.1061044
> >>
...
> > 
> > ## The perturbed matrix, A' = A + dA
> > 
> >> A_prime <- A + dA
> >> A_prime
> >            [,1]      [,2]       [,3]
> > [1,] 12.000310 -50.99927   4.000287
> > [2,]  6.000988 167.00072 -67.999136
> > [3,] -3.999051  24.00051 -40.999968
> > 
> > ## Compute eigenvalues for A and A_prime
> > 
> >> L <- eigen(A)$values
> >> L_prime <- eigen(A_prime)$values
> > 
> > ## Compute change in eigenvalues using hum9's theorem:
> >> B <- CT %*% dA %*% C
> >> B
> >                [,1]          [,2]          [,3]
> > [1,]  0.0002396215 -0.0010482326  8.106709e-04
> > [2,] -0.0004493261  0.0009161406 -1.123436e-03
> > [3,]  0.0005559052 -0.0004037818  4.391405e-05
> > 
> > ## Compare diagonal elements of above matrix
> > ## with actual change in eigenvalues
> >> L_prime - L
> > [1] 0.0004830448 0.0003955184 0.0001790197
> >>
> 
> [C]^T * [C] should be [I], but is:
> 
> 1.00000006788172 -0.31094703528148 -0.41952024562898
> -0.31094703528148 0.99999998984406 -0.12850432115031
> -0.41952024562898 -0.12850432115031 1.00000004505621
> 
> Looks that eigenvectors are not orthonormal.
> ...

Yes, you're quite correct. In order to ensure the eigenvectors are orthogonal, we must start with a real symmetric matrix. Let's check the theorem in the case of a real symmetric matrix:

> A <- matrix(c(12, -51, 4, -51, 167, -68, 4, -68, -41),
+ nrow=3, byrow=TRUE)
> A
     [,1] [,2] [,3]
[1,]   12  -51    4
[2,]  -51  167  -68
[3,]    4  -68  -41

## The eigenvectors of A; each column
## is a separate eigenvector.
> C = eigen(A)$vectors
> C
           [,1]       [,2]      [,3]
[1,]  0.2569501  0.9526354 0.1626730
[2,] -0.9291492  0.1972167 0.3127096
[3,]  0.2658164 -0.2314982 0.9358153

## Transpose of eigenvector matrix:
> CT = t(C)
> CT
          [,1]       [,2]       [,3]
[1,] 0.2569501 -0.9291492  0.2658164
[2,] 0.9526354  0.1972167 -0.2314982
[3,] 0.1626730  0.3127096  0.9358153

## Check orthogonality of eigenvectors
## of A. Zeroing elements which have
## abs value < 1e-15 will leave us with
## the identity matrix, below. Yes, we
## have orthogonal eigenvectors!
> CT %*% C
             [,1]         [,2]         [,3]
[1,] 1.000000e+00 1.595946e-16 1.665335e-16
[2,] 1.595946e-16 1.000000e+00 1.110223e-16
[3,] 1.665335e-16 1.110223e-16 1.000000e+00
> 

## Create a perturbation matrix, dA:
> p <- runif(length(A))/1000
> dA <- matrix(p, nrow=3, byrow=TRUE)
> dA
             [,1]         [,2]         [,3]
[1,] 0.0002611730 0.0001076001 0.0004993786
[2,] 0.0002825590 0.0004750525 0.0003822412
[3,] 0.0008472763 0.0005271824 0.0007238125

## The perturbed matrix, A' = A + dA:
> A_prime <- A + dA
> A_prime
           [,1]      [,2]       [,3]
[1,]  12.000261 -50.99989   4.000499
[2,] -50.999717 167.00048 -67.999618
[3,]   4.000847 -67.99947 -40.999276

## Eigenvalues of the original matrix, A:
> L <- eigen(A)$values
> L
[1] 200.5575446   0.4698347 -63.0273793

## Eigenvalues of perturbed matrix, A':
> L_prime <- eigen(A_prime)$values
> L_prime
[1] 200.5577974 -63.0262011   0.4698638

## In comparing the eigenvalues of A and A',
## we need to be careful about ordering them.
> L <- sort(L)
> L
[1] -63.0273793   0.4698347 200.5575446
> L_prime <- sort(L_prime)
> L_prime
[1] -63.0262011   0.4698638 200.5577974
> dL_computed <- L_prime - L

## Actual change in eigenvalues between
## A and A' = A + dA:
> dL_computed
[1] 1.178222e-03 2.908929e-05 2.527263e-04

## First-order change in eigenvalues computed 
## with hum9's theorem; the changes are on
## the diagonal elements of the matrix below,
## which requires only the perturbations and
## the (orthogonal) eigenvectors of the original
## matrix:
> CT %*% dA %*% C
              [,1]          [,2]          [,3]
[1,]  2.527263e-04 -1.760979e-05 -0.0001129474
[2,]  6.090657e-05  2.908421e-05  0.0003997496
[3,] -1.419097e-04  8.067480e-04  0.0011782275
> 
---------------------------

Yes, the diagonal elements of C^T * dA * C match
the computed change in eigenvalues, for this real
symmetric matrix, A.


Krishna

[toc] | [prev] | [next] | [standalone]

#60802

On 16.06.2017 05:06, krishna.myneni@ccreweb.org wrote:
> On Thursday, June 15, 2017 at 9:57:56 AM UTC-5, hum9 wrote:
>> On 15.06.2017 13:30, krishna...org wrote:
> ...
>>>
>>> That's a clearer statement, but a definition of "small" is needed. Let's try a numerical experiment, using R, below. The results seem to show that your theorem doesn't work very well, when the perturbations are at the 1e-3 level.
>>>
> ...
>>> ## Make a square matrix, [A]
>>>> A <- matrix(c(12, -51, 4, 6, 167, -68, -4, 24, -41),
>>>     nrow = 3, byrow = TRUE)
>>>> A
>>>        [,1] [,2] [,3]
>>> [1,]   12  -51    4
>>> [2,]    6  167  -68
>>> [3,]   -4   24  -41
>>>
>>>
>>> ## Obtain the eigenvectors of A as the matrix C.
>>>> C <- eigen(A)$vectors
>>>> C
>>>              [,1]      [,2]       [,3]
>>> [1,]  0.3281474 0.2547577 -0.9905263
>>> [2,] -0.9368814 0.3027934  0.0871754
>>> [3,] -0.1207170 0.9183761  0.1061044
>>>>
> ...
>>>
>>> ## The perturbed matrix, A' = A + dA
>>>
>>>> A_prime <- A + dA
>>>> A_prime
>>>             [,1]      [,2]       [,3]
>>> [1,] 12.000310 -50.99927   4.000287
>>> [2,]  6.000988 167.00072 -67.999136
>>> [3,] -3.999051  24.00051 -40.999968
>>>
>>> ## Compute eigenvalues for A and A_prime
>>>
>>>> L <- eigen(A)$values
>>>> L_prime <- eigen(A_prime)$values
>>>
>>> ## Compute change in eigenvalues using hum9's theorem:
>>>> B <- CT %*% dA %*% C
>>>> B
>>>                 [,1]          [,2]          [,3]
>>> [1,]  0.0002396215 -0.0010482326  8.106709e-04
>>> [2,] -0.0004493261  0.0009161406 -1.123436e-03
>>> [3,]  0.0005559052 -0.0004037818  4.391405e-05
>>>
>>> ## Compare diagonal elements of above matrix
>>> ## with actual change in eigenvalues
>>>> L_prime - L
>>> [1] 0.0004830448 0.0003955184 0.0001790197
>>>>
>>
>> [C]^T * [C] should be [I], but is:
>>
>> 1.00000006788172 -0.31094703528148 -0.41952024562898
>> -0.31094703528148 0.99999998984406 -0.12850432115031
>> -0.41952024562898 -0.12850432115031 1.00000004505621
>>
>> Looks that eigenvectors are not orthonormal.
>> ...
> 
> Yes, you're quite correct. In order to ensure the eigenvectors are orthogonal, we must start with a real symmetric matrix. Let's check the theorem in the case of a real symmetric matrix:
> 
>> A <- matrix(c(12, -51, 4, -51, 167, -68, 4, -68, -41),
> + nrow=3, byrow=TRUE)
>> A
>       [,1] [,2] [,3]
> [1,]   12  -51    4
> [2,]  -51  167  -68
> [3,]    4  -68  -41
> 
> ## The eigenvectors of A; each column
> ## is a separate eigenvector.
>> C = eigen(A)$vectors
>> C
>             [,1]       [,2]      [,3]
> [1,]  0.2569501  0.9526354 0.1626730
> [2,] -0.9291492  0.1972167 0.3127096
> [3,]  0.2658164 -0.2314982 0.9358153
> 
> ## Transpose of eigenvector matrix:
>> CT = t(C)
>> CT
>            [,1]       [,2]       [,3]
> [1,] 0.2569501 -0.9291492  0.2658164
> [2,] 0.9526354  0.1972167 -0.2314982
> [3,] 0.1626730  0.3127096  0.9358153
> 
> ## Check orthogonality of eigenvectors
> ## of A. Zeroing elements which have
> ## abs value < 1e-15 will leave us with
> ## the identity matrix, below. Yes, we
> ## have orthogonal eigenvectors!
>> CT %*% C
>               [,1]         [,2]         [,3]
> [1,] 1.000000e+00 1.595946e-16 1.665335e-16
> [2,] 1.595946e-16 1.000000e+00 1.110223e-16
> [3,] 1.665335e-16 1.110223e-16 1.000000e+00
>>
> 
> ## Create a perturbation matrix, dA:
>> p <- runif(length(A))/1000
>> dA <- matrix(p, nrow=3, byrow=TRUE)
>> dA
>               [,1]         [,2]         [,3]
> [1,] 0.0002611730 0.0001076001 0.0004993786
> [2,] 0.0002825590 0.0004750525 0.0003822412
> [3,] 0.0008472763 0.0005271824 0.0007238125
> 
> ## The perturbed matrix, A' = A + dA:
>> A_prime <- A + dA
>> A_prime
>             [,1]      [,2]       [,3]
> [1,]  12.000261 -50.99989   4.000499
> [2,] -50.999717 167.00048 -67.999618
> [3,]   4.000847 -67.99947 -40.999276
> 
> ## Eigenvalues of the original matrix, A:
>> L <- eigen(A)$values
>> L
> [1] 200.5575446   0.4698347 -63.0273793
> 
> ## Eigenvalues of perturbed matrix, A':
>> L_prime <- eigen(A_prime)$values
>> L_prime
> [1] 200.5577974 -63.0262011   0.4698638
> 
> ## In comparing the eigenvalues of A and A',
> ## we need to be careful about ordering them.
>> L <- sort(L)
>> L
> [1] -63.0273793   0.4698347 200.5575446
>> L_prime <- sort(L_prime)
>> L_prime
> [1] -63.0262011   0.4698638 200.5577974
>> dL_computed <- L_prime - L
> 
> ## Actual change in eigenvalues between
> ## A and A' = A + dA:
>> dL_computed
> [1] 1.178222e-03 2.908929e-05 2.527263e-04
> 
> ## First-order change in eigenvalues computed
> ## with hum9's theorem; the changes are on
> ## the diagonal elements of the matrix below,
> ## which requires only the perturbations and
> ## the (orthogonal) eigenvectors of the original
> ## matrix:
>> CT %*% dA %*% C
>                [,1]          [,2]          [,3]
> [1,]  2.527263e-04 -1.760979e-05 -0.0001129474
> [2,]  6.090657e-05  2.908421e-05  0.0003997496
> [3,] -1.419097e-04  8.067480e-04  0.0011782275
>>
> ---------------------------
> 
> Yes, the diagonal elements of C^T * dA * C match
> the computed change in eigenvalues, for this real
> symmetric matrix, A.
> 
> 
> Krishna
> 

Thank you for support and for inspiration! :)

I reformulated that theorem so it looks now:

If real NxN matrix [A] has N orthonormal eigenvectors,
with matrix [C] of orthonormal eigenvectors as columns.

In the case of small perturbation of [A], for every
unperturbed eigenvector {x} corresponds an perturbed
eigenvector {x+dx} collinear with {x}.

Changing basis by [C] to a small perturbed version of [A]
will result an slight perturbed diagonal matrix with small
perturbed eigenvalues on its diagonal.

Have a nice day,
hum9

[toc] | [prev] | [next] | [standalone]

#60810

On Friday, June 16, 2017 at 2:21:21 AM UTC-5, hum9 wrote:
...
> Thank you for support and for inspiration! :)
> 

Glad to be of help.

> I reformulated that theorem so it looks now:
> 
> If real NxN matrix [A] has N orthonormal eigenvectors,
> with matrix [C] of orthonormal eigenvectors as columns.
> 
> In the case of small perturbation of [A], for every
> unperturbed eigenvector {x} corresponds an perturbed
> eigenvector {x+dx} collinear with {x}.
> 
> Changing basis by [C] to a small perturbed version of [A]
> will result an slight perturbed diagonal matrix with small
> perturbed eigenvalues on its diagonal.
> 

I think you're not stating the main result, which is the expression for the first-order correction to the eigenvalues. I would state the result in the following way:

Thm: Given a real matrix A, for which the eigenvectors
are orthogonal, a matrix C containing the eigenvectors
as columns, and a perturbation matrix dA, the 
first-order corrections to the eigenvalues of the
perturbed matrix, A' = A + dA, are given by the
diagonal elements of the similarity transformation,

   B = C^T dA C

(end of theorem)

This restatement of the result is missing the geometric
interpretation you suggest above, though, and that's
essential to understanding the result, but not to stating
the result. It's not clear to me that the transformation,

  e' = C e

where e is an eigenvector of A, gives a vector e',
which is collinear with e.

Also, I think stating the result in this manner makes it
easy to generalize to the case of complex matrices. In
any case, you should check the literature to see if this
result is already well-known to mathematicians.

Krishna

[toc] | [prev] | [next] | [standalone]

#60818

On 16.06.2017 13:38, krishna.myneni@ccreweb.org wrote:
> On Friday, June 16, 2017 at 2:21:21 AM UTC-5, hum9 wrote:
> ...
>> Thank you for support and for inspiration! :)
>>
> 
> Glad to be of help.
> 
>> I reformulated that theorem so it looks now:
>>
>> If real NxN matrix [A] has N orthonormal eigenvectors,
>> with matrix [C] of orthonormal eigenvectors as columns.
>>
>> In the case of small perturbation of [A], for every
>> unperturbed eigenvector {x} corresponds an perturbed
>> eigenvector {x+dx} collinear with {x}.
>>
>> Changing basis by [C] to a small perturbed version of [A]
>> will result an slight perturbed diagonal matrix with small
>> perturbed eigenvalues on its diagonal.
>>
> 
> I think you're not stating the main result, which is the expression for the first-order correction to the eigenvalues. I would state the result in the following way:
> 
> Thm: Given a real matrix A, for which the eigenvectors
> are orthogonal, a matrix C containing the eigenvectors
> as columns, and a perturbation matrix dA, the
> first-order corrections to the eigenvalues of the
> perturbed matrix, A' = A + dA, are given by the
> diagonal elements of the similarity transformation,
> 
>     B = C^T dA C
> 
> (end of theorem)
> 
> This restatement of the result is missing the geometric
> interpretation you suggest above, though, and that's
> essential to understanding the result, but not to stating
> the result. It's not clear to me that the transformation,
> 
>    e' = C e
> 
> where e is an eigenvector of A, gives a vector e',
> which is collinear with e.
> 
> Also, I think stating the result in this manner makes it
> easy to generalize to the case of complex matrices. In
> any case, you should check the literature to see if this
> result is already well-known to mathematicians.
> 
> Krishna
> 

Inside proof, it was demonstrated that left side of equality 4) is zero:

for every k=1..N

     ([A] -Lk*[I]) * {dx}k = 0

That means delta-eigenvector {dx}k is also an eigenvector of [A]!

Corresponding to eigenvalue Lk is eigenvector {x} also {dx}k,
so {x}k is parallel to {dx}k.

That means for an eigenvalue Lk of [A], perturbed eigenvector
{x+dx}k is collinear with unperturbed eigenvector of {x}k :)

hum9

[toc] | [prev] | [next] | [standalone]

#60821

On Friday, June 16, 2017 at 7:58:17 AM UTC-5, hum9 wrote:
...
> Inside proof, it was demonstrated that left side of equality 4) is zero:
> 
> for every k=1..N
> 
>      ([A] -Lk*[I]) * {dx}k = 0
> 
> That means delta-eigenvector {dx}k is also an eigenvector of [A]!
> 

The eigenvectors of A' = A + dA cannot be exact eigenvectors of A, in general. You are assuming this as an approximation to obtain the first-order corrections to the eigenvalues of A. Please look up perturbation theory, and, in particular, first-order perturbation theory.

Krishna

[toc] | [prev] | [next] | [standalone]

#60822

On 16.06.2017 16:04, krishna.myneni@ccreweb.org wrote:
> On Friday, June 16, 2017 at 7:58:17 AM UTC-5, hum9 wrote:
> ...
>> Inside proof, it was demonstrated that left side of equality 4) is zero:
>>
>> for every k=1..N
>>
>>       ([A] -Lk*[I]) * {dx}k = 0
>>
>> That means delta-eigenvector {dx}k is also an eigenvector of [A]!
>>
> 
> The eigenvectors of A' = A + dA cannot be exact eigenvectors of A, in general. You are assuming this as an approximation to obtain the first-order corrections to the eigenvalues of A. Please look up perturbation theory, and, in particular, first-order perturbation theory.
> 
> Krishna
> 

For small changes [dA], in first-order approximation, I have found that
eigenvectors of [A+dA] differs by eigenvector of [A] only by magnitude,
not by direction.

hum9

[toc] | [prev] | [next] | [standalone]

#60811

On Friday, June 16, 2017 at 9:21:21 AM UTC+2, hum9 wrote:
> On 16.06.2017 05:06, krishna.myneni@ccreweb.org wrote:
> > On Thursday, June 15, 2017 at 9:57:56 AM UTC-5, hum9 wrote:
> >> On 15.06.2017 13:30, krishna...org wrote:
...
> I reformulated that theorem so it looks now:
> 
> If real NxN matrix [A] has N orthonormal eigenvectors,
> with matrix [C] of orthonormal eigenvectors as columns.
> 
> In the case of small perturbation of [A], for every
> unperturbed eigenvector {x} corresponds an perturbed
> eigenvector {x+dx} collinear with {x}.
> 
> Changing basis by [C] to a small perturbed version of [A]
> will result an slight perturbed diagonal matrix with small
> perturbed eigenvalues on its diagonal.

I don't understand your formulation, sorry.
A = X D X^-1, where X is matrix of eigenvectors and D the 
diagonal matrix of A's eigenvalues.
Isn't it obvious that a (very) small perturbation of A
will result in (very) small perturbations of D and X? 
I don't even need orthonormality for that.

Interpreted physically, in my domain, it means that if
the system matrix of a 'reasonable' application (like 
an electronic low-pass filter) changes slightly, the 
eigenvectors and eigenvalues won't change dramatically.
I would even expect collinearity of the eigenvectors. 

The theorem would be more useful if you could extend it 
and say something about the relationship between small 
changes in A and the resulting changes in D and X.

-marcel

[toc] | [prev] | [next] | [standalone]

Page 1 of 2  [1] 2  Next page →

Back to top | Article view | comp.lang.forth

csiph-web