aparler/RosettaCodeData
1
0
Fork 0
Code Issues 12 Pull requests 3 Projects Releases Packages Wiki Activity Actions

tasks a-s

Browse source
This commit is contained in:
Ingy döt Net 2013-04-10 23:57:08 -07:00
parent 47bf37c096
commit b83f433714
12433 changed files with 156208 additions and 123 deletions
Download patch file Download diff file Expand all files Collapse all files

95
Task/QR-decomposition/0DESCRIPTION Normal file
View file

@ -0,0 +1,95 @@
Any rectangular <math>m \times n</math> matrix <math>\mathit A</math> can be decomposed to a product of a orthogonal matrix <math>\mathit Q</math> and a upper (right) triangular matrix <math>\mathit R</math>, as described in [[wp:QR decomposition|QR decomposition]].
'''Task'''
Demonstrate the QR decomposition on the example matrix from the [[wp:QR_decomposition#Example_2|Wikipedia article]]:
::<math>A = \begin{pmatrix}
12 & -51 & 4 \\
6 & 167 & -68 \\
-4 & 24 & -41 \end{pmatrix}</math>
and the usage for linear least squares problems on the example from [[Polynomial_regression]]. The method of [[wp: Householder transformation|Householder reflections]] should be used:
'''Method'''
Multiplying a given vector <math>\mathit a</math>, for example the first column of matrix <math>\mathit A</math>, with the Householder matrix <math>\mathit H</math>, which is given as
::<math>H = I - \frac {2} {u^T u} u u^T</math>
reflects <math>\mathit a</math> about a plane given by its normal vector <math>\mathit u</math>. When the normal vector of the plane <math>\mathit u</math> is given as
::<math>u = a - \|a\|_2 \; e_1</math>
then the transformation reflects <math>\mathit a</math> onto the first standard basis vector
::<math>e_1 = [1 \; 0 \; 0 \; ...]^T</math>
which means that all entries but the first become zero. To avoid numerical cancellation errors, we should take the opposite sign of <math>a_1</math>:
::<math>u = a + \textrm{sign}(a_1)\|a\|_2 \; e_1</math>
and normalize with respect to the first element:
::<math>v = \frac{u}{u_1}</math>
The equation for <math>H</math> thus becomes:
::<math>H = I - \frac {2} {v^T v} v v^T</math>
or, in another form
::<math>H = I - \beta v v^T</math>
with
::<math>\beta = \frac {2} {v^T v}</math>
Applying <math>\mathit H</math> on <math>\mathit a</math> then gives
::<math>H \; a = -\textrm{sign}(a_1) \; \|a\|_2 \; e_1</math>
and applying <math>\mathit H</math> on the matrix <math>\mathit A</math> zeroes all subdiagonal elements of the first column:
::<math>H_1 \; A = \begin{pmatrix}
r_{11} & r_{12} & r_{13} \\
0 & * & * \\
0 & * & * \end{pmatrix}</math>
In the second step, the second column of <math>\mathit A</math>, we want to zero all elements but the first two, which means that we have to calculate <math>\mathit H</math> with the first column of the ''submatrix'' (denoted *), not on the whole second column of <math>\mathit A</math>.
To get <math>H_2</math>, we then embed the new <math>\mathit H</math> into an <math>m \times n</math> identity:
::<math>H_2 = \begin{pmatrix}
1 & 0 & 0 \\
0 & H & \\
0 & & \end{pmatrix}</math>
This is how we can, column by column, remove all subdiagonal elements of <math>\mathit A</math> and thus transform it into <math>\mathit R</math>.
::<math>H_n \; ... \; H_3 H_2 H_1 A = R</math>
The product of all the Householder matrices <math>\mathit H</math>, for every column, in reverse order, will then yield the orthogonal matrix <math>\mathit Q</math>.
::<math>H_1 H_2 H_3 \; ... \; H_n = Q</math>
The QR decomposition should then be used to solve linear least squares ([[Multiple regression]]) problems <math>\mathit A x = b</math> by solving
::<math>R \; x = Q^T \; b</math>
When <math>\mathit R</math> is not square, i.e. <math>m > n</math> we have to cut off the <math>\mathit m - n</math> zero padded bottom rows.
::<math>R =
\begin{pmatrix}
R_1 \\
0 \end{pmatrix}</math>
and the same for the RHS:
::<math>Q^T \; b =
\begin{pmatrix}
q_1 \\
q_2 \end{pmatrix}</math>
Finally, solve the square upper triangular system by back substitution:
::<math>R_1 \; x = q_1</math>

4
Task/QR-decomposition/1META.yaml Normal file
View file

@ -0,0 +1,4 @@
---
category:
- Mathematics
note: Matrices

94
Task/QR-decomposition/Ada/qr-decomposition.ada Normal file
View file

@ -0,0 +1,94 @@
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;

55
Task/QR-decomposition/Axiom/qr-decomposition-1.axiom Normal file
View file

@ -0,0 +1,55 @@
)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)

