Merge remote-tracking branch 'origin/dev' into refactor/ndalgebra

# Conflicts:
#	CHANGELOG.md
#	examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt
#	kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt
#	kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt
This commit is contained in:
Alexander Nozik 2021-01-28 20:07:43 +03:00
commit 45866b500f
9 changed files with 79 additions and 33 deletions

View File

@ -32,9 +32,10 @@
- Use `Point<Double>` instead of specialized type in `kmath-for-real` - Use `Point<Double>` instead of specialized type in `kmath-for-real`
- Optimized dot product for buffer matrices moved to `kmath-for-real` - Optimized dot product for buffer matrices moved to `kmath-for-real`
- EjmlMatrix context is an object - EjmlMatrix context is an object
- Matrix LUP `inverse` renamed to `inverseWithLUP` - Matrix LUP `inverse` renamed to `inverseWithLup`
- `NumericAlgebra` moved outside of regular algebra chain (`Ring` no longer implements it). - `NumericAlgebra` moved outside of regular algebra chain (`Ring` no longer implements it).
- Features moved to NDStructure and became transparent. - Features moved to NDStructure and became transparent.
- Capitalization of LUP in many names changed to Lup.
- Refactored `NDStructure` algebra to be more simple, preferring under-the-hood conversion to explicit NDStructure types - Refactored `NDStructure` algebra to be more simple, preferring under-the-hood conversion to explicit NDStructure types
### Deprecated ### Deprecated

View File

@ -8,6 +8,7 @@ import kscience.kmath.commons.linear.inverse
import kscience.kmath.ejml.EjmlMatrixContext import kscience.kmath.ejml.EjmlMatrixContext
import kscience.kmath.ejml.inverse import kscience.kmath.ejml.inverse
import kscience.kmath.operations.invoke import kscience.kmath.operations.invoke
import kscience.kmath.structures.Matrix
import org.openjdk.jmh.annotations.Scope import org.openjdk.jmh.annotations.Scope
import org.openjdk.jmh.annotations.State import org.openjdk.jmh.annotations.State
import kotlin.random.Random import kotlin.random.Random
@ -25,10 +26,8 @@ class LinearAlgebraBenchmark {
} }
@Benchmark @Benchmark
fun kmathLUPInversion() { fun kmathLupInversion() {
MatrixContext.real{ MatrixContext.real.inverseWithLup(matrix)
inverseWithLUP(matrix)
}
} }
@Benchmark @Benchmark

View File

@ -2,6 +2,8 @@ package kscience.kmath.commons.linear
import kscience.kmath.linear.* import kscience.kmath.linear.*
import kscience.kmath.misc.UnstableKMathAPI import kscience.kmath.misc.UnstableKMathAPI
import kscience.kmath.structures.Matrix
import kscience.kmath.structures.RealBuffer
import org.apache.commons.math3.linear.* import org.apache.commons.math3.linear.*
import kotlin.reflect.KClass import kotlin.reflect.KClass
import kotlin.reflect.cast import kotlin.reflect.cast
@ -13,8 +15,40 @@ public inline class CMMatrix(public val origin: RealMatrix) : Matrix<Double> {
@UnstableKMathAPI @UnstableKMathAPI
override fun <T : Any> getFeature(type: KClass<T>): T? = when (type) { override fun <T : Any> getFeature(type: KClass<T>): T? = when (type) {
DiagonalFeature::class -> if (origin is DiagonalMatrix) DiagonalFeature else null DiagonalFeature::class -> if (origin is DiagonalMatrix) DiagonalFeature else null
DeterminantFeature::class, LupDecompositionFeature::class -> object :
DeterminantFeature<Double>,
LupDecompositionFeature<Double> {
private val lup by lazy { LUDecomposition(origin) }
override val determinant: Double by lazy { lup.determinant }
override val l: Matrix<Double> by lazy { CMMatrix(lup.l) + LFeature }
override val u: Matrix<Double> by lazy { CMMatrix(lup.u) + UFeature }
override val p: Matrix<Double> by lazy { CMMatrix(lup.p) }
}
CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Double> {
override val l: Matrix<Double> by lazy {
val cholesky = CholeskyDecomposition(origin)
CMMatrix(cholesky.l) + LFeature
}
}
QRDecompositionFeature::class -> object : QRDecompositionFeature<Double> {
private val qr by lazy { QRDecomposition(origin) }
override val q: Matrix<Double> by lazy { CMMatrix(qr.q) + OrthogonalFeature }
override val r: Matrix<Double> by lazy { CMMatrix(qr.r) + UFeature }
}
SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature<Double> {
private val sv by lazy { SingularValueDecomposition(origin) }
override val u: Matrix<Double> by lazy { CMMatrix(sv.u) }
override val s: Matrix<Double> by lazy { CMMatrix(sv.s) }
override val v: Matrix<Double> by lazy { CMMatrix(sv.v) }
override val singularValues: Point<Double> by lazy { RealBuffer(sv.singularValues) }
}
else -> null else -> null
}?.let { type.cast(it) } }?.let(type::cast)
public override operator fun get(i: Int, j: Int): Double = origin.getEntry(i, j) public override operator fun get(i: Int, j: Int): Double = origin.getEntry(i, j)
} }

