From b46e8c5fe23bb64fa6868363848e2e13e0c1edfb Mon Sep 17 00:00:00 2001 From: Roland Grinis Date: Wed, 14 Apr 2021 22:13:54 +0100 Subject: [PATCH] LU and det refactored --- .../core/DoubleLinearOpsTensorAlgebra.kt | 64 ++++++------------- .../kscience/kmath/tensors/core/linutils.kt | 43 +++++++++++-- .../core/TestDoubleLinearOpsAlgebra.kt | 8 +-- 3 files changed, 59 insertions(+), 56 deletions(-) diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt index a6bd5d5f7..e72948c84 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt @@ -10,43 +10,15 @@ public class DoubleLinearOpsTensorAlgebra : LinearOpsTensorAlgebra, DoubleTensorAlgebra() { - override fun DoubleTensor.inv(): DoubleTensor = invLU() + override fun DoubleTensor.inv(): DoubleTensor = invLU(1e-9) - override fun DoubleTensor.det(): DoubleTensor = detLU() + override fun DoubleTensor.det(): DoubleTensor = detLU(1e-9) - internal fun DoubleTensor.luForDet(forDet: Boolean = false): Pair { - checkSquareMatrix(shape) + public fun DoubleTensor.lu(epsilon: Double): Pair = + computeLU(this, epsilon) ?: + throw RuntimeException("Tensor contains matrices which are singular at precision $epsilon") - val luTensor = copy() - - 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 } - ) - - for ((lu, pivots) in luTensor.matrixSequence().zip(pivotsTensor.vectorSequence())) - try { - luHelper(lu.as2D(), pivots.as1D(), m) - } catch (e: RuntimeException) { - if (forDet) { - lu.as2D()[intArrayOf(0, 0)] = 0.0 - } else { - throw IllegalStateException("LUP decomposition can't be performed") - } - } - - - return Pair(luTensor, pivotsTensor) - } - - override fun DoubleTensor.lu(): Pair { - return luForDet(false) - } + override fun DoubleTensor.lu(): Pair = lu(1e-9) override fun luPivot( luTensor: DoubleTensor, @@ -79,9 +51,7 @@ public class DoubleLinearOpsTensorAlgebra : public fun DoubleTensor.cholesky(epsilon: Double): DoubleTensor { checkSquareMatrix(shape) - checkPositiveDefinite(this) - //checkPositiveDefinite(this, epsilon) - + checkPositiveDefinite(this, epsilon) val n = shape.last() val lTensor = zeroesLike() @@ -139,15 +109,19 @@ public class DoubleLinearOpsTensorAlgebra : val shp = s.shape + intArrayOf(1) val utv = u.transpose() dot v val n = s.shape.last() - for( matrix in utv.matrixSequence()) - cleanSymHelper(matrix.as2D(),n) + for (matrix in utv.matrixSequence()) + cleanSymHelper(matrix.as2D(), n) val eig = (utv dot s.view(shp)).view(s.shape) return Pair(eig, v) } - public fun DoubleTensor.detLU(): DoubleTensor { - val (luTensor, pivotsTensor) = luForDet(forDet = true) + public fun DoubleTensor.detLU(epsilon: Double = 1e-9): DoubleTensor { + + checkSquareMatrix(this.shape) + val luTensor = this.copy() + val pivotsTensor = this.setUpPivots() + val n = shape.size val detTensorShape = IntArray(n - 1) { i -> shape[i] } @@ -160,15 +134,15 @@ public class DoubleLinearOpsTensorAlgebra : ) luTensor.matrixSequence().zip(pivotsTensor.vectorSequence()).forEachIndexed { index, (lu, pivots) -> - resBuffer[index] = luMatrixDet(lu.as2D(), pivots.as1D()) + resBuffer[index] = if (luHelper(lu.as2D(), pivots.as1D(), epsilon)) + 0.0 else luMatrixDet(lu.as2D(), pivots.as1D()) } return detTensor } - public fun DoubleTensor.invLU(): DoubleTensor { - //TODO("Andrei the det is non-zero") - val (luTensor, pivotsTensor) = lu() + public fun DoubleTensor.invLU(epsilon: Double = 1e-9): DoubleTensor { + val (luTensor, pivotsTensor) = lu(epsilon) val invTensor = luTensor.zeroesLike() val seq = luTensor.matrixSequence().zip(pivotsTensor.vectorSequence()).zip(invTensor.matrixSequence()) diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/linutils.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/linutils.kt index f2c9d8c76..ee159b614 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/linutils.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/linutils.kt @@ -61,7 +61,13 @@ internal inline fun dotHelper( } } -internal inline fun luHelper(lu: MutableStructure2D, pivots: MutableStructure1D, m: Int) { +internal inline fun luHelper( + lu: MutableStructure2D, + pivots: MutableStructure1D, + epsilon: Double): Boolean { + + val m = lu.rowNum + for (row in 0..m) pivots[row] = row for (i in 0 until m) { @@ -69,16 +75,15 @@ internal inline fun luHelper(lu: MutableStructure2D, pivots: MutableStru var maxInd = i for (k in i until m) { - val absA = kotlin.math.abs(lu[k, i]) + val absA = abs(lu[k, i]) if (absA > maxVal) { maxVal = absA maxInd = k } } - if (abs(maxVal) < 1e-9) { - throw RuntimeException() - } + if (abs(maxVal) < epsilon) + return true // matrix is singular if (maxInd != i) { @@ -103,6 +108,34 @@ internal inline fun luHelper(lu: MutableStructure2D, pivots: MutableStru } } } + return false +} + +internal inline fun BufferedTensor.setUpPivots(): IntTensor { + val n = this.shape.size + val m = this.shape.last() + val pivotsShape = IntArray(n - 1) { i -> this.shape[i] } + pivotsShape[n - 2] = m + 1 + + return IntTensor( + pivotsShape, + IntArray(pivotsShape.reduce(Int::times)) { 0 } + ) +} + +internal inline fun DoubleLinearOpsTensorAlgebra.computeLU( + tensor: DoubleTensor, + epsilon: Double): Pair? { + + checkSquareMatrix(tensor.shape) + val luTensor = tensor.copy() + val pivotsTensor = tensor.setUpPivots() + + for ((lu, pivots) in luTensor.matrixSequence().zip(pivotsTensor.vectorSequence())) + if(luHelper(lu.as2D(), pivots.as1D(), epsilon)) + return null + + return Pair(luTensor, pivotsTensor) } internal inline fun pivInit( diff --git a/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt b/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt index 75ff12355..d19d0b6f6 100644 --- a/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt +++ b/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt @@ -1,8 +1,6 @@ package space.kscience.kmath.tensors.core -import space.kscience.kmath.structures.toList import kotlin.math.abs -import kotlin.test.Ignore import kotlin.test.Test import kotlin.test.assertEquals import kotlin.test.assertTrue @@ -45,7 +43,7 @@ class TestDoubleLinearOpsTensorAlgebra { ) ) - assertTrue { abs(m.det().value() - expectedValue) < 1e-5} + assertTrue { abs(m.det().value() - expectedValue) < 1e-5 } } @Test @@ -57,7 +55,7 @@ class TestDoubleLinearOpsTensorAlgebra { ) ) - assertTrue { abs(m.det().value() - expectedValue) < 1e-5} + assertTrue { abs(m.det().value() - expectedValue) < 1e-5 } } @Test @@ -144,11 +142,9 @@ class TestDoubleLinearOpsTensorAlgebra { val sigma = (tensor dot tensor.transpose()) + diagonalEmbedding( fromArray(intArrayOf(2, 5), DoubleArray(10) { 0.1 }) ) - //checkPositiveDefinite(sigma) sigma must be positive definite val low = sigma.cholesky() val sigmChol = low dot low.transpose() assertTrue(sigma.eq(sigmChol)) - } @Test