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

This commit is contained in:
Iaroslav Postovalov 2021-01-16 15:51:36 +07:00
parent fed17a17a9
commit 7fdd001a77
No known key found for this signature in database
GPG Key ID: 46E15E4A31B3BCD7
7 changed files with 217 additions and 68 deletions

View File

@ -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>
/**

View File

@ -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 {

View File

@ -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 {
/**
* The lower triangular matrix in this decomposition. It may have [LFeature].
*/
public interface LUPDecompositionFeature<T : Any> : MatrixFeature {
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

View File

@ -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>

View File

@ -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()) }
},
object : DeterminantFeature<Double> {
override val determinant: Double by lazy(origin::determinant)
},
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 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) }
},
object : QRDecompositionFeature<Double> {
private val qr by lazy {
DecompositionFactory_DDRM.qr().apply { decompose(origin.ddrm.copy()) }
}
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 {
val ludecompositionF64 = DecompositionFactory_DDRM.lu(origin.numRows(), origin.numCols())
.also { it.decompose(origin.ddrm.copy()) }
Triple(
EjmlMatrix(SimpleMatrix(ludecompositionF64.getRowPivot(null))),
EjmlMatrix(SimpleMatrix(ludecompositionF64.getLower(null))),
EjmlMatrix(SimpleMatrix(ludecompositionF64.getUpper(null))),
)
DecompositionFactory_DDRM.lu(origin.numRows(), origin.numCols()).apply { decompose(origin.ddrm.copy()) }
}
override val l: FeaturedMatrix<Double>
get() = lup.second
override val u: FeaturedMatrix<Double>
get() = lup.third
override val p: FeaturedMatrix<Double>
get() = lup.first
override val l: FeaturedMatrix<Double> by lazy {
EjmlMatrix(SimpleMatrix(lup.getLower(null)), setOf(LFeature))
}
) union features.orEmpty()
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)

View File

@ -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())

View File

@ -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()) }