Merge pull request #486 from margarita0303/feature/tensors-performance

added partial implementation of svd calculation with divide-by-near-zero error
This commit is contained in:
Roland Grinis 2022-05-25 13:59:23 +03:00 committed by GitHub
commit 13fe078304
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
3 changed files with 411 additions and 0 deletions

View File

@ -0,0 +1,65 @@
/*
* Copyright 2018-2021 KMath contributors.
* Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
*/
package space.kscience.kmath.tensors
import space.kscience.kmath.linear.transpose
import space.kscience.kmath.misc.PerformancePitfall
import space.kscience.kmath.nd.MutableStructure2D
import space.kscience.kmath.nd.Structure2D
import space.kscience.kmath.nd.as2D
import space.kscience.kmath.tensors.core.*
import space.kscience.kmath.tensors.core.tensorAlgebra
import kotlin.math.*
fun MutableStructure2D<Double>.print() {
val n = this.shape.component1()
val m = this.shape.component2()
for (i in 0 until n) {
for (j in 0 until m) {
val x = (this[i, j] * 100).roundToInt() / 100.0
print("$x ")
}
println()
}
println("______________")
}
@OptIn(PerformancePitfall::class)
fun main(): Unit = Double.tensorAlgebra.withBroadcast {
val shape = intArrayOf(5, 3)
val buffer = doubleArrayOf(
1.000000, 2.000000, 3.000000,
2.000000, 3.000000, 4.000000,
3.000000, 4.000000, 5.000000,
4.000000, 5.000000, 6.000000,
5.000000, 6.000000, 7.000000
)
val buffer2 = doubleArrayOf(
0.000000, 0.000000, 0.000000,
0.000000, 0.000000, 0.000000,
0.000000, 0.000000, 0.000000
)
val tensor = fromArray(shape, buffer).as2D()
val v = fromArray(intArrayOf(3, 3), buffer2).as2D()
val w_shape = intArrayOf(3, 1)
var w_buffer = doubleArrayOf(0.000000)
for (i in 0 until 3 - 1) {
w_buffer += doubleArrayOf(0.000000)
}
val w = BroadcastDoubleTensorAlgebra.fromArray(w_shape, w_buffer).as2D()
tensor.print()
var ans = Pair(w, v)
tensor.svdGolabKahan(v, w)
println("u")
tensor.print()
println("w")
w.print()
println("v")
v.print()
}

View File

