FAQFAQ   SearchSearch   MemberlistMemberlist   UsergroupsUsergroups 
 ProfileProfile   PreferencesPreferences   Log in to check your private messagesLog in to check your private messages   Log inLog in 
Forum index » Science and Technology » Math » num-analysis
Matrix functions via EVD decomposition
Post new topic   Reply to topic Page 1 of 1 [8 Posts] View previous topic :: View next topic
Author Message
dan1112
science forum beginner


Joined: 15 Jul 2006
Posts: 3

PostPosted: Tue Jul 18, 2006 6:10 am    Post subject: Re: Matrix functions via EVD decomposition Reply with quote

NEVER use a Jordan Decomposition - it can't be computed stabley. You
can get similar information from the Schur factoristion (which is an
orthogonal decomposition and can therefore be computed stabley). The
Jordan Decomposition is only really useful for theory. I'm sure you'll
find that if you look a bit harder at Golub and van Loan
Jeremy Watts wrote:
Quote:
"~Glynne" <glynnec2002@yahoo.com> wrote in message
news:1152942674.413923.309690@35g2000cwc.googlegroups.com...
I was looking for a way to numerically evaluate functions of a real
matrix. The EigenvalueDecomposition class in the JAva MAtrix (JAMA)
package contained almost everything that was needed.

The decomposition produces matrices V and D such that the columns of V
contain the eigenvectors of A, while D contains the eigenvalues, i.e.
A.V = V.D

The eigenvalue matrix D is block diagonal with the real eigenvalues in
1-by-1 blocks and any complex eigenvalues, x+iy, in 2-by-2 blocks, [x,
y; -y, x].

Therefore, I should be able to compute f(A) using
f(A) = V.f(D)/V


Two computational issues arose:
1) how to calculate f(D)
2) how to deal with ill-conditioning of V


The first issue can be dealt with by representing the 2-by-2 diagonal
blocks as complex numbers and requiring that f(z) be defined for
complex arguments:
f(x + iy) = g(x,y) + i h(x,y)

This result means that a 2-by-2 block, [x, y; -y, x], is transformed by
f() into a new 2-by-2 block, [g, h; -h, g]


Handling an ill-conditioned V is trickier. I started by considering a
perturbation matrix E such that:
1) the EVD of (A+E) is not ill-conditioned
2) ||E|| << ||A||
3) E.A = A.E

yes you're likely to get an ill conditioned V if your eigenvalues are close
together, or a singular V of course
if you have repeat eigenvalues, in which case you cannot form the eigen
decomposition.
have you thought about using a jordan decomposition instead of an eigen?

for the jordan decomposition then A = SJS^-1 where J is a Jordan block
matrix, and S a canonical basis
a function then would be f(A) = Sf(J)S^-1 . it is a little bit more
difficult to find f(J) than f(D) of an eigen
decomposition but not too much (see golub von loan for details on finding
matrix functions using the Jordan
decomposition)

Utilizing such a perturbation results in a nice, simple Taylor
expansion
f(A+E) = f(A)/0! + E.{D}f(A)/1! + E^2.{D^2}f(A)/2! + O(E^3)


[[Beyond this point, I'm not completely satisfied with my reasoning,
although my Java application works just fine.]]


Finding a matrix which commutes with A was too difficult, so I
considered what would happen to the Taylor expansion if a small, random
matrix were used in two independent evaluations:
f(A+ E) = f(A)/0! + {D}f(E; A)/1! + {D^2}f(E; A)/2! + O(E^3)
f(A+2E) = f(A)/0! + {D}f(2E;A)/1! + {D^2}f(2E;A)/2! + O(E^3)


The non-commutivity of A with E makes the evaluation of {D^k}f(E;A)
terms very messy, but a scalar multiplier can be pulled clear of the
mess, i.e.

{D^k}f(s*E;A) = (s^k)*{D^k}f(E;A)


