Rosetta Code

QR decomposition: Difference between revisions

Page Discussion

QR decomposition (view source)

Revision as of 20:12, 26 March 2022

33,989 bytes added ,  2 years ago
m
→‎{{header|Phix}}: syntax coloured, improved output formatting a touch
(Realize in F#)
m (→‎{{header|Phix}}: syntax coloured, improved output formatting a touch)
Line 2,413:
 
=={{header|Phix}}==
using matrix_mul() from [[Matrix_multiplication#Phix]]
and matrix_transpose() from [[Matrix_transposition#Phix]]
<lang Phix>-- demo/rosettacode/QRdecomposition.exw
<!--<lang Phix>(phixonline)-->
function vtranspose(sequence v)
<span style="color: #000080;font-style:italic;">-- demo/rosettacode/QRdecomposition.exw</span>
-- transpose a vector of length m into an mx1 matrix,
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
-- eg {1,2,3} -> {{1},{2},{3}}
for i=1 to length(v) do v[i] = {v[i]} end for
<span style="color: #008080;">function</span> <span style="color: #000000;">matrix_mul</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">)</span>
return v
<span style="color: #004080;">integer</span> <span style="color: #000000;">arows</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">~</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">acols</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">~</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span>
end function
<span style="color: #000000;">brows</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">~</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">bcols</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">~</span><span style="color: #000000;">b</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span>
 
<span style="color: #008080;">if</span> <span style="color: #000000;">acols</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">brows</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">0</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
function mat_col(sequence a, integer col)
<span style="color: #004080;">sequence</span> <span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">bcols</span><span style="color: #0000FF;">),</span><span style="color: #000000;">arows</span><span style="color: #0000FF;">)</span>
sequence res = repeat(0,length(a))
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">arows</span> <span style="color: #008080;">do</span>
for i=col to length(a) do
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">bcols</span> <span style="color: #008080;">do</span>
res[i] = a[i,col]
<span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">acols</span> <span style="color: #008080;">do</span>
end for
<span style="color: #000000;">c</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">b</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span>
return res
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
end function
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
 
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
function mat_norm(sequence a)
<span style="color: #008080;">return</span> <span style="color: #000000;">c</span>
atom res = 0
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
for i=1 to length(a) do
res += a[i]*a[i]
<span style="color: #008080;">function</span> <span style="color: #000000;">vtranspose</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span>
end for
<span style="color: #000080;font-style:italic;">-- transpose a vector of length m into an mx1 matrix,
res = sqrt(res)
-- eg {1,2,3} -&gt; <nowiki>{{</nowiki>1},{2},{3<nowiki>}}</nowiki></span>
return res
<span style="color: #004080;">integer</span> <span style="color: #000000;">l</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">v</span><span style="color: #0000FF;">)</span>
end function
<span style="color: #004080;">sequence</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">l</span><span style="color: #0000FF;">)</span>
 
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">l</span> <span style="color: #008080;">do</span> <span style="color: #000000;">res</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">v</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]}</span> <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
function mat_ident(integer n)
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span>
sequence res = repeat(repeat(0,n),n)
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
for i=1 to n do
res[i,i] = 1
<span style="color: #008080;">function</span> <span style="color: #000000;">mat_col</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">col</span><span style="color: #0000FF;">)</span>
end for
<span style="color: #004080;">integer</span> <span style="color: #000000;">la</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
return res
<span style="color: #004080;">sequence</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">la</span><span style="color: #0000FF;">)</span>
end function
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">col</span> <span style="color: #008080;">to</span> <span style="color: #000000;">la</span> <span style="color: #008080;">do</span>
 
<span style="color: #000000;">res</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">col</span><span style="color: #0000FF;">]</span>
function QRHouseholder(sequence a)
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
integer columns = length(a[1]),
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span>
rows = length(a),
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
m = max(columns,rows),
n = min(rows,columns)
<span style="color: #008080;">function</span> <span style="color: #000000;">mat_norm</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
sequence q, I = mat_ident(m), Q = I, u, v
<span style="color: #004080;">atom</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
 
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
--
<span style="color: #000000;">res</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span>
-- Programming note: The code of this main loop was not as easily
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
-- written as the first glance might suggest. Explicitly setting
<span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">)</span>
-- to 0 any a[i,j] [etc] that should be 0 but have inadvertently
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span>
-- gotten set to +/-1e-15 or thereabouts may be advisable. The
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
-- commented-out code was retrieved from a backup and should be
-- treated as an example and not be trusted (iirc, it made no
<span style="color: #008080;">function</span> <span style="color: #000000;">mat_ident</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
-- difference to the test cases used, so I deleted it, and then
<span style="color: #004080;">sequence</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">),</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
-- had second thoughts a few days later).
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span>
--
<span style="color: #000000;">res</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span>
for j=1 to min(m-1,n) do
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
u = mat_col(a,j)
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span>
u[j] -= mat_norm(u)
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
v = sq_div(u,mat_norm(u))
q = sq_sub(I,sq_mul(2,matrix_mul(vtranspose(v),{v})))
<span style="color: #008080;">function</span> <span style="color: #000000;">QRHouseholder</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
a = matrix_mul(q,a)
<span style="color: #004080;">integer</span> <span style="color: #000000;">cols</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]),</span>
-- for row=j+1 to length(a) do
<span style="color: #000000;">rows</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">),</span>
-- a[row][j] = 0
<span style="color: #000000;">m</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">max</span><span style="color: #0000FF;">(</span><span style="color: #000000;">cols</span><span style="color: #0000FF;">,</span><span style="color: #000000;">rows</span><span style="color: #0000FF;">),</span>
-- end for
<span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">min</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rows</span><span style="color: #0000FF;">,</span><span style="color: #000000;">cols</span><span style="color: #0000FF;">)</span>
Q = matrix_mul(Q,q)
<span style="color: #004080;">sequence</span> <span style="color: #000000;">q</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">I</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">mat_ident</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">Q</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">I</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">u</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">v</span>
end for
<span style="color: #000080;font-style:italic;">--
 
-- Programming note: The code of this main loop was not as easily
-- Get the upper triangular matrix R.
-- written as the first glance might suggest. Explicitly setting
sequence R = repeat(repeat(0,n),m)
-- to 0 any a[i,j] [etc] that should be 0 but have inadvertently
for i=1 to n do -- (logically 1 to m(>=n), but no need)
-- gotten set to +/-1e-15 or thereabouts may be advisable. The
for j=i to n do
-- commented-out code was retrieved from a backup and should be
R[i,j] = a[i,j]
-- treated as an example and not be trusted (iirc, it made no
end for
-- difference to the test cases used, so I deleted it, and then
end for
-- had second thoughts about it a few days later).
--</span>
return {Q,R}
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">min</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
end function
<span style="color: #000000;">u</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">mat_col</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">j</span><span style="color: #0000FF;">)</span>
 
<span style="color: #000000;">u</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">-=</span> <span style="color: #000000;">mat_norm</span><span style="color: #0000FF;">(</span><span style="color: #000000;">u</span><span style="color: #0000FF;">)</span>
sequence a = {{12, -51, 4},
<span style="color: #000000;">v</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_div</span><span style="color: #0000FF;">(</span><span style="color: #000000;">u</span><span style="color: #0000FF;">,</span><span style="color: #000000;">mat_norm</span><span style="color: #0000FF;">(</span><span style="color: #000000;">u</span><span style="color: #0000FF;">))</span>
{ 6, 167, -68},
<span style="color: #000000;">q</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">I</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">matrix_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">vtranspose</span><span style="color: #0000FF;">(</span><span style="color: #000000;">v</span><span style="color: #0000FF;">),{</span><span style="color: #000000;">v</span><span style="color: #0000FF;">})))</span>
{-4, 24, -41}},
<span style="color: #000000;">a</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">matrix_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">q</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
{q,r} = QRHouseholder(a)
<span style="color: #000080;font-style:italic;">-- for row=j+1 to length(a) do
 
?"A" -- pp(a,{pp_Nest,1})[row][j] = 0
?"Q" -- end pp(q,{pp_Nest,1})for</span>
<span style="color: #000000;">Q</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">matrix_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">Q</span><span style="color: #0000FF;">,</span><span style="color: #000000;">q</span><span style="color: #0000FF;">)</span>
?"R" pp(r,{pp_Nest,1})
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
?"Q * R" pp(matrix_mul(q,r),{pp_Nest,1})</lang>
<span style="color: #000080;font-style:italic;">-- Get the upper triangular matrix R.</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">R</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">),</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span> <span style="color: #000080;font-style:italic;">-- (logically 1 to m(&gt;=n), but no need)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">i</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">R</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">Q</span><span style="color: #0000FF;">,</span><span style="color: #000000;">R</span><span style="color: #0000FF;">}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">a</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{</span><span style="color: #000000;">12</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">51</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #000000;">6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">167</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">68</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{-</span><span style="color: #000000;">4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">24</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">41</span><span style="color: #0000FF;">}}</span>
<span style="color: #004080;">sequence</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">q</span><span style="color: #0000FF;">,</span><span style="color: #000000;">r</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">QRHouseholder</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">ppOpt</span><span style="color: #0000FF;">({</span><span style="color: #004600;">pp_Nest</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #004600;">pp_IntFmt</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%4d"</span><span style="color: #0000FF;">,</span><span style="color: #004600;">pp_FltFmt</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%4g"</span><span style="color: #0000FF;">,</span><span style="color: #004600;">pp_IntCh</span><span style="color: #0000FF;">,</span><span style="color: #004600;">false</span><span style="color: #0000FF;">})</span>
<span style="color: #0000FF;">?</span><span style="color: #008000;">"A"</span> <span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
<span style="color: #0000FF;">?</span><span style="color: #008000;">"Q"</span> <span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">q</span><span style="color: #0000FF;">)</span>
<span style="color: #0000FF;">?</span><span style="color: #008000;">"R"</span> <span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">)</span>
<span style="color: #0000FF;">?</span><span style="color: #008000;">"Q * R"</span> <span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">matrix_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">q</span><span style="color: #0000FF;">,</span><span style="color: #000000;">r</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">matrix_transpose</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">mat</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">rows</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mat</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">cols</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mat</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">])</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">rows</span><span style="color: #0000FF;">),</span><span style="color: #000000;">cols</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">rows</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">c</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">cols</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">res</span><span style="color: #0000FF;">[</span><span style="color: #000000;">c</span><span style="color: #0000FF;">][</span><span style="color: #000000;">r</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">mat</span><span style="color: #0000FF;">[</span><span style="color: #000000;">r</span><span style="color: #0000FF;">][</span><span style="color: #000000;">c</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #000080;font-style:italic;">--?"Q * Q'" pp(matrix_mul(q,matrix_transpose(q))) -- (~1e-16s)</span>
<span style="color: #0000FF;">?</span><span style="color: #008000;">"Q * Q`"</span> <span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_round</span><span style="color: #0000FF;">(</span><span style="color: #000000;">matrix_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">q</span><span style="color: #0000FF;">,</span><span style="color: #000000;">matrix_transpose</span><span style="color: #0000FF;">(</span><span style="color: #000000;">q</span><span style="color: #0000FF;">)),</span><span style="color: #000000;">1e15</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">least_squares</span><span style="color: #0000FF;">()</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">x</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">10</span><span style="color: #0000FF;">},</span>
<span style="color: #000000;">y</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">17</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">34</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">57</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">86</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">121</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">162</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">209</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">262</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">321</span><span style="color: #0000FF;">},</span>
<span style="color: #000000;">a</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">3</span><span style="color: #0000FF;">),</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">3</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span><span style="color: #000000;">j</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">q</span><span style="color: #0000FF;">,</span><span style="color: #000000;">r</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">QRHouseholder</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">matrix_transpose</span><span style="color: #0000FF;">(</span><span style="color: #000000;">q</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">b</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">matrix_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">vtranspose</span><span style="color: #0000FF;">(</span><span style="color: #000000;">y</span><span style="color: #0000FF;">)),</span>
<span style="color: #000000;">z</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">3</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">3</span> <span style="color: #008080;">to</span> <span style="color: #000000;">1</span> <span style="color: #008080;">by</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">k</span><span style="color: #0000FF;"><</span><span style="color: #000000;">3</span> <span style="color: #008080;">then</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">3</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">s</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">z</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000000;">z</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">][</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]-</span><span style="color: #000000;">s</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">r</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Least-squares solution:\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">-- printf(1," %v\n",{z}) -- {1.0,2.0.3,0}
-- printf(1," %v\n",{sq_sub(z,{1,2,3})}) -- (+/- ~1e-14s)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" %v\n"</span><span style="color: #0000FF;">,{</span><span style="color: #7060A8;">sq_round</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1e13</span><span style="color: #0000FF;">)})</span> <span style="color: #000080;font-style:italic;">-- {1,2,3}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #000000;">least_squares</span><span style="color: #0000FF;">()</span>
<!--</lang>-->
{{out}}
<pre>
"A"
{{12  12, -51,4   4},
 {6   6,167 167, -68},
 {  -4,24  24, -41}}
"Q"
{{0.857142857185714,-0.39428571433943,0.331428571433143},
 {0.428571428642857,0.902857142990286,-0.034285714290343},
 {-0.28571428572857,0.171428571417143,0.942857142994286}}
"R"
{{14  14,21  21, -14},
 {0   0,175 175, -70},
 {0   0,0   0, -35}}
"Q Q * R R"
{{12  12, -51,4   4},
 {6   6,167 167, -68},
 {  -4,24  24, -41}}
"Q * Q`"
</pre>
{{   1,   0,   0},
using matrix_transpose from [[Matrix_transposition#Phix]]
 {   0,   1,   0},
<lang Phix>procedure least_squares()
 {   0,   0,   1}}
sequence x = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10},
Least-squares solution:
y = {1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321},
 {1,2,3}
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>
 
Petelomax
7,818

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/189739"