Skip to content

Commit

Permalink
Add sparse matrix builder
Browse files Browse the repository at this point in the history
  • Loading branch information
altavir committed Jan 26, 2025
1 parent 676cf5f commit c70a7c5
Show file tree
Hide file tree
Showing 22 changed files with 316 additions and 94 deletions.
1 change: 1 addition & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@

### Added
- Convenient matrix builders for rows, columns, vstacks and hstacks
- Spreadsheet matrix builder

### Changed

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -9,8 +9,8 @@ import kotlinx.benchmark.Benchmark
import kotlinx.benchmark.Blackhole
import kotlinx.benchmark.Scope
import kotlinx.benchmark.State
import space.kscience.kmath.linear.MatrixBuilder
import space.kscience.kmath.linear.linearSpace
import space.kscience.kmath.linear.matrix
import space.kscience.kmath.linear.symmetric
import space.kscience.kmath.operations.Float64Field
import space.kscience.kmath.tensors.core.symEigJacobi
Expand All @@ -24,7 +24,7 @@ internal class TensorAlgebraBenchmark {
private val random = Random(12224)
private const val dim = 30

private val matrix = Float64Field.linearSpace.matrix(dim, dim).symmetric { _, _ -> random.nextDouble() }
private val matrix = Float64Field.linearSpace.MatrixBuilder(dim, dim).symmetric { _, _ -> random.nextDouble() }
}

@Benchmark
Expand Down
2 changes: 1 addition & 1 deletion build.gradle.kts
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ allprojects {
}

group = "space.kscience"
version = "0.4.1"
version = "0.4.2-dev"
}