Letting e = E/||E||, we can say something like
f(A+ E) = f(A) + ||E||.{D}f(e;A) + ||E||^2.{D^2}f(e;A)/2 + ...
f(A+2E) = f(A) + ||2E||.{D}f(e;A) + ||2E||^2.{D^2}f(e;A)/2 + ...

without worrying over the commutation details of the individual terms.
The point is that the terms in each expansion are identical up to a
scale factor.

Subtraction yields the approximation
f(A) = 2f(A+E) - f(A+2E) + O(||E||^2)

If E is chosen such that
||E|| = sqrt(eps)*||A||

then the overall approximation should be accurate to O(eps).


Does anyone know of a nice expression for f(A+E) when A and E don't
commute?

Can this idea be extended to higher orders by evaluating f(A+E),
f(A+2E), f(A+3E), ... f(A+kE) and then combining them in such a way
that all terms cancel except for {D^k}f(E;A)/k! and f(A)?

~Glynne
Back to top
~Glynne
science forum beginner


Joined: 21 Nov 2005
Posts: 18

PostPosted: Mon Jul 17, 2006 8:21 pm    Post subject: Re: Matrix functions via EVD decomposition Reply with quote

~Glynne wrote:
Quote:
Carl Barron wrote:
~Glynne <glynnec2002@yahoo.com> wrote:

I was looking for a way to numerically evaluate functions of a real
matrix. The EigenvalueDecomposition class in the JAva MAtrix (JAMA)
package contained almost everything that was needed.

The decomposition produces matrices V and D such that the columns of V
contain the eigenvectors of A, while D contains the eigenvalues, i.e.
A.V = V.D

The eigenvalue matrix D is block diagonal with the real eigenvalues in
1-by-1 blocks and any complex eigenvalues, x+iy, in 2-by-2 blocks, [x,
y; -y, x].

Therefore, I should be able to compute f(A) using
f(A) = V.f(D)/V


Two computational issues arose:
1) how to calculate f(D)
2) how to deal with ill-conditioning of V


The first issue can be dealt with by representing the 2-by-2 diagonal
blocks as complex numbers and requiring that f(z) be defined for
complex arguments:
f(x + iy) = g(x,y) + i h(x,y)

This result means that a 2-by-2 block, [x, y; -y, x], is transformed by
f() into a new 2-by-2 block, [g, h; -h, g]


Handling an ill-conditioned V is trickier. I started by considering a
perturbation matrix E such that:
1) the EVD of (A+E) is not ill-conditioned
2) ||E|| << ||A||
3) E.A = A.E

Utilizing such a perturbation results in a nice, simple Taylor
expansion
f(A+E) = f(A)/0! + E.{D}f(A)/1! + E^2.{D^2}f(A)/2! + O(E^3)


[[Beyond this point, I'm not completely satisfied with my reasoning,
although my Java application works just fine.]]


Finding a matrix which commutes with A was too difficult, so I
considered what would happen to the Taylor expansion if a small, random
matrix were used in two independent evaluations:
f(A+ E) = f(A)/0! + {D}f(E; A)/1! + {D^2}f(E; A)/2! + O(E^3)
f(A+2E) = f(A)/0! + {D}f(2E;A)/1! + {D^2}f(2E;A)/2! + O(E^3)


The non-commutivity of A with E makes the evaluation of {D^k}f(E;A)
terms very messy, but a scalar multiplier can be pulled clear of the
mess, i.e.

{D^k}f(s*E;A) = (s^k)*{D^k}f(E;A)


Letting e = E/||E||, we can say something like
f(A+ E) = f(A) + ||E||.{D}f(e;A) + ||E||^2.{D^2}f(e;A)/2 + ...
f(A+2E) = f(A) + ||2E||.{D}f(e;A) + ||2E||^2.{D^2}f(e;A)/2 + ...

without worrying over the commutation details of the individual terms.
The point is that the terms in each expansion are identical up to a
scale factor.

