From 96a301ef16fb6d7d5f74f752889cb0c4782bad6d Mon Sep 17 00:00:00 2001 From: edoh Date: Sun, 19 Aug 2018 21:09:31 +0200 Subject: [PATCH 1/7] solve linear equation method of Jordan-Gauss triangularisation of the matrix check if all diagonal values are not null get the matrix solution of the equation --- src/ops/linalg_solve.ts | 101 +++++++++++++++++++++++++++++++++++ src/ops/linalg_solve_test.ts | 42 +++++++++++++++ src/ops/ops.ts | 1 + 3 files changed, 144 insertions(+) create mode 100644 src/ops/linalg_solve.ts create mode 100644 src/ops/linalg_solve_test.ts diff --git a/src/ops/linalg_solve.ts b/src/ops/linalg_solve.ts new file mode 100644 index 0000000000..d9ce642e52 --- /dev/null +++ b/src/ops/linalg_solve.ts @@ -0,0 +1,101 @@ +/** + * @license + * Copyright 2018 Google LLC. 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. + * ============================================================================= + */ + +/** + * Linear algebra resolution. + */ + +import {split, tensor1d} from '..'; +import {ENV} from '../environment'; +import {Tensor} from '../tensor'; +import {convertToTensor} from '../tensor_util_env'; +import {TensorLike} from '../types'; +import {assert} from '../util'; + +import {stack, unstack} from './array_ops'; +import {op} from './operation'; + +function diagonalMul(x: Tensor): Tensor { + assert(x.rank === 2, 'Input is not of rank 2'); + const [r, c] = x.shape; + assert(r === c, 'Input is not a square matrix'); + let mul = x.slice([0, 0], [1, 1]).as1D().asScalar(); + for (let i = 0; i < r; i++) { + mul = x.slice([i, i], [1, 1]).as1D().asScalar(); + } + return mul; +} + +function solve_(a: T|TensorLike, b: T|TensorLike): Tensor { + const $a = convertToTensor(a, 'a', 'solve'); + const $b = convertToTensor(b, 'b', 'solve'); + const [r, c] = $a.shape; + const [r2, c2] = $b.shape; + assert(r === r2, 'Second dimension size does not match'); + return ENV.engine.tidy(() => { + let inv: Tensor = $a.concat($b, 1); + const rows = Array.from({length: r}, (v, i) => i); + for (let i = 0; i < r; i++) { + inv = ENV.engine.tidy(() => { + for (let j = i + 1; j < r; j++) { + const elt = inv.slice([j, i], [1, 1]).as1D().asScalar(); + const pivot = inv.slice([i, i], [1, 1]).as1D().asScalar(); + if (elt.dataSync()[0] !== 0) { + const newrow = + inv.gather(tensor1d([i], 'int32')) + .sub(inv.gather(tensor1d([j], 'int32')).div(elt).mul(pivot)) + .as1D(); + const sli = + inv.gather(tensor1d(rows.filter(e => e !== j), 'int32')); + const arr: Tensor[] = []; + if (j === 0) { + arr.push(newrow); + } + unstack(sli).forEach((t, ind) => { + if (ind !== j) { + arr.push(t); + } else { + arr.push(newrow); + arr.push(t); + } + }); + if (j === r - 1) { + arr.push(newrow); + } + inv = stack(arr); + } + } + return inv; // inv is an upper triangular matrix + }); + } + const determinant = diagonalMul(split(inv, [c, c2], 1)[0]).dataSync(); + assert(determinant[0] !== 0, 'Input matrix is not inversible'); + const trian = unstack(inv); + const len = trian.length; + trian[len - 1] = + trian[len - 1].div(trian[len - 1].slice(r - 1, 1).asScalar()); + for (let i = r - 2; i > -1; i--) { + for (let j = r - 1; j > i; j--) { + trian[i] = trian[i].sub(trian[j].mul(trian[i].slice(j, 1).asScalar())); + } + trian[i] = trian[i].div(trian[i].slice(i, 1).asScalar()); + } + return split(stack(trian), [c, c2], 1)[1]; + }); +} + +export const solve = op({solve_}); diff --git a/src/ops/linalg_solve_test.ts b/src/ops/linalg_solve_test.ts new file mode 100644 index 0000000000..a64042cd35 --- /dev/null +++ b/src/ops/linalg_solve_test.ts @@ -0,0 +1,42 @@ +/** + * @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 tf from '../index'; +import {describeWithFlags} from '../jasmine_util'; +import {ALL_ENVS, expectArraysClose} from '../test_util'; + +describeWithFlags('solve_linear', ALL_ENVS, () => { + it('solve equation with a the matrix eye', () => { + const a = tf.eye(3); + const b = tf.tensor2d([1, 2, 3], [3, 1]); + const x = tf.solve(a, b); + expect(() => expectArraysClose(x, [1, 2, 3])); + }); + + it('solve equation with a the matrix eye times 2', () => { + const a = tf.eye(3).mul(tf.scalar(2)); + const b = tf.tensor2d([1, 2, 3], [3, 1]); + const x = tf.solve(a, b); + expect(() => expectArraysClose(x, [1, 2, 3])); + }); + + it('should throw error if a is not inversible', () => { + const a = tf.ones([3, 3]); + const b = tf.tensor2d([1, 2, 3], [3, 1]); + expect(() => tf.solve(a, b)).toThrowError('Input matrix is not inversible'); + }); +}); diff --git a/src/ops/ops.ts b/src/ops/ops.ts index 2c9be6f47f..3f88d4f812 100644 --- a/src/ops/ops.ts +++ b/src/ops/ops.ts @@ -39,6 +39,7 @@ export * from './lstm'; export * from './moving_average'; export * from './strided_slice'; export * from './topk'; +export * from './linalg_solve'; export {op} from './operation'; From ad328364fb2fca0182565bbf21045897ec84a6cc Mon Sep 17 00:00:00 2001 From: edoh Date: Tue, 21 Aug 2018 00:17:31 +0200 Subject: [PATCH 2/7] invert matrix --- src/ops/linalg_solve.ts | 108 ++++++++++++++++++++++------------- src/ops/linalg_solve_test.ts | 11 ++++ 2 files changed, 78 insertions(+), 41 deletions(-) diff --git a/src/ops/linalg_solve.ts b/src/ops/linalg_solve.ts index d9ce642e52..26c10a1f4f 100644 --- a/src/ops/linalg_solve.ts +++ b/src/ops/linalg_solve.ts @@ -19,16 +19,65 @@ * Linear algebra resolution. */ -import {split, tensor1d} from '..'; +import {scalar, split, tensor1d} from '..'; import {ENV} from '../environment'; -import {Tensor} from '../tensor'; +import {Scalar, Tensor} from '../tensor'; import {convertToTensor} from '../tensor_util_env'; import {TensorLike} from '../types'; import {assert} from '../util'; -import {stack, unstack} from './array_ops'; +import {eye, stack, unstack} from './array_ops'; import {op} from './operation'; +function gaussJordanTriangular( + $a: Tensor, $b: Tensor): {upperM: Tensor, det: Scalar} { + console.log('gauss'); + const [r, c] = $a.shape; + const [r2, c2] = $b.shape; + assert(r === r2, 'Second dimension size does not match'); + let inv: Tensor = $a.concat($b, 1); + const rows = Array.from({length: r}, (v, i) => i); + let coef = scalar(1); + for (let i = 0; i < r; i++) { + ({inv, coef} = ENV.engine.tidy(() => { + for (let j = i + 1; j < r; j++) { + const elt = inv.slice([j, i], [1, 1]).as1D().asScalar(); + const pivot = inv.slice([i, i], [1, 1]).as1D().asScalar(); + if (elt.dataSync()[0] !== 0) { + const factor = pivot.div(elt); + coef = coef.mul(factor).mul(scalar(-1)); + const newrow = + inv.gather(tensor1d([i], 'int32')) + .sub(inv.gather(tensor1d([j], 'int32')).mul(factor)) + .as1D(); + const sli = inv.gather(tensor1d(rows.filter(e => e !== j), 'int32')); + const arr: Tensor[] = []; + if (j === 0) { + arr.push(newrow); + } + unstack(sli).forEach((t, ind) => { + if (ind !== j) { + arr.push(t); + } else { + arr.push(newrow); + arr.push(t); + } + }); + if (j === r - 1) { + arr.push(newrow); + } + inv = stack(arr); + } + } + // the first c colomns of inv is an upper triangular matrix + return {inv, coef}; + })); + } + const determinant = + diagonalMul(split(inv, [c, c2], 1)[0]).mul(coef).asScalar(); + return {upperM: inv, det: determinant}; +} + function diagonalMul(x: Tensor): Tensor { assert(x.rank === 2, 'Input is not of rank 2'); const [r, c] = x.shape; @@ -47,44 +96,10 @@ function solve_(a: T|TensorLike, b: T|TensorLike): Tensor { const [r2, c2] = $b.shape; assert(r === r2, 'Second dimension size does not match'); return ENV.engine.tidy(() => { - let inv: Tensor = $a.concat($b, 1); - const rows = Array.from({length: r}, (v, i) => i); - for (let i = 0; i < r; i++) { - inv = ENV.engine.tidy(() => { - for (let j = i + 1; j < r; j++) { - const elt = inv.slice([j, i], [1, 1]).as1D().asScalar(); - const pivot = inv.slice([i, i], [1, 1]).as1D().asScalar(); - if (elt.dataSync()[0] !== 0) { - const newrow = - inv.gather(tensor1d([i], 'int32')) - .sub(inv.gather(tensor1d([j], 'int32')).div(elt).mul(pivot)) - .as1D(); - const sli = - inv.gather(tensor1d(rows.filter(e => e !== j), 'int32')); - const arr: Tensor[] = []; - if (j === 0) { - arr.push(newrow); - } - unstack(sli).forEach((t, ind) => { - if (ind !== j) { - arr.push(t); - } else { - arr.push(newrow); - arr.push(t); - } - }); - if (j === r - 1) { - arr.push(newrow); - } - inv = stack(arr); - } - } - return inv; // inv is an upper triangular matrix - }); - } - const determinant = diagonalMul(split(inv, [c, c2], 1)[0]).dataSync(); - assert(determinant[0] !== 0, 'Input matrix is not inversible'); - const trian = unstack(inv); + const {upperM, det} = gaussJordanTriangular($a, $b); + console.log('determinant', det); + assert(det.dataSync()[0] !== 0, 'Input matrix is not inversible'); + const trian = unstack(upperM); const len = trian.length; trian[len - 1] = trian[len - 1].div(trian[len - 1].slice(r - 1, 1).asScalar()); @@ -98,4 +113,15 @@ function solve_(a: T|TensorLike, b: T|TensorLike): Tensor { }); } +function invertMatrix_(x: T): T { + const $x = convertToTensor(x, 'x', 'invertMatrix'); + assert($x.rank === 2, 'Input is not of rank 2'); + const [r, c] = $x.shape; + assert(r === c, 'Input is not a square matrix'); + return solve($x, eye(r) as T) as T; + // assert(det.dataSync()[0] !== 0, 'Input matrix is not inversible'); + // return split(upperM, 2, 1)[1] as T; +} + export const solve = op({solve_}); +export const invertMatrix = op({invertMatrix_}); diff --git a/src/ops/linalg_solve_test.ts b/src/ops/linalg_solve_test.ts index a64042cd35..6cab2e25a0 100644 --- a/src/ops/linalg_solve_test.ts +++ b/src/ops/linalg_solve_test.ts @@ -18,6 +18,7 @@ import * as tf from '../index'; import {describeWithFlags} from '../jasmine_util'; import {ALL_ENVS, expectArraysClose} from '../test_util'; +import {invertMatrix} from './linalg_solve'; describeWithFlags('solve_linear', ALL_ENVS, () => { it('solve equation with a the matrix eye', () => { @@ -40,3 +41,13 @@ describeWithFlags('solve_linear', ALL_ENVS, () => { expect(() => tf.solve(a, b)).toThrowError('Input matrix is not inversible'); }); }); + +describeWithFlags('invert_matrix', ALL_ENVS, () => { + it('solve equation with a the matrix eye', () => { + const m = tf.tensor2d([1, 3, 3, 1, 4, 3, 1, 3, 4], [3, 3]); + const inv = invertMatrix(m); + m.print(); + inv.print(); + expect(() => expectArraysClose(inv, [7, -3, -3, -1, 1, 0, -1, 0, 1])); + }); +}); From a2bcda39e236a6f9e12dc62d7409f525e83f2bfc Mon Sep 17 00:00:00 2001 From: edoh Date: Tue, 21 Aug 2018 20:08:05 +0200 Subject: [PATCH 3/7] compute matrix adjoint --- src/ops/linalg_solve.ts | 84 +++++++++++++++++++++++++----------- src/ops/linalg_solve_test.ts | 20 +++++++-- 2 files changed, 75 insertions(+), 29 deletions(-) diff --git a/src/ops/linalg_solve.ts b/src/ops/linalg_solve.ts index 26c10a1f4f..aec433c8fa 100644 --- a/src/ops/linalg_solve.ts +++ b/src/ops/linalg_solve.ts @@ -21,21 +21,19 @@ import {scalar, split, tensor1d} from '..'; import {ENV} from '../environment'; -import {Scalar, Tensor} from '../tensor'; -import {convertToTensor} from '../tensor_util_env'; -import {TensorLike} from '../types'; +import {Scalar, Tensor, Tensor2D} from '../tensor'; import {assert} from '../util'; import {eye, stack, unstack} from './array_ops'; import {op} from './operation'; function gaussJordanTriangular( - $a: Tensor, $b: Tensor): {upperM: Tensor, det: Scalar} { + a: Tensor2D, b: Tensor2D): {upperM: Tensor2D, det: Scalar} { console.log('gauss'); - const [r, c] = $a.shape; - const [r2, c2] = $b.shape; + const [r, c] = a.shape; + const [r2, c2] = b.shape; assert(r === r2, 'Second dimension size does not match'); - let inv: Tensor = $a.concat($b, 1); + let inv: Tensor = a.concat(b, 1); const rows = Array.from({length: r}, (v, i) => i); let coef = scalar(1); for (let i = 0; i < r; i++) { @@ -73,14 +71,17 @@ function gaussJordanTriangular( return {inv, coef}; })); } + console.log('coef', coef.dataSync()[0]); const determinant = - diagonalMul(split(inv, [c, c2], 1)[0]).mul(coef).asScalar(); - return {upperM: inv, det: determinant}; + diagonalMul(split(inv, [c, c2], 1)[0] as Tensor2D).div(coef).asScalar(); + console.log('determinant', determinant.dataSync()[0]); + return {upperM: inv as Tensor2D, det: determinant}; } -function diagonalMul(x: Tensor): Tensor { - assert(x.rank === 2, 'Input is not of rank 2'); +function diagonalMul(x: Tensor2D): Scalar { const [r, c] = x.shape; + console.log('diag'); + x.print(); assert(r === c, 'Input is not a square matrix'); let mul = x.slice([0, 0], [1, 1]).as1D().asScalar(); for (let i = 0; i < r; i++) { @@ -89,14 +90,12 @@ function diagonalMul(x: Tensor): Tensor { return mul; } -function solve_(a: T|TensorLike, b: T|TensorLike): Tensor { - const $a = convertToTensor(a, 'a', 'solve'); - const $b = convertToTensor(b, 'b', 'solve'); - const [r, c] = $a.shape; - const [r2, c2] = $b.shape; +function solve_(a: Tensor2D, b: Tensor2D): Tensor2D { + const [r, c] = a.shape; + const [r2, c2] = b.shape; assert(r === r2, 'Second dimension size does not match'); return ENV.engine.tidy(() => { - const {upperM, det} = gaussJordanTriangular($a, $b); + const {upperM, det} = gaussJordanTriangular(a, b); console.log('determinant', det); assert(det.dataSync()[0] !== 0, 'Input matrix is not inversible'); const trian = unstack(upperM); @@ -109,19 +108,54 @@ function solve_(a: T|TensorLike, b: T|TensorLike): Tensor { } trian[i] = trian[i].div(trian[i].slice(i, 1).asScalar()); } - return split(stack(trian), [c, c2], 1)[1]; + return split(stack(trian), [c, c2], 1)[1] as Tensor2D; }); } -function invertMatrix_(x: T): T { - const $x = convertToTensor(x, 'x', 'invertMatrix'); - assert($x.rank === 2, 'Input is not of rank 2'); - const [r, c] = $x.shape; +function invertMatrix_(x: Tensor2D): Tensor2D { + const [r, c] = x.shape; + assert(r === c, 'Input is not a square matrix'); + return solve(x, eye(r)); +} + +function det_(m: Tensor2D): Scalar { + return gaussJordanTriangular(m, eye(m.shape[0]) as Tensor2D).det; +} + +function adjointM_(m: Tensor2D): Tensor2D { + /* const [r, c] = m.shape; assert(r === c, 'Input is not a square matrix'); - return solve($x, eye(r) as T) as T; - // assert(det.dataSync()[0] !== 0, 'Input matrix is not inversible'); - // return split(upperM, 2, 1)[1] as T; + const rows = Array.from({length: r}, (v, i) => i); + const dets: Tensor1D[] = []; + for (let i = 0; i < r; i++) { + for (let j = 0; j < r; j++) { + const mSub = m.gather(tensor1d(rows.filter(e => e !== i), 'int32')); + mSub.print(); + let sli; + if (j === 0) { + sli = mSub.slice([0, 1], [r - 1, r - 1]); + } else if (j === r - 1) { + sli = mSub.slice([0, 0], [r - 1, r - 1]); + } else { + const [a, b, c] = split(mSub, [j, 1, r - (j + 1)], 1); + a.print(); + c.print(); + b.dispose(); + sli = concat([a, c], 1); + sli.print(); + } + dets.push( + scalar(Math.pow(-1, (i + j))) + .mul(gaussJordanTriangular(sli, eye(sli.shape[0]) as Tensor2D) + .det.as1D())); + } + } + const adjM = transpose(concat(dets).reshape([r, r])); + return adjM as Tensor2D; */ + return invertMatrix(m).mul(det(m)); } export const solve = op({solve_}); export const invertMatrix = op({invertMatrix_}); +export const adjM = op({adjointM_}); +export const det = op({det_}); diff --git a/src/ops/linalg_solve_test.ts b/src/ops/linalg_solve_test.ts index 6cab2e25a0..f3067abdfe 100644 --- a/src/ops/linalg_solve_test.ts +++ b/src/ops/linalg_solve_test.ts @@ -18,7 +18,8 @@ import * as tf from '../index'; import {describeWithFlags} from '../jasmine_util'; import {ALL_ENVS, expectArraysClose} from '../test_util'; -import {invertMatrix} from './linalg_solve'; + +import {adjM, invertMatrix} from './linalg_solve'; describeWithFlags('solve_linear', ALL_ENVS, () => { it('solve equation with a the matrix eye', () => { @@ -29,21 +30,21 @@ describeWithFlags('solve_linear', ALL_ENVS, () => { }); it('solve equation with a the matrix eye times 2', () => { - const a = tf.eye(3).mul(tf.scalar(2)); + const a = tf.eye(3).mul(tf.scalar(2)) as tf.Tensor2D; const b = tf.tensor2d([1, 2, 3], [3, 1]); const x = tf.solve(a, b); expect(() => expectArraysClose(x, [1, 2, 3])); }); it('should throw error if a is not inversible', () => { - const a = tf.ones([3, 3]); + const a = tf.ones([3, 3]) as tf.Tensor2D; const b = tf.tensor2d([1, 2, 3], [3, 1]); expect(() => tf.solve(a, b)).toThrowError('Input matrix is not inversible'); }); }); describeWithFlags('invert_matrix', ALL_ENVS, () => { - it('solve equation with a the matrix eye', () => { + it('invert a matrix', () => { const m = tf.tensor2d([1, 3, 3, 1, 4, 3, 1, 3, 4], [3, 3]); const inv = invertMatrix(m); m.print(); @@ -51,3 +52,14 @@ describeWithFlags('invert_matrix', ALL_ENVS, () => { expect(() => expectArraysClose(inv, [7, -3, -3, -1, 1, 0, -1, 0, 1])); }); }); + +describeWithFlags('adjoint_matrix', ALL_ENVS, () => { + it('computes the adjoint of a matrix', () => { + const m = tf.tensor2d([1, 3, 3, 1, 4, 3, 1, 3, 4], [3, 3]); + const a = adjM(m); + m.print(); + console.log('adjm'); + a.print(); + expect(() => expectArraysClose(a, [7, -3, -3, -1, 1, 0, -1, 0, 1])); + }); +}); From aa0381a1a9bbfe03b4e588ed9bc79f418116c31f Mon Sep 17 00:00:00 2001 From: edoh Date: Tue, 21 Aug 2018 21:05:24 +0200 Subject: [PATCH 4/7] clean code remove comments and console.log --- src/ops/linalg_solve.ts | 34 ---------------------------------- src/ops/linalg_solve_test.ts | 1 - 2 files changed, 35 deletions(-) diff --git a/src/ops/linalg_solve.ts b/src/ops/linalg_solve.ts index aec433c8fa..a742b1bc86 100644 --- a/src/ops/linalg_solve.ts +++ b/src/ops/linalg_solve.ts @@ -29,7 +29,6 @@ import {op} from './operation'; function gaussJordanTriangular( a: Tensor2D, b: Tensor2D): {upperM: Tensor2D, det: Scalar} { - console.log('gauss'); const [r, c] = a.shape; const [r2, c2] = b.shape; assert(r === r2, 'Second dimension size does not match'); @@ -71,16 +70,13 @@ function gaussJordanTriangular( return {inv, coef}; })); } - console.log('coef', coef.dataSync()[0]); const determinant = diagonalMul(split(inv, [c, c2], 1)[0] as Tensor2D).div(coef).asScalar(); - console.log('determinant', determinant.dataSync()[0]); return {upperM: inv as Tensor2D, det: determinant}; } function diagonalMul(x: Tensor2D): Scalar { const [r, c] = x.shape; - console.log('diag'); x.print(); assert(r === c, 'Input is not a square matrix'); let mul = x.slice([0, 0], [1, 1]).as1D().asScalar(); @@ -96,7 +92,6 @@ function solve_(a: Tensor2D, b: Tensor2D): Tensor2D { assert(r === r2, 'Second dimension size does not match'); return ENV.engine.tidy(() => { const {upperM, det} = gaussJordanTriangular(a, b); - console.log('determinant', det); assert(det.dataSync()[0] !== 0, 'Input matrix is not inversible'); const trian = unstack(upperM); const len = trian.length; @@ -123,35 +118,6 @@ function det_(m: Tensor2D): Scalar { } function adjointM_(m: Tensor2D): Tensor2D { - /* const [r, c] = m.shape; - assert(r === c, 'Input is not a square matrix'); - const rows = Array.from({length: r}, (v, i) => i); - const dets: Tensor1D[] = []; - for (let i = 0; i < r; i++) { - for (let j = 0; j < r; j++) { - const mSub = m.gather(tensor1d(rows.filter(e => e !== i), 'int32')); - mSub.print(); - let sli; - if (j === 0) { - sli = mSub.slice([0, 1], [r - 1, r - 1]); - } else if (j === r - 1) { - sli = mSub.slice([0, 0], [r - 1, r - 1]); - } else { - const [a, b, c] = split(mSub, [j, 1, r - (j + 1)], 1); - a.print(); - c.print(); - b.dispose(); - sli = concat([a, c], 1); - sli.print(); - } - dets.push( - scalar(Math.pow(-1, (i + j))) - .mul(gaussJordanTriangular(sli, eye(sli.shape[0]) as Tensor2D) - .det.as1D())); - } - } - const adjM = transpose(concat(dets).reshape([r, r])); - return adjM as Tensor2D; */ return invertMatrix(m).mul(det(m)); } diff --git a/src/ops/linalg_solve_test.ts b/src/ops/linalg_solve_test.ts index f3067abdfe..c24ef38f64 100644 --- a/src/ops/linalg_solve_test.ts +++ b/src/ops/linalg_solve_test.ts @@ -58,7 +58,6 @@ describeWithFlags('adjoint_matrix', ALL_ENVS, () => { const m = tf.tensor2d([1, 3, 3, 1, 4, 3, 1, 3, 4], [3, 3]); const a = adjM(m); m.print(); - console.log('adjm'); a.print(); expect(() => expectArraysClose(a, [7, -3, -3, -1, 1, 0, -1, 0, 1])); }); From 06d3cdf4dd46f01f5caefadcb94e895be258046f Mon Sep 17 00:00:00 2001 From: edoh Date: Tue, 21 Aug 2018 22:35:27 +0200 Subject: [PATCH 5/7] array input or not --- src/ops/linalg_solve.ts | 30 ++++++++++++++++++++++++++---- 1 file changed, 26 insertions(+), 4 deletions(-) diff --git a/src/ops/linalg_solve.ts b/src/ops/linalg_solve.ts index a742b1bc86..b6a3396b4f 100644 --- a/src/ops/linalg_solve.ts +++ b/src/ops/linalg_solve.ts @@ -86,11 +86,33 @@ function diagonalMul(x: Tensor2D): Scalar { return mul; } -function solve_(a: Tensor2D, b: Tensor2D): Tensor2D { - const [r, c] = a.shape; - const [r2, c2] = b.shape; - assert(r === r2, 'Second dimension size does not match'); +// function solve_(a: Tensor2D): Tensor; +function solve_( + a: Tensor2D[]|Tensor2D, b: Tensor2D[]|Tensor2D, adjoint = false): Tensor2D| + Tensor2D[] { + if (Array.isArray(a) || Array.isArray(b)) { + assert( + (a as Tensor2D[]).length === (b as Tensor2D[]).length, + 'Second dimension size does not match'); + const sol: Tensor2D[] = []; + (a as Tensor2D[]).forEach((m, i) => { + sol.push(solve_unique_equation(m, (b as Tensor2D[])[i], adjoint)); + }); + return sol; + } else { + return solve_unique_equation(a, b, adjoint); + } +} + +function solve_unique_equation( + a: Tensor2D, b: Tensor2D, adjoint = false): Tensor2D { return ENV.engine.tidy(() => { + const [r, c] = a.shape; + const [r2, c2] = b.shape; + assert(r === r2, 'Second dimension size does not match'); + if (adjoint) { + a = adjM(a); + } const {upperM, det} = gaussJordanTriangular(a, b); assert(det.dataSync()[0] !== 0, 'Input matrix is not inversible'); const trian = unstack(upperM); From 628d585e4c0c320619be65e63a180262b04b747e Mon Sep 17 00:00:00 2001 From: edoh Date: Wed, 22 Aug 2018 00:03:14 +0200 Subject: [PATCH 6/7] solve linear equation Using Jordan-Gauss elimination to solve linear equation AX = B. Compute the determinant of a matrix Compute the adjoint of a matrix FEATURE solve linear equation --- src/ops/linalg_solve.ts | 65 ++++++++++++++++++++++++++++++------ src/ops/linalg_solve_test.ts | 19 +++++++---- 2 files changed, 66 insertions(+), 18 deletions(-) diff --git a/src/ops/linalg_solve.ts b/src/ops/linalg_solve.ts index b6a3396b4f..3a8e413193 100644 --- a/src/ops/linalg_solve.ts +++ b/src/ops/linalg_solve.ts @@ -27,6 +27,17 @@ import {assert} from '../util'; import {eye, stack, unstack} from './array_ops'; import {op} from './operation'; +/** + * + * @param a is a square matrix M (Tensor2d with shape `[r, c]` such that `r === + * c`) + * @param b is a Tensor2d with shape `[r2, c2]` + * @desc `r === r2` + * @returns a matrix of shape `[r, c+c2]` after a [jordan-gauss + * elimination](https://en.wikipedia.org/wiki/Gaussian_elimination) on the + * matrix given by the concatenation `[a, b]`. The first r or c columns is an + * upper triangular matrix + */ function gaussJordanTriangular( a: Tensor2D, b: Tensor2D): {upperM: Tensor2D, det: Scalar} { const [r, c] = a.shape; @@ -75,21 +86,37 @@ function gaussJordanTriangular( return {upperM: inv as Tensor2D, det: determinant}; } -function diagonalMul(x: Tensor2D): Scalar { - const [r, c] = x.shape; - x.print(); +/** + * + * @param m Tensor2d or matrix + * @returns the product of the diagonal elements of @param m as a `tf.scalar` + */ +function diagonalMul(m: Tensor2D): Scalar { + const [r, c] = m.shape; assert(r === c, 'Input is not a square matrix'); - let mul = x.slice([0, 0], [1, 1]).as1D().asScalar(); + let mul = m.slice([0, 0], [1, 1]).as1D().asScalar(); for (let i = 0; i < r; i++) { - mul = x.slice([i, i], [1, 1]).as1D().asScalar(); + mul = m.slice([i, i], [1, 1]).as1D().asScalar(); } return mul; } -// function solve_(a: Tensor2D): Tensor; +/** + * + * @param a is a unique square matrix M or a tensor of shape `[..., M, M]` whose + * inner-most 2 dimensions form square matrices + * @param b is a unique matrix M or a tensor of shape `[..., M, M]` + * @param adjoint is a boolean. + * If adjoint is false then each output matrix satisfies `a * output = b` + * (respectively `a[..., :, :] * output[..., :, :] = b[..., :, :]` if the inputs + * are arrays of matrixes) . If adjoint is true then each output matrix + * satisfies `adjoint(a) * output = b` (respectively `adjoint(a[..., :, :]) * + * output[..., :, :] = b[..., :, + * :]`). + */ function solve_( - a: Tensor2D[]|Tensor2D, b: Tensor2D[]|Tensor2D, adjoint = false): Tensor2D| - Tensor2D[] { + a: Tensor2D[]|Tensor2D, b: Tensor2D[]|Tensor2D, + adjoint = false): Tensor2D[]|Tensor2D { if (Array.isArray(a) || Array.isArray(b)) { assert( (a as Tensor2D[]).length === (b as Tensor2D[]).length, @@ -104,6 +131,7 @@ function solve_( } } +// helper to the solve equation function solve_unique_equation( a: Tensor2D, b: Tensor2D, adjoint = false): Tensor2D { return ENV.engine.tidy(() => { @@ -129,16 +157,31 @@ function solve_unique_equation( }); } -function invertMatrix_(x: Tensor2D): Tensor2D { - const [r, c] = x.shape; +/** + * + * @param x square matrix to invert + * @returns the invert matrix of @param m if inversible + */ +function invertMatrix_(m: Tensor2D): Tensor2D { + const [r, c] = m.shape; assert(r === c, 'Input is not a square matrix'); - return solve(x, eye(r)); + return solve(m, eye(r)) as Tensor2D; } +/** + * + * @param m Tensor2d or matrix + * @returns the determinant of @param m as a `tf.scalar` + */ function det_(m: Tensor2D): Scalar { return gaussJordanTriangular(m, eye(m.shape[0]) as Tensor2D).det; } +/** + * + * @param m Tensor2d or matrix + * @returns the adjoint of @param m if inversible + */ function adjointM_(m: Tensor2D): Tensor2D { return invertMatrix(m).mul(det(m)); } diff --git a/src/ops/linalg_solve_test.ts b/src/ops/linalg_solve_test.ts index c24ef38f64..08875ce583 100644 --- a/src/ops/linalg_solve_test.ts +++ b/src/ops/linalg_solve_test.ts @@ -16,23 +16,32 @@ */ import * as tf from '../index'; +import {concat} from '../index'; import {describeWithFlags} from '../jasmine_util'; import {ALL_ENVS, expectArraysClose} from '../test_util'; import {adjM, invertMatrix} from './linalg_solve'; describeWithFlags('solve_linear', ALL_ENVS, () => { - it('solve equation with a the matrix eye', () => { + it('solve equation with a: the matrix eye', () => { const a = tf.eye(3); const b = tf.tensor2d([1, 2, 3], [3, 1]); - const x = tf.solve(a, b); + const x = tf.solve(a, b) as tf.Tensor2D; + expect(() => expectArraysClose(x, [1, 2, 3])); + }); + + it('solve equation with a: an array of matrix eye', () => { + const a = [tf.eye(3), tf.eye(3), tf.eye(3)]; + const b = Array.from({length: 3}, (v, i) => tf.tensor2d([1, 2, 3], [3, 1])); + const x = concat(tf.solve(a, b) as tf.Tensor2D[]); + (tf.solve(a, b) as tf.Tensor2D[]).forEach(e => e.print()); expect(() => expectArraysClose(x, [1, 2, 3])); }); it('solve equation with a the matrix eye times 2', () => { const a = tf.eye(3).mul(tf.scalar(2)) as tf.Tensor2D; const b = tf.tensor2d([1, 2, 3], [3, 1]); - const x = tf.solve(a, b); + const x = tf.solve(a, b) as tf.Tensor2D; expect(() => expectArraysClose(x, [1, 2, 3])); }); @@ -47,8 +56,6 @@ describeWithFlags('invert_matrix', ALL_ENVS, () => { it('invert a matrix', () => { const m = tf.tensor2d([1, 3, 3, 1, 4, 3, 1, 3, 4], [3, 3]); const inv = invertMatrix(m); - m.print(); - inv.print(); expect(() => expectArraysClose(inv, [7, -3, -3, -1, 1, 0, -1, 0, 1])); }); }); @@ -57,8 +64,6 @@ describeWithFlags('adjoint_matrix', ALL_ENVS, () => { it('computes the adjoint of a matrix', () => { const m = tf.tensor2d([1, 3, 3, 1, 4, 3, 1, 3, 4], [3, 3]); const a = adjM(m); - m.print(); - a.print(); expect(() => expectArraysClose(a, [7, -3, -3, -1, 1, 0, -1, 0, 1])); }); }); From ad6e06161c573cfe6d2433cb65dae244268e4795 Mon Sep 17 00:00:00 2001 From: Jason Shin Date: Sat, 1 Dec 2018 18:04:26 +1100 Subject: [PATCH 7/7] fix: circular dependency issue --- src/ops/linalg_solve.ts | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/ops/linalg_solve.ts b/src/ops/linalg_solve.ts index 3a8e413193..d660c26fd4 100644 --- a/src/ops/linalg_solve.ts +++ b/src/ops/linalg_solve.ts @@ -19,7 +19,8 @@ * Linear algebra resolution. */ -import {scalar, split, tensor1d} from '..'; +import { scalar, tensor1d } from './tensor_ops'; +import { split } from './array_ops'; import {ENV} from '../environment'; import {Scalar, Tensor, Tensor2D} from '../tensor'; import {assert} from '../util';