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 +32 -35
2 changed files with 32 additions and 35 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files
Showing only changes of commit 9a46dd8966 - Show all commits

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

@ -833,8 +833,33 @@ public open class DoubleTensorAlgebra :
return qTensor to rTensor return qTensor to rTensor
} }
override fun StructureND<Double>.svd(): Triple<DoubleTensor, DoubleTensor, DoubleTensor> = override fun StructureND<Double>.svd(): Triple<DoubleTensor, DoubleTensor, DoubleTensor> {
svd(epsilon = 1e-10) 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(m, n))
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().svdHelper(uTensors[index].as2D(), sTensorVectors[index], vTensors[index].as2D())
}
return Triple(uTensor.transpose(), sTensor, vTensor)
}
/** /**
* Singular Value Decomposition. * Singular Value Decomposition.
@ -849,7 +874,7 @@ public open class DoubleTensorAlgebra :
* i.e., the precision with which the cosine approaches 1 in an iterative algorithm. * i.e., the precision with which the cosine approaches 1 in an iterative algorithm.
* @return a triple `Triple(U, S, V)`. * @return a triple `Triple(U, S, V)`.
*/ */
public fun StructureND<Double>.svd(epsilon: Double): Triple<DoubleTensor, DoubleTensor, DoubleTensor> { public fun StructureND<Double>.svdPowerMethod(epsilon: Double = 1e-10): Triple<DoubleTensor, DoubleTensor, DoubleTensor> {
val size = tensor.dimension val size = tensor.dimension
val commonShape = tensor.shape.sliceArray(0 until size - 2) val commonShape = tensor.shape.sliceArray(0 until size - 2)
val (n, m) = tensor.shape.sliceArray(size - 2 until size) val (n, m) = tensor.shape.sliceArray(size - 2 until size)
@ -876,40 +901,12 @@ public open class DoubleTensorAlgebra :
.slice(matrix.bufferStart until matrix.bufferStart + matrixSize) .slice(matrix.bufferStart until matrix.bufferStart + matrixSize)
.toDoubleArray() .toDoubleArray()
) )
svdHelper(curMatrix, usv, m, n, epsilon) svdPowerMethodHelper(curMatrix, usv, m, n, epsilon)
} }
return Triple(uTensor.transpose(), sTensor, vTensor.transpose()) return Triple(uTensor.transpose(), sTensor, vTensor.transpose())
} }
public fun StructureND<Double>.svdGolabKahan(): Triple<DoubleTensor, DoubleTensor, DoubleTensor> {
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(m, n))
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().svdGolabKahanHelper(uTensors[index].as2D(), sTensorVectors[index], vTensors[index].as2D())
}
return Triple(uTensor.transpose(), sTensor, vTensor)
}
override fun StructureND<Double>.symEig(): Pair<DoubleTensor, DoubleTensor> = symEigJacobi(maxIteration = 50, epsilon = 1e-15) override fun StructureND<Double>.symEig(): Pair<DoubleTensor, DoubleTensor> = symEigJacobi(maxIteration = 50, epsilon = 1e-15)
/** /**
@ -935,7 +932,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 shp = s.shape + intArrayOf(1)
val utv = u.transpose() dot v val utv = u.transpose() dot v
val n = s.shape.last() val n = s.shape.last()

4
kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt
View File

@ -300,7 +300,7 @@ internal fun DoubleTensorAlgebra.svd1d(a: DoubleTensor, epsilon: Double = 1e-10)
} }
} }
internal fun DoubleTensorAlgebra.svdHelper( internal fun DoubleTensorAlgebra.svdPowerMethodHelper(
matrix: DoubleTensor, matrix: DoubleTensor,
USV: Triple<BufferedTensor<Double>, BufferedTensor<Double>, BufferedTensor<Double>>, USV: Triple<BufferedTensor<Double>, BufferedTensor<Double>, BufferedTensor<Double>>,
m: Int, n: Int, epsilon: Double, m: Int, n: Int, epsilon: Double,
@ -372,7 +372,7 @@ private fun SIGN(a: Double, b: Double): Double {
else else
return -abs(a) return -abs(a)
} }
internal fun MutableStructure2D<Double>.svdGolabKahanHelper(u: MutableStructure2D<Double>, w: BufferedTensor<Double>, v: MutableStructure2D<Double>) { internal fun MutableStructure2D<Double>.svdHelper(u: MutableStructure2D<Double>, w: BufferedTensor<Double>, v: MutableStructure2D<Double>) {
val shape = this.shape val shape = this.shape
val m = shape.component1() val m = shape.component1()
val n = shape.component2() val n = shape.component2()