user
Any rectangular 

m
×
n
{\displaystyle {\displaystyle m\times n}}
 matrix 

A
{\displaystyle {\displaystyle {\mathit {A}}}}
 can be decomposed to a product of an orthogonal matrix 

Q
{\displaystyle {\displaystyle {\mathit {Q}}}}
 and an upper (right) triangular matrix 

R
{\displaystyle {\displaystyle {\mathit {R}}}}
, as described in QR decomposition.

**Task**

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

:   :   A
        =


        (

        12

        −
        51

        4


        6

        167

        −
        68


        −
        4

        24

        −
        41
        )
        {\displaystyle {\displaystyle 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 

a
{\displaystyle {\displaystyle {\mathit {a}}}}
, for example the first column of matrix 

A
{\displaystyle {\displaystyle {\mathit {A}}}}
, with the Householder matrix 

H
{\displaystyle {\displaystyle {\mathit {H}}}}
, which is given as

:   :   H
        =
        I
        −


        2


        u

        T
        u
        u

        u

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

reflects 

a
{\displaystyle {\displaystyle {\mathit {a}}}}
 about a plane given by its normal vector 

u
{\displaystyle {\displaystyle {\mathit {u}}}}
. When the normal vector of the plane 

u
{\displaystyle {\displaystyle {\mathit {u}}}}
 is given as

:   :   u
        =
        a
        −
        ‖
        a

        ‖

        2


        e

        1
        {\displaystyle {\displaystyle u=a-\|a\|\_{2}\;e\_{1}}}

then the transformation reflects 

a
{\displaystyle {\displaystyle {\mathit {a}}}}
 onto the first standard basis vector

:   :   e

        1
        =
        [
        1

        0

        0

        .
        .
        .

        ]

        T
        {\displaystyle {\displaystyle 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
{\displaystyle {\displaystyle a\_{1}}}
:

:   :   u
        =
        a
        +


        sign
        (

        a

        1
        )
        ‖
        a

        ‖

        2


        e

        1
        {\displaystyle {\displaystyle u=a+{\textrm {sign}}(a\_{1})\|a\|\_{2}\;e\_{1}}}

and normalize with respect to the first element:

:   :   v
        =


        u

        u

        1
        {\displaystyle {\displaystyle v={\frac {u}{u\_{1}}}}}

The equation for 

H
{\displaystyle {\displaystyle H}}
 thus becomes:

:   :   H
        =
        I
        −


        2


        v

        T
        v
        v

        v

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

or, in another form

:   :   H
        =
        I
        −
        β
        v

        v

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

with

:   :   β
        =


        2


        v

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

Applying 

H
{\displaystyle {\displaystyle {\mathit {H}}}}
 on 

a
{\displaystyle {\displaystyle {\mathit {a}}}}
 then gives

:   :   H

        a
        =
        −


        sign
        (

        a

        1
        )

        ‖
        a

        ‖

        2


        e

        1
        {\displaystyle {\displaystyle H\;a=-{\textrm {sign}}(a\_{1})\;\|a\|\_{2}\;e\_{1}}}

and applying 

H
{\displaystyle {\displaystyle {\mathit {H}}}}
 on the matrix 

A
{\displaystyle {\displaystyle {\mathit {A}}}}
 zeroes all subdiagonal elements of the first column:

:   :   H

        1

        A
        =


        (

        r

        11


        r

        12


        r

        13


        0

        ∗

        ∗


        0

        ∗

        ∗
        )
        {\displaystyle {\displaystyle H\_{1}\;A={\begin{pmatrix}r\_{11}&r\_{12}&r\_{13}\\0&\*&\*\\0&\*&\*\end{pmatrix}}}}

In the second step, the second column of 

A
{\displaystyle {\displaystyle {\mathit {A}}}}
, we want to zero all elements but the first two, which means that we have to calculate 

H
{\displaystyle {\displaystyle {\mathit {H}}}}
 with the first column of the *submatrix* (denoted \*), not on the whole second column of 

A
{\displaystyle {\displaystyle {\mathit {A}}}}
.

To get 

H

2
{\displaystyle {\displaystyle H\_{2}}}
, we then embed the new 

H
{\displaystyle {\displaystyle {\mathit {H}}}}
 into an 

m
×
n
{\displaystyle {\displaystyle m\times n}}
 identity:

:   :   H

        2
        =


        (

        1

        0

        0


        0

        H

        0

        )
        {\displaystyle {\displaystyle 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 

A
{\displaystyle {\displaystyle {\mathit {A}}}}
 and thus transform it into 

R
{\displaystyle {\displaystyle {\mathit {R}}}}
.

:   :   H

        n

        .
        .
        .


        H

        3

        H

        2

        H

        1
        A
        =
        R
        {\displaystyle {\displaystyle H\_{n}\;...\;H\_{3}H\_{2}H\_{1}A=R}}

The product of all the Householder matrices 

H
{\displaystyle {\displaystyle {\mathit {H}}}}
, for every column, in reverse order, will then yield the orthogonal matrix 

Q
{\displaystyle {\displaystyle {\mathit {Q}}}}
.

:   :   H

        1

        H

        2

        H

        3

        .
        .
        .


        H

        n
        =
        Q
        {\displaystyle {\displaystyle H\_{1}H\_{2}H\_{3}\;...\;H\_{n}=Q}}

The QR decomposition should then be used to solve linear least squares (Multiple regression) problems 

A
x
=
b
{\displaystyle {\displaystyle {\mathit {A}}x=b}}
 by solving

:   :   R

        x
        =

        Q

        T

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

When 

R
{\displaystyle {\displaystyle {\mathit {R}}}}
 is not square, i.e. 

m
>
n
{\displaystyle {\displaystyle m>n}}
 we have to cut off the 

m
−
n
{\displaystyle {\displaystyle {\mathit {m}}-n}}
 zero padded bottom rows.

:   :   R
        =


        (

        R

        1


        0
        )
        {\displaystyle {\displaystyle R={\begin{pmatrix}R\_{1}\\0\end{pmatrix}}}}

and the same for the RHS:

:   :   Q

        T

        b
        =


        (

        q

        1

        q

        2
        )
        {\displaystyle {\displaystyle 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
        {\displaystyle {\displaystyle R\_{1}\;x=q\_{1}}}
assistant
(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


```