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

SVD implementation #271

Merged
AlyaNovikova merged 2 commits from feature/tensor-algebra into feature/tensor-algebra 2021-04-06 11:00:27 +03:00
Conversation 0 Commits 2 Files Changed 4 +182 -2
4 changed files with 182 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

77
kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt
View File

@ -3,6 +3,8 @@ package space.kscience.kmath.tensors.core
import space.kscience.kmath.tensors.LinearOpsTensorAlgebra
import space.kscience.kmath.nd.as1D
import space.kscience.kmath.nd.as2D
import kotlin.math.abs
import kotlin.math.min
public class DoubleLinearOpsTensorAlgebra :
LinearOpsTensorAlgebra<Double, DoubleTensor, IntTensor>,
@ -89,8 +91,81 @@ 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<DoubleTensor, DoubleTensor, DoubleTensor> {
TODO("ALYA")
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 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())))) {
val res = ArrayList<Triple<Double, DoubleTensor, DoubleTensor>>(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())
}
return Triple(resU, resS, resV.transpose(size - 2, size - 1))
}
override fun DoubleTensor.symEig(eigenvectors: Boolean): Pair<DoubleTensor, DoubleTensor> {

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

@ -283,7 +283,47 @@ public open class DoubleTensorAlgebra : TensorPartialDivisionAlgebra<Double, Dou
}
override fun diagonalEmbedding(diagonalEntries: DoubleTensor, offset: Int, dim1: Int, dim2: Int): DoubleTensor {
TODO("Alya")
val n = diagonalEntries.shape.size
if (dim1 == dim2) {
throw RuntimeException("Diagonal dimensions cannot be identical $dim1, $dim2")
}
if (dim1 > n || dim2 > n) {
throw RuntimeException("Dimension out of range")
}
var lessDim = dim1
var greaterDim = dim2
var realOffset = offset
if (lessDim > greaterDim) {
realOffset *= -1
lessDim = greaterDim.also {greaterDim = lessDim}
}
val resShape = diagonalEntries.shape.slice(0 until lessDim).toIntArray() +
intArrayOf(diagonalEntries.shape[n - 1] + abs(realOffset)) +
diagonalEntries.shape.slice(lessDim until greaterDim - 1).toIntArray() +
intArrayOf(diagonalEntries.shape[n - 1] + abs(realOffset)) +
diagonalEntries.shape.slice(greaterDim - 1 until n - 1).toIntArray()
val resTensor = zeros(resShape)
for (i in 0 until diagonalEntries.linearStructure.size) {
val multiIndex = diagonalEntries.linearStructure.index(i)
var offset1 = 0
var offset2 = abs(realOffset)
if (realOffset < 0) {
offset1 = offset2.also {offset2 = offset1}
}
val diagonalMultiIndex = multiIndex.slice(0 until lessDim).toIntArray() +
intArrayOf(multiIndex[n - 1] + offset1) +
multiIndex.slice(lessDim until greaterDim - 1).toIntArray() +
intArrayOf(multiIndex[n - 1] + offset2) +
multiIndex.slice(greaterDim - 1 until n - 1).toIntArray()
resTensor[diagonalMultiIndex] = diagonalEntries[multiIndex]
}
return resTensor
}

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

@ -122,4 +122,36 @@ class TestDoubleLinearOpsTensorAlgebra {
assertTrue { p.dot(tensor).buffer.array().epsEqual(l.dot(u).buffer.array()) }
}
@Test
fun svd1d() = DoubleLinearOpsTensorAlgebra {
val tensor2 = fromArray(intArrayOf(2, 3), doubleArrayOf(1.0, 2.0, 3.0, 4.0, 5.0, 6.0))
val res = svd1d(tensor2)
assertTrue(res.shape contentEquals intArrayOf(2))
assertTrue { abs(abs(res.buffer.array()[res.bufferStart]) - 0.386) < 0.01}
assertTrue { abs(abs(res.buffer.array()[res.bufferStart + 1]) - 0.922) < 0.01}
}
@Test
fun svd() = DoubleLinearOpsTensorAlgebra {
val epsilon = 1e-10
fun test_tensor(tensor: DoubleTensor) {
val svd = tensor.svd()
val tensorSVD = svd.first
.dot(
diagonalEmbedding(svd.second, 0, 0, 1)
.dot(svd.third.transpose(0, 1))
)
for ((x1, x2) in tensor.buffer.array() zip tensorSVD.buffer.array()) {
assertTrue { abs(x1 - x2) < epsilon }
}
}
test_tensor(fromArray(intArrayOf(2, 3), doubleArrayOf(1.0, 2.0, 3.0, 4.0, 5.0, 6.0)))
test_tensor(fromArray(intArrayOf(2, 2), doubleArrayOf(-1.0, 0.0, 239.0, 238.0)))
}
}

33
kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleTensorAlgebra.kt
View File

@ -115,6 +115,39 @@ class TestDoubleTensorAlgebra {
assertTrue(tensor4.dot(tensor5).shape contentEquals intArrayOf(5, 4, 2, 8, 3, 8, 5))
}
@Test
fun diagonalEmbedding() = DoubleTensorAlgebra {
val tensor1 = fromArray(intArrayOf(3), doubleArrayOf(10.0, 20.0, 30.0))
val tensor2 = fromArray(intArrayOf(2, 3), doubleArrayOf(1.0, 2.0, 3.0, 4.0, 5.0, 6.0))
val tensor3 = zeros(intArrayOf(2, 3, 4, 5))
assertTrue(diagonalEmbedding(tensor3, 0, 3, 4).shape contentEquals
intArrayOf(2, 3, 4, 5, 5))
assertTrue(diagonalEmbedding(tensor3, 1, 3, 4).shape contentEquals
intArrayOf(2, 3, 4, 6, 6))
assertTrue(diagonalEmbedding(tensor3, 2, 0, 3).shape contentEquals
intArrayOf(7, 2, 3, 7, 4))
val diagonal1 = diagonalEmbedding(tensor1, 0, 1, 0)
assertTrue(diagonal1.shape contentEquals intArrayOf(3, 3))
assertTrue(diagonal1.buffer.array() contentEquals
doubleArrayOf(10.0, 0.0, 0.0, 0.0, 20.0, 0.0, 0.0, 0.0, 30.0))
val diagonal1_offset = diagonalEmbedding(tensor1, 1, 1, 0)
assertTrue(diagonal1_offset.shape contentEquals intArrayOf(4, 4))
assertTrue(diagonal1_offset.buffer.array() contentEquals
doubleArrayOf(0.0, 0.0, 0.0, 0.0, 10.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, 0.0, 0.0, 0.0, 30.0, 0.0))
val diagonal2 = diagonalEmbedding(tensor2, 1, 0, 2)
assertTrue(diagonal2.shape contentEquals intArrayOf(4, 2, 4))
assertTrue(diagonal2.buffer.array() contentEquals
doubleArrayOf(
0.0, 1.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0,
0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 5.0, 0.0,
0.0, 0.0, 0.0, 3.0, 0.0, 0.0, 0.0, 6.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0))
}
@Test
fun testContentEqual() = DoubleTensorAlgebra {
//TODO()