View File

@ -152,7 +152,7 @@ public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext
public fun MatrixContext<Double, Matrix<Double>>.lup(matrix: Matrix<Double>): LupDecomposition<Double> = public fun MatrixContext<Double, Matrix<Double>>.lup(matrix: Matrix<Double>): LupDecomposition<Double> =
lup(Buffer.Companion::real, RealField, matrix) { it < 1e-11 } lup(Buffer.Companion::real, RealField, matrix) { it < 1e-11 }
public fun <T : Any> LupDecomposition<T>.solveWithLUP( public fun <T : Any> LupDecomposition<T>.solveWithLup(
factory: MutableBufferFactory<T>, factory: MutableBufferFactory<T>,
matrix: Matrix<T>, matrix: Matrix<T>,
): Matrix<T> { ): Matrix<T> {
@ -200,14 +200,14 @@ public fun <T : Any> LupDecomposition<T>.solveWithLUP(
} }
} }
public inline fun <reified T : Any> LupDecomposition<T>.solveWithLUP(matrix: Matrix<T>): Matrix<T> = public inline fun <reified T : Any> LupDecomposition<T>.solveWithLup(matrix: Matrix<T>): Matrix<T> =
solveWithLUP(MutableBuffer.Companion::auto, matrix) solveWithLup(MutableBuffer.Companion::auto, matrix)
/** /**
* Solve a linear equation **a*x = b** using LUP decomposition * Solves a system of linear equations *ax = b** using LUP decomposition.
*/ */
@OptIn(UnstableKMathAPI::class) @OptIn(UnstableKMathAPI::class)
public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, Matrix<T>>.solveWithLUP( public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, Matrix<T>>.solveWithLup(
a: Matrix<T>, a: Matrix<T>,
b: Matrix<T>, b: Matrix<T>,
noinline bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::auto, noinline bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::auto,
@ -215,26 +215,26 @@ public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext
): Matrix<T> { ): Matrix<T> {
// Use existing decomposition if it is provided by matrix // Use existing decomposition if it is provided by matrix
val decomposition = a.getFeature() ?: lup(bufferFactory, elementContext, a, checkSingular) val decomposition = a.getFeature() ?: lup(bufferFactory, elementContext, a, checkSingular)
return decomposition.solveWithLUP(bufferFactory, b) return decomposition.solveWithLup(bufferFactory, b)
} }
public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, Matrix<T>>.inverseWithLUP( public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, Matrix<T>>.inverseWithLup(
matrix: Matrix<T>, matrix: Matrix<T>,
noinline bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::auto, noinline bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::auto,
noinline checkSingular: (T) -> Boolean, noinline checkSingular: (T) -> Boolean,
): Matrix<T> = solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum), bufferFactory, checkSingular) ): Matrix<T> = solveWithLup(matrix, one(matrix.rowNum, matrix.colNum), bufferFactory, checkSingular)
@OptIn(UnstableKMathAPI::class) @OptIn(UnstableKMathAPI::class)
public fun RealMatrixContext.solveWithLUP(a: Matrix<Double>, b: Matrix<Double>): Matrix<Double> { public fun RealMatrixContext.solveWithLup(a: Matrix<Double>, b: Matrix<Double>): Matrix<Double> {
// Use existing decomposition if it is provided by matrix // Use existing decomposition if it is provided by matrix
val bufferFactory: MutableBufferFactory<Double> = MutableBuffer.Companion::real val bufferFactory: MutableBufferFactory<Double> = MutableBuffer.Companion::real
val decomposition: LupDecomposition<Double> = a.getFeature() ?: lup(bufferFactory, RealField, a) { it < 1e-11 } val decomposition: LupDecomposition<Double> = a.getFeature() ?: lup(bufferFactory, RealField, a) { it < 1e-11 }
return decomposition.solveWithLUP(bufferFactory, b) return decomposition.solveWithLup(bufferFactory, b)
} }
/** /**
* Inverses a square matrix using LUP decomposition. Non square matrix will throw a error. * Inverses a square matrix using LUP decomposition. Non square matrix will throw a error.
*/ */
public fun RealMatrixContext.inverseWithLUP(matrix: Matrix<Double>): Matrix<Double> = public fun RealMatrixContext.inverseWithLup(matrix: Matrix<Double>): Matrix<Double> =
solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum)) solveWithLup(matrix, one(matrix.rowNum, matrix.colNum))

