Rosetta Code

QR decomposition: Difference between revisions

Page Discussion

QR decomposition (view source)

Revision as of 12:03, 29 January 2024

46,054 bytes added ,  3 months ago
m
→‎{{header|Wren}}: Minor tidy
m (syntax highlighting fixup automation)
m (→‎{{header|Wren}}: Minor tidy)
 
(8 intermediate revisions by 3 users not shown)
Line 203:
-0.0000 175.0000 -70.0000
-0.0000 0.0000 35.0000</pre>
 
=={{header|ATS}}==
 
Perhaps not every template function that was written below is actually used. Much of what is below amounts to a little library for working with matrices. To treat blocks and transposes as matrices themselves, I use a trick employed in some Scheme implementations of matrices: indices are mapped by closures, and the closures can nested.
 
<syntaxhighlight lang="ats">
%{^
#include <math.h>
#include <float.h>
%}
 
#include "share/atspre_staload.hats"
 
macdef NAN = g0f2f ($extval (float, "NAN"))
macdef Zero = g0i2f 0
macdef One = g0i2f 1
macdef Two = g0i2f 2
 
(* g0float_sqrt is available in the ats2-xprelude package, but let us
quickly add it here, with implementations for the g0float types
included in the prelude. *)
extern fn {tk : tkind} g0float_sqrt : g0float tk -<> g0float tk
overload sqrt with g0float_sqrt
implement g0float_sqrt<fltknd> x = $extfcall (float, "sqrtf", x)
implement g0float_sqrt<dblknd> x = $extfcall (double, "sqrt", x)
implement g0float_sqrt<ldblknd> x = $extfcall (ldouble, "sqrtl", x)
 
(* Similarly for g0float_copysign. *)
extern fn {tk : tkind}
g0float_copysign : (g0float tk, g0float tk) -<> g0float tk
overload copysign with g0float_copysign
implement
g0float_copysign<fltknd> (x, y) =
$extfcall (float, "copysignf", x, y)
implement
g0float_copysign<dblknd> (x, y) =
$extfcall (double, "copysign", x, y)
implement
g0float_copysign<ldblknd> (x, y) =
$extfcall (ldouble, "copysignl", x, y)
 
(*------------------------------------------------------------------*)
 
typedef Matrix_Index_Map (m1 : int, n1 : int, m0 : int, n0 : int) =
{i1, j1 : pos | i1 <= m1; j1 <= n1}
(int i1, int j1) -<cloref0>
[i0, j0 : pos | i0 <= m0; j0 <= n0]
@(int i0, int j0)
 
datatype Real_Matrix (tk : tkind,
m1 : int, n1 : int,
m0 : int, n0 : int) =
| Real_Matrix of (matrixref (g0float tk, m0, n0),
int m1, int n1, int m0, int n0,
Matrix_Index_Map (m1, n1, m0, n0))
typedef Real_Matrix (tk : tkind, m1 : int, n1 : int) =
[m0, n0 : pos] Real_Matrix (tk, m1, n1, m0, n0)
typedef Real_Vector (tk : tkind, m1 : int, n1 : int) =
[m1 == 1 || n1 == 1] Real_Matrix (tk, m1, n1)
typedef Real_Row (tk : tkind, n1 : int) = Real_Vector (tk, 1, n1)
typedef Real_Column (tk : tkind, m1 : int) = Real_Vector (tk, m1, 1)
 
extern fn {tk : tkind}
Real_Matrix_make_elt :
{m0, n0 : pos}
(int m0, int n0, g0float tk) -< !wrt >
Real_Matrix (tk, m0, n0, m0, n0)
 
extern fn {tk : tkind}
Real_Matrix_copy :
{m1, n1 : pos}
Real_Matrix (tk, m1, n1) -< !refwrt > Real_Matrix (tk, m1, n1)
 
extern fn {tk : tkind}
Real_Matrix_copy_to :
{m1, n1 : pos}
(Real_Matrix (tk, m1, n1), (* destination *)
Real_Matrix (tk, m1, n1)) -< !refwrt >
void
 
extern fn {tk : tkind}
Real_Matrix_fill_with_elt :
{m1, n1 : pos}
(Real_Matrix (tk, m1, n1), g0float tk) -< !refwrt > void
 
extern fn {}
Real_Matrix_dimension :
{tk : tkind}
{m1, n1 : pos}
Real_Matrix (tk, m1, n1) -<> @(int m1, int n1)
 
extern fn {tk : tkind}
Real_Matrix_get_at :
{m1, n1 : pos}
{i1, j1 : pos | i1 <= m1; j1 <= n1}
(Real_Matrix (tk, m1, n1), int i1, int j1) -< !ref > g0float tk
 
extern fn {tk : tkind}
Real_Matrix_set_at :
{m1, n1 : pos}
{i1, j1 : pos | i1 <= m1; j1 <= n1}
(Real_Matrix (tk, m1, n1), int i1, int j1, g0float tk) -< !refwrt >
void
 
extern fn {}
Real_Matrix_transpose :
(* This is transposed INDEXING. It does NOT copy the data. *)
{tk : tkind}
{m1, n1 : pos}
{m0, n0 : pos}
Real_Matrix (tk, m1, n1, m0, n0) -<>
Real_Matrix (tk, n1, m1, m0, n0)
 
extern fn {}
Real_Matrix_block :
(* This is block (submatrix) INDEXING. It does NOT copy the data. *)
{tk : tkind}
{p0, p1 : pos | p0 <= p1}
{q0, q1 : pos | q0 <= q1}
{m1, n1 : pos | p1 <= m1; q1 <= n1}
{m0, n0 : pos}
(Real_Matrix (tk, m1, n1, m0, n0),
int p0, int p1, int q0, int q1) -<>
Real_Matrix (tk, p1 - p0 + 1, q1 - q0 + 1, m0, n0)
 
extern fn {tk : tkind}
Real_Matrix_unit_matrix :
{m : pos}
int m -< !refwrt > Real_Matrix (tk, m, m)
 
extern fn {tk : tkind}
Real_Matrix_unit_matrix_to :
{m : pos}
Real_Matrix (tk, m, m) -< !refwrt > void
 
