QR decomposition: Difference between revisions
Simplify code; add pivots
(Changed householder: Vector(R)->Matrix(R)) |
(Simplify code; add pivots) |
||
Line 205:
=={{header|Axiom}}==
{{trans|Common Lisp}}
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
Line 246:
(m,n) := (nrows a, ncols a)
qm := scalarMatrix(m,1)
rm
for i in 1..(if m=n then n-1 else n) repeat
x := column(subMatrix(rm,i,m,i,
h := scalarMatrix(m,1)
setsubMatrix!(h,i,i,householder x)
Line 257:
dc := qr a
n := ncols(dc.r)
solveUpperTriangular(subMatrix(dc.r,1,n,1,n),
polyfit(x,y,n) ==
a := new(#x,n+1,0@R)$Matrix(R)
Line 270 ⟶ 269:
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
Line 291 ⟶ 290:
[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. This 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>
=={{header|C}}==
| |||