QR decomposition: Difference between revisions

Content added Content deleted

Inline

Revision as of 17:50, 19 November 2012

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.

Task

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. The method of Householder reflections should be used:

Method

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}

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
     u := a + length(a)*signValue(a(1))*unitVector(m) 
     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 default QR decomposition in R.

BBC BASIC

Works with: BBC BASIC for Windows

Makes heavy use of BBC BASIC's matrix arithmetic. <lang bbcbasic> *FLOAT 64

     @% = &2040A
     INSTALL @lib$+"ARRAYLIB"
     
     REM Test matrix for QR decomposition:
     DIM A(2,2)
     A() = 12, -51,   4, \
     \      6, 167, -68, \
     \     -4,  24, -41
     
     REM Do the QR decomposition:
     DIM Q(2,2), R(2,2)
     PROCqrdecompose(A(), Q(), R())
     PRINT "Q:"
     PRINT Q(0,0), Q(0,1), Q(0,2)
     PRINT Q(1,0), Q(1,1), Q(1,2)
     PRINT Q(2,0), Q(2,1), Q(2,2)
     PRINT "R:"
     PRINT R(0,0), R(0,1), R(0,2)
     PRINT R(1,0), R(1,1), R(1,2)
     PRINT R(2,0), R(2,1), R(2,2)
     
     REM Test data for least-squares solution:
     DIM x(10) : x() = 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
     DIM y(10) : y() = 1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321
     
     REM Do the least-squares solution:
     DIM a(10,2), q(10,10), r(10,2), t(10,10), b(10), z(2)
     FOR i% = 0 TO 10
       FOR j% = 0 TO 2
         a(i%,j%) = x(i%) ^ j%
       NEXT
     NEXT
     PROCqrdecompose(a(), q(), r())
     PROC_transpose(q(),t())
     b() = t() . y()
     FOR k% = 2 TO 0 STEP -1
       s = 0
       IF k% < 2 THEN
         FOR j% = k%+1 TO 2
           s += r(k%,j%) * z(j%)
         NEXT
       ENDIF
       z(k%) = (b(k%) - s) / r(k%,k%)
     NEXT k%
     PRINT '"Least-squares solution:"
     PRINT z(0), z(1), z(2)
     END
     
     DEF PROCqrdecompose(A(), Q(), R())
     LOCAL i%, k%, m%, n%, H()
     m% = DIM(A(),1) : n% = DIM(A(),2)
     DIM H(m%,m%)
     FOR i% = 0 TO m% : Q(i%,i%) = 1 : NEXT
     WHILE n%
       PROCqrstep(n%, k%, A(), H())
       A() = H() . A()
       Q() = Q() . H()
       k% += 1
       m% -= 1
       n% -= 1
     ENDWHILE
     R() = A()
     ENDPROC
     
     DEF PROCqrstep(n%, k%, A(), H())
     LOCAL a(), h(), i%, j%
     DIM a(n%,0), h(n%,n%)
     FOR i% = 0 TO n% : a(i%,0) = A(i%+k%,k%) : NEXT
     PROChouseholder(h(), a())
     H() = 0  : H(0,0) = 1
     FOR i% = 0 TO n%
       FOR j% = 0 TO n%
         H(i%+k%,j%+k%) = h(i%,j%)
       NEXT
     NEXT
     ENDPROC
     
     REM Create the Householder matrix for the supplied column vector:
     DEF PROChouseholder(H(), a())
     LOCAL e(), u(), v(), vt(), vvt(), I(), d()
     LOCAL i%, n% : n% = DIM(a(),1)
     REM Create the scaled standard basis vector e():
     DIM e(n%,0) : e(0,0) = SGN(a(0,0)) * MOD(a())
     REM Create the normal vector u():
     DIM u(n%,0) : u() = a() + e()
     REM Normalise with respect to the first element:
     DIM v(n%,0) : v() = u() / u(0,0)
     REM Get the transpose of v() and its dot product with v():
     DIM vt(0,n%), d(0) : PROC_transpose(v(), vt()) : d() = vt() . v()
     REM Get the product of v() and vt():
     DIM vvt(n%,n%) : vvt() = v() . vt()
     REM Create an identity matrix I():
     DIM I(n%,n%) : FOR i% = 0 TO n% : I(i%,i%) = 1 : NEXT
     REM Create the Householder matrix H() = I - 2/vt()v() v()vt():
     vvt() *= 2 / d(0) : H() = I() - vvt()
     ENDPROC</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