Subtraction yields the approximation
f(A) = 2f(A+E) - f(A+2E) + O(||E||^2)

If E is chosen such that
||E|| = sqrt(eps)*||A||

then the overall approximation should be accurate to O(eps).


Does anyone know of a nice expression for f(A+E) when A and E don't
commute?

Can this idea be extended to higher orders by evaluating f(A+E),
f(A+2E), f(A+3E), ... f(A+kE) and then combining them in such a way
that all terms cancel except for {D^k}f(E;A)/k! and f(A)?

~Glynne
Good luck! consider A below it really has one eignevector and three
identical eigenvalues 2, given the V produce from TNT's JAMA below
V is essentially |0|*|1 1 1| , it should be singular there is
|0| only one eigenvector (0,0,1)^T
|1|

A =
2 0 0
1 2 0
0 1 2

V =
0 0 3.15544e-30
0 -1.77636e-15 -1.77636e-15
1 1 1

D =
2 0 0
0 2 0
0 0 2

AV =
0 0 6.31089e-30
0 -3.55271e-15 -3.55271e-15
2 2 2

VD =
0 0 6.31089e-30
0 -3.55271e-15 -3.55271e-15
2 2 2

Yes, I realize that your A matrix (which is a 3x3 Jordan block) cannot
be evaluated by the simple formula
f(A) = V.f(D)/V

But this is precisely the point of my post. Sticking with your A
matrix. It has an EVD decomposition for which
cond(V) = infinity.

However, let's pick a random perturbation E matrix:

0.42 0.71 0.31
0.01 0.42 0.40
0.54 0.43 0.12

And scale it by 10^(-8)

Now the EVD of (A+E) produces a new V, for which
cond(V) = 1.85 x 10^(-8)

while the EVD of (A+2E) produces
cond(V) = 0.33 x 10^(-8)

These condition numbers are large, but sufficient for IEEE double
precision calculations. For IEEE single precision scaling E by 10^(-4)
would be more appropriate.

Now I CAN calculate f(A+E) by the simple formula
f(A+E) = V.f(D)/V

...of course, the V & D matrices used are not those of the original A
matrix, but those of the perturbed matrix (A+E).

Similarly, I can calculate f(A+2E), f(A+3E), .... f(A+kE).

So what is the best way to combine these functional values to
accurately approximate f(A) -- which cannot be evaluated from the EVD
of A?

~Glynne

Pardon the typos.
For (A+E) and (A+2E), cond(V) equals
1.85 x 10^8
0.33 x 10^8
respectively.
Back to top
~Glynne
science forum beginner


Joined: 21 Nov 2005
Posts: 18

PostPosted: Sat Jul 15, 2006 10:14 pm    Post subject: Re: Matrix functions via EVD decomposition Reply with quote

dan wrote:
Quote:
Using an eigenvalue decomposition is fine for nearly normal matrices,
in which case f(D) = diag(f(lambda_i), i=1..n). If the eigenbasis is
poorly conditioned, or if the matrix is not diagonalisable, calculate a
schur decomposition. There is a recurrence relationship known as
Parlett's recurrence (which is O(n^3) and has some problems which you
can look up). For other methods look up Golub and van Loan's book
"MAtrix Computions" which has a chapter devoted to matrix functions.
If the matrix is sparse, that's a whole other problem. The best known
way to compute f(A)b is by using Krylov subspace methods, however, this
is still a major area of research in numerical linear algebra.

To sumerise - if the eigenbasis is ill conditioned, don't use an
eigenvalue decomposition (even if you need to write your own code). On
a somewhat unrealated topic, why are you using Java - it's note really
the best language for scientific computing.

The reason that I'm using Java is because, by utilizing two packages --
JAMA and BeanShell -- you get an interactive system which has roughly
the same functionality as MATLAB....for free.

~Glynne
Back to top
~Glynne
science forum beginner


