# Problem with QRDecomposition [closed]

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?

• 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. – J. M.'s technical difficulties Nov 30 '19 at 5:36
• Even after 300000 times iteration (more than what was suggested), the results are not right. – AYBRXQD Nov 30 '19 at 5:42
• In MatLab after 30 times of iterations, I obtain very good results. – AYBRXQD Nov 30 '19 at 5:43
• @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. – J. M.'s technical difficulties Nov 30 '19 at 5:57
• It would appear that documentation was not consulted here. – Daniel Lichtblau Nov 30 '19 at 15:44

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.)

• 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. – AYBRXQD Nov 30 '19 at 6:06
• Correct, and this is in the documentation. – J. M.'s technical difficulties Nov 30 '19 at 6:10