From a2a7ddcddaa400ccb7b6533fe4021a9583fbbf22 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Thu, 17 Jan 2019 19:26:17 +0300 Subject: [PATCH] Commons-math implementation for linear algebra --- benchmarks/build.gradle | 1 + .../kmath/linear/LinearAlgebraBenchmark.kt | 48 +++++++++++ .../scientifik/kmath/linear/CMMatrix.kt | 83 +++++++++++++++++++ .../kmath/linear/LUPDecomposition.kt | 28 +++---- .../scientifik/kmath/linear/LinearAlgrebra.kt | 5 +- .../kotlin/scientifik/kmath/linear/Matrix.kt | 4 +- 6 files changed, 149 insertions(+), 20 deletions(-) create mode 100644 benchmarks/src/main/kotlin/scientifik/kmath/linear/LinearAlgebraBenchmark.kt create mode 100644 kmath-commons/src/main/kotlin/scientifik/kmath/linear/CMMatrix.kt diff --git a/benchmarks/build.gradle b/benchmarks/build.gradle index a9450a974..f732e2db6 100644 --- a/benchmarks/build.gradle +++ b/benchmarks/build.gradle @@ -7,6 +7,7 @@ plugins { dependencies { compile project(":kmath-core") compile project(":kmath-coroutines") + compile project(":kmath-commons") //jmh project(':kmath-core') } diff --git a/benchmarks/src/main/kotlin/scientifik/kmath/linear/LinearAlgebraBenchmark.kt b/benchmarks/src/main/kotlin/scientifik/kmath/linear/LinearAlgebraBenchmark.kt new file mode 100644 index 000000000..1a12d25b1 --- /dev/null +++ b/benchmarks/src/main/kotlin/scientifik/kmath/linear/LinearAlgebraBenchmark.kt @@ -0,0 +1,48 @@ +package scientifik.kmath.linear + +import org.apache.commons.math3.linear.Array2DRowRealMatrix +import kotlin.random.Random +import kotlin.system.measureTimeMillis + +fun main() { + val random = Random(12224) + val dim = 100 + //creating invertible matrix + val u = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 } + val l = Matrix.real(dim, dim) { i, j -> if (i >= j) random.nextDouble() else 0.0 } + val matrix = l dot u + + val n = 500 // iterations + + val solver = LUSolver.real + + repeat(50) { + val res = solver.inverse(matrix) + } + + val inverseTime = measureTimeMillis { + repeat(n) { + val res = solver.inverse(matrix) + } + } + + println("[kmath] Inversion of $n matrices $dim x $dim finished in $inverseTime millis") + + //commons + + val cmSolver = CMSolver + + val commonsMatrix = Array2DRowRealMatrix(dim, dim) + matrix.elements().forEach { (index, value) -> commonsMatrix.setEntry(index[0], index[1], value) } + + val commonsTime = measureTimeMillis { + val cm = matrix.toCM() + repeat(n) { + //overhead on coversion could be mitigated + val res = cmSolver.inverse(cm) + } + } + + + println("[commons-math] Inversion of $n matrices $dim x $dim finished in $commonsTime millis") +} \ No newline at end of file diff --git a/kmath-commons/src/main/kotlin/scientifik/kmath/linear/CMMatrix.kt b/kmath-commons/src/main/kotlin/scientifik/kmath/linear/CMMatrix.kt new file mode 100644 index 000000000..a3f206c54 --- /dev/null +++ b/kmath-commons/src/main/kotlin/scientifik/kmath/linear/CMMatrix.kt @@ -0,0 +1,83 @@ +package scientifik.kmath.linear + +import org.apache.commons.math3.linear.* +import org.apache.commons.math3.linear.RealMatrix +import org.apache.commons.math3.linear.RealVector +import scientifik.kmath.operations.RealField +import scientifik.kmath.structures.DoubleBuffer + +inline class CMMatrix(val origin: RealMatrix) : Matrix { + override val rowNum: Int get() = origin.rowDimension + override val colNum: Int get() = origin.columnDimension + + override val features: Set get() = emptySet() + + override fun get(i: Int, j: Int): Double = origin.getEntry(i, j) +} + +fun Matrix.toCM(): CMMatrix = if (this is CMMatrix) { + this +} else { + //TODO add feature analysis + val array = Array(rowNum) { i -> DoubleArray(colNum) { j -> get(i, j) } } + CMMatrix(Array2DRowRealMatrix(array)) +} + +fun RealMatrix.toMatrix() = CMMatrix(this) + +inline class CMVector(val origin: RealVector) : Point { + override val size: Int get() = origin.dimension + + override fun get(index: Int): Double = origin.getEntry(index) + + override fun iterator(): Iterator = origin.toArray().iterator() +} + +fun Point.toCM(): CMVector = if (this is CMVector) { + this +} else { + val array = DoubleArray(size) { this[it] } + CMVector(ArrayRealVector(array)) +} + +fun RealVector.toPoint() = DoubleBuffer(this.toArray()) + +object CMMatrixContext : MatrixContext { + override val elementContext: RealField get() = RealField + + override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> Double): Matrix { + val array = Array(rows) { i -> DoubleArray(columns) { j -> initializer(i, j) } } + return CMMatrix(Array2DRowRealMatrix(array)) + } + + override fun point(size: Int, initializer: (Int) -> Double): Point { + val array = DoubleArray(size, initializer) + return CMVector(ArrayRealVector(array)) + } + + override fun Matrix.dot(other: Matrix): Matrix = this.toCM().dot(other.toCM()) +} + +operator fun CMMatrix.plus(other: CMMatrix): CMMatrix = CMMatrix(this.origin.add(other.origin)) +operator fun CMMatrix.minus(other: CMMatrix): CMMatrix = CMMatrix(this.origin.subtract(other.origin)) + +infix fun CMMatrix.dot(other: CMMatrix): CMMatrix = CMMatrix(this.origin.multiply(other.origin)) + +object CMSolver : LinearSolver { + + override fun solve(a: Matrix, b: Matrix): Matrix { + val decomposition = LUDecomposition(a.toCM().origin) + return decomposition.solver.solve(b.toCM().origin).toMatrix() + } + + override fun solve(a: Matrix, b: Point): Point { + val decomposition = LUDecomposition(a.toCM().origin) + return decomposition.solver.solve(b.toCM().origin).toPoint() + } + + override fun inverse(a: Matrix): Matrix { + val decomposition = LUDecomposition(a.toCM().origin) + return decomposition.solver.inverse.toMatrix() + } + +} \ No newline at end of file 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 dbc45618d..79973bd37 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUPDecomposition.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUPDecomposition.kt @@ -2,13 +2,16 @@ package scientifik.kmath.linear import scientifik.kmath.operations.Field import scientifik.kmath.operations.Ring -import scientifik.kmath.structures.* +import scientifik.kmath.structures.MutableBuffer import scientifik.kmath.structures.MutableBuffer.Companion.boxing +import scientifik.kmath.structures.MutableBufferFactory +import scientifik.kmath.structures.NDStructure +import scientifik.kmath.structures.get class LUPDecomposition>( private val elementContext: Ring, - private val lu: NDStructure, + internal val lu: NDStructure, val pivot: IntArray, private val even: Boolean ) : DeterminantFeature { @@ -62,19 +65,12 @@ class LUPDecomposition>( class LUSolver, F : Field>( - override val context: MatrixContext, + val context: MatrixContext, val bufferFactory: MutableBufferFactory = ::boxing, val singularityCheck: (T) -> Boolean ) : LinearSolver { - /** - * In-place transformation for [MutableNDStructure], using given transformation for each element - */ - private operator fun MutableNDStructure.set(i: Int, j: Int, value: T) { - this[intArrayOf(i, j)] = value - } - private fun abs(value: T) = if (value > context.elementContext.zero) value else with(context.elementContext) { -value } @@ -180,19 +176,19 @@ class LUSolver, F : Field>( for (col in 0 until a.rowNum) { for (i in col + 1 until a.rowNum) { for (j in 0 until b.colNum) { - bp[i, j] -= bp[col, j] * l[i, col] + bp[i, j] -= bp[col, j] * lu[i, col] } } } // Solve UX = Y for (col in a.rowNum - 1 downTo 0) { - for(j in 0 until b.colNum){ - bp[col,j]/= u[col,col] + for (j in 0 until b.colNum) { + bp[col, j] /= lu[col, col] } for (i in 0 until col) { for (j in 0 until b.colNum) { - bp[i, j] -= bp[col, j] * u[i, col] + bp[i, j] -= bp[col, j] * lu[i, col] } } } @@ -202,7 +198,9 @@ class LUSolver, F : Field>( } } + override fun inverse(a: Matrix): Matrix = solve(a, context.one(a.rowNum, a.colNum)) + companion object { val real = LUSolver(MatrixContext.real, MutableBuffer.Companion::auto) { it < 1e-11 } } -} +} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt index fb542c16b..4c56b1e63 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt @@ -12,12 +12,9 @@ import scientifik.kmath.structures.asSequence * A group of methods to resolve equation A dot X = B, where A and B are matrices or vectors */ interface LinearSolver> { - val context: MatrixContext - - fun solve(a: Matrix, b: Matrix): Matrix fun solve(a: Matrix, b: Point): Point = solve(a, b.toMatrix()).toVector() - fun inverse(a: Matrix): Matrix = solve(a, context.one(a.rowNum, a.colNum)) + fun inverse(a: Matrix): Matrix } /** diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt index 6cd5e985f..7c98502f3 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt @@ -199,4 +199,6 @@ fun > Matrix.transpose(): Matrix { this.rowNum, setOf(TransposedFeature(this)) ) { i, j -> get(j, i) } -} \ No newline at end of file +} + +infix fun Matrix.dot(other: Matrix): Matrix = with(MatrixContext.real) { dot(other) } \ No newline at end of file