extern fn {tk : tkind}
Real_Matrix_matrix_sum :
{m, n : pos}
(Real_Matrix (tk, m, n), Real_Matrix (tk, m, n)) -< !refwrt >
Real_Matrix (tk, m, n)
 
extern fn {tk : tkind}
Real_Matrix_matrix_sum_to :
{m, n : pos}
(Real_Matrix (tk, m, n), (* destination*)
Real_Matrix (tk, m, n),
Real_Matrix (tk, m, n)) -< !refwrt >
void
 
extern fn {tk : tkind}
Real_Matrix_matrix_difference :
{m, n : pos}
(Real_Matrix (tk, m, n), Real_Matrix (tk, m, n)) -< !refwrt >
Real_Matrix (tk, m, n)
 
extern fn {tk : tkind}
Real_Matrix_matrix_difference_to :
{m, n : pos}
(Real_Matrix (tk, m, n), (* destination*)
Real_Matrix (tk, m, n),
Real_Matrix (tk, m, n)) -< !refwrt >
void
 
extern fn {tk : tkind}
Real_Matrix_matrix_product :
{m, n, p : pos}
(Real_Matrix (tk, m, n), Real_Matrix (tk, n, p)) -< !refwrt >
Real_Matrix (tk, m, p)
 
extern fn {tk : tkind}
Real_Matrix_matrix_product_to :
{m, n, p : pos}
(Real_Matrix (tk, m, p), (* destination*)
Real_Matrix (tk, m, n),
Real_Matrix (tk, n, p)) -< !refwrt >
void
 
extern fn {tk : tkind}
Real_Matrix_scalar_product :
{m, n : pos}
(Real_Matrix (tk, m, n), g0float tk) -< !refwrt >
Real_Matrix (tk, m, n)
 
extern fn {tk : tkind}
Real_Matrix_scalar_product_2 :
{m, n : pos}
(g0float tk, Real_Matrix (tk, m, n)) -< !refwrt >
Real_Matrix (tk, m, n)
 
extern fn {tk : tkind}
Real_Matrix_scalar_product_to :
{m, n : pos}
(Real_Matrix (tk, m, n), (* destination*)
Real_Matrix (tk, m, n), g0float tk) -< !refwrt > void
 
extern fn {tk : tkind}
Real_Vector_l2norm_squared :
{m, n : pos}
Real_Vector (tk, m, n) -< !ref > g0float tk
 
extern fn {tk : tkind}
Real_Matrix_QR_decomposition :
{m, n : pos}
Real_Matrix (tk, m, n) -< !refwrt >
@(Real_Matrix (tk, m, m), Real_Matrix (tk, m, n))
 
extern fn {tk : tkind}
Real_Matrix_least_squares_solution :
(* This can solve p problems at once. Use p=1 to solve just Ax=b. *)
{m, n, p : pos | n <= m}
(Real_Matrix (tk, m, n), Real_Matrix (tk, m, p)) -< !refwrt >
Real_Matrix (tk, n, p)
 
extern fn {tk : tkind}
Real_Matrix_fprint :
{m, n : pos}
(FILEref, Real_Matrix (tk, m, n)) -<1> void
 
overload copy with Real_Matrix_copy
overload copy_to with Real_Matrix_copy_to
overload fill_with_elt with Real_Matrix_fill_with_elt
overload dimension with Real_Matrix_dimension
overload [] with Real_Matrix_get_at
overload [] with Real_Matrix_set_at
overload transpose with Real_Matrix_transpose
overload block with Real_Matrix_block
overload unit_matrix with Real_Matrix_unit_matrix
overload unit_matrix_to with Real_Matrix_unit_matrix_to
overload matrix_sum with Real_Matrix_matrix_sum
overload matrix_sum_to with Real_Matrix_matrix_sum_to
overload matrix_difference with Real_Matrix_matrix_difference
overload matrix_difference_to with Real_Matrix_matrix_difference_to
overload matrix_product with Real_Matrix_matrix_product
overload matrix_product_to with Real_Matrix_matrix_product_to
overload scalar_product with Real_Matrix_scalar_product
overload scalar_product with Real_Matrix_scalar_product_2
overload scalar_product_to with Real_Matrix_scalar_product_to
overload + with matrix_sum
overload - with matrix_difference
overload * with matrix_product
overload * with scalar_product
 
(* Overload for a Real_Matrix_l2norm_squared, if we decided to have
one, would be given precedence 0. *)
overload l2norm_squared with Real_Vector_l2norm_squared of 1
 
overload QR_decomposition with Real_Matrix_QR_decomposition
overload least_squares_solution with
Real_Matrix_least_squares_solution
 
(*------------------------------------------------------------------*)
 
implement {tk}
Real_Matrix_make_elt (m0, n0, elt) =
Real_Matrix (matrixref_make_elt<g0float tk> (i2sz m0, i2sz n0, elt),
m0, n0, m0, n0, lam (i1, j1) => @(i1, j1))
 
implement {}
Real_Matrix_dimension A =
case+ A of Real_Matrix (_, m1, n1, _, _, _) => @(m1, n1)
 
implement {tk}
Real_Matrix_get_at (A, i1, j1) =
let
val+ Real_Matrix (storage, _, _, _, n0, index_map) = A
val @(i0, j0) = index_map (i1, j1)
in
matrixref_get_at<g0float tk> (storage, pred i0, n0, pred j0)
end
 
implement {tk}
Real_Matrix_set_at (A, i1, j1, x) =
let
val+ Real_Matrix (storage, _, _, _, n0, index_map) = A
val @(i0, j0) = index_map (i1, j1)
in
matrixref_set_at<g0float tk> (storage, pred i0, n0, pred j0, x)
end
 
implement {}
Real_Matrix_transpose A =
let
val+ Real_Matrix (storage, m1, n1, m0, n0, index_map) = A
in
Real_Matrix (storage, n1, m1, m0, n0,
lam (i1, j1) => index_map (j1, i1))
end
 
