QR decomposition

QR decomposition
You are encouraged to solve this task according to the task description, using any language you may know.

Any rectangular $m\times n$ matrix ${\mathit {A}}$ can be decomposed to a product of an orthogonal matrix ${\mathit {Q}}$ and an upper (right) triangular matrix ${\mathit {R}}$, as described in QR decomposition.

Task

Demonstrate the QR decomposition on the example matrix from the Wikipedia article:

$A={\begin{pmatrix}12&-51&4\\6&167&-68\\-4&24&-41\end{pmatrix}}$

and the usage for linear least squares problems on the example from Polynomial regression. The method of Householder reflections should be used:

Method

Multiplying a given vector ${\mathit {a}}$, for example the first column of matrix ${\mathit {A}}$, with the Householder matrix ${\mathit {H}}$, which is given as

$H=I-{\frac {2}{u^{T}u}}uu^{T}$

reflects ${\mathit {a}}$ about a plane given by its normal vector ${\mathit {u}}$. When the normal vector of the plane ${\mathit {u}}$ is given as

$u=a-\|a\|_{2}\;e_{1}$

then the transformation reflects ${\mathit {a}}$ onto the first standard basis vector

$e_{1}=[1\;0\;0\;...]^{T}$

which means that all entries but the first become zero. To avoid numerical cancellation errors, we should take the opposite sign of $a_{1}$:

$u=a+{\textrm {sign}}(a_{1})\|a\|_{2}\;e_{1}$

and normalize with respect to the first element:

$v={\frac {u}{u_{1}}}$

The equation for $H$ thus becomes:

$H=I-{\frac {2}{v^{T}v}}vv^{T}$

or, in another form

$H=I-\beta vv^{T}$

with

$\beta ={\frac {2}{v^{T}v}}$

Applying ${\mathit {H}}$ on ${\mathit {a}}$ then gives

$H\;a=-{\textrm {sign}}(a_{1})\;\|a\|_{2}\;e_{1}$

and applying ${\mathit {H}}$ on the matrix ${\mathit {A}}$ zeroes all subdiagonal elements of the first column:

$H_{1}\;A={\begin{pmatrix}r_{11}&r_{12}&r_{13}\\0&*&*\\0&*&*\end{pmatrix}}$

In the second step, the second column of ${\mathit {A}}$, we want to zero all elements but the first two, which means that we have to calculate ${\mathit {H}}$ with the first column of the submatrix (denoted *), not on the whole second column of ${\mathit {A}}$.

To get $H_{2}$, we then embed the new ${\mathit {H}}$ into an $m\times n$ identity:

$H_{2}={\begin{pmatrix}1&0&0\\0&H&\\0&&\end{pmatrix}}$

This is how we can, column by column, remove all subdiagonal elements of ${\mathit {A}}$ and thus transform it into ${\mathit {R}}$.

$H_{n}\;...\;H_{3}H_{2}H_{1}A=R$

The product of all the Householder matrices ${\mathit {H}}$, for every column, in reverse order, will then yield the orthogonal matrix ${\mathit {Q}}$.

$H_{1}H_{2}H_{3}\;...\;H_{n}=Q$

The QR decomposition should then be used to solve linear least squares (Multiple regression) problems ${\mathit {A}}x=b$ by solving

$R\;x=Q^{T}\;b$

When ${\mathit {R}}$ is not square, i.e. $m>n$ we have to cut off the ${\mathit {m}}-n$ zero padded bottom rows.

$R={\begin{pmatrix}R_{1}\\0\end{pmatrix}}$

and the same for the RHS:

$Q^{T}\;b={\begin{pmatrix}q_{1}\\q_{2}\end{pmatrix}}$

Finally, solve the square upper triangular system by back substitution:

$R_{1}\;x=q_{1}$

Ada

Output matches that of Matlab solution, not tested with other matrices.