5
Task/QR-decomposition/Axiom/qr-decomposition-2.axiom Normal file
View file

@ -0,0 +1,5 @@
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)

20
Task/QR-decomposition/Axiom/qr-decomposition-3.axiom Normal file
View file

@ -0,0 +1,20 @@
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)

97
Task/QR-decomposition/BBC-BASIC/qr-decomposition.bbc Normal file
View file

@ -0,0 +1,97 @@
*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

192
Task/QR-decomposition/C/qr-decomposition.c Normal file
View file

@ -0,0 +1,192 @@
#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;
}

54
Task/QR-decomposition/Common-Lisp/qr-decomposition-1.lisp Normal file
View file

@ -0,0 +1,54 @@
(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))

34
Task/QR-decomposition/Common-Lisp/qr-decomposition-2.lisp Normal file
View file

@ -0,0 +1,34 @@
(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)))

9
Task/QR-decomposition/Common-Lisp/qr-decomposition-3.lisp Normal file
View file

@ -0,0 +1,9 @@
(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))

29
Task/QR-decomposition/Common-Lisp/qr-decomposition-4.lisp Normal file
View file

@ -0,0 +1,29 @@
(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))

6
Task/QR-decomposition/Common-Lisp/qr-decomposition-5.lisp Normal file
View file

@ -0,0 +1,6 @@
;; 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))

205
Task/QR-decomposition/D/qr-decomposition.d Normal file
View file

@ -0,0 +1,205 @@
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 (immutable r, const 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 mat.all!(r => r.length == mat[0].length);
}
T[][] matMul(T)(in T[][] a, in T[][] b) /*pure nothrow*/
in {
assert(a.isRectangular && b.isRectangular &&
a[0].length == b.length);
} body {
auto result = new T[][](a.length, b[0].length);
auto aux = new T[b.length];
foreach (immutable j; 0 .. b[0].length) {
foreach (immutable k; 0 .. b.length)
aux[k] = b[k][j];
foreach (immutable 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 (immutable nr, const row; m)
foreach (immutable nc, immutable 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 (immutable 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 (immutable 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 (immutable i, const arow; A)
C[i][] = arow[]; // some wasted copies
foreach (immutable i, const 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]);
immutable e = makeUnitVector!T(m);
immutable u = msum(a, pmul(e, norm(a) * s));
immutable v = pdiv(u, u[0][0]);
immutable 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);
auto Q = matId!T(m);
// Work on n columns of A.
foreach (immutable i; 0 .. (m == n ? n-1 : n)) {
// Select the i-th submatrix. For i=0 this means the original
// matrix A.
immutable B = slice2D(A, i, m-1, i, n-1);
// Take the first column of the current submatrix B.
immutable x = mcol(B, 0);
// Create the Householder matrix for the column and embed it
// into an mxm identity.
immutable 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 (immutable k; 0 .. n) {
T tot = 0;
foreach (immutable 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) {
immutable size_t m = cols(x);
auto A = new Unqual!T[][](m, n + 1);
foreach (immutable i, row; A)
foreach (immutable 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]]);
immutable 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));
}

106
Task/QR-decomposition/Go/qr-decomposition.go Normal file
View file

@ -0,0 +1,106 @@
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
}

16
Task/QR-decomposition/J/qr-decomposition-1.j Normal file
View file

@ -0,0 +1,16 @@
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.
)

6
Task/QR-decomposition/J/qr-decomposition-2.j Normal file
View file

@ -0,0 +1,6 @@
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│
└────────────────────┴──────────┘

4
Task/QR-decomposition/J/qr-decomposition-3.j Normal file
View file

@ -0,0 +1,4 @@
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

4
Task/QR-decomposition/MATLAB/qr-decomposition.m Normal file
View file

@ -0,0 +1,4 @@
A = [12 -51 4
6 167 -68
-4 24 -41];
[Q,R]=qr(A)

11
Task/QR-decomposition/Mathematica/qr-decomposition.math Normal file
View file