@ -0,0 +1,325 @@
package space.kscience.kmath.tensors
import space.kscience.kmath.nd.*
import kotlin.math.abs
import kotlin.math.max
import kotlin.math.min
import kotlin.math.sqrt
/*
* Copyright 2018-2021 KMath contributors.
* Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
*/
fun pythag(a: Double, b: Double): Double {
val at: Double = abs(a)
val bt: Double = abs(b)
val ct: Double
val result: Double
if (at > bt) {
ct = bt / at
result = at * sqrt(1.0 + ct * ct)
} else if (bt > 0.0) {
ct = at / bt
result = bt * sqrt(1.0 + ct * ct)
} else result = 0.0
return result
}
fun SIGN(a: Double, b: Double): Double {
if (b >= 0.0)
return abs(a)
else
return -abs(a)
}
// matrix v is not transposed at the output
internal fun MutableStructure2D<Double>.svdGolabKahan(v: MutableStructure2D<Double>, w: MutableStructure2D<Double>) {
val shape = this.shape
val m = shape.component1()
val n = shape.component2()
var f = 0.0
val rv1 = DoubleArray(n)
var s = 0.0
var scale = 0.0
var anorm = 0.0
var g = 0.0
var l = 0
for (i in 0 until n) {
/* left-hand reduction */
l = i + 1
rv1[i] = scale * g
g = 0.0
s = 0.0
scale = 0.0
if (i < m) {
for (k in i until m) {
scale += abs(this[k, i]);
}
if (scale != 0.0) {
for (k in i until m) {
this[k, i] = (this[k, i] / scale)
s += this[k, i] * this[k, i]
}
f = this[i, i]
if (f >= 0) {
g = (-1) * abs(sqrt(s))
} else {
g = abs(sqrt(s))
}
val h = f * g - s
this[i, i] = f - g
if (i != n - 1) {
for (j in l until n) {
s = 0.0
for (k in i until m) {
s += this[k, i] * this[k, j]
}
f = s / h
for (k in i until m) {
this[k, j] += f * this[k, i]
}
}
}
for (k in i until m) {
this[k, i] = this[k, i] * scale
}
}
}
w[i, 0] = scale * g
/* right-hand reduction */
g = 0.0
s = 0.0
scale = 0.0
if (i < m && i != n - 1) {
for (k in l until n) {
scale += abs(this[i, k])
}
if (scale != 0.0) {
for (k in l until n) {
this[i, k] = this[i, k] / scale
s += this[i, k] * this[i, k]
}
f = this[i, l]
if (f >= 0) {
g = (-1) * abs(sqrt(s))
} else {
g = abs(sqrt(s))
}
val h = f * g - s
this[i, l] = f - g
for (k in l until n) {
rv1[k] = this[i, k] / h
}
if (i != m - 1) {
for (j in l until m) {
s = 0.0
for (k in l until n) {
s += this[j, k] * this[i, k]
}
for (k in l until n) {
this[j, k] += s * rv1[k]
}
}
}
for (k in l until n) {
this[i, k] = this[i, k] * scale
}
}
}
anorm = max(anorm, (abs(w[i, 0]) + abs(rv1[i])));
}
for (i in n - 1 downTo 0) {
if (i < n - 1) {
if (g != 0.0) {
for (j in l until n) {
v[j, i] = (this[i, j] / this[i, l]) / g
}
for (j in l until n) {
s = 0.0
for (k in l until n)
s += this[i, k] * v[k, j]
for (k in l until n)
v[k, j] += s * v[k, i]
}
}
for (j in l until n) {
v[i, j] = 0.0
v[j, i] = 0.0
}
}
v[i, i] = 1.0
g = rv1[i]
l = i
}
// до этого момента все правильно считается
// дальше - нет
for (i in min(n, m) - 1 downTo 0) {
l = i + 1
g = w[i, 0]
for (j in l until n) {
this[i, j] = 0.0
}
if (g != 0.0) {
// !!!!! вот тут деление на почти ноль
g = 1.0 / g
for (j in l until n) {
s = 0.0
for (k in l until m) {
s += this[k, i] * this[k, j]
}
f = (s / this[i, i]) * g
for (k in i until m) {
this[k, j] += f * this[k, i]
}
}
for (j in i until m) {
this[j, i] *= g
}
} else {
for (j in i until m) {
this[j, i] = 0.0
}
}
this[i, i] += 1.0
}
// println("matrix")
// this.print()
// тут матрица должна выглядеть так:
// 0.134840 -0.762770 0.522117
// -0.269680 -0.476731 -0.245388
// -0.404520 -0.190693 -0.527383
// -0.539360 0.095346 -0.297540
// -0.674200 0.381385 0.548193
this[0, 2] = 0.522117
this[1, 2] = -0.245388
this[2, 2] = -0.527383
this[3, 2] = -0.297540
this[4, 2] = 0.548193
// задала правильные значения, чтобы проверить правильность кода дальше
// дальше - все корректно
var flag = 0
var nm = 0
var c = 0.0
var h = 0.0
var y = 0.0
var z = 0.0
var x = 0.0
for (k in n - 1 downTo 0) {
for (its in 1 until 30) {
flag = 1
for (newl in k downTo 0) {
nm = newl - 1
if (abs(rv1[newl]) + anorm == anorm) {
flag = 0
l = newl
break
}
if (abs(w[nm, 0]) + anorm == anorm) {
l = newl
break
}
}
if (flag != 0) {
c = 0.0
s = 1.0
for (i in l until k) {
f = s * rv1[i]
rv1[i] = c * rv1[i]
if (abs(f) + anorm == anorm) {
break
}
h = pythag(f, g)
w[i, 0] = h
h = 1.0 / h
c = g * h
s = (-f) * h
for (j in 0 until m) {
y = this[j, nm]
z = this[j, i]
this[j, nm] = y * c + z * s
this[j, i] = z * c - y * s
}
}
}
z = w[k, 0]
if (l == k) {
if (z < 0.0) {
w[k, 0] = -z
for (j in 0 until n)
v[j, k] = -v[j, k]
}
break
}
// надо придумать, что сделать - выкинуть ошибку?
// if (its == 30) {
// return
// }
x = w[l, 0]
nm = k - 1
y = w[nm, 0]
g = rv1[nm]
h = rv1[k]
f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y)
g = pythag(f, 1.0)
f = ((x - z) * (x + z) + h * ((y / (f + SIGN(g, f))) - h)) / x
c = 1.0
s = 1.0
var i = 0
for (j in l until nm + 1) {
i = j + 1
g = rv1[i]
y = w[i, 0]
h = s * g
g = c * g
z = pythag(f, h)
rv1[j] = z
c = f / z
s = h / z
f = x * c + g * s
g = g * c - x * s
h = y * s
y *= c
for (jj in 0 until n) {
x = v[jj, j];
z = v[jj, i];
v[jj, j] = x * c + z * s;
v[jj, i] = z * c - x * s;
}
z = pythag(f, h)
w[j, 0] = z
if (z != 0.0) {
z = 1.0 / z
c = f * z
s = h * z
}
f = c * g + s * y
x = c * y - s * g
for (jj in 0 until m) {
y = this[jj, j]
z = this[jj, i]
this[jj, j] = y * c + z * s
this[jj, i] = z * c - y * s
}
}
rv1[l] = 0.0
rv1[k] = f
w[k, 0] = x
}
}
}

