From fea53af0eeaa3522a5754bfe83688f010a37a5f4 Mon Sep 17 00:00:00 2001
From: Andrei Kislitsyn <kisandy@yandex.ru>
Date: Wed, 24 Mar 2021 14:00:47 +0300
Subject: [PATCH 1/3] fix
---
.../space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt | 2 +-
.../kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt | 3 +--
2 files changed, 2 insertions(+), 3 deletions(-)
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 c4fa3e0a8..e412af7c6 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
@@ -14,7 +14,7 @@ public interface LinearOpsTensorAlgebra<T, TensorType : TensorStructure<T>, Inde
public fun TensorType.qr(): TensorType
//https://pytorch.org/docs/stable/generated/torch.lu.html
- public fun TensorType.lu(tol: T): Pair<TensorType, IndexTensorType>
+ public fun TensorType.lu(): Pair<TensorType, IndexTensorType>
//https://pytorch.org/docs/stable/generated/torch.lu_unpack.html
public fun luPivot(luTensor: TensorType, pivotsTensor: IndexTensorType): Triple<TensorType, TensorType, TensorType>
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 e0abc49b7..d090ce79c 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
@@ -11,7 +11,7 @@ public class DoubleLinearOpsTensorAlgebra :
TODO("Not yet implemented")
}
- override fun DoubleTensor.lu(tol: Double): Pair<DoubleTensor, IntTensor> {
+ override fun DoubleTensor.lu(): Pair<DoubleTensor, IntTensor> {
checkSquareMatrix(shape)
@@ -138,7 +138,6 @@ public class DoubleLinearOpsTensorAlgebra :
TODO("Not yet implemented")
}
-
override fun DoubleTensor.svd(): Triple<DoubleTensor, DoubleTensor, DoubleTensor> {
TODO("Not yet implemented")
}
--
2.34.1
From 206bcfc909d56a859ec783a30cc75a8ade178786 Mon Sep 17 00:00:00 2001
From: Andrei Kislitsyn <kisandy@yandex.ru>
Date: Wed, 24 Mar 2021 18:08:36 +0300
Subject: [PATCH 2/3] 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<T, TensorType : TensorStructure<T>, 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<TensorType, IndexTensorType>
//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<Double>, pivots: Structure1D<Int>): 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<Double>,
+ pivots: Structure1D<Int>,
+ invMatrix : MutableStructure2D<Double>
+ ): 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 <R> 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<Double, Dou
}
override fun DoubleTensor.det(): DoubleTensor {
- TODO("ANDREI")
+ TODO()
+ /*
+ checkSquareMatrix(shape)
+
+ val n = shape.size
+ val m = shape.last()
+
+ val detTensorShape = IntArray(n - 1) { i -> 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<Double, Dou
}
public fun DoubleTensor.contentEquals(other: DoubleTensor, eqFunction: (Double, Double) -> 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<Double, Dou
public fun DoubleTensor.eq(other: DoubleTensor, eqFunction: (Double, Double) -> 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) }
+ }
+}
--
2.34.1
From ab8137000146c52f10b3213c960b03917a0ada4f Mon Sep 17 00:00:00 2001
From: Andrei Kislitsyn <kisandy@yandex.ru>
Date: Wed, 24 Mar 2021 18:42:41 +0300
Subject: [PATCH 3/3] fixes
---
.../kmath/tensors/LinearOpsTensorAlgebra.kt | 2 +-
.../core/DoubleLinearOpsTensorAlgebra.kt | 22 +++++++------------
2 files changed, 9 insertions(+), 15 deletions(-)
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 f551d524a..e412af7c6 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<T, TensorType : TensorStructure<T>, Inde
//https://pytorch.org/docs/stable/linalg.html#torch.linalg.qr
public fun TensorType.qr(): TensorType
- //htt ps://pytorch.org/docs/stable/generated/torch.lu.html
+ //https://pytorch.org/docs/stable/generated/torch.lu.html
public fun TensorType.lu(): Pair<TensorType, IndexTensorType>
//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 c8c4c413f..0b9d824f0 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
@@ -10,9 +10,9 @@ public class DoubleLinearOpsTensorAlgebra :
LinearOpsTensorAlgebra<Double, DoubleTensor, IntTensor>,
DoubleTensorAlgebra() {
- override fun DoubleTensor.inv(): DoubleTensor {
- TODO("ANDREI")
- }
+ override fun DoubleTensor.inv(): DoubleTensor = invLU()
+
+ override fun DoubleTensor.det(): DoubleTensor = detLU()
override fun DoubleTensor.lu(): Pair<DoubleTensor, IntTensor> {
@@ -156,15 +156,11 @@ public class DoubleLinearOpsTensorAlgebra :
private fun luMatrixDet(lu: Structure2D<Double>, pivots: Structure1D<Int>): 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
+ return (0 until m).asSequence().map { lu[it, it] }.fold(sign) { left, right -> left * right }
}
public fun DoubleTensor.detLU(): DoubleTensor {
- val (luTensor, pivotsTensor) = this.lu()
+ val (luTensor, pivotsTensor) = lu()
val n = shape.size
val detTensorShape = IntArray(n - 1) { i -> shape[i] }
@@ -211,14 +207,12 @@ public class DoubleLinearOpsTensorAlgebra :
}
public fun DoubleTensor.invLU(): DoubleTensor {
- val (luTensor, pivotsTensor) = this.lu()
+ val (luTensor, pivotsTensor) = lu()
val n = shape.size
val invTensor = luTensor.zeroesLike()
- for (
- (luP, invMatrix) in
- luTensor.matrixSequence().zip(pivotsTensor.vectorSequence()).zip(invTensor.matrixSequence())
- ) {
+ val seq = luTensor.matrixSequence().zip(pivotsTensor.vectorSequence()).zip(invTensor.matrixSequence())
+ for ((luP, invMatrix) in seq) {
val (lu, pivots) = luP
luMatrixInv(lu, pivots, invMatrix)
}
--
2.34.1