with Ada.Text_IO; use Ada.Text_IO;
with Ada.Numerics.Real_Arrays; use Ada.Numerics.Real_Arrays;
with Ada.Numerics.Generic_Elementary_Functions;
procedure QR is

   procedure Show (mat : Real_Matrix) is
      package FIO is new Ada.Text_IO.Float_IO (Float);
   begin
      for row in mat'Range (1) loop
         for col in mat'Range (2) loop
            FIO.Put (mat (row, col), Exp => 0, Aft => 4, Fore => 5);
         end loop;
         New_Line;
      end loop;
   end Show;

   function GetCol (mat : Real_Matrix; n : Integer) return Real_Matrix is
      column : Real_Matrix (mat'Range (1), 1 .. 1);
   begin
      for row in mat'Range (1) loop
         column (row, 1) := mat (row, n);
      end loop;
      return column;
   end GetCol;

   function Mag (mat : Real_Matrix) return Float is
      sum : Real_Matrix := Transpose (mat) * mat;
      package Math is new Ada.Numerics.Generic_Elementary_Functions
         (Float);
   begin
      return Math.Sqrt (sum (1, 1));
   end Mag;

   function eVect (col : Real_Matrix; n : Integer) return Real_Matrix is
      vect : Real_Matrix (col'Range (1), 1 .. 1);
   begin
      for row in col'Range (1) loop
         if row /= n then vect (row, 1) := 0.0;
         else vect (row, 1) := 1.0; end if;
      end loop;
      return vect;
   end eVect;

   function Identity (n : Integer) return Real_Matrix is
      mat : Real_Matrix (1 .. n, 1 .. n) := (1 .. n => (others => 0.0));
   begin
      for i in Integer range 1 .. n loop mat (i, i) := 1.0; end loop;
      return mat;
   end Identity;

   function Chop (mat : Real_Matrix; n : Integer) return Real_Matrix is
      small : Real_Matrix (n .. mat'Length (1), n .. mat'Length (2));
   begin
      for row in small'Range (1) loop
         for col in small'Range (2) loop
            small (row, col) := mat (row, col);
         end loop;
      end loop;
      return small;
   end Chop;

   function H_n (inmat : Real_Matrix; n : Integer)
      return Real_Matrix is
      mat : Real_Matrix := Chop (inmat, n);
      col : Real_Matrix := GetCol (mat, n);
      colT : Real_Matrix (1 .. 1, mat'Range (1));
      H : Real_Matrix := Identity (mat'Length (1));
      Hall : Real_Matrix := Identity (inmat'Length (1));
   begin
      col := col - Mag (col) * eVect (col, n);
      col := col / Mag (col);
      colT := Transpose (col);
      H := H - 2.0 * (col * colT);
      for row in H'Range (1) loop
         for col in H'Range (2) loop
            Hall (n - 1 + row, n - 1 + col) := H (row, col);
         end loop;
      end loop;
      return Hall;
   end H_n;

   A : constant Real_Matrix (1 .. 3, 1 .. 3) := (
      (12.0, -51.0, 4.0),
      (6.0, 167.0, -68.0),
      (-4.0, 24.0, -41.0));
   Q1, Q2, Q3, Q, R: Real_Matrix (1 .. 3, 1 .. 3);
begin
   Q1 := H_n (A, 1);
   Q2 := H_n (Q1 * A, 2);
   Q3 := H_n (Q2 * Q1* A, 3);
   Q := Transpose (Q1) * Transpose (Q2) * TransPose(Q3);
   R := Q3 * Q2 * Q1 * A;
   Put_Line ("Q:"); Show (Q);
   Put_Line ("R:"); Show (R);
end QR;

Output:

Q:
    0.8571   -0.3943   -0.3314
    0.4286    0.9029    0.0343
   -0.2857    0.1714   -0.9429
R:
   14.0000   21.0000  -14.0000
   -0.0000  175.0000  -70.0000
   -0.0000    0.0000   35.0000

ATS

Perhaps not every template function that was written below is actually used. Much of what is below amounts to a little library for working with matrices. To treat blocks and transposes as matrices themselves, I use a trick employed in some Scheme implementations of matrices: indices are mapped by closures, and the closures can nested.

%{^
#include <math.h>
#include <float.h>
%}

#include "share/atspre_staload.hats"

macdef NAN = g0f2f ($extval (float, "NAN"))
macdef Zero = g0i2f 0
macdef One = g0i2f 1
macdef Two = g0i2f 2

(* g0float_sqrt is available in the ats2-xprelude package, but let us
   quickly add it here, with implementations for the g0float types
   included in the prelude. *)
extern fn {tk : tkind} g0float_sqrt : g0float tk -<> g0float tk
overload sqrt with g0float_sqrt
implement g0float_sqrt<fltknd> x = $extfcall (float, "sqrtf", x)
implement g0float_sqrt<dblknd> x = $extfcall (double, "sqrt", x)
implement g0float_sqrt<ldblknd> x = $extfcall (ldouble, "sqrtl", x)

(* Similarly for g0float_copysign. *)
extern fn {tk : tkind}
g0float_copysign : (g0float tk, g0float tk) -<> g0float tk
overload copysign with g0float_copysign
implement
g0float_copysign<fltknd> (x, y) =
  $extfcall (float, "copysignf", x, y)
implement
g0float_copysign<dblknd> (x, y) =
  $extfcall (double, "copysign", x, y)
implement
g0float_copysign<ldblknd> (x, y) =
  $extfcall (ldouble, "copysignl", x, y)

(*------------------------------------------------------------------*)

typedef Matrix_Index_Map (m1 : int, n1 : int, m0 : int, n0 : int) =
  {i1, j1 : pos | i1 <= m1; j1 <= n1}
  (int i1, int j1) -<cloref0>
  [i0, j0 : pos | i0 <= m0; j0 <= n0]
    @(int i0, int j0)

datatype Real_Matrix (tk : tkind,
                      m1 : int, n1 : int,
                      m0 : int, n0 : int) =
| Real_Matrix of (matrixref (g0float tk, m0, n0),
                  int m1, int n1, int m0, int n0,
                  Matrix_Index_Map (m1, n1, m0, n0))
typedef Real_Matrix (tk : tkind, m1 : int, n1 : int) =
  [m0, n0 : pos] Real_Matrix (tk, m1, n1, m0, n0)
typedef Real_Vector (tk : tkind, m1 : int, n1 : int) =
  [m1 == 1 || n1 == 1] Real_Matrix (tk, m1, n1)
typedef Real_Row (tk : tkind, n1 : int) = Real_Vector (tk, 1, n1)
typedef Real_Column (tk : tkind, m1 : int) = Real_Vector (tk, m1, 1)

extern fn {tk : tkind}
Real_Matrix_make_elt :
  {m0, n0 : pos}
  (int m0, int n0, g0float tk) -< !wrt >
    Real_Matrix (tk, m0, n0, m0, n0)

extern fn {tk : tkind}
Real_Matrix_copy :
  {m1, n1 : pos}
  Real_Matrix (tk, m1, n1) -< !refwrt > Real_Matrix (tk, m1, n1)

extern fn {tk : tkind}
Real_Matrix_copy_to :
  {m1, n1 : pos}
  (Real_Matrix (tk, m1, n1),    (* destination *)
   Real_Matrix (tk, m1, n1)) -< !refwrt >
    void

extern fn {tk : tkind}
Real_Matrix_fill_with_elt :
  {m1, n1 : pos}
  (Real_Matrix (tk, m1, n1), g0float tk) -< !refwrt > void

extern fn {}
Real_Matrix_dimension :
  {tk : tkind}
  {m1, n1 : pos}
  Real_Matrix (tk, m1, n1) -<> @(int m1, int n1)

extern fn {tk : tkind}
Real_Matrix_get_at :
  {m1, n1 : pos}
  {i1, j1 : pos | i1 <= m1; j1 <= n1}
  (Real_Matrix (tk, m1, n1), int i1, int j1) -< !ref > g0float tk

extern fn {tk : tkind}
Real_Matrix_set_at :
  {m1, n1 : pos}
  {i1, j1 : pos | i1 <= m1; j1 <= n1}
  (Real_Matrix (tk, m1, n1), int i1, int j1, g0float tk) -< !refwrt >
    void

extern fn {}
Real_Matrix_transpose :
  (* This is transposed INDEXING. It does NOT copy the data. *)
  {tk : tkind}
  {m1, n1 : pos}
  {m0, n0 : pos}
  Real_Matrix (tk, m1, n1, m0, n0) -<>
    Real_Matrix (tk, n1, m1, m0, n0)

extern fn {}
Real_Matrix_block :
  (* This is block (submatrix) INDEXING. It does NOT copy the data. *)
  {tk : tkind}
  {p0, p1 : pos | p0 <= p1}
  {q0, q1 : pos | q0 <= q1}
  {m1, n1 : pos | p1 <= m1; q1 <= n1}
  {m0, n0 : pos}
  (Real_Matrix (tk, m1, n1, m0, n0),
   int p0, int p1, int q0, int q1) -<>
    Real_Matrix (tk, p1 - p0 + 1, q1 - q0 + 1, m0, n0)

extern fn {tk : tkind}
Real_Matrix_unit_matrix :
  {m : pos}
  int m -< !refwrt > Real_Matrix (tk, m, m)

extern fn {tk : tkind}
Real_Matrix_unit_matrix_to :
  {m : pos}
  Real_Matrix (tk, m, m) -< !refwrt > void

extern fn {tk : tkind}
Real_Matrix_matrix_sum :
  {m, n : pos}
  (Real_Matrix (tk, m, n), Real_Matrix (tk, m, n)) -< !refwrt >
    Real_Matrix (tk, m, n)

extern fn {tk : tkind}
Real_Matrix_matrix_sum_to :
  {m, n : pos}
  (Real_Matrix (tk, m, n),      (* destination*)
   Real_Matrix (tk, m, n),
   Real_Matrix (tk, m, n)) -< !refwrt >
    void

extern fn {tk : tkind}
Real_Matrix_matrix_difference :
  {m, n : pos}
  (Real_Matrix (tk, m, n), Real_Matrix (tk, m, n)) -< !refwrt >
    Real_Matrix (tk, m, n)

extern fn {tk : tkind}
Real_Matrix_matrix_difference_to :
  {m, n : pos}
  (Real_Matrix (tk, m, n),      (* destination*)
   Real_Matrix (tk, m, n),
   Real_Matrix (tk, m, n)) -< !refwrt >
    void

extern fn {tk : tkind}
Real_Matrix_matrix_product :
  {m, n, p : pos}
  (Real_Matrix (tk, m, n), Real_Matrix (tk, n, p)) -< !refwrt >
    Real_Matrix (tk, m, p)

extern fn {tk : tkind}
Real_Matrix_matrix_product_to :
  {m, n, p : pos}
  (Real_Matrix (tk, m, p),      (* destination*)
   Real_Matrix (tk, m, n),
   Real_Matrix (tk, n, p)) -< !refwrt >
    void

extern fn {tk : tkind}
Real_Matrix_scalar_product :
  {m, n : pos}
  (Real_Matrix (tk, m, n), g0float tk) -< !refwrt >
    Real_Matrix (tk, m, n)

extern fn {tk : tkind}
Real_Matrix_scalar_product_2 :
  {m, n : pos}
  (g0float tk, Real_Matrix (tk, m, n)) -< !refwrt >
    Real_Matrix (tk, m, n)

extern fn {tk : tkind}
Real_Matrix_scalar_product_to :
  {m, n : pos}
  (Real_Matrix (tk, m, n),      (* destination*)
   Real_Matrix (tk, m, n), g0float tk) -< !refwrt > void

extern fn {tk : tkind}
Real_Vector_l2norm_squared :
  {m, n : pos}
  Real_Vector (tk, m, n) -< !ref > g0float tk

extern fn {tk : tkind}
Real_Matrix_QR_decomposition :
  {m, n : pos}
  Real_Matrix (tk, m, n) -< !refwrt >
    @(Real_Matrix (tk, m, m), Real_Matrix (tk, m, n))

extern fn {tk : tkind}
Real_Matrix_least_squares_solution :
  (* This can solve p problems at once. Use p=1 to solve just Ax=b. *)
  {m, n, p : pos | n <= m}
  (Real_Matrix (tk, m, n), Real_Matrix (tk, m, p)) -< !refwrt >
    Real_Matrix (tk, n, p)

extern fn {tk : tkind}
Real_Matrix_fprint :
  {m, n : pos}
  (FILEref, Real_Matrix (tk, m, n)) -<1> void

overload copy with Real_Matrix_copy
overload copy_to with Real_Matrix_copy_to
overload fill_with_elt with Real_Matrix_fill_with_elt
overload dimension with Real_Matrix_dimension
overload [] with Real_Matrix_get_at
overload [] with Real_Matrix_set_at
overload transpose with Real_Matrix_transpose
overload block with Real_Matrix_block
overload unit_matrix with Real_Matrix_unit_matrix
overload unit_matrix_to with Real_Matrix_unit_matrix_to
overload matrix_sum with Real_Matrix_matrix_sum
overload matrix_sum_to with Real_Matrix_matrix_sum_to
overload matrix_difference with Real_Matrix_matrix_difference
overload matrix_difference_to with Real_Matrix_matrix_difference_to
overload matrix_product with Real_Matrix_matrix_product
overload matrix_product_to with Real_Matrix_matrix_product_to
overload scalar_product with Real_Matrix_scalar_product
overload scalar_product with Real_Matrix_scalar_product_2
overload scalar_product_to with Real_Matrix_scalar_product_to
overload + with matrix_sum
overload - with matrix_difference
overload * with matrix_product
overload * with scalar_product

(* Overload for a Real_Matrix_l2norm_squared, if we decided to have
   one, would be given precedence 0. *)
overload l2norm_squared with Real_Vector_l2norm_squared of 1

overload QR_decomposition with Real_Matrix_QR_decomposition
overload least_squares_solution with
  Real_Matrix_least_squares_solution

(*------------------------------------------------------------------*)

implement {tk}
Real_Matrix_make_elt (m0, n0, elt) =
  Real_Matrix (matrixref_make_elt<g0float tk> (i2sz m0, i2sz n0, elt),
               m0, n0, m0, n0, lam (i1, j1) => @(i1, j1))

implement {}
Real_Matrix_dimension A =
  case+ A of Real_Matrix (_, m1, n1, _, _, _) => @(m1, n1)

implement {tk}
Real_Matrix_get_at (A, i1, j1) =
  let
    val+ Real_Matrix (storage, _, _, _, n0, index_map) = A
    val @(i0, j0) = index_map (i1, j1)
  in
    matrixref_get_at<g0float tk> (storage, pred i0, n0, pred j0)
  end

implement {tk}
Real_Matrix_set_at (A, i1, j1, x) =
  let
    val+ Real_Matrix (storage, _, _, _, n0, index_map) = A
    val @(i0, j0) = index_map (i1, j1)
  in
    matrixref_set_at<g0float tk> (storage, pred i0, n0, pred j0, x)
  end

implement {}
Real_Matrix_transpose A =
  let
    val+ Real_Matrix (storage, m1, n1, m0, n0, index_map) = A
  in
    Real_Matrix (storage, n1, m1, m0, n0,
                 lam (i1, j1) => index_map (j1, i1))
  end

implement {}
Real_Matrix_block (A, p0, p1, q0, q1) =
  let
    val+ Real_Matrix (storage, m1, n1, m0, n0, index_map) = A
  in
    Real_Matrix (storage, succ (p1 - p0), succ (q1 - q0), m0, n0,
                 lam (i1, j1) =>
                  index_map (p0 + pred i1, q0 + pred j1))
  end

implement {tk}
Real_Matrix_copy A =
  let
    val @(m1, n1) = dimension A
    val C = Real_Matrix_make_elt<tk> (m1, n1, A[1, 1])
    val () = copy_to<tk> (C, A)
  in
    C
  end

implement {tk}
Real_Matrix_copy_to (Dst, Src) =
  let
    val @(m1, n1) = dimension Src
    prval [m1 : int] EQINT () = eqint_make_gint m1
    prval [n1 : int] EQINT () = eqint_make_gint n1

    var i : intGte 1
  in
    for* {i : pos | i <= m1 + 1} .<(m1 + 1) - i>.
         (i : int i) =>
      (i := 1; i <> succ m1; i := succ i)
        let
          var j : intGte 1
        in
          for* {j : pos | j <= n1 + 1} .<(n1 + 1) - j>.
               (j : int j) =>
            (j := 1; j <> succ n1; j := succ j)
              Dst[i, j] := Src[i, j]
        end
  end

implement {tk}
Real_Matrix_fill_with_elt (A, elt) =
  let
    val @(m1, n1) = dimension A
    prval [m1 : int] EQINT () = eqint_make_gint m1
    prval [n1 : int] EQINT () = eqint_make_gint n1

    var i : intGte 1
  in
    for* {i : pos | i <= m1 + 1} .<(m1 + 1) - i>.
         (i : int i) =>
      (i := 1; i <> succ m1; i := succ i)
        let
          var j : intGte 1
        in
          for* {j : pos | j <= n1 + 1} .<(n1 + 1) - j>.
               (j : int j) =>
            (j := 1; j <> succ n1; j := succ j)
              A[i, j] := elt
        end
  end

implement {tk}
Real_Matrix_unit_matrix {m} m =
  let
    val A = Real_Matrix_make_elt<tk> (m, m, Zero)
    var i : intGte 1
  in
    for* {i : pos | i <= m + 1} .<(m + 1) - i>.
         (i : int i) =>
      (i := 1; i <> succ m; i := succ i)
        A[i, i] := One;
    A
  end

implement {tk}
Real_Matrix_unit_matrix_to A =
  let
    val @(m, _) = dimension A
    prval [m : int] EQINT () = eqint_make_gint m

    var i : intGte 1
  in
    for* {i : pos | i <= m + 1} .<(m + 1) - i>.
         (i : int i) =>
      (i := 1; i <> succ m; i := succ i)
        let
          var j : intGte 1
        in
          for* {j : pos | j <= m + 1} .<(m + 1) - j>.
               (j : int j) =>
               (j := 1; j <> succ m; j := succ j)
            A[i, j] := (if i = j then One else Zero)
        end
  end

implement {tk}
Real_Matrix_matrix_sum (A, B) =
  let
    val @(m, n) = dimension A
    val C = Real_Matrix_make_elt<tk> (m, n, NAN)
    val () = matrix_sum_to<tk> (C, A, B)
  in
    C
  end

implement {tk}
Real_Matrix_matrix_sum_to (C, A, B) =
  let
    val @(m, n) = dimension A
    prval [m : int] EQINT () = eqint_make_gint m
    prval [n : int] EQINT () = eqint_make_gint n

    var i : intGte 1
  in
    for* {i : pos | i <= m + 1} .<(m + 1) - i>.
         (i : int i) =>
      (i := 1; i <> succ m; i := succ i)
        let
          var j : intGte 1
        in
          for* {j : pos | j <= n + 1} .<(n + 1) - j>.
               (j : int j) =>
            (j := 1; j <> succ n; j := succ j)
              C[i, j] := A[i, j] + B[i, j]
        end
  end

implement {tk}
Real_Matrix_matrix_difference (A, B) =
  let
    val @(m, n) = dimension A
    val C = Real_Matrix_make_elt<tk> (m, n, NAN)
    val () = matrix_difference_to<tk> (C, A, B)
  in
    C
  end

implement {tk}
Real_Matrix_matrix_difference_to (C, A, B) =
  let
    val @(m, n) = dimension A
    prval [m : int] EQINT () = eqint_make_gint m
    prval [n : int] EQINT () = eqint_make_gint n

    var i : intGte 1
  in
    for* {i : pos | i <= m + 1} .<(m + 1) - i>.
         (i : int i) =>
      (i := 1; i <> succ m; i := succ i)
        let
          var j : intGte 1
        in
          for* {j : pos | j <= n + 1} .<(n + 1) - j>.
               (j : int j) =>
            (j := 1; j <> succ n; j := succ j)
              C[i, j] := A[i, j] - B[i, j]
        end
  end

implement {tk}
Real_Matrix_matrix_product (A, B) =
  let
    val @(m, n) = dimension A and @(_, p) = dimension B
    val C = Real_Matrix_make_elt<tk> (m, p, NAN)
    val () = matrix_product_to<tk> (C, A, B)
  in
    C
  end

implement {tk}
Real_Matrix_matrix_product_to (C, A, B) =
  let
    val @(m, n) = dimension A and @(_, p) = dimension B
    prval [m : int] EQINT () = eqint_make_gint m
    prval [n : int] EQINT () = eqint_make_gint n
    prval [p : int] EQINT () = eqint_make_gint p

    var i : intGte 1
  in
    for* {i : pos | i <= m + 1} .<(m + 1) - i>.
         (i : int i) =>
      (i := 1; i <> succ m; i := succ i)
        let
          var k : intGte 1
        in
          for* {k : pos | k <= p + 1} .<(p + 1) - k>.
               (k : int k) =>
            (k := 1; k <> succ p; k := succ k)
              let
                var j : intGte 1
              in
                C[i, k] := A[i, 1] * B[1, k];
                for* {j : pos | j <= n + 1} .<(n + 1) - j>.
                     (j : int j) =>
                  (j := 2; j <> succ n; j := succ j)
                    C[i, k] :=
                      C[i, k] + (A[i, j] * B[j, k])
              end
        end
  end

implement {tk}
Real_Matrix_scalar_product (A, r) =
  let
    val @(m, n) = dimension A
    val C = Real_Matrix_make_elt<tk> (m, n, NAN)
    val () = scalar_product_to<tk> (C, A, r)
  in
    C
  end

implement {tk}
Real_Matrix_scalar_product_2 (r, A) =
  Real_Matrix_scalar_product<tk> (A, r)

implement {tk}
Real_Matrix_scalar_product_to (C, A, r) =
  let
    val @(m, n) = dimension A
    prval [m : int] EQINT () = eqint_make_gint m
    prval [n : int] EQINT () = eqint_make_gint n

    var i : intGte 1
  in
    for* {i : pos | i <= m + 1} .<(m + 1) - i>.
         (i : int i) =>
      (i := 1; i <> succ m; i := succ i)
        let
          var j : intGte 1
        in
          for* {j : pos | j <= n + 1} .<(n + 1) - j>.
               (j : int j) =>
            (j := 1; j <> succ n; j := succ j)
              C[i, j] := A[i, j] * r
        end
  end

implement {tk}
Real_Vector_l2norm_squared v =
  $effmask_wrt
  let
    val @(m, n) = dimension v
    prval [m : int] EQINT () = eqint_make_gint m
    prval [n : int] EQINT () = eqint_make_gint n
  in
    if n = 1 then
      let
        var sum : g0float tk
        var i : intGte 1
        val v11 = v[1, 1]
      in
        sum := v11 * v11;
        for* {i : pos | i <= m + 1} .<(m + 1) - i>.
             (i : int i) =>
          (i := 2; i <> succ m; i := succ i)
            let
              val vi1 = v[i, 1]
            in
              sum := sum + (vi1 * vi1)
            end;
        sum
      end
    else
      let
        var sum : g0float tk
        var j : intGte 1
        val v11 = v[1, 1]
      in
        sum := v11 * v11;
        for* {j : pos | j <= n + 1} .<(n + 1) - j>.
             (j : int j) =>
          (j := 2; j <> succ n; j := succ j)
            let
              val v1j = v[1, j]
            in
              sum := sum + (v1j * v1j)
            end;
        sum
      end
  end

implement {tk}
Real_Matrix_QR_decomposition A =
  (* Some of what follows does needless allocation and work, but
     making this code more efficient would be a project of its own!
     Also, one would likely want to implement pivot selection. See,
     for instance, Businger, P., Golub, G.H. Linear least squares
     solutions by householder transformations. Numer. Math. 7, 269–276
     (1965). https://doi.org/10.1007/BF01436084
     (https://web.archive.org/web/20230514003458/https://pages.stat.wisc.edu/~bwu62/771/businger1965.pdf)

     Note that I follow
     https://en.wikipedia.org/w/index.php?title=QR_decomposition&oldid=1152640697#Using_Householder_reflections
     more closely than I do what is stated in the task description at
     the time of this writing (13 May 2023). The presentation there
     seems simpler to me, and I prefer seeing a norm used to normalize
     the u vector. *)
  let
    val @(m, n) = dimension A
    prval [m : int] EQINT () = eqint_make_gint m
    prval [n : int] EQINT () = eqint_make_gint n

    stadef min_mn = min (m, n)
    val min_mn : int min_mn = min (m, n)

    var Q : Real_Matrix (tk, m, m) = unit_matrix<tk> m
    val R : Real_Matrix (tk, m, n) = copy A

    (* I_mm is a unit matrix of the maximum size used. Smaller unit
       matrices will be had by the "identity" function, and unit
       column vectors by the "unit_column" function. *)
    val I_mm : Real_Matrix (tk, m, m) = unit_matrix<tk> m
    fn
    identity {p : pos | p <= m}
             (p : int p) :<> Real_Matrix (tk, p, p) =
      block (I_mm, 1, p, 1, p)
    fn
    unit_column {p, j : pos | j <= p; p <= m}
                (p    : int p,
                 j    : int j) :<> Real_Column (tk, p) =
      block (I_mm, 1, p, j, j)

    var k : intGte 1
  in
    for* {k : pos | k <= min_mn} .<min_mn - k>.
         (k : int k) =>
      (k := 1; k <> min_mn; k := succ k)
        let
          val x = block (R, k, m, k, k)
          val sigma = l2norm_squared x

          (* Choose the sign of alpha to increase the magnitude of the
             pivot. *)
          val alpha = copysign (sqrt sigma, ~x[1, 1])

          val e1 = unit_column (succ (m - k), 1)
          val u = x - (alpha * e1)
          val v = u * (One / sqrt (l2norm_squared u))
          val I = identity (succ (m - k))
          val H = I - (Two * v * transpose v)

          (* Update R, using block operations. *)
          val () = fill_with_elt<tk> (x, Zero)
          val () = x[1, 1] := alpha
          val R_ = block (R, k, m, succ k, n)
          val Tmp = H * R_
          val () = copy_to (R_, Tmp)

          (* Update Q. *)
          val Tmp = unit_matrix m
          val Tmp_ = block (Tmp, k, m, k, m)
          val () = copy_to (Tmp_, H)
          val () = Q := Q * Tmp
        in
        end;
    @(Q, R)
  end

implement {tk}
Real_Matrix_least_squares_solution (A, B) =
  let
    (* I use this algorithm for the back substitutions:
    https://algowiki-project.org/algowiki/en/index.php?title=Backward_substitution&oldid=10412#Approaches_and_features_of_implementing_the_back_substitution_algorithm_in_parallel
    *)

    val @(m, n) = dimension A and @(_, p) = dimension B
    prval [m : int] EQINT () = eqint_make_gint m
    prval [n : int] EQINT () = eqint_make_gint n
    prval [p : int] EQINT () = eqint_make_gint p

    val @(Q, R) = QR_decomposition<tk> A

    (* X is initialized for back substitutions. *)
    val X = block (transpose Q * B, 1, n, 1, p)
    and R = block (R, 1, n, 1, n)

    var k : intGte 1
  in
    (* Complete the back substitutions. *)
    for* {k : pos | k <= p + 1} .<(p + 1) - k>.
         (k : int k) =>
      (k := 1; k <> succ p; k := succ k)
        let
          val x = block (X, 1, n, k, k)
          var j : intGte 0
        in
          for* {j : nat | 0 <= j; j <= n} .<j>.
               (j : int j) =>
            (j := n; j <> 0; j := pred j)
              let
                var i : intGte 1
              in
                x[j, 1] := x[j, 1] / R[j, j];
                for* {i : pos | i <= j} .<j - i>.
                     (i : int i) =>
                  (i := 1; i <> j; i := succ i)
                    x[i, 1] := x[i, 1] - (R[i, j] * x[j, 1])
              end
        end;
    X
  end

implement {tk}
Real_Matrix_fprint {m, n} (outf, A) =
  let
    val @(m, n) = dimension A
    var i : intGte 1
  in
    for* {i : pos | i <= m + 1} .<(m + 1) - i>.
         (i : int i) =>
      (i := 1; i <> succ m; i := succ i)
        let
          var j : intGte 1
        in
          for* {j : pos | j <= n + 1} .<(n + 1) - j>.
               (j : int j) =>
            (j := 1; j <> succ n; j := succ j)
              let
                typedef FILEstar = $extype"FILE *"
                extern castfn FILEref2star : FILEref -<> FILEstar
                val _ = $extfcall (int, "fprintf", FILEref2star outf,
                                   "%16.6g", A[i, j])
              in
              end;
          fprintln! (outf)
        end
  end

(*------------------------------------------------------------------*)

implement
main0 () =
  let
    stadef fltknd = dblknd
    macdef i2flt = g0int2float<intknd,dblknd>

    val A = Real_Matrix_make_elt<fltknd> (3, 3, NAN)
    val () =
      begin
        A[1, 1] := i2flt 12;
        A[2, 1] := i2flt 6;
        A[3, 1] := i2flt ~4;

        A[1, 2] := i2flt ~51;
        A[2, 2] := i2flt 167;
        A[3, 2] := i2flt 24;

        A[1, 3] := i2flt 4;
        A[2, 3] := i2flt ~68;
        A[3, 3] := i2flt ~41
      end

    val @(Q, R) = QR_decomposition<fltknd> A

    (* Example of least-squares solution. (Copied from the BBC BASIC
       or Common Lisp entry, whichever you prefer to think it copied
       from.) *)
    val x = $list (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
    and y = $list (1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321)
    val X = Real_Matrix_make_elt<fltknd> (11, 3, NAN)
    and Y = Real_Matrix_make_elt<fltknd> (11, 1, NAN)
    val () =
      let
        var i : intGte 1
      in
        for* {i : pos | i <= 12} .<12 - i>.
             (i : int i) =>
          (i := 1; i <> 12; i := succ i)
            let
              val xi = x[pred i] : int
              and yi = y[pred i] : int
            in
              X[i, 1] := g0i2f (xi ** 0);
              X[i, 2] := g0i2f (xi ** 1);
              X[i, 3] := g0i2f (xi ** 2);
              Y[i, 1] := g0i2f yi
            end
      end
    val solution = least_squares_solution (X, Y)
  in
    println! ("A :");
    Real_Matrix_fprint (stdout_ref, A);
    println! ();
    println! ("Q :");
    Real_Matrix_fprint (stdout_ref, Q);
    println! ();
    println! ("R :");
    Real_Matrix_fprint (stdout_ref, R);
    println! ();
    println! ("Q * R :");
    Real_Matrix_fprint (stdout_ref, Q * R);
    println! ();
    println! ("least squares A in Ax=b :");
    Real_Matrix_fprint (stdout_ref, X);
    println! ();
    println! ("least squares b in Ax=b :");
    Real_Matrix_fprint (stdout_ref, Y);
    println! ();
    println! ("least squares solution :");
    Real_Matrix_fprint (stdout_ref, solution)
  end

(*------------------------------------------------------------------*)

Output:

$ patscc -std=gnu2x -g -O2 -DATS_MEMALLOC_GCBDW qr_decomposition_task.dats -lgc -lm && ./a.out
A :
              12             -51               4
               6             167             -68
              -4              24             -41

Q :
       -0.857143        0.394286        0.331429
       -0.428571       -0.902857      -0.0342857
        0.285714       -0.171429        0.942857

R :
             -14             -21              14
               0            -175              70
               0               0             -35

Q * R :
              12             -51               4
               6             167             -68
              -4              24             -41

least squares A in Ax=b :
               1               0               0
               1               1               1
               1               2               4
               1               3               9
               1               4              16
               1               5              25
               1               6              36
               1               7              49
               1               8              64
               1               9              81
               1              10             100

least squares b in Ax=b :
               1
               6
              17
              34
              57
              86
             121
             162
             209
             262
             321

least squares solution :
               1
               2
               3

Axiom

The following provides a generic QR decomposition for arbitrary precision floats, double floats and exact calculations:

)abbrev package TESTP TestPackage
TestPackage(R:Join(Field,RadicalCategory)): with
    unitVector: NonNegativeInteger -> Vector(R)
    "/": (Vector(R),R) -> Vector(R)
    "^": (Vector(R),NonNegativeInteger) -> Vector(R)
    solveUpperTriangular: (Matrix(R),Vector(R)) -> Vector(R)
    signValue: R -> R
    householder: Vector(R) -> Matrix(R)
    qr: Matrix(R) -> Record(q:Matrix(R),r:Matrix(R))
    lsqr: (Matrix(R),Vector(R)) -> Vector(R)
    polyfit: (Vector(R),Vector(R),NonNegativeInteger) -> Vector(R)
  == add
    unitVector(dim) ==
      out := new(dim,0@R)$Vector(R)
      out(1) := 1@R
      out
    v:Vector(R) / a:R == map((vi:R):R +-> vi/a, v)$Vector(R)
    v:Vector(R) ^ n:NonNegativeInteger == map((vi:R):R +-> vi^n, v)$Vector(R)
    solveUpperTriangular(r,b) ==
      n := ncols r
      x := new(n,0@R)$Vector(R)
      for k in n..1 by -1 repeat
        index := min(n,k+1)
	x(k) := (b(k)-reduce("+",subMatrix(r,k,k,index,n)*x.(index..n)))/r(k,k)
      x
    signValue(r) ==
      R has (sign: R -> Integer) => coerce(sign(r)$R)$R
      zero? r => r
      if sqrt(r*r) = r then 1 else -1
    householder(a) ==
      m := #a
      u := a + length(a)*signValue(a(1))*unitVector(m) 
      v := u/u(1) 
      beta := (1+1)/dot(v,v)
      scalarMatrix(m,1) - beta*transpose(outerProduct(v,v))
    qr(a) ==
      (m,n) := (nrows a, ncols a)
      qm := scalarMatrix(m,1)
      rm := copy a
      for i in 1..(if m=n then n-1 else n) repeat
        x := column(subMatrix(rm,i,m,i,i),1)
	h := scalarMatrix(m,1)
	setsubMatrix!(h,i,i,householder x)
	qm := qm*h
	rm := h*rm
      [qm,rm]
    lsqr(a,b) ==
      dc := qr a
      n := ncols(dc.r)
      solveUpperTriangular(subMatrix(dc.r,1,n,1,n),transpose(dc.q)*b)
    polyfit(x,y,n) ==
      a := new(#x,n+1,0@R)$Matrix(R)
      for j in 0..n repeat
        setColumn!(a,j+1,x^j)
      lsqr(a,y)

This can be called using:

m := matrix [[12, -51, 4], [6, 167, -68], [-4, 24, -41]];
qr m
x := vector [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10];
y := vector [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321];
polyfit(x, y, 2)

With output in exact form:

qr m

             +  6    69     58 +
             |- -   ---    --- |
             |  7   175    175 |
             |                 |    +- 14  - 21    14 +
             |  3    158     6 |    |                 |
         [q= |- -  - ---  - ---|,r= | 0    - 175   70 |]
             |  7    175    175|    |                 |
             |                 |    + 0      0    - 35+
             | 2      6    33  |
             | -   - --    --  |
             + 7     35    35  +

          Type: Record(q: Matrix(AlgebraicNumber),r: Matrix(AlgebraicNumber))

polyfit(x, y, 2)

   [1,2,3]
                              Type: Vector(AlgebraicNumber)

The calculations are comparable to those from the default QR decomposition in R.

BBC BASIC

Works with: BBC BASIC for Windows

Makes heavy use of BBC BASIC's matrix arithmetic.

      *FLOAT 64
      @% = &2040A
      INSTALL @lib$+"ARRAYLIB"
      
      REM Test matrix for QR decomposition:
      DIM A(2,2)
      A() = 12, -51,   4, \
      \      6, 167, -68, \
      \     -4,  24, -41
      
      REM Do the QR decomposition:
      DIM Q(2,2), R(2,2)
      PROCqrdecompose(A(), Q(), R())
      PRINT "Q:"
      PRINT Q(0,0), Q(0,1), Q(0,2)
      PRINT Q(1,0), Q(1,1), Q(1,2)
      PRINT Q(2,0), Q(2,1), Q(2,2)
      PRINT "R:"
      PRINT R(0,0), R(0,1), R(0,2)
      PRINT R(1,0), R(1,1), R(1,2)
      PRINT R(2,0), R(2,1), R(2,2)
      
      REM Test data for least-squares solution:
      DIM x(10) : x() = 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
      DIM y(10) : y() = 1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321
      
      REM Do the least-squares solution:
      DIM a(10,2), q(10,10), r(10,2), t(10,10), b(10), z(2)
      FOR i% = 0 TO 10
        FOR j% = 0 TO 2
          a(i%,j%) = x(i%) ^ j%
        NEXT
      NEXT
      PROCqrdecompose(a(), q(), r())
      PROC_transpose(q(),t())
      b() = t() . y()
      FOR k% = 2 TO 0 STEP -1
        s = 0
        IF k% < 2 THEN
          FOR j% = k%+1 TO 2
            s += r(k%,j%) * z(j%)
          NEXT
        ENDIF
        z(k%) = (b(k%) - s) / r(k%,k%)
      NEXT k%
      PRINT '"Least-squares solution:"
      PRINT z(0), z(1), z(2)
      END
      
      DEF PROCqrdecompose(A(), Q(), R())
      LOCAL i%, k%, m%, n%, H()
      m% = DIM(A(),1) : n% = DIM(A(),2)
      DIM H(m%,m%)
      FOR i% = 0 TO m% : Q(i%,i%) = 1 : NEXT
      WHILE n%
        PROCqrstep(n%, k%, A(), H())
        A() = H() . A()
        Q() = Q() . H()
        k% += 1
        m% -= 1
        n% -= 1
      ENDWHILE
      R() = A()
      ENDPROC
      
      DEF PROCqrstep(n%, k%, A(), H())
      LOCAL a(), h(), i%, j%
      DIM a(n%,0), h(n%,n%)
      FOR i% = 0 TO n% : a(i%,0) = A(i%+k%,k%) : NEXT
      PROChouseholder(h(), a())
      H() = 0  : H(0,0) = 1
      FOR i% = 0 TO n%
        FOR j% = 0 TO n%
          H(i%+k%,j%+k%) = h(i%,j%)
        NEXT
      NEXT
      ENDPROC
      
      REM Create the Householder matrix for the supplied column vector:
      DEF PROChouseholder(H(), a())
      LOCAL e(), u(), v(), vt(), vvt(), I(), d()
      LOCAL i%, n% : n% = DIM(a(),1)
      REM Create the scaled standard basis vector e():
      DIM e(n%,0) : e(0,0) = SGN(a(0,0)) * MOD(a())
      REM Create the normal vector u():
      DIM u(n%,0) : u() = a() + e()
      REM Normalise with respect to the first element:
      DIM v(n%,0) : v() = u() / u(0,0)
      REM Get the transpose of v() and its dot product with v():
      DIM vt(0,n%), d(0) : PROC_transpose(v(), vt()) : d() = vt() . v()
      REM Get the product of v() and vt():
      DIM vvt(n%,n%) : vvt() = v() . vt()
      REM Create an identity matrix I():
      DIM I(n%,n%) : FOR i% = 0 TO n% : I(i%,i%) = 1 : NEXT
      REM Create the Householder matrix H() = I - 2/vt()v() v()vt():
      vvt() *= 2 / d(0) : H() = I() - vvt()
      ENDPROC

Output:

Q:
   -0.8571    0.3943    0.3314
   -0.4286   -0.9029   -0.0343
    0.2857   -0.1714    0.9429
R:
  -14.0000  -21.0000   14.0000
    0.0000 -175.0000   70.0000
    0.0000    0.0000  -35.0000

Least-squares solution:
    1.0000    2.0000    3.0000

C

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

typedef struct {
	int m, n;
	double ** v;
} mat_t, *mat;

mat matrix_new(int m, int n)
{
	mat x = malloc(sizeof(mat_t));
	x->v = malloc(sizeof(double*) * m);
	x->v[0] = calloc(sizeof(double), m * n);
	for (int i = 0; i < m; i++)
		x->v[i] = x->v[0] + n * i;
	x->m = m;
	x->n = n;
	return x;
}

void matrix_delete(mat m)
{
	free(m->v[0]);
	free(m->v);
	free(m);
}

void matrix_transpose(mat m)
{
	for (int i = 0; i < m->m; i++) {
		for (int j = 0; j < i; j++) {
			double t = m->v[i][j];
			m->v[i][j] = m->v[j][i];
			m->v[j][i] = t;
		}
	}
}

mat matrix_copy(int n, double a[][n], int m)
{
	mat x = matrix_new(m, n);
	for (int i = 0; i < m; i++)
		for (int j = 0; j < n; j++)
			x->v[i][j] = a[i][j];
	return x;
}

mat matrix_mul(mat x, mat y)
{
	if (x->n != y->m) return 0;
	mat r = matrix_new(x->m, y->n);
	for (int i = 0; i < x->m; i++)
		for (int j = 0; j < y->n; j++)
			for (int k = 0; k < x->n; k++)
				r->v[i][j] += x->v[i][k] * y->v[k][j];
	return r;
}

mat matrix_minor(mat x, int d)
{
	mat m = matrix_new(x->m, x->n);
	for (int i = 0; i < d; i++)
		m->v[i][i] = 1;
	for (int i = d; i < x->m; i++)
		for (int j = d; j < x->n; j++)
			m->v[i][j] = x->v[i][j];
	return m;
}

/* c = a + b * s */
double *vmadd(double a[], double b[], double s, double c[], int n)
{
	for (int i = 0; i < n; i++)
		c[i] = a[i] + s * b[i];
	return c;
}

/* m = I - v v^T */
mat vmul(double v[], int n)
{
	mat x = matrix_new(n, n);
	for (int i = 0; i < n; i++)
		for (int j = 0; j < n; j++)
			x->v[i][j] = -2 *  v[i] * v[j];
	for (int i = 0; i < n; i++)
		x->v[i][i] += 1;

	return x;
}

/* ||x|| */
double vnorm(double x[], int n)
{
	double sum = 0;
	for (int i = 0; i < n; i++) sum += x[i] * x[i];
	return sqrt(sum);
}

/* y = x / d */
double* vdiv(double x[], double d, double y[], int n)
{
	for (int i = 0; i < n; i++) y[i] = x[i] / d;
	return y;
}

/* take c-th column of m, put in v */
double* mcol(mat m, double *v, int c)
{
	for (int i = 0; i < m->m; i++)
		v[i] = m->v[i][c];
	return v;
}

void matrix_show(mat m)
{
	for(int i = 0; i < m->m; i++) {
		for (int j = 0; j < m->n; j++) {
			printf(" %8.3f", m->v[i][j]);
		}
		printf("\n");
	}
	printf("\n");
}

void householder(mat m, mat *R, mat *Q)
{
	mat q[m->m];
	mat z = m, z1;
	for (int k = 0; k < m->n && k < m->m - 1; k++) {
		double e[m->m], x[m->m], a;
		z1 = matrix_minor(z, k);
		if (z != m) matrix_delete(z);
		z = z1;

		mcol(z, x, k);
		a = vnorm(x, m->m);
		if (m->v[k][k] > 0) a = -a;

		for (int i = 0; i < m->m; i++)
			e[i] = (i == k) ? 1 : 0;

		vmadd(x, e, a, e, m->m);
		vdiv(e, vnorm(e, m->m), e, m->m);
		q[k] = vmul(e, m->m);
		z1 = matrix_mul(q[k], z);
		if (z != m) matrix_delete(z);
		z = z1;
	}
	matrix_delete(z);
	*Q = q[0];
	*R = matrix_mul(q[0], m);
	for (int i = 1; i < m->n && i < m->m - 1; i++) {
		z1 = matrix_mul(q[i], *Q);
		if (i > 1) matrix_delete(*Q);
		*Q = z1;
		matrix_delete(q[i]);
	}
	matrix_delete(q[0]);
	z = matrix_mul(*Q, m);
	matrix_delete(*R);
	*R = z;
	matrix_transpose(*Q);
}

double in[][3] = {
	{ 12, -51,   4},
	{  6, 167, -68},
	{ -4,  24, -41},
	{ -1, 1, 0},
	{ 2, 0, 3},
};

int main()
{
	mat R, Q;
	mat x = matrix_copy(3, in, 5);
	householder(x, &R, &Q);

	puts("Q"); matrix_show(Q);
	puts("R"); matrix_show(R);

	// to show their product is the input matrix
	mat m = matrix_mul(Q, R);
	puts("Q * R"); matrix_show(m);

	matrix_delete(x);
	matrix_delete(R);
	matrix_delete(Q);
	matrix_delete(m);
	return 0;
}

Output:

Q
    0.846   -0.391    0.343    0.082    0.078
    0.423    0.904   -0.029    0.026    0.045
   -0.282    0.170    0.933   -0.047   -0.137
   -0.071    0.014   -0.001    0.980   -0.184
    0.141   -0.017   -0.106   -0.171   -0.969

R
   14.177   20.667  -13.402
   -0.000  175.043  -70.080
    0.000    0.000  -35.202
   -0.000   -0.000   -0.000
    0.000    0.000   -0.000

Q * R
   12.000  -51.000    4.000
    6.000  167.000  -68.000
   -4.000   24.000  -41.000
   -1.000    1.000   -0.000
    2.000   -0.000    3.000

C#

Library: Math.Net

using System;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;


class Program
{

    static void Main(string[] args)
    {
        Matrix<double> A = DenseMatrix.OfArray(new double[,]
        {
                {  12,  -51,    4 },
                {   6,  167,  -68 },
                {  -4,   24,  -41 }
        });
        Console.WriteLine("A:");
        Console.WriteLine(A);
        var qr = A.QR();
        Console.WriteLine();
        Console.WriteLine("Q:");
        Console.WriteLine(qr.Q);
        Console.WriteLine();
        Console.WriteLine("R:");
        Console.WriteLine(qr.R);
    }
}

Output:

A:
DenseMatrix 3x3-Double
12  -51    4
 6  167  -68
-4   24  -41


Q:
DenseMatrix 3x3-Double
-0.857143   0.394286  -0.331429
-0.428571  -0.902857  0.0342857
 0.285714  -0.171429  -0.942857


R:
DenseMatrix 3x3-Double
-14   -21  14
  0  -175  70
  0     0  35

C++

/*
 * g++ -O3 -Wall --std=c++11 qr_standalone.cpp -o qr_standalone
 */
#include <cstdio>
#include <cstdlib>
#include <cstring> // for memset
#include <limits>
#include <iostream>
#include <vector>

#include <math.h>

class Vector;

class Matrix {

public:
  // default constructor (don't allocate)
  Matrix() : m(0), n(0), data(nullptr) {}
  
  // constructor with memory allocation, initialized to zero
  Matrix(int m_, int n_) : Matrix() {
    m = m_;
    n = n_;
    allocate(m_,n_);
  }

  // copy constructor
  Matrix(const Matrix& mat) : Matrix(mat.m,mat.n) {

    for (int i = 0; i < m; i++)
      for (int j = 0; j < n; j++)
	(*this)(i,j) = mat(i,j);
  }
  
  // constructor from array
  template<int rows, int cols>
  Matrix(double (&a)[rows][cols]) : Matrix(rows,cols) {

    for (int i = 0; i < m; i++)
      for (int j = 0; j < n; j++)
	(*this)(i,j) = a[i][j];
  }

  // destructor
  ~Matrix() {
    deallocate();
  }


  // access data operators
  double& operator() (int i, int j) {
    return data[i+m*j]; }
  double  operator() (int i, int j) const {
    return data[i+m*j]; }

  // operator assignment
  Matrix& operator=(const Matrix& source) {
    
    // self-assignment check
    if (this != &source) { 
      if ( (m*n) != (source.m * source.n) ) { // storage cannot be reused
	allocate(source.m,source.n);          // re-allocate storage
      }
      // storage can be used, copy data
      std::copy(source.data, source.data + source.m*source.n, data);
    }
    return *this;
  }
  
  // compute minor
  void compute_minor(const Matrix& mat, int d) {

    allocate(mat.m, mat.n);
    
    for (int i = 0; i < d; i++)
      (*this)(i,i) = 1.0;
    for (int i = d; i < mat.m; i++)
      for (int j = d; j < mat.n; j++)
	(*this)(i,j) = mat(i,j);
    
  }

  // Matrix multiplication
  // c = a * b
  // c will be re-allocated here
  void mult(const Matrix& a, const Matrix& b) {

    if (a.n != b.m) {
      std::cerr << "Matrix multiplication not possible, sizes don't match !\n";
      return;
    }

    // reallocate ourself if necessary i.e. current Matrix has not valid sizes
    if (a.m != m or b.n != n)
      allocate(a.m, b.n);

    memset(data,0,m*n*sizeof(double));
    
    for (int i = 0; i < a.m; i++)
      for (int j = 0; j < b.n; j++)
	for (int k = 0; k < a.n; k++)
	  (*this)(i,j) += a(i,k) * b(k,j);
    
  }

  void transpose() {
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < i; j++) {
	double t = (*this)(i,j);
	(*this)(i,j) = (*this)(j,i);
	(*this)(j,i) = t;
      }
    }
  }

  // take c-th column of m, put in v
  void extract_column(Vector& v, int c);  

  // memory allocation
  void allocate(int m_, int n_) {

    // if already allocated, memory is freed
    deallocate();
    
    // new sizes
    m = m_;
    n = n_;
    
    data = new double[m_*n_];
    memset(data,0,m_*n_*sizeof(double));

  } // allocate

  // memory free
  void deallocate() {

    if (data)
      delete[] data;

    data = nullptr;

  }    
  
  int m, n;
  
private:
  double* data;
  
}; // struct Matrix

// column vector
class Vector {

public:
  // default constructor (don't allocate)
  Vector() : size(0), data(nullptr) {}
  
  // constructor with memory allocation, initialized to zero
  Vector(int size_) : Vector() {
    size = size_;
    allocate(size_);
  }

  // destructor
  ~Vector() {
    deallocate();
  }

  // access data operators
  double& operator() (int i) {
    return data[i]; }
  double  operator() (int i) const {
    return data[i]; }

  // operator assignment
  Vector& operator=(const Vector& source) {
    
    // self-assignment check
    if (this != &source) { 
      if ( size != (source.size) ) {   // storage cannot be reused
	allocate(source.size);         // re-allocate storage
      }
      // storage can be used, copy data
      std::copy(source.data, source.data + source.size, data);
    }
    return *this;
  }

  // memory allocation
  void allocate(int size_) {

    deallocate();
    
    // new sizes
    size = size_;
    
    data = new double[size_];
    memset(data,0,size_*sizeof(double));

  } // allocate

  // memory free
  void deallocate() {

    if (data)
      delete[] data;

    data = nullptr;

  }    

  //   ||x||
  double norm() {
    double sum = 0;
    for (int i = 0; i < size; i++) sum += (*this)(i) * (*this)(i);
    return sqrt(sum);
  }

  // divide data by factor
  void rescale(double factor) {
    for (int i = 0; i < size; i++) (*this)(i) /= factor;
  }

  void rescale_unit() {
    double factor = norm();
    rescale(factor);
  }
  
  int size;
  
private:
  double* data;

}; // class Vector

// c = a + b * s
void vmadd(const Vector& a, const Vector& b, double s, Vector& c)
{
  if (c.size != a.size or c.size != b.size) {
    std::cerr << "[vmadd]: vector sizes don't match\n";
    return;
  }
  
  for (int i = 0; i < c.size; i++)
    c(i) = a(i) + s * b(i);
}

// mat = I - 2*v*v^T
// !!! m is allocated here !!!
void compute_householder_factor(Matrix& mat, const Vector& v)
{

  int n = v.size;
  mat.allocate(n,n);
  for (int i = 0; i < n; i++)
    for (int j = 0; j < n; j++)
      mat(i,j) = -2 *  v(i) * v(j);
  for (int i = 0; i < n; i++)
    mat(i,i) += 1;  
}

// take c-th column of a matrix, put results in Vector v
void Matrix::extract_column(Vector& v, int c) {
  if (m != v.size) {
    std::cerr << "[Matrix::extract_column]: Matrix and Vector sizes don't match\n";
    return;
  }
  
  for (int i = 0; i < m; i++)
    v(i) = (*this)(i,c);
}

void matrix_show(const Matrix&  m, const std::string& str="")
{
  std::cout << str << "\n";
  for(int i = 0; i < m.m; i++) {
    for (int j = 0; j < m.n; j++) {
      printf(" %8.3f", m(i,j));
    }
    printf("\n");
  }
  printf("\n");
}

// L2-norm ||A-B||^2
double matrix_compare(const Matrix& A, const Matrix& B) {
  // matrices must have same size
  if (A.m != B.m or  A.n != B.n)
    return std::numeric_limits<double>::max();

  double res=0;
  for(int i = 0; i < A.m; i++) {
    for (int j = 0; j < A.n; j++) {
      res += (A(i,j)-B(i,j)) * (A(i,j)-B(i,j));
    }
  }

  res /= A.m*A.n;
  return res;
}

void householder(Matrix& mat,
		 Matrix& R,
		 Matrix& Q)
{

  int m = mat.m;
  int n = mat.n;

  // array of factor Q1, Q2, ... Qm
  std::vector<Matrix> qv(m);

  // temp array
  Matrix z(mat);
  Matrix z1;
  
  for (int k = 0; k < n && k < m - 1; k++) {

    Vector e(m), x(m);
    double a;
    
    // compute minor
    z1.compute_minor(z, k);
    
    // extract k-th column into x
    z1.extract_column(x, k);
    
    a = x.norm();
    if (mat(k,k) > 0) a = -a;
    
    for (int i = 0; i < e.size; i++)
      e(i) = (i == k) ? 1 : 0;

    // e = x + a*e
    vmadd(x, e, a, e);

    // e = e / ||e||
    e.rescale_unit();
    
    // qv[k] = I - 2 *e*e^T
    compute_householder_factor(qv[k], e);

    // z = qv[k] * z1
    z.mult(qv[k], z1);

  }
  
  Q = qv[0];

  // after this loop, we will obtain Q (up to a transpose operation)
  for (int i = 1; i < n && i < m - 1; i++) {

    z1.mult(qv[i], Q);
    Q = z1;
    
  }
  
  R.mult(Q, mat);
  Q.transpose();
}
 
double in[][3] = {
  { 12, -51,   4},
  {  6, 167, -68},
  { -4,  24, -41},
  { -1,   1,   0},
  {  2,   0,   3},
};
 
int main()
{
  Matrix A(in); 
  Matrix Q, R;

  matrix_show(A,"A");  

  // compute QR decompostion
  householder(A, R, Q);
  
  matrix_show(Q,"Q");
  matrix_show(R,"R");
 
  // compare Q*R to the original matrix A
  Matrix A_check;
  A_check.mult(Q, R);

  // compute L2 norm ||A-A_check||^2
  double l2 = matrix_compare(A,A_check);

  // display Q*R
  matrix_show(A_check, l2 < 1e-12 ? "A == Q * R ? yes" : "A == Q * R ? no");
 
  return EXIT_SUCCESS;
}

Output:

A
   12.000  -51.000    4.000
    6.000  167.000  -68.000
   -4.000   24.000  -41.000
   -1.000    1.000    0.000
    2.000    0.000    3.000

Q
    0.846   -0.391    0.343    0.082    0.078
    0.423    0.904   -0.029    0.026    0.045
   -0.282    0.170    0.933   -0.047   -0.137
   -0.071    0.014   -0.001    0.980   -0.184
    0.141   -0.017   -0.106   -0.171   -0.969

R
   14.177   20.667  -13.402
   -0.000  175.043  -70.080
    0.000    0.000  -35.202
   -0.000   -0.000   -0.000
    0.000    0.000   -0.000

A == Q * R ? yes
   12.000  -51.000    4.000
    6.000  167.000  -68.000
   -4.000   24.000  -41.000
   -1.000    1.000   -0.000
    2.000   -0.000    3.000

With Polynomial Fitting

#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <string>
#include <vector>

class Matrix {
public:
	Matrix(const std::vector<std::vector<double>>& data) : data(data) {
		initialise();
	}

	Matrix(const Matrix& matrix) : data(matrix.data) {
		initialise();
	}

	Matrix(const uint64_t& row_count, const uint64_t& column_count) {
		data.assign(row_count, std::vector<double>(column_count, 0.0));
		initialise();
	}

	Matrix add(const Matrix& other) {
		if ( other.row_count != row_count || other.column_count != column_count ) {
			throw std::invalid_argument("Incompatible matrix dimensions.");
		}

		Matrix result(data);
		for ( int32_t i = 0; i < row_count; ++i ) {
			for ( int32_t j = 0; j < column_count; ++j ) {
				result.data[i][j] = data[i][j] + other.data[i][j];
			}
		}
		return result;
	}

	Matrix multiply(const Matrix& other) {
		if ( column_count != other.row_count ) {
			throw std::invalid_argument("Incompatible matrix dimensions.");
		}

		Matrix result(row_count, other.column_count);
		for ( int32_t i = 0; i < row_count; ++i ) {
			for ( int32_t j = 0; j < other.column_count; ++j ) {
				for ( int32_t k = 0; k < row_count; k++ ) {
					result.data[i][j] += data[i][k] * other.data[k][j];
				}
			}
		}
		return result;
	}

	Matrix transpose() {
		Matrix result(column_count, row_count);
		for ( int32_t i = 0; i < row_count; ++i ) {
			for ( int32_t j = 0; j < column_count; ++j ) {
				result.data[j][i] = data[i][j];
			}
		}
		return result;
	}

	Matrix minor(const int32_t& index) {
		Matrix result(row_count, column_count);
		for ( int32_t i = 0; i < index; ++i ) {
			result.set_entry(i, i, 1.0);
		}

		for ( int32_t i = index; i < row_count; ++i ) {
			for ( int32_t j = index; j < column_count; ++j ) {
				result.set_entry(i, j, data[i][j]);
			}
		}
		return result;
	}

	Matrix column(const int32_t& index) {
		Matrix result(row_count, 1);
		for ( int32_t i = 0; i < row_count; ++i ) {
			result.set_entry(i, 0, data[i][index]);
		}
		return result;
	}

	Matrix scalarMultiply(const double& value) {
		if ( column_count != 1 ) {
			throw std::invalid_argument("Incompatible matrix dimension.");
		}

		Matrix result(row_count, column_count);
		for ( int32_t i = 0; i < row_count; ++i ) {
			result.data[i][0] = data[i][0] * value;
		}
		return result;
	}

	Matrix unit() {
		if ( column_count != 1 ) {
			throw std::invalid_argument("Incompatible matrix dimensions.");
		}

		const double the_magnitude = magnitude();
		Matrix result(row_count, column_count);
		for ( int32_t i = 0; i < row_count; ++i ) {
			result.data[i][0] = data[i][0] / the_magnitude;
		}
		return result;
	}

	double magnitude() {
		if ( column_count != 1 ) {
			throw std::invalid_argument("Incompatible matrix dimensions.");
		}

		double norm = 0.0;
		for ( int32_t i = 0; i < row_count; ++i ) {
			norm += data[i][0] * data[i][0];
		}
		return std::sqrt(norm);
	}

	int32_t size() {
		if ( column_count != 1 ) {
			throw std::invalid_argument("Incompatible matrix dimensions.");
		}
		return row_count;
	}

	void display(const std::string& title) {
		std::cout << title << std::endl;
		for ( int32_t i = 0; i < row_count; ++i ) {
			for ( int32_t j = 0; j < column_count; ++j ) {
				std::cout << std::setw(9) << std::fixed << std::setprecision(4) << data[i][j];
			}
			std::cout << std::endl;
		}
		std::cout << std::endl;
	}

	double get_entry(const int32_t& row, const int32_t& col) {
		return data[row][col];
	}

	void set_entry(const int32_t& row, const int32_t& col, const double& value) {
		data[row][col] = value;
	}

	int32_t get_row_count() {
		return row_count;
	}

	int32_t get_column_count() {
		return column_count;
	}

private:
	void initialise() {
		row_count = data.size();
		column_count = data[0].size();
	}

	int32_t row_count;
	int32_t column_count;
	std::vector<std::vector<double>> data;
};

typedef std::pair<Matrix, Matrix> matrix_pair;

Matrix householder_factor(Matrix vector) {
	if ( vector.get_column_count() != 1 ) {
		throw std::invalid_argument("Incompatible matrix dimensions.");
	}

	const int32_t size = vector.size();
	Matrix result(size, size);
	for ( int32_t i = 0; i < size; ++i ) {
		for ( int32_t j = 0; j < size; ++j ) {
			result.set_entry(i, j, -2 * vector.get_entry(i, 0) * vector.get_entry(j, 0));
		}
	}

	for ( int32_t i = 0; i < size; ++i ) {
		result.set_entry(i, i, result.get_entry(i, i) + 1.0);
	}
	return result;
}

matrix_pair householder(Matrix matrix) {
	const int32_t row_count = matrix.get_row_count();
	const int32_t column_count = matrix.get_column_count();
	std::vector<Matrix> versions_of_Q;
	Matrix z(matrix);
	Matrix z1(row_count, column_count);

	for ( int32_t k = 0; k < column_count && k < row_count - 1; ++k ) {
		Matrix vectorE(row_count, 1);
		z1 = z.minor(k);
		Matrix vectorX = z1.column(k);
		double magnitudeX = vectorX.magnitude();
		if ( matrix.get_entry(k, k) > 0 ) {
			magnitudeX = -magnitudeX;
		}

		for ( int32_t i = 0; i < vectorE.size(); ++i ) {
			vectorE.set_entry(i, 0, ( i == k ) ? 1 : 0);
		}
		vectorE = vectorE.scalarMultiply(magnitudeX).add(vectorX).unit();
		versions_of_Q.emplace_back(householder_factor(vectorE));
		z = versions_of_Q[k].multiply(z1);
	}

	Matrix Q = versions_of_Q[0];
	for ( int32_t i = 1; i < column_count && i < row_count - 1; ++i ) {
		Q = versions_of_Q[i].multiply(Q);
	}

	Matrix R = Q.multiply(matrix);
	Q = Q.transpose();
	return matrix_pair(R, Q);
}

Matrix solve_upper_triangular(Matrix r, Matrix b) {
	const int32_t column_count = r.get_column_count();
	Matrix result(column_count, 1);

	for ( int32_t k = column_count - 1; k >= 0; --k ) {
		double total = 0.0;
		for ( int32_t j = k + 1; j < column_count; ++j ) {
			total += r.get_entry(k, j) * result.get_entry(j, 0);
		}
		result.set_entry(k, 0, ( b.get_entry(k, 0) - total ) / r.get_entry(k, k));
	}
	return result;
}

Matrix least_squares(Matrix vandermonde, Matrix b) {
	matrix_pair pair = householder(vandermonde);
	return solve_upper_triangular(pair.first, pair.second.transpose().multiply(b));
}

Matrix fit_polynomial(Matrix x, Matrix y, const int32_t& polynomial_degree) {
	Matrix vandermonde(x.get_column_count(), polynomial_degree + 1);
	for ( int32_t i = 0; i < x.get_column_count(); ++i ) {
		for ( int32_t j = 0; j < polynomial_degree + 1; ++j ) {
			vandermonde.set_entry(i, j, std::pow(x.get_entry(0, i), j));
		}
	}
	return least_squares(vandermonde, y.transpose());
}

int main() {
	const std::vector<std::vector<double>> data = { { 12.0, -51.0,   4.0 },
					  						        {  6.0, 167.0, -68.0 },
					  						        { -4.0,  24.0, -41.0 },
					  						        { -1.0,   1.0,   0.0 },
					  						        {  2.0,   0.0,   3.0 } };

	// Task 1
	Matrix A(data);
	A.display("Initial matrix A:");

	matrix_pair pair = householder(A);
	Matrix Q = pair.second;
	Matrix R = pair.first;

	Q.display("Matrix Q:");
	R.display("Matrix R:");

	Matrix result = Q.multiply(R);
	result.display("Matrix Q * R:");

	// Task 2
	Matrix x( std::vector<std::vector<double>>{ { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 } } );
	Matrix y( std::vector<std::vector<double>>{
		 { 1.0, 6.0, 17.0, 34.0, 57.0, 86.0, 121.0, 162.0, 209.0, 262.0, 321.0 } } );

	result = fit_polynomial(x, y, 2);
	result.display("Result of fitting polynomial:");
}

Output:

Initial matrix A:
  12.0000 -51.0000   4.0000
   6.0000 167.0000 -68.0000
  -4.0000  24.0000 -41.0000
  -1.0000   1.0000   0.0000
   2.0000   0.0000   3.0000

Matrix Q:
   0.8464  -0.3913   0.3431   0.0815   0.0781
   0.4232   0.9041  -0.0293   0.0258   0.0447
  -0.2821   0.1704   0.9329  -0.0474  -0.1374
  -0.0705   0.0140  -0.0011   0.9804  -0.1836
   0.1411  -0.0167  -0.1058  -0.1713  -0.9692

Matrix R:
  14.1774  20.6666 -13.4016
  -0.0000 175.0425 -70.0803
   0.0000   0.0000 -35.2015
  -0.0000  -0.0000  -0.0000
   0.0000   0.0000  -0.0000

Matrix Q * R:
  12.0000 -51.0000   4.0000
   6.0000 167.0000 -68.0000
  -4.0000  24.0000 -41.0000
  -1.0000   1.0000  -0.0000
   2.0000  -0.0000   3.0000

Result of fitting polynomial:
   1.0000
   2.0000
   3.0000

Common Lisp

Uses the routines m+, m-, .*, ./ from Element-wise_operations, mmul from Matrix multiplication, mtp from Matrix transposition.

Helper functions:

(defun sign (x)
  (if (zerop x)
      x
      (/ x (abs x))))

(defun norm (x)
  (let ((len (car (array-dimensions x))))
    (sqrt (loop for i from 0 to (1- len) sum (expt (aref x i 0) 2)))))

(defun make-unit-vector (dim)
  (let ((vec (make-array `(,dim ,1) :initial-element 0.0d0)))
    (setf (aref vec 0 0) 1.0d0)
    vec))

;; Return a nxn identity matrix.
(defun eye (n)
  (let ((I (make-array `(,n ,n) :initial-element 0)))
    (loop for j from 0 to (- n 1) do
          (setf (aref I j j) 1))
    I))

(defun array-range (A ma mb na nb)
  (let* ((mm (1+ (- mb ma)))
         (nn (1+ (- nb na)))
         (B (make-array `(,mm ,nn) :initial-element 0.0d0)))

    (loop for i from 0 to (1- mm) do
         (loop for j from 0 to (1- nn) do
              (setf (aref B i j)
                    (aref A (+ ma i) (+ na j)))))
    B))
 
(defun rows (A) (car  (array-dimensions A)))
(defun cols (A) (cadr (array-dimensions A)))
(defun mcol (A n) (array-range A 0 (1- (rows A)) n n))
(defun mrow (A n) (array-range A n n 0 (1- (cols A))))

(defun array-embed (A B row col)
  (let* ((ma (rows A))
         (na (cols A))
         (mb (rows B))
         (nb (cols B))
         (C  (make-array `(,ma ,na) :initial-element 0.0d0)))

    (loop for i from 0 to (1- ma) do
         (loop for j from 0 to (1- na) do
              (setf (aref C i j) (aref A i j))))

    (loop for i from 0 to (1- mb) do
         (loop for j from 0 to (1- nb) do
              (setf (aref C (+ row i) (+ col j))
                    (aref B i j))))

    C))

Main routines:

(defun make-householder (a)
  (let* ((m    (car (array-dimensions a)))
         (s    (sign (aref a 0 0)))
         (e    (make-unit-vector m))
         (u    (m+ a (.* (* (norm a) s) e)))
         (v    (./ u (aref u 0 0)))
         (beta (/ 2 (aref (mmul (mtp v) v) 0 0))))
    
    (m- (eye m)
        (.* beta (mmul v (mtp v))))))

(defun qr (A)
  (let* ((m (car  (array-dimensions A)))
         (n (cadr (array-dimensions A)))
         (Q (eye m)))

    ;; Work on n columns of A.
    (loop for i from 0 to (if (= m n) (- n 2) (- n 1)) do

         ;; Select the i-th submatrix. For i=0 this means the original matrix A.
         (let* ((B (array-range A i (1- m) i (1- n)))
                ;; Take the first column of the current submatrix B.
                (x (mcol B 0))
                ;; Create the Householder matrix for the column and embed it into an mxm identity.
                (H (array-embed (eye m) (make-householder x) i i)))

           ;; The product of all H matrices from the right hand side is the orthogonal matrix Q.
           (setf Q (mmul Q H))

           ;; The product of all H matrices with A from the LHS is the upper triangular matrix R.
           (setf A (mmul H A))))

    ;; Return Q and R.
    (values Q A)))

Example 1:

(qr #2A((12 -51 4) (6 167 -68) (-4 24 -41)))

#2A((-0.85  0.39  0.33)
    (-0.42 -0.90 -0.03)
    ( 0.28 -0.17  0.94))

#2A((-14.0  -21.0  14.0)
    (  0.0 -175.0  70.0)
    (  0.0    0.0 -35.0))

Example 2, Polynomial regression:

(defun polyfit (x y n)
  (let* ((m (cadr (array-dimensions x)))
         (A (make-array `(,m ,(+ n 1)) :initial-element 0)))
    (loop for i from 0 to (- m 1) do
          (loop for j from 0 to n do
                (setf (aref A i j)
                      (expt (aref x 0 i) j))))
    (lsqr A (mtp y))))

;; Solve a linear least squares problem by QR decomposition.
(defun lsqr (A b)
  (multiple-value-bind (Q R) (qr A)
    (let* ((n (cadr (array-dimensions R))))
      (solve-upper-triangular (array-range R                0 (- n 1) 0 (- n 1))
                              (array-range (mmul (mtp Q) b) 0 (- n 1) 0 0)))))

;; Solve an upper triangular system by back substitution.
(defun solve-upper-triangular (R b)
  (let* ((n (cadr (array-dimensions R)))
         (x (make-array `(,n 1) :initial-element 0.0d0)))

    (loop for k from (- n 1) downto 0
       do (setf (aref x k 0)
                (/ (- (aref b k 0)
                      (loop for j from (+ k 1) to (- n 1)
                         sum (* (aref R k j)
                                (aref x j 0))))
                   (aref R k k))))
    x))

;; Finally use the data:
(let ((x #2A((0 1 2 3 4 5 6 7 8 9 10)))
      (y #2A((1 6 17 34 57 86 121 162 209 262 321))))
    (polyfit x y 2))

#2A((0.999999966345088) (2.000000015144699) (2.99999999879804))

D

Translation of: Common Lisp

Uses the functions copied from Element-wise_operations, Matrix multiplication, and Matrix transposition.

import std.stdio, std.math, std.algorithm, std.traits,
       std.typecons, std.numeric, std.range, std.conv;

template elementwiseMat(string op) {
    T[][] elementwiseMat(T)(in T[][] A, in T B) pure nothrow {
        if (A.empty)
            return null;
        auto R = new typeof(return)(A.length, A[0].length);
        foreach (immutable r, const row; A)
            R[r][] = mixin("row[] " ~ op ~ "B");
        return R;
    }

    T[][] elementwiseMat(T, U)(in T[][] A, in U[][] B)
    pure nothrow if (is(Unqual!T == Unqual!U)) {
        assert(A.length == B.length);
        if (A.empty)
            return null;
        auto R = new typeof(return)(A.length, A[0].length);
        foreach (immutable r, const row; A) {
            assert(row.length == B[r].length);
            R[r][] = mixin("row[] " ~ op ~ "B[r][]");
        }
        return R;
    }
}

alias mSum = elementwiseMat!q{ + },
      mSub = elementwiseMat!q{ - },
      pMul = elementwiseMat!q{ * },
      pDiv = elementwiseMat!q{ / };

bool isRectangular(T)(in T[][] mat) pure nothrow {
    return mat.all!(r => r.length == mat[0].length);
}

T[][] matMul(T)(in T[][] a, in T[][] b) pure nothrow
in {
    assert(a.isRectangular && b.isRectangular &&
           a[0].length == b.length);
} body {
    auto result = new T[][](a.length, b[0].length);
    auto aux = new T[b.length];
    foreach (immutable j; 0 .. b[0].length) {
        foreach (immutable k; 0 .. b.length)
            aux[k] = b[k][j];
        foreach (immutable i; 0 .. a.length)
            result[i][j] = a[i].dotProduct(aux);
    }
    return result;
}

Unqual!T[][] transpose(T)(in T[][] m) pure nothrow {
    auto r = new Unqual!T[][](m[0].length, m.length);
    foreach (immutable nr, row; m)
        foreach (immutable nc, immutable c; row)
            r[nc][nr] = c;
    return r;
}

T norm(T)(in T[][] m) pure nothrow {
    return transversal(m, 0).map!q{ a ^^ 2 }.sum.sqrt;
}

Unqual!T[][] makeUnitVector(T)(in size_t dim) pure nothrow {
    auto result = new Unqual!T[][](dim, 1);
    foreach (row; result)
        row[] = 0;
    result[0][0] = 1;
    return result;
}

/// Return a nxn identity matrix.
Unqual!T[][] matId(T)(in size_t n) pure nothrow {
    auto Id = new Unqual!T[][](n, n);
    foreach (immutable r, row; Id) {
        row[] = 0;
        row[r] = 1;
    }
    return Id;
}

T[][] slice2D(T)(in T[][] A,
                 in size_t ma, in size_t mb,
                 in size_t na, in size_t nb) pure nothrow {
    auto B = new T[][](mb - ma + 1, nb - na + 1);
    foreach (immutable i, brow; B)
        brow[] = A[ma + i][na .. na + brow.length];
    return B;
}

size_t rows(T)(in T[][] A) pure nothrow { return A.length; }

size_t cols(T)(in T[][] A) pure nothrow {
    return A.length ? A[0].length : 0;
}

T[][] mcol(T)(in T[][] A, in size_t n) pure nothrow {
    return slice2D(A, 0, A.rows - 1, n, n);
}

T[][] matEmbed(T)(in T[][] A, in T[][] B,
                  in size_t row, in size_t col) pure nothrow {
    auto C = new T[][](rows(A), cols(A));
    foreach (immutable i, const arow; A)
        C[i][] = arow[]; // Some wasted copies.
    foreach (immutable i, const brow; B)
        C[row + i][col .. col + brow.length] = brow[];
    return C;
}

// Main routines ---------------

T[][] makeHouseholder(T)(in T[][] a) {
    immutable m = a.rows;
    immutable T s = a[0][0].sgn;
    immutable e = makeUnitVector!T(m);
    immutable u = mSum(a, pMul(e, a.norm * s));
    immutable v = pDiv(u, u[0][0]);
    immutable beta = 2.0 / v.transpose.matMul(v)[0][0];
    return mSub(matId!T(m), pMul(v.matMul(v.transpose), beta));
}

Tuple!(T[][],"Q", T[][],"R") QRdecomposition(T)(T[][] A) {
    immutable m = A.rows;
    immutable n = A.cols;
    auto Q = matId!T(m);

    // Work on n columns of A.
    foreach (immutable i; 0 .. (m == n ? n - 1 : n)) {
        // Select the i-th submatrix. For i=0 this means the original
        // matrix A.
        immutable B = slice2D(A, i, m - 1, i, n - 1);

        // Take the first column of the current submatrix B.
        immutable x = mcol(B, 0);

        // Create the Householder matrix for the column and embed it
        // into an mxm identity.
        immutable H = matEmbed(matId!T(m), x.makeHouseholder, i, i);

        // The product of all H matrices from the right hand side is
        // the orthogonal matrix Q.
        Q = Q.matMul(H);

        // The product of all H matrices with A from the LHS is the
        // upper triangular matrix R.
        A  = H.matMul(A);
    }

    // Return Q and R.
    return typeof(return)(Q, A);
}

// Polynomial regression ---------------

/// Solve an upper triangular system by back substitution.
T[][] solveUpperTriangular(T)(in T[][] R, in T[][] b) pure nothrow {
    immutable n = R.cols;
    auto x = new T[][](n, 1);

    foreach_reverse (immutable k; 0 .. n) {
        T tot = 0;
        foreach (immutable j; k + 1 .. n)
            tot += R[k][j] * x[j][0];
        x[k][0] = (b[k][0] - tot) / R[k][k];
    }

    return x;
}

/// Solve a linear least squares problem by QR decomposition.
T[][] lsqr(T)(T[][] A, in T[][] b) pure nothrow {
    const qr = A.QRdecomposition;
    immutable n = qr.R.cols;
    return solveUpperTriangular(
        slice2D(qr.R, 0, n - 1, 0, n - 1),
        slice2D(qr.Q.transpose.matMul(b), 0, n - 1, 0, 0));
}

T[][] polyFit(T)(in T[][] x, in T[][] y, in size_t n) pure nothrow {
    immutable size_t m = x.cols;
    auto A = new T[][](m, n + 1);
    foreach (immutable i, row; A)
        foreach (immutable j, ref item; row)
            item = x[0][i] ^^ j;
    return lsqr(A, y.transpose);
}

void main() {
    // immutable (Q, R) = QRdecomposition([[12.0, -51,   4],
    immutable qr = QRdecomposition([[12.0, -51,   4],
                                    [ 6.0, 167, -68],
                                    [-4.0,  24, -41]]);
    immutable form = "[%([%(%2.3f, %)]%|,\n %)]\n";
    writefln(form, qr.Q);
    writefln(form, qr.R);

    immutable x = [[0.0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]];
    immutable y = [[1.0, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]];
    polyFit(x, y, 2).writeln;
}

Output:

[[-0.857, 0.394, 0.331],
 [-0.429, -0.903, -0.034],
 [0.286, -0.171, 0.943]]

[[-14.000, -21.000, 14.000],
 [0.000, -175.000, 70.000],
 [0.000, -0.000, -35.000]]

[[1], [2], [3]]

F#

// QR decomposition. Nigel Galloway: January 11th., 2022
let n=[[12.0;-51.0;4.0];[6.0;167.0;-68.0];[-4.0;24.0;-41.0]]|>MathNet.Numerics.LinearAlgebra.MatrixExtensions.matrix
let g=n|>MathNet.Numerics.LinearAlgebra.Matrix.qr
printfn $"Matrix\n------\n%A{n}\nQ\n-\n%A{g.Q}\nR\n-\n%A{g.R}"

Output:

Matrix
------
DenseMatrix 3x3-Double
12  -51    4
 6  167  -68
-4   24  -41

Q
-
DenseMatrix 3x3-Double
-0.857143   0.394286  -0.331429
-0.428571  -0.902857  0.0342857
 0.285714  -0.171429  -0.942857

R
-
DenseMatrix 3x3-Double
-14   -21  14
  0  -175  70
  0     0  35

Fortran

Library: LAPACK

See the documentation for the DGEQRF and DORGQR routines. Here the example matrix is the magic square from Albrecht Dürer's Melencolia I.

program qrtask
    implicit none
    integer, parameter :: n = 4
    real(8) :: durer(n, n) = reshape(dble([ &
        16,  5,  9,  4, &
         3, 10,  6, 15, &
         2, 11,  7, 14, &
        13,  8, 12,  1  &
    ]), [n, n])
    real(8) :: q(n, n), r(n, n), qr(n, n), id(n, n), tau(n)
    integer, parameter :: lwork = 1024
    real(8) :: work(lwork)
    integer :: info, i, j
 
    q = durer
    call dgeqrf(n, n, q, n, tau, work, lwork, info)
 
    r = 0d0
    forall (i = 1:n, j = 1:n, j >= i) r(i, j) = q(i, j)
 
    call dorgqr(n, n, n, q, n, tau, work, lwork, info)
 
    qr = matmul(q, r)
    id = matmul(q, transpose(q))
 
    call show(4, durer, "A")
    call show(4, q, "Q")
    call show(4, r, "R")
    call show(4, qr, "Q*R")
    call show(4, id, "Q*Q'")
contains
    subroutine show(n, a, s)
        character(*) :: s
        integer :: n, i
        real(8) :: a(n, n)
 
        print *, s
        do i = 1, n
            print 1, a(i, :)
          1 format (*(f12.6,:,' '))
        end do
    end subroutine
end program

Output:

 A
   16.000000     3.000000     2.000000    13.000000
    5.000000    10.000000    11.000000     8.000000
    9.000000     6.000000     7.000000    12.000000
    4.000000    15.000000    14.000000     1.000000
 Q
   -0.822951     0.376971     0.361447    -0.223607
   -0.257172    -0.454102    -0.526929    -0.670820
   -0.462910    -0.060102    -0.576283     0.670820
   -0.205738    -0.805029     0.509510     0.223607
 R
  -19.442222   -10.904103   -10.595497   -18.516402
    0.000000   -15.846152   -15.932298    -0.258437
    0.000000     0.000000    -1.974168    -5.922505
    0.000000     0.000000     0.000000    -0.000000
 Q*R
   16.000000     3.000000     2.000000    13.000000
    5.000000    10.000000    11.000000     8.000000
    9.000000     6.000000     7.000000    12.000000
    4.000000    15.000000    14.000000     1.000000
 Q*Q'
    1.000000    -0.000000    -0.000000     0.000000
   -0.000000     1.000000     0.000000     0.000000
   -0.000000     0.000000     1.000000    -0.000000
    0.000000     0.000000    -0.000000     1.000000

Futhark

import "lib/github.com/diku-dk/linalg/linalg"

module linalg_f64 = mk_linalg f64

let eye (n: i32): [n][n]f64 =
  let arr = map (\ind -> let (i,j) = (ind/n,ind%n) in if (i==j) then 1.0 else 0.0) (iota (n*n))
  in unflatten n n arr

let norm v = linalg_f64.dotprod v v |> f64.sqrt

let qr [n] [m] (a: [m][n]f64): ([m][m]f64, [m][n]f64) =

  let make_householder [d] (x: [d]f64): [d][d]f64 =
    let div = if x[0] > 0 then x[0] + norm x else x[0] - norm x
    let v = map (/div) x
    let v[0] = 1
    let fac = 2.0 / linalg_f64.dotprod v v
    in map2 (map2 (-)) (eye d) (map (map (*fac)) (linalg_f64.outer v v))

  let step ((x,y):([m][m]f64,[m][n]f64)) (i:i32): ([m][m]f64,[m][n]f64) =
    let h = eye m
    let h[i:m,i:m] = make_householder y[i:m,i]
    let q': [m][m]f64 = linalg_f64.matmul x h
    let a': [m][n]f64 = linalg_f64.matmul h y
    in (q',a')

  let q = eye m
  in foldl step (q,a) (iota n)

entry main = qr [[12.0, -51.0, 4.0],[6.0, 167.0, -68.0],[-4.0, 24.0, -41.0]]

Output:

$ ./qr 
[[-0.857143f64, 0.394286f64, -0.331429f64], [-0.428571f64, -0.902857f64, 0.034286f64], [0.285714f64, -0.171429f64, -0.942857f64]]
[[-14.000000f64, -21.000000f64, 14.000000f64], [0.000000f64, -175.000000f64, 70.000000f64], [-0.000000f64, 0.000000f64, 35.000000f64]]

Go

Method of task description, library go.matrix

Translation of: Common Lisp

A fairly close port of the Common Lisp solution, this solution uses the go.matrix library for supporting functions. Note though, that go.matrix has QR decomposition, as shown in the Go solution to Polynomial regression. The solution there is coded more directly than by following the CL example here. Similarly, examination of the go.matrix QR source shows that it computes the decomposition more directly.

package main

import (
    "fmt"
    "math"

    "github.com/skelterjohn/go.matrix"
)

func sign(s float64) float64 {
    if s > 0 {
        return 1
    } else if s < 0 {
        return -1
    }
    return 0
}

func unitVector(n int) *matrix.DenseMatrix {
    vec := matrix.Zeros(n, 1)
    vec.Set(0, 0, 1)
    return vec
}

func householder(a *matrix.DenseMatrix) *matrix.DenseMatrix {
    m := a.Rows()
    s := sign(a.Get(0, 0))
    e := unitVector(m)
    u := matrix.Sum(a, matrix.Scaled(e, a.TwoNorm()*s))
    v := matrix.Scaled(u, 1/u.Get(0, 0))
    // (error checking skipped in this solution)
    prod, _ := v.Transpose().TimesDense(v)
    β := 2 / prod.Get(0, 0)

    prod, _ = v.TimesDense(v.Transpose())
    return matrix.Difference(matrix.Eye(m), matrix.Scaled(prod, β))
}

func qr(a *matrix.DenseMatrix) (q, r *matrix.DenseMatrix) {
    m := a.Rows()
    n := a.Cols()
    q = matrix.Eye(m)

    last := n - 1
    if m == n {
        last--
    }
    for i := 0; i <= last; i++ {
        // (copy is only for compatibility with an older version of gomatrix)
        b := a.GetMatrix(i, i, m-i, n-i).Copy()
        x := b.GetColVector(0)
        h := matrix.Eye(m)
        h.SetMatrix(i, i, householder(x))
        q, _ = q.TimesDense(h)
        a, _ = h.TimesDense(a)
    }
    return q, a
}

func main() {
    // task 1: show qr decomp of wp example
    a := matrix.MakeDenseMatrixStacked([][]float64{
        {12, -51, 4},
        {6, 167, -68},
        {-4, 24, -41}})
    q, r := qr(a)
    fmt.Println("q:\n", q)
    fmt.Println("r:\n", r)

    // task 2: use qr decomp for polynomial regression example
    x := matrix.MakeDenseMatrixStacked([][]float64{
        {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}})
    y := matrix.MakeDenseMatrixStacked([][]float64{
        {1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}})
    fmt.Println("\npolyfit:\n", polyfit(x, y, 2))
}

func polyfit(x, y *matrix.DenseMatrix, n int) *matrix.DenseMatrix {
    m := x.Cols()
    a := matrix.Zeros(m, n+1)
    for i := 0; i < m; i++ {
        for j := 0; j <= n; j++ {
            a.Set(i, j, math.Pow(x.Get(0, i), float64(j)))
        }
    }
    return lsqr(a, y.Transpose())
}

func lsqr(a, b *matrix.DenseMatrix) *matrix.DenseMatrix {
    q, r := qr(a)
    n := r.Cols()
    prod, _ := q.Transpose().TimesDense(b)
    return solveUT(r.GetMatrix(0, 0, n, n), prod.GetMatrix(0, 0, n, 1))
}

func solveUT(r, b *matrix.DenseMatrix) *matrix.DenseMatrix {
    n := r.Cols()
    x := matrix.Zeros(n, 1)
    for k := n - 1; k >= 0; k-- {
        sum := 0.
        for j := k + 1; j < n; j++ {
            sum += r.Get(k, j) * x.Get(j, 0)
        }
        x.Set(k, 0, (b.Get(k, 0)-sum)/r.Get(k, k))
    }
    return x
}

Output:

q:
 {-0.857143,  0.394286,  0.331429,
 -0.428571, -0.902857, -0.034286,
  0.285714, -0.171429,  0.942857}
r:
 { -14,  -21,   14,
    0, -175,   70,
    0,    0,  -35}

polyfit:
 {1,
 2,
 3}

Library QR, gonum/matrix

package main

import (
    "fmt"

    "github.com/gonum/matrix/mat64"
)

func main() {
    // task 1: show qr decomp of wp example
    a := mat64.NewDense(3, 3, []float64{
        12, -51, 4,
        6, 167, -68,
        -4, 24, -41,
    })
    var qr mat64.QR
    qr.Factorize(a)
    var q, r mat64.Dense
    q.QFromQR(&qr)
    r.RFromQR(&qr)
    fmt.Printf("q: %.3f\n\n", mat64.Formatted(&q, mat64.Prefix("   ")))
    fmt.Printf("r: %.3f\n\n", mat64.Formatted(&r, mat64.Prefix("   ")))

    // task 2: use qr decomp for polynomial regression example
    x := []float64{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}
    y := []float64{1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}
    a = Vandermonde(x, 2)
    b := mat64.NewDense(11, 1, y)
    qr.Factorize(a)
    var f mat64.Dense
    f.SolveQR(&qr, false, b)
    fmt.Printf("polyfit: %.3f\n",
        mat64.Formatted(&f, mat64.Prefix("         ")))
}

func Vandermonde(a []float64, degree int) *mat64.Dense {
    x := mat64.NewDense(len(a), degree+1, nil)
    for i := range a {
        for j, p := 0, 1.; j <= degree; j, p = j+1, p*a[i] {
            x.Set(i, j, p)
        }
    }
    return x
}

Output:

q: ⎡-0.857   0.394   0.331⎤
   ⎢-0.429  -0.903  -0.034⎥
   ⎣ 0.286  -0.171   0.943⎦

r: ⎡ -14.000   -21.000    14.000⎤
   ⎢   0.000  -175.000    70.000⎥
   ⎣   0.000     0.000   -35.000⎦

polyfit: ⎡1.000⎤
         ⎢2.000⎥
         ⎣3.000⎦

Haskell

Square matrices only; decompose A and solve Rx = q by back substitution

import Data.List
import Text.Printf (printf)

eps = 1e-6 :: Double

-- a matrix is represented as a list of columns
mmult :: Num a => [[a]] -> [[a]] -> [[a]] 
nth :: Num a => [[a]] -> Int -> Int -> a
mmult_num :: Num a => [[a]] -> a -> [[a]]
madd :: Num a => [[a]] -> [[a]] -> [[a]]
idMatrix :: Num a => Int -> Int -> [[a]]

adjustWithE :: [[Double]] -> Int -> [[Double]]

mmult a b = [ [ sum $ zipWith (*) ak bj | ak <- (transpose a) ] | bj <- b ]
nth mA i j = (mA !! j) !! i
mmult_num mA n = map (\c -> map (*n) c) mA
madd mA mB = zipWith (\c1 c2 -> zipWith (+) c1 c2) mA mB
idMatrix n m = [ [if (i==j) then 1 else 0 | i <- [1..n]] | j <- [1..m]]

adjustWithE mA n = let lA = length mA in
    (idMatrix n (n - lA)) ++ (map (\c -> (take (n - lA) (repeat 0.0)) ++ c ) mA)

-- auxiliary functions
sqsum :: Floating a => [a] -> a
norm :: Floating a => [a] -> a
epsilonize :: [[Double]] -> [[Double]]

sqsum a = foldl (\x y -> x + y*y) 0 a
norm a = sqrt $! sqsum a
epsilonize mA = map (\c -> map (\x -> if abs x <= eps then 0 else x) c) mA

-- Householder transformation; householder A = (Q, R)
uTransform :: [Double] -> [Double]
hMatrix :: [Double] -> Int -> Int -> [[Double]]
householder :: [[Double]] -> ([[Double]], [[Double]])

-- householder_rec Q R A
householder_rec :: [[Double]] -> [[Double]] -> Int -> ([[Double]], [[Double]])

uTransform a = let t = (head a) + (signum (head a))*(norm a) in
    1 : map (\x -> x/t) (tail a)

hMatrix a n i = let u = uTransform (drop i a) in
    madd
        (idMatrix (n-i) (n-i))
        (mmult_num
            (mmult [u] (transpose [u]))
            ((/) (-2) (sqsum u)))

householder_rec mQ mR 0 = (mQ, mR)
householder_rec mQ mR n = let mSize = length mR in
    let mH = adjustWithE (hMatrix (mR!!(mSize - n)) mSize (mSize - n)) mSize in
        householder_rec (mmult mQ mH) (mmult mH mR) (n - 1)

householder mA = let mSize = length mA in
    let (mQ, mR) = householder_rec (idMatrix mSize mSize) mA mSize in
        (epsilonize mQ, epsilonize mR)

backSubstitution :: [[Double]] -> [Double] -> [Double] -> [Double]
backSubstitution mR [] res = res
backSubstitution mR@(hR:tR) q@(h:t) res =
    let x = (h / (head hR)) in
        backSubstitution
            (map tail tR)
            (tail (zipWith (-) q (map (*x) hR)))
            (x : res)

showMatrix :: [[Double]] -> String
showMatrix mA =
    concat $ intersperse "\n"
        (map (\x -> unwords $ printf "%10.4f" <$> (x::[Double])) (transpose mA))

mY = [[12, 6, -4], [-51, 167, 24], [4, -68, -41]] :: [[Double]]
q = [21, 245, 35] :: [Double]
main = let (mQ, mR) = householder mY in
    putStrLn ("Q: \n" ++ showMatrix mQ) >>
    putStrLn ("R: \n" ++ showMatrix mR) >>
    putStrLn ("q: \n" ++ show q) >>
    putStrLn ("x: \n" ++ show (backSubstitution (reverse (map reverse mR)) (reverse q) []))

Output:

Q: 
   -0.8571     0.3943    -0.3314
   -0.4286    -0.9029     0.0343
    0.2857    -0.1714    -0.9429
R: 
  -14.0000   -21.0000    14.0000
    0.0000  -175.0000    70.0000
    0.0000     0.0000    35.0000
q: 
[21.0,245.0,35.0]
x: 
[1.0000000000000004,-0.9999999999999999,1.0]

QR decomposition with Numeric.LinearAlgebra

import Numeric.LinearAlgebra

a :: Matrix R
a = (3><3) 
  [ 12, -51,   4
  ,  6, 167, -68
  , -4,  24, -41]

main = do
  print $ qr a

Output:

((3><3)
 [ -0.8571428571428572,   0.3942857142857143,   0.33142857142857146
 , -0.4285714285714286,  -0.9028571428571428, -3.428571428571427e-2
 , 0.28571428571428575, -0.17142857142857137,    0.9428571428571428 ],(3><3)
 [ -14.0,               -21.0, 14.000000000000002
 ,   0.0, -175.00000000000003,  70.00000000000001
 ,   0.0, 

J

Solution (built-in):

   QR =: 128!:0

Solution (custom implementation):

   mp=: +/ . *  NB. matrix product
   h =: +@|:    NB. conjugate transpose

   QR=: 3 : 0
    n=.{:$A=.y
    if. 1>:n do.
     A ((% {.@,) ; ]) %:(h A) mp A
    else.
     m =.>.n%2
     A0=.m{."1 A
     A1=.m}."1 A
     'Q0 R0'=.QR A0
     'Q1 R1'=.QR A1 - Q0 mp T=.(h Q0) mp A1
     (Q0,.Q1);(R0,.T),(-n){."1 R1
    end.
   )

Example:

   QR 12 _51 4,6 167 _68,:_4 24 _41
+-----------------------------+----------+
| 0.857143 _0.394286 _0.331429|14  21 _14|
| 0.428571  0.902857 0.0342857| 0 175 _70|
|_0.285714  0.171429 _0.942857| 0   0  35|
+-----------------------------+----------+

Example (polynomial fitting using QR reduction):

   X=:i.# Y=:1 6 17 34 57 86 121 162 209 262 321
   'Q R'=: QR X ^/ i.3
   R %.~(|:Q)+/ .* Y
1 2 3

Notes:J offers a built-in QR decomposition function, 128!:0. If J did not offer this function as a built-in, it could be written in J along the lines of the second version, which is covered in an essay on the J wiki.

Java

JAMA

Using the JAMA library. Compile with: javac -cp Jama-1.0.3.jar Decompose.java.

import Jama.Matrix;
import Jama.QRDecomposition;

public class Decompose {
    public static void main(String[] args) {
        var matrix = new Matrix(new double[][] {
            {12, -51,   4},
            { 6, 167, -68},
            {-4,  24, -41},
        });

        var qr = new QRDecomposition(matrix);
        qr.getQ().print(10, 4);
        qr.getR().print(10, 4);
    }
}

Output:

     -0.8571      0.3943     -0.3314
     -0.4286     -0.9029      0.0343
      0.2857     -0.1714     -0.9429


    -14.0000    -21.0000     14.0000
      0.0000   -175.0000     70.0000
      0.0000      0.0000     35.0000

Colt

Using the Colt library. Compile with: javac -cp colt.jar Decompose.java.

import cern.colt.matrix.impl.DenseDoubleMatrix2D;
import cern.colt.matrix.linalg.QRDecomposition;

public class Decompose {
    public static void main(String[] args) {
        var a = new DenseDoubleMatrix2D(new double[][] {
            {12, -51,   4},
            { 6, 167, -68},
            {-4,  24, -41}
        });
        var qr = new QRDecomposition(a);
        System.out.println(qr.getQ());
        System.out.println();
        System.out.println(qr.getR());
    }
}

Output:

3 x 3 matrix
-0.857143  0.394286 -0.331429
-0.428571 -0.902857  0.034286
 0.285714 -0.171429 -0.942857

3 x 3 matrix
-14  -21 14
  0 -175 70
  0    0 35

Apache Commons Math

Using the Apache Commons Math library.

Compile with: javac -cp commons-math3-3.6.1.jar Decompose.java.

import java.util.Locale;

import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.QRDecomposition;
import org.apache.commons.math3.linear.RealMatrix;

public class Decompose {
    public static void main(String[] args) {
        var a = new Array2DRowRealMatrix(new double[][] {
            {12, -51,   4},
            { 6, 167, -68},
            {-4,  24, -41}
        });
                                                         
        var qr = new QRDecomposition(a);
        print(qr.getQ());
        System.out.println();
        print(qr.getR());
    }
    
    public static void print(RealMatrix a) {
        for (double[] u: a.getData()) {
            System.out.print("[ ");
            for (double x: u) {
                System.out.printf(Locale.ROOT, "%10.4f ", x);
            }
            System.out.println("]");
        }
    }
}

Output:

[    -0.8571     0.3943    -0.3314 ]
[    -0.4286    -0.9029     0.0343 ]
[     0.2857    -0.1714    -0.9429 ]

[   -14.0000   -21.0000    14.0000 ]
[     0.0000  -175.0000    70.0000 ]
[     0.0000     0.0000    35.0000 ]

la4j

Using the la4j library. Compile with: javac -cp la4j-0.6.0.jar Decompose.java.

import org.la4j.Matrix;
import org.la4j.decomposition.QRDecompositor;

public class Decompose {
    public static void main(String[] args) {
        var a = Matrix.from2DArray(new double[][] {
            {12, -51,   4},
            { 6, 167, -68},
            {-4,  24, -41},
        });
        
        Matrix[] qr = new QRDecompositor(a).decompose();
        System.out.println(qr[0]);
        System.out.println(qr[1]);
    }
}

Output:

-0,857  0,394 -0,331
-0,429 -0,903  0,034
 0,286 -0,171 -0,943

-14,000  -21,000 14,000
  0,000 -175,000 70,000
  0,000    0,000 35,000

Without external libraries

import java.util.ArrayList;
import java.util.List;

public final class QRDecomposition {

	public static void main(String[] aArgs) {
		final double[][] data = new double [][] { { 12.0, -51.0,   4.0 },
				  						          {  6.0, 167.0, -68.0 },
				  						          { -4.0,  24.0, -41.0 },
				  						          { -1.0,   1.0,   0.0 },
				  						          {  2.0,   0.0,   3.0 } };
		
		// Task 1		  						          
		Matrix A = new Matrix(data);				  						 
		A.display("Initial matrix A:");
		
		MatrixPair pair = householder(A);		
		Matrix Q = pair.q;
		Matrix R = pair.r;
		
		Q.display("Matrix Q:"); 
		R.display("Matrix R:");
		
		Matrix result = Q.multiply(R);
		result.display("Matrix Q * R:");
		
		// Task 2
		Matrix x = new Matrix ( new double[][] { { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 } } );
		Matrix y = new Matrix(
			new double[][] { { 1.0, 6.0, 17.0, 34.0, 57.0, 86.0, 121.0, 162.0, 209.0, 262.0, 321.0 } } );
		
		result = fitPolynomial(x, y, 2);
		result.display("Result of fitting polynomial:");
	}
	
	private static MatrixPair householder(Matrix aMatrix) {
		final int rowCount = aMatrix.getRowCount();
		final int columnCount = aMatrix.getColumnCount();
		List<Matrix> versionsOfQ = new ArrayList<Matrix>(rowCount);
		Matrix z = new Matrix(aMatrix);
		Matrix z1 = new Matrix(rowCount, columnCount);
		
		for ( int k = 0; k < columnCount && k < rowCount - 1; k++ ) {
			Matrix vectorE = new Matrix(rowCount, 1);		    
		    z1 = z.minor(k);
		    Matrix vectorX = z1.column(k);	    	  
		    double magnitudeX = vectorX.magnitude();
		    if ( aMatrix.getEntry(k, k) > 0 ) {
		    	magnitudeX = -magnitudeX;
		    }

		    for ( int i = 0; i < vectorE.size(); i++ ) {
		        vectorE.setEntry(i, 0, ( i == k ) ? 1 : 0);
		    }			
		    vectorE = vectorE.scalarMultiply(magnitudeX).add(vectorX).unit();	    
		    versionsOfQ.add(householderFactor(vectorE));
		    z = versionsOfQ.get(k).multiply(z1);
		}
		
		Matrix Q = versionsOfQ.get(0);
	    for ( int i = 1; i < columnCount && i < rowCount - 1; i++ ) {
	    	Q = versionsOfQ.get(i).multiply(Q);
	    }

	    Matrix R = Q.multiply(aMatrix);
	    Q = Q.transpose();
	    return new MatrixPair(R, Q);
	}
	
	public static Matrix householderFactor(Matrix aVector) {
    	if ( aVector.getColumnCount() != 1 ) {
        	throw new RuntimeException("Incompatible matrix dimensions.");
        } 
    	
    	final int size = aVector.size();
    	Matrix result = new Matrix(size, size);
	    for ( int i = 0; i < size; i++ ) {
	        for ( int j = 0; j < size; j++ ) {
	        	result.setEntry(i, j, -2 * aVector.getEntry(i, 0) * aVector.getEntry(j, 0));
	        }
	    }
	    
	    for ( int i = 0; i < size; i++ ) {
	    	result.setEntry(i, i, result.getEntry(i, i) + 1.0);
	    }
	    return result;
    }
	
	private static Matrix fitPolynomial(Matrix aX, Matrix aY, int aPolynomialDegree) {
	    Matrix vandermonde = new Matrix(aX.getColumnCount(), aPolynomialDegree + 1);
	    for ( int i = 0; i < aX.getColumnCount(); i++ ) {
	        for ( int j = 0; j < aPolynomialDegree + 1; j++ ) {
	            vandermonde.setEntry(i, j, Math.pow(aX.getEntry(0, i), j));
	        }
		}
	    return leastSquares(vandermonde, aY.transpose());
	}
	
	private static Matrix leastSquares(Matrix aVandermonde, Matrix aB) {
		MatrixPair pair = householder(aVandermonde);
		return solveUpperTriangular(pair.r, pair.q.transpose().multiply(aB));
	}
	
	private static Matrix solveUpperTriangular(Matrix aR, Matrix aB) {
		final int columnCount = aR.getColumnCount();
	    Matrix result = new Matrix(columnCount, 1);

	    for ( int k = columnCount - 1; k >= 0; k-- ) {
	        double total = 0.0;
	        for ( int j = k + 1; j < columnCount; j++ ) {
	            total += aR.getEntry(k, j) * result.getEntry(j, 0);
	        }
	        result.setEntry(k, 0, ( aB.getEntry(k, 0) - total ) / aR.getEntry(k, k));
	    }
	    return result;
	}
	
	private static record MatrixPair(Matrix r, Matrix q) {}

}

final class Matrix {   

    public Matrix(double[][] aData) {
        rowCount = aData.length;
        columnCount = aData[0].length;
        data = new double[rowCount][columnCount];
        for ( int i = 0; i < rowCount; i++ ) {
            for ( int j = 0; j < columnCount; j++ ) {
                data[i][j] = aData[i][j];
            }
        }
    }
    
    public Matrix(Matrix aMatrix) {
    	this(aMatrix.data);
    } 
    
    public Matrix(int aRowCount, int aColumnCount) {
        this( new double[aRowCount][aColumnCount] );
    }   

    public Matrix add(Matrix aOther) {
        if ( aOther.rowCount != rowCount || aOther.columnCount != columnCount ) {
        	throw new IllegalArgumentException("Incompatible matrix dimensions.");
        }
        
        Matrix result = new Matrix(data);
        for ( int i = 0; i < rowCount; i++ ) {
            for ( int j = 0; j < columnCount; j++ ) {
                result.data[i][j] = data[i][j] + aOther.data[i][j];
            }
        }
        return result;
    }

    public Matrix multiply(Matrix aOther) {
        if ( columnCount != aOther.rowCount ) {
        	throw new IllegalArgumentException("Incompatible matrix dimensions.");
        }
        
        Matrix result = new Matrix(rowCount, aOther.columnCount);
        for ( int i = 0; i < rowCount; i++ ) {
            for ( int j = 0; j < aOther.columnCount; j++ ) {
                for ( int k = 0; k < rowCount; k++ ) {
                    result.data[i][j] += data[i][k] * aOther.data[k][j];
                }
            }
        }
        return result;
    }
    
    public Matrix transpose() {
        Matrix result = new Matrix(columnCount, rowCount);
        for ( int i = 0; i < rowCount; i++ ) {
            for ( int j = 0; j < columnCount; j++ ) {
                result.data[j][i] = data[i][j];
            }
        }
        return result;
    }
        
    public Matrix minor(int aIndex) {
    	Matrix result = new Matrix(rowCount, columnCount);
        for ( int i = 0; i < aIndex; i++ ) {
            result.setEntry(i, i, 1.0);
        }
        
        for ( int i = aIndex; i < rowCount; i++ ) {
            for ( int j = aIndex; j < columnCount; j++ ) {
        	    result.setEntry(i, j, data[i][j]);
            }
        }
        return result;
    }	
	
	public Matrix column(int aIndex) {
		Matrix result = new Matrix(rowCount, 1);
		for ( int i = 0; i < rowCount; i++ ) {
		    result.setEntry(i, 0, data[i][aIndex]);
		}
		return result;		
	}
	
	public Matrix scalarMultiply(double aValue) {
		if ( columnCount != 1 ) { 			
			throw new IllegalArgumentException("Incompatible matrix dimension.");
		}
		
		Matrix result = new Matrix(rowCount, columnCount);
		for ( int i = 0; i < rowCount; i++ ) {
		    result.data[i][0] = data[i][0] * aValue;
	    }
		return result;
	}
	
	public Matrix unit() {
		if ( columnCount != 1 ) { 			
			throw new IllegalArgumentException("Incompatible matrix dimensions.");
		}
		
		final double magnitude = magnitude();
		Matrix result = new Matrix(rowCount, columnCount);
		for ( int i = 0; i < rowCount; i++ ) {
		    result.data[i][0] = data[i][0] / magnitude;
	    }
		return result;
	}
	
	public double magnitude() {
		if ( columnCount != 1 ) { 			
			throw new IllegalArgumentException("Incompatible matrix dimensions.");
		}
		
		double norm = 0.0;
	    for ( int i = 0; i < data.length; i++ ) {
		    norm += data[i][0] * data[i][0];
	    }
	    return Math.sqrt(norm);		
	}
	
	public int size() {
		if ( columnCount != 1 ) { 			
			throw new IllegalArgumentException("Incompatible matrix dimensions.");
		}
		return rowCount;
	}
	
	public void display(String aTitle) {
		System.out.println(aTitle);
		for ( int i = 0; i < rowCount; i++ ) {
		    for ( int j = 0; j < columnCount; j++ ) {
		    	System.out.print(String.format("%9.4f", data[i][j]));
		    }
		    System.out.println();
		}
		System.out.println();
	}
	
    public double getEntry(int aRow, int aColumn) {
    	return data[aRow][aColumn];
    }
    
    public void setEntry(int aRow, int aColumn, double aValue) {
    	data[aRow][aColumn] = aValue;
    }
    
    public int getRowCount() {
    	return rowCount;
    }
    
    public int getColumnCount() {
    	return columnCount;
    }    
    
    private final int rowCount;
    private final int columnCount;
    private final double[][] data;
    
}

Output:

Initial matrix A:
  12.0000 -51.0000   4.0000
   6.0000 167.0000 -68.0000
  -4.0000  24.0000 -41.0000
  -1.0000   1.0000   0.0000
   2.0000   0.0000   3.0000

Matrix Q:
   0.8464  -0.3913   0.3431   0.0815   0.0781
   0.4232   0.9041  -0.0293   0.0258   0.0447
  -0.2821   0.1704   0.9329  -0.0474  -0.1374
  -0.0705   0.0140  -0.0011   0.9804  -0.1836
   0.1411  -0.0167  -0.1058  -0.1713  -0.9692

Matrix R:
  14.1774  20.6666 -13.4016
  -0.0000 175.0425 -70.0803
   0.0000   0.0000 -35.2015
  -0.0000  -0.0000  -0.0000
   0.0000   0.0000  -0.0000

Matrix Q * R:
  12.0000 -51.0000   4.0000
   6.0000 167.0000 -68.0000
  -4.0000  24.0000 -41.0000
  -1.0000   1.0000  -0.0000
   2.0000  -0.0000   3.0000

Result of fitting polynomial:
   1.0000
   2.0000
   3.0000

jq

Adapted from Wren

Works with: jq

Also works with gojq, the Go implementation of jq.

General utilities

def sum(s): reduce s as $_ (0; . + $_);

# Sum of squares
def ss(s): sum(s|.*.);

# Create an m x n matrix
def matrix(m; n; init):
  if m == 0 then []
  elif m == 1 then [range(0;n) | init]
  elif m > 0 then
    matrix(1;n;init) as $row
    | [range(0;m) | $row ]
  else error("matrix\(m);_;_) invalid")
  end;

def dot_product(a; b):
  reduce range(0;a|length) as $i (0; . + (a[$i] * b[$i]) );

# A and B should both be numeric matrices, A being m by n, and B being n by p.
def multiply($A; $B):
  ($B[0]|length) as $p
  | ($B|transpose) as $BT
  | reduce range(0; $A|length) as $i
       ([];
       reduce range(0; $p) as $j 
         (.;
          .[$i][$j] = dot_product( $A[$i]; $BT[$j] ) ));

# $ndec decimal places
def round($ndec):
  def rpad: tostring | ($ndec - length) as $l | . + ("0" * $l);
  def abs: if . < 0 then -. else . end;
  pow(10; $ndec) as $p
  | round as $round
  | if $p * ((. - $round)|abs) < 0.1
    then ($round|tostring)  + "." + ($ndec * "0")
    else  . * $p | round / $p
    | tostring
    | capture("(?<left>[^.]*)[.](?<right>.*)")
    | .left + "." + (.right|rpad)
    end;

# pretty-print a 2-d matrix
def pp($ndec; $width):
  def pad(n): tostring | (n - length) * " " + .;
  def row: map(round($ndec) | pad($width)) | join(" ");
  reduce .[] as $row (""; . + "\n\($row|row)");

QR-Decomposition

def minor($x; $d):
   ($x|length) as $nr
   | ($x[0]|length) as $nc
   | reduce range(0; $d) as $i (matrix($nr;$nc;0); .[$i][$i] = 1)
   | reduce range($d; $nr) as $i (.;
       reduce range($d;$nc) as $j (.; .[$i][$j] = $x[$i][$j] ) );

def vmadd($a; $b; $s):
  reduce range (0; $a|length) as $i ([];
    .[$i] = $a[$i] + $s * $b[$i] );

def vmul($v):
  ($v|length) as $n
  | reduce range(0;$n) as $i (null;
      reduce range(0;$n) as $j (.; .[$i][$j] = -2 * $v[$i] * $v[$j] ))
  | reduce range(0;$n) as $i (.; .[$i][$i] += 1 );

def vnorm($x):
  sum($x[] | .*.) | sqrt;

def vdiv($x; $d):
  [range (0;$x|length) | $x[.] / $d];

def mcol($m; $c):
  [range (0;$m|length) | $m[.][$c]];

def householder($m):
    ($m|length) as $nr
    | ($m[0]|length) as $nc
    | { q: [], # $nr
        z: $m,
        k: 0 }
    | until( .k >= $nc or .k >= $nr-1;
        .z = minor(.z; .k)
        | .x = mcol(.z; .k)
        | .a = vnorm(.x)
        | if ($m[.k][.k] > 0) then .a = -.a else . end
        | .e = [range (0; $nr) as $i | if ($i == .k) then 1 else 0 end]
        | .e = vmadd(.x; .e; .a)
        | .e = vdiv(.e; vnorm(.e))
        | .q[.k] = vmul(.e)
        | .z = multiply(.q[.k]; .z)
        | .k += 1 )
    | .Q = .q[0]
    | .R = multiply(.q[0]; $m)
    | .i = 1
    | until (.i >= $nc or .i >= $nr-1;
        .Q = multiply(.q[.i]; .Q)
        | .i += 1 )
    | .R = multiply(.Q; $m)
    | .Q |= transpose
    | [.Q, .R] ;

def x: [
  [12, -51,   4],
  [ 6, 167, -68],
  [-4,  24, -41],
  [-1,   1,   0],
  [ 2,   0,   3]
];

def task:
  def pp: pp(3;8);

  # Assume $a and $b are conformal
  def ssd($a; $b):
    [$a[][]] as $a
    | [$b[][]] as $b
    | ss( range(0;$a|length) | $a[.] - $b[.] );

  householder(x) as [$Q, $R]
  | multiply($Q; $R) as $m
  | "Q:",      ($Q|pp),
   "\nR:",     ($R|pp),
   "\nQ * R:", ($m|pp),
   "\nSum of squared discrepancies: \(ssd(x; $m))" 
;

task

Output:

Q:

   0.846   -0.391    0.343    0.082    0.078
   0.423    0.904   -0.029    0.026    0.045
  -0.282    0.170    0.933   -0.047   -0.137
  -0.071    0.014   -0.001    0.980   -0.184
   0.141   -0.017   -0.106   -0.171   -0.969

R:

  14.177   20.667  -13.402
  -0.000  175.043  -70.080
   0.000    0.000  -35.202
  -0.000   -0.000   -0.000
   0.000    0.000   -0.000

Q * R:

  12.000  -51.000    4.000
   6.000  167.000  -68.000
  -4.000   24.000  -41.000
  -1.000    1.000   -0.000
   2.000   -0.000    3.000

Sum of squared discrepancies: 1.1675699109208862e-26

Julia

Built-in function

Q, R = qr([12 -51 4; 6 167 -68; -4 24 -41])

Output:

(
3x3 Array{Float64,2}:
 -0.857143   0.394286   0.331429 
 -0.428571  -0.902857  -0.0342857
  0.285714  -0.171429   0.942857 ,

3x3 Array{Float64,2}:
 -14.0   -21.0   14.0
   0.0  -175.0   70.0
   0.0     0.0  -35.0)

Maple

with(LinearAlgebra):
A:=<12,-51,4;6,167,-68;-4,24,-41>:
Q,R:=QRDecomposition(A):
Q;
R;

Output:

     [        -69      -58 ]
     [6/7     ---      --- ]
     [        175      175 ]
     [                     ]
     [        158          ]
     [3/7     ---     6/175]
     [        175          ]
     [                     ]
     [                 -33 ]
     [-2/7    6/35     --- ]
     [                 35  ]


       [14     21    -14]
       [                ]
       [ 0    175    -70]
       [                ]
       [ 0      0     35]

Mathematica/Wolfram Language

{q,r}=QRDecomposition[{{12, -51, 4}, {6, 167, -68}, {-4, 24, -41}}];
q//MatrixForm

-> 6/7 3/7 -(2/7)
-69/175 158/175 6/35
-58/175 6/175 -33/35

r//MatrixForm
-> 14 21 -14
   0  175 -70
   0  0  35

MATLAB / Octave

 A = [12 -51   4 
       6 167 -68
      -4  24 -41];
 [Q,R]=qr(A)

Output:

Q =

   0.857143  -0.394286  -0.331429
   0.428571   0.902857   0.034286
  -0.285714   0.171429  -0.942857

R =

    14    21   -14
     0   175   -70
     0     0    35

Maxima

load(lapack)$   /* This may hang up in wxMaxima, if this happens, use xMaxima or plain Maxima in a terminal */

a: matrix([12, -51,   4],
          [ 6, 167, -68],
          [-4,  24, -41])$

[q, r]: dgeqrf(a)$

mat_norm(q . r - a, 1);
4.2632564145606011E-14

/* Note: the lapack package is a lisp translation of the fortran lapack library */

For an exact or arbitrary precision solution:

load("linearalgebra")$
load("eigen")$
unitVector(n) := ematrix(n,1,1,1,1);
signValue(r) := block([s:sign(r)],
  if s='pos then 1 else if s='zero then 0 else -1);
householder(a) := block([m : length(a),u,v,beta],
  u : a + sqrt(a .  a)*signValue(a[1,1])*unitVector(m),
  v : u / u[1,1],
  beta : 2/(v . v),
  diagmatrix(m,1) - beta*transpose(v . transpose(v)));
getSubmatrix(obj,i1,j1,i2,j2) :=
genmatrix(lambda([i,j], obj[i+i1-1,j+j1-1]),i2-i1+1,j2-j1+1);
setSubmatrix(obj,i1,j1,subobj) := block([m,n],
  [m,n] : matrix_size(subobj),
  for i: 0 thru m-1 do
  (for j: 0 thru n-1 do
    obj[i1+i,j1+j] : subobj[i+1,j+1]));
qr(obj) := block([m,n,qm,rm,i],
  [m,n] : matrix_size(obj),
  qm : diagmatrix(m,1),
  rm : copymatrix(obj),
  for i: 1 thru (if m=n then n-1 else n) do
  block([x,h],
    x : getSubmatrix(rm,i,i,m,i),
    h : diagmatrix(m,1),
    setSubmatrix(h,i,i,householder(x)),
    qm : qm . h,
    rm : h . rm),
  [qm,rm]);
solveUpperTriangular(r,b) := block([n,x,index,k],
  n : second(matrix_size(r)),
  x : genmatrix(lambda([a, b], 0), n, 1),
  for k: n thru 1 step -1 do
  (index : min(n,k+1),
    x[k,1] : (b[k,1] - (getSubmatrix(r,k,index,k,n) . getSubmatrix(x,index,1,n,1)))/r[k,k]),
  x);
lsqr(a,b) := block([q,r,n],
  [q,r] : qr(a),
  n : second(matrix_size(r)),
  solveUpperTriangular(getSubmatrix(r,1,1,n,n), transpose(q) . b));
polyfit(x,y,n) := block([a,j],
  a : genmatrix(lambda([i,j], if j=1 then 1.0b0 else bfloat(x[i,1]^(j-1))),
    length(x),n+1),
  lsqr(a,y));

Then we have the examples:

(%i) [q,r] : qr(a);

                 [   6   69     58   ]
                 [ - -   ---    ---  ]
                 [   7   175    175  ]
                 [                   ]  [ - 14  - 21    14  ]
                 [   3    158     6  ]  [                   ]
(%o)            [[ - -  - ---  - --- ], [  0    - 175   70  ]]
                 [   7    175    175 ]  [                   ]
                 [                   ]  [  0      0    - 35 ]
                 [  2     6     33   ]
                 [  -   - --    --   ]
                 [  7     35    35   ]
(%i) mat_norm(q . r - a, 1);

(%o)                                   0
(%i) x : transpose(matrix([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]))$

(%i) y : transpose(matrix([1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]))$

(%i) fpprec : 30$

(%i) polyfit(x, y, 2);

                    [ 9.99999999999999999999999999996b-1 ]
                    [                                    ]
(%o)                [ 2.00000000000000000000000000002b0  ]
                    [                                    ]
                    [               3.0b0                ]

Nim

Translation of: Python

Library: arraymancer

The library “arraymancer” provides a function “qr” to get the QR decomposition. Using the Tensor type of “arraymancer” we propose here an implementation of this decomposition adapted from the Python version.

import math, strformat, strutils
import arraymancer

####################################################################################################
# First part: QR decomposition.

proc eye(n: Positive): Tensor[float] =
  ## Return the (n, n) identity matrix.
  result = newTensor[float](n.int, n.int)
  for i in 0..<n: result[i, i] = 1

proc norm(v: Tensor[float]): float =
  ## return the norm of a vector.
  assert v.shape.len == 1
  result = sqrt(dot(v, v)) * sgn(v[0]).toFloat

proc houseHolder(a: Tensor[float]): Tensor[float] =
  ## return the house holder of vector "a".
  var v = a / (a[0] + norm(a))
  v[0] = 1
  result = eye(a.shape[0]) - (2 / dot(v, v)) * (v.unsqueeze(1) * v.unsqueeze(0))

proc qrDecomposition(a: Tensor): tuple[q, r: Tensor] =
  ## Return the QR decomposition of matrix "a".
  assert a.shape.len == 2
  let m = a.shape[0]
  let n = a.shape[1]
  result.q = eye(m)
  result.r = a.clone
  for i in 0..<(n - ord(m == n)):
    var h = eye(m)
    h[i..^1, i..^1] = houseHolder(result.r[i..^1, i].squeeze(1))
    result.q = result.q * h
    result.r = h * result.r

####################################################################################################
# Second part: polynomial regression example.

proc lsqr(a, b: Tensor[float]): Tensor[float] =
  let (q, r) = a.qrDecomposition()
  let n = r.shape[1]
  result = solve(r[0..<n, _], (q.transpose() * b)[0..<n])

proc polyfit(x, y: Tensor[float]; n: int): Tensor[float] =
  var z = newTensor[float](x.shape[0], n + 1)
  var t = x.reshape(x.shape[0], 1)
  for i in 0..n: z[_, i] = t^.i.toFloat
  result = lsqr(z, y.transpose())

#———————————————————————————————————————————————————————————————————————————————————————————————————

proc printMatrix(a: Tensor) =
  var str: string
  for i in 0..<a.shape[0]:
    let start = str.len
    for j in 0..<a.shape[1]:
      str.addSep(" ", start)
      str.add &"{a[i, j]:8.3f}"
    str.add '\n'
  stdout.write str

proc printVector(a: Tensor) =
  var str: string
  for i in 0..<a.shape[0]:
    str.addSep(" ")
    str.add &"{a[i]:4.1f}"
  echo str


let mat = [[12, -51,   4],
           [ 6, 167, -68],
           [-4,  24, -41]].toTensor.astype(float)

let (q, r) = mat.qrDecomposition()
echo "Q:"
printMatrix q
echo "R:"
printMatrix r
echo()

let x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10].toTensor.astype(float)
let y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321].toTensor.astype(float)
echo "polyfit:"
printVector polyfit(x, y, 2)

Output:

Q:
  -0.857    0.394    0.331
  -0.429   -0.903   -0.034
   0.286   -0.171    0.943
R:
 -14.000  -21.000   14.000
   0.000 -175.000   70.000
   0.000   -0.000  -35.000

polyfit:
 1.0  2.0  3.0

PARI/GP

Works with: PARI/GP version 2.6.0 and above

matqr(M)

Perl

Letting the PDL module do all the work.

use strict;
use warnings;

use PDL;
use PDL::LinearAlgebra qw(mqr);

my $a = pdl(
      [12, -51,   4],
      [ 6, 167, -68],
      [-4,  24, -41],
      [-1,   1,   0],
      [ 2,   0,   3]
);

my ($q, $r) = mqr($a);
print $q, $r, $q x $r;

Output:

[
 [ -0.84641474   0.39129081  -0.34312406]
 [ -0.42320737  -0.90408727  0.029270162]
 [  0.28213825  -0.17042055  -0.93285599]
 [ 0.070534562 -0.014040652  0.001099372]
 [ -0.14106912  0.016655511   0.10577161]
]

[
 [-14.177447 -20.666627  13.401567]
 [         0 -175.04254  70.080307]
 [         0          0  35.201543]
]

[
 [           12           -51             4]
 [            6           167           -68]
 [           -4            24           -41]
 [           -1             1             0]
 [            2             0             3]
]

Phix

using matrix_mul() from Matrix_multiplication#Phix and matrix_transpose() from Matrix_transposition#Phix

-- demo/rosettacode/QRdecomposition.exw
with javascript_semantics

function matrix_mul(sequence a, b)
    integer arows = ~a, acols = ~a[1],
            brows = ~b, bcols = ~b[1]
    if acols!=brows then return 0 end if
    sequence c = repeat(repeat(0,bcols),arows)
    for i=1 to arows do
        for j=1 to bcols do
            for k=1 to acols do
                c[i][j] += a[i][k]*b[k][j]
            end for
        end for
    end for
    return c
end function

function vtranspose(sequence v)
-- transpose a vector of length m into an mx1 matrix, 
--                       eg {1,2,3} -> {{1},{2},{3}}
    integer l = length(v)
    sequence res = repeat(0,l)
    for i=1 to l do res[i] = {v[i]} end for
    return res
end function

function mat_col(sequence a, integer col)
    integer la = length(a)
    sequence res = repeat(0,la)
    for i=col to la do
        res[i] = a[i,col]
    end for
    return res
end function

function mat_norm(sequence a)
    atom res = 0
    for i=1 to length(a) do
        res += a[i]*a[i]
    end for
    res = sqrt(res)
    return res
end function

function mat_ident(integer n)
    sequence res = repeat(repeat(0,n),n)
    for i=1 to n do
        res[i,i] = 1
    end for
    return res
end function

function QRHouseholder(sequence a)
    integer cols = length(a[1]),
            rows = length(a),
            m = max(cols,rows),
            n = min(rows,cols)
    sequence q, I = mat_ident(m), Q = I, u, v
--
-- Programming note: The code of this main loop was not as easily
-- written as the first glance might suggest. Explicitly setting 
-- to 0 any a[i,j] [etc] that should be 0 but have inadvertently 
-- gotten set to +/-1e-15 or thereabouts may be advisable. The
-- commented-out code was retrieved from a backup and should be
-- treated as an example and not be trusted (iirc, it made no
-- difference to the test cases used, so I deleted it, and then
-- had second thoughts about it a few days later).
--
    for j=1 to min(m-1,n) do
        u = mat_col(a,j)
        u[j] -= mat_norm(u)
        v = sq_div(u,mat_norm(u))
        q = sq_sub(I,sq_mul(2,matrix_mul(vtranspose(v),{v})))
        a = matrix_mul(q,a)
--      for row=j+1 to length(a) do
--          a[row][j] = 0
--      end for
        Q = matrix_mul(Q,q)
    end for

    -- Get the upper triangular matrix R.
    sequence R = repeat(repeat(0,n),m)
    for i=1 to n do -- (logically 1 to m(>=n), but no need)
        for j=i to n do
            R[i,j] = a[i,j]
        end for
    end for
        
    return {Q,R}
end function

constant a = {{12, -51,   4},
              { 6, 167, -68},
              {-4,  24, -41}}

sequence {q,r} = QRHouseholder(a)

ppOpt({pp_Nest,1,pp_IntFmt,"%4d",pp_FltFmt,"%4g",pp_IntCh,false})
?"A"        pp(a)
?"Q"        pp(q)
?"R"        pp(r)
?"Q * R"    pp(matrix_mul(q,r))

function matrix_transpose(sequence mat)
    integer rows = length(mat),
            cols = length(mat[1])
    sequence res = repeat(repeat(0,rows),cols)
    for r=1 to rows do
        for c=1 to cols do
            res[c][r] = mat[r][c]
        end for
    end for
    return res
end function

--?"Q * Q'" pp(matrix_mul(q,matrix_transpose(q))) -- (~1e-16s)
?"Q * Q`"   pp(sq_round(matrix_mul(q,matrix_transpose(q)),1e15))

procedure least_squares()
    sequence x = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10},
             y = {1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321},
             a = repeat(repeat(0,3),length(x))
    for i=1 to length(x) do
        for j=1 to 3 do
            a[i,j] = power(x[i],j-1)
        end for
    end for
    {q,r} = QRHouseholder(a)
    sequence t = matrix_transpose(q),
             b = matrix_mul(t,vtranspose(y)),
             z = repeat(0,3)
    for k=3 to 1 by -1 do
        atom s = 0
        if k<3 then
            for j = k+1 to 3 do
                s += r[k,j]*z[j]
            end for
        end if
        z[k] = (b[k][1]-s)/r[k,k]
    end for
    printf(1,"Least-squares solution:\n")
--  printf(1," %v\n",{z})                   -- {1.0,2.0.3,0}
--  printf(1," %v\n",{sq_sub(z,{1,2,3})})   -- (+/- ~1e-14s)
    printf(1," %v\n",{sq_round(z,1e13)})    -- {1,2,3}
end procedure
least_squares()

Output:

"A"
{{  12, -51,   4},
 {   6, 167, -68},
 {  -4,  24, -41}}
"Q"
{{0.85714,-0.3943,0.33143},
 {0.42857,0.90286,-0.0343},
 {-0.2857,0.17143,0.94286}}
"R"
{{  14,  21, -14},
 {   0, 175, -70},
 {   0,   0, -35}}
"Q * R"
{{  12, -51,   4},
 {   6, 167, -68},
 {  -4,  24, -41}}
"Q * Q`"
{{   1,   0,   0},
 {   0,   1,   0},
 {   0,   0,   1}}
Least-squares solution:
 {1,2,3}

PowerShell

function qr([double[][]]$A) {
    $m,$n = $A.count, $A[0].count
    $pm,$pn = ($m-1), ($n-1)
    [double[][]]$Q = 0..($m-1) | foreach{$row = @(0) * $m; $row[$_] = 1; ,$row} 
    [double[][]]$R = $A | foreach{$row = $_; ,@(0..$pn | foreach{$row[$_]})}
    foreach ($h in 0..$pn) { 
        [double[]]$u = $R[$h..$pm] | foreach{$_[$h]} 
        [double]$nu = $u | foreach {[double]$sq = 0} {$sq += $_*$_} {[Math]::Sqrt($sq)} 
        $u[0] -= if ($u[0] -lt 0) {$nu} else {-$nu}
        [double]$nu = $u | foreach {$sq = 0} {$sq += $_*$_} {[Math]::Sqrt($sq)} 
        [double[]]$u = $u | foreach { $_/$nu}
        [double[][]]$v = 0..($u.Count - 1) | foreach{$i = $_; ,($u | foreach{2*$u[$i]*$_})}
        [double[][]]$CR = $R | foreach{$row = $_; ,@(0..$pn | foreach{$row[$_]})}
        [double[][]]$CQ = $Q | foreach{$row = $_; ,@(0..$pm | foreach{$row[$_]})}
        foreach ($i in  $h..$pm) {
            foreach ($j in  $h..$pn) {
                $R[$i][$j] -=  $h..$pm | foreach {[double]$sum = 0} {$sum += $v[$i-$h][$_-$h]*$CR[$_][$j]} {$sum}
            }
        }
        if (0 -eq $h)  {
            foreach ($i in  $h..$pm) {
                foreach ($j in  $h..$pm) {
                    $Q[$i][$j] -=  $h..$pm | foreach {$sum = 0} {$sum += $v[$i][$_]*$CQ[$_][$j]} {$sum}
                }
            }
        } else  {
            $p = $h-1
            foreach ($i in  $h..$pm) {
                foreach ($j in  0..$p) {
                    $Q[$i][$j] -=  $h..$pm | foreach {$sum = 0} {$sum += $v[$i-$h][$_-$h]*$CQ[$_][$j]} {$sum}
                }
                foreach ($j in  $h..$pm) {
                    $Q[$i][$j] -=  $h..$pm | foreach {$sum = 0} {$sum += $v[$i-$h][$_-$h]*$CQ[$_][$j]} {$sum}
                }
            }
        }
    }
    foreach ($i in  0..$pm) {
        foreach ($j in  $i..$pm) {$Q[$i][$j],$Q[$j][$i] = $Q[$j][$i],$Q[$i][$j]}
    }
    [PSCustomObject]@{"Q" = $Q; "R" = $R}
}

function leastsquares([Double[][]]$A,[Double[]]$y) {
    $QR = qr $A
    [Double[][]]$Q = $QR.Q
    [Double[][]]$R = $QR.R
    $m,$n = $A.count, $A[0].count
    [Double[]]$z = foreach ($j in  0..($m-1)) { 
            0..($m-1) | foreach {$sum = 0} {$sum += $Q[$_][$j]*$y[$_]} {$sum}
    }
    [Double[]]$x = @(0)*$n
    for ($i = $n-1; $i -ge 0; $i--) {
        for ($j = $i+1; $j -lt $n; $j++) {
            $z[$i] -= $x[$j]*$R[$i][$j]
        }
        $x[$i] = $z[$i]/$R[$i][$i]
    }
    $x
}

function polyfit([Double[]]$x,[Double[]]$y,$n) {
    $m = $x.Count 
    [Double[][]]$A = 0..($m-1) | foreach{$row = @(1) * ($n+1); ,$row} 
    for ($i = 0; $i -lt $m; $i++) {
        for ($j = $n-1; 0 -le $j; $j--) {
            $A[$i][$j] = $A[$i][$j+1]*$x[$i]
        }
    }
    leastsquares $A $y
}

function show($m) {$m | foreach {write-host "$_"}}

$A = @(@(12,-51,4), @(6,167,-68), @(-4,24,-41))
$x = @(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
$y = @(1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321)
$QR = qr $A
$ps = (polyfit $x $y 2)
"Q = "
show $QR.Q
"R = "
show $QR.R
"polyfit "
"X^2 X constant"
"$(polyfit $x $y 2)"

Output:

Q = 
-0.857142857142857 0.394285714285714 -0.331428571428571
-0.428571428571429 -0.902857142857143 0.0342857142857143
0.285714285714286 -0.171428571428571 -0.942857142857143
R = 
-14 -21 14
8.88178419700125E-16 -175 70
-4.44089209850063E-16 0 35
polyfit 
X^2 X constant
3 1.99999999999998 1.00000000000005

Python

Library: NumPy

Numpy has a qr function but here is a reimplementation to show construction and use of the Householder reflections.

#!/usr/bin/env python3

import numpy as np

def qr(A):
    m, n = A.shape
    Q = np.eye(m)
    for i in range(n - (m == n)):
        H = np.eye(m)
        H[i:, i:] = make_householder(A[i:, i])
        Q = np.dot(Q, H)
        A = np.dot(H, A)
    return Q, A

def make_householder(a):
    v = a / (a[0] + np.copysign(np.linalg.norm(a), a[0]))
    v[0] = 1
    H = np.eye(a.shape[0])
    H -= (2 / np.dot(v, v)) * np.dot(v[:, None], v[None, :])
    return H

# task 1: show qr decomp of wp example
a = np.array(((
    (12, -51,   4),
    ( 6, 167, -68),
    (-4,  24, -41),
)))

q, r = qr(a)
print('q:\n', q.round(6))
print('r:\n', r.round(6))

# task 2: use qr decomp for polynomial regression example
def polyfit(x, y, n):
    return lsqr(x[:, None]**np.arange(n + 1), y.T)

def lsqr(a, b):
    q, r = qr(a)
    _, n = r.shape
    return np.linalg.solve(r[:n, :], np.dot(q.T, b)[:n])

x = np.array((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10))
y = np.array((1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321))

print('\npolyfit:\n', polyfit(x, y, 2))

Output:

q:
 [[-0.857143  0.394286  0.331429]
 [-0.428571 -0.902857 -0.034286]
 [ 0.285714 -0.171429  0.942857]]
r:
 [[ -14.  -21.   14.]
 [   0. -175.   70.]
 [   0.    0.  -35.]]

polyfit:
 [ 1.  2.  3.]

R

# R has QR decomposition built-in (using LAPACK or LINPACK)

a <- matrix(c(12, -51, 4, 6, 167, -68, -4, 24, -41), nrow=3, ncol=3, byrow=T)
d <- qr(a)
qr.Q(d)
qr.R(d)

# now fitting a polynomial
x <- 0:10
y <- 3*x^2 + 2*x + 1

# using QR decomposition directly
a <- cbind(1, x, x^2)
qr.coef(qr(a), y)

# using least squares
a <- cbind(x, x^2)
lsfit(a, y)$coefficients

# using a linear model
xx <- x*x
m <- lm(y ~ x + xx)
coef(m)

Racket

Racket has QR-decomposition builtin:

> (require math)
> (matrix-qr (matrix [[12 -51   4]
                      [ 6 167 -68]
                      [-4  24 -41]]))
(array #[#[6/7 -69/175 -58/175] #[3/7 158/175 6/175] #[-2/7 6/35 -33/35]])
(array #[#[14 21 -14] #[0 175 -70] #[0 0 35]])

The builtin QR-decomposition uses the Gram-Schmidt algorithm.

Here is an implementation of the Householder method:

#lang racket
(require math/matrix math/array)
(define-values (T I col size)
  (values ; short names
   matrix-transpose identity-matrix matrix-col matrix-num-rows))

(define (scale c A) (matrix-scale A c))
(define (unit n i) (build-matrix n 1 (λ (j _) (if (= j i) 1 0))))

(define (H u)
  (matrix- (I (size u))
           (scale (/ 2 (matrix-dot u u))
                  (matrix* u (T u)))))

(define (normal a)
  (define a0 (matrix-ref a 0 0))
  (matrix- a (scale (* (sgn a0) (matrix-2norm a))
                    (unit (size a) 0))))

(define (QR A)
  (define n (size A))
  (for/fold ([Q (I n)] [R A]) ([i (- n 1)])
    (define Hi (H (normal (submatrix R (:: i n) (:: i (+ i 1))))))
    (define Hi* (if (= i 0) Hi (block-diagonal-matrix (list (I i) Hi)))) 
    (values (matrix* Q Hi*) (matrix* Hi* R))))

(QR (matrix [[12 -51   4]
               [ 6 167 -68]
               [-4  24 -41]]))

Output:

(array #[#[6/7 69/175 -58/175] 
             #[3/7 -158/175 6/175] 
             #[-2/7 -6/35 -33/35]])
(array #[#[14 21 -14] 
             #[0 -175 70] 
             #[0 0 35]])

Raku

(formerly Perl 6)

Works with: rakudo version 2018.06

# sub householder translated from https://codereview.stackexchange.com/questions/120978/householder-transformation

use v6;

sub identity(Int:D $m --> Array of Array) {
   my Array @M; 
   
   for 0 ..^ $m -> $i {
      @M.push: [0 xx $m];
      @M[$i; $i] = 1;
   }

   @M;
}

multi multiply(Array:D @A, @b where Array:D --> Array) {
   my @c;

   for ^@A X ^@b -> ($i, $j) {
      @c[$i] += @A[$i; $j] * @b[$j];
   }

   @c;
}

multi multiply(Array:D @A, Array:D @B --> Array of Array) {
   my Array @C;

   for ^@A X ^@B[0] -> ($i, $j) {
      @C[$i; $j] += @A[$i; $_] * @B[$_; $j] for ^@B;
   }

   @C;
}

sub transpose(Array:D @M --> Array of Array) {
   my ($rows, $cols) = (@M.elems, @M[0].elems);

   my Array @T;

   for ^$cols X ^$rows -> ($j, $i) {
      @T[$j; $i] = @M[$i; $j];
   }

   @T;
}

####################################################
# NOTE: @A gets overwritten and becomes @R, only need
# to return @Q.
####################################################
sub householder(Array:D @A --> Array) {
   my Int ($m, $n) = (@A.elems, @A[0].elems);
   my @v = 0 xx $m;
   my Array @Q = identity($m);

   for 0 ..^ $n -> $k {
      my Real $sum = 0;
      my Real $A0 = @A[$k; $k];
      my Int $sign = $A0 < 0 ?? -1 !! 1;

      for $k ..^ $m -> $i {
         $sum += @A[$i; $k] * @A[$i; $k];
      }

      my Real $sqr_sum = $sign * sqrt($sum);
      my Real $tmp = sqrt(2 * ($sum + $A0 * $sqr_sum));     
      @v[$k] = ($sqr_sum  + $A0) / $tmp;

      for ($k + 1) ..^ $m -> $i {
         @v[$i] = @A[$i; $k] / $tmp;
      }              

      for 0 ..^ $n -> $j {
         $sum = 0;

         for $k ..^ $m -> $i {
            $sum += @v[$i] * @A[$i; $j];
         }

         for $k ..^ $m -> $i {
            @A[$i; $j] -= 2 * @v[$i] * $sum;
         }
      }       

      for 0 ..^ $m -> $j {
         $sum = 0;

         for $k ..^ $m -> $i {
            $sum += @v[$i] * @Q[$i; $j];
         }

         for $k ..^ $m -> $i {
            @Q[$i; $j] -= 2 * @v[$i] * $sum;
         }
      }
   }

   @Q
}
 
sub dotp(@a where Array:D, @b where Array:D --> Real) {
   [+] @a >>*<< @b;
}

sub upper-solve(Array:D @U, @b where Array:D, Int:D $n --> Array) {
   my @y = 0 xx $n;

   @y[$n - 1] = @b[$n - 1] / @U[$n - 1; $n - 1];

   for reverse ^($n - 1) -> $i {
      @y[$i] = (@b[$i] - (dotp(@U[$i], @y))) / @U[$i; $i];
   }
   
   @y;
}

sub polyfit(@x where Array:D, @y where Array:D, Int:D $n) {
   my Int $m = @x.elems;
   my Array @V;

   # Vandermonde matrix
   for ^$m X (0 .. $n) -> ($i, $j) {
      @V[$i; $j] = @x[$i] ** $j
   }

   # least squares
   my $Q = householder(@V);
   my @b = multiply($Q, @y);

   return upper-solve(@V, @b, $n + 1);
}

sub print-mat(Array:D @M, Str:D $name) {
   my Int ($m, $n) = (@M.elems, @M[0].elems);
   print "\n$name:\n";

   for 0 ..^ $m -> $i {
      for 0 ..^ $n -> $j {
         print @M[$i; $j].fmt("%12.6f ");
      }

      print "\n";
   }
}

sub MAIN() {
   ############
   # 1st part #
   ############
   my Array @A = (
      [12, -51,   4],
      [ 6, 167, -68],
      [-4,  24, -41],
      [-1,   1,   0],
      [ 2,   0,   3]
   );

   print-mat(@A, 'A');
   my $Q = householder(@A);
   $Q = transpose($Q);
   print-mat($Q, 'Q');
   # after householder, @A is now @R
   print-mat(@A, 'R');
   print-mat(multiply($Q, @A), 'check Q x R = A');

   ############
   # 2nd part #
   ############
   my @x = [^11];
   my @y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321];

   my @coef = polyfit(@x, @y, 2);

   say 
      "\npolyfit:\n", 
      <constant X X^2>.fmt("%12s"),
      "\n",
      @coef.fmt("%12.6f");
}

output:

A:
   12.000000   -51.000000     4.000000 
    6.000000   167.000000   -68.000000 
   -4.000000    24.000000   -41.000000 
   -1.000000     1.000000     0.000000 
    2.000000     0.000000     3.000000 

Q:
   -0.846415     0.391291    -0.343124     0.066137    -0.091462 
   -0.423207    -0.904087     0.029270     0.017379    -0.048610 
    0.282138    -0.170421    -0.932856    -0.021942     0.143712 
    0.070535    -0.014041     0.001099     0.997401     0.004295 
   -0.141069     0.016656     0.105772     0.005856     0.984175 

R:
  -14.177447   -20.666627    13.401567 
   -0.000000  -175.042539    70.080307 
    0.000000     0.000000    35.201543 
   -0.000000     0.000000     0.000000 
    0.000000    -0.000000     0.000000 

check Q x R = A:
   12.000000   -51.000000     4.000000 
    6.000000   167.000000   -68.000000 
   -4.000000    24.000000   -41.000000 
   -1.000000     1.000000    -0.000000 
    2.000000    -0.000000     3.000000 

polyfit:
    constant            X          X^2
    1.000000     2.000000     3.000000

Rascal

This function applies the Gram Schmidt algorithm. Q is printed in the console, R can be printed or visualized.

import util::Math;
import Prelude;
import vis::Figure;
import vis::Render;

public rel[real,real,real] QRdecomposition(rel[real x, real y, real v] matrix){
	//orthogonalcolumns
	oc = domainR(matrix, {0.0});
	for (x <- sort(toList(domain(matrix)-{0.0}))){
		c = domainR(matrix, {x});
		o = domainR(oc, {x-1});
		
		for (n <- [1.0 .. x]){
			o = domainR(oc, {n-1});
			c = matrixSubtract(c, matrixMultiplybyN(o, matrixDotproduct(o, c)/matrixDotproduct(o, o)));
			}
			
		oc += c;
	}
	
	Q = {};
	//from orthogonal to orthonormal columns
	for (el <- oc){
		c = domainR(oc, {el[0]});
		Q += matrixNormalize({el}, c);
	}
	
	//from Q to R
	R= matrixMultiplication(matrixTranspose(Q), matrix);
	R= {<x,y,toReal(round(v))> | <x,y,v> <- R};
	
	println("Q:");
	iprintlnExp(Q);
	println();
	println("R:");
	return R;
}

//a function that takes the transpose of a matrix, see also Rosetta Code problem "Matrix transposition"
public rel[real, real, real] matrixTranspose(rel[real x, real y, real v] matrix){
	return {<y, x, v> | <x, y, v> <- matrix};
}

//a function to normalize an element of a matrix by the normalization of a column
public rel[real,real,real] matrixNormalize(rel[real x, real y, real v] element, rel[real x, real y, real v] column){
	normalized = 1.0/nroot((0.0 | it + v*v | <x,y,v> <- column), 2);
	return matrixMultiplybyN(element, normalized);
}

//a function that takes the dot product, see also Rosetta Code problem "Dot product"
public real matrixDotproduct(rel[real x, real y, real v] column1, rel[real x, real y, real v] column2){
	return (0.0 | it + v1*v2 | <x1,y1,v1> <- column1, <x2,y2,v2> <- column2, y1==y2);
}

//a function to subtract two columns
public rel[real,real,real] matrixSubtract(rel[real x, real y, real v] column1, rel[real x, real y, real v] column2){
	return {<x1,y1,v1-v2> | <x1,y1,v1> <- column1, <x2,y2,v2> <- column2, y1==y2};
}

//a function to multiply a column by a number
public rel[real,real,real] matrixMultiplybyN(rel[real x, real y, real v] column, real n){
	return {<x,y,v*n> | <x,y,v> <- column};
}

//a function to perform matrix multiplication, see also Rosetta Code problem "Matrix multiplication".
public rel[real, real, real] matrixMultiplication(rel[real x, real y, real v] matrix1, rel[real x, real y, real v] matrix2){
	if (max(matrix1.x) == max(matrix2.y)){
		p = {<x1,y1,x2,y2, v1*v2> | <x1,y1,v1> <- matrix1, <x2,y2,v2> <- matrix2};

		result = {};
		for (y <- matrix1.y){
			for (x <- matrix2.x){
				v = (0.0 | it + v | <x1, y1, x2, y2, v> <- p,  x==x2 && y==y1, x1==y2 && y2==x1);
				result += <x,y,v>;
			}
		}
		return result;
	}
	else throw "Matrix sizes do not match.";
} 	

// a function to visualize the result
public void displayMatrix(rel[real x, real y, real v] matrix){
	points = [box(text("<v>"), align(0.3333*(x+1),0.3333*(y+1)),shrink(0.25)) | <x,y,v> <- matrix];
	render(overlay([*points], aspectRatio(1.0)));
}

//a matrix, given by a relation of <x-coordinate, y-coordinate, value>.
public rel[real x, real y, real v] matrixA = {
<0.0,0.0,12.0>, <0.0,1.0, 6.0>, <0.0,2.0,-4.0>, 
<1.0,0.0,-51.0>, <1.0,1.0,167.0>, <1.0,2.0,24.0>, 
<2.0,0.0,4.0>, <2.0,1.0,-68.0>, <2.0,2.0,-41.0>
};

Example using visualization

rascal>displayMatrix(QRdecomposition(matrixA))

Q:
{
  <1.0,0.0,-0.394285714285714285714285714285714285714285714285714285714285714285713300>,
  <0.0,0.0,0.857142857142857142857142857142857142857142857142857142857142857142840>,
  <0.0,1.0,0.428571428571428571428571428571428571428571428571428571428571428571420>,
  <0.0,2.0,-0.285714285714285714285714285714285714285714285714285714285714285714280>,
  <2.0,0.0,-0.33142857142857142857142857142857142857142857142857142857142857142858800>,
  <1.0,2.0,0.171428571428571428571428571428571428571428571428571428571428571428571000>,
  <2.0,2.0,-0.94285714285714285714285714285714285714285714285714285714285714285719000>,
  <1.0,1.0,0.902857142857142857142857142857142857142857142857142857142857142857140600>,
  <2.0,1.0,0.03428571428571428571428571428571428571428571428571428571428571428571600>
}
See R in picture

SAS

/* See http://support.sas.com/documentation/cdl/en/imlug/63541/HTML/default/viewer.htm#imlug_langref_sect229.htm */

proc iml;
a={12 -51 4,6 167 -68,-4 24 -41};
print(a);
call qr(q,r,p,d,a);
print(q);
print(r);
quit;

/*
                  a

           12       -51         4
            6       167       -68
           -4        24       -41


                  q

    -0.857143 0.3942857 -0.331429
    -0.428571 -0.902857 0.0342857
    0.2857143 -0.171429 -0.942857


                  r

          -14       -21        14
            0      -175        70
            0         0        35
*/

Scala

Output:

Best seen running in your browser Scastie (remote JVM).

import java.io.{PrintWriter, StringWriter}

import Jama.{Matrix, QRDecomposition}

object QRDecomposition extends App {
  val matrix =
    new Matrix(
      Array[Array[Double]](Array(12, -51, 4),
        Array(6, 167, -68),
        Array(-4, 24, -41)))
  val d = new QRDecomposition(matrix)

  def toString(m: Matrix): String = {
    val sw = new StringWriter
    m.print(new PrintWriter(sw, true), 8, 6)
    sw.toString
  }

  print(toString(d.getQ))
  print(toString(d.getR))

}

SequenceL

Translation of: Go

import <Utilities/Math.sl>;
import <Utilities/Sequence.sl>;
import <Utilities/Conversion.sl>;

main :=
    let
        qrTest := [[12.0, -51.0,   4.0],
                   [ 6.0, 167.0, -68.0],
                   [-4.0,  24.0, -41.0]];
        
        qrResult := qr(qrTest);
        
        x := 1.0*(0 ... 10);
        y := 1.0*[1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321];
        
        regResult := polyfit(x, y, 2);
    in
        "q:\n" ++ delimit(delimit(floatToString(qrResult[1], 6), ','), '\n') ++ "\n\n" ++ 
        "r:\n" ++ delimit(delimit(floatToString(qrResult[2], 1), ','), '\n') ++ "\n\n" ++
        "polyfit:\n" ++ "[" ++ delimit(floatToString(regResult, 1), ',') ++ "]";

//---Polynomial Regression---

polyfit(x(1), y(1), n) :=
    let
        a[j] := x ^ j foreach j within 0 ... n;
    in  
        lsqr(transpose(a), transpose([y]));
    
lsqr(a(2), b(2)) :=
    let
        qrDecomp := qr(a);
        prod := mm(transpose(qrDecomp[1]), b);
    in
        solveUT(qrDecomp[2], prod);
        
solveUT(r(2), b(2)) := 
    let 
        n := size(r[1]);
    in
        solveUTHelper(r, b, n, duplicate(0.0, n)); 

solveUTHelper(r(2), b(2), k, x(1)) :=
    let
        n := size(r[1]);
        newX :=  setElementAt(x, k, (b[k][1] - sum(r[k][(k+1) ... n] * x[(k+1) ... n])) / r[k][k]);
    in
        x when k <= 0
    else
        solveUTHelper(r, b, k - 1, newX);

//---QR Decomposition---

qr(A(2)) := qrHelper(A, id(size(A)), 1);

qrHelper(A(2), Q(2), i) :=
    let
        m := size(A);
        n := size(A[1]);
        
        householder := makeHouseholder(A[i ... m, i]);
        
        H[j,k] := 
                householder[j - i + 1][k - i + 1] when j >= i and k >= i 
            else
                1.0 when j = k else 0.0
            foreach j within 1 ... m,
                    k within 1 ... m;
    in
        [Q,A] when i > (n - 1 when m = n else n)
    else
        qrHelper(mm(H, A), mm(Q, H), i + 1);
    

makeHouseholder(a(1)) := 
    let
        v := [1.0] ++ tail(a / (a[1] + sqrt(sum(a ^ 2)) * sign(a[1])));
        
        H := id(size(a)) - (2.0 / mm([v], transpose([v])))[1,1] * mm(transpose([v]), [v]);
    in
        H;

//---Utilities---

id(n)[i,j] := 1.0 when i = j else 0.0
              foreach i within 1 ... n,
                      j within 1 ... n;
                        
mm(A(2), B(2))[i,j] := sum( A[i] * transpose(B)[j] );

Output:

"q:
-0.857143,0.394286,0.331429
-0.428571,-0.902857,-0.034286
0.285714,-0.171429,0.942857

r:
-14.0,-21.0,14.0
-0.0,-175.0,70.0
0.0,0.0,-35.0

polyfit:
[1.0,2.0,3.0]"

SPAD

See QR_decomposition#Axiom in Axiom.

Standard ML

Translation of: Axiom

We first define a signature for a radical category joined with a field. We then define a functor with (a) structures to define operators and functions for Array and Array2, and (b) functions for the QR decomposition:

signature RADCATFIELD = sig
type real
val zero : real
val one : real
val + : real * real -> real
val - : real * real -> real
val * : real * real -> real
val / : real * real -> real
val sign : real -> real
val sqrt : real -> real
end
		      
functor QR(F: RADCATFIELD) = struct
structure A = struct
local
    open Array
in
fun unitVector n = tabulate (n, fn i => if i=0 then F.one else F.zero)
fun map f x = tabulate(length x, fn i => f(sub(x,i)))
fun map2 f (x, y) = tabulate(length x, fn i => f(sub(x,i),sub(y,i)))
val op + = map2 F.+
val op - = map2 F.-
val op * = map2 F.*
fun multc(c,x) = array(length x,c)*x
fun dot (x,y) = foldl F.+ F.zero (x*y)
fun outer f (x,y) =
    Array2.tabulate Array2.RowMajor (length x, length y,
				     fn (i,j) => f(sub(x,i),sub(y,j)))
fun copy x = map (fn x => x) x
fun fromVector v = tabulate(Vector.length v, fn i => Vector.sub(v,i))
fun slice(x,i,sz) =
    let	open ArraySlice
	val s = slice(x,i,sz)
    in Array.tabulate(length s, fn i => sub(s,i)) end
end
end
structure M = struct
local
    open Array2
in
fun map f x = tabulate RowMajor (nRows x, nCols x, fn (i,j) => f(sub(x,i,j)))
fun map2 f (x, y) =
    tabulate RowMajor (nRows x, nCols x, fn (i,j) => f(sub(x,i,j),sub(y,i,j)))
fun scalarMatrix(m, x) = tabulate RowMajor (m,m,fn (i,j) => if i=j then x else F.zero)
fun multc(c, x) = map (fn xij => F.*(c,xij)) x
val op + = map2 F.+
val op - = map2 F.-
fun column(x,i) = A.fromVector(Array2.column(x,i))
fun row(x,i) = A.fromVector(Array2.row(x,i))
fun x*y = tabulate RowMajor (nRows x, nCols y,
			     fn (i,j) => A.dot(row(x,i), column(y,j)))
fun multa(x,a) = Array.tabulate (nRows x, fn i => A.dot(row(x,i), a))
fun copy x = map (fn x => x) x
fun subMatrix(h, i1, i2, j1, j2) =
    tabulate RowMajor (Int.+(Int.-(i2,i1),1),
		       Int.+(Int.-(j2,j1),1),
		       fn (a,b) => sub(h,Int.+(i1,a),Int.+(j1,b)))
fun transpose m = tabulate RowMajor (nCols m,
				     nRows m,
				     fn (i,j) => sub(m,j,i))
fun updateSubMatrix(h,i,j,s) =
    tabulate RowMajor (nRows s, nCols s, fn (a,b) => update(h,Int.+(i,a),Int.+(j,b),sub(s,a,b)))
end
end
fun toList a =
    List.tabulate(Array2.nRows a, fn i => List.tabulate(Array2.nCols a, fn j => Array2.sub(a,i,j)))
fun householder a =
    let open Array
	val m = length a
	val len = F.sqrt(A.dot(a,a))
	val u = A.+(a, A.multc(F.*(len,F.sign(sub(a,0))), A.unitVector m))
	val v = A.multc(F./(F.one,sub(u,0)), u)
	val beta = F./(F.+(F.one,F.one),A.dot(v,v))
    in
	M.-(M.scalarMatrix(m,F.one), M.multc(beta,A.outer F.* (v,v)))
    end
fun qr mat =
    let open Array2
	val (m,n) = dimensions mat
	val upperIndex = if m=n then Int.-(n,1) else n
	fun loop(i,qm,rm) = if i=upperIndex then {q=qm,r=rm} else 
			    let val x = A.slice(A.fromVector(column(rm,i)),i,NONE)
				val h = M.scalarMatrix(m,F.one)
				val _ = M.updateSubMatrix(h,i,i,householder x)
			    in
				loop(Int.+(i,1), M.*(qm,h), M.*(h,rm))
			    end
    in
	loop(0, M.scalarMatrix(m,F.one), mat)
    end
fun solveUpperTriangular(r,b) =
    let open Array
	val n = Array2.nCols r
	val x = array(n, F.zero)
	fun loop k =
	    let val index = Int.min(Int.-(n,1),Int.+(k,1))
		val _ = update(x,k,
			       F./(F.-(sub(b,k),
				       A.dot(A.slice(x,index,NONE),
					     A.slice(M.row(r,k),index,NONE))),
				   Array2.sub(r,k,k)))
	    in
		if k=0 then x else loop(Int.-(k,1))
	    end
    in
	loop (Int.-(n,1))
    end
fun lsqr(a,b) =
    let val {q,r} = qr a
	val n = Array2.nCols r
    in
	solveUpperTriangular(M.subMatrix(r, 0, Int.-(n,1), 0, Int.-(n,1)),
			     M.multa(M.transpose(q), b))
    end
fun pow(x,1) = x
  | pow(x,n) = F.*(x,pow(x,Int.-(n,1)))  
fun polyfit(x,y,n) =
    let open Array2
	val a = tabulate RowMajor (Array.length x,
				   Int.+(n,1),
				   fn (i,j) => if j=0 then F.one else
					       pow(Array.sub(x,i),j))
    in
	lsqr(a,y)
    end
end

We can then show the examples:

structure RealRadicalCategoryField : RADCATFIELD = struct
open Real
val one = 1.0
val zero = 0.0
val sign = real o Real.sign
val sqrt = Real.Math.sqrt
end

structure Q = QR(RealRadicalCategoryField);

let
    val mat = Array2.fromList [[12.0, ~51.0, 4.0], [6.0, 167.0, ~68.0], [~4.0, 24.0, ~41.0]]
    val {q,r} = Q.qr(mat)
in
    {q=Q.toList q; r=Q.toList r}
end;
(* output *)
val it =
  {q=[[~0.857142857143,0.394285714286,0.331428571429],
      [~0.428571428571,~0.902857142857,~0.0342857142857],
      [0.285714285714,~0.171428571429,0.942857142857]],
   r=[[~14.0,~21.0,14.0],[5.97812397875E~18,~175.0,70.0],
      [4.47505280695E~16,0.0,~35.0]]} : {q:real list list, r:real list list}

let open Array
    val x = fromList [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
    val y = fromList [1.0, 6.0, 17.0, 34.0, 57.0, 86.0, 121.0, 162.0, 209.0, 262.0, 321.0]
in
    Q.polyfit(x, y, 2)
end;

(* output *)
val it = [|1.0,2.0,3.0|] : real array

Stata

See QR decomposition in Stata help.

mata
: qrd(a=(12,-51,4\6,167,-68\-4,24,-41),q=.,r=.)

: a
         1     2     3
    +-------------------+
  1 |   12   -51     4  |
  2 |    6   167   -68  |
  3 |   -4    24   -41  |
    +-------------------+

: q
                  1              2              3
    +----------------------------------------------+
  1 |  -.8571428571    .3942857143    .3314285714  |
  2 |  -.4285714286   -.9028571429   -.0342857143  |
  3 |   .2857142857   -.1714285714    .9428571429  |
    +----------------------------------------------+

: r
          1      2      3
    +----------------------+
  1 |   -14    -21     14  |
  2 |     0   -175     70  |
  3 |     0      0    -35  |
    +----------------------+

Tcl

Assuming the presence of the Tcl solutions to these tasks: Element-wise operations, Matrix multiplication, Matrix transposition

Translation of: Common Lisp

package require Tcl 8.5
namespace path {::tcl::mathfunc ::tcl::mathop}
proc sign x {expr {$x == 0 ? 0 : $x < 0 ? -1 : 1}}
proc norm vec {
    set s 0
    foreach x $vec {set s [expr {$s + $x**2}]}
    return [sqrt $s]
}
proc unitvec n {
    set v [lrepeat $n 0.0]
    lset v 0 1.0
    return $v
}
proc I n {
    set m [lrepeat $n [lrepeat $n 0.0]]
    for {set i 0} {$i < $n} {incr i} {lset m $i $i 1.0}
    return $m
}

proc arrayEmbed {A B row col} {
    # $A will be copied automatically; Tcl values are copy-on-write
    lassign [size $B] mb nb
    for {set i 0} {$i < $mb} {incr i} {
	for {set j 0} {$j < $nb} {incr j} {
	    lset A [expr {$row + $i}] [expr {$col + $j}] [lindex $B $i $j]
	}
    }
    return $A
}

# Unlike the Common Lisp version, here we use a specialist subcolumn
# extraction function: like that, there's a lot less intermediate memory allocation
# and the code is actually clearer.
proc subcolumn {A size column} {
    for {set i $column} {$i < $size} {incr i} {lappend x [lindex $A $i $column]}
    return $x
}

proc householder A {
    lassign [size $A] m
    set U [m+ $A [.* [unitvec $m] [expr {[norm $A] * [sign [lindex $A 0 0]]}]]]
    set V [./ $U [lindex $U 0 0]]
    set beta [expr {2.0 / [lindex [matrix_multiply [transpose $V] $V] 0 0]}]
    return [m- [I $m] [.* [matrix_multiply $V [transpose $V]] $beta]]
}

proc qrDecompose A {
    lassign [size $A] m n
    set Q [I $m]
    for {set i 0} {$i < ($m==$n ? $n-1 : $n)} {incr i} {
	# Construct the Householder matrix
	set H [arrayEmbed [I $m] [householder [subcolumn $A $n $i]] $i $i]
	# Apply to build the decomposition
	set Q [matrix_multiply $Q $H]
	set A [matrix_multiply $H $A]
    }
    return [list $Q $A]
}

Demonstrating:

set demo [qrDecompose {{12 -51 4} {6 167 -68} {-4 24 -41}}]
puts "==Q=="
print_matrix [lindex $demo 0] "%f"
puts "==R=="
print_matrix [lindex $demo 1] "%.1f"

Output:

==Q==
-0.857143  0.394286  0.331429 
-0.428571 -0.902857 -0.034286 
 0.285714 -0.171429  0.942857 
==R==
-14.0  -21.0  14.0 
  0.0 -175.0  70.0 
  0.0    0.0 -35.0 

VBA

Translation of: Phix

Option Base 1
Private Function vtranspose(v As Variant) As Variant
'-- transpose a vector of length m into an mx1 matrix,
'--                       eg {1,2,3} -> {1;2;3}
    vtranspose = WorksheetFunction.Transpose(v)
End Function
 
Private Function mat_col(a As Variant, col As Integer) As Variant
    Dim res() As Double
    ReDim res(UBound(a))
    For i = col To UBound(a)
        res(i) = a(i, col)
    Next i
    mat_col = res
End Function
 
Private Function mat_norm(a As Variant) As Double
    mat_norm = Sqr(WorksheetFunction.SumProduct(a, a))
End Function
 
Private Function mat_ident(n As Integer) As Variant
    mat_ident = WorksheetFunction.Munit(n)
End Function

Private Function sq_div(a As Variant, p As Double) As Variant
    Dim res() As Variant
    ReDim res(UBound(a))
    For i = 1 To UBound(a)
        res(i) = a(i) / p
    Next i
    sq_div = res
End Function

Private Function sq_mul(p As Double, a As Variant) As Variant
    Dim res() As Variant
    ReDim res(UBound(a), UBound(a, 2))
    For i = 1 To UBound(a)
        For j = 1 To UBound(a, 2)
            res(i, j) = p * a(i, j)
        Next j
    Next i
    sq_mul = res
End Function

Private Function sq_sub(x As Variant, y As Variant) As Variant
    Dim res() As Variant
    ReDim res(UBound(x), UBound(x, 2))
    For i = 1 To UBound(x)
        For j = 1 To UBound(x, 2)
            res(i, j) = x(i, j) - y(i, j)
        Next j
    Next i
    sq_sub = res
End Function

Private Function matrix_mul(x As Variant, y As Variant) As Variant
    matrix_mul = WorksheetFunction.MMult(x, y)
End Function

Private Function QRHouseholder(ByVal a As Variant) As Variant
    Dim columns As Integer: columns = UBound(a, 2)
    Dim rows As Integer: rows = UBound(a)
    Dim m As Integer: m = WorksheetFunction.Max(columns, rows)
    Dim n As Integer: n = WorksheetFunction.Min(rows, columns)
    I_ = mat_ident(m)
    Q_ = I_
    Dim q As Variant
    Dim u As Variant, v As Variant, j As Integer
    For j = 1 To WorksheetFunction.Min(m - 1, n)
        u = mat_col(a, j)
        u(j) = u(j) - mat_norm(u)
        v = sq_div(u, mat_norm(u))
        q = sq_sub(I_, sq_mul(2, matrix_mul(vtranspose(v), v)))
        a = matrix_mul(q, a)
        Q_ = matrix_mul(Q_, q)
    Next j
 
    '-- Get the upper triangular matrix R.
    Dim R() As Variant
    ReDim R(m, n)
    For i = 1 To m 'in Phix this is n
        For j = 1 To n 'in Phix this is i to n. starting at 1 to fill zeroes
            R(i, j) = a(i, j)
        Next j
    Next i
    Dim res(2) As Variant
    res(1) = Q_
    res(2) = R
    QRHouseholder = res
End Function

Private Sub pp(m As Variant)
    For i = 1 To UBound(m)
        For j = 1 To UBound(m, 2)
            Debug.Print Format(m(i, j), "0.#####"),
        Next j
        Debug.Print
    Next i
End Sub
Public Sub main()
    a = [{12, -51,   4; 6, 167, -68; -4,  24, -41;-1,1,0;2,0,3}]
    result = QRHouseholder(a)
    q = result(1)
    r_ = result(2)
    Debug.Print "A"
    pp a
    Debug.Print "Q"
    pp q
    Debug.Print "R"
    pp r_
    Debug.Print "Q * R"
    pp matrix_mul(q, r_)
End Sub

Output:

A
12,           -51,          4,            
6,            167,          -68,          
-4,           24,           -41,          
-1,           1,            0,            
2,            0,            3,            
Q
0,84641       -0,39129      -0,34312      0,06641       -0,09126      
0,42321       0,90409       0,02927       0,01752       -0,04856      
-0,28214      0,17042       -0,93286      -0,02237      0,14365       
-0,07053      0,01404       0,0011        0,99738       0,00728       
0,14107       -0,01666      0,10577       0,00291       0,98419       
R
14,17745      20,66663      -13,40157     
0,            175,04254     -70,08031     
0,            0,            35,20154      
0,            0,            0,            
0,            0,            0,            
Q * R
12,           -51,          4,            
6,            167,          -68,          
-4,           24,           -41,          
-1,           1,            0,            
2,            0,            3, 

Least squares

Public Sub least_squares()
    x = [{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}]
    y = [{1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}]
    Dim a() As Double
    ReDim a(UBound(x), 3)
    For i = 1 To UBound(x)
        For j = 1 To 3
            a(i, j) = x(i) ^ (j - 1)
        Next j
    Next i
    result = QRHouseholder(a)
    q = result(1)
    r_ = result(2)
    t = WorksheetFunction.Transpose(q)
    b = matrix_mul(t, vtranspose(y))
    Dim z(3) As Double
    For k = 3 To 1 Step -1
        Dim s As Double: s = 0
        If k < 3 Then
            For j = k + 1 To 3
                s = s + r_(k, j) * z(j)
            Next j
        End If
        z(k) = (b(k, 1) - s) / r_(k, k)
    Next k
    Debug.Print "Least-squares solution:",
    For i = 1 To 3
        Debug.Print Format(z(i), "0.#####"),
    Next i
End Sub

Output:

Least-squares solution:     1,            2,            3,  

Wren

Translation of: C

Library: Wren-matrix

Library: Wren-fmt

import "./matrix" for Matrix
import "./fmt" for Fmt

var minor = Fn.new { |x, d|
    var nr = x.numRows
    var nc = x.numCols
    var m = Matrix.new(nr, nc)
    for (i in 0...d) m[i, i] = 1
    for (i in d...nr) {
        for (j in d...nc) m[i, j] = x[i, j]
    }
    return m
}

var vmadd = Fn.new { |a, b, s|
    var n = a.count
    var c = List.filled(n, 0)
    for (i in 0...n) c[i] = a[i] + s * b[i]
    return c
}

var vmul = Fn.new { |v|
    var n = v.count
    var x = Matrix.new(n, n)
    for (i in 0...n) {
        for (j in 0...n) x[i, j] = -2 * v[i] * v[j]
    }
    for (i in 0...n) x[i, i] = x[i, i] + 1
    return x
}

var vnorm = Fn.new { |x|
    var n = x.count
    var sum = 0
    for (i in 0...n) sum = sum + x[i] * x[i]
    return sum.sqrt
}

var vdiv = Fn.new { |x, d|
    var n = x.count
    var y = List.filled(n, 0)
    for (i in 0...n) y[i] = x[i] / d
    return y
}

var mcol = Fn.new { |m, c|
    var n = m.numRows
    var v = List.filled(n, 0)
    for (i in 0...n) v[i] = m[i, c]
    return v
}

var householder = Fn.new { |m|
    var nr = m.numRows
    var nc = m.numCols
    var q = List.filled(nr, null)
    var z = m.copy()
    var k = 0
    while (k < nc && k < nr-1) {
       var e = List.filled(nr, 0)
       z = minor.call(z, k)
       var x = mcol.call(z, k)
       var a = vnorm.call(x)
       if (m[k, k] > 0) a = -a
       for (i in 0...nr) e[i] = (i == k) ? 1 : 0
       e = vmadd.call(x, e, a)
       e = vdiv.call(e, vnorm.call(e))
       q[k] = vmul.call(e)
       z = q[k] * z
       k = k + 1
    }
    var Q = q[0]
    var R = q[0] * m
    var i = 1
    while (i < nc && i < nr-1) {
        Q = q[i] * Q
        i = i + 1
    }
    R = Q * m
    Q = Q.transpose
    return [Q, R]
}

var inp = [
    [12, -51,   4],
	[ 6, 167, -68],
	[-4,  24, -41],
	[-1,   1,   0],
	[ 2,   0,   3]
]
var x = Matrix.new(inp)
var res = householder.call(x)
var Q = res[0]
var R = res[1]
var m = Q * R
System.print("Q:")
Fmt.mprint(Q, 8, 3)
System.print("\nR:")
Fmt.mprint(R, 8, 3)
System.print("\nQ * R:")
Fmt.mprint(m, 8, 3)

Output:

Q:
|   0.846   -0.391    0.343    0.082    0.078|
|   0.423    0.904   -0.029    0.026    0.045|
|  -0.282    0.170    0.933   -0.047   -0.137|
|  -0.071    0.014   -0.001    0.980   -0.184|
|   0.141   -0.017   -0.106   -0.171   -0.969|

R:
|  14.177   20.667  -13.402|
|  -0.000  175.043  -70.080|
|   0.000    0.000  -35.202|
|  -0.000   -0.000   -0.000|
|   0.000    0.000   -0.000|

Q * R:
|  12.000  -51.000    4.000|
|   6.000  167.000  -68.000|
|  -4.000   24.000  -41.000|
|  -1.000    1.000   -0.000|
|   2.000   -0.000    3.000|

zkl

var [const] GSL=Import("zklGSL");	// libGSL (GNU Scientific Library)
A:=GSL.Matrix(3,3).set(12.0, -51.0,   4.0,
			6.0, 167.0, -68.0,
			4.0,  24.0, -41.0);
Q,R:=A.QRDecomp();
println("Q:\n",Q.format());
println("R:\n",R.format());
println("Q*R:\n",(Q*R).format());

Output:

Q:
     -0.86,      0.47,     -0.22
     -0.43,     -0.88,     -0.20
     -0.29,     -0.08,      0.95
R:
    -14.00,    -34.71,     37.43
      0.00,   -172.80,     65.07
      0.00,      0.00,    -26.19
Q*R:
     12.00,    -51.00,      4.00
      6.00,    167.00,    -68.00
      4.00,     24.00,    -41.00