4
$\begingroup$

Note: From the comments submitted below by other StackExchange users it appears that this apparent bug was fixed sometime after version 7.0.1.

I asked a vague variant of this question a few days ago, but now I will give a concrete example.

Floating point matrices with repeated eigenvalues often arise as the solutions to optimization problems.

The Mathematica functions MatrixPower and Eigensystem appear to have sporadic difficulty dealing with such matrices. (It is not clear to me whether this is because of numerical instability when trying to determine whether such matrices have nontrivial Jordan blocks, or if there is some other problem.) I have had to write my own routines for matrix power (and more generally for the functional calculus) to handle such cases, but maybe someone knows a quick fix.

In the notebook

http://tinyurl.com/lhodpdv

I generated the following matrix by random search: (I'm using Mathematica 7.0.1 on 32-bit Windows)

 {{0.461964, -0.314388, -0.138647, -2.23115*10^-6}, 
 {-0.314388, 0.510134, -0.128432, -2.06676*10^-6}, 
 {-0.138647, -0.128432, 0.744719, -9.11454*10^-7}, 
{-2.23115*10^-6, -2.06676*10^-6, -9.11454*10^-7, 0.801359}}

This matrix was defined by taking the average of a matrix and its transpose, so it is exactly self-adjoint. (I searched for it by finding rare matrices for which my usual code for inverting the Eigensystem function acting on a Hermitian matrix failed.)

The list of eigenvalues

{0.801359, 0.801359, 0.801359, 0.1141}

has three repeated elements. All eigenvalues are positive, so Mathematica should have no trouble computing the unique positive semidefinite square root of this matrix.

However, using MatrixPower[%, 1/2], Mathematica reports that the square root is

{0.690223, -0.170008, -0.477225, 0.},
{-0.189862, 0.737704, -0.442063, 0.},
{-0.0837305, -0.0694507, 0.700234, 0.},
{-1.34742*10^-6, -1.11762*10^-6, -3.13724*10^-6, 0.895186}

This matrix is not even symmetric. Much worse, Mathematica computes the matrix product

(square root).(square root)

as 

{{0.548644, -0.209616, -0.588406, 0.}, 
{-0.234095, 0.607187, -0.545053, 0.}, 
{-0.103238, -0.0856309, 0.560987, 0.}, 
{-1.66133*10^-6, -1.37799*10^-6, -3.86814*10^-6, 0.801359}},

which is nowhere near the original matrix:

{{0.461964, -0.314388, -0.138647, -2.23115*10^-6}, 
{-0.314388, 0.510134, -0.128432, -2.06676*10^-6}, 
{-0.138647, -0.128432, 0.744719, -9.11454*10^-7}, 
{-2.23115*10^-6, -2.06676*10^-6, -9.11454*10^-7, 0.801359}}
Improve this question
$\endgroup$
21
  • $\begingroup$ If numerical stability is an issue the bug may not reproduce if you cut and past the numbers, although the linked notebook can be used to reproduce them. $\endgroup$
    – J Tyson
    Commented Jul 20, 2013 at 8:09
  • $\begingroup$ I just let Mathematica 9 run through more 4 million iterations of your code, without finding a example where it fails. Although it is not indicated in the docs, MatrixPower must have changed between 8 and 9, as the handling of non positive definite matrices is now different. Maybe it was a bug and it has been fixed? $\endgroup$
    – sebhofer
    Commented Jul 20, 2013 at 9:54
  • 1
    $\begingroup$ Rojo: I should add that taking a square root of a positive semidefinite matrix is a common operation, for example it is used when one computes the polar decomposition of an arbitrary matrix. More generally, if f:R->R is a function from the reals to the reals and A is a self-adjoint matrix, then one defines f(A) using the so-called "functional calculus" as follows: Take a spectral decomposition of A, with orthonormal eigenvectors. Then f(A) is computed by replacing each eigenvalue lamda_i by f(lamda_i). Note that if f=P is a polynomial then f(A) agrees with the usual notion of P(A). $\endgroup$
    – J Tyson
    Commented Jul 21, 2013 at 15:42
  • 1
    $\begingroup$ @Rojo no; in fact, I think it should be documented specifically if valid exponents are limited to the integers, but this doesn't appear anywhere in the documentation, and no message is produced if non-integers are used. The default position for all Mathematica functions seems to be that they can accept any type of input unless otherwise noted, and the only restriction given is that the matrix must be square. To lend additional weight to this argument, the problem is fixed in later versions, but the documentation isn't updated, and MatrixPower is listed as "last updated in 6". $\endgroup$
    – Oleksandr R.
    Commented Jul 28, 2013 at 12:43
  • 2
    $\begingroup$ @Rojo I understand your point of view. Although I mostly disagree, without getting into a detailed discussion, I think it's better to remove the tag. Although I'd personally still consider it a (very likely if not absolutely certain) bug, since the tag is meant to be used by community consensus, I can't justify adding it on my own initiative when others legitimately disagree. Incidentally, the matrix square root has been studied widely, and I agree with OP that this is likely just a numerical stability issue, perhaps due to special-casing this computation to avoid using the more general SVD. $\endgroup$
    – Oleksandr R.
    Commented Jul 29, 2013 at 4:08

1 Answer 1

Reset to default
2
$\begingroup$

I consider this question to be answered by the comments above that this bug in MATHEMATICA 7 was fixed in later versions.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.