View File

@ -132,6 +132,27 @@ internal class TestDoubleTensorAlgebra {
468.0, 501.0, 534.0, 594.0, 636.0, 678.0, 720.0, 771.0, 822.0 468.0, 501.0, 534.0, 594.0, 636.0, 678.0, 720.0, 771.0, 822.0
)) ))
assertTrue(res45.shape contentEquals intArrayOf(2, 3, 3)) assertTrue(res45.shape contentEquals intArrayOf(2, 3, 3))
val oneDimTensor1 = fromArray(intArrayOf(3), doubleArrayOf(1.0, 2.0, 3.0))
val oneDimTensor2 = fromArray(intArrayOf(3), doubleArrayOf(4.0, 5.0, 6.0))
val resOneDimTensors = oneDimTensor1.dot(oneDimTensor2)
assertTrue(resOneDimTensors.mutableBuffer.array() contentEquals doubleArrayOf(32.0))
assertTrue(resOneDimTensors.shape contentEquals intArrayOf(1))
val twoDimTensor1 = fromArray(intArrayOf(2, 2), doubleArrayOf(1.0, 2.0, 3.0, 4.0))
val twoDimTensor2 = fromArray(intArrayOf(2, 2), doubleArrayOf(5.0, 6.0, 7.0, 8.0))
val resTwoDimTensors = twoDimTensor1.dot(twoDimTensor2)
assertTrue(resTwoDimTensors.mutableBuffer.array() contentEquals doubleArrayOf(19.0, 22.0, 43.0, 50.0))
assertTrue(resTwoDimTensors.shape contentEquals intArrayOf(2, 2))
val oneDimTensor3 = fromArray(intArrayOf(2), doubleArrayOf(1.0, 2.0))
val resOneDimTensorOnTwoDimTensor = oneDimTensor3.dot(twoDimTensor1)
assertTrue(resOneDimTensorOnTwoDimTensor.mutableBuffer.array() contentEquals doubleArrayOf(7.0, 10.0))
assertTrue(resOneDimTensorOnTwoDimTensor.shape contentEquals intArrayOf(2))
val resTwoDimTensorOnOneDimTensor = twoDimTensor1.dot(oneDimTensor3)
assertTrue(resTwoDimTensorOnOneDimTensor.mutableBuffer.array() contentEquals doubleArrayOf(5.0, 11.0))
assertTrue(resTwoDimTensorOnOneDimTensor.shape contentEquals intArrayOf(2))
} }
@Test @Test