12
$\begingroup$

I'm attempting to calculate the QRDecomposition of the matrix:

a = {{1, 3}, {0, 5}, {2, -8}}

QRDecomposition[a]

The answer Mathematica gives me is:
- $\mathbf Q$ which is a $2\times3$ matrix
- $\mathbf R$ which is a $2\times2$ matrix

The actual answer that I got using a number of QR decomposition calculators online is a $3\times3$ matrix for $\mathbf Q$ and $3\times2$ matrix for $\mathbf R$. The TA of my class told me that Mathematica gives $\mathbf Q^\top$ instead of $\mathbf Q$ so I should just transpose the answer, but clearly $2\times3$ transposed is not $3\times3$. Any ideas?

$\endgroup$
  • 6
    $\begingroup$ QRDecomposition[] is computing what is called a "thin" or "economy" QR, where the orthonormal factor inherits the dimensions of the rectangular matrix. There are relations with this and "full QR" that you can use, however. Search around. $\endgroup$ – J. M. will be back soon Aug 1 '15 at 4:38
15
$\begingroup$

As I have previously noted, QRDecomposition[] is by default set to return the so-called "thin QR" or "economy QR" decomposition; this is often the form desired in applications, since the triangular factor does not have the unneeded zero rows. MATLAB's qr(), by contrast, returns the full QR decomposition by default, and the economy QR through an option setting.

Here, then, is a general routine for producing the full QR decomposition, adapted from this old linear algebra tutorial. I have elected to also extend this routine to the case where QR with column pivoting is wanted (sometimes a must in certain least squares applications):

Options[FullQRDecomposition] = Options[QRDecomposition];
FullQRDecomposition[mat_?MatrixQ, opts : OptionsPattern[FullQRDecomposition]] := 
      Module[{prec = Precision[mat], dec},  
             dec = QRDecomposition[mat, Sequence @@ FilterRules[{opts} ~Join~
                                   Options[FullQRDecomposition],
                                   Options[QRDecomposition]]];
             dec = MapAt[PadRight[#, Dimensions[mat], N[0, {prec, prec}]] &,
                         dec, 2];
             dec = MapAt[Orthogonalize[Join[#, NullSpace[#]]] &, dec, 1];
             dec]

Using the OP's example:

a = {{1., 3.}, {0., 5.}, {2., -8.}};

{qq, rr} = FullQRDecomposition[a];
Transpose[qq].rr - a // Chop
   {{0, 0}, {0, 0}, {0, 0}}

{qq, rr, piv} = FullQRDecomposition[a, Pivoting -> True];
Transpose[qq].rr - a.piv // Chop
   {{0, 0}, {0, 0}, {0, 0}}
$\endgroup$
  • 1
    $\begingroup$ P.S. As a point of contrast, SingularValueDecomposition[] returns the "full SVD" by default (although one can use the second argument to yield the "thin SVD") in the same manner as MATLAB's svd(). $\endgroup$ – J. M. will be back soon Aug 1 '15 at 13:20
  • $\begingroup$ Your qq does not appear to be a unitary matrix; i.e., the third row of qq, although orthogonal to the other two rows, is not normalized. Thus, Transpose[qq].qq does not yield the identity matrix. This is, of course, easy to fix. $\endgroup$ – bbgodfrey Aug 2 '15 at 12:30
  • $\begingroup$ Huh… what results are you getting? I did test qq for orthogonality with and without pivoting, and I obtain something that looks like the identity matrix after a Chop[]. $\endgroup$ – J. M. will be back soon Aug 2 '15 at 12:54
  • $\begingroup$ Ah, I see what you mean now. The discrepancy is due to NullSpace[] using SVD in the inexact case, but Gaussian elimination otherwise. Let me fix that… $\endgroup$ – J. M. will be back soon Aug 2 '15 at 13:43
11
$\begingroup$

The relationship between Q and R as computed by QRDecomposition and the "full QR" results (as described by J. M.} can be found in, for instance, Wikipedia. The following illustrates how to go from the Mathematica to the Wikipedia formulation. With a as defined in the question,

{q, r} = QRDecomposition[a]
(* {{{1/Sqrt[5], 0, 2/Sqrt[5]}, 
     {28/Sqrt[1605], 5 Sqrt[5/321], -(14/Sqrt[1605])}}, 
    {{Sqrt[5], -(13/Sqrt[5])}, {0, Sqrt[321/5]}}} *)

To obtain the R that you are seeking, add a third row of zeroes.

rr = Join[r, {{0, 0}}]
(* {{Sqrt[5], -(13/Sqrt[5])}, {0, Sqrt[321/5]}, {0, 0}} *)

To obtain the full Q, one must determine a third row, orthonormal to the first two, and append it.

Solve[{q.{q1, q2, q3} == 0, q1^2 + q2^2 + q3^2 == 1}, {q1, q2, q3}][[1]];
qq = Join[q, {{q1, q2, q3}}]/.%
(* {{1/Sqrt[5], 0, 2/Sqrt[5]}, 
    {28/Sqrt[1605], 5 Sqrt[5/321], -(14/Sqrt[1605])}, 
    {-(10/Sqrt[321]), 14/Sqrt[321], 5/Sqrt[321]}} *)

To prove that the decomposition is correct, recover a

Transpose[qq].rr
(* {{1, 3}, {0, 5}, {2, -8}} *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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