Re: eig accuracy after similarity transformation

[Thomas Shores]

Rolf Fabian wrote:

> Hi,
> please consider the following example
>
> % dimension
> :) N=10;
>
> % we create an upper tri matrix with all N eigenvalues =1
> :) A=diag( ones(N,1) )+triu(randn(N),1)  |- ans e.g.
> A =
> 1.000  -0.706  -0.552  -0.154  -0.553   1.004   0.925  -1.900   0.951  1.138
> 0.000   1.000  -0.664   0.366  -0.373   0.899   1.526   0.924  -0.720  1.803
> 0.000   0.000   1.000  -0.697   1.062   1.352   1.215   0.298   0.097 -1.577
> 0.000   0.000   0.000   1.000   0.268   0.614   1.267  -0.401  -1.363  1.022
> 0.000   0.000   0.000   0.000   1.000   0.792  -0.041   1.094   0.917  0.554
> 0.000   0.000   0.000   0.000   0.000   1.000  -0.928   0.608   0.140  0.233
> 0.000   0.000   0.000   0.000   0.000   0.000   1.000   0.518  -0.384 -0.712
> 0.000   0.000   0.000   0.000   0.000   0.000   0.000   1.000  -0.042  0.192
> 0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   1.000 -1.844
> 0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000  1.000
>
> % all diagonal elements (eigenval) are 1 exactly
> :) eig(A).'             |- ans =
>   1  1  1  1  1  1  1  1  1  1
>
> % let's create a random unitary transformation matrix RU
> :) [ RU,dummy ] = schur( randn(N) );
> :) max(max(abs(RU*RU'-eye(N)))) |- ans = 6.6613e-16
> %  RU*RU' == Identity, hence RU is definitely unitary
>
> % now we look at calculated eigenval of the
> % unitarily similarity transformed matrix A
> :) eig( RU'*A*RU )      |- ans =
>   1.01319 + 0.00000i
>   1.01051 + 0.00788i
>   1.01051 - 0.00788i
>   1.00368 + 0.01245i
>   1.00368 - 0.01245i
>   0.99570 + 0.01200i
>   0.99570 - 0.01200i
>   0.98972 + 0.00715i
>   0.98972 - 0.00715i
>   0.98758 + 0.00000i
>
> What's the origin of this extremely large differences (>1%)
> to the expected eigenvalues (all 1) ?
>
> If matrix A is constructed with different eigenvalues
> only, the analogous differences of eig( RU'*A*RU )
> are observed to be smaller by many orders of magnitude!
>
> Rolf
>
> -------------------------------------------------------------
> Octave is freely available under the terms of the GNU GPL.
>
> Octave's home on the web:  http://www.octave.org
> How to fund new projects:  http://www.octave.org/funding.html
> Subscription information:  http://www.octave.org/archive.html
> -------------------------------------------------------------

Your question is a classic in numerical linear algebra (NLA), and it has a
classic answer. You might consult Golub and Van Loan's text "Matrix
Computations" for more details.  In a nutshell here is the reason for the
problem: by the very act of disguising your A by the conjugation, you
introduced error into the problem , addmittedly very small of order  eps =
1e-15.  So look at it this way.  Suppose that the *only* error that occurs is
the eps you introduced, and thereafter you (or Octave) is going to do exact
arithmetic (of course, this is too optimistic) on the purturbed matrix A~.
There is a standard result in NLA that says that if you have a defective
eigenvalue (that is, its geometric multiplicity is less than its algebraic
multiplicity), say belonging to a Jordan block of order p, than the discrepency
between the exact eigenvalues of the purturbed matrix and the eigenvalues of
the original matrix is of order eps^(1/p).  Now p can be as large as N, the
size of the matrix. This means that you could expect a deviation of the order
of magnitude of 1e-1.5 in youy example.  This explains the error you saw.  It
also underscores the fact that the problem is not with Octave, but with the
more fundamental  and unavoidable fact that only finite precision arithmetic is
performed.  BTW, your example is most definitely defective: the only
nondefective upper triangular matrix with ones down the diagonal is the
identity matrix.  Hope this helps.

Tom Shores




-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------