diff --git a/docs/index.html b/docs/index.html index 4b788124..2ddbf557 100644 --- a/docs/index.html +++ b/docs/index.html @@ -843,6 +843,40 @@

ml-matrix

+ + + + + +
  • + NNMF + ▸ + + +
    • @@ -2014,7 +2048,7 @@

    ▸ - randInt(rows, columns, maxValue, rng) + randInt(rows, columns, maxValue, minValue, rng)
    +
    +
    + minValue (number + = 0) + Minimum value + +
    + +
    +
    rng (function diff --git a/examples/linearlyDependent.js b/examples/linearlyDependent.js new file mode 100644 index 00000000..46359dad --- /dev/null +++ b/examples/linearlyDependent.js @@ -0,0 +1,11 @@ +import Matrix from '../src'; + +const test = [ + [1, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0, 0, 0, 0, 0], + [1, 1, 0, 0, 0, 0, 0, 0, 0, 0], + [1, 1, 1, 0, 0, 0, 0, 0, 0, 0], + [1, 1, 1, 1, 0, 0, 0, 0, 0, 0], + [1, 1, 1, 1, 1, 0, 0, 0, 0, 0], + [5, 4, 3, 2, 1, 0, 0, 0, 0, 0] +]; diff --git a/examples/nnmf.js b/examples/nnmf.js new file mode 100644 index 00000000..c865cb84 --- /dev/null +++ b/examples/nnmf.js @@ -0,0 +1,21 @@ +import Matrix from '../src'; + +const test = [ + [1, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [1, 1, 0, 0, 0, 0, 0, 0, 0, 0], + [1, 1, 1, 0, 0, 0, 0, 0, 0, 0], + [1, 1, 1, 1, 0, 0, 0, 0, 0, 0], + [1, 1, 1, 1, 1, 0, 0, 0, 0, 0], + [1, 1, 1, 1, 1, 1, 0, 0, 0, 0], + [1, 1, 1, 1, 1, 1, 1, 0, 0, 0], + [1, 1, 1, 1, 1, 1, 1, 1, 0, 0], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1] +]; +const targets = [ + [1, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 1, 1, 0, 0, 0, 0, 0, 0, 0], + [10, 9, 8, 7, 6, 5, 4, 3, 2, 1], + [4, 4, 2, 2, 0, 0, 0, 0, 0, 0] +]; diff --git a/matrix.d.ts b/matrix.d.ts index 65ca5659..2fff963e 100644 --- a/matrix.d.ts +++ b/matrix.d.ts @@ -331,6 +331,10 @@ declare module 'ml-matrix' { readonly diagonal: number[]; readonly diagonalMatrix: Matrix; } + class NNMF{ + constructor(value: MaybeMatrix, options?: number); + error(): number + } export interface ISVDOptions { computeLeftSingularVectors?: boolean; computeRightSingularVectors?: boolean; diff --git a/package-lock.json b/package-lock.json index 6bbd6e6f..b75eec19 100644 --- a/package-lock.json +++ b/package-lock.json @@ -1700,6 +1700,12 @@ "integrity": "sha512-qzm/XxIbxm/FHyH341ZrbnMUpe+5Bocte9xkmFMzPMjRaZMcXww+MpBptFvtU+79L362nqiLhekCxCxDPaUMBQ==", "dev": true }, + "esm": { + "version": "3.0.84", + "resolved": "https://registry.npmjs.org/esm/-/esm-3.0.84.tgz", + "integrity": "sha512-SzSGoZc17S7P+12R9cg21Bdb7eybX25RnIeRZ80xZs+VZ3kdQKzqTp2k4hZJjR7p9l0186TTXSgrxzlMDBktlw==", + "dev": true + }, "espree": { "version": "5.0.0", "resolved": "https://registry.npmjs.org/espree/-/espree-5.0.0.tgz", diff --git a/src/__tests__/decompositions/nnmf.js b/src/__tests__/decompositions/nnmf.js new file mode 100644 index 00000000..babfd0c7 --- /dev/null +++ b/src/__tests__/decompositions/nnmf.js @@ -0,0 +1,111 @@ +import { toBeDeepCloseTo } from 'jest-matcher-deep-close-to'; + +import { Matrix, NNMF } from '../..'; + +expect.extend({ toBeDeepCloseTo }); + +function positivity(An) { + let positive = true; + for (let i = 0; i < An.m; i++) { + for (let j = 0; j < An.r; j++) { + if (An.X.get(i, j) < 0) { + positive = false; + } + } + } + for (let i = 0; i < An.r; i++) { + for (let j = 0; j < An.n; j++) { + if (An.X.get(i, j) < 0) { + positive = false; + } + } + } + return (positive); +} + +describe('Non-negative Matrix Factorization', () => { + it('Factorization test I version 1', () => { + /** + * Global minimum : + * X = [ + * [2, 4], + * [8, 16], + * [32, 64], + * [128, 256]] + * + * Y = [ + * [1, 2, 3, 4], + * [6, 7, 8, 9]] + */ + let A = new Matrix([ + [26, 32, 38, 44], + [104, 128, 152, 176], + [416, 512, 608, 704], + [1664, 2048, 2432, 2816] + ]); + + let nA = + new NNMF(A, 1, { targetRelativeError: 0.000001, maxIterations: 100, version: 1 }); + + expect(positivity(nA)).toStrictEqual(true); + expect(nA.error.max()).toBeLessThan(0.000001); + }); + it('Factorization test I version 2', () => { + /** + * Global minimum : + * X = [ + * [2, 4], + * [8, 16], + * [32, 64], + * [128, 256]] + * + * Y = [ + * [1, 2, 3, 4], + * [6, 7, 8, 9]] + */ + let A = new Matrix([ + [26, 32, 38, 44], + [104, 128, 152, 176], + [416, 512, 608, 704], + [1664, 2048, 2432, 2816] + ]); + + let nA = + new NNMF(A, 1, { targetRelativeError: 0.000001, maxIterations: 100 }); + + expect(positivity(nA)).toStrictEqual(true); + expect(nA.error.max()).toBeLessThan(0.000001); + }); + it('Factorization test II', () => { + let A = new Matrix([ + [1, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0, 0, 0, 0, 0], + [1, 1, 0, 0, 0, 0, 0, 0, 0, 0], + [1, 1, 1, 0, 0, 0, 0, 0, 0, 0], + [1, 1, 1, 1, 0, 0, 0, 0, 0, 0], + [1, 1, 1, 1, 1, 0, 0, 0, 0, 0], + [5, 4, 3, 2, 1, 0, 0, 0, 0, 0] + ]); + + let nA = new NNMF(A, 3, { targetRelativeError: 1, maxIterations: 100 }); + + expect(positivity(nA)).toStrictEqual(true); + expect(nA.error.max()).toBeLessThan(1); + }); + it('Comparation of the error with matlab', () => { + // Working with 10*10 matrix and 100*100 matrix + // Matlab as a max error of 0.8 for 10*10 matrix and 0.98 for 100*100 matrix + // However, we don't have the same factorization + // But we have the same position of zeros in X and Y matrix + let A = Matrix.zeros(10, 10); + for (let i = 0; i < 10; i++) { + for (let j = 0; j <= i; j++) { + A.set(i, j, 1); + } + } + let nA = new NNMF(A, 1, { targetRelativeError: 1, maxIterations: 10 }); + + expect(positivity(nA)).toStrictEqual(true); + expect(nA.error.max()).toBeLessThan(1); + }); +}); diff --git a/src/abstractMatrix.js b/src/abstractMatrix.js index afaf3d9f..8147d4f3 100644 --- a/src/abstractMatrix.js +++ b/src/abstractMatrix.js @@ -132,16 +132,18 @@ export default function AbstractMatrix(superCtor) { * @param {number} rows - Number of rows * @param {number} columns - Number of columns * @param {number} [maxValue=1000] - Maximum value + * @param {number} [minValue=0] - Minimum value * @param {function} [rng=Math.random] - Random number generator * @return {Matrix} The new matrix */ - static randInt(rows, columns, maxValue, rng) { + static randInt(rows, columns, maxValue, minValue, rng) { if (maxValue === undefined) maxValue = 1000; + if (minValue === undefined) minValue = 0; if (rng === undefined) rng = Math.random; var matrix = this.empty(rows, columns); for (var i = 0; i < rows; i++) { for (var j = 0; j < columns; j++) { - var value = Math.floor(rng() * maxValue); + var value = Math.floor(rng() * (maxValue - minValue) + minValue); matrix.set(i, j, value); } } diff --git a/src/dc/nnmf.js b/src/dc/nnmf.js new file mode 100644 index 00000000..41548290 --- /dev/null +++ b/src/dc/nnmf.js @@ -0,0 +1,185 @@ +import { Matrix, WrapperMatrix2D } from '../index'; + +/** + * @class NNMF + * @link + * http://papers.nips.cc/paper/1861-algorithms-for-non-negative-matrix-factorization.pdf + * @param {Matrix} A + * @param {number} r + * @param {object} [options={}] + * @param {number} [options.targetRelativeError=0.001] //Maybe better with derivated + * @param {number} [options.maxIterations=10000] + * @param {number} [options.version=2] + */ +export default class NNMF { + constructor(A, r, options = {}) { + const { targetRelativeError = 0.001, maxIterations = 1000, version = 2 } = + options; + A = WrapperMatrix2D.checkMatrix(A); + var m = A.rows; + var n = A.columns; + + this.targetRelativeError = targetRelativeError; + this.r = r; + this.A = A; + this.m = m; + this.n = n; + this.X = Matrix.rand(this.m, this.r); + this.Y = Matrix.rand(this.r, this.n); + + let condition = false; + let time = 0; + + + do { + if (version === 1) { + doNnmf1.call(this, maxIterations / 2); + } else { + doNnmf2.call(this, maxIterations / 2); + } + condition = false; + for (let i = 0; i < m; i++) { + for (let j = 0; j < n; j++) { + if (this.error.get(i, j) > targetRelativeError) { + condition = true; + } + } + } + time++; + } while (condition && (time < 2)); + } + + /** + * Compute the error + */ + get error() { + let error = new Matrix(this.m, this.n); + let A2 = this.X.mmul(this.Y); + for (let i = 0; i < this.m; i++) { + for (let j = 0; j < this.n; j++) { + if (this.A.get(i, j) === 0) { + error.set(i, j, Math.abs(A2.get(i, j) - this.A.get(i, j))); + } else { + error.set(i, j, Math.abs(A2.get(i, j) - this.A.get(i, j)) / + this.A.get(i, j)); + } + } + } + return (error); + } +} + +/** + * Do the NNMF of a matrix A into two matrix X and Y + * @param {number} [numberIterations=1000] + */ +function doNnmf1(numberIterations = 1000) { + let A2 = this.X.mmul(this.Y); + let numX = Matrix.empty(this.m, this.r); + let denumX = Matrix.empty(this.m, this.r); + + let numY = Matrix.empty(this.r, this.n); + let denumY = Matrix.empty(this.r, this.n); + + let temp1 = Matrix.empty(this.m, this.n); + let temp2 = Matrix.empty(this.n, this.m); + + for (let a = 0; a < numberIterations; a++) { + for (let i = 0; i < this.m; i++) { + for (let j = 0; j < this.n; j++) { + if (A2.get(i, j) === 0) { + temp1.set(i, j, 0); + temp2.set(i, j, 0); + } else { + temp1.set(i, j, Math.pow(A2.get(i, j), -2) * this.A.get(i, j)); + temp2.set(i, j, Math.pow(A2.get(i, j), -1)); + } + } + } + numX = temp1.mmul(this.Y.transpose()); + denumX = temp2.mmul(this.Y.transpose()); + for (let i = 0; i < this.m; i++) { + for (let j = 0; j < this.r; j++) { + numX.set(i, j, numX.get(i, j) + Number.EPSILON); + denumX.set(i, j, denumX.get(i, j) + Number.EPSILON); + } + } + for (let i = 0; i < this.m; i++) { + for (let j = 0; j < this.r; j++) { + this.X.set(i, j, + (this.X.get(i, j) * numX.get(i, j)) / denumX.get(i, j)); + } + } + A2 = this.X.mmul(this.Y); + for (let i = 0; i < this.m; i++) { + for (let j = 0; j < this.n; j++) { + A2.set(i, j, A2.get(i, j) + Number.EPSILON); + } + } + + for (let i = 0; i < this.m; i++) { + for (let j = 0; j < this.n; j++) { + temp1.set(i, j, Math.pow(A2.get(i, j), -2) * this.A.get(i, j)); + temp2.set(i, j, Math.pow(A2.get(i, j), -1)); + } + } + numY = this.X.transpose().mmul(temp1); + denumY = this.X.transpose().mmul(temp2); + for (let i = 0; i < this.r; i++) { + for (let j = 0; j < this.n; j++) { + numY.set(i, j, numY.get(i, j) + Number.EPSILON); + denumY.set(i, j, denumY.get(i, j) + Number.EPSILON); + } + } + for (let i = 0; i < this.r; i++) { + for (let j = 0; j < this.n; j++) { + this.Y.set(i, j, + (this.Y.get(i, j) * numY.get(i, j)) / denumY.get(i, j)); + } + } + A2 = this.X.mmul(this.Y); + for (let i = 0; i < this.m; i++) { + for (let j = 0; j < this.n; j++) { + A2.set(i, j, A2.get(i, j) + Number.EPSILON); + } + } + } +} + +function doNnmf2(numberIterations = 1000) { + let numX = Matrix.empty(this.m, this.r); + let denumX = Matrix.empty(this.m, this.r); + let numY = Matrix.empty(this.r, this.n); + let denumY = Matrix.empty(this.r, this.n); + + for (let a = 0; a < numberIterations; a++) { + numX = this.A.mmul(this.Y.transpose()); + denumX = (this.X.mmul(this.Y)).mmul(this.Y.transpose()); + for (let i = 0; i < this.m; i++) { + for (let j = 0; j < this.r; j++) { + numX.set(i, j, numX.get(i, j) + Number.EPSILON); + denumX.set(i, j, denumX.get(i, j) + Number.EPSILON); + } + } + for (let i = 0; i < this.m; i++) { + for (let j = 0; j < this.r; j++) { + this.X.set(i, j, + (this.X.get(i, j) * numX.get(i, j)) / denumX.get(i, j)); + } + } + numY = this.X.transpose().mmul(this.A); + denumY = (this.X.transpose().mmul(this.X)).mmul(this.Y); + for (let i = 0; i < this.r; i++) { + for (let j = 0; j < this.n; j++) { + numY.set(i, j, numY.get(i, j) + Number.EPSILON); + denumY.set(i, j, denumY.get(i, j) + Number.EPSILON); + } + } + for (let i = 0; i < this.r; i++) { + for (let j = 0; j < this.n; j++) { + this.Y.set(i, j, + (this.Y.get(i, j) * numY.get(i, j)) / denumY.get(i, j)); + } + } + } +} diff --git a/src/index.js b/src/index.js index 8ab4329a..a994423f 100644 --- a/src/index.js +++ b/src/index.js @@ -21,3 +21,4 @@ export { } from './dc/cholesky.js'; export { default as LuDecomposition, default as LU } from './dc/lu.js'; export { default as QrDecomposition, default as QR } from './dc/qr.js'; +export { default as NNMFDecomposition, default as NNMF } from './dc/nnmf.js';