View File

@ -37,6 +37,8 @@ public interface InverseMatrixFeature<T : Any> : MatrixFeature {
/** /**
* Matrices with this feature can compute their determinant. * Matrices with this feature can compute their determinant.
*
* @param T the type of matrices' items.
*/ */
public interface DeterminantFeature<T : Any> : MatrixFeature { public interface DeterminantFeature<T : Any> : MatrixFeature {
/** /**

View File

@ -8,7 +8,7 @@ class RealLUSolverTest {
@Test @Test
fun testInvertOne() { fun testInvertOne() {
val matrix = MatrixContext.real.one(2, 2) val matrix = MatrixContext.real.one(2, 2)
val inverted = MatrixContext.real.inverseWithLUP(matrix) val inverted = MatrixContext.real.inverseWithLup(matrix)
assertEquals(matrix, inverted) assertEquals(matrix, inverted)
} }
@ -36,7 +36,7 @@ class RealLUSolverTest {
1.0, 3.0 1.0, 3.0
) )
val inverted = MatrixContext.real.inverseWithLUP(matrix) val inverted = MatrixContext.real.inverseWithLup(matrix)
val expected = Matrix.square( val expected = Matrix.square(
0.375, -0.125, 0.375, -0.125,

View File

@ -15,21 +15,20 @@ import kotlin.reflect.cast
* @property origin the underlying [SimpleMatrix]. * @property origin the underlying [SimpleMatrix].
* @author Iaroslav Postovalov * @author Iaroslav Postovalov
*/ */
public class EjmlMatrix( public class EjmlMatrix(public val origin: SimpleMatrix) : Matrix<Double> {
public val origin: SimpleMatrix,
) : Matrix<Double> {
public override val rowNum: Int get() = origin.numRows() public override val rowNum: Int get() = origin.numRows()
public override val colNum: Int get() = origin.numCols() public override val colNum: Int get() = origin.numCols()
@UnstableKMathAPI @UnstableKMathAPI
override fun <T : Any> getFeature(type: KClass<T>): T? = when (type) { public override fun <T : Any> getFeature(type: KClass<T>): T? = when (type) {
InverseMatrixFeature::class -> object : InverseMatrixFeature<Double> { InverseMatrixFeature::class -> object : InverseMatrixFeature<Double> {
override val inverse: Matrix<Double> by lazy { EjmlMatrix(origin.invert()) } override val inverse: Matrix<Double> by lazy { EjmlMatrix(origin.invert()) }
} }
DeterminantFeature::class -> object : DeterminantFeature<Double> { DeterminantFeature::class -> object : DeterminantFeature<Double> {
override val determinant: Double by lazy(origin::determinant) override val determinant: Double by lazy(origin::determinant)
} }
SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature<Double> { SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature<Double> {
private val svd by lazy { private val svd by lazy {
DecompositionFactory_DDRM.svd(origin.numRows(), origin.numCols(), true, true, false) DecompositionFactory_DDRM.svd(origin.numRows(), origin.numCols(), true, true, false)
@ -41,14 +40,19 @@ public class EjmlMatrix(
override val v: Matrix<Double> by lazy { EjmlMatrix(SimpleMatrix(svd.getV(null, false))) } override val v: Matrix<Double> by lazy { EjmlMatrix(SimpleMatrix(svd.getV(null, false))) }
override val singularValues: Point<Double> by lazy { RealBuffer(svd.singularValues) } override val singularValues: Point<Double> by lazy { RealBuffer(svd.singularValues) }
} }
QRDecompositionFeature::class -> object : QRDecompositionFeature<Double> { QRDecompositionFeature::class -> object : QRDecompositionFeature<Double> {
private val qr by lazy { private val qr by lazy {
DecompositionFactory_DDRM.qr().apply { decompose(origin.ddrm.copy()) } DecompositionFactory_DDRM.qr().apply { decompose(origin.ddrm.copy()) }
} }
override val q: Matrix<Double> by lazy { EjmlMatrix(SimpleMatrix(qr.getQ(null, false))) } override val q: Matrix<Double> by lazy {
override val r: Matrix<Double> by lazy { EjmlMatrix(SimpleMatrix(qr.getR(null, false))) } EjmlMatrix(SimpleMatrix(qr.getQ(null, false))) + OrthogonalFeature
} }
override val r: Matrix<Double> by lazy { EjmlMatrix(SimpleMatrix(qr.getR(null, false))) + UFeature }
}
CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Double> { CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Double> {
override val l: Matrix<Double> by lazy { override val l: Matrix<Double> by lazy {
val cholesky = val cholesky =
@ -57,6 +61,7 @@ public class EjmlMatrix(
EjmlMatrix(SimpleMatrix(cholesky.getT(null))) + LFeature EjmlMatrix(SimpleMatrix(cholesky.getT(null))) + LFeature
} }
} }
LupDecompositionFeature::class -> object : LupDecompositionFeature<Double> { LupDecompositionFeature::class -> object : LupDecompositionFeature<Double> {
private val lup by lazy { private val lup by lazy {
DecompositionFactory_DDRM.lu(origin.numRows(), origin.numCols()).apply { decompose(origin.ddrm.copy()) } DecompositionFactory_DDRM.lu(origin.numRows(), origin.numCols()).apply { decompose(origin.ddrm.copy()) }
@ -72,8 +77,9 @@ public class EjmlMatrix(
override val p: Matrix<Double> by lazy { EjmlMatrix(SimpleMatrix(lup.getRowPivot(null))) } override val p: Matrix<Double> by lazy { EjmlMatrix(SimpleMatrix(lup.getRowPivot(null))) }
} }
else -> null else -> null
}?.let { type.cast(it) } }?.let(type::cast)
public override operator fun get(i: Int, j: Int): Double = origin[i, j] public override operator fun get(i: Int, j: Int): Double = origin[i, j]

View File

@ -1,8 +1,12 @@
package kscience.kmath.real package kscience.kmath.real
import kscience.kmath.linear.* import kscience.kmath.linear.MatrixContext
import kscience.kmath.linear.VirtualMatrix
import kscience.kmath.linear.inverseWithLup
import kscience.kmath.linear.real
import kscience.kmath.misc.UnstableKMathAPI import kscience.kmath.misc.UnstableKMathAPI
import kscience.kmath.structures.Buffer import kscience.kmath.structures.Buffer
import kscience.kmath.structures.Matrix
import kscience.kmath.structures.RealBuffer import kscience.kmath.structures.RealBuffer
import kscience.kmath.structures.asIterable import kscience.kmath.structures.asIterable
import kotlin.math.pow import kotlin.math.pow
@ -148,7 +152,7 @@ public inline fun RealMatrix.map(crossinline transform: (Double) -> Double): Rea
/** /**
* Inverse a square real matrix using LUP decomposition * Inverse a square real matrix using LUP decomposition
*/ */
public fun RealMatrix.inverseWithLUP(): RealMatrix = MatrixContext.real.inverseWithLUP(this) public fun RealMatrix.inverseWithLup(): RealMatrix = MatrixContext.real.inverseWithLup(this)
//extended operations //extended operations