QR decomposition: Difference between revisions
m
→{{header|Wren}}: Minor tidy
m (Changed signValue function) |
m (→{{header|Wren}}: Minor tidy) |
||
| (98 intermediate revisions by 44 users not shown) | |||
Line 1:
{{task|Matrices}}[[Category:Mathematics]]
Any rectangular <math>m \times n</math> matrix <math>\mathit A</math> can be decomposed to a product of
'''Task'''
Line 11:
-4 & 24 & -41 \end{pmatrix}</math>
and the usage for linear least squares problems on the example from [[
'''Method'''
Line 98:
=={{header|Ada}}==
Output matches that of Matlab solution, not tested with other matrices.
<syntaxhighlight lang="ada">
with Ada.Text_IO; use Ada.Text_IO;
with Ada.Numerics.Real_Arrays; use Ada.Numerics.Real_Arrays;
Line 110:
for col in mat'Range (2) loop
FIO.Put (mat (row, col), Exp => 0, Aft => 4, Fore => 5);
end loop
New_Line;
end loop;
end Show;
Line 192 ⟶ 193:
Put_Line ("Q:"); Show (Q);
Put_Line ("R:"); Show (R);
end QR;</
{{out}}
<pre>Q:
Line 202 ⟶ 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}}==
The following provides a generic QR decomposition for arbitrary precision floats, double floats and exact calculations:
<syntaxhighlight lang="axiom">)abbrev package TESTP TestPackage
TestPackage(R:Join(Field,RadicalCategory)): with
unitVector: NonNegativeInteger -> Vector(R)
"^": (Vector(R),NonNegativeInteger) -> Vector(R)
solveUpperTriangular: (Matrix(R),Vector(R)) -> Vector(R)
signValue: R -> R
householder:
qr: Matrix(R) -> Record(q:Matrix(R),r:Matrix(R))
lsqr: (Matrix(R),Vector(R)) -> Vector(R)
Line 223 ⟶ 1,071:
out(1) := 1@R
out
v:Vector(R) ^ n:NonNegativeInteger == map((vi:R):R +-> vi^n, v)$Vector(R)
solveUpperTriangular(r,b) ==
Line 238 ⟶ 1,085:
if sqrt(r*r) = r then 1 else -1
householder(a) ==
m :=
scalarMatrix(m,1) - beta*transpose(outerProduct(v,v))
qr(a) ==
(m,n) := (nrows a, ncols a)
qm := scalarMatrix(m,1)
rm
for i in 1..(if m=n then n-1 else n) repeat
h := scalarMatrix(m,1)
setsubMatrix!(h,i,i,householder x)
Line 261 ⟶ 1,104:
dc := qr a
n := ncols(dc.r)
solveUpperTriangular(subMatrix(dc.r,1,n,1,n),
polyfit(x,y,n) ==
a := new(#x,n+1,0@R)$Matrix(R)
for j in 0..n repeat
setColumn!(a,j+1,x^j)
lsqr(a,y)</
This can be called using:
<
qr m
x := vector [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10];
y := vector [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321];
polyfit(x, y, 2)</
With output in exact form:
<
+ 6 69 58 +
Line 294 ⟶ 1,136:
[1,2,3]
Type: Vector(AlgebraicNumber)</
The calculations are comparable to those from the default QR decomposition in R.
=={{header|BBC BASIC}}==
{{works with|BBC BASIC for Windows}}
Makes heavy use of BBC BASIC's matrix arithmetic.
<syntaxhighlight lang="bbcbasic"> *FLOAT 64
@% = &2040A
INSTALL @lib$+"ARRAYLIB"
REM Test matrix for QR decomposition:
DIM A(2,2)
A() = 12, -51, 4, \
\ 6, 167, -68, \
\ -4, 24, -41
REM Do the QR decomposition:
DIM Q(2,2), R(2,2)
PROCqrdecompose(A(), Q(), R())
PRINT "Q:"
PRINT Q(0,0), Q(0,1), Q(0,2)
PRINT Q(1,0), Q(1,1), Q(1,2)
PRINT Q(2,0), Q(2,1), Q(2,2)
PRINT "R:"
PRINT R(0,0), R(0,1), R(0,2)
PRINT R(1,0), R(1,1), R(1,2)
PRINT R(2,0), R(2,1), R(2,2)
REM Test data for least-squares solution:
DIM x(10) : x() = 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
DIM y(10) : y() = 1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321
REM Do the least-squares solution:
DIM a(10,2), q(10,10), r(10,2), t(10,10), b(10), z(2)
FOR i% = 0 TO 10
FOR j% = 0 TO 2
a(i%,j%) = x(i%) ^ j%
NEXT
NEXT
PROCqrdecompose(a(), q(), r())
PROC_transpose(q(),t())
b() = t() . y()
FOR k% = 2 TO 0 STEP -1
s = 0
IF k% < 2 THEN
FOR j% = k%+1 TO 2
s += r(k%,j%) * z(j%)
NEXT
ENDIF
z(k%) = (b(k%) - s) / r(k%,k%)
NEXT k%
PRINT '"Least-squares solution:"
PRINT z(0), z(1), z(2)
END
DEF PROCqrdecompose(A(), Q(), R())
LOCAL i%, k%, m%, n%, H()
m% = DIM(A(),1) : n% = DIM(A(),2)
DIM H(m%,m%)
FOR i% = 0 TO m% : Q(i%,i%) = 1 : NEXT
WHILE n%
PROCqrstep(n%, k%, A(), H())
A() = H() . A()
Q() = Q() . H()
k% += 1
m% -= 1
n% -= 1
ENDWHILE
R() = A()
ENDPROC
DEF PROCqrstep(n%, k%, A(), H())
LOCAL a(), h(), i%, j%
DIM a(n%,0), h(n%,n%)
FOR i% = 0 TO n% : a(i%,0) = A(i%+k%,k%) : NEXT
PROChouseholder(h(), a())
H() = 0 : H(0,0) = 1
FOR i% = 0 TO n%
FOR j% = 0 TO n%
H(i%+k%,j%+k%) = h(i%,j%)
NEXT
NEXT
ENDPROC
REM Create the Householder matrix for the supplied column vector:
DEF PROChouseholder(H(), a())
LOCAL e(), u(), v(), vt(), vvt(), I(), d()
LOCAL i%, n% : n% = DIM(a(),1)
REM Create the scaled standard basis vector e():
DIM e(n%,0) : e(0,0) = SGN(a(0,0)) * MOD(a())
REM Create the normal vector u():
DIM u(n%,0) : u() = a() + e()
REM Normalise with respect to the first element:
DIM v(n%,0) : v() = u() / u(0,0)
REM Get the transpose of v() and its dot product with v():
DIM vt(0,n%), d(0) : PROC_transpose(v(), vt()) : d() = vt() . v()
REM Get the product of v() and vt():
DIM vvt(n%,n%) : vvt() = v() . vt()
REM Create an identity matrix I():
DIM I(n%,n%) : FOR i% = 0 TO n% : I(i%,i%) = 1 : NEXT
REM Create the Householder matrix H() = I - 2/vt()v() v()vt():
vvt() *= 2 / d(0) : H() = I() - vvt()
ENDPROC</syntaxhighlight>
'''Output:'''
<pre>
Q:
-0.8571 0.3943 0.3314
-0.4286 -0.9029 -0.0343
0.2857 -0.1714 0.9429
R:
-14.0000 -21.0000 14.0000
0.0000 -175.0000 70.0000
0.0000 0.0000 -35.0000
Least-squares solution:
1.0000 2.0000 3.0000
</pre>
=={{header|C}}==
<
#include <stdlib.h>
#include <math.h>
Line 309 ⟶ 1,267:
{
mat x = malloc(sizeof(mat_t));
x->v = malloc(sizeof(double*) * m);
x->v[0] = calloc(sizeof(double), m * n);
for (int i = 0; i < m; i++)
Line 336 ⟶ 1,294:
}
mat matrix_copy(
{
mat x = matrix_new(m, n);
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
x->v[i][j] =
return x;
}
Line 415 ⟶ 1,373:
for(int i = 0; i < m->m; i++) {
for (int j = 0; j < m->n; j++) {
printf("
}
printf("\n");
Line 426 ⟶ 1,384:
mat q[m->m];
mat z = m, z1;
for (int k = 0; k < m->n && k < m->m - 1; k++) {
double e[m->m], x[m->m], a;
z1 = matrix_minor(z, k);
Line 449 ⟶ 1,407:
*Q = q[0];
*R = matrix_mul(q[0], m);
for (int i = 1; i < m->n && i < m->m - 1; i++) {
z1 = matrix_mul(q[i], *Q);
if (i > 1) matrix_delete(*Q);
Line 463 ⟶ 1,421:
double in[][3] = {
{ 12, -51, 4
{ 6, 167, -68
{ -4, 24, -41
{ -1, 1, 0},
{ 2, 0, 3},
};
Line 471 ⟶ 1,431:
{
mat R, Q;
mat x = matrix_copy(
householder(x, &R, &Q);
// to show their product is the input matrix
mat m = matrix_mul(Q, R);
puts("Q * R"); matrix_show(m);
matrix_delete(x);
matrix_delete(R);
matrix_delete(Q);
matrix_delete(m);
return 0;
}</syntaxhighlight>
{{out}}
<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.
-0.000 175.
-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
</pre>
=={{header|C sharp|C#}}==
{{libheader|Math.Net}}
<syntaxhighlight lang="csharp">using System;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;
class Program
{
static void Main(string[] args)
{
Matrix<double> A = DenseMatrix.OfArray(new double[,]
{
{ 12, -51, 4 },
{ 6, 167, -68 },
{ -4, 24, -41 }
});
Console.WriteLine("A:");
Console.WriteLine(A);
var qr = A.QR();
Console.WriteLine();
Console.WriteLine("Q:");
Console.WriteLine(qr.Q);
Console.WriteLine();
Console.WriteLine("R:");
Console.WriteLine(qr.R);
}
}</syntaxhighlight>
{{out}}
<pre>A:
DenseMatrix 3x3-Double
12 -51 4
6 167 -68
-4 24 -41
Q:
DenseMatrix 3x3-Double
-0.857143 0.394286 -0.331429
-0.428571 -0.902857 0.0342857
0.285714 -0.171429 -0.942857
R:
DenseMatrix 3x3-Double
-14 -21 14
0 -175 70
0 0 35</pre>
=={{header|C++}}==
<syntaxhighlight lang="cpp">/*
* g++ -O3 -Wall --std=c++11 qr_standalone.cpp -o qr_standalone
*/
#include <cstdio>
#include <cstdlib>
#include <cstring> // for memset
#include <limits>
#include <iostream>
#include <vector>
#include <math.h>
class Vector;
class Matrix {
public:
// default constructor (don't allocate)
Matrix() : m(0), n(0), data(nullptr) {}
// constructor with memory allocation, initialized to zero
Matrix(int m_, int n_) : Matrix() {
m = m_;
n = n_;
allocate(m_,n_);
}
// copy constructor
Matrix(const Matrix& mat) : Matrix(mat.m,mat.n) {
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
(*this)(i,j) = mat(i,j);
}
// constructor from array
template<int rows, int cols>
Matrix(double (&a)[rows][cols]) : Matrix(rows,cols) {
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
(*this)(i,j) = a[i][j];
}
// destructor
~Matrix() {
deallocate();
}
// access data operators
double& operator() (int i, int j) {
return data[i+m*j]; }
double operator() (int i, int j) const {
return data[i+m*j]; }
// operator assignment
Matrix& operator=(const Matrix& source) {
// self-assignment check
if (this != &source) {
if ( (m*n) != (source.m * source.n) ) { // storage cannot be reused
allocate(source.m,source.n); // re-allocate storage
}
// storage can be used, copy data
std::copy(source.data, source.data + source.m*source.n, data);
}
return *this;
}
// compute minor
void compute_minor(const Matrix& mat, int d) {
allocate(mat.m, mat.n);
for (int i = 0; i < d; i++)
(*this)(i,i) = 1.0;
for (int i = d; i < mat.m; i++)
for (int j = d; j < mat.n; j++)
(*this)(i,j) = mat(i,j);
}
// Matrix multiplication
// c = a * b
// c will be re-allocated here
void mult(const Matrix& a, const Matrix& b) {
if (a.n != b.m) {
std::cerr << "Matrix multiplication not possible, sizes don't match !\n";
return;
}
// reallocate ourself if necessary i.e. current Matrix has not valid sizes
if (a.m != m or b.n != n)
allocate(a.m, b.n);
memset(data,0,m*n*sizeof(double));
for (int i = 0; i < a.m; i++)
for (int j = 0; j < b.n; j++)
for (int k = 0; k < a.n; k++)
(*this)(i,j) += a(i,k) * b(k,j);
}
void transpose() {
for (int i = 0; i < m; i++) {
for (int j = 0; j < i; j++) {
double t = (*this)(i,j);
(*this)(i,j) = (*this)(j,i);
(*this)(j,i) = t;
}
}
}
// take c-th column of m, put in v
void extract_column(Vector& v, int c);
// memory allocation
void allocate(int m_, int n_) {
// if already allocated, memory is freed
deallocate();
// new sizes
m = m_;
n = n_;
data = new double[m_*n_];
memset(data,0,m_*n_*sizeof(double));
} // allocate
// memory free
void deallocate() {
if (data)
delete[] data;
data = nullptr;
}
int m, n;
private:
double* data;
}; // struct Matrix
// column vector
class Vector {
public:
// default constructor (don't allocate)
Vector() : size(0), data(nullptr) {}
// constructor with memory allocation, initialized to zero
Vector(int size_) : Vector() {
size = size_;
allocate(size_);
}
// destructor
~Vector() {
deallocate();
}
// access data operators
double& operator() (int i) {
return data[i]; }
double operator() (int i) const {
return data[i]; }
// operator assignment
Vector& operator=(const Vector& source) {
// self-assignment check
if (this != &source) {
if ( size != (source.size) ) { // storage cannot be reused
allocate(source.size); // re-allocate storage
}
// storage can be used, copy data
std::copy(source.data, source.data + source.size, data);
}
return *this;
}
// memory allocation
void allocate(int size_) {
deallocate();
// new sizes
size = size_;
data = new double[size_];
memset(data,0,size_*sizeof(double));
} // allocate
// memory free
void deallocate() {
if (data)
delete[] data;
data = nullptr;
}
// ||x||
double norm() {
double sum = 0;
for (int i = 0; i < size; i++) sum += (*this)(i) * (*this)(i);
return sqrt(sum);
}
// divide data by factor
void rescale(double factor) {
for (int i = 0; i < size; i++) (*this)(i) /= factor;
}
void rescale_unit() {
double factor = norm();
rescale(factor);
}
int size;
private:
double* data;
}; // class Vector
// c = a + b * s
void vmadd(const Vector& a, const Vector& b, double s, Vector& c)
{
if (c.size != a.size or c.size != b.size) {
std::cerr << "[vmadd]: vector sizes don't match\n";
return;
}
for (int i = 0; i < c.size; i++)
c(i) = a(i) + s * b(i);
}
// mat = I - 2*v*v^T
// !!! m is allocated here !!!
void compute_householder_factor(Matrix& mat, const Vector& v)
{
int n = v.size;
mat.allocate(n,n);
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
mat(i,j) = -2 * v(i) * v(j);
for (int i = 0; i < n; i++)
mat(i,i) += 1;
}
// take c-th column of a matrix, put results in Vector v
void Matrix::extract_column(Vector& v, int c) {
if (m != v.size) {
std::cerr << "[Matrix::extract_column]: Matrix and Vector sizes don't match\n";
return;
}
for (int i = 0; i < m; i++)
v(i) = (*this)(i,c);
}
void matrix_show(const Matrix& m, const std::string& str="")
{
std::cout << str << "\n";
for(int i = 0; i < m.m; i++) {
for (int j = 0; j < m.n; j++) {
printf(" %8.3f", m(i,j));
}
printf("\n");
}
printf("\n");
}
// L2-norm ||A-B||^2
double matrix_compare(const Matrix& A, const Matrix& B) {
// matrices must have same size
if (A.m != B.m or A.n != B.n)
return std::numeric_limits<double>::max();
double res=0;
for(int i = 0; i < A.m; i++) {
for (int j = 0; j < A.n; j++) {
res += (A(i,j)-B(i,j)) * (A(i,j)-B(i,j));
}
}
res /= A.m*A.n;
return res;
}
void householder(Matrix& mat,
Matrix& R,
Matrix& Q)
{
int m = mat.m;
int n = mat.n;
// array of factor Q1, Q2, ... Qm
std::vector<Matrix> qv(m);
// temp array
Matrix z(mat);
Matrix z1;
for (int k = 0; k < n && k < m - 1; k++) {
Vector e(m), x(m);
double a;
// compute minor
z1.compute_minor(z, k);
// extract k-th column into x
z1.extract_column(x, k);
a = x.norm();
if (mat(k,k) > 0) a = -a;
for (int i = 0; i < e.size; i++)
e(i) = (i == k) ? 1 : 0;
// e = x + a*e
vmadd(x, e, a, e);
// e = e / ||e||
e.rescale_unit();
// qv[k] = I - 2 *e*e^T
compute_householder_factor(qv[k], e);
// z = qv[k] * z1
z.mult(qv[k], z1);
}
Q = qv[0];
// after this loop, we will obtain Q (up to a transpose operation)
for (int i = 1; i < n && i < m - 1; i++) {
z1.mult(qv[i], Q);
Q = z1;
}
R.mult(Q, mat);
Q.transpose();
}
double in[][3] = {
{ 12, -51, 4},
{ 6, 167, -68},
{ -4, 24, -41},
{ -1, 1, 0},
{ 2, 0, 3},
};
int main()
{
Matrix A(in);
Matrix Q, R;
matrix_show(A,"A");
// compute QR decompostion
householder(A, R, Q);
matrix_show(Q,"Q");
matrix_show(R,"R");
// compare Q*R to the original matrix A
Matrix A_check;
A_check.mult(Q, R);
// compute L2 norm ||A-A_check||^2
double l2 = matrix_compare(A,A_check);
// display Q*R
matrix_show(A_check, l2 < 1e-12 ? "A == Q * R ? yes" : "A == Q * R ? no");
return EXIT_SUCCESS;
}
</syntaxhighlight>
{{out}}
<pre>
A
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
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
A == Q * R ? yes
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
</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>
=={{header|Common Lisp}}==
Line 495 ⟶ 2,276:
Helper functions:
<
(if (zerop x)
x
Line 549 ⟶ 2,330:
C))
</syntaxhighlight>
Main routines:
<
(defun make-householder (a)
(let* ((m (car (array-dimensions a)))
Line 587 ⟶ 2,368:
;; Return Q and R.
(values Q A)))
</syntaxhighlight>
Example 1:
<
#2A((-0.85 0.39 0.33)
Line 599 ⟶ 2,380:
#2A((-14.0 -21.0 14.0)
( 0.0 -175.0 70.0)
( 0.0 0.0 -35.0))</
Example 2, [[Polynomial regression]]:
<
(let* ((m (cadr (array-dimensions x)))
(A (make-array `(,m ,(+ n 1)) :initial-element 0)))
Line 631 ⟶ 2,412:
(aref x j 0))))
(aref R k k))))
x))</
<
(let ((x #2A((0 1 2 3 4 5 6 7 8 9 10)))
(y #2A((1 6 17 34 57 86 121 162 209 262 321))))
(polyfit x y 2))
#2A((0.999999966345088) (2.000000015144699) (2.99999999879804))</
=={{header|D}}==
{{trans|Common Lisp}}
Uses the functions copied from [[Element-wise_operations]], [[Matrix multiplication]], and [[Matrix transposition]].
<
std.typecons, std.numeric, std.range, std.conv;
R[r][] = mixin("row[] " ~ op ~ "B");
}
T[][] elementwiseMat(T, U)(in T[][] A, in U[][] B)
pure nothrow if (is(Unqual!T == Unqual!U)) {
assert(A.length == B.length);
if (A.empty)
return null;
auto R = new typeof(return)(A.length, A[0].length);
foreach (immutable r, const row; A) {
assert(row.length == B[r].length);
R[r][] = mixin("row[] " ~ op ~ "B[r][]");
}
return R;
}
alias mSum = elementwiseMat!q{ + },
pMul = elementwiseMat!q{ * },
pDiv = elementwiseMat!q{ / };
bool isRectangular(T)(in T[][]
return mat.all!(r => r.length == mat[0].length);
}
T[][] matMul(T)(in T[][] a, in T[][] b) pure nothrow
in {
assert
a[0].length == b.length);
} body {
auto result = new T[][](a.length, b[0].length);
auto aux = new T[b.length];
foreach (immutable j; 0 .. b[0].length) {
foreach (immutable k; 0 .. b.length)
aux[k] = b[k][j];
foreach (immutable i; 0 .. a.length)
result[i][j] =
}
return result;
}
Unqual!T[][] transpose(T)(in T[][] m) pure nothrow {
auto r = new Unqual!T[][](m[0].length, m.length);
foreach (immutable nr, row; m)
foreach (immutable nc, immutable c; row)
r[nc][nr] = c;
return r;
}
T norm(T)(in T[][]
return
}
Unqual!T[][] makeUnitVector(T)(in size_t dim) pure nothrow {
auto result = new Unqual!T[][](dim, 1);
foreach (row; result)
row[] = 0;
Line 727 ⟶ 2,497:
/// Return a nxn identity matrix.
Unqual!T[][] matId(T)(in size_t n) pure nothrow {
auto Id = new Unqual!T[][](n, n);
foreach (immutable r, row; Id) {
row[] = 0;
row[r] = 1;
Line 736 ⟶ 2,506:
}
auto B = new
foreach (immutable i, brow; B)
brow[] = A[ma + i][na .. na + brow.length];
return B;
Line 746 ⟶ 2,516:
size_t rows(T)(in T[][] A) pure nothrow { return A.length; }
size_t cols(T)(in T[][] A) pure nothrow {
return A.length ? A[0].length : 0;
Line 751 ⟶ 2,522:
T[][] mcol(T)(in T[][] A, in size_t n) pure nothrow {
return slice2D(A, 0,
}
Line 757 ⟶ 2,528:
in size_t row, in size_t col) pure nothrow {
auto C = new T[][](rows(A), cols(A));
foreach (immutable i, const arow; A)
C[i][] = arow[]; //
foreach (immutable i, const brow; B)
C[row + i][col .. col + brow.length] = brow[];
return C;
Line 766 ⟶ 2,537:
// Main routines ---------------
T[][] makeHouseholder(T)(in T[][] a) {
return
}
Tuple!(T[][],"Q", T[][],"R") QRdecomposition(T)(T[][] A) {
// Work on n columns of A.
foreach (immutable i; 0 .. (m == n ? n - 1 : n)) {
// Select the i-th submatrix. For i=0 this means the original
// matrix A.
// Take the first column of the current submatrix B.
// Create the Householder matrix for the column and embed it
// into an mxm identity.
// The product of all H matrices from the right hand side is
// the orthogonal matrix Q.
Q = Q.matMul(
// The product of all H matrices with A from the LHS is the
// upper triangular matrix R.
A = H.matMul(
}
Line 811 ⟶ 2,582:
/// Solve an upper triangular system by back substitution.
T[][] solveUpperTriangular(T)(in T[][] R, in T[][] b) pure nothrow {
auto x = new T[][](n, 1);
foreach_reverse (immutable k; 0 .. n) {
T tot = 0;
foreach (immutable j; k + 1 .. n)
tot += R[k][j] * x[j][0];
x[k][0] = (b[k][0] - tot) / R[k][k];
Line 825 ⟶ 2,596:
/// Solve a linear least squares problem by QR decomposition.
T[][] lsqr(T)(T[][] A, in T[][] b) pure nothrow {
const qr =
return solveUpperTriangular(
slice2D(qr.R, 0, n - 1, 0, n - 1),
slice2D
}
auto A = new
foreach (immutable i, row; A)
foreach (immutable j, ref item; row)
item = x[0][i] ^^ j;
return lsqr(A,
}
void main() {
immutable qr =
[
[-4.0, 24, -41]]);
writefln(form, qr.Q);
writefln(form, qr.R);
}</
{{out}}
[-0.
[
[0.000, -175.000, 70.000],
[0.000, -0.000, -35.000]]
[[1], [2], [3]]</pre>
=={{header|F_Sharp|F#}}==
<syntaxhighlight lang="fsharp">
// QR decomposition. Nigel Galloway: January 11th., 2022
let n=[[12.0;-51.0;4.0];[6.0;167.0;-68.0];[-4.0;24.0;-41.0]]|>MathNet.Numerics.LinearAlgebra.MatrixExtensions.matrix
let g=n|>MathNet.Numerics.LinearAlgebra.Matrix.qr
printfn $"Matrix\n------\n%A{n}\nQ\n-\n%A{g.Q}\nR\n-\n%A{g.R}"
</syntaxhighlight>
{{out}}
<pre>
Matrix
------
DenseMatrix 3x3-Double
12 -51 4
6 167 -68
-4 24 -41
Q
-
DenseMatrix 3x3-Double
-0.857143 0.394286 -0.331429
-0.428571 -0.902857 0.0342857
0.285714 -0.171429 -0.942857
R
-
DenseMatrix 3x3-Double
-14 -21 14
0 -175 70
0 0 35
</pre>
=={{header|Fortran}}==
{{libheader|LAPACK}}
See the documentation for the [http://www.netlib.org/lapack/lapack-3.1.1/html/dgeqrf.f.html DGEQRF] and [http://www.netlib.org/lapack/lapack-3.1.1/html/dorgqr.f.html DORGQR] routines. Here the example matrix is the magic square from Albrecht Dürer's ''[https://en.wikipedia.org/wiki/Melencolia_I Melencolia I]''.
<syntaxhighlight lang="fortran">program qrtask
implicit none
integer, parameter :: n = 4
real(8) :: durer(n, n) = reshape(dble([ &
16, 5, 9, 4, &
3, 10, 6, 15, &
2, 11, 7, 14, &
13, 8, 12, 1 &
]), [n, n])
real(8) :: q(n, n), r(n, n), qr(n, n), id(n, n), tau(n)
integer, parameter :: lwork = 1024
real(8) :: work(lwork)
integer :: info, i, j
q = durer
call dgeqrf(n, n, q, n, tau, work, lwork, info)
r = 0d0
forall (i = 1:n, j = 1:n, j >= i) r(i, j) = q(i, j)
call dorgqr(n, n, n, q, n, tau, work, lwork, info)
qr = matmul(q, r)
id = matmul(q, transpose(q))
call show(4, durer, "A")
call show(4, q, "Q")
call show(4, r, "R")
call show(4, qr, "Q*R")
call show(4, id, "Q*Q'")
contains
subroutine show(n, a, s)
character(*) :: s
integer :: n, i
real(8) :: a(n, n)
print *, s
do i = 1, n
print 1, a(i, :)
1 format (*(f12.6,:,' '))
end do
end subroutine
end program</syntaxhighlight>
{{out}}
<pre> A
16.000000 3.000000 2.000000 13.000000
5.000000 10.000000 11.000000 8.000000
9.000000 6.000000 7.000000 12.000000
4.000000 15.000000 14.000000 1.000000
Q
-0.822951 0.376971 0.361447 -0.223607
-0.257172 -0.454102 -0.526929 -0.670820
-0.462910 -0.060102 -0.576283 0.670820
-0.205738 -0.805029 0.509510 0.223607
R
-19.442222 -10.904103 -10.595497 -18.516402
0.000000 -15.846152 -15.932298 -0.258437
0.000000 0.000000 -1.974168 -5.922505
0.000000 0.000000 0.000000 -0.000000
Q*R
16.000000 3.000000 2.000000 13.000000
5.000000 10.000000 11.000000 8.000000
9.000000 6.000000 7.000000 12.000000
4.000000 15.000000 14.000000 1.000000
Q*Q'
1.000000 -0.000000 -0.000000 0.000000
-0.000000 1.000000 0.000000 0.000000
-0.000000 0.000000 1.000000 -0.000000
0.000000 0.000000 -0.000000 1.000000</pre>
=={{header|Futhark}}==
<syntaxhighlight lang="futhark">
import "lib/github.com/diku-dk/linalg/linalg"
module linalg_f64 = mk_linalg f64
let eye (n: i32): [n][n]f64 =
let arr = map (\ind -> let (i,j) = (ind/n,ind%n) in if (i==j) then 1.0 else 0.0) (iota (n*n))
in unflatten n n arr
let norm v = linalg_f64.dotprod v v |> f64.sqrt
let qr [n] [m] (a: [m][n]f64): ([m][m]f64, [m][n]f64) =
let make_householder [d] (x: [d]f64): [d][d]f64 =
let div = if x[0] > 0 then x[0] + norm x else x[0] - norm x
let v = map (/div) x
let v[0] = 1
let fac = 2.0 / linalg_f64.dotprod v v
in map2 (map2 (-)) (eye d) (map (map (*fac)) (linalg_f64.outer v v))
let step ((x,y):([m][m]f64,[m][n]f64)) (i:i32): ([m][m]f64,[m][n]f64) =
let h = eye m
let h[i:m,i:m] = make_householder y[i:m,i]
let q': [m][m]f64 = linalg_f64.matmul x h
let a': [m][n]f64 = linalg_f64.matmul h y
in (q',a')
let q = eye m
in foldl step (q,a) (iota n)
entry main = qr [[12.0, -51.0, 4.0],[6.0, 167.0, -68.0],[-4.0, 24.0, -41.0]]
</syntaxhighlight>
{{out}}
<pre>
$ ./qr
[[-0.857143f64, 0.394286f64, -0.331429f64], [-0.428571f64, -0.902857f64, 0.034286f64], [0.285714f64, -0.171429f64, -0.942857f64]]
[[-14.000000f64, -21.000000f64, 14.000000f64], [0.000000f64, -175.000000f64, 70.000000f64], [-0.000000f64, 0.000000f64, 35.000000f64]]
</pre>
=={{header|Go}}==
===Method of task description, library go.matrix===
{{trans|Common Lisp}}
A fairly close port of the Common Lisp solution, this solution uses the [http://github.com/skelterjohn/go.matrix go.matrix library] for supporting functions. Note though, that go.matrix has QR decomposition, as shown in the [[Polynomial_regression#Go|Go solution]] to Polynomial regression. The solution there is coded more directly than by following the CL example here. Similarly, examination of the go.matrix QR source shows that it computes the decomposition more directly.
<syntaxhighlight lang="go">package main
import (
"fmt"
"math"
"github.com/skelterjohn/go.matrix"
)
Line 970 ⟶ 2,895:
}
return x
}</
Output:
<pre>
Line 987 ⟶ 2,912:
3}
</pre>
===Library QR, gonum/matrix===
<syntaxhighlight lang="go">package main
import (
"fmt"
"github.com/gonum/matrix/mat64"
)
func main() {
// task 1: show qr decomp of wp example
a := mat64.NewDense(3, 3, []float64{
12, -51, 4,
6, 167, -68,
-4, 24, -41,
})
var qr mat64.QR
qr.Factorize(a)
var q, r mat64.Dense
q.QFromQR(&qr)
r.RFromQR(&qr)
fmt.Printf("q: %.3f\n\n", mat64.Formatted(&q, mat64.Prefix(" ")))
fmt.Printf("r: %.3f\n\n", mat64.Formatted(&r, mat64.Prefix(" ")))
// task 2: use qr decomp for polynomial regression example
x := []float64{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}
y := []float64{1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}
a = Vandermonde(x, 2)
b := mat64.NewDense(11, 1, y)
qr.Factorize(a)
var f mat64.Dense
f.SolveQR(&qr, false, b)
fmt.Printf("polyfit: %.3f\n",
mat64.Formatted(&f, mat64.Prefix(" ")))
}
func Vandermonde(a []float64, degree int) *mat64.Dense {
x := mat64.NewDense(len(a), degree+1, nil)
for i := range a {
for j, p := 0, 1.; j <= degree; j, p = j+1, p*a[i] {
x.Set(i, j, p)
}
}
return x
}</syntaxhighlight>
{{out}}
<pre>
q: ⎡-0.857 0.394 0.331⎤
⎢-0.429 -0.903 -0.034⎥
⎣ 0.286 -0.171 0.943⎦
r: ⎡ -14.000 -21.000 14.000⎤
⎢ 0.000 -175.000 70.000⎥
⎣ 0.000 0.000 -35.000⎦
polyfit: ⎡1.000⎤
⎢2.000⎥
⎣3.000⎦
</pre>
=={{header|Haskell}}==
Square matrices only; decompose A and solve Rx = q by back substitution
<syntaxhighlight lang="haskell">
import Data.List
import Text.Printf (printf)
eps = 1e-6 :: Double
-- a matrix is represented as a list of columns
mmult :: Num a => [[a]] -> [[a]] -> [[a]]
nth :: Num a => [[a]] -> Int -> Int -> a
mmult_num :: Num a => [[a]] -> a -> [[a]]
madd :: Num a => [[a]] -> [[a]] -> [[a]]
idMatrix :: Num a => Int -> Int -> [[a]]
adjustWithE :: [[Double]] -> Int -> [[Double]]
mmult a b = [ [ sum $ zipWith (*) ak bj | ak <- (transpose a) ] | bj <- b ]
nth mA i j = (mA !! j) !! i
mmult_num mA n = map (\c -> map (*n) c) mA
madd mA mB = zipWith (\c1 c2 -> zipWith (+) c1 c2) mA mB
idMatrix n m = [ [if (i==j) then 1 else 0 | i <- [1..n]] | j <- [1..m]]
adjustWithE mA n = let lA = length mA in
(idMatrix n (n - lA)) ++ (map (\c -> (take (n - lA) (repeat 0.0)) ++ c ) mA)
-- auxiliary functions
sqsum :: Floating a => [a] -> a
norm :: Floating a => [a] -> a
epsilonize :: [[Double]] -> [[Double]]
sqsum a = foldl (\x y -> x + y*y) 0 a
norm a = sqrt $! sqsum a
epsilonize mA = map (\c -> map (\x -> if abs x <= eps then 0 else x) c) mA
-- Householder transformation; householder A = (Q, R)
uTransform :: [Double] -> [Double]
hMatrix :: [Double] -> Int -> Int -> [[Double]]
householder :: [[Double]] -> ([[Double]], [[Double]])
-- householder_rec Q R A
householder_rec :: [[Double]] -> [[Double]] -> Int -> ([[Double]], [[Double]])
uTransform a = let t = (head a) + (signum (head a))*(norm a) in
1 : map (\x -> x/t) (tail a)
hMatrix a n i = let u = uTransform (drop i a) in
madd
(idMatrix (n-i) (n-i))
(mmult_num
(mmult [u] (transpose [u]))
((/) (-2) (sqsum u)))
householder_rec mQ mR 0 = (mQ, mR)
householder_rec mQ mR n = let mSize = length mR in
let mH = adjustWithE (hMatrix (mR!!(mSize - n)) mSize (mSize - n)) mSize in
householder_rec (mmult mQ mH) (mmult mH mR) (n - 1)
householder mA = let mSize = length mA in
let (mQ, mR) = householder_rec (idMatrix mSize mSize) mA mSize in
(epsilonize mQ, epsilonize mR)
backSubstitution :: [[Double]] -> [Double] -> [Double] -> [Double]
backSubstitution mR [] res = res
backSubstitution mR@(hR:tR) q@(h:t) res =
let x = (h / (head hR)) in
backSubstitution
(map tail tR)
(tail (zipWith (-) q (map (*x) hR)))
(x : res)
showMatrix :: [[Double]] -> String
showMatrix mA =
concat $ intersperse "\n"
(map (\x -> unwords $ printf "%10.4f" <$> (x::[Double])) (transpose mA))
mY = [[12, 6, -4], [-51, 167, 24], [4, -68, -41]] :: [[Double]]
q = [21, 245, 35] :: [Double]
main = let (mQ, mR) = householder mY in
putStrLn ("Q: \n" ++ showMatrix mQ) >>
putStrLn ("R: \n" ++ showMatrix mR) >>
putStrLn ("q: \n" ++ show q) >>
putStrLn ("x: \n" ++ show (backSubstitution (reverse (map reverse mR)) (reverse q) []))
</syntaxhighlight>
{{out}}
<pre>
Q:
-0.8571 0.3943 -0.3314
-0.4286 -0.9029 0.0343
0.2857 -0.1714 -0.9429
R:
-14.0000 -21.0000 14.0000
0.0000 -175.0000 70.0000
0.0000 0.0000 35.0000
q:
[21.0,245.0,35.0]
x:
[1.0000000000000004,-0.9999999999999999,1.0]
</pre>
===QR decomposition with Numeric.LinearAlgebra===
<syntaxhighlight lang="haskell">import Numeric.LinearAlgebra
a :: Matrix R
a = (3><3)
[ 12, -51, 4
, 6, 167, -68
, -4, 24, -41]
main = do
print $ qr a
</syntaxhighlight>
{{out}}
<pre>((3><3)
[ -0.8571428571428572, 0.3942857142857143, 0.33142857142857146
, -0.4285714285714286, -0.9028571428571428, -3.428571428571427e-2
, 0.28571428571428575, -0.17142857142857137, 0.9428571428571428 ],(3><3)
[ -14.0, -21.0, 14.000000000000002
, 0.0, -175.00000000000003, 70.00000000000001
, 0.0,
</pre>
=={{header|J}}==
'''Solution''' (built-in):<syntaxhighlight lang="j"> QR =: 128!:0</syntaxhighlight>
'''Solution''' (custom implementation): <syntaxhighlight lang="j"> mp=: +/ . * NB. matrix product
h =: +@|: NB. conjugate transpose
QR=: 3 : 0
n=.{:$A=.y
if. 1>:n do.
A ((% {.@,) ; ]) %:(h A) mp A
else.
m =.>.n%2
A0=.m{."1 A
A1=.m}."1 A
'Q0 R0'=.QR A0
'Q1 R1'=.QR A1 - Q0 mp T=.(h Q0) mp A1
(Q0,.Q1);(R0,.T),(-n){."1 R1
end.
)</syntaxhighlight>
'''Example''': <syntaxhighlight lang="j"> QR 12 _51 4,6 167 _68,:_4 24 _41
+-----------------------------+----------+
| 0.857143 _0.394286 _0.331429|14 21 _14|
| 0.428571 0.902857 0.0342857| 0 175 _70|
|_0.285714 0.171429 _0.942857| 0 0 35|
+-----------------------------+----------+</syntaxhighlight>
'''Example''' (polynomial fitting using QR reduction):<syntaxhighlight lang="j"> X=:i.# Y=:1 6 17 34 57 86 121 162 209 262 321
'Q R'=: QR X ^/ i.3
R %.~(|:Q)+/ .* Y
1 2 3</
'''Notes''':J offers a built-in QR decomposition function, <tt>128!:0</tt>. If J did not offer this function as a built-in, it could be written in J along the lines of the second version, which is covered in [[j:Essays/QR Decomposition|an essay on the J wiki]].
=={{header|Java}}==
=== JAMA ===
Using the [https://math.nist.gov/javanumerics/jama/ JAMA] library. Compile with: '''javac -cp Jama-1.0.3.jar Decompose.java'''.
<syntaxhighlight lang="java">import Jama.Matrix;
import Jama.QRDecomposition;
public class Decompose {
public static void main(String[] args) {
var matrix = new Matrix(new double[][] {
{12, -51, 4},
{ 6, 167, -68},
{-4, 24, -41},
});
var qr = new QRDecomposition(matrix);
qr.getQ().print(10, 4);
qr.getR().print(10, 4);
}
}</syntaxhighlight>
{{out}}
<pre> -0.8571 0.3943 -0.3314
-0.4286 -0.9029 0.0343
0.2857 -0.1714 -0.9429
-14.0000 -21.0000 14.0000
0.0000 -175.0000 70.0000
0.0000 0.0000 35.0000</pre>
=== Colt ===
Using the [https://dst.lbl.gov/ACSSoftware/colt/ Colt] library. Compile with: '''javac -cp colt.jar Decompose.java'''.
<syntaxhighlight lang="java">import cern.colt.matrix.impl.DenseDoubleMatrix2D;
import cern.colt.matrix.linalg.QRDecomposition;
public class Decompose {
public static void main(String[] args) {
var a = new DenseDoubleMatrix2D(new double[][] {
{12, -51, 4},
{ 6, 167, -68},
{-4, 24, -41}
});
var qr = new QRDecomposition(a);
System.out.println(qr.getQ());
System.out.println();
System.out.println(qr.getR());
}
}</syntaxhighlight>
{{out}}
<pre>3 x 3 matrix
-0.857143 0.394286 -0.331429
-0.428571 -0.902857 0.034286
0.285714 -0.171429 -0.942857
3 x 3 matrix
-14 -21 14
0 -175 70
0 0 35</pre>
=== Apache Commons Math ===
Using the Apache Commons [http://commons.apache.org/proper/commons-math/ Math] library.
Compile with: '''javac -cp commons-math3-3.6.1.jar Decompose.java'''.
<syntaxhighlight lang="java">import java.util.Locale;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.QRDecomposition;
import org.apache.commons.math3.linear.RealMatrix;
public class Decompose {
public static void main(String[] args) {
var a = new Array2DRowRealMatrix(new double[][] {
{12, -51, 4},
{ 6, 167, -68},
{-4, 24, -41}
});
var qr = new QRDecomposition(a);
print(qr.getQ());
System.out.println();
print(qr.getR());
}
public static void print(RealMatrix a) {
for (double[] u: a.getData()) {
System.out.print("[ ");
for (double x: u) {
System.out.printf(Locale.ROOT, "%10.4f ", x);
}
System.out.println("]");
}
}
}</syntaxhighlight>
{{out}}
<pre>[ -0.8571 0.3943 -0.3314 ]
[ -0.4286 -0.9029 0.0343 ]
[ 0.2857 -0.1714 -0.9429 ]
[ -14.0000 -21.0000 14.0000 ]
[ 0.0000 -175.0000 70.0000 ]
[ 0.0000 0.0000 35.0000 ]</pre>
=== la4j ===
Using the [http://la4j.org/ la4j] library. Compile with: '''javac -cp la4j-0.6.0.jar Decompose.java'''.
<syntaxhighlight lang="java">import org.la4j.Matrix;
import org.la4j.decomposition.QRDecompositor;
public class Decompose {
public static void main(String[] args) {
var a = Matrix.from2DArray(new double[][] {
{12, -51, 4},
{ 6, 167, -68},
{-4, 24, -41},
});
Matrix[] qr = new QRDecompositor(a).decompose();
System.out.println(qr[0]);
System.out.println(qr[1]);
}
}</syntaxhighlight>
{{out}}
<pre>-0,857 0,394 -0,331
-0,429 -0,903 0,034
0,286 -0,171 -0,943
-14,000 -21,000 14,000
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}}==
Built-in function
<syntaxhighlight lang="julia">Q, R = qr([12 -51 4; 6 167 -68; -4 24 -41])</syntaxhighlight>
{{out}}
<pre>
(
3x3 Array{Float64,2}:
-0.857143 0.394286 0.331429
-0.428571 -0.902857 -0.0342857
0.285714 -0.171429 0.942857 ,
3x3 Array{Float64,2}:
-14.0 -21.0 14.0
0.0 -175.0 70.0
0.0 0.0 -35.0)
</pre>
=={{header|Maple}}==
<syntaxhighlight lang="maple">with(LinearAlgebra):
A:=<12,-51,4;6,167,-68;-4,24,-41>:
Q,R:=QRDecomposition(A):
Q;
R;</syntaxhighlight>
{{out}}
<pre> [ -69 -58 ]
[6/7 --- --- ]
[ 175 175 ]
[ ]
[ 158 ]
[3/7 --- 6/175]
[ 175 ]
[ ]
[ -33 ]
[-2/7 6/35 --- ]
[ 35 ]
[14 21 -14]
[ ]
[ 0 175 -70]
[ ]
[ 0 0 35]</pre>
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<
q//MatrixForm
Line 1,035 ⟶ 3,812:
-> 14 21 -14
0 175 -70
0 0 35</
=={{header|MATLAB}} / {{header|Octave}}==
<
6 167 -68
-4 24 -41];
[Q,R]=qr(A) </
Output:
<pre>Q =
Line 1,054 ⟶ 3,831:
0 175 -70
0 0 35</pre>
=={{header|Maxima}}==
<syntaxhighlight lang="maxima">load(lapack)$ /* This may hang up in wxMaxima, if this happens, use xMaxima or plain Maxima in a terminal */
a: matrix([12, -51, 4],
[ 6, 167, -68],
[-4, 24, -41])$
[q, r]: dgeqrf(a)$
mat_norm(q . r - a, 1);
4.2632564145606011E-14
/* Note: the lapack package is a lisp translation of the fortran lapack library */</syntaxhighlight>
For an exact or arbitrary precision solution:<syntaxhighlight lang="maxima">load("linearalgebra")$
load("eigen")$
unitVector(n) := ematrix(n,1,1,1,1);
signValue(r) := block([s:sign(r)],
if s='pos then 1 else if s='zero then 0 else -1);
householder(a) := block([m : length(a),u,v,beta],
u : a + sqrt(a . a)*signValue(a[1,1])*unitVector(m),
v : u / u[1,1],
beta : 2/(v . v),
diagmatrix(m,1) - beta*transpose(v . transpose(v)));
getSubmatrix(obj,i1,j1,i2,j2) :=
genmatrix(lambda([i,j], obj[i+i1-1,j+j1-1]),i2-i1+1,j2-j1+1);
setSubmatrix(obj,i1,j1,subobj) := block([m,n],
[m,n] : matrix_size(subobj),
for i: 0 thru m-1 do
(for j: 0 thru n-1 do
obj[i1+i,j1+j] : subobj[i+1,j+1]));
qr(obj) := block([m,n,qm,rm,i],
[m,n] : matrix_size(obj),
qm : diagmatrix(m,1),
rm : copymatrix(obj),
for i: 1 thru (if m=n then n-1 else n) do
block([x,h],
x : getSubmatrix(rm,i,i,m,i),
h : diagmatrix(m,1),
setSubmatrix(h,i,i,householder(x)),
qm : qm . h,
rm : h . rm),
[qm,rm]);
solveUpperTriangular(r,b) := block([n,x,index,k],
n : second(matrix_size(r)),
x : genmatrix(lambda([a, b], 0), n, 1),
for k: n thru 1 step -1 do
(index : min(n,k+1),
x[k,1] : (b[k,1] - (getSubmatrix(r,k,index,k,n) . getSubmatrix(x,index,1,n,1)))/r[k,k]),
x);
lsqr(a,b) := block([q,r,n],
[q,r] : qr(a),
n : second(matrix_size(r)),
solveUpperTriangular(getSubmatrix(r,1,1,n,n), transpose(q) . b));
polyfit(x,y,n) := block([a,j],
a : genmatrix(lambda([i,j], if j=1 then 1.0b0 else bfloat(x[i,1]^(j-1))),
length(x),n+1),
lsqr(a,y));</syntaxhighlight>Then we have the examples:<syntaxhighlight lang="maxima">(%i) [q,r] : qr(a);
[ 6 69 58 ]
[ - - --- --- ]
[ 7 175 175 ]
[ ] [ - 14 - 21 14 ]
[ 3 158 6 ] [ ]
(%o) [[ - - - --- - --- ], [ 0 - 175 70 ]]
[ 7 175 175 ] [ ]
[ ] [ 0 0 - 35 ]
[ 2 6 33 ]
[ - - -- -- ]
[ 7 35 35 ]
(%i) mat_norm(q . r - a, 1);
(%o) 0
(%i) x : transpose(matrix([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]))$
(%i) y : transpose(matrix([1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]))$
(%i) fpprec : 30$
(%i) polyfit(x, y, 2);
[ 9.99999999999999999999999999996b-1 ]
[ ]
(%o) [ 2.00000000000000000000000000002b0 ]
[ ]
[ 3.0b0 ]</syntaxhighlight>
=={{header|Nim}}==
{{trans|Python}}
{{libheader|arraymancer}}
The library “arraymancer” provides a function “qr” to get the QR decomposition. Using the Tensor type of “arraymancer” we propose here an implementation of this decomposition adapted from the Python version.
<syntaxhighlight lang="nim">import math, strformat, strutils
import arraymancer
####################################################################################################
# First part: QR decomposition.
proc eye(n: Positive): Tensor[float] =
## Return the (n, n) identity matrix.
result = newTensor[float](n.int, n.int)
for i in 0..<n: result[i, i] = 1
proc norm(v: Tensor[float]): float =
## return the norm of a vector.
assert v.shape.len == 1
result = sqrt(dot(v, v)) * sgn(v[0]).toFloat
proc houseHolder(a: Tensor[float]): Tensor[float] =
## return the house holder of vector "a".
var v = a / (a[0] + norm(a))
v[0] = 1
result = eye(a.shape[0]) - (2 / dot(v, v)) * (v.unsqueeze(1) * v.unsqueeze(0))
proc qrDecomposition(a: Tensor): tuple[q, r: Tensor] =
## Return the QR decomposition of matrix "a".
assert a.shape.len == 2
let m = a.shape[0]
let n = a.shape[1]
result.q = eye(m)
result.r = a.clone
for i in 0..<(n - ord(m == n)):
var h = eye(m)
h[i..^1, i..^1] = houseHolder(result.r[i..^1, i].squeeze(1))
result.q = result.q * h
result.r = h * result.r
####################################################################################################
# Second part: polynomial regression example.
proc lsqr(a, b: Tensor[float]): Tensor[float] =
let (q, r) = a.qrDecomposition()
let n = r.shape[1]
result = solve(r[0..<n, _], (q.transpose() * b)[0..<n])
proc polyfit(x, y: Tensor[float]; n: int): Tensor[float] =
var z = newTensor[float](x.shape[0], n + 1)
var t = x.reshape(x.shape[0], 1)
for i in 0..n: z[_, i] = t^.i.toFloat
result = lsqr(z, y.transpose())
#———————————————————————————————————————————————————————————————————————————————————————————————————
proc printMatrix(a: Tensor) =
var str: string
for i in 0..<a.shape[0]:
let start = str.len
for j in 0..<a.shape[1]:
str.addSep(" ", start)
str.add &"{a[i, j]:8.3f}"
str.add '\n'
stdout.write str
proc printVector(a: Tensor) =
var str: string
for i in 0..<a.shape[0]:
str.addSep(" ")
str.add &"{a[i]:4.1f}"
echo str
let mat = [[12, -51, 4],
[ 6, 167, -68],
[-4, 24, -41]].toTensor.astype(float)
let (q, r) = mat.qrDecomposition()
echo "Q:"
printMatrix q
echo "R:"
printMatrix r
echo()
let x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10].toTensor.astype(float)
let y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321].toTensor.astype(float)
echo "polyfit:"
printVector polyfit(x, y, 2)</syntaxhighlight>
{{out}}
<pre>Q:
-0.857 0.394 0.331
-0.429 -0.903 -0.034
0.286 -0.171 0.943
R:
-14.000 -21.000 14.000
0.000 -175.000 70.000
0.000 -0.000 -35.000
polyfit:
1.0 2.0 3.0</pre>
=={{header|PARI/GP}}==
{{works with|PARI/GP|2.6.0 and above}}
<syntaxhighlight lang="parigp">matqr(M)</syntaxhighlight>
=={{header|Perl}}==
Letting the <code>PDL</code> module do all the work.
<syntaxhighlight lang="perl">use strict;
use warnings;
use PDL;
use PDL::LinearAlgebra qw(mqr);
my $a = pdl(
[12, -51, 4],
[ 6, 167, -68],
[-4, 24, -41],
[-1, 1, 0],
[ 2, 0, 3]
);
my ($q, $r) = mqr($a);
print $q, $r, $q x $r;</syntaxhighlight>
{{out}}
<pre>[
[ -0.84641474 0.39129081 -0.34312406]
[ -0.42320737 -0.90408727 0.029270162]
[ 0.28213825 -0.17042055 -0.93285599]
[ 0.070534562 -0.014040652 0.001099372]
[ -0.14106912 0.016655511 0.10577161]
]
[
[-14.177447 -20.666627 13.401567]
[ 0 -175.04254 70.080307]
[ 0 0 35.201543]
]
[
[ 12 -51 4]
[ 6 167 -68]
[ -4 24 -41]
[ -1 1 0]
[ 2 0 3]
]</pre>
=={{header|Phix}}==
using matrix_mul() from [[Matrix_multiplication#Phix]]
and matrix_transpose() from [[Matrix_transposition#Phix]]
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #000080;font-style:italic;">-- demo/rosettacode/QRdecomposition.exw</span>
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<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>
<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>
<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>
<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>
<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>
<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>
<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>
<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>
<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;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">c</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<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>
<span style="color: #000080;font-style:italic;">-- transpose a vector of length m into an mx1 matrix,
-- eg {1,2,3} -> <nowiki>{{</nowiki>1},{2},{3<nowiki>}}</nowiki></span>
<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>
<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>
<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: #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>
<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>
<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>
<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>
<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: #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>
<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>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<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>
<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: #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>
<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>
<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>
<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: #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>
<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>
<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>
<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>
<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>
<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>
<span style="color: #000080;font-style:italic;">--
-- Programming note: The code of this main loop was not as easily
-- written as the first glance might suggest. Explicitly setting
-- to 0 any a[i,j] [etc] that should be 0 but have inadvertently
-- gotten set to +/-1e-15 or thereabouts may be advisable. The
-- commented-out code was retrieved from a backup and should be
-- treated as an example and not be trusted (iirc, it made no
-- difference to the test cases used, so I deleted it, and then
-- had second thoughts about it a few days later).
--</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: #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>
<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>
<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>
<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>
<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>
<span style="color: #000080;font-style:italic;">-- for row=j+1 to length(a) do
-- a[row][j] = 0
-- end 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>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<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(>=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>
<!--</syntaxhighlight>-->
{{out}}
<pre>
"A"
{{ 12, -51, 4},
{ 6, 167, -68},
{ -4, 24, -41}}
"Q"
{{0.85714,-0.3943,0.33143},
{0.42857,0.90286,-0.0343},
{-0.2857,0.17143,0.94286}}
"R"
{{ 14, 21, -14},
{ 0, 175, -70},
{ 0, 0, -35}}
"Q * R"
{{ 12, -51, 4},
{ 6, 167, -68},
{ -4, 24, -41}}
"Q * Q`"
{{ 1, 0, 0},
{ 0, 1, 0},
{ 0, 0, 1}}
Least-squares solution:
{1,2,3}
</pre>
=={{header|PowerShell}}==
<syntaxhighlight lang="powershell">
function qr([double[][]]$A) {
$m,$n = $A.count, $A[0].count
$pm,$pn = ($m-1), ($n-1)
[double[][]]$Q = 0..($m-1) | foreach{$row = @(0) * $m; $row[$_] = 1; ,$row}
[double[][]]$R = $A | foreach{$row = $_; ,@(0..$pn | foreach{$row[$_]})}
foreach ($h in 0..$pn) {
[double[]]$u = $R[$h..$pm] | foreach{$_[$h]}
[double]$nu = $u | foreach {[double]$sq = 0} {$sq += $_*$_} {[Math]::Sqrt($sq)}
$u[0] -= if ($u[0] -lt 0) {$nu} else {-$nu}
[double]$nu = $u | foreach {$sq = 0} {$sq += $_*$_} {[Math]::Sqrt($sq)}
[double[]]$u = $u | foreach { $_/$nu}
[double[][]]$v = 0..($u.Count - 1) | foreach{$i = $_; ,($u | foreach{2*$u[$i]*$_})}
[double[][]]$CR = $R | foreach{$row = $_; ,@(0..$pn | foreach{$row[$_]})}
[double[][]]$CQ = $Q | foreach{$row = $_; ,@(0..$pm | foreach{$row[$_]})}
foreach ($i in $h..$pm) {
foreach ($j in $h..$pn) {
$R[$i][$j] -= $h..$pm | foreach {[double]$sum = 0} {$sum += $v[$i-$h][$_-$h]*$CR[$_][$j]} {$sum}
}
}
if (0 -eq $h) {
foreach ($i in $h..$pm) {
foreach ($j in $h..$pm) {
$Q[$i][$j] -= $h..$pm | foreach {$sum = 0} {$sum += $v[$i][$_]*$CQ[$_][$j]} {$sum}
}
}
} else {
$p = $h-1
foreach ($i in $h..$pm) {
foreach ($j in 0..$p) {
$Q[$i][$j] -= $h..$pm | foreach {$sum = 0} {$sum += $v[$i-$h][$_-$h]*$CQ[$_][$j]} {$sum}
}
foreach ($j in $h..$pm) {
$Q[$i][$j] -= $h..$pm | foreach {$sum = 0} {$sum += $v[$i-$h][$_-$h]*$CQ[$_][$j]} {$sum}
}
}
}
}
foreach ($i in 0..$pm) {
foreach ($j in $i..$pm) {$Q[$i][$j],$Q[$j][$i] = $Q[$j][$i],$Q[$i][$j]}
}
[PSCustomObject]@{"Q" = $Q; "R" = $R}
}
function leastsquares([Double[][]]$A,[Double[]]$y) {
$QR = qr $A
[Double[][]]$Q = $QR.Q
[Double[][]]$R = $QR.R
$m,$n = $A.count, $A[0].count
[Double[]]$z = foreach ($j in 0..($m-1)) {
0..($m-1) | foreach {$sum = 0} {$sum += $Q[$_][$j]*$y[$_]} {$sum}
}
[Double[]]$x = @(0)*$n
for ($i = $n-1; $i -ge 0; $i--) {
for ($j = $i+1; $j -lt $n; $j++) {
$z[$i] -= $x[$j]*$R[$i][$j]
}
$x[$i] = $z[$i]/$R[$i][$i]
}
$x
}
function polyfit([Double[]]$x,[Double[]]$y,$n) {
$m = $x.Count
[Double[][]]$A = 0..($m-1) | foreach{$row = @(1) * ($n+1); ,$row}
for ($i = 0; $i -lt $m; $i++) {
for ($j = $n-1; 0 -le $j; $j--) {
$A[$i][$j] = $A[$i][$j+1]*$x[$i]
}
}
leastsquares $A $y
}
function show($m) {$m | foreach {write-host "$_"}}
$A = @(@(12,-51,4), @(6,167,-68), @(-4,24,-41))
$x = @(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
$y = @(1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321)
$QR = qr $A
$ps = (polyfit $x $y 2)
"Q = "
show $QR.Q
"R = "
show $QR.R
"polyfit "
"X^2 X constant"
"$(polyfit $x $y 2)"
</syntaxhighlight>
{{out}}
<pre>
Q =
-0.857142857142857 0.394285714285714 -0.331428571428571
-0.428571428571429 -0.902857142857143 0.0342857142857143
0.285714285714286 -0.171428571428571 -0.942857142857143
R =
-14 -21 14
8.88178419700125E-16 -175 70
-4.44089209850063E-16 0 35
polyfit
X^2 X constant
3 1.99999999999998 1.00000000000005
</pre>
=={{header|Python}}==
{{libheader|NumPy}}
Numpy has a qr function but here is a reimplementation to show construction and use of the Householder reflections.
<syntaxhighlight lang="python">#!/usr/bin/env python3
import numpy as np
def qr(A):
m, n = A.shape
Q = np.eye(m)
for i in range(n - (m == n)):
H = np.eye(m)
H[i:, i:] = make_householder(A[i:, i])
Q = np.dot(Q, H)
A = np.dot(H, A)
return Q, A
def make_householder(a):
v = a / (a[0] + np.copysign(np.linalg.norm(a), a[0]))
v[0] = 1
H = np.eye(a.shape[0])
H -= (2 / np.dot(v, v)) * np.dot(v[:, None], v[None, :])
return H
# task 1: show qr decomp of wp example
a = np.array(((
(12, -51, 4),
( 6, 167, -68),
(-4, 24, -41),
)))
q, r = qr(a)
print('q:\n', q.round(6))
print('r:\n', r.round(6))
# task 2: use qr decomp for polynomial regression example
def polyfit(x, y, n):
return lsqr(x[:, None]**np.arange(n + 1), y.T)
def lsqr(a, b):
q, r = qr(a)
_, n = r.shape
return np.linalg.solve(r[:n, :], np.dot(q.T, b)[:n])
x = np.array((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10))
y = np.array((1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321))
print('\npolyfit:\n', polyfit(x, y, 2))</syntaxhighlight>
{{out}}
<pre>
q:
[[-0.857143 0.394286 0.331429]
[-0.428571 -0.902857 -0.034286]
[ 0.285714 -0.171429 0.942857]]
r:
[[ -14. -21. 14.]
[ 0. -175. 70.]
[ 0. 0. -35.]]
polyfit:
[ 1. 2. 3.]
</pre>
=={{header|R}}==
<
a <- matrix(c(12, -51, 4, 6, 167, -68, -4, 24, -41), nrow=3, ncol=3, byrow=T)
Line 1,078 ⟶ 4,433:
xx <- x*x
m <- lm(y ~ x + xx)
coef(m)</
=={{header|Racket}}==
Racket has QR-decomposition builtin:
<syntaxhighlight lang="racket">
> (require math)
> (matrix-qr (matrix [[12 -51 4]
[ 6 167 -68]
[-4 24 -41]]))
(array #[#[6/7 -69/175 -58/175] #[3/7 158/175 6/175] #[-2/7 6/35 -33/35]])
(array #[#[14 21 -14] #[0 175 -70] #[0 0 35]])
</syntaxhighlight>
The builtin QR-decomposition uses the Gram-Schmidt algorithm.
Here is an implementation of the Householder method:
<syntaxhighlight lang="racket">
#lang racket
(require math/matrix math/array)
(define-values (T I col size)
(values ; short names
matrix-transpose identity-matrix matrix-col matrix-num-rows))
(define (scale c A) (matrix-scale A c))
(define (unit n i) (build-matrix n 1 (λ (j _) (if (= j i) 1 0))))
(define (H u)
(matrix- (I (size u))
(scale (/ 2 (matrix-dot u u))
(matrix* u (T u)))))
(define (normal a)
(define a0 (matrix-ref a 0 0))
(matrix- a (scale (* (sgn a0) (matrix-2norm a))
(unit (size a) 0))))
(define (QR A)
(define n (size A))
(for/fold ([Q (I n)] [R A]) ([i (- n 1)])
(define Hi (H (normal (submatrix R (:: i n) (:: i (+ i 1))))))
(define Hi* (if (= i 0) Hi (block-diagonal-matrix (list (I i) Hi))))
(values (matrix* Q Hi*) (matrix* Hi* R))))
(QR (matrix [[12 -51 4]
[ 6 167 -68]
[-4 24 -41]]))
</syntaxhighlight>
Output:
<syntaxhighlight lang="racket">
(array #[#[6/7 69/175 -58/175]
#[3/7 -158/175 6/175]
#[-2/7 -6/35 -33/35]])
(array #[#[14 21 -14]
#[0 -175 70]
#[0 0 35]])
</syntaxhighlight>
=={{header|Raku}}==
(formerly Perl 6)
{{Works with|rakudo|2018.06}}
<syntaxhighlight lang="raku" line># sub householder translated from https://codereview.stackexchange.com/questions/120978/householder-transformation
use v6;
sub identity(Int:D $m --> Array of Array) {
my Array @M;
for 0 ..^ $m -> $i {
@M.push: [0 xx $m];
@M[$i; $i] = 1;
}
@M;
}
multi multiply(Array:D @A, @b where Array:D --> Array) {
my @c;
for ^@A X ^@b -> ($i, $j) {
@c[$i] += @A[$i; $j] * @b[$j];
}
@c;
}
multi multiply(Array:D @A, Array:D @B --> Array of Array) {
my Array @C;
for ^@A X ^@B[0] -> ($i, $j) {
@C[$i; $j] += @A[$i; $_] * @B[$_; $j] for ^@B;
}
@C;
}
sub transpose(Array:D @M --> Array of Array) {
my ($rows, $cols) = (@M.elems, @M[0].elems);
my Array @T;
for ^$cols X ^$rows -> ($j, $i) {
@T[$j; $i] = @M[$i; $j];
}
@T;
}
####################################################
# NOTE: @A gets overwritten and becomes @R, only need
# to return @Q.
####################################################
sub householder(Array:D @A --> Array) {
my Int ($m, $n) = (@A.elems, @A[0].elems);
my @v = 0 xx $m;
my Array @Q = identity($m);
for 0 ..^ $n -> $k {
my Real $sum = 0;
my Real $A0 = @A[$k; $k];
my Int $sign = $A0 < 0 ?? -1 !! 1;
for $k ..^ $m -> $i {
$sum += @A[$i; $k] * @A[$i; $k];
}
my Real $sqr_sum = $sign * sqrt($sum);
my Real $tmp = sqrt(2 * ($sum + $A0 * $sqr_sum));
@v[$k] = ($sqr_sum + $A0) / $tmp;
for ($k + 1) ..^ $m -> $i {
@v[$i] = @A[$i; $k] / $tmp;
}
for 0 ..^ $n -> $j {
$sum = 0;
for $k ..^ $m -> $i {
$sum += @v[$i] * @A[$i; $j];
}
for $k ..^ $m -> $i {
@A[$i; $j] -= 2 * @v[$i] * $sum;
}
}
for 0 ..^ $m -> $j {
$sum = 0;
for $k ..^ $m -> $i {
$sum += @v[$i] * @Q[$i; $j];
}
for $k ..^ $m -> $i {
@Q[$i; $j] -= 2 * @v[$i] * $sum;
}
}
}
@Q
}
sub dotp(@a where Array:D, @b where Array:D --> Real) {
[+] @a >>*<< @b;
}
sub upper-solve(Array:D @U, @b where Array:D, Int:D $n --> Array) {
my @y = 0 xx $n;
@y[$n - 1] = @b[$n - 1] / @U[$n - 1; $n - 1];
for reverse ^($n - 1) -> $i {
@y[$i] = (@b[$i] - (dotp(@U[$i], @y))) / @U[$i; $i];
}
@y;
}
sub polyfit(@x where Array:D, @y where Array:D, Int:D $n) {
my Int $m = @x.elems;
my Array @V;
# Vandermonde matrix
for ^$m X (0 .. $n) -> ($i, $j) {
@V[$i; $j] = @x[$i] ** $j
}
# least squares
my $Q = householder(@V);
my @b = multiply($Q, @y);
return upper-solve(@V, @b, $n + 1);
}
sub print-mat(Array:D @M, Str:D $name) {
my Int ($m, $n) = (@M.elems, @M[0].elems);
print "\n$name:\n";
for 0 ..^ $m -> $i {
for 0 ..^ $n -> $j {
print @M[$i; $j].fmt("%12.6f ");
}
print "\n";
}
}
sub MAIN() {
############
# 1st part #
############
my Array @A = (
[12, -51, 4],
[ 6, 167, -68],
[-4, 24, -41],
[-1, 1, 0],
[ 2, 0, 3]
);
print-mat(@A, 'A');
my $Q = householder(@A);
$Q = transpose($Q);
print-mat($Q, 'Q');
# after householder, @A is now @R
print-mat(@A, 'R');
print-mat(multiply($Q, @A), 'check Q x R = A');
############
# 2nd part #
############
my @x = [^11];
my @y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321];
my @coef = polyfit(@x, @y, 2);
say
"\npolyfit:\n",
<constant X X^2>.fmt("%12s"),
"\n",
@coef.fmt("%12.6f");
}</syntaxhighlight>
output:
<pre>
A:
12.000000 -51.000000 4.000000
6.000000 167.000000 -68.000000
-4.000000 24.000000 -41.000000
-1.000000 1.000000 0.000000
2.000000 0.000000 3.000000
Q:
-0.846415 0.391291 -0.343124 0.066137 -0.091462
-0.423207 -0.904087 0.029270 0.017379 -0.048610
0.282138 -0.170421 -0.932856 -0.021942 0.143712
0.070535 -0.014041 0.001099 0.997401 0.004295
-0.141069 0.016656 0.105772 0.005856 0.984175
R:
-14.177447 -20.666627 13.401567
-0.000000 -175.042539 70.080307
0.000000 0.000000 35.201543
-0.000000 0.000000 0.000000
0.000000 -0.000000 0.000000
check Q x R = A:
12.000000 -51.000000 4.000000
6.000000 167.000000 -68.000000
-4.000000 24.000000 -41.000000
-1.000000 1.000000 -0.000000
2.000000 -0.000000 3.000000
polyfit:
constant X X^2
1.000000 2.000000 3.000000
</pre>
=={{header|Rascal}}==
[[File:Qrresult.jpeg||200px|thumb|right]]
This function applies the Gram Schmidt algorithm. Q is printed in the console, R can be printed or visualized.
<syntaxhighlight lang="rascal">import util::Math;
import Prelude;
import vis::Figure;
import vis::Render;
public rel[real,real,real] QRdecomposition(rel[real x, real y, real v] matrix){
//orthogonalcolumns
oc = domainR(matrix, {0.0});
for (x <- sort(toList(domain(matrix)-{0.0}))){
c = domainR(matrix, {x});
o = domainR(oc, {x-1});
for (n <- [1.0 .. x]){
o = domainR(oc, {n-1});
c = matrixSubtract(c, matrixMultiplybyN(o, matrixDotproduct(o, c)/matrixDotproduct(o, o)));
}
oc += c;
}
Q = {};
//from orthogonal to orthonormal columns
for (el <- oc){
c = domainR(oc, {el[0]});
Q += matrixNormalize({el}, c);
}
//from Q to R
R= matrixMultiplication(matrixTranspose(Q), matrix);
R= {<x,y,toReal(round(v))> | <x,y,v> <- R};
println("Q:");
iprintlnExp(Q);
println();
println("R:");
return R;
}
//a function that takes the transpose of a matrix, see also Rosetta Code problem "Matrix transposition"
public rel[real, real, real] matrixTranspose(rel[real x, real y, real v] matrix){
return {<y, x, v> | <x, y, v> <- matrix};
}
//a function to normalize an element of a matrix by the normalization of a column
public rel[real,real,real] matrixNormalize(rel[real x, real y, real v] element, rel[real x, real y, real v] column){
normalized = 1.0/nroot((0.0 | it + v*v | <x,y,v> <- column), 2);
return matrixMultiplybyN(element, normalized);
}
//a function that takes the dot product, see also Rosetta Code problem "Dot product"
public real matrixDotproduct(rel[real x, real y, real v] column1, rel[real x, real y, real v] column2){
return (0.0 | it + v1*v2 | <x1,y1,v1> <- column1, <x2,y2,v2> <- column2, y1==y2);
}
//a function to subtract two columns
public rel[real,real,real] matrixSubtract(rel[real x, real y, real v] column1, rel[real x, real y, real v] column2){
return {<x1,y1,v1-v2> | <x1,y1,v1> <- column1, <x2,y2,v2> <- column2, y1==y2};
}
//a function to multiply a column by a number
public rel[real,real,real] matrixMultiplybyN(rel[real x, real y, real v] column, real n){
return {<x,y,v*n> | <x,y,v> <- column};
}
//a function to perform matrix multiplication, see also Rosetta Code problem "Matrix multiplication".
public rel[real, real, real] matrixMultiplication(rel[real x, real y, real v] matrix1, rel[real x, real y, real v] matrix2){
if (max(matrix1.x) == max(matrix2.y)){
p = {<x1,y1,x2,y2, v1*v2> | <x1,y1,v1> <- matrix1, <x2,y2,v2> <- matrix2};
result = {};
for (y <- matrix1.y){
for (x <- matrix2.x){
v = (0.0 | it + v | <x1, y1, x2, y2, v> <- p, x==x2 && y==y1, x1==y2 && y2==x1);
result += <x,y,v>;
}
}
return result;
}
else throw "Matrix sizes do not match.";
}
// a function to visualize the result
public void displayMatrix(rel[real x, real y, real v] matrix){
points = [box(text("<v>"), align(0.3333*(x+1),0.3333*(y+1)),shrink(0.25)) | <x,y,v> <- matrix];
render(overlay([*points], aspectRatio(1.0)));
}
//a matrix, given by a relation of <x-coordinate, y-coordinate, value>.
public rel[real x, real y, real v] matrixA = {
<0.0,0.0,12.0>, <0.0,1.0, 6.0>, <0.0,2.0,-4.0>,
<1.0,0.0,-51.0>, <1.0,1.0,167.0>, <1.0,2.0,24.0>,
<2.0,0.0,4.0>, <2.0,1.0,-68.0>, <2.0,2.0,-41.0>
};</syntaxhighlight>
Example using visualization
<pre>rascal>displayMatrix(QRdecomposition(matrixA))
Q:
{
<1.0,0.0,-0.394285714285714285714285714285714285714285714285714285714285714285713300>,
<0.0,0.0,0.857142857142857142857142857142857142857142857142857142857142857142840>,
<0.0,1.0,0.428571428571428571428571428571428571428571428571428571428571428571420>,
<0.0,2.0,-0.285714285714285714285714285714285714285714285714285714285714285714280>,
<2.0,0.0,-0.33142857142857142857142857142857142857142857142857142857142857142858800>,
<1.0,2.0,0.171428571428571428571428571428571428571428571428571428571428571428571000>,
<2.0,2.0,-0.94285714285714285714285714285714285714285714285714285714285714285719000>,
<1.0,1.0,0.902857142857142857142857142857142857142857142857142857142857142857140600>,
<2.0,1.0,0.03428571428571428571428571428571428571428571428571428571428571428571600>
}
See R in picture</pre>
=={{header|SAS}}==
<syntaxhighlight lang="sas">/* See http://support.sas.com/documentation/cdl/en/imlug/63541/HTML/default/viewer.htm#imlug_langref_sect229.htm */
proc iml;
a={12 -51 4,6 167 -68,-4 24 -41};
print(a);
call qr(q,r,p,d,a);
print(q);
print(r);
quit;
/*
a
12 -51 4
6 167 -68
-4 24 -41
q
-0.857143 0.3942857 -0.331429
-0.428571 -0.902857 0.0342857
0.2857143 -0.171429 -0.942857
r
-14 -21 14
0 -175 70
0 0 35
*/</syntaxhighlight>
=={{header|Scala}}==
{{Out}}Best seen running in your browser [https://scastie.scala-lang.org/NMueO16uQl6oivliBKZHew Scastie (remote JVM)].
<syntaxhighlight lang="scala">import java.io.{PrintWriter, StringWriter}
import Jama.{Matrix, QRDecomposition}
object QRDecomposition extends App {
val matrix =
new Matrix(
Array[Array[Double]](Array(12, -51, 4),
Array(6, 167, -68),
Array(-4, 24, -41)))
val d = new QRDecomposition(matrix)
def toString(m: Matrix): String = {
val sw = new StringWriter
m.print(new PrintWriter(sw, true), 8, 6)
sw.toString
}
print(toString(d.getQ))
print(toString(d.getR))
}</syntaxhighlight>
=={{header|SequenceL}}==
{{trans|Go}}
<syntaxhighlight lang="sequencel">import <Utilities/Math.sl>;
import <Utilities/Sequence.sl>;
import <Utilities/Conversion.sl>;
main :=
let
qrTest := [[12.0, -51.0, 4.0],
[ 6.0, 167.0, -68.0],
[-4.0, 24.0, -41.0]];
qrResult := qr(qrTest);
x := 1.0*(0 ... 10);
y := 1.0*[1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321];
regResult := polyfit(x, y, 2);
in
"q:\n" ++ delimit(delimit(floatToString(qrResult[1], 6), ','), '\n') ++ "\n\n" ++
"r:\n" ++ delimit(delimit(floatToString(qrResult[2], 1), ','), '\n') ++ "\n\n" ++
"polyfit:\n" ++ "[" ++ delimit(floatToString(regResult, 1), ',') ++ "]";
//---Polynomial Regression---
polyfit(x(1), y(1), n) :=
let
a[j] := x ^ j foreach j within 0 ... n;
in
lsqr(transpose(a), transpose([y]));
lsqr(a(2), b(2)) :=
let
qrDecomp := qr(a);
prod := mm(transpose(qrDecomp[1]), b);
in
solveUT(qrDecomp[2], prod);
solveUT(r(2), b(2)) :=
let
n := size(r[1]);
in
solveUTHelper(r, b, n, duplicate(0.0, n));
solveUTHelper(r(2), b(2), k, x(1)) :=
let
n := size(r[1]);
newX := setElementAt(x, k, (b[k][1] - sum(r[k][(k+1) ... n] * x[(k+1) ... n])) / r[k][k]);
in
x when k <= 0
else
solveUTHelper(r, b, k - 1, newX);
//---QR Decomposition---
qr(A(2)) := qrHelper(A, id(size(A)), 1);
qrHelper(A(2), Q(2), i) :=
let
m := size(A);
n := size(A[1]);
householder := makeHouseholder(A[i ... m, i]);
H[j,k] :=
householder[j - i + 1][k - i + 1] when j >= i and k >= i
else
1.0 when j = k else 0.0
foreach j within 1 ... m,
k within 1 ... m;
in
[Q,A] when i > (n - 1 when m = n else n)
else
qrHelper(mm(H, A), mm(Q, H), i + 1);
makeHouseholder(a(1)) :=
let
v := [1.0] ++ tail(a / (a[1] + sqrt(sum(a ^ 2)) * sign(a[1])));
H := id(size(a)) - (2.0 / mm([v], transpose([v])))[1,1] * mm(transpose([v]), [v]);
in
H;
//---Utilities---
id(n)[i,j] := 1.0 when i = j else 0.0
foreach i within 1 ... n,
j within 1 ... n;
mm(A(2), B(2))[i,j] := sum( A[i] * transpose(B)[j] );</syntaxhighlight>
{{out}}
<pre>
"q:
-0.857143,0.394286,0.331429
-0.428571,-0.902857,-0.034286
0.285714,-0.171429,0.942857
r:
-14.0,-21.0,14.0
-0.0,-175.0,70.0
0.0,0.0,-35.0
polyfit:
[1.0,2.0,3.0]"
</pre>
=={{header|SPAD}}==
See [[QR_decomposition#Axiom]] in Axiom.
=={{header|Standard ML}}==
{{trans|Axiom}}
We first define a signature for a radical category joined with a field. We then define a functor with (a) structures to define operators and functions for Array and Array2, and (b) functions for the QR decomposition:
<syntaxhighlight lang="sml">signature RADCATFIELD = sig
type real
val zero : real
val one : real
val + : real * real -> real
val - : real * real -> real
val * : real * real -> real
val / : real * real -> real
val sign : real -> real
val sqrt : real -> real
end
functor QR(F: RADCATFIELD) = struct
structure A = struct
local
open Array
in
fun unitVector n = tabulate (n, fn i => if i=0 then F.one else F.zero)
fun map f x = tabulate(length x, fn i => f(sub(x,i)))
fun map2 f (x, y) = tabulate(length x, fn i => f(sub(x,i),sub(y,i)))
val op + = map2 F.+
val op - = map2 F.-
val op * = map2 F.*
fun multc(c,x) = array(length x,c)*x
fun dot (x,y) = foldl F.+ F.zero (x*y)
fun outer f (x,y) =
Array2.tabulate Array2.RowMajor (length x, length y,
fn (i,j) => f(sub(x,i),sub(y,j)))
fun copy x = map (fn x => x) x
fun fromVector v = tabulate(Vector.length v, fn i => Vector.sub(v,i))
fun slice(x,i,sz) =
let open ArraySlice
val s = slice(x,i,sz)
in Array.tabulate(length s, fn i => sub(s,i)) end
end
end
structure M = struct
local
open Array2
in
fun map f x = tabulate RowMajor (nRows x, nCols x, fn (i,j) => f(sub(x,i,j)))
fun map2 f (x, y) =
tabulate RowMajor (nRows x, nCols x, fn (i,j) => f(sub(x,i,j),sub(y,i,j)))
fun scalarMatrix(m, x) = tabulate RowMajor (m,m,fn (i,j) => if i=j then x else F.zero)
fun multc(c, x) = map (fn xij => F.*(c,xij)) x
val op + = map2 F.+
val op - = map2 F.-
fun column(x,i) = A.fromVector(Array2.column(x,i))
fun row(x,i) = A.fromVector(Array2.row(x,i))
fun x*y = tabulate RowMajor (nRows x, nCols y,
fn (i,j) => A.dot(row(x,i), column(y,j)))
fun multa(x,a) = Array.tabulate (nRows x, fn i => A.dot(row(x,i), a))
fun copy x = map (fn x => x) x
fun subMatrix(h, i1, i2, j1, j2) =
tabulate RowMajor (Int.+(Int.-(i2,i1),1),
Int.+(Int.-(j2,j1),1),
fn (a,b) => sub(h,Int.+(i1,a),Int.+(j1,b)))
fun transpose m = tabulate RowMajor (nCols m,
nRows m,
fn (i,j) => sub(m,j,i))
fun updateSubMatrix(h,i,j,s) =
tabulate RowMajor (nRows s, nCols s, fn (a,b) => update(h,Int.+(i,a),Int.+(j,b),sub(s,a,b)))
end
end
fun toList a =
List.tabulate(Array2.nRows a, fn i => List.tabulate(Array2.nCols a, fn j => Array2.sub(a,i,j)))
fun householder a =
let open Array
val m = length a
val len = F.sqrt(A.dot(a,a))
val u = A.+(a, A.multc(F.*(len,F.sign(sub(a,0))), A.unitVector m))
val v = A.multc(F./(F.one,sub(u,0)), u)
val beta = F./(F.+(F.one,F.one),A.dot(v,v))
in
M.-(M.scalarMatrix(m,F.one), M.multc(beta,A.outer F.* (v,v)))
end
fun qr mat =
let open Array2
val (m,n) = dimensions mat
val upperIndex = if m=n then Int.-(n,1) else n
fun loop(i,qm,rm) = if i=upperIndex then {q=qm,r=rm} else
let val x = A.slice(A.fromVector(column(rm,i)),i,NONE)
val h = M.scalarMatrix(m,F.one)
val _ = M.updateSubMatrix(h,i,i,householder x)
in
loop(Int.+(i,1), M.*(qm,h), M.*(h,rm))
end
in
loop(0, M.scalarMatrix(m,F.one), mat)
end
fun solveUpperTriangular(r,b) =
let open Array
val n = Array2.nCols r
val x = array(n, F.zero)
fun loop k =
let val index = Int.min(Int.-(n,1),Int.+(k,1))
val _ = update(x,k,
F./(F.-(sub(b,k),
A.dot(A.slice(x,index,NONE),
A.slice(M.row(r,k),index,NONE))),
Array2.sub(r,k,k)))
in
if k=0 then x else loop(Int.-(k,1))
end
in
loop (Int.-(n,1))
end
fun lsqr(a,b) =
let val {q,r} = qr a
val n = Array2.nCols r
in
solveUpperTriangular(M.subMatrix(r, 0, Int.-(n,1), 0, Int.-(n,1)),
M.multa(M.transpose(q), b))
end
fun pow(x,1) = x
| pow(x,n) = F.*(x,pow(x,Int.-(n,1)))
fun polyfit(x,y,n) =
let open Array2
val a = tabulate RowMajor (Array.length x,
Int.+(n,1),
fn (i,j) => if j=0 then F.one else
pow(Array.sub(x,i),j))
in
lsqr(a,y)
end
end</syntaxhighlight>
We can then show the examples:<syntaxhighlight lang="sml">structure RealRadicalCategoryField : RADCATFIELD = struct
open Real
val one = 1.0
val zero = 0.0
val sign = real o Real.sign
val sqrt = Real.Math.sqrt
end
structure Q = QR(RealRadicalCategoryField);
let
val mat = Array2.fromList [[12.0, ~51.0, 4.0], [6.0, 167.0, ~68.0], [~4.0, 24.0, ~41.0]]
val {q,r} = Q.qr(mat)
in
{q=Q.toList q; r=Q.toList r}
end;
(* output *)
val it =
{q=[[~0.857142857143,0.394285714286,0.331428571429],
[~0.428571428571,~0.902857142857,~0.0342857142857],
[0.285714285714,~0.171428571429,0.942857142857]],
r=[[~14.0,~21.0,14.0],[5.97812397875E~18,~175.0,70.0],
[4.47505280695E~16,0.0,~35.0]]} : {q:real list list, r:real list list}
let open Array
val x = fromList [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
val y = fromList [1.0, 6.0, 17.0, 34.0, 57.0, 86.0, 121.0, 162.0, 209.0, 262.0, 321.0]
in
Q.polyfit(x, y, 2)
end;
(* output *)
val it = [|1.0,2.0,3.0|] : real array</syntaxhighlight>
=={{header|Stata}}==
See [http://www.stata.com/help.cgi?mf_qrd QR decomposition] in Stata help.
<syntaxhighlight lang="stata">mata
: qrd(a=(12,-51,4\6,167,-68\-4,24,-41),q=.,r=.)
: a
1 2 3
+-------------------+
1 | 12 -51 4 |
2 | 6 167 -68 |
3 | -4 24 -41 |
+-------------------+
: q
1 2 3
+----------------------------------------------+
1 | -.8571428571 .3942857143 .3314285714 |
2 | -.4285714286 -.9028571429 -.0342857143 |
3 | .2857142857 -.1714285714 .9428571429 |
+----------------------------------------------+
: r
1 2 3
+----------------------+
1 | -14 -21 14 |
2 | 0 -175 70 |
3 | 0 0 -35 |
+----------------------+</syntaxhighlight>
=={{header|Tcl}}==
Assuming the presence of the Tcl solutions to these tasks: [[Element-wise operations]], [[Matrix multiplication]], [[Matrix transposition]]
{{trans|Common Lisp}}
<
namespace path {::tcl::mathfunc ::tcl::mathop}
proc sign x {expr {$x == 0 ? 0 : $x < 0 ? -1 : 1}}
Line 1,140 ⟶ 5,248:
}
return [list $Q $A]
}</
Demonstrating:
<
puts "==Q=="
print_matrix [lindex $demo 0] "%f"
puts "==R=="
print_matrix [lindex $demo 1] "%.1f"</
Output:
<pre>
Line 1,157 ⟶ 5,265:
0.0 -175.0 70.0
0.0 0.0 -35.0
</pre>
=={{header|VBA}}==
{{trans|Phix}}<syntaxhighlight lang="vb">Option Base 1
Private Function vtranspose(v As Variant) As Variant
'-- transpose a vector of length m into an mx1 matrix,
'-- eg {1,2,3} -> {1;2;3}
vtranspose = WorksheetFunction.Transpose(v)
End Function
Private Function mat_col(a As Variant, col As Integer) As Variant
Dim res() As Double
ReDim res(UBound(a))
For i = col To UBound(a)
res(i) = a(i, col)
Next i
mat_col = res
End Function
Private Function mat_norm(a As Variant) As Double
mat_norm = Sqr(WorksheetFunction.SumProduct(a, a))
End Function
Private Function mat_ident(n As Integer) As Variant
mat_ident = WorksheetFunction.Munit(n)
End Function
Private Function sq_div(a As Variant, p As Double) As Variant
Dim res() As Variant
ReDim res(UBound(a))
For i = 1 To UBound(a)
res(i) = a(i) / p
Next i
sq_div = res
End Function
Private Function sq_mul(p As Double, a As Variant) As Variant
Dim res() As Variant
ReDim res(UBound(a), UBound(a, 2))
For i = 1 To UBound(a)
For j = 1 To UBound(a, 2)
res(i, j) = p * a(i, j)
Next j
Next i
sq_mul = res
End Function
Private Function sq_sub(x As Variant, y As Variant) As Variant
Dim res() As Variant
ReDim res(UBound(x), UBound(x, 2))
For i = 1 To UBound(x)
For j = 1 To UBound(x, 2)
res(i, j) = x(i, j) - y(i, j)
Next j
Next i
sq_sub = res
End Function
Private Function matrix_mul(x As Variant, y As Variant) As Variant
matrix_mul = WorksheetFunction.MMult(x, y)
End Function
Private Function QRHouseholder(ByVal a As Variant) As Variant
Dim columns As Integer: columns = UBound(a, 2)
Dim rows As Integer: rows = UBound(a)
Dim m As Integer: m = WorksheetFunction.Max(columns, rows)
Dim n As Integer: n = WorksheetFunction.Min(rows, columns)
I_ = mat_ident(m)
Q_ = I_
Dim q As Variant
Dim u As Variant, v As Variant, j As Integer
For j = 1 To WorksheetFunction.Min(m - 1, n)
u = mat_col(a, j)
u(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)
Next j
'-- Get the upper triangular matrix R.
Dim R() As Variant
ReDim R(m, n)
For i = 1 To m 'in Phix this is n
For j = 1 To n 'in Phix this is i to n. starting at 1 to fill zeroes
R(i, j) = a(i, j)
Next j
Next i
Dim res(2) As Variant
res(1) = Q_
res(2) = R
QRHouseholder = res
End Function
Private Sub pp(m As Variant)
For i = 1 To UBound(m)
For j = 1 To UBound(m, 2)
Debug.Print Format(m(i, j), "0.#####"),
Next j
Debug.Print
Next i
End Sub
Public Sub main()
a = [{12, -51, 4; 6, 167, -68; -4, 24, -41;-1,1,0;2,0,3}]
result = QRHouseholder(a)
q = result(1)
r_ = result(2)
Debug.Print "A"
pp a
Debug.Print "Q"
pp q
Debug.Print "R"
pp r_
Debug.Print "Q * R"
pp matrix_mul(q, r_)
End Sub</syntaxhighlight>{{out}}
<pre>A
12, -51, 4,
6, 167, -68,
-4, 24, -41,
-1, 1, 0,
2, 0, 3,
Q
0,84641 -0,39129 -0,34312 0,06641 -0,09126
0,42321 0,90409 0,02927 0,01752 -0,04856
-0,28214 0,17042 -0,93286 -0,02237 0,14365
-0,07053 0,01404 0,0011 0,99738 0,00728
0,14107 -0,01666 0,10577 0,00291 0,98419
R
14,17745 20,66663 -13,40157
0, 175,04254 -70,08031
0, 0, 35,20154
0, 0, 0,
0, 0, 0,
Q * R
12, -51, 4,
6, 167, -68,
-4, 24, -41,
-1, 1, 0,
2, 0, 3, </pre>
Least squares
<syntaxhighlight lang="vb">Public Sub least_squares()
x = [{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}]
y = [{1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}]
Dim a() As Double
ReDim a(UBound(x), 3)
For i = 1 To UBound(x)
For j = 1 To 3
a(i, j) = x(i) ^ (j - 1)
Next j
Next i
result = QRHouseholder(a)
q = result(1)
r_ = result(2)
t = WorksheetFunction.Transpose(q)
b = matrix_mul(t, vtranspose(y))
Dim z(3) As Double
For k = 3 To 1 Step -1
Dim s As Double: s = 0
If k < 3 Then
For j = k + 1 To 3
s = s + r_(k, j) * z(j)
Next j
End If
z(k) = (b(k, 1) - s) / r_(k, k)
Next k
Debug.Print "Least-squares solution:",
For i = 1 To 3
Debug.Print Format(z(i), "0.#####"),
Next i
End Sub</syntaxhighlight>{{out}}
<pre>Least-squares solution: 1, 2, 3, </pre>
=={{header|Wren}}==
{{trans|C}}
{{libheader|Wren-matrix}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="wren">import "./matrix" for Matrix
import "./fmt" for Fmt
var minor = Fn.new { |x, d|
var nr = x.numRows
var nc = x.numCols
var m = Matrix.new(nr, nc)
for (i in 0...d) m[i, i] = 1
for (i in d...nr) {
for (j in d...nc) m[i, j] = x[i, j]
}
return m
}
var vmadd = Fn.new { |a, b, s|
var n = a.count
var c = List.filled(n, 0)
for (i in 0...n) c[i] = a[i] + s * b[i]
return c
}
var vmul = Fn.new { |v|
var n = v.count
var x = Matrix.new(n, n)
for (i in 0...n) {
for (j in 0...n) x[i, j] = -2 * v[i] * v[j]
}
for (i in 0...n) x[i, i] = x[i, i] + 1
return x
}
var vnorm = Fn.new { |x|
var n = x.count
var sum = 0
for (i in 0...n) sum = sum + x[i] * x[i]
return sum.sqrt
}
var vdiv = Fn.new { |x, d|
var n = x.count
var y = List.filled(n, 0)
for (i in 0...n) y[i] = x[i] / d
return y
}
var mcol = Fn.new { |m, c|
var n = m.numRows
var v = List.filled(n, 0)
for (i in 0...n) v[i] = m[i, c]
return v
}
var householder = Fn.new { |m|
var nr = m.numRows
var nc = m.numCols
var q = List.filled(nr, null)
var z = m.copy()
var k = 0
while (k < nc && k < nr-1) {
var e = List.filled(nr, 0)
z = minor.call(z, k)
var x = mcol.call(z, k)
var a = vnorm.call(x)
if (m[k, k] > 0) a = -a
for (i in 0...nr) e[i] = (i == k) ? 1 : 0
e = vmadd.call(x, e, a)
e = vdiv.call(e, vnorm.call(e))
q[k] = vmul.call(e)
z = q[k] * z
k = k + 1
}
var Q = q[0]
var R = q[0] * m
var i = 1
while (i < nc && i < nr-1) {
Q = q[i] * Q
i = i + 1
}
R = Q * m
Q = Q.transpose
return [Q, R]
}
var inp = [
[12, -51, 4],
[ 6, 167, -68],
[-4, 24, -41],
[-1, 1, 0],
[ 2, 0, 3]
]
var x = Matrix.new(inp)
var res = householder.call(x)
var Q = res[0]
var R = res[1]
var m = Q * R
System.print("Q:")
Fmt.mprint(Q, 8, 3)
System.print("\nR:")
Fmt.mprint(R, 8, 3)
System.print("\nQ * R:")
Fmt.mprint(m, 8, 3)</syntaxhighlight>
{{out}}
<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|
</pre>
=={{header|zkl}}==
<syntaxhighlight lang="zkl">var [const] GSL=Import("zklGSL"); // libGSL (GNU Scientific Library)
A:=GSL.Matrix(3,3).set(12.0, -51.0, 4.0,
6.0, 167.0, -68.0,
4.0, 24.0, -41.0);
Q,R:=A.QRDecomp();
println("Q:\n",Q.format());
println("R:\n",R.format());
println("Q*R:\n",(Q*R).format());</syntaxhighlight>
{{out}}
<pre>
Q:
-0.86, 0.47, -0.22
-0.43, -0.88, -0.20
-0.29, -0.08, 0.95
R:
-14.00, -34.71, 37.43
0.00, -172.80, 65.07
0.00, 0.00, -26.19
Q*R:
12.00, -51.00, 4.00
6.00, 167.00, -68.00
4.00, 24.00, -41.00
</pre>
| |||