Author 
Message 
dan1112 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
1by1 blocks and any complex eigenvalues, x+iy, in 2by2 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 illconditioning of V
The first issue can be dealt with by representing the 2by2 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 2by2 block, [x, y; y, x], is transformed by
f() into a new 2by2 block, [g, h; h, g]
Handling an illconditioned V is trickier. I started by considering a
perturbation matrix E such that:
1) the EVD of (A+E) is not illconditioned
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 noncommutivity 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: 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
1by1 blocks and any complex eigenvalues, x+iy, in 2by2 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 illconditioning of V
The first issue can be dealt with by representing the 2by2 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 2by2 block, [x, y; y, x], is transformed by
f() into a new 2by2 block, [g, h; h, g]
Handling an illconditioned V is trickier. I started by considering a
perturbation matrix E such that:
1) the EVD of (A+E) is not illconditioned
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 noncommutivity 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.15544e30
0 1.77636e15 1.77636e15
1 1 1
D =
2 0 0
0 2 0
0 0 2
AV =
0 0 6.31089e30
0 3.55271e15 3.55271e15
2 2 2
VD =
0 0 6.31089e30
0 3.55271e15 3.55271e15
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

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: 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
1by1 blocks and any complex eigenvalues, x+iy, in 2by2 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 illconditioning of V
The first issue can be dealt with by representing the 2by2 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 2by2 block, [x, y; y, x], is transformed by
f() into a new 2by2 block, [g, h; h, g]
Handling an illconditioned V is trickier. I started by considering a
perturbation matrix E such that:
1) the EVD of (A+E) is not illconditioned
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 noncommutivity 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.15544e30
0 1.77636e15 1.77636e15
1 1 1
D =
2 0 0
0 2 0
0 0 2
AV =
0 0 6.31089e30
0 3.55271e15 3.55271e15
2 2 2
VD =
0 0 6.31089e30
0 3.55271e15 3.55271e15
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

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
1by1 blocks and any complex eigenvalues, x+iy, in 2by2 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 illconditioning of V
The first issue can be dealt with by representing the 2by2 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 2by2 block, [x, y; y, x], is transformed by
f() into a new 2by2 block, [g, h; h, g]
Handling an illconditioned V is trickier. I started by considering a
perturbation matrix E such that:
1) the EVD of (A+E) is not illconditioned
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 noncommutivity 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

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 


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
1by1 blocks and any complex eigenvalues, x+iy, in 2by2 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 illconditioning of V
The first issue can be dealt with by representing the 2by2 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 2by2 block, [x, y; y, x], is transformed by
f() into a new 2by2 block, [g, h; h, g]
Handling an illconditioned V is trickier. I started by considering a
perturbation matrix E such that:
1) the EVD of (A+E) is not illconditioned
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 noncommutivity 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.15544e30
0 1.77636e15 1.77636e15
1 1 1
D =
2 0 0
0 2 0
0 0 2
AV =
0 0 6.31089e30
0 3.55271e15 3.55271e15
2 2 2
VD =
0 0 6.31089e30
0 3.55271e15 3.55271e15
2 2 2 

Back to top 


~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
1by1 blocks and any complex eigenvalues, x+iy, in 2by2 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 illconditioning of V
The first issue can be dealt with by representing the 2by2 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 2by2 block, [x, y; y, x], is transformed by
f() into a new 2by2 block, [g, h; h, g]
Handling an illconditioned V is trickier. I started by considering a
perturbation matrix E such that:
1) the EVD of (A+E) is not illconditioned
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 noncommutivity 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 Oct 23, 2017 11:20 am  All times are GMT