implement {}
Real_Matrix_block (A, p0, p1, q0, q1) =
let
val+ Real_Matrix (storage, m1, n1, m0, n0, index_map) = A
in
Real_Matrix (storage, succ (p1 - p0), succ (q1 - q0), m0, n0,
lam (i1, j1) =>
index_map (p0 + pred i1, q0 + pred j1))
end
 
implement {tk}
Real_Matrix_copy A =
let
val @(m1, n1) = dimension A
val C = Real_Matrix_make_elt<tk> (m1, n1, A[1, 1])
val () = copy_to<tk> (C, A)
in
C
end
 
implement {tk}
Real_Matrix_copy_to (Dst, Src) =
let
val @(m1, n1) = dimension Src
prval [m1 : int] EQINT () = eqint_make_gint m1
prval [n1 : int] EQINT () = eqint_make_gint n1
 
var i : intGte 1
in
for* {i : pos | i <= m1 + 1} .<(m1 + 1) - i>.
(i : int i) =>
(i := 1; i <> succ m1; i := succ i)
let
var j : intGte 1
in
for* {j : pos | j <= n1 + 1} .<(n1 + 1) - j>.
(j : int j) =>
(j := 1; j <> succ n1; j := succ j)
Dst[i, j] := Src[i, j]
end
end
 
implement {tk}
Real_Matrix_fill_with_elt (A, elt) =
let
val @(m1, n1) = dimension A
prval [m1 : int] EQINT () = eqint_make_gint m1
prval [n1 : int] EQINT () = eqint_make_gint n1
 
var i : intGte 1
in
for* {i : pos | i <= m1 + 1} .<(m1 + 1) - i>.
(i : int i) =>
(i := 1; i <> succ m1; i := succ i)
let
var j : intGte 1
in
for* {j : pos | j <= n1 + 1} .<(n1 + 1) - j>.
(j : int j) =>
(j := 1; j <> succ n1; j := succ j)
A[i, j] := elt
end
end
 
implement {tk}
Real_Matrix_unit_matrix {m} m =
let
val A = Real_Matrix_make_elt<tk> (m, m, Zero)
var i : intGte 1
in
for* {i : pos | i <= m + 1} .<(m + 1) - i>.
(i : int i) =>
(i := 1; i <> succ m; i := succ i)
A[i, i] := One;
A
end
 
implement {tk}
Real_Matrix_unit_matrix_to A =
let
val @(m, _) = dimension A
prval [m : int] EQINT () = eqint_make_gint m
 
var i : intGte 1
in
for* {i : pos | i <= m + 1} .<(m + 1) - i>.
(i : int i) =>
(i := 1; i <> succ m; i := succ i)
let
var j : intGte 1
in
for* {j : pos | j <= m + 1} .<(m + 1) - j>.
(j : int j) =>
(j := 1; j <> succ m; j := succ j)
A[i, j] := (if i = j then One else Zero)
end
end
 
implement {tk}
Real_Matrix_matrix_sum (A, B) =
let
val @(m, n) = dimension A
val C = Real_Matrix_make_elt<tk> (m, n, NAN)
val () = matrix_sum_to<tk> (C, A, B)
in
C
end
 
implement {tk}
Real_Matrix_matrix_sum_to (C, A, B) =
let
val @(m, n) = dimension A
prval [m : int] EQINT () = eqint_make_gint m
prval [n : int] EQINT () = eqint_make_gint n
 
var i : intGte 1
in
for* {i : pos | i <= m + 1} .<(m + 1) - i>.
(i : int i) =>
(i := 1; i <> succ m; i := succ i)
let
var j : intGte 1
in
for* {j : pos | j <= n + 1} .<(n + 1) - j>.
(j : int j) =>
(j := 1; j <> succ n; j := succ j)
C[i, j] := A[i, j] + B[i, j]
end
end
 
implement {tk}
Real_Matrix_matrix_difference (A, B) =
let
val @(m, n) = dimension A
val C = Real_Matrix_make_elt<tk> (m, n, NAN)
val () = matrix_difference_to<tk> (C, A, B)
in
C
end
 
implement {tk}
Real_Matrix_matrix_difference_to (C, A, B) =
let
val @(m, n) = dimension A
prval [m : int] EQINT () = eqint_make_gint m
prval [n : int] EQINT () = eqint_make_gint n
 
var i : intGte 1
in
for* {i : pos | i <= m + 1} .<(m + 1) - i>.
(i : int i) =>
(i := 1; i <> succ m; i := succ i)
let
var j : intGte 1
in
for* {j : pos | j <= n + 1} .<(n + 1) - j>.
(j : int j) =>
(j := 1; j <> succ n; j := succ j)
C[i, j] := A[i, j] - B[i, j]
end
end
 
implement {tk}
Real_Matrix_matrix_product (A, B) =
let
val @(m, n) = dimension A and @(_, p) = dimension B
val C = Real_Matrix_make_elt<tk> (m, p, NAN)
val () = matrix_product_to<tk> (C, A, B)
in
C
end
 
implement {tk}
Real_Matrix_matrix_product_to (C, A, B) =
let
val @(m, n) = dimension A and @(_, p) = dimension B
prval [m : int] EQINT () = eqint_make_gint m
prval [n : int] EQINT () = eqint_make_gint n
prval [p : int] EQINT () = eqint_make_gint p
 
var i : intGte 1
in
for* {i : pos | i <= m + 1} .<(m + 1) - i>.
(i : int i) =>
(i := 1; i <> succ m; i := succ i)
let
var k : intGte 1
in
for* {k : pos | k <= p + 1} .<(p + 1) - k>.
(k : int k) =>
(k := 1; k <> succ p; k := succ k)
let
var j : intGte 1
in
C[i, k] := A[i, 1] * B[1, k];
for* {j : pos | j <= n + 1} .<(n + 1) - j>.
(j : int j) =>
(j := 2; j <> succ n; j := succ j)
C[i, k] :=
C[i, k] + (A[i, j] * B[j, k])
end
end
end
 
implement {tk}
Real_Matrix_scalar_product (A, r) =
let
val @(m, n) = dimension A
val C = Real_Matrix_make_elt<tk> (m, n, NAN)
val () = scalar_product_to<tk> (C, A, r)
in
C
end
 
