KMP library for tensors #300

Merged
grinisrit merged 215 commits from feature/tensor-algebra into dev 2021-05-08 09:48:04 +03:00
8 changed files with 51 additions and 33 deletions
Showing only changes of commit 21b5d45b96 - Show all commits

View File

@ -69,7 +69,7 @@ fun main () {
val n = l.shape[0] val n = l.shape[0]
val x = zeros(intArrayOf(n)) val x = zeros(intArrayOf(n))
for (i in 0 until n){ for (i in 0 until n){
x[intArrayOf(i)] = (b[intArrayOf(i)] - l[i].dot(x).valueOrNull()!!) / l[intArrayOf(i, i)] x[intArrayOf(i)] = (b[intArrayOf(i)] - l[i].dot(x).value()) / l[intArrayOf(i, i)]
} }
return x return x
} }

View File

@ -60,7 +60,7 @@ fun main() {
require(yTrue.shape contentEquals yPred.shape) require(yTrue.shape contentEquals yPred.shape)
val diff = yTrue - yPred val diff = yTrue - yPred
return diff.dot(diff).sqrt().valueOrNull()!! return diff.dot(diff).sqrt().value()
} }
println("MSE: ${mse(alpha, alphaOLS)}") println("MSE: ${mse(alpha, alphaOLS)}")

View File

