lu inv and det complete + tests
This commit is contained in:
parent
5c0674f1f5
commit
206bcfc909
@ -13,7 +13,7 @@ public interface LinearOpsTensorAlgebra<T, TensorType : TensorStructure<T>, Inde
|
|||||||
//https://pytorch.org/docs/stable/linalg.html#torch.linalg.qr
|
//https://pytorch.org/docs/stable/linalg.html#torch.linalg.qr
|
||||||
public fun TensorType.qr(): TensorType
|
public fun TensorType.qr(): TensorType
|
||||||
|
|
||||||
//https://pytorch.org/docs/stable/generated/torch.lu.html
|
//htt ps://pytorch.org/docs/stable/generated/torch.lu.html
|
||||||
public fun TensorType.lu(): Pair<TensorType, IndexTensorType>
|
public fun TensorType.lu(): Pair<TensorType, IndexTensorType>
|
||||||
|
|
||||||
//https://pytorch.org/docs/stable/generated/torch.lu_unpack.html
|
//https://pytorch.org/docs/stable/generated/torch.lu_unpack.html
|
||||||
|
@ -1,5 +1,8 @@
|
|||||||
package space.kscience.kmath.tensors.core
|
package space.kscience.kmath.tensors.core
|
||||||
|
|
||||||
|
import space.kscience.kmath.nd.MutableStructure2D
|
||||||
|
import space.kscience.kmath.nd.Structure1D
|
||||||
|
import space.kscience.kmath.nd.Structure2D
|
||||||
import space.kscience.kmath.tensors.LinearOpsTensorAlgebra
|
import space.kscience.kmath.tensors.LinearOpsTensorAlgebra
|
||||||
import kotlin.math.sqrt
|
import kotlin.math.sqrt
|
||||||
|
|
||||||
@ -20,6 +23,8 @@ public class DoubleLinearOpsTensorAlgebra :
|
|||||||
val n = shape.size
|
val n = shape.size
|
||||||
val m = shape.last()
|
val m = shape.last()
|
||||||
val pivotsShape = IntArray(n - 1) { i -> shape[i] }
|
val pivotsShape = IntArray(n - 1) { i -> shape[i] }
|
||||||
|
pivotsShape[n - 2] = m + 1
|
||||||
|
|
||||||
val pivotsTensor = IntTensor(
|
val pivotsTensor = IntTensor(
|
||||||
pivotsShape,
|
pivotsShape,
|
||||||
IntArray(pivotsShape.reduce(Int::times)) { 0 }
|
IntArray(pivotsShape.reduce(Int::times)) { 0 }
|
||||||
@ -54,6 +59,8 @@ public class DoubleLinearOpsTensorAlgebra :
|
|||||||
lu[maxInd, k] = tmp
|
lu[maxInd, k] = tmp
|
||||||
}
|
}
|
||||||
|
|
||||||
|
pivots[m] += 1
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
for (j in i + 1 until m) {
|
for (j in i + 1 until m) {
|
||||||
@ -146,6 +153,78 @@ public class DoubleLinearOpsTensorAlgebra :
|
|||||||
TODO("ANDREI")
|
TODO("ANDREI")
|
||||||
}
|
}
|
||||||
|
|
||||||
|
private fun luMatrixDet(lu: Structure2D<Double>, pivots: Structure1D<Int>): Double {
|
||||||
|
val m = lu.shape[0]
|
||||||
|
val sign = if((pivots[m] - m) % 2 == 0) 1.0 else -1.0
|
||||||
|
var det = sign
|
||||||
|
for (i in 0 until m){
|
||||||
|
det *= lu[i, i]
|
||||||
|
}
|
||||||
|
return det
|
||||||
|
}
|
||||||
|
|
||||||
|
public fun DoubleTensor.detLU(): DoubleTensor {
|
||||||
|
val (luTensor, pivotsTensor) = this.lu()
|
||||||
|
val n = shape.size
|
||||||
|
|
||||||
|
val detTensorShape = IntArray(n - 1) { i -> shape[i] }
|
||||||
|
detTensorShape[n - 2] = 1
|
||||||
|
val resBuffer = DoubleArray(detTensorShape.reduce(Int::times)) { 0.0 }
|
||||||
|
|
||||||
|
val detTensor = DoubleTensor(
|
||||||
|
detTensorShape,
|
||||||
|
resBuffer
|
||||||
|
)
|
||||||
|
|
||||||
|
luTensor.matrixSequence().zip(pivotsTensor.vectorSequence()).forEachIndexed { index, (luMatrix, pivots) ->
|
||||||
|
resBuffer[index] = luMatrixDet(luMatrix, pivots)
|
||||||
|
}
|
||||||
|
|
||||||
|
return detTensor
|
||||||
|
}
|
||||||
|
|
||||||
|
private fun luMatrixInv(
|
||||||
|
lu: Structure2D<Double>,
|
||||||
|
pivots: Structure1D<Int>,
|
||||||
|
invMatrix : MutableStructure2D<Double>
|
||||||
|
): Unit {
|
||||||
|
val m = lu.shape[0]
|
||||||
|
|
||||||
|
for (j in 0 until m) {
|
||||||
|
for (i in 0 until m) {
|
||||||
|
if (pivots[i] == j){
|
||||||
|
invMatrix[i, j] = 1.0
|
||||||
|
}
|
||||||
|
|
||||||
|
for (k in 0 until i){
|
||||||
|
invMatrix[i, j] -= lu[i, k] * invMatrix[k, j]
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
for (i in m - 1 downTo 0) {
|
||||||
|
for (k in i + 1 until m) {
|
||||||
|
invMatrix[i, j] -= lu[i, k] * invMatrix[k, j]
|
||||||
|
}
|
||||||
|
invMatrix[i, j] /= lu[i, i]
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
public fun DoubleTensor.invLU(): DoubleTensor {
|
||||||
|
val (luTensor, pivotsTensor) = this.lu()
|
||||||
|
val n = shape.size
|
||||||
|
val invTensor = luTensor.zeroesLike()
|
||||||
|
|
||||||
|
for (
|
||||||
|
(luP, invMatrix) in
|
||||||
|
luTensor.matrixSequence().zip(pivotsTensor.vectorSequence()).zip(invTensor.matrixSequence())
|
||||||
|
) {
|
||||||
|
val (lu, pivots) = luP
|
||||||
|
luMatrixInv(lu, pivots, invMatrix)
|
||||||
|
}
|
||||||
|
|
||||||
|
return invTensor
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
public inline fun <R> DoubleLinearOpsTensorAlgebra(block: DoubleLinearOpsTensorAlgebra.() -> R): R =
|
public inline fun <R> DoubleLinearOpsTensorAlgebra(block: DoubleLinearOpsTensorAlgebra.() -> R): R =
|
||||||
|
@ -1,5 +1,8 @@
|
|||||||
package space.kscience.kmath.tensors.core
|
package space.kscience.kmath.tensors.core
|
||||||
|
|
||||||
|
import space.kscience.kmath.linear.Matrix
|
||||||
|
import space.kscience.kmath.nd.MutableStructure2D
|
||||||
|
import space.kscience.kmath.nd.Structure2D
|
||||||
import space.kscience.kmath.tensors.TensorPartialDivisionAlgebra
|
import space.kscience.kmath.tensors.TensorPartialDivisionAlgebra
|
||||||
import kotlin.math.abs
|
import kotlin.math.abs
|
||||||
|
|
||||||
@ -262,7 +265,31 @@ public open class DoubleTensorAlgebra : TensorPartialDivisionAlgebra<Double, Dou
|
|||||||
}
|
}
|
||||||
|
|
||||||
override fun DoubleTensor.det(): DoubleTensor {
|
override fun DoubleTensor.det(): DoubleTensor {
|
||||||
TODO("ANDREI")
|
TODO()
|
||||||
|
/*
|
||||||
|
checkSquareMatrix(shape)
|
||||||
|
|
||||||
|
val n = shape.size
|
||||||
|
val m = shape.last()
|
||||||
|
|
||||||
|
val detTensorShape = IntArray(n - 1) { i -> shape[i] }
|
||||||
|
detTensorShape[n - 1] = 1
|
||||||
|
val resBuffer = DoubleArray(detTensorShape.reduce(Int::times)) { 0.0 }
|
||||||
|
|
||||||
|
val detTensor = DoubleTensor(
|
||||||
|
detTensorShape,
|
||||||
|
resBuffer
|
||||||
|
)
|
||||||
|
|
||||||
|
this.matrixSequence().forEachIndexed{i, matrix ->
|
||||||
|
// todo need Matrix determinant algo
|
||||||
|
// todo resBuffer[i] = matrix.det()
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
return detTensor
|
||||||
|
|
||||||
|
*/
|
||||||
}
|
}
|
||||||
|
|
||||||
override fun DoubleTensor.square(): DoubleTensor {
|
override fun DoubleTensor.square(): DoubleTensor {
|
||||||
@ -294,7 +321,7 @@ public open class DoubleTensorAlgebra : TensorPartialDivisionAlgebra<Double, Dou
|
|||||||
}
|
}
|
||||||
|
|
||||||
public fun DoubleTensor.contentEquals(other: DoubleTensor, eqFunction: (Double, Double) -> Boolean): Boolean {
|
public fun DoubleTensor.contentEquals(other: DoubleTensor, eqFunction: (Double, Double) -> Boolean): Boolean {
|
||||||
if (!(this.shape contentEquals other.shape)){
|
if (!(this.shape contentEquals other.shape)) {
|
||||||
return false
|
return false
|
||||||
}
|
}
|
||||||
return this.eq(other, eqFunction)
|
return this.eq(other, eqFunction)
|
||||||
@ -303,10 +330,10 @@ public open class DoubleTensorAlgebra : TensorPartialDivisionAlgebra<Double, Dou
|
|||||||
public fun DoubleTensor.eq(other: DoubleTensor, eqFunction: (Double, Double) -> Boolean): Boolean {
|
public fun DoubleTensor.eq(other: DoubleTensor, eqFunction: (Double, Double) -> Boolean): Boolean {
|
||||||
// todo broadcasting checking
|
// todo broadcasting checking
|
||||||
val n = this.linearStructure.size
|
val n = this.linearStructure.size
|
||||||
if (n != other.linearStructure.size){
|
if (n != other.linearStructure.size) {
|
||||||
return false
|
return false
|
||||||
}
|
}
|
||||||
for (i in 0 until n){
|
for (i in 0 until n) {
|
||||||
if (!eqFunction(this.buffer[this.bufferStart + i], other.buffer[other.bufferStart + i])) {
|
if (!eqFunction(this.buffer[this.bufferStart + i], other.buffer[other.bufferStart + i])) {
|
||||||
return false
|
return false
|
||||||
}
|
}
|
||||||
|
@ -3,7 +3,6 @@ package space.kscience.kmath.tensors.core
|
|||||||
import kotlin.math.abs
|
import kotlin.math.abs
|
||||||
import kotlin.math.exp
|
import kotlin.math.exp
|
||||||
import kotlin.test.Test
|
import kotlin.test.Test
|
||||||
import kotlin.test.assertEquals
|
|
||||||
import kotlin.test.assertTrue
|
import kotlin.test.assertTrue
|
||||||
|
|
||||||
class TestDoubleAnalyticTensorAlgebra {
|
class TestDoubleAnalyticTensorAlgebra {
|
||||||
@ -16,9 +15,9 @@ class TestDoubleAnalyticTensorAlgebra {
|
|||||||
return this.map(transform).toDoubleArray()
|
return this.map(transform).toDoubleArray()
|
||||||
}
|
}
|
||||||
|
|
||||||
fun DoubleArray.deltaEqual(other: DoubleArray, delta: Double = 1e-5): Boolean {
|
fun DoubleArray.epsEqual(other: DoubleArray, eps: Double = 1e-5): Boolean {
|
||||||
for ((elem1, elem2) in this.asSequence().zip(other.asSequence())) {
|
for ((elem1, elem2) in this.asSequence().zip(other.asSequence())) {
|
||||||
if (abs(elem1 - elem2) > delta) {
|
if (abs(elem1 - elem2) > eps) {
|
||||||
return false
|
return false
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
@ -29,7 +28,7 @@ class TestDoubleAnalyticTensorAlgebra {
|
|||||||
fun testExp() = DoubleAnalyticTensorAlgebra {
|
fun testExp() = DoubleAnalyticTensorAlgebra {
|
||||||
tensor.exp().let {
|
tensor.exp().let {
|
||||||
assertTrue { shape contentEquals it.shape }
|
assertTrue { shape contentEquals it.shape }
|
||||||
assertTrue { buffer.fmap(::exp).deltaEqual(it.buffer.array())}
|
assertTrue { buffer.fmap(::exp).epsEqual(it.buffer.array())}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
@ -0,0 +1,73 @@
|
|||||||
|
package space.kscience.kmath.tensors.core
|
||||||
|
|
||||||
|
import kotlin.math.abs
|
||||||
|
import kotlin.math.exp
|
||||||
|
import kotlin.test.Test
|
||||||
|
import kotlin.test.assertEquals
|
||||||
|
import kotlin.test.assertTrue
|
||||||
|
|
||||||
|
class TestDoubleLinearOpsTensorAlgebra {
|
||||||
|
|
||||||
|
private val eps = 1e-5
|
||||||
|
|
||||||
|
private fun Double.epsEqual(other: Double): Boolean {
|
||||||
|
return abs(this - other) < eps
|
||||||
|
}
|
||||||
|
|
||||||
|
fun DoubleArray.epsEqual(other: DoubleArray, eps: Double = 1e-5): Boolean {
|
||||||
|
for ((elem1, elem2) in this.asSequence().zip(other.asSequence())) {
|
||||||
|
if (abs(elem1 - elem2) > eps) {
|
||||||
|
return false
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return true
|
||||||
|
}
|
||||||
|
|
||||||
|
@Test
|
||||||
|
fun testDetLU() = DoubleLinearOpsTensorAlgebra {
|
||||||
|
val tensor = fromArray(
|
||||||
|
intArrayOf(2, 2, 2),
|
||||||
|
doubleArrayOf(
|
||||||
|
1.0, 3.0,
|
||||||
|
1.0, 2.0,
|
||||||
|
1.5, 1.0,
|
||||||
|
10.0, 2.0
|
||||||
|
)
|
||||||
|
)
|
||||||
|
|
||||||
|
val expectedShape = intArrayOf(2, 1)
|
||||||
|
val expectedBuffer = doubleArrayOf(
|
||||||
|
-1.0,
|
||||||
|
-7.0
|
||||||
|
)
|
||||||
|
val detTensor = tensor.detLU()
|
||||||
|
|
||||||
|
assertTrue { detTensor.shape contentEquals expectedShape }
|
||||||
|
assertTrue { detTensor.buffer.array().epsEqual(expectedBuffer) }
|
||||||
|
}
|
||||||
|
|
||||||
|
@Test
|
||||||
|
fun testInvLU() = DoubleLinearOpsTensorAlgebra {
|
||||||
|
val tensor = fromArray(
|
||||||
|
intArrayOf(2, 2, 2),
|
||||||
|
doubleArrayOf(
|
||||||
|
1.0, 0.0,
|
||||||
|
0.0, 2.0,
|
||||||
|
1.0, 1.0,
|
||||||
|
1.0, 0.0
|
||||||
|
)
|
||||||
|
)
|
||||||
|
|
||||||
|
val expectedShape = intArrayOf(2, 2, 2)
|
||||||
|
val expectedBuffer = doubleArrayOf(
|
||||||
|
1.0, 0.0,
|
||||||
|
0.0, 0.5,
|
||||||
|
0.0, 1.0,
|
||||||
|
1.0, -1.0
|
||||||
|
)
|
||||||
|
|
||||||
|
val invTensor = tensor.invLU()
|
||||||
|
assertTrue { invTensor.shape contentEquals expectedShape }
|
||||||
|
assertTrue { invTensor.buffer.array().epsEqual(expectedBuffer) }
|
||||||
|
}
|
||||||
|
}
|
Loading…
Reference in New Issue
Block a user