implement {tk}
Real_Matrix_scalar_product_2 (r, A) =
Real_Matrix_scalar_product<tk> (A, r)
 
implement {tk}
Real_Matrix_scalar_product_to (C, A, r) =
let
val @(m, n) = dimension A
prval [m : int] EQINT () = eqint_make_gint m
prval [n : int] EQINT () = eqint_make_gint n
 
var i : intGte 1
in
for* {i : pos | i <= m + 1} .<(m + 1) - i>.
(i : int i) =>
(i := 1; i <> succ m; i := succ i)
let
var j : intGte 1
in
for* {j : pos | j <= n + 1} .<(n + 1) - j>.
(j : int j) =>
(j := 1; j <> succ n; j := succ j)
C[i, j] := A[i, j] * r
end
end
 
implement {tk}
Real_Vector_l2norm_squared v =
$effmask_wrt
let
val @(m, n) = dimension v
prval [m : int] EQINT () = eqint_make_gint m
prval [n : int] EQINT () = eqint_make_gint n
in
if n = 1 then
let
var sum : g0float tk
var i : intGte 1
val v11 = v[1, 1]
in
sum := v11 * v11;
for* {i : pos | i <= m + 1} .<(m + 1) - i>.
(i : int i) =>
(i := 2; i <> succ m; i := succ i)
let
val vi1 = v[i, 1]
in
sum := sum + (vi1 * vi1)
end;
sum
end
else
let
var sum : g0float tk
var j : intGte 1
val v11 = v[1, 1]
in
sum := v11 * v11;
for* {j : pos | j <= n + 1} .<(n + 1) - j>.
(j : int j) =>
(j := 2; j <> succ n; j := succ j)
let
val v1j = v[1, j]
in
sum := sum + (v1j * v1j)
end;
sum
end
end
 
