QR decomposition using Gram Schmidt [closed]

Question

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?

Answer 1

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

Answer 2

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