From bbc012d8cd120b33b8a0b7ba54ae7f92a2ca67cd Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Sat, 20 Apr 2019 11:43:30 +0300 Subject: [PATCH] Optimizing decomposition performance --- .../kmath/linear/LinearAlgebraBenchmark.kt | 2 +- build.gradle.kts | 5 +- .../scientifik/kmath/linear/BufferMatrix.kt | 13 + .../kmath/linear/LUPDecomposition.kt | 335 ++++++++---------- .../scientifik/kmath/linear/MatrixContext.kt | 3 +- .../kmath/structures/BufferAccessor2D.kt | 45 +++ .../scientifik/kmath/structures/Buffers.kt | 2 +- 7 files changed, 205 insertions(+), 200 deletions(-) create mode 100644 kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferAccessor2D.kt diff --git a/benchmarks/src/main/kotlin/scientifik/kmath/linear/LinearAlgebraBenchmark.kt b/benchmarks/src/main/kotlin/scientifik/kmath/linear/LinearAlgebraBenchmark.kt index c2d74379e..c397e9bd0 100644 --- a/benchmarks/src/main/kotlin/scientifik/kmath/linear/LinearAlgebraBenchmark.kt +++ b/benchmarks/src/main/kotlin/scientifik/kmath/linear/LinearAlgebraBenchmark.kt @@ -9,7 +9,7 @@ import kotlin.system.measureTimeMillis @ExperimentalContracts fun main() { - val random = Random(12224) + val random = Random(1224) val dim = 100 //creating invertible matrix val u = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 } diff --git a/build.gradle.kts b/build.gradle.kts index 3f25a622c..bff826b37 100644 --- a/build.gradle.kts +++ b/build.gradle.kts @@ -1,9 +1,12 @@ import com.moowork.gradle.node.NodeExtension import com.moowork.gradle.node.npm.NpmTask import com.moowork.gradle.node.task.NodeTask +import org.jetbrains.kotlin.gradle.dsl.KotlinMultiplatformExtension +import org.jetbrains.kotlin.gradle.tasks.Kotlin2JsCompile +import org.jetbrains.kotlin.gradle.tasks.KotlinCompile buildscript { - val kotlinVersion: String by rootProject.extra("1.3.21") + val kotlinVersion: String by rootProject.extra("1.3.30") val ioVersion: String by rootProject.extra("0.1.5") val coroutinesVersion: String by rootProject.extra("1.1.1") val atomicfuVersion: String by rootProject.extra("0.12.1") diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/BufferMatrix.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/BufferMatrix.kt index 739c170d4..34d9a6883 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/BufferMatrix.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/BufferMatrix.kt @@ -1,5 +1,6 @@ package scientifik.kmath.linear +import scientifik.kmath.operations.RealField import scientifik.kmath.operations.Ring import scientifik.kmath.structures.* @@ -23,6 +24,18 @@ class BufferMatrixContext>( } } +object RealMatrixContext : GenericMatrixContext { + + override val elementContext = RealField + + override inline fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> Double): Matrix { + val buffer = DoubleBuffer(rows * columns) { offset -> initializer(offset / columns, offset % columns) } + return BufferMatrix(rows, columns, buffer) + } + + override inline fun point(size: Int, initializer: (Int) -> Double): Point = DoubleBuffer(size,initializer) +} + class BufferMatrix( override val rowNum: Int, override val colNum: Int, diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUPDecomposition.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUPDecomposition.kt index 4bc463333..9c4382e38 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUPDecomposition.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUPDecomposition.kt @@ -3,22 +3,22 @@ package scientifik.kmath.linear import scientifik.kmath.operations.Field import scientifik.kmath.operations.RealField import scientifik.kmath.operations.Ring -import scientifik.kmath.structures.* -import kotlin.contracts.ExperimentalContracts -import kotlin.contracts.InvocationKind -import kotlin.contracts.contract -import kotlin.reflect.KClass +import scientifik.kmath.structures.BufferAccessor2D +import scientifik.kmath.structures.Matrix +import scientifik.kmath.structures.Structure2D /** * Common implementation of [LUPDecompositionFeature] */ class LUPDecomposition( - private val elementContext: Ring, + val context: GenericMatrixContext>, val lu: Structure2D, val pivot: IntArray, private val even: Boolean ) : LUPDecompositionFeature, DeterminantFeature { + val elementContext get() = context.elementContext + /** * Returns the matrix L of the decomposition. * @@ -66,102 +66,14 @@ class LUPDecomposition( } -internal open class BufferAccessor( - val type: KClass, - val field: Field, - val rowNum: Int, - val colNum: Int -) { - open operator fun MutableBuffer.get(i: Int, j: Int) = get(i + colNum * j) - open operator fun MutableBuffer.set(i: Int, j: Int, value: T) { - set(i + colNum * j, value) - } - - fun create(init: (i: Int, j: Int) -> T) = - MutableBuffer.auto(type, rowNum * colNum) { offset -> init(offset / colNum, offset % colNum) } - - fun create(mat: Structure2D) = create { i, j -> mat[i, j] } - - //TODO optimize wrapper - fun MutableBuffer.collect(): Structure2D = - NDStructure.auto(type, rowNum, colNum) { (i, j) -> get(i, j) }.as2D() - - open fun MutableBuffer.innerProduct(row: Int, col: Int, max: Int): T { - var sum = field.zero - field.run { - for (i in 0 until max) { - sum += get(row, i) * get(i, col) - } - } - return sum - } - - open fun MutableBuffer.divideInPlace(i: Int, j: Int, factor: T) { - field.run { set(i, j, get(i, j) / factor) } - } - - open fun MutableBuffer.subtractInPlace(i: Int, j: Int, lu: MutableBuffer, col: Int) { - field.run { - set(i, j, get(i, j) - get(col, j) * lu[i, col]) - } - } -} - -/** - * Specialized LU operations for Doubles - */ -private class RealBufferAccessor(rowNum: Int, colNum: Int) : - BufferAccessor(Double::class, RealField, rowNum, colNum) { - override fun MutableBuffer.get(i: Int, j: Int) = (this as DoubleBuffer).array[i + colNum * j] - override fun MutableBuffer.set(i: Int, j: Int, value: Double) { - (this as DoubleBuffer).array[i + colNum * j] = value - } - - override fun MutableBuffer.innerProduct(row: Int, col: Int, max: Int): Double { - var sum = 0.0 - for (i in 0 until max) { - sum += get(row, i) * get(i, col) - } - return sum - } - - override fun MutableBuffer.divideInPlace(i: Int, j: Int, factor: Double) { - set(i, j, get(i, j) / factor) - } - - override fun MutableBuffer.subtractInPlace(i: Int, j: Int, lu: MutableBuffer, col: Int) { - set(i, j, get(i, j) - get(col, j) * lu[i, col]) - } -} - -@ExperimentalContracts -private inline fun , F : Field> GenericMatrixContext.withAccessor( - type: KClass, - rowNum: Int, - colNum: Int, - block: BufferAccessor.() -> Unit -) { - contract { - callsInPlace(block, InvocationKind.EXACTLY_ONCE) - } - if (elementContext == RealField) { - @Suppress("UNCHECKED_CAST") - RealBufferAccessor(rowNum, colNum) as BufferAccessor - } else { - BufferAccessor(type, elementContext, rowNum, colNum) - }.run(block) -} - -private fun , F : Field> GenericMatrixContext.abs(value: T) = +fun , F : Field> GenericMatrixContext.abs(value: T) = if (value > elementContext.zero) value else with(elementContext) { -value } /** * Create a lup decomposition of generic matrix */ -@ExperimentalContracts -fun , F : Field> GenericMatrixContext.lup( - type: KClass, +inline fun , F : Field> GenericMatrixContext.lup( matrix: Matrix, checkSingular: (T) -> Boolean ): LUPDecomposition { @@ -169,133 +81,166 @@ fun , F : Field> GenericMatrixContext.lup( error("LU decomposition supports only square matrices") } - val m = matrix.colNum val pivot = IntArray(matrix.rowNum) - withAccessor(type, matrix.rowNum, matrix.colNum) { + //TODO just waits for KEEP-176 + BufferAccessor2D(T::class, matrix.rowNum, matrix.colNum).run { + elementContext.run { - val lu = create(matrix) + val lu = create(matrix) - // Initialize permutation array and parity - for (row in 0 until m) { - pivot[row] = row - } - var even = true - - // Loop over columns - for (col in 0 until m) { - - // upper - for (row in 0 until col) { - val sum = lu.innerProduct(row, col, row) - lu[row, col] = field.run { lu[row, col] - sum } + // Initialize permutation array and parity + for (row in 0 until m) { + pivot[row] = row } + var even = true - // lower - val max = (col until m).maxBy { row -> - val sum = lu.innerProduct(row, col, col) - lu[row, col] = field.run { lu[row, col] - sum } - abs(sum) - } ?: col - - // Singularity check - if (checkSingular(lu[max, col])) { - error("Singular matrix") + // Initialize permutation array and parity + for (row in 0 until m) { + pivot[row] = row } + var singular = false - // Pivot if necessary - if (max != col) { - for (i in 0 until m) { - lu[max, i] = lu[col, i] - lu[col, i] = lu[max, i] + // Loop over columns + for (col in 0 until m) { + + // upper + for (row in 0 until col) { + val luRow = lu.row(row) + var sum = luRow[col] + for (i in 0 until row) { + sum -= luRow[i] * lu[i, col] + } + luRow[col] = sum + } + + // lower + var max = col // permutation row + var largest = -one + for (row in col until m) { + val luRow = lu.row(row) + var sum = luRow[col] + for (i in 0 until col) { + sum -= luRow[i] * lu[i, col] + } + luRow[col] = sum + + // maintain best permutation choice + if (abs(sum) > largest) { + largest = abs(sum) + max = row + } + } + + // Singularity check + if (checkSingular(abs(lu[max, col]))) { + error("The matrix is singular") + } + + // Pivot if necessary + if (max != col) { + val luMax = lu.row(max) + val luCol = lu.row(col) + for (i in 0 until m) { + val tmp = luMax[i] + luMax[i] = luCol[i] + luCol[i] = tmp + } + val temp = pivot[max] + pivot[max] = pivot[col] + pivot[col] = temp + even = !even + } + + // Divide the lower elements by the "winning" diagonal elt. + val luDiag = lu[col, col] + for (row in col + 1 until m) { + lu[row, col] /= luDiag } - val temp = pivot[max] - pivot[max] = pivot[col] - pivot[col] = temp - even = !even } - // Divide the lower elements by the "winning" diagonal elt. - val luDiag = lu[col, col] - for (row in col + 1 until m) { - lu.divideInPlace(row, col, luDiag) - //lu[row, col] = lu[row, col] / luDiag - } + return LUPDecomposition(this@lup, lu.collect(), pivot, even) + } + } +} + +fun GenericMatrixContext.lup(matrix: Matrix) = lup(matrix) { it < 1e-11 } + +inline fun LUPDecomposition.solve(matrix: Matrix): Matrix { + + if (matrix.rowNum != pivot.size) { + error("Matrix dimension mismatch. Expected ${pivot.size}, but got ${matrix.colNum}") + } + + BufferAccessor2D(T::class, matrix.rowNum, matrix.colNum).run { + elementContext.run { + + val lu = create{i,j-> this@solve.lu[i,j]} + + // Apply permutations to b + val bp = create { i, j -> zero } + for (row in 0 until pivot.size) { + val bpRow = bp.row(row) + val pRow = pivot[row] + for (col in 0 until matrix.colNum) { + bpRow[col] = matrix[pRow, col] + } + } + + // Solve LY = b + for (col in 0 until pivot.size) { + val bpCol = bp.row(col) + for (i in col + 1 until pivot.size) { + val bpI = bp.row(i) + val luICol = lu[i, col] + for (j in 0 until matrix.colNum) { + bpI[j] -= bpCol[j] * luICol + } + } + } + + // Solve UX = Y + for (col in pivot.size - 1 downTo 0) { + val bpCol = bp.row(col) + val luDiag = lu[col, col] + for (j in 0 until matrix.colNum) { + bpCol[j] /= luDiag + } + for (i in 0 until col) { + val bpI = bp.row(i) + val luICol = lu[i, col] + for (j in 0 until matrix.colNum) { + bpI[j] -= bpCol[j] * luICol + } + } + } + return context.produce(pivot.size, matrix.colNum) { i, j -> bp[i, j] } } - return LUPDecomposition(elementContext, lu.collect(), pivot, even) } } -@ExperimentalContracts -fun GenericMatrixContext.lup(matrix: Matrix) = lup(Double::class, matrix) { it < 1e-11 } /** * Solve a linear equation **a*x = b** */ -@ExperimentalContracts -fun , F : Field> GenericMatrixContext.solve( - type: KClass, +inline fun , F : Field> GenericMatrixContext.solve( a: Matrix, b: Matrix, - checkSingular: (T) -> Boolean + crossinline checkSingular: (T) -> Boolean ): Matrix { - if (b.rowNum != a.colNum) { - error("Matrix dimension mismatch. Expected ${a.rowNum}, but got ${b.colNum}") - } - // Use existing decomposition if it is provided by matrix - val decomposition = a.getFeature() ?: lup(type, a, checkSingular) - - withAccessor(type, a.rowNum, a.colNum) { - - val lu = create(decomposition.lu) - - // Apply permutations to b - val bp = create { i, j -> - b[decomposition.pivot[i], j] - } - - // Solve LY = b - for (col in 0 until a.rowNum) { - for (i in col + 1 until a.rowNum) { - for (j in 0 until b.colNum) { - bp.subtractInPlace(i, j, lu, col) - //bp[i, j] -= bp[col, j] * lu[i, col] - } - } - } - - // Solve UX = Y - for (col in a.rowNum - 1 downTo 0) { - val luDiag = lu[col, col] - for (j in 0 until b.colNum) { - bp.divideInPlace(col, j, luDiag) - //bp[col, j] /= lu[col, col] - } - for (i in 0 until col) { - for (j in 0 until b.colNum) { - bp.subtractInPlace(i, j, lu, col) - //bp[i, j] -= bp[col, j] * lu[i, col] - } - } - } - - return produce(a.rowNum, a.colNum) { i, j -> bp[i, j] } - } + val decomposition = a.getFeature() ?: lup(a, checkSingular) + return decomposition.solve(b) } -@ExperimentalContracts -fun GenericMatrixContext.solve(a: Matrix, b: Matrix) = - solve(Double::class, a, b) { it < 1e-11 } +fun RealMatrixContext.solve(a: Matrix, b: Matrix) = + solve(a, b) { it < 1e-11 } -@ExperimentalContracts inline fun , F : Field> GenericMatrixContext.inverse( matrix: Matrix, noinline checkSingular: (T) -> Boolean -) = - solve(T::class, matrix, one(matrix.rowNum, matrix.colNum), checkSingular) +) = solve(matrix, one(matrix.rowNum, matrix.colNum), checkSingular) -@ExperimentalContracts -fun GenericMatrixContext.inverse(matrix: Matrix) = - inverse(matrix) { it < 1e-11 } \ No newline at end of file +fun RealMatrixContext.inverse(matrix: Matrix) = + solve(matrix, one(matrix.rowNum, matrix.colNum)) { it < 1e-11 } \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixContext.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixContext.kt index f97c55cc4..7797fdadf 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixContext.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixContext.kt @@ -1,6 +1,5 @@ package scientifik.kmath.linear -import scientifik.kmath.operations.RealField import scientifik.kmath.operations.Ring import scientifik.kmath.operations.SpaceOperations import scientifik.kmath.operations.sum @@ -30,7 +29,7 @@ interface MatrixContext : SpaceOperations> { /** * Non-boxing double matrix */ - val real = BufferMatrixContext(RealField, Buffer.Companion::auto) + val real = RealMatrixContext /** * A structured matrix with custom buffer diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferAccessor2D.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferAccessor2D.kt new file mode 100644 index 000000000..ac0727f7c --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferAccessor2D.kt @@ -0,0 +1,45 @@ +package scientifik.kmath.structures + +import kotlin.reflect.KClass + +/** + * A context that allows to operate on a [MutableBuffer] as on 2d array + */ +class BufferAccessor2D(val type: KClass, val rowNum: Int, val colNum: Int) { + + inline operator fun Buffer.get(i: Int, j: Int) = get(i + colNum * j) + + inline operator fun MutableBuffer.set(i: Int, j: Int, value: T) { + set(i + colNum * j, value) + } + + inline fun create(init: (i: Int, j: Int) -> T) = + MutableBuffer.auto(type, rowNum * colNum) { offset -> init(offset / colNum, offset % colNum) } + + fun create(mat: Structure2D) = create { i, j -> mat[i, j] } + + //TODO optimize wrapper + fun MutableBuffer.collect(): Structure2D = + NDStructure.auto(type, rowNum, colNum) { (i, j) -> get(i, j) }.as2D() + + + inner class Row(val buffer: MutableBuffer, val rowIndex: Int) : MutableBuffer { + override val size: Int get() = colNum + + override fun get(index: Int): T = buffer[rowIndex, index] + + override fun set(index: Int, value: T) { + buffer[rowIndex, index] = value + } + + override fun copy(): MutableBuffer = MutableBuffer.auto(type, colNum) { get(it) } + + override fun iterator(): Iterator = (0 until colNum).map(::get).iterator() + + } + + /** + * Get row + */ + fun MutableBuffer.row(i: Int) = Row(this, i) +} diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt index 51970ebae..c8d091017 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt @@ -84,7 +84,7 @@ interface MutableBuffer : Buffer { MutableListBuffer(MutableList(size, initializer)) @Suppress("UNCHECKED_CAST") - inline fun auto(type: KClass, size: Int, initializer: (Int) -> T): MutableBuffer { + inline fun auto(type: KClass, size: Int, initializer: (Int) -> T): MutableBuffer { return when (type) { Double::class -> DoubleBuffer(DoubleArray(size) { initializer(it) as Double }) as MutableBuffer Short::class -> ShortBuffer(ShortArray(size) { initializer(it) as Short }) as MutableBuffer