Refactor singular value computation in Matrix class b to return both the maximum singular value and vector. (#1041)

Author: ishii-norimi

lib/util/matrix.js

@@ -4147,35 +4147,6 @@ 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>, number[], Matrix<number>]} Unitary matrix and singular values
@@ -4288,6 +4259,35 @@ export default class Matrix {
 		}
 	}
 
+	/**
+	 * Returns the maximum singular value and vector.
+	 * @param {number} [maxIteration] Maximum iteration
+	 * @returns {[Matrix<number>, number, Matrix<number>]} Maximum singular value and vector
+	 */
+	svdPowerIteration(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).toScaler()
+			v2.div(v2.norm())
+			u = this.dot(v2)
+			u.div(u.norm())
+			const e = Math.abs(s - ps)
+			if (e < tol || isNaN(e)) {
+				return [u, s, v2]
+			}
+			v = v2
+			ps = s
+		}
+		throw new MatrixException('svdPowerIteration not converged.', this)
+	}
+
 	/**
 	 * Returns a cholesky decomposition.
 	 * @returns {Matrix<number>} Cholesky decomposition matrix

tests/lib/util/matrix.test.js

@@ -1,3 +1,4 @@
+import { expect } from 'vitest'
 import Matrix from '../../../lib/util/matrix.js'
 import Tensor from '../../../lib/util/tensor.js'
 
@@ -6252,28 +6253,6 @@ 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],
@@ -6503,6 +6482,39 @@ describe('Matrix', () => {
 		})
 	})
 
+	describe('svdPowerIteration', () => {
+		test.each([
+			[10, 6],
+			[6, 10],
+		])('size [%i, %i]', (r, c) => {
+			const mat = Matrix.randn(r, c)
+			const [u, sv, v] = mat.svdPowerIteration()
+			expect(sv).toBeGreaterThanOrEqual(0)
+			expect(u.norm()).toBeCloseTo(1)
+			expect(v.norm()).toBeCloseTo(1)
+
+			const mv = mat.dot(v)
+			for (let i = 0; i < r; i++) {
+				expect(mv.at(i, 0)).toBeCloseTo(sv * u.at(i, 0))
+			}
+			const mu = mat.tDot(u)
+			for (let i = 0; i < c; i++) {
+				expect(mu.at(i, 0)).toBeCloseTo(sv * v.at(i, 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.svdPowerIteration(1)).toThrow('svdPowerIteration not converged.')
+		})
+	})
+
 	describe('cholesky', () => {
 		test('success', () => {
 			const n = 10