Joined: 21 Nov 2005
Posts: 18

PostPosted: Sat Jul 15, 2006 10:08 pm    Post subject: Re: Matrix functions via EVD decomposition Reply with quote

Carl Barron wrote:
Quote:
~Glynne <glynnec2002@yahoo.com> wrote:

I was looking for a way to numerically evaluate functions of a real
matrix. The EigenvalueDecomposition class in the JAva MAtrix (JAMA)
package contained almost everything that was needed.

The decomposition produces matrices V and D such that the columns of V
contain the eigenvectors of A, while D contains the eigenvalues, i.e.
A.V = V.D

The eigenvalue matrix D is block diagonal with the real eigenvalues in
1-by-1 blocks and any complex eigenvalues, x+iy, in 2-by-2 blocks, [x,
y; -y, x].

Therefore, I should be able to compute f(A) using
f(A) = V.f(D)/V


Two computational issues arose:
1) how to calculate f(D)
2) how to deal with ill-conditioning of V


The first issue can be dealt with by representing the 2-by-2 diagonal
blocks as complex numbers and requiring that f(z) be defined for
complex arguments:
f(x + iy) = g(x,y) + i h(x,y)

This result means that a 2-by-2 block, [x, y; -y, x], is transformed by
f() into a new 2-by-2 block, [g, h; -h, g]


Handling an ill-conditioned V is trickier. I started by considering a
perturbation matrix E such that:
1) the EVD of (A+E) is not ill-conditioned
2) ||E|| << ||A||
3) E.A = A.E

Utilizing such a perturbation results in a nice, simple Taylor
expansion
f(A+E) = f(A)/0! + E.{D}f(A)/1! + E^2.{D^2}f(A)/2! + O(E^3)


[[Beyond this point, I'm not completely satisfied with my reasoning,
although my Java application works just fine.]]


Finding a matrix which commutes with A was too difficult, so I
considered what would happen to the Taylor expansion if a small, random
matrix were used in two independent evaluations:
f(A+ E) = f(A)/0! + {D}f(E; A)/1! + {D^2}f(E; A)/2! + O(E^3)
f(A+2E) = f(A)/0! + {D}f(2E;A)/1! + {D^2}f(2E;A)/2! + O(E^3)


The non-commutivity of A with E makes the evaluation of {D^k}f(E;A)
terms very messy, but a scalar multiplier can be pulled clear of the
mess, i.e.

{D^k}f(s*E;A) = (s^k)*{D^k}f(E;A)


Letting e = E/||E||, we can say something like
f(A+ E) = f(A) + ||E||.{D}f(e;A) + ||E||^2.{D^2}f(e;A)/2 + ...
f(A+2E) = f(A) + ||2E||.{D}f(e;A) + ||2E||^2.{D^2}f(e;A)/2 + ...

without worrying over the commutation details of the individual terms.
The point is that the terms in each expansion are identical up to a
scale factor.

Subtraction yields the approximation
f(A) = 2f(A+E) - f(A+2E) + O(||E||^2)

If E is chosen such that
||E|| = sqrt(eps)*||A||

then the overall approximation should be accurate to O(eps).


Does anyone know of a nice expression for f(A+E) when A and E don't
commute?

Can this idea be extended to higher orders by evaluating f(A+E),
f(A+2E), f(A+3E), ... f(A+kE) and then combining them in such a way
that all terms cancel except for {D^k}f(E;A)/k! and f(A)?

~Glynne
Good luck! consider A below it really has one eignevector and three
identical eigenvalues 2, given the V produce from TNT's JAMA below
V is essentially |0|*|1 1 1| , it should be singular there is
|0| only one eigenvector (0,0,1)^T
|1|

A =
2 0 0
1 2 0
0 1 2

V =
0 0 3.15544e-30
0 -1.77636e-15 -1.77636e-15
1 1 1

