|
|
| Author |
Message |
~Glynne
science forum beginner
Joined: 21 Nov 2005
Posts: 18
|
 Posted: Sat Jul 15, 2006 5:51 am Post subject:
Matrix functions via EVD decomposition
|
|
|
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 |
|
 |
Carl Barron
science forum beginner
Joined: 12 Jun 2005
Posts: 27
|
| Posted: Sat Jul 15, 2006 9:23 am Post subject:
Re: Matrix functions via EVD decomposition
|
|
|
~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 |
|
|
dan
science forum beginner
Joined: 15 Jul 2006
Posts: 3
|
| Posted: Sat Jul 15, 2006 1:57 pm Post subject:
Re: Matrix functions via EVD decomposition
|
|
|
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 |
|
|
Jeremy Watts
science forum Guru Wannabe
Joined: 24 Mar 2005
Posts: 239
|
| Posted: Sat Jul 15, 2006 6:02 pm Post subject:
Re: Matrix functions via EVD decomposition
|
|
|
"~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 |
|
|
~Glynne
science forum beginner
Joined: 21 Nov 2005
Posts: 18
|
| Posted: Sat Jul 15, 2006 10:08 pm Post subject:
Re: Matrix functions via EVD decomposition
|
|
|
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 |
|
|
~Glynne
science forum beginner
Joined: 21 Nov 2005
Posts: 18
|
| Posted: Sat Jul 15, 2006 10:14 pm Post subject:
Re: Matrix functions via EVD decomposition
|
|
|
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
|
| Posted: Mon Jul 17, 2006 8:21 pm Post subject:
Re: Matrix functions via EVD decomposition
|
|
|
~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 |
|
|
dan
science forum beginner
Joined: 15 Jul 2006
Posts: 3
|
| Posted: Tue Jul 18, 2006 6:10 am Post subject:
Re: Matrix functions via EVD decomposition
|
|
|
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 |
|
|
Google
|
|
| Back to top |
|
|
|
|
The time now is Mon Dec 01, 2008 5:34 pm | All times are GMT
|
|
Personal Injury Attorney Los Angeles | Mortgage Calculator | Debt Help | Mobile Phone deals | Credit Card
|
|
Copyright © 2004-2005 DeniX Solutions SRL
|
|
Other DeniX Solutions sites:
Electronics forum |
Medicine forum |
Unix/Linux blog |
Unix/Linux documentation |
Unix/Linux forums
|
Powered by phpBB © 2001, 2005 phpBB Group
|
|
FAQ
Search
Memberlist
Usergroups
Register
Profile
Log in to check your private messages
Log in
Forum index
»
Science and Technology
»
Math
»
num-analysis
