QR decomposition: Difference between revisions
Content added Content deleted
m (→{{header|J}}: missing word) |
(Updated D entry) |
||
Line 778:
template elementwiseMat(string op) {
T[][] elementwiseMat(T
pure nothrow if (is(U == T) || 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)
assert(row.length == B[r].length);▼
R[r][] = mixin("row[] " ~ op ~ "B[r][]");▼
}▼
T[][] elementwiseMat(T, U)(in T[][] A, in U[][] B)
if (A.empty)
return null;
auto R = new typeof(return)(A.length, A[0].length);
foreach (immutable r, const row; A) {
return R;
}
}
alias
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,
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:
}
auto B = new
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,
}
Line 885 ⟶ 888:
T[][] makeHouseholder(T)(in T[][] a) {
immutable
immutable T s =
immutable e = makeUnitVector!T(m);
immutable u =
immutable v =
immutable beta = 2.0 / v.transpose.matMul(
return
}
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),
// 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 929 ⟶ 932:
/// Solve an upper triangular system by back substitution.
T[][] solveUpperTriangular(T)(in T[][] R, in T[][] b) pure nothrow {
immutable
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 =
immutable
return solveUpperTriangular(
slice2D(qr.R, 0, n - 1, 0, n - 1),
slice2D
}
immutable size_t m = x.cols;
foreach (immutable i, row; A)
foreach (immutable j, ref item; row)
item = x[0][i] ^^ j;
return lsqr(A,
}
void main() {
//
[ 6.0, 167, -68],
[-4.0, 24, -41]]);
immutable
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]];
}</lang>
{{out}}
|