Least-squares solution:
    1.0000    2.0000    3.0000

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(int n;double a[][n], 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] = 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(" %8.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->n && 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->n && 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}, { -1, 1, 0}, { 2, 0, 3}, };

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

puts("Q"); matrix_show(Q); puts("R"); matrix_show(R);

// to show their product is the input matrix mat m = matrix_mul(Q, R); puts("Q * R"); matrix_show(m);

matrix_delete(x); matrix_delete(R); matrix_delete(Q); matrix_delete(m); return 0; }</lang>

Output:

Q
    0.846   -0.391    0.343    0.082    0.078
    0.423    0.904   -0.029    0.026    0.045
   -0.282    0.170    0.933   -0.047   -0.137
   -0.071    0.014   -0.001    0.980   -0.184
    0.141   -0.017   -0.106   -0.171   -0.969

R
   14.177   20.667  -13.402
   -0.000  175.043  -70.080
    0.000    0.000  -35.202
   -0.000   -0.000   -0.000
    0.000    0.000   -0.000

Q * R
   12.000  -51.000    4.000
    6.000  167.000  -68.000
   -4.000   24.000  -41.000
   -1.000    1.000   -0.000
    2.000   -0.000    3.000

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 /*nothrow*/ if (is(U == T) || is(U == T[][])) {

   static if (is(U == T[][]))
       assert(A.length == B.length);
   if (A.empty)
       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 /*nothrow*/ {

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

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

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

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

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

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

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

}

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

   //return matrix.all!(r => r.length == mat[0].length)();
   foreach (row; mat)
       if (row.length != mat[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] = a[i].dotProduct(aux);
   }
   return result;

}

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[][] m) /*pure nothrow*/ {

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

}

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)(in T[][] a) {

   immutable size_t m = rows(a);
   immutable 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) {

   immutable m = rows(A);
   immutable 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 {

   immutable 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, in T[][] b) {

   const qr = QRdecomposition(A);
   immutable 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 (Q, R) = QRdecomposition(
   const qr = QRdecomposition([[12.0, -51,   4],
                               [ 6.0, 167, -68],
                               [-4.0,  24, -41]]);
   enum string form = "[%([%(%2.3f, %)]%|,\n %)]\n";
   writefln(form, qr.Q);
   writefln(form, qr.R);

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

}</lang>

Output:

[[-0.857, 0.394, 0.331],
 [-0.429, -0.903, -0.034],
 [0.286, -0.171, 0.943]]

[[-14.000, -21.000, 14.000],
 [0.000, -175.000, 70.000],
 [0.000, -0.000, -35.000]]

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

Go

Translation of: Common Lisp

Library: Gomatrix

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:

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

Maxima

<lang maxima>load(lapack)$ /* This may hang up in wxMaxima, if this happens, use xMaxima or plain MAxima in a terminal */

a: matrix([12, -51, 4],

         [ 6, 167, -68],
         [-4,  24, -41])$

[q, r]: dgeqrf(a)$

mat_norm(q . r - a, 1); 4.2632564145606011E-14

/* Note: the lapack package is a lisp translation of the fortran lapack library */</lang>

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>

Rascal

This function applies the Gram Schmidt algorithm. Q is printed in the console, R can be printed or visualized.

<lang Rascal>import util::Math; import Prelude; import vis::Figure; import vis::Render;

public rel[real,real,real] QRdecomposition(rel[real x, real y, real v] matrix){ //orthogonalcolumns oc = domainR(matrix, {0.0}); for (x <- sort(toList(domain(matrix)-{0.0}))){ c = domainR(matrix, {x}); o = domainR(oc, {x-1});

for (n <- [1.0 .. x]){ o = domainR(oc, {n-1}); c = matrixSubtract(c, matrixMultiplybyN(o, matrixDotproduct(o, c)/matrixDotproduct(o, o))); }

oc += c; }

Q = {}; //from orthogonal to orthonormal columns for (el <- oc){ c = domainR(oc, {el[0]}); Q += matrixNormalize({el}, c); }

//from Q to R R= matrixMultiplication(matrixTranspose(Q), matrix); R= {<x,y,toReal(round(v))> | <x,y,v> <- R};

println("Q:"); iprintlnExp(Q); println(); println("R:"); return R; }

//a function that takes the transpose of a matrix, see also Rosetta Code problem "Matrix transposition" public rel[real, real, real] matrixTranspose(rel[real x, real y, real v] matrix){ return {<y, x, v> | <x, y, v> <- matrix}; }

//a function to normalize an element of a matrix by the normalization of a column public rel[real,real,real] matrixNormalize(rel[real x, real y, real v] element, rel[real x, real y, real v] column){ normalized = 1.0/nroot((0.0 | it + v*v | <x,y,v> <- column), 2); return matrixMultiplybyN(element, normalized); }

//a function that takes the dot product, see also Rosetta Code problem "Dot product" public real matrixDotproduct(rel[real x, real y, real v] column1, rel[real x, real y, real v] column2){ return (0.0 | it + v1*v2 | <x1,y1,v1> <- column1, <x2,y2,v2> <- column2, y1==y2); }

//a function to subtract two columns public rel[real,real,real] matrixSubtract(rel[real x, real y, real v] column1, rel[real x, real y, real v] column2){ return {<x1,y1,v1-v2> | <x1,y1,v1> <- column1, <x2,y2,v2> <- column2, y1==y2}; }

//a function to multiply a column by a number public rel[real,real,real] matrixMultiplybyN(rel[real x, real y, real v] column, real n){ return {<x,y,v*n> | <x,y,v> <- column}; }

//a function to perform matrix multiplication, see also Rosetta Code problem "Matrix multiplication". public rel[real, real, real] matrixMultiplication(rel[real x, real y, real v] matrix1, rel[real x, real y, real v] matrix2){ if (max(matrix1.x) == max(matrix2.y)){ p = {<x1,y1,x2,y2, v1*v2> | <x1,y1,v1> <- matrix1, <x2,y2,v2> <- matrix2};

result = {}; for (y <- matrix1.y){ for (x <- matrix2.x){ v = (0.0 | it + v | <x1, y1, x2, y2, v> <- p, x==x2 && y==y1, x1==y2 && y2==x1); result += <x,y,v>; } } return result; } else throw "Matrix sizes do not match."; }

