QR decomposition: Difference between revisions
Content added Content deleted
(Scala contribution added.) |
|||
Line 1,897: | Line 1,897: | ||
{{works with|PARI/GP|2.6.0 and above}} |
{{works with|PARI/GP|2.6.0 and above}} |
||
<lang parigp>matqr(M)</lang> |
<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}}== |
=={{header|Python}}== |