Merge pull request #180 from mipt-npm/commandertvis/matrix-decompositions

Update KDoc comments for Matrix classes, improve MatrixFeature API, implement new features with EJML matrix, delete inversion API from EJML in favor of InverseMatrixFeature, override point by EJML matrix

kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt

@@ -7,10 +7,16 @@ import kscience.kmath.structures.asBuffer
 import kotlin.math.sqrt
 
 /**
- * A 2d structure plus optional matrix-specific features
+ * A [Matrix] that holds [MatrixFeature] objects.
+ *
+ * @param T the type of items.
  */
 public interface FeaturedMatrix<T : Any> : Matrix<T> {
-    override val shape: IntArray get() = intArrayOf(rowNum, colNum)
+    public override val shape: IntArray get() = intArrayOf(rowNum, colNum)
+
+    /**
+     * The set of features this matrix possesses.
+     */
     public val features: Set<MatrixFeature>
 
     /**

kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt

@@ -4,16 +4,15 @@ import kscience.kmath.operations.*
 import kscience.kmath.structures.*
 
 /**
- * Common implementation of [LUPDecompositionFeature]
+ * Common implementation of [LupDecompositionFeature].
  */
 public class LUPDecomposition<T : Any>(
     public val context: MatrixContext<T, FeaturedMatrix<T>>,
     public val elementContext: Field<T>,
-    public val lu: Structure2D<T>,
+    public val lu: Matrix<T>,
     public val pivot: IntArray,
     private val even: Boolean,
-) : LUPDecompositionFeature<T>, DeterminantFeature<T> {
-
+) : LupDecompositionFeature<T>, DeterminantFeature<T> {
     /**
      * Returns the matrix L of the decomposition.
      *
@@ -151,7 +150,10 @@ public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext
 public fun MatrixContext<Double, FeaturedMatrix<Double>>.lup(matrix: Matrix<Double>): LUPDecomposition<Double> =
     lup(Buffer.Companion::real, RealField, matrix) { it < 1e-11 }
 
-public fun <T : Any> LUPDecomposition<T>.solveWithLUP(factory: MutableBufferFactory<T>, matrix: Matrix<T>): FeaturedMatrix<T> {
+public fun <T : Any> LUPDecomposition<T>.solveWithLUP(
+    factory: MutableBufferFactory<T>,
+    matrix: Matrix<T>
+): FeaturedMatrix<T> {
     require(matrix.rowNum == pivot.size) { "Matrix dimension mismatch. Expected ${pivot.size}, but got ${matrix.colNum}" }
 
     BufferAccessor2D(matrix.rowNum, matrix.colNum, factory).run {

kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixFeatures.kt

@@ -1,62 +1,154 @@
 package kscience.kmath.linear
 
 /**
- * A marker interface representing some matrix feature like diagonal, sparse, zero, etc. Features used to optimize matrix
- * operations performance in some cases.
+ * A marker interface representing some properties of matrices or additional transformations of them. Features are used
+ * to optimize matrix operations performance in some cases or retrieve the APIs.
  */
 public interface MatrixFeature
 
 /**
- * The matrix with this feature is considered to have only diagonal non-null elements
+ * Matrices with this feature are considered to have only diagonal non-null elements.
  */
 public object DiagonalFeature : MatrixFeature
 
 /**
- * Matrix with this feature has all zero elements
+ * Matrices with this feature have all zero elements.
  */
 public object ZeroFeature : MatrixFeature
 
 /**
- * Matrix with this feature have unit elements on diagonal and zero elements in all other places
+ * Matrices with this feature have unit elements on diagonal and zero elements in all other places.
  */
 public object UnitFeature : MatrixFeature
 
 /**
- * Inverted matrix feature
+ * Matrices with this feature can be inverted: [inverse] = `a`<sup>-1</sup> where `a` is the owning matrix.
+ *
+ * @param T the type of matrices' items.
  */
 public interface InverseMatrixFeature<T : Any> : MatrixFeature {
+    /**
+     * The inverse matrix of the matrix that owns this feature.
+     */
     public val inverse: FeaturedMatrix<T>
 }
 
 /**
- * A determinant container
+ * Matrices with this feature can compute their determinant.
  */
 public interface DeterminantFeature<T : Any> : MatrixFeature {
+    /**
+     * The determinant of the matrix that owns this feature.
+     */
     public val determinant: T
 }
 
+/**
+ * Produces a [DeterminantFeature] where the [DeterminantFeature.determinant] is [determinant].
+ *
+ * @param determinant the value of determinant.
+ * @return a new [DeterminantFeature].
+ */
 @Suppress("FunctionName")
 public fun <T : Any> DeterminantFeature(determinant: T): DeterminantFeature<T> = object : DeterminantFeature<T> {
     override val determinant: T = determinant
 }
 
 /**
- * Lower triangular matrix
+ * Matrices with this feature are lower triangular ones.
  */
 public object LFeature : MatrixFeature
 
 /**
- * Upper triangular feature
+ * Matrices with this feature are upper triangular ones.
  */
 public object UFeature : MatrixFeature
 
 /**
- * TODO add documentation
+ * Matrices with this feature support LU factorization with partial pivoting: *[p] &middot; a = [l] &middot; [u]* where
+ * *a* is the owning matrix.
+ *
+ * @param T the type of matrices' items.
  */
-public interface LUPDecompositionFeature<T : Any> : MatrixFeature {
+public interface LupDecompositionFeature<T : Any> : MatrixFeature {
+    /**
+     * The lower triangular matrix in this decomposition. It may have [LFeature].
+     */
     public val l: FeaturedMatrix<T>
+
+    /**
+     * The upper triangular matrix in this decomposition. It may have [UFeature].
+     */
     public val u: FeaturedMatrix<T>
+
+    /**
+     * The permutation matrix in this decomposition.
+     */
     public val p: FeaturedMatrix<T>
 }
 
+/**
+ * Matrices with this feature are orthogonal ones: *a &middot; a<sup>T</sup> = u* where *a* is the owning matrix, *u*
+ * is the unit matrix ([UnitFeature]).
+ */
+public object OrthogonalFeature : MatrixFeature
+
+/**
+ * Matrices with this feature support QR factorization: *a = [q] &middot; [r]* where *a* is the owning matrix.
+ *
+ * @param T the type of matrices' items.
+ */
+public interface QRDecompositionFeature<T : Any> : MatrixFeature {
+    /**
+     * The orthogonal matrix in this decomposition. It may have [OrthogonalFeature].
+     */
+    public val q: FeaturedMatrix<T>
+
+    /**
+     * The upper triangular matrix in this decomposition. It may have [UFeature].
+     */
+    public val r: FeaturedMatrix<T>
+}
+
+/**
+ * Matrices with this feature support Cholesky factorization: *a = [l] &middot; [l]<sup>H</sup>* where *a* is the
+ * owning matrix.
+ *
+ * @param T the type of matrices' items.
+ */
+public interface CholeskyDecompositionFeature<T : Any> : MatrixFeature {
+    /**
+     * The triangular matrix in this decomposition. It may have either [UFeature] or [LFeature].
+     */
+    public val l: FeaturedMatrix<T>
+}
+
+/**
+ * Matrices with this feature support SVD: *a = [u] &middot; [s] &middot; [v]<sup>H</sup>* where *a* is the owning
+ * matrix.
+ *
+ * @param T the type of matrices' items.
+ */
+public interface SingularValueDecompositionFeature<T : Any> : MatrixFeature {
+    /**
+     * The matrix in this decomposition. It is unitary, and it consists from left singular vectors.
+     */
+    public val u: FeaturedMatrix<T>
+
+    /**
+     * The matrix in this decomposition. Its main diagonal elements are singular values.
+     */
+    public val s: FeaturedMatrix<T>
+
+    /**
+     * The matrix in this decomposition. It is unitary, and it consists from right singular vectors.
+     */
+    public val v: FeaturedMatrix<T>
+
+    /**
+     * The buffer of singular values of this SVD.
+     */
+    public val singularValues: Point<T>
+}
+
 //TODO add sparse matrix feature

kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Structure2D.kt

@@ -1,12 +1,40 @@
 package kscience.kmath.structures
 
 /**
- * A structure that is guaranteed to be two-dimensional
+ * A structure that is guaranteed to be two-dimensional.
+ *
+ * @param T the type of items.
  */
 public interface Structure2D<T> : NDStructure<T> {
+    /**
+     * The number of rows in this structure.
+     */
     public val rowNum: Int get() = shape[0]
+
+    /**
+     * The number of columns in this structure.
+     */
     public val colNum: Int get() = shape[1]
 
+    /**
+     * The buffer of rows of this structure. It gets elements from the structure dynamically.
+     */
+    public val rows: Buffer<Buffer<T>>
+        get() = VirtualBuffer(rowNum) { i -> VirtualBuffer(colNum) { j -> get(i, j) } }
+
+    /**
+     * The buffer of columns of this structure. It gets elements from the structure dynamically.
+     */
+    public val columns: Buffer<Buffer<T>>
+        get() = VirtualBuffer(colNum) { j -> VirtualBuffer(rowNum) { i -> get(i, j) } }
+
+    /**
+     * Retrieves an element from the structure by two indices.
+     *
+     * @param i the first index.
+     * @param j the second index.
+     * @return an element.
+     */
     public operator fun get(i: Int, j: Int): T
 
     override operator fun get(index: IntArray): T {
@@ -14,15 +42,9 @@ public interface Structure2D<T> : NDStructure<T> {
         return get(index[0], index[1])
     }
 
-    public val rows: Buffer<Buffer<T>>
-        get() = VirtualBuffer(rowNum) { i -> VirtualBuffer(colNum) { j -> get(i, j) } }
-
-    public val columns: Buffer<Buffer<T>>
-        get() = VirtualBuffer(colNum) { j -> VirtualBuffer(rowNum) { i -> get(i, j) } }
-
     override fun elements(): Sequence<Pair<IntArray, T>> = sequence {
-        for (i in (0 until rowNum))
-            for (j in (0 until colNum)) yield(intArrayOf(i, j) to get(i, j))
+        for (i in 0 until rowNum)
+            for (j in 0 until colNum) yield(intArrayOf(i, j) to get(i, j))
     }
 
     public companion object
@@ -47,4 +69,9 @@ public fun <T> NDStructure<T>.as2D(): Structure2D<T> = if (shape.size == 2)
 else
     error("Can't create 2d-structure from ${shape.size}d-structure")
 
+/**
+ * Alias for [Structure2D] with more familiar name.
+ *
+ * @param T the type of items.
+ */
 public typealias Matrix<T> = Structure2D<T>

kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrix.kt

@@ -1,12 +1,10 @@
 package kscience.kmath.ejml
 
+import kscience.kmath.linear.*
+import kscience.kmath.structures.NDStructure
+import kscience.kmath.structures.RealBuffer
 import org.ejml.dense.row.factory.DecompositionFactory_DDRM
 import org.ejml.simple.SimpleMatrix
-import kscience.kmath.linear.DeterminantFeature
-import kscience.kmath.linear.FeaturedMatrix
-import kscience.kmath.linear.LUPDecompositionFeature
-import kscience.kmath.linear.MatrixFeature
-import kscience.kmath.structures.NDStructure
 
 /**
  * Represents featured matrix over EJML [SimpleMatrix].
@@ -14,42 +12,71 @@ import kscience.kmath.structures.NDStructure
  * @property origin the underlying [SimpleMatrix].
  * @author Iaroslav Postovalov
  */
-public class EjmlMatrix(public val origin: SimpleMatrix, features: Set<MatrixFeature>? = null) : FeaturedMatrix<Double> {
+public class EjmlMatrix(public val origin: SimpleMatrix, features: Set<MatrixFeature> = emptySet()) :
+    FeaturedMatrix<Double> {
     public override val rowNum: Int
         get() = origin.numRows()
 
     public override val colNum: Int
         get() = origin.numCols()
 
-    public override val shape: IntArray
-        get() = intArrayOf(origin.numRows(), origin.numCols())
+    public override val shape: IntArray by lazy { intArrayOf(rowNum, colNum) }
 
-    public override val features: Set<MatrixFeature> = setOf(
-        object : LUPDecompositionFeature<Double>, DeterminantFeature<Double> {
-            override val determinant: Double
-                get() = origin.determinant()
+    public override val features: Set<MatrixFeature> = hashSetOf(
+        object : InverseMatrixFeature<Double> {
+            override val inverse: FeaturedMatrix<Double> by lazy { EjmlMatrix(origin.invert()) }
+        },
 
-            private val lup by lazy {
-                val ludecompositionF64 = DecompositionFactory_DDRM.lu(origin.numRows(), origin.numCols())
-                    .also { it.decompose(origin.ddrm.copy()) }
+        object : DeterminantFeature<Double> {
+            override val determinant: Double by lazy(origin::determinant)
+        },
 
-                Triple(
-                    EjmlMatrix(SimpleMatrix(ludecompositionF64.getRowPivot(null))),
-                    EjmlMatrix(SimpleMatrix(ludecompositionF64.getLower(null))),
-                    EjmlMatrix(SimpleMatrix(ludecompositionF64.getUpper(null))),
-                )
+        object : SingularValueDecompositionFeature<Double> {
+            private val svd by lazy {
+                DecompositionFactory_DDRM.svd(origin.numRows(), origin.numCols(), true, true, false)
+                    .apply { decompose(origin.ddrm.copy()) }
             }
 
-            override val l: FeaturedMatrix<Double>
-                get() = lup.second
+            override val u: FeaturedMatrix<Double> by lazy { EjmlMatrix(SimpleMatrix(svd.getU(null, false))) }
+            override val s: FeaturedMatrix<Double> by lazy { EjmlMatrix(SimpleMatrix(svd.getW(null))) }
+            override val v: FeaturedMatrix<Double> by lazy { EjmlMatrix(SimpleMatrix(svd.getV(null, false))) }
+            override val singularValues: Point<Double> by lazy { RealBuffer(svd.singularValues) }
+        },
 
-            override val u: FeaturedMatrix<Double>
-                get() = lup.third
+        object : QRDecompositionFeature<Double> {
+            private val qr by lazy {
+                DecompositionFactory_DDRM.qr().apply { decompose(origin.ddrm.copy()) }
+            }
 
-            override val p: FeaturedMatrix<Double>
-                get() = lup.first
-        }
-    ) union features.orEmpty()
+            override val q: FeaturedMatrix<Double> by lazy { EjmlMatrix(SimpleMatrix(qr.getQ(null, false))) }
+            override val r: FeaturedMatrix<Double> by lazy { EjmlMatrix(SimpleMatrix(qr.getR(null, false))) }
+        },
+
+        object : CholeskyDecompositionFeature<Double> {
+            override val l: FeaturedMatrix<Double> by lazy {
+                val cholesky =
+                    DecompositionFactory_DDRM.chol(rowNum, true).apply { decompose(origin.ddrm.copy()) }
+
+                EjmlMatrix(SimpleMatrix(cholesky.getT(null)), setOf(LFeature))
+            }
+        },
+
+        object : LupDecompositionFeature<Double> {
+            private val lup by lazy {
+                DecompositionFactory_DDRM.lu(origin.numRows(), origin.numCols()).apply { decompose(origin.ddrm.copy()) }
+            }
+
+            override val l: FeaturedMatrix<Double> by lazy {
+                EjmlMatrix(SimpleMatrix(lup.getLower(null)), setOf(LFeature))
+            }
+
+            override val u: FeaturedMatrix<Double> by lazy {
+                EjmlMatrix(SimpleMatrix(lup.getUpper(null)), setOf(UFeature))
+            }
+
+            override val p: FeaturedMatrix<Double> by lazy { EjmlMatrix(SimpleMatrix(lup.getRowPivot(null))) }
+        },
+    ) union features
 
     public override fun suggestFeature(vararg features: MatrixFeature): EjmlMatrix =
         EjmlMatrix(origin, this.features + features)

kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt

@@ -17,7 +17,6 @@ public fun Matrix<Double>.toEjml(): EjmlMatrix =
  * @author Iaroslav Postovalov
  */
 public object EjmlMatrixContext : MatrixContext<Double, EjmlMatrix> {
-
     /**
      * Converts this vector to EJML one.
      */
@@ -33,6 +32,11 @@ public object EjmlMatrixContext : MatrixContext<Double, EjmlMatrix> {
             }
         })
 
+    override fun point(size: Int, initializer: (Int) -> Double): Point<Double> =
+        EjmlVector(SimpleMatrix(size, 1).also {
+            (0 until it.numRows()).forEach { row -> it[row, 0] = initializer(row) }
+        })
+
     public override fun Matrix<Double>.dot(other: Matrix<Double>): EjmlMatrix =
         EjmlMatrix(toEjml().origin.mult(other.toEjml().origin))
 
@@ -73,12 +77,3 @@ public fun EjmlMatrixContext.solve(a: Matrix<Double>, b: Matrix<Double>): EjmlMa
  */
 public fun EjmlMatrixContext.solve(a: Matrix<Double>, b: Point<Double>): EjmlVector =
     EjmlVector(a.toEjml().origin.solve(b.toEjml().origin))
-
-/**
- * Returns the inverse of given matrix: b = a^(-1).
- *
- * @param a the matrix.
- * @return the inverse of this matrix.
- * @author Iaroslav Postovalov
- */
-public fun EjmlMatrixContext.inverse(a: Matrix<Double>): EjmlMatrix = EjmlMatrix(a.toEjml().origin.invert())

kmath-ejml/src/test/kotlin/kscience/kmath/ejml/EjmlMatrixTest.kt

@@ -1,7 +1,7 @@
 package kscience.kmath.ejml
 
 import kscience.kmath.linear.DeterminantFeature
-import kscience.kmath.linear.LUPDecompositionFeature
+import kscience.kmath.linear.LupDecompositionFeature
 import kscience.kmath.linear.MatrixFeature
 import kscience.kmath.linear.getFeature
 import org.ejml.dense.row.factory.DecompositionFactory_DDRM
@@ -44,7 +44,7 @@ internal class EjmlMatrixTest {
         val w = EjmlMatrix(m)
         val det = w.getFeature<DeterminantFeature<Double>>() ?: fail()
         assertEquals(m.determinant(), det.determinant)
-        val lup = w.getFeature<LUPDecompositionFeature<Double>>() ?: fail()
+        val lup = w.getFeature<LupDecompositionFeature<Double>>() ?: fail()
 
         val ludecompositionF64 = DecompositionFactory_DDRM.lu(m.numRows(), m.numCols())
             .also { it.decompose(m.ddrm.copy()) }