From 206bcfc909d56a859ec783a30cc75a8ade178786 Mon Sep 17 00:00:00 2001 From: Andrei Kislitsyn Date: Wed, 24 Mar 2021 18:08:36 +0300 Subject: [PATCH] lu inv and det complete + tests --- .../kmath/tensors/LinearOpsTensorAlgebra.kt | 2 +- .../core/DoubleLinearOpsTensorAlgebra.kt | 79 +++++++++++++++++++ .../kmath/tensors/core/DoubleTensorAlgebra.kt | 35 +++++++- .../core/TestDoubleAnalyticTensorAlgebra.kt | 7 +- .../core/TestDoubleLinearOpsAlgebra.kt | 73 +++++++++++++++++ 5 files changed, 187 insertions(+), 9 deletions(-) create mode 100644 kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt index e412af7c6..f551d524a 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt @@ -13,7 +13,7 @@ public interface LinearOpsTensorAlgebra, Inde //https://pytorch.org/docs/stable/linalg.html#torch.linalg.qr public fun TensorType.qr(): TensorType - //https://pytorch.org/docs/stable/generated/torch.lu.html + //htt ps://pytorch.org/docs/stable/generated/torch.lu.html public fun TensorType.lu(): Pair //https://pytorch.org/docs/stable/generated/torch.lu_unpack.html diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt index 7779ad029..c8c4c413f 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt @@ -1,5 +1,8 @@ package space.kscience.kmath.tensors.core +import space.kscience.kmath.nd.MutableStructure2D +import space.kscience.kmath.nd.Structure1D +import space.kscience.kmath.nd.Structure2D import space.kscience.kmath.tensors.LinearOpsTensorAlgebra import kotlin.math.sqrt @@ -20,6 +23,8 @@ public class DoubleLinearOpsTensorAlgebra : val n = shape.size val m = shape.last() val pivotsShape = IntArray(n - 1) { i -> shape[i] } + pivotsShape[n - 2] = m + 1 + val pivotsTensor = IntTensor( pivotsShape, IntArray(pivotsShape.reduce(Int::times)) { 0 } @@ -54,6 +59,8 @@ public class DoubleLinearOpsTensorAlgebra : lu[maxInd, k] = tmp } + pivots[m] += 1 + } for (j in i + 1 until m) { @@ -146,6 +153,78 @@ public class DoubleLinearOpsTensorAlgebra : TODO("ANDREI") } + private fun luMatrixDet(lu: Structure2D, pivots: Structure1D): Double { + val m = lu.shape[0] + val sign = if((pivots[m] - m) % 2 == 0) 1.0 else -1.0 + var det = sign + for (i in 0 until m){ + det *= lu[i, i] + } + return det + } + + public fun DoubleTensor.detLU(): DoubleTensor { + val (luTensor, pivotsTensor) = this.lu() + val n = shape.size + + val detTensorShape = IntArray(n - 1) { i -> shape[i] } + detTensorShape[n - 2] = 1 + val resBuffer = DoubleArray(detTensorShape.reduce(Int::times)) { 0.0 } + + val detTensor = DoubleTensor( + detTensorShape, + resBuffer + ) + + luTensor.matrixSequence().zip(pivotsTensor.vectorSequence()).forEachIndexed { index, (luMatrix, pivots) -> + resBuffer[index] = luMatrixDet(luMatrix, pivots) + } + + return detTensor + } + + private fun luMatrixInv( + lu: Structure2D, + pivots: Structure1D, + invMatrix : MutableStructure2D + ): Unit { + val m = lu.shape[0] + + for (j in 0 until m) { + for (i in 0 until m) { + if (pivots[i] == j){ + invMatrix[i, j] = 1.0 + } + + for (k in 0 until i){ + invMatrix[i, j] -= lu[i, k] * invMatrix[k, j] + } + } + + for (i in m - 1 downTo 0) { + for (k in i + 1 until m) { + invMatrix[i, j] -= lu[i, k] * invMatrix[k, j] + } + invMatrix[i, j] /= lu[i, i] + } + } + } + + public fun DoubleTensor.invLU(): DoubleTensor { + val (luTensor, pivotsTensor) = this.lu() + val n = shape.size + val invTensor = luTensor.zeroesLike() + + for ( + (luP, invMatrix) in + luTensor.matrixSequence().zip(pivotsTensor.vectorSequence()).zip(invTensor.matrixSequence()) + ) { + val (lu, pivots) = luP + luMatrixInv(lu, pivots, invMatrix) + } + + return invTensor + } } public inline fun DoubleLinearOpsTensorAlgebra(block: DoubleLinearOpsTensorAlgebra.() -> R): R = diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt index c67687a09..ff65ca027 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt @@ -1,5 +1,8 @@ package space.kscience.kmath.tensors.core +import space.kscience.kmath.linear.Matrix +import space.kscience.kmath.nd.MutableStructure2D +import space.kscience.kmath.nd.Structure2D import space.kscience.kmath.tensors.TensorPartialDivisionAlgebra import kotlin.math.abs @@ -262,7 +265,31 @@ public open class DoubleTensorAlgebra : TensorPartialDivisionAlgebra shape[i] } + detTensorShape[n - 1] = 1 + val resBuffer = DoubleArray(detTensorShape.reduce(Int::times)) { 0.0 } + + val detTensor = DoubleTensor( + detTensorShape, + resBuffer + ) + + this.matrixSequence().forEachIndexed{i, matrix -> + // todo need Matrix determinant algo + // todo resBuffer[i] = matrix.det() + } + + + return detTensor + + */ } override fun DoubleTensor.square(): DoubleTensor { @@ -294,7 +321,7 @@ public open class DoubleTensorAlgebra : TensorPartialDivisionAlgebra Boolean): Boolean { - if (!(this.shape contentEquals other.shape)){ + if (!(this.shape contentEquals other.shape)) { return false } return this.eq(other, eqFunction) @@ -303,10 +330,10 @@ public open class DoubleTensorAlgebra : TensorPartialDivisionAlgebra Boolean): Boolean { // todo broadcasting checking val n = this.linearStructure.size - if (n != other.linearStructure.size){ + if (n != other.linearStructure.size) { return false } - for (i in 0 until n){ + for (i in 0 until n) { if (!eqFunction(this.buffer[this.bufferStart + i], other.buffer[other.bufferStart + i])) { return false } diff --git a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleAnalyticTensorAlgebra.kt b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleAnalyticTensorAlgebra.kt index 8028ce175..4dcdb7848 100644 --- a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleAnalyticTensorAlgebra.kt +++ b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleAnalyticTensorAlgebra.kt @@ -3,7 +3,6 @@ package space.kscience.kmath.tensors.core import kotlin.math.abs import kotlin.math.exp import kotlin.test.Test -import kotlin.test.assertEquals import kotlin.test.assertTrue class TestDoubleAnalyticTensorAlgebra { @@ -16,9 +15,9 @@ class TestDoubleAnalyticTensorAlgebra { return this.map(transform).toDoubleArray() } - fun DoubleArray.deltaEqual(other: DoubleArray, delta: Double = 1e-5): Boolean { + fun DoubleArray.epsEqual(other: DoubleArray, eps: Double = 1e-5): Boolean { for ((elem1, elem2) in this.asSequence().zip(other.asSequence())) { - if (abs(elem1 - elem2) > delta) { + if (abs(elem1 - elem2) > eps) { return false } } @@ -29,7 +28,7 @@ class TestDoubleAnalyticTensorAlgebra { fun testExp() = DoubleAnalyticTensorAlgebra { tensor.exp().let { assertTrue { shape contentEquals it.shape } - assertTrue { buffer.fmap(::exp).deltaEqual(it.buffer.array())} + assertTrue { buffer.fmap(::exp).epsEqual(it.buffer.array())} } } } \ No newline at end of file diff --git a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt new file mode 100644 index 000000000..2da79382f --- /dev/null +++ b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt @@ -0,0 +1,73 @@ +package space.kscience.kmath.tensors.core + +import kotlin.math.abs +import kotlin.math.exp +import kotlin.test.Test +import kotlin.test.assertEquals +import kotlin.test.assertTrue + +class TestDoubleLinearOpsTensorAlgebra { + + private val eps = 1e-5 + + private fun Double.epsEqual(other: Double): Boolean { + return abs(this - other) < eps + } + + fun DoubleArray.epsEqual(other: DoubleArray, eps: Double = 1e-5): Boolean { + for ((elem1, elem2) in this.asSequence().zip(other.asSequence())) { + if (abs(elem1 - elem2) > eps) { + return false + } + } + return true + } + + @Test + fun testDetLU() = DoubleLinearOpsTensorAlgebra { + val tensor = fromArray( + intArrayOf(2, 2, 2), + doubleArrayOf( + 1.0, 3.0, + 1.0, 2.0, + 1.5, 1.0, + 10.0, 2.0 + ) + ) + + val expectedShape = intArrayOf(2, 1) + val expectedBuffer = doubleArrayOf( + -1.0, + -7.0 + ) + val detTensor = tensor.detLU() + + assertTrue { detTensor.shape contentEquals expectedShape } + assertTrue { detTensor.buffer.array().epsEqual(expectedBuffer) } + } + + @Test + fun testInvLU() = DoubleLinearOpsTensorAlgebra { + val tensor = fromArray( + intArrayOf(2, 2, 2), + doubleArrayOf( + 1.0, 0.0, + 0.0, 2.0, + 1.0, 1.0, + 1.0, 0.0 + ) + ) + + val expectedShape = intArrayOf(2, 2, 2) + val expectedBuffer = doubleArrayOf( + 1.0, 0.0, + 0.0, 0.5, + 0.0, 1.0, + 1.0, -1.0 + ) + + val invTensor = tensor.invLU() + assertTrue { invTensor.shape contentEquals expectedShape } + assertTrue { invTensor.buffer.array().epsEqual(expectedBuffer) } + } +}