kscience/kmath
9
0
Fork 1
Code Issues Pull Requests 6 Packages Projects Releases 16 Wiki Activity

Golub-Kahan SVD algorithm for KMP tensors #499

Closed
grinisrit wants to merge 64 commits from dev into dev
pull from: dev
merge into: kscience:dev
Conversation 1 Commits 64 Files Changed 9 +54 -48
3 changed files with 54 additions and 48 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files
Showing only changes of commit 9cc9d959c4 - Show all commits

28
kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BroadcastDoubleTensorAlgebra.kt
View File

@ -27,7 +27,7 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() {
val newThis = broadcast[0] val newThis = broadcast[0]
val newOther = broadcast[1] val newOther = broadcast[1]
val resBuffer = DoubleArray(newThis.indices.linearSize) { i -> 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) return DoubleTensor(newThis.shape, resBuffer)
} }
@ -35,8 +35,8 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() {
override fun Tensor<Double>.plusAssign(arg: StructureND<Double>) { override fun Tensor<Double>.plusAssign(arg: StructureND<Double>) {
val newOther = broadcastTo(arg.tensor, tensor.shape) val newOther = broadcastTo(arg.tensor, tensor.shape)
for (i in 0 until tensor.indices.linearSize) { for (i in 0 until tensor.indices.linearSize) {
tensor.mutableBuffer.array()[tensor.bufferStart + i] += tensor.mutableBuffer[tensor.bufferStart + i] +=
newOther.mutableBuffer.array()[tensor.bufferStart + i] newOther.mutableBuffer[tensor.bufferStart + i]
} }
} }
@ -45,7 +45,7 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() {
val newThis = broadcast[0] val newThis = broadcast[0]
val newOther = broadcast[1] val newOther = broadcast[1]
val resBuffer = DoubleArray(newThis.indices.linearSize) { i -> 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) return DoubleTensor(newThis.shape, resBuffer)
} }
@ -53,8 +53,8 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() {
override fun Tensor<Double>.minusAssign(arg: StructureND<Double>) { override fun Tensor<Double>.minusAssign(arg: StructureND<Double>) {
val newOther = broadcastTo(arg.tensor, tensor.shape) val newOther = broadcastTo(arg.tensor, tensor.shape)
for (i in 0 until tensor.indices.linearSize) { for (i in 0 until tensor.indices.linearSize) {
tensor.mutableBuffer.array()[tensor.bufferStart + i] -= tensor.mutableBuffer[tensor.bufferStart + i] -=
newOther.mutableBuffer.array()[tensor.bufferStart + i] newOther.mutableBuffer[tensor.bufferStart + i]
} }
} }
@ -63,8 +63,8 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() {
val newThis = broadcast[0] val newThis = broadcast[0]
val newOther = broadcast[1] val newOther = broadcast[1]
val resBuffer = DoubleArray(newThis.indices.linearSize) { i -> val resBuffer = DoubleArray(newThis.indices.linearSize) { i ->
newThis.mutableBuffer.array()[newThis.bufferStart + i] * newThis.mutableBuffer[newThis.bufferStart + i] *
newOther.mutableBuffer.array()[newOther.bufferStart + i] newOther.mutableBuffer[newOther.bufferStart + i]
} }
return DoubleTensor(newThis.shape, resBuffer) return DoubleTensor(newThis.shape, resBuffer)
} }
@ -72,8 +72,8 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() {
override fun Tensor<Double>.timesAssign(arg: StructureND<Double>) { override fun Tensor<Double>.timesAssign(arg: StructureND<Double>) {
val newOther = broadcastTo(arg.tensor, tensor.shape) val newOther = broadcastTo(arg.tensor, tensor.shape)
for (i in 0 until tensor.indices.linearSize) { for (i in 0 until tensor.indices.linearSize) {
tensor.mutableBuffer.array()[tensor.bufferStart + i] *= tensor.mutableBuffer[tensor.bufferStart + i] *=
newOther.mutableBuffer.array()[tensor.bufferStart + i] newOther.mutableBuffer[tensor.bufferStart + i]
} }
} }
@ -82,8 +82,8 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() {
val newThis = broadcast[0] val newThis = broadcast[0]
val newOther = broadcast[1] val newOther = broadcast[1]
val resBuffer = DoubleArray(newThis.indices.linearSize) { i -> val resBuffer = DoubleArray(newThis.indices.linearSize) { i ->
newThis.mutableBuffer.array()[newOther.bufferStart + i] / newThis.mutableBuffer[newOther.bufferStart + i] /
newOther.mutableBuffer.array()[newOther.bufferStart + i] newOther.mutableBuffer[newOther.bufferStart + i]
} }
return DoubleTensor(newThis.shape, resBuffer) return DoubleTensor(newThis.shape, resBuffer)
} }
@ -91,8 +91,8 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() {
override fun Tensor<Double>.divAssign(arg: StructureND<Double>) { override fun Tensor<Double>.divAssign(arg: StructureND<Double>) {
val newOther = broadcastTo(arg.tensor, tensor.shape) val newOther = broadcastTo(arg.tensor, tensor.shape)
for (i in 0 until tensor.indices.linearSize) { for (i in 0 until tensor.indices.linearSize) {
tensor.mutableBuffer.array()[tensor.bufferStart + i] /= tensor.mutableBuffer[tensor.bufferStart + i] /=
newOther.mutableBuffer.array()[tensor.bufferStart + i] newOther.mutableBuffer[tensor.bufferStart + i]
} }
} }
} }

