QR decomposition: Difference between revisions

Content added Content deleted

Inline

Revision as of 16:30, 12 July 2011

Any rectangular Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle m \times n} matrix Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathit A} can be decomposed to a product of a orthogonal matrix Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathit Q} and a upper (right) triangular matrix Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathit R} , as described in QR decomposition.

There are several methods to compute Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathit Q} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathit R} , but in this task, the method of Householder reflections should be used. Multiplying a given vector Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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}}}

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

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle H \; a = -\textrm{sign}(a_1) \; \|a\|_2 \; e_1}

and applying Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathit H} on the matrix Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathit A} zeroes all subdiagonal elements of the first column:

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle H_1 \; A = \begin{pmatrix} r_{11} & r_{12} & r_{13} \\ 0 & * & * \\ 0 & * & * \end{pmatrix}}

In the second step, the second column of Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathit A} , we want to zero all elements but the first two, which means that we have to calculate Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathit H} with the first column of the submatrix (denoted *), not on the whole second column of Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathit A} .

To get Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle H_2} , we then embed the new Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathit H} into an Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle m \times n} identity:

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathit A} and thus transform it into Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathit R} .

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle H_n \; ... \; H_3 H_2 H_1 A = R}

The product of all the Householder matrices Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathit H} , for every column, in reverse order, will then yield the orthogonal matrix Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathit Q} .

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathit A x = b} by solving

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle R \; x = Q^T \; b}

When Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathit R} is not square, i.e. Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle m > n} we have to cut off the Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathit m - n} zero padded bottom rows.

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle R = \begin{pmatrix} R_1 \\ 0 \end{pmatrix}}

and the same for the RHS:

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle Q^T \; b = \begin{pmatrix} q_1 \\ q_2 \end{pmatrix}}

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

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle R_1 \; x = q_1}

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

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\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.

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 (mmul (mtp v) v))))
   
   (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

This is not complete. First part of the translation of the Common Lisp code. <lang d>import std.stdio, std.math, std.algorithm, std.traits;

T norm(T)(T[] array) if (isFloatingPoint!T) {

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

}

T[][] makeUnitVector(T)(size_t dim) if (isFloatingPoint!T) {

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

}

T[][] arrayRange(T)(T[][] A, size_t ma, size_t mb,

                            size_t na, size_t nb) {
   const size_t mm = 1 + (mb - ma);
   const size_t nn = 1 + (nb - na);
   auto B = new T[][](nn, mm);

   foreach (i; 0 .. mm)
       foreach (j; 0 .. nn)
           B[i][j] = A[ma + i][na + j];
   return B;

}

size_t rows(T)(T[][] A) { return A.length; } size_t cols(T)(T[][] A) { return A.length ? A[0].length : 0; } T[][] mcol(T)(T[][] A, size_t n) {

   return arrayRange(A, 0, rows(A) - 1, n, n);

} T[][] mrow(T)(T[][] A, size_t n) {

   return arrayRange(A, n, n, 0, cols(A) -1);

}

template msum(T) { alias elementwiseMat!(q{ + }, T, T[][]) msum; } template msub(T) { alias elementwiseMat!(q{ - }, T, T[][]) msub; } template smul(T) { alias elementwiseMat!(q{ * }, T, T) smul; } template sdiv(T) { alias elementwiseMat!(q{ / }, T, T) sdiv; }

T[][] arrayEmbed(T)(T[][] A, T[][] B, size_t row, size_t col) {

   const ma = rows(A);
   const na = cols(A);
   const mb = rows(B);
   const nb = cols(B);
   auto C = new T[][](ma, na);

   foreach (i; 0 .. ma)
       foreach (j; 0 .. na)
           C[i][j] = A[i][j];

   foreach (i; 0 .. mb)
       foreach (j; 0 .. nb)
           C[row + i][col + j] = B[i][j];

   return C;

}</lang>

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>