QR decomposition: Difference between revisions
→{{header|jq}}
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
|||
Line 2,101:
0,000 -175,000 70,000
0,000 0,000 35,000</pre>
=={{header|jq}}==
'''Adapted from [[#Wren|Wren]]'''
{{works with|jq}}
'''Also works with gojq, the Go implementation of jq.'''
'''General utilities'''
<syntaxhighlight lang=jq>
def sum(s): reduce s as $_ (0; . + $_);
# Sum of squares
def ss(s): sum(s|.*.);
# Create an m x n matrix
def matrix(m; n; init):
if m == 0 then []
elif m == 1 then [range(0;n) | init]
elif m > 0 then
matrix(1;n;init) as $row
| [range(0;m) | $row ]
else error("matrix\(m);_;_) invalid")
end;
def dot_product(a; b):
reduce range(0;a|length) as $i (0; . + (a[$i] * b[$i]) );
# A and B should both be numeric matrices, A being m by n, and B being n by p.
def multiply($A; $B):
($B[0]|length) as $p
| ($B|transpose) as $BT
| reduce range(0; $A|length) as $i
([];
reduce range(0; $p) as $j
(.;
.[$i][$j] = dot_product( $A[$i]; $BT[$j] ) ));
# $ndec decimal places
def round($ndec):
def rpad: tostring | ($ndec - length) as $l | . + ("0" * $l);
def abs: if . < 0 then -. else . end;
pow(10; $ndec) as $p
| round as $round
| if $p * ((. - $round)|abs) < 0.1
then ($round|tostring) + "." + ($ndec * "0")
else . * $p | round / $p
| tostring
| capture("(?<left>[^.]*)[.](?<right>.*)")
| .left + "." + (.right|rpad)
end;
# pretty-print a 2-d matrix
def pp($ndec; $width):
def pad(n): tostring | (n - length) * " " + .;
def row: map(round($ndec) | pad($width)) | join(" ");
reduce .[] as $row (""; . + "\n\($row|row)");
</syntaxhighlight>
'''QR-Decomposition'''
<syntaxhighlight lang=jq>
def minor($x; $d):
($x|length) as $nr
| ($x[0]|length) as $nc
| reduce range(0; $d) as $i (matrix($nr;$nc;0); .[$i][$i] = 1)
| reduce range($d; $nr) as $i (.;
reduce range($d;$nc) as $j (.; .[$i][$j] = $x[$i][$j] ) );
def vmadd($a; $b; $s):
reduce range (0; $a|length) as $i ([];
.[$i] = $a[$i] + $s * $b[$i] );
def vmul($v):
($v|length) as $n
| reduce range(0;$n) as $i (null;
reduce range(0;$n) as $j (.; .[$i][$j] = -2 * $v[$i] * $v[$j] ))
| reduce range(0;$n) as $i (.; .[$i][$i] += 1 );
def vnorm($x):
sum($x[] | .*.) | sqrt;
def vdiv($x; $d):
[range (0;$x|length) | $x[.] / $d];
def mcol($m; $c):
[range (0;$m|length) | $m[.][$c]];
def householder($m):
($m|length) as $nr
| ($m[0]|length) as $nc
| { q: [], # $nr
z: $m,
k: 0 }
| until( .k >= $nc or .k >= $nr-1;
.z = minor(.z; .k)
| .x = mcol(.z; .k)
| .a = vnorm(.x)
| if ($m[.k][.k] > 0) then .a = -.a else . end
| .e = [range (0; $nr) as $i | if ($i == .k) then 1 else 0 end]
| .e = vmadd(.x; .e; .a)
| .e = vdiv(.e; vnorm(.e))
| .q[.k] = vmul(.e)
| .z = multiply(.q[.k]; .z)
| .k += 1 )
| .Q = .q[0]
| .R = multiply(.q[0]; $m)
| .i = 1
| until (.i >= $nc or .i >= $nr-1;
.Q = multiply(.q[.i]; .Q)
| .i += 1 )
| .R = multiply(.Q; $m)
| .Q |= transpose
| [.Q, .R] ;
def x: [
[12, -51, 4],
[ 6, 167, -68],
[-4, 24, -41],
[-1, 1, 0],
[ 2, 0, 3]
];
def task:
def pp: pp(3;8);
# Assume $a and $b are conformal
def ssd($a; $b):
[$a[][]] as $a
| [$b[][]] as $b
| ss( range(0;$a|length) | $a[.] - $b[.] );
householder(x) as [$Q, $R]
| multiply($Q; $R) as $m
| "Q:", ($Q|pp),
"\nR:", ($R|pp),
"\nQ * R:", ($m|pp),
"\nSum of squared discrepancies: \(ssd(x; $m))"
;
task
</syntaxhighlight>
{{output}}
<pre>
Q:
0.846 -0.391 0.343 0.082 0.078
0.423 0.904 -0.029 0.026 0.045
-0.282 0.170 0.933 -0.047 -0.137
-0.071 0.014 -0.001 0.980 -0.184
0.141 -0.017 -0.106 -0.171 -0.969
R:
14.177 20.667 -13.402
-0.000 175.043 -70.080
0.000 0.000 -35.202
-0.000 -0.000 -0.000
0.000 0.000 -0.000
Q * R:
12.000 -51.000 4.000
6.000 167.000 -68.000
-4.000 24.000 -41.000
-1.000 1.000 -0.000
2.000 -0.000 3.000
Sum of squared discrepancies: 1.1675699109208862e-26
</pre>
=={{header|Julia}}==
| |||