From 4336788a6ba846af6a63be1352d1bec1c1f39b56 Mon Sep 17 00:00:00 2001 From: Roland Grinis Date: Tue, 6 Apr 2021 09:00:13 +0100 Subject: [PATCH] Moving Alya's SVD implementation to linutils --- .../core/DoubleLinearOpsTensorAlgebra.kt | 66 +--------------- .../kscience/kmath/tensors/core/linutils.kt | 76 ++++++++++++++++++- .../core/TestDoubleLinearOpsAlgebra.kt | 4 +- 3 files changed, 78 insertions(+), 68 deletions(-) diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt index fbd9da18e..7db812d9d 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 @@ -91,29 +91,6 @@ public class DoubleLinearOpsTensorAlgebra : return qTensor to rTensor } - internal fun svd1d(a: DoubleTensor, epsilon: Double = 1e-10): DoubleTensor { - val (n, m) = a.shape - var v: DoubleTensor - val b: DoubleTensor - if (n > m) { - b = a.transpose(0, 1).dot(a) - v = DoubleTensor(intArrayOf(m), getRandomNormals(m, 0)) - } else { - b = a.dot(a.transpose(0, 1)) - v = DoubleTensor(intArrayOf(n), getRandomNormals(n, 0)) - } - - var lastV: DoubleTensor - while (true) { - lastV = v - v = b.dot(lastV) - val norm = DoubleAnalyticTensorAlgebra { (v dot v).sqrt().value() } - v = v.times(1.0 / norm) - if (abs(v.dot(lastV).value()) > 1 - epsilon) { - return v - } - } - } override fun DoubleTensor.svd(): Triple { val size = this.shape.size @@ -124,47 +101,8 @@ public class DoubleLinearOpsTensorAlgebra : val resV = zeros(commonShape + intArrayOf(min(n, m), m)) for ((matrix, USV) in this.matrixSequence() - .zip(resU.matrixSequence().zip(resS.vectorSequence().zip(resV.matrixSequence())))) { - val res = ArrayList>(0) - val (matrixU, SV) = USV - val (matrixS, matrixV) = SV - - for (k in 0 until min(n, m)) { - var a = matrix.asTensor().copy() - for ((singularValue, u, v) in res.slice(0 until k)) { - val outerProduct = DoubleArray(u.shape[0] * v.shape[0]) - for (i in 0 until u.shape[0]) { - for (j in 0 until v.shape[0]) { - outerProduct[i * v.shape[0] + j] = u[i].value() * v[j].value() - } - } - a = a - singularValue.times(DoubleTensor(intArrayOf(u.shape[0], v.shape[0]), outerProduct)) - } - var v: DoubleTensor - var u: DoubleTensor - var norm: Double - if (n > m) { - v = svd1d(a) - u = matrix.asTensor().dot(v) - norm = DoubleAnalyticTensorAlgebra { (u dot u).sqrt().value() } - u = u.times(1.0 / norm) - } else { - u = svd1d(a) - v = matrix.asTensor().transpose(0, 1).dot(u) - norm = DoubleAnalyticTensorAlgebra { (v dot v).sqrt().value() } - v = v.times(1.0 / norm) - } - - res.add(Triple(norm, u, v)) - } - - 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()) - } + .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)) } 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 f6cf71b07..e1d5c1819 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 @@ -4,6 +4,8 @@ import space.kscience.kmath.nd.MutableStructure1D import space.kscience.kmath.nd.MutableStructure2D import space.kscience.kmath.nd.as1D import space.kscience.kmath.nd.as2D +import kotlin.math.abs +import kotlin.math.min import kotlin.math.sqrt @@ -195,8 +197,8 @@ internal inline fun DoubleLinearOpsTensorAlgebra.qrHelper( checkSquareMatrix(matrix.shape) val n = matrix.shape[0] val qM = q.as2D() - val matrixT = matrix.transpose(0,1) - val qT = q.transpose(0,1) + val matrixT = matrix.transpose(0, 1) + val qT = q.transpose(0, 1) for (j in 0 until n) { val v = matrixT[j] @@ -216,3 +218,73 @@ internal inline fun DoubleLinearOpsTensorAlgebra.qrHelper( } } } + +internal inline fun DoubleLinearOpsTensorAlgebra.svd1d(a: DoubleTensor, epsilon: Double = 1e-10): DoubleTensor { + val (n, m) = a.shape + var v: DoubleTensor + val b: DoubleTensor + if (n > m) { + b = a.transpose(0, 1).dot(a) + v = DoubleTensor(intArrayOf(m), getRandomNormals(m, 0)) + } else { + b = a.dot(a.transpose(0, 1)) + v = DoubleTensor(intArrayOf(n), getRandomNormals(n, 0)) + } + + var lastV: DoubleTensor + while (true) { + lastV = v + v = b.dot(lastV) + val norm = DoubleAnalyticTensorAlgebra { (v dot v).sqrt().value() } + v = v.times(1.0 / norm) + if (abs(v.dot(lastV).value()) > 1 - epsilon) { + return v + } + } +} + +internal inline fun DoubleLinearOpsTensorAlgebra.svdHelper( + matrix: DoubleTensor, + USV: Pair, Pair, BufferedTensor>>, + m: Int, n: Int +): Unit { + val res = ArrayList>(0) + val (matrixU, SV) = USV + val (matrixS, matrixV) = SV + + for (k in 0 until min(n, m)) { + var a = matrix.copy() + for ((singularValue, u, v) in res.slice(0 until k)) { + val outerProduct = DoubleArray(u.shape[0] * v.shape[0]) + for (i in 0 until u.shape[0]) { + for (j in 0 until v.shape[0]) { + outerProduct[i * v.shape[0] + j] = u[i].value() * v[j].value() + } + } + a = a - singularValue.times(DoubleTensor(intArrayOf(u.shape[0], v.shape[0]), outerProduct)) + } + var v: DoubleTensor + var u: DoubleTensor + var norm: Double + if (n > m) { + v = svd1d(a) + u = matrix.dot(v) + norm = DoubleAnalyticTensorAlgebra { (u dot u).sqrt().value() } + u = u.times(1.0 / norm) + } else { + u = svd1d(a) + v = matrix.transpose(0, 1).dot(u) + norm = DoubleAnalyticTensorAlgebra { (v dot v).sqrt().value() } + v = v.times(1.0 / norm) + } + + res.add(Triple(norm, u, v)) + } + + 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()) +} 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 843707153..60f7f4a97 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 @@ -124,7 +124,7 @@ class TestDoubleLinearOpsTensorAlgebra { } @Test - fun svd1d() = DoubleLinearOpsTensorAlgebra { + fun testSVD1D() = DoubleLinearOpsTensorAlgebra { val tensor2 = fromArray(intArrayOf(2, 3), doubleArrayOf(1.0, 2.0, 3.0, 4.0, 5.0, 6.0)) val res = svd1d(tensor2) @@ -135,7 +135,7 @@ class TestDoubleLinearOpsTensorAlgebra { } @Test - fun svd() = DoubleLinearOpsTensorAlgebra { + fun testSVD() = DoubleLinearOpsTensorAlgebra { val epsilon = 1e-10 fun test_tensor(tensor: DoubleTensor) { val svd = tensor.svd()