QR decomposition: Difference between revisions
Content added Content deleted
Line 2,666: | Line 2,666: | ||
See [[QR_decomposition#Axiom]] in Axiom. |
See [[QR_decomposition#Axiom]] in Axiom. |
||
=={{header|Standard ML}}== |
|||
We first define a signature for a radical category joined with a field and then a functor: |
|||
<lang sml>signature RADCATFIELD = sig |
|||
type real |
|||
val zero : real |
|||
val one : real |
|||
val + : real * real -> real |
|||
val - : real * real -> real |
|||
val * : real * real -> real |
|||
val / : real * real -> real |
|||
val toString : real -> string |
|||
val sign : real -> real |
|||
val sqrt : real -> real |
|||
end |
|||
functor QR(F: RADCATFIELD) = struct |
|||
structure A = struct |
|||
local |
|||
open Array |
|||
in |
|||
fun unitVector n = tabulate (n, fn i => if i=0 then F.one else F.zero) |
|||
fun map f x = tabulate(length x, fn i => f(sub(x,i))) |
|||
fun map2 f (x, y) = tabulate(length x, fn i => f(sub(x,i),sub(y,i))) |
|||
val op + = map2 F.+ |
|||
val op - = map2 F.- |
|||
val op * = map2 F.* |
|||
fun multc(c,x) = array(length x,c)*x |
|||
fun dot (x,y) = foldl F.+ F.zero (x*y) |
|||
fun outer f (x,y) = |
|||
Array2.tabulate Array2.RowMajor (length x, length y, |
|||
fn (i,j) => f(sub(x,i),sub(y,j))) |
|||
fun copy x = map (fn x => x) x |
|||
fun fromVector v = tabulate(Vector.length v, fn i => Vector.sub(v,i)) |
|||
fun slice(x,i,sz) = |
|||
let open ArraySlice |
|||
val s = slice(x,i,sz) |
|||
in Array.tabulate(length s, fn i => sub(s,i)) end |
|||
end |
|||
end |
|||
structure M = struct |
|||
local |
|||
open Array2 |
|||
in |
|||
fun map f x = tabulate RowMajor (nRows x, nCols x, fn (i,j) => f(sub(x,i,j))) |
|||
fun map2 f (x, y) = |
|||
tabulate RowMajor (nRows x, nCols x, fn (i,j) => f(sub(x,i,j),sub(y,i,j))) |
|||
fun scalarMatrix(m, x) = tabulate RowMajor (m,m,fn (i,j) => if i=j then x else F.zero) |
|||
fun multc(c, x) = map (fn xij => F.*(c,xij)) x |
|||
val op + = map2 F.+ |
|||
val op - = map2 F.- |
|||
fun column(x,i) = A.fromVector(Array2.column(x,i)) |
|||
fun row(x,i) = A.fromVector(Array2.row(x,i)) |
|||
fun x*y = tabulate RowMajor (nRows x, nCols y, |
|||
fn (i,j) => A.dot(row(x,i), column(y,j))) |
|||
fun multa(x,a) = Array.tabulate (nRows x, fn i => A.dot(row(x,i), a)) |
|||
fun copy x = map (fn x => x) x |
|||
fun subMatrix(h, i1, i2, j1, j2) = |
|||
tabulate RowMajor (Int.+(Int.-(i2,i1),1), |
|||
Int.+(Int.-(j2,j1),1), |
|||
fn (a,b) => sub(h,Int.+(i1,a),Int.+(j1,b))) |
|||
fun transpose m = tabulate RowMajor (nCols m, |
|||
nRows m, |
|||
fn (i,j) => sub(m,j,i)) |
|||
fun updateSubMatrix(h,i,j,s) = |
|||
tabulate RowMajor (nRows s, nCols s, fn (a,b) => update(h,Int.+(i,a),Int.+(j,b),sub(s,a,b))) |
|||
fun print a = |
|||
let |
|||
val (m,n) = dimensions a |
|||
fun loop(i,j,x) = |
|||
(if i=0 andalso j=0 then TextIO.print "Array2.fromList([" else (); |
|||
if j=0 andalso Int.>(i,0) then TextIO.print "," else (); |
|||
if j=0 then TextIO.print "[" else (); |
|||
if Int.>(j,0) then TextIO.print "," else (); |
|||
(TextIO.print o F.toString) x; |
|||
if j=Int.-(n,1) then TextIO.print "]" else (); |
|||
if i=Int.-(m,1) andalso j=Int.-(n,1) then TextIO.print "])\n" else ()) |
|||
in |
|||
appi RowMajor loop {base=a,nrows=SOME m,ncols=SOME n,row=0,col=0} |
|||
end |
|||
end |
|||
end |
|||
fun householder a = |
|||
let open Array |
|||
val m = length a |
|||
val len = F.sqrt(A.dot(a,a)) |
|||
val u = A.+(a, A.multc(F.*(len,F.sign(sub(a,0))), A.unitVector m)) |
|||
val v = A.multc(F./(F.one,sub(u,0)), u) |
|||
val beta = F./(F.+(F.one,F.one),A.dot(v,v)) |
|||
in |
|||
M.-(M.scalarMatrix(m,F.one), M.multc(beta,A.outer F.* (v,v))) |
|||
end |
|||
fun qr mat = |
|||
let open Array2 |
|||
val (m,n) = dimensions mat |
|||
val upperIndex = if m=n then Int.-(n,1) else n |
|||
fun loop(i,qm,rm) = if i=upperIndex then {q=qm,r=rm} else |
|||
let val x = A.slice(A.fromVector(column(rm,i)),i,NONE) |
|||
val h = M.scalarMatrix(m,F.one) |
|||
val _ = M.updateSubMatrix(h,i,i,householder x) |
|||
in |
|||
loop(Int.+(i,1), M.*(qm,h), M.*(h,rm)) |
|||
end |
|||
in |
|||
loop(0, M.scalarMatrix(m,F.one), mat) |
|||
end |
|||
fun solveUpperTriangular(r,b) = |
|||
let open Array |
|||
val n = Array2.nCols r |
|||
val x = array(n, F.zero) |
|||
fun loop k = |
|||
let val index = Int.min(Int.-(n,1),Int.+(k,1)) |
|||
val _ = update(x,k, |
|||
F./(F.-(sub(b,k), |
|||
A.dot(A.slice(x,index,NONE), |
|||
A.slice(M.row(r,k),index,NONE))), |
|||
Array2.sub(r,k,k))) |
|||
in |
|||
if k=0 then x else loop(Int.-(k,1)) |
|||
end |
|||
in |
|||
loop (Int.-(n,1)) |
|||
end |
|||
fun lsqr(a,b) = |
|||
let val {q,r} = qr a |
|||
val n = Array2.nCols r |
|||
in |
|||
solveUpperTriangular(M.subMatrix(r, 0, Int.-(n,1), 0, Int.-(n,1)), |
|||
M.multa(M.transpose(q), b)) |
|||
end |
|||
fun pow(x,1) = x |
|||
| pow(x,n) = F.*(x,pow(x,Int.-(n,1))) |
|||
fun polyfit(x,y,n) = |
|||
let open Array2 |
|||
val a = tabulate RowMajor (Array.length x, |
|||
Int.+(n,1), |
|||
fn (i,j) => if j=0 then F.one else |
|||
pow(Array.sub(x,i),j)) |
|||
in |
|||
lsqr(a,y) |
|||
end |
|||
end</lang> |
|||
We can then show the examples:<lang sml>structure RealField : RADCATFIELD = struct |
|||
open Real |
|||
val one = 1.0 |
|||
val zero = 0.0 |
|||
val sign = real o Real.sign |
|||
val sqrt = Real.Math.sqrt |
|||
end |
|||
val mat = Array2.fromList [[12.0, ~51.0, 4.0], [6.0, 167.0, ~68.0], [~4.0, 24.0, ~41.0]]; |
|||
structure Q = QR(RealField); |
|||
structure M = Q.M; |
|||
let val {q,r} = Q.qr(mat) in (M.print q; M.print r) end; |
|||
val fit = |
|||
let open Array |
|||
val x = fromList [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0] |
|||
val y = fromList [1.0, 6.0, 17.0, 34.0, 57.0, 86.0, 121.0, 162.0, 209.0, 262.0, 321.0] |
|||
in |
|||
Q.polyfit(x, y, 2) |
|||
end |
|||
Array2.fromList([[~0.857142857143,0.394285714286,0.331428571429], |
|||
[~0.428571428571,~0.902857142857,~0.0342857142857], |
|||
[0.285714285714,~0.171428571429,0.942857142857]]) |
|||
Array2.fromList([[~14.0,~21.0,14.0], |
|||
[5.97812397875E~18,~175.0,70.0], |
|||
[4.47505280695E~16,0.0,~35.0]]) |
|||
val fit = [|1.0,2.0,3.0|] : real array</lang> |
|||
=={{header|Stata}}== |
=={{header|Stata}}== |
||
See [http://www.stata.com/help.cgi?mf_qrd QR decomposition] in Stata help. |
See [http://www.stata.com/help.cgi?mf_qrd QR decomposition] in Stata help. |