Jump to content

QR decomposition: Difference between revisions

QR decomposition (view source)

Revision as of 23:41, 22 March 2012

77 bytes added ,  12 years ago
Updated D entry
(Updated D entry)
Line 639:
std.typecons, std.numeric, std.range, std.conv;
 
T[][] elementwiseMat(string op, T, U)(in T[][] A, in U B) pure
pure /*nothrow*/ if (is(U == T) || is(U == T[][])) {
static if (is(U == T[][]))
assert(A.length == B.length);
if (!A.lengthempty)
return null;
auto R = new typeof(return)(A.length, A[0].length);
Line 658:
}
 
T[][] msum(T)(in T[][] A, in T[][] B) pure /*nothrow*/ {
return elementwiseMat!(q{ + }, T, T[][])(A, B);
}
T[][] msub(T)(in T[][] A, in T[][] B) pure /*nothrow*/ {
return elementwiseMat!(q{ - }, T, T[][])(A, B);
}
T[][] pmul(T)(in T[][] A, in T x) pure /*nothrow*/ {
return elementwiseMat!(q{ * }, T, T)(A, x);
}
T[][] pdiv(T)(in T[][] A, in T x) pure /*nothrow*/ {
return elementwiseMat!(q{ / }, T, T)(A, x);
}
 
bool isRectangular(T)(in T[][] matrixmat) pure nothrow {
//return matrix.all!(r => r.length == mat[0].length)();
foreach (row; matrix)
ifforeach (row.length !=; matrix[0].lengthmat)
if (row.length != mat[0].length)
return false;
return true;
Line 689 ⟶ 690:
aux[k] = b[k][j];
foreach (i; 0 .. a.length)
result[i][j] = dotProduct(a[i], .dotProduct(aux);
}
return result;
}
 
string prettyPrint(T)(in T[][] A) {
return "[" ~ array(map!text(A)).join(",\n ") ~ "]";
}
 
Line 706 ⟶ 703:
}
 
T norm(T)(in T[][] arraym) /*pure nothrow*/ {
return sqrt(reduce!q{ a + b ^^ 2 }(cast(T)0, transversal(m, 0)).sqrt();
transversal(array, 0)));
}
 
Line 759 ⟶ 755:
// Main routines ---------------
 
T[][] makeHouseholder(T)(in T[][] a) {
constimmutable size_t m = rows(a);
constimmutable T s = sgn(a[0][0]);
T[][] e = makeUnitVector!T(m);
T[][] u = msum(a, pmul(e, norm(a) * s));
Line 770 ⟶ 766:
 
Tuple!(T[][],"Q", T[][],"R") QRdecomposition(T)(T[][] A) {
constimmutable m = rows(A);
constimmutable n = cols(A);
T[][] Q = matId!T(m);
 
Line 804 ⟶ 800:
/// Solve an upper triangular system by back substitution.
T[][] solveUpperTriangular(T)(in T[][] R, in T[][] b) pure nothrow {
constimmutable size_t n = cols(R);
auto x = new T[][](n, 1);
 
Line 818 ⟶ 814:
 
/// Solve a linear least squares problem by QR decomposition.
T[][] lsqr(T)(T[][] A, in T[][] b) {
const qr = QRdecomposition(A);
constimmutable size_t n = cols(qr.R);
return solveUpperTriangular(
slice2D(qr.R, 0, n-1, 0, n-1),
Line 836 ⟶ 832:
 
void main() {
// const (Q, R) = QRdecomposition(
const qr = QRdecomposition([[12.0, -51, 4],
[ 6.0, 167, -68],
[-4.0, 24, -41]]);
enum string form = "[%([%(%2.3f, %)]%|,\n %)]\n";
writeln(prettyPrint(qr.Q));
writeln(prettyPrintwritefln(form, qr.R), "\n"Q);
writefln(form, qr.R);
 
constimmutable x = [[0.0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]];
constimmutable y = [[1.0, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]];
writeln(polyFit(x, y, 2));
}</lang>
Output:<pre>[[-0.857143857, 0.394286394, 0.331429331],
[-0.428571429, -0.902857903, -0.0342857034],
[0.285714286, -0.171429171, 0.942857943]]
 
[[-14, -21, 14],
[2.26656e[-1614.000, -17521.000, 7014.000],
[40.72101e-16000, -7175.42462e-16000, -35]70.000],
[0.000, -0.000, -35.000]]
 
[[1], [2], [3]]</pre>
 
=={{header|Go}}==
{{trans|Common Lisp}}
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/189657"
Cookies help us deliver our services. By using our services, you agree to our use of cookies.
More information