QR decomposition: Difference between revisions
Content added Content deleted
Line 203: | Line 203: | ||
-0.0000 175.0000 -70.0000 |
-0.0000 175.0000 -70.0000 |
||
-0.0000 0.0000 35.0000</pre> |
-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. |
|||
<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}}== |
=={{header|Axiom}}== |