QR decomposition
Any rectangular matrix can be decomposed to a product of a orthogonal matrix and a upper (right) triangular matrix , as described in QR decomposition.
You are encouraged to solve this task according to the task description, using any language you may know.
Task
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. The method of Householder reflections should be used:
Method
Multiplying a given vector , for example the first column of matrix , with the Householder matrix , which is given as
reflects about a plane given by its normal vector . When the normal vector of the plane is given as
then the transformation reflects onto the first standard basis vector
which means that all entries but the first become zero. To avoid numerical cancellation errors, we should take the opposite sign of :
and normalize with respect to the first element:
The equation for thus becomes:
or, in another form
with
Applying on then gives
and applying on the matrix zeroes all subdiagonal elements of the first column:
In the second step, the second column of , we want to zero all elements but the first two, which means that we have to calculate with the first column of the submatrix (denoted *), not on the whole second column of .
To get , we then embed the new into an identity:
This is how we can, column by column, remove all subdiagonal elements of and thus transform it into .
The product of all the Householder matrices , for every column, in reverse order, will then yield the orthogonal matrix .
The QR decomposition should then be used to solve linear least squares (Multiple regression) problems by solving
When is not square, i.e. we have to cut off the zero padded bottom rows.
and the same for the RHS:
Finally, solve the square upper triangular system by back substitution:
Ada
Output matches that of Matlab solution, not tested with other matrices. <lang Ada> 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;</lang>
- 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
Axiom
The following provides a generic QR decomposition for arbitrary precision floats, double floats and exact calculations: <lang Axiom>)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 s := signValue(a(1)) e := unitVector(m) u := a + length(a)*s*e 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)</lang>
This can be called using: <lang Axiom>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)</lang> With output in exact form: <lang Axiom>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)</lang>
The calculations are comparable to those from the QR decomposition in R, which by default uses LINPACK. R's LAPACK implementation uses column pivoting. Pivots can be implemented using: <lang Axiom>)abbrev package TESTP2 TestPackage2 TestPackage2(R:Join(Field,RadicalCategory)): with
positionMax: (Vector(R),Integer) -> Integer ++ positionMax(v,i) finds the maximum in a vector v from position i qrPivot: Matrix(R) -> Record(q:Matrix(R),r:Matrix(R),perm:Vector(Integer)) ++ qrPivot(m) finds the QR decomposition with pivots from matrix m == add import TestPackage(R) positionMax(v,i) == vsubset := vector [v(j) for j in i..#v] vmax : R := reduce((x:R,y:R):R +-> if signValue(x-y)=-1 then y else x,v) i+position(vmax,v)-1 qrPivot(a) == -- http://www.netlib.org/lapack/lawnspdf/lawn114.pdf, Figure 1 (m,n) := (nrows a, ncols a) perms : Vector Integer := vector [j for j in 1..n] colnorms : Vector(R) := vector [dot(column(a,j),column(a,j)) for j in 1..n] qm := scalarMatrix(m,1) rm : Matrix(R) := copy a for j in 1..(if m=n then n-1 else n) repeat p := positionMax(colnorms,j)
if colnorms(p)=0@R then return [qm,rm,perms] if j ~= p then swap!(perms,p,j) swapColumns!(rm,p,j) swap!(colnorms,p,j) h := householder column(subMatrix(rm,j,m,j,j),1) setsubMatrix!(qm,1,j,subMatrix(qm,1,m,j,n)*h) setsubMatrix!(rm,j,1,h*subMatrix(rm,j,m,1,n)) for k in (j+1)..n repeat colnorms(k) := colnorms(k) - rm(j,k)^2
[qm,rm,perms]</lang>
Called using:<lang Axiom>qrPivot m</lang>With output:<lang Axiom> + +---+ +-+ +---+ +-+ +---+ +---+ +
| 51\|634 3957\|2 \|317 29\|2 \|317 \|634 | |- -------- - -------------- ------------------ | | 4438 110950 55475 | | | | +---+ +-+ +---+ +-+ +---+ +---+| [q= |167\|634 1401\|2 \|317 3\|2 \|317 \|634 |, |--------- - -------------- - -----------------| | 4438 110950 55475 | | | | +---+ +-+ +---+ +-+ +---+ +---+ | | 12\|634 134\|2 \|317 33\|2 \|317 \|634 | | -------- ------------- ------------------ | + 2219 11095 22190 + + +---+ +---+ + | +---+ 21\|634 896\|634 | |7\|634 -------- - --------- | | 317 317 | | | | +-+ +---+ +-+ +---+ | r= | 175\|2 \|317 70\|2 \|317 |, perm= [2,1,3]] | 0 - ------------- ------------ | | 317 317 | | | | +-+ +---+ +---+| | 35\|2 \|317 \|634 | | 0 0 - ------------------| + 634 +
Type: Record(q: Matrix(AlgebraicNumber),r: Matrix(AlgebraicNumber),perm: Vector(Integer))</lang>
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
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]]
Go
A fairly close port of the Common Lisp solution, this solution uses the gomatrix for supporting functions. Note though, that gomatrix 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 Gomatrix QR source shows that it computes the decomposition more directly. <lang go>package main
import (
"code.google.com/p/gomatrix/matrix" "fmt" "math"
)
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
}</lang> 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}
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>
Mathematica
<lang Mathematica>{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</lang>
MATLAB / Octave
<lang Matlab> A = [12 -51 4
6 167 -68 -4 24 -41]; [Q,R]=qr(A) </lang>
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
R
<lang 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)</lang>
Tcl
Assuming the presence of the Tcl solutions to these tasks: Element-wise operations, Matrix multiplication, Matrix transposition
<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; 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]
}</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