From 7a72a0b979e04ba928c31e01058ff561cebcc69f Mon Sep 17 00:00:00 2001 From: Ivan Kylchik Date: Sun, 13 Feb 2022 13:16:35 +0300 Subject: [PATCH 1/5] Implement Jacobi algorithm to find eigenvalues --- .../kmath/tensors/core/DoubleTensorAlgebra.kt | 98 ++++++++++++++++++- 1 file changed, 97 insertions(+), 1 deletion(-) 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 ebe7a10b6..418cf16b9 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 @@ -895,7 +895,7 @@ public open class DoubleTensorAlgebra : * and when the cosine approaches 1 in the SVD algorithm. * @return a pair `eigenvalues to eigenvectors`. */ - public fun StructureND.symEig(epsilon: Double): Pair { + public fun StructureND.symEigSvd(epsilon: Double): Pair { checkSymmetric(tensor, epsilon) fun MutableStructure2D.cleanSym(n: Int) { @@ -922,6 +922,102 @@ public open class DoubleTensorAlgebra : return eig to v } + // TODO + // 1. Cyclic method + // 2. Sort eigenvalues + public fun StructureND.symEig(epsilon: Double): Pair { + checkSymmetric(tensor, epsilon) + + val ii = tensor.minusIndex(-2) + val jj = tensor.minusIndex(-1) + val n = tensor.numElements + + val size = this.dimension + val commonShape = this.shape.sliceArray(0 until size - 2) + intArrayOf(1, 1) + + var d = this.copy() + var s = diagonalEmbedding(ones(this.shape.sliceArray(0 until size - 1))) + + do { + // 1. Find max element by abs value that is not on diagonal + val buffer = MutableBuffer.boxing(commonShape.reduce(Int::times)) { Triple(0.0, 0, 0) } + val maxOffDiagonalElements = BufferedTensor(commonShape, buffer, 0) + for (offset in 0 until n) { + val multiIndex = d.linearStructure.index(offset) + if (multiIndex[ii] != multiIndex[jj]) { + val value = d.mutableBuffer.array()[offset] + + val commonIndex = multiIndex.sliceArray(0 until size - 2) + intArrayOf(0, 0) + if (abs(value) > maxOffDiagonalElements[commonIndex].first) { + maxOffDiagonalElements[commonIndex] = Triple(abs(value), multiIndex[ii], multiIndex[jj]) + } + } + } + + // 2. Evaluate "rotation" angle + val angles = zeros(commonShape) + for (offset in 0 until maxOffDiagonalElements.numElements) { + val (_, i, j) = maxOffDiagonalElements.mutableBuffer[offset] + val multiIndex = maxOffDiagonalElements.linearStructure.index(offset) + + val dIJ = d.mutableBuffer[d.linearStructure.offset(multiIndex.also { it[ii] = i; it[jj] = j })] + val dII = d.mutableBuffer[d.linearStructure.offset(multiIndex.also { it[ii] = i; it[jj] = i })] + val dJJ = d.mutableBuffer[d.linearStructure.offset(multiIndex.also { it[ii] = j; it[jj] = j })] + + angles.mutableBuffer.array()[offset] = if (dII == dJJ) { + if (dIJ > 0) PI / 4 else -PI / 4 + } else { + 0.5 * atan(2 * dIJ / (dJJ - dII)) + } + } + + // 3. Build rotation tensor `s1` + val s1 = diagonalEmbedding(ones(this.shape.sliceArray(0 until size - 1))) + for (offset in 0 until n) { + val multiIndex = d.linearStructure.index(offset) + + val commonIndex = multiIndex.sliceArray(0 until size - 2) + intArrayOf(0, 0) + val (_, i, j) = maxOffDiagonalElements[commonIndex] + val angleValue = angles[commonIndex] + s1.mutableBuffer.array()[offset] = when { + multiIndex[ii] == i && multiIndex[jj] == i || multiIndex[ii] == j && multiIndex[jj] == j -> cos(angleValue) + multiIndex[ii] == i && multiIndex[jj] == j -> sin(angleValue) + multiIndex[ii] == j && multiIndex[jj] == i -> -sin(angleValue) + else -> s1.mutableBuffer.array()[offset] + } + } + + // 4. Evaluate new tensor + d = (s1.transpose() dot d) dot s1 + s = s dot s1 + if (d.isDiagonal(epsilon)) break + } while(true) + + val eigenvalues = zeros(d.shape.sliceArray(0 until size - 1)) + for (offset in 0 until n) { + val multiIndex = d.linearStructure.index(offset) + if (multiIndex[ii] == multiIndex[jj]) { + eigenvalues[multiIndex.sliceArray(0 until size - 1)] = d.mutableBuffer.array()[offset] + } + } + + return eigenvalues to s + } + + public fun StructureND.isDiagonal(epsilon: Double = 1e-9): Boolean { + val ii = tensor.minusIndex(-2) + val jj = tensor.minusIndex(-1) + + for (offset in 0 until tensor.numElements) { + val multiIndex = tensor.linearStructure.index(offset) + if (multiIndex[ii] != multiIndex[jj] && abs(tensor.mutableBuffer.array()[offset]) > epsilon) { + return false + } + } + + return true + } + /** * Computes the determinant of a square matrix input, or of each square matrix in a batched input * using LU factorization algorithm. From 7aff774bc166d537c40a7958418ac9a40403d74e Mon Sep 17 00:00:00 2001 From: Ivan Kylchik Date: Sun, 13 Feb 2022 21:49:06 +0300 Subject: [PATCH 2/5] Improve Jacobi algorithm readability by extracting some logic into helper fun --- .../kmath/tensors/core/DoubleTensorAlgebra.kt | 128 +++++++++--------- 1 file changed, 66 insertions(+), 62 deletions(-) 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 418cf16b9..7fa6e885d 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 @@ -10,6 +10,7 @@ package space.kscience.kmath.tensors.core import space.kscience.kmath.misc.PerformancePitfall import space.kscience.kmath.nd.MutableStructure2D +import space.kscience.kmath.nd.Structure2D import space.kscience.kmath.nd.StructureND import space.kscience.kmath.nd.as1D import space.kscience.kmath.nd.as2D @@ -885,7 +886,7 @@ public open class DoubleTensorAlgebra : return Triple(uTensor.transpose(), sTensor, vTensor.transpose()) } - override fun StructureND.symEig(): Pair = symEig(epsilon = 1e-15) + override fun StructureND.symEig(): Pair = symEigJacobi(epsilon = 1e-10) /** * Returns eigenvalues and eigenvectors of a real symmetric matrix input or a batch of real symmetric matrices, @@ -922,96 +923,99 @@ public open class DoubleTensorAlgebra : return eig to v } - // TODO - // 1. Cyclic method - // 2. Sort eigenvalues - public fun StructureND.symEig(epsilon: Double): Pair { + public fun StructureND.symEigJacobi(epsilon: Double): Pair { checkSymmetric(tensor, epsilon) - val ii = tensor.minusIndex(-2) - val jj = tensor.minusIndex(-1) - val n = tensor.numElements - val size = this.dimension - val commonShape = this.shape.sliceArray(0 until size - 2) + intArrayOf(1, 1) + val s = zeros(this.shape) + val eigenvalues = zeros(this.shape.sliceArray(0 until size - 1)) - var d = this.copy() - var s = diagonalEmbedding(ones(this.shape.sliceArray(0 until size - 1))) + var eigenvalueStart = 0 + var eigenvectorStart = 0 + for (matrix in tensor.matrixSequence()) { + matrix.as2D().jacobiHelper(eigenvalues, s, eigenvalueStart, eigenvectorStart, epsilon) + eigenvalueStart += this.shape.last() + eigenvectorStart += this.shape.last() * this.shape.last() + } + // TODO sort eigenvalues + return eigenvalues to s + } + + private fun MutableStructure2D.jacobiHelper( + eigenvalues: DoubleTensor, + eigenvectors: DoubleTensor, + eigenvalueStart: Int, + eigenvectorStart: Int, + epsilon: Double + ) { + var d = this + var s = eye(this.shape[0]) + + // TODO implement cyclic method do { // 1. Find max element by abs value that is not on diagonal - val buffer = MutableBuffer.boxing(commonShape.reduce(Int::times)) { Triple(0.0, 0, 0) } - val maxOffDiagonalElements = BufferedTensor(commonShape, buffer, 0) - for (offset in 0 until n) { - val multiIndex = d.linearStructure.index(offset) - if (multiIndex[ii] != multiIndex[jj]) { - val value = d.mutableBuffer.array()[offset] - - val commonIndex = multiIndex.sliceArray(0 until size - 2) + intArrayOf(0, 0) - if (abs(value) > maxOffDiagonalElements[commonIndex].first) { - maxOffDiagonalElements[commonIndex] = Triple(abs(value), multiIndex[ii], multiIndex[jj]) + var maxOffDiagonalElement = 0.0 + var maxElementIndex = Pair(0, 0) + for (i in 0 until this.rowNum) { + for (j in 0 until this.colNum) { + if (i == j) continue + if (abs(d[i, j]) > maxOffDiagonalElement) { + maxOffDiagonalElement = abs(d[i, j]) + maxElementIndex = i to j } } } // 2. Evaluate "rotation" angle - val angles = zeros(commonShape) - for (offset in 0 until maxOffDiagonalElements.numElements) { - val (_, i, j) = maxOffDiagonalElements.mutableBuffer[offset] - val multiIndex = maxOffDiagonalElements.linearStructure.index(offset) + val dIJ = d[maxElementIndex.first, maxElementIndex.second] + val dII = d[maxElementIndex.first, maxElementIndex.first] + val dJJ = d[maxElementIndex.second, maxElementIndex.second] - val dIJ = d.mutableBuffer[d.linearStructure.offset(multiIndex.also { it[ii] = i; it[jj] = j })] - val dII = d.mutableBuffer[d.linearStructure.offset(multiIndex.also { it[ii] = i; it[jj] = i })] - val dJJ = d.mutableBuffer[d.linearStructure.offset(multiIndex.also { it[ii] = j; it[jj] = j })] - - angles.mutableBuffer.array()[offset] = if (dII == dJJ) { - if (dIJ > 0) PI / 4 else -PI / 4 - } else { - 0.5 * atan(2 * dIJ / (dJJ - dII)) - } + val angle = if (dII == dJJ) { + if (dIJ > 0) PI / 4 else -PI / 4 + } else { + 0.5 * atan(2 * dIJ / (dJJ - dII)) } // 3. Build rotation tensor `s1` - val s1 = diagonalEmbedding(ones(this.shape.sliceArray(0 until size - 1))) - for (offset in 0 until n) { - val multiIndex = d.linearStructure.index(offset) - - val commonIndex = multiIndex.sliceArray(0 until size - 2) + intArrayOf(0, 0) - val (_, i, j) = maxOffDiagonalElements[commonIndex] - val angleValue = angles[commonIndex] - s1.mutableBuffer.array()[offset] = when { - multiIndex[ii] == i && multiIndex[jj] == i || multiIndex[ii] == j && multiIndex[jj] == j -> cos(angleValue) - multiIndex[ii] == i && multiIndex[jj] == j -> sin(angleValue) - multiIndex[ii] == j && multiIndex[jj] == i -> -sin(angleValue) - else -> s1.mutableBuffer.array()[offset] + val s1 = eye(this.rowNum) + for (i in 0 until this.rowNum) { + for (j in 0 until this.colNum) { + s1.mutableBuffer.array()[i * this.rowNum + j] = when { + maxElementIndex.first == i && maxElementIndex.first == j -> cos(angle) + maxElementIndex.second == i && maxElementIndex.second == j -> cos(angle) + maxElementIndex.first == i && maxElementIndex.second == j -> sin(angle) + maxElementIndex.first == j && maxElementIndex.second == i -> -sin(angle) + else -> s1.mutableBuffer.array()[i * this.rowNum + j] + } } } // 4. Evaluate new tensor - d = (s1.transpose() dot d) dot s1 + d = ((s1.transpose() dot d) dot s1).as2D() s = s dot s1 if (d.isDiagonal(epsilon)) break } while(true) - val eigenvalues = zeros(d.shape.sliceArray(0 until size - 1)) - for (offset in 0 until n) { - val multiIndex = d.linearStructure.index(offset) - if (multiIndex[ii] == multiIndex[jj]) { - eigenvalues[multiIndex.sliceArray(0 until size - 1)] = d.mutableBuffer.array()[offset] + // 5. Copy result + for (i in 0 until this.rowNum) { + for (j in 0 until this.colNum) { + eigenvectors.mutableBuffer.array()[eigenvectorStart + i * this.rowNum + j] = s.mutableBuffer.array()[i * this.rowNum + j] } } - return eigenvalues to s + for (i in 0 until this.rowNum) { + eigenvalues.mutableBuffer.array()[eigenvalueStart + i] = d[i, i] + } } - public fun StructureND.isDiagonal(epsilon: Double = 1e-9): Boolean { - val ii = tensor.minusIndex(-2) - val jj = tensor.minusIndex(-1) - - for (offset in 0 until tensor.numElements) { - val multiIndex = tensor.linearStructure.index(offset) - if (multiIndex[ii] != multiIndex[jj] && abs(tensor.mutableBuffer.array()[offset]) > epsilon) { - return false + public fun Structure2D.isDiagonal(epsilon: Double = 1e-9): Boolean { + for (i in 0 until this.rowNum) { + for (j in 0 until this.colNum) { + if (i != j && abs(this[i, j]) > epsilon) { + return false + } } } From b13765ec1928e0f2b2eb6165a8a5a0e1a057b08f Mon Sep 17 00:00:00 2001 From: Ivan Kylchik Date: Sun, 20 Feb 2022 02:21:52 +0300 Subject: [PATCH 3/5] Implement much faster Jacobi algorithm --- benchmarks/build.gradle.kts | 5 + .../benchmarks/TensorAlgebraBenchmark.kt | 37 ++++ .../kmath/tensors/core/DoubleTensorAlgebra.kt | 182 ++++++++++-------- 3 files changed, 149 insertions(+), 75 deletions(-) create mode 100644 benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/TensorAlgebraBenchmark.kt diff --git a/benchmarks/build.gradle.kts b/benchmarks/build.gradle.kts index e46f4a9b4..f8d39b9c5 100644 --- a/benchmarks/build.gradle.kts +++ b/benchmarks/build.gradle.kts @@ -124,6 +124,11 @@ benchmark { include("JafamaBenchmark") } + configurations.register("tensorAlgebra") { + commonConfiguration() + include("TensorAlgebraBenchmark") + } + configurations.register("viktor") { commonConfiguration() include("ViktorBenchmark") diff --git a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/TensorAlgebraBenchmark.kt b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/TensorAlgebraBenchmark.kt new file mode 100644 index 000000000..38e064e53 --- /dev/null +++ b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/TensorAlgebraBenchmark.kt @@ -0,0 +1,37 @@ +/* + * 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 space.kscience.kmath.linear.linearSpace +import space.kscience.kmath.linear.matrix +import space.kscience.kmath.linear.symmetric +import space.kscience.kmath.operations.DoubleField +import space.kscience.kmath.tensors.core.tensorAlgebra +import kotlin.random.Random + +@State(Scope.Benchmark) +internal class TensorAlgebraBenchmark { + companion object { + private val random = Random(12224) + private const val dim = 30 + + private val matrix = DoubleField.linearSpace.matrix(dim, dim).symmetric { _, _ -> random.nextDouble() } + } + + @Benchmark + fun tensorSymEigSvd(blackhole: Blackhole) = with(Double.tensorAlgebra) { + blackhole.consume(matrix.symEigSvd(1e-10)) + } + + @Benchmark + fun tensorSymEigJacobi(blackhole: Blackhole) = with(Double.tensorAlgebra) { + blackhole.consume(matrix.symEigJacobi(50, 1e-10)) + } +} \ No newline at end of file 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 7fa6e885d..675ff4191 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 @@ -9,11 +9,7 @@ package space.kscience.kmath.tensors.core import space.kscience.kmath.misc.PerformancePitfall -import space.kscience.kmath.nd.MutableStructure2D -import space.kscience.kmath.nd.Structure2D -import space.kscience.kmath.nd.StructureND -import space.kscience.kmath.nd.as1D -import space.kscience.kmath.nd.as2D +import space.kscience.kmath.nd.* import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.structures.MutableBuffer import space.kscience.kmath.structures.indices @@ -886,7 +882,7 @@ public open class DoubleTensorAlgebra : return Triple(uTensor.transpose(), sTensor, vTensor.transpose()) } - override fun StructureND.symEig(): Pair = symEigJacobi(epsilon = 1e-10) + override fun StructureND.symEig(): Pair = symEigJacobi(maxIteration = 50, epsilon = 1e-15) /** * Returns eigenvalues and eigenvectors of a real symmetric matrix input or a batch of real symmetric matrices, @@ -923,103 +919,139 @@ public open class DoubleTensorAlgebra : return eig to v } - public fun StructureND.symEigJacobi(epsilon: Double): Pair { + public fun StructureND.symEigJacobi(maxIteration: Int, epsilon: Double): Pair { checkSymmetric(tensor, epsilon) val size = this.dimension - val s = zeros(this.shape) + val eigenvectors = zeros(this.shape) val eigenvalues = zeros(this.shape.sliceArray(0 until size - 1)) var eigenvalueStart = 0 var eigenvectorStart = 0 for (matrix in tensor.matrixSequence()) { - matrix.as2D().jacobiHelper(eigenvalues, s, eigenvalueStart, eigenvectorStart, epsilon) + 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] + } + } + + for (i in 0 until matrix2D.rowNum) { + eigenvalues.mutableBuffer.array()[eigenvalueStart + i] = d[i] + } + eigenvalueStart += this.shape.last() eigenvectorStart += this.shape.last() * this.shape.last() } - // TODO sort eigenvalues - return eigenvalues to s + return eigenvalues to eigenvectors } private fun MutableStructure2D.jacobiHelper( - eigenvalues: DoubleTensor, - eigenvectors: DoubleTensor, - eigenvalueStart: Int, - eigenvectorStart: Int, + maxIteration: Int, epsilon: Double - ) { - var d = this - var s = eye(this.shape[0]) + ): Pair, Structure2D> { + val A_ = this.copy().as2D() + val V = eye(this.shape[0]).as2D() + val D = DoubleTensor(intArrayOf(this.shape[0]), (0 until this.rowNum).map { this[it, it] }.toDoubleArray()).as1D() + val B = DoubleTensor(intArrayOf(this.shape[0]), (0 until this.rowNum).map { this[it, it] }.toDoubleArray()).as1D() + val Z = zeros(intArrayOf(this.shape[0])).as1D() - // TODO implement cyclic method - do { - // 1. Find max element by abs value that is not on diagonal + fun maxOffDiagonal(matrix: MutableStructure2D): Double { var maxOffDiagonalElement = 0.0 - var maxElementIndex = Pair(0, 0) - for (i in 0 until this.rowNum) { - for (j in 0 until this.colNum) { - if (i == j) continue - if (abs(d[i, j]) > maxOffDiagonalElement) { - maxOffDiagonalElement = abs(d[i, j]) - maxElementIndex = i to j + for (i in 0 until matrix.rowNum - 1) { + for (j in i + 1 until matrix.colNum) { + maxOffDiagonalElement = max(maxOffDiagonalElement, abs(matrix[i, j])) + } + } + return maxOffDiagonalElement + } + + fun rotate(a: MutableStructure2D, s: Double, tau: Double, i: Int, j: Int, k: Int, l: Int) { + val g = a[i, j] + val h = a[k, l] + a[i, j] = g - s * (h + g * tau) + a[k, l] = h + s * (g - h * tau) + } + + fun jacobiIteration( + a: MutableStructure2D, + v: MutableStructure2D, + d: MutableStructure1D, + z: MutableStructure1D, + ) { + for (ip in 0 until a.rowNum - 1) { + for (iq in ip + 1 until a.colNum) { + val g = 100.0 * abs(a[ip, iq]) + + if (g <= epsilon * abs(d[ip]) && g <= epsilon * abs(d[iq])) { + a[ip, iq] = 0.0 + continue + } + + var h = d[iq] - d[ip] + val t = when { + g <= epsilon * abs(h) -> (a[ip, iq]) / h + else -> { + val theta = 0.5 * h / (a[ip, iq]) + val denominator = abs(theta) + sqrt(1.0 + theta * theta) + if (theta < 0.0) -1.0 / denominator else 1.0 / denominator + } + } + + val c = 1.0 / sqrt(1 + t * t) + val s = t * c + val tau = s / (1.0 + c) + h = t * a[ip, iq] + z[ip] -= h + z[iq] += h + d[ip] -= h + d[iq] += h + a[ip, iq] = 0.0 + + for (j in 0 until ip) { + rotate(a, s, tau, j, ip, j, iq) + } + for (j in (ip + 1) until iq) { + rotate(a, s, tau, ip, j, j, iq) + } + for (j in (iq + 1) until a.rowNum) { + rotate(a, s, tau, ip, j, iq, j) + } + for (j in 0 until a.rowNum) { + rotate(v, s, tau, j, ip, j, iq) } } } + } - // 2. Evaluate "rotation" angle - val dIJ = d[maxElementIndex.first, maxElementIndex.second] - val dII = d[maxElementIndex.first, maxElementIndex.first] - val dJJ = d[maxElementIndex.second, maxElementIndex.second] - - val angle = if (dII == dJJ) { - if (dIJ > 0) PI / 4 else -PI / 4 - } else { - 0.5 * atan(2 * dIJ / (dJJ - dII)) - } - - // 3. Build rotation tensor `s1` - val s1 = eye(this.rowNum) - for (i in 0 until this.rowNum) { - for (j in 0 until this.colNum) { - s1.mutableBuffer.array()[i * this.rowNum + j] = when { - maxElementIndex.first == i && maxElementIndex.first == j -> cos(angle) - maxElementIndex.second == i && maxElementIndex.second == j -> cos(angle) - maxElementIndex.first == i && maxElementIndex.second == j -> sin(angle) - maxElementIndex.first == j && maxElementIndex.second == i -> -sin(angle) - else -> s1.mutableBuffer.array()[i * this.rowNum + j] - } - } - } - - // 4. Evaluate new tensor - d = ((s1.transpose() dot d) dot s1).as2D() - s = s dot s1 - if (d.isDiagonal(epsilon)) break - } while(true) - - // 5. Copy result - for (i in 0 until this.rowNum) { - for (j in 0 until this.colNum) { - eigenvectors.mutableBuffer.array()[eigenvectorStart + i * this.rowNum + j] = s.mutableBuffer.array()[i * this.rowNum + j] + fun updateDiagonal( + d: MutableStructure1D, + z: MutableStructure1D, + b: MutableStructure1D, + ) { + for (ip in 0 until d.size) { + b[ip] += z[ip] + d[ip] = b[ip] + z[ip] = 0.0 } } - for (i in 0 until this.rowNum) { - eigenvalues.mutableBuffer.array()[eigenvalueStart + i] = d[i, i] - } - } - - public fun Structure2D.isDiagonal(epsilon: Double = 1e-9): Boolean { - for (i in 0 until this.rowNum) { - for (j in 0 until this.colNum) { - if (i != j && abs(this[i, j]) > epsilon) { - return false - } + var sm = maxOffDiagonal(A_) + for (iteration in 0 until maxIteration) { + if (sm < epsilon) { + break } + + jacobiIteration(A_, V, D, Z) + updateDiagonal(D, Z, B) + sm = maxOffDiagonal(A_) } - return true + // TODO sort eigenvalues + return D to V } /** From dda6602ed4c2b8b9d7eea017b2786697e6b02ec8 Mon Sep 17 00:00:00 2001 From: Ivan Kylchik Date: Sun, 20 Feb 2022 02:44:27 +0300 Subject: [PATCH 4/5] Replace complex access to tensor with faster access to buffer in Jacobi algorithm --- .../kmath/tensors/core/DoubleTensorAlgebra.kt | 42 ++++++++++++------- 1 file changed, 26 insertions(+), 16 deletions(-) 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 675ff4191..e9dc34748 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 @@ -953,23 +953,33 @@ public open class DoubleTensorAlgebra : maxIteration: Int, epsilon: Double ): Pair, Structure2D> { - val A_ = this.copy().as2D() - val V = eye(this.shape[0]).as2D() - val D = DoubleTensor(intArrayOf(this.shape[0]), (0 until this.rowNum).map { this[it, it] }.toDoubleArray()).as1D() - val B = DoubleTensor(intArrayOf(this.shape[0]), (0 until this.rowNum).map { this[it, it] }.toDoubleArray()).as1D() - val Z = zeros(intArrayOf(this.shape[0])).as1D() + val n = this.shape[0] + val A_ = this.copy() + val V = eye(n) + val D = DoubleTensor(intArrayOf(n), (0 until this.rowNum).map { this[it, it] }.toDoubleArray()).as1D() + val B = DoubleTensor(intArrayOf(n), (0 until this.rowNum).map { this[it, it] }.toDoubleArray()).as1D() + val Z = zeros(intArrayOf(n)).as1D() - fun maxOffDiagonal(matrix: MutableStructure2D): Double { + // 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] + } + + operator fun BufferedTensor.set(i: Int, j: Int, value: Double) { + this.mutableBuffer.array()[bufferStart + i * this.shape[0] + j] = value + } + + fun maxOffDiagonal(matrix: BufferedTensor): Double { var maxOffDiagonalElement = 0.0 - for (i in 0 until matrix.rowNum - 1) { - for (j in i + 1 until matrix.colNum) { + for (i in 0 until n - 1) { + for (j in i + 1 until n) { maxOffDiagonalElement = max(maxOffDiagonalElement, abs(matrix[i, j])) } } return maxOffDiagonalElement } - fun rotate(a: MutableStructure2D, s: Double, tau: Double, i: Int, j: Int, k: Int, l: Int) { + fun rotate(a: BufferedTensor, s: Double, tau: Double, i: Int, j: Int, k: Int, l: Int) { val g = a[i, j] val h = a[k, l] a[i, j] = g - s * (h + g * tau) @@ -977,13 +987,13 @@ public open class DoubleTensorAlgebra : } fun jacobiIteration( - a: MutableStructure2D, - v: MutableStructure2D, + a: BufferedTensor, + v: BufferedTensor, d: MutableStructure1D, z: MutableStructure1D, ) { - for (ip in 0 until a.rowNum - 1) { - for (iq in ip + 1 until a.colNum) { + for (ip in 0 until n - 1) { + for (iq in ip + 1 until n) { val g = 100.0 * abs(a[ip, iq]) if (g <= epsilon * abs(d[ip]) && g <= epsilon * abs(d[iq])) { @@ -1017,10 +1027,10 @@ public open class DoubleTensorAlgebra : for (j in (ip + 1) until iq) { rotate(a, s, tau, ip, j, j, iq) } - for (j in (iq + 1) until a.rowNum) { + for (j in (iq + 1) until n) { rotate(a, s, tau, ip, j, iq, j) } - for (j in 0 until a.rowNum) { + for (j in 0 until n) { rotate(v, s, tau, j, ip, j, iq) } } @@ -1051,7 +1061,7 @@ public open class DoubleTensorAlgebra : } // TODO sort eigenvalues - return D to V + return D to V.as2D() } /** From a621dd7c5b97fa3753ec3fe7f55b238b914fd544 Mon Sep 17 00:00:00 2001 From: Ivan Kylchik Date: Sun, 20 Feb 2022 02:55:37 +0300 Subject: [PATCH 5/5] Drop duplicate test from `DorBenchmark` --- .../kscience/kmath/benchmarks/DotBenchmark.kt | 14 ++------------ 1 file changed, 2 insertions(+), 12 deletions(-) diff --git a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt index 4a5cd4aa2..16fd544a8 100644 --- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt +++ b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt @@ -15,10 +15,8 @@ import space.kscience.kmath.linear.invoke import space.kscience.kmath.linear.linearSpace import space.kscience.kmath.multik.multikAlgebra import space.kscience.kmath.operations.DoubleField -import space.kscience.kmath.operations.invoke import space.kscience.kmath.structures.Buffer import space.kscience.kmath.tensorflow.produceWithTF -import space.kscience.kmath.tensors.core.DoubleTensorAlgebra import space.kscience.kmath.tensors.core.tensorAlgebra import kotlin.random.Random @@ -36,9 +34,6 @@ internal class DotBenchmark { random.nextDouble() } - val tensor1 = DoubleTensorAlgebra.randomNormal(shape = intArrayOf(dim, dim), 12224) - val tensor2 = DoubleTensorAlgebra.randomNormal(shape = intArrayOf(dim, dim), 12225) - val cmMatrix1 = CMLinearSpace { matrix1.toCM() } val cmMatrix2 = CMLinearSpace { matrix2.toCM() } @@ -48,10 +43,10 @@ internal class DotBenchmark { @Benchmark - fun tfDot(blackhole: Blackhole){ + fun tfDot(blackhole: Blackhole) { blackhole.consume( DoubleField.produceWithTF { - tensor1 dot tensor2 + matrix1 dot matrix1 } ) } @@ -95,9 +90,4 @@ internal class DotBenchmark { fun doubleDot(blackhole: Blackhole) = with(DoubleField.linearSpace) { blackhole.consume(matrix1 dot matrix2) } - - @Benchmark - fun doubleTensorDot(blackhole: Blackhole) = DoubleTensorAlgebra.invoke { - blackhole.consume(tensor1 dot tensor2) - } }