67
kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt
View File

@ -89,7 +89,7 @@ public open class DoubleTensorAlgebra :
} }
override fun StructureND<Double>.valueOrNull(): Double? = if (tensor.shape contentEquals intArrayOf(1)) override fun StructureND<Double>.valueOrNull(): Double? = if (tensor.shape contentEquals intArrayOf(1))
tensor.mutableBuffer.array()[tensor.bufferStart] else null tensor.mutableBuffer[tensor.bufferStart] else null
override fun StructureND<Double>.value(): Double = valueOrNull() override fun StructureND<Double>.value(): Double = valueOrNull()
?: throw IllegalArgumentException("The tensor shape is $shape, but value method is allowed only for shape [1]") ?: 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<Double>): DoubleTensor { override fun Double.plus(arg: StructureND<Double>): DoubleTensor {
val resBuffer = DoubleArray(arg.tensor.numElements) { i -> 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) return DoubleTensor(arg.shape, resBuffer)
} }
@ -218,35 +218,35 @@ public open class DoubleTensorAlgebra :
override fun StructureND<Double>.plus(arg: StructureND<Double>): DoubleTensor { override fun StructureND<Double>.plus(arg: StructureND<Double>): DoubleTensor {
checkShapesCompatible(tensor, arg.tensor) checkShapesCompatible(tensor, arg.tensor)
val resBuffer = DoubleArray(tensor.numElements) { i -> 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) return DoubleTensor(tensor.shape, resBuffer)
} }
override fun Tensor<Double>.plusAssign(value: Double) { override fun Tensor<Double>.plusAssign(value: Double) {
for (i in 0 until tensor.numElements) { for (i in 0 until tensor.numElements) {
tensor.mutableBuffer.array()[tensor.bufferStart + i] += value tensor.mutableBuffer[tensor.bufferStart + i] += value
} }
} }
override fun Tensor<Double>.plusAssign(arg: StructureND<Double>) { override fun Tensor<Double>.plusAssign(arg: StructureND<Double>) {
checkShapesCompatible(tensor, arg.tensor) checkShapesCompatible(tensor, arg.tensor)
for (i in 0 until tensor.numElements) { for (i in 0 until tensor.numElements) {
tensor.mutableBuffer.array()[tensor.bufferStart + i] += tensor.mutableBuffer[tensor.bufferStart + i] +=
arg.tensor.mutableBuffer.array()[tensor.bufferStart + i] arg.tensor.mutableBuffer[tensor.bufferStart + i]
} }
} }
override fun Double.minus(arg: StructureND<Double>): DoubleTensor { override fun Double.minus(arg: StructureND<Double>): DoubleTensor {
val resBuffer = DoubleArray(arg.tensor.numElements) { i -> 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) return DoubleTensor(arg.shape, resBuffer)
} }
override fun StructureND<Double>.minus(arg: Double): DoubleTensor { override fun StructureND<Double>.minus(arg: Double): DoubleTensor {
val resBuffer = DoubleArray(tensor.numElements) { i -> val resBuffer = DoubleArray(tensor.numElements) { i ->
tensor.mutableBuffer.array()[tensor.bufferStart + i] - arg tensor.mutableBuffer[tensor.bufferStart + i] - arg
} }
return DoubleTensor(tensor.shape, resBuffer) return DoubleTensor(tensor.shape, resBuffer)
} }
@ -254,28 +254,28 @@ public open class DoubleTensorAlgebra :
override fun StructureND<Double>.minus(arg: StructureND<Double>): DoubleTensor { override fun StructureND<Double>.minus(arg: StructureND<Double>): DoubleTensor {
checkShapesCompatible(tensor, arg) checkShapesCompatible(tensor, arg)
val resBuffer = DoubleArray(tensor.numElements) { i -> 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) return DoubleTensor(tensor.shape, resBuffer)
} }
override fun Tensor<Double>.minusAssign(value: Double) { override fun Tensor<Double>.minusAssign(value: Double) {
for (i in 0 until tensor.numElements) { for (i in 0 until tensor.numElements) {
tensor.mutableBuffer.array()[tensor.bufferStart + i] -= value tensor.mutableBuffer[tensor.bufferStart + i] -= value
} }
} }
override fun Tensor<Double>.minusAssign(arg: StructureND<Double>) { override fun Tensor<Double>.minusAssign(arg: StructureND<Double>) {
checkShapesCompatible(tensor, arg) checkShapesCompatible(tensor, arg)
for (i in 0 until tensor.numElements) { for (i in 0 until tensor.numElements) {
tensor.mutableBuffer.array()[tensor.bufferStart + i] -= tensor.mutableBuffer[tensor.bufferStart + i] -=
arg.tensor.mutableBuffer.array()[tensor.bufferStart + i] arg.tensor.mutableBuffer[tensor.bufferStart + i]
} }
} }
override fun Double.times(arg: StructureND<Double>): DoubleTensor { override fun Double.times(arg: StructureND<Double>): DoubleTensor {
val resBuffer = DoubleArray(arg.tensor.numElements) { i -> 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) return DoubleTensor(arg.shape, resBuffer)
} }
@ -285,36 +285,36 @@ public open class DoubleTensorAlgebra :
override fun StructureND<Double>.times(arg: StructureND<Double>): DoubleTensor { override fun StructureND<Double>.times(arg: StructureND<Double>): DoubleTensor {
checkShapesCompatible(tensor, arg) checkShapesCompatible(tensor, arg)
val resBuffer = DoubleArray(tensor.numElements) { i -> val resBuffer = DoubleArray(tensor.numElements) { i ->
tensor.mutableBuffer.array()[tensor.bufferStart + i] * tensor.mutableBuffer[tensor.bufferStart + i] *
arg.tensor.mutableBuffer.array()[arg.tensor.bufferStart + i] arg.tensor.mutableBuffer[arg.tensor.bufferStart + i]
} }
return DoubleTensor(tensor.shape, resBuffer) return DoubleTensor(tensor.shape, resBuffer)
} }
override fun Tensor<Double>.timesAssign(value: Double) { override fun Tensor<Double>.timesAssign(value: Double) {
for (i in 0 until tensor.numElements) { for (i in 0 until tensor.numElements) {
tensor.mutableBuffer.array()[tensor.bufferStart + i] *= value tensor.mutableBuffer[tensor.bufferStart + i] *= value
} }
} }
override fun Tensor<Double>.timesAssign(arg: StructureND<Double>) { override fun Tensor<Double>.timesAssign(arg: StructureND<Double>) {
checkShapesCompatible(tensor, arg) checkShapesCompatible(tensor, arg)
for (i in 0 until tensor.numElements) { for (i in 0 until tensor.numElements) {
tensor.mutableBuffer.array()[tensor.bufferStart + i] *= tensor.mutableBuffer[tensor.bufferStart + i] *=
arg.tensor.mutableBuffer.array()[tensor.bufferStart + i] arg.tensor.mutableBuffer[tensor.bufferStart + i]
} }
} }
override fun Double.div(arg: StructureND<Double>): DoubleTensor { override fun Double.div(arg: StructureND<Double>): DoubleTensor {
val resBuffer = DoubleArray(arg.tensor.numElements) { i -> 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) return DoubleTensor(arg.shape, resBuffer)
} }
override fun StructureND<Double>.div(arg: Double): DoubleTensor { override fun StructureND<Double>.div(arg: Double): DoubleTensor {
val resBuffer = DoubleArray(tensor.numElements) { i -> val resBuffer = DoubleArray(tensor.numElements) { i ->
tensor.mutableBuffer.array()[tensor.bufferStart + i] / arg tensor.mutableBuffer[tensor.bufferStart + i] / arg
} }
return DoubleTensor(shape, resBuffer) return DoubleTensor(shape, resBuffer)
} }
@ -322,29 +322,29 @@ public open class DoubleTensorAlgebra :
override fun StructureND<Double>.div(arg: StructureND<Double>): DoubleTensor { override fun StructureND<Double>.div(arg: StructureND<Double>): DoubleTensor {
checkShapesCompatible(tensor, arg) checkShapesCompatible(tensor, arg)
val resBuffer = DoubleArray(tensor.numElements) { i -> val resBuffer = DoubleArray(tensor.numElements) { i ->
tensor.mutableBuffer.array()[arg.tensor.bufferStart + i] / tensor.mutableBuffer[arg.tensor.bufferStart + i] /
arg.tensor.mutableBuffer.array()[arg.tensor.bufferStart + i] arg.tensor.mutableBuffer[arg.tensor.bufferStart + i]
} }
return DoubleTensor(tensor.shape, resBuffer) return DoubleTensor(tensor.shape, resBuffer)
} }
override fun Tensor<Double>.divAssign(value: Double) { override fun Tensor<Double>.divAssign(value: Double) {
for (i in 0 until tensor.numElements) { for (i in 0 until tensor.numElements) {
tensor.mutableBuffer.array()[tensor.bufferStart + i] /= value tensor.mutableBuffer[tensor.bufferStart + i] /= value
} }
} }
override fun Tensor<Double>.divAssign(arg: StructureND<Double>) { override fun Tensor<Double>.divAssign(arg: StructureND<Double>) {
checkShapesCompatible(tensor, arg) checkShapesCompatible(tensor, arg)
for (i in 0 until tensor.numElements) { for (i in 0 until tensor.numElements) {
tensor.mutableBuffer.array()[tensor.bufferStart + i] /= tensor.mutableBuffer[tensor.bufferStart + i] /=
arg.tensor.mutableBuffer.array()[tensor.bufferStart + i] arg.tensor.mutableBuffer[tensor.bufferStart + i]
} }
} }
override fun StructureND<Double>.unaryMinus(): DoubleTensor { override fun StructureND<Double>.unaryMinus(): DoubleTensor {
val resBuffer = DoubleArray(tensor.numElements) { i -> val resBuffer = DoubleArray(tensor.numElements) { i ->
tensor.mutableBuffer.array()[tensor.bufferStart + i].unaryMinus() tensor.mutableBuffer[tensor.bufferStart + i].unaryMinus()
} }
return DoubleTensor(tensor.shape, resBuffer) return DoubleTensor(tensor.shape, resBuffer)
} }
@ -367,8 +367,8 @@ public open class DoubleTensorAlgebra :
newMultiIndex[ii] = newMultiIndex[jj].also { newMultiIndex[jj] = newMultiIndex[ii] } newMultiIndex[ii] = newMultiIndex[jj].also { newMultiIndex[jj] = newMultiIndex[ii] }
val linearIndex = resTensor.indices.offset(newMultiIndex) val linearIndex = resTensor.indices.offset(newMultiIndex)
resTensor.mutableBuffer.array()[linearIndex] = resTensor.mutableBuffer[linearIndex] =
tensor.mutableBuffer.array()[tensor.bufferStart + offset] tensor.mutableBuffer[tensor.bufferStart + offset]
} }
return resTensor return resTensor
} }
@ -958,18 +958,21 @@ public open class DoubleTensorAlgebra :
var eigenvalueStart = 0 var eigenvalueStart = 0
var eigenvectorStart = 0 var eigenvectorStart = 0
val eigenvaluesBuffer = eigenvalues.mutableBuffer
val eigenvectorsBuffer = eigenvectors.mutableBuffer
for (matrix in tensor.matrixSequence()) { for (matrix in tensor.matrixSequence()) {
val matrix2D = matrix.as2D() val matrix2D = matrix.as2D()
val (d, v) = matrix2D.jacobiHelper(maxIteration, epsilon) val (d, v) = matrix2D.jacobiHelper(maxIteration, epsilon)
for (i in 0 until matrix2D.rowNum) { for (i in 0 until matrix2D.rowNum) {
for (j in 0 until matrix2D.colNum) { 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) { for (i in 0 until matrix2D.rowNum) {
eigenvalues.mutableBuffer.array()[eigenvalueStart + i] = d[i] eigenvaluesBuffer[eigenvalueStart + i] = d[i]
} }
eigenvalueStart += this.shape.last() eigenvalueStart += this.shape.last()
@ -992,11 +995,11 @@ public open class DoubleTensorAlgebra :
// assume that buffered tensor is square matrix // assume that buffered tensor is square matrix
operator fun BufferedTensor<Double>.get(i: Int, j: Int): Double { operator fun BufferedTensor<Double>.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<Double>.set(i: Int, j: Int, value: Double) { operator fun BufferedTensor<Double>.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>): Double { fun maxOffDiagonal(matrix: BufferedTensor<Double>): Double {

7
kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt
View File

@ -156,9 +156,12 @@ internal class TestDoubleLinearOpsTensorAlgebra {
val res = svd1d(tensor2) val res = svd1d(tensor2)
val resBuffer = res.mutableBuffer
val resStart = res.bufferStart
assertTrue(res.shape contentEquals intArrayOf(2)) assertTrue(res.shape contentEquals intArrayOf(2))
assertTrue { abs(abs(res.mutableBuffer.array()[res.bufferStart]) - 0.386) < 0.01 } assertTrue { abs(abs(resBuffer[resStart]) - 0.386) < 0.01 }
assertTrue { abs(abs(res.mutableBuffer.array()[res.bufferStart + 1]) - 0.922) < 0.01 } assertTrue { abs(abs(resBuffer[resStart + 1]) - 0.922) < 0.01 }
} }
@Test @Test