# QR decomposition using Gram Schmidt [closed]

First of all I don't want to use the reflection method, only Gram-Schmidt.

So here's the program for Gram-Schmidt:

GS[A_] := Module[{u = {}, col, e = {}, ae, t},
col = Transpose[A];(* list with the column vectors of A *)
u = Append[u, col[[1]]];
e = AppendTo[e, u[[1]]/Norm[u[[1]]]];
For[i = 2, i <= Length[A], i++, ae = 0; t = 1;
While[t <= i - 1, ae = ae - (col[[i]].e[[t]])*e[[t]]; t++];
u = AppendTo[u, col[[i]] + ae];
e = Append[e, u[[i]]/Norm[u[[i]]]]]; {u, e}]


With this simple matrix:

B = {{1, 3, 1}, {2, 2, 1}, {3, 2, 3}};


it gives

{{{1, 2, 3}, {29/14, 1/7, -(11/14)}, {14/69, -(49/69), 28/69}}, {{1/
Sqrt[14], Sqrt[2/7], 3/Sqrt[14]}, {29/Sqrt[966], Sqrt[2/
483], -(11/Sqrt[966])}, {2/Sqrt[69], -(7/Sqrt[69]), 4/Sqrt[69]}}}


Now the code for QR decomposition is:

QR[A_] := Module[{Ak, i, Q, R = {}, a = {}, col, k},
Ak = FullSimplify[GS[A]];
col = Transpose[A];(* list with the column vectors of A *)
Q = MatrixForm[Transpose[Ak[[2]]]];
For[i = 1, i <= Length[Ak[[2]]], i++, a = {}; k = 1;
While[k <= Length[col], a = AppendTo[a, col[[k]]. Ak[[2]][[i]]];
k++];
R = AppendTo[R, a]];
R = UpperTriangularize[R]; {Q, R}]


and it works for a square matrix but it doesn't work for a non square matrix, unless I change my GS code to i< Length[A]. What do I have to do for this code work in both cases?

• You should write : k = 1 instead of k == 1. The latter is a comparison. Dec 5, 2020 at 15:21
• Thank you so much Daniel Huber! Dec 5, 2020 at 15:59
• Tip 1: Don't use MatrixForm the way you do. Omit it. You can use it later, like MatrixForm /@ QR[B], to make the output pretty. Dec 5, 2020 at 19:44
• Tip 2: Length[A] in GS[] should be Length[col]. Then it works on B = {{1, 3, 1}, {2, 2, 1}, {3, 2, 3}, {4, 5, 6}}. Dec 5, 2020 at 19:46

These methods are all built in:

B = {{1, 3, 1}, {2, 2, 1}, {3, 2, 3}};

Orthogonalize[B, Method -> "GramSchmidt"]
(*    {{1/Sqrt[11], 3/Sqrt[11], 1/Sqrt[11]},
{13/(3 Sqrt[22]), -5/(3 Sqrt[22]), Sqrt[2/11]/3},
{-1/(3 Sqrt[2]), -1/(3 Sqrt[2]), 2 Sqrt[2]/3}}    *)

QRDecomposition[B, Pivoting -> False]
(*    {{{1/Sqrt[14], Sqrt[2/7], 3/Sqrt[14]},
{29/Sqrt[966], Sqrt[2/483], -11/Sqrt[966]},
{2/Sqrt[69], -7/Sqrt[69], 4/Sqrt[69]}},
{{Sqrt[14], 5 Sqrt[2/7] + 3/Sqrt[14], 6 Sqrt[2/7]},
{0, Sqrt[69/14], Sqrt[2/483] - 11 Sqrt[3/322] + 29/Sqrt[966]},
{0, 0, 7/Sqrt[69]}}}                                              *)

• Thanks but it´s for a project so I have to compute them... Dec 5, 2020 at 19:05

Your code seems basically to work after fixing a few typos (see the comments above), so I'll share my take on your problem:

B = {{1, 3, 1}, {2, 2, 1}, {3, 2, 3}, {4, 5, 6}};
mat = Transpose@ B;
qmat = Transpose@
Fold[ (* G-S *)
Append[#1, #2 - #1.#2.#1 // Normalize] &,
{Normalize@ First@ mat},
Take[mat, {2, UpTo@ Length@B}]
];
rmat = Transpose[qmat].B;

OrthogonalMatrixQ@ qmat
UpperTriangularMatrixQ@ rmat
B == qmat.rmat // Simplify

(*
True  <-- Q orthogonal (code assumes B has real entries)
True  <-- R upper-triangular
True  <-- B == Q.R
*)


If you're studying numerical linear algebra, you can explore the precision loss:

SeedRandom[1]; (* results change slightly with a different seed *)
wp = 32;  (* initial precision 32 digits *)
B = RandomReal[10, {10, 50}, WorkingPrecision -> wp];
mat = Transpose@B;
qmat = Transpose@
Fold[ (* G-S *)
Append[#1, #2 - #1.#2.#1 // Normalize] &,
{Normalize@First@mat},
Take[mat, {2, UpTo@ Length@B}]
];
rmat = Transpose[qmat].B;

OrthogonalMatrixQ@qmat
UpperTriangularMatrixQ@rmat
B == qmat.rmat
qmat.rmat // Precision
B - SetPrecision[qmat.rmat, wp] // Abs // Max
(*
True  <-- Q orthogonal (code assumes B has real entries)
True  <-- R upper-triangular
True  <-- B == Q.R
14.7916               <-- Precision (up to 17.2 digits lost)
5.29635583648*10^-20  <-- Accuracy (19 dec. pl.)
*)