D =
2 0 0
0 2 0
0 0 2

AV =
0 0 6.31089e-30
0 -3.55271e-15 -3.55271e-15
2 2 2

VD =
0 0 6.31089e-30
0 -3.55271e-15 -3.55271e-15
2 2 2

Yes, I realize that your A matrix (which is a 3x3 Jordan block) cannot
be evaluated by the simple formula
f(A) = V.f(D)/V

But this is precisely the point of my post. Sticking with your A
matrix. It has an EVD decomposition for which
cond(V) = infinity.

However, let's pick a random perturbation E matrix:

0.42 0.71 0.31
0.01 0.42 0.40
0.54 0.43 0.12

And scale it by 10^(-8)

Now the EVD of (A+E) produces a new V, for which
cond(V) = 1.85 x 10^(-8)

while the EVD of (A+2E) produces
cond(V) = 0.33 x 10^(-8)

These condition numbers are large, but sufficient for IEEE double
precision calculations. For IEEE single precision scaling E by 10^(-4)
would be more appropriate.

Now I CAN calculate f(A+E) by the simple formula
f(A+E) = V.f(D)/V

....of course, the V & D matrices used are not those of the original A
matrix, but those of the perturbed matrix (A+E).

Similarly, I can calculate f(A+2E), f(A+3E), .... f(A+kE).

So what is the best way to combine these functional values to
accurately approximate f(A) -- which cannot be evaluated from the EVD
of A?

~Glynne
Back to top
Jeremy Watts
science forum Guru Wannabe


Joined: 24 Mar 2005
Posts: 239

PostPosted: Sat Jul 15, 2006 6:02 pm    Post subject: Re: Matrix functions via EVD decomposition Reply with quote

"~Glynne" <glynnec2002@yahoo.com> wrote in message
news:1152942674.413923.309690@35g2000cwc.googlegroups.com...
Quote:
I was looking for a way to numerically evaluate functions of a real
matrix. The EigenvalueDecomposition class in the JAva MAtrix (JAMA)
package contained almost everything that was needed.

The decomposition produces matrices V and D such that the columns of V
contain the eigenvectors of A, while D contains the eigenvalues, i.e.
A.V = V.D

The eigenvalue matrix D is block diagonal with the real eigenvalues in
1-by-1 blocks and any complex eigenvalues, x+iy, in 2-by-2 blocks, [x,
y; -y, x].

Therefore, I should be able to compute f(A) using
f(A) = V.f(D)/V


Two computational issues arose:
1) how to calculate f(D)
2) how to deal with ill-conditioning of V


The first issue can be dealt with by representing the 2-by-2 diagonal
blocks as complex numbers and requiring that f(z) be defined for
complex arguments:
f(x + iy) = g(x,y) + i h(x,y)

This result means that a 2-by-2 block, [x, y; -y, x], is transformed by
f() into a new 2-by-2 block, [g, h; -h, g]


Handling an ill-conditioned V is trickier. I started by considering a
perturbation matrix E such that:
1) the EVD of (A+E) is not ill-conditioned
2) ||E|| << ||A||
3) E.A = A.E

yes you're likely to get an ill conditioned V if your eigenvalues are close
together, or a singular V of course
if you have repeat eigenvalues, in which case you cannot form the eigen
decomposition.
have you thought about using a jordan decomposition instead of an eigen?

for the jordan decomposition then A = SJS^-1 where J is a Jordan block
matrix, and S a canonical basis
a function then would be f(A) = Sf(J)S^-1 . it is a little bit more
difficult to find f(J) than f(D) of an eigen
decomposition but not too much (see golub von loan for details on finding
matrix functions using the Jordan
decomposition)
Quote:

Utilizing such a perturbation results in a nice, simple Taylor
expansion
f(A+E) = f(A)/0! + E.{D}f(A)/1! + E^2.{D^2}f(A)/2! + O(E^3)


