From 2f2307fff71a16df12eaf835ebb4c56aad6314c8 Mon Sep 17 00:00:00 2001 From: viktorpergjoka Date: Mon, 4 Dec 2017 23:18:01 +0100 Subject: [PATCH 1/2] Added linalg with cholesky decomposition --- src/math/linalg.ts | 78 +++++++++++++++++++++++++++++++++++++ src/math/linalg_test.ts | 86 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 164 insertions(+) create mode 100644 src/math/linalg.ts create mode 100644 src/math/linalg_test.ts diff --git a/src/math/linalg.ts b/src/math/linalg.ts new file mode 100644 index 0000000000..ab40f25df8 --- /dev/null +++ b/src/math/linalg.ts @@ -0,0 +1,78 @@ +/** + * @license + * Copyright 2017 Google Inc. All Rights Reserved. + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * ============================================================================= + */ + +import {Array2D} from "./ndarray"; + +/** + * Checks if a matrix is a square matrix + * @param {Array2D} a Matrix + * @returns {boolean} + */ +export function isSquare(a: Array2D): boolean { + return a.shape[0] === a.shape[1]; +} + +/** + * Checks if a matrix is symmetric + * @param {Array2D} a Matrix + * @returns {boolean} + */ +export function isSymmetric(a: Array2D): boolean { + if (!isSquare(a)) { + return false; + } + const n: number = a.shape[0]; + for (let i = 0; i < n; i++) { + for (let j = 0; j < i; j++) { + if (a.get(i, j) !== a.get(j, i)) { + return false; + } + } + } + return true; +} + +/** + * Return the Cholesky decomposition + * @param {Array2D} a Matrix + * @returns {Array2D} l lower triangular matrix + */ +export function cholesky(a: Array2D): Array2D { + if (!isSquare(a)) { + throw new Error('Cholesky error: matrix is not square'); + } + if (!isSymmetric(a)) { + throw new Error('Cholesky error: matrix is not symmetric'); + } + const n: number = a.shape[0]; + const l: Array2D = Array2D.zeros(a.shape, a.dtype); + for (let i = 0; i < n; i++) { + for (let j = 0; j < (i + 1); j++) { + let sum = 0; + for (let k = 0; k < j; k++) { + sum += l.get(i, k) * l.get(j, k); + } + const value = i === j ? Math.sqrt(a.get(i, i) - sum) || 0 : + (1.0 / l.get(j, j) * (a.get(i, j) - sum)); + l.set(value, i, j); + } + if (l.get(i, i) <= 0) { + throw new Error('Cholesky error: matrix is not positive definite'); + } + } + return l; +} diff --git a/src/math/linalg_test.ts b/src/math/linalg_test.ts new file mode 100644 index 0000000000..a2d9bea512 --- /dev/null +++ b/src/math/linalg_test.ts @@ -0,0 +1,86 @@ +/** + * @license + * Copyright 2017 Google Inc. All Rights Reserved. + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * ============================================================================= + */ + +import * as test_util from "../test_util"; +import {Array2D} from "./ndarray"; +import * as linalg from "./linalg"; + +function squareMatrixEquals(a: Array2D, b: Array2D): boolean { + const n = a.shape[0]; + for (let i = 0; i < n; i++) { + for (let j = 0; j < n; j++) { + if (a.get(i, j) !== b.get(i, j)) { + return false; + } + } + } + return true; +} + +test_util.describeCustom('linalg', () => { + + it('Cholesky decompose', () => { + const t1 = Array2D.new([3, 3], + [1, 2, 3, 2, 5, 7, 3, 7, 26]); + const chol1: Array2D = linalg.cholesky(t1); + const exp1 = Array2D.new([3, 3], [1, 0, 0, 2, 1, 0, 3, 1, 4]); + expect(squareMatrixEquals(chol1, exp1)).toEqual(true); + + const t2 = Array2D.new([3, 3], [4, 12, -16, 12, 37, -43, -16, -43, 98]); + const chol2 = linalg.cholesky(t2); + const exp2 = Array2D.new([3, 3], [2, 0, 0, 6, 1, 0, -8, 5, 3]); + expect(squareMatrixEquals(chol2, exp2)).toEqual(true); + + const t3 = Array2D.new([3, 3], [25, 15, -5, 15, 18, 0, -5, 0, 11]); + const chol3 = linalg.cholesky(t3); + const exp3 = Array2D.new([3, 3], [5, 0, 0, 3, 3, 0, -1, 1, 3]); + expect(squareMatrixEquals(chol3, exp3)).toEqual(true); + + const t4 = Array2D.new([4, 4], [18, 22, 54, 42, 22, 70, 86, 62, + 54, 86, 174, 134, 42, 62, 134, 106]); + const chol4 = linalg.cholesky(t4); + const exp4 = Array2D.new([4, 4], + new Float32Array([4.24264, 0, 0, 0, 5.18545, 6.56591, + 0, 0, 12.72792, 3.04604, 1.64974, 0, 9.89949, 1.62455, + 1.84971, 1.39262])); + test_util.expectArraysClose(chol4.dataSync(), exp4.dataSync()); + t4.dispose(); + exp4.dispose(); + }); + + it('Cholesky throw error for non square matrix', () => { + const t = Array2D.new([2, 1], [1, 2]); + expect(() => linalg.cholesky(t)) + .toThrowError('Cholesky error: matrix is not square'); + }); + + it('Cholesky throw error for non symmetric matrix', () => { + const t = Array2D.new([3, 3], [2, 6, 0, 8, 3, -11, 1, -1, 4]); + expect(() => linalg.cholesky(t)) + .toThrowError('Cholesky error: matrix is not symmetric'); + }); + + it('Cholesky throw error for not positive definite matrix', () => { + const t1 = Array2D.new([2, 2], [1, 2, 2, 1]); + expect(() => linalg.cholesky(t1)) + .toThrowError('Cholesky error: matrix is not positive definite'); + + const t2 = Array2D.new([3, 3], [-3, 0, 0, 0, -2, 0, 0, 0, 1]); + expect(() => linalg.cholesky(t2)) + .toThrowError('Cholesky error: matrix is not positive definite'); + }); +}); \ No newline at end of file From 15f8e0eefc4ac1b05bc38dacc3f212aef49ecb0b Mon Sep 17 00:00:00 2001 From: viktorpergjoka Date: Thu, 7 Dec 2017 20:42:24 +0100 Subject: [PATCH 2/2] added missing disposes --- src/math/linalg_test.ts | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/math/linalg_test.ts b/src/math/linalg_test.ts index a2d9bea512..7db1248883 100644 --- a/src/math/linalg_test.ts +++ b/src/math/linalg_test.ts @@ -58,6 +58,13 @@ test_util.describeCustom('linalg', () => { 0, 0, 12.72792, 3.04604, 1.64974, 0, 9.89949, 1.62455, 1.84971, 1.39262])); test_util.expectArraysClose(chol4.dataSync(), exp4.dataSync()); + + t1.dispose(); + exp1.dispose(); + t2.dispose(); + exp2.dispose(); + t3.dispose(); + exp3.dispose(); t4.dispose(); exp4.dispose(); }); @@ -66,12 +73,14 @@ test_util.describeCustom('linalg', () => { const t = Array2D.new([2, 1], [1, 2]); expect(() => linalg.cholesky(t)) .toThrowError('Cholesky error: matrix is not square'); + t.dispose(); }); it('Cholesky throw error for non symmetric matrix', () => { const t = Array2D.new([3, 3], [2, 6, 0, 8, 3, -11, 1, -1, 4]); expect(() => linalg.cholesky(t)) .toThrowError('Cholesky error: matrix is not symmetric'); + t.dispose(); }); it('Cholesky throw error for not positive definite matrix', () => { @@ -82,5 +91,7 @@ test_util.describeCustom('linalg', () => { const t2 = Array2D.new([3, 3], [-3, 0, 0, 0, -2, 0, 0, 0, 1]); expect(() => linalg.cholesky(t2)) .toThrowError('Cholesky error: matrix is not positive definite'); + t1.dispose(); + t2.dispose(); }); }); \ No newline at end of file