1
$\begingroup$

I was using Mathematica for the QR decomposition method. But I got strange results. I wanted to find eigenvalues of a matrix, say

A = {{1, 1, 1, 1, 1},
     {16, 8, 4, 2, 1},
     {81, 27, 9, 3, 1},
     {256, 64, 16, 4, 1},
     {625, 125, 25, 5, 1}};

To obtain eigenvalues using the QR decomposition method, I decompose the matrix to {q, r} matrices, then in an iterative manner, I calculate A = r.q. After enough iterative steps, the diagonal elements of A would be the eigenvalues.

For example, let me do it 30 times, so I can use the following code

For[i = 1, i <= 30, i++,
  Clear[qq, rr];
  {qq, rr} = QRDecomposition[N[A]];
  A = rr.qq;];

TableForm[A]

gives

   -13.8505           279.402         -319.744           106.778       -539.532
    0.827919        23.8641         -21.5706           9.68315       -15.3848
    2.21191*10^-30  6.3774*10^-29    1.05619          -1.3251         12.3686
    6.79367*10^-38  1.2065*10^-36   -0.954018         -0.472083      -4.47852
    2.2827*10^-57   3.83204*10^-56  -1.42879*10^-20   -7.2897*10^-21  0.290785

But

Eigenvalues[N[A]]

gives

{45.9604, -28.9437, 6.79539, -0.849619, 0.0374986}

which is different than the diagonal elements obtained in the iterative process. I did the same process in MATLAB and I obtained exactly the same results as I got using Eigenvalues. When I checked, I find that Mathematica calculates the QR decomposition differently than MATLAB. For instance, q matrix calculated by Mathematica is

   -0.00146967    -0.0235147  -0.119043   -0.376235   -0.918543
   -0.0485342   -0.285599   -0.617372   -0.642007   0.350366
    0.421646    0.708267    0.236505    -0.491504   0.151863
   -0.848679    0.133537    0.387736    -0.324186   0.0804755
   -0.315597    0.631194    -0.631194   0.315597    -0.0631194

but MATLAB calculated the same matrix as

  -0.0015   -0.0485    0.4216   -0.8487   -0.3156
  -0.0235   -0.2856    0.7083    0.1335    0.6312
  -0.1190   -0.6174    0.2365    0.3877   -0.6312
  -0.3762   -0.6420   -0.4915   -0.3242    0.3156
  -0.9185    0.3504    0.1519    0.0805   -0.0631

I'm sure that MATLAB is doing the computation correctly because it derives the same eigenvalues as Mathematica's Eigenvalues did.

I have this question: Does my results mean that QRDecomposition is working wrong?

$\endgroup$
  • 1
    $\begingroup$ Do not use the bugs tag unless other people have confirmed that what you encountered is a bug. (In this case, the error comes from a less-than-careful reading of the documentation.) I'll write something more detailed in a little bit. $\endgroup$ – J. M.'s technical difficulties Nov 30 '19 at 5:36
  • $\begingroup$ Even after 300000 times iteration (more than what was suggested), the results are not right. $\endgroup$ – AYBRXQD Nov 30 '19 at 5:42
  • $\begingroup$ In MatLab after 30 times of iterations, I obtain very good results. $\endgroup$ – AYBRXQD Nov 30 '19 at 5:43
  • $\begingroup$ @bill, I can't speak for the convergence of your proposal, but the OP's original implementation is manifestly not a similarity transformation, which is why the expected results are not seen. $\endgroup$ – J. M.'s technical difficulties Nov 30 '19 at 5:57
  • 2
    $\begingroup$ It would appear that documentation was not consulted here. $\endgroup$ – Daniel Lichtblau Nov 30 '19 at 15:44
4
$\begingroup$

See if you can spot the difference between Mathematica and MATLAB:

{qq, rr} = QRDecomposition[N[{{1, 2, 3},
                              {1, 4, 9},
                              {1, 8, 27}}]]
   {{{-0.57735, -0.57735, -0.57735},
     {0.617213, 0.154303, -0.771517},
     {0.534522, -0.801784, 0.267261}},
    {{-1.73205, -8.0829, -22.5167},
     {0., -4.32049, -17.5906},
     {0., 0., 1.60357}}}

Transpose[qq].rr
   {{1., 2., 3.}, {1., 4., 9.}, {1., 8., 27.}}

mat = [1 2 3; 1 4 9; 1 8 27];

[qq, rr] = qr(mat)

qq =

   -0.5774    0.6172    0.5345
   -0.5774    0.1543   -0.8018
   -0.5774   -0.7715    0.2673


rr =

   -1.7321   -8.0829  -22.5167
         0   -4.3205  -17.5906
         0         0    1.6036

qq*rr

ans =

    1.0000    2.0000    3.0000
    1.0000    4.0000    9.0000
    1.0000    8.0000   27.0000

With that in mind, here is a compact way to perform the QR iteration in Mathematica:

mat = N[{{1, 1, 1, 1, 1},
         {16, 8, 4, 2, 1},
         {81, 27, 9, 3, 1},
         {256, 64, 16, 4, 1},
         {625, 125, 25, 5, 1}}];

Nest[With[{qr = QRDecomposition[#]},
          Apply[Dot, Reverse[MapAt[Transpose, qr, 1]]]] &, mat, 30]
   {{45.9605, 333.072, -514.756, 307.72, -66.7936},
    {-0.0000205211, -28.9438, 62.9179, -47.4505, 12.6878},
    {-8.60691*10^-24, -1.6602*10^-17, 6.79539, -9.37454, 3.58764},
    {-1.30163*10^-50, -2.94067*10^-44, 1.93237*10^-26, -0.849619, 0.647286},
    {-3.31603*10^-91, -8.34873*10^-85, 6.99777*10^-67, -5.11405*10^-41, 0.0374986}}

Sort[Diagonal[%], LessEqual]
   {-28.9438, -0.849619, 0.0374986, 6.79539, 45.9605}

Sort[Eigenvalues[mat], LessEqual]
   {-28.9437, -0.849619, 0.0374986, 6.79539, 45.9604}

(Use NestList[] instead of Nest[] if you want to see the individual iterates.)

| improve this answer | |
$\endgroup$
  • 2
    $\begingroup$ Dear @J. M., I think I find the solution, q matrix in the Mathematica must be transformed to obtain correct results. Please check this point. When you use {qq, rr} = QRDecomposition[N[A]] , you can not obtain A by A = qq.rr. But instead you can obtain A by A = Transpose[qq].rr. Please check this point if I'm wrong. $\endgroup$ – AYBRXQD Nov 30 '19 at 6:06
  • 2
    $\begingroup$ Correct, and this is in the documentation. $\endgroup$ – J. M.'s technical difficulties Nov 30 '19 at 6:10

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