@ -15,13 +15,21 @@ import space.kscience.kmath.operations.Algebra
*/ */
public interface TensorAlgebra<T>: Algebra<Tensor<T>> { public interface TensorAlgebra<T>: Algebra<Tensor<T>> {
/**
*
* Returns a single tensor value of unit dimension if tensor shape equals to [1].
*
* @return a nullable value of a potentially scalar tensor.
*/
public fun Tensor<T>.valueOrNull(): T?
/** /**
* *
* Returns a single tensor value of unit dimension. The tensor shape must be equal to [1]. * Returns a single tensor value of unit dimension. The tensor shape must be equal to [1].
* *
* @return the value of a scalar tensor. * @return the value of a scalar tensor.
*/ */
public fun Tensor<T>.valueOrNull(): T? public fun Tensor<T>.value(): T
/** /**
* Each element of the tensor [other] is added to this value. * Each element of the tensor [other] is added to this value.

View File

@ -35,10 +35,11 @@ public open class DoubleTensorAlgebra :
public companion object : DoubleTensorAlgebra() public companion object : DoubleTensorAlgebra()
override fun Tensor<Double>.valueOrNull(): Double? = if(tensor.shape contentEquals intArrayOf(1)) { override fun Tensor<Double>.valueOrNull(): Double? = if (tensor.shape contentEquals intArrayOf(1))
// Inconsistent value for tensor of with this shape tensor.mutableBuffer.array()[tensor.bufferStart] else null
tensor.mutableBuffer.array()[tensor.bufferStart]
} else null override fun Tensor<Double>.value(): Double =
valueOrNull() ?: throw IllegalArgumentException("Inconsistent value for tensor of with $shape shape")
/** /**
* Constructs a tensor with the specified shape and data. * Constructs a tensor with the specified shape and data.
@ -62,8 +63,10 @@ public open class DoubleTensorAlgebra :
* @return tensor with the [shape] shape and data generated by the [initializer]. * @return tensor with the [shape] shape and data generated by the [initializer].
*/ */
public fun produce(shape: IntArray, initializer: (IntArray) -> Double): DoubleTensor = public fun produce(shape: IntArray, initializer: (IntArray) -> Double): DoubleTensor =
fromArray(shape, fromArray(
TensorLinearStructure(shape).indices().map(initializer).toMutableList().toDoubleArray()) shape,
TensorLinearStructure(shape).indices().map(initializer).toMutableList().toDoubleArray()
)
override operator fun Tensor<Double>.get(i: Int): DoubleTensor { override operator fun Tensor<Double>.get(i: Int): DoubleTensor {
val lastShape = tensor.shape.drop(1).toIntArray() val lastShape = tensor.shape.drop(1).toIntArray()
@ -621,7 +624,7 @@ public open class DoubleTensorAlgebra :
keepDim keepDim
) )
private fun cov(x: DoubleTensor, y:DoubleTensor): Double{ private fun cov(x: DoubleTensor, y: DoubleTensor): Double {
val n = x.shape[0] val n = x.shape[0]
return ((x - x.mean()) * (y - y.mean())).mean() * n / (n - 1) return ((x - x.mean()) * (y - y.mean())).mean() * n / (n - 1)
} }
@ -633,10 +636,10 @@ public open class DoubleTensorAlgebra :
check(tensors.all { it.shape contentEquals intArrayOf(m) }) { "Tensors must have same shapes" } check(tensors.all { it.shape contentEquals intArrayOf(m) }) { "Tensors must have same shapes" }
val resTensor = DoubleTensor( val resTensor = DoubleTensor(
intArrayOf(n, n), intArrayOf(n, n),
DoubleArray(n * n) {0.0} DoubleArray(n * n) { 0.0 }
) )
for (i in 0 until n){ for (i in 0 until n) {
for (j in 0 until n){ for (j in 0 until n) {
resTensor[intArrayOf(i, j)] = cov(tensors[i].tensor, tensors[j].tensor) resTensor[intArrayOf(i, j)] = cov(tensors[i].tensor, tensors[j].tensor)
} }
} }
@ -779,8 +782,10 @@ public open class DoubleTensorAlgebra :
val qTensor = zeroesLike() val qTensor = zeroesLike()
val rTensor = zeroesLike() val rTensor = zeroesLike()
tensor.matrixSequence() tensor.matrixSequence()
.zip((qTensor.matrixSequence() .zip(
.zip(rTensor.matrixSequence()))).forEach { (matrix, qr) -> (qTensor.matrixSequence()
.zip(rTensor.matrixSequence()))
).forEach { (matrix, qr) ->
val (q, r) = qr val (q, r) = qr
qrHelper(matrix.asTensor(), q.asTensor(), r.as2D()) qrHelper(matrix.asTensor(), q.asTensor(), r.as2D())
} }
@ -812,9 +817,13 @@ public open class DoubleTensorAlgebra :
val vTensor = zeros(commonShape + intArrayOf(min(n, m), m)) val vTensor = zeros(commonShape + intArrayOf(min(n, m), m))
tensor.matrixSequence() tensor.matrixSequence()
.zip(uTensor.matrixSequence() .zip(
.zip(sTensor.vectorSequence() uTensor.matrixSequence()
.zip(vTensor.matrixSequence()))).forEach { (matrix, USV) -> .zip(
sTensor.vectorSequence()
.zip(vTensor.matrixSequence())
)
).forEach { (matrix, USV) ->
val matrixSize = matrix.shape.reduce { acc, i -> acc * i } val matrixSize = matrix.shape.reduce { acc, i -> acc * i }
val curMatrix = DoubleTensor( val curMatrix = DoubleTensor(
matrix.shape, matrix.shape,
@ -918,10 +927,11 @@ public open class DoubleTensorAlgebra :
* @return triple of P, L and U tensors * @return triple of P, L and U tensors
*/ */
public fun Tensor<Double>.lu(epsilon: Double = 1e-9): Triple<DoubleTensor, DoubleTensor, DoubleTensor> { public fun Tensor<Double>.lu(epsilon: Double = 1e-9): Triple<DoubleTensor, DoubleTensor, DoubleTensor> {
val (lu, pivots) = this.luFactor(epsilon) val (lu, pivots) = tensor.luFactor(epsilon)
return luPivot(lu, pivots) return luPivot(lu, pivots)
} }
override fun Tensor<Double>.lu(): Triple<DoubleTensor, DoubleTensor, DoubleTensor> = lu(1e-9) override fun Tensor<Double>.lu(): Triple<DoubleTensor, DoubleTensor, DoubleTensor> = lu(1e-9)
} }

View File

@ -58,7 +58,7 @@ internal fun DoubleTensorAlgebra.checkSymmetric(
internal fun DoubleTensorAlgebra.checkPositiveDefinite(tensor: DoubleTensor, epsilon: Double = 1e-6) { internal fun DoubleTensorAlgebra.checkPositiveDefinite(tensor: DoubleTensor, epsilon: Double = 1e-6) {
checkSymmetric(tensor, epsilon) checkSymmetric(tensor, epsilon)
for (mat in tensor.matrixSequence()) for (mat in tensor.matrixSequence())
check(mat.asTensor().detLU().valueOrNull()!! > 0.0) { check(mat.asTensor().detLU().value() > 0.0) {
"Tensor contains matrices which are not positive definite ${mat.asTensor().detLU().valueOrNull()!!}" "Tensor contains matrices which are not positive definite ${mat.asTensor().detLU().value()}"
} }
} }

View File

@ -240,14 +240,14 @@ internal fun DoubleTensorAlgebra.qrHelper(
val vv = v.as1D() val vv = v.as1D()
if (j > 0) { if (j > 0) {
for (i in 0 until j) { for (i in 0 until j) {
r[i, j] = (qT[i] dot matrixT[j]).valueOrNull()!! r[i, j] = (qT[i] dot matrixT[j]).value()
for (k in 0 until n) { for (k in 0 until n) {
val qTi = qT[i].as1D() val qTi = qT[i].as1D()
vv[k] = vv[k] - r[i, j] * qTi[k] vv[k] = vv[k] - r[i, j] * qTi[k]
} }
} }
} }
r[j, j] = DoubleTensorAlgebra { (v dot v).sqrt().valueOrNull()!! } r[j, j] = DoubleTensorAlgebra { (v dot v).sqrt().value() }
for (i in 0 until n) { for (i in 0 until n) {
qM[i, j] = vv[i] / r[j, j] qM[i, j] = vv[i] / r[j, j]
} }
@ -270,9 +270,9 @@ internal fun DoubleTensorAlgebra.svd1d(a: DoubleTensor, epsilon: Double = 1e-10)
while (true) { while (true) {
lastV = v lastV = v
v = b.dot(lastV) v = b.dot(lastV)
val norm = DoubleTensorAlgebra { (v dot v).sqrt().valueOrNull()!! } val norm = DoubleTensorAlgebra { (v dot v).sqrt().value() }
v = v.times(1.0 / norm) v = v.times(1.0 / norm)
if (abs(v.dot(lastV).valueOrNull()!!) > 1 - epsilon) { if (abs(v.dot(lastV).value()) > 1 - epsilon) {
return v return v
} }
} }
@ -293,7 +293,7 @@ internal fun DoubleTensorAlgebra.svdHelper(
val outerProduct = DoubleArray(u.shape[0] * v.shape[0]) val outerProduct = DoubleArray(u.shape[0] * v.shape[0])
for (i in 0 until u.shape[0]) { for (i in 0 until u.shape[0]) {
for (j in 0 until v.shape[0]) { for (j in 0 until v.shape[0]) {
outerProduct[i * v.shape[0] + j] = u[i].valueOrNull()!! * v[j].valueOrNull()!! 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)) a = a - singularValue.times(DoubleTensor(intArrayOf(u.shape[0], v.shape[0]), outerProduct))
@ -304,12 +304,12 @@ internal fun DoubleTensorAlgebra.svdHelper(
if (n > m) { if (n > m) {
v = svd1d(a, epsilon) v = svd1d(a, epsilon)
u = matrix.dot(v) u = matrix.dot(v)
norm = DoubleTensorAlgebra { (u dot u).sqrt().valueOrNull()!! } norm = DoubleTensorAlgebra { (u dot u).sqrt().value() }
u = u.times(1.0 / norm) u = u.times(1.0 / norm)
} else { } else {
u = svd1d(a, epsilon) u = svd1d(a, epsilon)
v = matrix.transpose(0, 1).dot(u) v = matrix.transpose(0, 1).dot(u)
norm = DoubleTensorAlgebra { (v dot v).sqrt().valueOrNull()!! } norm = DoubleTensorAlgebra { (v dot v).sqrt().value() }
v = v.times(1.0 / norm) v = v.times(1.0 / norm)
} }

View File

@ -46,7 +46,7 @@ internal class TestDoubleLinearOpsTensorAlgebra {
) )
) )
assertTrue { abs(m.det().valueOrNull()!! - expectedValue) < 1e-5 } assertTrue { abs(m.det().value() - expectedValue) < 1e-5 }
} }
@Test @Test
@ -58,7 +58,7 @@ internal class TestDoubleLinearOpsTensorAlgebra {
) )
) )
assertTrue { abs(m.det().valueOrNull()!! - expectedValue) < 1e-5 } assertTrue { abs(m.det().value() - expectedValue) < 1e-5 }
} }
@Test @Test
@ -90,7 +90,7 @@ internal class TestDoubleLinearOpsTensorAlgebra {
fun testScalarProduct() = DoubleTensorAlgebra { fun testScalarProduct() = DoubleTensorAlgebra {
val a = fromArray(intArrayOf(3), doubleArrayOf(1.8, 2.5, 6.8)) val a = fromArray(intArrayOf(3), doubleArrayOf(1.8, 2.5, 6.8))
val b = fromArray(intArrayOf(3), doubleArrayOf(5.5, 2.6, 6.4)) val b = fromArray(intArrayOf(3), doubleArrayOf(5.5, 2.6, 6.4))
assertEquals(a.dot(b).valueOrNull()!!, 59.92) assertEquals(a.dot(b).value(), 59.92)
} }
@Test @Test

View File

@ -21,7 +21,7 @@ internal class TestDoubleTensor {
fun testValue() = DoubleTensorAlgebra { fun testValue() = DoubleTensorAlgebra {
val value = 12.5 val value = 12.5
val tensor = fromArray(intArrayOf(1), doubleArrayOf(value)) val tensor = fromArray(intArrayOf(1), doubleArrayOf(value))
assertEquals(tensor.valueOrNull()!!, value) assertEquals(tensor.value(), value)
} }
@Test @Test