QR decomposition

Any rectangular $m\times n$ matrix ${\mathit {A}}$ can be decomposed to a product of a orthogonal matrix ${\mathit {Q}}$ and a 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}

C

<lang 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(void * a, int m, int n) { mat x = matrix_new(m, n); for (int i = 0; i < m; i++) for (int j = 0; j < n; j++) x->v[i][j] = ((double(*)[n])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(" %5.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->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->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 }, };

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

printf("Q\n"); matrix_show(Q); printf("R\n"); matrix_show(R);

matrix_delete(x); matrix_delete(R); matrix_delete(Q); return 0; }</lang>Output:<lang>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</lang>

Common Lisp

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

Helper functions: <lang lisp>(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))

</lang>

Main routines: <lang lisp> (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)))

</lang>

Example 1:

<lang lisp>(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))</lang>

Example 2, Polynomial regression:

<lang lisp>(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))</lang>

<lang lisp>;; 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))</lang>

D

Translation of: Common Lisp

Uses the functions copied from Element-wise_operations, Matrix multiplication, and Matrix transposition. <lang d>import std.stdio, std.math, std.algorithm, std.traits,

      std.typecons, std.numeric, std.range, std.conv;

T[][] elementwiseMat(string op, T, U)(in T[][] A, in U B) pure if (is(U == T) || is(U == T[][])) {

   static if (is(U == T[][]))
       assert(A.length == B.length);
   if (!A.length)
       return null;
   auto R = new typeof(return)(A.length, A[0].length);

   foreach (r, row; A)
       static if (is(U == T)) {
           R[r][] = mixin("row[] " ~ op ~ "B");
       } else {
           assert(row.length == B[r].length);
           R[r][] = mixin("row[] " ~ op ~ "B[r][]");
       }

   return R;

}

T[][] msum(T)(in T[][] A, in T[][] B) pure {

   return elementwiseMat!(q{ + }, T, T[][])(A, B);

} T[][] msub(T)(in T[][] A, in T[][] B) pure {

   return elementwiseMat!(q{ - }, T, T[][])(A, B);

} T[][] pmul(T)(in T[][] A, in T x) pure {

   return elementwiseMat!(q{ * }, T, T)(A, x);

} T[][] pdiv(T)(in T[][] A, in T x) pure {

   return elementwiseMat!(q{ / }, T, T)(A, x);

}

bool isRectangular(T)(in T[][] matrix) pure nothrow {

   foreach (row; matrix)
       if (row.length != matrix[0].length)
           return false;
   return true;

}

T[][] matMul(T)(in T[][] a, in T[][] b) pure nothrow in {

   assert(isRectangular(a) && isRectangular(b) &&
          a[0].length == b.length);

} body {

   auto result = new T[][](a.length, b[0].length);
   auto aux = new T[b.length];
   foreach (j; 0 .. b[0].length) {
       foreach (k; 0 .. b.length)
           aux[k] = b[k][j];
       foreach (i; 0 .. a.length)
           result[i][j] = dotProduct(a[i], aux);
   }
   return result;

}

string prettyPrint(T)(in T[][] A) {

   return "[" ~ array(map!text(A)).join(",\n ") ~ "]";

}

Unqual!T[][] transpose(T)(in T[][] m) pure nothrow {

   auto r = new Unqual!T[][](m[0].length, m.length);
   foreach (nr, row; m)
       foreach (nc, c; row)
           r[nc][nr] = c;
   return r;

}

T norm(T)(in T[][] array) {

   return sqrt(reduce!q{ a + b ^^ 2 }(cast(T)0,
                                      transversal(array, 0)));

}

T[][] makeUnitVector(T)(in size_t dim) pure nothrow {

   auto result = new T[][](dim, 1);
   foreach (row; result)
       row[] = 0;
   result[0][0] = 1;
   return result;

}

/// Return a nxn identity matrix. T[][] matId(T)(in size_t n) pure nothrow {

   auto Id = new T[][](n, n);
   foreach (r, row; Id) {
       row[] = 0;
       row[r] = 1;
   }
   return Id;

}

Unqual!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 Unqual!T[][](mb - ma + 1, nb - na + 1);
   foreach (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, rows(A)-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 (i, arow; A)
       C[i][] = arow[]; // some wasted copies
   foreach (i, brow; B)
       C[row + i][col .. col + brow.length] = brow[];
   return C;

}

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

T[][] makeHouseholder(T)(T[][] a) {

   const size_t m = rows(a);
   const T s = sgn(a[0][0]);
   T[][] e = makeUnitVector!T(m);
   T[][] u = msum(a, pmul(e, norm(a) * s));
   T[][] v = pdiv(u, u[0][0]);
   T beta = 2.0 / matMul(transpose(v), v)[0][0];
   return msub(matId!T(m), pmul(matMul(v, transpose(v)), beta));

}

Tuple!(T[][],"Q", T[][],"R") QRdecomposition(T)(T[][] A) {

   const m = rows(A);
   const n = cols(A);
   T[][] Q = matId!T(m);

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

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

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

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

       // The product of all H matrices with A from the LHS is the
       // upper triangular matrix R.
       A  = matMul(H, 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 {

   const size_t n = cols(R);
   auto x = new T[][](n, 1);

   foreach_reverse (k; 0 .. n) {
       T tot = 0;
       foreach (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, T[][] b) {

   const qr = QRdecomposition(A);
   const size_t n = cols(qr.R);
   return solveUpperTriangular(
       slice2D(qr.R, 0, n-1, 0, n-1),
       slice2D(matMul(transpose(qr.Q), b), 0, n-1, 0, 0));

}

Unqual!T[][] polyFit(T)(in T[][] x, in T[][] y, in size_t n) {

   const size_t m = cols(x);
   auto A = new Unqual!T[][](m, n+1);
   foreach (i, row; A)
       foreach (j, ref item; row)
           item = x[0][i] ^^ j;
   return lsqr(A, transpose(y));

}

void main() {

   const qr = QRdecomposition([[12.0, -51,   4],
                               [ 6.0, 167, -68],
                               [-4.0,  24, -41]]);
   writeln(prettyPrint(qr.Q));
   writeln(prettyPrint(qr.R), "\n");

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

}</lang>

Output:

[[-0.857143, 0.394286, 0.331429],
 [-0.428571, -0.902857, -0.0342857],
 [0.285714, -0.171429, 0.942857]]
[[-14, -21, 14],
 [2.26656e-16, -175, 70],
 [4.72101e-16, -7.42462e-16, -35]]

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

J

From j:Essays/QR Decomposition

<lang j>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.

)</lang>

Example use:

Polynomial fitting using QR reduction:

<lang j> 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</lang>

Tcl

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

Translation of: Common Lisp

<lang tcl>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
   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

}

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]

}</lang> Demonstrating: <lang tcl>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"</lang> 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