# Why is Mathematica eating a row from QRDecomposition

I'm attempting to calculate the QRDecomposition of the following 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?

• 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. Aug 1, 2015 at 4:38

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}}

• 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(). Aug 1, 2015 at 13:20
• 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. Aug 2, 2015 at 12:30
• 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[]. Aug 2, 2015 at 12:54
• 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… Aug 2, 2015 at 13:43

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}} *)