[[Beyond this point, I'm not completely satisfied with my reasoning,
although my Java application works just fine.]]


Finding a matrix which commutes with A was too difficult, so I
considered what would happen to the Taylor expansion if a small, random
matrix were used in two independent evaluations:
f(A+ E) = f(A)/0! + {D}f(E; A)/1! + {D^2}f(E; A)/2! + O(E^3)
f(A+2E) = f(A)/0! + {D}f(2E;A)/1! + {D^2}f(2E;A)/2! + O(E^3)


The non-commutivity of A with E makes the evaluation of {D^k}f(E;A)
terms very messy, but a scalar multiplier can be pulled clear of the
mess, i.e.

{D^k}f(s*E;A) = (s^k)*{D^k}f(E;A)


Letting e = E/||E||, we can say something like
f(A+ E) = f(A) + ||E||.{D}f(e;A) + ||E||^2.{D^2}f(e;A)/2 + ...
f(A+2E) = f(A) + ||2E||.{D}f(e;A) + ||2E||^2.{D^2}f(e;A)/2 + ...

without worrying over the commutation details of the individual terms.
The point is that the terms in each expansion are identical up to a
scale factor.

Subtraction yields the approximation
f(A) = 2f(A+E) - f(A+2E) + O(||E||^2)

If E is chosen such that
||E|| = sqrt(eps)*||A||

then the overall approximation should be accurate to O(eps).


Does anyone know of a nice expression for f(A+E) when A and E don't
commute?

Can this idea be extended to higher orders by evaluating f(A+E),
f(A+2E), f(A+3E), ... f(A+kE) and then combining them in such a way
that all terms cancel except for {D^k}f(E;A)/k! and f(A)?

~Glynne
Back to top
dan1112
science forum beginner


Joined: 15 Jul 2006
Posts: 3

PostPosted: Sat Jul 15, 2006 1:57 pm    Post subject: Re: Matrix functions via EVD decomposition Reply with quote

Using an eigenvalue decomposition is fine for nearly normal matrices,
in which case f(D) = diag(f(lambda_i), i=1..n). If the eigenbasis is
poorly conditioned, or if the matrix is not diagonalisable, calculate a
schur decomposition. There is a recurrence relationship known as
Parlett's recurrence (which is O(n^3) and has some problems which you
can look up). For other methods look up Golub and van Loan's book
"MAtrix Computions" which has a chapter devoted to matrix functions.
If the matrix is sparse, that's a whole other problem. The best known
way to compute f(A)b is by using Krylov subspace methods, however, this
is still a major area of research in numerical linear algebra.

To sumerise - if the eigenbasis is ill conditioned, don't use an
eigenvalue decomposition (even if you need to write your own code). On
a somewhat unrealated topic, why are you using Java - it's note really
the best language for scientific computing.
Back to top
Carl Barron
science forum beginner


Joined: 12 Jun 2005
Posts: 27

PostPosted: Sat Jul 15, 2006 9:23 am    Post subject: Re: Matrix functions via EVD decomposition Reply with quote

~Glynne <glynnec2002@yahoo.com> wrote:

Quote:
I was looking for a way to numerically evaluate functions of a real
matrix. The EigenvalueDecomposition class in the JAva MAtrix (JAMA)
package contained almost everything that was needed.

The decomposition produces matrices V and D such that the columns of V
contain the eigenvectors of A, while D contains the eigenvalues, i.e.
A.V = V.D

The eigenvalue matrix D is block diagonal with the real eigenvalues in
1-by-1 blocks and any complex eigenvalues, x+iy, in 2-by-2 blocks, [x,
y; -y, x].

Therefore, I should be able to compute f(A) using
f(A) = V.f(D)/V


Two computational issues arose:
1) how to calculate f(D)
2) how to deal with ill-conditioning of V


The first issue can be dealt with by representing the 2-by-2 diagonal
blocks as complex numbers and requiring that f(z) be defined for
complex arguments:
f(x + iy) = g(x,y) + i h(x,y)

