Basic matrix inversion

kmath-commons/src/main/kotlin/scientifik/kmath/commons/CommonsMatrix.kt

@@ -1,3 +0,0 @@
-package scientifik.kmath.commons
-
-//val solver: DecompositionSolver

kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt

@@ -1,44 +1,12 @@
 package scientifik.kmath.linear
 
-import scientifik.kmath.operations.Field
 import scientifik.kmath.structures.MutableNDArray
 import scientifik.kmath.structures.NDArray
 import scientifik.kmath.structures.NDArrays
+import kotlin.math.absoluteValue
 
 /**
- * Calculates the LUP-decomposition of a square matrix.
+ * Implementation copier from Apache common-maths
- *
- * The LUP-decomposition of a matrix A consists of three matrices L, U and
- * P that satisfy: PA = LU. L is lower triangular (with unit
- * diagonal terms), U is upper triangular and P is a permutation matrix. All
- * matrices are mm.
- *
- * As shown by the presence of the P matrix, this decomposition is
- * implemented using partial pivoting.
- *
- * This class is based on the class with similar name from the
- * [JAMA](http://math.nist.gov/javanumerics/jama/) library.
- *
- *  * a [getP][.getP] method has been added,
- *  * the `det` method has been renamed as [   getDeterminant][.getDeterminant],
- *  * the `getDoublePivot` method has been removed (but the int based
- * [getPivot][.getPivot] method has been kept),
- *  * the `solve` and `isNonSingular` methods have been replaced
- * by a [getSolver][.getSolver] method and the equivalent methods
- * provided by the returned [DecompositionSolver].
- *
- *
- * @see [MathWorld](http://mathworld.wolfram.com/LUDecomposition.html)
- *
- * @see [Wikipedia](http://en.wikipedia.org/wiki/LU_decomposition)
- *
- * @since 2.0 (changed to concrete class in 3.0)
- *
- * @param matrix The matrix to decompose.
- * @param singularityThreshold threshold (based on partial row norm)
- * under which a matrix is considered singular
- * @throws NonSquareMatrixException if matrix is not square
-
  */
 abstract class LUDecomposition<T : Comparable<T>>(val matrix: Matrix<T>) {
 
@@ -51,7 +19,7 @@ abstract class LUDecomposition<T : Comparable<T>>(val matrix: Matrix<T>) {
     private var even: Boolean = false
 
     init {
-        val pair = matrix.context.field.calculateLU()
+        val pair = calculateLU()
         lu = pair.first
         pivot = pair.second
     }
@@ -134,72 +102,76 @@ abstract class LUDecomposition<T : Comparable<T>>(val matrix: Matrix<T>) {
 
     abstract fun isSingular(value: T): Boolean
 
-    private fun Field<T>.calculateLU(): Pair<NDArray<T>, IntArray> {
+    private fun abs(value: T) = if (value > matrix.context.field.zero) value else with(matrix.context.field) { -value }
+
+    private fun calculateLU(): Pair<NDArray<T>, IntArray> {
         if (matrix.rows != matrix.columns) {
             error("LU decomposition supports only square matrices")
         }
 
-        fun T.abs() = if (this > zero) this else -this
-
         val m = matrix.columns
         val pivot = IntArray(matrix.rows)
         //TODO fix performance
         val lu: MutableNDArray<T> = NDArrays.createMutable(matrix.context.field, listOf(matrix.rows, matrix.columns)) { index -> matrix[index[0], index[1]] }
 
-        // Initialize permutation array and parity
-        for (row in 0 until m) {
-            pivot[row] = row
-        }
-        even = true
 
-        // Loop over columns
+        with(matrix.context.field) {
-        for (col in 0 until m) {
+            // Initialize permutation array and parity
-
+            for (row in 0 until m) {
-            // upper
+                pivot[row] = row
-            for (row in 0 until col) {
-                var sum = lu[row, col]
-                for (i in 0 until row) {
-                    sum -= lu[row, i] * lu[i, col]
-                }
-                lu[row, col] = sum
             }
+            even = true
 
-            // lower
+            // Loop over columns
-            val max = (col until m).maxBy { row ->
+            for (col in 0 until m) {
-                var sum = lu[row, col]
+
-                for (i in 0 until col) {
+                // upper
-                    sum -= lu[row, i] * lu[i, col]
+                for (row in 0 until col) {
+                    var sum = lu[row, col]
+                    for (i in 0 until row) {
+                        sum -= lu[row, i] * lu[i, col]
+                    }
+                    lu[row, col] = sum
                 }
-                //luRow[col] = sum
-                lu[row, col] = sum
 
-                sum.abs()
+                // lower
-            } ?: col
+                val max = (col until m).maxBy { row ->
+                    var sum = lu[row, col]
+                    for (i in 0 until col) {
+                        sum -= lu[row, i] * lu[i, col]
+                    }
+                    //luRow[col] = sum
+                    lu[row, col] = sum
 
-            // Singularity check
+                    abs(sum)
-            if (isSingular(lu[max, col].abs())) {
+                } ?: col
-                error("Singular matrix")
-            }
 
-            // Pivot if necessary
+                // Singularity check
-            if (max != col) {
+                if (isSingular(lu[max, col])) {
-                //var tmp = zero
+                    error("Singular matrix")
-                //val luMax = lu[max]
-                //val luCol = lu[col]
-                for (i in 0 until m) {
-                    lu[max, i] = lu[col, i]
-                    lu[col, i] = lu[max, i]
                 }
-                val temp = pivot[max]
-                pivot[max] = pivot[col]
-                pivot[col] = temp
-                even = !even
-            }
 
-            // Divide the lower elements by the "winning" diagonal elt.
+                // Pivot if necessary
-            val luDiag = lu[col, col]
+                if (max != col) {
-            for (row in col + 1 until m) {
+                    //var tmp = zero
-                lu[row, col] /= luDiag
+                    //val luMax = lu[max]
+                    //val luCol = lu[col]
+                    for (i in 0 until m) {
+                        lu[max, i] = lu[col, i]
+                        lu[col, i] = lu[max, i]
+                    }
+                    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] = lu[row, col] / luDiag
+//                    lu[row, col] /= luDiag
+                }
             }
         }
         return Pair(lu, pivot)
@@ -214,21 +186,22 @@ abstract class LUDecomposition<T : Comparable<T>>(val matrix: Matrix<T>) {
         return pivot.copyOf()
     }
 
+}
+
+class RealLUDecomposition(matrix: Matrix<Double>, private val singularityThreshold: Double = DEFAULT_TOO_SMALL) : LUDecomposition<Double>(matrix) {
+    override fun isSingular(value: Double): Boolean {
+        return value.absoluteValue < singularityThreshold
+    }
+
     companion object {
         /** Default bound to determine effective singularity in LU decomposition.  */
         private const val DEFAULT_TOO_SMALL = 1e-11
     }
 }
 
-class RealLUDecomposition(matrix: Matrix<Double>, private val singularityThreshold: Double = 1e-11) : LUDecomposition<Double>(matrix) {
-    override fun isSingular(value: Double): Boolean {
-        return value < singularityThreshold
-    }
-}
-
 
 /** Specialized solver.  */
-class RealLUSolver : LinearSolver<Double> {
+object RealLUSolver : LinearSolver<Double> {
 
 //
 //    /** {@inheritDoc}  */

kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt

@@ -14,31 +14,31 @@ import scientifik.kmath.structures.realNDFieldFactory
  * @param T type of individual element of the vector or matrix
  * @param V the type of vector space element
  */
-abstract class LinearSpace<T : Any, V : Matrix<T>>(val rows: Int, val columns: Int, val field: Field<T>) : Space<V> {
+abstract class MatrixSpace<T : Any>(val rows: Int, val columns: Int, val field: Field<T>) : Space<Matrix<T>> {
 
     /**
      * Produce the element of this space
      */
-    abstract fun produce(initializer: (Int, Int) -> T): V
+    abstract fun produce(initializer: (Int, Int) -> T): Matrix<T>
 
     /**
-     * Produce new linear space with given dimensions. The space produced could be raised from cache since [LinearSpace] does not have mutable elements
+     * Produce new matrix space with given dimensions. The space produced could be raised from cache since [MatrixSpace] does not have mutable elements
      */
-    abstract fun produceSpace(rows: Int, columns: Int): LinearSpace<T, V>
+    abstract fun produceSpace(rows: Int, columns: Int): MatrixSpace<T>
 
-    override val zero: V by lazy {
+    override val zero: Matrix<T> by lazy {
         produce { _, _ -> field.zero }
     }
 
-    val one: V by lazy {
+//    val one: Matrix<T> by lazy {
-        produce { i, j -> if (i == j) field.one else field.zero }
+//        produce { i, j -> if (i == j) field.one else field.zero }
-    }
+//    }
 
-    override fun add(a: V, b: V): V {
+    override fun add(a: Matrix<T>, b: Matrix<T>): Matrix<T> {
         return produce { i, j -> with(field) { a[i, j] + b[i, j] } }
     }
 
-    override fun multiply(a: V, k: Double): V {
+    override fun multiply(a: Matrix<T>, k: Double): Matrix<T> {
         //TODO it is possible to implement scalable linear elements which normed values and adjustable scale to save memory and processing poser
         return produce { i, j -> with(field) { a[i, j] * k } }
     }
@@ -46,29 +46,23 @@ abstract class LinearSpace<T : Any, V : Matrix<T>>(val rows: Int, val columns: I
     /**
      * Dot product. Throws exception on dimension mismatch
      */
-    fun multiply(a: V, b: V): V {
+    fun multiply(a: Matrix<T>, b: Matrix<T>): Matrix<T> {
         if (a.rows != b.columns) {
             //TODO replace by specific exception
             error("Dimension mismatch in linear structure dot product: [${a.rows},${a.columns}]*[${b.rows},${b.columns}]")
         }
         return produceSpace(a.rows, b.columns).produce { i, j ->
-            (0..a.columns).asSequence().map { k -> field.multiply(a[i, k], b[k, j]) }.reduce { first, second -> field.add(first, second) }
+            (0 until a.columns).asSequence().map { k -> field.multiply(a[i, k], b[k, j]) }.reduce { first, second -> field.add(first, second) }
         }
     }
-
-    infix fun V.dot(b: V): V = multiply(this, b)
 }
 
-/**
+infix fun <T : Any> Matrix<T>.dot(b: Matrix<T>): Matrix<T> = this.context.multiply(this, b)
- * A specialized [LinearSpace] which works with vectors
- */
-abstract class VectorSpace<T : Any, V : Vector<T>>(size: Int, field: Field<T>) : LinearSpace<T, V>(size, 1, field)
 
 /**
  * A matrix-like structure
  */
-interface Matrix<T : Any> {
+interface Matrix<T : Any> : SpaceElement<Matrix<T>, MatrixSpace<T>> {
-    val context: LinearSpace<T, out Matrix<T>>
     /**
      * Number of rows
      */
@@ -83,34 +77,58 @@ interface Matrix<T : Any> {
      */
     operator fun get(i: Int, j: Int): T
 
+    override val self: Matrix<T>
+        get() = this
+
     fun transpose(): Matrix<T> {
         return object : Matrix<T> {
-            override val context: LinearSpace<T, out Matrix<T>> = this@Matrix.context
+            override val context: MatrixSpace<T> = this@Matrix.context
             override val rows: Int = this@Matrix.columns
             override val columns: Int = this@Matrix.rows
             override fun get(i: Int, j: Int): T = this@Matrix[j, i]
         }
     }
+
+    companion object {
+        fun <T : Any> one(rows: Int, columns: Int, field: Field<T>): Matrix<T> {
+            return matrix(rows, columns, field) { i, j -> if (i == j) field.one else field.zero }
+        }
+    }
 }
 
-interface Vector<T : Any> : Matrix<T> {
-    override val context: VectorSpace<T, Vector<T>>
-    override val columns: Int
-        get() = 1
 
-    operator fun get(i: Int) = get(i, 0)
+/**
+ * A linear space for vectors
+ */
+abstract class VectorSpace<T : Any>(val size: Int, val field: Field<T>) : Space<Vector<T>> {
+
+    abstract fun produce(initializer: (Int) -> T): Vector<T>
+
+    override val zero: Vector<T> by lazy { produce { field.zero } }
+
+    override fun add(a: Vector<T>, b: Vector<T>): Vector<T> = produce { with(field) { a[it] + b[it] } }
+
+    override fun multiply(a: Vector<T>, k: Double): Vector<T> = produce { with(field) { a[it] * k } }
+}
+
+
+interface Vector<T : Any> : SpaceElement<Vector<T>, VectorSpace<T>> {
+    val size: Int
+        get() = context.size
+
+    operator fun get(i: Int): T
 }
 
 
 /**
  * NDArray-based implementation of vector space. By default uses slow [SimpleNDField], but could be overridden with custom [NDField] factory.
  */
-class ArraySpace<T : Any>(
+class ArrayMatrixSpace<T : Any>(
         rows: Int,
         columns: Int,
         field: Field<T>,
         val ndFactory: NDFieldFactory<T> = createFactory(field)
-) : LinearSpace<T, Matrix<T>>(rows, columns, field) {
+) : MatrixSpace<T>(rows, columns, field) {
 
     val ndField by lazy {
         ndFactory(listOf(rows, columns))
@@ -118,8 +136,8 @@ class ArraySpace<T : Any>(
 
     override fun produce(initializer: (Int, Int) -> T): Matrix<T> = ArrayMatrix(this, initializer)
 
-    override fun produceSpace(rows: Int, columns: Int): ArraySpace<T> {
+    override fun produceSpace(rows: Int, columns: Int): ArrayMatrixSpace<T> {
-        return ArraySpace(rows, columns, field, ndFactory)
+        return ArrayMatrixSpace(rows, columns, field, ndFactory)
     }
 }
 
@@ -127,26 +145,20 @@ class ArrayVectorSpace<T : Any>(
         size: Int,
         field: Field<T>,
         val ndFactory: NDFieldFactory<T> = createFactory(field)
-) : VectorSpace<T, Vector<T>>(size, field) {
+) : VectorSpace<T>(size, field) {
     val ndField by lazy {
         ndFactory(listOf(size))
     }
 
-    override fun produce(initializer: (Int, Int) -> T): Vector<T> = produceVector { i -> initializer(i, 0) }
+    override fun produce(initializer: (Int) -> T): Vector<T> = ArrayVector(this, initializer)
-
-    fun produceVector(initializer: (Int) -> T): Vector<T> = ArrayVector(this, initializer)
-
-    override fun produceSpace(rows: Int, columns: Int): LinearSpace<T, Vector<T>> {
-        TODO("not implemented") //To change body of created functions use File | Settings | File Templates.
-    }
 }
 
 /**
- * Member of [ArraySpace] which wraps 2-D array
+ * Member of [ArrayMatrixSpace] which wraps 2-D array
  */
-class ArrayMatrix<T : Any> internal constructor(override val context: ArraySpace<T>, val array: NDArray<T>) : Matrix<T>, SpaceElement<Matrix<T>, ArraySpace<T>> {
+class ArrayMatrix<T : Any> internal constructor(override val context: ArrayMatrixSpace<T>, val array: NDArray<T>) : Matrix<T> {
 
-    constructor(context: ArraySpace<T>, initializer: (Int, Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0], list[1]) })
+    constructor(context: ArrayMatrixSpace<T>, initializer: (Int, Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0], list[1]) })
 
     override val rows: Int get() = context.rows
 
@@ -160,32 +172,21 @@ class ArrayMatrix<T : Any> internal constructor(override val context: ArraySpace
 }
 
 
-class ArrayVector<T : Any> internal constructor(override val context: ArrayVectorSpace<T>, val array: NDArray<T>) : Vector<T>, SpaceElement<Vector<T>, ArrayVectorSpace<T>> {
+class ArrayVector<T : Any> internal constructor(override val context: ArrayVectorSpace<T>, val array: NDArray<T>) : Vector<T> {
 
     constructor(context: ArrayVectorSpace<T>, initializer: (Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0]) })
 
     init {
-        if (context.columns != 1) {
+        if (context.size != array.shape[0]) {
-            error("Vector must have single column")
-        }
-        if (context.rows != array.shape[0]) {
             error("Array dimension mismatch")
         }
     }
 
-    //private val array = context.ndField.produce { list -> initializer(list[0]) }
+    override fun get(i: Int): T {
-
-
-    override val rows: Int get() = context.rows
-
-    override val columns: Int = 1
-
-    override fun get(i: Int, j: Int): T {
         return array[i]
     }
 
     override val self: ArrayVector<T> get() = this
-
 }
 
 /**
@@ -193,8 +194,8 @@ class ArrayVector<T : Any> internal constructor(override val context: ArrayVecto
  */
 interface LinearSolver<T : Any> {
     fun solve(a: Matrix<T>, b: Matrix<T>): Matrix<T>
-    fun solve(a: Matrix<T>, b: Vector<T>): Vector<T> = solve(a, b as Matrix<T>).toVector()
+    fun solve(a: Matrix<T>, b: Vector<T>): Vector<T> = solve(a, b.toMatrix()).toVector()
-    fun inverse(a: Matrix<T>): Matrix<T> = solve(a, a.context.one)
+    fun inverse(a: Matrix<T>): Matrix<T> = solve(a, Matrix.one(a.rows, a.columns, a.context.field))
 }
 
 /**
@@ -220,13 +221,13 @@ fun DoubleArray.asVector() = realVector(this.size) { this[it] }
  * Create [ArrayMatrix] with custom field
  */
 fun <T : Any> matrix(rows: Int, columns: Int, field: Field<T>, initializer: (Int, Int) -> T) =
-        ArrayMatrix(ArraySpace(rows, columns, field), initializer)
+        ArrayMatrix(ArrayMatrixSpace(rows, columns, field), initializer)
 
 /**
  * Create [ArrayMatrix] of doubles.
  */
 fun realMatrix(rows: Int, columns: Int, initializer: (Int, Int) -> Double) =
-        ArrayMatrix(ArraySpace(rows, columns, DoubleField, realNDFieldFactory), initializer)
+        ArrayMatrix(ArrayMatrixSpace(rows, columns, DoubleField, realNDFieldFactory), initializer)
 
 
 /**
@@ -234,17 +235,28 @@ fun realMatrix(rows: Int, columns: Int, initializer: (Int, Int) -> Double) =
  */
 fun <T : Any> Matrix<T>.toVector(): Vector<T> {
     return when {
-        this is Vector -> return this
         this.columns == 1 -> {
-            if (this is ArrayMatrix) {
+//            if (this is ArrayMatrix) {
-                //Reuse existing underlying array
+//                //Reuse existing underlying array
-                ArrayVector(ArrayVectorSpace(rows, context.field, context.ndFactory), array)
+//                ArrayVector(ArrayVectorSpace(rows, context.field, context.ndFactory), array)
-            } else {
+//            } else {
-                //Generic vector
+//                //Generic vector
-                vector(rows, context.field) { get(it, 0) }
+//                vector(rows, context.field) { get(it, 0) }
-            }
+//            }
+            vector(rows, context.field) { get(it, 0) }
         }
         else -> error("Can't convert matrix with more than one column to vector")
     }
 }
 
+fun <T : Any> Vector<T>.toMatrix(): Matrix<T> {
+//    return if (this is ArrayVector) {
+//        //Reuse existing underlying array
+//        ArrayMatrix(ArrayMatrixSpace(size, 1, context.field, context.ndFactory), array)
+//    } else {
+//        //Generic vector
+//        matrix(size, 1, context.field) { i, j -> get(i) }
+//    }
+    return matrix(size, 1, context.field) { i, j -> get(i) }
+}
+

kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt

@@ -7,6 +7,11 @@ package scientifik.kmath.operations
  * @param S the type of mathematical context for this element
  */
 interface MathElement<T, S> {
+    /**
+     * Self value. Needed for static type checking.
+     */
+    val self: T
+
     /**
      * The context this element belongs to
      */
@@ -54,12 +59,6 @@ interface Space<T> {
  * @param S the type of space
  */
 interface SpaceElement<T, S : Space<T>> : MathElement<T, S> {
-
-    /**
-     * Self value. Needed for static type checking. Needed to avoid type erasure on JVM.
-     */
-    val self: T
-
     operator fun plus(b: T): T = context.add(self, b)
     operator fun minus(b: T): T = context.add(self, context.multiply(b, -1.0))
     operator fun times(k: Number): T = context.multiply(self, k.toDouble())

kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/ArrayMatrixTest.kt

@@ -1,6 +1,5 @@
-package scientifik.kmath.structures
+package scientifik.kmath.linear
 
-import scientifik.kmath.linear.realVector
 import kotlin.test.Test
 import kotlin.test.assertEquals
 
@@ -11,17 +10,25 @@ class ArrayMatrixTest {
         val vector1 = realVector(5) { it.toDouble() }
         val vector2 = realVector(5) { 5 - it.toDouble() }
         val sum = vector1 + vector2
-        assertEquals(5.0, sum[2, 0])
+        assertEquals(5.0, sum[2])
     }
 
+    @Test
+    fun testVectorToMatrix() {
+        val vector = realVector(5) { it.toDouble() }
+        val matrix = vector.toMatrix()
+        assertEquals(4.0, matrix[4, 0])
+    }
+
+
     @Test
     fun testDot() {
         val vector1 = realVector(5) { it.toDouble() }
         val vector2 = realVector(5) { 5 - it.toDouble() }
-        val product = with(vector1.context) {
+        val product = vector1.toMatrix() dot (vector2.toMatrix().transpose())
-            vector1 dot (vector2.transpose())
-        }
 
-        assertEquals(10.0, product[1, 0])
+
+        assertEquals(5.0, product[1, 0])
+        assertEquals(6.0, product[2, 2])
     }
 }

kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/RealLUSolverTest.kt

@@ -0,0 +1,14 @@
+package scientifik.kmath.linear
+
+import scientifik.kmath.operations.DoubleField
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+class RealLUSolverTest {
+    @Test
+    fun testInvert() {
+        val matrix = Matrix.one(2, 2, DoubleField)
+        val inverted = RealLUSolver.inverse(matrix)
+        assertEquals(1.0, inverted[0, 0])
+    }
+}

settings.gradle

@@ -10,5 +10,4 @@ enableFeaturePreview('GRADLE_METADATA')
 
 rootProject.name = 'kmath'
 include ':kmath-core'
-include ':kmath-commons'