diff --git a/dimred/ISOMAP.js b/dimred/ISOMAP.js index 8630fc6..76abdeb 100755 --- a/dimred/ISOMAP.js +++ b/dimred/ISOMAP.js @@ -69,9 +69,11 @@ export class ISOMAP extends DR { for (let i = 0; i < rows; ++i) { for (let j = 0; j < rows; ++j) { + let min_val = G.entry(i, j); for (let k = 0; k < rows; ++k) { - G.set_entry(i, j, Math.min(G.entry(i, j), G.entry(i, k) + G.entry(k, j))); + min_val = Math.min(min_val, G.entry(i, k) + G.entry(k, j)); } + G.set_entry(i, j, min_val); } } diff --git a/dimred/LDA.js b/dimred/LDA.js index cd04b44..4bb18bb 100755 --- a/dimred/LDA.js +++ b/dimred/LDA.js @@ -70,7 +70,7 @@ export class LDA extends DR { const v = V_mean.row(unique_labels[label].id); const m = new Matrix(cols, 1, (j) => v[j] - X_mean); const N = unique_labels[label].count; - S_b = S_b.add(m.dot(m.transpose()).mult(N)); + S_b = S_b.add(m.dotTrans(m).mult(N)); } // scatter_within @@ -81,7 +81,7 @@ export class LDA extends DR { const R = unique_labels[label].rows; for (let i = 0, n = unique_labels[label].count; i < n; ++i) { const row_v = new Matrix(cols, 1, (j, _) => R[i][j] - m.entry(j, 0)); - S_w = S_w.add(row_v.dot(row_v.transpose())); + S_w = S_w.add(row_v.dotTrans(row_v)); } } diff --git a/dimred/LLE.js b/dimred/LLE.js index 06a9654..de322c3 100755 --- a/dimred/LLE.js +++ b/dimred/LLE.js @@ -49,11 +49,11 @@ export class LLE extends DR { for (let row = 0; row < rows; ++row) { const nN_row = nN[row]; const Z = new Matrix(neighbors, cols, (i, j) => X.entry(nN_row[i].j, j) - X.entry(row, j)); - const C = Z.dot(Z.T); + const C = Z.dotTrans(Z); if (neighbors > cols) { const C_trace = neumair_sum(C.diag) / 1000; for (let j = 0; j < neighbors; ++j) { - C.set_entry(j, j, C.entry(j, j) + C_trace); + C.add_entry(j, j, C_trace); } } // reconstruct; @@ -66,7 +66,7 @@ export class LLE extends DR { // comp embedding const I = new Matrix(rows, rows, "identity"); const IW = I.sub(W); - const M = IW.T.dot(IW); + const M = IW.transDot(IW); const { eigenvectors: V } = simultaneous_poweriteration(M.T.inverse(), d + 1, eig_args); this.Y = Matrix.from(V.slice(1, 1 + d)).T; diff --git a/dimred/LSP.js b/dimred/LSP.js index eb78376..a41fe39 100755 --- a/dimred/LSP.js +++ b/dimred/LSP.js @@ -85,10 +85,9 @@ export class LSP extends DR { transform() { this.check_init(); const A = this._A; - const AT = A.T; const b = this._b; - const ATA = AT.dot(A); - const ATb = AT.dot(b); + const ATA = A.transDot(A); + const ATb = A.transDot(b); this.Y = Matrix.solve_CG(ATA, ATb, this._randomizer); return this.projection; } diff --git a/dimred/LTSA.js b/dimred/LTSA.js index 67b0fb9..12219d2 100755 --- a/dimred/LTSA.js +++ b/dimred/LTSA.js @@ -55,17 +55,16 @@ export class LTSA extends DR { // center X_i X_i = X_i.dot(O); // correlation matrix - const C = X_i.dot(X_i.transpose()); + const C = X_i.dotTrans(X_i); const { eigenvectors: g } = simultaneous_poweriteration(C, d, eig_args); //g.push(linspace(0, k).map(_ => 1 / Math.sqrt(k + 1))); const G_i_t = Matrix.from(g); // 2. Constructing alignment matrix - const W_i = G_i_t.transpose() - .dot(G_i_t) + const W_i = G_i_t.transDot(G_i_t) .add(1 / Math.sqrt(neighbors + 1)); for (let i = 0; i < neighbors + 1; ++i) { for (let j = 0; j < neighbors + 1; ++j) { - B.set_entry(I_i[i], I_i[j], B.entry(I_i[i], I_i[j]) - (i === j ? 1 : 0) + W_i.entry(i, j)); + B.add_entry(I_i[i], I_i[j], W_i.entry(i, j) - (i === j ? 1 : 0)); } } } diff --git a/dimred/MLLE.js b/dimred/MLLE.js index 09f798c..1862125 100755 --- a/dimred/MLLE.js +++ b/dimred/MLLE.js @@ -34,7 +34,7 @@ export class MLLE{ let X_i = Matrix.from(I_i.map(n => X.row(n))); X_i = X_i.sub(x_i); //X_i = X_i.dot(new Matrix(X_i._cols, X_i._cols, "center")) - let C_i = X_i.dot(X_i.transpose()) // k by k + let C_i = X_i.dotTrans(X_i); // k by k let gamma = neumair_sum(C_i.diag) / 1000; for (let j = 0; j < k; ++j) { @@ -91,9 +91,9 @@ export class MLLE{ H_i = H_i.sub(h.mult(2).outer(h)); let W_i = V_i.sub(V_i.dot(h).dot(h.T).mult(2)).add(w_i.mult(1 - alpha_i)) */ - let W_i = V_i.sub(V_i.dot(h).dot(h.T).mult(2)).add(w_i.mult(1 - alpha_i).dot(ones.T)) - - W_i = W_i.dot(W_i.T) + let W_i = V_i.sub(V_i.dot(h).dotTrans(h).mult(2)).add(w_i.mult(1 - alpha_i).dotTrans(ones)) + + W_i = W_i.dotTrans(W_i); for (let i = 0; i < k + 1; ++i) { for (let j = 0; j < s_i; ++j) { Phi.set_entry(I_i[i], I_i[j], Phi.entry(I_i[i], I_i[j]) - (i === j ? 1 : 0 ) + W_i.entry(i, j)); diff --git a/dimred/OAP.js b/dimred/OAP.js index bd28701..bae12f8 100755 --- a/dimred/OAP.js +++ b/dimred/OAP.js @@ -220,7 +220,7 @@ export class OAP { return i < j ? -d_X.entry(i, j) / d_Y.entry(i, j) : ratio.entry(j, i); }] for (let i = 0; i < N; ++i) { - ratio.set_entry(i, i, ratio.entry(i, i) - neumair_sum(ratio.row(i))); + ratio.sub_entry(i, i, neumair_sum(ratio.row(i))); } const mds_Y = ratio.dot(Y).divide(N); @@ -232,7 +232,7 @@ export class OAP { const dm = M(mds_Y.row(i)); const h_i = h[i]; for (let d = 0; d < dim; ++d) { - Y.set_entry(i, d, Y.entry(i, d) + step_size * (diff_Y.entry(i, d) + w * 2 * (m - h_i) * dm)); + Y.add_entry(i, d, step_size * (diff_Y.entry(i, d) + w * 2 * (m - h_i) * dm)); } } diff --git a/dimred/PCA.js b/dimred/PCA.js index 9d314d1..a5071d4 100755 --- a/dimred/PCA.js +++ b/dimred/PCA.js @@ -57,9 +57,8 @@ export class PCA extends DR { } const { d, eig_args } = this._parameters; const X = this.X; - const means = Matrix.from(X.meanCols); - const X_cent = X.sub(means); - const C = X_cent.transpose().dot(X_cent); + const X_cent = X.sub(X.meanCols); + const C = X_cent.transDot(X_cent); const { eigenvectors: V } = simultaneous_poweriteration(C, d, eig_args); this.V = Matrix.from(V).transpose(); return this.V; diff --git a/dimred/TSNE.js b/dimred/TSNE.js index 210d0b5..3352043 100755 --- a/dimred/TSNE.js +++ b/dimred/TSNE.js @@ -1,5 +1,5 @@ import { Matrix } from "../matrix/index.js"; -import { euclidean } from "../metrics/index.js"; +import { euclidean_squared } from "../metrics/index.js"; import { DR } from "./DR.js"; /** @@ -18,15 +18,15 @@ export class TSNE extends DR { * @param {Number} [parameters.perplexity = 50] - perplexity. * @param {Number} [parameters.epsilon = 10] - learning parameter. * @param {Number} [parameters.d = 2] - the dimensionality of the projection. - * @param {Function|"precomputed"} [parameters.metric = euclidean] - the metric which defines the distance between two points. + * @param {Function|"precomputed"} [parameters.metric = euclidean_squared] - the metric which defines the distance between two points. * @param {Number} [parameters.seed = 1212] - the seed for the random number generator. * @returns {TSNE} */ constructor(X, parameters) { - super(X, { perplexity: 50, epsilon: 10, d: 2, metric: euclidean, seed: 1212 }, parameters); + super(X, { perplexity: 50, epsilon: 10, d: 2, metric: euclidean_squared, seed: 1212 }, parameters); [this._N, this._D] = this.X.shape; this._iter = 0; - this.Y = new Matrix(this._N, this.parameter("d"), () => this._randomizer.random); + this.Y = new Matrix(this._N, this.parameter("d"), () => this._randomizer.gauss_random() * 1e-4); return this; } @@ -56,66 +56,62 @@ export class TSNE extends DR { } } - const P = new Matrix(N, N, "zeros"); + const P = new Matrix(N, N, 0); - this._ystep = new Matrix(N, D, "zeros"); + this._ystep = new Matrix(N, D, 0); this._gains = new Matrix(N, D, 1); // search for fitting sigma - let prow = new Float64Array(N) const tol = 1e-4; const maxtries = 50; for (let i = 0; i < N; ++i) { + const dist_i = Delta.row(i); + const prow = P.row(i); let betamin = -Infinity; let betamax = Infinity; let beta = 1; + let cnt = maxtries; let done = false; + let psum; - let num = 0; - while (!done) { - let psum = 0; + while (!done && cnt--) { + // compute entropy and kernel row with beta precision + psum = 0; + let dp_sum = 0; for (let j = 0; j < N; ++j) { - let pj = Math.exp(-Delta.entry(i, j) * beta); - if (i === j) pj = 0; + const dist = dist_i[j]; + const pj = (i !== j) ? Math.exp(-dist * beta) : 0; + dp_sum += dist * pj; prow[j] = pj; psum += pj; } - let Hhere = 0; - for (let j = 0; j < N; ++j) { - let pj = psum === 0 ? 0 : prow[j] / psum; - prow[j] = pj; - if (pj > 1e-7) { - Hhere -= pj * Math.log(pj); - } - } - if (Hhere > Htarget) { + // compute entropy + const H = psum > 0 ? Math.log(psum) + beta * dp_sum / psum : 0; + if (H > Htarget) { betamin = beta; beta = betamax === Infinity ? beta * 2 : (beta + betamax) / 2; } else { betamax = beta; beta = betamin === -Infinity ? beta / 2 : (beta + betamin) / 2; } - ++num; - if (Math.abs(Hhere - Htarget) < tol) done = true; - if (num >= maxtries) done = true; + done = Math.abs(H - Htarget) < tol; } - + // normalize p for (let j = 0; j < N; ++j) { - P.set_entry(i, j, prow[j]); + prow[j] /= psum; } } - //compute probabilities - const Pout = new Matrix(N, N, "zeros"); + // compute probabilities const N2 = N * 2; for (let i = 0; i < N; ++i) { for (let j = i; j < N; ++j) { const p = Math.max((P.entry(i, j) + P.entry(j, i)) / N2, 1e-100); - Pout.set_entry(i, j, p); - Pout.set_entry(j, i, p); + P.set_entry(i, j, p); + P.set_entry(j, i, p); } } - this._P = Pout; + this._P = P; return this; } @@ -195,7 +191,7 @@ export class TSNE extends DR { for (let j = 0; j < N; ++j) { const premult = 4 * (pmul * P.entry(i, j) - Q.entry(i, j)) * Qu.entry(i, j); for (let d = 0; d < dim; ++d) { - grad.set_entry(i, d, grad.entry(i, d) + premult * (Y.entry(i, d) - Y.entry(j, d))); + grad.add_entry(i, d, premult * (Y.entry(i, d) - Y.entry(j, d))); } } } @@ -216,14 +212,14 @@ export class TSNE extends DR { const newsid = momval * sid - epsilon * newgain * gid; ystep.set_entry(i, d, newsid); - Y.set_entry(i, d, Y.entry(i, d) + newsid); + Y.add_entry(i, d, newsid); ymean[d] += Y.entry(i, d); } } for (let i = 0; i < N; ++i) { - for (let d = 0; d < 2; ++d) { - Y.set_entry(i, d, Y.entry(i, d) - ymean[d] / N); + for (let d = 0; d < dim; ++d) { + Y.sub_entry(i, d, ymean[d] / N); } } diff --git a/dimred/TriMap.js b/dimred/TriMap.js index 65ae919..c4003c5 100755 --- a/dimred/TriMap.js +++ b/dimred/TriMap.js @@ -40,11 +40,11 @@ export class TriMap extends DR { init(pca = null, knn = null) { const X = this.X; const N = X.shape[0]; - const { d, metric, c } = this._parameters; + const { c, d, metric, seed } = this._parameters; this.n_inliers = 2 * c; this.n_outliers = 1 * c; this.n_random = 1 * c; - this.Y = pca || new PCA(X, d).transform(); + this.Y = pca || new PCA(X, { d, seed }).transform(); this.knn = knn || new BallTree(X.to2dArray, metric); const { triplets, weights } = this._generate_triplets(this.n_inliers, this.n_outliers, this.n_random); this.triplets = triplets; @@ -156,7 +156,7 @@ export class TriMap extends DR { const triplets = new Matrix(N * n_inliers * n_outliers, 3); for (let i = 0; i < N; ++i) { let n_i = i * n_inliers * n_outliers; - const sort_indices = this.__argsort(P.row(i).map((d) => -d)); + const sort_indices = this.__argsort(P.row(i)); for (let j = 0; j < n_inliers; ++j) { let n_j = j * n_outliers; const sim = nbrs.entry(i, sort_indices[j]); @@ -179,11 +179,7 @@ export class TriMap extends DR { * @param {Array} A */ __argsort(A) { - return A.map((d, i) => { - return { d: d, i: i }; - }) - .sort((a, b) => a.d - b.d) - .map((d) => d.i); + return linspace(0, A.length - 1).sort((i, j) => A[j] - A[i]); } /** @@ -314,9 +310,9 @@ export class TriMap extends DR { for (let d = 0; d < dim; ++d) { const gs = y_ij[d] * d_ik * w; const go = y_ik[d] * d_ij * w; - grad.set_entry(i, d, grad.entry(i, d) + gs - go); - grad.set_entry(j, d, grad.entry(j, d) - gs); - grad.set_entry(k, d, grad.entry(k, d) + go); + grad.add_entry(i, d, gs - go); + grad.sub_entry(j, d, gs); + grad.add_entry(k, d, go); } } return { grad, loss, n_viol }; diff --git a/dimred/UMAP.js b/dimred/UMAP.js index d808f99..ed7c497 100755 --- a/dimred/UMAP.js +++ b/dimred/UMAP.js @@ -86,9 +86,11 @@ export class UMAP extends DR { */ _compute_membership_strengths(distances, sigmas, rhos) { for (let i = 0, n = distances.length; i < n; ++i) { - for (let j = 0, m = distances[i].length; j < m; ++j) { - const v = distances[i][j].value - rhos[i]; - distances[i][j].value = v > 0 ? Math.exp(-v / sigmas[i]) : 1; + const rho = rhos[i]; + const curr_dist = distances[i]; + for (let j = 0, m = curr_dist.length; j < m; ++j) { + const v = curr_dist[j].value - rho; + curr_dist[j].value = v > 0 ? Math.exp(-v / sigmas[i]) : 1.0; } } return distances; @@ -123,32 +125,33 @@ export class UMAP extends DR { } } + const index = Math.floor(local_connectivity); + const interpolation = local_connectivity - index; for (let i = 0; i < N; ++i) { let lo = 0; let hi = Infinity; let mid = 1; + let rho = 0; const search_result = distances[i]; const non_zero_dist = search_result.filter((d) => d.value > 0); const non_zero_dist_length = non_zero_dist.length; if (non_zero_dist_length >= local_connectivity) { - const index = Math.floor(local_connectivity); - const interpolation = local_connectivity - index; if (index > 0) { - rhos.push(non_zero_dist[index - 1]); + rho = non_zero_dist[index - 1].value; if (interpolation > SMOOTH_K_TOLERANCE) { - rhos[i].value += interpolation * (non_zero_dist[index].value - non_zero_dist[index - 1]); + rho += interpolation * (non_zero_dist[index].value - non_zero_dist[index - 1].value); } } else { - rhos[i].value = interpolation * non_zero_dist[0].value; + rho = interpolation * non_zero_dist[0].value; } } else if (non_zero_dist_length > 0) { - rhos[i] = non_zero_dist[non_zero_dist_length - 1].value; + rho = non_zero_dist[non_zero_dist_length - 1].value; } for (let x = 0; x < n_iter; ++x) { let psum = 0; for (let j = 0; j < k; ++j) { - const d = search_result[j].value - rhos[i]; + const d = search_result[j].value - rho; psum += d > 0 ? Math.exp(-(d / mid)) : 1; } if (Math.abs(psum - target) < SMOOTH_K_TOLERANCE) { @@ -164,20 +167,21 @@ export class UMAP extends DR { } } } - sigmas[i] = mid; - const mean_ithd = search_result.reduce((a, b) => a + b.value, 0) / search_result.length; //let mean_d = null; - if (rhos[i] > 0) { - if (sigmas[i] < MIN_K_DIST_SCALE * mean_ithd) { - sigmas[i] = MIN_K_DIST_SCALE * mean_ithd; + if (rho > 0) { + const mean_ithd = search_result.reduce((a, b) => a + b.value, 0) / search_result.length; + if (mid < MIN_K_DIST_SCALE * mean_ithd) { + mid = MIN_K_DIST_SCALE * mean_ithd; } } else { const mean_d = distances.reduce((acc, res) => acc + res.reduce((a, b) => a + b.value, 0) / res.length); - if (sigmas[i] > MIN_K_DIST_SCALE * mean_d) { - sigmas[i] = MIN_K_DIST_SCALE * mean_d; + if (mid < MIN_K_DIST_SCALE * mean_d) { + mid = MIN_K_DIST_SCALE * mean_d; } } + rhos[i] = rho; + sigmas[i] = mid; } return { distances: distances, @@ -223,9 +227,11 @@ export class UMAP extends DR { _make_epochs_per_sample(n_epochs) { const weights = this._weights; const result = new Float32Array(weights.length).fill(-1); - const weights_max = max(weights); - const n_samples = weights.map((w) => n_epochs * (w / weights_max)); - for (let i = 0; i < result.length; ++i) if (n_samples[i] > 0) result[i] = Math.round(n_epochs / n_samples[i]); + const weight_scl = n_epochs / max(weights); + weights.forEach((w, i) => { + const sample = w * weight_scl; + if (sample > 0) result[i] = Math.round(n_epochs / sample); + }) return result; } @@ -350,18 +356,13 @@ export class UMAP extends DR { const current = head_embedding.row(j); const other = tail_embedding.row(k); const dist = euclidean_squared(current, other); - let grad_coeff = 0; if (dist > 0) { - grad_coeff = (-2 * a * b * Math.pow(dist, b - 1)) / (a * Math.pow(dist, b) + 1); - } - for (let d = 0; d < dim; ++d) { - const grad_d = clip(grad_coeff * (current[d] - other[d])) * alpha; - const c = current[d] + grad_d; - const o = other[d] - grad_d; - current[d] = c; - other[d] = o; - head_embedding.set_entry(j, d, c); - tail_embedding.set_entry(k, d, o); + const grad_coeff = (-2 * a * b * Math.pow(dist, b - 1)) / (a * Math.pow(dist, b) + 1); + for (let d = 0; d < dim; ++d) { + const grad_d = clip(grad_coeff * (current[d] - other[d])) * alpha; + current[d] += grad_d; + other[d] -= grad_d; + } } epoch_of_next_sample[i] += epochs_per_sample[i]; const n_neg_samples = (this._iter - epoch_of_next_negative_sample[i]) / epochs_per_negative_sample[i]; @@ -369,21 +370,16 @@ export class UMAP extends DR { const k = randomizer.random_int % tail_length; const other = tail_embedding.row(tail[k]); const dist = euclidean_squared(current, other); - let grad_coeff = 0; if (dist > 0) { - grad_coeff = (2 * _repulsion_strength * b) / ((0.01 + dist) * (a * Math.pow(dist, b) + 1)); + const grad_coeff = (2 * _repulsion_strength * b) / ((0.01 + dist) * (a * Math.pow(dist, b) + 1)); + for (let d = 0; d < dim; ++d) { + const grad_d = clip(grad_coeff * (current[d] - other[d])) * alpha; + current[d] += grad_d; + other[d] -= grad_d; + } } else if (j === k) { continue; } - for (let d = 0; d < dim; ++d) { - const grad_d = clip(grad_coeff * (current[d] - other[d])) * alpha; - const c = current[d] + grad_d; - const o = other[d] - grad_d; - current[d] = c; - other[d] = o; - head_embedding.set_entry(j, d, c); - tail_embedding.set_entry(tail[k], d, o); - } } epoch_of_next_negative_sample[i] += n_neg_samples * epochs_per_negative_sample[i]; } diff --git a/linear_algebra/poweriteration_m.js b/linear_algebra/poweriteration_m.js index 0e76c00..81bc97e 100755 --- a/linear_algebra/poweriteration_m.js +++ b/linear_algebra/poweriteration_m.js @@ -12,7 +12,7 @@ import { Randomizer } from "../util/index.js"; export default function(data, x0, beta, max_iter=20, seed) { let randomizer = new Randomizer(seed); let [ n, d ] = data.shape; - let A = data.T.dot(data).divide(n) + let A = data.transDot(data).divide(n) if (x0 === null) x0 = new Matrix(d, 1, () => randomizer.random) x0 = x0.divide(norm(x0)); let x = x0.clone() diff --git a/linear_algebra/qr.js b/linear_algebra/qr.js index 9e2bc4c..c3c5d7e 100755 --- a/linear_algebra/qr.js +++ b/linear_algebra/qr.js @@ -20,8 +20,10 @@ export default function (A) { for (let i = 0; i < j; ++i) { const q = Q.col(i); const q_dot_v = neumair_sum(q.map((q_, k) => q_ * v[k])); + for (let k = 0; k < rows; ++k) { + v[k] -= q_dot_v * q[k]; + } R.set_entry(i, j, q_dot_v); - v = v.map((v_, k) => v_ - q_dot_v * q[k]); } const v_norm = norm(v, euclidean); for (let k = 0; k < rows; ++k) { diff --git a/linear_algebra/qr_givens.js b/linear_algebra/qr_givens.js index 711b3fe..696c152 100755 --- a/linear_algebra/qr_givens.js +++ b/linear_algebra/qr_givens.js @@ -22,7 +22,7 @@ export default function (A) { Gji.set_entry(i - 1, i, -s); Gji.set_entry(i, i - 1, s); Gji.set_entry(i, i, c); - R = Gji.T.dot(R); + R = Gji.transDot(R); Q = Q.dot(Gji); } } diff --git a/linear_algebra/sapi.js b/linear_algebra/sapi.js index 2b377d1..54bf7f2 100755 --- a/linear_algebra/sapi.js +++ b/linear_algebra/sapi.js @@ -26,8 +26,8 @@ export default function (A, max_iterations = 100, batch_size = 10, beta = 0.05, r = r_next.divide(r_next_norm); } - let u = r.transpose().dot(A).dot(r); - let l = r.transpose().dot(r); + let u = r.transDot(A).dot(r); + let l = r.transDot(r); let lambda = u.divide(l).entry(0, 0); return { eigenvector: r.transpose().to2dArray[0], diff --git a/linear_algebra/simultaneous_poweriteration.js b/linear_algebra/simultaneous_poweriteration.js index 671212b..87e49eb 100755 --- a/linear_algebra/simultaneous_poweriteration.js +++ b/linear_algebra/simultaneous_poweriteration.js @@ -22,7 +22,7 @@ export default function (A, k = 2, {seed = 1212, max_iterations = 100, qr = qr_g const n = A.shape[0]; let { Q, R } = qr(new Matrix(n, k, () => (randomizer.random - .5) * 2)); while (max_iterations--) { - const oldQ = Q.clone(); + const oldQ = Q; const Z = A.dot(Q); const QR = qr(Z); Q = QR.Q; diff --git a/linear_algebra/svrg.js b/linear_algebra/svrg.js index d025439..5dd0e62 100755 --- a/linear_algebra/svrg.js +++ b/linear_algebra/svrg.js @@ -20,15 +20,15 @@ export default function(data, x, beta, epoch=20, m=10, s=1, seed) { x = new Matrix(d, 1, () => randomizer.random) x = x.divide(norm(x)); let x0 = x.clone(); - let A = data.T.dot(data).divide(n); + let A = data.transDot(data).divide(n); let x_tilde = x.clone(); for (let t = 0; t < epoch; ++t) { const gx = A.dot(x_tilde); for (let i = 0; i < m; ++i) { - const ang = x.T.dot(x_tilde).entry(0, 0); + const ang = x.transDot(x_tilde).entry(0, 0); const sample = Matrix.from(Randomizer.choice(data, s)); - const sampleT_dot_sample = sample.T.dot(sample) + const sampleT_dot_sample = sample.transDot(sample) const x_tmp = x.clone(); const XTXx = sampleT_dot_sample .dot(x.divide(s)); diff --git a/matrix/Matrix.js b/matrix/Matrix.js index bbfc1f0..7b72b4d 100755 --- a/matrix/Matrix.js +++ b/matrix/Matrix.js @@ -91,11 +91,11 @@ export class Matrix { static from(A, type = "row") { if (A instanceof Matrix) { return A.clone(); - } else if (Array.isArray(A) || A instanceof Float64Array) { + } else if (Matrix.isArray(A)) { let m = A.length; if (m === 0) throw new Error("Array is empty"); // 1d - if (!Array.isArray(A[0]) && !(A[0] instanceof Float64Array)) { + if (!Matrix.isArray(A[0])) { if (type === "row") { return new Matrix(1, m, (_, j) => A[j]); } else if (type === "col") { @@ -106,7 +106,7 @@ export class Matrix { throw new Error("1d array has NaN entries"); } // 2d - } else if (Array.isArray(A[0]) || A[0] instanceof Float64Array) { + } else { let n = A[0].length; for (let row = 0; row < m; ++row) { if (A[row].length !== n) { @@ -164,7 +164,7 @@ export class Matrix { */ set_row(row, values) { const cols = this._cols; - if ((Array.isArray(values) || values instanceof Float64Array) && values.length === cols) { + if (Matrix.isArray(values) && values.length === cols) { const offset = row * cols; for (let col = 0; col < cols; ++col) { this.values[offset + col] = values[col]; @@ -180,6 +180,22 @@ export class Matrix { return this; } + /** + * Swaps the rows {@link row1} and {@link row2} of the Matrix. + * @param {Number} row1 + * @param {Number} row2 + * @returns {Matrix} + */ + swap_rows(row1, row2) { + const cols = this._cols; + const data = this.values; + for (let i = row1 * cols, j = row2 * cols, col = 0; col < cols; ++col, ++i, ++j) { + const t = data[i]; + data[i] = data[j]; + data[j] = t; + } + } + /** * Returns the {@link col}^th column from the Matrix. * @param {Number} col @@ -215,6 +231,30 @@ export class Matrix { return this; } + /** + * Adds a given {@link value} to the {@link col}^th entry from the {@link row}^th row of the Matrix. + * @param {int} row + * @param {int} col + * @param {float64} value + * @returns {Matrix} + */ + add_entry(row, col, value) { + this.values[row * this._cols + col] += value; + return this; + } + + /** + * Subtracts a given {@link value} from the {@link col}^th entry from the {@link row}^th row of the Matrix. + * @param {int} row + * @param {int} col + * @param {float64} value + * @returns {Matrix} + */ + sub_entry(row, col, value) { + this.values[row * this._cols + col] -= value; + return this; + } + /** * Returns a new transposed Matrix. * @returns {Matrix} @@ -239,96 +279,174 @@ export class Matrix { inverse() { const rows = this._rows; const cols = this._cols; - let B = new Matrix(rows, 2 * cols, (i, j) => { - if (j >= cols) { - return i === j - cols ? 1 : 0; - } else { - return this.entry(i, j); - } - }); - let h = 0; - let k = 0; - while (h < rows && k < cols) { - var i_max = 0; - let max_val = -Infinity; - for (let i = h; i < rows; ++i) { - let val = Math.abs(B.entry(i, k)); + const A = this.clone(); + const B = new Matrix(rows, cols, 'I'); + + // foreach column + for (let col = 0; col < cols; ++col) { + // Search for maximum in this column (pivot) + let max_idx = col; + let max_val = Math.abs(A.entry(col, col)); + for (let row = col + 1; row < rows; ++row) { + const val = Math.abs(A.entry(row, col)); if (max_val < val) { - i_max = i; + max_idx = row; max_val = val; } } - if (B.entry(i_max, k) == 0) { - k++; - } else { - // swap rows - for (let j = 0; j < 2 * cols; ++j) { - let h_val = B.entry(h, j); - let i_val = B.entry(i_max, j); - B.set_entry(h, j, h_val); - B.set_entry(i_max, j, i_val); - } - for (let i = h + 1; i < rows; ++i) { - let f = B.entry(i, k) / B.entry(h, k); - B.set_entry(i, k, 0); - for (let j = k + 1; j < 2 * cols; ++j) { - B.set_entry(i, j, B.entry(i, j) - B.entry(h, j) * f); + if (max_val === 0) { + throw new Error('Cannot compute inverse of Matrix, determinant is zero'); + } + // Swap maximum row with current row + if (max_idx !== col) { + A.swap_rows(col, max_idx); + B.swap_rows(col, max_idx); + } + + // eliminate non-zero values on the other rows at column c + const A_col = A.row(col); + const B_col = B.row(col); + for (let row = 0; row < rows; ++row) { + if (row !== col) { + // eliminate value at column c and row r + const A_row = A.row(row); + const B_row = B.row(row); + if (A_row[col] !== 0) { + const f = A_row[col] / A_col[col]; + // sub (f * row c) from row r to eliminate the value at column c + for (let s = col; s < cols; ++s) { + A_row[s] -= (f * A_col[s]); + } + for (let s = 0; s < cols; ++s) { + B_row[s] -= (f * B_col[s]); + } + } + } else { + // normalize value at Acc to 1 (diagonal): + // divide each value of row r=c by the value at Acc + const f = A_col[col]; + for (let s = col; s < cols; ++s) { + A_col[s] /= f; + } + for (let s = 0; s < cols; ++s) { + B_col[s] /= f; } } - h++; - k++; } } + return B; + } - for (let row = 0; row < rows; ++row) { - let f = B.entry(row, row); - for (let col = row; col < 2 * cols; ++col) { - B.set_entry(row, col, B.entry(row, col) / f); + /** + * Returns the dot product. If {@link B} is an Array or Float64Array then an Array gets returned. If {@link B} is a Matrix then a Matrix gets returned. + * @param {(Matrix|Array|Float64Array)} B the right side + * @returns {(Matrix|Array)} + */ + dot(B) { + if (B instanceof Matrix) { + let A = this; + const [rows_A, cols_A] = A.shape; + const [rows_B, cols_B] = B.shape; + if (cols_A !== rows_B) { + throw new Error(`A.dot(B): A is a ${A.shape.join(" ⨯ ")}-Matrix, B is a ${B.shape.join(" ⨯ ")}-Matrix: + A has ${cols_A} cols and B ${rows_B} rows. + Must be equal!`); + } + const C = new Matrix(rows_A, cols_B, (row, col) => { + const A_i = A.row(row); + const B_val = B.values; + let sum = 0; + for (let i = 0, j = col; i < cols_A; ++i, j += cols_B) { + sum += A_i[i] * B_val[j]; + } + return sum; + }); + return C; + } else if (Matrix.isArray(B)) { + let rows = this._rows; + if (B.length !== rows) { + throw new Error(`A.dot(B): A has ${rows} cols and B has ${B.length} rows. Must be equal!`); + } + let C = new Array(rows); + for (let row = 0; row < rows; ++row) { + C[row] = neumair_sum(this.row(row).map((e) => e * B[row])); } + return C; + } else { + throw new Error(`B must be Matrix or Array`); } + } - for (let row = rows - 1; row >= 0; --row) { - let B_row_row = B.entry(row, row); - for (let i = 0; i < row; i++) { - let B_i_row = B.entry(i, row); - let f = B_i_row / B_row_row; - for (let j = i; j < 2 * cols; ++j) { - let B_i_j = B.entry(i, j); - let B_row_j = B.entry(row, j); - B_i_j = B_i_j - B_row_j * f; - B.set_entry(i, j, B_i_j); + /** + * Transposes the current matrix and returns the dot product with {@link B}. + * If {@link B} is an Array or Float64Array then an Array gets returned. + * If {@link B} is a Matrix then a Matrix gets returned. + * @param {(Matrix|Array|Float64Array)} B the right side + * @returns {(Matrix|Array)} + */ + transDot(B) { + if (B instanceof Matrix) { + let A = this; + const [cols_A, rows_A] = A.shape; // transpose matrix + const [rows_B, cols_B] = B.shape; + if (cols_A !== rows_B) { + throw new Error(`A.dot(B): A is a ${[rows_A, cols_A].join(" ⨯ ")}-Matrix, B is a ${B.shape.join(" ⨯ ")}-Matrix: + A has ${cols_A} cols and B ${rows_B} rows, which must be equal!`); + } + // let B = new Matrix(this._cols, this._rows, (row, col) => this.entry(col, row)); + // this.values[row * this._cols + col]; + const C = new Matrix(rows_A, cols_B, (row, col) => { + const A_val = A.values; + const B_val = B.values; + let sum = 0; + for (let i = 0, j = row, k = col; i < cols_A; ++i, j += rows_A, k += cols_B) { + sum += A_val[j] * B_val[k]; } + return sum; + }); + return C; + } else if (Matrix.isArray(B)) { + let rows = this._cols; + if (B.length !== rows) { + throw new Error(`A.dot(B): A has ${rows} cols and B has ${B.length} rows. Must be equal!`); + } + let C = new Array(rows); + for (let row = 0; row < rows; ++row) { + C[row] = neumair_sum(this.col(row).map((e) => e * B[row])); } + return C; + } else { + throw new Error(`B must be Matrix or Array`); } - - return new Matrix(rows, cols, (i, j) => B.entry(i, j + cols)); } /** - * Returns the dot product. If {@link B} is an Array or Float64Array then an Array gets returned. If {@link B} is a Matrix then a Matrix gets returned. + * Returns the dot product with the transposed version of {@link B}. + * If {@link B} is an Array or Float64Array then an Array gets returned. + * If {@link B} is a Matrix then a Matrix gets returned. * @param {(Matrix|Array|Float64Array)} B the right side * @returns {(Matrix|Array)} */ - dot(B) { + dotTrans(B) { if (B instanceof Matrix) { let A = this; - if (A.shape[1] !== B.shape[0]) { - throw new Error(`A.dot(B): A is a ${A.shape.join(" ⨯ ")}-Matrix, B is a ${B.shape.join(" ⨯ ")}-Matrix: - A has ${A.shape[1]} cols and B ${B.shape[0]} rows. - Must be equal!`); + const [rows_A, cols_A] = A.shape; + const [cols_B, rows_B] = B.shape; + if (cols_A !== rows_B) { + throw new Error(`A.dot(B): A is a ${A.shape.join(" ⨯ ")}-Matrix, B is a ${[rows_B, cols_B].join(" ⨯ ")}-Matrix: + A has ${cols_A} cols and B ${rows_B} rows, which must be equal!`); } - let I = A.shape[1]; - let C = new Matrix(A.shape[0], B.shape[1], (row, col) => { + const C = new Matrix(rows_A, cols_B, (row, col) => { const A_i = A.row(row); - const B_i = B.col(col); + const B_i = B.row(col); let sum = 0; - for (let i = 0; i < I; ++i) { + for (let i = 0; i < cols_A; ++i) { sum += A_i[i] * B_i[i]; } return sum; }); return C; - } else if (Array.isArray(B) || B instanceof Float64Array) { + } else if (Matrix.isArray(B)) { let rows = this._rows; if (B.length !== rows) { throw new Error(`A.dot(B): A has ${rows} cols and B has ${B.length} rows. Must be equal!`); @@ -420,15 +538,10 @@ export class Matrix { * @returns {Matrix} */ set_block(offset_row, offset_col, B) { - let [rows, cols] = B.shape; + const rows = Math.min(this._rows - offset_row, B.shape[0]); + const cols = Math.min(this._cols - offset_col, B.shape[1]); for (let row = 0; row < rows; ++row) { - if (row > this._rows) { - continue; - } for (let col = 0; col < cols; ++col) { - if (col > this._cols) { - continue; - } this.set_entry(row + offset_row, col + offset_col, B.entry(row, col)); } } @@ -458,7 +571,7 @@ export class Matrix { end_col = end_col ?? cols; if (end_row <= start_row || end_col <= start_col) { throw new Error(` - end_row must be greater than start_row, and + end_row must be greater than start_row, and end_col must be greater than start_col, but end_row = ${end_row}, start_row = ${start_row}, end_col = ${end_col}, and start_col = ${start_col}!`); } @@ -503,10 +616,8 @@ export class Matrix { _apply_array(f, v) { const data = this.values; const [rows, cols] = this.shape; - for (let row = 0; row < rows; ++row) { - const offset = row * cols; - for (let col = 0; col < cols; ++col) { - const i = offset + col; + for (let i = 0, row = 0; row < rows; ++row) { + for (let col = 0; col < cols; ++col, ++i) { data[i] = f(data[i], v(row, col)); } } @@ -520,68 +631,66 @@ export class Matrix { _apply_colwise_array(values, f) { const data = this.values; const [rows, cols] = this.shape; - for (let row = 0; row < rows; ++row) { - const offset = row * cols; - for (let col = 0; col < cols; ++col) { - const i = offset + col; - data[i] = f(data[i], values[row]); + for (let i = 0, row = 0; row < rows; ++row) { + const val = values[row]; + for (let col = 0; col < cols; ++col, ++i) { + data[i] = f(data[i], val); } } return this; } _apply(value, f) { - let data = this.values; + const data = this.values; + const [rows, cols] = this.shape; if (value instanceof Matrix) { - let [value_rows, value_cols] = value.shape; - let [rows, cols] = this.shape; + const values = value.values; + const [value_rows, value_cols] = value.shape; if (value_rows === 1) { if (cols !== value_cols) { throw new Error(`cols !== value_cols`); } - for (let row = 0; row < rows; ++row) { - for (let col = 0; col < cols; ++col) { - data[row * cols + col] = f(data[row * cols + col], value.entry(0, col)); + for (let i = 0, row = 0; row < rows; ++row) { + for (let col = 0; col < cols; ++col, ++i) { + data[i] = f(data[i], values[col]); } } } else if (value_cols === 1) { if (rows !== value_rows) { throw new Error(`rows !== value_rows`); } - for (let row = 0; row < rows; ++row) { - for (let col = 0; col < cols; ++col) { - data[row * cols + col] = f(data[row * cols + col], value.entry(row, 0)); + for (let i = 0, row = 0; row < rows; ++row) { + const v = values[row]; + for (let col = 0; col < cols; ++col, ++i) { + data[i] = f(data[i], v); } } } else if (rows == value_rows && cols == value_cols) { - for (let row = 0; row < rows; ++row) { - for (let col = 0; col < cols; ++col) { - data[row * cols + col] = f(data[row * cols + col], value.entry(row, col)); - } + for (let i = 0, n = rows * cols; i < n; ++i) { + data[i] = f(data[i], values[i]); } } else { throw new Error(`error`); } - } else if (Array.isArray(value)) { - let rows = this._rows; - let cols = this._cols; + } else if (Matrix.isArray(value)) { if (value.length === rows) { - for (let row = 0; row < rows; ++row) { - for (let col = 0; col < cols; ++col) { - data[row * cols + col] = f(data[row * cols + col], value[row]); + for (let i = 0, row = 0; row < rows; ++row) { + const v = value[row]; + for (let col = 0; col < cols; ++col, ++i) { + data[i] = f(data[i], v); } } } else if (value.length === cols) { - for (let row = 0; row < rows; ++row) { - for (let col = 0; col < cols; ++col) { - data[row * cols + col] = f(data[row * cols + col], value[col]); + for (let i = 0, row = 0; row < rows; ++row) { + for (let col = 0; col < cols; ++col, ++i) { + data[i] = f(data[i], value[col]); } } } else { throw new Error(`error`); } - } else { - for (let i = 0, n = this._rows * this._cols; i < n; ++i) { + } else { // scalar value + for (let i = 0, n = rows * cols; i < n; ++i) { data[i] = f(data[i], value); } } @@ -693,9 +802,9 @@ export class Matrix { this._rows = rows; this._cols = cols; this._data = new Float64Array(rows * cols); - for (let row = 0; row < rows; ++row) { - for (let col = 0; col < cols; ++col) { - this._data[row * cols + col] = value(row, col); + for (let i = 0, row = 0; row < rows; ++row) { + for (let col = 0; col < cols; ++col, ++i) { + this._data[i] = value(row, col); } } return this; @@ -760,7 +869,7 @@ export class Matrix { } /** - * Returns the sum oof all entries of the Matrix. + * Returns the entries of the Matrix. * @returns {Float64Array} */ get values() { @@ -777,12 +886,12 @@ export class Matrix { const rows = this._rows; const cols = this._cols; const result = Float64Array.from({ length: rows }); - for (let row = 0; row < rows; ++row) { - result[row] = 0; - for (let col = 0; col < cols; ++col) { - result[row] += data[row * cols + col]; + for (let i = 0, row = 0; row < rows; ++row) { + let sum = 0; + for (let col = 0; col < cols; ++col, ++i) { + sum += data[i]; } - result[row] /= cols; + result[row] = sum / cols; } return result; } @@ -796,11 +905,11 @@ export class Matrix { const cols = this._cols; const result = Float64Array.from({ length: cols }); for (let col = 0; col < cols; ++col) { - result[col] = 0; - for (let row = 0; row < rows; ++row) { - result[col] += data[row * cols + col]; + let sum = 0; + for (let i = col, row = 0; row < rows; ++row, i += cols) { + sum += data[i]; } - result[col] /= rows; + result[col] = sum / rows; } return result; } @@ -827,10 +936,10 @@ export class Matrix { let d = r.clone(); do { const z = A.dot(d); - const alpha = r.T.dot(r).entry(0, 0) / d.T.dot(z).entry(0, 0); + const alpha = r.transDot(r).entry(0, 0) / d.transDot(z).entry(0, 0); x = x.add(d.mult(alpha)); const r_next = r.sub(z.mult(alpha)); - const beta = r_next.T.dot(r_next).entry(0, 0) / r.T.dot(r).entry(0, 0); + const beta = r_next.transDot(r_next).entry(0, 0) / r.transDot(r).entry(0, 0); d = r_next.add(d.mult(beta)); r = r_next; } while (Math.abs(r.mean) > tol); @@ -853,7 +962,7 @@ export class Matrix { // forward for (let row = 0; row < rows; ++row) { for (let col = 0; col < row - 1; ++col) { - x.set_entry(0, row, x.entry(0, row) - L.entry(row, col) * x.entry(1, col)); + x.sub_entry(0, row, L.entry(row, col) * x.entry(1, col)); } x.set_entry(0, row, x.entry(0, row) / L.entry(row, row)); } @@ -861,7 +970,7 @@ export class Matrix { // backward for (let row = rows - 1; row >= 0; --row) { for (let col = rows - 1; col > row; --col) { - x.set_entry(0, row, x.entry(0, row) - U.entry(row, col) * x.entry(0, col)); + x.sub_entry(0, row, U.entry(row, col) * x.entry(0, col)); } x.set_entry(0, row, x.entry(0, row) / U.entry(row, row)); } @@ -926,9 +1035,8 @@ export class Matrix { * @returns {{U: Matrix, Sigma: Matrix, V: Matrix}} */ static SVD(M, k = 2) { - const MT = M.T; - let MtM = MT.dot(M); - let MMt = M.dot(MT); + let MtM = M.transDot(M); + let MMt = M.dotTrans(M); let { eigenvectors: V, eigenvalues: Sigma } = simultaneous_poweriteration(MtM, k); let { eigenvectors: U } = simultaneous_poweriteration(MMt, k); return { U: U, Sigma: Sigma.map((sigma) => Math.sqrt(sigma)), V: V }; @@ -943,4 +1051,8 @@ export class Matrix { console.log(U,V,B) return { U: U, "Sigma": B, V: V }; */ } + + static isArray(A) { + return Array.isArray(A) || A instanceof Float64Array || A instanceof Float32Array; + } } diff --git a/metrics/euclidean_squared.js b/metrics/euclidean_squared.js index 47a9354..e7b4d12 100755 --- a/metrics/euclidean_squared.js +++ b/metrics/euclidean_squared.js @@ -1,4 +1,3 @@ -import { neumair_sum } from "../numerical/index.js"; /** * Computes the squared euclidean distance (l₂²) between a and b. * @memberof module:metrics @@ -10,12 +9,10 @@ import { neumair_sum } from "../numerical/index.js"; export default function (a, b) { if (a.length != b.length) return undefined; const n = a.length; - const s = new Float64Array(n); + let sum = 0; for (let i = 0; i < n; ++i) { - const x = a[i]; - const y = b[i]; - const x_y = x - y; - s[i] = x_y * x_y; + const a_b = a[i] - b[i]; + sum += a_b * a_b; } - return neumair_sum(s); + return sum; } diff --git a/test/linear_algebra.js b/test/linear_algebra.js index 75fd4e4..03bf37d 100644 --- a/test/linear_algebra.js +++ b/test/linear_algebra.js @@ -2,6 +2,7 @@ import * as druid from "./test_index.js"; import * as assert from "assert"; +const eps = 0.000001; describe("eig", () => { const N = 20; const R = new druid.Randomizer(1212); @@ -11,6 +12,50 @@ describe("eig", () => { it("qr", () => { assert.ok(druid.qr(M)); assert.ok(druid.qr_householder(M)); + const { Q: QM, R: RM } = druid.qr(M); + checkDecomposition(M, QM, RM); + // const { Q: QMh, R: RMh } = druid.qr_householder(M); + // checkDecomposition(M, QMh, RMh); + + const M1 = druid.Matrix.from([[15, 42], [20, 81]]); + const { Q: QM1, R: RM1 } = druid.qr(M1); + approxEqual(QM1.values, Float64Array.from([0.6, -0.8, 0.8, 0.6])); + approxEqual(RM1.values, Float64Array.from([25, 90, 0, 15])); + checkDecomposition(M1, QM1, RM1); + + const M2 = druid.Matrix.from([[12, -51, 4], [6, 167, -68], [-4, 24, -41]]); + const { Q: QM2, R: RM2 } = druid.qr(M2); + approxEqual(QM2.values, Float64Array.from([ + 6/7, -69/175, -58/175, + 3/7, 158/175, 6/175, + -2/7, 6/35, -33/35 + ])); + approxEqual(RM2.values, Float64Array.from([14, 21, -14, 0, 175, -70, 0, 0, 35])); + checkDecomposition(M2, QM2, RM2); + + const M3 = druid.Matrix.from([ + [1, -1, 4], + [1, 4, -2], + [1, 4, 2], + [1, -1, 0] + ]); + const { Q: QM3, R: RM3 } = druid.qr(M3); + approxEqual(QM3.values, Float64Array.from([ + 0.5, -0.5, 0.5, + 0.5, 0.5, -0.5, + 0.5, 0.5, 0.5, + 0.5, -0.5, -0.5 + ])); + approxEqual(RM3.values, Float64Array.from([2, 3, 2, 0, 5, -2, 0, 0, 4])); + checkDecomposition(M3, QM3, RM3); + + const M4 = druid.Matrix.from([ + [7.507, 9.868, 5.057], + [4.482, 2.536, 9.744], + [6.527, 1.094, 3.321] + ]); + const { Q: QM4, R: RM4 } = druid.qr(M4); + checkDecomposition(M4, QM4, RM4); }); it("simultanious poweriteration", () => { assert.ok(druid.simultaneous_poweriteration(M, 2, {max_iterations: 100, seed: 1212, qr: druid.qr})); @@ -18,9 +63,55 @@ describe("eig", () => { let eigs = druid.simultaneous_poweriteration(M, 2, {max_iterations: 100, seed: 1212, qr: druid.qr}); assert.equal(eigs.eigenvectors.length, 2) assert.equal(eigs.eigenvectors[0].length, N) + checkEigs(M, eigs.eigenvalues, eigs.eigenvectors); eigs = druid.simultaneous_poweriteration(M, 6, {max_iterations: 100, seed: 1212, qr: druid.qr}); assert.equal(eigs.eigenvectors.length, 6) assert.equal(eigs.eigenvectors[0].length, N) + checkEigs(M, eigs.eigenvalues, eigs.eigenvectors); + + const A = druid.Matrix.from([[1, 0.1], [0.1, 1]]); + const { eigenvalues: A_val } = druid.simultaneous_poweriteration(A, 2); + approxEqual(A_val, Float64Array.from([1.1, 0.9])); + + // const B = druid.Matrix.from([[3, -1, -1], [-12, 0, 5], [4, -2, -1]]); + const B = druid.Matrix.from([[13, -4, 2], [-4, 11, -2], [2, -2, 8]]); + const { eigenvalues: B_val, eigenvectors: B_vec } = druid.simultaneous_poweriteration(B, 3); + approxEqual(B_val, Float64Array.from([17, 8, 7])); + checkEigs(B, B_val, B_vec); }).timeout(10000); }); + +function checkDecomposition(A, Q, R) { + approxEqual(Q.dot(R).values, A.values); + + // R is upper triangular + const [rows, cols] = A.shape; + for (let i = 0; i < cols; ++i) { + for (let j = 0; j < i; ++j) { + assert.ok(Math.abs(R.entry(i, j)) < eps); + } + } + + // All elements on leading diagonal of R are positive + for (let i = 0; i < Math.min(rows, cols); ++i) { + assert.ok(R.entry(i, i) >= 0); + } +} + +function checkEigs(A, vals, vecs) { + vals.forEach((val, i) => { + const proj1 = A.dot(druid.Matrix.from(vecs[i], 'col')); + const proj2 = vecs[i].map(d => d * val); + approxEqual(proj1.values, proj2, 0.005); + }); +} + +function approxEqual(a, b, e = eps) { + // assert.deepEqual(a, b); + const N = a.length; + assert.ok(a.length === b.length); + for (let i = 0; i < N; ++i) { + assert.ok(Math.abs(a[i] - b[i]) < e, a[i] + ' != ' + b[i]); + } +} \ No newline at end of file diff --git a/test/matrix.js b/test/matrix.js index 0a1328c..b7bbf3b 100644 --- a/test/matrix.js +++ b/test/matrix.js @@ -27,6 +27,58 @@ describe("Matrix", () => { } assert.deepEqual(A.dot(B).shape, [5, 2]); assert.throws(() => B.dot(A), Error, "error thrown"); + + const D = new druid.Matrix(2, 3); + D._data = Float64Array.from([0, 1, 2, 3, 4, 5]); + + assert.ok(D.dot(D.T)); + assert.ok(D.dotTrans(D)); + const D_dot_DT = Float64Array.from([5, 14, 14, 50]); + assert.deepEqual(D.dot(D.T).values, D_dot_DT); + assert.deepEqual(D.dotTrans(D).values, D_dot_DT); + + assert.ok(D.T.dot(D)); + assert.ok(D.transDot(D)); + const DT_dot_D = Float64Array.from([9, 12, 15, 12, 17, 22, 15, 22, 29]); + assert.deepEqual(D.T.dot(D).values, DT_dot_D); + assert.deepEqual(D.transDot(D).values, DT_dot_D); + }); + + it("Matrix inversion", () => { + const A = new druid.Matrix(2, 2); + A._data = Float64Array.from([2, 3, 4, 7]); + assert.ok(A.inverse()); + assert.deepEqual(A.inverse().values, Float64Array.from([3.5, -1.5, -2, 1])); + + const B = new druid.Matrix(3, 3); + B._data = Float64Array.from([1, 4, 7, 3, 0, 5, -1, 9, 11]); + deepEqual(B.inverse().values, Float64Array.from([5.625, -2.375, -2.5, 4.75, -2.25, -2, -3.375, 1.625, 1.5])); + + const C = new druid.Matrix(3, 3); + C._data = Float64Array.from([2, -1, 0, -1, 2, -1, 0, -1, 2]); + deepEqual(C.inverse().values, Float64Array.from([3/4, 1/2, 1/4, 1/2, 1, 1/2, 1/4, 1/2, 3/4])); + + const D = new druid.Matrix(3, 3); + D._data = Float64Array.from([1, 0, 0, 0, 0, 1, 0, 1, 0]); + assert.deepEqual(D.inverse().values, Float64Array.from([1, 0, 0, 0, 0, 1, 0, 1, 0])); + + const E = new druid.Matrix(3, 3); + E._data = Float64Array.from([1, 0, 0, 0, -1, -1, 0, 0, 1]); + assert.deepEqual(E.inverse().values, Float64Array.from([1, 0, 0, 0, -1, -1, 0, 0, 1])); + + const F = new druid.Matrix(4, 4); + F._data = Float64Array.from([ + 0, 1, 0, 788, + -1, 0, 0, 692, + 0, 0, 1, 0, + 0, 0, 0, 1 + ]); + assert.deepEqual(F.inverse().values, Float64Array.from([ + 0, -1, 0, 692, + 1, 0, 0, -788, + 0, 0, 1, 0, + 0, 0, 0, 1 + ])); }); it("LU decomposition", () => { @@ -51,6 +103,14 @@ describe("norm", () => { it("norm", () => { const v = druid.normalize(Float64Array.from({length: 100}, () => Math.random() * 100 - 50)); assert.equal(Math.abs(druid.norm(v) - 1) < 1e-12, 1) - + }) -}) \ No newline at end of file +}) + +function deepEqual(a, b) { + assert.ok(a.length === b.length); + const N = a.length; + for (let i = 0; i < N; ++i) { + assert.ok(Math.abs(a[i] - b[i]) < 0.0001, (a + ' ~= ' + b)); + } +} \ No newline at end of file diff --git a/util/randomizer.js b/util/randomizer.js index eb75218..c6387f0 100755 --- a/util/randomizer.js +++ b/util/randomizer.js @@ -95,6 +95,21 @@ export class Randomizer { return y >>> 0; } + gauss_random() { + let x, y, r; + if (this._val != null) { + x = this._val, this._val = null; + return x; + } else do { + x = 2 * this.random - 1; + y = 2 * this.random - 1; + r = x * x + y * y; + } while (!r || r > 1); + const c = Math.sqrt(-2 * Math.log(r) / r); + this._val = y * c; // cache this for next function call for efficiency + return x * c; + } + /** * Returns samples from an input Matrix or Array. * @param {Matrix|Array|Float64Array} A - The input Matrix or Array.