This result means that a 2-by-2 block, [x, y; -y, x], is transformed by
f() into a new 2-by-2 block, [g, h; -h, g]


Handling an ill-conditioned V is trickier. I started by considering a
perturbation matrix E such that:
1) the EVD of (A+E) is not ill-conditioned
2) ||E|| << ||A||
3) E.A = A.E

Utilizing such a perturbation results in a nice, simple Taylor
expansion
f(A+E) = f(A)/0! + E.{D}f(A)/1! + E^2.{D^2}f(A)/2! + O(E^3)


[[Beyond this point, I'm not completely satisfied with my reasoning,
although my Java application works just fine.]]


Finding a matrix which commutes with A was too difficult, so I
considered what would happen to the Taylor expansion if a small, random
matrix were used in two independent evaluations:
f(A+ E) = f(A)/0! + {D}f(E; A)/1! + {D^2}f(E; A)/2! + O(E^3)
f(A+2E) = f(A)/0! + {D}f(2E;A)/1! + {D^2}f(2E;A)/2! + O(E^3)


The non-commutivity of A with E makes the evaluation of {D^k}f(E;A)
terms very messy, but a scalar multiplier can be pulled clear of the
mess, i.e.

{D^k}f(s*E;A) = (s^k)*{D^k}f(E;A)


Letting e = E/||E||, we can say something like
f(A+ E) = f(A) + ||E||.{D}f(e;A) + ||E||^2.{D^2}f(e;A)/2 + ...
f(A+2E) = f(A) + ||2E||.{D}f(e;A) + ||2E||^2.{D^2}f(e;A)/2 + ...

without worrying over the commutation details of the individual terms.
The point is that the terms in each expansion are identical up to a
scale factor.

Subtraction yields the approximation
f(A) = 2f(A+E) - f(A+2E) + O(||E||^2)

If E is chosen such that
||E|| = sqrt(eps)*||A||

then the overall approximation should be accurate to O(eps).


Does anyone know of a nice expression for f(A+E) when A and E don't
commute?

Can this idea be extended to higher orders by evaluating f(A+E),
f(A+2E), f(A+3E), ... f(A+kE) and then combining them in such a way
that all terms cancel except for {D^k}f(E;A)/k! and f(A)?

~Glynne
Good luck! consider A below it really has one eignevector and three

identical eigenvalues 2, given the V produce from TNT's JAMA below
V is essentially |0|*|1 1 1| , it should be singular there is
|0| only one eigenvector (0,0,1)^T
|1|

A =
2 0 0
1 2 0
0 1 2

V =
0 0 3.15544e-30
0 -1.77636e-15 -1.77636e-15
1 1 1

D =
2 0 0
0 2 0
0 0 2

AV =
0 0 6.31089e-30
0 -3.55271e-15 -3.55271e-15
2 2 2

VD =
0 0 6.31089e-30
0 -3.55271e-15 -3.55271e-15
2 2 2
Back to top
~Glynne
science forum beginner


Joined: 21 Nov 2005
Posts: 18

PostPosted: Sat Jul 15, 2006 5:51 am    Post subject: Matrix functions via EVD decomposition Reply with quote

I was looking for a way to numerically evaluate functions of a real
matrix. The EigenvalueDecomposition class in the JAva MAtrix (JAMA)
package contained almost everything that was needed.

The decomposition produces matrices V and D such that the columns of V
contain the eigenvectors of A, while D contains the eigenvalues, i.e.
A.V = V.D

The eigenvalue matrix D is block diagonal with the real eigenvalues in
1-by-1 blocks and any complex eigenvalues, x+iy, in 2-by-2 blocks, [x,
y; -y, x].

Therefore, I should be able to compute f(A) using
f(A) = V.f(D)/V


Two computational issues arose:
1) how to calculate f(D)
2) how to deal with ill-conditioning of V


