From 7fdd001a77e14aa745d629d528971b979a2f019a Mon Sep 17 00:00:00 2001 From: Iaroslav Postovalov Date: Sat, 16 Jan 2021 15:51:36 +0700 Subject: [PATCH] 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 --- .../kscience/kmath/linear/FeaturedMatrix.kt | 10 +- .../kscience/kmath/linear/LUPDecomposition.kt | 14 ++- .../kscience/kmath/linear/MatrixFeatures.kt | 114 ++++++++++++++++-- .../kscience/kmath/structures/Structure2D.kt | 45 +++++-- .../kotlin/kscience/kmath/ejml/EjmlMatrix.kt | 83 ++++++++----- .../kscience/kmath/ejml/EjmlMatrixContext.kt | 15 +-- .../kscience/kmath/ejml/EjmlMatrixTest.kt | 4 +- 7 files changed, 217 insertions(+), 68 deletions(-) diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt index 68272203c..119f5d844 100644 --- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt +++ b/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 : Matrix { - 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 /** 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 bf2a9f59e..75acaf8a1 100644 --- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt +++ b/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( public val context: MatrixContext>, public val elementContext: Field, - public val lu: Structure2D, + public val lu: Matrix, public val pivot: IntArray, private val even: Boolean, -) : LUPDecompositionFeature, DeterminantFeature { - +) : LupDecompositionFeature, DeterminantFeature { /** * Returns the matrix L of the decomposition. * @@ -151,7 +150,10 @@ public inline fun , F : Field> GenericMatrixContext public fun MatrixContext>.lup(matrix: Matrix): LUPDecomposition = lup(Buffer.Companion::real, RealField, matrix) { it < 1e-11 } -public fun LUPDecomposition.solveWithLUP(factory: MutableBufferFactory, matrix: Matrix): FeaturedMatrix { +public fun LUPDecomposition.solveWithLUP( + factory: MutableBufferFactory, + matrix: Matrix +): FeaturedMatrix { require(matrix.rowNum == pivot.size) { "Matrix dimension mismatch. Expected ${pivot.size}, but got ${matrix.colNum}" } BufferAccessor2D(matrix.rowNum, matrix.colNum, factory).run { @@ -217,4 +219,4 @@ public inline fun , F : Field> GenericMatrixContext matrix: Matrix, noinline bufferFactory: MutableBufferFactory = MutableBuffer.Companion::auto, noinline checkSingular: (T) -> Boolean, -): FeaturedMatrix = solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum), bufferFactory, checkSingular) \ No newline at end of file +): FeaturedMatrix = solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum), bufferFactory, checkSingular) diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixFeatures.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixFeatures.kt index a82032e50..767b58eba 100644 --- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixFeatures.kt +++ b/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`-1 where `a` is the owning matrix. + * + * @param T the type of matrices' items. */ public interface InverseMatrixFeature : MatrixFeature { + /** + * The inverse matrix of the matrix that owns this feature. + */ public val inverse: FeaturedMatrix } /** - * A determinant container + * Matrices with this feature can compute their determinant. */ public interface DeterminantFeature : 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 DeterminantFeature(determinant: T): DeterminantFeature = object : DeterminantFeature { 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] · a = [l] · [u]* where + * *a* is the owning matrix. + * + * @param T the type of matrices' items. */ -public interface LUPDecompositionFeature : MatrixFeature { +public interface LupDecompositionFeature : MatrixFeature { + /** + * The lower triangular matrix in this decomposition. It may have [LFeature]. + */ public val l: FeaturedMatrix + + /** + * The upper triangular matrix in this decomposition. It may have [UFeature]. + */ public val u: FeaturedMatrix + + /** + * The permutation matrix in this decomposition. + */ public val p: FeaturedMatrix } +/** + * Matrices with this feature are orthogonal ones: *a · aT = 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] · [r]* where *a* is the owning matrix. + * + * @param T the type of matrices' items. + */ +public interface QRDecompositionFeature : MatrixFeature { + /** + * The orthogonal matrix in this decomposition. It may have [OrthogonalFeature]. + */ + public val q: FeaturedMatrix + + /** + * The upper triangular matrix in this decomposition. It may have [UFeature]. + */ + public val r: FeaturedMatrix +} + +/** + * Matrices with this feature support Cholesky factorization: *a = [l] · [l]H* where *a* is the + * owning matrix. + * + * @param T the type of matrices' items. + */ +public interface CholeskyDecompositionFeature : MatrixFeature { + /** + * The triangular matrix in this decomposition. It may have either [UFeature] or [LFeature]. + */ + public val l: FeaturedMatrix +} + +/** + * Matrices with this feature support SVD: *a = [u] · [s] · [v]H* where *a* is the owning + * matrix. + * + * @param T the type of matrices' items. + */ +public interface SingularValueDecompositionFeature : MatrixFeature { + /** + * The matrix in this decomposition. It is unitary, and it consists from left singular vectors. + */ + public val u: FeaturedMatrix + + /** + * The matrix in this decomposition. Its main diagonal elements are singular values. + */ + public val s: FeaturedMatrix + + /** + * The matrix in this decomposition. It is unitary, and it consists from right singular vectors. + */ + public val v: FeaturedMatrix + + /** + * The buffer of singular values of this SVD. + */ + public val singularValues: Point +} + //TODO add sparse matrix feature diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Structure2D.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Structure2D.kt index 25fdf3f3d..bac7d3389 100644 --- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Structure2D.kt +++ b/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 : NDStructure { + /** + * 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> + 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> + 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 : NDStructure { return get(index[0], index[1]) } - public val rows: Buffer> - get() = VirtualBuffer(rowNum) { i -> VirtualBuffer(colNum) { j -> get(i, j) } } - - public val columns: Buffer> - get() = VirtualBuffer(colNum) { j -> VirtualBuffer(rowNum) { i -> get(i, j) } } - override fun elements(): Sequence> = 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 NDStructure.as2D(): Structure2D = 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 = Structure2D diff --git a/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrix.kt b/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrix.kt index ed6b1571e..5b7d0a01b 100644 --- a/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrix.kt +++ b/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? = null) : FeaturedMatrix { +public class EjmlMatrix(public val origin: SimpleMatrix, features: Set = emptySet()) : + FeaturedMatrix { 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 = setOf( - object : LUPDecompositionFeature, DeterminantFeature { - override val determinant: Double - get() = origin.determinant() + public override val features: Set = hashSetOf( + object : InverseMatrixFeature { + override val inverse: FeaturedMatrix 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 { + 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 { + private val svd by lazy { + DecompositionFactory_DDRM.svd(origin.numRows(), origin.numCols(), true, true, false) + .apply { decompose(origin.ddrm.copy()) } } - override val l: FeaturedMatrix - get() = lup.second + override val u: FeaturedMatrix by lazy { EjmlMatrix(SimpleMatrix(svd.getU(null, false))) } + override val s: FeaturedMatrix by lazy { EjmlMatrix(SimpleMatrix(svd.getW(null))) } + override val v: FeaturedMatrix by lazy { EjmlMatrix(SimpleMatrix(svd.getV(null, false))) } + override val singularValues: Point by lazy { RealBuffer(svd.singularValues) } + }, - override val u: FeaturedMatrix - get() = lup.third + object : QRDecompositionFeature { + private val qr by lazy { + DecompositionFactory_DDRM.qr().apply { decompose(origin.ddrm.copy()) } + } - override val p: FeaturedMatrix - get() = lup.first - } - ) union features.orEmpty() + override val q: FeaturedMatrix by lazy { EjmlMatrix(SimpleMatrix(qr.getQ(null, false))) } + override val r: FeaturedMatrix by lazy { EjmlMatrix(SimpleMatrix(qr.getR(null, false))) } + }, + + object : CholeskyDecompositionFeature { + override val l: FeaturedMatrix by lazy { + val cholesky = + DecompositionFactory_DDRM.chol(rowNum, true).apply { decompose(origin.ddrm.copy()) } + + EjmlMatrix(SimpleMatrix(cholesky.getT(null)), setOf(LFeature)) + } + }, + + object : LupDecompositionFeature { + private val lup by lazy { + DecompositionFactory_DDRM.lu(origin.numRows(), origin.numCols()).apply { decompose(origin.ddrm.copy()) } + } + + override val l: FeaturedMatrix by lazy { + EjmlMatrix(SimpleMatrix(lup.getLower(null)), setOf(LFeature)) + } + + override val u: FeaturedMatrix by lazy { + EjmlMatrix(SimpleMatrix(lup.getUpper(null)), setOf(UFeature)) + } + + override val p: FeaturedMatrix by lazy { EjmlMatrix(SimpleMatrix(lup.getRowPivot(null))) } + }, + ) union features public override fun suggestFeature(vararg features: MatrixFeature): EjmlMatrix = EjmlMatrix(origin, this.features + features) diff --git a/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt b/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt index 31792e39c..f8791b72c 100644 --- a/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt +++ b/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt @@ -17,7 +17,6 @@ public fun Matrix.toEjml(): EjmlMatrix = * @author Iaroslav Postovalov */ public object EjmlMatrixContext : MatrixContext { - /** * Converts this vector to EJML one. */ @@ -33,6 +32,11 @@ public object EjmlMatrixContext : MatrixContext { } }) + override fun point(size: Int, initializer: (Int) -> Double): Point = + EjmlVector(SimpleMatrix(size, 1).also { + (0 until it.numRows()).forEach { row -> it[row, 0] = initializer(row) } + }) + public override fun Matrix.dot(other: Matrix): EjmlMatrix = EjmlMatrix(toEjml().origin.mult(other.toEjml().origin)) @@ -73,12 +77,3 @@ public fun EjmlMatrixContext.solve(a: Matrix, b: Matrix): EjmlMa */ public fun EjmlMatrixContext.solve(a: Matrix, b: Point): 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): EjmlMatrix = EjmlMatrix(a.toEjml().origin.invert()) diff --git a/kmath-ejml/src/test/kotlin/kscience/kmath/ejml/EjmlMatrixTest.kt b/kmath-ejml/src/test/kotlin/kscience/kmath/ejml/EjmlMatrixTest.kt index e0f15be83..70b82a3cb 100644 --- a/kmath-ejml/src/test/kotlin/kscience/kmath/ejml/EjmlMatrixTest.kt +++ b/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>() ?: fail() assertEquals(m.determinant(), det.determinant) - val lup = w.getFeature>() ?: fail() + val lup = w.getFeature>() ?: fail() val ludecompositionF64 = DecompositionFactory_DDRM.lu(m.numRows(), m.numCols()) .also { it.decompose(m.ddrm.copy()) }