@ -0,0 +1,11 @@
{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

12
Task/QR-decomposition/Maxima/qr-decomposition-1.maxima Normal file
View file

@ -0,0 +1,12 @@
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 */

44
Task/QR-decomposition/Maxima/qr-decomposition-2.maxima Normal file
View file

@ -0,0 +1,44 @@
load("linearalgebra")$
load("eigen")$
unitVector(n) := ematrix(n,1,1,1,1);
signValue(r) := block([s:sign(r)],
if s='pos then 1 else if s='zero then 0 else -1);
householder(a) := block([m : length(a),u,v,beta],
u : a + sqrt(a . a)*signValue(a[1,1])*unitVector(m),
v : u / u[1,1],
beta : 2/(v . v),
diagmatrix(m,1) - beta*transpose(v . transpose(v)));
getSubmatrix(obj,i1,j1,i2,j2) :=
genmatrix(lambda([i,j], obj[i+i1-1,j+j1-1]),i2-i1+1,j2-j1+1);
setSubmatrix(obj,i1,j1,subobj) := block([m,n],
[m,n] : matrix_size(subobj),
for i: 0 thru m-1 do
(for j: 0 thru n-1 do
obj[i1+i,j1+j] : subobj[i+1,j+1]));
qr(obj) := block([m,n,qm,rm,i],
[m,n] : matrix_size(obj),
qm : diagmatrix(m,1),
rm : copymatrix(obj),
for i: 1 thru (if m=n then n-1 else n) do
block([x,h],
x : getSubmatrix(rm,i,i,m,i),
h : diagmatrix(m,1),
setSubmatrix(h,i,i,householder(x)),
qm : qm . h,
rm : h . rm),
[qm,rm]);
solveUpperTriangular(r,b) := block([n,x,index,k],
n : second(matrix_size(r)),
x : genmatrix(lambda([a, b], 0), n, 1),
for k: n thru 1 step -1 do
(index : min(n,k+1),
x[k,1] : (b[k,1] - (getSubmatrix(r,k,index,k,n) . getSubmatrix(x,index,1,n,1)))/r[k,k]),
x);
lsqr(a,b) := block([q,r,n],
[q,r] : qr(a),
n : second(matrix_size(r)),
solveUpperTriangular(getSubmatrix(r,1,1,n,n), transpose(q) . b));
polyfit(x,y,n) := block([a,j],
a : genmatrix(lambda([i,j], if j=1 then 1.0b0 else bfloat(x[i,1]^(j-1))),
length(x),n+1),
lsqr(a,y));

29
Task/QR-decomposition/Maxima/qr-decomposition-3.maxima Normal file
View file

@ -0,0 +1,29 @@
(%i) [q,r] : qr(a);
[ 6 69 58 ]
[ - - --- --- ]
[ 7 175 175 ]
[ ] [ - 14 - 21 14 ]
[ 3 158 6 ] [ ]
(%o) [[ - - - --- - --- ], [ 0 - 175 70 ]]
[ 7 175 175 ] [ ]
[ ] [ 0 0 - 35 ]
[ 2 6 33 ]
[ - - -- -- ]
[ 7 35 35 ]
(%i) mat_norm(q . r - a, 1);
(%o) 0
(%i) x : transpose(matrix([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]))$
(%i) y : transpose(matrix([1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]))$
(%i) fpprec : 30$
(%i) polyfit(x, y, 2);
[ 9.99999999999999999999999999996b-1 ]
[ ]
(%o) [ 2.00000000000000000000000000002b0 ]
[ ]
[ 3.0b0 ]

1
Task/QR-decomposition/PARI-GP/qr-decomposition.pari Normal file
View file

@ -0,0 +1 @@
matqr(M)

23
Task/QR-decomposition/R/qr-decomposition.r Normal file
View file

@ -0,0 +1,23 @@
# 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)

93
Task/QR-decomposition/Rascal/qr-decomposition.rascal Normal file
View file

@ -0,0 +1,93 @@
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>
};

32
Task/QR-decomposition/SAS/qr-decomposition.sas Normal file
View file

@ -0,0 +1,32 @@
/* See http://support.sas.com/documentation/cdl/en/imlug/63541/HTML/default/viewer.htm#imlug_langref_sect229.htm */
proc iml;
a={12 -51 4,6 167 -68,-4 24 -41};
print(a);
call qr(q,r,p,d,a);
print(q);
print(r);
quit;
/*
a
12 -51 4
6 167 -68
-4 24 -41
q
-0.857143 0.3942857 -0.331429
-0.428571 -0.902857 0.0342857
0.2857143 -0.171429 -0.942857
r
-14 -21 14
0 -175 70
0 0 35
*/

58
Task/QR-decomposition/Tcl/qr-decomposition-1.tcl Normal file
View file

@ -0,0 +1,58 @@
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]
}

5
Task/QR-decomposition/Tcl/qr-decomposition-2.tcl Normal file
View file

@ -0,0 +1,5 @@
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"