// a function to visualize the result public void displayMatrix(rel[real x, real y, real v] matrix){ points = [box(text("<v>"), align(0.3333*(x+1),0.3333*(y+1)),shrink(0.25)) | <x,y,v> <- matrix]; render(overlay([*points], aspectRatio(1.0))); }

//a matrix, given by a relation of <x-coordinate, y-coordinate, value>. public rel[real x, real y, real v] matrixA = { <0.0,0.0,12.0>, <0.0,1.0, 6.0>, <0.0,2.0,-4.0>, <1.0,0.0,-51.0>, <1.0,1.0,167.0>, <1.0,2.0,24.0>, <2.0,0.0,4.0>, <2.0,1.0,-68.0>, <2.0,2.0,-41.0> };</lang>

Example using visualization

rascal>displayMatrix(QRdecomposition(matrixA))

Q:
{
  <1.0,0.0,-0.394285714285714285714285714285714285714285714285714285714285714285713300>,
  <0.0,0.0,0.857142857142857142857142857142857142857142857142857142857142857142840>,
  <0.0,1.0,0.428571428571428571428571428571428571428571428571428571428571428571420>,
  <0.0,2.0,-0.285714285714285714285714285714285714285714285714285714285714285714280>,
  <2.0,0.0,-0.33142857142857142857142857142857142857142857142857142857142857142858800>,
  <1.0,2.0,0.171428571428571428571428571428571428571428571428571428571428571428571000>,
  <2.0,2.0,-0.94285714285714285714285714285714285714285714285714285714285714285719000>,
  <1.0,1.0,0.902857142857142857142857142857142857142857142857142857142857142857140600>,
  <2.0,1.0,0.03428571428571428571428571428571428571428571428571428571428571428571600>
}
See R in picture

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; 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