From 47a425bbba79992cb1f78cfeecb0f5a3e1a97dd8 Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Mon, 25 Apr 2022 14:51:25 +0200
Subject: [PATCH 01/22] Bugfix UMAP and small optimizations
---
dimred/UMAP.js | 38 +++++++++++++++++---------------------
1 file changed, 17 insertions(+), 21 deletions(-)
diff --git a/dimred/UMAP.js b/dimred/UMAP.js
index d808f99..77cd057 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 curr_dist = distances[i];
+ const rho = rhos[i].value;
+ 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;
@@ -137,7 +139,7 @@ export class UMAP extends DR {
if (index > 0) {
rhos.push(non_zero_dist[index - 1]);
if (interpolation > SMOOTH_K_TOLERANCE) {
- rhos[i].value += interpolation * (non_zero_dist[index].value - non_zero_dist[index - 1]);
+ rhos[i].value += interpolation * (non_zero_dist[index].value - non_zero_dist[index - 1].value);
}
} else {
rhos[i].value = interpolation * non_zero_dist[0].value;
@@ -148,7 +150,7 @@ export class UMAP extends DR {
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 - rhos[i].value;
psum += d > 0 ? Math.exp(-(d / mid)) : 1;
}
if (Math.abs(psum - target) < SMOOTH_K_TOLERANCE) {
@@ -168,7 +170,7 @@ export class UMAP extends DR {
const mean_ithd = search_result.reduce((a, b) => a + b.value, 0) / search_result.length;
//let mean_d = null;
- if (rhos[i] > 0) {
+ if (rhos[i].value > 0) {
if (sigmas[i] < MIN_K_DIST_SCALE * mean_ithd) {
sigmas[i] = MIN_K_DIST_SCALE * mean_ithd;
}
@@ -223,9 +225,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;
}
@@ -356,12 +360,8 @@ export class UMAP extends DR {
}
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);
+ 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];
@@ -377,12 +377,8 @@ export class UMAP extends DR {
}
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);
+ current[d] += grad_d;
+ other[d] -= grad_d;
}
}
epoch_of_next_negative_sample[i] += n_neg_samples * epochs_per_negative_sample[i];
From a6fb9742fd1ed408b3988593613d810b9b7c43ce Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Mon, 25 Apr 2022 17:15:28 +0200
Subject: [PATCH 02/22] Rewrite Matrix.dot to avoid copying column to new
array; bugfix in set_block
---
matrix/Matrix.js | 28 ++++++++++++----------------
1 file changed, 12 insertions(+), 16 deletions(-)
diff --git a/matrix/Matrix.js b/matrix/Matrix.js
index bbfc1f0..f483075 100755
--- a/matrix/Matrix.js
+++ b/matrix/Matrix.js
@@ -312,18 +312,19 @@ export class Matrix {
dot(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.
+ 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!`);
}
- 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_val = B.values;
let sum = 0;
- for (let i = 0; i < I; ++i) {
- sum += A_i[i] * B_i[i];
+ for (let i = 0, j = col; i < cols_A; ++i, j += cols_B) {
+ sum += A_i[i] * B_val[j];
}
return sum;
});
@@ -420,15 +421,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 +454,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}!`);
}
From ec354b09c8d551d0cf1a83d6181ce35b7f93804a Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Fri, 13 May 2022 11:40:19 +0200
Subject: [PATCH 03/22] UMAP: Skip updating current / other, when grad_coeff is
0 (has no effect anyways)
---
dimred/UMAP.js | 26 ++++++++++++--------------
1 file changed, 12 insertions(+), 14 deletions(-)
diff --git a/dimred/UMAP.js b/dimred/UMAP.js
index 77cd057..7d0e2f5 100755
--- a/dimred/UMAP.js
+++ b/dimred/UMAP.js
@@ -354,14 +354,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;
- current[d] += grad_d;
- other[d] -= grad_d;
+ 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,17 +368,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;
- current[d] += grad_d;
- other[d] -= grad_d;
- }
}
epoch_of_next_negative_sample[i] += n_neg_samples * epochs_per_negative_sample[i];
}
From b65ea88ba4cd169f0aa819ca13d64e59cfd1fec9 Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Tue, 17 May 2022 12:54:18 +0200
Subject: [PATCH 04/22] Bugfix in TSNE
---
dimred/TSNE.js | 13 +++++--------
matrix/Matrix.js | 24 ++++++++++++++++++++++++
2 files changed, 29 insertions(+), 8 deletions(-)
diff --git a/dimred/TSNE.js b/dimred/TSNE.js
index 210d0b5..6c50d8d 100755
--- a/dimred/TSNE.js
+++ b/dimred/TSNE.js
@@ -62,7 +62,7 @@ export class TSNE extends DR {
this._gains = new Matrix(N, D, 1);
// search for fitting sigma
- let prow = new Float64Array(N)
+ let prow = new Float64Array(N);
const tol = 1e-4;
const maxtries = 50;
for (let i = 0; i < N; ++i) {
@@ -99,10 +99,7 @@ export class TSNE extends DR {
if (Math.abs(Hhere - Htarget) < tol) done = true;
if (num >= maxtries) done = true;
}
-
- for (let j = 0; j < N; ++j) {
- P.set_entry(i, j, prow[j]);
- }
+ P.set_row(i, prow);
}
//compute probabilities
@@ -216,14 +213,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/matrix/Matrix.js b/matrix/Matrix.js
index f483075..c19874c 100755
--- a/matrix/Matrix.js
+++ b/matrix/Matrix.js
@@ -215,6 +215,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}
From dce969cf71822ea5db9acacb359a1f84834fa5e0 Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Thu, 19 May 2022 13:50:07 +0200
Subject: [PATCH 05/22] Matrix class can now also deal with Float32Arrays
---
matrix/Matrix.js | 16 ++++++++++------
1 file changed, 10 insertions(+), 6 deletions(-)
diff --git a/matrix/Matrix.js b/matrix/Matrix.js
index c19874c..65211d6 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];
@@ -353,7 +353,7 @@ export class Matrix {
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!`);
@@ -780,7 +780,7 @@ export class Matrix {
}
/**
- * Returns the sum oof all entries of the Matrix.
+ * Returns the entries of the Matrix.
* @returns {Float64Array}
*/
get values() {
@@ -963,4 +963,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;
+ }
}
From 8a11f1d0fbeaf0d81b3671f6984093f810245818 Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Mon, 23 May 2022 15:37:39 +0200
Subject: [PATCH 06/22] Small optimizations in UMAP
---
dimred/UMAP.js | 32 +++++++++++++++++---------------
1 file changed, 17 insertions(+), 15 deletions(-)
diff --git a/dimred/UMAP.js b/dimred/UMAP.js
index 7d0e2f5..ed7c497 100755
--- a/dimred/UMAP.js
+++ b/dimred/UMAP.js
@@ -86,8 +86,8 @@ export class UMAP extends DR {
*/
_compute_membership_strengths(distances, sigmas, rhos) {
for (let i = 0, n = distances.length; i < n; ++i) {
+ const rho = rhos[i];
const curr_dist = distances[i];
- const rho = rhos[i].value;
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;
@@ -125,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].value);
+ 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].value;
+ const d = search_result[j].value - rho;
psum += d > 0 ? Math.exp(-(d / mid)) : 1;
}
if (Math.abs(psum - target) < SMOOTH_K_TOLERANCE) {
@@ -166,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].value > 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,
From 7e274d6d82feece4e8d362d6f80fb476b700c715 Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Tue, 24 May 2022 14:29:08 +0200
Subject: [PATCH 07/22] Speed up matrix operations
---
dimred/PCA.js | 3 +-
matrix/Matrix.js | 86 +++++++++++++++++++++++-------------------------
2 files changed, 42 insertions(+), 47 deletions(-)
diff --git a/dimred/PCA.js b/dimred/PCA.js
index 9d314d1..ba16f62 100755
--- a/dimred/PCA.js
+++ b/dimred/PCA.js
@@ -57,8 +57,7 @@ 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 X_cent = X.sub(X.meanCols);
const C = X_cent.transpose().dot(X_cent);
const { eigenvectors: V } = simultaneous_poweriteration(C, d, eig_args);
this.V = Matrix.from(V).transpose();
diff --git a/matrix/Matrix.js b/matrix/Matrix.js
index 65211d6..41726d5 100755
--- a/matrix/Matrix.js
+++ b/matrix/Matrix.js
@@ -523,10 +523,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));
}
}
@@ -540,68 +538,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 arr = 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], arr[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 = arr[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], arr[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);
}
}
@@ -713,9 +709,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;
@@ -797,12 +793,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;
}
@@ -816,11 +812,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;
}
From 59bcf7fc596eb2d2bc0e3ba0d62ab987a9e813e8 Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Tue, 24 May 2022 15:34:49 +0200
Subject: [PATCH 08/22] Bugfix when directly accessing matrix values
---
matrix/Matrix.js | 8 ++++----
1 file changed, 4 insertions(+), 4 deletions(-)
diff --git a/matrix/Matrix.js b/matrix/Matrix.js
index 41726d5..f278825 100755
--- a/matrix/Matrix.js
+++ b/matrix/Matrix.js
@@ -551,7 +551,7 @@ export class Matrix {
const data = this.values;
const [rows, cols] = this.shape;
if (value instanceof Matrix) {
- const arr = value.values();
+ const values = value.values;
const [value_rows, value_cols] = value.shape;
if (value_rows === 1) {
if (cols !== value_cols) {
@@ -559,7 +559,7 @@ export class Matrix {
}
for (let i = 0, row = 0; row < rows; ++row) {
for (let col = 0; col < cols; ++col, ++i) {
- data[i] = f(data[i], arr[col]);
+ data[i] = f(data[i], values[col]);
}
}
} else if (value_cols === 1) {
@@ -567,14 +567,14 @@ export class Matrix {
throw new Error(`rows !== value_rows`);
}
for (let i = 0, row = 0; row < rows; ++row) {
- const v = arr[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 i = 0, n = rows * cols; i < n; ++i) {
- data[i] = f(data[i], arr[i]);
+ data[i] = f(data[i], values[i]);
}
} else {
throw new Error(`error`);
From 2ffebc4f8a51345ae59743201a1e49ecdada415c Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Wed, 25 May 2022 09:45:44 +0200
Subject: [PATCH 09/22] Add convenience function that performs
Matrix.transpose().dot()
---
dimred/LLE.js | 2 +-
dimred/PCA.js | 2 +-
linear_algebra/poweriteration_m.js | 2 +-
linear_algebra/qr_givens.js | 2 +-
linear_algebra/sapi.js | 4 +--
linear_algebra/svrg.js | 6 ++--
matrix/Matrix.js | 45 ++++++++++++++++++++++++++++--
7 files changed, 52 insertions(+), 11 deletions(-)
diff --git a/dimred/LLE.js b/dimred/LLE.js
index 06a9654..8b2ca9f 100755
--- a/dimred/LLE.js
+++ b/dimred/LLE.js
@@ -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/PCA.js b/dimred/PCA.js
index ba16f62..a5071d4 100755
--- a/dimred/PCA.js
+++ b/dimred/PCA.js
@@ -58,7 +58,7 @@ export class PCA extends DR {
const { d, eig_args } = this._parameters;
const X = this.X;
const X_cent = X.sub(X.meanCols);
- const C = X_cent.transpose().dot(X_cent);
+ 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/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_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/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 f278825..dda36d1 100755
--- a/matrix/Matrix.js
+++ b/matrix/Matrix.js
@@ -368,6 +368,47 @@ export class Matrix {
}
}
+ /**
+ * 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)}
+ */
+ 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`);
+ }
+ }
+
/**
* Computes the outer product from {@link this} and {@link B}.
* @param {Matrix} B
@@ -843,10 +884,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);
From 3922f7a21f49aa946d0f10133d579bba71453cba Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Wed, 25 May 2022 10:46:44 +0200
Subject: [PATCH 10/22] Add convenience function that performs
Matrix.dot(B.transposed())
---
dimred/LDA.js | 4 ++--
dimred/LSP.js | 5 ++---
dimred/LTSA.js | 5 ++---
dimred/MLLE.js | 8 ++++----
matrix/Matrix.js | 45 ++++++++++++++++++++++++++++++++++++++++++++-
5 files changed, 54 insertions(+), 13 deletions(-)
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/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..05ae436 100755
--- a/dimred/LTSA.js
+++ b/dimred/LTSA.js
@@ -55,13 +55,12 @@ 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) {
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/matrix/Matrix.js b/matrix/Matrix.js
index dda36d1..1794ee6 100755
--- a/matrix/Matrix.js
+++ b/matrix/Matrix.js
@@ -369,7 +369,9 @@ export class Matrix {
}
/**
- * 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.
+ * 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)}
*/
@@ -409,6 +411,47 @@ export class Matrix {
}
}
+ /**
+ * 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)}
+ */
+ dotTrans(B) {
+ if (B instanceof Matrix) {
+ let A = this;
+ 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!`);
+ }
+ 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 * rows_B; i < cols_A; ++i, ++j) {
+ 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`);
+ }
+ }
+
/**
* Computes the outer product from {@link this} and {@link B}.
* @param {Matrix} B
From 246df159d433382818ca8955e2b1f44af293b3a5 Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Fri, 27 May 2022 13:08:37 +0200
Subject: [PATCH 11/22] Bugfix in TriMap: parameters passed to PCA, argsort
---
dimred/TriMap.js | 12 ++++--------
1 file changed, 4 insertions(+), 8 deletions(-)
diff --git a/dimred/TriMap.js b/dimred/TriMap.js
index 65ae919..b82b260 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]);
}
/**
From a6c62026d2594fb5fa50661e15aa632a997f330c Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Fri, 27 May 2022 14:49:59 +0200
Subject: [PATCH 12/22] Add test cases for matrix methods
---
matrix/Matrix.js | 6 +++---
test/matrix.js | 24 +++++++++++++++++++++++-
2 files changed, 26 insertions(+), 4 deletions(-)
diff --git a/matrix/Matrix.js b/matrix/Matrix.js
index 1794ee6..eb82434 100755
--- a/matrix/Matrix.js
+++ b/matrix/Matrix.js
@@ -429,10 +429,10 @@ export class Matrix {
}
const C = new Matrix(rows_A, cols_B, (row, col) => {
const A_i = A.row(row);
- const B_val = B.values;
+ const B_i = B.row(col);
let sum = 0;
- for (let i = 0, j = col * rows_B; i < cols_A; ++i, ++j) {
- sum += A_i[i] * B_val[j];
+ for (let i = 0; i < cols_A; ++i) {
+ sum += A_i[i] * B_i[i];
}
return sum;
});
diff --git a/test/matrix.js b/test/matrix.js
index 0a1328c..7655db4 100644
--- a/test/matrix.js
+++ b/test/matrix.js
@@ -27,6 +27,28 @@ 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]));
});
it("LU decomposition", () => {
@@ -51,6 +73,6 @@ 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
From 02600ccc1fba269eaeabef8dc5d6747729b716d5 Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Mon, 30 May 2022 12:01:39 +0200
Subject: [PATCH 13/22] Rewrite matrix inversion since it failed in some cases
---
matrix/Matrix.js | 115 +++++++++++++++++++++++++----------------------
test/matrix.js | 40 ++++++++++++++++-
2 files changed, 101 insertions(+), 54 deletions(-)
diff --git a/matrix/Matrix.js b/matrix/Matrix.js
index eb82434..e50a847 100755
--- a/matrix/Matrix.js
+++ b/matrix/Matrix.js
@@ -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
@@ -263,69 +279,62 @@ 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);
- }
- }
- h++;
- k++;
+ if (max_val === 0) {
+ throw new Error('Cannot compute inverse of Matrix, determinant is zero');
}
- }
-
- 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);
+ // Swap maximum row with current row
+ if (max_idx !== col) {
+ A.swap_rows(col, max_idx);
+ B.swap_rows(col, max_idx);
}
- }
- 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);
+ // 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) {
+ const A_row = A.row(row);
+ const B_row = B.row(row);
+ if (row !== col) {
+ // eliminate value at column c and row r
+ 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,
+ // divide each value on row r with the value at Acc
+ const f = A_col[col]
+ for (let s = col; s < cols; ++s) {
+ A_row[s] /= f;
+ }
+ for (let s = 0; s < cols; ++s) {
+ B_row[s] /= f;
+ }
}
}
}
-
- return new Matrix(rows, cols, (i, j) => B.entry(i, j + cols));
+ return B;
}
/**
diff --git a/test/matrix.js b/test/matrix.js
index 7655db4..b7bbf3b 100644
--- a/test/matrix.js
+++ b/test/matrix.js
@@ -49,6 +49,36 @@ describe("Matrix", () => {
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", () => {
@@ -75,4 +105,12 @@ describe("norm", () => {
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
From 052e83ddec54ced49b558ea598bd6f1d5cce11dd Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Tue, 31 May 2022 11:49:15 +0200
Subject: [PATCH 14/22] Small changes to matrix inversion
---
matrix/Matrix.js | 10 +++++-----
1 file changed, 5 insertions(+), 5 deletions(-)
diff --git a/matrix/Matrix.js b/matrix/Matrix.js
index e50a847..16e6e43 100755
--- a/matrix/Matrix.js
+++ b/matrix/Matrix.js
@@ -307,10 +307,10 @@ export class Matrix {
const A_col = A.row(col);
const B_col = B.row(col);
for (let row = 0; row < rows; ++row) {
- const A_row = A.row(row);
- const B_row = B.row(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
@@ -324,12 +324,12 @@ export class Matrix {
} else {
// normalize value at Acc to 1,
// divide each value on row r with the value at Acc
- const f = A_col[col]
+ const f = A_col[col];
for (let s = col; s < cols; ++s) {
- A_row[s] /= f;
+ A_col[s] /= f;
}
for (let s = 0; s < cols; ++s) {
- B_row[s] /= f;
+ B_col[s] /= f;
}
}
}
From 9a7df9bef1c9f5e72de59dfda9e30bab3ace745d Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Tue, 31 May 2022 14:47:47 +0200
Subject: [PATCH 15/22] Use convenience functions to change matrix elements
(add, sub)
---
dimred/ISOMAP.js | 4 +++-
dimred/LLE.js | 2 +-
dimred/LTSA.js | 2 +-
dimred/TSNE.js | 2 +-
dimred/TriMap.js | 6 +++---
matrix/Matrix.js | 8 ++++----
6 files changed, 13 insertions(+), 11 deletions(-)
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/LLE.js b/dimred/LLE.js
index 8b2ca9f..cfa64a0 100755
--- a/dimred/LLE.js
+++ b/dimred/LLE.js
@@ -53,7 +53,7 @@ export class LLE extends DR {
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;
diff --git a/dimred/LTSA.js b/dimred/LTSA.js
index 05ae436..12219d2 100755
--- a/dimred/LTSA.js
+++ b/dimred/LTSA.js
@@ -64,7 +64,7 @@ export class LTSA extends DR {
.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/TSNE.js b/dimred/TSNE.js
index 6c50d8d..0093384 100755
--- a/dimred/TSNE.js
+++ b/dimred/TSNE.js
@@ -192,7 +192,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)));
}
}
}
diff --git a/dimred/TriMap.js b/dimred/TriMap.js
index b82b260..c4003c5 100755
--- a/dimred/TriMap.js
+++ b/dimred/TriMap.js
@@ -310,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/matrix/Matrix.js b/matrix/Matrix.js
index 16e6e43..158088e 100755
--- a/matrix/Matrix.js
+++ b/matrix/Matrix.js
@@ -322,8 +322,8 @@ export class Matrix {
}
}
} else {
- // normalize value at Acc to 1,
- // divide each value on row r with the value at Acc
+ // 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;
@@ -962,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));
}
@@ -970,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));
}
From 9674722fa98dad897e59961dcbc169ffd489d2c0 Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Thu, 2 Jun 2022 16:44:25 +0200
Subject: [PATCH 16/22] Add tests for QR method
---
test/linear_algebra.js | 54 ++++++++++++++++++++++++++++++++++++++++++
1 file changed, 54 insertions(+)
diff --git a/test/linear_algebra.js b/test/linear_algebra.js
index 75fd4e4..f293cc2 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.00000001;
describe("eig", () => {
const N = 20;
const R = new druid.Randomizer(1212);
@@ -11,6 +12,34 @@ 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 A = druid.Matrix.from([[12, -51, 4], [6, 167, -68], [-4, 24, -41]]);
+ const { Q, R } = druid.qr(A);
+ approxEqual(Q.values, Float64Array.from([
+ 6/7, -69/175, -58/175,
+ 3/7, 158/175, 6/175,
+ -2/7, 6/35, -33/35
+ ]));
+ approxEqual(R.values, Float64Array.from([14, 21, -14, 0, 175, -70, 0, 0, 35]));
+
+ const B = druid.Matrix.from([
+ [1, -1, 4],
+ [1, 4, -2],
+ [1, 4, 2],
+ [1, -1, 0]
+ ]);
+ const { Q: QB, R: RB } = druid.qr(B);
+ approxEqual(QB.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(RB.values, Float64Array.from([2, 3, 2, 0, 5, -2, 0, 0, 4]));
});
it("simultanious poweriteration", () => {
assert.ok(druid.simultaneous_poweriteration(M, 2, {max_iterations: 100, seed: 1212, qr: druid.qr}));
@@ -24,3 +53,28 @@ describe("eig", () => {
assert.equal(eigs.eigenvectors[0].length, N)
}).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 < rows; ++i) {
+ for (let j = 0; j < i && j < cols; ++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 approxEqual(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]) < eps);
+ }
+}
\ No newline at end of file
From e4e79a4b8a260509f70dc855831e166d8bd54058 Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Tue, 7 Jun 2022 13:49:22 +0200
Subject: [PATCH 17/22] Small changes to code
---
dimred/LLE.js | 2 +-
dimred/OAP.js | 4 ++--
linear_algebra/simultaneous_poweriteration.js | 2 +-
matrix/Matrix.js | 5 ++---
4 files changed, 6 insertions(+), 7 deletions(-)
diff --git a/dimred/LLE.js b/dimred/LLE.js
index cfa64a0..de322c3 100755
--- a/dimred/LLE.js
+++ b/dimred/LLE.js
@@ -49,7 +49,7 @@ 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) {
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/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/matrix/Matrix.js b/matrix/Matrix.js
index 158088e..7b72b4d 100755
--- a/matrix/Matrix.js
+++ b/matrix/Matrix.js
@@ -1035,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 };
From b322a3d97d8e15960563af1f2e14409378de0bbc Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Wed, 8 Jun 2022 15:35:53 +0200
Subject: [PATCH 18/22] Add test cases for QR decomposition and eigenvalues
---
linear_algebra/qr.js | 4 ++-
test/linear_algebra.js | 60 +++++++++++++++++++++++++++++++-----------
2 files changed, 47 insertions(+), 17 deletions(-)
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/test/linear_algebra.js b/test/linear_algebra.js
index f293cc2..0082878 100644
--- a/test/linear_algebra.js
+++ b/test/linear_algebra.js
@@ -2,7 +2,7 @@
import * as druid from "./test_index.js";
import * as assert from "assert";
-const eps = 0.00000001;
+const eps = 0.0000001;
describe("eig", () => {
const N = 20;
const R = new druid.Randomizer(1212);
@@ -14,32 +14,48 @@ describe("eig", () => {
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 { Q: QMh, R: RMh } = druid.qr_householder(M);
+ // checkDecomposition(M, QMh, RMh);
- const A = druid.Matrix.from([[12, -51, 4], [6, 167, -68], [-4, 24, -41]]);
- const { Q, R } = druid.qr(A);
- approxEqual(Q.values, Float64Array.from([
+ 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(R.values, Float64Array.from([14, 21, -14, 0, 175, -70, 0, 0, 35]));
+ approxEqual(RM2.values, Float64Array.from([14, 21, -14, 0, 175, -70, 0, 0, 35]));
+ checkDecomposition(M2, QM2, RM2);
- const B = druid.Matrix.from([
+ const M3 = druid.Matrix.from([
[1, -1, 4],
[1, 4, -2],
[1, 4, 2],
[1, -1, 0]
]);
- const { Q: QB, R: RB } = druid.qr(B);
- approxEqual(QB.values, Float64Array.from([
+ 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(RB.values, Float64Array.from([2, 3, 2, 0, 5, -2, 0, 0, 4]));
+ 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}));
@@ -51,6 +67,17 @@ describe("eig", () => {
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)
+
+ 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 { eigenvalues: B_val, eigenvectors: B_vec } = druid.simultaneous_poweriteration(B, 3);
+ // const B_val2 = druid.Matrix.from(B_vec);
+ // // approxEqual(B_val, Float64Array.from([3, 1, 0]));
+ // approxEqual(B.dot(B_vec), [0, 0]);
+
}).timeout(10000);
});
@@ -59,19 +86,20 @@ function checkDecomposition(A, Q, R) {
// R is upper triangular
const [rows, cols] = A.shape;
- for (let i = 0; i < rows; ++i) {
- for (let j = 0; j < i && j < cols; ++j) {
+ 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);
- // }
+ for (let i = 0; i < Math.min(rows, cols); ++i) {
+ assert.ok(R.entry(i, i) >= 0);
+ }
}
function approxEqual(a, b) {
+ // assert.deepEqual(a, b);
const N = a.length;
assert.ok(a.length === b.length);
for (let i = 0; i < N; ++i) {
From 2e2ed6ab94520509a0556ff1ece5b6b04baa89bd Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Thu, 23 Jun 2022 10:13:58 +0200
Subject: [PATCH 19/22] Test cases for eigenvalue decomposition
---
test/linear_algebra.js | 25 +++++++++++++++++--------
1 file changed, 17 insertions(+), 8 deletions(-)
diff --git a/test/linear_algebra.js b/test/linear_algebra.js
index 0082878..03bf37d 100644
--- a/test/linear_algebra.js
+++ b/test/linear_algebra.js
@@ -2,7 +2,7 @@
import * as druid from "./test_index.js";
import * as assert from "assert";
-const eps = 0.0000001;
+const eps = 0.000001;
describe("eig", () => {
const N = 20;
const R = new druid.Randomizer(1212);
@@ -63,21 +63,22 @@ 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 { eigenvalues: B_val, eigenvectors: B_vec } = druid.simultaneous_poweriteration(B, 3);
- // const B_val2 = druid.Matrix.from(B_vec);
- // // approxEqual(B_val, Float64Array.from([3, 1, 0]));
- // approxEqual(B.dot(B_vec), [0, 0]);
-
+ 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);
});
@@ -98,11 +99,19 @@ function checkDecomposition(A, Q, R) {
}
}
-function approxEqual(a, b) {
+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]) < eps);
+ assert.ok(Math.abs(a[i] - b[i]) < e, a[i] + ' != ' + b[i]);
}
}
\ No newline at end of file
From 7397cea33b00c28840f40eaa5a88f6cd59865582 Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Tue, 19 Jul 2022 13:12:46 +0200
Subject: [PATCH 20/22] Initialize TSNE with random Gaussian distribution
---
dimred/TSNE.js | 2 +-
util/randomizer.js | 15 +++++++++++++++
2 files changed, 16 insertions(+), 1 deletion(-)
diff --git a/dimred/TSNE.js b/dimred/TSNE.js
index 0093384..5d42fb2 100755
--- a/dimred/TSNE.js
+++ b/dimred/TSNE.js
@@ -26,7 +26,7 @@ export class TSNE extends DR {
super(X, { perplexity: 50, epsilon: 10, d: 2, metric: euclidean, 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;
}
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.
From 677d0c2980029957a829fecbb294c68b2d6db67f Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Tue, 19 Jul 2022 16:33:50 +0200
Subject: [PATCH 21/22] Optimization to TSNE
---
dimred/TSNE.js | 51 +++++++++++++++++++++++++-------------------------
1 file changed, 25 insertions(+), 26 deletions(-)
diff --git a/dimred/TSNE.js b/dimred/TSNE.js
index 5d42fb2..0a01a02 100755
--- a/dimred/TSNE.js
+++ b/dimred/TSNE.js
@@ -56,63 +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) {
+ prow[j] /= psum;
}
- P.set_row(i, prow);
}
- //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;
}
From 225ffadf53db912af8205c6ea463434fa9fcd358 Mon Sep 17 00:00:00 2001
From: jkehrer
Date: Wed, 20 Jul 2022 16:10:35 +0200
Subject: [PATCH 22/22] TSNE now uses euclidean_squared as default metric
---
dimred/TSNE.js | 6 +++---
metrics/euclidean_squared.js | 11 ++++-------
2 files changed, 7 insertions(+), 10 deletions(-)
diff --git a/dimred/TSNE.js b/dimred/TSNE.js
index 0a01a02..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,12 +18,12 @@ 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.gauss_random() * 1e-4);
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 (l22) 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;
}