saehm · saehm · Aug 16, 2022 · Apr 25, 2022 · Apr 25, 2022 · May 13, 2022
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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 };