Anonymous user
QR decomposition: Difference between revisions
D entry updated
(GP) |
(D entry updated) |
||
Line 784:
auto R = new typeof(return)(A.length, A[0].length);
foreach (immutable r, const row; A)
static if (is(U == T)) {
R[r][] = mixin("row[] " ~ op ~ "B");
Line 808:
}
bool isRectangular(T)(in T[][] mat) /*pure nothrow*/ {
}
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] = a[i].dotProduct(aux);
}
Line 834 ⟶ 830:
Unqual!T[][] transpose(T)(in T[][] m) pure nothrow {
auto r = new Unqual!T[][](m[0].length, m.length);
foreach (immutable nr, const row; m)
foreach (immutable nc, immutable c; row)
r[nc][nr] = c;
return r;
}
T norm(T)(in T[][] m)
return reduce!q{ a + b ^^ 2 }(cast(T)0, transversal(m, 0)).sqrt
}
Line 855 ⟶ 851:
T[][] matId(T)(in size_t n) pure nothrow {
auto Id = new T[][](n, n);
foreach (immutable r, row; Id) {
row[] = 0;
row[r] = 1;
Line 866 ⟶ 862:
in size_t na, in size_t nb) pure nothrow {
auto B = new Unqual!T[][](mb - ma + 1, nb - na + 1);
foreach (immutable i, brow; B)
brow[] = A[ma + i][na .. na + brow.length];
return B;
Line 872 ⟶ 868:
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 883 ⟶ 880:
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[]; // some wasted copies
foreach (immutable i, const brow; B)
C[row + i][col .. col + brow.length] = brow[];
return C;
Line 895 ⟶ 892:
immutable size_t m = rows(a);
immutable T s = sgn(a[0][0]);
return msub(matId!T(m), pmul(matMul(v, transpose(v)), beta));
}
Line 905 ⟶ 902:
immutable m = rows(A);
immutable n = cols(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
Line 940 ⟶ 937:
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 960 ⟶ 957:
Unqual!T[][] polyFit(T)(in T[][] x, in T[][] y, in size_t n) {
auto A = new Unqual!T[][](m, n + 1);
foreach (immutable i, row; A)
foreach (immutable j, ref item; row)
item = x[0][i] ^^ j;
return lsqr(A, transpose(y));
Line 973 ⟶ 970:
[ 6.0, 167, -68],
[-4.0, 24, -41]]);
writefln(form, qr.Q);
writefln(form, qr.R);
Line 981 ⟶ 978:
writeln(polyFit(x, y, 2));
}</lang>
{{out}}
[-0.429, -0.903, -0.034],
[0.286, -0.171, 0.943]]
| |||