dependencies{
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -6,22 +6,21 @@
package space.kscience.kmath.operations

import space.kscience.kmath.commons.linear.CMLinearSpace
import space.kscience.kmath.linear.matrix
import space.kscience.kmath.linear.MatrixBuilder
import space.kscience.kmath.linear.fill
import space.kscience.kmath.nd.Float64BufferND
import space.kscience.kmath.nd.Structure2D
import space.kscience.kmath.nd.mutableStructureND
import space.kscience.kmath.nd.ndAlgebra
import space.kscience.kmath.structures.Float64
import space.kscience.kmath.viktor.viktorAlgebra
import kotlin.collections.component1
import kotlin.collections.component2

fun main() {
val viktorStructure = Float64Field.viktorAlgebra.mutableStructureND(2, 2) { (i, j) ->
if (i == j) 2.0 else 0.0
}

val cmMatrix: Structure2D<Float64> = CMLinearSpace.matrix(2, 2)(0.0, 1.0, 0.0, 3.0)
val cmMatrix: Structure2D<Float64> = CMLinearSpace.MatrixBuilder(2, 2).fill(0.0, 1.0, 0.0, 3.0)

val res: Float64BufferND = Float64Field.ndAlgebra {
exp(viktorStructure) + 2.0 * cmMatrix
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -5,70 +5,59 @@

package space.kscience.kmath.linear

import space.kscience.attributes.FlagAttribute
import space.kscience.attributes.Attributes
import space.kscience.attributes.SafeType
import space.kscience.attributes.WithType
import space.kscience.kmath.UnstableKMathAPI
import space.kscience.kmath.operations.Ring
import space.kscience.kmath.structures.BufferAccessor2D
import space.kscience.kmath.structures.MutableBufferFactory

public class MatrixBuilder<T : Any, out A : Ring<T>>(
/**
* A builder for matrix with fixed size
*/
@UnstableKMathAPI
public class MatrixBuilder<T, out A : Ring<T>>(
public val linearSpace: LinearSpace<T, A>,
public val rows: Int,
public val columns: Int,
public val rowNum: Int,
public val colNum: Int,
) : WithType<T> {

override val type: SafeType<T> get() = linearSpace.type

public operator fun invoke(vararg elements: T): Matrix<T> {
require(rows * columns == elements.size) { "The number of elements ${elements.size} is not equal $rows * $columns" }
return linearSpace.buildMatrix(rows, columns) { i, j -> elements[i * columns + j] }
}
}

//TODO add specific matrix builder functions like diagonal, etc
@UnstableKMathAPI
public fun <T, A : Ring<T>> MatrixBuilder<T, A>.sparse(): SparseMatrix<T> =
SparseMatrix(rowNum, colNum, linearSpace.elementAlgebra.zero)

@UnstableKMathAPI
public fun <T, A : Ring<T>> MatrixBuilder<T, A>.fill(vararg elements: T): Matrix<T> {
require(rowNum * colNum == elements.size) { "The number of elements ${elements.size} is not equal $rowNum * $colNum" }
return linearSpace.buildMatrix(rowNum, colNum) { i, j -> elements[i * colNum + j] }
}

/**
* Create a matrix builder with given number of rows and columns
*/
@UnstableKMathAPI
public fun <T : Any, A : Ring<T>> LinearSpace<T, A>.matrix(rows: Int, columns: Int): MatrixBuilder<T, A> =
public fun <T, A : Ring<T>> LinearSpace<T, A>.MatrixBuilder(rows: Int, columns: Int): MatrixBuilder<T, A> =
MatrixBuilder(this, rows, columns)

@UnstableKMathAPI
public fun <T : Any> LinearSpace<T, Ring<T>>.vector(vararg elements: T): Point<T> {
return buildVector(elements.size) { elements[it] }
}

public inline fun <T : Any> LinearSpace<T, Ring<T>>.row(
size: Int,
crossinline builder: (Int) -> T,
): Matrix<T> = buildMatrix(1, size) { _, j -> builder(j) }

public fun <T : Any> LinearSpace<T, Ring<T>>.row(vararg values: T): Matrix<T> = row(values.size, values::get)

public inline fun <T : Any> LinearSpace<T, Ring<T>>.column(
size: Int,
crossinline builder: (Int) -> T,
): Matrix<T> = buildMatrix(size, 1) { i, _ -> builder(i) }

public fun <T : Any> LinearSpace<T, Ring<T>>.column(vararg values: T): Matrix<T> = column(values.size, values::get)

public object Symmetric : MatrixAttribute<Unit>, FlagAttribute

/**
* Naive implementation of a symmetric matrix builder, that adds a [Symmetric] tag. The resulting matrix contains
* full `size^2` number of elements, but caches elements during calls to save [builder] calls. [builder] is always called in the
* upper triangle region meaning that `i <= j`
* Naive implementation of a symmetric matrix builder, that adds a [Symmetric] tag.
* The resulting matrix contains full `size^2` number of elements,
* but caches elements during calls to save [builder] calls.
* Always called in the upper triangle region meaning that `i <= j`
*/
public fun <T : Any, A : Ring<T>> MatrixBuilder<T, A>.symmetric(
builder: (i: Int, j: Int) -> T,
@UnstableKMathAPI
public fun <T, A : Ring<T>> MatrixBuilder<T, A>.symmetric(
builder: A.(i: Int, j: Int) -> T,
): Matrix<T> {
require(columns == rows) { "In order to build symmetric matrix, number of rows $rows should be equal to number of columns $columns" }
return with(BufferAccessor2D<T?>(rows, rows, MutableBufferFactory(type))) {
require(colNum == rowNum) { "In order to build symmetric matrix, number of rows $rowNum should be equal to number of columns $colNum" }
return with(BufferAccessor2D<T?>(rowNum, rowNum, MutableBufferFactory(type))) {
val cache = HashMap<IntArray, T>()
linearSpace.buildMatrix(rows, rows) { i, j ->
linearSpace.buildMatrix(this@symmetric.rowNum, this@symmetric.rowNum) { i, j ->
val index = intArrayOf(i, j)
val cached = cache[index]
if (cached == null) {
Expand All @@ -81,4 +70,50 @@ public fun <T : Any, A : Ring<T>> MatrixBuilder<T, A>.symmetric(
}
}.withAttribute(Symmetric)
}
}
}

/**
* Create a diagonal matrix with given factory.
*/
@UnstableKMathAPI
public fun <T, A : Ring<T>> MatrixBuilder<T, A>.diagonal(
builder: A.(Int) -> T
): Matrix<T> = with(linearSpace.elementAlgebra) {
require(colNum == rowNum) { "In order to build symmetric matrix, number of rows $rowNum should be equal to number of columns $colNum" }
return VirtualMatrix(rowNum, colNum, attributes = Attributes(IsDiagonal)) { i, j ->
check(i in 0 until rowNum) { "$i out of bounds: 0..<$rowNum" }
check(j in 0 until colNum) { "$j out of bounds: 0..<$colNum" }
if (i == j) {
builder(i)
} else {
zero
}
}
}

/**
* Create a diagonal matrix from elements
*/
@UnstableKMathAPI
public fun <T> MatrixBuilder<T, Ring<T>>.diagonal(vararg elements: T): Matrix<T> {
require(colNum == rowNum) { "In order to build symmetric matrix, number of rows $rowNum should be equal to number of columns $colNum" }
return return VirtualMatrix(rowNum, colNum, attributes = Attributes(IsDiagonal)) { i, j ->
check(i in 0 until rowNum) { "$i out of bounds: 0..<$rowNum" }
check(j in 0 until colNum) { "$j out of bounds: 0..<$colNum" }
if (i == j) {
elements[i]
} else {
linearSpace.elementAlgebra.zero
}
}
}

/**
* Create a lazily evaluated virtual matrix with a given size
*/
@UnstableKMathAPI
public fun <T : Any> MatrixBuilder<T, *>.virtual(
attributes: Attributes = Attributes.EMPTY,
generator: (i: Int, j: Int) -> T,
): VirtualMatrix<T> = VirtualMatrix(rowNum, colNum, attributes, generator)

Original file line number Diff line number Diff line change
@@ -0,0 +1,69 @@
/*
* Copyright 2018-2025 KMath contributors.
* Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
*/

package space.kscience.kmath.linear

import space.kscience.attributes.Attributes
import space.kscience.kmath.PerformancePitfall
import space.kscience.kmath.UnstableKMathAPI
import space.kscience.kmath.operations.Ring

/**
* Mutable sparse matrix that stores values only for non-zero cells ([DOK format](https://en.wikipedia.org/wiki/Sparse_matrix#Dictionary_of_keys_(DOK))).
*
* [SparseMatrix] is ineffective, but does not depend on particular [LinearSpace]
*
* Using this class is almost always a [PerformancePitfall]. It should be used only for special cases like manual matrix building.
*/
@UnstableKMathAPI
public class SparseMatrix<T>(
override val rowNum: Int,
override val colNum: Int,
private val zero: T,
cells: Map<Pair<Int, Int>, T> = emptyMap(),
override val attributes: Attributes = Attributes.EMPTY,
) : MutableMatrix<T> {

private val cells = cells.toMutableMap()

override fun get(i: Int, j: Int): T {
if (i !in 0 until rowNum) throw IndexOutOfBoundsException("Row index $i out of row range 0..$rowNum")
if (j !in 0 until colNum) throw IndexOutOfBoundsException("Column index $j out of column range 0..$colNum")
return cells[i to j] ?: zero
}

override fun set(i: Int, j: Int, value: T) {
require(i in 0 until rowNum) { "Row index $i is out of bounds: 0..<$rowNum" }
require(j in 0 until colNum) { "Row index $j is out of bounds: 0..<$colNum" }
val coordinates = i to j
if (cells[coordinates] != null || value != zero) {
cells[coordinates] = value
}
}

override fun set(index: IntArray, value: T) {
require(index.size == 2) { "Index array must contain two elements." }
set(index[0], index[1], value)
}
}

@UnstableKMathAPI
public fun <T> SparseMatrix<T>.fill(vararg elements: T): SparseMatrix<T> {
require(rowNum * colNum == elements.size) { "The number of elements ${elements.size} is not equal $rowNum * $colNum" }
for (i in 0 until rowNum) {
for (j in 0 until colNum) {
set(i, j, elements[i * rowNum + j])
}
}
return this
}

@UnstableKMathAPI
public fun <T> LinearSpace<T, Ring<T>>.sparse(
rows: Int,
columns: Int,
attributes: Attributes = Attributes.EMPTY,
builder: SparseMatrix<T>.() -> Unit = {}
): SparseMatrix<T> = SparseMatrix(rows, columns, elementAlgebra.zero, attributes = attributes).apply(builder)
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,6 @@
package space.kscience.kmath.linear

import space.kscience.attributes.Attributes
import space.kscience.kmath.nd.ShapeND


/**
Expand All @@ -20,13 +19,5 @@ public class VirtualMatrix<out T>(
override val attributes: Attributes = Attributes.EMPTY,
public val generator: (i: Int, j: Int) -> T,
) : Matrix<T> {

override val shape: ShapeND get() = ShapeND(rowNum, colNum)

override operator fun get(i: Int, j: Int): T = generator(i, j)
}

public fun <T : Any> MatrixBuilder<T, *>.virtual(
attributes: Attributes = Attributes.EMPTY,
generator: (i: Int, j: Int) -> T,
): VirtualMatrix<T> = VirtualMatrix(rows, columns, attributes, generator)
}
Original file line number Diff line number Diff line change
Expand Up @@ -23,10 +23,17 @@ public interface MatrixScope<T> : WithType<T>
*/
public interface MatrixAttribute<T> : StructureAttribute<T>

/**
* Matrices with this feature are symmetric, meaning `matrix[i,j] == matrix[j,i]`
*/
public interface Symmetric : MatrixAttribute<Unit>, FlagAttribute{
public companion object: Symmetric
}

/**
* Matrices with this feature are considered to have only diagonal non-zero elements.
*/
public interface IsDiagonal : MatrixAttribute<Unit>, FlagAttribute {
public interface IsDiagonal : Symmetric {
public companion object : IsDiagonal
}

Expand Down
Loading

0 comments on commit c70a7c5

Please sign in to comment.