From b288704528ca54ace552369ca068cd1bf1983992 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Thu, 7 Jan 2021 18:07:00 +0300 Subject: [PATCH] Optimize RealMatrix dot operation --- ...iplicationBenchmark.kt => DotBenchmark.kt} | 24 ++++-- .../kmath/structures/typeSafeDimensions.kt | 5 +- .../kscience/kmath/commons/linear/CMMatrix.kt | 1 + .../kscience/kmath/linear/BufferMatrix.kt | 24 +----- .../kscience/kmath/linear/LUPDecomposition.kt | 11 +-- .../kscience/kmath/linear/MatrixContext.kt | 14 ++-- .../kmath/linear/RealMatrixContext.kt | 84 +++++++++++++++++++ .../kscience/kmath/dimensions/Wrappers.kt | 38 ++++----- .../kscience/dimensions/DMatrixContextTest.kt | 1 + .../kotlin/kscience/kmath/real/RealMatrix.kt | 11 +-- 10 files changed, 135 insertions(+), 78 deletions(-) rename examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/{MultiplicationBenchmark.kt => DotBenchmark.kt} (73%) create mode 100644 kmath-core/src/commonMain/kotlin/kscience/kmath/linear/RealMatrixContext.kt diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/MultiplicationBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/DotBenchmark.kt similarity index 73% rename from examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/MultiplicationBenchmark.kt rename to examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/DotBenchmark.kt index 9d2b02245..8823e86db 100644 --- a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/MultiplicationBenchmark.kt +++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/DotBenchmark.kt @@ -2,19 +2,22 @@ package kscience.kmath.benchmarks import kotlinx.benchmark.Benchmark import kscience.kmath.commons.linear.CMMatrixContext -import kscience.kmath.commons.linear.CMMatrixContext.dot import kscience.kmath.commons.linear.toCM import kscience.kmath.ejml.EjmlMatrixContext import kscience.kmath.ejml.toEjml +import kscience.kmath.linear.BufferMatrixContext +import kscience.kmath.linear.RealMatrixContext import kscience.kmath.linear.real +import kscience.kmath.operations.RealField import kscience.kmath.operations.invoke +import kscience.kmath.structures.Buffer import kscience.kmath.structures.Matrix import org.openjdk.jmh.annotations.Scope import org.openjdk.jmh.annotations.State import kotlin.random.Random @State(Scope.Benchmark) -class MultiplicationBenchmark { +class DotBenchmark { companion object { val random = Random(12224) val dim = 1000 @@ -32,14 +35,14 @@ class MultiplicationBenchmark { @Benchmark fun commonsMathMultiplication() { - CMMatrixContext.invoke { + CMMatrixContext { cmMatrix1 dot cmMatrix2 } } @Benchmark fun ejmlMultiplication() { - EjmlMatrixContext.invoke { + EjmlMatrixContext { ejmlMatrix1 dot ejmlMatrix2 } } @@ -48,13 +51,22 @@ class MultiplicationBenchmark { fun ejmlMultiplicationwithConversion() { val ejmlMatrix1 = matrix1.toEjml() val ejmlMatrix2 = matrix2.toEjml() - EjmlMatrixContext.invoke { + EjmlMatrixContext { ejmlMatrix1 dot ejmlMatrix2 } } @Benchmark fun bufferedMultiplication() { - matrix1 dot matrix2 + BufferMatrixContext(RealField, Buffer.Companion::real).invoke{ + matrix1 dot matrix2 + } + } + + @Benchmark + fun realMultiplication(){ + RealMatrixContext { + matrix1 dot matrix2 + } } } \ No newline at end of file diff --git a/examples/src/main/kotlin/kscience/kmath/structures/typeSafeDimensions.kt b/examples/src/main/kotlin/kscience/kmath/structures/typeSafeDimensions.kt index 987eea16f..96684f7dc 100644 --- a/examples/src/main/kotlin/kscience/kmath/structures/typeSafeDimensions.kt +++ b/examples/src/main/kotlin/kscience/kmath/structures/typeSafeDimensions.kt @@ -4,9 +4,8 @@ import kscience.kmath.dimensions.D2 import kscience.kmath.dimensions.D3 import kscience.kmath.dimensions.DMatrixContext import kscience.kmath.dimensions.Dimension -import kscience.kmath.operations.RealField -private fun DMatrixContext.simple() { +private fun DMatrixContext.simple() { val m1 = produce { i, j -> (i + j).toDouble() } val m2 = produce { i, j -> (i + j).toDouble() } @@ -18,7 +17,7 @@ private object D5 : Dimension { override val dim: UInt = 5u } -private fun DMatrixContext.custom() { +private fun DMatrixContext.custom() { val m1 = produce { i, j -> (i + j).toDouble() } val m2 = produce { i, j -> (i - j).toDouble() } val m3 = produce { i, j -> (i - j).toDouble() } diff --git a/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt b/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt index 712927400..49888f8d6 100644 --- a/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt +++ b/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt @@ -29,6 +29,7 @@ public class CMMatrix(public val origin: RealMatrix, features: Set.toCM(): CMMatrix = if (this is CMMatrix) { this } else { diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt index 8b50bbe33..402161207 100644 --- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt +++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt @@ -1,8 +1,10 @@ package kscience.kmath.linear -import kscience.kmath.operations.RealField import kscience.kmath.operations.Ring -import kscience.kmath.structures.* +import kscience.kmath.structures.Buffer +import kscience.kmath.structures.BufferFactory +import kscience.kmath.structures.NDStructure +import kscience.kmath.structures.asSequence /** * Basic implementation of Matrix space based on [NDStructure] @@ -21,24 +23,6 @@ public class BufferMatrixContext>( public companion object } -@Suppress("OVERRIDE_BY_INLINE") -public object RealMatrixContext : GenericMatrixContext> { - public override val elementContext: RealField - get() = RealField - - public override inline fun produce( - rows: Int, - columns: Int, - initializer: (i: Int, j: Int) -> Double, - ): BufferMatrix { - val buffer = RealBuffer(rows * columns) { offset -> initializer(offset / columns, offset % columns) } - return BufferMatrix(rows, columns, buffer) - } - - public override inline fun point(size: Int, initializer: (Int) -> Double): Point = - RealBuffer(size, initializer) -} - public class BufferMatrix( public override val rowNum: Int, public override val colNum: Int, diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt index 099fa1909..bf2a9f59e 100644 --- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt +++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt @@ -213,17 +213,8 @@ public inline fun , F : Field> GenericMatrixContext return decomposition.solveWithLUP(bufferFactory, b) } -public fun RealMatrixContext.solveWithLUP(a: Matrix, b: Matrix): FeaturedMatrix = - solveWithLUP(a, b) { it < 1e-11 } - public inline fun , F : Field> GenericMatrixContext>.inverseWithLUP( matrix: Matrix, noinline bufferFactory: MutableBufferFactory = MutableBuffer.Companion::auto, noinline checkSingular: (T) -> Boolean, -): FeaturedMatrix = solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum), bufferFactory, checkSingular) - -/** - * Inverses a square matrix using LUP decomposition. Non square matrix will throw a error. - */ -public fun RealMatrixContext.inverseWithLUP(matrix: Matrix): FeaturedMatrix = - solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum), Buffer.Companion::real) { it < 1e-11 } +): FeaturedMatrix = solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum), bufferFactory, checkSingular) \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt index d9dc57b0f..9bc79e12b 100644 --- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt +++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt @@ -18,6 +18,11 @@ public interface MatrixContext> : SpaceOperations T): M + /** + * Produce a point compatible with matrix space (and possibly optimized for it) + */ + public fun point(size: Int, initializer: (Int) -> T): Point = Buffer.boxing(size, initializer) + @Suppress("UNCHECKED_CAST") public override fun binaryOperation(operation: String, left: Matrix, right: Matrix): M = when (operation) { "dot" -> left dot right @@ -61,10 +66,6 @@ public interface MatrixContext> : SpaceOperations): M = m * this public companion object { - /** - * Non-boxing double matrix - */ - public val real: RealMatrixContext = RealMatrixContext /** * A structured matrix with custom buffer @@ -88,11 +89,6 @@ public interface GenericMatrixContext, out M : Matrix> : */ public val elementContext: R - /** - * Produce a point compatible with matrix space - */ - public fun point(size: Int, initializer: (Int) -> T): Point - public override infix fun Matrix.dot(other: Matrix): M { //TODO add typed error require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" } diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/RealMatrixContext.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/RealMatrixContext.kt new file mode 100644 index 000000000..772b20f3b --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/RealMatrixContext.kt @@ -0,0 +1,84 @@ +package kscience.kmath.linear + +import kscience.kmath.operations.RealField +import kscience.kmath.structures.Matrix +import kscience.kmath.structures.MutableBuffer +import kscience.kmath.structures.MutableBufferFactory +import kscience.kmath.structures.RealBuffer + +@Suppress("OVERRIDE_BY_INLINE") +public object RealMatrixContext : MatrixContext> { + + public override inline fun produce( + rows: Int, + columns: Int, + initializer: (i: Int, j: Int) -> Double, + ): BufferMatrix { + val buffer = RealBuffer(rows * columns) { offset -> initializer(offset / columns, offset % columns) } + return BufferMatrix(rows, columns, buffer) + } + + private fun Matrix.wrap(): BufferMatrix = if (this is BufferMatrix) this else { + produce(rowNum, colNum) { i, j -> get(i, j) } + } + + public fun one(rows: Int, columns: Int): FeaturedMatrix = VirtualMatrix(rows, columns, DiagonalFeature) { i, j -> + if (i == j) 1.0 else 0.0 + } + + public override infix fun Matrix.dot(other: Matrix): BufferMatrix { + require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" } + return produce(rowNum, other.colNum) { i, j -> + var res = 0.0 + for (l in 0 until colNum) { + res += get(i, l) * other.get(l, j) + } + res + } + } + + public override infix fun Matrix.dot(vector: Point): Point { + require(colNum == vector.size) { "Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})" } + return RealBuffer(rowNum) { i -> + var res = 0.0 + for (j in 0 until colNum) { + res += get(i, j) * vector[j] + } + res + } + } + + override fun add(a: Matrix, b: Matrix): BufferMatrix { + require(a.rowNum == b.rowNum) { "Row number mismatch in matrix addition. Left side: ${a.rowNum}, right side: ${b.rowNum}" } + require(a.colNum == b.colNum) { "Column number mismatch in matrix addition. Left side: ${a.colNum}, right side: ${b.colNum}" } + return produce(a.rowNum, a.colNum) { i, j -> + a[i, j] + b[i, j] + } + } + + override fun Matrix.times(value: Double): BufferMatrix = + produce(rowNum, colNum) { i, j -> get(i, j) * value } + + + override fun multiply(a: Matrix, k: Number): BufferMatrix = + produce(a.rowNum, a.colNum) { i, j -> a.get(i, j) * k.toDouble() } +} + + +/** + * Partially optimized real-valued matrix + */ +public val MatrixContext.Companion.real: RealMatrixContext get() = RealMatrixContext + +public fun RealMatrixContext.solveWithLUP(a: Matrix, b: Matrix): FeaturedMatrix { + // Use existing decomposition if it is provided by matrix + val bufferFactory: MutableBufferFactory = MutableBuffer.Companion::real + val decomposition = a.getFeature() ?: lup(bufferFactory, RealField, a) { it < 1e-11 } + return decomposition.solveWithLUP(bufferFactory, b) +} + +/** + * Inverses a square matrix using LUP decomposition. Non square matrix will throw a error. + */ +public fun RealMatrixContext.inverseWithLUP(matrix: Matrix): FeaturedMatrix = + solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum)) diff --git a/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt b/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt index 68a5dc262..0422d11b2 100644 --- a/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt +++ b/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt @@ -1,11 +1,6 @@ package kscience.kmath.dimensions -import kscience.kmath.linear.GenericMatrixContext -import kscience.kmath.linear.MatrixContext -import kscience.kmath.linear.Point -import kscience.kmath.linear.transpose -import kscience.kmath.operations.RealField -import kscience.kmath.operations.Ring +import kscience.kmath.linear.* import kscience.kmath.operations.invoke import kscience.kmath.structures.Matrix import kscience.kmath.structures.Structure2D @@ -42,7 +37,7 @@ public interface DMatrix : Structure2D { * An inline wrapper for a Matrix */ public inline class DMatrixWrapper( - private val structure: Structure2D + private val structure: Structure2D, ) : DMatrix { override val shape: IntArray get() = structure.shape override operator fun get(i: Int, j: Int): T = structure[i, j] @@ -81,7 +76,7 @@ public inline class DPointWrapper(public val point: Point) /** * Basic operations on dimension-safe matrices. Operates on [Matrix] */ -public inline class DMatrixContext>(public val context: GenericMatrixContext>) { +public inline class DMatrixContext(public val context: MatrixContext>) { public inline fun Matrix.coerce(): DMatrix { require(rowNum == Dimension.dim().toInt()) { "Row number mismatch: expected ${Dimension.dim()} but found $rowNum" @@ -115,7 +110,7 @@ public inline class DMatrixContext>(public val context: Ge } public inline infix fun DMatrix.dot( - other: DMatrix + other: DMatrix, ): DMatrix = context { this@dot dot other }.coerce() public inline infix fun DMatrix.dot(vector: DPoint): DPoint = @@ -139,18 +134,19 @@ public inline class DMatrixContext>(public val context: Ge public inline fun DMatrix.transpose(): DMatrix = context { (this@transpose as Matrix).transpose() }.coerce() - /** - * A square unit matrix - */ - public inline fun one(): DMatrix = produce { i, j -> - if (i == j) context.elementContext.one else context.elementContext.zero - } - - public inline fun zero(): DMatrix = produce { _, _ -> - context.elementContext.zero - } - public companion object { - public val real: DMatrixContext = DMatrixContext(MatrixContext.real) + public val real: DMatrixContext = DMatrixContext(MatrixContext.real) } } + + +/** + * A square unit matrix + */ +public inline fun DMatrixContext.one(): DMatrix = produce { i, j -> + if (i == j) 1.0 else 0.0 +} + +public inline fun DMatrixContext.zero(): DMatrix = produce { _, _ -> + 0.0 +} \ No newline at end of file diff --git a/kmath-dimensions/src/commonTest/kotlin/kscience/dimensions/DMatrixContextTest.kt b/kmath-dimensions/src/commonTest/kotlin/kscience/dimensions/DMatrixContextTest.kt index f44b16753..5b330fcce 100644 --- a/kmath-dimensions/src/commonTest/kotlin/kscience/dimensions/DMatrixContextTest.kt +++ b/kmath-dimensions/src/commonTest/kotlin/kscience/dimensions/DMatrixContextTest.kt @@ -3,6 +3,7 @@ package kscience.dimensions import kscience.kmath.dimensions.D2 import kscience.kmath.dimensions.D3 import kscience.kmath.dimensions.DMatrixContext +import kscience.kmath.dimensions.one import kotlin.test.Test internal class DMatrixContextTest { diff --git a/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt index e8ad835e5..772abfbed 100644 --- a/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt +++ b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt @@ -1,13 +1,7 @@ package kscience.kmath.real -import kscience.kmath.linear.FeaturedMatrix -import kscience.kmath.linear.MatrixContext -import kscience.kmath.linear.RealMatrixContext.elementContext -import kscience.kmath.linear.VirtualMatrix -import kscience.kmath.linear.inverseWithLUP +import kscience.kmath.linear.* import kscience.kmath.misc.UnstableKMathAPI -import kscience.kmath.operations.invoke -import kscience.kmath.operations.sum import kscience.kmath.structures.Buffer import kscience.kmath.structures.RealBuffer import kscience.kmath.structures.asIterable @@ -122,8 +116,7 @@ public fun RealMatrix.extractColumn(columnIndex: Int): RealMatrix = extractColumns(columnIndex..columnIndex) public fun RealMatrix.sumByColumn(): RealBuffer = RealBuffer(colNum) { j -> - val column = columns[j] - elementContext { sum(column.asIterable()) } + columns[j].asIterable().sum() } public fun RealMatrix.minByColumn(): RealBuffer = RealBuffer(colNum) { j ->