Fix bugs in svd #275
@ -96,14 +96,20 @@ public class DoubleLinearOpsTensorAlgebra :
|
|||||||
val size = this.shape.size
|
val size = this.shape.size
|
||||||
val commonShape = this.shape.sliceArray(0 until size - 2)
|
val commonShape = this.shape.sliceArray(0 until size - 2)
|
||||||
val (n, m) = this.shape.sliceArray(size - 2 until size)
|
val (n, m) = this.shape.sliceArray(size - 2 until size)
|
||||||
val resU = zeros(commonShape + intArrayOf(n, min(n, m)))
|
val resU = zeros(commonShape + intArrayOf(min(n, m), n))
|
||||||
val resS = zeros(commonShape + intArrayOf(min(n, m)))
|
val resS = zeros(commonShape + intArrayOf(min(n, m)))
|
||||||
val resV = zeros(commonShape + intArrayOf(min(n, m), m))
|
val resV = zeros(commonShape + intArrayOf(min(n, m), m))
|
||||||
|
|
||||||
for ((matrix, USV) in this.matrixSequence()
|
for ((matrix, USV) in this.matrixSequence()
|
||||||
.zip(resU.matrixSequence().zip(resS.vectorSequence().zip(resV.matrixSequence()))))
|
.zip(resU.matrixSequence().zip(resS.vectorSequence().zip(resV.matrixSequence())))) {
|
||||||
svdHelper(matrix.asTensor(), USV, m, n)
|
val size = matrix.shape.reduce { acc, i -> acc * i }
|
||||||
return Triple(resU, resS, resV.transpose(size - 2, size - 1))
|
val curMatrix = DoubleTensor(
|
||||||
|
matrix.shape,
|
||||||
|
matrix.buffer.array().slice(matrix.bufferStart until matrix.bufferStart + size).toDoubleArray()
|
||||||
|
)
|
||||||
|
svdHelper(curMatrix, USV, m, n)
|
||||||
|
}
|
||||||
|
return Triple(resU.transpose(size - 2, size - 1), resS, resV)
|
||||||
}
|
}
|
||||||
|
|
||||||
override fun DoubleTensor.symEig(eigenvectors: Boolean): Pair<DoubleTensor, DoubleTensor> {
|
override fun DoubleTensor.symEig(eigenvectors: Boolean): Pair<DoubleTensor, DoubleTensor> {
|
||||||
|
@ -225,10 +225,10 @@ internal inline fun DoubleLinearOpsTensorAlgebra.svd1d(a: DoubleTensor, epsilon:
|
|||||||
val b: DoubleTensor
|
val b: DoubleTensor
|
||||||
if (n > m) {
|
if (n > m) {
|
||||||
b = a.transpose(0, 1).dot(a)
|
b = a.transpose(0, 1).dot(a)
|
||||||
v = DoubleTensor(intArrayOf(m), getRandomNormals(m, 0))
|
v = DoubleTensor(intArrayOf(m), getRandomUnitVector(m, 0))
|
||||||
} else {
|
} else {
|
||||||
b = a.dot(a.transpose(0, 1))
|
b = a.dot(a.transpose(0, 1))
|
||||||
v = DoubleTensor(intArrayOf(n), getRandomNormals(n, 0))
|
v = DoubleTensor(intArrayOf(n), getRandomUnitVector(n, 0))
|
||||||
}
|
}
|
||||||
|
|
||||||
var lastV: DoubleTensor
|
var lastV: DoubleTensor
|
||||||
@ -284,7 +284,13 @@ internal inline fun DoubleLinearOpsTensorAlgebra.svdHelper(
|
|||||||
val s = res.map { it.first }.toDoubleArray()
|
val s = res.map { it.first }.toDoubleArray()
|
||||||
val uBuffer = res.map { it.second }.flatMap { it.buffer.array().toList() }.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()
|
val vBuffer = res.map { it.third }.flatMap { it.buffer.array().toList() }.toDoubleArray()
|
||||||
uBuffer.copyInto(matrixU.buffer.array())
|
for (i in uBuffer.indices) {
|
||||||
s.copyInto(matrixS.buffer.array())
|
matrixU.buffer.array()[matrixU.bufferStart + i] = uBuffer[i]
|
||||||
vBuffer.copyInto(matrixV.buffer.array())
|
}
|
||||||
|
for (i in s.indices) {
|
||||||
|
matrixS.buffer.array()[matrixS.bufferStart + i] = s[i]
|
||||||
|
}
|
||||||
|
for (i in vBuffer.indices) {
|
||||||
|
matrixV.buffer.array()[matrixV.bufferStart + i] = vBuffer[i]
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
@ -3,6 +3,7 @@ package space.kscience.kmath.tensors.core
|
|||||||
import space.kscience.kmath.samplers.GaussianSampler
|
import space.kscience.kmath.samplers.GaussianSampler
|
||||||
import space.kscience.kmath.stat.RandomGenerator
|
import space.kscience.kmath.stat.RandomGenerator
|
||||||
import space.kscience.kmath.structures.*
|
import space.kscience.kmath.structures.*
|
||||||
|
import kotlin.math.sqrt
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Returns a reference to [IntArray] containing all of the elements of this [Buffer].
|
* Returns a reference to [IntArray] containing all of the elements of this [Buffer].
|
||||||
@ -42,6 +43,12 @@ internal inline fun getRandomNormals(n: Int, seed: Long): DoubleArray {
|
|||||||
return distribution.sample(generator).nextBufferBlocking(n).toDoubleArray()
|
return distribution.sample(generator).nextBufferBlocking(n).toDoubleArray()
|
||||||
}
|
}
|
||||||
|
|
||||||
|
internal inline fun getRandomUnitVector(n: Int, seed: Long): DoubleArray {
|
||||||
|
val unnorm = getRandomNormals(n, seed)
|
||||||
|
val norm = sqrt(unnorm.map { it * it }.sum())
|
||||||
|
return unnorm.map { it / norm }.toDoubleArray()
|
||||||
|
}
|
||||||
|
|
||||||
internal inline fun minusIndexFrom(n: Int, i: Int) : Int = if (i >= 0) i else {
|
internal inline fun minusIndexFrom(n: Int, i: Int) : Int = if (i >= 0) i else {
|
||||||
val ii = n + i
|
val ii = n + i
|
||||||
check(ii >= 0) {
|
check(ii >= 0) {
|
||||||
|
@ -126,9 +126,9 @@ class TestDoubleLinearOpsTensorAlgebra {
|
|||||||
testSVDFor(fromArray(intArrayOf(2, 2), doubleArrayOf(-1.0, 0.0, 239.0, 238.0)))
|
testSVDFor(fromArray(intArrayOf(2, 2), doubleArrayOf(-1.0, 0.0, 239.0, 238.0)))
|
||||||
}
|
}
|
||||||
|
|
||||||
@Test @Ignore
|
@Test
|
||||||
fun testBatchedSVD() = DoubleLinearOpsTensorAlgebra {
|
fun testBatchedSVD() = DoubleLinearOpsTensorAlgebra {
|
||||||
val tensor = randNormal(intArrayOf(7, 5, 3), 0)
|
val tensor = randNormal(intArrayOf(1, 15, 4, 7, 5, 3), 0)
|
||||||
val (tensorU, tensorS, tensorV) = tensor.svd()
|
val (tensorU, tensorS, tensorV) = tensor.svd()
|
||||||
val tensorSVD = tensorU dot (diagonalEmbedding(tensorS) dot tensorV)
|
val tensorSVD = tensorU dot (diagonalEmbedding(tensorS) dot tensorV)
|
||||||
assertTrue(tensor.eq(tensorSVD))
|
assertTrue(tensor.eq(tensorSVD))
|
||||||
|
Loading…
Reference in New Issue
Block a user