QR decomposition: Difference between revisions

From Rosetta Code
Content added Content deleted
(Complete D version)
(Added tag to D entry)
Line 436: Line 436:


=={{header|D}}==
=={{header|D}}==
{{trans|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,
<lang d>import std.stdio, std.math, std.algorithm, std.traits,
std.typecons, std.numeric, std.range, std.conv;
std.typecons, std.numeric, std.range, std.conv;

Revision as of 10:28, 14 July 2011

QR decomposition is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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 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} , with the Householder 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 H} , which is given as

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 = I - \frac {2} {u^T u} u u^T}

reflects 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} about a plane given by its normal 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 u} . When the normal vector of the plane 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 u} is given as

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 u = a - \|a\|_2 \; e_1}

then the transformation reflects 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} onto the first standard basis 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 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 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_1} :

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 u = a + \textrm{sign}(a_1)\|a\|_2 \; e_1}

and normalize with respect to the first element:

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 v = \frac{u}{u_1}}

The equation for 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} thus becomes:

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 = I - \frac {2} {v^T v} v v^T}

or, in another form

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 = I - \beta v v^T}

with

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 \beta = \frac {2} {v^T v}}

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 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} 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:

and the usage for linear least squares problems on the example from Polynomial_regression.

C

<lang C>#include <stdio.h>

  1. include <stdlib.h>
  2. 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)))

  1. 2A((-0.85 0.39 0.33)
   (-0.42 -0.90 -0.03)
   ( 0.28 -0.17  0.94))
  1. 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))
  1. 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:

<lang j> QR x:12 _51 4,6 167 _68,:_4 24 _41 ┌────────────────────┬──────────┐ │ 6r7 _69r175 _58r175│14 21 _14│ │ 3r7 158r175 6r175│ 0 175 _70│ │_2r7 6r35 _33r35│ 0 0 35│ └────────────────────┴──────────┘</lang>

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>

Retrieved from "https://rosettacode.org/wiki/QR_decomposition?oldid=189625"