Jump to content

QR decomposition: Difference between revisions

QR decomposition (view source)

Revision as of 15:40, 6 May 2014

106 bytes added ,  10 years ago
Updated D entry
m (→‎{{header|J}}: missing word)
(Updated D entry)
Line 778:
 
template elementwiseMat(string op) {
T[][] elementwiseMat(T, U)(in T[][] A, in UT B) pure nothrow {
pure nothrow if (is(U == T) || is(U == T[][])) {▼
static if (is(U == T[][]))
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)
staticR[r][] if= mixin(is(U"row[] ==" T~ op ~ "B")) {;
return R[r][] = mixin("row[] " ~ op ~ "B");
} else {
assert(row.length == B[r].length);▼
R[r][] = mixin("row[] " ~ op ~ "B[r][]");▼
}▼
 
T[][] elementwiseMat(T, U)(in T[][] A, in U[][] B)
▲ pure nothrow if (is(UUnqual!T == T) || is(Unqual!U == T[][])) {
▲ 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 msummSum = elementwiseMat!q{ + },
msubmSub = elementwiseMat!q{ - },
pmulpMul = elementwiseMat!q{ * },
pdivpDiv = elementwiseMat!q{ / };
 
bool isRectangular(T)(in T[][] mat) pure nothrow {
Line 825 ⟶ 828:
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;
Line 835 ⟶ 838:
}
 
Unqual!T[][] makeUnitVector(T)(in size_t dim) pure nothrow {
auto result = new Unqual!T[][](dim, 1);
foreach (row; result)
row[] = 0;
Line 844 ⟶ 847:
 
/// 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;
Line 853 ⟶ 856:
}
 
Unqual!T[][] slice2D(T)(in T[][] A,
in size_t ma, in size_t mb,
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];
Line 869 ⟶ 872:
 
T[][] mcol(T)(in T[][] A, in size_t n) pure nothrow {
return slice2D(A, 0, rows(A).rows - 1, n, n);
}
 
Line 885 ⟶ 888:
 
T[][] makeHouseholder(T)(in T[][] a) {
immutable size_t m = a.rows(a);
immutable T s = sgn(a[0][0]).sgn;
immutable e = makeUnitVector!T(m);
immutable u = msummSum(a, pmulpMul(e, a.norm(a) * s));
immutable v = pdivpDiv(u, u[0][0]);
immutable beta = 2.0 / v.transpose.matMul(transpose(v), v)[0][0];
return msubmSub(matId!T(m), pmulpMul(v.matMul(v, .transpose(v)), beta));
}
 
Line 900 ⟶ 903:
 
// 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.
immutable B = slice2D(A, i, m - 1, i, n - 1);
 
// Take the first column of the current submatrix B.
Line 910 ⟶ 913:
// Create the Householder matrix for the column and embed it
// into an mxm identity.
immutable H = matEmbed(matId!T(m), makeHouseholder(x).makeHouseholder, i, i);
 
// The product of all H matrices from the right hand side is
// the orthogonal matrix Q.
Q = Q.matMul(Q, H);
 
// The product of all H matrices with A from the LHS is the
// upper triangular matrix R.
A = H.matMul(H, A);
}
 
Line 929 ⟶ 932:
/// Solve an upper triangular system by back substitution.
T[][] solveUpperTriangular(T)(in T[][] R, in T[][] b) pure nothrow {
immutable size_t n = R.cols(R);
auto x = new T[][](n, 1);
 
Line 944 ⟶ 947:
/// Solve a linear least squares problem by QR decomposition.
T[][] lsqr(T)(T[][] A, in T[][] b) pure nothrow {
const qr = QRdecomposition(A).QRdecomposition;
immutable size_t n = cols(qr.R).cols;
return solveUpperTriangular(
slice2D(qr.R, 0, n - 1, 0, n - 1),
slice2D(matMul(transpose(qr.Q), .transpose.matMul(b), 0, n - 1, 0, 0));
}
 
Unqual!T[][] polyFit(T)(in T[][] x, in T[][] y, in size_t n) pure nothrow {
immutable size_t m = x.cols;
pure nothrow {
immutableauto size_t mA = colsnew T[][](xm, n + 1);
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).transpose);
}
 
void main() {
// constimmutable (Q, R) = QRdecomposition([[12.0, -51, 4],
constimmutable qr = QRdecomposition([[12.0, -51, 4],
[ 6.0, 167, -68],
[-4.0, 24, -41]]);
immutable string form = "[%([%(%2.3f, %)]%|,\n %)]\n";
writefln(form, qr.Q);
writefln(form, qr.R);
Line 972 ⟶ 974:
immutable x = [[0.0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]];
immutable y = [[1.0, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]];
writeln(polyFit(x, y, 2)).writeln;
}</lang>
{{out}}
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/189681"
Cookies help us deliver our services. By using our services, you agree to our use of cookies.
More information