The first issue can be dealt with by representing the 2-by-2 diagonal
blocks as complex numbers and requiring that f(z) be defined for
complex arguments:
f(x + iy) = g(x,y) + i h(x,y)

This result means that a 2-by-2 block, [x, y; -y, x], is transformed by
f() into a new 2-by-2 block, [g, h; -h, g]


Handling an ill-conditioned V is trickier. I started by considering a
perturbation matrix E such that:
1) the EVD of (A+E) is not ill-conditioned
2) ||E|| << ||A||
3) E.A = A.E

Utilizing such a perturbation results in a nice, simple Taylor
expansion
f(A+E) = f(A)/0! + E.{D}f(A)/1! + E^2.{D^2}f(A)/2! + O(E^3)


[[Beyond this point, I'm not completely satisfied with my reasoning,
although my Java application works just fine.]]


Finding a matrix which commutes with A was too difficult, so I
considered what would happen to the Taylor expansion if a small, random
matrix were used in two independent evaluations:
f(A+ E) = f(A)/0! + {D}f(E; A)/1! + {D^2}f(E; A)/2! + O(E^3)
f(A+2E) = f(A)/0! + {D}f(2E;A)/1! + {D^2}f(2E;A)/2! + O(E^3)


The non-commutivity of A with E makes the evaluation of {D^k}f(E;A)
terms very messy, but a scalar multiplier can be pulled clear of the
mess, i.e.

{D^k}f(s*E;A) = (s^k)*{D^k}f(E;A)


Letting e = E/||E||, we can say something like
f(A+ E) = f(A) + ||E||.{D}f(e;A) + ||E||^2.{D^2}f(e;A)/2 + ...
f(A+2E) = f(A) + ||2E||.{D}f(e;A) + ||2E||^2.{D^2}f(e;A)/2 + ...

without worrying over the commutation details of the individual terms.
The point is that the terms in each expansion are identical up to a
scale factor.

Subtraction yields the approximation
f(A) = 2f(A+E) - f(A+2E) + O(||E||^2)

If E is chosen such that
||E|| = sqrt(eps)*||A||

then the overall approximation should be accurate to O(eps).


Does anyone know of a nice expression for f(A+E) when A and E don't
commute?

Can this idea be extended to higher orders by evaluating f(A+E),
f(A+2E), f(A+3E), ... f(A+kE) and then combining them in such a way
that all terms cancel except for {D^k}f(E;A)/k! and f(A)?

~Glynne
Back to top
Google

Back to top
Display posts from previous:   
Post new topic   Reply to topic Page 1 of 1 [8 Posts] View previous topic :: View next topic
The time now is Mon Oct 23, 2017 11:20 am | All times are GMT
Forum index » Science and Technology » Math » num-analysis
Jump to:  

Similar Topics
Topic Author Forum Replies Last Post
No new posts Diagonalizable matrix aline Math 0 Wed Nov 29, 2006 3:08 am
No new posts Generating function for Mathieu functions cosmicstring@gmail.com Math 1 Fri Jul 21, 2006 8:39 am
No new posts Entire functions, polynomial bounds david petry Math 2 Thu Jul 20, 2006 11:09 pm
No new posts sign of the determinant of an augmented matrix? Mark Math 4 Thu Jul 20, 2006 1:30 am
No new posts spectrum of a symmetric tridiagonal random matrix pf.buonsante@gmail.com Math 0 Wed Jul 19, 2006 9:45 am

Copyright © 2004-2005 DeniX Solutions SRL
Other DeniX Solutions sites: Electronics forum |  Medicine forum |  Unix/Linux blog |  Unix/Linux documentation |  Unix/Linux forums  |  send newsletters
 


Powered by phpBB © 2001, 2005 phpBB Group
[ Time: 0.0323s ][ Queries: 20 (0.0047s) ][ GZIP on - Debug on ]