Add balancing method in matrix class (#724)

Author: ishii-norimi

lib/util/matrix.js

@@ -3420,6 +3420,60 @@ export default class Matrix {
 		return h
 	}
 
+	/**
+	 * Returns a doubly stochastic matrix.
+	 *
+	 * @returns {[number[], Matrix, number[]]} Doubly stochastic matrix
+	 */
+	balancing() {
+		if (!this.isSquare()) {
+			throw new MatrixException('Doubly stochastic matrix only defined for square matrix.', this)
+		}
+		if (this._value.some(v => v <= 0)) {
+			throw new MatrixException('Doubly stochastic matrix only calculate for non negative matrix.', this)
+		}
+		if (this.rows === 1) {
+			return [[this._value[0]], new Matrix(1, 1, 1), [1]]
+		}
+		return this.balancingSinkhornKnopp()
+	}
+
+	/**
+	 * Returns a doubly stochastic matrix by Sinkhorn-Knopp algorithm.
+	 *
+	 * @returns {[number[], Matrix, number[]]} Doubly stochastic matrix
+	 */
+	balancingSinkhornKnopp() {
+		if (!this.isSquare()) {
+			throw new MatrixException('Doubly stochastic matrix only defined for square matrix.', this)
+		}
+		if (this._value.some(v => v <= 0)) {
+			throw new MatrixException('Doubly stochastic matrix only calculate for non negative matrix.', this)
+		}
+		const n = this.rows
+		const d1 = Array(n).fill(1)
+		const d2 = Array(n).fill(1)
+		const a = this.copy()
+		let maxCount = 1.0e4
+		while (maxCount-- > 0) {
+			const s1 = a.sum(1)
+			a.div(s1)
+			for (let i = 0; i < n; i++) {
+				d1[i] *= s1.at(i, 0)
+			}
+			const s2 = a.sum(0)
+			a.div(s2)
+			for (let i = 0; i < n; i++) {
+				d2[i] *= s2.at(0, i)
+			}
+			const e = s1.reduce((s, v) => s + Math.abs(v - 1), 0) + s2.reduce((s, v) => s + Math.abs(v - 1), 0)
+			if (e < 1.0e-8) {
+				break
+			}
+		}
+		return [d1, a, d2]
+	}
+
 	/**
 	 * Returns a LU decomposition.
 	 *

tests/lib/util/matrix.test.js

@@ -5405,6 +5405,82 @@ describe('Matrix', () => {
 		})
 	})
 
+	describe('balancing', () => {
+		test.each([1, 2, 3, 4])('square %d', n => {
+			const mat = Matrix.random(n, n, 0.1, 10)
+			const [d1, a, d2] = mat.balancing()
+
+			expect(a.sizes).toEqual([n, n])
+			expect(d1).toHaveLength(n)
+			expect(d2).toHaveLength(n)
+			const s0 = a.sum(0).value
+			const s1 = a.sum(1).value
+			for (let i = 0; i < n; i++) {
+				expect(s0[i]).toBeCloseTo(1)
+				expect(s1[i]).toBeCloseTo(1)
+			}
+
+			const dad = Matrix.diag(d1).dot(a).dot(Matrix.diag(d2))
+			for (let i = 0; i < n; i++) {
+				for (let j = 0; j < n; j++) {
+					expect(dad.at(i, j)).toBeCloseTo(mat.at(i, j))
+				}
+			}
+		})
+
+		test('fail neg value', () => {
+			const mat = Matrix.randn(3, 3)
+			mat.set(0, 0, -0.1)
+			expect(() => mat.balancing()).toThrow('Doubly stochastic matrix only calculate for non negative matrix.')
+		})
+
+		test.each([
+			[2, 3],
+			[3, 2],
+		])('fail(%i, %i)', (r, c) => {
+			const mat = Matrix.randn(r, c)
+			expect(() => mat.balancing()).toThrow('Doubly stochastic matrix only defined for square matrix.')
+		})
+	})
+
+	describe('balancingSinkhornKnopp', () => {
+		test.each([1, 2, 3, 4])('square %d', n => {
+			const mat = Matrix.random(n, n, 0.1, 10)
+			const [d1, a, d2] = mat.balancingSinkhornKnopp()
+
+			expect(a.sizes).toEqual([n, n])
+			expect(d1).toHaveLength(n)
+			expect(d2).toHaveLength(n)
+			const s0 = a.sum(0).value
+			const s1 = a.sum(1).value
+			for (let i = 0; i < n; i++) {
+				expect(s0[i]).toBeCloseTo(1)
+				expect(s1[i]).toBeCloseTo(1)
+			}
+
+			const dad = Matrix.diag(d1).dot(a).dot(Matrix.diag(d2))
+			for (let i = 0; i < n; i++) {
+				for (let j = 0; j < n; j++) {
+					expect(dad.at(i, j)).toBeCloseTo(mat.at(i, j))
+				}
+			}
+		})
+
+		test('fail neg value', () => {
+			const mat = Matrix.randn(3, 3)
+			mat.set(0, 0, -0.1)
+			expect(() => mat.balancingSinkhornKnopp()).toThrow('Doubly stochastic matrix only calculate for non negative matrix.')
+		})
+
+		test.each([
+			[2, 3],
+			[3, 2],
+		])('fail(%i, %i)', (r, c) => {
+			const mat = Matrix.randn(r, c)
+			expect(() => mat.balancingSinkhornKnopp()).toThrow('Doubly stochastic matrix only defined for square matrix.')
+		})
+	})
+
 	describe('lu', () => {
 		test.each([0, 1, 2, 3, 5])('success %i', n => {
 			const mat = Matrix.randn(n, n)
@@ -6577,7 +6653,7 @@ describe('Matrix', () => {
 			expect(() => mat.eigenValuesQR()).toThrow('Eigen values only define square matrix.')
 		})
 
-		test('iteration not converged', () => {
+		test.only('iteration not converged', () => {
 			const mat = new Matrix(3, 3, [
 				[-0.3, -0.4, 1.7],
 				[-0.2, -1.8, -0.8],