implement {tk}
Real_Matrix_QR_decomposition A =
(* Some of what follows does needless allocation and work, but
making this code more efficient would be a project of its own!
Also, one would likely want to implement pivot selection. See,
for instance, Businger, P., Golub, G.H. Linear least squares
solutions by householder transformations. Numer. Math. 7, 269–276
(1965). https://doi.org/10.1007/BF01436084
(https://web.archive.org/web/20230514003458/https://pages.stat.wisc.edu/~bwu62/771/businger1965.pdf)
 
Note that I follow
https://en.wikipedia.org/w/index.php?title=QR_decomposition&oldid=1152640697#Using_Householder_reflections
more closely than I do what is stated in the task description at
the time of this writing (13 May 2023). The presentation there
seems simpler to me, and I prefer seeing a norm used to normalize
the u vector. *)
let
val @(m, n) = dimension A
prval [m : int] EQINT () = eqint_make_gint m
prval [n : int] EQINT () = eqint_make_gint n
 
stadef min_mn = min (m, n)
val min_mn : int min_mn = min (m, n)
 
var Q : Real_Matrix (tk, m, m) = unit_matrix<tk> m
val R : Real_Matrix (tk, m, n) = copy A
 
(* I_mm is a unit matrix of the maximum size used. Smaller unit
matrices will be had by the "identity" function, and unit
column vectors by the "unit_column" function. *)
val I_mm : Real_Matrix (tk, m, m) = unit_matrix<tk> m
fn
identity {p : pos | p <= m}
(p : int p) :<> Real_Matrix (tk, p, p) =
block (I_mm, 1, p, 1, p)
fn
unit_column {p, j : pos | j <= p; p <= m}
(p : int p,
j : int j) :<> Real_Column (tk, p) =
block (I_mm, 1, p, j, j)
 
var k : intGte 1
in
for* {k : pos | k <= min_mn} .<min_mn - k>.
(k : int k) =>
(k := 1; k <> min_mn; k := succ k)
let
val x = block (R, k, m, k, k)
val sigma = l2norm_squared x
 
(* Choose the sign of alpha to increase the magnitude of the
pivot. *)
val alpha = copysign (sqrt sigma, ~x[1, 1])
 
val e1 = unit_column (succ (m - k), 1)
val u = x - (alpha * e1)
val v = u * (One / sqrt (l2norm_squared u))
val I = identity (succ (m - k))
val H = I - (Two * v * transpose v)
 
(* Update R, using block operations. *)
val () = fill_with_elt<tk> (x, Zero)
val () = x[1, 1] := alpha
val R_ = block (R, k, m, succ k, n)
val Tmp = H * R_
val () = copy_to (R_, Tmp)
 
(* Update Q. *)
val Tmp = unit_matrix m
val Tmp_ = block (Tmp, k, m, k, m)
val () = copy_to (Tmp_, H)
val () = Q := Q * Tmp
in
end;
@(Q, R)
end
 
implement {tk}
Real_Matrix_least_squares_solution (A, B) =
let
(* I use this algorithm for the back substitutions:
https://algowiki-project.org/algowiki/en/index.php?title=Backward_substitution&oldid=10412#Approaches_and_features_of_implementing_the_back_substitution_algorithm_in_parallel
*)
 
val @(m, n) = dimension A and @(_, p) = dimension B
prval [m : int] EQINT () = eqint_make_gint m
prval [n : int] EQINT () = eqint_make_gint n
prval [p : int] EQINT () = eqint_make_gint p
 
val @(Q, R) = QR_decomposition<tk> A
 
(* X is initialized for back substitutions. *)
val X = block (transpose Q * B, 1, n, 1, p)
and R = block (R, 1, n, 1, n)
 
var k : intGte 1
in
(* Complete the back substitutions. *)
for* {k : pos | k <= p + 1} .<(p + 1) - k>.
(k : int k) =>
(k := 1; k <> succ p; k := succ k)
let
val x = block (X, 1, n, k, k)
var j : intGte 0
in
for* {j : nat | 0 <= j; j <= n} .<j>.
(j : int j) =>
(j := n; j <> 0; j := pred j)
let
var i : intGte 1
in
x[j, 1] := x[j, 1] / R[j, j];
for* {i : pos | i <= j} .<j - i>.
(i : int i) =>
(i := 1; i <> j; i := succ i)
x[i, 1] := x[i, 1] - (R[i, j] * x[j, 1])
end
end;
X
end
 
implement {tk}
Real_Matrix_fprint {m, n} (outf, A) =
let
val @(m, n) = dimension A
var i : intGte 1
in
for* {i : pos | i <= m + 1} .<(m + 1) - i>.
(i : int i) =>
(i := 1; i <> succ m; i := succ i)
let
var j : intGte 1
in
for* {j : pos | j <= n + 1} .<(n + 1) - j>.
(j : int j) =>
(j := 1; j <> succ n; j := succ j)
let
typedef FILEstar = $extype"FILE *"
extern castfn FILEref2star : FILEref -<> FILEstar
val _ = $extfcall (int, "fprintf", FILEref2star outf,
"%16.6g", A[i, j])
in
end;
fprintln! (outf)
end
end
 
(*------------------------------------------------------------------*)
 
implement
main0 () =
let
stadef fltknd = dblknd
macdef i2flt = g0int2float<intknd,dblknd>
 
val A = Real_Matrix_make_elt<fltknd> (3, 3, NAN)
val () =
begin
A[1, 1] := i2flt 12;
A[2, 1] := i2flt 6;
A[3, 1] := i2flt ~4;
 
A[1, 2] := i2flt ~51;
A[2, 2] := i2flt 167;
A[3, 2] := i2flt 24;
 
A[1, 3] := i2flt 4;
A[2, 3] := i2flt ~68;
A[3, 3] := i2flt ~41
end
 
val @(Q, R) = QR_decomposition<fltknd> A
 
(* Example of least-squares solution. (Copied from the BBC BASIC
or Common Lisp entry, whichever you prefer to think it copied
from.) *)
val x = $list (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
and y = $list (1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321)
val X = Real_Matrix_make_elt<fltknd> (11, 3, NAN)
and Y = Real_Matrix_make_elt<fltknd> (11, 1, NAN)
val () =
let
var i : intGte 1
in
for* {i : pos | i <= 12} .<12 - i>.
(i : int i) =>
(i := 1; i <> 12; i := succ i)
let
val xi = x[pred i] : int
and yi = y[pred i] : int
in
X[i, 1] := g0i2f (xi ** 0);
X[i, 2] := g0i2f (xi ** 1);
X[i, 3] := g0i2f (xi ** 2);
Y[i, 1] := g0i2f yi
end
end
val solution = least_squares_solution (X, Y)
in
println! ("A :");
Real_Matrix_fprint (stdout_ref, A);
println! ();
println! ("Q :");
Real_Matrix_fprint (stdout_ref, Q);
println! ();
println! ("R :");
Real_Matrix_fprint (stdout_ref, R);
println! ();
println! ("Q * R :");
Real_Matrix_fprint (stdout_ref, Q * R);
println! ();
println! ("least squares A in Ax=b :");
Real_Matrix_fprint (stdout_ref, X);
println! ();
println! ("least squares b in Ax=b :");
Real_Matrix_fprint (stdout_ref, Y);
println! ();
println! ("least squares solution :");
Real_Matrix_fprint (stdout_ref, solution)
end
 
(*------------------------------------------------------------------*)
</syntaxhighlight>
 
{{out}}
<pre>$ patscc -std=gnu2x -g -O2 -DATS_MEMALLOC_GCBDW qr_decomposition_task.dats -lgc -lm && ./a.out
A :
12 -51 4
6 167 -68
-4 24 -41
 
Q :
-0.857143 0.394286 0.331429
-0.428571 -0.902857 -0.0342857
0.285714 -0.171429 0.942857
 
R :
-14 -21 14
0 -175 70
0 0 -35
 
Q * R :
12 -51 4
6 167 -68
-4 24 -41
 
least squares A in Ax=b :
1 0 0
1 1 1
1 2 4
1 3 9
1 4 16
1 5 25
1 6 36
1 7 49
1 8 64
1 9 81
1 10 100
 
least squares b in Ax=b :
1
6
17
34
57
86
121
162
209
262
321
 
least squares solution :
1
2
3
</pre>
 
=={{header|Axiom}}==
Line 1,102 ⟶ 1,951:
-1.000 1.000 -0.000
2.000 -0.000 3.000
</pre>
 
===With Polynomial Fitting===
<syntaxhighlight lang="c++">
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <string>
#include <vector>
 
class Matrix {
public:
Matrix(const std::vector<std::vector<double>>& data) : data(data) {
initialise();
}
 
Matrix(const Matrix& matrix) : data(matrix.data) {
initialise();
}
 
Matrix(const uint64_t& row_count, const uint64_t& column_count) {
data.assign(row_count, std::vector<double>(column_count, 0.0));
initialise();
}
 
Matrix add(const Matrix& other) {
if ( other.row_count != row_count || other.column_count != column_count ) {
throw std::invalid_argument("Incompatible matrix dimensions.");
}
 
Matrix result(data);
for ( int32_t i = 0; i < row_count; ++i ) {
for ( int32_t j = 0; j < column_count; ++j ) {
result.data[i][j] = data[i][j] + other.data[i][j];
}
}
return result;
}
 
Matrix multiply(const Matrix& other) {
if ( column_count != other.row_count ) {
throw std::invalid_argument("Incompatible matrix dimensions.");
}
 
Matrix result(row_count, other.column_count);
for ( int32_t i = 0; i < row_count; ++i ) {
for ( int32_t j = 0; j < other.column_count; ++j ) {
for ( int32_t k = 0; k < row_count; k++ ) {
result.data[i][j] += data[i][k] * other.data[k][j];
}
}
}
return result;
}
 
Matrix transpose() {
Matrix result(column_count, row_count);
for ( int32_t i = 0; i < row_count; ++i ) {
for ( int32_t j = 0; j < column_count; ++j ) {
result.data[j][i] = data[i][j];
}
}
return result;
}
 
Matrix minor(const int32_t& index) {
Matrix result(row_count, column_count);
for ( int32_t i = 0; i < index; ++i ) {
result.set_entry(i, i, 1.0);
}
 
for ( int32_t i = index; i < row_count; ++i ) {
for ( int32_t j = index; j < column_count; ++j ) {
result.set_entry(i, j, data[i][j]);
}
}
return result;
}
 
Matrix column(const int32_t& index) {
Matrix result(row_count, 1);
for ( int32_t i = 0; i < row_count; ++i ) {
result.set_entry(i, 0, data[i][index]);
}
return result;
}
 
Matrix scalarMultiply(const double& value) {
if ( column_count != 1 ) {
throw std::invalid_argument("Incompatible matrix dimension.");
}
 
Matrix result(row_count, column_count);
for ( int32_t i = 0; i < row_count; ++i ) {
result.data[i][0] = data[i][0] * value;
}
return result;
}
 
Matrix unit() {
if ( column_count != 1 ) {
throw std::invalid_argument("Incompatible matrix dimensions.");
}
 
const double the_magnitude = magnitude();
Matrix result(row_count, column_count);
for ( int32_t i = 0; i < row_count; ++i ) {
result.data[i][0] = data[i][0] / the_magnitude;
}
return result;
}
 
double magnitude() {
if ( column_count != 1 ) {
throw std::invalid_argument("Incompatible matrix dimensions.");
}
 
double norm = 0.0;
for ( int32_t i = 0; i < row_count; ++i ) {
norm += data[i][0] * data[i][0];
}
return std::sqrt(norm);
}
 
int32_t size() {
if ( column_count != 1 ) {
throw std::invalid_argument("Incompatible matrix dimensions.");
}
return row_count;
}
 
void display(const std::string& title) {
std::cout << title << std::endl;
for ( int32_t i = 0; i < row_count; ++i ) {
for ( int32_t j = 0; j < column_count; ++j ) {
std::cout << std::setw(9) << std::fixed << std::setprecision(4) << data[i][j];
}
std::cout << std::endl;
}
std::cout << std::endl;
}
 
double get_entry(const int32_t& row, const int32_t& col) {
return data[row][col];
}
 
void set_entry(const int32_t& row, const int32_t& col, const double& value) {
data[row][col] = value;
}
 
int32_t get_row_count() {
return row_count;
}
 
int32_t get_column_count() {
return column_count;
}
 
private:
void initialise() {
row_count = data.size();
column_count = data[0].size();
}
 
int32_t row_count;
int32_t column_count;
std::vector<std::vector<double>> data;
};
 
typedef std::pair<Matrix, Matrix> matrix_pair;
 
Matrix householder_factor(Matrix vector) {
if ( vector.get_column_count() != 1 ) {
throw std::invalid_argument("Incompatible matrix dimensions.");
}
 
const int32_t size = vector.size();
Matrix result(size, size);
for ( int32_t i = 0; i < size; ++i ) {
for ( int32_t j = 0; j < size; ++j ) {
result.set_entry(i, j, -2 * vector.get_entry(i, 0) * vector.get_entry(j, 0));
}
}
 
for ( int32_t i = 0; i < size; ++i ) {
result.set_entry(i, i, result.get_entry(i, i) + 1.0);
}
return result;
}
 
matrix_pair householder(Matrix matrix) {
const int32_t row_count = matrix.get_row_count();
const int32_t column_count = matrix.get_column_count();
std::vector<Matrix> versions_of_Q;
Matrix z(matrix);
Matrix z1(row_count, column_count);
 
for ( int32_t k = 0; k < column_count && k < row_count - 1; ++k ) {
Matrix vectorE(row_count, 1);
z1 = z.minor(k);
Matrix vectorX = z1.column(k);
double magnitudeX = vectorX.magnitude();
if ( matrix.get_entry(k, k) > 0 ) {
magnitudeX = -magnitudeX;
}
 
for ( int32_t i = 0; i < vectorE.size(); ++i ) {
vectorE.set_entry(i, 0, ( i == k ) ? 1 : 0);
}
vectorE = vectorE.scalarMultiply(magnitudeX).add(vectorX).unit();
versions_of_Q.emplace_back(householder_factor(vectorE));
z = versions_of_Q[k].multiply(z1);
}
 
Matrix Q = versions_of_Q[0];
for ( int32_t i = 1; i < column_count && i < row_count - 1; ++i ) {
Q = versions_of_Q[i].multiply(Q);
}
 
Matrix R = Q.multiply(matrix);
Q = Q.transpose();
return matrix_pair(R, Q);
}
 
Matrix solve_upper_triangular(Matrix r, Matrix b) {
const int32_t column_count = r.get_column_count();
Matrix result(column_count, 1);
 
for ( int32_t k = column_count - 1; k >= 0; --k ) {
double total = 0.0;
for ( int32_t j = k + 1; j < column_count; ++j ) {
total += r.get_entry(k, j) * result.get_entry(j, 0);
}
result.set_entry(k, 0, ( b.get_entry(k, 0) - total ) / r.get_entry(k, k));
}
return result;
}
 
Matrix least_squares(Matrix vandermonde, Matrix b) {
matrix_pair pair = householder(vandermonde);
return solve_upper_triangular(pair.first, pair.second.transpose().multiply(b));
}
 
Matrix fit_polynomial(Matrix x, Matrix y, const int32_t& polynomial_degree) {
Matrix vandermonde(x.get_column_count(), polynomial_degree + 1);
for ( int32_t i = 0; i < x.get_column_count(); ++i ) {
for ( int32_t j = 0; j < polynomial_degree + 1; ++j ) {
vandermonde.set_entry(i, j, std::pow(x.get_entry(0, i), j));
}
}
return least_squares(vandermonde, y.transpose());
}
 
int main() {
const std::vector<std::vector<double>> data = { { 12.0, -51.0, 4.0 },
{ 6.0, 167.0, -68.0 },
{ -4.0, 24.0, -41.0 },
{ -1.0, 1.0, 0.0 },
{ 2.0, 0.0, 3.0 } };
 
// Task 1
Matrix A(data);
A.display("Initial matrix A:");
 
matrix_pair pair = householder(A);
Matrix Q = pair.second;
Matrix R = pair.first;
 
Q.display("Matrix Q:");
R.display("Matrix R:");
 
Matrix result = Q.multiply(R);
result.display("Matrix Q * R:");
 
// Task 2
Matrix x( std::vector<std::vector<double>>{ { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 } } );
Matrix y( std::vector<std::vector<double>>{
{ 1.0, 6.0, 17.0, 34.0, 57.0, 86.0, 121.0, 162.0, 209.0, 262.0, 321.0 } } );
 
result = fit_polynomial(x, y, 2);
result.display("Result of fitting polynomial:");
}
</syntaxhighlight>
{{ out }}
<pre>
Initial matrix A:
12.0000 -51.0000 4.0000
6.0000 167.0000 -68.0000
-4.0000 24.0000 -41.0000
-1.0000 1.0000 0.0000
2.0000 0.0000 3.0000
 
Matrix Q:
0.8464 -0.3913 0.3431 0.0815 0.0781
0.4232 0.9041 -0.0293 0.0258 0.0447
-0.2821 0.1704 0.9329 -0.0474 -0.1374
-0.0705 0.0140 -0.0011 0.9804 -0.1836
0.1411 -0.0167 -0.1058 -0.1713 -0.9692
 
Matrix R:
14.1774 20.6666 -13.4016
-0.0000 175.0425 -70.0803
0.0000 0.0000 -35.2015
-0.0000 -0.0000 -0.0000
0.0000 0.0000 -0.0000
 
Matrix Q * R:
12.0000 -51.0000 4.0000
6.0000 167.0000 -68.0000
-4.0000 24.0000 -41.0000
-1.0000 1.0000 -0.0000
2.0000 -0.0000 3.0000
 
Result of fitting polynomial:
1.0000
2.0000
3.0000
</pre>
 
Line 2,101 ⟶ 3,269:
0,000 -175,000 70,000
0,000 0,000 35,000</pre>
 
===Without external libraries===
<syntaxhighlight lang="java">
import java.util.ArrayList;
import java.util.List;
 
public final class QRDecomposition {
 
public static void main(String[] aArgs) {
final double[][] data = new double [][] { { 12.0, -51.0, 4.0 },
{ 6.0, 167.0, -68.0 },
{ -4.0, 24.0, -41.0 },
{ -1.0, 1.0, 0.0 },
{ 2.0, 0.0, 3.0 } };
// Task 1
Matrix A = new Matrix(data);
A.display("Initial matrix A:");
MatrixPair pair = householder(A);
Matrix Q = pair.q;
Matrix R = pair.r;
Q.display("Matrix Q:");
R.display("Matrix R:");
Matrix result = Q.multiply(R);
result.display("Matrix Q * R:");
// Task 2
Matrix x = new Matrix ( new double[][] { { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 } } );
Matrix y = new Matrix(
new double[][] { { 1.0, 6.0, 17.0, 34.0, 57.0, 86.0, 121.0, 162.0, 209.0, 262.0, 321.0 } } );
result = fitPolynomial(x, y, 2);
result.display("Result of fitting polynomial:");
}
private static MatrixPair householder(Matrix aMatrix) {
final int rowCount = aMatrix.getRowCount();
final int columnCount = aMatrix.getColumnCount();
List<Matrix> versionsOfQ = new ArrayList<Matrix>(rowCount);
Matrix z = new Matrix(aMatrix);
Matrix z1 = new Matrix(rowCount, columnCount);
for ( int k = 0; k < columnCount && k < rowCount - 1; k++ ) {
Matrix vectorE = new Matrix(rowCount, 1);
z1 = z.minor(k);
Matrix vectorX = z1.column(k);
double magnitudeX = vectorX.magnitude();
if ( aMatrix.getEntry(k, k) > 0 ) {
magnitudeX = -magnitudeX;
}
 
for ( int i = 0; i < vectorE.size(); i++ ) {
vectorE.setEntry(i, 0, ( i == k ) ? 1 : 0);
}
vectorE = vectorE.scalarMultiply(magnitudeX).add(vectorX).unit();
versionsOfQ.add(householderFactor(vectorE));
z = versionsOfQ.get(k).multiply(z1);
}
Matrix Q = versionsOfQ.get(0);
for ( int i = 1; i < columnCount && i < rowCount - 1; i++ ) {
Q = versionsOfQ.get(i).multiply(Q);
}
 
Matrix R = Q.multiply(aMatrix);
Q = Q.transpose();
return new MatrixPair(R, Q);
}
public static Matrix householderFactor(Matrix aVector) {
if ( aVector.getColumnCount() != 1 ) {
throw new RuntimeException("Incompatible matrix dimensions.");
}
final int size = aVector.size();
Matrix result = new Matrix(size, size);
for ( int i = 0; i < size; i++ ) {
for ( int j = 0; j < size; j++ ) {
result.setEntry(i, j, -2 * aVector.getEntry(i, 0) * aVector.getEntry(j, 0));
}
}
for ( int i = 0; i < size; i++ ) {
result.setEntry(i, i, result.getEntry(i, i) + 1.0);
}
return result;
}
private static Matrix fitPolynomial(Matrix aX, Matrix aY, int aPolynomialDegree) {
Matrix vandermonde = new Matrix(aX.getColumnCount(), aPolynomialDegree + 1);
for ( int i = 0; i < aX.getColumnCount(); i++ ) {
for ( int j = 0; j < aPolynomialDegree + 1; j++ ) {
vandermonde.setEntry(i, j, Math.pow(aX.getEntry(0, i), j));
}
}
return leastSquares(vandermonde, aY.transpose());
}
private static Matrix leastSquares(Matrix aVandermonde, Matrix aB) {
MatrixPair pair = householder(aVandermonde);
return solveUpperTriangular(pair.r, pair.q.transpose().multiply(aB));
}
private static Matrix solveUpperTriangular(Matrix aR, Matrix aB) {
final int columnCount = aR.getColumnCount();
Matrix result = new Matrix(columnCount, 1);
 
for ( int k = columnCount - 1; k >= 0; k-- ) {
double total = 0.0;
for ( int j = k + 1; j < columnCount; j++ ) {
total += aR.getEntry(k, j) * result.getEntry(j, 0);
}
result.setEntry(k, 0, ( aB.getEntry(k, 0) - total ) / aR.getEntry(k, k));
}
return result;
}
private static record MatrixPair(Matrix r, Matrix q) {}
 
}
 
final class Matrix {
 
public Matrix(double[][] aData) {
rowCount = aData.length;
columnCount = aData[0].length;
data = new double[rowCount][columnCount];
for ( int i = 0; i < rowCount; i++ ) {
for ( int j = 0; j < columnCount; j++ ) {
data[i][j] = aData[i][j];
}
}
}
public Matrix(Matrix aMatrix) {
this(aMatrix.data);
}
public Matrix(int aRowCount, int aColumnCount) {
this( new double[aRowCount][aColumnCount] );
}
 
public Matrix add(Matrix aOther) {
if ( aOther.rowCount != rowCount || aOther.columnCount != columnCount ) {
throw new IllegalArgumentException("Incompatible matrix dimensions.");
}
Matrix result = new Matrix(data);
for ( int i = 0; i < rowCount; i++ ) {
for ( int j = 0; j < columnCount; j++ ) {
result.data[i][j] = data[i][j] + aOther.data[i][j];
}
}
return result;
}
 
public Matrix multiply(Matrix aOther) {
if ( columnCount != aOther.rowCount ) {
throw new IllegalArgumentException("Incompatible matrix dimensions.");
}
Matrix result = new Matrix(rowCount, aOther.columnCount);
for ( int i = 0; i < rowCount; i++ ) {
for ( int j = 0; j < aOther.columnCount; j++ ) {
for ( int k = 0; k < rowCount; k++ ) {
result.data[i][j] += data[i][k] * aOther.data[k][j];
}
}
}
return result;
}
public Matrix transpose() {
Matrix result = new Matrix(columnCount, rowCount);
for ( int i = 0; i < rowCount; i++ ) {
for ( int j = 0; j < columnCount; j++ ) {
result.data[j][i] = data[i][j];
}
}
return result;
}
public Matrix minor(int aIndex) {
Matrix result = new Matrix(rowCount, columnCount);
for ( int i = 0; i < aIndex; i++ ) {
result.setEntry(i, i, 1.0);
}
for ( int i = aIndex; i < rowCount; i++ ) {
for ( int j = aIndex; j < columnCount; j++ ) {
result.setEntry(i, j, data[i][j]);
}
}
return result;
}
public Matrix column(int aIndex) {
Matrix result = new Matrix(rowCount, 1);
for ( int i = 0; i < rowCount; i++ ) {
result.setEntry(i, 0, data[i][aIndex]);
}
return result;
}
public Matrix scalarMultiply(double aValue) {
if ( columnCount != 1 ) {
throw new IllegalArgumentException("Incompatible matrix dimension.");
}
Matrix result = new Matrix(rowCount, columnCount);
for ( int i = 0; i < rowCount; i++ ) {
result.data[i][0] = data[i][0] * aValue;
}
return result;
}
public Matrix unit() {
if ( columnCount != 1 ) {
throw new IllegalArgumentException("Incompatible matrix dimensions.");
}
final double magnitude = magnitude();
Matrix result = new Matrix(rowCount, columnCount);
for ( int i = 0; i < rowCount; i++ ) {
result.data[i][0] = data[i][0] / magnitude;
}
return result;
}
public double magnitude() {
if ( columnCount != 1 ) {
throw new IllegalArgumentException("Incompatible matrix dimensions.");
}
double norm = 0.0;
for ( int i = 0; i < data.length; i++ ) {
norm += data[i][0] * data[i][0];
}
return Math.sqrt(norm);
}
public int size() {
if ( columnCount != 1 ) {
throw new IllegalArgumentException("Incompatible matrix dimensions.");
}
return rowCount;
}
public void display(String aTitle) {
System.out.println(aTitle);
for ( int i = 0; i < rowCount; i++ ) {
for ( int j = 0; j < columnCount; j++ ) {
System.out.print(String.format("%9.4f", data[i][j]));
}
System.out.println();
}
System.out.println();
}
public double getEntry(int aRow, int aColumn) {
return data[aRow][aColumn];
}
public void setEntry(int aRow, int aColumn, double aValue) {
data[aRow][aColumn] = aValue;
}
public int getRowCount() {
return rowCount;
}
public int getColumnCount() {
return columnCount;
}
private final int rowCount;
private final int columnCount;
private final double[][] data;
}
</syntaxhighlight>
{{ out }}
<pre>
Initial matrix A:
12.0000 -51.0000 4.0000
6.0000 167.0000 -68.0000
-4.0000 24.0000 -41.0000
-1.0000 1.0000 0.0000
2.0000 0.0000 3.0000
 
Matrix Q:
0.8464 -0.3913 0.3431 0.0815 0.0781
0.4232 0.9041 -0.0293 0.0258 0.0447
-0.2821 0.1704 0.9329 -0.0474 -0.1374
-0.0705 0.0140 -0.0011 0.9804 -0.1836
0.1411 -0.0167 -0.1058 -0.1713 -0.9692
 
Matrix R:
14.1774 20.6666 -13.4016
-0.0000 175.0425 -70.0803
0.0000 0.0000 -35.2015
-0.0000 -0.0000 -0.0000
0.0000 0.0000 -0.0000
 
Matrix Q * R:
12.0000 -51.0000 4.0000
6.0000 167.0000 -68.0000
-4.0000 24.0000 -41.0000
-1.0000 1.0000 -0.0000
2.0000 -0.0000 3.0000
 
Result of fitting polynomial:
1.0000
2.0000
3.0000
</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}}==
Line 3,789 ⟶ 5,442:
{{libheader|Wren-matrix}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="ecmascriptwren">import "./matrix" for Matrix
import "./fmt" for Fmt
 
var minor = Fn.new { |x, d|
PureFox
9,479

edits

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