Add some singular / eigen value calculate method (#904)

Author: ishii-norimi

lib/util/matrix.js

@@ -3614,6 +3614,35 @@ export default class Matrix {
 		return ev
 	}
 
+	/**
+	 * Returns the maximum singular value.
+	 * @param {number} [maxIteration] Maximum iteration
+	 * @returns {number} Maximum singular value
+	 */
+	singularValuePowerIteration(maxIteration = 1.0e4) {
+		const tol = 1.0e-15
+		let v = Matrix.randn(this.cols, 1)
+		let u = Matrix.randn(this.rows, 1)
+		v.div(v.norm())
+		u.div(u.norm())
+		let ps = Infinity
+		let maxCount = maxIteration
+		while (maxCount-- > 0) {
+			const v2 = this.tDot(u)
+			const s = v.tDot(v2).at(0, 0)
+			v2.div(v2.norm())
+			u = this.dot(v2)
+			u.div(u.norm())
+			const e = Math.abs(s - ps)
+			if (e < tol || isNaN(e)) {
+				return s
+			}
+			v = v2
+			ps = s
+		}
+		throw new MatrixException('singularValuePowerIteration not converged.', this)
+	}
+
 	/**
 	 * Returns a singular value decomposition.
 	 * @returns {[Matrix, number[], Matrix]} Unitary matrix and singular values
@@ -4390,6 +4419,41 @@ export default class Matrix {
 		throw new MatrixException('eigenPowerIteration not converged.', this)
 	}
 
+	/**
+	 * Returns the k highest eigenvalues and its eigenvectors.
+	 * @param {number} k Number of calculated eigenvalues
+	 * @param {number} [maxIteration] Maximum iteration
+	 * @returns {[number[], Matrix]} The k highest eigenvalues and its eigenvectors. Eigenvectors correspond to each column of the matrix.
+	 */
+	eigenSimultaneousIteration(k, maxIteration = 1.0e4) {
+		if (!this.isSymmetric(1.0e-15)) {
+			throw new MatrixException('Simultaneous iteration method can only use symmetric matrix.', this)
+		}
+
+		const n = this.rows
+		const tol = 1.0e-8
+		let px = Matrix.ones(n, k)
+		px.div(px.norm())
+		let pl = new Matrix(1, k, Infinity)
+		const r0 = this.row(0)
+		let maxCount = maxIteration
+		while (maxCount-- > 0) {
+			const x = this.dot(px)
+			const [q] = x.qr()
+			q.resize(q.rows, k)
+			let ev = r0.dot(q)
+			ev.div(q.row(0))
+			pl.sub(ev)
+			const e = pl.reduce((s, v) => s + Math.abs(v), 0) / k
+			if (e < tol || isNaN(e)) {
+				return [ev.value, q]
+			}
+			px = q
+			pl = ev
+		}
+		throw new MatrixException('eigenSimultaneousIteration not converged.', this)
+	}
+
 	/**
 	 * Returns the nearest eigenvalue and its eigenvector to the specified value.
 	 * @param {number} [ev] Target value

tests/lib/util/matrix.test.js

@@ -5923,6 +5923,28 @@ describe('Matrix', () => {
 		})
 	})
 
+	describe('singularValuePowerIteration', () => {
+		test.each([
+			[10, 6],
+			[6, 10],
+		])('size [%i, %i]', (r, c) => {
+			const mat = Matrix.randn(r, c)
+			const sv = mat.singularValuePowerIteration()
+			expect(sv).toBeGreaterThanOrEqual(0)
+
+			const [ev] = mat.tDot(mat).eigenPowerIteration()
+			expect(sv ** 2).toBeCloseTo(ev)
+		})
+
+		test('iteration not converged', () => {
+			const mat = new Matrix(2, 2, [
+				[-1, -2],
+				[2, -2],
+			])
+			expect(() => mat.singularValuePowerIteration(1)).toThrow('singularValuePowerIteration not converged.')
+		})
+	})
+
 	describe('svd', () => {
 		test.each([
 			[10, 10],
@@ -7018,6 +7040,67 @@ describe('Matrix', () => {
 		})
 	})
 
+	describe('eigenSimultaneousIteration', () => {
+		test.each([
+			[1, 1],
+			[2, 2],
+			[5, 3],
+		])('symmetric %i,%i', (n, k) => {
+			const mat = Matrix.randn(n, n).gram()
+			const [eigvalues, eigvectors] = mat.eigenSimultaneousIteration(k)
+
+			for (let t = 0; t < k; t++) {
+				const cmat = mat.copy()
+				for (let i = 0; i < n; i++) {
+					cmat.subAt(i, i, eigvalues[t])
+				}
+				expect(cmat.det()).toBeCloseTo(0)
+			}
+
+			const x = mat.dot(eigvectors)
+			const y = Matrix.mult(eigvectors, new Matrix(1, k, eigvalues))
+			for (let t = 0; t < k; t++) {
+				for (let i = 0; i < n; i++) {
+					expect(x.at(i, t)).toBeCloseTo(y.at(i, t))
+				}
+			}
+			const eye = eigvectors.tDot(eigvectors)
+			for (let t = 0; t < k; t++) {
+				expect(eye.at(t, t)).toBeCloseTo(1)
+			}
+		})
+
+		test('non symmetric', () => {
+			const mat = new Matrix(4, 4, [
+				[16, -1, 1, 2],
+				[2, 12, 1, -1],
+				[1, 3, -24, 2],
+				[4, -2, 1, 20],
+			])
+			expect(() => mat.eigenSimultaneousIteration(1)).toThrow(
+				'Simultaneous iteration method can only use symmetric matrix.'
+			)
+		})
+
+		test.each([
+			[2, 3],
+			[3, 2],
+		])('fail(%i, %i)', (r, c) => {
+			const mat = Matrix.randn(r, c)
+			expect(() => mat.eigenSimultaneousIteration(1)).toThrow(
+				'Simultaneous iteration method can only use symmetric matrix.'
+			)
+		})
+
+		test('iteration not converged', () => {
+			const mat = new Matrix(2, 2, [
+				[-1, 2],
+				[2, -2],
+			])
+			expect(() => mat.eigenSimultaneousIteration(1, 1)).toThrow('eigenSimultaneousIteration not converged.')
+		})
+	})
+
 	describe('eigenInverseIteration', () => {
 		test.each([1, 2, 5])('symmetric %i', n => {
 			const mat = Matrix.randn(n, n).gram()