Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]
Groups > comp.lang.forth > #60780
| From | hum9 <obi65537@gmail.com> |
|---|---|
| Newsgroups | comp.lang.forth |
| Subject | Re: fast matrix determinant |
| Date | 2017-06-15 17:57 +0300 |
| Organization | A noiseless patient Spider |
| Message-ID | <ohu72p$nob$1@dont-email.me> (permalink) |
| References | (2 earlier) <adf343c7-467e-4e6c-b4e1-5d72f3d9ece3@googlegroups.com> <b6974c55-fda0-4102-8630-71522e32e842@googlegroups.com> <ohrj4a$jra$1@dont-email.me> <ohtb25$t0m$1@dont-email.me> <ef709c75-138b-4e9c-ad5b-08d40b24914f@googlegroups.com> |
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
Back to comp.lang.forth | Previous | Next — Previous in thread | Next in thread | Find similar | Unroll thread
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
csiph-web