QR decomposition: Difference between revisions
→{{header|Phix}}
(Scala contribution added.) |
|||
Line 1,897:
{{works with|PARI/GP|2.6.0 and above}}
<lang parigp>matqr(M)</lang>
=={{header|Phix}}==
using matrix_mul from [[Matrix_multiplication#Phix]]
<lang Phix>-- demo/rosettacode/QRdecomposition.exw
function vtranspose(sequence v)
-- transpose a vector of length m into an mx1 matrix,
-- eg {1,2,3} -> {{1},{2},{3}}
for i=1 to length(v) do v[i] = {v[i]} end for
return v
end function
function mat_col(sequence a, integer col)
sequence res = repeat(0,length(a))
for i=col to length(a) do
res[i] = a[i,col]
end for
return res
end function
function mat_norm(sequence a)
atom res = 0
for i=1 to length(a) do
res += a[i]*a[i]
end for
res = sqrt(res)
return res
end function
function mat_ident(integer n)
sequence res = repeat(repeat(0,n),n)
for i=1 to n do
res[i,i] = 1
end for
return res
end function
function QRHouseholder(sequence a)
integer columns = length(a[1]),
rows = length(a),
m = max(columns,rows),
n = min(rows,columns)
sequence q, I = mat_ident(m), Q = I, u, v
for j=1 to min(m-1,n) do
u = mat_col(a,j)
u[j] -= mat_norm(u)
v = sq_div(u,mat_norm(u))
q = sq_sub(I,sq_mul(2,matrix_mul(vtranspose(v),{v})))
a = matrix_mul(q,a)
Q = matrix_mul(Q,q)
end for
-- Get the upper triangular matrix R.
sequence R = repeat(repeat(0,n),m)
for i=1 to n do
for j=i to n do
R[i,j] = a[i,j]
end for
end for
return {Q,R}
end function
sequence a = {{12, -51, 4},
{ 6, 167, -68},
{-4, 24, -41}},
{q,r} = QRHouseholder(a)
?"A" pp(a,{pp_Nest,1})
?"Q" pp(q,{pp_Nest,1})
?"R" pp(r,{pp_Nest,1})
?"Q * R" pp(matrix_mul(q,r),{pp_Nest,1})</lang>
{{out}}
<pre>
"A"
{{12,-51,4},
{6,167,-68},
{-4,24,-41}}
"Q"
{{0.8571428571,-0.3942857143,0.3314285714},
{0.4285714286,0.9028571429,-0.03428571429},
{-0.2857142857,0.1714285714,0.9428571429}}
"R"
{{14,21,-14},
{0,175,-70},
{0,0,-35}}
"Q * R"
{{12,-51,4},
{6,167,-68},
{-4,24,-41}}
</pre>
using matrix_transpose from [[Matrix_transposition#Phix]]
<lang Phix>procedure least_squares()
sequence x = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10},
y = {1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321},
a = repeat(repeat(0,3),length(x))
for i=1 to length(x) do
for j=1 to 3 do
a[i,j] = power(x[i],j-1)
end for
end for
{q,r} = QRHouseholder(a)
sequence t = matrix_transpose(q),
b = matrix_mul(t,vtranspose(y)),
z = repeat(0,3)
for k=3 to 1 by -1 do
atom s = 0
if k<3 then
for j = k+1 to 3 do
s += r[k,j]*z[j]
end for
end if
z[k] = (b[k][1]-s)/r[k,k]
end for
?{"Least-squares solution:",z}
end procedure
least_squares()</lang>
{{out}}
<pre>
{"Least-squares solution:",{1.0,2.0,3.0}}
</pre>
=={{header|Python}}==
| |||