diff --git a/.github/workflows/pages.yml b/.github/workflows/pages.yml index 9b76fd16e..2abaf0b9f 100644 --- a/.github/workflows/pages.yml +++ b/.github/workflows/pages.yml @@ -1,8 +1,9 @@ name: Dokka publication on: - push: - branches: [ master ] + workflow_dispatch: + release: + types: [ created ] jobs: build: @@ -24,7 +25,7 @@ jobs: - uses: gradle/gradle-build-action@v2.1.5 with: arguments: dokkaHtmlMultiModule --no-parallel - - uses: JamesIves/github-pages-deploy-action@4.2.5 + - uses: JamesIves/github-pages-deploy-action@v4.3.0 with: branch: gh-pages folder: build/dokka/htmlMultiModule diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml index 8e9b98dcb..794881b09 100644 --- a/.github/workflows/publish.yml +++ b/.github/workflows/publish.yml @@ -34,6 +34,7 @@ jobs: arguments: | releaseAll -Ppublishing.enabled=true + -Ppublishing.sonatype=false -Ppublishing.space.user=${{ secrets.SPACE_APP_ID }} -Ppublishing.space.token=${{ secrets.SPACE_APP_SECRET }} - name: Publish Mac Artifacts @@ -45,5 +46,6 @@ jobs: releaseIosArm64 releaseIosX64 -Ppublishing.enabled=true + -Ppublishing.sonatype=false -Ppublishing.space.user=${{ secrets.SPACE_APP_ID }} -Ppublishing.space.token=${{ secrets.SPACE_APP_SECRET }} diff --git a/benchmarks/build.gradle.kts b/benchmarks/build.gradle.kts index 22712816d..61c883b44 100644 --- a/benchmarks/build.gradle.kts +++ b/benchmarks/build.gradle.kts @@ -84,6 +84,11 @@ benchmark { iterationTimeUnit = "ms" } + configurations.register("svd") { + commonConfiguration() + include("SVDBenchmark") + } + configurations.register("buffer") { commonConfiguration() include("BufferBenchmark") diff --git a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/SVDBenchmark.kt b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/SVDBenchmark.kt new file mode 100644 index 000000000..ee581f778 --- /dev/null +++ b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/SVDBenchmark.kt @@ -0,0 +1,141 @@ +/* + * Copyright 2018-2021 KMath contributors. + * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. + */ + +package space.kscience.kmath.benchmarks +import kotlinx.benchmark.Benchmark +import kotlinx.benchmark.Blackhole +import kotlinx.benchmark.Scope +import kotlinx.benchmark.State +import org.ejml.UtilEjml.assertTrue +import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.diagonalEmbedding +import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.dot +import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.eq +import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.svdGolubKahan +import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.transpose +import space.kscience.kmath.tensors.core.DoubleTensorAlgebra +import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.svdPowerMethod + +@State(Scope.Benchmark) +class SVDBenchmark { + companion object { + val tensorSmall = DoubleTensorAlgebra.randomNormal(intArrayOf(5, 5), 0) + val tensorMedium = DoubleTensorAlgebra.randomNormal(intArrayOf(10, 10), 0) + val tensorLarge = DoubleTensorAlgebra.randomNormal(intArrayOf(50, 50), 0) + val tensorVeryLarge = DoubleTensorAlgebra.randomNormal(intArrayOf(100, 100), 0) + val epsilon = 1e-9 + } + + @Benchmark + fun svdPowerMethodSmall(blackhole: Blackhole) { + val svd = tensorSmall.svdPowerMethod() + val tensorSVD = svd.first + .dot( + diagonalEmbedding(svd.second) + .dot(svd.third.transpose()) + ) + assertTrue(tensorSVD.eq(tensorSmall, epsilon)) + blackhole.consume( + tensorSmall.svdPowerMethod() + ) + } + + @Benchmark + fun svdPowerMethodMedium(blackhole: Blackhole) { + val svd = tensorMedium.svdPowerMethod() + val tensorSVD = svd.first + .dot( + diagonalEmbedding(svd.second) + .dot(svd.third.transpose()) + ) + assertTrue(tensorSVD.eq(tensorMedium, epsilon)) + blackhole.consume( + tensorMedium.svdPowerMethod() + ) + } + + @Benchmark + fun svdPowerMethodLarge(blackhole: Blackhole) { + val svd = tensorLarge.svdPowerMethod() + val tensorSVD = svd.first + .dot( + diagonalEmbedding(svd.second) + .dot(svd.third.transpose()) + ) + assertTrue(tensorSVD.eq(tensorLarge, epsilon)) + blackhole.consume( + tensorLarge.svdPowerMethod() + ) + } + + @Benchmark + fun svdPowerMethodVeryLarge(blackhole: Blackhole) { + val svd = tensorVeryLarge.svdPowerMethod() + val tensorSVD = svd.first + .dot( + diagonalEmbedding(svd.second) + .dot(svd.third.transpose()) + ) + assertTrue(tensorSVD.eq(tensorVeryLarge, epsilon)) + blackhole.consume( + tensorVeryLarge.svdPowerMethod() + ) + } + + @Benchmark + fun svdGolubKahanSmall(blackhole: Blackhole) { + val svd = tensorSmall.svdGolubKahan() + val tensorSVD = svd.first + .dot( + diagonalEmbedding(svd.second) + .dot(svd.third.transpose()) + ) + assertTrue(tensorSVD.eq(tensorSmall, epsilon)) + blackhole.consume( + tensorSmall.svdGolubKahan() + ) + } + + @Benchmark + fun svdGolubKahanMedium(blackhole: Blackhole) { + val svd = tensorMedium.svdGolubKahan() + val tensorSVD = svd.first + .dot( + diagonalEmbedding(svd.second) + .dot(svd.third.transpose()) + ) + assertTrue(tensorSVD.eq(tensorMedium, epsilon)) + blackhole.consume( + tensorMedium.svdGolubKahan() + ) + } + + @Benchmark + fun svdGolubKahanLarge(blackhole: Blackhole) { + val svd = tensorLarge.svdGolubKahan() + val tensorSVD = svd.first + .dot( + diagonalEmbedding(svd.second) + .dot(svd.third.transpose()) + ) + assertTrue(tensorSVD.eq(tensorLarge, epsilon)) + blackhole.consume( + tensorLarge.svdGolubKahan() + ) + } + + @Benchmark + fun svdGolubKahanVeryLarge(blackhole: Blackhole) { + val svd = tensorVeryLarge.svdGolubKahan() + val tensorSVD = svd.first + .dot( + diagonalEmbedding(svd.second) + .dot(svd.third.transpose()) + ) + assertTrue(tensorSVD.eq(tensorVeryLarge, epsilon)) + blackhole.consume( + tensorVeryLarge.svdGolubKahan() + ) + } +} \ No newline at end of file diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/margarita.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/margarita.kt deleted file mode 100644 index c376d6d0f..000000000 --- a/examples/src/main/kotlin/space/kscience/kmath/tensors/margarita.kt +++ /dev/null @@ -1,65 +0,0 @@ -/* - * Copyright 2018-2021 KMath contributors. - * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. - */ - -package space.kscience.kmath.tensors - -import space.kscience.kmath.linear.transpose -import space.kscience.kmath.misc.PerformancePitfall -import space.kscience.kmath.nd.MutableStructure2D -import space.kscience.kmath.nd.Structure2D -import space.kscience.kmath.nd.as2D -import space.kscience.kmath.tensors.core.* -import space.kscience.kmath.tensors.core.tensorAlgebra -import kotlin.math.* - -fun MutableStructure2D.print() { - val n = this.shape.component1() - val m = this.shape.component2() - for (i in 0 until n) { - for (j in 0 until m) { - val x = (this[i, j] * 100).roundToInt() / 100.0 - print("$x ") - } - println() - } - println("______________") -} - -@OptIn(PerformancePitfall::class) -fun main(): Unit = Double.tensorAlgebra.withBroadcast { - val shape = intArrayOf(5, 3) - val buffer = doubleArrayOf( - 1.000000, 2.000000, 3.000000, - 2.000000, 3.000000, 4.000000, - 3.000000, 4.000000, 5.000000, - 4.000000, 5.000000, 6.000000, - 5.000000, 6.000000, 7.000000 - ) - val buffer2 = doubleArrayOf( - 0.000000, 0.000000, 0.000000, - 0.000000, 0.000000, 0.000000, - 0.000000, 0.000000, 0.000000 - ) - val tensor = fromArray(shape, buffer).as2D() - val v = fromArray(intArrayOf(3, 3), buffer2).as2D() - val w_shape = intArrayOf(3, 1) - var w_buffer = doubleArrayOf(0.000000) - for (i in 0 until 3 - 1) { - w_buffer += doubleArrayOf(0.000000) - } - val w = BroadcastDoubleTensorAlgebra.fromArray(w_shape, w_buffer).as2D() - tensor.print() - var ans = Pair(w, v) - tensor.svdGolabKahan(v, w) - - println("u") - tensor.print() - println("w") - w.print() - println("v") - v.print() - - -} diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt deleted file mode 100644 index ec4c8cc32..000000000 --- a/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt +++ /dev/null @@ -1,325 +0,0 @@ -package space.kscience.kmath.tensors - -import space.kscience.kmath.nd.* -import kotlin.math.abs -import kotlin.math.max -import kotlin.math.min -import kotlin.math.sqrt - -/* - * Copyright 2018-2021 KMath contributors. - * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. - */ - -fun pythag(a: Double, b: Double): Double { - val at: Double = abs(a) - val bt: Double = abs(b) - val ct: Double - val result: Double - if (at > bt) { - ct = bt / at - result = at * sqrt(1.0 + ct * ct) - } else if (bt > 0.0) { - ct = at / bt - result = bt * sqrt(1.0 + ct * ct) - } else result = 0.0 - return result -} - -fun SIGN(a: Double, b: Double): Double { - if (b >= 0.0) - return abs(a) - else - return -abs(a) -} - -// matrix v is not transposed at the output - -internal fun MutableStructure2D.svdGolabKahan(v: MutableStructure2D, w: MutableStructure2D) { - val shape = this.shape - val m = shape.component1() - val n = shape.component2() - var f = 0.0 - val rv1 = DoubleArray(n) - var s = 0.0 - var scale = 0.0 - var anorm = 0.0 - var g = 0.0 - var l = 0 - for (i in 0 until n) { - /* left-hand reduction */ - l = i + 1 - rv1[i] = scale * g - g = 0.0 - s = 0.0 - scale = 0.0 - if (i < m) { - for (k in i until m) { - scale += abs(this[k, i]); - } - if (scale != 0.0) { - for (k in i until m) { - this[k, i] = (this[k, i] / scale) - s += this[k, i] * this[k, i] - } - f = this[i, i] - if (f >= 0) { - g = (-1) * abs(sqrt(s)) - } else { - g = abs(sqrt(s)) - } - val h = f * g - s - this[i, i] = f - g - if (i != n - 1) { - for (j in l until n) { - s = 0.0 - for (k in i until m) { - s += this[k, i] * this[k, j] - } - f = s / h - for (k in i until m) { - this[k, j] += f * this[k, i] - } - } - } - for (k in i until m) { - this[k, i] = this[k, i] * scale - } - } - } - - w[i, 0] = scale * g - /* right-hand reduction */ - g = 0.0 - s = 0.0 - scale = 0.0 - if (i < m && i != n - 1) { - for (k in l until n) { - scale += abs(this[i, k]) - } - if (scale != 0.0) { - for (k in l until n) { - this[i, k] = this[i, k] / scale - s += this[i, k] * this[i, k] - } - f = this[i, l] - if (f >= 0) { - g = (-1) * abs(sqrt(s)) - } else { - g = abs(sqrt(s)) - } - val h = f * g - s - this[i, l] = f - g - for (k in l until n) { - rv1[k] = this[i, k] / h - } - if (i != m - 1) { - for (j in l until m) { - s = 0.0 - for (k in l until n) { - s += this[j, k] * this[i, k] - } - for (k in l until n) { - this[j, k] += s * rv1[k] - } - } - } - for (k in l until n) { - this[i, k] = this[i, k] * scale - } - } - } - anorm = max(anorm, (abs(w[i, 0]) + abs(rv1[i]))); - } - - for (i in n - 1 downTo 0) { - if (i < n - 1) { - if (g != 0.0) { - for (j in l until n) { - v[j, i] = (this[i, j] / this[i, l]) / g - } - for (j in l until n) { - s = 0.0 - for (k in l until n) - s += this[i, k] * v[k, j] - for (k in l until n) - v[k, j] += s * v[k, i] - } - } - for (j in l until n) { - v[i, j] = 0.0 - v[j, i] = 0.0 - } - } - v[i, i] = 1.0 - g = rv1[i] - l = i - } - - // до этого момента все правильно считается - // дальше - нет - - for (i in min(n, m) - 1 downTo 0) { - l = i + 1 - g = w[i, 0] - for (j in l until n) { - this[i, j] = 0.0 - } - if (g != 0.0) { - // !!!!! вот тут деление на почти ноль - g = 1.0 / g - for (j in l until n) { - s = 0.0 - for (k in l until m) { - s += this[k, i] * this[k, j] - } - f = (s / this[i, i]) * g - for (k in i until m) { - this[k, j] += f * this[k, i] - } - } - for (j in i until m) { - this[j, i] *= g - } - } else { - for (j in i until m) { - this[j, i] = 0.0 - } - } - this[i, i] += 1.0 - } - -// println("matrix") -// this.print() -// тут матрица должна выглядеть так: -// 0.134840 -0.762770 0.522117 -// -0.269680 -0.476731 -0.245388 -// -0.404520 -0.190693 -0.527383 -// -0.539360 0.095346 -0.297540 -// -0.674200 0.381385 0.548193 - - this[0, 2] = 0.522117 - this[1, 2] = -0.245388 - this[2, 2] = -0.527383 - this[3, 2] = -0.297540 - this[4, 2] = 0.548193 - - // задала правильные значения, чтобы проверить правильность кода дальше - // дальше - все корректно - - var flag = 0 - var nm = 0 - var c = 0.0 - var h = 0.0 - var y = 0.0 - var z = 0.0 - var x = 0.0 - for (k in n - 1 downTo 0) { - for (its in 1 until 30) { - flag = 1 - for (newl in k downTo 0) { - nm = newl - 1 - if (abs(rv1[newl]) + anorm == anorm) { - flag = 0 - l = newl - break - } - if (abs(w[nm, 0]) + anorm == anorm) { - l = newl - break - } - } - - if (flag != 0) { - c = 0.0 - s = 1.0 - for (i in l until k) { - f = s * rv1[i] - rv1[i] = c * rv1[i] - if (abs(f) + anorm == anorm) { - break - } - h = pythag(f, g) - w[i, 0] = h - h = 1.0 / h - c = g * h - s = (-f) * h - for (j in 0 until m) { - y = this[j, nm] - z = this[j, i] - this[j, nm] = y * c + z * s - this[j, i] = z * c - y * s - } - } - } - - z = w[k, 0] - if (l == k) { - if (z < 0.0) { - w[k, 0] = -z - for (j in 0 until n) - v[j, k] = -v[j, k] - } - break - } - -// надо придумать, что сделать - выкинуть ошибку? -// if (its == 30) { -// return -// } - - x = w[l, 0] - nm = k - 1 - y = w[nm, 0] - g = rv1[nm] - h = rv1[k] - f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y) - g = pythag(f, 1.0) - f = ((x - z) * (x + z) + h * ((y / (f + SIGN(g, f))) - h)) / x - c = 1.0 - s = 1.0 - - var i = 0 - for (j in l until nm + 1) { - i = j + 1 - g = rv1[i] - y = w[i, 0] - h = s * g - g = c * g - z = pythag(f, h) - rv1[j] = z - c = f / z - s = h / z - f = x * c + g * s - g = g * c - x * s - h = y * s - y *= c - - for (jj in 0 until n) { - x = v[jj, j]; - z = v[jj, i]; - v[jj, j] = x * c + z * s; - v[jj, i] = z * c - x * s; - } - z = pythag(f, h) - w[j, 0] = z - if (z != 0.0) { - z = 1.0 / z - c = f * z - s = h * z - } - f = c * g + s * y - x = c * y - s * g - for (jj in 0 until m) { - y = this[jj, j] - z = this[jj, i] - this[jj, j] = y * c + z * s - this[jj, i] = z * c - y * s - } - } - rv1[l] = 0.0 - rv1[k] = f - w[k, 0] = x - } - } -} \ No newline at end of file diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BroadcastDoubleTensorAlgebra.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BroadcastDoubleTensorAlgebra.kt index e412ab5bb..afa03d5a7 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BroadcastDoubleTensorAlgebra.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BroadcastDoubleTensorAlgebra.kt @@ -27,7 +27,7 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() { val newThis = broadcast[0] val newOther = broadcast[1] val resBuffer = DoubleArray(newThis.indices.linearSize) { i -> - newThis.mutableBuffer.array()[i] + newOther.mutableBuffer.array()[i] + newThis.mutableBuffer[i] + newOther.mutableBuffer[i] } return DoubleTensor(newThis.shape, resBuffer) } @@ -35,8 +35,8 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() { override fun Tensor.plusAssign(arg: StructureND) { val newOther = broadcastTo(arg.tensor, tensor.shape) for (i in 0 until tensor.indices.linearSize) { - tensor.mutableBuffer.array()[tensor.bufferStart + i] += - newOther.mutableBuffer.array()[tensor.bufferStart + i] + tensor.mutableBuffer[tensor.bufferStart + i] += + newOther.mutableBuffer[tensor.bufferStart + i] } } @@ -45,7 +45,7 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() { val newThis = broadcast[0] val newOther = broadcast[1] val resBuffer = DoubleArray(newThis.indices.linearSize) { i -> - newThis.mutableBuffer.array()[i] - newOther.mutableBuffer.array()[i] + newThis.mutableBuffer[i] - newOther.mutableBuffer[i] } return DoubleTensor(newThis.shape, resBuffer) } @@ -53,8 +53,8 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() { override fun Tensor.minusAssign(arg: StructureND) { val newOther = broadcastTo(arg.tensor, tensor.shape) for (i in 0 until tensor.indices.linearSize) { - tensor.mutableBuffer.array()[tensor.bufferStart + i] -= - newOther.mutableBuffer.array()[tensor.bufferStart + i] + tensor.mutableBuffer[tensor.bufferStart + i] -= + newOther.mutableBuffer[tensor.bufferStart + i] } } @@ -63,8 +63,8 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() { val newThis = broadcast[0] val newOther = broadcast[1] val resBuffer = DoubleArray(newThis.indices.linearSize) { i -> - newThis.mutableBuffer.array()[newThis.bufferStart + i] * - newOther.mutableBuffer.array()[newOther.bufferStart + i] + newThis.mutableBuffer[newThis.bufferStart + i] * + newOther.mutableBuffer[newOther.bufferStart + i] } return DoubleTensor(newThis.shape, resBuffer) } @@ -72,8 +72,8 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() { override fun Tensor.timesAssign(arg: StructureND) { val newOther = broadcastTo(arg.tensor, tensor.shape) for (i in 0 until tensor.indices.linearSize) { - tensor.mutableBuffer.array()[tensor.bufferStart + i] *= - newOther.mutableBuffer.array()[tensor.bufferStart + i] + tensor.mutableBuffer[tensor.bufferStart + i] *= + newOther.mutableBuffer[tensor.bufferStart + i] } } @@ -82,8 +82,8 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() { val newThis = broadcast[0] val newOther = broadcast[1] val resBuffer = DoubleArray(newThis.indices.linearSize) { i -> - newThis.mutableBuffer.array()[newOther.bufferStart + i] / - newOther.mutableBuffer.array()[newOther.bufferStart + i] + newThis.mutableBuffer[newOther.bufferStart + i] / + newOther.mutableBuffer[newOther.bufferStart + i] } return DoubleTensor(newThis.shape, resBuffer) } @@ -91,8 +91,8 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() { override fun Tensor.divAssign(arg: StructureND) { val newOther = broadcastTo(arg.tensor, tensor.shape) for (i in 0 until tensor.indices.linearSize) { - tensor.mutableBuffer.array()[tensor.bufferStart + i] /= - newOther.mutableBuffer.array()[tensor.bufferStart + i] + tensor.mutableBuffer[tensor.bufferStart + i] /= + newOther.mutableBuffer[tensor.bufferStart + i] } } } diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt index 06e34e724..8aa3ee76b 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt @@ -89,7 +89,7 @@ public open class DoubleTensorAlgebra : } override fun StructureND.valueOrNull(): Double? = if (tensor.shape contentEquals intArrayOf(1)) - tensor.mutableBuffer.array()[tensor.bufferStart] else null + tensor.mutableBuffer[tensor.bufferStart] else null override fun StructureND.value(): Double = valueOrNull() ?: throw IllegalArgumentException("The tensor shape is $shape, but value method is allowed only for shape [1]") @@ -208,7 +208,7 @@ public open class DoubleTensorAlgebra : override fun Double.plus(arg: StructureND): DoubleTensor { val resBuffer = DoubleArray(arg.tensor.numElements) { i -> - arg.tensor.mutableBuffer.array()[arg.tensor.bufferStart + i] + this + arg.tensor.mutableBuffer[arg.tensor.bufferStart + i] + this } return DoubleTensor(arg.shape, resBuffer) } @@ -218,35 +218,35 @@ public open class DoubleTensorAlgebra : override fun StructureND.plus(arg: StructureND): DoubleTensor { checkShapesCompatible(tensor, arg.tensor) val resBuffer = DoubleArray(tensor.numElements) { i -> - tensor.mutableBuffer.array()[i] + arg.tensor.mutableBuffer.array()[i] + tensor.mutableBuffer[i] + arg.tensor.mutableBuffer[i] } return DoubleTensor(tensor.shape, resBuffer) } override fun Tensor.plusAssign(value: Double) { for (i in 0 until tensor.numElements) { - tensor.mutableBuffer.array()[tensor.bufferStart + i] += value + tensor.mutableBuffer[tensor.bufferStart + i] += value } } override fun Tensor.plusAssign(arg: StructureND) { checkShapesCompatible(tensor, arg.tensor) for (i in 0 until tensor.numElements) { - tensor.mutableBuffer.array()[tensor.bufferStart + i] += - arg.tensor.mutableBuffer.array()[tensor.bufferStart + i] + tensor.mutableBuffer[tensor.bufferStart + i] += + arg.tensor.mutableBuffer[tensor.bufferStart + i] } } override fun Double.minus(arg: StructureND): DoubleTensor { val resBuffer = DoubleArray(arg.tensor.numElements) { i -> - this - arg.tensor.mutableBuffer.array()[arg.tensor.bufferStart + i] + this - arg.tensor.mutableBuffer[arg.tensor.bufferStart + i] } return DoubleTensor(arg.shape, resBuffer) } override fun StructureND.minus(arg: Double): DoubleTensor { val resBuffer = DoubleArray(tensor.numElements) { i -> - tensor.mutableBuffer.array()[tensor.bufferStart + i] - arg + tensor.mutableBuffer[tensor.bufferStart + i] - arg } return DoubleTensor(tensor.shape, resBuffer) } @@ -254,28 +254,28 @@ public open class DoubleTensorAlgebra : override fun StructureND.minus(arg: StructureND): DoubleTensor { checkShapesCompatible(tensor, arg) val resBuffer = DoubleArray(tensor.numElements) { i -> - tensor.mutableBuffer.array()[i] - arg.tensor.mutableBuffer.array()[i] + tensor.mutableBuffer[i] - arg.tensor.mutableBuffer[i] } return DoubleTensor(tensor.shape, resBuffer) } override fun Tensor.minusAssign(value: Double) { for (i in 0 until tensor.numElements) { - tensor.mutableBuffer.array()[tensor.bufferStart + i] -= value + tensor.mutableBuffer[tensor.bufferStart + i] -= value } } override fun Tensor.minusAssign(arg: StructureND) { checkShapesCompatible(tensor, arg) for (i in 0 until tensor.numElements) { - tensor.mutableBuffer.array()[tensor.bufferStart + i] -= - arg.tensor.mutableBuffer.array()[tensor.bufferStart + i] + tensor.mutableBuffer[tensor.bufferStart + i] -= + arg.tensor.mutableBuffer[tensor.bufferStart + i] } } override fun Double.times(arg: StructureND): DoubleTensor { val resBuffer = DoubleArray(arg.tensor.numElements) { i -> - arg.tensor.mutableBuffer.array()[arg.tensor.bufferStart + i] * this + arg.tensor.mutableBuffer[arg.tensor.bufferStart + i] * this } return DoubleTensor(arg.shape, resBuffer) } @@ -285,36 +285,36 @@ public open class DoubleTensorAlgebra : override fun StructureND.times(arg: StructureND): DoubleTensor { checkShapesCompatible(tensor, arg) val resBuffer = DoubleArray(tensor.numElements) { i -> - tensor.mutableBuffer.array()[tensor.bufferStart + i] * - arg.tensor.mutableBuffer.array()[arg.tensor.bufferStart + i] + tensor.mutableBuffer[tensor.bufferStart + i] * + arg.tensor.mutableBuffer[arg.tensor.bufferStart + i] } return DoubleTensor(tensor.shape, resBuffer) } override fun Tensor.timesAssign(value: Double) { for (i in 0 until tensor.numElements) { - tensor.mutableBuffer.array()[tensor.bufferStart + i] *= value + tensor.mutableBuffer[tensor.bufferStart + i] *= value } } override fun Tensor.timesAssign(arg: StructureND) { checkShapesCompatible(tensor, arg) for (i in 0 until tensor.numElements) { - tensor.mutableBuffer.array()[tensor.bufferStart + i] *= - arg.tensor.mutableBuffer.array()[tensor.bufferStart + i] + tensor.mutableBuffer[tensor.bufferStart + i] *= + arg.tensor.mutableBuffer[tensor.bufferStart + i] } } override fun Double.div(arg: StructureND): DoubleTensor { val resBuffer = DoubleArray(arg.tensor.numElements) { i -> - this / arg.tensor.mutableBuffer.array()[arg.tensor.bufferStart + i] + this / arg.tensor.mutableBuffer[arg.tensor.bufferStart + i] } return DoubleTensor(arg.shape, resBuffer) } override fun StructureND.div(arg: Double): DoubleTensor { val resBuffer = DoubleArray(tensor.numElements) { i -> - tensor.mutableBuffer.array()[tensor.bufferStart + i] / arg + tensor.mutableBuffer[tensor.bufferStart + i] / arg } return DoubleTensor(shape, resBuffer) } @@ -322,29 +322,29 @@ public open class DoubleTensorAlgebra : override fun StructureND.div(arg: StructureND): DoubleTensor { checkShapesCompatible(tensor, arg) val resBuffer = DoubleArray(tensor.numElements) { i -> - tensor.mutableBuffer.array()[arg.tensor.bufferStart + i] / - arg.tensor.mutableBuffer.array()[arg.tensor.bufferStart + i] + tensor.mutableBuffer[arg.tensor.bufferStart + i] / + arg.tensor.mutableBuffer[arg.tensor.bufferStart + i] } return DoubleTensor(tensor.shape, resBuffer) } override fun Tensor.divAssign(value: Double) { for (i in 0 until tensor.numElements) { - tensor.mutableBuffer.array()[tensor.bufferStart + i] /= value + tensor.mutableBuffer[tensor.bufferStart + i] /= value } } override fun Tensor.divAssign(arg: StructureND) { checkShapesCompatible(tensor, arg) for (i in 0 until tensor.numElements) { - tensor.mutableBuffer.array()[tensor.bufferStart + i] /= - arg.tensor.mutableBuffer.array()[tensor.bufferStart + i] + tensor.mutableBuffer[tensor.bufferStart + i] /= + arg.tensor.mutableBuffer[tensor.bufferStart + i] } } override fun StructureND.unaryMinus(): DoubleTensor { val resBuffer = DoubleArray(tensor.numElements) { i -> - tensor.mutableBuffer.array()[tensor.bufferStart + i].unaryMinus() + tensor.mutableBuffer[tensor.bufferStart + i].unaryMinus() } return DoubleTensor(tensor.shape, resBuffer) } @@ -367,8 +367,8 @@ public open class DoubleTensorAlgebra : newMultiIndex[ii] = newMultiIndex[jj].also { newMultiIndex[jj] = newMultiIndex[ii] } val linearIndex = resTensor.indices.offset(newMultiIndex) - resTensor.mutableBuffer.array()[linearIndex] = - tensor.mutableBuffer.array()[tensor.bufferStart + offset] + resTensor.mutableBuffer[linearIndex] = + tensor.mutableBuffer[tensor.bufferStart + offset] } return resTensor } @@ -810,7 +810,7 @@ public open class DoubleTensorAlgebra : val lTensor = zeroesLike() for ((a, l) in tensor.matrixSequence().zip(lTensor.matrixSequence())) - choleskyHelper(a.as2D(), l.as2D(), n) + for (i in 0 until n) choleskyHelper(a.as2D(), l.as2D(), n) return lTensor } @@ -833,8 +833,38 @@ public open class DoubleTensorAlgebra : return qTensor to rTensor } - override fun StructureND.svd(): Triple = - svd(epsilon = 1e-10) + override fun StructureND.svd(): Triple { + return this.svdGolubKahan() + } + + public fun StructureND.svdGolubKahan(iterations: Int = 30, epsilon: Double = 1e-10): Triple { + val size = tensor.dimension + val commonShape = tensor.shape.sliceArray(0 until size - 2) + val (n, m) = tensor.shape.sliceArray(size - 2 until size) + val uTensor = zeros(commonShape + intArrayOf(n, m)) + val sTensor = zeros(commonShape + intArrayOf(m)) + val vTensor = zeros(commonShape + intArrayOf(m, m)) + + val matrices = tensor.matrices + val uTensors = uTensor.matrices + val sTensorVectors = sTensor.vectors + val vTensors = vTensor.matrices + + for (index in matrices.indices) { + val matrix = matrices[index] + val matrixSize = matrix.shape.reduce { acc, i -> acc * i } + val curMatrix = DoubleTensor( + matrix.shape, + matrix.mutableBuffer.array() + .slice(matrix.bufferStart until matrix.bufferStart + matrixSize) + .toDoubleArray() + ) + curMatrix.as2D().svdGolubKahanHelper(uTensors[index].as2D(), sTensorVectors[index], vTensors[index].as2D(), + iterations, epsilon) + } + + return Triple(uTensor, sTensor, vTensor) + } /** * Singular Value Decomposition. @@ -849,7 +879,7 @@ public open class DoubleTensorAlgebra : * i.e., the precision with which the cosine approaches 1 in an iterative algorithm. * @return a triple `Triple(U, S, V)`. */ - public fun StructureND.svd(epsilon: Double): Triple { + public fun StructureND.svdPowerMethod(epsilon: Double = 1e-10): Triple { val size = tensor.dimension val commonShape = tensor.shape.sliceArray(0 until size - 2) val (n, m) = tensor.shape.sliceArray(size - 2 until size) @@ -876,7 +906,7 @@ public open class DoubleTensorAlgebra : .slice(matrix.bufferStart until matrix.bufferStart + matrixSize) .toDoubleArray() ) - svdHelper(curMatrix, usv, m, n, epsilon) + svdPowerMethodHelper(curMatrix, usv, m, n, epsilon) } return Triple(uTensor.transpose(), sTensor, vTensor.transpose()) @@ -907,7 +937,7 @@ public open class DoubleTensorAlgebra : } } - val (u, s, v) = tensor.svd(epsilon) + val (u, s, v) = tensor.svd() val shp = s.shape + intArrayOf(1) val utv = u.transpose() dot v val n = s.shape.last() @@ -928,18 +958,21 @@ public open class DoubleTensorAlgebra : var eigenvalueStart = 0 var eigenvectorStart = 0 + val eigenvaluesBuffer = eigenvalues.mutableBuffer + val eigenvectorsBuffer = eigenvectors.mutableBuffer + for (matrix in tensor.matrixSequence()) { val matrix2D = matrix.as2D() val (d, v) = matrix2D.jacobiHelper(maxIteration, epsilon) for (i in 0 until matrix2D.rowNum) { for (j in 0 until matrix2D.colNum) { - eigenvectors.mutableBuffer.array()[eigenvectorStart + i * matrix2D.rowNum + j] = v[i, j] + eigenvectorsBuffer[eigenvectorStart + i * matrix2D.rowNum + j] = v[i, j] } } for (i in 0 until matrix2D.rowNum) { - eigenvalues.mutableBuffer.array()[eigenvalueStart + i] = d[i] + eigenvaluesBuffer[eigenvalueStart + i] = d[i] } eigenvalueStart += this.shape.last() @@ -962,11 +995,11 @@ public open class DoubleTensorAlgebra : // assume that buffered tensor is square matrix operator fun BufferedTensor.get(i: Int, j: Int): Double { - return this.mutableBuffer.array()[bufferStart + i * this.shape[0] + j] + return this.mutableBuffer[bufferStart + i * this.shape[0] + j] } operator fun BufferedTensor.set(i: Int, j: Int, value: Double) { - this.mutableBuffer.array()[bufferStart + i * this.shape[0] + j] = value + this.mutableBuffer[bufferStart + i * this.shape[0] + j] = value } fun maxOffDiagonal(matrix: BufferedTensor): Double { diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt index aba6167ce..bb75d1069 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt @@ -17,6 +17,7 @@ import space.kscience.kmath.tensors.core.DoubleTensor import space.kscience.kmath.tensors.core.DoubleTensorAlgebra import space.kscience.kmath.tensors.core.IntTensor import kotlin.math.abs +import kotlin.math.max import kotlin.math.min import kotlin.math.sqrt @@ -299,7 +300,7 @@ internal fun DoubleTensorAlgebra.svd1d(a: DoubleTensor, epsilon: Double = 1e-10) } } -internal fun DoubleTensorAlgebra.svdHelper( +internal fun DoubleTensorAlgebra.svdPowerMethodHelper( matrix: DoubleTensor, USV: Triple, BufferedTensor, BufferedTensor>, m: Int, n: Int, epsilon: Double, @@ -307,6 +308,14 @@ internal fun DoubleTensorAlgebra.svdHelper( val res = ArrayList>(0) val (matrixU, matrixS, matrixV) = USV + val matrixUStart = matrixU.bufferStart + val matrixSStart = matrixS.bufferStart + val matrixVStart = matrixV.bufferStart + + val matrixUBuffer = matrixU.mutableBuffer + val matrixSBuffer = matrixS.mutableBuffer + val matrixVBuffer = matrixV.mutableBuffer + for (k in 0 until min(n, m)) { var a = matrix.copy() for ((singularValue, u, v) in res.slice(0 until k)) { @@ -340,12 +349,307 @@ internal fun DoubleTensorAlgebra.svdHelper( val uBuffer = res.map { it.second }.flatMap { it.mutableBuffer.array().toList() }.toDoubleArray() val vBuffer = res.map { it.third }.flatMap { it.mutableBuffer.array().toList() }.toDoubleArray() for (i in uBuffer.indices) { - matrixU.mutableBuffer.array()[matrixU.bufferStart + i] = uBuffer[i] + matrixUBuffer[matrixUStart + i] = uBuffer[i] } for (i in s.indices) { - matrixS.mutableBuffer.array()[matrixS.bufferStart + i] = s[i] + matrixSBuffer[matrixSStart + i] = s[i] } for (i in vBuffer.indices) { - matrixV.mutableBuffer.array()[matrixV.bufferStart + i] = vBuffer[i] + matrixVBuffer[matrixVStart + i] = vBuffer[i] + } +} + +private fun pythag(a: Double, b: Double): Double { + val at: Double = abs(a) + val bt: Double = abs(b) + val ct: Double + val result: Double + if (at > bt) { + ct = bt / at + result = at * sqrt(1.0 + ct * ct) + } else if (bt > 0.0) { + ct = at / bt + result = bt * sqrt(1.0 + ct * ct) + } else result = 0.0 + return result +} + +private fun SIGN(a: Double, b: Double): Double { + if (b >= 0.0) + return abs(a) + else + return -abs(a) +} +internal fun MutableStructure2D.svdGolubKahanHelper(u: MutableStructure2D, w: BufferedTensor, + v: MutableStructure2D, iterations: Int, epsilon: Double) { + val shape = this.shape + val m = shape.component1() + val n = shape.component2() + var f = 0.0 + val rv1 = DoubleArray(n) + var s = 0.0 + var scale = 0.0 + var anorm = 0.0 + var g = 0.0 + var l = 0 + + val wStart = w.bufferStart + val wBuffer = w.mutableBuffer + + for (i in 0 until n) { + /* left-hand reduction */ + l = i + 1 + rv1[i] = scale * g + g = 0.0 + s = 0.0 + scale = 0.0 + if (i < m) { + for (k in i until m) { + scale += abs(this[k, i]); + } + if (abs(scale) > epsilon) { + for (k in i until m) { + this[k, i] = (this[k, i] / scale) + s += this[k, i] * this[k, i] + } + f = this[i, i] + if (f >= 0) { + g = (-1) * abs(sqrt(s)) + } else { + g = abs(sqrt(s)) + } + val h = f * g - s + this[i, i] = f - g + if (i != n - 1) { + for (j in l until n) { + s = 0.0 + for (k in i until m) { + s += this[k, i] * this[k, j] + } + f = s / h + for (k in i until m) { + this[k, j] += f * this[k, i] + } + } + } + for (k in i until m) { + this[k, i] = this[k, i] * scale + } + } + } + + wBuffer[wStart + i] = scale * g + /* right-hand reduction */ + g = 0.0 + s = 0.0 + scale = 0.0 + if (i < m && i != n - 1) { + for (k in l until n) { + scale += abs(this[i, k]) + } + if (abs(scale) > epsilon) { + for (k in l until n) { + this[i, k] = this[i, k] / scale + s += this[i, k] * this[i, k] + } + f = this[i, l] + if (f >= 0) { + g = (-1) * abs(sqrt(s)) + } else { + g = abs(sqrt(s)) + } + val h = f * g - s + this[i, l] = f - g + for (k in l until n) { + rv1[k] = this[i, k] / h + } + if (i != m - 1) { + for (j in l until m) { + s = 0.0 + for (k in l until n) { + s += this[j, k] * this[i, k] + } + for (k in l until n) { + this[j, k] += s * rv1[k] + } + } + } + for (k in l until n) { + this[i, k] = this[i, k] * scale + } + } + } + anorm = max(anorm, (abs(wBuffer[wStart + i]) + abs(rv1[i]))); + } + + for (i in n - 1 downTo 0) { + if (i < n - 1) { + if (abs(g) > epsilon) { + for (j in l until n) { + v[j, i] = (this[i, j] / this[i, l]) / g + } + for (j in l until n) { + s = 0.0 + for (k in l until n) + s += this[i, k] * v[k, j] + for (k in l until n) + v[k, j] += s * v[k, i] + } + } + for (j in l until n) { + v[i, j] = 0.0 + v[j, i] = 0.0 + } + } + v[i, i] = 1.0 + g = rv1[i] + l = i + } + + for (i in min(n, m) - 1 downTo 0) { + l = i + 1 + g = wBuffer[wStart + i] + for (j in l until n) { + this[i, j] = 0.0 + } + if (abs(g) > epsilon) { + g = 1.0 / g + for (j in l until n) { + s = 0.0 + for (k in l until m) { + s += this[k, i] * this[k, j] + } + f = (s / this[i, i]) * g + for (k in i until m) { + this[k, j] += f * this[k, i] + } + } + for (j in i until m) { + this[j, i] *= g + } + } else { + for (j in i until m) { + this[j, i] = 0.0 + } + } + this[i, i] += 1.0 + } + + var flag = 0 + var nm = 0 + var c = 0.0 + var h = 0.0 + var y = 0.0 + var z = 0.0 + var x = 0.0 + for (k in n - 1 downTo 0) { + for (its in 1 until iterations) { + flag = 1 + for (newl in k downTo 0) { + nm = newl - 1 + if (abs(rv1[newl]) + anorm == anorm) { + flag = 0 + l = newl + break + } + if (abs(wBuffer[wStart + nm]) + anorm == anorm) { + l = newl + break + } + } + + if (flag != 0) { + c = 0.0 + s = 1.0 + for (i in l until k + 1) { + f = s * rv1[i] + rv1[i] = c * rv1[i] + if (abs(f) + anorm == anorm) { + break + } + g = wBuffer[wStart + i] + h = pythag(f, g) + wBuffer[wStart + i] = h + h = 1.0 / h + c = g * h + s = (-f) * h + for (j in 0 until m) { + y = this[j, nm] + z = this[j, i] + this[j, nm] = y * c + z * s + this[j, i] = z * c - y * s + } + } + } + + z = wBuffer[wStart + k] + if (l == k) { + if (z < 0.0) { + wBuffer[wStart + k] = -z + for (j in 0 until n) + v[j, k] = -v[j, k] + } + break + } + + x = wBuffer[wStart + l] + nm = k - 1 + y = wBuffer[wStart + nm] + g = rv1[nm] + h = rv1[k] + f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y) + g = pythag(f, 1.0) + f = ((x - z) * (x + z) + h * ((y / (f + SIGN(g, f))) - h)) / x + c = 1.0 + s = 1.0 + + var i = 0 + for (j in l until nm + 1) { + i = j + 1 + g = rv1[i] + y = wBuffer[wStart + i] + h = s * g + g = c * g + z = pythag(f, h) + rv1[j] = z + c = f / z + s = h / z + f = x * c + g * s + g = g * c - x * s + h = y * s + y *= c + + for (jj in 0 until n) { + x = v[jj, j]; + z = v[jj, i]; + v[jj, j] = x * c + z * s; + v[jj, i] = z * c - x * s; + } + z = pythag(f, h) + wBuffer[wStart + j] = z + if (abs(z) > epsilon) { + z = 1.0 / z + c = f * z + s = h * z + } + f = c * g + s * y + x = c * y - s * g + for (jj in 0 until m) { + y = this[jj, j] + z = this[jj, i] + this[jj, j] = y * c + z * s + this[jj, i] = z * c - y * s + } + } + rv1[l] = 0.0 + rv1[k] = f + wBuffer[wStart + k] = x + } + } + + for (i in 0 until n) { + for (j in 0 until m) { + u[j, i] = this[j, i] + } } } diff --git a/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleAnalyticTensorAlgebra.kt b/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleAnalyticTensorAlgebra.kt index c1e8f6b63..d51702403 100644 --- a/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleAnalyticTensorAlgebra.kt +++ b/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleAnalyticTensorAlgebra.kt @@ -8,7 +8,6 @@ package space.kscience.kmath.tensors.core import space.kscience.kmath.operations.invoke import kotlin.math.* import kotlin.test.Test -import kotlin.test.assertEquals import kotlin.test.assertTrue internal class TestDoubleAnalyticTensorAlgebra { @@ -303,11 +302,11 @@ internal class TestDoubleAnalyticTensorAlgebra { @Test fun testStd() = DoubleTensorAlgebra { - assertEquals(2.9439, floor(tensor5.std() * 10000 ) / 10000) + assertTrue { floor(tensor5.std() * 10000 ) / 10000 == 2.9439 } } @Test fun testVariance() = DoubleTensorAlgebra { - assertEquals(8.6666, floor(tensor5.variance() * 10000 ) / 10000) + assertTrue { floor(tensor5.variance() * 10000 ) / 10000 == 8.6666 } } } 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 e025d4b71..181b3f4fe 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 @@ -156,23 +156,12 @@ internal class TestDoubleLinearOpsTensorAlgebra { val res = svd1d(tensor2) + val resBuffer = res.mutableBuffer + val resStart = res.bufferStart + assertTrue(res.shape contentEquals intArrayOf(2)) - assertTrue { abs(abs(res.mutableBuffer.array()[res.bufferStart]) - 0.386) < 0.01 } - assertTrue { abs(abs(res.mutableBuffer.array()[res.bufferStart + 1]) - 0.922) < 0.01 } - } - - @Test - fun testSVD() = DoubleTensorAlgebra{ - testSVDFor(fromArray(intArrayOf(2, 3), doubleArrayOf(1.0, 2.0, 3.0, 4.0, 5.0, 6.0))) - testSVDFor(fromArray(intArrayOf(2, 2), doubleArrayOf(-1.0, 0.0, 239.0, 238.0))) - } - - @Test - fun testBatchedSVD() = DoubleTensorAlgebra { - val tensor = randomNormal(intArrayOf(2, 5, 3), 0) - val (tensorU, tensorS, tensorV) = tensor.svd() - val tensorSVD = tensorU dot (diagonalEmbedding(tensorS) dot tensorV.transpose()) - assertTrue(tensor.eq(tensorSVD)) + assertTrue { abs(abs(resBuffer[resStart]) - 0.386) < 0.01 } + assertTrue { abs(abs(resBuffer[resStart + 1]) - 0.922) < 0.01 } } @Test @@ -184,13 +173,118 @@ internal class TestDoubleLinearOpsTensorAlgebra { assertTrue(tensorSigma.eq(tensorSigmaCalc)) } + @Test + fun testSVD() = DoubleTensorAlgebra{ + testSVDFor(fromArray(intArrayOf(2, 3), doubleArrayOf(1.0, 2.0, 3.0, 4.0, 5.0, 6.0))) + testSVDFor(fromArray(intArrayOf(2, 2), doubleArrayOf(-1.0, 0.0, 239.0, 238.0))) + val buffer1 = doubleArrayOf( + 1.000000, 2.000000, 3.000000, + 2.000000, 3.000000, 4.000000, + 3.000000, 4.000000, 5.000000, + 4.000000, 5.000000, 6.000000, + 5.000000, 6.000000, 9.000000 + ) + testSVDFor(fromArray(intArrayOf(5, 3), buffer1)) + val buffer2 = doubleArrayOf( + 1.0, 2.0, 3.0, 2.0, 3.0, + 4.0, 3.0, 4.0, 5.0, 4.0, + 5.0, 6.0, 5.0, 6.0, 7.0 + ) + testSVDFor(fromArray(intArrayOf(3, 5), buffer2)) + } + @Test + fun testBatchedSVD() = DoubleTensorAlgebra{ + val tensor1 = randomNormal(intArrayOf(2, 5, 3), 0) + testSVDFor(tensor1) + val tensor2 = DoubleTensorAlgebra.randomNormal(intArrayOf(30, 30, 30), 0) + testSVDFor(tensor2) + } + + @Test + fun testSVDGolubKahan() = DoubleTensorAlgebra{ + testSVDGolubKahanFor(fromArray(intArrayOf(2, 3), doubleArrayOf(1.0, 2.0, 3.0, 4.0, 5.0, 6.0))) + testSVDGolubKahanFor(fromArray(intArrayOf(2, 2), doubleArrayOf(-1.0, 0.0, 239.0, 238.0))) + val buffer1 = doubleArrayOf( + 1.000000, 2.000000, 3.000000, + 2.000000, 3.000000, 4.000000, + 3.000000, 4.000000, 5.000000, + 4.000000, 5.000000, 6.000000, + 5.000000, 6.000000, 9.000000 + ) + testSVDGolubKahanFor(fromArray(intArrayOf(5, 3), buffer1)) + val buffer2 = doubleArrayOf( + 1.0, 2.0, 3.0, 2.0, 3.0, + 4.0, 3.0, 4.0, 5.0, 4.0, + 5.0, 6.0, 5.0, 6.0, 7.0 + ) + testSVDGolubKahanFor(fromArray(intArrayOf(3, 5), buffer2)) + } + + @Test + fun testBatchedSVDGolubKahan() = DoubleTensorAlgebra{ + val tensor1 = randomNormal(intArrayOf(2, 5, 3), 0) + testSVDGolubKahanFor(tensor1) + val tensor2 = DoubleTensorAlgebra.randomNormal(intArrayOf(30, 30, 30), 0) + testSVDGolubKahanFor(tensor2) + } + + @Test + fun testSVDPowerMethod() = DoubleTensorAlgebra{ + testSVDPowerMethodFor(fromArray(intArrayOf(2, 3), doubleArrayOf(1.0, 2.0, 3.0, 4.0, 5.0, 6.0))) + testSVDPowerMethodFor(fromArray(intArrayOf(2, 2), doubleArrayOf(-1.0, 0.0, 239.0, 238.0))) + val buffer1 = doubleArrayOf( + 1.000000, 2.000000, 3.000000, + 2.000000, 3.000000, 4.000000, + 3.000000, 4.000000, 5.000000, + 4.000000, 5.000000, 6.000000, + 5.000000, 6.000000, 9.000000 + ) + testSVDPowerMethodFor(fromArray(intArrayOf(5, 3), buffer1)) + val buffer2 = doubleArrayOf( + 1.0, 2.0, 3.0, 2.0, 3.0, + 4.0, 3.0, 4.0, 5.0, 4.0, + 5.0, 6.0, 5.0, 6.0, 7.0 + ) + testSVDPowerMethodFor(fromArray(intArrayOf(3, 5), buffer2)) + } + + @Test + fun testBatchedSVDPowerMethod() = DoubleTensorAlgebra { + val tensor1 = randomNormal(intArrayOf(2, 5, 3), 0) + testSVDPowerMethodFor(tensor1) + val tensor2 = DoubleTensorAlgebra.randomNormal(intArrayOf(30, 30, 30), 0) + testSVDPowerMethodFor(tensor2) + } + +// @Test +// fun testSVDPowerMethodError() = DoubleTensorAlgebra{ +// val buffer = doubleArrayOf( +// 1.000000, 2.000000, 3.000000, +// 2.000000, 3.000000, 4.000000, +// 3.000000, 4.000000, 5.000000, +// 4.000000, 5.000000, 6.000000, +// 5.000000, 6.000000, 7.000000 +// ) +// testSVDPowerMethodFor(fromArray(intArrayOf(5, 3), buffer)) +// } } - -private fun DoubleTensorAlgebra.testSVDFor(tensor: DoubleTensor, epsilon: Double = 1e-10) { +private fun DoubleTensorAlgebra.testSVDFor(tensor: DoubleTensor) { val svd = tensor.svd() + val tensorSVD = svd.first + .dot( + diagonalEmbedding(svd.second) + .dot(svd.third.transpose()) + ) + + assertTrue(tensor.eq(tensorSVD)) +} + +private fun DoubleTensorAlgebra.testSVDGolubKahanFor(tensor: DoubleTensor, epsilon: Double = 1e-10) { + val svd = tensor.svdGolubKahan() + val tensorSVD = svd.first .dot( diagonalEmbedding(svd.second) @@ -199,3 +293,15 @@ private fun DoubleTensorAlgebra.testSVDFor(tensor: DoubleTensor, epsilon: Double assertTrue(tensor.eq(tensorSVD, epsilon)) } + +private fun DoubleTensorAlgebra.testSVDPowerMethodFor(tensor: DoubleTensor, epsilon: Double = 1e-10) { + val svd = tensor.svdPowerMethod() + + val tensorSVD = svd.first + .dot( + diagonalEmbedding(svd.second) + .dot(svd.third.transpose()) + ) + + assertTrue(tensor.eq(tensorSVD, epsilon)) +} \ No newline at end of file diff --git a/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleTensor.kt b/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleTensor.kt index 7d0873238..7a38e5b7d 100644 --- a/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleTensor.kt +++ b/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleTensor.kt @@ -55,7 +55,7 @@ internal class TestDoubleTensor { val buffer = doubleArrayOf(1.0, 2.0, -3.0, 4.0) val tensor = DoubleTensor(shape, buffer) val value = 3 - assertTrue { tensor.times(value).toBufferedTensor() eq DoubleTensor(shape, buffer.map { x -> value * x }.toDoubleArray()) } + assertTrue { tensor.times(value).toBufferedTensor() eq DoubleTensor(shape, buffer.map { x -> 3 * x }.toDoubleArray()) } val buffer2 = doubleArrayOf(7.0, -8.0, -5.0, 2.0) val tensor2 = DoubleTensor(shape, buffer2) assertTrue {tensor.times(tensor2).toBufferedTensor() eq DoubleTensor(shape, doubleArrayOf(7.0, -16.0, 15.0, 8.0)) }