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 7db812d9d..83a41b80d 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 @@ -96,14 +96,20 @@ public class DoubleLinearOpsTensorAlgebra : val size = this.shape.size val commonShape = this.shape.sliceArray(0 until size - 2) val (n, m) = this.shape.sliceArray(size - 2 until size) - val resU = zeros(commonShape + intArrayOf(n, min(n, m))) + val resU = zeros(commonShape + intArrayOf(min(n, m), n)) val resS = zeros(commonShape + intArrayOf(min(n, m))) val resV = zeros(commonShape + intArrayOf(min(n, m), m)) for ((matrix, USV) in this.matrixSequence() - .zip(resU.matrixSequence().zip(resS.vectorSequence().zip(resV.matrixSequence())))) - svdHelper(matrix.asTensor(), USV, m, n) - return Triple(resU, resS, resV.transpose(size - 2, size - 1)) + .zip(resU.matrixSequence().zip(resS.vectorSequence().zip(resV.matrixSequence())))) { + val size = matrix.shape.reduce { acc, i -> acc * i } + val curMatrix = DoubleTensor( + matrix.shape, + matrix.buffer.array().slice(matrix.bufferStart until matrix.bufferStart + size).toDoubleArray() + ) + svdHelper(curMatrix, USV, m, n) + } + return Triple(resU.transpose(size - 2, size - 1), resS, resV) } override fun DoubleTensor.symEig(eigenvectors: Boolean): Pair { 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 e1d5c1819..1bdfee2d5 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 @@ -225,10 +225,10 @@ internal inline fun DoubleLinearOpsTensorAlgebra.svd1d(a: DoubleTensor, epsilon: val b: DoubleTensor if (n > m) { b = a.transpose(0, 1).dot(a) - v = DoubleTensor(intArrayOf(m), getRandomNormals(m, 0)) + v = DoubleTensor(intArrayOf(m), getRandomUnitVector(m, 0)) } else { b = a.dot(a.transpose(0, 1)) - v = DoubleTensor(intArrayOf(n), getRandomNormals(n, 0)) + v = DoubleTensor(intArrayOf(n), getRandomUnitVector(n, 0)) } var lastV: DoubleTensor @@ -284,7 +284,13 @@ internal inline fun DoubleLinearOpsTensorAlgebra.svdHelper( val s = res.map { it.first }.toDoubleArray() val uBuffer = res.map { it.second }.flatMap { it.buffer.array().toList() }.toDoubleArray() val vBuffer = res.map { it.third }.flatMap { it.buffer.array().toList() }.toDoubleArray() - uBuffer.copyInto(matrixU.buffer.array()) - s.copyInto(matrixS.buffer.array()) - vBuffer.copyInto(matrixV.buffer.array()) + for (i in uBuffer.indices) { + matrixU.buffer.array()[matrixU.bufferStart + i] = uBuffer[i] + } + for (i in s.indices) { + matrixS.buffer.array()[matrixS.bufferStart + i] = s[i] + } + for (i in vBuffer.indices) { + matrixV.buffer.array()[matrixV.bufferStart + i] = vBuffer[i] + } } diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/utils.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/utils.kt index c99f9d156..392abd1c2 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/utils.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/utils.kt @@ -3,6 +3,7 @@ package space.kscience.kmath.tensors.core import space.kscience.kmath.samplers.GaussianSampler import space.kscience.kmath.stat.RandomGenerator import space.kscience.kmath.structures.* +import kotlin.math.sqrt /** * Returns a reference to [IntArray] containing all of the elements of this [Buffer]. @@ -42,6 +43,12 @@ internal inline fun getRandomNormals(n: Int, seed: Long): DoubleArray { return distribution.sample(generator).nextBufferBlocking(n).toDoubleArray() } +internal inline fun getRandomUnitVector(n: Int, seed: Long): DoubleArray { + val unnorm = getRandomNormals(n, seed) + val norm = sqrt(unnorm.map { it * it }.sum()) + return unnorm.map { it / norm }.toDoubleArray() +} + internal inline fun minusIndexFrom(n: Int, i: Int) : Int = if (i >= 0) i else { val ii = n + i check(ii >= 0) { 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 56f9332f6..35390968e 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 @@ -126,9 +126,9 @@ class TestDoubleLinearOpsTensorAlgebra { testSVDFor(fromArray(intArrayOf(2, 2), doubleArrayOf(-1.0, 0.0, 239.0, 238.0))) } - @Test @Ignore + @Test fun testBatchedSVD() = DoubleLinearOpsTensorAlgebra { - val tensor = randNormal(intArrayOf(7, 5, 3), 0) + val tensor = randNormal(intArrayOf(1, 15, 4, 7, 5, 3), 0) val (tensorU, tensorS, tensorV) = tensor.svd() val tensorSVD = tensorU dot (diagonalEmbedding(tensorS) dot tensorV) assertTrue(tensor.eq(tensorSVD))